李紹武，張弛，楊學斌，李文善

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

基于有限體積法的平面二維水流數學模型的改進

李紹武，張弛，楊學斌，李文善

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

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

水流數學模型;底坡補償項;有限體積法;潛堤;淹沒度;動邊界

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

1 平面二維淺水波數學模型的控制方程

1.1 控制方程

近岸潮波運動二維淺水波方程包括質量守恒方程和動量守恒方程。

質量守恒方程

x方向動量守恒方程

y方向動量守恒方程

式中：H為全水深；U、V分別為垂向平均流速的x、y方向分量；g為重力加速度；Zb為底高程；C為謝才系數，可由曼寧公式計算，其中n為曼寧糙率系數；f為科氏參數，f=2Ωsin?，Ω為地球的自轉頻率，Ω=2π/86 164，?為當地緯度；νh為水平方向上的紊動粘滯系數，由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 控制方程的數值求解方法

2.1 有限體積法計算格式

根據有限體積法的求解思路［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）在單元體上進行積分，并應用格林公式，得到

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

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

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

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

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

2.2 淹沒度概念

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

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

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

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

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

情形一：當η≤Z1時，αs=0，Hs=0。

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

情形三：當Z2≤η≤Z3時，淹沒度平均水深為

情形四：當η≥Z3時

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

3 水流模型的驗證

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

3.1 恒定水位算例

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

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

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

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

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

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

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

3.4 導堤算例

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

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

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

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

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

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

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

圖8 流速計算結果Fig.8Computational results of velocity

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

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

4 結語

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

［1］曹祖德，王運洪.水動力學泥沙數值模擬［M］.天津：天津大學出版社，1994.

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

［3］譚維炎，胡四一.淺水流動計算中一階有限體積法Osher格式的實現(xiàn)［J］.水科學進展，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］胡四一，譚維炎.無結構網格上二維淺水流動的數值模擬［J］.水科學進展，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］李紹武，盧麗鋒.河口準三維涌潮數學模型研究［J］.水動力學研究與進展，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］楊學斌.波流聯(lián)合作用下二維水流泥沙數學模型研究［D］.天津：天津大學，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

國家自然科學基金（51379143）

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

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

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