湍黏系数对浅滩海域三维风暴潮的影响

熊 伟,朱志夏,董 佳,齐庆辉

(1.江苏省交通规划设计院股份有限公司,江苏 南京 210014;2.江苏省水运工程技术研究中心,江苏 南京 210014)

湍黏系数对浅滩海域三维风暴潮的影响

熊 伟1,2,朱志夏1,2,董 佳1,2,齐庆辉1,2

(1.江苏省交通规划设计院股份有限公司,江苏 南京 210014;2.江苏省水运工程技术研究中心,江苏 南京 210014)

利用WRF大气模式和基于有限体积法的三维海洋数值模型FVCOM,建立了连云港及其附近浅滩海域的三维风暴潮流数学模型。运用k-ε和MY-2.5两种不同的湍流模型,对“韦帕”台风作用下的三维风暴潮流进行了数值模拟,并对计算得到的三维风暴潮流流场结果和垂向湍流黏性系数结果进行了对比分析。计算分析结果表明,两种不同湍流模型计算所得的风暴潮增水差别很小,可以忽略;但在风生流较强的时间段内,MY-2.5湍流模型计算所得的水平流速流向沿垂向分布并不一致,而k-ε湍流模型的计算结果较好,与实测分层潮流资料更为符合。研究结果表明,垂向湍流黏性系数沿水深的垂向分布对浅滩海域风暴潮流水平流速的垂向结构至关重要,建议选用k-ε湍流模型计算垂向湍流黏性系数和湍流扩散系数。

湍黏系数;浅滩海域;三维数值模型;风暴潮;垂向分布

风暴潮是指强烈的大气扰动(强风和气压骤变)引起的海面异常升高或降低的现象[1]。风暴潮位若遇到天文大潮,就会发生大的潮灾。而伴随强风产生的波浪和水流流速很大,使得浅滩海岸泥沙的启动变得异常容易,其运动也变得较为复杂,极易形成港口航道工程的淤积,如具有浅滩淤泥质海岸性质的连云港海域。在大风浪作用下,水流和泥沙具有明显的分层现象,因此加强风暴潮特别是三维风暴潮的数值模拟研究也显得尤为必要。

20世纪末,风暴潮数值模拟研究大多采用二维数值模型计算增减水[2-4]。随着计算机的发展,近年来三维风暴潮数值模式得到了较大发展,出现了许多海洋数值模型如POM,ECOM和FVCOM等[5]。从物理角度看,三维模型可以用来描述风暴潮的水流垂向结构,并且三维模型在近底处的底部摩阻假设更为真实[6]。强风天气下的水流运动、泥沙浓度分层现象较为明显。因此,对风暴潮作用下的三维水流进行准确的模拟,是研究港口航道工程骤淤问题的前提条件。李大鸣等[7]对渤海湾的风暴潮进行了三维数值模拟,曹丛华等[8]利用FVCOM模型在胶州湾进行了三维风暴潮漫滩数值模拟研究,但对三维风暴潮研究都只进行了水位的验证,缺乏三维流速流向的验证分析;张娜等[9]利用SWAN和MOHID联合模型对连云港深水航道进行了风暴潮骤淤研究,取得了较好的成果,但未对三维水流结构进行详细分析研究。陈永平等[10]讨论了垂向紊动黏性系数对三维潮流数值计算结果的影响,结果表明垂向紊动黏性系数越大,流速分层越不明显,流速的绝对值也越小,垂向紊动黏性系数的梯度变化则对水平流速的垂向分布起着决定性的作用。本文利用基于有限体积法的FVCOM模型和已有的三维潮流资料,对连云港海域在“韦帕”台风作用下的三维风暴潮进行数值模拟,着重研究风暴潮作用下的水流垂向结构和不同湍流模型的应用对比情况。

1 数学模型的建立

“韦帕”(WIPHA)台风于2007年9月16日在西北太平洋洋面上生成,在浙江省苍南县霞关镇登陆,登陆时台风中心附近最大风速达45 m/s,台风中心向北移动,过境江苏省连云港市,对该海域造成一定的影响。有关研究单位在连云港徐圩海域对此次台风进行了连续现场观测,获得了较为详细的波浪潮流泥沙等观测资料。本文建立台风模型和三维潮流模型,对该时段的三维风暴潮流进行数值模拟研究。

1.1 大气模式

WRF(Weather Research Forecast)大气模式为新一代中尺度预报模式与同化系统,具有可移植、易维护、可扩充、高效率等特性,可用于数值天气预报的研究、大气-海洋模式的耦合等。目前,WRF模式包括2个动力框架ARW(the Advanced Research WRF)和NMM(the Nonhydrostatic Mesoscale Model),分别由美国国家大气研究中心(NCAR)和美国国家环境预报中心(NCEP)主要开发并维护更新。本文使用的是WRF-ARW动力框架,WRF-ARW系统主要由4部分组成:WRF预处理系统(WPS)、WRF数据同化系统(WRFDA)、ARW求解、程序后处理以及可视化工具[11]。

WRF模式是一个完全可压的非静力模式,控制方程都写为通量的形式。垂直方向采用地形跟随的质量坐标,水平方向采用Arakawa C网格[12]。

1.2 风暴潮数值模型

FVCOM采用非结构化的三角形网格单元,能更精确拟合复杂不规则的岸界,如处理连云港东西连岛海域岸线;采用有限体积法,能够在整个计算区域内更好地保证各物理量的守恒;采用垂向σ坐标,可以体现不规则变化的海底地形,并在浅水区域具有更高的垂向分辨率。另外,该数值模型程序代码的并行化计算,大大缩短了数值计算求解的时间。其具体控制方程和数值离散求解模式详见文献[5]。

适当地描述和计算湍流,对于三维水流的模拟、预测泥沙输运和海床变化具有重要的影响, L.O.Amoudry等人详细研究了湍流模型对泥沙分层和输运的影响[13]。本文在研究过程中发现,选用不同的湍流模型,得到的三维水流计算结果差异较大。在此,对湍流模型作简要介绍。FVCOM海洋数值模型提供了多样性的湍流封闭模型,用于计算垂向涡黏性系数和垂向扩散系数。FVCOM自带了改进过的MY-2.5湍流封闭模型,此外,它通过耦合GOTM湍流模式,也包含了k-ε湍流封闭模型供研究人员选择。

1.2.1 MY-2.5湍流模型 目前比较常用的MY-2.5湍流模型,它包括湍动能方程和混合长度方程,可表示如下:

式中:~W=1+E2l2/(κL)2为固壁迫近函数,κ=0.4为卡曼常数;E1,E3和B1为模型闭合常数;Km,Kh和Kq分别为垂向涡黏系数、垂向扩散系数和垂向湍流扩散系数;q2/2为湍动能;l为湍流混合长度。MY-2.5模型中,表面及底面的边界条件分别为:

式中:η和H分别为潮位和水深。

1.2.2 k-ε湍流模型 GOTM湍流模式实现了众多的从简单的Richardson参数化到复杂的雷诺应力等湍流模型,其中k-ε湍流模型可表示为:

式中:Km为垂向涡黏系数;k为湍流动能;ε为湍流耗散;经验常数cμ,c1,c2,σk,σε的取值分别为0.09,1.44, 1.92,1.00,1.30。表面及底面边界条件分别为:

2 数学模型的验证

2.1 韦帕台风验证

为了给潮流计算提供更高精度的风场数据,减小边界对计算结果的影响,并模拟出台风条件下的台风路径。本次研究中,“韦帕”台风的风场模型计算区域采用大模型、中模型与小模型三重网格双向嵌套技术,即大模型的计算结果直接影响小模型,且小模型的计算结果也会对大模型进行反馈。模型的计算范围与网格精度如表1所示,模型垂直方向分为34层。模拟的台风路径与实际路径对比见图1,风速、风向与连云港大西山海洋站实测数据对比见图2。由图2可见,建立的WRF大气模式能较好地反映连云港海域风场状况。

图1 模拟台风路径与实际路径对比Fig.1 Comparison between simulated typhoon path and actual path

表1 WRF模型计算范围与网格精度Tab.1 Range and accuracy of grid in WRF model

图2 台风风速风向对比验证Fig.2 Validation of typhoon speed and direction

2.2 三维风暴潮验证

本次风暴潮数值模拟研究采用二维大模型、三维小模型双层嵌套的方法计算。二维大模型主要为三维小模型的计算提供合理的边界条件。二维、三维模型的计算范围如图3所示。模型均采用非结构化的三角形网格,网格最小尺度为30 m,位于连云港主港区和主航道;网格最大尺度为2 000 m,位于计算域的离岸开边界处。测站相对位置如图4所示,1#高程为-3.0 m,2#高程为-5.0 m(高程基准为当地理论基面),但1#点没有流速验证资料,因此采用2#的实测潮位流速资料进行验证。

图3 双重嵌套模型计算范围Fig.3 Calculated range of nested models

图4 计算验证点位置Fig.4 Position of measured points

计算时间为2007-09-17T0:00/2007-09-21T20:00,整个过程共计117 h。三维模型的计算时间步长内模取1.0 s,外模为内模的5倍;垂向均匀分为10层。由于连云港海域基本为淤泥质海岸性质,摩擦阻力系数相对较小,底部粗糙厚度取0.5 mm。风场数据源于WRF模型计算得到的结果。需要说明的是,连云港海域水深较浅,水流垂向结构较为复杂。因此,本次数值模拟研究中湍流模型选MY-2.5和k-ε两种模型分别计算,计算验证结果如图5和6所示。

图5 2#观测点潮位验证Fig.5 Verification of tidal level for 2#measured point

图6 2#观测点流速流向验证Fig.6 Verification of velocity and direction of current for 2#measured point

从潮位计算结果可知,两种湍流模型都能够很好地模拟风暴潮位的过程,潮位过程差别不大,这说明湍流黏性系数数值变化对风暴潮位的影响微弱;从分层水流计算结果可知,相对于MY-2.5模型来说,k-ε模型的计算结果和实测资料更相符。很显然,在其他基本参数一致的情况下,三维风暴潮流结果差异显著的原因是采用的湍流模型不同造成的。下面通过MY-2.5模型和k-ε模型的计算,对流速过程和垂向涡黏性系数等进行对比分析。

3 三维水流结构与湍黏系数

选取2个特征性的代表时刻进行对比分析。风前时刻为2007-09-18T06:00,模拟风速为0.94 m/s;大风时刻为09-20T08:00,模拟风速为16.31 m/s。2个时刻的表层与底层流场分别如图7和8所示。可以看出,风前时刻流场保持着正常天气下的潮流涨落特性,但是在风中时刻,受大风影响,k-ε模型此刻近岸水体表底层流向基本一致,全变为沿岸流的性质,这与实测资料基本吻合。这说明在浅滩区域,浅滩海域的水体在垂向上作为一个整体对台风风暴作出响应。而大风时刻MY-2.5模型的计算结果中,表底层流向差异较大,表现为表层是沿岸流,底层近乎为离岸流,这与实测情况不符。

图7 大风时刻MY-2.5模型计算的表、底层流场Fig.7 Surface and bottom flow field calculated by MY-2.5 model at strong wind moment

图8 大风时刻k-ε模型计算的表、底层流场Fig.8 Surface and bottom flow field calculated by k-ε model during strong wind

计算输出2个测点在风前、风中两个不同时刻的垂向涡黏性系数Km和垂向扩散系数Kh(见图9)。由于Km与Kh具有相关性,仅对Km作详细比较,Kh类同。

从不同水深点的Km值比较来看,风前时刻与大风时刻,-5.0与-3.0 m的Km值垂向分布趋势基本一致,不同的是-5.0 m的湍流黏性系数值比-3.0 m稍大,这与两者水流流速大小有关。另外,k-ε模型计算得到表底面处的Km为0,这与模型的表底面边界条件式(8)和(9)有关。湍流模型中没有考虑波浪的影响,因此表底面的Km计算值与实际情况存在差别。

从不同时刻来看,如图9(a)和(b)所示,风前时刻的两种模型计算的湍流黏性系数分布具有些许相似,都呈现中间层大、表底层小的特征,但是Km的最大值分布位置不同,MY-2.5模型偏底层,而k-ε模型偏表层。图9(c)和(d)则说明两种模型计算结果差异很大,MY-2.5模型的Km表现出了间断性不连续现象,且Km最大值分布在表层,数值达1.0左右;而k-ε模型的结果依旧与风前时刻相似,但Km的最大值在垂向上的位置分布向下移,更趋于水深中间层,数值也大了许多,-3.0和-5.0 m处的Km的最大值分别为0.014和0.019 m2/s。需要强调的是,MY-2.5模型的Km表现出了间断性不连续现象,造成了三维风暴潮流在垂向流速分布上不一致,且间断层与水流垂向流速流向的突变层均位于0.6H层。为此,本文尝试在MY-2.5模型计算中将垂向分层层数加大为20,但Km值依旧在0.6H层出现间断性分布。这说明,在浅滩海域的三维风暴潮流数值模拟计算中,垂向湍流黏性系数这一参数的计算至关重要,它关系到风暴潮流水平流速的垂向结构,与文献[10]中的结论一致。另外,由于Km与Kh具有相关性,而垂向扩散系数Kh值也关系到风暴潮流作用下泥沙、污染物等物质的输运扩散;因此应选用适当的湍流模型计算垂向涡黏性系数和扩散系数。

图9 风前时刻和大风时刻的测站垂向湍流黏性系数和扩散系数分布Fig.9 Vertical distribution of turbulent viscosity coefficients and diffusion coefficients before and during strong wind

4 结 语

本文对“韦帕”台风作用下连云港浅滩海域的三维风暴潮流进行数值模拟和研究,采用k-ε模型得到了比较好的计算结果,与实测风暴潮流资料符合较好,并得出以下结论:

(1)在浅滩海域的三维风暴潮流数值模拟计算中,垂向湍流黏性系数这一参数的计算和取值至关重要,它关系到风暴潮流的垂向结构;应选用适当的湍流模型,建议采用k-ε湍流模型进行计算。

(2)在大风作用达到一定时间之后,浅滩海域的水体在垂向上作为一个整体对台风风暴作出响应。

[1]冯士筰.风暴潮导论[M].北京:科学出版社,1982:1-2.(FENG Shi-zuo.Storm surge introduction[M].Beijing:Science Press,1982:1-2.(in Chinese))

[2]李树华,陈文广,梁善任,等.广西沿海8007号台风暴潮数值模拟及台风暴潮某些特性的分析[J].广西科学,1994,1 (2):14-20.(LI Shu-hua,CHEN Weng-guang,LIANG Shan-ren.Numerical simulation and certain characteristics analysis on the storm surge No.8007 along the coast of Guangxi[J].Guangxi Sciences,1994,1(2):14-20.(in Chinese))

[3]余晓军,吕开华.9618号台风风暴潮特征及其数值模拟[J].海洋通报,1998,17(6):7-11.(YU Xiao-jun,LV Kai-hua. Characteristics and numerical simulation of the storm surge caused by typhone No.9618[J].Marine Science Bulletin,1998,17 (6):7-11.(in Chinese))

[4]柴扉,汪景庸.东海风暴潮与天文潮的非线性相互作用[J].青岛海洋大学学报,1990,20(3):56-62.(CAI Fei,WANG Jing-yong.The nonlinear interaction of surge and tide in the East China Sea[J].Journal of Ocean University of Qingdao,1990, 20(3):56-62.(in Chinese))

[5]熊伟.非结构网格下的三维潮流泥沙数值模拟[D].大连:大连理工大学,2012:3-4.(XIONG Wei.3D numericalsimulation of tidal currents and sediment based on unstructured mesh model[D].Dalian:Dalian University of Technology,2012: 3-4.(in Chinese))

[6]黄华,朱建荣,吴辉.长江口与杭州湾风暴潮三维数值模拟[J].华东师范大学学报:自然科学版,2007(4):9-19. (HUANG Hua,ZHU Jian-rong,WU Hui.3-D numerical simulation of storm surge in the Changjiang estuary and the Hangzhou bay[J].Journal of East China Normal University(Natural Science),2007(4):9-19.(in Chinese))

[7]李大鸣,范玉海,徐亚男,等.风暴潮三维数值计算模式的研究及在渤海湾的应用[J].海洋科学,2012,36(7):7-13. (LI Da-ming,FAN Yu-hai,XU Ya-nan.Research on 3-D numerical calculation model of storm surge and application in the Bohai Bay[J].Marine Sciences,2012,36(7):7-13.(in Chinese))

[8]曹丛华,白涛,高松,等.胶州湾高分辨率三维风暴潮漫滩数值模拟[J].海洋科学,2013,37(2):118-125.(CAO Cong-hua,BAI Tao,GAO Song,et al.High resolution 3D storm surge and inundation numerical model used in the Jiaozhou bay [J].Marine Sciences,2013,37(2):118-125.(in Chinese))

[9]张娜,杨华,严冰,等.SWAN和MOHID联合模型在连云港港30万吨级航道建设中的应用[J].水运工程,2012(2): 128-133.(ZHANG Na,YANG Hua,YAN Bing,et al.Application of SWN/MOHID-based three-dimensional numerical model in 300 000 DWT deep-water channel of Lianyungang port[J].Port&Waterway Engineering,2012(2):128-133.(in Chinese))

[10]陈永平,刘家驹,喻国华.潮流数值模拟中紊动粘性系数的研究[J].河海大学学报:自然科学版,2002,30(1):39-43.(CHEN Yong-ping,LIU Jia-ju,YU Guo-hua.A study of eddy viscosity coefficient in numerical tidal simulation[J]. Journal of Hohai University(Natural Sciences),2002,30(1):39-43.(in Chinese))

[11]WANG W,BRUYÈRE C,DUDA M,et al.ARW Version 3 modeling system user′s guide:mesoscale and microscale meteorology division[M].Boulder:National Center for Atmosphere Research,2011.

[12]王勇.基于WRF模式的酒泉风电基地近地层风场特征模拟研究[D].兰州:兰州大学,2012:25-26.(WANG Yong. Simulation and analysis of wind characteristics near surface layer in Jiuquan with WRF[D].Lanzhou:Lanzhou University, 2012:25-26.(in Chinese))

[13]AMOUDRY L O,SOUZA A J.Impact of sediment-induced stratification and turbulence closures on sediment transport and morphological modeling[J].Continental Shelf Research,2011,31(9):912-928.

Effects of turbulent viscosity coefficient on 3-D storm surge within shallow seas

XIONG Wei1,2,ZHU Zhi-xia1,2,DONG Jia1,2,QI Qing-hui1,2
(1.Jiangsu Province Communications Planning and Design Institute Limited Company,Nanjing 210014,China; 2.Jiangsu Province Waterway Engineering Technology Research Center,Nanjing 210014,China)

Using a WRF atmospheric model and a 3-D ocean numerical model FVCOM based on finite volume method,a three-dimensional numerical model of storm surge within Lianyungang shallow waters is developed.The 3-D storm surge during typhoon WIPHA is simulated by use of different turbulence models such as k-ε and MY-2.5.The simulation results of the 3-D storm surge flow and vertical turbulence viscosity coefficient are analyzed in the study.The analysis results show that,compared with the poor results from MY-2.5 model,the vertical velocity and direction distribution of the storm surge calculated by k-ε model has a good agreement with the measured data,and the values of the vertical turbulence viscosity coefficient are crucial for the vertical structure of the storm surge,therefore a right turbulence model is preferred.The k-ε turbulence model is suggested for simulation of the 3-D storm surge in the future.

viscosity coefficients;shallow seas;3D numerical simulation;storm surge;vertical distribution

P731.23

A

1009-640X(2014)06-0045-07

2014-04-26

国家高新技术研究发展计划(863计划)资助项目(2012AA112509)

熊 伟(1988-),男,湖北荆州人,硕士,主要从事海岸近海水动力与泥沙研究。E-mail:xw_7816@126.com通信作者:朱志夏(E-mail:zhixiazhu@sina.com)