基于CC-CV混合有限體積法的干濕邊界處理技術及其在淺水問題中的應用

趙旭東，孫家文,2，孫昭晨*，梁書秀

(1.大連理工大學 海岸和近海工程國家重點實驗室，大連 116024；2.國家海洋環境監測中心，大連 116023)

近年來，海岸水動力數值模型在解決洪水演進、潮波漫灘等近岸淺水問題中得到了廣泛的應用，此類問題難點在于復雜干濕邊界條件處理。其中，干濕邊界位置底摩阻源項處理方法對數值格式穩定性和計算結果的準確性具有較大的影響，宋利祥[1]提出的一種半隱式處理方法，可以有效的保證數值格式的穩定性；潘存鴻[2]等人提出了水位床底方程法，呂彪[3]等人提出了特殊的底坡源項處理技術，在一定程度上克服了靜水平衡問題求解時出現的不和諧現象。但求解復雜地形條件下水體流動問題時，如果僅對控制方程中的源項進行簡單處理，容易導致數值計算過程中的質量不守恒現象，尤其是針對封閉水域及淺水問題時，對計算結果的影響更加明顯，黃玉新[4]和柏祿海[5]等人分別提出了不同的零質量誤差邊界處理方法。

數值格式和諧性又稱為C特性，在計算復雜地形下靜水問題時，計算域內各個位置處水位靜止，流速始終為零的性質。本文結合淺水模型中干濕邊界處理技術最新研究成果，對其加以改進，以滿足數值格式的和諧性，通過修改負水深節點及其相鄰節點的流出水體體積消除計算過程中出現的質量不守恒問題。將本文模型應用于2個典型算例，驗證了數值格式在解決靜水平衡問題和非平底潰壩問題中的格式穩定性、和諧性和質量守恒性。

1 二維水動力數值模型

1.1 控制方程

守恒型的二維水動力數值模型控制方程如下

(1)

(2)

(3)

其中

(4)

(5)

控制方程中η為水位；D為水深；u,v分別為x,y方向的垂向平均流速；ρ0為水體密度；g為重力加速度；AM為水平渦粘系數；(τsx，τsy)和(τbx,τby)分別為x和y方向的表面切應力和底部切應力。

本文模型在水平空間將計算域剖分為三角形網格，采用單元中心(CC)和單元節點(CV)混合格式有限體積法進行空間離散(圖1)，矢量數據保存在控制體網格單元中心，水位等矢量數據保存在網格節點位置處。

圖1 三角形網格空間離散示意圖Fig.1 Schematic diagram of triangular grid spatial discretization

1.2 底摩阻項修正

在地形變化復雜、起伏較大的天然水域，容易出現水位或流速突變的情況，此時底摩阻項對數值模型的穩定性、計算效率具有較大的影響[6]。尤其是在干濕界面處，水深趨近于零，如果仍然采用曼寧系數直接計算底摩阻，由于底摩阻計算公式中水深為分母項，由此會產生不符合實際情況的摩阻值，從而引起數值格式的不穩定。針對這一問題，常用的處理方法是當水深小于某個值后，將底摩阻給定為常數，但這種處理方法對數值的計算精度影響較大，Liang等人[7]提出分離隱式方法處理底摩阻項，通過分裂算子法求解，其處理方法在一定程度上解決底摩阻項處理不當引起的數值格式失穩問題，但計算過程相對復雜。本文提出了一個相對簡化的處理方法，從底摩阻項作為阻力項的特性出發，對其進行修正，基本過程如下。

第一步：考慮底摩阻作用情況下求解二維動量方程(2)和(3)，計算流速的中間變量(u*,v*)

(6)

上式中Ωi為控制體在水平面的投影面積；△t為時間步間隔；R為控制方程整理后除了時間離散項和底摩阻源項以外其余各項移至等號右側后，再進行空間積分的結果；Sf為底摩阻源項空間積分的計算結果。

(7)

根據兩步計算得到的流速方向關系，修正底摩阻源項得到第n+1時間步的流速

(8)

如上式所示，在考慮底摩阻作用情況下，如果n與n+1時刻的流速方向相同，則說明底摩阻在計算過程中，始終起到阻力作用，沒有導致流向發生改變；在考慮底摩阻作用情況下，如果n與n+1時刻的流速方向相反，而在不考慮底摩阻作用情況下，如果n與n+1時刻的流速方向相同，流向發生變化的主要原因是由于底摩阻作用引起的，這一現象不符合底摩阻作為阻力的物理性質，為此，可以認為在這一時間步計算過程中，在底摩阻的作用下，流速為零，不再導致流向發生變化；在不考慮底摩阻作用情況下，如果n與n+1時刻的流速方向相反，說明在不考慮底摩阻所用，依然會導致流向發生變化，而在計算過程中，底摩阻的作用方向發生了變化，為此在這一計算過程中可以不考慮底摩阻的作用。

2 干濕邊界處理方法

2.1 干濕界面位置水位梯度項修正

在非結構化網格模型中，通常有兩種底高程存儲方式：(1)將底高程變量保存在網格形心位置處，在控制單元內底部高程一致；(2)將底高程數據保存節點位置處，在控制單元內呈線性分布。第一種方法僅為一階精度，而第二種方法則具有二階精度[8]。為了能夠更加準確的模擬復雜地形，提高計算精度和數值格式穩定性，本文模型采用第二種方案，利用存儲在網格節點的底高程數據，構建連續變化的地形，即斜底三角形網格模型。

圖2 水位梯度修正示意圖Fig.2 Schematic diagram of surface gradient modification

在求解連續性方程的過程中，采用單元節點格式的有限體積法，水位變量保存在三角形網格的節點處。設水位η=D+zb(zb為底高程)，在靜水平衡問題計算過程中，必須保持水面靜止狀態，而在干濕邊界處的水面狀態會如圖 2-a所示。

(9)

(10)

2.2 基于CC-CV混合有限體積格式的物質守恒干濕邊界捕捉方法

在數值格式中每一個控制體內可供流出的水體體積是有限的，但計算過程中由于時間步長的關系，如不做修正處理，則在干濕邊界位置處一個時間步內部分控制單元出現負水深。為了保證數值格式的穩定性，在大多數的淺水模型僅利用最小水深值對水位進行了修正

η=max(η,zb+dmin)

(11)

該方法在修正水位的過程中，無形中增加了計算域內總的水體體積，進而破壞了數值模型質量守恒性。針對這一問題，通常采用兩種方法處理，一是根據計算域內所有控制單元內剩余水體體積和流出水體速度，利用自適應方法計算最大時間步長以保證節點控制單元剛好干出，該方法雖然可以保證水體質量的守恒，但是由于需要實時計算時間步長，且計算得到的時間步長相對較小，嚴重影響了數值模型的計算效率；二是黃玉新[4]等人提出了水深修正法，針對在求解過程中出現負水深(Di)的單元，令這些單元的水深和速度為零，相當于向網格單元輸入了額外的水體體積△VA=|DAΩA|，后續計算需要減少周圍網格單元體積，保證水體質量的守恒，然而進行水深修正后容易出現新的負水深單元，黃玉新[4]等人基于CC格式的有限體積法，利用迭代修正的方法避免了這一問題的出現。

針對質量不守恒的問題，本文提出了針對CC-CV混合格式有限體積法的體積修正方法，網格間拓撲關系如圖 3所示。對節點控制單元的水體進行了迭代修正，直到消除所有負水深節點控制單元。

首先，在求解連續性方程的過程中，檢查節點控制單元是否出現負水深節點。

圖3 零質量誤差計算過程網格拓撲關系示意圖Fig.3 Schematic diagram of the grid topological relationship for zero mass error computation process

然后，針對出現負體積的節點，搜索與其相鄰的各個節點，以周圍節點水體質量變化為入流的節點控制體的剩余水體體積作為加權系數，對于水體質量為出流的節點控制單元，不參與加權計算。周圍各受體節點控制單元的水體調整后的入流水體體積為

(12)

最后，在新出現的負水深節點，重復以上步驟，直到無負水深節點出現。

3 模型驗證

基于上文中提出的干濕邊界處理方法和底摩阻項處理方法，建立了能夠保證數值格式和諧性、質量守恒性和數值格式穩定的二維水動力模型，并將模型應用于求解靜水平衡問題和非平底潰壩問題中。為了更直觀有效的證明本文提出的處理方法在數值計算過程中的優越性，在本文模型中將計算結果轉化為MIKE21的結果文件格式，并利用其后處理工具實現了結果的可視化，分別對比了改進前后的流場分布和水體質量誤差計算結果。

3.1 靜水平衡問題

在靜水平衡問題中為長1 m、寬1 m的矩形封閉水池，水池四周為封閉的固邊界，邊界法向流速為0 m/s；底邊界Manning系數n=0.001。水池底高程zb通過式(13)計算獲得

zb(x,y)=max{0,0.25-5[(x-0.5)2]+(y-0.5)2]}

(13)

初始條件

η=d+zb=0.2 m

(14)

u=v=0 m/s

(15)

圖4 靜水平衡問題改進前后計算流場對比Fig.4 Comparison of flow filed result of still water problem

從改進前后數值模型在模擬靜水平衡試驗時得到的流場計算結果對比可以看出：在改進前，水體在初始靜止且無任何外力驅動條件下，由于干濕邊界節點處水位為max(η,zb+dmin)，引起了干濕界面處水面梯度發生變化，從而導致在計算域內出現了不符合實際情況的擾動；而在本文模型改進后的計算結果中，水體始終保持靜止狀態，不會產生因為數值格式引起的錯誤流速，有效的保證了數值格式的和諧性，能夠正確模擬部分淹沒的非平底地形問題。

3.2 非平底潰壩問題

國內外諸多學者，如黃玉新[9]、柏祿海[10]、Song[11]和Brufau[12]等人均模擬了非平底潰壩問題，以用于檢驗水動力數值模型在干濕邊界位置處的處理能力和封閉水池內的質量守恒性。

在該算例中，模型計算域為70 m×30 m的矩形水池，四周為固壁邊界，在距離左側邊界30 m和47.5 m位置處分別有3個凸起的擋水建筑物，計算域內的底高程通過式(16)計算獲得

(16)

初始條件：流速u=v=0 m/s，初始水位為如式(17)所示，水體總體積為900 m3。

(17)

為了更好地驗證本文模型在解決非平底潰壩問題時的優越性，分別于不同的數值模型計算結果進行了對比，圖 5為本文模型和Song等改進的淺水模型在不同時刻時模擬得到的自由水面分布。從計算結果對比可以看出2個模型模擬獲得的水體傳播過程基本一致：在t=6 s時刻，水體已經淹沒左側的兩個小的凸起建筑物，水體前沿在大的凸起建筑物的阻礙作用下左側出現增水現象；t=12 s時刻，水體繞過大的凸起建筑物后繼續向下游傳播，傳播過程中大的凸起建筑物始終有大部分區域未被淹沒；在t=30 s時刻，水體到達右側固邊界，開始向回傳播；在t=300 s時刻，在各個阻力項的共同作用下，計算域內水體基本達到靜止狀態。由此可見本文模型能夠正確反映非平底潰壩問題中水體演進的過程。

5-a 本文模型

5-b 文獻[11]模型

6-a 本文模型

6-b MIKE 21模型

圖7 水體質量相對誤差對比Fig.7 Comparison of relative error of water mass

圖6為本文模型和MIKE 21的流速計算結果對比，從非平底潰壩實驗的地形設計可以看出，原始地形和水位呈對稱布置，因此計算結果也應具有對稱性，通過對比不同時刻本文模型與MIKE的流場計算結果對比可以看出，本文模型的流場對稱性更好，更加符合實際流場情況；在圖 7的對比中可以看出，本文模型在質量守恒性方面明顯優于MIKE21的計算結果，本文模型計算過程中產生的質量誤差約為MIKE21計算誤差的60%左右，能夠較好保證數值計算過程中質量守恒性。

4 結論

本文基于二維非結構網格水動力數值模型的基礎上提出了一種改進的干濕界面處理方法和零質量誤差修正方法，并將模型應用于靜水平衡問題和非平底潰壩問題。通過計算結果對比表明：本文模型改進的干濕邊界處理方法在靜水平衡問題中能夠得到和諧的計算結果；在非平底地形潰壩問題應用中，通過與其他模型計算結果的對比表明，本文模型可以很好地模擬在復雜地形下的潰壩問題中水體演進的過程，并能夠有效地保證水體質量的守恒性。可見本文模型具有較好地數值格式和諧性、穩定性和質量守恒性。

[1] 宋利祥. 潰壩洪水數學模型及水動力學特性研究[D]. 武漢：華中科技大學, 2012.

[2] 潘存鴻. 三角形網格下求解二維淺水方程的和諧Godunov格式[J]. 水科學進展, 2007(2): 204-209.

PAN C H. Well-balanced Godunov-type scheme for 2D shallow water flow with triangle mesh [J]. ADVANCES IN WATER SCIENCE, 2007(2): 204-209.

[3] 呂彪, 金生, 艾叢芳. 非結構化網格下求解二維淺水方程的和諧Roe格式[J]. 水利水運工程學報, 2010(2): 39-44.

LV B, JIN S, AI C F.Well-balanced roe-type scheme for 2D shallow water flow using unstructured grids[J]. HYDRO-SCIENCE AND ENGINEERING, 2010(2): 39-44.

[4] 黃玉新, 張寧川. 零質量誤差的干濕邊界技術及其在淺水模型中的應用[J]. 水動力學研究與進展:A輯, 2013(6): 745-753.

HUANG Y X, ZHANG N C. An efficient wetting-drying front technique with zero mass error and its application in shallow water model [J]. CHINESE JOURNAL OF HYDRODYNAMIC, 2013(6): 745-753.

[5] 柏祿海, 金生. 無質量誤差的具有復雜地形的淺水問題[J]. 水動力學研究與進展, 2008, 23(5): 515-521.

BAI L H, JIN S. Numerical simulation of zero mass error using wetting-drying conditions in shallow flow over complex topography[J]. CHINESE JOURNAL OF HYDRODYNAMIC, 2008, 23(5): 515-521.

[6] 呂彪. 基于非結構化網格的具有自由表面水波流動數值模擬研究[D]. 大連:大連理工大學, 2010.

[7] Liang Q, Marche F. Numerical resolution of well-balanced shallow water equations with complex source terms[J]. Advances in water resources, 2009, 32(6): 873-884.

[8] Begnudelli L, Sanders B F. Unstructured grid finite-volume algorithm for shallow-water flow and scalar transport with wetting and drying[J]. Journal of hydraulic engineering, 2006, 132(4): 371-384.

[9] 黃玉新. 多功能淺水模型的建立及其應用研究[D]. 大連:大連理工大學, 2013.

[10] 柏祿海. 淺水方程高分辨率算法的研究[D]. 大連:大連理工大學, 2013.

[11] Song L, Zhou J, Li Q, et al. An unstructured finite volume model for dam-break floods with wet/dry fronts over complex topography[J]. International Journal for Numerical Methods in Fluids, 2011, 67(8): 960-980.

[12] Brufau P, Vázquez-Cendón M E, García-Navarro P. A numerical model for the flooding and drying of irregular domains[J]. International Journal for Numerical Methods in Fluids, 2002, 39(3): 247-275.