MPS方法模拟三维圆柱形液舱晃荡问题

田鑫 ，万德成*

1上海交通大学船舶海洋与建筑工程学院，上海200240

2上海交通大学海洋工程国家重点实验室，上海200240

3高新船舶与深海开发装备协同创新中心，上海200240

0 引 言

液舱晃荡是指在外部激励作用下，舱内装载的部分液体所产生的波动及其与舱壁结构相互作用的现象。随着世界航运业的发展以及对能源需求的不断提高，相继开发应用了大型原油货轮、液化石油气船以及液化天然气船等载液货船。这些容积巨大的船装载着巨量的液体货物，其内部产生晃荡现象往往会对船舶产生巨大危害，因此晃荡现象成为船舶与海洋工程研究的热点问题。

船舶液舱的形状多种多样，以往对液舱晃荡的研究主要集中在矩形液舱及薄膜型液舱上，例如，张书谊和段文洋［1］就曾采用CFD软件Fluent对二维矩形液舱不同舱内水深、不同激振频率时的横荡进行过数值计算和对比分析。对于其他形状液舱的晃荡，国内外的研究则较少。Kobayashi等［2］通过实验方法，对水平圆柱形液舱的晃荡进行了研究，并对载液率、幅值、液舱长度等因素进行了分析探讨。Faltinsen等［3］对球形液舱中的液舱晃荡现象予以了研究，并对其中的线性晃荡进行了理论分析。Gavrilyuk等［4］对直立环形圆柱的晃荡进行了理论及实验研究，分析了环形挡板的减晃作用。刘戈等［5-6］通过实验，研究了独立C型液舱在纵摇激励下的晃荡，讨论并分析了其晃荡时自由液面的波形特征，并对该型液舱的频域共振特性进行了研究，发现0.5倍理论固有频率为晃荡起始频率，终止频率会随着液位的升高而下降。

本文将基于课题组自主研发的无网格粒子法求解器MPSGPU-SJTU，对横荡激励下圆柱形液舱中的晃荡现象进行研究。首先，模拟50%充水率下，激励频率为一阶固有频率时的晃荡现象，通过对比实验结果，验证模拟的准确性。然后，对不同激励频率下的晃荡现象进行计算，对比不同激励频率下的液舱受力与流场情况，以分析激励频率对圆柱形液舱晃荡的影响。

1 数值方法

MPSGPU-SJTU求解器是在本课题组之前开发的CPU并行求解器MLParticle-SJTU［7-9］的基础上，引入 GPU（Graphics Process Unit）加速技术开发的新一代求解器，其数值方法为经张雨新［10］改进的移动粒子半隐式（Moving Particle Semi-implicit，MPS）方法，GPU加速基于CUDA平台实现。

1.1 控制方程

在MPS方法中，控制方程包括连续性方程和N-S方程，对于不可压缩流体，其形式如下：

式中：ρ为流体密度；P为压力；V为速度向量；Fm为质量力，一般为重力；ν为流体的运动粘性系数；t为时间。

1.2 核函数

在粒子法中，控制方程将被写成粒子形式，两个粒子间的相互影响是通过核函数来实现的，其与传统的核函数不同。MPSGPU-SJTU求解器采用的核函数为

式中：r为两个粒子间的距离；rε为粒子作用域的半径。

1.3 梯度模型

MPSGPU-SJTU求解器采用的动量守恒的梯度模型为

式中：D为空间维数；ri，rj分别为粒子i及其邻居粒子j的坐标；Pi，Pj分别为粒子i和j的压力；n0为初始粒子数密度。

1.4 拉普拉斯算子模型

在MPS方法中，拉普拉斯算子模型如下式所示：

式中，φi和φj为物理量φ在i，j处的值。

其中，

系数λ的引入是为了使数值结果与扩散方程的解析解相一致。

1.5 压力泊松方程

MPSGPU-SJTU求解器采用的压力泊松方程是由Tanaka 和 Masunaga［11］提出的混合源项法（Mixed source term method）。该方法结合了速度散度和粒子数密度，后来被Lee等［12］写成了更为合理的表达形式：

1.6 自由面判断方法

在MPS方法中，对自由液面的准确判断对计算的精度和稳定性来说十分重要。在本文采用的自由面判断方法中，首先定义矢量

式中，F为粒子数密度的不确定度。计算F的模，当粒子满足时，即被判定为自由面粒子，其中α=0.9，为一参数。

1.7 GPU加速技术

GPU具有强大的计算能力，在并行计算中，应用GPU能获得显著的效果。NVIDIA公司推出的CUDA技术为GPU编程提供了方便的平台，MPSGPU-SJTU求解器即是基于此平台而开发。本研究对MPS方法的计算流程进行了合理的设计，实现了MPS方法主要计算过程的GPU加速。计算流程如图1所示。

图1 MPSGPU-SJTU求解器计算流程图Fig.1 The flow chart of MPSGPU-SJTU solver

2 模型验证

本文构造的数值模型是参考Kobayashi等［2］的模型而建立，其物理模型如图2所示。液舱为圆柱体，圆柱底面直径0.47 m，长0.94 m，载液率为50%。

图2 几何模型Fig.2 Geometric model

液舱横荡激励的运动方程为

式中：Ax为激励幅值，Ax=0.015 m；f为激励频率，Hz。

构造的数值模型如图3所示。流体粒子数为79 258，总粒子数为137 380，粒子间距0.01 m。水的密度ρ=1 000 kg/m3，运动 粘性系 数ν=10-6m2/s，重力加速度g=9.81 m/s2，时间步长dt=0.000 5 s。

图3 数值模型Fig.3 Numerical model

本节对f=1.2 Hz、横荡激励幅值Ax=0.015 m的晃荡进行了模拟。图4给出了横荡激励作用下，实验与模拟中液舱所受水平半径方向合力Fy的时历曲线对比。从图中可以看出，数值模拟给出的合力曲线与实验结果吻合较好，两者的变化趋势、周期峰值基本一致。可见，采用求解器MPSGPU-SJTU能够很好地预测液舱所受合力。

图4 液舱所受水平半径方向合力的时历曲线Fig.4 Time history curves of horizontal radial direction resultant force applied to tank

3 激励频率对晃荡的影响

本节分别对激励频率f=1.0，1.1，1.2，1.3，1.4 Hz，幅值Ax=0.015 m的横荡激励进行了数值模拟。图5给出了不同激励频率下晃荡流体的运动情况。

图5 不同激励频率下的流场Fig.5 Flow field at different frequencies

由图中可以看出，激励频率对晃荡流体的流动特征具有显著影响。流场随频率变化的过程可分为以下4个阶段：

1）如图5（a）所示，当激励频率f=1.0 Hz时，晃荡的幅度很小，自由面非常光滑，流体运动比较平稳。

2）随着激励频率的增加，晃荡波的波陡迅速增大，反应十分剧烈。当f=1.1和1.2 Hz时（图5（b）和图5（c））出现了冲顶现象，晃荡波沿壁面爬升至舱顶后下落，形成一股射流猛烈拍击在自由面上，从而发生了波面破碎和液体飞溅现象。整个晃荡过程非常剧烈，并伴随着一些复杂的流动现象，具有较强的非线性。其中当f=1.2 Hz时晃荡波的运动更加剧烈，其沿舱壁爬升的高度更高，并且流体动能更大，一部分液体甚至飞溅到了另一侧舱壁。

3）随着激励频率的继续增大，晃荡波的波高和波陡迅速减小。当f=1.3 Hz时（图5（d）），晃荡虽然仍较剧烈，但晃荡波的动能已不足以支撑其爬升至舱顶。波浪在砰击壁面并爬升很小的高度后，波面翻卷冲击自由面，发生了轻微的破波现象。在这一过程中，液体的运动幅度是随着激励频率的增加而减小的。

4）当f=1.4 Hz时，液舱运动较快，但晃荡波的波幅继续减小，流动变得平缓，不再出现波浪的翻卷破碎现象。

图6所示为不同激励频率下的液舱受力时历曲线。从整体上看，其呈现出与激励频率一致的周期性。但从细节上可以看到，液舱受力曲线的形状和峰值也随着不同激励频率而发生变化。当激励频率f=1.0 Hz时，受力曲线为单峰特征（图6（a）），受力的峰值相对较小，其大小约为100 N，此时晃荡现象较为平稳，没有发生拍击现象。当激励频率增大至1.1或1.2 Hz时，晃荡产生共振现象，受力曲线出现了双峰现象。其中第1个波峰达到了250 N甚至是300 N以上，这是流体直接冲击壁面形成；而第2个受力峰值大小约为100～150 N，是由爬升的流体回落，对底部流体产生拍击作用所致。随着激励频率继续增大，受力峰值开始下降，当f=1.3 Hz时，受力曲线仍然存在双峰特征，但2个峰值已经非常接近，此时，流场仍然存在明显的非线性流动。当激励频率增大至1.4 Hz时，只在一些周期内出现了第2个受力峰值，这是由舱壁附近微弱的破波现象所致。

图7给出了受力峰值随激励频率的变化规律。从图中可以看出，随着激励频率的增加，受力峰值是先增大后减小。

当f＜1.1 Hz时，受力峰值与激励频率呈正相关。当f=1.0 Hz时，流体的运动幅度较小，晃荡较为平缓，液舱的受力幅值在100 N左右。随着激励频率的增大，受力峰值增大，在1.1～1.2 Hz时达到峰值，约为220 N，此时晃荡最剧烈，说明此时处于共振频率。当f＞1.2 Hz时，受力峰值与激励频率呈负相关。随着激励频率的增大，液舱所受水平半径方向合力的峰值减小，晃荡幅度减小。当f=1.3 Hz时，受力的幅值降至120 N左右。当f=1.4 Hz时，受力幅值降至90 N。

图6 不同激励频率下的液舱受力Fig.6 Forces acting on tanks at different excitation frequencies

图7 液舱受力峰值随激励频率的变化Fig.7 Variation of peak force in tank with excitation frequency

4 结 论

本文采用课题组自主开发的MPSGPU-SJTU求解器，研究了圆柱形液舱在横荡运动中激励频率对液舱受力的影响，主要得到以下结论：

1）MPSGPU-SJTU求解器能够较好地模拟圆柱形液舱内的液体晃荡现象，能很好地捕捉晃荡现象中的液体飞溅和波浪翻卷破碎等现象，并且也可以较为准确地计算液舱受力。

2）不同激励频率下，晃荡的波形区别很大。当f=1.1与1.2 Hz时，有冲顶现象出现；当f=1.3 Hz时，有轻微的破波现象出现；当f=1.1和1.4 Hz时，流动较平缓。

3）当激励频率在固有频率附近时，晃荡最为激烈，液舱的受力幅值最大，但随着激励频率远离固有频率，晃荡逐渐减弱。

4）在船舶运营中，要避免在液舱固有频率附近产生激励振动。

本文研究了激励频率对圆柱形液舱晃荡的影响，给出了液舱整体受力的结果，可为圆柱形液舱的设计应用提供有效依据。在实际液舱晃荡问题中，由于很多破坏是由于局部应力过大所造成，所以在今后的研究中，会计算液舱的局部压力，对液舱晃荡中的现象与载荷进行更为细致的研究。