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

(1.大連理工大學(xué) 海岸和近海工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，大連 116024；2.國(guó)家海洋環(huán)境監(jiān)測(cè)中心，大連 116023)

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

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

1 二維水動(dòng)力數(shù)值模型

1.1 控制方程

守恒型的二維水動(dòng)力數(shù)值模型控制方程如下

(1)

(2)

(3)

其中

(4)

(5)

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

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

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

1.2 底摩阻項(xiàng)修正

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

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

(6)

上式中Ωi為控制體在水平面的投影面積；△t為時(shí)間步間隔；R為控制方程整理后除了時(shí)間離散項(xiàng)和底摩阻源項(xiàng)以外其余各項(xiàng)移至等號(hào)右側(cè)后，再進(jìn)行空間積分的結(jié)果；Sf為底摩阻源項(xiàng)空間積分的計(jì)算結(jié)果。

(7)

根據(jù)兩步計(jì)算得到的流速方向關(guān)系，修正底摩阻源項(xiàng)得到第n+1時(shí)間步的流速

(8)

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

2 干濕邊界處理方法

2.1 干濕界面位置水位梯度項(xiàng)修正

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

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

在求解連續(xù)性方程的過程中，采用單元節(jié)點(diǎn)格式的有限體積法，水位變量保存在三角形網(wǎng)格的節(jié)點(diǎn)處。設(shè)水位η=D+zb(zb為底高程)，在靜水平衡問題計(jì)算過程中，必須保持水面靜止?fàn)顟B(tài)，而在干濕邊界處的水面狀態(tài)會(huì)如圖 2-a所示。

(9)

(10)

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

在數(shù)值格式中每一個(gè)控制體內(nèi)可供流出的水體體積是有限的，但計(jì)算過程中由于時(shí)間步長(zhǎng)的關(guān)系，如不做修正處理，則在干濕邊界位置處一個(gè)時(shí)間步內(nèi)部分控制單元出現(xiàn)負(fù)水深。為了保證數(shù)值格式的穩(wěn)定性，在大多數(shù)的淺水模型僅利用最小水深值對(duì)水位進(jìn)行了修正

η=max(η,zb+dmin)

(11)

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

針對(duì)質(zhì)量不守恒的問題，本文提出了針對(duì)CC-CV混合格式有限體積法的體積修正方法，網(wǎng)格間拓?fù)潢P(guān)系如圖 3所示。對(duì)節(jié)點(diǎn)控制單元的水體進(jìn)行了迭代修正，直到消除所有負(fù)水深節(jié)點(diǎn)控制單元。

首先，在求解連續(xù)性方程的過程中，檢查節(jié)點(diǎn)控制單元是否出現(xiàn)負(fù)水深節(jié)點(diǎn)。

圖3 零質(zhì)量誤差計(jì)算過程網(wǎng)格拓?fù)潢P(guān)系示意圖Fig.3 Schematic diagram of the grid topological relationship for zero mass error computation process

然后，針對(duì)出現(xiàn)負(fù)體積的節(jié)點(diǎn)，搜索與其相鄰的各個(gè)節(jié)點(diǎn)，以周圍節(jié)點(diǎn)水體質(zhì)量變化為入流的節(jié)點(diǎn)控制體的剩余水體體積作為加權(quán)系數(shù)，對(duì)于水體質(zhì)量為出流的節(jié)點(diǎn)控制單元，不參與加權(quán)計(jì)算。周圍各受體節(jié)點(diǎn)控制單元的水體調(diào)整后的入流水體體積為

(12)

最后，在新出現(xiàn)的負(fù)水深節(jié)點(diǎn)，重復(fù)以上步驟，直到無負(fù)水深節(jié)點(diǎn)出現(xiàn)。

3 模型驗(yàn)證

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

3.1 靜水平衡問題

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

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 靜水平衡問題改進(jìn)前后計(jì)算流場(chǎng)對(duì)比Fig.4 Comparison of flow filed result of still water problem

從改進(jìn)前后數(shù)值模型在模擬靜水平衡試驗(yàn)時(shí)得到的流場(chǎng)計(jì)算結(jié)果對(duì)比可以看出：在改進(jìn)前，水體在初始靜止且無任何外力驅(qū)動(dòng)條件下，由于干濕邊界節(jié)點(diǎn)處水位為max(η,zb+dmin)，引起了干濕界面處水面梯度發(fā)生變化，從而導(dǎo)致在計(jì)算域內(nèi)出現(xiàn)了不符合實(shí)際情況的擾動(dòng)；而在本文模型改進(jìn)后的計(jì)算結(jié)果中，水體始終保持靜止?fàn)顟B(tài)，不會(huì)產(chǎn)生因?yàn)閿?shù)值格式引起的錯(cuò)誤流速，有效的保證了數(shù)值格式的和諧性，能夠正確模擬部分淹沒的非平底地形問題。

3.2 非平底潰壩問題

國(guó)內(nèi)外諸多學(xué)者，如黃玉新[9]、柏祿海[10]、Song[11]和Brufau[12]等人均模擬了非平底潰壩問題，以用于檢驗(yàn)水動(dòng)力數(shù)值模型在干濕邊界位置處的處理能力和封閉水池內(nèi)的質(zhì)量守恒性。

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

(16)

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

(17)

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

5-a 本文模型

5-b 文獻(xiàn)[11]模型

6-a 本文模型

6-b MIKE 21模型

圖7 水體質(zhì)量相對(duì)誤差對(duì)比Fig.7 Comparison of relative error of water mass

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

4 結(jié)論

本文基于二維非結(jié)構(gòu)網(wǎng)格水動(dòng)力數(shù)值模型的基礎(chǔ)上提出了一種改進(jìn)的干濕界面處理方法和零質(zhì)量誤差修正方法，并將模型應(yīng)用于靜水平衡問題和非平底潰壩問題。通過計(jì)算結(jié)果對(duì)比表明：本文模型改進(jìn)的干濕邊界處理方法在靜水平衡問題中能夠得到和諧的計(jì)算結(jié)果；在非平底地形潰壩問題應(yīng)用中，通過與其他模型計(jì)算結(jié)果的對(duì)比表明，本文模型可以很好地模擬在復(fù)雜地形下的潰壩問題中水體演進(jìn)的過程，并能夠有效地保證水體質(zhì)量的守恒性。可見本文模型具有較好地?cái)?shù)值格式和諧性、穩(wěn)定性和質(zhì)量守恒性。

[1] 宋利祥. 潰壩洪水?dāng)?shù)學(xué)模型及水動(dòng)力學(xué)特性研究[D]. 武漢：華中科技大學(xué), 2012.

[2] 潘存鴻. 三角形網(wǎng)格下求解二維淺水方程的和諧Godunov格式[J]. 水科學(xué)進(jìn)展, 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] 呂彪, 金生, 艾叢芳. 非結(jié)構(gòu)化網(wǎng)格下求解二維淺水方程的和諧Roe格式[J]. 水利水運(yùn)工程學(xué)報(bào), 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] 黃玉新, 張寧川. 零質(zhì)量誤差的干濕邊界技術(shù)及其在淺水模型中的應(yīng)用[J]. 水動(dòng)力學(xué)研究與進(jìn)展: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] 柏祿海, 金生. 無質(zhì)量誤差的具有復(fù)雜地形的淺水問題[J]. 水動(dòng)力學(xué)研究與進(jìn)展, 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] 呂彪. 基于非結(jié)構(gòu)化網(wǎng)格的具有自由表面水波流動(dòng)數(shù)值模擬研究[D]. 大連:大連理工大學(xué), 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] 黃玉新. 多功能淺水模型的建立及其應(yīng)用研究[D]. 大連:大連理工大學(xué), 2013.

[10] 柏祿海. 淺水方程高分辨率算法的研究[D]. 大連:大連理工大學(xué), 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.