速度边界造波法与推板造波法的差异化分析

2021-12-30程友良焦慎俐

海洋技术学报 2021年5期

程友良，许强，焦慎俐

（ 1.华北电力大学动力工程系，河北保定 071003；2.华北电力大学河北省低碳高效发电技术重点实验室，河北保定 071003； 3.华北电力大学保定市低碳高效发电技术重点实验室，河北保定 071003）

数值波浪水池是用于模拟自由表面波、水动力和浮体运动的数值工具的总称[1]。在海洋工程领域，物理模型和数值模拟是两种主要的实验方法。物理模型可以在实验室中准确的再现波的传播和变形过程，但这种方法成本较高，而且受到实验设备设置的限制。随着计算机技术的飞速发展，数值模拟已成为研究波浪及其对海洋结构物作用的有效工具。数值模拟的成本更低，且能够得到更多难以通过物理实验测量的细节[2]。BOO S Y[3]将线性规则波和不规则波的仿真结果与理论值进行了比较，线性波仿真结果与理论值吻合良好，但对非线性波的模拟要进一步研究。RYU S等[4]为研究非线性波浪开发了一种具有完全非线性自由表面边界条件的数值波浪水池。此后，数值波浪水池的优化一直伴随着数值波浪水池的发展。CHEN Y P 等[5]提出了一种基于变换的N-S方程求解的全非线性数值波浪水池，数值计算结果与线性或非线性解析解吻合较好，验证了该数值波浪水池的准确性。ZHANG X T等[6]研究了在完全非线性数值波浪水池中的波传播问题，解决了每个时间步长的潜在流动边界值问题。随着数值波浪水池可用性以及功能性快速增长，其精确性越来越重要。BUNNIK T等[7]对基于势流理论和VOF（Volume of Fluid）理论的数值波浪水池准确性进行评估。LIU Z等[8]建立了基于Fluent的不可压缩粘性流体两相VOF数值波浪水池，以模拟波浪的产生和传播，且波浪运动诱发了气相的运动。LAL A等[9]使用ANSYS CFX模拟了襟翼类型造波机，并对各种湍流模型进行了比较，模拟采用层流模型、k-ε模型和SST（Shear-Stress Transport）切应力输运模型，结果显示这3种方法没有显著差异。GOMES M N等[10]采用了基于VOF模型的二维数值波浪水池，模拟了两种不同形式的规则重力波，表明VOF是一种能够正确模拟水与空气之间相互作用的模型。无网格方法也适用于数值波浪水池。SENTURK U等[11]研究了在自由表面上受到非线性边界条件的二维数值波浪水池数值模拟，采用局部无网格RBF（Radial Basis Function）方法。这项研究表明了无网格方法的灵活性，特别是对于移动边界问题，其中求解过程中的每个时间步仅需要更新节点位置而不是重新整理整个域。随着计算机能力的进步，波浪模拟精确性的提高，好的消波方式能进一步提高数值波浪水池的性能。KIM M W等[12]评估了用于渐进水波的各种人工阻尼方案。将五种类型的阻尼方案应用于计算域的末端，并比较了它们的阻尼能力，提出了在一定范围的阻尼区域上与适当的阻尼形状相关联的最优数值阻尼方案。但反射通常是无法避免的，PERIC R等[13]提出了一种通过数值流动模拟来量化这种影响，使数值波浪水池的波反射最小化。3年后，PERIC R等[14]进一步提出了一种解析预测反射系数的理论，该理论可用于在进行仿真之前优化选择源项参数，与规则自由表面波的有限体积流动模拟结果进行了对比验证，结果表明，该方法具有令人满意的实用精度。随着造波与消波方法的优化发展，造波的其他方面也引起人们重视。GALERA-CALERO L等[15]对采用VOF的k-ε、k-ω和LES（Large Eddy Simulation）模型分别进行了分析和比较，在物理波浪水池上进行了实测，并与理论进行了比较，验证了各种湍流模型的正确性。UDDIN M N等[16]提出了一种二维数值波浪水池来计算实验波浪水槽下壁的静压变化，得出对于小型二维实验波浪水槽，非定常无粘流体法可以较准确地预测表面压力分布。

数值波浪水池模型有多种，但国内对不同种类造波方法的讨论极少，本文对速度边界造波法和推板造波法的波幅历时曲线、历程曲线、压力分布图和u、v速度分布图进行了分析，得出的结论对海洋工程具有重要意义。

1 波浪要素

1.1 求解波浪要素

设周期T=1.94 s，静水深d=1 m，波高H=0.1 m。由弥散关系对波长进行迭代求解，利用MATLAB的求解程序求得波长L=5 m。

1.2 验证波浪要素

1.2.1 二阶stokes波的理论验证

根据梅沃特图[17]采用两个无因次独立参数。

该波浪要素符合理论要求。

1.2.2 模拟条件验证

（1）二阶stokes波只适用于中、深水深的规则波。对于浅水波，；对于中、深水波，。该波浪要素d/L=1/5＞1/20，符合条件。

（2）波函数方程中的二次项必须要小于线性项，从而确保收敛，没有波浪破碎。要满足以上条件必须要满足式（4）和式（5）[18]。

式中，A为振幅；k为波数；Ur为努塞尔数。

该波浪要素Ur=0.005＜1/3，A=0.05＜0.26，符合条件。

（3）波浪破碎极限

当波浪破碎时，没有任何波浪理论是有效的。最有效的波动特征值在水深较大时为波陡H/L，在水深较小时为相对波高H/d。在深水中波浪的破碎极限如下。

式中，H0为深水波波高；L0为深水波波长。

该波浪要素未超过破碎极限，符合Stokes二阶波。

2 速度边界造波法

2.1 二阶Stokes波

流体运动时，必须遵循的另一规律为动量守恒定理。在研究液体波浪运动时，一般认为流体是无粘性的理想流体，此时可以忽略流体的粘性效应，只考虑压强而不必考虑切向应力。

Stokes二阶波的速度势函数见式（7）。

Stokes二阶波的波面方程见式（8）。

Stokes波的水质点运动速度见式（9）和式（10）。

式中，u、v分别为x、y方向的水质点速度；为波频。

2.2 控制方程

连续性方程见式（11）。

Navier-Stokes方程（N-S方程）见式（12）至式（14）。

式中，ρ为液体密度（kg/m3）；g为重力加速（m/s）；t为时间（s）；p为压强（Pa）；μ为动力粘度（Pa·s）。

2.3 VOF模型

2.3.1 VOF机理

VOF模型中界面跟踪是通过求解相连续方程实现的，界面的位置是通过求出体积分量中急剧变化的点来确定的。VOF模型可以对两种或更多种不相溶的流体进行建模。典型的应用包括射流破裂的预测，液体中大气泡的运动，溃坝后液体的运动以及任何液—气界面的稳定或瞬态跟踪。

VOF模型中，不同的流体组份共用着一套动量方程，通过引进相体积分数这一变量，实现对每一个计算单元相界面的追踪。在每个控制体积内，所有相体积分数总和为1。因此，在任何确定单元中的变量和属性要么纯粹代表其中一个相，要么代表各相混合，这取决于体积分数值。故在单元中，若第q相流体体积分数为αq，那么可能存在以下三种情况：

（1）αq=0：单元里不存在第q相流体。

（2）αq=1：单元里充满了第q相流体。

（3）0＜αq＜1：单元里包含了第q相流体和一相或者其他多相流体的界面。

基于αq的值，为计算域中每个计算单元分配合适的属性与变量。

2.3.2 VOF方程

跟踪相之间的界面是通过求解一相或者多相的容积比率的连续方程来完成的。对第q相，有

式中，S为质量源项；V为速度矢量。在默认情况下该方程右端源项为0，但当你给每一项制定常数或用户定义的质量源，则右端不为0。初相的体积分数的计算基于公式（16）。

2.3.3 VOF使用限制

（1）VOF只能应用于压力求解器，在密度基求解器中无法使用。

（2）每一计算单元必须充满一种或多种流体，VOF模型不允许区域中不存在任何流体的情况。

（3）仅有一相可定义为可压缩理想气体，但对于使用UDF定义的可压缩液体则无限制。

（4）流向周期运动（指定质量流率或指定压力降）无法与VOF一起使用。

（5）二阶隐式时间步格式无法与VOF显示格式一起使用。

（6）当利用并行计算进行粒子追踪时，在共享内存选项被激活的情况下，DPM模型无法与VOF模型一起使用[19]。

2.4 消波

2.4.1 阻尼消波

在消波区的动量方程添加阻尼源项。

式中，μ（x）为消波系数，根据公式（19）确定。

式中，a为常数，根据经验选取；x0为消波区起始位置。

2.4.2 拉伸网格

在消波区增大网格尺寸，通过计算过程中的数值耗散来达到消波的目的。

2.4.3 混合消波

阻尼消波与拉伸网格方法相结合的混合网格消波方法。本文中所用消波方法都为此法。

2.5 模型建立

数值波浪水池如图1所示，长40 m，高2 m，静水深为1 m，20～40 m处为消波区。

图1 速度边界造波水池模型

2.5.1 网格划分

网格划分为表1的3种方式，网格划分方法参考了相关文献[20-21]。Δx为x方向单个单元格长度；加密处Δy为静水面±H区域内y方向单个单元格长度，x方向不变；非加密区Δy为加密区之外区域y方向单个单元格长度，x方向不变。此外，消波区网格在x方向数量取150，最小为Δx，并以1.1的网格尺寸比率增加。

表1 网格划分方式

网格的全局示意图如图2所示。为了降低压力传递过程中的数值耗散，需要对波面附近的网格进行加密，加密区（图3）。

图2 全局网格

图3 局部网格

2.5.2 模拟设置

该模型采用表2所示的设置。

表2 模型设置

在消波区添加编译好的阻尼源项，和拉伸后的网格一起达到消波的目的。文献[18]中提到使用层流模型与常用的模拟结果没有显著差异，故选用层流模型。在距进口5 m、10 m和15 m的波面处设置浪高仪，监测波幅的变化。

2.5.3 独立性验证

（1）网格独立性验证

在进行数值计算过程中，为了降低网格密度对计算结果的影响，通常要对多套疏密程度不同的网格系统进行计算，并比较不同网格系统下的计算结果。

从图4可得出，3种网格模拟的波幅随时间变化的趋势相同，随着模拟时间的增加，模拟结果与理论波幅的误差逐渐增大。这是因为在消波区的混合消波方法虽然起到很好的消波效果，但并不能完全消除反射。故随着模拟时间的增加，波浪的反射会越来越大，影响模拟的正确性。因此为了确保结果的准确性，需要控制模拟波浪的时间。经过多次模拟，该数值波浪水池模型下的模拟时间控制在30～50 s。由图5可得波幅波峰和波谷的变化趋势。

图4 波幅历时曲线

图5（a）中3种网格下模拟的波幅与理论波幅在波峰相一致。在波谷：Mesh3与理论波幅近乎一致；Mesh2与Mesh3在10～16 s相一致，但在16～20 s会与Mesh3拉开距离。相比其他两种网格，Mesh1下波谷处一直与理论结果相差较大。

图5（b）中3种网格下模拟的波幅与理论波幅较符合。波峰：Mesh1与Mesh2下的波峰的波幅与理论波幅的波峰相一致；Mesh3下波峰的波幅比理论波幅的波峰要高。波谷：Mesh3下的波谷的波幅与理论波幅几乎一致，Mesh2与Mesh1下波谷的波幅比理论波幅的波谷高，且Mesh1下最高。

图5 波峰与波谷的局部对比图

故Mesh1的结果与理论波幅相差较大，Mesh2与Mesh3较为接近。但Mesh2在波峰表现更好，Mesh3在波谷表现更佳。

为进一步验证，引入波幅的相对误差wc。

式中，y为模拟波幅下的纵坐标；y理论为理论波幅下的纵坐标。

由图6可以看出：1）在波幅稳定后的35 s内，相对误差符合wc1＞wc2＞wc3；2）在35 s后，Mesh1的相对误差最大，Mesh2和Mesh3的相对误差相似；3）由图中的最大相对误差值小于0.6%可验证该模型造波情况良好。

图6 相对误差

根据以上结论和计算资源考虑，选取Mesh2的网格划分方式。

（2）网格时间独立性验证

对于瞬态计算，其结果不仅依赖于计算网格，还依赖于时间步长。在进行瞬态计算过程中，时间步长的选择固然是按照所需的时间分辨率来进行，但依照CFL条件，瞬态计算对于时间步长也存在明确的要求。

时间独立性验证方法与网格独立性方法类似，取不同大小的时间步长进行瞬态计算，考察相同时间点上的计算结果，分析其误差，最终确定最大时间步长用于计算[19]。

如图7中局部放大图所示，在波峰处，时间步长t=0.004 s下波幅较理论波幅下有较大的增大，t=0.01 s与t=0.02 s下波幅较吻合，且与理论波幅相比增大较小。在波谷处，3种时间步长下波幅与理论波长的差距：0.02＞0.01＞0.004，但t=0.01 s与t=0.004 s差距接近，都接近理论波长的波谷。故将对比情况与计算资源结合考虑，最终选择Mesh2网格，时间步长为0.01 s，计算时间为50 s。

图7 x=5 m处波幅历时曲线

2.5.4 模拟结果

为了验证模拟结果的正确性，需要将模拟结果进行多方面论证。

在计算前监测距进口5 m、10 m和15 m处波幅的历时曲线（图8），并与理论历时曲线进行对比。波幅随位置的变化并未发生波长方向的偏移，波幅方向与理论波幅吻合良好。故从历时曲线的结果可以验证模拟结果波幅的正确性。

图8 5 m、10 m和15 m处波幅历时曲线

消波区加入阻尼项的阻尼系数过小会导致消波不完全，影响波浪模拟的准确性，而通过加长消波区域的长度，则会增大计算机资源消耗。但阻尼系数过大，消波效果仿佛在20 m处放置了一堵墙，因此造成的反射同样会影响波浪模拟的正确性。如图9和图10所示，波浪在10 s左右到达20 m处的消波区，且在20～35 m之间消波效果良好,在35～40 m趋于平稳。整个区域的波浪情况在20 s时基本稳定，20～50 s消波区域保持稳定，且在整个模拟时间段内消波区末端水平面一直保持水平状态。故验证了波浪模拟的准确性和消波区良好的消波效果。

图9 5～25 s波幅历程曲线

图10 30～50 s波幅历程曲线

如图11所示，压力分布图随时间的变化过程。在波浪传递过程中，压力随着波浪水池深度的增加而增大，且压力数值与理论计算值相符。图中虚线为该时刻理论波幅的状态，通过与压力分布图对比得出，静压等值线与理论波幅相符。且静压强必然是沿着作用面的内法线法向，即可验证流动状态的波形与理论波形相符合。

图11 10～50 s压力分布图

一个周期内x方向水质点的运动速度u（图12）：波峰处最大；在波幅与静水面的交界处为零；在交界处运动到波谷处速度增至最大，方向为负；在波谷运动到下一个交界处速度再次降为0；最终回到波峰时速度再次回到正向最大，以此循环往复。

图12 x方向10～50 s 速度分布图

一个周期内y方向水质点的运动速度v（图13）：波峰处为零；在波幅与静水面的交界处达到正向最大；在交界处运动到波谷处降为0；在波谷运动到下一个交界处速度增至最大，方向为负；最终回到波峰时速度再次降为0，以此循环往复。

图13 y方向10～50 s速度分布图

以上模拟结果的分析皆符合理论上一个周期内水质点运动的速度情况。

经过以上论证，该模型的精确性良好。

3 推板造波

如图14所示，数值波浪水池与速度边界模型相同，但进口、出口边界条件及消波区的设置不同。

图14 推板造波水池模型

3.1 推板造波机理

推板造波方法模仿的是实验中的推板式造波机，通过动边界实现。控制方程如下：

式中，xp为推板位移；S为推板运动的冲程。

3.2 模拟设置

在模拟过程中发现推板造波方法的一些问题，对设置进行了一些调整。

（1）进口边界采用动网格layering方法。为预防计算中负体积的计算，需要根据分裂因子α及合并因子对初始高度h0进行设置，即可在边界运动中控制相邻边界的网格的分裂及合并[22]。

（2）阻尼项的适用范围从4L变为L，网格拉伸区域仍为4L。因为在模拟过程中发现适用边界速度法的设置会发生库朗数过大的错误。库朗数表达式如式（25）所示[19]：

在相同网格数量时减小时间步长，不改变阻尼项适用长度，在25 s左右发现计算不会收敛，经过调整松弛因子发现波幅会突然增大。最终将阻尼项适用长度改为L且网格拉伸区长度不变，得到的结果符合预期且收敛。

（3）时间步长的改变。原时间步长在推板造波方法上会发生库朗数过大，故降低时间步长，并对其进行时间独立性验证。

如图15所示，在波谷时y0.0025＞y＞y0.005＞y0.01（y为理论波幅），且y0.0025和y0.005与理论波幅相差较小，y0.01与理论波幅相差最大；在波峰时，y＞y0.0025＞y0.005＞y0.01，且y0.0025与y0.005几乎一致，y0.01与理论波幅的相差最大；y0.01波幅虽然也较好，但在30 s左右会发生不收敛的情况。故选择0.005 s为时间步长。