0.   引　言

远洋大型集装箱船具有大开口结构，总体结构刚度较柔，在航行于深远海航线时，极有可能遭遇超过其自身极限承载能力的大波浪，导致船体结构崩溃断裂，从而带来巨大的生命和财产损失。在1969—1994年，全世界船舶遭遇极端海况而导致的海难事故达到22起，超过542人因而丧生^[1]。近些年，船舶遭遇波浪发生的崩溃断裂事故仍时有发生，如2013年日本“舒适”号集装箱船遭遇波浪发生的船体断裂事故^[2]和2021年乌克兰Arvin号货船在土耳其海域遭遇波浪发生的断裂事故。因此，为了避免船舶结构崩溃事故的发生，需要对集装箱船在物面大变形、甲板上浪、砰击等非线性波浪载荷作用下的结构动态响应及崩溃模式予以准确评估。并且，对于此强非线性流固耦合问题，国内外船舶研究学者已将其作为流固耦合的研究热点。

最早，学界将基于势流理论的波浪载荷求解器与结构非线性求解器结合，采用“两步法”的求解思路，只将波浪载荷单向映射到流结构模型中，以此计算船舶结构崩溃响应。Masaoka^[3]提出了一种结合非线性有限元与Ursell–Tasai法的数值方法，该方法考虑作用在船体切片上的非线性流体力，施加在船体非线性梁上计算船体的非线性弯矩。Iijima^[4]提出了水弹塑性的概念（称之为的水弹塑性方法），将切片理论与非线性梁的理论结合，对设计的水弹塑性试验进行数值仿真，计算船舯转角大变形。Liu^[5]采用非线性切片理论计算了集装箱船在畸形波作用下船首的最大砰击载荷，并将砰击载荷加入到船体梁上得到其非线性结构响应。刘维勤等^[6-7]将三维面元法求解器与非线性有限元求解器结合，实现了对任意三维船舶结构崩溃响应的模拟，与水弹屈曲试验结果对比结果良好，此外，还将三维势流理论的面元法与非线性有限元进行联合仿真，模拟了某集装箱船在波浪下的结构崩溃响应。

近年来，有学者将黏流理论CFD求解器与有限元求解器结合，构建双向流固耦合路径，开展船舶结构响应计算。Takami^[8]采用CFD–FEA方法计算船体结构在波浪中的总体和局部弹性响应，并将计算结果与非线性切片法、面元法和模型试验结果进行了对比。陈月 ^[9]基于STAR-CCM+和ABAQUS分析了砰击载荷作用下的船体响应。Jiao等^[10]基于STAR-CCM+和ABAQUS计算了S-175集装箱船的甲板上浪和砰击现象。Liu等^[11]基于STAR-CCM+和ABAQUS采用单向耦合和双向耦合方法分别分析了模型尺度的近岸结构水弹塑性，发现相较于单向耦合结果，双向耦合结果与试验结果对比更为吻合和准确。目前，CFD−线性有限元双向耦合方法多应用于船舶结构水弹性问题，尚未用于实船结构非线性问题的研究，故有望启用非线性有限元求解器与CFD求解器，构建双向流固耦合路径，实现对波浪下的实船结构崩溃的模拟。

基于势流理论的“两步法”仅能施行单向耦合，不仅无法充分考虑流体的黏性及自由面破碎、上浪、砰击等非线性波浪载荷，且也无法考虑动态大变形对波浪载荷的影响，用于评估船舶在波浪下的结构崩溃问题时仍存在不足。因此，本文将基于CFD求解器与显式动力学非线性有限元求解器，采用双向耦合方法对4600 TEU集装箱船在迎浪规则波下的弹性、弹塑性和后极限强度工况进行数值仿真计算，以分析其弹塑性响应及结构动态崩溃模式。

1.   基本理论和数值方法

本文研究基于CFD理论与显式动力学非线性有限元理论，采用双向耦合方法求解集装箱船在波浪下的结构崩溃响应。因此，以下将依次介绍CFD、显式动态非线性有限元和双向流固耦合相关的数值理论和方法路径。

1.1   计算流体力学方法

CFD求解器基于有限体积法求解N–S方程，其将流场划分为大小不一的控制体，控制体中的流体流动通过连续性方程和动量方程对其予以控制。同时，考虑到计算精度和成本，采用将速度与压力分解为瞬时脉动值和平均值的雷诺时间平均方程（RANS）对N–S方程进行求解，如式（1）所示。

式中：$\bar{{{\boldsymbol{u}}}_{{i}}}和\bar{{{\boldsymbol{u}}}_{{j}}}=(\bar{{{u}}},\bar{v},\bar{w})$为笛卡尔直角坐标系中x，y，z方向上的雷诺平均速度分量；${{\boldsymbol{x}}}_{i}和{{\boldsymbol{x}}}_j=(x,y,{\textit{z}})$为笛卡尔直角坐标值；t为时间；$\rho$为流体密度；$\bar p$为时间的平均压力；$\mu$为流体的动力黏性系数；$\rho {{\boldsymbol{u}}'_i}{{\boldsymbol{u}}'_j}$为雷诺应力项，是一个二阶张量。

随着表现湍流脉动的雷诺应力项的引入，采用湍流模型来封闭方程。本文采用${\rm{SST}}{\text{ }}k - \omega$湍流模型与RANS方程相结合来求解整个流体域，并利用流体体积法（VOF）对自由液面进行捕捉。${\rm{SST}}{\text{ }}k - \omega$模型的湍流动能和耗散率的输运方程描述如下：

式中：k为湍动能；$\tilde {{G_k}}$为湍动能的生成项；${\mu _{\rm{t}}}$为湍流黏度；$\omega$为湍流耗散率；${G_\omega }$为比耗散率生成项；${D_\omega }$为$\omega$的交叉扩散项；${Y_k}$为k的湍动能耗散项；${Y_\omega }$为$\omega$的湍动能耗散项；${S_k}$和${S_\omega }$为自定义源项；$\sigma _k$是k的普朗特数；$\sigma _\omega$是$\omega$的普朗特数。

1.2   显式动力学非线性有限元方法

结构系统的通用运动学方程如下：

式中：M，C和K 分别是物体的质量矩阵、阻尼矩阵和刚度矩阵； P为外载荷；u，$\dot {\boldsymbol{u}}$和$\ddot {\boldsymbol{u}}$分别是某个时刻的物体位移、速度和加速度。

本文研究采用显式方法求解船体结构的动态非线性响应。显式动力学非线性有限元在增量步开始时，程序求解动力学平衡方程，表示为用节点质量矩阵M乘以节点加速度等于节点所施加的外力P与I单元内力间的差值，如式（5）所示。

在当前增量步开始时的t时刻，计算加速度：

随后，对加速度在时间上采用中心差分法进行积分。在计算速度的变化时，假定加速度为常数，速度和位移只与前一个时刻的加速度和位移有关，从而不断地在时间上以“显式”的方式向前推进速度和位移，最终完成运动学方程的求解。

1.3   双向流固耦合方法

结构在波浪下的崩溃过程中会出现极强的结构非线性和波浪载荷非线性特性，属于强非线性的流固耦合问题，此类问题在一个耦合迭代时间步内需要多次交换物理量才能获得准确的计算结果。本文通过CFD平台STAR-CCM+和大型非线性有限元平台ABAQUS之间的协同仿真功能，在流固耦合交界面交换压力和节点位移数据以实现双向隐式耦合仿真。从图1所示耦合过程可以发现，在一个耦合迭代时间步内，STAR-CCM+和ABAQUS多次交换波浪载荷F与结构位移及变形U，每次交换都会应用重叠网格技术和网格“变形”功能对水动力模型的网格予以更新，并获得新的船体空间位置X和波浪载荷F。经过多次耦合迭代交换，波浪载荷F、结构位移和变形U、船体表面压力P和船体空间位置X不断趋于收敛，最终在一个耦合迭代时间步结束时求解出稳定的结构响应和流场分布。

图  1  双向流固耦合迭代时间步内的数据交换

Figure  1.  Data exchange of fluid-solid bi-directional coupling in iteration time step

下载: 全尺寸图片 幻灯片

由于在流固耦合过程中需要在流固耦合边界上进行数据传递，因此流固耦合交界面需要满足如下基本守恒定律：

式中：${\tau _{\rm{f}}}$表示流体域的应力；${\tau _{\rm{s}}}$表示结构域的应力；n表示方向矢量；${d_{\rm{f}}}$表示流体域的位移；${d_{\rm{s}}}$表示结构域的位移。

同时，由于在流固耦合边界上的数据传递过程中有限元模型和CFD模型的流固耦合交界面存在许多不同尺寸大小的网格，因此通过径向基函数完成网格插值问题，其公式如下：

式中：$\sigma (x)$为关于x的函数；N为插值点个数；${\alpha _i}$为待求的加权系数；$\varphi$为径向基函数；模$\left\| {x - \left. {{x_i}} \right\|} \right.$为插值点x到${x_i}$的欧氏距离；$p(x)$为低价多项式的插值函数。

2.   4600 TEU集装箱船流固耦合数值模型

建立CFD水动力模型和有限元模型，基于CFD−非线性有限元双向流固耦合方法对一艘远洋航行的4600 TEU集装箱船在波浪下的结构崩溃响应进行计算。其中，4600 TEU集装箱船结构形式采用混合骨架式、双层底和双舷侧设计，全船主要布置有3个大开口货舱，分布于驾驶室和艏楼之间。表1为4600 TEU集装箱船主尺度参数，图2为其结构几何模型。

表  1  4600 TEU集装箱船主尺度

Table  1.  Main dimensions of the 4600 TEU container ship

参数	数值		参数	数值
总长/m	260		型深/m	19.6
垂线间长/m	247		设计吃水/m	11.0
型宽/m	37.3		型排水体积/ m3	64647.7

下载: 导出CSV 

| 显示表格

图  2  4600 TEU集装箱船几何模型

Figure  2.  Geometric model of the 4600 TEU container ship

下载: 全尺寸图片 幻灯片

2.1   水动力模型

建立4600 TEU集装箱船的CFD水动力模型，如图3所示。采用压力耦合方程组的半隐式方法（Simple算法）求解N–S方程，其中，对流项采用迎风差分格式离散，湍流模型采用k–ε湍流模型求解，空气和水之间的自由液面采用VOF方法捕捉。最后，将计算得到的波浪压力结果作为载荷输入到结构计算中。

图  3  4600 TEU集装箱船CFD水动力模型

Figure  3.  CFD hydrodynamic model of the 4600 TEU container ship

下载: 全尺寸图片 幻灯片

计算域包含数值波浪背景域和重叠区域，考虑节约计算成本，应用对称边界条件。首先，在STAR-CCM+上建立半个波浪水池背景域（长1750 m、宽500 m和高750 m），坐标范围相对于船长L约为−3.85L<x<2.88L，0<y<1.92L和−1.88L<z<L。然后，建立用于模拟刚体运动以及与背景区域进行流体数据传输的半个重叠区域（长350 m、宽50 m和高75 m），坐标范围相对于船长约为−0.19L<x<1.15L，0<y<0.19L和−0.096L<z<0.19L。计算域的坐标系原点位于船体0号肋位基线处，x轴正向沿船长方向指向船首，y轴正向沿船宽方向指向左舷，z轴正向沿型深方向。最终，4600 TEU集装箱船的计算域和边界示意图如图3（a）所示。其中，重叠网格交界面处的边界类型为重叠网格，船体模型边界类型为壁面，对称剖面边界类型为对称平面，入口、右侧和底部边界类型为速度入口，顶部和出口边界类型为压力出口，压力出口边界处还设有2倍船长的消波区域，采用消波阻尼方法消除壁面反射波和船体运动兴波对流场的影响。

本文在研究时采用了切割体网格单元生成器和棱柱层网格生成器对计算域进行网格划分，其中船体表面尺寸大小为0.5 m，与结构有限元模型网格尺寸大小相近，以避免映射过程中在流固耦合交界面处发生网格插值错误。壁面附近棱柱层为8层（y+ > 300），船体重叠网格区域最大尺寸为1.5 m，背景域最大尺寸为7.5 m，并在重叠区域与背景域之间划分有一套最大尺寸为4 m的过渡网格，以此来消除粗细网格差异产生的通量误差。最后，根据ITTC^[12]规则波自由液面上波长不少于40个网格加密和波高上不少于20个网格加密的建议，对自由液面波长方向采用波长的1/80进行网格划分，垂直方向采用波高的1/20进行网格划分，以保证波浪的稳定性以及液面破碎、上浪和砰击等自由面非线性现象的捕捉，最终的划分网格如图3（b）所示，总计527万个单元，其中重叠区域为202万个单元。

2.2   结构非线性有限元模型

结合前处理软件HYPERMESH和大型结构非线性计算平台ABAQUS建立4600 TEU集装箱船结构一半有限元模型，如图4所示。为确保模型在计算过程中不会报告由于刚度矩阵的奇异性而产生的错误，在船尾选择一排节点约束船体沿x轴运动来防止船体在波浪作用下发生飘移，同时在船体中纵剖面亦使用对称边界条件，边界条件如表2所示。启用显式动力学非线性有限元计算船舶结构的时域响应，划分的有限元网格尺寸为0.5 m，以确保水动力模型船体表面网格的匹配和对船体结构屈曲行为的模拟，最终的有限元模型共包含有216 485个单元。

图  4  4600 TEU集装箱船有限元模型

Figure  4.  Finite element model of the 4600 TEU container ship

下载: 全尺寸图片 幻灯片

表  2  边界条件

Table  2.  Boundary conditions

位置	Ux	Uy	Uz	URx	URy	URz
船尾	约束	释放	释放	释放	释放	释放
中纵剖面	释放	约束	释放	约束	释放	约束

下载: 导出CSV 

| 显示表格

本文主要采用基于黏流求解器的双向耦合方法对船舶的崩溃问题进行方法研究，聚焦于瞬时湿表面非线性下的结构崩溃问题，不过大增加波高。并且，为了数值模拟计算的效率，在赋予板材厚度时，在保证船体结构刚度分布与实际船体结构一致条件下对其进行了一定的强度削弱处理。本文船体结构的材料为Q235钢，本构模型采用理想弹塑性材料，具体材料参数见表3。

表  3  船体结构材料参数

Table  3.  Material parameters of hull structure

密度/（t·m−3）	弹性模量/MPa	泊松比	屈服强度/MPa	许用应力/MPa
7850	2.06×105	0.33	235	211

下载: 导出CSV 

| 显示表格

2.3   耦合计算设置

在STAR-CCM+和ABAQUS中设定外船体表面为流固耦合交界面。流场的压力数据和结构响应的位移场数据通过该流固耦合交界面进行实时交换，并完成双向流固耦合的求解。选择Surface to Surface的方式进行数据映射。压力场中包含流场动压力和壁面剪切力两个部分，通过流固耦合交界面从流场映射至结构求解器的船体表面；水压力作用下的位移场在ABAQUS中进行求解，并通过流固耦合交界面传递回流场，流场更新后完成一次完整的双向耦合计算。选用隐式耦合求解算法，设定一个时间步内交换5次数据，且每次交换经多次内部迭代，流体计算器在一个计算时间步内进行10次内部迭代。选用恒定的耦合时间步0.005 s。协同仿真采用20核40线程的Intel Xeon E5-2660服务器并行计算，一个工况约需计算15 天。

2.4   计算工况设置

制定计算工况用于研究4600 TEU集装箱船在不同波浪参数下的结构崩溃响应特征。计算工况如表4所示，所有工况采用的波浪皆为一阶线性波，波长均取1倍船长247 m，波高从2.5 m开始每间隔2.5 m等间距增长到15 m，保证船体结构在波浪载荷作用下发生弹性变形、弹塑性变形、极限强度变形和后极限强度变形。本文在船舯崩溃区域自甲板至船底处选取4个典型节点，分别位于船舯主甲板、舷侧外板、船舭外板和船底板，并依次命名为节点1、节点2、节点3和节点4，如图5所示，用于输出时域弹塑性应力曲线，判断其是否崩溃。

表  4  全船崩溃响应数值仿真工况

Table  4.  Numerical simulation conditions for the collapse response of whole ship hull

工况	波长/船长	波高/m		工况	波长/船长	波高/m
H_1	1	2.5		H_4	1	10.0
H_2	1	5.0		H_5	1	12.5
H_3	1	7.5		H_6	1	15.0

下载: 导出CSV 

| 显示表格

图  5  典型应力节点的位置

Figure  5.  Positions of typical stress nodes

下载: 全尺寸图片 幻灯片

3.   结果分析

3.1   CFD网格无关性验证

首先，针对工况H_1的计算进行CFD网格无关性验证。设置了a（粗）、b（中）、c（细）3套网格，基础尺寸分别为12.5，7.5，5 m。网格数量分别为284万，527万，998 万。生成的3套网格的升沉、纵摇及节点2的应力结果如图6所示。可以看出a（粗）网格计算幅值较小，比细网格的升沉幅值约小9.1%，纵摇幅值约小7.6%，应力幅值约小7.9%；中网格和细网格的计算幅值非常接近，升沉幅值偏差约为3.3%，纵摇偏差约为2.1%，应力偏差约为2.5%。考虑到数值计算的效率问题，本文选用b（中）网格进行计算。

图  6  网格无关性验证

Figure  6.  Grid independence verification

下载: 全尺寸图片 幻灯片

3.2   典型节点时域应力曲线

采用本文提出的CFD−非线性有限元双向流固耦合方法对4600 TEU集装箱船进行双向流固耦合计算，得到4600 TEU集装箱船在如表4所示的6个工况下的弹性响应及弹塑性响应。最终典型节点沿船长方向的纵向应力结果如图7所示，其中正值代表受压，负值代表受拉。从图中可以看出，整体应力曲线未出现明显的水弹性振动响应，一方面是因为波激振动现象通常在船体一阶总振动频率和波浪遭遇频率成倍数时才会产生，即便是集装箱船的波激振动问题，也不是任何频率下都会产生明显的波激振动现象；另一方面，本文计算时船舶静置在波浪下未给予其航速，其砰击力较小，未能产生明显颤振现象，并且本文关注的重点是弹塑性响应及崩溃问题，计算选取的工况有意避开了水弹性振动响应工况。

图  7  各典型节点在各工况下的时域弹塑性应力曲线

Figure  7.  Time-domain elastoplastic stress curves of typical nodes under various working conditions

下载: 全尺寸图片 幻灯片

图7（a）和图7（b）分别展示了工况H_1（波高2.5 m）和H_2（波高5.0 m）的应力曲线，从图中可以看出，双向耦合过程中波浪载荷较小，典型节点的应力曲线皆呈弹性响应；图7（c）和图7（d）分别展示了工况H_3（波高7.5 m）和H_4（波高10.0 m）的应力曲线，从图中可以看出，工况H_3的节点1在14.96 s首次进入塑性，节点2在26.76 s首次进入塑性，节点3和节点4未进入塑性，工况H_4的节点1在14.01 s首次进入塑性，节点2在28.26 s首次进入塑性，节点3和节点4未进入塑性，这表明随着波高的增大，船舶结构进入塑性的速率增快，且随着局部结构的塑性产生，其刚度开始下降，抵抗变形的能力开始下降，在对应时刻的中拱中垂变形也进一步增大；图7（e）和图7（f）分别展示了工况H_5（波高12.5 m）和H_6（波高15.0 m）的应力曲线，从图中可以看出，随着波高的增大，船舶结构进入塑性的速率增快，大部分节点均已进入了塑性，且已进入结构崩溃状态，应力曲线开始丧失一定的周期性，呈现不规则形状，究其原因，一方面是因为其崩溃过程中承载能力不断下降甚至丧失，另一方面是结构崩溃过程中产生的复杂面外大变形对流场的波浪载荷进一步造成了影响。

3.3   典型工况流固耦合崩溃模式分析

在双向流固耦合过程中，船体结构响应、船体浮态和船体水压力在任一时刻的对应关系皆可通过水动力模型及结构模型计算云图得到直观体现。为此本文选取工况H_5作为流固耦合崩溃模式分析的典型工况，展示其在中垂时刻下的流固耦合崩溃模式，如图8～图10所示。其中，图8为4600 TEU集装箱船在典型工况H_5的波浪自由面分布图。从图中可以看出，波浪波谷接近船舯剖面时，水动力模型通过流固耦合交界面完成了有限元模型产生的纵摇、升沉及中垂变形的网格映射，产生了明显的崩溃变形及纵倾现象，同时艏部出现明显的上浪现象，艉部出现砰击现象；图9为集装箱船的湿表面压力分布图，动态波浪使船舶产生了很大的压力分布，并在双向耦合过程中映射至了结构有限元模型中，两者压力分布完全一致，波浪载荷映射良好；图10为集装箱船结构应力与变形云图，可见，集装箱船的上部至中部结构此时已广泛进入塑性，主甲板、舷侧板、甲板纵骨和舷侧纵骨等构件在波浪载荷作用下发生了明显的屈服及屈曲失稳，主甲板和舷侧板等板材出现明显的褶皱，甲板纵骨和舷侧纵骨等骨材发生严重的侧向失稳，伴随塑性变形区域由舷侧外板延伸至船底板，表明工况H_5已使4600 TEU集装箱船广泛进入结构崩溃，船体结构彻底丧失了承载能力。

图  8  12.5 m波高工况波浪自由面分布

Figure  8.  Wave free surface distribution (wave height is 12.5 m)

下载: 全尺寸图片 幻灯片

图  9  12.5 m波高工况波浪压力云图

Figure  9.  Wave pressure contours (wave height is 12.5 m)

下载: 全尺寸图片 幻灯片

图  10  12.5 m波高工况船体结构应力与变形云图

Figure  10.  Stress contours of hull structure and deformation (wave height is 12.5 m)

下载: 全尺寸图片 幻灯片

3.4   波浪载荷对整体和局部变形的影响

为直观体现船体整体变形响应及局部变形响应与波浪载荷大小的对应关系，本文提取6个典型工况船舯转角变形结果的最大值和典型节点应力结果的最大值进行对比分析，结果如图11所示。由图11（a）可见，船舯最大转角伴随波高的增大而增大，且当波高未超过10 m时，由于波浪载荷较小，船舯结构处在弹性和弹塑性阶段，最大转角变形较小，而波高超过10 m时，波浪载荷超过船舶承载能力，船舯结构大规模进入塑性，承载能力极度下降，船舯最大转角变形值迅速增大，船体结构发生了明显的崩溃破坏。由图11（b）可见，局部结构随着波高的增大，应力不断增加，并且甲板（节点1）和舷侧外板（节点2）相较于船舭外板（节点3）和船底板（节点4），因远离中和轴其进入塑性更快。

图  11  局部最大变形和整体最大变形随波高的变化

Figure  11.  Variation of local maximum deformation and global maximum deformation with wave height

下载: 全尺寸图片 幻灯片

4.   结　论

本文建立4600 TEU集装箱船的CFD水动力模型和结构非线性有限元模型，选取系列典型波浪工况，采用CFD−非线性有限元双向耦合方法对集装箱船进行结构弹塑性响应仿真，研究典型节点应力曲线和流固耦合崩溃模式，得到如下结论：

1） CFD和非线性有限元方法构建的双向流固耦合方法可以考虑甲板上浪、瞬时湿表面、船体大变形等非线性波浪载荷的影响，实现流体载荷和结构位移及变形的双向流固耦合映射，是一种准确评估船舶在波浪下的结构响应及动态崩溃过程的新方法。

2） 在波浪高度较小时，船体结构只发生线弹性变形或小部分结构塑性变形，结构仍处于安全状态，超过一定波高时船舶由部分结构塑性进入大部分结构的广泛塑性，发生结构崩溃现象，且波高越大，结构进入塑性的速率越快。

3） 4600 TEU集装箱船为大开口船，中和轴远离甲板区域，在波浪下的结构崩溃模式为中垂破坏，表现为甲板和上部船体结构先开始出现屈服现象及屈曲失稳，并在塑性累计效应的影响下逐渐由舷侧外板向船底板蔓延，直至完全失去承载能力，整体出现纵向大变形。

本文仅对4600 TEU集装箱船在0航速下的结构崩溃响应进行了评估，对航行状态下砰击颤振和波激振动对崩溃响应的影响有待接下来进一步研究。