航海模拟器中破冰船破冰阻力数值计算

王晓雪, 张秀凤, 刘兆春, 孟耀

(1.浙江师范大学 工学院,浙江 金华 321004;2.大连海事大学 航海学院,辽宁 大连 116026)

随着北极航道的开通,极地船舶航行安全极为重要。由于北极海域环境恶劣,海冰成为影响极地船舶安全航行的最主要威胁。冰区航海模拟器在船舶驾驶员培训中起到重要作用,其关键技术之一是冰区船舶运动数学模型。而冰区船舶运动数学模型的建立需要确定船舶与冰碰撞时的受力。因此,研究船舶与冰碰撞时的受力具有十分重要的意义。

对破冰阻力的预报常采用理论方法、数值方法和试验方法。试验方法最简单直接,但是不足之处是海冰与海洋结构物之间的相互作用极其复杂,实船观测的数据有限,且这些数据受到船型和海域的影响,会受到实验条件和测试水平的限制,不能被广泛应用。而经验公式与数值模拟的应用范围更为广泛,实用性强,可在大量假设条件下以及近似情况的基础上进行。所以学者们一般都采用理论与经验公式相结合的方式来研究。

Lindqvist[1]将总阻力划分为3个部分:破冰阻力、水下阻力和速度相关的阻力,根据试验数据提出半经验公式对冰阻力进行预报。Riska方法[2]是在一系列经验数据基础上,根据基础理论相结合推导出来的公式。其中经验系数是根据波罗的海的不同船型的多组实船试验得到的。Jeong方法[3]是基于模型试验数据,对Spencer所建立的近似估算方法进行了验证并推广。王超[4]利用经验方法对连续破冰过程中船舶的冰阻力进行估算,对平整冰中船舶操纵性能进行了研究。

与理论方法和试验方法相比,数值模拟方法作是目前最主要的研究方法,Luo等[5]考虑了船-水-冰的相互作用,利用CFD-DEM耦合数值方法对冰区航道航行的船舶阻力进行了研究,并与实船试验数据作了对比。Zhang等[6]通过拖曳水池试验和计算流体力学软件,确定了船体斜拖MMG模型中的水动力和力矩,得到了船体4种前进速度下的水动力导数,最后总结了水动力导数与前进速度之间的关系。Su等[7]采用数值方法对船舶在平整冰中的操纵性作了相应预报,并针对在船舶与平整冰的接触过程中分析了利用船体重力将海冰压碎的过程。Su等[8]研究了船舶在冰区航行的推力,研究结果表明船舶的推力在冰区中相较于开阔水域中有一定的减额,并对船体将平整冰压碎的过程进行了受力分析。Raed等[9]研究了船与冰相互作用过程的数值模拟,使用新的封闭式求解方式描述了海冰破裂过程。Wang等[10]采用有限元方法对航行于碎冰区的散货船进行了数值分析,模拟了4种速度和3种碎冰密度下的冰区船舶受力,分析了航行速度和碎冰浓度对冰载荷的影响。张健等[11]并对封冻在层冰中的船舶受到层冰挤压的损伤变形进行了分析,比较了不同结构形式的外板对船舶所受冰载荷的影响。刘为民等[12]应用LS-DYNA流固耦合方法对船舶在平整冰条件下航行的冰载荷进行了数值模拟,分析了速度和冰厚对破冰阻力的影响,发现在时域上各个方向的冰载荷都具有一定周期性,且破冰阻力随航速和冰厚的增加基本呈线性增加的趋势。

综上,目前大多学者利用数值模拟方法对船-冰碰撞过程开展研究,主要分析了破冰过程中不同船速、船艏形状、冰厚等对破冰能力的影响,对船舶不同角度与平整冰的碰撞以及船舶在冰区旋回航行的破冰能力研究的文献很少。本文基于ANSYS/LS-DYNA软件,对破冰船与平整冰的碰撞进行数值求解,并研究了船舶与平整冰在不同速度、不同碰撞角度和旋回角度以及不同的平整冰厚度发生碰撞时的动态响应特性,定性分析了不同工况下船冰碰撞的船舶的碰撞力、能量损耗等规律,并且对船速与冰厚2个因素进行了定量分析,拟合出公式,进一步为破冰船在冰区的破冰作业提供一定的理论参考价值。

1 船冰碰撞数值仿真理论

数值模拟法作为目前船冰碰撞问题的最主要研究方法,能够有效地计算出船舶在发生碰撞过程中船舶的破坏变形、碰撞力的变化和能量变化等各种响应。ANSYS/LS-DYNA是现阶段最优秀的显性动力分析的有限元软件,非常适合求解结构的非线性高速碰撞等动态冲击问题。

1.1 有限元分析方法

应用有限元法求解船冰碰撞非线性动力学问题时,系统的一般方程描述为:

(1)

当采用中心差分时间积分显式求解加速度向量时,则计算结构系统各节点在第n个时间步结束时刻tn的加速度矢量可描述为:

a(tn)=M-1[P(tn)-Fint(tn)]

(2)

式中Fint为内力矢量,其计算为:

(3)

式中:BTσdΩ、Fhg、Fcontact分别是当前时刻单元应力场等效节点力、沙漏阻力和接触力矢量。

在系统中的节点速度和位移矢量为:

(4)

时间步和时间点定义为:

Δtn-1=(tn-tn-1),Δtn=(tn+1-tn)

(5)

(6)

下一个模型构型为:

xt+Δt=x0+ut+Δt

(7)

1.2 LS_DYNA时间步控制

LS_DYNA软件显式有限元方法采用中心差分法离散化时间来求解运动微分方程。LS_DYNA软件会通过检查所有单元来自动计算所需时间步长:

Δt=0.9l/c

(8)

式中c为材料声速,计算公式为:

(9)

l为实体单元特征长度,计算公式为:

l=V/Amax

(10)

式中:Amax是单元的最大面积;V为单元的体积。

2 船冰碰撞数值模拟

2.1 有限元模型建立

“雪龙2号”极地考察船是我国第1艘自主建造的极地科学考察破冰船,是全球第1艘采用船艏、船艉双向破冰技术的极地考察船。本文以“雪龙2号”主尺度为原型,绘制一艘用于破冰作业的仿“雪龙2号”的船舶模型。只考虑船艏破冰[13],因此对船艉做流线型简化,如图1所示。

为了降低对计算机性能的要求,在保证精度的同时,缩减仿真计算运行时间,验证试验方法可行性。本文参照“雪龙2号”船艏几何形状,建立了简化的楔形体模型,进行船舶破冰的仿真试验。利用楔形体简化模型,通过改变多项影响因子,对船冰碰撞进行多次模拟,得到大量数据。根据数据对各项参数进行灵敏度分析,总结各项参数影响下的冰阻力变化规律。图2、3为船-冰的有限元模型。楔形体模型与平整冰的主要几何参数如表1所示。

表1 楔形体与平整冰的主要参数Table 1 Main dimension of wedge and level ice m

图2 破冰船-平整冰有限元模型Fig.2 Finite element model of icebreaker and level ice

图3 楔形体与平整冰有限元模型Fig.3 Finite element model of wedge and level ice

在船体破冰模拟中,船体有限元模型采用壳单元,网格采用三角形网格,尺寸长度50 cm(为避免网格畸形,防止在仿真计算过程中产生负体积,定义最小网格尺寸不小于10 cm);平整冰设定为无限冰区,为了提高计算效率,故平整冰采用中心加密方式,即平整冰采用沿船舶中纵剖面方向为中心线,依次向两侧稀疏的方式。其中最小网格尺寸为10 cm×10 cm×20 cm,如图4所示。

图4 “雪龙2号”与平整冰的网格划分Fig.4 The meshing of “Xuelong 2” and level ice

2.2 模型材料选取

1)破冰船材料选取。

本文采用塑性动态模型模拟破冰船船体材料,具体的塑性动态模型材料的参数设置如表2所示。

表2 塑性动态模型材料的参数设置Table 2 Parameter setting of plastic dynamic model material

2)平整冰材料选取。

平整冰材料模型的选择一直是船-冰碰撞仿真的重点与难点。如表3所示,本文选取各向同性的弹塑性断裂模型[14]模拟冰材料,失效准则为Von Mises屈服准则,破坏模式为最大塑性应变模式,分离模式为恒定最小压力模式。

表3 冰的材料模型参数Table 3 Material model parameters of level ice

2.3 船-冰碰撞模型数值计算

破冰过程主要分为破碎、翻转、滑行。设定船舶初始速度为1.5 m/s,与1 m厚的平整冰发生接触作用。如图5、6所示,平整冰首先被船艏撞出一个缺口,随着破冰船的持续破冰,缺口逐渐增大,直至达到冰层的失效应变或截断压力,平整冰破碎断裂。“雪龙2号”船舶破冰时,海冰的失效主要表现为弯曲破坏,如图7所示;楔形体破冰时,海冰破损主要为挤压破环,如图8所示。

图6 楔形体破冰应力应变效果Fig.6 Stress-strain effect of wedge ice breaking

图7 弯曲破坏Fig.7 Bending failure

图8 挤压破坏Fig.8 Extrusion failure

为检验仿真结果和数据的可靠性,本文采用破冰船Terry Fox号的船模水池试验[15-16]数据进行对比。由于Terry Fox船模与本文数值模拟所采用船模尺寸不同,需对两者所得出的结果进行无因次化处理。根据相似原理和量纲分析法,定义阻力的无量纲系数为:

(11)

(12)

(13)

式中:ρice、hice分别是冰的密度和厚度;Δρ是冰和水的密度差;BWL是水线宽;T是吃水;Rb、Rc、Rbr分别是冰浮阻力、冰滑阻力和破冰阻力;Cb、Cc、Cbr分别是冰浮阻力系数、冰滑阻力系数和破冰阻力系数。

总阻力系数的表达式为:

CIt=Cbr+Cc+Cb

(14)

绘制破冰船Terry Fox号的破冰阻力系数与冰厚傅汝德数Fri的关系如图9所示。由于本文在数值计算中没有考虑冰清阻力和冰滑阻力,在Fri=0.5时,误差最大为22%,具有较好的精度,说明了本文数值计算的可靠性。

图9 破冰阻力系数与冰厚傅汝德数关系Fig.9 Ice breaking resistance coefficient and Fri

如图10所示,破冰船与平整冰碰撞过程中,平整冰处于“接触-挤压-变形-失效-应力释放-再接触”状态。船艏与冰层首先单点接触,由于船艏接触点突然受到冰阻碍,造成暂时停顿,船体受到的冰阻力此时产生较大的峰值,时历第1个峰值后,船舶继续前行破冰,破冰力逐渐稳定在一个范围内;与此同时,冰层断裂过程中,冰层因受到船艏冲击在接触处产生微裂纹,船艏对冰层施力时,冰层内部的受力并不均匀,大部分会集中在微裂纹处。冰阻力到达峰值,造成进一步开裂,裂纹扩展是一个正反馈过程,即裂纹越长,应力集中效应越明显,进一步开裂也越容易。

图10 冰阻力Fig.10 Ice resistance

对该时间段内的冰载荷进行数据统计得到船体受到的平均阻力为6.358×106N。将数值模拟得到的破冰阻力与Lindqvist[1]经验公式及Jeong[2]经验公式进行对比,如图11所示,2种经验公式计算得到的冰阻力值随时间变化的总体趋势与数值模拟结果相符,Lindqvist经验公式的阻力平均值4.489×106N及Jeong经验公式的阻力平均值1.230×106N,与数值模拟得到的平均破冰阻力6.358×106N为同一个数量级,且误差值在20%以内。故本文对船舶破冰阻力预报数值模拟研究具有一定可靠性。

图11 经验公式计算的冰阻力Fig.11 Ice resistance calculated by empirical formula

3 船舶破冰阻力参数敏感性分析

破冰船进行破冰作业时,影响船-冰碰撞的因素有很多,本文选取不同平整冰厚度、不同船舶初速度以及不同船冰碰撞角度、旋回角度等易于量化的参数来分析对船舶破冰阻力的影响。

3.1 不同平整冰厚度

本文选取了平整冰厚度为0.5、0.8、1.0、1.2、1.5 m的5个不同的工况,破冰船以初速度为1.5 m/s向平整冰域航行破冰。破开冰的长度、动能、破冰阻力变化曲线分别如图12～14所示,破冰船进行破冰作业的同一时间内,随着平整冰厚度的增加,破冰船破开平整冰的长度渐渐变小,动能迅速下降,随着破冰船与平整冰的不断接触,冰力呈现出锯齿状的非线性特征波动,并伴随有多个峰值出现。

图12 不同平整冰厚度下破冰长度Fig.12 Ice breaking length under different ice thickness

图13 不同冰层厚度下动能Fig.13 Ship kinetic energy under different ice thickness

图14 不同平整冰厚度下碰撞力Fig.14 Force under different ice thickness

如图15所示,数值模拟得到的平均破冰阻力与由经验公式(见文献[1-3])计算得到的平均阻力值进行对比,数据结果如表4所示,与Lindqvist经验公式计算所得平均值相对误差约为10.85%,与Riska[3]方法的相对误差约为28.21%,与Jeong公式的相对误差约为8.15%。

表4 不同冰厚的平均破冰阻力值Table 4 Average ice breaking resistance under different ice thicknesses 106 N

图15 不同冰厚的平均破冰阻力Fig.15 Average ice breaking resistance under different ice thickness

由于本文的仿真试验忽略了破碎浮冰在滑行阶段对船体造成的摩擦阻力,仅考虑了海冰环境干扰力,所以数值模拟得到的值比经验公式计算所得值偏小,但误差仍在一定范围内,故本文对船舶破冰阻力数值模拟预报具有一定可靠性。

3.2 不同碰撞速度

破冰船以不同初速度1、1.5、2、2.5 m/s冲撞冰层时,破冰阻力和动能时历曲线如图16、 17所示,由于在船-冰碰撞过程中,不断发生“碰撞-平整冰破碎-碰撞-平整冰破碎”这一现象,碰撞力曲线呈现出了一定程度内的波动。

图16 不同速度下冰阻力Fig.16 Ice force at different speeds

图17 不同速度下动能Fig.17 Kinetic energy at different speeds

3.3 参数灵敏度

将上述得到的有关不同冰厚与不同速度得到的冰阻力如表5所示的数据,为利于冰区航海模拟器数学模型的实时解算,采用矩形域的最小二乘法曲面拟合方法,得到适用于“雪龙2号”船舶运动预报的阻力表达式为:

表5 不同工况下船-冰碰撞计算冰载荷Table 5 Ice load in ship ice collision under different working conditions 106 N

(15)

3.4 不同碰撞角度

船舶在无限冰区航行时,不管船艏对着哪个方向都是直航破冰,但在要进入冰区航行时,破冰船船艏与所要接触的平整冰区所呈角度不同,对碰撞的影响也不同。如图18所示,破冰船船身与平整冰边缘成θ角,船速V为2 m/s破冰。

图18 船舶斜向破冰Fig.18 Oblique ship ice breaking

改变角度θ,取值为60°、65°、70°、75°、80°、85°。如图19所示,破冰船1/4船身进入平整冰区域的碰撞力时间历程曲线。船舶冲撞式破冰,撞击力过大,导致平整冰破损缺口过大,船舶行进时暂时与冰层无接触,致使有冰力为0的状态,而且船舶由于船体左右舷两侧受到的冰载荷大小差值过大,碰撞力大幅度波动。

图19 船舶斜向破冰的碰撞力Fig.19 Force of ship oblique ice breaking

在破冰船以θ角在冰区航行时,将船冰碰撞过程中最大的碰撞力作为破冰船以θ角航行时的碰撞力。如表6所示,船冰碰撞角度为90°时,碰撞力最大,破冰船与平整冰成直角撞击为最危险工况,由于船舶设计时,一般考虑最危险工况,故本文其他工况船舶与平整冰碰撞的角度均为90°。

表6 不同破冰角度的碰撞力Table 6 Force of different ice breaking angles

3.5 不同旋回角度

如图20所示,本文设计破冰船船身与平整冰成直角,给予船体前进速度Vy和横向速度Vx,使破冰船以合速度2 m/s旋回破冰。

图20 船舶旋回破冰Fig.20 Ship turning ice breaking

但由于破冰船前期以合速度∑V斜向上与平整冰碰撞,如图20的虚线船身所示,导致船艏右侧受到的冰力远大于左侧,致使破冰船向左转向旋回,直至船舶停止。这一过程不易界定船舶的破冰阻力,故将该过程的平均阻力作为船舶以α角在冰区航行时的破冰阻力。

船舶旋回破冰作业时,改变角度α,取值为5°、10°、15°、20°、25°、30°、10 s内的碰撞力时间历程曲线如图21所示,可见碰撞力随着船舶逐渐进入冰区,船舶与平整冰接触碰撞的面积逐渐增大而增大。在前期船舶向右前方破冰,碰撞力缓慢增加;在约5～6 s之后,船舶开始在冰区旋回破冰,这个阶段的碰撞力时间历程曲线与图19的碰撞力状况相似,运动趋势一致。

图21 不同破冰角度的碰撞力Fig.21 Force at different ice breaking angles

船舶在冰区旋回破冰的这7个工况的碰撞力值如表7所示,将船舶旋回破冰的碰撞力与直航连续式破冰的碰撞力对比,发现船舶在冰区旋回破冰比直航连续式破冰受到的碰撞力更大。碰撞力随着旋回角度的增大,也小幅度的增大。

表7 不同破冰角度的碰撞力Table 7 Force at different ice breaking angles

由于冰区船舶操纵模拟器行为真实感的实时性的要求,依据上述有限元软件数值计算出的破冰阻力值,进行了适合于航海模拟器实时性要求的处理:1)拟合出适合模拟器用的较高精度的实时估算的阻力公式;2)建立了适用于运动数学模型实时解算的破冰阻力离线数据库,利用插值方法,实时获得冰阻力,并在冰区船舶操纵模拟器中得到了应用。

4 结论

1)利用Terry Fox号的船模水池阻力试验结果与“雪龙2号”船舶的数值计算的无量纲化结果比较,验证了数值计算的可靠性,并将仿真得到的阻力平均值与3种经验公式计算的平均冰阻力值进行了验证对比,进一步说明了数值计算的可信度;

2)利用“雪龙2号”简化后的楔形体模型,通过LS_DYNA软件对楔形体模型以不同的碰撞速度、碰撞角度、旋回角度以及平整冰厚度,进行多次船-冰碰撞数值计算,利用矩形域的最小二乘法定量拟合了航速、冰厚的冰阻力计算公式,适用于船舶低速破冰的阻力值计算,本文拟合出的估算公式可应用到航海模拟器中破冰船开辟航道的作业模拟中。

3)相比经验公式方法,数值计算在可视化、船舶初始环境设定等问题上有很大优势,今后将利用破冰船实船试验数据进行仿真验证,以提高数学模型的仿真精度。