基于WENO格式的高精度高分辨臺(tái)風(fēng)風(fēng)暴潮數(shù)值模式

王如云，汪 天，吳楚敏，羌丹丹，周 鈞，張 鑫，張 彬，占飛

基于WENO格式的高精度高分辨臺(tái)風(fēng)風(fēng)暴潮數(shù)值模式

王如云1，汪 天2，吳楚敏1，羌丹丹3，周 鈞4，張 鑫1，張 彬1，占飛1

（1.河海大學(xué)海洋學(xué)院，江蘇南京210098；2.河海大學(xué)港航學(xué)院，江蘇南京210098；3.南通市生產(chǎn)力促進(jìn)中心，江蘇南通226019；4.河海大學(xué)水文水資源學(xué)院，江蘇南京210098)

考慮到臺(tái)風(fēng)風(fēng)暴潮在近岸淺水地區(qū)的非線性效應(yīng)，基于無(wú)結(jié)構(gòu)網(wǎng)格，通過(guò)采用有限體積法和高精度高分辨率的WENO數(shù)值格式對(duì)二維淺水方程進(jìn)行空間離散，并利用三階的Runge-Kutta格式進(jìn)行時(shí)間離散，最后利用Rogers方法解決復(fù)雜海底地形造成的通量梯度項(xiàng)與源項(xiàng)數(shù)值離散后的不平衡問(wèn)題，從而建立了二維臺(tái)風(fēng)風(fēng)暴潮數(shù)值模式。模式中的風(fēng)場(chǎng)和氣壓場(chǎng)分別采用宮崎正衛(wèi)風(fēng)場(chǎng)模式和藤田氣壓場(chǎng)模式。最后通過(guò)對(duì)江蘇沿海的風(fēng)暴增水的模擬和驗(yàn)證，表明了該數(shù)值模式對(duì)臺(tái)風(fēng)風(fēng)暴潮模擬的有效性和可行性。

WENO格式；無(wú)結(jié)構(gòu)網(wǎng)格；臺(tái)風(fēng)風(fēng)暴潮；數(shù)值模式；高精度高分辨

1 引言

臺(tái)風(fēng)風(fēng)暴潮是由臺(tái)風(fēng)引起的海水水位異常升降的現(xiàn)象[1]，對(duì)其進(jìn)行準(zhǔn)確、快速數(shù)值模擬方法方面的研究，一直是人們關(guān)注的問(wèn)題。從數(shù)值計(jì)算的網(wǎng)格來(lái)講，為了適應(yīng)復(fù)雜岸邊界形態(tài)，人們從最初采用結(jié)構(gòu)網(wǎng)格[2]，逐漸采用嵌套網(wǎng)格[3]、曲線坐標(biāo)網(wǎng)格[4]，目前為止，不少研究工作基于適應(yīng)性更好、便于整個(gè)流場(chǎng)設(shè)計(jì)統(tǒng)一計(jì)算格式而采用無(wú)結(jié)構(gòu)網(wǎng)格[5]。從數(shù)值計(jì)算的方程離散方法來(lái)講，主要是有限差分法[2]、有限元法[6]、有限體積法[5]。另外，在已有的水動(dòng)力學(xué)數(shù)值模型基礎(chǔ)上，加上臺(tái)風(fēng)場(chǎng)開(kāi)展風(fēng)暴潮數(shù)值模擬的工作也有不少成果，如李云川等[7]運(yùn)用非靜力MM5的二重網(wǎng)格雙向嵌套數(shù)值模式，模擬了2003年渤海灣的風(fēng)暴潮過(guò)程。Choi等[8]通過(guò)耦合一個(gè)三代的波浪模型、WAM模型和二維的風(fēng)暴潮模型，建立了耦合潮汐、風(fēng)暴潮和波浪的模式，并將其應(yīng)用到黃海海域中去。Kim等[9]將沿深度積分的非線性淺水波方程和SWAM模型相結(jié)合，用2006年的艾云尼臺(tái)風(fēng)作為算例，分析了風(fēng)暴潮和風(fēng)浪對(duì)天文潮的影響。Xie等[10]運(yùn)用三維風(fēng)暴潮模式，針對(duì)5種不同類型的數(shù)值實(shí)驗(yàn)來(lái)研究不對(duì)稱風(fēng)場(chǎng)對(duì)風(fēng)暴潮和漫灘的影響。

臺(tái)風(fēng)風(fēng)暴潮在近岸淺水地區(qū)的非線性效應(yīng)，可能會(huì)產(chǎn)生陡度較大的波面，甚至出現(xiàn)間斷波現(xiàn)象，這些問(wèn)題是上述數(shù)值模式進(jìn)一步發(fā)展的瓶頸。為解決這些問(wèn)題，2006年，Wang等[5]基于無(wú)結(jié)構(gòu)網(wǎng)格、有限體積法、Roe平均法，嘗試建立了高分辨率的風(fēng)暴潮二維數(shù)值預(yù)報(bào)模式，并對(duì)珠江口的風(fēng)暴潮進(jìn)行了模擬。考慮到Roe格式的參向量、間斷解展開(kāi)、線性化替代矩陣等，都是在最簡(jiǎn)單的左右常數(shù)態(tài)的假定下進(jìn)行的，這種左右常數(shù)態(tài)的假定，從根本上限制了方法的精度。1994年，Liu等[11]提出并構(gòu)造了一維空間上的三階有限體積加權(quán)本質(zhì)無(wú)震蕩(WENO)格式，它是由網(wǎng)格平均形式的ENO格式發(fā)展起來(lái)的。其后，Jiang等[12]在1996年構(gòu)造了多維空間上的三階和五階有限差分WENO格式，提出了現(xiàn)在被廣泛采用的光滑因子和非線性權(quán)[13-14]。本文基于無(wú)結(jié)構(gòu)網(wǎng)格，并利用WENO格式設(shè)計(jì)思想，設(shè)計(jì)建立了可模擬間斷波現(xiàn)象的高精度高分辨率的二維臺(tái)風(fēng)風(fēng)暴潮數(shù)值模式。運(yùn)用該數(shù)值模式對(duì)江蘇沿海地區(qū)臺(tái)風(fēng)風(fēng)暴潮引起的純風(fēng)暴增水進(jìn)行數(shù)值模擬，并對(duì)7910號(hào)臺(tái)風(fēng)的增水過(guò)程進(jìn)行驗(yàn)證，取得了比較滿意的結(jié)果。

2 基本方程和初邊界條件

2.1 二維淺水風(fēng)暴潮波控制方程

假設(shè)水體是均質(zhì)不可壓的，并且水流紊動(dòng)應(yīng)力滿足Boussinesq假定和靜水壓強(qiáng)假定，則二維淺水風(fēng)暴潮波控制方程組的向量形式為：

式中：U為守恒量向量；F、G為通量向量；S為源項(xiàng)向量。各項(xiàng)分別為：

式中：u、v分別為x、y方向的流速；D=h+ζ為總水深，其中h為靜水深，ζ為波高；g為重力加速度；f=2Ωsinφ為柯氏力參數(shù)，其中Ω=7.29× 10-5rad/s為地球自轉(zhuǎn)角速度，φ為計(jì)算區(qū)域的地理緯度；Zb為海底的底面高程；ρ為海水密度，這里取為常量；Pa為計(jì)算區(qū)域的氣壓；τax、τay分別為水體表面沿x、y方向的風(fēng)應(yīng)力分量；(τbx,τby)=ρCb(u,v)為底摩擦效應(yīng)項(xiàng)，其中為底摩擦系數(shù)，n為曼寧系數(shù)。

2.2 氣壓場(chǎng)模式

采用藤田氣壓場(chǎng)模式：

式中：ΔP=P∞-P0，其中P∞為臺(tái)風(fēng)外圍氣壓，這里取為標(biāo)準(zhǔn)大氣壓1.013 25×105Pa；P0為臺(tái)風(fēng)中心氣壓；為計(jì)算區(qū)域任意點(diǎn)與臺(tái)風(fēng)中心的距離；r0為臺(tái)風(fēng)最大風(fēng)速半徑。

2.3 風(fēng)場(chǎng)模式

采用宮崎正衛(wèi)風(fēng)場(chǎng)模式：

此模式由兩部分組成：一是由臺(tái)風(fēng)中心移動(dòng)而產(chǎn)生的風(fēng)速；二是由臺(tái)風(fēng)氣壓梯度產(chǎn)生的對(duì)稱梯度風(fēng)速。Wx、Wy分別為沿x、y軸的風(fēng)速分量；C1為風(fēng)場(chǎng)系數(shù)，取值范圍為（4/7，6/7）；C2為梯度風(fēng)系數(shù)，弱臺(tái)風(fēng)一般取0.6，較強(qiáng)臺(tái)風(fēng)取值范圍為0.7～0.8，也可由實(shí)測(cè)值確定；β是梯度風(fēng)與海面風(fēng)的訂正角，這里取值為30°；V0x、V0y為臺(tái)風(fēng)中心移動(dòng)速度在x、y方向上的分量；ρa(bǔ)為空氣密度；為參數(shù)；不同時(shí)刻的V0x、V0y、P0可根據(jù)臺(tái)風(fēng)年鑒資料線性插值得到。

2.4 初始條件與邊界條件

初始條件：u=0；v=0；ζ=ζ0。

邊界條件：對(duì)于水邊界可分為急流開(kāi)邊界和緩流開(kāi)邊界。邊界線的法向量記為(nx,ny)，邊界單元上的流速設(shè)為(u,v)，則法向流速和切向流速分別為un=(u,v)·(nx,ny)、uτ=(u,v)·(-ny,nx)。設(shè)與邊界鄰接的內(nèi)單元（位于邊界左側(cè)）上的水深、法向和切向速度分別為DL、un,L、uτ,L，則邊界處右單元（水邊界單元）的狀態(tài)可由表1中公式計(jì)算。

表1中DR=DL+ζR，其中ζR是邊界處臺(tái)風(fēng)氣壓降引起的海面靜壓升

陸地邊界采用鏡像法：DR=DL，un,R=-un,L，uτ,R=uτ,L。

表1 開(kāi)邊界條件

3 有限體積模型的建立

3.1 空間離散

3.1.1 有限體積離散

對(duì)雙曲守恒型方程組（1），令H=(F,G)，則式（1）可表示為：

采用無(wú)結(jié)構(gòu)三角形網(wǎng)格單元對(duì)計(jì)算區(qū)域進(jìn)行離散，將單一的網(wǎng)格單元作為控制元，并在每個(gè)單元的中心點(diǎn)設(shè)定物理變量。將第i個(gè)三角形控制元記為Ωi，在Ωi上對(duì)式（2）進(jìn)行積分，利用Green公式將單元積分化為線積分，得：

式中：dΩ為面積分微元，dl為線積分微元，Lk(k=1,2,3)表示Ωi的第k條邊，nk=(nkx,nky)= (cosθ,sinθ)為L(zhǎng)k的單位外法向量。

對(duì)式（6）中的線積分采用Gauss積分公式可得控制方程（6）的有限體積半離散化格式：

3.1.2 WENO重構(gòu)

UL、UR的計(jì)算方法是有限體積模式空間離散中的另一個(gè)重要環(huán)節(jié)，它決定著數(shù)值模式的空間精度和分辨率，這里采用高精度的WENO格式對(duì)UL、UR進(jìn)行重構(gòu)。

（1）二次多項(xiàng)式的構(gòu)造

在Gauss積分點(diǎn)對(duì)點(diǎn)值進(jìn)行高階重構(gòu)。首先構(gòu)造一個(gè)二次多項(xiàng)式P(x,y)，使在當(dāng)前單元Ω(0)上嚴(yán)格成立：

而在其余單元上使下式在最小二乘意義下成立：

因此在某一固定Gauss點(diǎn)(xG,yG)上二次多項(xiàng)式可表示成如下形式：

式中：m為大模板單元總數(shù)，這里取m=10，gl為每個(gè)單元對(duì)應(yīng)的常系數(shù)。

（2）一次多項(xiàng)式的構(gòu)造

構(gòu)造一次多項(xiàng)式需將大模板S分解成若干個(gè)小模板，然后在每個(gè)小模板上構(gòu)造一次多項(xiàng)式。

通過(guò)求解上式，可得重構(gòu)的線性多項(xiàng)式R(x,y)，并可達(dá)到二次多項(xiàng)式的精度，因此UL(Gj,t)= R(xGj,yGj)，UR(Gj,t)=可對(duì)精確解U(Gj,t)的近似達(dá)三階精度[13]。其中R表示以Ω0為中心單元構(gòu)建的模板的一次函數(shù)；記與當(dāng)前單元Ω0相鄰且邊界含(xGj,yGj)的單元為Ωn，則R?表示以Ωn為中心單元構(gòu)建的模板的一次函數(shù)。

3.2 時(shí)間離散

二維淺水波控制方程（1）經(jīng)過(guò)空間離散后，可表示為：

式中：L(U)是空間離散化算子。

采用三步Runge-Kutta時(shí)間離散方法對(duì)式（13）進(jìn)行離散來(lái)獲得三階精度的時(shí)間離散格式：

在進(jìn)行時(shí)間離散時(shí)，時(shí)間步長(zhǎng)Δt可用下式確定：

式中：lj為三角形網(wǎng)格單元的內(nèi)接圓直徑；CFL條件數(shù)v取值為0.2。

4 通量梯度項(xiàng)與源項(xiàng)的平衡

在有限體積格式的離散過(guò)程中，為將控制方程組寫成守恒型方程組，把表面梯度項(xiàng)人為地分成通量梯度項(xiàng)和源項(xiàng)，但會(huì)因此產(chǎn)生數(shù)值的不平衡。因此必須對(duì)通量梯度項(xiàng)和源項(xiàng)進(jìn)行平衡。采用Rogers提出的方法[15]，令U=Ueq+～U，其中Ueq是平衡量，使模擬能夠平衡的穩(wěn)定條件。

圖1 江蘇省沿海地形圖

5 臺(tái)風(fēng)風(fēng)暴潮數(shù)值模式在江蘇沿岸中的應(yīng)用與結(jié)果分析

圖2 三港站7910號(hào)臺(tái)風(fēng)風(fēng)暴潮計(jì)算增水與實(shí)測(cè)增水對(duì)比（1979年8月22日12時(shí)起）

為驗(yàn)證所建立的臺(tái)風(fēng)風(fēng)暴潮數(shù)值模式的可靠性，將該模式應(yīng)用到江蘇沿海臺(tái)風(fēng)風(fēng)暴潮的模擬中。選取29.988°～36.236°N和119.190°～123.965°E范圍之間的江蘇沿海區(qū)域（見(jiàn)圖1）。采用無(wú)結(jié)構(gòu)三角網(wǎng)格進(jìn)行剖分，為能夠充分反映研究區(qū)域及附近水域的變化情況，同時(shí)又不會(huì)帶來(lái)太大的計(jì)算量，對(duì)研究海域的網(wǎng)格進(jìn)行局部加密，越靠近研究區(qū)域網(wǎng)格尺寸越小，網(wǎng)格單元數(shù)為12 239個(gè)，節(jié)點(diǎn)數(shù)為6 990個(gè)。取本海區(qū)多年月平均氣壓值的均值1 013 hPa，風(fēng)場(chǎng)中β取30°。對(duì)7910號(hào)臺(tái)風(fēng)，選取r0=120km進(jìn)行計(jì)算。利用連云港港、呂四港和蘆潮港的觀測(cè)資料，進(jìn)行實(shí)測(cè)值與計(jì)算值的對(duì)比分析。

7910號(hào)臺(tái)風(fēng)計(jì)算的起始時(shí)間為1979年8月22日12時(shí)（北京時(shí)，下同），計(jì)算結(jié)束時(shí)間為1979年8月26日00時(shí)。圖2為連云港站、呂四港站和蘆潮港站計(jì)算值與實(shí)測(cè)值的對(duì)比圖。

由于這里僅僅是模擬純風(fēng)暴增水，沒(méi)有考慮天文潮等的影響，因而計(jì)算值所呈現(xiàn)出的僅是一個(gè)增水過(guò)程。另一方面，7910號(hào)臺(tái)風(fēng)風(fēng)暴潮期間的實(shí)測(cè)風(fēng)暴增水值沒(méi)有完全去除掉天文潮的影響，呈現(xiàn)出周期性的變化，但是增水趨勢(shì)比較明顯。由上述3個(gè)測(cè)站的實(shí)測(cè)增水與計(jì)算增水對(duì)比圖可以看出，7910號(hào)臺(tái)風(fēng)風(fēng)暴潮期間，3個(gè)測(cè)站的計(jì)算增水趨勢(shì)與實(shí)測(cè)增水趨勢(shì)基本吻合，而且計(jì)算所得的增水曲線基本都在實(shí)測(cè)值的高低值之間，由此可以看出，若不考慮天文潮的影響，計(jì)算所得的風(fēng)暴增水與實(shí)際的風(fēng)暴增水還是比較接近的。圖3展示了1979年8月24日00時(shí)，7910號(hào)臺(tái)風(fēng)旋轉(zhuǎn)風(fēng)場(chǎng)北側(cè)引起的流場(chǎng)運(yùn)動(dòng)。由于此刻該海域處于臺(tái)風(fēng)7910的北側(cè)，向岸的大風(fēng)把海水推向岸邊，且由于近岸處水深變淺，從而水流速度增大。綜上所述，我們建立的風(fēng)暴增水計(jì)算模式可以較準(zhǔn)確地模擬風(fēng)暴增水和流場(chǎng)。本文建立的模型可望進(jìn)一步應(yīng)用在涌潮和風(fēng)暴潮的相互作用等方面的研究中。

圖3 1979年8月24日00時(shí)7910號(hào)臺(tái)風(fēng)引起的局部海域流場(chǎng)圖

[1]馮士筰.風(fēng)暴潮導(dǎo)論[M].北京:科學(xué)出版社,1982:1-28.

[2]王喜年,尹慶江,張保明.中國(guó)海臺(tái)風(fēng)風(fēng)暴潮預(yù)報(bào)模式的研究與應(yīng)用[J].水科學(xué)進(jìn)展,1991,2(1):1-10.

[3]于福江,張占海.一個(gè)東海嵌套網(wǎng)格臺(tái)風(fēng)暴潮數(shù)值預(yù)報(bào)模式的研制與應(yīng)用[J].海洋學(xué)報(bào)(中文版),2002,24(4):23-33.

[4]Aschariyaphotha N,Wongwises P,Humphries U W,et al.Study of storm surge due to typhoon Linda(1997)in the gulf of Thailand using a three dimensional ocean model[J].Applied Mathematics and Computation,2011,217(21):8640-8654.

[5]Wang R Y,Li W,Chen Y D,et al.Numerical forecasting model for storm surge based upon unstructured meshes and finite volume method[C]//Third Chinese-German Joint Symposium on Coastal and Ocean Engineering.Tainan:National Cheng Kung University, 2006.

[6]李燕初,蔡文理.風(fēng)暴潮的有限元分析[J].海洋學(xué)報(bào)(中文版), 1982,4(4):404-414.

[7]李云川,張迎新,王福俠,等.2003年10月風(fēng)暴潮的形成及數(shù)值模擬分析[J].氣象,2005,31(11):15-18.

[8]Choi B H,Eum H M,Woo S B.Modeling of coupled tide-wavesurge process in the Yellow Sea[J].Ocean Engineering,2003,30 (6):739-759.

[9]Kim S Y,Yasuda T,Mase H.Numerical analysis of effects of tidal variations on storm surges and waves[J].Applied Ocean Research, 2008,30(4):311-322.

[10]Xie L,Liu H Q,Liu B,et al.A numerical study of the effect of hurricane wind asymmetry on storm surge and inundation[J]. Ocean Modelling,2011,36(1-2):71-79.

[11]Liu X D,Osher S,Chan T.Weighted essentially non-oscillatory schemes[J].Journal of Computational Physics,1994,115(1):200-212.

[12]Jiang G S,Shu C W.Efficient implementation of weighted ENO schemes[J].Journal of Computational Physics,1996,126(1):202-228.

[13]Montarnal P,Shu C W.Real gas computation using an energy relaxation method and high-order WENO schemes[J].Journal of Computational Physics,1999,148(1):59-80.

[14]Balsara D S,Shu C W.Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high orderof accuracy[J].Journal of Computational Physics,2000,160(2): 405-452.

[15]Rogers B D,Borthwick A G L,Taylor P H.Mathematical balancing of flux gradient and source terms prior to using Roe’s approximateRiemannsolver[J].JournalofComputational Physics,2003,192(2):422-451.

A high-order and high-resolution numerical simulation of storm surge based on the WENO method

WANG Ru-yun1，WANG Tian1，WU Chu-min1，QIANG Dan-dan2，ZHOU Jun3，ZHANG Xin1，ZHANG Bin1，ZHAN Fei1
（1.College of Harbor,Coastal and Offshore Engineering,Hohai University,Nanjing 210098 China; 2.Nantong Productivity Promotion Center,Nantong 226019 China; 3.College of Hydrology and Water Resources,Hohai University,Nanjing 210098 China）

Considering the non-linear effects of typhoon storm surge in coastal areas,based on the unstructured meshes,the numerical model of two-dimensional typhoon storm surge was established by using the finite volume method,high-order high-resolution WENO scheme to complete the space discretization of the two-dimensional typhoon storm surge equation,the third-order Runge-Kutta scheme for the time discretization,and the Rogers method to solve the problem of the imbalance between the flux gradient and the discrete source terms,which was caused by the complex submarine topography.Fujita formula and Veno Takeo formula are used to simulate pressure and wind in the model,respectively.At last,through the simulation and verification of the typhoon storm surge along Jiangsu coastal areas,it’s proved that the simulation of typhoon storm surge by this numerical model is validity and feasibility.

WENO scheme;unstructured meshes;typhoon storm surge;numerical simulation;high-precision and high-resolution

P731.23

A

1003-0239（2017）02-0021-06

10.11737/j.issn.1003-0239.2017.02.003

2016-05-16；

2016-08-01。

中國(guó)江蘇省水利科技重點(diǎn)項(xiàng)目（2010500312）；中央高校基本科研業(yè)務(wù)費(fèi)專項(xiàng)資金（2014B06314）；江蘇省普通高校研究生科研創(chuàng)新計(jì)劃項(xiàng)目（CXZZ12_0224）。

王如云（1963-），男，教授，博士，主要從事物理海洋學(xué)研究。E-mail：wangry@hhu.edu.cn