基于WENO格式的高精度高分辨台风暴潮数值模拟(投稿海洋会议)

1基于WENO格式的高精度高分辨台风暴潮数值模式王如云1，羌丹丹1，周钧2，张彬1，张鑫1，占飞1(1.河海大学港口海岸与近海工程学院，江苏南京，210098；2.河海大学水文水资源学院，江苏南京，210098)摘要：考虑到台风暴潮在近岸浅水地区的非线性效应，可能会产生陡度较大的波面，甚至出现间断波现象，基于无结构网格，并利用WENO格式设计思想，设计建立了高精度高分辨率的二维台风暴潮数值模式。模式中所用的风场和气压场，分别采用宫崎正卫风场模式和藤田气压场模式。采用Rogers方法解决复杂海底地形造成的通量梯度项与源项数值离散后的不平衡问题。通过对一维溃坝波的模拟，验证了模式对间断波进行模拟的可行性。通过对江苏沿海的台风暴潮增水的模拟和验证，表明该模式对台风暴潮进行数值模拟是可行的。关键词：WENO格式；无结构网格；台风暴潮；数值模式；高精度高分辨中图分类号：P731.23台风暴潮是自然界中一种严重的，由台风引起的海水水位异常升降的灾害现象。若与天文高潮相叠加，往往会使其影响所及的海域潮位暴涨，甚至海水漫溢陆地，吞没码头、城镇、村庄，酿成巨灾[1]，因此对风暴潮的研究具有实际意义。上世纪90年代，王喜年等[2]建立了以五个子域重复覆盖整个中国海岸海域的台风风暴潮模式(FBM模式)。李云川等[3]运用非静力MM5的二重网格双向嵌套数值模式，模拟了2003年渤海湾的风暴潮过程。B.H.Choi等[4]通过耦合一个三代的波浪模型、WAM模型和二维的风暴潮模型，建立了耦合潮汐、风暴潮和波浪的模式，并将其应用到黄海海域中去。SooYoulKim等[5]将沿深度积分的非线性浅水波方程和SWAM模型相结合，用2006年的艾云尼台风作为算例，分析了风暴潮和风浪对天文潮的影响。LianXie等[6]运用三维风暴潮模式，针对五种不同类型的数值实验来研究不对称风场对风暴潮和漫滩的影响。另外，根据实际模拟地域的特征，一般可建立直角坐标系、极坐标系或球坐标系[5,7,8]下的控制方程。台风暴潮在近岸浅水地区的非线性效应，可能会产生陡度较大的波面，甚至出现间断波现象，这些问题是上述数值模式得到进一步发展的瓶颈。为了解决这些问题，2006年，王如云等[9]基于有限体积法，建立了高分辨率的风暴潮二维数值预报模式，并对珠江口的风基金项目：中国江苏省水利科技重点项目(2010500312),自然科学基金(40906048)作者简介：王如云(1963-),教授,男,安徽芜湖人,从事计算物理学研究,E-mail:[e-mail]暴潮进行了模拟。本文基于无结构网格，并利用WENO[10]格式设计思想，设计建立了可模2拟间断波现象的高精度高分辨率的二维台风暴潮数值模式。运用该数值模式对江苏沿海地区台风暴潮引起的纯台风暴潮增水进行数值模拟，并对7910号台风的增水过程进行验证，取得了比较满意的结果。1基本控制方程假设水体是均质不可压的，并且水流紊动应力满足Boussinesq假定和静水压强假定，则平面二维风暴潮运动的守恒型方程为()()txyUFUGUS（1）其中U为守恒量向量；F、G为通量向量；S为源项向量。,,TDuDvDU22,2,TuDuDgDuvDF22,,2TvDuvDvDgDG0,,TaybybaaxbxbaZPZPDDDfvgDDfugDxxyyS式中u、v分别为x和y方向的流速分量，Dh为总水深，其中为波高，h为静水深，2sinf为柯氏力参数，这里57.2910rads为地球自转角速度，为计算区域的地理纬度；bZ为海底的底面高程；为海水密度，这里取为常量；ax，ay为水体表面沿x轴和y轴的风应力分量；aP为计算区域的气压；bx，by为沿x和y方向的底摩擦效应项，其表达式为：vuCbbybx,,，31222hvugnCb其中，n为曼宁系数。风场采用宫崎正卫风场模式:31002522041,,exp115102xyxyarfPWWCVVCZfr00sin,coscos,sinxxyy此模式由两部分组成：一是由台风中心移动而产生的风速；二是由台风气压梯度产生的对称梯度风速。,xyWW分别为沿x轴和y轴的风速分量；1C为风场系数，取值范围为(4/7,6/7)；2C为梯度风系数，弱台风一般取0.6，稍强或强台风的取值范围为(0.7,0.8)，3也可由实测值确定，是梯度风与海面风的订正角；00,xyVV为台风中心移动速度的x和y方向分量；a为空气密度；1/2201rZr；0PPP，其中P为无穷远处的气压，0P为台风中心的气压；不同时刻的000,,xyVVP可根据台风年鉴中的资料线性插值得到。给定风场以后，海面上任一点的风应力为：22,,axayawxyxyC其中，a为空气密度，220.00110.07wxyCWW为风应力系数，xW，yW为沿x轴和y轴方向的风速（水面以上10m的风速）。气压场采用藤田公式：200()1arPPPPr其中0P为台风中心气压；P为无穷远处气压，即台风外围气压；r为计算区域任意点与台风中心的距离；0r为台风最大风速半径，它的选取将直接影响风场模式中的风速轮廓线，亦即风场的真实性。因此，最大风速半径对风场的模拟起着关键作用。计算最大风速半径的方法有很多，通常根据具体台风中心附近的实测气压值或实测风速反求，也有采用七级大风半径的十分之一作为最大风速半径。由于本文缺少相应的气压实测资料，且目前主要着眼于该模式对台风暴潮数值模拟的可行性研究上,因此本文选取固定的0120rkm。初始条件：0uv。边界条件：对于水边界可分为急流开边界和缓流开边界。设边界处的已知状态为LD、,nLu、,Lu，其中L表示边界处的左单元（已知单元），(,)(,)nxyuuvnn，(,)(,)yxuuvnn，分别表示单元法向流速和切向流速。未知单元的状态如表1所示：表1.开边界条件Table1.Openboundarycondition缓流急流水边界条件,,2()nRnLLRuugDD,,nRnLuu,,RLuu,,RLuu附：表中RLRDD，其中R是边界处台风气压降引起的海面静压升高400.009911/RPrr。陆地边界采用镜像法[11]：RLDD，,,nRnLuu，,,RLuu。2构建有限体积模式：对双曲守恒型方程组（1），令(,)HFG，yx,为梯度算子，则可表示为：tUHS(2)本文采用三角形网格单元对计算区域进行离散，将单一的网格单元作为控制元，物理变量配置在每个单元的重心。将第i个三角形控制元记为i，在i上对式(2)进行积分，利用Green公式将单元积分化为线积分，得：31iikiikLkdddlddtUHnS(3)其中d为面积分微元，dl为线积分微元，ikL3,2,1k表示i的第k条边，,cos,siniixiynnn为ikL的单位外法向量（见图1）。图1三角形网格控制元及其外法向示意图Map1.Triangulargridcontrolelementandschematicdiagramoftheouternormalvector对式（3）中的线积分，利用Gauss积分公式和Lax-Friedrichs数值通量，可以得到控制方程(3)的有限体积半离散化格式：3111ˆ,qiikijijRijLikikjidULFSdtUUn(4)其中Lax-Friedrichs(LF)通量格式如下：inＸＹiixniyn51ˆ,,,,,2LFjRjLLjRjRjLjFttttUUnHUGHUGnUGUG这里，q是单元i上的高斯点个数，ij为高斯点对应的权，jG为第j个高斯点。3WENO格式LU、RU的重构是有限体积模式空间离散中的另一个重要环节，它决定着数值模式的空间精度和分辨率，本文采用高精度的WENO格式对LU、RU进行重构。（1）二次多项式的构造在Gauss积分点对点值进行三阶重构，需要构造一个二次多项式(,)Pxy，使得该多项式在当前三角单元0上满足：0001(,)PxydxdyU(5)记0的重心为00(,)xy，定义变量00xxh和00yyh，其中1200h。可记二次多项式为：22123456(,)Pxybbbbbb(6)为利用最小二乘法定出系数126,,,bbb，需要当前三角形单元0及其相邻单元信息。以当前单元0为中心的大模板S还包括：三个相邻单元,,ijk，i的两个相邻单元,abii，j的两个相邻单元,abjj，k的两个相邻单元,abkk。为便于表述将大模板0{;,,;,;,;,}abababijkiijjkkS，按顺序编号对应为：(1)(2)(3)(4)(5)(6)(7)(8)(9)(10){,,,,,,,,,}S我们要求二次多项式(,)Pxy在当前单元(1)上使下式严格成立：(1)0(1)1(,)PxydxdyU(7)而在其余单元上使下式在最小二乘意义下成立：()1()1(,)2,3,,10lllPxydxdyUl(8)记6(),1()11llled(),2()1llled(),3()1llled()2,4()1llled()2,5()1llled(),6()1llled其中1,2,,10l。由式(7)可得：0121,231,341,451,561,6bUbebebebebe将上式代入式(6)得：021,231,3(,)()()PxyUbebe2241,451,561,6()()()bebebe(9)将式(9)代入式(7)得：02,21,23,31,3()()llUbeebee4,41,45,51,56,61,61()()()llllbeebeebeeU2,3,,10l构造函数：1002362,21,23,31,321(,,...,)[()()2lllbbbUbeebee214,41,45,51,56,61,6()()()]llllbeebeebeeU根据最小二乘法原理，需使所构造函数最小，分别对236,,,bbb求导，并令其等于零可得：10101010,1,12,2,13,5,16,110222210101010,1,22,2,23,5,26,21022221010,1,52,2,5322()()()()()()()()()()(lllllllllllllllllllllllllllllllaabaabaabaUUaabaabaabaUUaabaaba1010,5,56,51022)()lllllabaUU其中：,,11,1ljljjaee，2,3,,101,2,,5lj求解上式可得到236,,,bbb，然后可由式(9)求出1b。实际上，126,,,bbb均是019,,,UUU的函数，因此在某一固定Gauss点(,)GGxy上二次多项式可表示成如下形式：710(,)mGGlllPxygU(10)其中m为大模板单元总数，这里取=10m，lg为每个单元对应的常系数。（2