0 引　言

1992年，国际海事组织（International Maritime Organization，IMO）修订了《国际防止船舶造成污染公约》（MARPOL公约）^[1]，要求油船采用双壳船体结构。到2006年，实施《IACS双壳油船共同结构规范》^[2]，自此，全球远洋油船逐步禁止使用单壳结构，双壳油船逐渐占据主流。双壳油船破损后，外界水体会通过船体外壳破裂处涌入船体内，轻则影响油船正常行驶，严重的会导致倾覆。虽然双层壳设计在一定程度上大大提高了油船的破损稳性，但仍有必要对双壳油船的破损进水过程及运动响应进行研究。

对破损船舶进水过程及运动响应的研究主要集中在理论研究、模型实验和数值模拟3个方面。在理论研究方面，卜淑霞等^[3]提出了一种基于势流理论、利用修正的伯努利方程模拟船舶破损进水的方法，并以一艘ITTC驳船标模为例，研究了不同破口位置对破损船体运动响应的影响。Zhang等^[4]应用改进的MPS方法，对破损舱室的瞬态进水过程进行了数值模拟，分析了破口位置和舱内挡板对破损舱室进水过程和运动特性的影响。Gao等^[5]提出了一种基于势流理论的耐波性求解器与基于VOF模型的Navier-Stokes（NS）求解器耦合的数值模拟方法，并模拟了一艘破损客滚船的横摇衰减和在规则横浪下的横摇运动。在模型实验研究方面：Begovic等^[6]在波浪水池中对完好和破损两种状态的DTMB-5415军舰模型在不同浪向下的运动响应特性进行了研究。Siddiqui等^[7]在实验水池中对完好和破损两种状态的箱型体模型在横浪下进行了自由漂浮试验，研究了波浪周期、初始装载条件、破口尺寸等对受损模型运动的影响。在数值模拟研究方面，陈波^[8]运用STAR-CCM+软件，研究了一艘油船在不同浪向下破损后的油水交换过程以及船体横摇、垂荡等运动的变化。顾乾乾^[9]以装有货油的双层壳箱体和单层壳箱体为对象，数值模拟了破口大小和位置不同时原油的泄漏过程，分析了双层壳对舱室油品泄漏速度和泄漏量的影响。厉京达等^[10]以ITTC标准驳船模型为研究对象，分析了空气压缩效应对破损舱室内部压强、进水速度等因素的影响。Liu等^[11]研究了舱室不同破口位置、破口尺寸以及舱室内部不同挡板布置形式对破损舱室进水过程的影响。Ming等^[12]研究了破损船舶在横浪海况下，波浪前进方向正对破口和背对破口对破损舱室进水过程的影响，发现波浪前进方向背对破口时，破损舱室的进水速度更快，进水量更大。Zhang等^[13]研究了破损舱室内部横纵舱壁在不同开孔布置形式下对破损舱室进水过程的影响。李月萌等^[14]基于Fluent对一艘舷侧破损的箱型船进行了时域模拟，研究了该船破损后的进水过程。王世取等^[15]对一艘破损舰船在有航速和无航速下的进水过程进行了研究，研究发现，相较于无航速状态，破损舰船在有航速状态下的破口进水速度更慢，但进水量更多。董良志等^[16]运用STAR-CCM+软件，对一艘破损邮轮单舱进水和多舱进水工况进行了模拟，并对破损进水后邮轮的横摇、纵摇以及垂荡运动进行了对比分析。于博文等^[17] 运用该软件研究了一艘破损渔船在不同海况下的瞬时进水过程，分析了渔船破损进水后的运动响应和压强分布。

综上所述，目前对破损船舶进水过程及运动响应的研究多数都集中在单壳船舶，或简化船体双壳结构为单壳等效模型，对双壳船舶破损的研究较少，鲜有人对双壳船单壳破损情况进行研究，更鲜有对单、双壳破损后船体性能的对比分析。但在实际工程应用中，单、双壳破损情况下船体的进水过程和运动响应具有差异性，对制定不同抗沉决策具有重要影响。因此，有必要对双壳船舶单、双壳破损进水过程及运动响应进行研究和对比分析。

本文以一艘113300 DWT双壳油船为研究对象，采用（Volume of Fluid, VOF）流体体积法捕捉破损舱室进水过程，采用（Dynamic Fluid Body Interaction, DFBI）模型与重叠网格技术模拟船体在水面上的运动，运用STAR-CCM+软件分别对该船右舷第三货舱段的单舱单壳破损和单舱双壳破损工况进行数值模拟，分别探究了这两种破损工况在静水、尾随浪和横浪3种海况下的进水过程和运动响应，并将结果进行了综合对比分析，期望该研究结果能为双壳油船的抗沉决策提供一定的参考。

1 数值模型 1.1 RANS方程

RANS方法是流体力学中最常用的数值模拟方法之一，其基本原理是通过对流体的运动方程进行求解，得到流体的速度、压力和温度等物理量的分布情况，该方法的雷诺时均连续性方程和动量方程^[18]如下。

雷诺时均连续性方程：

$ \frac{{\partial \rho }}{{\partial t}} + \frac{\partial }{{\partial {x_i}}}\left( {\rho {{\overline u }_i}} \right) = 0。$

(1)

动量方程：

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

(2)

式中：$\rho $为流体的密度；t为时间；${\overline u _i}$、${\overline u _j}$为流体的时均速度；$u_i^{'}$、$u_j^{'}$为流体的脉动速度；$\overline p $为流体的平均压力；$\mu $为流体的动力粘度系数；${\overline F _i}$为体积力。

1.2 湍流模型

本文使用Jones等提出的双方程模型的标准k-epsilon模型^[19]。k-epsilon模型的湍流控制方程为：

${\begin{split} \frac{{\partial \left( {\rho k} \right)}}{{\partial t}}+ &\frac{{\partial \left( {\rho ku_i} \right)}}{{\partial {x_i}}} = \frac{\partial }{{\partial {x_j}}}\left[ \left( {\mu + \frac{{{\mu _t}}}{{{\sigma _k}}}} \right)\frac{{\partial k}}{{\partial {x_j}}} \right]+ \\ &{G_k}+ {G_b}-\rho \varepsilon -{Y_M} + {S_k}，\\ \frac{{\partial \left( {\rho \varepsilon } \right)}}{{\partial t}} + & \frac{{\partial \left( {\rho \varepsilon {u_i}} \right)}}{{\partial {x_i}}}=\frac{\partial }{{\partial {x_j}}}\left[ {\left( {\mu + \frac{{{\mu _t}}}{{{\sigma _\varepsilon }}}} \right)\frac{{\partial \varepsilon }}{{\partial {x_j}}}} \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}}$

(3)

式中：${G_k}$为速度梯度产生的湍动能项；${G_b}$为浮力产生的湍动能项；$i$,$j$=1、2、3分别表示$x$、$y$、$z$方向；${Y_M}$为脉动扩张项；$ {C_{1\varepsilon }} $、$ {C_{2\varepsilon }} $、$ {C_{3\varepsilon }} $为经验常数；$ {\sigma _k} $、$ {\sigma _\varepsilon } $分别为与湍动能$k$和耗散率相对应的Prandtl数；${S_k}$和${S_\varepsilon }$为用户自定义的源项。

1.3 VOF法

VOF法^[20]一般用在两相或多相流动的流体问题中，通过计算每个网格单元中各相流体所占的体积分数来捕捉流体的交界面，对于体积分数$ \alpha $的传输方程为：

$ \frac{{\partial \alpha }}{{\partial t}} + \vec u \cdot \nabla \alpha = 0 \alpha =\left\{\begin{array}{l}0\text{，}{\rm air}，\\ 1\text{，}{\rm water}。\end{array}\right. $

(4)

式中：$ \vec u $为速度矢量；$ t $为时间。

对于单个网格单元，若$ \alpha $=0，则该网格单元充满空气；若0＜$ \alpha $＜1，则该网格处于水和空气的交界面上；若$ \alpha $=1，则该网格单元充满水。

在VOF方法中，对于任意单元的流体密度和粘性系数是通过体积分数来计算的。流场中任意位置网格单元中的流体密度$\rho $和粘性系数$\mu $可表示为：

$ \left\{ \begin{gathered} \rho = \alpha {\rho _0} + (1 - \alpha ){\rho _1}，\\ \mu = \alpha {\mu _0} + (1 - \alpha ){\mu _1}。\\ \end{gathered} \right. $

(5)

式中：$ {\rho _0} $、$ {\rho _1} $分别为网格单元中2种流体的密度；$ {\mu _0} $、$ {\mu _1} $分别为网格单元中2种流体的粘性系数。

1.4 数值模拟方法有效性验证

通过模拟文献[21]的浮体淹没实验来验证本文的数值模拟方法的有效性。浮体破口位于底部中央，浮体上方设有6个相同大小的通风口。实验在一个长方形水池中进行，水池的长和宽均为500 mm，水深为400 mm。浮体用细绳通过定滑轮水平悬挂在水面上方，浮体底部与水面相接触，随之剪断细绳，浮体下落，通过高速摄影机捕捉浮体的淹没过程。浮体模型参数如表1所示。建立好的破损浮体模型如图1所示。

通过数值模拟监测破损浮体的进水过程和垂向速度，将所得破损浮体的进水过程图和垂向速度曲线与文献[21]中模型实验所得结果进行对比，破损浮体同一时刻进水过程如图2所示，破损浮体垂向速度曲线如图3所示。破损浮体的进水过程和自由液面位置基本吻合，数值模拟所得浮体垂向速度的变化趋势与模型试验的基本一致，证明了本文采用数值模拟方法的有效性。

2 数值计算模型 2.1 建立船体及舱室模型

以一艘113300 DWT双壳油船为例进行研究，通过SolidWorks建模软件建立船体三维模型（主尺度要素见表2）。坐标原点取在艉垂线与基线的交点处，以船长方向为x轴，向艏为正；沿船宽方向为y轴，向左舷为正，型深方向为z轴，向上为正。建模比例为1∶1，模型只建立裸船体，不考虑螺旋桨和上层建筑等附体的影响。船体模型如图4(a)所示。

其次建立舱室模型，如图4(b)所示。选择该船右舷第三货舱段的双壳舱室建立三维模型。其后舱壁距离艉垂线139.03 m，破口中心点距艉垂线149.84 m。该船设计的双壳宽度和双层底高均为2.40 m。舱室内外破口均设置为圆形，考虑外破口通常大于内破口的情况，取内破口半径为2.80 m，外破口半径为4 m，破口中心点距基线10.29 m、距舱室后舱壁10.81 m。考虑到空气压缩性对进水过程的影响，在主甲板上设置2个圆形通风口，通风口直径均为0.50 m。

2.2 计算域的建立与网格划分

将上述船体和舱室模型导入到STAR-CCM+软件，并构建破损船舶水动力计算域，如图5所示，计算域尺寸为1420 m×740 m×800 m，其中水深为500 m，空气部分为300 m。计算域划分为背景区域、重叠区域和波浪区域。计算域边界条件设置如下：计算域前侧、顶部和底部边界为速度入口，后侧边界为压力出口，左右两侧边界为对称平面。舱室破口和甲板通风口设置为交界面，船体表面设置为光滑壁面。

生成计算域网格时，选择切割体网格生成器来生成网格，选择自动表面修复和表面重构来保证生成网格质量。背景区域的网格基础尺寸为10 m，为更好地捕捉气液流动和船体运动，将重叠网格区域和波浪区域的网格进行细化，细化后的尺寸为基础尺寸的25%。由于舱室内部、舱室破口以及甲板通风口处的流体运动更为剧烈，将舱室内部的网格基础尺寸设置为1 m，舱室破口、甲板通风口的尺寸为舱室内部网格尺寸的6.25%。计算域网格截面和船体表面网格细节如图6所示。

2.3 网格无关性分析

选择合适的网格密度是高效运用CFD进行数值模拟的关键，在此对在静水中113300 DWT双壳油船单壳破损工况进行网格无关性验证，分别选取了3套不同疏密程度的网格对破损双壳油船在静水中的破损进水进行数值模拟，网格数量分别为147.28万、321.13万和687.42万，模拟进水时长为60 s，分析网格尺寸对计算精度的影响。3套网格尺寸差异见表3。

由图7可以看出，对于破损舱室的进水量变化，网格1与网格2、3的数值结果略有差别，网格2和网格3的数值模拟结果较为接近，表明321.12万网格已能满足计算精度的要求。可见网格尺寸过大对计算结果略有影响，但网格尺寸过小会使计算时间呈指数增长，对计算资源消耗过大，为兼顾计算精度和求解效率，选择网格2进行计算。

2.4 仿真工况

运用STAR-CCM+对右舷第三货舱段的单舱单壳破损（以下简称“单壳破损”）工况和单舱双壳破损（以下简称“双壳破损”）工况进行了数值模拟（单、双壳破损工况均为第三类破损舱室进水）。模拟海况分别为静水、尾随浪和横浪，其中横浪的方向取正对舱室破口。共模拟了6种工况（见表4），经过试算，单工况模拟时长设置为200 s。模拟时，假定所有破损舱室均为空舱，监测破损舱室的进水过程及船体的六自由度运动。

模拟使用一台双路Intel® Xeon® E5—2696 v3处理器（36核心72线程），内存为64 G的图形工作站，模拟时调用40线程，内存占用约27 G，模拟6种工况总消耗物理时间约为613 h（约为25天）。

3 模拟结果分析 3.1 不同海况下单、双壳破损进水过程对比分析

图8(a)为单壳破损时不同海况下同一时刻的进水过程对比图，以横浪为例，破损船体在横浪的作用下，外界水体在t=0.02 s开始涌入压载水舱内，冲击到船体内壳后向四周扩散。在t=7.78 s时，外界水体不断流入压载水舱底部。在t=23.39 s时，压载水舱内自由液面进一步升高。在t=187.97 s时压载水舱内的进水量达到最大值2786.06 t。

图8(b)为双壳破损时不同海况下同一时刻的进水过程对比图，以横浪为例，破损船体在横浪的作用下，外界水体在t=0.02 s时开始涌入破损舱室；在t=13.89 s时，外界水体在水压的作用下形成水柱，部分跨过双层舷侧涌入货油舱内部，并对货油舱底部形成冲击，部分涌入水体在冲击船体内壳后涌入压载水舱内；在t=25.73 s时，货油舱内自由液面到达货油舱破口底部，此时压载水舱内的进水速度开始加快，货油舱内水体在外部水体的冲击下形成涡流；在t=193.46 s时，货油舱及压载水舱的进水量达到最大值9782.31 t。船体随波浪起伏，仍存在周期性的内外液交换。

对比不同海况下单、双壳破损进水量和进水量速率历时曲线（见图9）和不同海况下单、双壳破损进水数据（见表5）可知：

1） 单壳破损时，在静水、尾随浪和横浪海况下，破损压载水舱分别在t=197.43、192.45、187.97 s时的进水量达到最大值，最大值分别为2524.41 、2621.93 、2786.06 t。可见，横浪下的进水量最大，进水速度也最快，尾随浪次之，静水下最小。尾随浪和横浪下的较静水下的最大进水量分别增加了3.86%和10.36%。

2） 单壳破损时，不同海况对进水量波动和进水量速率波动幅度的影响也较大，其中，横浪下的进水量最大波动幅度（437.23 t）是尾随浪下（301.34 t）的1.45倍，进水量速率最大波动幅度（333.70 t/s）是尾随浪下（239.46 t/s）的1.39倍。

3） 双壳破损时，在静水、尾随浪和横浪海况下，分别在t=185.73、174.41、193.46 s时刻的进水量达到最大值，最大值分别为9856.34、9705.12、9782.31 t。尾随浪和横浪下的较静水下的最大进水量分别减少了1.53%和0.75%。由此可见，不同海况对该双壳破损工况的最大进水量值影响不大。尾随浪和横浪对最大进水量不只有增加的作用（如单壳破损时），也有减少的作用（如双壳破损时），具体是增加还是减少的作用，可能因海况和不同破损工况而异。

4） 双壳破损时，不同海况对该工况的进水量和进水量速率波动幅度影响也较大，横浪下进水量最大波动幅度（868.34 t）是尾随浪下（524.85 t）的1.65倍，进水量速率最大波动幅度（330.15 t/s）是尾随浪下（223.37 t/s）的1.48倍。可见，横浪对该单、双壳破损工况的进水量波动幅度和进水量速率波动幅度的影响都比尾随浪大。结合单、双壳破损数据可知，单、双壳破损工况下的进水量速率最大波动幅度值相差较小。

3.2 不同海况下单、双壳破损船体六自由度变化对比分析

图10为不同海况下单、双壳破损船体六自由度变化历时曲线。表6为（计算时域200 s内）不同海况下单、双壳破损船体六自由度变化数据。

由图10和表6可知：

1） 无论是单壳破损还是双壳破损，尾随浪对最大横摇角的影响都不大，且与静水下的相差较小，横浪对最大横摇角的影响较大，其中，横浪对单壳破损时最大横摇角的影响（静水的1.66倍）远大于对双壳破损时（静水的1.10倍）的影响。双壳破损时，由于进水量更大，无论哪种海况下都比单壳破损时的最大横摇角大。横浪对单、双壳破损时的最大横摇波动幅度的影响最大，静水下次之，尾随浪下最小。

2） 对比单、双壳破损时产生的艏摇趋势可以看出，单壳破损时尾随浪下的艏摇角增长最快，在横浪下，艏摇角增大到9.99 °时开始出现减小的趋势。在单壳破损时，无论哪种海况下的艏摇角增长速度都大于双壳破损时的，且两者艏摇方向相反，这是由于船体破损位置更靠近船艏，在单壳破损时，外界水体在水压力的作用下跨过船体外壳破口，对船体内壳产生冲击，导致船艏逐渐向左舷偏转；在双壳破损时，由于破损进水后的船体整体向右舷倾斜，在尾随浪和横浪的推动下，导致船艏逐渐向右舷偏转。

3） 单壳破损时，横浪下产生的最大纵摇角最大（1.42°），其余5种工况的最大纵摇角值都不超过1.40°。

4） 单、双壳破损时，尾随浪下的纵荡值增长最快，呈现向船艏方向运动的趋势，且单壳破损时比双壳破损时的纵荡值增长更快，这是由于单壳破损工况比双壳破损工况进水量更少导致的。横浪下和双壳破损静水下的纵荡呈现向船艉方向运动的趋势。

5） 单、双壳破损时，横浪下的横荡值增长最快，呈现向左舷方向运动的趋势，且单壳破损时比双壳破损时的横荡值增长更快。单、双壳破损在静水和尾随浪下的横荡方向相反，单壳破损时，呈现向左舷偏移的趋势，双壳破损时，呈现向右舷方向运动的趋势。

6） 横浪对单、双壳破损时的垂荡影响也较大，分别可达1.61 m和1.41 m。其余4种工况中双壳破损时比单壳破损时的垂荡值更大。最大垂荡波动幅度也是在横浪下最大，且双壳破损时比单壳破损更大。

4 结　语

本文对不同海况下一艘113300 DWT双壳油船的右舷第三货舱段单舱单壳破损和单舱双壳破损2种工况进行了数值模拟，结果表明：

1）单壳破损时，尾随浪和横浪下的最大进水量分别比静水下增加了3.86%和10.36%。在双壳破损时，尾随浪和横浪下的最大进水量与静水下的相差较小（不超过1.53%）。横浪对单壳破损和双壳破损工况的进水量波动幅度和进水量速率波动幅度的影响最大。

2）无论是单、双壳破损，横浪对最大横摇角值的影响较大，尾随浪对最大横摇角值的影响不大，与静水下的相差较小。横浪对最大横摇波动幅度的影响最大，静水下次之，尾随浪下最小。在单壳破损时，无论哪种海况下的艏摇角增长速度都大于双壳破损时的，且单、双壳破损时的艏摇方向相反。

3）单、双壳破损时的船体纵荡受尾随浪影响较大；横荡受横浪影响较大，单、双壳破损在静水和尾随浪下的横荡方向相反。单壳破损时比双壳破损时的纵荡和横荡增长更快。船体垂荡受横浪影响较大，单、双壳破损时的最大垂荡值达到了1.61 m和1.41 m。