李紹武，張弛，楊學(xué)斌，李文善

（天津大學(xué)水利工程仿真與安全國家重點實驗室，天津300072）

基于有限體積法的平面二維水流數(shù)學(xué)模型的改進(jìn)

李紹武，張弛，楊學(xué)斌，李文善

（天津大學(xué)水利工程仿真與安全國家重點實驗室，天津300072）

在基于有限體積法的近岸水流數(shù)學(xué)模型中，單元交界面數(shù)值通量的計算精度、潛堤淹沒和出露交替的模擬以及動邊界技術(shù)都是至關(guān)重要的。首先通過加入數(shù)值通量底坡補(bǔ)償項對Osher格式的數(shù)值通量計算格式進(jìn)行了改進(jìn)；又提出單元淹沒度的概念，提高了模型在模擬潛堤出露和潮間帶干濕變化中的適應(yīng)性。利用典型算例對模型進(jìn)行了驗證，得到了滿意的結(jié)果。

水流數(shù)學(xué)模型;底坡補(bǔ)償項;有限體積法;潛堤;淹沒度;動邊界

基于有限體積法的水流數(shù)學(xué)模型的精度在很大程度上取決于單元交界面數(shù)值通量計算格式的精度。有限體積法計算格式大都將主要變量定義在單元中心，由于變量定義點的錯位，難免造成界面數(shù)值通量計算產(chǎn)生一定數(shù)值誤差，這種誤差在時間步進(jìn)中不斷積累，將影響計算精度和格式穩(wěn)定性。另一方面，將水流數(shù)學(xué)模型運(yùn)用于潮汐河口或緩坡海域時，需要對活動邊界進(jìn)行處理。目前常用的干濕網(wǎng)格法、水邊線步進(jìn)法、窄逢法和井點法等［1-2］都有一定局限性，如干濕網(wǎng)格法和水邊線步進(jìn)法中邊界的干濕變化以一個網(wǎng)格為單位更替，精度取決于網(wǎng)格大??；窄逢法和井點法雖然精度較高，但數(shù)值穩(wěn)定性稍差。數(shù)學(xué)模型在近岸應(yīng)用中的另一個問題是潛堤出露給建模帶來的困難。結(jié)合有限體積數(shù)值求解方法［3-5］，首先提出一種交界面數(shù)值通量計算精度的改進(jìn)方法，導(dǎo)出了相關(guān)計算公式；然后，針對活動邊界和潛堤出露問題，提出了淹沒度的概念，即將所有單元（包括內(nèi)單元、水陸交界單元和陸單元）都視為具有不同淹沒度的網(wǎng)格單元?？刂品匠淘谒袉卧线M(jìn)行統(tǒng)一離散，從而使變量在水、陸單元之間可以連續(xù)過渡，計算更為簡便，精確度有一定提高。

1 平面二維淺水波數(shù)學(xué)模型的控制方程

1.1 控制方程

近岸潮波運(yùn)動二維淺水波方程包括質(zhì)量守恒方程和動量守恒方程。

質(zhì)量守恒方程

x方向動量守恒方程

y方向動量守恒方程

式中：H為全水深；U、V分別為垂向平均流速的x、y方向分量；g為重力加速度；Zb為底高程；C為謝才系數(shù)，可由曼寧公式計算，其中n為曼寧糙率系數(shù)；f為科氏參數(shù)，f=2Ωsin?，Ω為地球的自轉(zhuǎn)頻率，Ω=2π/86 164，?為當(dāng)?shù)鼐暥?；νh為水平方向上的紊動粘滯系數(shù)，由Smagorinsky亞格子紊動模型得到，即

式中：cs=0.1～0.2；Δs為離散單元面積。

1.2 初始條件

初始水位按η=η0(x,y)給定，初始速度按0給定。

1.3 邊界條件

固邊界采用“不穿透”邊界條件，即V×n=0，其中V為流速矢量，n為沿固邊界外法線方向的單位矢量。開邊界處給定水位過程η=ηb(x,y,t)。

2 控制方程的數(shù)值求解方法

2.1 有限體積法計算格式

根據(jù)有限體積法的求解思路［5］，將控制方程寫為向量形式

式中：U=H(1,U,V)T為守恒物理量向量，F(xiàn)=(HU,HU2+gH2/2,HUV)T為x方向上的通量，G=(HU,HUV,HV2+ gH2/2)T為y方向上的通量為右端項，其中的Fb、Ff、Fc和Ft分別代表底坡項、底摩阻項、科氏力項和紊動摻混項。

對式（5）在單元體上進(jìn)行積分，并應(yīng)用格林公式，得到

式中：U?為轉(zhuǎn)換到交界面局部坐標(biāo)下的狀態(tài)變量。

定義單元第j條邊的邊長為Lj，U?和V?分別為局部坐標(biāo)系下狀態(tài)變量的分量，h（l）為單元邊上的水深變化函數(shù)，則單元數(shù)值通量計算公式為

假設(shè)h（l）沿單元邊線性變化，且單元邊兩端點的底高程差為ΔZj，hˉsj為第j條邊的平均水深，則可以導(dǎo)出數(shù)值通量計算公式為

式（6）采用時間向前差分即可進(jìn)行顯式求解。

圖1 有效單元示意圖Fig.1 Sketch of an effective element

2.2 淹沒度概念

假設(shè)任意多邊形單元面積為A，被水淹沒的部分面積為As（圖1），定義淹沒度αs=As/A［6］，該量表示單元中淹沒部分與全部面積的比值，可以表征單元的干濕狀態(tài)。

顯然，當(dāng)αs=0時，網(wǎng)格處于全干狀態(tài)；當(dāng)αs=1時，As=A，網(wǎng)格處于全淹狀態(tài)；當(dāng)0＜αs＜1時，單元處于半干半淹狀態(tài)。

對于任意淹沒度（0＜αs≤1）的單元，單元水體體積為

式中：Vw為單元中被水淹沒部分的水體積。

假設(shè)某單元3個節(jié)點的底高程從小到大依次為Z1，Z2，Z3，根據(jù)三角單元的淹沒情況，可能出現(xiàn)以下4種情形。

情形一：當(dāng)η≤Z1時，αs=0，Hs=0。

情形二：當(dāng)Z1≤η≤Z2時，單元僅一個節(jié)點被淹沒，淹沒部分仍然為三角形，淹沒度此時，平均水深

情形三：當(dāng)Z2≤η≤Z3時，淹沒度平均水深為

情形四：當(dāng)η≥Z3時

圖2 計算區(qū)域水位及地形縱剖面圖Fig.2Longitudinal profile of computational domain

3 水流模型的驗證

以下通過具體算例討論模型的適用性。

3.1 恒定水位算例

為了考察通量補(bǔ)償項的影響效果，構(gòu)造尺寸為500 m－250 m的非平底模型（圖2和圖3），底部糙率為0.02，計算區(qū)域上下邊界為固邊界，左右為開邊界，水位保持10 m不變。計算穩(wěn)定后，沿x軸方向剖面水位變化如圖4所示，結(jié)果表明，缺少補(bǔ)償項將會在底坡度變化處引起顯著誤差，且誤差有振蕩，最大達(dá)到0.003 m，而在兩側(cè)的平底處誤差為0，而這正是引入底坡補(bǔ)償項的效果。

圖3 計算區(qū)域地形平面及網(wǎng)格剖分圖Fig.3Topography and gridding of computational domain

圖4 有無補(bǔ)償項計算結(jié)果對比Fig.4Comparison of numerical results of water surface elevation for the cases with and without ompensational term

3.2 均勻底坡上的圓形淺灘算例

計算區(qū)域為一5000m×2500m的矩形斜底水槽，底部糙率為0.02，中心設(shè)一球冠，半徑為673 m，其底高程由以下方程給出（圖5）。

模型上下兩側(cè)設(shè)置為固邊界，左側(cè)為岸邊界，右側(cè)為開邊界。開邊界按正弦波給出水位過程，平均水位為0 m，振幅為3 m，周期為12 h25 min。

圖6給出漲落潮各特征時刻的流速和淹沒度計算結(jié)果。低潮時球冠部分露灘，高潮時淺灘完全淹沒。淹沒度計算結(jié)果可以較好反映出球冠及潮灘的淹沒和露灘過程。

3.4 導(dǎo)堤算例

在實際工程中，常常會遇到導(dǎo)堤在高潮時淹沒，而低潮時出露的情況，這種交替變化給建模帶來困難。此外，由于導(dǎo)堤外形狹長，若采用加密網(wǎng)格的方法勢必導(dǎo)致局部網(wǎng)格尺度極小，網(wǎng)格數(shù)量激增。同時，根據(jù)克朗條件，計算時間步長也需相應(yīng)減小，從而導(dǎo)致計算時間劇增。

圖5 均勻底坡上圓形淺灘地形Fig.5Topography of circular shoal on sloping bottom

圖6 圓形淺灘流速及淹沒度計算結(jié)果Fig.6Computational results of velocity and submergence on sloping bottom with a circular shoal

考慮到導(dǎo)堤形狀窄長的特點，若將導(dǎo)堤軸線所在位置確定為網(wǎng)格邊線來進(jìn)行網(wǎng)格剖分，對與導(dǎo)堤軸線相連的單元通過引入淹沒度的方法，來反映其淹沒和出露的特點，則可以避免在導(dǎo)堤所在位置進(jìn)行過分加密的問題。

構(gòu)造5000m×2500m的矩形平底水槽（圖7），底部糙率為0.02，在區(qū)域中設(shè)置一頂高程為5 m的弧形堤，其軸線所在的圓弧半徑為1350m，暫不考慮科氏力作用。開邊界同樣按正弦潮波給定。網(wǎng)格剖分時保證網(wǎng)格邊線通過潛堤。

圖8為漲落潮過程中各特征時刻的流速計算結(jié)果?？梢钥闯?，導(dǎo)堤淹沒后，對流場仍有較大影響，這與實際是相符的。由于潛堤的影響，堤頭處流速變大，潛堤彎曲形狀的使得漲潮時堤頭流速明顯增大。圖9為水位計算結(jié)果，圖10為淹沒度計算結(jié)果，可以看出，隨著漲落潮的交替，導(dǎo)堤所在位置淹沒度也隨之變化。

圖7 潛堤模型地形及網(wǎng)格剖分Fig.7Topography and meshes of the submerged dike model

圖8 流速計算結(jié)果Fig.8Computational results of velocity

圖9 水位計算結(jié)果Fig.9Computational results of water surface elevation

圖10 淹沒度計算結(jié)果Fig.10Computational results of submergence index

4 結(jié)語

對基于有限體積法的數(shù)值通量計算方法進(jìn)行了改進(jìn)，導(dǎo)出了底坡影響下的數(shù)值通量計算格式，使得該數(shù)值模型對較大底坡的地形具有更強(qiáng)的適應(yīng)性。提出了淹沒度概念，并給出了任意網(wǎng)格淹沒度的有限體積離散方法，使得單元可以在全干全濕兩種狀態(tài)下進(jìn)行平滑過度，實現(xiàn)了水陸邊界的連續(xù)性推移的過程，且計算耗時無明顯增加。典型算例計算結(jié)果表明，模型在處理底坡變化劇烈的淺灘、潛堤等復(fù)雜地形時具有良好的性能。

［1］曹祖德，王運(yùn)洪.水動力學(xué)泥沙數(shù)值模擬［M］.天津：天津大學(xué)出版社，1994.

［2］陶建華.水波的數(shù)值模擬［M］.天津：天津大學(xué)出版社，2005.

［3］譚維炎，胡四一.淺水流動計算中一階有限體積法Osher格式的實現(xiàn)［J］.水科學(xué)進(jìn)展，1994，5（4）：262-270. TAN W Y，HU S Y.Implementation of First?order Finite?volume Osher Scheme in Shallow?water Flow Computation［J］.Advances in Water Science，1994，5（4）：262-270.

［4］胡四一，譚維炎.無結(jié)構(gòu)網(wǎng)格上二維淺水流動的數(shù)值模擬［J］.水科學(xué)進(jìn)展，1995，6（1）：1-9. HU S Y，TAN W Y.Numerical Modelling of Two?Dimensional Shallow Water Flows on Unstructured Grids［J］.Advances in Water Science，1995，6（1）：1-9.

［5］李紹武，盧麗鋒.河口準(zhǔn)三維涌潮數(shù)學(xué)模型研究［J］.水動力學(xué)研究與進(jìn)展，2004，19（4）：407-415. LI S W，LU L F.A Quasi?3D Numerical Model of Estuarine Tidal Bore［J］.Journal of Hydrodynamics，2004，19（4）：407-415.

［6］楊學(xué)斌.波流聯(lián)合作用下二維水流泥沙數(shù)學(xué)模型研究［D］.天津：天津大學(xué)，2008.

Improvements on a plane 2D numerical flow model based on FVM

LI Shao?wu,ZHANG Chi,YANG Xue?bin,LI Wen?shan
（State Key Laboratory of Hydraulic Engineering Simulation and Safety,Tianjin University,Tianjin 300072, China）

Precision of the numerical flux through mesh interface,treatment for dike emerging?submerging al?ternating and techniques of movable boundary have important significance for the numerical current models based on FVM.Compensational term of bottom slope in calculating numerical flux was added into the model in order to im?prove the interface numerical flux that was calculated by the Osher scheme.A concept of submergence index was proposed to enhance the adaptability of the model in simulation of dike emerging?submerging and movable bound?ary in the intertidal zone.A series of ideal computational examples were simulated to verify the properties of the model.Finally,satisfactory results were obtained.

numerical current model;compensation term of bottom slope;finite volume method;submerged dike;submergence index;movable boundary

TV 143；O 242.1

A

1005-8443（2014）05-0475-06

中俄將合作在俄羅斯共建大型海港

2013-10-22；

2013-12-12

國家自然科學(xué)基金（51379143）

李紹武（1962-），男，山東省人，教授，主要從事海岸動力學(xué)及海岸工程研究工作。

Biography：LI Shao?wu（1962-）,male,professor.

據(jù)俄羅斯最大港口運(yùn)營商“俄羅斯蘇瑪集團(tuán)”在其2014年9月舉行的莫斯科專場推介會上披露，俄中兩國企業(yè)將合作建設(shè)年吞吐量可達(dá)到6 000萬t的扎魯比諾大型萬能海港。該港位于俄遠(yuǎn)東濱海邊疆區(qū)東南部，距離中國邊境18 km，是俄遠(yuǎn)東地區(qū)的天然不凍港，有鐵路、公路與俄內(nèi)陸和中國吉林省琿春市相連。據(jù)悉，早在2014年5月舉行的上海亞信峰會上，吉林省與蘇瑪集團(tuán)簽訂了合作建設(shè)扎魯比諾萬能海港的框架協(xié)議。（殷缶，梅深）