黏性数值波浪水槽动边界造波技术研究

蒋学炼，赵 悦，杨伟超，谷汉斌，朱福明

（1.天津城建大学 天津市土木建筑结构防护与加固重点实验室，天津300384；2.宁波大学 土木与环境工程学院，浙江 宁波315211；3.天津市北洋水运水利勘察设计研究院有限公司 水运设计室，天津300452）

随着计算流体力学的发展，数值波浪水槽在海岸及近海工程的研究中得到广泛应用，其关键要求之一是数值造波的稳定性和可靠性.常用的数值造波技术包括边界造波和源函数造波.边界造波即根据目标波形指定波浪水槽入流边界的速度和波面的时间历程[1]，这一方法存在两个问题：一是水槽出流边界或域内建筑物会反射入射波，在造波边界形成二次反射，目前所提出的辐射边界方法只适用于小振幅弱反射波，长时间模拟存在误差累积问题[2]；二是在有限水深中波浪的非线性会导致周期平均净质量流，长时间模拟后会出现计算域总水体体积减小、水位下降问题[3].源函数造波是在流体连续性方程右端增加质量源项，各种波形的质量源项可根据波面方程积分得到，由于入射波和反射波分别从左右边界流出，可有效消除二次反射问题，但相对波高较大时存在波形不稳定的问题[4].

本文模仿物理模型水槽的推板造波方法，设置数值波浪水槽的一端边界为刚性推板，以其运动为扰动源驱使计算域内的水体运动，形成所需的波形.以下各节将分别介绍数值波浪水槽、造波理论和多种波形的实现方法.

1 数值波浪水槽

数值波浪水槽布置如图1.左端为刚性造波板（动边界），底部为固边界，流体与固体交界面均为无滑移边界，法向压力梯度为零，右侧设置海绵消浪层，端部为辐射边界，顶部为自由液面，采用流体体积法（volume of fluid，VOF）追踪.坐标原点x=0置于造波板平衡位置，顺浪向为正，z=0置于静水面，向上为正.

图1 二维数值波浪水槽布置图

数值波浪水槽中的液体运动采用描述不可压缩黏性流动的雷诺时均纳维尔-斯托克斯方程控制（Reynolds-averaged-Navier-Stokes equations，RANS），为了简化固体动边界的处理，避免追踪移动物体而动态剖分网格所消耗的巨大计算资源[5]，采用部分单元法（partial cell technique，PCT）描述网格单元性质，其张量形式为[6]

式中：i（i=1，2）和j（j=1，2）分别代表水平向和垂向；t为时间；ui为i向时均速度；ρ为液体密度；p为液体压力；gi为i向重力加速度；μ为动力黏滞系数；-为雷诺应力项，采用k-ε紊流模型求解；β为部分单元法中定义的开度函数，代表单元中流体面积与总面积之比.通过引入开度函数，整个计算域可分解为流体单元（β=1）、固体单元（β=0）和固液混合单元（0<β<1），单元中的变量取值设为开度函数与单元变量计算值的乘积.

采用交错网格有限差分法配合两步映射法求解控制方程（1）和（2），其中的对流项采用中心差分和迎风格式混合离散，压力梯度项和应力梯度项采用中心差分离散.这一数值模式已应用于众多波浪-建筑物相互作用问题中[7]，详细求解过程可参考文献[8].

2 造波理论基础

在实验室中，以变速装置推动造波板做往复强迫运动（见图2），带动水体产生波浪，其运动-速度传递方程为[9]

式中：XB（t）为t时刻推板位置；uB（t）为推板运动速度；为推板与水体接触界面处的水质点垂线平均水平速度.

图2 推板造波示意图

在数值波浪水槽中，设置左端边界为刚性推板模拟造波板运动（见图1）.在一个微小的时间步长Δt内，可以认为刚性推板的运动速度近似保持不变，因此，动边界与水体交界面单元内的控制方程（1）和（2）可改写为

可以看出，式（3）和式（4）的物理意义一致，均表示推板和水质点交界处两者的运动速度相等，以满足不可压缩流体连续性的要求.式（5）则表明，刚性动边界的运动通过对流和扩散传递给水体，对水体的压力和黏性切应力产生影响，进而引起水体波动.因此，如能计算得到对应不同波形的动边界速度时程uB（t），即可实现数值造波.

3 多波形造波技术及案例分析

3.1 一阶线性波

一阶线性波适用于较深水域（d/Lh≥0.5），其速度势为

式中：A为波浪振幅；ω为圆频率；k1、k3n分别对应行进波和驻波的波数；Cn为系数.

式（6）第一项表示沿x正向传播的行进波，即目标波形.第二项表示驻波系列，是由于造波板与交界面处的水质点运动轨迹不一致而产生的局部扰动，顺波形传播方向呈指数衰减，可估算出一阶驻波的振幅在距离推板2倍水深处衰减96%，3倍水深处衰减99%.因此，数值模拟的试验段应放置在距离推板一定距离以外，以消除驻波的影响.

假定刚性动边界以正弦形式做强迫运动

式中：SB/2为动边界的运动幅值.

将式（6）和式（7）代入式（3），可得

式（8）两边同乘cosh[k1（z+h）]，并沿水深（z=-h～0）积分得到波浪振幅的表达式

将式（6）代入一阶自由液面边界条件[10]，可得波面方程

在远离动边界的位置，式（10）右端第二项的驻波系列将耗散，右端第一项应与一阶线性波的波面方程相等，可得波高为

将式（9）代入式（11），并考虑弥散关系ω2=gk1tanh（k1h），可得波浪要素与动边界运动幅值之间的传递函数

对式（7）求导，可得到动边界的速度时程计算式

设置下列数值案例进行验证：水槽长L=5.0 m，水深h=0.5 m，波高H=0.04 m，周期T=0.7 s，波长Lh=0.76 m，海绵消浪层L2=1.0 m.计算网格取dx=0.02 m，dz=0.01 m，时间步长Δt=0.02 s.采样位置设在x=2.5 m.图3给出了模拟波形与理论波形的对比，其中理论波形由描述.相对误差计算式为，其中ηm为模拟波面高程，ηt为理论波面高程.

图3 一阶线性波模拟波形与理论波形对比

由图3a可以看出，仅采用辐射边界难以完全透射波能，模拟波形有一定程度的雍高.而图3b显示，联合运用海绵消浪层和辐射边界条件有效减小了二次反射.基于人工衰减消波的思想，在海绵消浪层内，控制方程（2）修正为[11]

式（14）右侧最后一项为人工阻尼项，Cd为阻尼系数，表示通过增加海绵层内的附加阻尼力逐渐吸收波能.为保持波形稳定，后续案例的出流边界前均设置了海绵消浪层.

3.2 二阶斯托克斯波

二阶斯托克斯波适用于有限水深区域（0.5＞d/Lh＞0.05），理论波形如下

对应的动边界运动方程为

对时间求导可得动边界的速度时程

图4 二阶斯托克斯波模拟波形与理论波形对比

3.3 孤立波

孤立波是浅水波浪运动的极限状态，波面全部在静水面上，波长无限大，传播过程中波形保持不变，水质点只朝波浪传播方向运动，常用于研究近岸区的波浪破碎和海啸波.基于连续性条件可推导出实验室造孤立波的推板运动-速度传递方程[9]

积分式（18），可得到造孤立波对应的动边界位移方程

Goring和Raichlen[12]建议，物理模型试验中推板的单向运动持续时间由式0.999确定，即，并推荐Newton-Raphson切线法迭代求解式（18）位移时程XB（t），再进行数值微分得到速度时程uB（t）.但在数值模拟时发现，按照t=（0，tf）的时间区间造波，只能出现孤立波的右半部分.为此，本文采用四阶Runge-Kutta法对式（18）在时间区间t=（-tf，tf）进行数值积分，直接获得各时刻的速度时程uB（t），再将时间区间平移至t=（0，2tf）获得完整的孤立波形.

设计两个算例验证孤立波的生成效果：①单孤立波，水深h=0.3 m，波高H1=0.06 m；②双孤立波，水深h=0.3 m，波高H1=0.06 m，波高H2=0.03 m.其他参数取值相同：水槽长L=10.0 m，海绵消浪层L2=2.0 m，计算网格取dx=0.02 m，dz=0.02 m，时间步长Δt=0.02 s.波形的采样位置设在x=6.0 m.图5为数值模拟结果，其中理论波形由式（19）描述.可见，两种孤立波的数值波形与理论波形吻合良好，证明了上述数值处理方法的有效性.

图5 孤立波模拟波形与理论波形对比

4 结论

（1）模拟实验室模型水槽推板造波的方法，在数值波浪水槽中设置刚性动边界，以其强迫运动为扰动源驱使计算域内的水体运动，形成所需的波形.基于实验室造波理论实现了一阶线性波、二阶斯托克斯波和孤立波的数值模拟，与理论波形比较，数值波形吻合良好.

（2）模仿实验室水槽末端设置孔隙介质消减波能的做法，基于人工衰减消波的思想在数值波浪水槽出流边界前增加海绵消浪层，与辐射边界条件结合吸收透射波能，数值结果显示这一处理可有效减小二次反射的影响.

（4）在浅水孤立波的生成中，不同于迭代求解位移时程再进行数值微分得到速度时程的传统做法，本文采用四阶Runge-Kutta法数值积分运动-速度传递方程直接获得速度时程，更为便捷.

致谢：衷心感谢四川大学的林鹏智教授为本文工作提供了RANS-VOF源代码.