规则波与流相互作用的数值模拟与不确定度分析

姚顺， 马宁， 丁俊杰， 顾解忡

(1.上海交通大学 海洋工程国家重点实验室，上海 200240；2. 上海交通大学 船舶海洋与建筑工程学院，上海 200240)

在实际海洋环境中，波流相互作用是一种非常普遍的现象，流的存在会使波浪的波高、波长等运动学特性发生显著变化。依照波浪理论，叠加规则波可以得到不规则波等实际波浪，可见研究规则波与流相互作用可为实际波浪与流相互作用研究打下基础，对模拟实际海域具有重要意义。

研究人员通常结合数值模拟和试验的方法对波浪及波流相互作用问题进行研究。胡捍红[1]选择四阶非线性薛定谔方程模拟了过渡水波在传播过程的非线性演变过程；贾岩等[2]基于波浪模式SWAN和海流模式POM研究了台风中的波流相互作用；Gao等[3]和Liu等[4]分别利用商业软件开发数值波浪水池，研究了规则波与水平圆柱相互作用问题。通过试验能够直观地验证模拟结果的准确性，Fernando等[5]、Qi等[6]开展试验研究了波浪与顺逆流的非线性相互作用；Soltanpour等[7]和Wei等[8]通过试验分别研究了波流相互作用下波浪在泥床上的耗散和桥梁基础的水动力响应特性；Jiang等[9]利用CFD软件模拟了不规则波与淹没式防波堤相互作用，并将模拟结果与实验进行了对比。然而，对比模拟结果与试验结果验证数值模型的有效性从本质上说是定性的。在第22届国际船模拖曳水池会议(ITTC)阻力委员会上推出了定量评估数值模拟结果的具体方法，即不确定度分析方法[10]，至此相关学者开始关注如何定量评估计算流体力学数值模拟结果的可靠性。通常，数值模型的数值误差只能估计得到，Richardson外推法和网格收敛指数[11]因此被广泛应用于估计误差和数值不确定性。基于此Wu等[12]对船模阻力进行了不确定度分析，通过分析结果发现有效控制近船模表面第1层网格高度并选择合适湍流模型可大幅度降低船模阻力数值试验的不确定度。Zhu等[13]和Deng等[14]分别对船模横摇运动和小水线面双体船在波浪中的纵向运动进行了数值模拟与不确定度分析，给出了减小船舶运动参数数值模拟不确定度的有效建议。Silva等[15]、柏君励等[16]则对规则波和聚焦波的波参数进行了不确定分析。然而，文献[1-16]没有研究流对波浪模型不确定度的影响，这一问题的难点在于开展波流相互作用试验时很难保证造波精度。随着丁俊杰等[17]对上海交通大学风洞循环水槽消波装置的改进，该水槽能很好的满足试验要求。本文对均匀流作用下规则波的生成与演化进行了物理试验，然后基于RANS方程建立数值波浪水池进行波流相互作用下规则波特性研究，并依照ITTC规程开展计算结果的不确定度分析。给出了规则波的运动学特性随流速变化的一般规律，并讨论了顺流对波流相互作用数值模型不确定度的影响情况。

1 数值模型及验证

1.1 试验安排

本文的试验在上海交通大学风洞循环水槽中进行，该循环水槽装置的工作原理如图1所示。装置整体长24.6 m，宽4.5 m，高8.5 m，主要分为工作段、上游整流段、下游整流段、底部回流段以及动力段等部分。其中动力段能使循环水槽中的水体沿顺时针方向进行循环流动。工作段长为8.0 m，宽为3.0 m，高1.95 m，工作段的水深根据循环水槽工作状态的不同可小幅度调节，水深范围为1.60～1.63 m。安装于工作段前方的摇摆式造波机，可以生成规则波、JONSWAP谱不规则波与聚焦波等多类型波浪。实际造波时，根据试验需要设置造波机角度振幅和周期，即可造出目标波浪。当在静水中造波时水槽水深固定为1.63 m，而单独造流或同时造波和造流时水深固定为1.6 m。循环水槽能够实现长时间稳定的造流，稳定造流的流速范围为0.1～3.0 m/s。

试验的目的是获取无流及有流情况下规则波自由液面处波高时历曲线，验证波流相互作用数值模型有效性，并分析流对规则波运动学特性的影响。试验布置如图2所示，其中消波装置距造波端6 m。为满足试验需要，在水槽工作段共布置5个浪高仪，1#、2#和3#浪高仪布置在消波装置前侧，分别距造波端3.2、3.6、4 m。为了减小前方浪高仪对后方浪高仪的影响，1#、2#和3#浪高仪分别距水槽左端1、2、1.5 m；4#和5#浪高仪平行布置在消波装置后侧，距造波端6.7 m，且分别距左侧槽壁1、2 m，用于测试消波装置的性能。实验开始时，造波机生成的前5～10个波并不稳定，记录数据时将前10 s时间段内的波浪忽略。所有浪高仪量程为30 cm，测量频率为100 Hz，精度为0.5%。

图2 试验布置示意Fig.2 The arrangement of wave experiments

1.2 数值模型

本文利用计算流体力学商业软件STARCCM+建立数值模型[18]，控制方程为Navier-Stokes方程：

(1)

式中：ui为xi方向的速度分量；ρ为流体的密度；p为压强；μ为动力粘度；gi=g是重力加速度。本文关注湍流引起的流场在时间上的平均变化，主要采用雷诺平均的Navier-Stokes(RANS)方程：

(2)

(3)

式中：μt是湍流黏度；τij是Kronecker函数；k是湍动能。可采用Renormalization-Group(RNG)k-ε模型求解湍流黏度和湍动能，从而得到雷诺应力。

本文依照上海交通大学风洞循环水槽构建二维数值波浪水池。数值水池如图3所示，主要分为工作段与消波段。由于只考虑波浪沿x方向的传播，数值水池的宽设为0.02 m，小于最小网格y方向的尺寸，采样点位置与试验中浪高仪位置对应。

图3 数值波浪水池Fig.3 Numerical wave flume

网格划分时自由液面作为液相与气相交界面，造波过程中会发生明显的起伏运动。通常对自由液面处的网格进行多层加密，从而减少造波误差。计算域x、y、z方向对应的最小网格尺寸分别为0.01、0.1、0.002 m。主要采集x=4 m自由液面处波高信息作为输出数据，采样频率为100 Hz。

造波方法采用STARCCM+中造波模块的源项造波法；消波时，在尾部2 m添加阻尼项模拟阻尼消波段，并在阻尼消波段前加入一种多层变角度开孔折弯板透水消波装置[19]一同构成消波区。该消波板能够保证较高的透水效率的同时实现高效的消波。

采用STARCCM+库函数CURRENT设定水体速度造流。数值模型空间离散采用二阶迎风格式，压力和速度的耦合求解采用SIMPLE算法。自由液面捕捉采用流体体积(VOF)法。入口AC边界选择速度入口，出口BD边界为压力出口，底部AB边界为无滑移固壁边界而顶部边界CD为压力入口边界，两侧的边界选用对称边界。

按照Stokes波浪理论生成规则波，Stokes波浪理论与线性波类似，波浪运动也是势运动。对于弱非线性问题，可用摄动法进行求解，先假设势函数φ和波面曲线函数η都是某个小参数εs的幂级数，则：

(4)

(5)

当n=2，5时，代入自由表面条件泰勒展开式中，即可得到弱非线性的二阶和五阶Stokes波面方程和势函数。本文数值模拟的研究对象为Stokes五阶波，在软件中选择波浪类型并输入指定波高和周期，即可得到目标波浪。试验及模拟工况如表1所示，其中C为流速，H为波高，L为波长，T为周期，EXP和NUM分别表示对该工况进行试验或者数值模拟。每个工况重复3组，每组试验或者模拟时长为50 s。

表1 规则波试验与数值模拟工况

1.3 数值模型验证

图4 工况R1数值模拟结果的验证Fig.4 Verification of simulation results of case R1

表2 数值模拟与试验数据的平均波高对比Table 2 Experimental and numerical results of wave height

2 结果及不确定度分析

2.1 结果分析

图5给出了无流及顺流情况下规则波自由表面分布情况。从图5可知，C=0 m/s时规则波的自由表面比较光滑；而在顺流(C=0.3 m/s)作用下，规则波波浪自由表面变得粗糙。

数值模拟得到的工况R1～R5波高时历曲线如图6所示。图6表明，随着均匀顺流流速增加，规则波波高降低，波峰和波谷相比于C=0 m/s时更平坦，逆流对规则波的影响与顺流工况的结果相反。

图5 数值水池中波浪自由表面情况Fig.5 The water surface of waves in the numerical wave flume

图6 波高时历数值模拟结果Fig.6 Simulation results of the water surface elevations

从图6也发现当逆流流速达到-0.3 m/s时，规则波传播变得缓慢，落后于无流及顺流工况；逆流流速过大时(C=-0.6 m/s)，规则波波高几乎为0，这是由于当逆流流速过大时，规则波在传播时会在某个位置被阻隔。图7给出了C=-0.6 m/s时波浪被阻隔的情况，被阻隔后规则波波面趋于平稳。

图7 波浪阻隔现象Fig.7 The wave blocking

2.2 不确定度分析

之前对数值模型的验证从本质上来说是定性的，无法定量评估数值模拟的误差。本节拟采用ITTC推荐的不确定度分析方法，定量评估数值模拟的误差。ITTC推荐的CFD数值模型的不确定度分析包括验证和确认2个部分，其中评估数值不确定度的过程称为验证，评估模型不确定度的过程称为确认，具体过程在ITTC规程中有详细介绍[6]。

评估数值不确定度时，需要对数值模型的参数进行收敛性分析。通常考虑的参数包括迭代次数、网格尺寸、时间步长和其他参数，本文规则波与流相互作用属于非稳态问题，网格尺寸和时间步长引起的误差比迭代次数等其他参数引起的误差要大，因此本节的不确定度分析将网格尺寸和时间步长作为主要因素加以考虑。

首先对工况R1中的网格参数进行收敛性研究，时间步长为0.001 s。用3套网格分别进行数值模拟后，得到3组波高时历曲线，如图8所示。

图8 不同网格对应的波高时历曲线Fig.8 Water surface elevations simulated with different grids

相邻网格参数波高时历之差的平均L2范数为：

(6)

(7)

其中ηG1i、ηG2i、ηG3i分别表示在细网格、中等网格和粗网格上模拟得到的波高时间序列中的第i个波高值，全局收敛因子为：

(8)

(9)

(10)

网格收敛因子：

(11)

式中pG,est是极限阶数的估计值，一般取2。由网格尺寸导致的数值模拟误差为：

(12)

当CG远小于或远大于1时，网格尺寸引起的数值模拟不确定度为：

(13)

当|1-CG|≥0.25，网格尺寸引起的数值模拟不确定度为：

(14)

接下来在细网格Grid1上进行时间步长收敛性分析，得到3个时间步长对应的各个波高时历曲线，如图9所示。

图9 不同时间步长对应的波高时历Fig.9 Water surface elevations simulated with different time steps

相邻时间步长参数的波高时历之差的平均L2范数为：

(15)

(16)

其中ηT1i、ηT2i、ηT3i分别表示时间步长0.001 s、0.001 41 s、0.002 s时数值模拟得到的波高时间序列中的第i个波高值,故全局收敛因子：

(17)

(18)

(19)

pT,est同样取2，时间步长修正因子：

(20)

由时间步长导致的数值模拟误差为：

(21)

时间步长导致的数值模拟不确定度：

(22)

综上，在x=4 m自由液面处t∈[15,25]s时间段内波高时历的数值模拟不确定度为：

(23)

数值模拟误差估计值为：

(24)

表3 数值模拟验证结果Table 3 Numerical results of verification

表3给出，R1和R22个工况数值模拟的网格不确定度都大于时间步长不确定度，可知数值模拟无流及顺流情况下的规则波运动时，数值模拟结果对网格尺寸依赖程度大于时间步长，实际模拟时计算资源有限情况下应尽量满足网格精度要求。对比R1和R2工况数据可知，顺流工况R2数值模拟的网格不确定度是无流工况R1的3.46倍，时间步长不确定度是无流时的1.60倍，可见相比于无流工况R1，顺流提高了规则波时历曲线数值模拟结果对网格尺寸和时间步长的依赖程度。并且，顺流时规则波的平均波高小于无流工况，考虑顺流与规则波相互作用问题时，可适当减少自由液面加密区的高度，但需要注意选择合适的网格尺寸和时间步长，减小数值模拟带来的误差。

接下来根据ITTC规程，对数值模拟结果进行确认，工况R1中x=4 m自由液面处模拟的波高时历曲线确认不确定度UV为：

(25)

(26)

式中ηEXP1,x=4 m和ηNUM1,x=4 m分别表示工况R1在x=4 m自由液面处波高时历曲线的试验和模拟值。

表4 数值模拟确认结果Table 4 Numerical results of validation

3 结论

1)均匀顺流使规则波波高降低，传播速度加快，波峰和波谷与无流工况相比更平坦；逆流对规则波的影响与顺流的影响相反，逆流流速过大时，规则波会被阻隔。

2)无流和顺流工况规则波波高时历曲线的比较误差都小于对应的确认不确定度，数值模拟结果得到确认，构建的数值波浪水池具有可靠性。

3)无流及顺流情况下规则波的数值模拟结果对网格尺寸依赖程度大于时间步长，实际模拟时计算资源有限情况下应尽量满足网格精度要求。

4)相比无流工况，顺流加剧了规则波波高时历曲线计算结果对网格尺寸和时间步长的依赖程度。

目前在循环水槽中进行逆流与规则波相互作用的试验时很难控制造波精度，下一步将从改善造波板水力传递函数等方面提高逆流中的造波质量，开展逆流与规则波相互作用的不确定度分析。