多周期稳定数值造波的参数设计方法

吕磊，陈作钢

1 上海交通大学 海洋工程国家重点实验室，上海 200240

2 高新船舶与深海开发装备协同创新中心，上海 200240

3 上海交通大学 船舶海洋与建筑工程学院，上海 200240

0 引 言

在船舶与海洋工程领域，利用计算流体力学（CFD）方法研究波浪中结构物的各类响应一直是重要的研究方向。然而，在进行CFD计算时，数值模拟收敛需要一定的物理时间，但数值波在计算过程中会发生波幅衰减和相位偏差的问题，影响结构物响应的计算精度。因此，实现多周期稳定的数值造波显得尤为重要。

目前，较为常用的数值造波方法为有限体积法（finite volume method，FVM），通过求解雷诺平均（Reynolds averaged Navier-Stokes，RANS）数值方程，并结合流体体积（volume of fluid，VOF）多相流模型来模拟波浪传播过程。对于此类数值造波问题，国内学者进行了广泛研究。例如：杨云涛等[1]基于不可压缩黏性流体的Navier-Stokes方程，采用质量源造波法（mass source-based wave generation）建立了三维无反射数值波浪水池，克服了边界造波法无法消除二次反射波的困难，并成功模拟了波浪与浮体相互作用的问题。王东旭等[2]选择两相不可压缩黏性流求解器 InterFoam作为平台进行了数值波浪水槽的开发，使用动边界造波法和质量源造波法进行了规则波和不规则波的模拟计算，结果表明，所开发的数值波浪水槽具有较高的精度与良好的造波能力。余阳等[3]基于推板造波理论和摇板造波理论，在Open-FOAM平台上采用重叠网格技术建立黏性数值波浪水槽，对比研究了此数值模型分别嵌入层流和湍流两种模型后的计算精度及效率，并对三维数值造波问题进行了进一步研究，验证了在三维两相流体域中求解运动物体与流场交互的可靠性和正确性。上述研究从各方面探讨了数值造波问题，并发现数值波在传播过程中普遍存在波幅衰减和相位偏移的现象。总之，包括网格设置在内的多个参数共同影响了数值造波的精度，而对于数值造波误差的评估及造成误差的关键参数仍缺乏定量的系统性研究。

因此，本文将基于STAR-CCM+商业求解器进行数值模拟，以五阶Stokes波作为研究对象，首先将CFD计算结果与理论解进行对比，通过对一维波动方程的离散形式的分析，得到影响造波精度的几个关键参数，进一步研究这些参数对造波效果的影响。然后，通过算例计算，验证内迭代次数的选择、网格分辨率的设置及库朗数等对波幅衰减和相位偏移的影响，并合理选择参数来显著提高多周期造波精度，最后，总结得出一套合理的参数设计方法，基于此方法，以KCS(KRISO container ship)船型为例，模拟计算其在迎浪中阻力和运动响应，再通过与实验数据的对比，以验证本文参数设计方法的可行性。

1 数学模型

1.1 控制方程

计算的控制方程为Navier-Stokes方程，包括连续性方程和动量守恒方程，如式（1）和式（2）所示。式中：V为控制体；S为控制体的边界；v为流体的速度矢量；ui速度矢量的笛卡尔分量；n为垂直于边界S并指向边界外的单位矢量；t为时间；p为压力；ρ为流体的密度；τij为黏性应力张量的分量；ii和ij分别为xi方向和xj方向的单位矢量；g为体积力，在此即为重力；qi为一个可选的动量源项，将在之后的消波设置中具体说明。

1.2 湍流模型

在对多周期的波浪进行数值模拟时，随着计算物理时间的增加，自由液面处的湍流黏度急剧升高，此时会对波幅产生削减作用，且不同湍流模型对波浪传播产生的负面影响也不尽相同。因此，本文使用标准的低雷诺数K-Epsilon模型，将其中的湍动能和湍流耗散率设为非常低的初始值，这会有助于改善湍流黏度对波浪传播的不利影响[4]。本文所有算例也均采用该K-Epsilon模型，根据文献[4] 湍动能和湍流耗散率的初始值分别设为：1×10−5J/kg和1×10−4m2/s3，与默认值相比分别小了2和3个数量级。

1.3 自由液面捕捉

使用VOF法追踪流体自由液面的位置。该方法使用流体体积比ai表 示第i相流体在网格单元中所占的体积分数。每个网格单元中的体积分数ai之和必须等于1。体积分数的输运方程为

式中，qai为体积分数ai的一个体积源。

1.4 造 波

造波方法采用边界输入法。定义纵向坐标x对应于波浪传播方向，垂向坐标z为铅垂向上的方向。根据五阶Stokes波的解析解，在计算域的边界分别设定水质点x方向（波传播方向）速度ux，z方向速度uz和波面瞬时高度η，如下所示。

式中：ω为圆频率；c=coshkd，其中，d为水深，k为波数；n为Stokes波阶数。其他符号含义见文献[5]。

1.5 消 波

在对波浪传播进行数值模拟时，因计算域范围有限，在各边界处常会发生波浪反射的情况，影响造波精度，因此，在研究时使用了阻尼消波法消除边界处的反射波，以使模拟计算环境更接近于真实情况。该方法在动量守恒方程 (2) 中将动量源项qi定义如下：

式中：fl为线性阻尼系数；f2为二次阻尼系数;uz为垂向坐标z方向的速度分量；eκ−1/e1−1为弯曲项，其中κ为弯曲指数；xsd和xed分别为阻尼层开始和结束时的x坐标（波传播方向），一般xed位于计算域末端，故阻尼消波区的长度为xd=|xed−xsd|；nd控制阻尼消波区的波形衰减形状，一般取nd=2。

为实现理想的消波效果，xd,fl和f2需根据波浪参数进行精确调整[6]。推荐只使用线性消波项（f2=0），阻尼消波区的长度设为2倍波长（xd=2λ），线性阻尼常数f1=Ψ1ω，其中系数Ψ1=π。

1.6 离散格式及其他设置

本文使用有限体积法（finite volume method,FVM）将控制方程离散化。所有数值差分、微分和积分都是基于二阶近似。每个代数方程包含来自网格中心及所有相邻网格中心的未知值，获得的耦合方程组被线性化，并由隐式非定常求解器求解。计算时，使用SIMPLE方法获取压力值并校正速度。压力的松弛因子为0.4，其他所有变量的松弛因子均为0.9。使用高分辨率接触面捕捉（high-resolution interface capturing，HRIC）方法提高自由液面的捕捉精度。同时，激活梯度平滑，减少在大网格长高比下的波面振荡。

2 误差分析

数值模拟的误差来源主要有：模型误差，即数学模型理论解与物理情况实际值间的误差；离散误差，即离散方程的解与数学模型理论解间的误差；迭代误差，即迭代算法的解与离散方程解间的误差。

本研究只考虑实际的数值造波结果与理论解间的差异，故无需考虑模型误差。而通过对线性方程组的求解和利用SIMPLE法进行足够次数的迭代，迭代误差可远远小于离散误差[7]，因而，离散误差是主要误差来源。为便于讨论，本文参考一维波动方程来分析影响离散误差的具体参数。波面控制方程为：

若以显式一阶迎风格式对波面方程(9)进行有限差分，离散误差则为式(10)等号右端项之和。

式中：φ为速度势；c为波速；σ为库朗数（Courant–Friedrichs–Lewy number），σ=cΔt/Δx，其中，Δx为网格尺寸，Δt为时间步长。由式(10)可知，波浪传播的离散误差主要由2个无量纲数影响：

1） 库朗数σ，即单个时间步内波浪传播的网格数。

2）与波浪尺度相比的几何分辨率Δx/λ（其中λ为波长）。当Δx越小，几何分辨率越高，单个网格内的梯度就会越小，离散误差也越小[8]。

离散误差的具体形式会因具体方程及所用的离散格式而不同，但无论哪种形式，误差的大小都与库朗数和几何分辨率有关。对于二维造波问题来说，库朗数定义为：

式中，uz,max为z方向的速度最大值。几何分辨率则与长度 方 向x和高度 方 向z的 尺度比Δx/λ和Δz/H有关（其中H为波高，λ 为波长）。为控制网格几何分辨率的大小，若定义波面附近加密区z方向的网格数为Nz，网格长高比为r，则有Δz=H/Nz，Δx=rΔz=rH/Nz。

当网格尺寸Δz和Δx确定后，可通过改变时间步长Δt来确定库朗数大小，具体如式（12）所示。

综上，后续研究将围绕最小加密区z方向的网格数Nz、网格长高比r以及库朗数σ这些参数展开。造波误差则由波幅误差em和相位偏差ep来衡量[9]。定义波幅的理论解为E=|E|ejθ，数值解为E′=|E′|ejθ（其中θ 为极坐标系旋转角，θˆ为其估计值），则坐标点(xN，zN)(N表示数据点)的波幅百分比误差em和相位偏差ep(单位：rad)为：

在后处理分析误差时，分别对各位置波幅误差em在波周期 (0～10)T，(10～20)T，(20～30)T和(30～40)T之间的时刻进行时间平均，相位偏差ep则主要关注每个位置在t=10，20，30，40T时刻的相位偏差ep。

3 算例设置

3.1 计算域设置

计算模拟二维深水波情况，计算域尺寸根据所对应的波浪波长λ进行调整，计算域如图1所示。为研究波浪在空间上传播可能存在的衰减，计算域长度Lx=10λ。为模拟深水条件，消除上、下边界对流场的影响，计算域高度Lz=6λ，其中水深d=4λ，其他为空气相。坐标原点设在计算域入口边界静水液面高度处，纵向坐标x对应于波浪传播方向，垂向坐标z对应于铅垂向上的方向。

图1 计算域设置示意图Fig.1 Schematic view of computational domain setting

入口边界（x=0）设为速度入口，数值波采用五阶Stokes波形式；出口边界（x=Lx）设为压力出口，并向前设置xd=2λ长度的消波区；底部边界（z=−4λ）设为滑移壁面，对应底部不可穿透条件；顶部边界（z=2λ）设为速度入口，对应于开放的顶部条件。

3.2 网格设置

为了能够精确地捕捉数值波的交界面，需对波面附近网格进行局部加密。这里，选择波高方向网格数Nz和网格长高比r=Δx/Δz作为网格设计策略的参数。如图2所示，在自由液面附近沿波浪传播方向共设置5个网格加密区。

图2 网格加密区示意图Fig.2 Schematic view of the zones by grid refinement

1) 加密区①，高度h1=1.25H，使加密区覆盖整个波高范围。加密区内使用各向异性网格尺寸，高度方向的网格尺寸Δz=1/Nz；长度方向的网格尺寸Δx取Δz的r倍。Nz和r的取值将在下文说明。

2) 加密区②和③的高度h2,h3分别为2.5H和5H，网格尺寸分别为加密区①的2和4倍，并保持长高比r不变，以实现网格尺寸的平滑过渡。

根据线性波的理论解，波面以下z=−0.5λ处质点运动速度小于波面处的5%，z=−λ处质点运动速度小于波面处的0.5%。因此，加密区④和⑤的高度h4，h5分别为0.5λ和λ，同样保持网格长高比r不变，网格尺寸分别为加密区①的8倍和16倍。

3.3 波高监测设置

所有算例模拟计算了各波形从物理时间t=0到t=40T的传播过程。为监测波幅及周期在时间和空间上的变化，在x=1,2,···,8λ的位置设置了8个探针。对于x>8λ的区域，因属于消波区，故不再进行监测。

4 各参数对造波精度的影响

根据上述分析，针对可能影响造波精度的各重要参数进行了研究，参数包括：影响迭代误差的内迭代次数n以及影响离散误差的几何分辨率和库朗数σ。其中，几何分辨率由加密区①波高方向的网格数Nz和网格长高比r来控制，库朗数σ则通过改变时间步长来控制。误差主要关注x=8λ位置的探针在t=40T时刻的波幅误差em和相位偏差ep，即每组算例中em和ep的最大值。作为对照算例A0，其参数设置如表1所示，其他算例则在A0算例的基础上修改相关参数得到。

表1 对照算例A0的参数设置Table 1 Parameter setting of Case A0 for comparison

4.1 内迭代次数

为研究内迭代次数n对迭代误差的影响，在对照算例A0的基础上修改了4个不同的内迭代次数。分别设为5,8,10和20次，对应算例编号为B1～B4，如表2所示。表中，加黑数字表示本组算例中的变量参数（以下同）。误差计算结果如图3所示。

表2 B1～B4算例参数设置Table 2 Parameter setting for Case B1-B4

图3 不同内迭代次数对应的误差变化Fig.3 Error variation corresponding to different inner iteration times

由图3计算结果可知，波幅误差em和相位偏差ep随内迭代次数n的增加而降低；当内迭代次数达到15次以上时，波幅误差em降至3.18%，相位偏差ep降至0.17 rad，这两个误差已达到较低水平，且随着内迭代步数的继续增大，并无明显改善，因此，选择15次作为以下算例的内迭代次数。

4.2 库朗数

在对照算例A0的基础上修改库朗数取值，并研究其对造波精度的影响。设置从0.2～0.5的4个不同库朗数取值，算例编号对应C1～C4，如表3所示。误差计算结果如图4所示。

表3 C1～C4算例参数设置Table 3 Parameter setting for Case C1-C4

图4 不同库朗数对应的误差变化Fig.4 Error variation corresponding to different Courant numbers

由图4计算结果可知，波幅误差em和相位偏差ep都随库朗数的增加而急剧增大。在库朗数为0.5时，数值波波幅与理论值相比衰减了33%，相位偏差可达2.31 rad，对离散误差有极大影响。当波面处的库朗数为0.1时，波幅误差em可控制在5%以内，相位偏差ep约为0.2 rad。继续调低库朗数虽可以进一步降低误差，但却会使时间步过小，计算成本将成倍增加。鉴此，综合考虑计算精度和成本后，本文使用0.1作为后续算例的库朗数较为合理。

4.3 波高方向网格数

波高方向网格数Nz影响沿波高方向网格的几何分辨率，进而对离散误差造成影响。设置4个不同的Nz值算例进行误差影响的验证。具体参数设置和误差结果见表4和图5。

表4 波高方向网格数影响算例参数设置Table 4 Parameter setting of cases with different mesh number in the wave height

图5 波高方向网格数对应的误差变化Fig.5 Error variation corresponding to the mesh number in the wave height direction

由表4和图5所示计算结果可知，随着波高方向网格数Nz的增加，波幅误差em和相位偏差ep都呈明显的下降趋势。当Nz为20时，波幅误差em降至3.5%以下，相位偏差ep降至0.2 rad以下，造波精度在可以接受的范围内。算例D3和D4相较于A0，虽然两种误差都有一定的降低，但过大的网格数会造成计算成本的急剧增大。鉴此，综合考虑计算精度及计算成本，后续计算选择20作为波高方向的网格数值。

4.4 网格长高比

波高方向网格数Nz确定后，网格长高比r决定了波传播方向的网格分辨率，这是影响离散误差的另一个重要因素。在试算中发现，网格长高比r对误差的影响在不同波陡下有很大差异。为研究其具体影响，在保持波长不变的前提下，进行了0.01～0.07共7组不同波陡下的数值计算，每组波陡下完成3～4个不同长高比网格的算例，共计26个。各算例的具体参数设置和误差结果如表5和图6所示。

表5 不同网格长高比算例参数设置Table 5 Parameter setting for cases with different mesh aspect ratios

波幅误差em和相位偏差ep的对比分别在t=20T和t=40T这2个时刻及x=8λ位置处展开。由误差结果分析可知，随着网格长高比r的增加，网格分辨率的降低使得em和ep都逐渐增大，且变化趋势基本相同；不同波陡的波浪对于网格长高比r的敏感性不同，波陡较小的算例（例如，E1和E2两组）在使用较大的r来划分网格时仍能保持较低的误差，而对于大波陡的算例（例如，E6和E7两组），则需要非常小的r才能满足误差要求。鉴此，综合考虑计算精度及计算成本，对于波陡在0.01～0.06之间的波浪，本文通过对图6及图7所示的结果进行拟合可得到网格长高比，具体由下式计算：

图6 em与ep随网格长高比的变化（t=20T，x=8λ）Fig.6 Variation of em and ep with grid aspect ratio at t=20T and x=8λ

图7 em与ep随网格长高比的变化（t= 40T，x=8λ）Fig.7 Variation of em and ep with grid aspect ratio at t= 40T and x=8λ

根据式（15）计算的r可以满足：t=20T时刻的波幅误差em约5%，相位偏差ep约0.2 rad；t=40T时刻波幅误差em约10%, 相位偏差ep约0.4 rad。

而对于波陡H/λ≥0.06的波浪而言，即使使用网格长高比为1的网格进行计算，仍然不能满足上述误差要求。对于此类大波陡波浪，可能需要改进自由液面捕捉方法或者使用其他创新方法等才能够得到较为理想的结果。

4.5 参数设置的总结

经过上述数值计算及结果对比可知，为实现多周期稳定数值造波，推荐采用以下设置：

1） 选用标准低雷诺数K-Epsilon湍流模型，湍动能和湍流耗散率的初始值分别为1×10−5J/kg和1×10−4m2/s3，以降低湍流黏度对波面衰减的影响；

2） 在计算域下游边界设置阻尼消波区，阻尼消波区的长度设为两倍波长（xd=2λ）；仅使用线性消波项（f2=0），线性阻尼常数f1=Ψ1ω；

3） 使用HRIC方法，并激活梯度平滑选项；

4） 内迭代次数设为15次或以上，自由液面处最大库朗数维持在0.1以下，并据此调整时间步长；

5） 对自由液面附近进行网格加密，加密区①波高方向的网格数Nz设为20，网格长高比r使用式(15)计算，其他加密区网格尺寸参考4.3节。

计算结果表明，使用上述参数设置方法，可以实现多周期的稳定数值造波，波幅误差em和相位偏差ep满足工程应用的精度要求。图8所示为算例E4.4在x=8λ位置波高随时间的变化曲线，图9所示为t=40T时刻波高随空间位置的变化曲线，其中波高、时刻和位置都进行了无量纲化处理。由图8和图9可以看出，实际波幅线（黑色）与理论波幅线（红色）基本重叠，实际波幅与理论波幅在时间上吻合较好，空间上波形分布均匀，但随着传播距离的增加存在一定的衰减，可见，消波区的波形消波效果良好。

图8 算例E4.4在x=8λ位置波高随时间变化Fig.8 Variation of wave height with time at x=8λ in Case E4.4

图9 算例E4.4在t=40T时刻波高随空间位置的变化Fig.9 Variation of wave height with spatial position at t=40T in Case E4.4

5 数值造波方法应用

以上算例都是针对模型尺度下的二维数值造波研究，为验证本文所提参数设计方法在三维数值造波问题上的适用性，计算了KCS船型在迎浪中的阻力和运动响应，并与工程流体动力学（EFD）试验值[10]进行对比。

为与实验条件保持一致，波长设为λ=0.5LBP=115 m，波陡H/λ=1/60，弗劳德数Fr=0.26。由于船体左右对称，为减少计算量，仅对半船的流场进行模拟，并使用重叠网格方法及DFBI (displaying dynamic fluid body interaction)模型实现船体迎浪时的升沉与纵倾运动。

背景区域及重叠域网格采用切割体非结构网格形式，其中，背景区域网格基准尺寸取为4.0H，重叠区域网格基准尺寸取为0.2H，并在背景区域设置多层嵌套的加密区以保证体网格尺寸的过渡平滑；船体边界层采用棱柱层网格来捕捉，边界层共划分有16层，增长率为1.2。自由液面处则根据4.3节所述参数进行加密，y轴方向网格大小取为4.0H。最终，得到的背景区域体网格总数为430万，重叠区域体网格总数为1 100万，图10所示为中纵剖面船体附近的网格示意图。

图10 船体附近网格Fig.10 Grids near the hull

总共计算20个波周期，在中纵剖面船体前方一倍LBP处设置探针，监测得到的波高变化曲线如图11所示，由图11可知，在整个计算时间内来波的波幅都无明显衰减。计算结果主要关注波浪增阻系数σaw、升沉和纵倾的传递函数（transfer function，TF），其计算公式分别如下：

图11 船前一倍LBP处波高的变化Fig.11 Variation of wave height at a distance of onefold-LBP ahead the ship

式中：Fx,wave和Fx,calm为船体在波浪和静水中的纵向受力的平均值；BWL为水线面宽度；LBP为垂线间长；x31，x51和ζI1分别为升沉、纵摇和波高瞬时幅值的一阶傅里叶系数；波数k=2π/λ。其中，升沉和纵摇以船体重心为参考点。

图12所示为自由液面在初始和计算一段时间后的自由液面波形对比图，由自由液面波高图可清晰地捕捉到开尔文船行波。增阻系数、升沉和纵倾的传递函数的计算结果见表6。

图12 自由液面波形对比Fig.12 Comparison of wave pattern for free surface

表6 CFD与EFD计算结果对比Table 6 Comparison of calculation results by CFD and EFD

由表6的的计算结果可知，CFD计算结果与EFD结果相比均偏低，误差约为4%～10%，这进一步验证了本文所提造波参数设计方法可以推广到三维数值造波的相关问题之中。

6 结 论

本文分析了一维波动方程的离散形式，得到了影响数值造波精度的几个关键参数及其对造波效果的具体影响。通过设置不同参数对深水二维水池进行了波浪模拟，对每组算例的计算误差进行评估，综合考虑计算精度和成本，确定了各参数的最佳值。而且，提出了包括湍流模型、消波方法、内迭代次数、库朗数及网格设置的一整套参数设计方法。通过分析计算，得到如下结论：

1） 数值造波误差主要包括离散误差和迭代误差。离散误差是影响数值造波精度的主要误差因素，与之相比的迭代误差可通过增加内迭代次数来降低到远小于离散误差的范围；

2） 影响离散误差大小的主要因素为库朗数和网格分辨率。其中，库朗数对误差影响较大，需确保自由液面库朗数在0.1以下，网格分辨率则需通过调整波高方向网格数及网格长高比来实现；

3） 波陡H/λ对数值造波的结果影响很大。在H/λ≤0.06的波陡范围内，使用4.5节参数设计方法可使得20个波周期和40个波周期内波幅误差分别约为5%和10%，相位偏差分别约为0.2和0.4 rad，并可推广到三维造波问题中；

4） 对于五阶Stokes波以外的其他规则的非线性波或线性波，计算结果表明，上述参数设计方法均可一定程度上减少造波误差，使造波结果满足精度要求；

5） 对于波陡H/λ≥0.06的大波陡波浪，即使使用网格长高比r=1的网格计算仍有较大误差，若为提高计算精度而进一步增加网格分辨率，则计算成本将变得过高，不具有实际应用的可行性。对于此类大波陡数值造波问题，仍需通过进一步研究和创新方法来解决。