杜建廷,华 锋,王道龙,江兴杰,江志辉
(1.海洋环境科学与数值模拟国家海洋局重点实验室,山东 青岛266061;2.国家海洋局 第一海洋研究所,山东 青岛266061)
MASNUM-WAM(Key Lab of Marine Science and Numerical Modeling-WAM),其前身为 LAGFDWAM[1-2](Laboratory of Geophysical Fluid Dynamics-WAM),是当今世界上最先进的第三代海浪数值模式之一,自1991年研发至今一直应用于中国近海海浪的数值模拟以及海洋工程中,其在高风速条件下模拟海浪的能力受到国内外研究人员的广泛认可。
MASNUM-WAM的破碎耗散源函数采用Yuan等改进的参数化方案[3],非线性波-波相互作用源函数参照Hasselmann等的参数化方案[4-5],在数值模式中首次考虑了波流相互作用源函数[1],数值计算格式上采用复杂特征线嵌入计算格式[6]。杨永增等在LAGFD-WAM区域海浪数值模式基础上建立了球坐标系下的全球海浪数值模式[7]。华锋等为LAGFD-WAM区域海浪数值模式设计了一种新的特征线数值计算格式[8],使其能够更好地进行浅水海浪数值模拟。
但是,大量的数值试验和工程计算结果表明,MASNUM-WAM在模拟低风速条件下的海浪时,其有效波高较实际偏低。在人类涉足的海域,除台风、寒潮等恶劣天气之外,海面大多处于低风速条件下,因此在海上平台设计、海岸建筑设计、海浪能估计等工程应用中,往往需要统计一般条件下的海浪。这也在一定程度上限制了MASNUM-WAM的应用。本文基于量纲分析方法改进了模式中的谱增长限制,使其更好地应用于低风速条件下海浪的数值模拟。
第三代海浪数值模式普遍采用谱能平衡方程[1]求解海浪谱,即。其中,为局地变化项与平流变化项之和;为波动群速度;为背景流速;Stot为各源函数之和,深水波在不考虑波-流相互作用的前提下,由风输入项、波-波非线性相互作用项以及破碎耗散项三部分组成:Stot=Sin+Snl+Sdis。在MASNUM-WAM中,能量的平流变化项采用特征线嵌入格式计算,而局地变化项与源函数项采用半隐半显格式处理。采用半隐半显格式数值积分平衡方程时,源函数对能谱的Jacobian导数矩阵只保留对角线项不能保证数值稳定性,海浪模式中通常采用谱的增长限制来抑制不稳定过程。
最初,在 WAM 模式中谱的增长限制采用与时间步长无关的形式:|ΔF|max=6.4×10-7g2f-5[9]。Tolman根据数值模拟结果指出,这种形式对时间步长具有很强的依赖性[10]。为解除这种依赖性,在 WAM Cylce4中采用了比例与时间步长的限制条件:|ΔF|max=6.4×10-7g2f-5Δt/τ(τ=1 200s)[11]。该模式应用于欧洲中期天气预报中心(European Centre For Medium Range Weather Forecasts)的业务化海浪预报系统当中。对于网格步长Δx>50km的情况,该模式取得了较好的预报结果。但是在高分辨率下,模式结果与有限风区海浪的成长规律不相符。为此,Hersbach和Janssen[12-13]基于特征尺度的不变性对增长限制进行了改进。改进后的增长限制更少地限制了海浪谱的成长,尤其是在高频谱段的成长,使得WAM在高分辨率下模拟海浪的能力有了很大的提高。改进后的增长限制为,其中,为无量纲PM谱峰频率,fc为截断频率。Yuan等[1]在LAGFD-WAM海浪数值模式中采用Phillips增长限制[14],。其中,为波数谱;θ为风向与波向的夹角;k为波数;c为波速;u*为摩擦速度;p,q为常数,通常p=0.025,q=2。
为消除其对时间步长的依赖性,在MASNUM-WAM中采用了比例与时间步长的限制条件:
其中E=E(k),为方向平均的波数谱;α为常数,α=9.1×10-6p(p=0.025);u*为摩擦速度;k为波数;Δt为时间步长。
当设定模式的时间步长和波数的划分方案,谱的增长限制只是摩擦速度的函数。而在MASNUMWAM中,摩擦速度仅仅是海面10m风速的函数u*=f(u10)。对于模式中的任意一个计算点,波谱在该点的增长限制只与风速有关,即ΔEmax=g(u10),而同样的风速在不同的波谱下所产生的谱能变化(ΔEs)是不同的。因此,为研究谱增长限制与风速以及与波谱之间的关系,采用Tsagareli的复合谱[15],方向分布采用Babanin and Soloviev方案[16],给出风速分别为5m/s和30m/s,无因次风区分别为χ=1×105,χ=1×101和χ=1×103的谱增长限制与源函数导致的谱能增量之间的关系。
图1 不同风速和风区下谱增长限制对沿风速方向的谱能增量的影响Fig.1 Effects of the growth spectrum limit on the spectrum energy increment in the wind direction under the conditions of different wind speeds and fetches
当无因次风区较小时,增长限制的绝对值小于源函数导致的谱能增量的绝对值,从而削减了谱能的增长(图1a和图1c),在低风速下的削减明显大于高风速的;而当无因次分区较大时,增长限制对谱能的增长影响较小(图1b和图1d)。在同一风速下,不同的波数谱对应的谱能增量不同,这就要求谱增长限制也要做出相应的调整。这些都表明以往采用的谱增长限制随风速或者随谱形的变化不合理是导致低海况下模拟有效波高较实际偏低的一个重要因素。由此推测,要改进谱增长限制可以在原谱增长限制的基础上降低u*的指数或者加入一个表征谱形的物理量。
MASNUM-WAM中的各个方程都是建立在以摩擦速度为基本量的量纲分析基础之上的。量纲分析需要用到以下3个无量纲量:(无量纲波数谱的谱能量),(无量纲波数)(无量纲时间)。模式中所用到的方程,必须满足量纲齐次原则。
将MASNUM-WAM中采用的增长限制式(1)写成无量纲量与特征量乘积的形式:
消去特征量得到无量纲形式的方程:
由此可知原谱增长限制中的参数α是一个隐含量纲的常数,其量纲为,该量纲与波数的量纲一致。
Yuan等最初在LAGFD-WAM海浪数值模式中采用的是Phillips谱增长限制,其方向平均的形式[1]:
其中β为待定系数。将其写成无量纲量与特征量的乘积形式:
消去特征量得到无量纲形式的方程:
当右式加入Δt项时,为满足量纲齐次原则,需在无量纲形式方程中加入无量纲时间,即
还原量纲:
或
参考式(1)中k-4关系,将式(9)中的k-3改为k-4,同时根据我们数值试验结果,加入表征谱峰位置的量
其中,
选取不用的时间步长和网格划分,经过大量数值试验得出,β取值为1.0×10-6时,可以得到稳定的结果。
将改进后的谱增长限制应用于MASNUM-WAM,采用风浪的有限风区成长关系作为参考标准对改进后的模式结果进行检验。检验内容包括无因次谱能随无因次风区的成长关系以及无因次谱峰波数随无因次风区的成长关系。其中,无因次谱能为,无因次谱峰波数为,无因次风区为。作为参考标准的成长关系包括Jonswap成长关系,其无因次谱能、无因次谱峰波数随无因次风区的成长关系为ε=1.6×10-7χ,Kp=3.5χ-0.33[17];Donelan等成长关系,其无因次谱能、无因次谱峰波数随无因次风区的成长关系为ε=8.415×10-7χ0.76,Kp=1.85χ-0.23[18];Young成长关系,其无因次谱能、无因次谱峰波数随无因次风区的成长关系为ε=7.5×10-7χ0.8,Kp=2.0χ-0.25[19]。我们共进行4组对比试验,其中 EXP1,EXP2采用的波数划分范围为0.007 1~0.688 8m-1,EXP3,EXP4采用的波数划分范围为0.001 0~4.390 9m-1,EXP1,EXP3采用原谱增长限制,EXP2,EXP4采用改进后的谱增长限制。
图2 改进谱增长限制前后模式结果的有限风区成长关系对比(波数划分范围:0.007 1~0.688 8m-1)Fig.2 Comparison of fetch-limited growth relationships of the model results before and after the improvement of growth spectrum limit(wavenumber range:0.007 1~0.688 8m-1)
图2为EXP1与EXP2的模式结果检验。波数划分范围采用MASNUM-WAM中通常采取的波数划分方案,即0.007 1~0.688 8m-1。对比可知,当风速为15~45m/s时,改进谱增长限制后的模式结果与参考值符合程度较高,数据也相对集中,而改进之前的模式结果与参考值符合程度较低,数据也相对离散。当风速为5m/s时,改进后的模式结果,其无因次谱能在小风区大于改进前,而在大风区则差别不大,二者均小于参考值。这是因为在低风速下,当风区较小时,谱能增量集中于波数较大的区域(图1a),采用0.688 8m-1为最大波数划分范围严重限制了海浪谱的成长,使得模式结果偏小。在高风速下,当风区较大时,谱能增量集中于波数较小的区域(图1d),采用0.007 1m-1为最小波数划分范围也会在一定程度上限制海浪谱的成长。基于上述2个原因,为了提高模式的准确度,EXP3与EXP4在保证数值稳定性的基础上扩大波数划分范围为0.001 0~4.390 9m-1。
图3为EXP3与EXP4的模式结果检验。对比可知,与参考值相比,在不同风速下改进后的模式结果比改进前均有不同程度的提高。在低风速下主要表现在小风区,在高风速下主要表现在大风区。同时,由于波数划分范围大于EXP1与EXP2,因此在极端条件下EXP3与EXP4的模式结果与参考值更接近。
图3 改进谱增长限制前后模式结果的有限风区成长关系对比(波数划分范围:0.001 0~4.390 9m-1)Fig.3 Comparison of fetch-limited growth relationships of the model results before and after the improvement of growth spectrum limit(wavenumber range:0.001 0~4.390 9m-1)
第三代海浪数值模式中普遍引入谱增长限制。谱增长限制作用于谱能增量,防止其在一个时间步长内超过合理的范围,保证了模式计算的稳定性。过分宽松的谱增长限制会导致模式计算溢出,而过紧的谱增长限制则会限制谱峰波数较大、谱宽度较宽的风浪谱的成长,导致海浪在低风速下和成长初期成长缓慢。
在实际应用中,大量的数值试验和工程计算结果表明,MASNUM-WAM在低风速下模拟海浪有效波高偏低的重要原因之一来自谱增长限制。基于量纲齐次原则,在原谱增长限制的基础上加入了平均波数,更符合风浪谱的成长规律,在保证模式计算稳定性的前提下更少地限制了风浪谱的成长。
采用风浪的有限风区成长关系作为参考标准对改进后的模式进行检验。检验结果表明,改进后的模式结果在低风速和小风区情况下更接近理论值,在高风速和大风区等极端条件下也与理论值更为接近。同时,当风速的变化范围较大时,应在保证计算稳定的前提下适当扩大波数划分范围。目前上述结论尚处于试验阶段,距离实际应用尚有一定差距,需要进一步采用实测资料进行模式的检验和参数的调整。
(References):
[1]YUAN Y L,HUA F,PAN Z D,et al.LAGFD-WAM numerical wave model-I.Basic physical model[J].Acta Oceanologica Sinica,1991,10(4):483-488.
[2]YUAN Y L,PAN Z D,HUA F,et al.LAGFD-WAM numerical wave model-II.Characteristics inlaid scheme and its application[J].Acta Oceanologica Sinica,1991,11(1):13-23.
[3]YUAN Y L,HUA F,PAN Z D,et al.Dissipation source function and an improvement to LAGFD-WAM model[J].Acta Oceanologica Sinica,1992,11(4):471-481.
[4]HASSELMANN S,HASSELMANN K.Computations and parameterizations of the nonlinear energy transfer in gravity-wave spectrum-Part I:A new method for efficient computations of the exact nonlinear transfer integral[J].Journal of Physical Oceanography,1985,15:1369-1377.
[5]HASSELMANN S,HASSELMANN K.Computations and parameterizations of the nonlinear energy transfer in gravity-wave spectrum-Part II:Parameterizations of the nonlinear energy transfer for application in wave models[J].Journal of Physical Oceanography,1985,15:1378-1391.
[6]PAN Z D,SUN L T,HUA F,et al.LAGFD-Ⅱ.regional numerical wave model and its application-Ⅱ.Characteristics inlaid computational scheme[J].Oceanologica et Limnologia Sinica,1992,23(5):459-467.潘增弟,孙乐涛,华锋,等.LAGFD-Ⅱ区域性海浪数值模式及其应用-Ⅱ.特征线嵌入网格计算方法[J].海洋与湖沼,1992,23(5):459-467.
[7]YANG Y Z,QIAO F L,ZHAO W,et a1.MASNUM ocean wave numerical model in spherical coordinates and its application[J].Acta Oceanologica Sinica,2005,27(2):l-7.杨永增,乔方利,赵伟,等.球坐标系下 MASNUM海浪数值模式的建立及其应用[J].海洋学报,2005,27(2):1-7.
[8]HUA F,WANG D L,YUAN Y L,et al.Characteristics computational scheme of numerical wave model under complex topography conditions[J].Advances in Marine Science,2005,23(3):272-280.华锋,王道龙,袁业立,等.复杂地形下海浪数值模式的特征线计算格式[J].海洋科学进展,2005,23(3):272-280.
[9]WAMDI Group.The WAM model-A third generation ocean wave prediction model[J].Journal of Physical Oceanography,1988,18:1775-1810.
[10]TOLMAN H L.Effects of numerics on the physics in a third-generation wind-wave model[J].Journal of Physical Oceanography,1992,22:1095-1111.
[11]KOMEN G J,CAVALERI L,DONELAN M,et al.Dynamics and modelling of ocean waves[M].Cambridge:Cambridge University Press,1994.
[12]HERSBACH H,JANSSEN P A E M.Improvement of the short-fetch behavior in the Wave Ocean Model(WAM)[J].Journal of Atmospheric and Oceanic Technology,1999,16:884-892.
[13]TOLMAN H L.Limiters in third-generation wind wave models[J].The Global Atmosphere and Ocean System,2001,8(1):67-83.
[14]PHILLIPS O M.Spectral and statistical properties of the equilibrium range of wind-generated gravity waves[J].Journal of Fluid Mechanics,1985,156:505-531.
[15]TSAGARELI K.Numerical investigation of wind input and spectral dissipation in evolution of wind waves[D].Adelaide,SA,Australia:University of Adelaide,2008.
[16]BABANIN A V,SOLOVIEV Y P.Variability of directional spectra of wind-generated waves,studied by means of wave staff arrays[J].Marine and Freshwater Research,1998,49(2):89-101.
[17]HASSELMANN K,BARNETT T P,BOUWS E,et al.Measurements of wind-wave growth and swell decay during the Joint North Sea Wave Project(JONSWAP)[R].[S.l.]:Deutsches Hydrographisches Institute.1973:7-94.
[18]DONELAN M A,HAMILTON J,HUI W H.Directional spectra of wind-generated waves[J].Philosphical Transaction of the Royal Society of London,1985,315(1534):509-562.
[19]YOUNG I R.Wind generated ocean waves[M].Amsterdam:Elsevier Science,1999.