0 引　言

随着各国对清洁能源的需求不断增加，全球范围内天然气运输量也在快速增长。与气态形式相比，LNG能量密度高、易于储存，非常适合陆上或海上长距离运输。LNG船液货舱通常可分为2种类型：薄膜型和独立型^[1]。独立型液货舱又可以分为A型、B型和C型^[2]。B型液货舱的货物维护系统主要包括液货舱和支撑结构，其特点在于液货舱具有部分次屏蔽且不作为压力容器。液货舱和船体通过弹性支撑结构连接。船体运动时，支撑可将作用力传递给液货舱，引起液货舱运动，进而激发舱内液体的晃荡。此时，晃荡载荷不仅会威胁液货舱的结构安全，也会使支撑结构受到的作用力增大。因此，在设计阶段，必须考虑晃荡载荷对支撑结构作用力的影响。

随着计算机性能的提高和相关理论的发展，越来越多的数值方法和求解技术被应用于液舱晃荡的研究^[3]，如有限元法（Finite Element，FE）^[4]，流体体积法（Volume of Fluid，VOF）^[5]，任意拉格朗日-欧拉法（Arbitrary Lagrangian-Eulerian，ALE）^[6]等。

本文以1艘LNG船B型独立液货舱为研究对象，建立二维模型，基于ALE方法和体积模量缩减技术，运用有限元软件Abaqus对液舱晃荡进行分析，得到横摇状态下支撑结构作用力的时间响应。讨论了船体横摇幅值、货物状态以及液货舱挡板布置对支撑结构作用力的影响。

1 研究对象

B型液货舱按形状可分为棱柱形（SPB型）和球形（MOSS型），如图1所示^[7]。本文的研究对象是棱柱形液货舱，主要由耐低温的平面板材制成。

支撑结构用于连接船体和液货舱。根据其安装位置和功能的不同，支撑结构可分为垂向支撑、防横摇支撑、防纵摇支撑以及止浮支撑等，如图2所示。垂向支撑和止浮支撑用于支持液货舱，防止其与船体发生碰撞，可提供y方向的作用力，且只受压不受拉。防横摇支撑和防纵摇支撑用于减轻由船体运动引起的液货舱晃动，可提供x或z方向的作用力。

2 数值仿真 2.1 数值方法

本文利用有限元软件Abaqus的显式动态分析模块（Abaqus/Explicit，Dynamic），采用任意拉格朗日-欧拉自适应网格技术（Arbitrary Lagrangian-Eulerian adaptive meshing technique）和体积模量缩减技术（Bulk modulus reduction technique），对B型独立液货舱的晃荡问题进行研究。

使用有限元法对液舱晃荡进行分析时，若流体单元变形过大则极易导致计算不收敛。ALE自适应网格技术可以有效解决这一问题，它结合了拉格朗日方法（网格点与物质点重合，网格之间没有物质的转移）和欧拉法（网格点固定，物质可以在网格间自由转移）的优点，能够在不改变网格原有拓扑结构的情况下保证网格质量，增大计算结果的收敛性和准确性。ALE方法的物质导数描述为^[8]：

$ \frac{{\partial f\left( {{X_i},t} \right)}}{{\partial t}} = \frac{{\partial f\left( {{x_i},t} \right)}}{{\partial t}} + \frac{{\partial f\left( {{x_i},t} \right)}}{{\partial {x_i}}}{w_i} {\text{。}}$

(1)

其中： ${X_i}$ 为拉格朗日坐标； ${x_i}$ 为欧拉坐标； ${w_i}$ 为物质速度 $v$ 和网格速度 $u$ 之间的相对速度。

该方法的质量守恒、动量守恒和能量守恒控制方程分别为：

$ \frac{{\partial \rho }}{{\partial t}} = - \rho \frac{{\partial {v_i}}}{{\partial {x_i}}} - {w_i}\frac{{\partial \rho }}{{\partial {x_i}}}{\text{，}} $

(2)

$ \rho \frac{{\partial {v_i}}}{{\partial t}} = {\sigma _{ij,j}} + \rho {b_i} - \rho {w_i}\frac{{\partial {v_i}}}{{\partial {x_j}}} {\text{，}}$

(3)

$ \rho \frac{{\partial E}}{{\partial t}} = {\sigma _{ij,j}}{v_{i,j}} + \rho {b_i}{v_i} - \rho {w_j}\frac{{\partial E}}{{\partial {x_j}}}{\text{。}} $

(4)

其中： $\rho $ 为流体密度； ${\sigma _{ij}}$ 为应力张量； $b$ 为单位质量体积力； $E$ 为单位质量总能量。

本文将LNG假定为无粘性不可压缩流体，在Abaqus中使用Mie-Grüneisen状态方程对其建模。一种常见的Mie-Grüneisen状态方程表达为^[9-10]：

$ p - {p_{ref}} = \mathit{\Gamma} \rho \left( {e - {e_{ref}}} \right){\text{。}} $

(5)

其中： $p$ 为压力； $\rho $ 为密度； $e$ 为比内能； ${p_{ref}}$ 和 ${e_{ref}}$ 为仅与密度 $\rho $ 有关的压力和比内能的参考值。 $\mathit{\Gamma} $ 为Grüneisen系数，有如下表达：

$ \Gamma = \frac{1}{\rho }{\left. {\frac{{\partial p}}{{\partial e}}} \right|_\rho }{\text{。}} $

(6)

Abaqus引入了线性 ${U_S} - {U_P}$ Hugoniot 关系，可将 ${p_{ref}}$ 改写为^[11-13]：

$ {p_{ref}} = \frac{{{\rho _{ref}}c_0^2\eta }}{{{{\left( {1 - s\eta } \right)}^2}}}{\text{。}} $

(7)

其中： ${\rho _{ref}}$ 为材料密度； $\eta = 1 - {\rho \mathord{\left/ {\vphantom {\rho {{\rho _{ref}}}}} \right.} {{\rho _{ref}}}}$ 是名义体积压缩应变； ${c_0}$ 为材料中声速； $s$ 为材料常数， ${c_0}$ 和 $s$ 用于表达冲击波速度 ${U_S}$ 和粒子速度 ${U_P}$ 之间的线性关系，即

$ {U_S} = {c_0} + s{U_P}{\text{。}} $

(8)

在小应变情况下，有下式成立：

$ K = {\rho _{ref}}c_0^2{\text{，}} $

(9)

其中， $K$ 为材料的体积模量。

数值分析中采用体积模量缩减技术能够在保证结果精度的情况下缩短计算时间。有限元分析时，稳态时间增量步长 $\Delta t$ 约等于通过网格中任意单元的膨胀波的最小传播时间，即

$ \Delta t \approx \frac{{{L_{\min }}}}{{{c_d}}}{\text{。}} $

(10)

其中： ${L_{\min }}$ 为网格中单元最小尺寸； ${c_d}$ 为膨胀波速度，由下式计算^[14]：

$ {c_d} = \sqrt {\frac{{E\left( {1 - \mu } \right)}}{{\rho \left( {1 + \mu } \right)\left( {1 - 2\mu } \right)}}}{\text{。}} $

(11)

考虑到体积模量的表达式为：

$ K = \frac{E}{{3\left( {1 - 2\mu } \right)}}{\text{，}} $

(12)

将式(12)代入式(11)，得到：

$ {c_d} = \sqrt {K\frac{{3\left( {1 - \mu } \right)}}{{\rho \left( {1 + \mu } \right)}}}{\text{。}} $

(13)

可以看出，若减小材料的体积模量 $K$ ，则降低了膨胀波速度 ${c_d}$ ，从而增大了稳态时间增量步长 $\Delta t$ ，进而缩短了计算时间。对于本文研究的晃荡问题，可以将液体的体积模量减小2～3个数量级，且不会影响液体的流动行为^[11]。Abaqus中，根据式(9)，在液体密度一定的情况下，可以通过设定不同的声速值实现体积模量的缩减。

2.2 方法验证

为验证上述数值方法的准确性，选取参考文献[15]中的二维矩形容器晃荡算例进行分析，如图3所示。

该矩形容器在水平方向外力作用下发生受迫振动，其水平加速度可表示为：

$ {a_x}\left( t \right) = - {X_0}{\omega ^2}{\rm sin} \omega t{\text{。}} $

(14)

其中：加速度幅值X₀=0.002 m；角速度ω=5.5 rad/s；重力加速度g=9.8 m/s²。容器内液体为水，密度 ${\rho _{water}}$ =983.2 kg/m³，真实体积模量 ${K_{water}} = 2.07 \times {10^9}\;\rm{Pa}$ 。t=0时液体完全静止。计算得到自由液面最右端A点的液面高度随时间的变化。为考察体积模量缩减对晃荡分析的影响，选取4个声速不同的模型进行计算，将结果与文献数据进行对比，如表1和图4所示。

可以看出，随着体积模量的缩减，计算时间也在减少。将体积模量缩减为真实值的1/‰，即 $ {K_{water}} =$ $ 2.07 \times {10^6}\;\rm{Pa}$ ，得到的结果依然与文献数据吻合较好。因此，采用ALE自适应网格技术和体积模量缩减技术能够对液舱晃荡问题进行有效分析，并在保证结果精度的同时缩短计算时间。

2.3 有限元模型

本文以1艘典型LNG船B型独立液货舱为例，建立横剖面二维模型。液货舱宽11.0 m，高6.0 m，位于中线处的纵向水密舱壁将其分为2个体积相等的舱室。液货舱主体结构，如外板、纵向舱壁、桁材等，均由刚体单元（R2D2）建模。模型包含5个支撑结构，中纵舱壁底部有1个防横摇支撑，液货舱底部关于中纵舱壁对称分布有4个垂向支撑。为方便计算，将这些支撑简化为弹性连接，由连接单元（CONN2D2）建模。防横摇支撑具有x方向刚度，垂向支撑具有y方向受压刚度，刚度系数取值见表2。有限元模型如图5所示。其中，VS1～VS4表示垂向支撑，AR1表示防横摇支撑。

液货舱装载率为40%，LNG由平面应变单元（CPE4R）建模。对自由表面附近的网格进行加密处理。不考虑液体和液货舱之间的摩擦力。LNG物理属性如表3所示^[16]。

2.4 计算工况

从运动状态、货物状态和液货舱内有无挡板3个方面考虑，选择7种计算工况，如表4所示。

液货舱按如下正弦规律做横摇运动：

$ \varphi = - {\varphi _0}\sin \omega t{\text{。}} $

(15)

其中，横摇幅值为φ₀，角速度ω=0.748 rad/s，相应的运动周期T=8.4 s，为获得稳定结果，计算时长选为5个周期。

为验证连接单元的准确性与可靠性，计算LC1工况，比较支撑结构作用力与液货舱自身重力。为分析不同横摇幅值对支撑结构作用力的影响，分别对φ₀=5°和φ₀=8°两种情况下的液舱晃荡进行计算。为分析液体和固体货物对支撑结构作用力的影响，计算LC4和LC7工况，货物状态为液体表示装载LNG，货物状态为固体表示装载与LNG具有相同密度和体积的固体，计算时当作刚体处理。为分析液货舱内挡板对支撑结构作用力的影响，分别计算有挡板和无挡板2种情况。

3 计算结果

对于LC1工况，液货舱静止时，仅有4个垂向支撑受到y方向的作用力。将计算得到的连接单元作用力与液货舱重力进行对比，结果如表5所示。

4个垂向支撑的总支持力与液货舱重力的误差小于1%，说明连接单元能够较好的模拟支撑结构，传递作用力。

图6和图7分别对比了无挡板和有挡板的液货舱模型在不同横摇幅值下支撑结构的受力情况。可以看出，支撑结构的作用力随着船体横摇呈现出周期性变化。对于垂向支撑VS，液货舱静止时平均作用力 ${F_{VS0}} = - 31\;538.5\; {\rm N}$ ，横摇运动时作用力以 ${F_{VS0}}$ 为平衡位置，在其附近波动，周期与横摇运动相同。远离液货舱中心线的VS1的作用力变化范围明显大于靠近中心线的VS2。对于防横摇支撑AR，液货舱静止时平均作用力 ${F_{AR0}} = 0$ ，横摇运动时作用力以 ${F_{AR0}}$ 为平衡位置，波动规律与垂向支撑类似。对于所有支撑，横摇幅值越大，支撑结构作用力的变化范围越大，作用力绝对值的最大值越大。

图8和图9分别对比了横摇幅值不同的液货舱模型在无挡板、有挡板以及装载固体货物时支撑结构的受力情况。可以看出，无论液货舱内有挡板还是没有挡板，装载液体还是固体，支撑结构作用力均具有类似的周期性变化规律。对于所有支撑，液货舱无挡板时支撑结构作用力变化范围最大，有挡板时次之，装载固体货物时最小。这说明，由于晃荡效应的存在，相对于固体货物，液体货物使支撑结构承受了更大的作用力。同时，在液货舱设置挡板能够有效降低作用力变化范围。

4 结　语

本文利用有限元法，采用任意拉格朗日-欧拉自适应网格技术和体积模量缩减技术，对LNG船B型独立液货舱进行晃荡分析，并探讨了晃荡对支撑结构作用力的影响。通过与文献数据对比，验证了数值方法的可靠性；建立有二维限元模型，根据横摇幅值、液货舱布置以及货物形态选取7个工况进行计算；比较不同工况下支撑结构作用力随时间的变化，得到以下结论：

1）船体做周期性横摇运动时，支撑结构作用力也呈现出周期性变化，且在平衡位置附近波动；

2）船体横摇幅值越大，支撑结构作用力的变化范围越大；

3）由于液舱晃荡，装载液体货物时支撑结构作用力的变化范围大于装载相同密度和体积的固体货物时的变化范围；

4）液货舱内设置挡板能够降低支撑结构作用力的变化范围。