刘文夫，薛鸿祥，唐文勇

（1.上海交通大学 海洋工程国家重点实验室，上海 200240；2.高新船舶与深海开发装备协同创新中心，上海 200240）

计及弹性支撑效应的独立液舱晃荡数值分析研究

刘文夫，薛鸿祥，唐文勇

（1.上海交通大学 海洋工程国家重点实验室，上海 200240；2.高新船舶与深海开发装备协同创新中心，上海 200240）

文章针对弹性支撑结构作用下的液舱晃荡问题，采用基于Hilber-Hughes-Taylor（HHT）格式的隐式直接积分法求解液舱运动，应用VOF法求解流体晃荡，并结合部分单元参数的概念处理棱形液舱边界，建立了液舱运动与流体晃荡双向耦合迭代算法。以独立液舱为研究对象，通过耦合求解弹性支撑液舱晃荡，并与非弹性支撑液舱晃荡进行比较，分析了在不同刚度的弹性支撑结构作用下液舱液体运动及晃荡载荷的变化规律。该文建立的分析方法为含弹性支撑的独立液舱晃荡载荷预报提供了一种快速、有效的分析手段。

弹性支撑；独立液舱；液舱晃荡；VOF方法

0 引 言

随着世界经济对LNG需求量的增长，各种型式的LNG船和海上浮式LNG装置概念得到广泛关注。目前世界范围内采用的LNG液舱主要分为薄膜型和自支撑式两种[1]，其中自支撑式液舱完全由自身支撑，不构成船体的一部分，是完全独立的液舱（以下简称独立液舱），包括日本石川岛播磨重工株式会社（IHI）研制的SPB棱形液舱（Self-supporting,Prismatic-shape IMO type B）和Moss Maritime研制的Moss球形液舱。独立液舱具有无装载限制、可靠性高、易于经营和维护等优点。

晃荡问题一直以来是大型LNG液舱设计的关键技术，Faltinsen等[2-3]曾基于线性势流理论在晃荡幅值较小的情况下对液舱晃荡进行了研究，然而在剧烈晃荡的情况下，液面伴随有波面破碎、飞溅和合并等强非线性运动现象[4]，理论分析方法已很难适用，目前主要通过模型试验和数值模拟进行研究。解决液体晃荡的关键问题是如何描述自由液面的变化，因此一系列捕捉或追踪自由液面运动的方法和技术得到了发展。其中有基于网格的MAC（Marker and Cell）法、VOF（Volume of Fluid）法、Ls（Levelset）法、CIP（Constrained Interpolation Profile）法，基于无网格的SPH（Smoothed Particle Hydrodynamics）法、MPS（Moving Particle Semifinal Implicit）法、MLPG（Meshless Local Petrov-Galerkin）法等[5]。朱仁庆[6]采用VOF法对矩形液舱内液体晃荡进行了数值模拟。沈猛等[7]用改进的VOF法和混合自由表面边界速度条件对棱形液舱液体晃荡进行了分析。Nagashima[8]应用有限元法并结合Level-set法研究了不同装载下液舱的晃荡并与理论和试验结果进行了对比。Cao等[9]应用SPH法对矩形液舱晃荡进行了研究。

对于油船、薄膜型LNG船而言，液舱作为主船体结构不可分割的一部分，其晃荡研究一般将液舱假定为刚性结构，将船体运动激励直接作用在刚性壁上，本文称之为非弹性支撑液舱。独立液舱的支撑结构相对于刚性液舱边界而言属于弹性体结构，船体运动激励需要通过弹性支撑结构传递到液舱内部，从而诱发舱室内部的液体晃荡运动，同时，作用于液舱边界的晃荡载荷又将反馈到弹性支撑结构上。因此，独立液舱的运动将是弹性支撑结构对外部激励的传递以及液舱晃荡共同作用的结果，其晃荡问题是考虑弹性支撑结构作用与液舱内流体冲击作用相互耦合的动力学问题。

本文针对带有弹性支撑的液舱晃荡问题，参考独立B型液舱结构形式，建立了弹性支撑下液舱作横摇运动的动力学方程，引入基于Hilber-Hughes-Taylor（HHT）[10]格式的隐式直接积分法求解液舱运动方程，运用VOF法[11]求解流体方程，并引入部分单元参数[12]处理棱形液舱边界，建立了液舱运动与流体晃荡双向耦合迭代算法。通过耦合计算研究弹性效应对液舱运动、晃荡液面变化以及晃荡载荷的影响规律。本文建立的方法适用于一般考虑弹性支撑的液舱晃荡问题。

1 数学方程和数值方法

1.1 液舱运动方程

带有弹性支撑的液舱晃荡模型如图1所示，液舱横摇运动动力学方程为如下形式：

式中：I0为体系的总惯性矩，Cθ为液舱转动阻尼，Kθ为液舱转动刚度，θ、、¨别为液舱的转角、角速度和角加速度，Mex为外部激励力矩，Mfluid为液体晃荡运动对液舱的作用力矩，Mg为液舱重力引起的力矩。

对船体施加正弦位移激励，将对船体的激励转化为对液舱的力矩：

式中：φ0为位移激励幅值；ωex为位移激励圆频率；β为相位差，β=arctan（ Cθωex/Kθ）。

流体作用在液舱上的力矩为流体在液舱壁以及内部构件压力积分的结果，其中压力由流体方程计算得到，液舱重力作用力矩为

图1 液舱及弹性支撑示意图Fig.1 Liquid tank with elastic support for sloshing analysis

（1）式为二阶非线性微分方程，运用t时刻液舱运动瞬时状态量通过流体方程求出通过（3）式求出Mg（）θ，代入（1）式，则原二阶非线性微分方程转换成t时刻二阶线性微分方程，采用HHT法进行离散求解，HHT法以Newmark法为基础，具有绝对稳定性，并且在保证二阶精度的基础上对数值结果引入数值阻尼，从而有效抑制了数值计算过程中尤其是初始状态的数值震荡。

1.2 流体晃荡控制方程

本文运用VOF法来计算流体晃荡，其通过引入体积分数F来追踪自由表面的变化，代表流体的体积与整个网格的体积比值。同时引入单元参数的概念，考虑体积通度λ以后，假设液体不可压缩，以随液舱运动的动参考系为坐标系统，流体相控制方程中连续性方程、动量方程和体积输运方程可表示为：

外部力包括重力和由于流体控制域的非匀速运动而引起的惯性力：

对于粘性流体，其在固壁上满足不可滑移条件和固壁不可穿透条件：

在自由表面上，流体必须同时满足运动学边界条件和动力学边界条件：

式中：un和ut分别为自由表面上的法向速度和切向速度，p0为舱室内蒸汽压力。

采用交错网格，压力P，体积分数F定义在网格中心，速度定义在网格的边界中心，采用有限差分法，均匀网格离散。对压力梯度▽p和扩散项采用中心差分离散，对于强非线性项▽·u→u→采用迎风格式与中心差分格式相结合的混合差分格式，保证稳定性和精度。

自由表面的重构是VOF中的重要部分，本文采用计算精度较高的Youngs法对自由表面进行重构[13]。计算中采用预测—矫正法解决不可压缩流动的压力—速度耦合问题[14-15]，在每个时间步内，首先忽略压力梯度项，通过动量方程求出预测速度，进而得到压力泊松方程，使用逐次超松弛迭代对该方程进行求解，最后用压力梯度对预测速度进行校正。

1.3 液舱运动与流体晃荡耦合计算

弹性支撑下液舱液体晃荡问题是双向运动与流体耦合问题。文中采用迭代耦合的方法进行求解，其主要思想是流体方程和运动方程分别按顺序相互迭代求解，在每一步将各自得到的结果提供给对方，直到耦合系统的解达到收敛，迭代停止。具体步骤如下，从液舱运动计算模块开始，将计算得到的角位移、角速度和角加速度传递给流体计算模块，流体计算模块由传递来的数据计算得到流体对液舱的力矩，将力矩传递给液舱运动计算模块进行下一步的计算，依次循环。

2 数值计算

2.1 物理模型

本文参考某一型LNG船独立B型棱形液舱的结构形式，取一段典型液舱截面作为研究对象，将液舱和船体之间的支撑结构等效成对应刚度的弹簧，编制二维弹性支撑下液舱运动与流体晃荡耦合计算程序，数值计算模型如图2所示。取液舱30%充装率，液面上P1和P2点为压力监测点，液体晃荡一阶共振周期为5.276 s，船体横摇运动固有周期为10.63 s，计算取结构化网格，网格尺寸为120×80，液体密度为500 kg/m3，其他主要参数见表1，其中外部激励作用在船体上，相当于弹簧的支座做强迫简谐位移运动。

表1 数值模拟主要参数Tab.1 Parameters of numericalmodel

图2 数值计算模型Fig.2 Numericalmodel for the simulation of sloshing

针对上述模型，将弹性支撑作用效果等效成作用于液舱上的横摇刚度，同理对应有液舱的等效横摇阻尼。在保证无因次阻尼系数不变的情况下改变等效横摇刚度大小，设定三个典型工况，如表2所示。其中，工况1为弹性支撑刚度较弱时的工况，工况2为系统接近共振时的工况，工况3为弹性支撑刚度较强时的工况。

2.2 弹性效应对液舱运动的影响

以接近共振状态的工况2为例，图3-5为液舱在非弹性支撑液舱充装30%液体、弹性支撑液舱充装30%液体（即为工况2）和弹性支撑空舱三种情况下液舱横摇角度、角速度和角加速度时间历程，其中上述三种情况外部激励参数均与表1中参数相同。通过分析发现，工况2条件下弹性支撑液舱相比非弹性支撑下液舱运动幅值显著增大，并产生相位差。弹性支撑下充装30%液体液舱运动幅值约为25.97°，弹性支撑下空舱的运动幅值约为34.38°，液舱中的液体晃荡抑制了液舱中的运动。图6为弹性支撑下装有30%液体液舱晃荡过程中激励矩、重力矩和液体对液舱作用矩的时间历程，由于相位差的原因，重力矩和液体作用矩减弱了船体激励矩对液舱的作用，从而抑制了液舱的运动，这正是利用液体晃荡进行减震的原理。

表2 工况设定Tab.2 Parameters of design conditions

图3 液舱横摇角度时间历程对比Fig.3 Time history of tank roll angle

图4 液舱横摇角速度时间历程对比Fig.4 Time history of tank roll angular velocity

图5 液舱横摇角加速度时间历程对比Fig.5 Time history of tank roll angular acceleration

图6 各成分弯矩Fig.6 Time histories of differentmoments

图7给出了不同弹性刚度条件下液舱运动时间历程，工况1船体激励通过弹性支撑传递到液舱后比非弹性支撑液舱的运动要小且呈包络曲线状，该工况下弹性支撑不足以完全传递船体运动激励，因此液舱的运动较小；工况2的参数设定使液舱的固有频率与激励频率相同，此时系统达到共振状态，液舱运动剧烈；工况3较工况2弹性支撑刚度更大，船体激励通过弹性支撑传递到液舱后比非弹性支撑液舱的运动要大，但小于工况2液舱的运动。

图7 不同弹性刚度下以及非弹性支撑下液舱运动时程曲线Fig.7 Time histories of tank rollmotion in different conditions

2.3 弹性效应对晃荡液面和载荷的影响

图8-9给出了工况2在不同时刻下非弹性支撑液舱和弹性支撑液舱液面运动模拟对比，由图8（a）和图9（a）可以看到非弹性支撑液舱的晃荡在t/T处于7.26～7.44时间段液面处于下降阶段，由较高液位恢复到接近于水平液位。图8（b）和图9（b）为弹性支撑下晃荡在不同时刻液面运动模拟，其中在t/T=7.26时刻左舱液面呈现明显的非线性波面特征，右舱液面接近静水液面；在t/T=7.44时刻左舱液面呈线性特性，右舱处于冲击作用阶段。由此可知，弹性效应作用下液面运动与非弹性支撑液面具有明显的区别，且左右舱晃荡呈现明显相位差。

图8 在t/T=7.26时刻液面运动模拟Fig.8 Liquid surfacemotion simulation at t/T=7.26

图9 在t/T=7.44时刻液面运动模拟Fig.9 Liquid surfacemotion simulation at t/T=7.44

比较不同弹性刚度以及非弹性支撑下液舱晃荡监测点压力时间历程，如图10所示，可以看到监测点P1在工况1的压力小于非弹性支撑下的压力，而工况2则相反，监测点P1压力较非弹性支撑下的压力大了一倍之多，也就是说弹性支撑的作用大大加强了液体晃荡对结构冲击的效应；工况3中监测点P1压力虽然依然比非弹性支撑下的压力大，但已远小于工况2的压力大小，由此可以看到弹性效应对晃荡载荷的影响规律与对液舱运动的影响规律基本一致。

以监测点P1为参考点，改变弹性支撑的刚度，并假设弹性支撑的变形量不受限制，即船体和液舱的相对运动角度可以足够大，计算对应刚度下P1点的晃荡压力，取晃荡处于稳定状态下的极值点，绘成曲线，如图11所示。可以看到，在初始阶段弹性支撑刚度较小时，有弹性支撑晃荡载荷比非弹性支撑晃荡载荷还要小，随着弹性支撑刚度的变大，液舱晃荡载荷呈逐渐增大趋势并逐渐会大于非弹性支撑下的晃荡载荷，直到系统达到共振状态，此时晃荡最为剧烈，随着刚度继续增大，晃荡压力逐渐减弱，当弹性支撑刚度足够大时，弹性支撑将船体激励完全传递到液舱上，此时独立液舱的晃荡与非弹性支撑液舱的晃荡趋于一致。进一步分析可以看到，当弹性支撑刚度趋于0时，外部激励将无法传递到液舱内部，此时将没有晃荡发生；当弹性支撑刚度无穷大时，此时的船体和液舱可以看成一个整体，液舱的运动将与船体运动保持一致；在系统共振点附近的一定范围内，晃荡载荷变化极其敏感，并且大于非弹性支撑下晃荡载荷的大小，最大可达非弹性支撑下晃荡载荷2倍大小，因此在设计中弹性支撑的刚度的选取应该避免这个范围。需要说明的是对于实际情况，一方面弹性支撑的变形是会受限制的，船体和液舱的相对运动角度不会很大；另一方面目前弹性支撑的选材和布置密度所形成的等效刚度一般位于图11曲线的后半程，此时弹性支撑下的液舱晃荡载荷相比于非弹性支撑下液舱晃荡载荷将略微增大。

图10 监测点P1压力时间历程对比 Fig.10 Pressure comparison of P1

图11 晃荡载荷随弹性支撑刚度变化曲线Fig.11 Sloshing load varied with the stiffness

3 结 论

本文针对带有弹性支撑的液舱晃荡问题，采用基于Hilber-Hughes-Taylor格式的隐式直接积分法求解液舱运动方程，应用含部分单元参数的VOF法求解流体方程，建立了弹性支撑条件下液舱运动与液体晃荡的耦合求解算法。通过耦合求解弹性支撑液舱晃荡，并与非弹性支撑液舱晃荡进行比较，得到如下结论：

（1）弹性支撑与非弹性支撑下液舱运动有明显区别，并随着弹性刚度的变化而变化，当激励频率达到液舱共振频率时，液舱运动最为剧烈；

（2）利用作用于液舱的外力矩与内部液体晃荡作用矩产生的相位差，舱内液体对液舱运动具有一定的抑制作用；

（3）弹性支撑下液面运动与非弹性支撑下液面运动相比有明显的区别，左右舱室晃荡会产生明显相位差；

（4）当弹性支撑刚度较小时，弹性支撑结构不足以传递船体的运动激励，液舱运动幅度较小，晃荡运动较平缓，晃荡载荷较小；随着弹性支撑刚度增大，液舱的自振频率接近船体激励频率时，液舱运动和舱内液体晃荡运动最为剧烈，此时晃荡载荷最大；当弹性支撑刚度继续增大时，液舱运动逐渐减小，晃荡剧烈程度逐渐降低，当弹性支撑刚度增加到足够大时，液舱运动与船体运动同步，此时独立液舱的晃荡趋近于非弹性支撑下液舱的晃荡。

本文建立的分析方法为含弹性支撑的独立液舱晃荡载荷预报提供了一种快速、有效的分析手段，同时也可为弹性支撑模块的材料选型和设计提供一定的参考。

[1]朱小松,谢 彬,喻西崇.LNG/LPG液舱晃荡研究进展综述[J].中国造船,2013,54(1):230-234. Zhu X S,Xie B,Yu X C.Research progress of liquid sloshing in LNG/LPG tanks[J].Shipbuilding of China,2013,54(1): 230-234.

[2]Faltinsen O M,Rognebakke O F,Timokha A N.Resonance three-dimensional nonlinear sloshing in a square-base basin [J].Journal of Fluid Mechanics,2003,487:1-42.

[3]Faltinsen O M,Rognebakke O F,Timokha A N.Transient and steady-state amplitudes of resonant three-dimensional sloshing in a square base tank with a finite fluid depth[J].Physics of Fluids,2006,18(1):1-14.

[4]Ibrahim R A,Pilipchuk V N,Ikeda T.Recent advances in liquid sloshing dynamics[J].Applied Mechanics Reviews, 2001,54(2):133-199.

[5]朱仁庆,马海潇,缪泉明,郑文涛.LNG船液舱晃荡压强预报[J].船舶力学,2013,17(1-2):42-48. Zhu R Q,Ma H X,Miao Q M,ZhengW T.Prediction of pressure induced by liquid sloshing for LNG carrier[J].Journal of Ship Mechanics,2013,17(1-2):42-48.

[6]朱仁庆,吴有生.液舱内流体晃荡特性数值模拟[J].中国造船,2002,43(2):15-28. Zhu R Q,Wu Y S.A numerical study on sloshing phenomena in a liquid tank[J].Shipbuilding of China,2002,43(2): 15-28.

[7]沈 猛,王 刚,唐文勇.基于改进VOF法的棱形液舱液体晃荡分析[J].中国造船,2009,50(1):1-9. Shen M,Wang G,TangW Y.Liquid sloshing analysis on prismatic tanks based on improved VOFmethod[J].Shipbuilding of China,2009,50(1):1-9.

[8]Nagashima T.Sloshing analysis of a liquid storage container using level set X-FEM[J].Communications in Numerical Methods in Engineering,2009,25(4):357-379.

[9]Cao X Y,Ming F R,Zhang A M.Sloshing in a rectangular tank based on SPH simulation[J].Applied Ocean Research, 2014,47:241-254.

[10]Hilber H M,Hughes T JR,Taylor R L.Improved numerical dissipation for time integration algorithms in structural dynamics[J].Earthquake Engineering and Structural Dynamics,1997,5:283-292.

[11]Hirt CW,Nichols B D.Volume of fluid(VOF)method of free boundaries[J].Journal of Computational Physics,1981,39: 201-225.

[12]Ingram D M,Causon D M,Mingham CG.Developments in Cartesian cut cellmethods[J].Mathematics and Computers in Simulation,2003,61:561-572.

[13]金 晶,薛鸿祥,唐文勇,张圣坤.大型补给舰船液舱晃荡载荷计算方法[J].中国舰船研究,2012,7(6):50-56. Jin J,Xue H X,Tang W Y,Zhang SK.Calculation procedure of sloshing loads for large-scale depot ships[J].Chinese Journal of Ship Research,2012,7(6):50-56.

[14]Griffith B E.An accurate and efficientmethod for the incompressible Navier Stokes equations using the projection method as a preconditioner[J].Journal of Computational Physics,2009,228:7565-7595.

[15]Brown D L,Cortez R,Minion M L.Accurate projectionmethods for the incompressible Navier Stokes equations[J].Journal of Computational Physics,2001,168(2):464-499.

Liquid sloshing analysis for independent tank considering the effect of elastic supports

LIUWen-fu,XUE Hong-xiang,TANGWen-yong
(1.State Key Laboratory of Ocean Engineering,Shanghai Jiao Tong University,Shanghai200240,China；2.Collaborative Innovation Center for Advanced Ship and Deep-Sea Exploration,Shanghai200240,China)

The independent tank is selected to investigate the sloshing characteristicswith the elastic support.A coupled algorithmis proposed involving the Hilber-Hughes-Taylor implicitmethod for tank motion and the VOFmethod for liquid sloshing.By conducting coupled calculation and comparing these numerical resultswith those obtained by ignoring the elastic impact,the connection between the stiffness of the elastic support and the liquid motion aswell as the sloshing load is analyzed.An efficient and effectivemethod is finally presented for predicting the sloshing load of independent tank considering the elastic support effect.

elastic support;independent tank;sloshing;VOFmethod

U661.32

：Adoi:10.3969/j.issn.1007-7294.2017.01.004

2016-06-20

工信部高技术船舶专项"新型液化天然气船液货围护系统预先研究"资助

刘文夫（1990-），男，硕士研究生；

薛鸿祥（1981-），男，博士，副教授，通讯作者，E-mail：hongxiangxue@sjtu.edu.cn。

1007-7294（2017）01-0023-08