动态边界粒子SPH法在养殖工船横摇液舱晃荡模拟中的应用比较

赵新颖，黄温赟，黄文超，管延敏

(1 中国水产科学研究院渔业机械仪器研究所，农业农村部远洋渔船与装备重点实验室，上海 200092;2 青岛海洋科学与技术国家试点实验室深蓝渔业工程联合实验室，山东 青岛 266237;3 江苏科技大学 船舶与海洋工程学院，江苏 镇江 212003)

近年来，作为渔业转型升级和绿色发展的重要方向，工船养殖受到广泛关注，养殖水舱是养殖工船的核心部分，与载液船舶相似，养殖水舱在波浪作用下产生的液舱晃荡现象，对船体性能以及液舱舱壁结构强度产生影响，进而影响养殖安全。为了对此类现象进行研究，有限差分法[1]、有限体积法[2]、有限元法[3]、边界元法[4]、无网格法(MPS[5]、SPH[6]、CPM[7])等诸多方法相继被用于晃荡现象数值预报。剧烈的液体晃荡会产生的自由液面翻卷、破碎等现象的模拟给传统有网格方法带来一定的难度，光滑粒子流体动力学(SPH)法作为无网格法的一种，摆脱了传统方法对网格、单元的依赖，避免了拉格朗日方法中的网格扭曲及网格重构问题。该方法通过粒子的运动模拟流场流动现象，可以实现对自由液面的精确追踪[1]，在处理大变形[9]、动边界[10]、自由液面翻卷、破碎[11-12]等问题时具有明显的优势，因此SPH已成为研究液舱晃荡问题的热门方法。如李大鸣等[13]对中、高液位液体晃荡及防荡隔板对液体晃荡现象的抑制作用进行了数值模拟，Shao等[14]提出了一种改进SPH方法用于模拟液体晃荡现象，Gotoh等[15]运用不可压缩SPH法模拟了剧烈液体晃荡现象，曹雪雁等[16-17]运用SPH法模拟了矩形箱体液体晃荡现象以及船舶破舱进水特性。

SPH方法中，壁面边界的处理是提高SPH的精度和稳定性的重要因素，一直是SPH研究关注的热点[18]。位于固壁边界上或边界附近的粒子由于其核函数被边界截断，使得该区域的变量梯度计算不准确，产生截断误差从而引起数值结果畸变，无法保证二阶精度[19]。近年来，不同的壁面边界的处理方法不断被提出，传统的镜像对称边界法[20]受镜像对称算法的限制，一般只应用于规则的平面或者直角边界，在复杂边界问题中很少被使用；Monaghan[21]提出的排斥力法方法简单，但因守恒性差，无法精确模拟粒子在边界附近的运动，影响了壁面载荷计算精度；Libersky等[22]引入的镜像粒子法能够对内部粒子进行准确计算，但可能会产生局部流体粒子穿透边界现象；Marrone等[12]在镜像粒子的基础上提出的固定虚粒子法，满足了壁面不可穿透条件，在处理复杂边界问题时，依然存在虚粒子布置的困难；动态边界粒子由Crespo等[24]详细介绍，并在Liu等[25]的研究中被进一步使用，该方法虚粒子满足与流体粒子相同的连续性和状态方程，可以与流体粒子在同一循环内简单地进行计算，可以节省大量的计算时间。

为准确模拟养殖工船养殖舱自由液面变化以及计算舱壁受力情况，保证养殖工船结构强度和养殖安全，本研究采用动态边界粒子SPH法对养殖舱横摇激励下液体晃荡现象展开模拟，并对动态边界粒子处理的两种方法展开比较研究，数值结果表明本文采取的方法精确有效，可以为养殖水舱晃荡现象模拟提供数值参考。

1 SPH法数值模型

1.1 控制方程及其离散形式

对于拉格朗日形式下的弱可压缩流体，其连续性方程和动量方程可以表示为：

(1)

(2)

式中：p、ρ、u、μ分别表示流体的压力、密度、速度和黏性；g为重力加速度。流体的弱可压缩性采用Monaghan等[21]提出的人工压缩法，如式(3)。

(3)

利用SPH法的基本原理，连续性方程和动量方程可以离散为：

(4)

(5)

式中：pi、ρi、ui分别是第i个粒子的压力、密度和速度；mj、ρj为粒子i支持域中相邻粒子j的质量和密度。ij为人工黏性项[26]，其表达式为：

(6)

1.2 核函数

核函数对数值模拟的精度、效率及稳定性，目前应用较为广泛的核函数有B-样条核函数、3次样条核函数、5次样条核函数和高斯核函数等[27]。Dehnen等[28]研究发现，随着光滑长度和周围粒子数量的增加，传统的核函数受到聚集不稳定因素粒子容易产生结对现象，推荐了Wendland核函数[29]。本研究在计算过程中采用Wendland的C2核函数：

(7)

式中；对于二维问题αd=7/4πh2。

1.3 时间积分

本研究采用预估校正法求解控制方程，计算粒子位置、密度和速度。

在预估步有：

(8)

在校正步有：

(9)

(10)

2 壁面边界处理的动态边界条件

2.1 动态边界条件基本原理

动态边界条件是在边界处布置一组动态边界粒子(图1)，这些边界粒子跟流体粒子一样参与连续性方程计算，不参与动量方程求解，当流体粒子靠近边界与边界粒子的距离减小至2倍光滑长度以下时，受流体粒子影响，边界粒子密度增加导致压力增大，边界粒子压力的变化又反作用于流体粒子动量方程，引起流体粒子速度的变化，从而防止流体粒子穿透边界。动态边界粒子位置可以保持固定不动，也可以通过某种方式使其运动，本研究针对这两种方式展开比较分析。

图1 动态边界条件示意图Fig.1 Schematic diagram of dynamic particle condition

2.2 动态边界粒子运动

假定壁面平移速度为u0，绕轴心O旋转的角速度为ω，对于任一边界粒子，其运动速度uw为

uw=u0+ω×(rw-rO)

(11)

式中：rw、rO分别为边界粒子、旋转轴心的空间位置。

在预估步，边界粒子的位置为：

(12)

在校正步，边界粒子的位置为：

(13)

在每一预估步及校正步，先进行边界粒子速度、位置的更新，再求解流体粒子连续性方程、动量方程。

2.3 动态边界粒子固定不动

以物体壁面随体坐标系为参考坐标系，保持壁面边界不动，流场所受重力方向发生变化，对于二维问题，有如下公式：

(14)

式中：gx、gy为重力g的两个分量；θ为物体绕轴心O旋转的角位移。

任一流场粒子a所受质量力为：

(15)

式中：第一项为角加速度项，第二项为向心加速度项，第三项为科氏加速度项。

将公式(15)代入公式(5)，动量方程变为：

(16)

3 数值算例

3.1 轻微液体晃荡算例分析

通常来讲，水平激励的幅值小于15%的箱体内液面半径时为小振幅激励，即轻微晃荡；反之，则为大振幅激励，相应地为剧烈晃荡。针对轻微液体晃荡，分别采用两种动态边界粒子模型对其进行数值模拟。计算模型定义如图2所示，矩形箱体的长度L=0.92 m，高度H=0.62 m，水深h=0.465 m。不考虑空气影响，水箱里的水密度为ρ=1 000 kg/m3。对箱体施加横摇激励使其绕箱体中心(轴x=L/2,y=H/2)转动，横摇角位移θr=θ0sinωrt，式中，θ0为最大角位移，ωr为横摇圆频率。

图2 矩形箱体坐标系设置示意图Fig.2 Schematic diagram of sloshing tank and the co-ordinate system

在水箱右壁面距离底部0.17 m处设置一个压力监测点A(x=0.92 m,y=0.17 m)，并通过与文献[31]的实验值进行比较对该算例进行验证。图3为在θ0=4°、ωr=2.0 rad/s工况下，分别采用两种动态边界粒子运动方法A点压力随时间的变化曲线以及与实验结果的对比。

图3 点A处压强比较Fig.3 Pressure comparison at point A

由图3可见，两种方法得出的计算结果趋势基本一致，但相对于边界粒子运动法，边界粒子固定不动法计算值与文献[31]更为接近。

为了进一步比较两种方法对计算精度的影响，对作用在箱体上的水动力和旋转绕轴的力矩进行了对比，如图4所示。由图4可见，边界粒子运动法在初始计算时会产生较大数值振荡，在t=1.0 s左右计算趋于稳定，两种计算方法得到的水动力和力矩的变化趋势基本一致，两种方法计算值大小略有差异。

图4 作用于箱体上的力和力矩Fig.4 Time history of hydrodynamic force and moment

图5为两种方法在模拟过程中几个典型时刻压力场的对比，左边为边界粒子运动的压力场，右边为边界粒子不动的压力场。总体来看两种方法的模拟结果接近，均能很好地捕捉自由液面形状。

图5 计算结果比较Fig.5 Comparison of the numerical results between moving boundary particles(left)and fixed particles(right)

3.2 剧烈液体晃荡算例分析

定义如图6所示的矩形箱体，箱体宽L=0.90 m，高H=0.508 m，h=0.093 m对箱体施加横摇激励使其绕箱体中心(轴x=0,y=0)做简谐振动，根据箱体自振频率表达式：

图6 矩形箱体坐标系设置示意图Fig.6 Schematic diagram of sloshing tank and the co-ordinate system

(17)

得ω0=1.919 rad/s，选取ωr=1.631 rad/s(ωr=0.85ω0)，横摇角位移随时间的变化曲线如图7所示。建立绝对坐标系，位于初始时刻箱体的下壁面中点，坐标系不随箱体振动而运动。二维箱体初始水深h=0.093 m，内部区域采用矩形布点，粒子间距Δx=Δy=0.002 m，CFL常数设置为0.2。

图7 横摇角位移曲线Fig.7 Time history of roll angle

为了验证剧烈晃荡工况下模型的精度，在距离箱体底部0.093 m处的左侧壁面设置一个压力监测点，坐标值为B(x=-0.45 m,y=0.009 3 m)，用来检测该点水压随时间的变化。图8为采用两种边界粒子法计算出的B点处压强随时间变化的曲线以及与文献[32]试验值的比较。由图中可见两种方法计算出的压力峰值的时刻与试验值是一致的，三条压力曲线的趋势变化也基本一致，算例表明对于水压值的计算，两种方法均具有良好的计算精度。

图8 点B处压强比较Fig.8 Pressure comparison at point B

图9为两种边界粒子法自由液面高度变化模拟结果与实验值的比较，可见两种方法均能捕捉到舱壁抨击、液体飞溅(t=3.45 s、t=4.35 s)以及自由液面翻卷(t=3.95 s、t=4.75 s)现象，适用于剧烈液体晃荡现象数值模拟。

图9 边界粒子运动(图左)、边界粒子固定不动(图中)计算结果与试验值(图右)比较Fig.9 Comparison of the numerical results of moving boundary particles(left)and fixed particles(middle) with experiments(right)

图11为相同时间步长不同粒子间距下计算结果比较，可见随粒子间距的变小两种方法计算值与试验值均误差变小，边界粒子固定法在粒子间距变化时，计算结果变化较小，相对于边界粒子运动法，边界粒子固定法具有更好的收敛性。

图10 边界粒子运动(图左)、边界粒子固定不动(图右)不同时间步长下计算结果比较Fig.10 Comparison of the numerical results of moving boundary particles(left)and fixed particles(right) by different time steps

图11 边界粒子运动(图左)、边界粒子固定不动(图右)不同粒子间距下计算结果比较Fig.11 Comparison of the numerical results of moving boundary particles(left)and fixed particles(right)by different particle spacings

4 结论

以养殖工船养殖水舱晃荡现象为研究对象，运用SPH法分别对横摇激励下轻微液体晃荡现象和剧烈液体晃荡现象展开模拟研究，并引入动态边界粒子法解决SPH法固壁边界上或边界附近的粒子截断误差影响精度问题，实现了对SPH法的改进。研究发现，边界粒子运动法、边界粒子固定不动法均能有效地进行液舱晃荡现象数值模拟；相对于边界粒子运动法，边界粒子固定不动法具有更好的算法稳定性。本研究方法适用于液体晃荡现象数值模拟，可以为养殖工船养殖水舱晃荡模拟提供数值方法参考。

□