在船舶与海洋工程领域，基于求解雷诺平均纳维-斯托克斯方程(RANS方程)的数值水槽已被广泛应用于不规则波与结构物相互作用研究．然而该类黏性数值水槽计算精度受多种因素影响[1-3]，例如理论模型、数值方法、网格尺寸和计算时间步长等．因此，如何确保不规则波在长时间模拟中的精度一直是黏性数值水槽研究的关键问题[4]．文献[5]采用摇板造波法在数值水槽内模拟不规则波，在单个有义波高范围设置15个网格，有义波高误差为12%．文献[6]在采用RNG k-ε模型的数值水槽内，通过速度造波入口边界进行造波，分析了网格尺寸、时间步长等对不规则波模拟精度的影响．该研究认为对于不规则波的模拟，须要在单个有义波高范围内设置12～20个网格，时间步长为平均周期的1/2 000，才可得到与目标值符合较好的模拟结果．文献[7]采用RNG k-ε模型对不规则波进行模拟，在自由液面附近垂直方向上的网格尺寸设置为有义波高的1/20，有义波高误差为4%～8%．文献[8]采用不同的湍流模型对不规则波进行模拟研究，发现SST k-ω模型的结果与目标值误差最大．文献[9]也在研究中发现SST k-ω模型得到的不规则波模拟结果与目标值有明显差别．以上研究说明湍流模型对不规则波的模拟精度有较大影响，精度合理的不规则波模拟结果大多采用RNG k-ε模型，而采用SST k-ω模型得到的结果会有明显误差．为了提高波浪的模拟精度，目前已有研究者针对湍流模型提出了改进措施，然而大部分都以规则波为研究对象[10]．文献[11-12]在采用SST k-ω模型模拟大波陡规则波过程中发现湍流黏度增长异常，导致波浪衰减．为此，该研究提出了一种可将波面附近流场限制为层流的浮力修正SST k-ω模型，明显抑制了湍流动能的过度增大．文献[13]在采用SST k-ω模型模拟大波陡规则波过程中也遇到了湍流黏度异常增长的问题．随着波陡增大，水槽内湍流黏度异常增长的范围和幅度越大，波浪衰减越明显．研究采用将流体域划分为层流区和湍流区的方案以抑制湍流黏度的伪增长，然而该方案仅在有限的模拟时长段内起效．文献[14]指出传统的两方程湍流模型用于规则波模拟中，湍流黏度会随时间变化出现指数增长的现象．该研究对多种湍流模型提出了改进的湍流黏度计算公式，明显抑制了湍流黏度的非物理增长，提高了波浪的长时间模拟精度．在黏性数值水槽中，湍流黏度的伪增长会使得不规则波的传播特性发生改变，这也导致船舶与海洋结构物的运动预报发生偏差．虽然文献[14]提出的改进湍流黏度计算公式已被用于波浪与结构物相互作用的研究，然而主要还是以规则波作为研究对象．目前，在不规则波模拟研究中，湍流黏度对模拟精度的影响往往被忽略．因此，研究不规则波模拟中湍流黏度伪增长现象及相应的抑制方法对数值波浪水槽发展有重要意义．为此，本研究分别采用RNG k-ε模型、标准k-ω模型、BSL k-ω模型、SST k-ω模型及SST k-ω湍流黏度改进模型对不同海况的不规则波进行模拟对比，分析了湍流黏度在不规则波模拟中的变化及其对波浪模拟精度的影响．1 理论方法1.1　流体控制方程本研究以单向不规则波为研究对象，流体运动可简化为二维问题．水和空气被视为不可压缩流体，采用流体体积法(VOF)对自由液面进行捕捉．流体运动满足连续性方程、RANS方程和流体体积分数运输方程，控制方程形式为∂ui/∂xi=0；(1)    ∂(ρui)/∂t+∂(ρuiuj)/∂xj=-∂p/∂xi+μ[∂2ui/(∂xj∂xj)]-∂τij/∂xj+Fi+Si;(2)∂α/∂t+∂(αui)/∂xi=0，(3)式中：下标i和j取值1和2分别为水平方向和垂直方向分量；x为位置坐标；t为时间；u为流体速度；p为压力；F为体积力；S为附加源项；τ为雷诺应力；ρ和μ分别为流体的有效密度和有效黏度；α为流体体积分数．本研究采用黏性流求解器ANSYS-FLUENT，分别采用RNG k-ε模型[15]、标准k-ω模型[16]、BSL k-ω模型[17]和SST k-ω模型[17]闭合上述控制方程，湍流模型及参数均为求解器默认设置．此外，针对SST k-ω模型出现的湍流黏度伪增长问题，采用了文献[14]提出的湍流黏度计算改进公式，以分析湍流黏度伪增长抑制的有效性．湍流黏度改进公式为    μt=ρk/{ωmax[1/α*,S*F2/(α1ω),λ2βp0/(β*α2pΩ)]},(4)式中：k为湍流动能；ω为比耗散率；α*=1；α1=0.31；λ2=0.05；β*=0.09；α2=0.52；S*为平均应变率张量．与文献[17]提出的湍流黏度计算公式相比，改进后的计算公式在分母中多了一个应力限制项(λ2 βp0/(β*α2 pΩ))，其中β，p0，pΩ的详细计算公式参见文献[14]．1.2　造波和消波方法不规则波的生成基于线性叠加原理[18]，波高和流体速度表达式分别为η=∑i=1MAicos((Kix1-Ψit)+Φi);(5)    u1=∑i=1M(gAiKi/Ψi)[coshKi((x2+d))/cosh(Kid)]cos((Kix1-Ψit)+Φi);(6)    u2=∑i=1M(gAiKi/Ψi)[sinh(Ki(x2+d))/cosh(Kid)]sin((Kix1-Ψit)+Φi),(7)式中：d为水深；Ai，Ki，Ψi和Φi分别为第i个组成波的波幅、波数、圆频率和初相位；M为组成波的个数，取M=200[19]．成分波的波幅通过海浪谱密度函数(Sw)求得，即Ai=(2SwΔΨi)1/2．为了保证不规则波模拟精度及稳定性，在水槽上游区域内添加造波动量源项，表达式为    Si=aρ(ut,i-uc,i){{exp[(x1,e-x1)/(x1,e-x1,s)]-1}/(e-1)},(8)式中：a为经验系数，取a=6[20]；ut为流场速度理论值，由式(6)和(7)确定；uc为流场速度计算值；x1，s和x1，e分别为造波源项作用区域在水平方向上的起点和终点坐标(x1，s＜x1＜x1，e)．为了防止波浪在水槽末端产生反射，在水槽下游区域内添加消波动量源项，表达式为    S2=-ρ(b1u2+b2u2u2)⋅[(x2-x2,b)/(x2,f-x2,b)]{{exp[(x1-x1,s')/(x1,e'-x1,s')]2-1}/(e-1)},(9)式中：b1和b2分别为线性和非线性阻尼系数，取值分别为10和50[20]；x1,s'和x1,e'分别为消波源作用区域在水平方向上的起点和终点坐标(x1,s'＜x1＜x1,e')；x2，f和x2，b分别为波面和水槽底部在垂直方向上的坐标(x2，b＜x2＜x2，f)．1.3　数值方法本研究采用ANSYS-FLUENT对控制方程进行求解，时间项采用一阶隐式格式，对流项采用二阶迎风格式，扩散项采用中心差分格式，流体体积分数计算采用几何重构算法，速度和压力耦合采用压力耦合方程组的半隐式算法(SIMPLE)实现．此外，通过用户自定义函数功能开发并编译上述改进湍流黏度计算公式、造波和消波动量源项到求解器中．2 结果分析2.1　研究工况及数值水槽本研究海浪谱选用ITTC双参数谱，海况等级为3～6级，各海况的有义波高(H1/3)和平均周期(T0)见表1．所建立的数值水槽的长度为16倍平均波长(λ0)，高度为7 m，水深为5 m，如图1所示．水槽左边界为造波速度入口，右边界为开渠流压力出口，上边界和下边界为固壁．分别将距离入口边界4λ0和出口边界8λ0的区域设为动量源造波区和消波区[20]，水槽中剩下的区域为工作区．在工作区布置三个浪高仪(标记为A号、B号和C号)，用于监测波面变化．浪高仪的位置分别为距离入口边界5λ0，6λ0和7λ0．此外，在工作区中央位置静水面以下0.1 m处，布置一个湍流黏度变化监测点Ptv．10.13245/j.hust.241212.T001表1不规则波海况海况实际尺度模型尺度H1/3 /mT0 /sH1/3 /mT0 /s3级0.8803.9000.0170.5464级2.1005.4000.0410.7565级4.0007.0000.0780.9806级7.0008.7000.1371.21810.13245/j.hust.241212.F001图1数值波浪水槽示意图2.2　网格和时间步长依赖性分析以5级海况不规则波为对象，分析网格布置和时间步长对模拟精度的影响．为了排除湍流模型的影响，不规则波模拟均采用层流模型．数值水槽采用结构化网格布置形式，如图2所示．为了保证波浪传播精度，在造波区和工作区的波面范围内，网格须要沿垂直方向和水平方向进行加密；在波面以外区域，网格垂直方向尺寸以1.2倍增长率朝水槽顶部和底部逐渐变大．此外，消波区内的网格水平方向尺寸以1.05倍增长率朝出口边界逐渐变大，以利用数值耗散加快波浪衰减．在波面加密区内，通过对比单个平均波长内的网格数量(nx)和单个有义波高内的网格数量(nz)对造波精度的影响，确定合理的网格尺寸．考虑到网格单元长宽比的影响，选取不同的nx和nz进行搭配，共生成5套网格方案(见表2)，∆x和∆z分别代表网格水平和垂直方向的尺寸．总模拟时长为250 s，网格依赖性分析中的计算时间步长∆t=T0/1 000．10.13245/j.hust.241212.F002图2数值水槽网格布置示意图10.13245/j.hust.241212.T002表2不同网格尺寸的计算误差对比(Δt=T0/1 000，层流模型)nx×nz(∆x×∆z)/(mm×mm)网格总数耗时/hA号/%B号/%C号/%平均/%H1/3T0H1/3T0H1/3T0H1/3T040×425×2028 33676-4.35.5-4.55.8-5.66.7-4.86.040×625×1332 25689-2.83.0-2.93.0-4.23.4-3.33.140×825×1041 088104-2.63.0-2.73.1-4.13.6-3.13.220×650×1320 26768-5.88.4-5.88.8-7.09.6-6.28.960×616×1050 722120-2.52.9-2.62.9-3.53.6-2.93.1表2汇总了通过上跨零点法得到的不同网格方案的有义波高和平均周期误差．对比表中的数据，可以发现不同位置处有义波高的模拟结果均小于目标值，而平均周期的模拟结果均大于目标值．此外，距离入口边界越远，两个波浪特征参数的误差也会越大．当nx=40时，nz从4增加至6，有义波高和平均周期的平均误差分别减小1%和3%，计算耗时增加17%；当nz继续增加至8时，相较于nz=6的结果，模拟精度没有明显改善，然而计算耗时增加了17%．因此，nz设置为6即可．当nz=6时，nx从20增加至40后，有义波高和平均周期的平均误差分别减小3%和6%，有较为明显的改善．然而当nx=60时，模拟精度相较于nx=40的结果提升较小，计算耗时却增长了35%．综上所述，在兼顾计算效率的同时，采用nx=40和nz=6的网格划分方案能得到较高的模拟精度．选取三种时间步长计算方案分析时间步长对不规则波数值模拟精度的影响，∆t分别为T0/750，T0/1 000和T0/1 250．表3汇总了不同时间步长方案得到的有义波高和平均周期误差，对比表中的数据，可以发现T0/1 000和T0/1 250得到的模拟结果几乎相同，有义波高和平均周期的平均误差相较于∆t=T0/750的结果分别有2%和3%的提升．然而采用T0/1 250的时间步长比采用T0/1 000的计算耗时增加了40%，因此将时间步长设置为T0/1 000即可．图3展示了采用自相关函数法[18]得到的计算谱与目标谱的对比(Sw为谱密度，f为频率)，从图上可以看到：基于nx=40，nz=6和∆t=T0/1 000的设置，在水槽中部位置得到的计算谱与目标谱整体符合较好，谱值仅在1.2～1.8 Hz内略低于目标值．上述分析表明该网格划分方案和时间步长能够满足本研究不规则波模拟精度要求．10.13245/j.hust.241212.T003表3不同时间步长的计算误差对比(nx=40，nz=6，层流模型)Δt/T0耗时/hA号/%B号/%C号/%平均/%H1/3T0H1/3T0H1/3T0H1/3T01/75065-4.95.7-4.95.9-6.07.1-5.36.21/1 00089-2.83.0-2.93.0-4.23.4-3.33.11/1 250124-2.62.8-2.73.0-4.23.4-3.23.110.13245/j.hust.241212.F003图3层流模型得到的计算谱与目标谱对比(B号浪高仪，nx=40，nz=6，∆t=T0/1 000)2.3　湍流黏度影响分析分别采用RNG k-ε模型、标准k-ω模型、BSL k-ω模型、SST k-ω模型及SST k-ω湍流黏度改进模型对表1中四种海况的不规则波进行模拟，分析湍流黏度对不规则波模拟精度的影响．表4统计了采用不同湍流模型得到的工作区内有义波高和平均周期的平均误差，对比表中的数据可以发现：有义波高和平均周期的误差均随海况等级的增大而增大；采用RNG k-ε模型得到的不同海况的模拟结果相差不大，H1/3和T0的误差均在4%以内，而采用k-ω模型得到的不同海况的误差均在10%以上；对比不同k-ω模型的结果，可以发现采用BSL k-ω模型和SST k-ω模型得到的结果非常相近，而采用标准k-ω模型得到的结果误差略低于另外两种k-ω模型结果的误差．值得注意的是，在采用湍流黏度改进模型后，SST k-ω模型得到的H1/3和T0误差降低约3/4．此外，采用SST k-ω湍流黏度改进模型得到的模拟结果略优于采用RNG k-ε模型得到的结果．10.13245/j.hust.241212.T004表4不同湍流模型得到的计算平均误差统计海况RNG k-ε标准k-ωBSL k-ωSST k-ω改进SST k-ωH1/3T0H1/3T0H1/3T0H1/3T0H1/3T03级-3.33.1-10.310.4-11.511.3-11.411.2-2.82.74级-3.53.4-12.410.8-13.611.8-13.212.0-2.92.85级-3.93.6-15.811.5-17.912.6-17.512.3-3.53.36级-3.93.8-16.111.9-18.413.2-17.813.0-3.53.4%图4为四种海况下不同湍流模型得到的计算谱与目标谱的对比．在相同海况下，采用RNG k-ε模型和SST k-ω湍流黏度改进模型得到的计算谱相似，并且在主要能量段与目标谱符合较好，谱峰值的误差均在3%以内．采用三种传统k-ω模型得到的计算谱比较相似，仅在低频区内与目标谱符合，在高频区域与目标谱有较为明显的区别，谱峰值的误差约为10%～20%．通过以上对比可以发现：基于三种传统k-ω模型的不规则波模拟精度均较差，而基于RNG k-ε模型和SST k-ω湍流黏度改进模型的模拟精度较高．10.13245/j.hust.241212.F004图4不同湍流模型得到的计算谱与目标谱对比(B号浪高仪)对不同湍流模型得到的湍流黏度进行对比，分析其对不规则波模拟的影响．图5为采用不同湍流模型模拟在Ptv处湍流黏度变化的时历对比，可以看出：不同湍流模型的湍流黏度在模拟过程的前期并没有明显区别，随着模拟时间的持续增加，三种传统k-ω模型的湍流黏度均出现了过度增长现象，且增幅超出了合理的范围．这是由于在波浪模拟过程中，三种传统k-ω模型在波面下方的近势流区10.13245/j.hust.241212.F005图5不同湍流模型得到的湍流黏度时历对比(Ptv处)(p0≥pΩ)内会变得不稳定，湍流动能呈指数性增长[14]，从而导致了湍流黏度的过度增长．RNG k-ε模型的湍流黏度一般会在波浪破碎后变得不稳定[14]，而本研究在不同海况的模拟过程中均未出现波浪破碎现象，因此该模型相对稳定，其湍流黏度的增幅较小，且整体保持在合理的范围内．而对于SST k-ω湍流黏度改进模型，在湍流黏度计算公式中加入新的应力限制项(λ2 βp0/(β*α2 pΩ))后(式(4))，当波面下方近势流区内湍流动能发生不稳定变化时，该限制项影响权重会增大，从而抑制湍流黏度的变化[14]，使得湍流模型保持稳定．在不同海况下，不同湍流模型的湍流黏度增长情况也不相同．从图5还可以看到：虽然三种传统k-ω模型的湍流黏度的变化趋势相似，但是标准k-ω模型的增长幅度小于另外两种k-ω模型．表5统计了采用不同湍流模型在Ptv处湍流黏度的最大值，对比表中的数据可以发现：随着海况等级的提升，湍流黏度的最大值也在上升；对于BSL k-ω模型和SST k-ω模型，当海况从3级提升至6级时，最大湍流黏度均增大了55 Pa∙s；对于标准k-ω模型，最大湍流黏度增大达到46 Pa∙s；相较于三种传统k-ω模型，RNG k-ε模型湍流黏度最大值的变化幅度较小，不同海况之间仅相差0.7 Pa∙s；与RNG k-ε模型相似，SST k-ω湍流黏度改进模型的最大湍流黏度变化幅度仅为0.5 Pa∙s．10.13245/j.hust.241212.T005表5不同湍流模型得到的最大湍流黏度对比海况RNG k-ε标准k-ωBSL k-ωSST k-ω改进SSTk-ω3级0.060.650.931.020.044级0.292.523.373.310.225级0.568.3014.3714.240.316级0.7346.355.7555.700.53Pa∙s图6为采用不同湍流模型模拟不同海况下，水槽内湍流黏度随时间的变化情况．对于RNG k-ε模型，湍流黏度主要在波面附近发生变化，并且变化幅度相对较小，因此该模型在不同海况下计算的波浪特征参数误差区别不大．对于三种传统k-ω模型，湍流黏度在造波区内的增长更为明显．在3级海况下，造波区内的湍流黏度最大可达10 Pa∙s，已超出合理范围．在6级海况下，造波区内的湍流黏度最大甚至达到100 Pa∙s．此外，三种传统k-ω模型的湍流黏度过度增长的主要影响区域变化相似．随模拟时间的增加，湍流黏度过度增长的主要影响区域会逐步扩大至工作区．此外，主要影响区域的范围也会随海况等级的增大而增大．对于SST k-ω湍流黏度改进模型，湍流黏度主要在消波区内增长，并且增长幅度相对较小，不会影响造波区和工作区的波浪传播．10.13245/j.hust.241212.F006图6不同湍流模型得到的水槽内湍流黏度随时间变化图7为不同湍流模型得到的波高时历对比，从图7可以看到：采用三种传统k-ω模型得到的波高在80 s后出现了明显衰减；采用RNG k-ε模型和SST k-ω湍流黏度改进模型得到的结果几乎重合，由于这两种湍流模型的湍流黏度变化均在合理范围内，波浪能维持长时间稳定传播且无明显衰减，模拟精度较高．综上所述，当采用k-ω模型模拟不规则波时，湍流黏度的过度增长是影响波浪模拟精度的主要因素．湍流黏度改进公式的使用可抑制湍流黏度伪增长，改善波浪模拟精度．10.13245/j.hust.241212.F007图7不同湍流模型得到的波高时历对比(B号浪高仪)3 结论本研究基于RANS方法建立二维数值波浪水槽，分析了湍流黏度对不规则波模拟精度的影响．分别采用RNG k-ε模型、标准k-ω模型、BSL k-ω模型、SST k-ω模型及SST k-ω湍流黏度改进模型对不同海况的不规则波进行数值模拟，分析了湍流黏度变化对波浪模拟精度的影响．本研究在不规则波模拟过程中对k-ω模型的伪增长问题进行了分析，并以SST k-ω模型为例证明了湍流黏度计算改进公式对抑制湍流黏度过度增长的有效性．研究表明这些公式可整合至其他湍流模型中．研究结论总结如下．a．对于RNG k-ε模型，湍流黏度变化区域主要在波面附近，并且随时间变化的增长幅度较小，整体增长幅度在合理范围内，计算波浪谱与目标谱符合较好，波浪特征参数误差在4%以内．b．对于传统k-ω模型，湍流黏度变化区域主要分布在造波区．随着模拟时间的增加，湍流黏度均会过度增长．随着模拟时长和海况等级增加，湍流黏度伪增长区域还会扩大至工作区，导致波浪出现严重衰减，计算波浪谱与目标谱误差明显，波浪特征参数误差均在10%以上．c．文献[14]提出的湍流黏度计算改进公式实施方便，在SST k-ω模型中使用可明显抑制造波区和工作区内的湍流黏度伪增长现象，波浪模拟精度略优于RNG k-ε模型．