王加夏,周天九,刘 昆,陈 月

(江苏科技大学 船舶与海洋工程学院, 镇江 212003)

船舶在航行过程中，会发生摇荡运动，这时剧烈的波浪会对船体发生猛烈冲击，这一现象被称为砰击[1-2]．波浪砰击过程涉及波浪强非线性、瞬时效应的物理特性，也涉及水弹性效应及流固耦合等复杂的分析技术，一直是海洋工程领域研究的热点课题之一．严重的砰击不仅可能会对冲击区域产生巨大的冲击压力，造成局部结构破坏;而且，频繁的波浪砰击也会诱导船体出现二节点的颤振现象，造成船体垂向弯矩极限值的增加，影响结构的极限强度．

文献[3]中首先对砰击问题进行了理论研究，利用动量定理得出了经典的砰击压力的计算公式．在此基础上，计入自由液面升高现象，将结构入水砰击过程简化为平板入水的模型，得到了冲击理论模型[4]．后续学者在不同方面进行了大量的修正．其中较为突出并且应用较广的分别为广义Wagner模型[5]和修正的Logvinovich模型[6]．文献[7]中对一艘豪华游轮进行了尾部砰击荷载和振动响应的模型试验研究，提出了一种用试验测量尾砰击压力场的方法．文献[8]中在拖曳水池中对某舰船的分段模型进行了尾砰击响应试验研究，指出在尾随浪、零航速、波长船长比λ/L=1.0左右时观察到了严重的尾砰击现象．文献[9]中针对静水和波浪场景下的楔形体入水砰击进行试验研究，文献[10]中针对平台立柱的入水砰击特性进行研究．文献[11-13]中研究了砰击场景下多船体甲板结构的动态响应问题，通过水弹性梁模型来模拟甲板，并通过初始的结构惯性阶段和随后的自由振动阶段模拟砰击场景，结构应变和位移结果表明理论与试验结果吻合较好．文献[14]中针对不同船级社关于滚装船外飘砰击校核公式进行对比分析研究．文献[15]中利用STAR-CCM+(STAR)对不同二维和三维刚性和弹性楔体的入水冲击特性进行了研究，认为STAR可以较为准确预测楔形体的位移和速度的时间历程．

文中使用基于CFD求解器的STAR与有限元软件Abaqus进行交互耦合计算，对船体模型[16]在试验水池的砰击现象进行数值仿真，分析船体发生砰击时的自由液面形态分布状况、砰击压力、结构变形以及应力应变分布状况．

1 波浪砰击的流固耦合数值模型

1.1 流体控制方程

STAR软件基于FVM法将控制方程的积分形式进行时间和空间维度的离散，从而得到可求解的代数方程组．首先计算域被细分为有限数量的邻接控制体，这些控制体可以是任意多面体形状．离散的控制方程在计算过程中还需要用到面积积分、体积积分以及时间和空间导数．

对于粘性的三维流动，假定流动由RANS方程控制，其中湍流效应包括涡流模型和粘性模型．此时就需要求解连续性方程、动量方程和两个湍流特性方程．模型选取K-Epsilon湍流模型来模拟湍流对流体的影响，其中积分形式的质量守恒和动量守恒控制方程可以表示为[15]:

(1)

(2)

式中：ρ为流体密度；v为流体速度矢量；vb为控制体表面的速度；n为垂直于控制体表面的单位向量法线，其中S和V分别为控制体的面积和体积；T为指应力张量(表示速度梯度和涡粘性)；p为压力，I为单位张量．

由于VOF模型假定在控制体中各个不相融流体都具有相同的速度、压力和温度场．因此，在控制体中只需要求解单相流场的动量、质量和能量守恒的基本控制方程即可．描述体积分数αi的输运的守恒方程为:

(3)

式中：sαi为第i相流场的源或汇；Dρi/Dt为相密度ρi的物质导数．

1.2 结构响应方程

对于研究三维舰船在砰击载荷下的动态响应时，需要利用有限单元法构造考虑弹性结构的响应方程，将三维模型进行网格离散，得到结构的有限元模型，是进行有限元计算的基本步骤．三维弹性动力方程的基本方程[17]是在得到船体结构的有限元模型的基础上，进行三维实体动力分析的有限元求解步骤:

(4)

1.3 流固耦合

船体结构波浪砰击问题属于强非线性的流固耦合问题，文中通过STAR与Abaqus之间的协同耦合功能，自动在流固交界面处进行数据的交互传递实现双向的流固耦合计算．其中STAR求解流体方程，Abaqus求解结构方程，属于分区式流固耦合．采用STAR和Abaqus之间通过隐式的耦合方法进行耦合计算，主要是因为隐式耦合方法适用于结构和流体之间的强耦合应用，它允许STAR和Abaqus在每一时间步内交换不止一次数据，不停迭代直到满足某些收敛准则．这种耦合方案更稳定，二阶精度更高．

由于波浪砰击现象属于强耦合现象，因此在数值模型中，选择隐式强耦合算法．文中依单个耦合步骤为例，针对隐式耦合时间步的流程进行介绍，具体如图1．在STAR中设置运动规范为“变形”，在Abaqus设置模型的弹性模量、板厚、初始速度、重力、时间步长和耦合面等参量．具体计算流程如下:首先STAR开始流场载荷计算，然后将压力和剪切力传递给FEM求解器，之后Abaqus根据耦合界面的载荷进行结构响应分析并将得到的位移和变形量传递给STAR，STAR再根据传递的位移量通过“变形”功能实现网格的移动，这样完成一次协同耦合进行计算更新，重复直至达到最大物理时间，tn为计算时间步．

图1 隐式算法的单个耦合时间步的流程Fig.1 Procedure for a single coupling step with implicit solvers

2 模型设置

2.1 模型尺寸

文中建立的波浪砰击载荷下船体结构动态响应的数值仿真模型和文献[16]中的铝制模型船试验作为对比，文献模型的主尺度参数如表1，船体的有限元模型如图2．

表1 主尺度参数Table 1 Principal dimension

图2 铝船的有限元模型Fig.2 FEM model of the aluminum ship

2.2 计算域设置及网格划分情况

坐标系原点在船尾垂线与船底平面的交点处，X轴沿船长方向指向船首，Z轴垂直向上，并满足右手直角坐标系．计算域如图3，其中背景区域(background)坐标为[-8，-4，-4]，[24， 4， 4]，重叠区域(overset)坐标为[-1.5，-1.2，-0.5]，[5， 1.2， 0.8]，网格的总数为1 178 788个．模拟收敛精度为98.779%．计算域背景区域[18]和重叠区域的边界条件设置可参见表2．

图3 计算域Fig.3 Computational domain

表2 边界条件Table 2 Boundary condition

整体计算域在Y=0处剖面网格划分如图4，最大网格尺寸为0.2 m，加密区域网格尺寸为0.1 m，切割体网格单元为各项同性．另外，为捕捉自由液面处波浪形态和存在的波浪翻卷及破碎现象，对船体表面邻近区域的网格进行加密，如图4，水面附近网格Z方向加密，Z方向尺寸0.025 m，重叠区域网格尺寸为0.1 m，船表面尺寸为0.05 m，棱柱层厚度为0.025 m．

图4 液面处网格加密区域Fig.4 Section meshing at the free surface domain

文中在Abaqus中设置试验模型的板厚、弹性模量等材料属性，模型使用5 mm厚度的铝板，并对结构施加初始速度和重 力，导出inp文件．在STAR的“协同仿真”中确定导入位移和速度，导出压力和壁面剪切应力，读取Abaqus中模型的inp文件，实现STAR与Abaqus之间的数据交换．设置区域中重叠区域(overset)的运动规范为“变形”．

文中基于建立的波浪砰击耦合数值模型，研究了规则波迎浪情况下铝制模型船的砰击场景及其船体结构响应．选取船体无航速(工况1)和有航速(工况2)进行数值仿真，分析结构运动及砰击响应特性．选取的两组工况如表3．

表3 仿真选取的两种工况Table 3 Two working conditions selected by simulation

为保证船体波浪砰击耦合分析过程中数值模型的正确性，将数值结果和文献结果在砰击载荷进行对比，具体如图5，t为试验及仿真计算时刻，图中给出了数值和试验关于支杆与船体连接处的垂向载荷的对比，从图中可以看出，仿真结果与文献值整体曲线趋势一致，吻合度较高，仿真载荷幅值稍小于文献值，这可能是由于波浪试验是在水域有限的水池中进行的，壁面对波的反射可能会造成试验结果的偏大，而且工况1的整体误差稍大于工况2的情况，这可能是由于工况1的波高比工况2更大，导致文献中试验结果与仿真之间的差值也越大，另一方面，数值计算误差也对精度有一定的影响．

图5 垂向力的试验值和仿真结果的对比Fig.5 Comparison of experimental and numerical values of heave force

根据图5的载荷对比曲线可知，文中所提出的波浪砰击下船体耦合响应数值仿真方法是有效的．

3 模型船体砰击压力

3.1 砰击时刻自由液面形态和模型压力云图

图6为铝制模型船在工况1时发生砰击的前后时刻情况，当船舶漂浮在波浪上时，由于波浪的周期性起伏运动，船体也随之出现升沉和纵摇现象，图中主要给出了不同时刻模型自由液面形态及船底表面砰击压力云图．图6(a)为t=1.48 s时自由液面形态和模型表面压力云图．从图中可知，该时刻波峰位于模型中部偏船尾位置，模型首部位于波谷处，首部表面完全出水，压力集中于船底偏尾部区域，首部压力最小．由图6(b)可知，t=2.06 s时模型尾部即将进入波谷，船首位于波浪波峰处，船首底部倾斜处发生波浪冲击，模型船首表面与波浪水面完全接触，最大压力集中在首部．该工况时波浪表面未发生大的表面变形甚至波浪破碎现象．

图6 工况1船体各时刻自由液面形态和表面压力云图Fig.6 Free surface elevation and the pressure distribution on the ship model in case 1

图7为铝制模型船在工况2时发生砰击前后时刻的模型自由液面形态及船底表面砰击压力云图．

图7 工况2船体各时刻自由液面形态和表面压力云图Fig.7 Free surface elevation and the pressure distribution on the ship model in case 2

由图7(a)可知，在t=1.25 s时波峰处于模型后方近船尾处，模型首部位于波谷处，首部表面完全出水，压力主要集中在船底偏尾部，而船首压力较小．由图7(b)可知，在t=1.56 s时刻波峰运动到船首处，船尾处离开上一个波峰位置即将进入波谷，模型首部倾斜处发生相对波浪冲击，波浪表面沿模型表面出现射流攀升现象，最大压力出现在模型首部倾斜位置．

对比图6、7可知，在工况2时，由于模型具有向前运动的航速，使得船体与波浪之间存在较大的相对速度，造成模型在较小波高的情况下，船首倾斜处与波浪水面交界处出现了液面的猛烈冲击现象，甚至出现波浪破碎射流分离的强非线性现象，造成船体首部出现较大的砰击载荷．

3.2 测点砰击压力和速度对比曲线

如图8，在船底布置一系列测点以探测船体在波浪中受到的砰击压力及结构变形位移．P2位于船首底部转折处，坐标为(3.4， 0， 0);P4点位于船尾底部转折处，坐标为(0.4， 0， 0)．

图8 测点位置示意Fig.8 Schematic map of measuring point position

图9为工况1时不同测点的压力和速度曲线历程曲线，其中实线为砰击压力曲线，虚线为速度曲线，左侧Y轴为压力值，右侧Y轴为速度值．

图9 工况1测点压力和速度对比曲线Fig.9 Correlation curve of pressure and velocity at measuring point in case 1

从图中可知，压力p和速度v曲线均呈现周期性变化，反映了砰击载荷的频繁出现的特性．并且图中可知P2测点(船首底部)的压力幅值最大为1.64 kPa，P4测点(船尾底部)幅值为0.90 kPa，船体首尾处均相继出现砰击现象．图中压力曲线峰值稍有一定的衰减，可能是由于湍流模型存在粘度以及阻尼造成波高衰减造成的，而且P2测点压力曲线呈现正弦函数形式，与规则波的输入相符，此时砰击压力曲线与速度曲线峰值向对应，表明了砰击载荷与砰击速度的相关性，其余测点的压力幅值与速度幅值存在一定的相位差，且测点的压力曲线均较为规则．

图10为工况2时测点的压力和速度曲线对比，对比图10、9，可以发现两种工况下相应各测点整体曲线特征相似，工况2时船首区域测点压力峰值为3.294 kPa，大于工况1中对应测点峰值，这主要是因为工况2中增加了船舶航速，使得砰击现象较工况1更加明显．工况2中P4测点均出现负压载荷，这是由测点所处位置的流场载荷的往复变化以及结构的变形效应共同作用下的结果．P2测点位于船艏底部，周围流场呈规律性变化，而P4测点位于船艉，且周围流场由于受到波浪破碎的原因，导致压力分布特性变化较大．

图10 工况2测点压力和速度对比曲线Fig.10 Correlation curve of pressure and velocity at measuring point in case 2

3.3 船体结构响应

图11为工况1时不同时刻的模型结构应力(单位:MPa)和应变云图.图12为工况2时不同时刻的模型结构应力和应变云图．

图11 工况1模型应力和应变云图Fig.11 Stress and strain nephogram of case 1

图11为工况1下t=1.48 s时模型应力和应变云图，由图6可知，t=1.48 s时刻船体中后部处于波峰位置，船首处完全出水，船体出现中拱现象，对应图11中船体底部及舷侧出现大的应力和应变集中，而船首处应力、应变较小．此时的对应的最大应力值为2.2 Mpa，最大应变值为3.053×10-5．

图12为工况2时船体受波浪砰击前后时刻应力、应变云图，对应图12、7可知，工况2中t=1.56 s时船首处发生波浪砰击现象，此时船首底部与波浪发生猛烈砰击，出现波浪破碎的情况，因此造成船首应力、应变增大，出现相应的应力、应变集中区，船体尾部区域的应力明显较小．此时的对应的最大应力值为0.683 MPa，最大应变值为7.34×10-6．

图12 工况2模型应力和应变云图Fig.12 Stress and strain nephogram of case 2

4 结论

基于文献中的模型试验进行了仿真分析，通过对比试验测量值和数值仿真结果，检验文中建立的波浪砰击载荷下船体结构动态响应的数值仿真模型的可行性．分析船体的发生砰击时的自由液面形态分布状况、砰击压力、结构变形以及应力应变分布状况．

(1) 所建立的计及结构变形效应的波浪砰击流固耦合仿真方法与文献中试验结果整体趋势吻合较好，只是幅值相对试验稍小，这可能是由于数值模拟与实船试验场景的差异如忽略了水池实验的壁面效应以及数值误差的影响造成的．因此文中建立的数值模型对于模拟船舶在波浪中的运动及砰击现象具有较大的适用性和可行性．

(2) 当船体在波浪中运动时，船首区域由于波浪的冲击会发生砰击现象，且当船体有航速时自由液面出现了波浪表面破碎等强非线性的砰击现象，并分析了在不同航行阶段中船体结构应力和应变的变化特征，观察到船底应力变化，表明船首底部是波浪砰击问题中需要重要关注的区域．