非静压方程与波作用谱模型的波浪传播嵌套模拟

邹国良，张庆河

(天津大学水利仿真与安全国家重点实验室，天津300072)

选取合理的数学模型来模拟不同范围和水深的波浪运动，对准确描述波浪传播变形、获得建筑物附近的设计波浪要素有重要的工程意义。工程中在描述波浪从大范围到小范围传播变形时通常将基于相位平均的波作用谱模型与基于求解相位的沿水深积分模型进行嵌套计算，即在大范围海域(以下简称大模型)采用计算效率较高的波作用谱模型来模拟波浪传播变形，在工程关心的小范围区域如港池附近(以下简称小模型)则采用求解相位的模型计算。小模型的边界条件由大模型提供。工程中常用的做法是由大模型推算出嵌套边界处代表波浪要素以作为小模型的嵌套边界输入条件。这与嵌套边界处的真实波况往往是不一致的，由此得到的小模型计算结果与实际情况必然会有出入。Tozer等［1-5］提出了更为合理的嵌套方法，即将波作用谱模型与浅水方程模型进行空间上谱嵌套，研究结果表明，嵌套边界条件采用沿边界变化的实际波浪谱时，复杂地形中的波浪传播变形模拟才能具有较高的精度。

近年来，越来越多的学者开始通过求解含非静压项的非线性浅水方程来模拟波浪的传播变形［6-11］。非静压模型在垂向往往只需要2～3层［6-8］即可较好地描述波浪的强非线性和强色散性，使得近岸波浪变形模拟的精度和计算效率提高。已有的针对求解相位模型与求解相位平均模型的嵌套工作主要是将波作用谱模型与平面二维的水深积分模型［1-4］及静压假定的浅水方程模型［5］进行嵌套，尚未见到波作用谱模型与三维波浪模型及非静压模型嵌套的文献。因此，为了更高效、准确地描述波浪从大尺度范围传播至中小尺度范围的传播变形，并将嵌套模型拓展到三维，本文将波作用谱模型和非静压浅水方程波浪模型进行嵌套，嵌套边界条件采用由波作用谱模型提供的空间变化谱作为非静压方程的波浪输入条件。

1 嵌套数学模型

本文的波作用谱大模型模拟采用目前应用非常广泛的第3代风浪模型SWAN模型，求解相位的小模型模拟采用考虑非静压项的非线性浅水方程模型SWASH。下面针对两种模型做简要的介绍。

1.1 波作用谱模型

SWAN模型采用波作用谱方程描述风浪生成及其在近岸区的演化过程。关于SWAN模型的控制方程详见文献［13］以及 SWAN 技术手册［14］，这里不再赘述。

SWAN模型计算域的边界类型分为陆地和水边界，其中陆地边界为波能吸收边界，水边界可根据现场实测和计算结果给定。针对多向不规则波模拟，这里在迎浪水边界采用JONSWAP频谱［15］和光易型方向分布函数的乘积作为输入的方向谱，即

式中:Hs为有效波波高;Tp和fp分别为谱峰周期和谱峰频率;f为组成波的频率;γ为谱峰因子;ms为方向集中度;θ0为波浪传播的主波向。

1.2 非静压浅水方程模型

忽略黏性影响，同时假定水体密度恒定为ρ0，SWASH模型的控制方程可表示为

式中:u、v和w分别为沿x、y和z方向的流速;η为自由面;d为静水深;h=η+d为总水深;g为重力加速度;压力p分解成静水压力g(η-z)和非静水压力q。自由面和底部的运动学边界条件分别为

自由面可通过由连续性方程(4)沿水深积分，并结合式(8)获得的自由面方程求解:

式中:Qx、Qy分别为x和y方向的流量，可表示为

忽略风应力和表面张力，自由面处的边界条件为

底部边界条件为

固边界采用自由滑移的边界条件，以平行于y轴的边界为例:

入流边界通过指定速度分布并结合弱反射边界条件［16］吸收二次反射波。在模拟不规则波时，可将方向谱按照频率分为M份，方向分为N份，共MN个谐波叠加，其速度边界为

出流边界采用改进的海绵层吸收边界［17-19］，海绵层内每个时间步内水平向速度采用显式阻尼消波，阻尼系数可表示为

2 数值方法

2.1 波作用谱模型的数值方法

SWAN模型方程求解采用有限差分格式进行离散。在空间上采用一阶迎风格式离散;在谱方向上采用迎风格式与中心差分相结合的混合格式离散，并取等间距方向步长Δθ;在频率上采用符合对数频率分布的等间距相对频率步长Δf/f。时间离散采用欧拉全隐格式，具体的离散形式可参照文献［14］。

2.2 非静压浅水方程模型的数值方法

非静压浅水方程模型SWASH基于交错网格系统对沿垂向各层积分的控制方程采用具有二阶精度的显式有限差分格式进行空间离散。动量方程中的对流项采用具有守恒特性的对流项离散格式进行离散［19］。为了减少垂向的分层数同时保证模型能较好的描述短周期波的色散特性，动量方程中的非静压梯度项采用紧致Keller-box格式进行离散，并将压力项定义在单元的边中心(如图1)，从而使表层自由面的压力边界自动满足，而不需要进行假定［7］。对于时间的积分，连续性方程和动量方程均采用二阶蛙跳格式进行时间积分。此外，模型还采用了干-湿处理技术来描述动边界。对于SWASH模型控制方程的详细数值离散过程可参照文献［6，18］

图1 SWASH交错网格变量布置Fig.1 Arrangement of the unknown with a staggered grid system in SWASH

模型在求解压力速度耦合方程时采用分步法。在静压步仅考虑非静压梯度项中的显式部分，求解新时刻的过渡流速;在非静压步求解局部质量守恒推导出的压力泊松方程获得非静压项修正项，并对过渡流速进行修正，最终通过自由面方程更新自由面。

2.3 嵌套方法

本文对2种波浪模型进行的嵌套计算并不考虑由SWASH模型所计算的反射波对SWAN模型的影响。嵌套过程需要解决2个问题:1)采用非静压浅水方程进行多向不规则波模拟时，为了使模拟的方向谱在空间上较为均匀，同时为了满足平稳性要求，方向谱按照频率分割时的数量一般要大于波作用谱模型所要求的数量;2)SWASH模型所需要的网格尺寸要远小于SWAN模型。为此，在嵌套时先采用3次样条插值函数重新插值SWAN模型计算出嵌套边界上的每一个方向上的频率谱［2〛，再将新插值的网格尺寸较粗节点处的方向谱线性插值到网格尺寸相对较小的SWASH模型的边界处。

SWASH模型嵌套边界处的方向谱采用傅里叶逆变换将其转换成多组谐波，由此可获得各组成谐波的波浪参数，如振幅、相位、圆频率等，进而根据式(15)可获得SWASH模型的输入条件。

3 数值计算结果与分析

为了说明所提出的嵌套方法的合理性，设计了两组多向不规则波在理想地形上传播的算例进行数值模拟试验。两组算例中的波浪条件以及水深地形均不考虑破碎情况的发生。

3.1 等水深波浪传播算例

算例1中大模型的计算地形为1 km×1 km，水深为-5m，网格尺寸为Δx=Δy=20m，计算地形如图2所示。多向不规则波的有效波高Hs为1 m，谱峰周期为8 s，谱峰因子为3.3，频率分布范围为0.08～1.0 Hz，按照频率分为28等份，波浪入射边界为南边界，主波向为90°，方向分布范围为75°～105°，按照方向分为21份，方向谱集中度取2。小模型计算域为400 m×350 m，南边界(图2中bnd)为嵌套边界，北边界设置100 m宽的海绵层进行消波，东西侧边界为自由滑移边界。模型计算网格尺寸为Δx=Δy=2 m，时间步长Δt=0.1 s，总计算时间为45 min。

图2 等水深计算地形Fig.2 Com putation bathymetry of hydroisobath

图3比较了由SWAN模型以及嵌套后SWASH模型计算出点P1～P3位置处的能量谱。表1为计算出的有效波高和波谱平均周期(Tm01)对比结果。其中，SWASH模型计算出时序列的水位通过快速傅里叶变化进行了能量谱估计，采样时间间隔为0.1 s，采样时长15 min。根据计算结果可看出，由采用谱作为嵌套边界输入条件的浅水方程模型计算出的能量谱、波高和波周期与用SWAN计算出的结果基本一致。

图3 P1～P3处SWAN模型和SWSAH模型计算谱比较Fig.3 Comparison of SWAN and SWASH com puted energy density spectra at P1～P3

表1 P1～P3处SWAN和SWASH计算波浪参数比较Table 1 Comparison of SWAN and SWASH computed wave parameters at P1～P3

3.2 斜坡上的波浪传播算例

为了进一步体现嵌套模型在模拟波浪发生折射变形时的准确性，本文设置了水深以1/400坡度从-10 m变化至-5 m的地形，如图4所示。为了消除SWAN计算时侧边界的影响从而得到合理的嵌套边界波浪要素，大模型计算范围为10 km×2 km，网格尺寸为Δx=Δy=20 m。多向不规则波的有效波高、谱峰周期以及方向谱等相关参数同3.1节，波浪入射边界为南边界，入射主波向为45°，方向分布为0°～90°，方向集中度为15。

图4 斜坡水深地形Fig.4 Sketch of linear beach bathymetry

图5 两种嵌套模式计算域Fig.5 Plan view of the com putation domain for two types nestingmode

模型嵌套采用2种模式:1)嵌套边界取为南边界bnd1和西边界bnd2，北、东边界为150 m宽的海绵层消波，如图5(a)所示;2)采用工程中常见做法即波浪正向入射嵌套边界，即取嵌套边界bnd3垂直于入射波向，其他边界均设置为150 m宽的海绵层，如图5(b)所示。2种模式的计算范围均为1.7 km×1 km，网格尺寸为Δx=1 m，Δy=1 m，垂向等分为2层，时间步长Δt=0.025 s，总计算时间为50 min。

图6为大模型提供的沿嵌套边界的波高分布图以及平均波向图。由图可看出:采用第1种嵌套模式时，波高和波向沿边界bnd1的变化不大，而沿边界bnd2变化较大;采用第2种嵌套模式时，波高沿嵌套边界bnd3的变化较明显，而平均波向沿边界bnd3基本在48°左右。

图6 沿嵌套边界的有效波高和平均波向分布Fig.6 Significant wave height and mean direction along nesting boundaries

图7、8比较了2种嵌套模式条件下大模型以及小模型嵌套计算出的沿同一位置S1～S2以及S3～S4的有效波高分布以及能量谱的对比结果(图中MODE1表示嵌套模式1，MODE2为嵌套模式2)。

图7 2种嵌套模式下SWAN和SWASH模型计算出有效波高分布对比Fig.7 Comparison of SWAN and SWASH computed wave significant wave height for two types of nestingmode

图8 S1～S4处SWAN模型和SWSAH模型计算谱比较Fig.8 Comparison of SWAN-and SWASH-computed energy density spectra at S1～S4

根据上述波浪在等水深地形传播和在斜坡上进行浅水折射变形的算例结果，嵌套后的非静压模型计算波高、周期与波作用谱模型接近，且两者计算出的谱比较吻合。这说明所建立的嵌套模型可以合理描述波浪从大范围区域向小范围区域的传播，所建立的嵌套模型可为工程区小范围波浪模拟(如港内波高分布等)提供一种新的计算模式。

表2则给出相应的计算波参数对比结果.根据图表结果可知，由两种嵌套模式分别计算出的有效波高沿程分布与采用波作用谱模型计算出的有效波高均基本一致，能量谱也基本接近，部分位置处的波高变动主要是因为波浪在斜坡上传播会产生部分反射所造成。

表2 S1～S4处SWAN和嵌套模式计算出的波浪参数比较Table 2 Comparison of SWAN and nesting mode computed wave parameters at S1～S4

4 结束语

通过波作用谱模型提供给非静压模型嵌套边界处的波浪谱参数，并对其在空间和频率方向上进行插值使其满足求解相位模型的边界条件，较好地实现了非静压浅水方程模型与波作用谱模型的嵌套。波浪在等水深地形传播和斜向入射波在斜坡上传播的2个算例结果表明，嵌套模型可以获得较为准确的波高、周期以及波浪谱，可以合理描述波浪从大范围区域向小范围区域的传播，在解决类似于港内波高分布等工程波浪计算问题时，可以获得更高精度的结果。非静压波浪模型与谱作用模型的嵌套有着十分广阔的应用前景。

［1］TOZER N P，DURAND N.Improved coastal wave modelling using coupled spectral wave transformation and wave disturbance models［C］//Proceedings of 28th International Conference on Coastal Engineering.Wales，UK，2002:448-458.

［2］OKIK，SAKAIT.Coupling of phase-averaging model and phase-resolvingmodel by using double connection boundaries［C］//Proceedings of 31st International Conference on Coastal Engineering.Hamburg，Germany，2008:302-313.

［3］CHENW，DEMIRBILEK Z，LIN L.Coupling phase-resolving nearshore wave model with phase-averaged spectral wave models in coastal applications［C］//Proceedings of 11st International Conference on Estuarine and CoastalModeling.Seattle，Washington，USA，2009:582-600.

［4］MCCUBEM.Coupled wave action and shallow-watermodelling for random wave runup on a slop［J］.Journal of Hydraulic Research，2011，49(4):515-522.

［5］RUSU E，SOARES C G.Modeling waves in open coastal areas and harbourswith phase-resolving and phase-averaged models［J］.Journal of Coastal Research，2013，29(1):1309-1325.

［6］ZIJLEMA M，STELLINGG S，SMITP.SWASH:An operational public domain code for simulating wave fields and rapidly varied flows in coastalwaters［J］.Coastal Engineering，2011，58(10):992-1012.

［7］ZIJLEMA M，STELLING G S.Further experiences with computing non-hydrostatic free-surface flows involvingwater waves［J］.International Journal for Numerical Methods in Fluids，2005，48(2):169-197.

［8］ZIJLEMA M，STELLING G S.Efficient computation of surf zone waves using the nonlinear shallow water equations with non-hydrostatic pressure［J］.Coastal Engineering，2008，55:780-790.

［9］AIC，JIN S，LV B.A new fully non-hydrostatic 3D free surface flowmodel forwaterwavemotions［J］.International Journal for Numerical Methods in Fluids，2011，66(11):1354-1370.

［10］STELLING G S，ZIJLEMA M.Numerical modeling of wave propagation，breaking and run-up on a beach［J］.Computational Science and Engineering，2010，71:373-401.

［11］SMIT P，ZIJLEMA M，STELLING G S.Depth-induced wave breaking in a non-hydrostatic，near-shorewavemodel［J］.Coastal Engineering，2013，76:1-16.

［12］BOOIJN，RISR C，HOLTHUIJSEN LH.A third generation wavemodel for coastal regions，part1:model description and validation［J］.Journal of Geophysical Research，1999，104(C4):7649-7666.

［13］HOLTHUIJSEN L H，A HERMAN，BOOIJN.Phase-decoupled refraction-diffraction for spectralwavemodels［J］.Coastal Engineering，2003，49:291-305.

［14］SWAN Team.SWAN Usermanual(version 40.85)［M］.Delft:Delft University of Technology，2011.

［15］GODA Y.A comparative review on the functional forms of directionalwave spectrum［J］.Coastal Engineering，1999，41(1):1-20.

［16］BLAYO E，DEBREU L.Revisiting open boundary conditions from the point of view of characteristic variables［J］.Ocean Modelling，2005，9:231-252.

［17］ZOU Guoliang，ZHANG Qinghe.Improvement of absorbing boundary conditions for non-hydrostatic wave-flow model SWSAH［J］.Applied Mechanics and Materials，2013(353/354/355/356):2676-2682.

［18］邹国良，张庆河.非静压波浪模型无反射造波［J］.海洋工程，2012，30(4):55-61.ZOU Guoliang，ZHANG Qinghe.Wave generation without re-reflection for non-hydrostatic wavemodel［J］.The Ocean Engineering，2012，30(4):55-61.

［19］STELLING G S，DUINMMEIJER SP A.A staggered conservative scheme for every Froude number in rapidly varied shallow water flows［J］.International Journal for Numerical Methods in Fluids，2003，43:1-23.