基于WENO格式的高精度高分辨台风风暴潮数值模式

王如云，汪 天，吴楚敏，羌丹丹，周 钧，张 鑫，张 彬，占飞

基于WENO格式的高精度高分辨台风风暴潮数值模式

王如云1，汪 天2，吴楚敏1，羌丹丹3，周 钧4，张 鑫1，张 彬1，占飞1

（1.河海大学海洋学院，江苏南京210098；2.河海大学港航学院，江苏南京210098；3.南通市生产力促进中心，江苏南通226019；4.河海大学水文水资源学院，江苏南京210098)

考虑到台风风暴潮在近岸浅水地区的非线性效应，基于无结构网格，通过采用有限体积法和高精度高分辨率的WENO数值格式对二维浅水方程进行空间离散，并利用三阶的Runge-Kutta格式进行时间离散，最后利用Rogers方法解决复杂海底地形造成的通量梯度项与源项数值离散后的不平衡问题，从而建立了二维台风风暴潮数值模式。模式中的风场和气压场分别采用宫崎正卫风场模式和藤田气压场模式。最后通过对江苏沿海的风暴增水的模拟和验证，表明了该数值模式对台风风暴潮模拟的有效性和可行性。

WENO格式；无结构网格；台风风暴潮；数值模式；高精度高分辨

1 引言

台风风暴潮是由台风引起的海水水位异常升降的现象[1]，对其进行准确、快速数值模拟方法方面的研究，一直是人们关注的问题。从数值计算的网格来讲，为了适应复杂岸边界形态，人们从最初采用结构网格[2]，逐渐采用嵌套网格[3]、曲线坐标网格[4]，目前为止，不少研究工作基于适应性更好、便于整个流场设计统一计算格式而采用无结构网格[5]。从数值计算的方程离散方法来讲，主要是有限差分法[2]、有限元法[6]、有限体积法[5]。另外，在已有的水动力学数值模型基础上，加上台风场开展风暴潮数值模拟的工作也有不少成果，如李云川等[7]运用非静力MM5的二重网格双向嵌套数值模式，模拟了2003年渤海湾的风暴潮过程。Choi等[8]通过耦合一个三代的波浪模型、WAM模型和二维的风暴潮模型，建立了耦合潮汐、风暴潮和波浪的模式，并将其应用到黄海海域中去。Kim等[9]将沿深度积分的非线性浅水波方程和SWAM模型相结合，用2006年的艾云尼台风作为算例，分析了风暴潮和风浪对天文潮的影响。Xie等[10]运用三维风暴潮模式，针对5种不同类型的数值实验来研究不对称风场对风暴潮和漫滩的影响。

台风风暴潮在近岸浅水地区的非线性效应，可能会产生陡度较大的波面，甚至出现间断波现象，这些问题是上述数值模式进一步发展的瓶颈。为解决这些问题，2006年，Wang等[5]基于无结构网格、有限体积法、Roe平均法，尝试建立了高分辨率的风暴潮二维数值预报模式，并对珠江口的风暴潮进行了模拟。考虑到Roe格式的参向量、间断解展开、线性化替代矩阵等，都是在最简单的左右常数态的假定下进行的，这种左右常数态的假定，从根本上限制了方法的精度。1994年，Liu等[11]提出并构造了一维空间上的三阶有限体积加权本质无震荡(WENO)格式，它是由网格平均形式的ENO格式发展起来的。其后，Jiang等[12]在1996年构造了多维空间上的三阶和五阶有限差分WENO格式，提出了现在被广泛采用的光滑因子和非线性权[13-14]。本文基于无结构网格，并利用WENO格式设计思想，设计建立了可模拟间断波现象的高精度高分辨率的二维台风风暴潮数值模式。运用该数值模式对江苏沿海地区台风风暴潮引起的纯风暴增水进行数值模拟，并对7910号台风的增水过程进行验证，取得了比较满意的结果。

2 基本方程和初边界条件

2.1 二维浅水风暴潮波控制方程

假设水体是均质不可压的，并且水流紊动应力满足Boussinesq假定和静水压强假定，则二维浅水风暴潮波控制方程组的向量形式为：

式中：U为守恒量向量；F、G为通量向量；S为源项向量。各项分别为：

式中：u、v分别为x、y方向的流速；D=h+ζ为总水深，其中h为静水深，ζ为波高；g为重力加速度；f=2Ωsinφ为柯氏力参数，其中Ω=7.29× 10-5rad/s为地球自转角速度，φ为计算区域的地理纬度；Zb为海底的底面高程；ρ为海水密度，这里取为常量；Pa为计算区域的气压；τax、τay分别为水体表面沿x、y方向的风应力分量；(τbx,τby)=ρCb(u,v)为底摩擦效应项，其中为底摩擦系数，n为曼宁系数。

2.2 气压场模式

采用藤田气压场模式：

式中：ΔP=P∞-P0，其中P∞为台风外围气压，这里取为标准大气压1.013 25×105Pa；P0为台风中心气压；为计算区域任意点与台风中心的距离；r0为台风最大风速半径。

2.3 风场模式

采用宫崎正卫风场模式：

此模式由两部分组成：一是由台风中心移动而产生的风速；二是由台风气压梯度产生的对称梯度风速。Wx、Wy分别为沿x、y轴的风速分量；C1为风场系数，取值范围为（4/7，6/7）；C2为梯度风系数，弱台风一般取0.6，较强台风取值范围为0.7～0.8，也可由实测值确定；β是梯度风与海面风的订正角，这里取值为30°；V0x、V0y为台风中心移动速度在x、y方向上的分量；ρa为空气密度；为参数；不同时刻的V0x、V0y、P0可根据台风年鉴资料线性插值得到。

2.4 初始条件与边界条件

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

边界条件：对于水边界可分为急流开边界和缓流开边界。边界线的法向量记为(nx,ny)，边界单元上的流速设为(u,v)，则法向流速和切向流速分别为un=(u,v)·(nx,ny)、uτ=(u,v)·(-ny,nx)。设与边界邻接的内单元（位于边界左侧）上的水深、法向和切向速度分别为DL、un,L、uτ,L，则边界处右单元（水边界单元）的状态可由表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）可表示为：

采用无结构三角形网格单元对计算区域进行离散，将单一的网格单元作为控制元，并在每个单元的中心点设定物理变量。将第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重构

UL、UR的计算方法是有限体积模式空间离散中的另一个重要环节，它决定着数值模式的空间精度和分辨率，这里采用高精度的WENO格式对UL、UR进行重构。

（1）二次多项式的构造

在Gauss积分点对点值进行高阶重构。首先构造一个二次多项式P(x,y)，使在当前单元Ω(0)上严格成立：

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

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

式中：m为大模板单元总数，这里取m=10，gl为每个单元对应的常系数。

（2）一次多项式的构造

构造一次多项式需将大模板S分解成若干个小模板，然后在每个小模板上构造一次多项式。

通过求解上式，可得重构的线性多项式R(x,y)，并可达到二次多项式的精度，因此UL(Gj,t)= R(xGj,yGj)，UR(Gj,t)=可对精确解U(Gj,t)的近似达三阶精度[13]。其中R表示以Ω0为中心单元构建的模板的一次函数；记与当前单元Ω0相邻且边界含(xGj,yGj)的单元为Ωn，则Rˆ表示以Ωn为中心单元构建的模板的一次函数。

3.2 时间离散

二维浅水波控制方程（1）经过空间离散后，可表示为：

式中：L(U)是空间离散化算子。

采用三步Runge-Kutta时间离散方法对式（13）进行离散来获得三阶精度的时间离散格式：

在进行时间离散时，时间步长Δt可用下式确定：

式中：lj为三角形网格单元的内接圆直径；CFL条件数v取值为0.2。

4 通量梯度项与源项的平衡

在有限体积格式的离散过程中，为将控制方程组写成守恒型方程组，把表面梯度项人为地分成通量梯度项和源项，但会因此产生数值的不平衡。因此必须对通量梯度项和源项进行平衡。采用Rogers提出的方法[15]，令U=Ueq+～U，其中Ueq是平衡量，使模拟能够平衡的稳定条件。

图1 江苏省沿海地形图

5 台风风暴潮数值模式在江苏沿岸中的应用与结果分析

图2 三港站7910号台风风暴潮计算增水与实测增水对比（1979年8月22日12时起）

为验证所建立的台风风暴潮数值模式的可靠性，将该模式应用到江苏沿海台风风暴潮的模拟中。选取29.988°～36.236°N和119.190°～123.965°E范围之间的江苏沿海区域（见图1）。采用无结构三角网格进行剖分，为能够充分反映研究区域及附近水域的变化情况，同时又不会带来太大的计算量，对研究海域的网格进行局部加密，越靠近研究区域网格尺寸越小，网格单元数为12 239个，节点数为6 990个。取本海区多年月平均气压值的均值1 013 hPa，风场中β取30°。对7910号台风，选取r0=120km进行计算。利用连云港港、吕四港和芦潮港的观测资料，进行实测值与计算值的对比分析。

7910号台风计算的起始时间为1979年8月22日12时（北京时，下同），计算结束时间为1979年8月26日00时。图2为连云港站、吕四港站和芦潮港站计算值与实测值的对比图。

由于这里仅仅是模拟纯风暴增水，没有考虑天文潮等的影响，因而计算值所呈现出的仅是一个增水过程。另一方面，7910号台风风暴潮期间的实测风暴增水值没有完全去除掉天文潮的影响，呈现出周期性的变化，但是增水趋势比较明显。由上述3个测站的实测增水与计算增水对比图可以看出，7910号台风风暴潮期间，3个测站的计算增水趋势与实测增水趋势基本吻合，而且计算所得的增水曲线基本都在实测值的高低值之间，由此可以看出，若不考虑天文潮的影响，计算所得的风暴增水与实际的风暴增水还是比较接近的。图3展示了1979年8月24日00时，7910号台风旋转风场北侧引起的流场运动。由于此刻该海域处于台风7910的北侧，向岸的大风把海水推向岸边，且由于近岸处水深变浅，从而水流速度增大。综上所述，我们建立的风暴增水计算模式可以较准确地模拟风暴增水和流场。本文建立的模型可望进一步应用在涌潮和风暴潮的相互作用等方面的研究中。

图3 1979年8月24日00时7910号台风引起的局部海域流场图

[1]冯士筰.风暴潮导论[M].北京:科学出版社,1982:1-28.

[2]王喜年,尹庆江,张保明.中国海台风风暴潮预报模式的研究与应用[J].水科学进展,1991,2(1):1-10.

[3]于福江,张占海.一个东海嵌套网格台风暴潮数值预报模式的研制与应用[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月风暴潮的形成及数值模拟分析[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）；中央高校基本科研业务费专项资金（2014B06314）；江苏省普通高校研究生科研创新计划项目（CXZZ12_0224）。

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