0 引　言

近年来，我国天然气进口量已突破千亿立方米，单一的管道运输已经无法满足需求，而低温液化天然气的高效运输方式受到了更多关注。因此，对于LNG船的需求也不断上升，发展LNG船逐渐成为了航运的一大趋势^[1]。

作为典型的液货运输设备，LNG船液舱晃荡是设计和运维中不可忽略的重要问题，一直以来都受到学术界与工程界的广泛关注^{[2 − 4]}。尤其是薄膜型液舱，由于其特殊的外形和设计，在部分装载时容易产生强烈的晃荡冲击载荷，对舱壁的结构安全构成威胁。一旦结构发生破坏，泄漏的天然气可能会导致灾难性的事故^{[5 − 6]}。因此，液舱围护系统结构的安全设计和评估是LNG运输船设计的关键环节。

当前CFD方法和有限元方法在预报波浪冲击载荷及其引起的船体结构响应方面应用越来越广泛。You等^[7]采用CIP的方法，研究了孤立波作用于水下板的水动力特性。石晓等^[8]采用CFD方法，建立了气相可压缩水相不可压缩的波浪数值计算模型，模拟了波浪对建筑物冲击作用等强非线性流体运动问题。张友林等^[9]采用无网格粒子法MPS模拟了孤立波对平板的冲击过程，并计算了作用在平板上的波浪载荷。王佳东等^[10]采用CIP的方法，建立了二维高精度粘性流数值水槽模型，研究了大幅波浪对板式结构物的冲击破坏作用以及结构物与波浪的强非线性耦合作用。

在薄膜型液舱的结构分析中，Pillon等^[11]对薄膜型液舱进行静力分析，评估其力学性能。Arswendy等^[12]建立了NO96型液货船的有限元模型，研究材料性能和静力响应，验证了有限元方法的合理性。刘婷婷等^[13]采用有限元数值仿真方法系统地分析了准静力加载下围护系统可能的失效模式，得到了不同失效模式下的危险区域及对应的极限载荷，为液舱围护系统结构安全设计提供参考。骆松等^[14]建立了围护系统单元截面模型，确定围护系统结构的关键失效位置，研究了不同晃荡冲击载荷作用下的结构响应情况。徐家晨等^[15]对围护系统设计进行了有价值的探讨。

本文采用CFD方法及动网格技术建立数值水槽，基于UDF编写波浪场初始化程序；采用有限元法建立Mark Ⅲ型液舱围护系统的数值模型，然后基于Ansys System Coupling模块构建水槽波浪冲击围护系统结构的流固耦合力学模型。具体研究了围护系统受到的波浪砰击载荷时空分布特性及围护系统在砰击载荷下的可能失效模式，从而为Mark Ⅲ型液舱围护系统后续开展直线水槽波浪砰击试验奠定研究基础。

1 数值方法 1.1 流场数学模型

液舱内的晃荡液体满足连续性方程以及雷诺平均纳维-斯托克斯方程（Reynolds-Averaged Navier-Stokes Equations，RANS）。

连续性方程：

$ \frac{{\partial {u_i}}}{{\partial {x_i}}} = 0 。$

(1)

RANS：

$ \rho \frac{{\partial {{\bar u}_i}}}{{\partial t}} + \rho {\bar u_j}\frac{{\partial {{\bar u}_i}}}{{\partial {x_j}}} = \rho {\bar F_i} - \frac{{\partial \bar p}}{{\partial {x_i}}} + \frac{\partial }{{\partial {x_i}}}\left( {\mu \frac{{\partial {{\bar u}_i}}}{{\partial {x_j}}} - \rho \overline {u_i^{'}u_j^{'}} } \right)。$

(2)

式中：u_i为i方向上的速度；u_i 为平均速度；u_i'为脉动速度；p为微元所受总压强；F_i为微元所受质量力；μ为动力黏度；ρ为密度。

本文研究波浪砰击作用下围护系统的结构响应，属于典型的湍流问题，针对这类问题，一般采用k-ε模型、k-ω模型或LES模型，本文采用标准的k-ε模型计算，输运方程为：

$ \begin{split}\frac{\partial }{{\partial t}}\left( {\rho k} \right) + &\frac{\partial }{{\partial {x_i}}}\left( {\rho k{u_i}} \right) = \frac{\partial }{{\partial {x_i}}}\left[ {\left( {\mu + \frac{{{\mu _i}}}{{{\sigma _k}}}} \right)\frac{{\partial k}}{{\partial {x_i}}}} \right] +\\ &{G_k} + {G_B} + \rho \varepsilon - {Y_M} + {S_k} ，\end{split}$

(3)

$ \begin{split}&\frac{\partial }{{\partial t}}\left( {\rho \varepsilon } \right) + \frac{\partial }{{\partial {x_i}}}\left( {\rho \varepsilon {u_i}} \right) = \frac{\partial }{{\partial {x_i}}}\left[ {\left( {\mu + \frac{{{\mu _i}}}{{{\sigma _\varepsilon }}}} \right)\frac{{\partial \varepsilon }}{{\partial {x_i}}}} \right] + \\ &{C_{1\varepsilon }}\frac{\varepsilon }{k}\left( {{G_k} + {C_{3\varepsilon }}{G_b}} \right) - {C_{2\varepsilon }}\rho \frac{{{\varepsilon ^2}}}{k} + {S_\varepsilon } 。\end{split}$

(4)

式中：G_k和G_b分别是平均速度梯度和浮力引起的派生项；C_1ε、C_2ε、C_3ε、σ_k和σ_ε和为经验常数；μ_i为湍动粘度。

1.2 自由液面捕捉技术

为了更好地模拟流体在计算过程中的各种非线性现象，本文采用VOF法来处理自由液面的追踪问题。VOF法是通过研究网格单元内流体和网格体积分数比来追踪自由液面的变化，其控制方程如下：

$ \frac{\partial a_w}{\partial t}+\boldsymbol{U}\cdot\nabla a_w=0，$

(5)

$ {a_w} + {a_a} = 1。$

(6)

式中：α_w和α_a为水和空气的体积分数；U为速度矢量，且当α_w=1时，单元内都为水，α_a=1时，单元内都为空气，0<α_w或α_a<1，则表示单元内是多相流体的混合，即为交界面。

1.3 结构场数学模型

假定结构为线弹性体，结构动力学控制方程的有限元表达式为：

$ {\boldsymbol{M\ddot U + C\dot U + KU = Q}}{\text{(}}t{\text{)}} 。$

(7)

式中：$ {\boldsymbol{\ddot U}} $、$ {\boldsymbol{\dot U}} $和$ {\boldsymbol{U}} $分别为系统的加速度向量、速度向量和节点位移向量；$ {\boldsymbol{M}} $、$ {\boldsymbol{C}} $和$ {\boldsymbol{K}} $分别为系统的质量矩阵、阻尼矩阵和刚度矩阵；$ {\boldsymbol Q}(t) $为系统受到的外部结点载荷向量。

本文采用Newmark法对式(7)求解。将二阶微分方程组降阶为线性方程组：

$ \begin{split} &\left[ {{\boldsymbol{K}} + \frac{1}{{\alpha \Delta {t^2}}}{\boldsymbol{M}} + \frac{\delta }{{\alpha \Delta t}}{\boldsymbol{C}}} \right]{{\boldsymbol{U}}_{{\boldsymbol{n + }}{\text{1}}}} = {{\boldsymbol{Q}}_{{\boldsymbol{n + }}{\text{1}}}} + \\ &{\boldsymbol{M}}\left[ {\frac{1}{{\alpha \Delta {t^2}}}{{\boldsymbol{U}}_{\boldsymbol{n}}} + \frac{1}{{\alpha \Delta t}}{{{\boldsymbol{\dot U}}}_{\boldsymbol{n}}} + \left( {\frac{1}{{2\alpha }} - 1} \right){{{\boldsymbol{\ddot U}}}_{\boldsymbol{n}}}} \right] +\\ &{\boldsymbol{C}}\left[ {\frac{\delta }{{\alpha \Delta }}{{\boldsymbol{U}}_{\boldsymbol{n}}} + \left( {\frac{\delta }{\alpha } - 1} \right){{{\boldsymbol{\dot U}}}_{\boldsymbol{n}}} + \left( {\frac{\delta }{{2\alpha }} - 1} \right)\Delta t{{{\boldsymbol{\ddot U}}}_{\boldsymbol{n}}}} \right]。\end{split}$

(8)

式中：$ \alpha $和$ \delta $为积分参数，与积分精度和算法稳定性要求有关。当$ \alpha $和$ \delta $取不同数值时，将得到不同的积分方案。

1.4 流固耦合边界条件

在Ansys双向流固耦合理论中，位移连续性和力守恒是确保流体和固体相互作用正确模拟的2个关键原则，其控制方程如下：

位移连续性方程：

$ {u_f} = {u_s} 。$

(9)

力守恒方程：

$ {t_f} = {t_s} 。$

(10)

式中：$ {u_f} $和$ {u_s} $分别为流体在流固界面上的位移向量和固体在流固界面上的位移向量；$ {t_f} $和$ {t_s} $分别为流体在流固界面上的应力向量和固体在流固界面上的应力向量。

2 Mark Ⅲ型LNG船围护系统结构

Mark Ⅲ型是薄膜型 LNG 船围护系统结构的一种主要类型，其结构如图1所示。其隔热层是由2层加强聚氨酯泡沫层和中间次屏壁薄膜层组成。次层隔热层由整块聚氨酯泡沫预制，通过后层合板和树脂绳安装在船舶内壳。首层隔热层和顶层合板被预制插槽用于安装首层不锈钢薄膜层。为了提高隔热层力学性能，在聚氨酯泡沫内加人了玻璃纤维，将次屏壁薄膜层制成由玻璃纤维和铝片组成的3层膜结构，同时为了克服金属层的温度应变，将首层屏障薄膜外形设计为波纹。

组成Mark Ⅲ型隔热层的主要材料是层合板、聚氨酯泡沫及树脂绳。在本文有限元模拟中假定层合板和聚氨酯泡沫遵循正交各向异性弹性规律，树脂绳为各向同性弹性材料。有限元模型尺寸如图2所示，材料参数如表1所示。

3 数值计算模型 3.1 数值水槽计算模型

图3为本文计算模型示意图，其中水槽高度H为10.6 m，长度B为20 m，围护系统距水槽底部h₁为 2.8 m，围护系统模型高h₂为2.38 m，围护系统的厚度b为0.421 m。

通过孤立波来砰击围护系统波纹面，孤立波波形通过UDF初始化加载。孤立波波面方程如下：

$ \eta = D{\text{sec}}{{\text{h}}^2}\left[ {\sqrt {\frac{{3D}}{{4d}}} \left( {\frac{x}{d} - \frac{{ct}}{d}} \right)} \right] ,$

(11)

$ c = \sqrt {g\left( {d + h} \right)}。$

(12)

式中：x为波浪传播方向上任一点；D为孤立波波高；d为水深；c为孤立波波速；g为重力加速度。波高为2.8 m，水深设为4 m，重力加速度设为9.81 m/s²，流体密度设为998 kg/m³。

如图4所示，流体计算域采用CFD类型的网格对其进行划分，通过multizone方法划分出部分三棱柱网格，用其进行过渡，其余部分采用六面体网格，这样不仅可以提高网格质量降低网格量，还提高了计算效率。

整体流体计算域网格共划分单元数700502个，节点数445839个，且正交质量在0.2以上。流体模块，压力速度耦合选SIMPLE算法，压力选择Body Force Weighted方法离散。时间步取0.0001 s，每个时间步迭代30次。

如图5所示，对于流体计算域，其中EH为自由液面，液面以上部分为空气，液面以下部分为水。上边界AB、右边界BG：设为压力出口条件；下边界CD、右边界DF：设为固壁条件；左边界AC：除波浪速度入口以外，其余均设为固壁条件；IJ面设为流固耦合面。由于围护系统结构在流场中运动，会引起流体网格的变形甚至畸变，从而导致流场的计算结果精度降低甚至计算中断，因此流体域数值模型需要使用动网格技术，本文采用的是光顺模型。

3.2 围护系统有限元模型

如图6所示，固体计算域采用Mechanical类型的网格进行划分，共划分单元数136094个，节点154981个，且平均单元质量在0.6以上。

固体计算域将主屏壁层下表面与层合板上表面设置为绑定接触，次屏壁层与主、次绝热层之间也设置绑定接触，注意接触面方向设置。根据围护系统结构的安装方式，在船底板的底部施加固定约束。

4 结果与讨论 4.1 刚性与弹性围护系统的砰击载荷分析

本文研究了2组工况，一组是刚性围护系统，分析孤立波砰击围护系统表面的冲击载荷时空分布特性；另一组是弹性围护系统，分析弹性效应对砰击压力的影响以及在该砰击压力下结构的变形及可能的失效情况。

在围护系统波纹面，即流固耦合面选取6个关键位置，来监测波浪冲击围护系统结构的压力，具体监测位置如图7所示。

图8为P₁、P₂、P₃位置受到的波浪作用到围护系统表面的压力时历。由图8(a)可知，无论围护系统是刚性或弹性结构，流体载荷均具有冲击特性，而且刚性围护系统受到的冲击载荷提前于弹性系统。但是对于如图8(b)所示的波节相交的角隅附近位置，刚性围护系统受到的流体载荷表现为明显的冲击特征，而弹性围护系统受到的流体载荷则没有冲击特征。

由图8(c)可知，在该位置无论围护系统是刚性或弹性结构，流体载荷均没有明显的冲击特征。这是因为波浪垂直砰击自由液面位置附近的围护系统后，沿着舱壁爬升到该位置附近，液面运动方向几乎与围护系统表面平行，故没有产生明显的冲击效应。

图9为P₄位置受到的波浪作用到围护系统表面的压力时历。由图9(a)可知，围护系统刚性和弹性工况下，P₄位置波浪作用于围护系统壁面引起明显的冲击载荷，且均表现出双峰特点。液体猛烈砰击竖直壁面，产生持续时间短、幅值大的主峰；次峰可能由于波浪跨越波纹节引起的压力波动所致，其峰值较小、持续随时间较长。

图9(a)表明弹性围护系统的冲击压力要小于刚性情况。这是因为在外界激励条件相同的情况下，弹性围护系统结构能够吸收波浪冲击过程中的能量，从而产生运动和形变，把液体的动能转化为弹性结构的弹性势能，并且在这个过程中由于摩擦力等阻尼因素，一部分机械能将被转化成内部损耗能，从而减小冲击压力。

由图9(b)可知，相较于弹性结构，刚性结构的冲击载荷轮廓更加光滑，且主峰对应的时刻提前于弹性围护系统情况。弹性围护系统对应的砰击载荷轮廓中除了明显的主峰外，还有较多波动的次峰，这由于波浪冲击引起的冲击压力作用在波纹板上导致其振动，而板的振动又反过来影响冲击压力，从而导致冲击压力的高频振动，这也说明流固耦合效应对砰击载荷特性的影响是明显的。

图10为P₅位置受到的波浪作用到围护系统表面的压力时历。由图10(a)可知，刚性与弹性围护系统的砰击载荷均有明显的冲击特征，且具有次峰；刚性情况下的砰击压力幅值大于弹性情况；在上升时间方面，两者基本一样。由图10(b)可知，刚性围护系统的冲击载荷轮廓较为光滑，而弹性情况的冲击载荷在上升阶段有多个压力波动，这种压力波动是由于结构的振动引起。它们的次峰远低于主峰值，这是由于波浪爬升波纹节引起的压力波动造成的。

图11为P₆位置受到的波浪作用到围护系统表面的压力时历。由图11(a)可知，刚性与弹性围护系统的砰击载荷均有明显的冲击特征，且具有次峰；刚性情况下的砰击压力幅值大于弹性情况。由图11(b)可知，相较于弹性结构，刚性结构的冲击载荷轮廓更加光滑，且主峰对应的时刻提前于弹性围护系统情况。弹性围护系统对应的砰击载荷轮廓中除了明显的主峰外，还有较多的波动的次峰，这是由于波浪冲击引起的冲击压力作用在波纹板上导致其振动，而板的振动又反过来影响冲击压力，从而导致冲击压力的高频振动，这也说明流固耦合效应对砰击载荷特性的影响是明显的。结合图9、图10可知，最大冲击压力均发生在波纹节中间位置且迎浪方向的平区附近。

图12为刚性和弹性围护系统在相同时刻自由液面对比；图13为弹性围护系统在t = 0.024 s时刻下的自由液面；图14为刚性和弹性围护系统在相同时刻压力对比。由图12、图13可以看出，自由液面猛烈地砰击竖直舱壁的局部位置，发生了飞溅及破碎现象，但是围护系统无论是刚性还是弹性，自由液面运动没有观察到明显差别。结合图11，在出现最大砰击压力时刻附近，即在0.012～0.020 s自由液面沿壁面爬升，在很小的局部区域液体接触到壁面，这些位置是产生大的冲击压力的区域，而更多的区域是液体包围着气体，气体作用于壁面。

4.2 Mark Ⅲ围护系统结构响应

根据对砰击载荷的分析，在t = 0.0172 s左右时，水槽内冲击压力达到最大，约为0.6 MPa。在该时刻，围护系统的结构响应云图如图15、图16所示。

分析可知，在砰击压力达到最大时，围护系统的最大应力响应发生在靠近波纹节附近，是由几何突变引起的应力集中现象，应力可达604 MPa。由图15(b)的围护系统结构变形剖面图可以看出，围护系统结构的变形主要在载荷冲击方向，最大变形出现在顶层合板上（最大值为3.4 mm），主要变形在顶层合板及首层聚氨酯泡沫隔热层上。这应该引起足够的重视，因为隔热层损坏，引起LNG相变、内压增大，进而加剧结构的失效。

由于波浪冲击载荷作用时间比较短，导致结构响应还没有传递至下层层合板时其作用过程就己经完成，因此动力作用下结构的应力及位移最大值主要分布在主屏壁层及首层层合板中。通过模拟可以看出，此时在波浪冲击作用下，顶层合板沿垂向方向发生较大的变形，这对于主屏壁层的安全性是极为不利的。

5 结　语

1）在波节相交的角隅附近位置，刚性围护系统受到的流体载荷表现为明显的冲击特征，而弹性围护系统受到的流体载荷则没有冲击特征。

2）流固耦合效应对砰击载荷特性的影响是明显的。在波纹板中间位置且迎浪方向的平区附近，刚性与弹性围护系统的砰击载荷均有明显的冲击特征，刚性围护系统的冲击载荷轮廓较为光滑，而弹性围护系统对应的砰击载荷轮廓中存在较多波动的次峰。

3）在给定波浪条件下，水槽内的冲击压力达到最大时，弹性围护系统结构的变形主要发生在顶层合板及首层聚氨酯泡沫隔热层上（最大值为3.4 mm），而最大应力响应发生在靠近波纹节附近，应力可达604 MPa。