考虑船-冰-水耦合的极地航行舰船载荷响应分析

丛 滨，慕 鹏，宋星晨，李陈峰

(1. 中国人民解放军92941 部队，辽宁 葫芦岛 125000；2. 哈尔滨工程大学，黑龙江 哈尔滨 150001)

0 引 言

与开阔水域所遭受的风浪流载荷作用不同，极地地区冰情复杂，存在许多大小形状不定的浮冰，舰船在极地地区还要承受冰载荷的碰撞与挤压作用，对船体结构造成极大的威胁。因此对船舶在极地海域航行过程中可能遭受的载荷成分进行分析，对船舶在冰区航行操纵性与安全性有重要意义。

目前国内外学者开展了大量的船-冰碰撞数值模拟分析，随着计算机技术的发展，形成了以网格单元建模的有限元方法(FEM)、无网格化建模的离散元方法(DEM)、光滑粒子流体动力学方法(SPH)及近场动力学方法(PD)[1]。蒋昱妍[2]基于粘聚单元法进行了层冰-结构物碰撞的数值仿真，对海冰失效模式和碰撞冰阻力进行了研究，但忽略了海水流固耦合效应的影响。Gagnon[3]基于流体动力学原理，建立了油船右舷舷侧与冰山碰撞的数值模型，讨论了碰撞过程中右舷结构上的负载与压力分布。Jou[4]对破冰船压溃式破冰方式进行了研究，通过DEM 方法对海冰及其裂缝进行建模，预测了不同形状破冰船船首与浮冰接触时可能产生的浮冰破坏模式。Khayyer[5]提出了1 种增强的ISPH 流体模型和新的半隐式ISPH-SPH 耦合方法，并针对典型流固耦合问题进行了模拟，做了严格的验证和性能研究，是SPH 方法的巨大进步。薛彦卓[6]将PD方法用于研究海洋结构物与冰的相互作用问题，如基于近场动力学进行了冰的3 点弯曲试验模拟、船舶与水平层冰相互作用的数值模拟等。上述文献对于海冰-海洋结构物相互作用的研究多集中分析船舶所受冰载荷作用，对浮冰的运动状态进行了约束，部分忽略了水域的作用，未对水阻力及水域作用下浮冰-碰撞的运动效果进行分析。

为了研究船-冰-水三相耦合时的船体载荷响应情况，本文根据船在冰缘区浮冰漂浮特性和船舶航行特性，在Ls-dyna 中通过S-ALE 流固耦合方法进行船-冰-水相互作用机理研究。在此基础上，研究不同尺寸浮冰与舰船相互作用的水阻力与冰载荷响应规律。

1 数值方法原理

1.1 S-ALE 流固耦合方法

船舶与浮冰发生碰撞，由于体积与质量差异，浮冰运动姿态的改变更为明显，因此船-冰-水的流固耦合模拟效果更多体现在浮冰运动动态的真实性与船体受载的准确性方面。本文在研究过程中对流体与结构的耦合问题描述，采用结构化-任意拉格朗日-欧拉(structured arbitrary Lagrangian-Euler algorithm，SALE)方法。

S-ALE 方法是在ALE 方法上发展的更为高阶的流固耦合模拟方法。ALE 方法描述的流固耦合是指发生在流体与固体介质交界处的协调关系和平衡，采用拉格朗日变量对固态结构进行运动边界的描述、采用欧拉变量对流体的运动特性进行描述，可以有效跟踪物质结构边界运动和防止网格产生畸变。S-ALE 方法在原有算法和逻辑的基础上，进行了规律性简化和自适应网格生成，在计算时间上比ALE 方法减少20%～40%。用光滑delta 近似函数通过添加“源项”、分布节点力和插值速度来表示流场和结构物的交互作用，实现结构与流场的作用及反作用模拟，优化了ALE 算法由于流体-结构之间接触定义产生的流体渗漏问题和压力分布的连续性问题[7-8]。一般不可压N-S 方程和力源项的分布函数为：

式中：ρ为流体密度；u为流场速度；p为流场压力；t为时间；µ为动态粘度系数。F(x,t)为流场受到边界的所有体积力，其中x为浸入边界位移；f(s,t)为浸入边界产生的单位力；X(s,t)为Dirac delta 函数；δ为光滑化的Dirac 函数。

1.2 罚函数接触算法

船冰碰撞过程中存在大量接触行为，如浮冰与船舶初始行为阶段的碰撞，首柱舷侧两边受到冰体的挤压行为，因此需要通过接触理论定义求解。船体与浮冰均为弹塑性结构，接触设置选择为面对面接触，接触类型选择为侵蚀接触，接触算法为罚函数法[9]。罚函数法从判断物体接触是否穿透的角度进行分析，将要接触的2 个面设置为主面与从面，在接触进行的各时间段内，自动检测主从面是否发生穿透。未穿透时，接触继续进行；当发生穿透时，在主面与从面之间引入垂直于主面法向的接触力，作用方式类似于弹簧，如图1 所示[10]。这一接触力F可由下式得出：

式中：k为主面刚度；δ为穿透深度值。

主面刚度k与接触面刚度的惩罚因子相关，在Lsdyna 中，惩罚因子的设置关乎数值模拟精度[11]。k的计算公式如下：

式中：pf为惩罚因子；V为主体积；K为接触单元体积弹性模量；A为主从面接触段面积。

图1 罚函数法作用原理Fig. 1 Principle of penalty function method

2 船-冰-水耦合模型

以某船首部为研究对象，如图2 所示。首部钢材为907A 高强钢，名义屈服强度为390 MPa，材料模型选择各向同性塑性模型，材料参数如表1 所示。

图2 船舶首部模型Fig. 2 Model of ship bow

表1 船体材料Tab. 1 Ship material parameters

船-冰-水耦合数值模型由船体、浮冰、海水、空气四部分构成。在Ls-dyna 中，S-ALE 方法对于流体域PART 的设置通过网格PART 和材料PART 的组合进行定义。同时，水及空气流体材料的物理特性通过材料与状态方程结合来描述的，表达式如下：

式中：C0～C6为多项式拟合参数，ς为体积参数，E0为初始内能。

浮冰采用SOLID164 实体单元进行建模，模型尺寸为10m×10m×1m，网格尺寸为0.2m×0.2m×0.2m。冰材料模型选择遵循Mohr-Coulomb 屈服准则[12]的各向同性弹塑性本构模型。根据国内外学者相关试验研究和海冰数据统计，本文海冰材料参数如表2 所示。

表2 空气、水状态方程参数Tab. 2 The state equation parameters of air and water

船-冰-水耦合数值模型如图3 所示

图3 船-冰-水耦合模型Fig. 3 Ship-ice-water coupling model

3 数值模拟结果

3.1 船-冰-水耦合结果分析

船-冰水相互作用过程，水阻力与冰载荷是船体载荷响应的主要成分。首先对敞水环境下水阻力进行研究，对船体施加纵荡方向的速度模拟船舶在恒定功率下匀速前进。航速控制分别为2 kn，5 kn，8 kn，10 kn，15 kn。在Ls-dyna 可以实现船舶前进过程中的兴波现象模拟，如图4 所示。船舶首部型破开水域前进，将水向两边推开，自由液面发生改变，波高出现区域船身压力增大，波谷区域船身压力降低，形成兴波现象。随着船舶继续前行，在船体水线下部分出现扰流，最终在尾部形成尾波现象，从首波、船身两侧的V-形扩散波和尾波轨迹与开尔文船行波系一致[13]。

图4 舰船航行数值模拟结果Fig. 4 Numerical simulation results of ship navigation

图5 开尔文波波峰轨迹Fig. 5 Kelvin wave crest trajectory

表3 海冰力学参数Tab. 3 Mechanical parameters of sea ice

图6 为15 kn 航速下船体外板的应力云图，水线下部分应力随着水深而逐渐变大，船身表面最大应力出现在水线以下首柱破波区域和球鼻艏处。

图6 水阻力作用下外板应力云图Fig. 6 Hull stress diagram under water resistance

图7 为不同航速下敞水阻力时历曲线，对数据进行滤波处理来减小振荡幅度。可以看出，在0～2 s内，船体在推力作用下，导致静水扰动，水阻力值周期振荡提升，3 s 后当水域趋于稳定后，船体表面受力趋于固定值，逐渐收敛。对收敛后水阻力值与航速关系进行分析，如图8 所示，当航速从2 kn 逐渐增加至15 kn，水阻力值从3.5×105N 增加至7.5×105N，本船未进行船体模型水阻力试验，将结果与S60 系列船型通过三因次法换算得到的船体总阻力对比，阻力值量级相同。并且在Fr值相同的情况下，水总阻力预报值与船体航速存在正相关的趋势。

图7 不同航速下水阻力时历曲线Fig. 7 Time history curves of water resistance

图8 水阻力稳定值随航速变化Fig. 8 Variation of water resistance with speed

在确定船-水耦合下的水阻力规律后，对冰水耦合下的水阻力、冰阻力作用机理进行研究。在多相耦合时，船体速度不再以恒定功率下航行，给船舶施加10 kn 的初速度，拥有初始动能的船首与浮冰发生冲撞，浮冰的存在导致水阻力与冰载荷作用机理之间存在相互影响。

图9 为船-冰-水三相耦合下的作用过程，船与冰接触发生碰撞。在碰撞力作用下，产生排水效应，在浮冰周围形成环形波浪，同时浮冰浮态发生改变，在波浪中呈现动态平衡效果。

图9 船-冰-水相互作用数值模拟Fig. 9 Numerical simulation of ship-ice-water interaction

对碰撞时间内的冰载荷时历曲线进行分析，浮冰与船体发生碰撞时，首柱处受到冰载荷的碰撞与挤压作用，并分为与船体运动方向相反的纵向力、沿船宽方向对称分布的横向力及沿Z轴向上的垂向力。时历曲线如图10 所示，而冰单元在碰撞后发生局部破碎失效导致峰值出现多次卸载现象。从时历曲线可发现，纵向力是冰阻力的主要组成部分，并且在时间上先于横向力出现。因为浮冰挤压作用是关于船线对称，横向破冰力表现为沿0 附近的振荡状态。船体在与浮冰接触过程中，船体自身重力对浮冰存在压溃作用，在冰阻力上表现为持续存在的垂向冰载荷。

图10 冰载荷时历曲线Fig. 10 Time curve of ice load

将冰-水耦合作用下的冰载荷与无水域状态时约束Z向自由度的冰载荷进行对比，如图11 所示。水域浮力的作用下弱化了船-冰碰撞时的接触，使得峰值数据较无水状态下，减小了32%，整体载荷趋势也有所下降。因此计算船冰碰撞时，水域的影响效果是不能忽视的。

图11 冰阻力时历曲线对比Fig. 11 Comparison of ice load time curve

船-冰-水耦合作用下的水阻力时历曲线如图12 所示，呈逐渐衰减直至消失趋势。这是因为水阻力在船舶与浮冰发生碰撞后，船体动能一部分转化为浮冰内能与动能，转化为水域内能部分减少。同时随着船体速度下降，击水作用力与排水速度降低，导致船体表面水阻力逐渐衰减。

图12 水阻力时历曲线Fig. 12 Time curve of water resistance

3.2 冰型尺寸影响分析

浮冰厚度与尺寸是规则浮冰与船舶首部碰撞的重要冰型参数，在船-冰-水相互作用机理研究的基础上，改变浮冰厚度或浮冰尺寸，展开初速度为10 kn 下的碰撞仿真分析。工况设置如表4 所示，水域及碰撞参数设置不作改变。

表4 船-冰碰撞计算工况Tab. 4 Ship-ice collision calculation conditions

分别以冰厚和浮冰尺寸为变形分析case1～case7 工况船体所受冰阻力与水阻力，各冰厚和浮冰尺寸冰阻力时历曲线对比如图13 所示。

图13 不同冰型参数下冰阻力时历曲线Fig. 13 Ice load time history curves under different ice type parameters

浮冰厚度越小，冰单元破坏程度越大，失效单元数越多，峰值卸载频率越高，峰值与冰厚变化成正相关关系；冰阻力在浮冰边长达到20 m 与50 m 时出现畸变，船-冰作用时间随浮冰尺寸的增加呈现增加趋势。当浮冰尺寸为50 m 时，冰阻力持续存在，这是因为船舶本身动能不足以使巨型浮冰产生明显速度变化，反之船舶受到浮冰质量的阻力影响。其初速度下降更为明显，二者之间一直处于接触状态，挤压与碰撞冰载荷也一直存在。不同冰厚下的冰载荷峰值如表5 所示，不同浮冰尺寸下的载荷峰值如表6 所示。

表5 不同冰厚下的冰载荷峰值Tab. 5 The peak ice load under different ice thickness

表6 不同浮冰尺寸下的载荷峰值Tab. 6 The peak ice load under different floating ice sizes

水阻力与船舶航速关系密切，水阻力值几乎不随浮冰参数发生改变，趋势与图12 基本相同，仅随船舶速度衰减而逐渐减小。从冰载荷在总载荷值中占比发现，浮冰厚度为0.5m 时，冰载荷峰值为2.11MN。当浮冰厚度增加1 倍、1.4 倍和2 倍时，峰值分别增加29%、65%与138%。浮冰边长增加给浮冰质量的增加是二次方级的，因此，冰载荷增长速度存在跃迁趋势，且当浮冰尺寸大于型宽时，增长趋势减缓，达到冰载荷阈值。假设浮冰尺寸无限增加，即可视为覆盖在海面的层冰，不产生移动速度，船体初始动能对浮冰的作用效果几乎完全转化为浮冰破碎失效的变形能。而浮冰越厚，尺寸越大，水域对冰载荷的弱化效果越弱，冰载荷在总载荷中占比越高。

图14 冰阻力峰值随冰型参数变化Fig. 14 The peak value of ice load varies with ice parameters.

4 结 语

为了研究舰船在航行时遭遇浮冰，载荷响应情况，本文基于S-ALE 流固耦合方法对冰区海域浮冰漂浮和船-冰-水相互作用进行数值模拟，并对不同冰型参数下载荷影响进行了分析。主要结论归纳如下：

1）通过对船-冰-水相互作用过程中船行波以及船-冰碰撞现象的准确模拟，验证了S-ALE 流固耦合方法处理船-冰问题的可行性。

2）船舶敞水域航行过程中，当航速达到稳定后，会产生开尔文船行波，并且水阻力趋于稳定值，与航速存在正相关关系。同时将船-冰-水相互作用的冰载荷结果与无水状态下约束浮冰Z向自由度的浮冰相比，冰载荷峰值减弱了32%，证明了在考虑船-冰作用时水域的重要性。

3）水阻力仅与航速相关，而冰载荷受到水域外部因素及冰厚、浮冰尺寸等自身因素的影响。水域在船冰碰撞过程对冰载荷起到了削弱作用，而这一削弱效果随浮冰尺寸、冰厚的增加而减弱。浮冰尺寸、冰厚越大，冰载荷在舰船航行总载荷中占比越大，且尺寸因素影响大于浮冰厚度。