基于WENO格式的高精度高分辨臺風風暴潮數(shù)值模式

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

基于WENO格式的高精度高分辨臺風風暴潮數(shù)值模式

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

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

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

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

1 引言

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

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

2 基本方程和初邊界條件

2.1 二維淺水風暴潮波控制方程

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

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

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

2.2 氣壓場模式

采用藤田氣壓場模式：

式中：ΔP=P∞-P0，其中P∞為臺風外圍氣壓，這里取為標準大氣壓1.013 25×105Pa；P0為臺風中心氣壓；為計算區(qū)域任意點與臺風中心的距離；r0為臺風最大風速半徑。

2.3 風場模式

采用宮崎正衛(wèi)風場模式：

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

2.4 初始條件與邊界條件

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

邊界條件：對于水邊界可分為急流開邊界和緩流開邊界。邊界線的法向量記為(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中公式計算。

表1中DR=DL+ζR，其中ζR是邊界處臺風氣壓降引起的海面靜壓升

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

表1 開邊界條件

3 有限體積模型的建立

3.1 空間離散

3.1.1 有限體積離散

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

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

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

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

3.1.2 WENO重構(gòu)

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

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

在Gauss積分點對點值進行高階重構(gòu)。首先構(gòu)造一個二次多項式P(x,y)，使在當前單元Ω(0)上嚴格成立：

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

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

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

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

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

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

3.2 時間離散

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

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

采用三步Runge-Kutta時間離散方法對式（13）進行離散來獲得三階精度的時間離散格式：

在進行時間離散時，時間步長Δt可用下式確定：

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

4 通量梯度項與源項的平衡

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

圖1 江蘇省沿海地形圖

5 臺風風暴潮數(shù)值模式在江蘇沿岸中的應用與結(jié)果分析

圖2 三港站7910號臺風風暴潮計算增水與實測增水對比（1979年8月22日12時起）

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

7910號臺風計算的起始時間為1979年8月22日12時（北京時，下同），計算結(jié)束時間為1979年8月26日00時。圖2為連云港站、呂四港站和蘆潮港站計算值與實測值的對比圖。

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

圖3 1979年8月24日00時7910號臺風引起的局部海域流場圖

[1]馮士筰.風暴潮導論[M].北京:科學出版社,1982:1-28.

[2]王喜年,尹慶江,張保明.中國海臺風風暴潮預報模式的研究與應用[J].水科學進展,1991,2(1):1-10.

[3]于福江,張占海.一個東海嵌套網(wǎng)格臺風暴潮數(shù)值預報模式的研制與應用[J].海洋學報(中文版),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]李燕初,蔡文理.風暴潮的有限元分析[J].海洋學報(中文版), 1982,4(4):404-414.

[7]李云川,張迎新,王福俠,等.2003年10月風暴潮的形成及數(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。

中國江蘇省水利科技重點項目（2010500312）；中央高?；究蒲袠I(yè)務(wù)費專項資金（2014B06314）；江蘇省普通高校研究生科研創(chuàng)新計劃項目（CXZZ12_0224）。

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