SPH求解器的开发及其在液舱晃荡数值模拟中的应用

尹梓炜 陈晓明 赵成璧

摘 要：本文提出一种改良的SPH固壁边界条件处理方法，并基于此开发适用于模拟液舱晃荡的数值求解器。采取该求解器对液舱晃荡现象进行数值模拟，结果与相关文献结果相当吻合，验证了其可靠性和准确性。在此基础上，通过计算分析横荡状态下不同隔板高度液舱晃荡时所产生的壁面压力和液体重心高度变化趋势，探讨了隔板高度对液舱晃荡的影响规律和物理机制，并为防荡隔板的设计提出合理建议。

关键词：SPH；固壁边界条件；液舱晃荡；数值模拟

中图分类号：U661.74 文献标识码：A

Abstract： A modified treatment for rigid boundary condition of SPH method is proposed in this paper， and a numerical solver that suitable for sloshing simulation is developed based on this SPH method. The phenomenon of liquid sloshing in a tank is simulated by using this solver， and the numerical simulation results are in good agreement with the results from the related literature， which verifies the reliability and accuracy of the model. On the basis of that， through the calculation and analysis of the wall pressure and the trend of liquid gravity center height， which generated by sloshing with different baffle lengths under the sway condition， this paper preliminarily discusses the effect and physical mechanism of baffle length on sloshing and proposes reasonable suggestions about the design of anti-sloshing baffle.

Key words： SPH； Rigid boundary condition； Sloshing； Numerical simulation

1 引言

液艙晃荡是指在部分充液的容器中，其内部液体在外部激励下产生运动的现象。液舱晃荡是船舶航行中遇到的常见现象，具有较强的随机性和非线性。液体晃荡产生的瞬时压力可对舱壁结构造成损伤，同时晃荡产生的附加力矩会影响船舶稳性，危害船舶安全，因此液舱防晃是船舶工程领域的基础性课题，具有重要的工程应用价值。

液舱晃荡问题的研究方法主要分为实验、理论分析和数值仿真。随着计算机科学的发展，数值分析已经成为研究液舱晃荡问题的重要手段。传统的数值方法如MAC、VOF等是基于网格对流场进行求解，在处理液舱晃荡这类带自由液面的强非线性运动问题时存在困难。光滑粒子流体动力学方法（SPH），最早由Lucy[1]和Monaghan[2]等分别独立提出用于解决天体物理学问题，后经过改进引入到水动力学中。SPH作为一种无网格的拉格朗日方法，它克服了无网格方法在处理流体自由面大变形问题上的弱点，近年来逐渐被应用于液舱晃荡问题的研究。

较早期的有Iglesias[3][4]等利用SPH方法先后研究了减摇水舱晃荡特性和不同频率下液舱晃荡力矩幅值等问题，后来崔岩[5]等运用SPH方法模拟分析了二维矩形水槽在激励频率接近一阶固有频率的纵荡过程中的晃荡现象。在基于SPH方法的液舱晃荡数值模拟中，边界条件的处理十分关键，而上述文献均采用Lennard-Jones斥力模型实施固壁边界条件。该模型对靠近边界的粒子施加强斥力，这容易导致粒子的初始分布被破坏，且时间步长受到严格的限制。同时为了防止粒子穿透边界，特别是在边界形状突变处，往往使用多层边界粒子，降低了计算效率，且会导致边界附近产生明显的压力振荡。为改善该问题，Shao[6]等提出耦合动力学边界的SPH方法并研究了二维液舱晃荡问题，得到较为精确的流场压力。Chen[7]等在此基础上修正其人工斥力的表达式，模拟了二维矩形容器在液舱晃荡下的舱壁砰击压力并进行了实验对比，获得了让人满意的结果。

基于上述研究成果，为了获得更为精确的液舱晃荡模拟结果，本文在提出一种改良的固壁边界条件处理方法基础上，自主开发了适用于液舱晃荡的SPH数值求解器，并进一步研究了隔板高度对晃荡行为的影响，从而扩展了SPH方法在液舱晃荡中的应用。

2 SPH求解器的开发

2.1 SPH基本方法和理论

SPH方法中[8]，计算域被离散为一系列具有相互作用的粒子，某时刻粒子的物理量可由一定光滑半径内的周围粒子通过核函数（如三次样条函数等）估计得到。基于SPH的核近似和粒子近似思想，考虑重力下控制流体运动的Navier-Stokes方程可离散为如下形式：

压力通过下列状态方程求得：

2.2 自适应边界粒子法向斥力模型

Akinic[10]提出了一种流固耦合边界方法，该方法引入边界粒子质量函数，考虑边界粒子的影响对流体粒子的密度和受力进行修正。本文基于该方法的思想并对其进行改进，只考虑法向斥力，提出一种改良的自适应边界粒子法向斥力模型，得到考虑边界粒子b的贡献的SPH控制方程如下：

由上式可知，由于核函数随距离增加而单调递减的性质，边界粒子分布密集的地方不会产生过大斥力，分布稀疏的地方也不会产生过小的斥力，因此相对于传统SPH边界处理方法具有自适应性。且该方法斥力大小又与流体粒子自身压力大小有关，只在法向施加排斥力，减少边界粒子对流场的扰动，可以使粒子在边界处受力更均匀。

以简单的溃坝模型为算例。从图1的速度矢量对比可以看出：传统的Lennard-Jones斥力模型中溃壩水头前沿的粒子运动十分不均匀，明显是由于产生的瞬时斥力过大导致；而本文方法中的底部边界粒子产生足够而均匀的斥力，有效避免粒子在固壁边界发生非物理穿透的现象，同时使粒子运动平滑有序，证明本文方法的效果明显优于传统的斥力模型。

2.3 算法流程

SPH求解器的算法流程如图2。

其中：邻域粒子搜索方法有全配对搜索法、链表搜索法和树形搜索法等方法；时间积分的方法可采用蛙跳法、预测校正法和龙格库塔法等方法。本文采用链表搜索法进行邻域搜索，采用蛙跳法进行时间积分，时间步长满足CFL条件，具体见参考文献[8]。

3 基于SPH求解器的液舱晃荡数值模拟

3.1 典型液舱晃荡实验

Chen[7]的实验是研究液舱晃荡的典型实验，并给出了较准确的SPH模拟结果。为了验证本文液舱晃荡SPH方法数值求解器的可靠性和精确性，本文模拟了Chen液舱晃荡算例，并与其实验结果和模拟结果进行比较。

计算模型中矩形液舱长和高为L=H=1 m、充液深度D=0.3 m、横摇频率3.81rad/s、振幅5°，容器底部中心处为轴心。计算采用流体粒子数12 870、时间步长1×10-4 s进行模拟，流体压缩率控制在0.1%附近。

图3和图4分别为矩形液舱在激励频率3.81 rad/s的横摇下晃荡的四个典型时刻的chen实验结果和本文模拟结果。从图中可以看出，本文的SPH模型有效地模拟了自由液面形态，同时获得了较精确的均匀压力分布。

在本文中，基于SPH的核近似和粒子近似的思想，边界粒子b的压强通过周围流体粒子根据

插值估算。取距底部0.2 m高度处舱壁作为压力监测点。

图5所示为压力监测点在计算时段2～10 s的压力时间历程。由图5可知，本文SPH方法的压力计算结果与Chen的数值模拟结果的变化趋势基本一致，同时压力的峰值则与Chen的实验结果更为接近，证明本文SPH方法在液舱晃荡的数值模拟应用具有更高的可靠性和精确性。

3.2 防荡隔板高度对液舱晃荡影响

为进一步研究不同液舱结构对晃荡行为的影响，本文在与3.1算例尺寸相同的矩形液舱中添加防荡隔板，模拟了不同高度的I型单隔板在横荡频率接近液舱一阶模态的横荡工况下的防荡效果对比。

根据线性理论，由 可得矩形液舱晃荡的一阶固有频率为ω1=4.76 rad/s，对液舱施加简谐激励x=Asin（ωt），取ω=1.166ω1、振幅A=0.05 m。在横荡工况下，为达到防荡效果最大化，在液舱中央竖向布置I型隔板，隔板高度l分别取0.25 D、0.5 D和0.75 D。同时为了分析方便，取距离舱壁底部0.3 m处侧壁作为压力监测点。

图6（1）描述了不同单隔板尺寸工况下舱壁压力计算监测点在计算时段2～10 s的压力时间历程。图6（2）描述了不同单隔板尺寸工况下在计算时段0～10 s液体重心高度的时间历程。由图6（1）可知，随着隔板高度的增加，整体上压力峰值呈减少趋势且双峰特性减弱，当隔板高度为0.75 D时最大压力峰值能降至无隔板状态下的24.5%。同时，当隔板增加到一定高度后，压力峰值对隔板高度变化的敏感度降低。由图6（2）可知，隔板高度和液体重心高度呈明显的正相关，其中当隔板高度为0.75D时，重心高度波动最大幅值可降至无隔板状态时的19.2%。

为了进一步研究液舱晃荡时隔板高度对压力峰值影响的物理机理，我们对不同隔板高度下的液舱流场速度进行对比分析。

图7反映了t=2.26 s时刻不同隔板高度的流场速度矢量和舱壁压力。可以看出：随着隔板高度的增加，隔板附近流场受到干扰，出现了紊乱并逐渐产生了涡。此时，涡街阻尼效应对较大波高的形成产生了削弱作用，有效抑制了液舱晃荡时的液面剧烈波动，从而有效地降低舱壁晃荡时的舱壁压力。

不过，结合图6（1）可知，当隔板到达一定高度后，其对舱壁压力峰值的影响效果将减弱；同时图7表明，随着隔板高度的增加，隔板的涡街振动响应不断增强，隔板的液体冲击荷载作用不断增大，这会导致隔板端部应力集中，因此从安全性和经济型考虑，隔板高度应适中，建议隔板在具体设计时可采用本文数值模拟方法确定隔板高度尺寸。

4 结论

本文利用自主开发改进的SPH求解器研究了二维矩形液舱晃荡问题，取得较好的应用效果。由数值模拟结果可以看出：

（1）本文改良的固壁边界条件处理方法，能有效阻止粒子穿透边界并使粒子在边界处的运动均匀有序。以此为基础自主开发的SPH求解器具有较高的可靠性和精确性，能很好地捕捉液舱晃荡时的自由液面状态，获得较为准确的流场压力值，是对液舱晃荡现象有效的数值模拟手段，具有较好的应用前景；

（2）本文进一步的研究结果表明，防荡隔板能明显降低液舱晃荡的幅度和舱壁压力峰值，但从安全性和经济性的角度出发，隔板高度应适中，不宜过高。

此外，SPH方法中在增加算法稳定性、提高计算精度和效率等方面尚有许多需要完善之处，本文的研究为进一步开发高效通用的SPH求解器以及拓展其应用领域打下了基础。

参考文献

[1]Lucy L B. A numerical approach to the testing of the fission hypothesis[J]. Astronomical Journal， 1977， （82）.

[2]Gingold R A， Monaghan J J. Smoothed particle hydrodynamics： theory and application to non-spherical stars[J]. Monthly Notices of the Royal Astronomical Society， 1977， 181（3）.

[3]Souto Iglesias A， Pérez Rojas L， Zamora Rodríguez R. Simulation of anti-roll tanks and sloshing type problems with smoothed particle hydrodynamics[J]. Ocean Engineering， 2004， 31（8）.

[4]Souto-Iglesias A， Delorme L， Pérez-Rojas L， et al. Liquid moment amplitude assessment in sloshing type problems with smooth particle hydrodynamics[J]. Ocean Engineering， 2006， 33（11）.

[5]崔巖， 吴卫， 刘桦. SPH方法模拟二维矩形水槽晃荡过程[C]// 全国水动力学研讨会. 2006.

[6]Shao J R， Li H Q， Liu G R， et al. An improved SPH method for modeling liquid sloshing dynamics[J]. Computers & Structures， 2012， s 100–101（6）.

[7]Chen Z， Zong Z， Li H T， et al. An investigation into the pressure on solid walls in 2D sloshing using SPH method[J]. Ocean Engineering， 2013， 59（2）.

[8]G.R.Liu， M.B.Liu， Liu，等. 光滑粒子流体动力学：一种无网格粒子法[M]. 湖南大学出版社， 2005.

[9]Monaghan J J， Gingold R A. Shock simulation by the particle method SPH[J]. Journal of Computational Physics， 1983， 52（2）.

[10]Akinci N， Ihmsen M， Akinci G， et al. Versatile rigid-fluid coupling for incompressible SPH[J]. Acm Transactions on Graphics， 2012， 31（4）.