航行体出水破冰的多场耦合效应与相似律1)

岳军政 吴先前,2) 黄晨光

* (中国科学院力学研究所,中国科学院流固耦合系统力学重点实验室,北京 100190)

† (中国科学院合肥物质科学研究院,合肥 230031)

引言

高速航行体出水破冰具有重大的工程应用背景.与常规水下发射过程相比,航行体需要穿越水面覆盖冰层或穿越含冰块结构的冰水混合区域,是一个涉及多相多介质动力耦合的复杂过程.该过程中,航行体的受载环境及弹道稳定性与在纯水介质出水过程有很大差别,例如冰载荷的增加,冰层受到冲击后复杂的断裂破碎行为,都将直接导致航行体和冰面的受力分布发生剧烈变化,增大出水冰发射的不确定性.因此,开展航行体出水破冰的动力学特征与载荷作用机理研究,具有重要的应用价值和科学意义.

目前冰与结构相互作用的研究主要集中在冰的力学性能、船−冰碰撞方面,而结构物穿冰作用的研究还较少.在冰的力学性能方面,温度和应变率是两个重要的影响因素[1-6],例如,Arakawa 和Maeno[1]的研究表明,在准静态单轴压缩条件下,0～ −100 °C时,多晶冰的抗压强度随着温度降低而升高,但是温度在 −100～ −175 °C 时,抗压强度并没有明显变化.Wu 和Prakash[7]通过分离式霍普金森压杆(split Hopkinson pressure bar,SHPB)研究了冰的动态压缩强度,实验结果表明,在 −15～ −125 °C 范围内降低温度时,冰的强度从32 MPa 增大到112 MPa,而在−125～ −173 °C 范围内继续降低温度时,强度大约为110～ 120 MPa,并无明显增强.在应变率效应方面,Jones[8]对温度 −11 °C 时应变率10−1～ 10 s−1范围内多晶冰力学性能的研究表明,冰的强度随应变率的增大而增大,并建立了压缩强度与应变率的近似关系式.季顺迎等[9]则结合Hibler 黏塑性本构模型和K-V 体黏弹性模型,建立了海冰的黏弹塑性本构模型,兼顾小应变率下的弹性以及大应变率下的塑性,以模拟海冰的相互作用.Gao 等[10]通过建立新的失效准则,结合“Tsai-Wu”屈服模型建立了海冰的弹塑性本构模型,以模拟船−冰山碰撞过程.在冰的破碎失效方面,静态和准静态压缩条件下,冰的塑性应变可超过0.1,翼型裂纹萌生、扩展及导致的劈裂破坏是主要的失效机制.随着应变率的增大,冰表现出由韧转脆的现象[11].在更高应变率条件下,冰的失效模式以大块的劈裂转向了碎片和碎粒,Shazly 等[12]采用SHPB 对单晶和多晶冰在应变率60～ 1400 s−1范围内的力学性能进行了研究,实验中单晶和多晶冰均碎裂成很多细小颗粒,并且冰失效后具有较高的残余强度.Carney 等[13]发展了高应变率下冰的失效模型.冰载荷是船冰碰撞过程中的研究重点,因此国内外学者对冰破坏的物理过程开展了大量的试验和数值模拟研究[14-17].例如,Sodhi[18]对冰层进行了小尺度压痕试验,发现低压痕速度下冰发生延性变形,主要是蠕变,而高速下发生连续脆性压碎,并且断裂前的冰变形大多是弹性的.郭春雨等[19]对冰区船舶在碎冰中的航行过程进行数值模拟,获得的航行阻力和实验结果在定性上保持一致.吴炜和黄焱[20]通过开展冰水池模型试验,研究了海冰漂移和船速对冰载荷分布的影响.程文霞[21]则通过在近场动力学模型中考虑爆炸破冰中的气泡载荷作用,对涉及的载荷形式进行了分析.

目前关于航行体出水破冰的研究鲜有报道,因此,本文基于量纲分析,系统地开展航行体出水破冰数值模拟,研究冰−水−结构多介质耦合作用机理,以获得航行体的运动特性及动载荷特征,为复杂冰区环境下航行体的高速出水发射提供理论依据.

1 量纲分析

为了系统地开展数值模拟研究,首先对航行体出水破冰多介质耦合过程(见图1 示意图)进行量纲分析.对于航行体的速度v、过载a以及头部应力σ,控制参数主要来自以下5 个方面.

图1 航行体出水破冰示意图Fig.1 Schematic diagram of ice break by vehicle launched underwater

(1) 时间t;

(2) 冰的参数:密度ρi,厚度hi,弹性模量Ei,泊松比υi,极限应变εimax,特征尺寸li;

(3) 水的参数:密度ρw,声速cw;

(4) 航行体与冰的初始距离s;

(5) 航行体参数:密度ρp,弹性模量Ep,泊松比υp,长度lp,直径dp,初始速度v0.

1.1 航行体速度v

v应当是上述5 组控制参数的函数,即

取ρp，Ep和dp作为基本量,式(1)可化为无量纲关系

若材料不变,并且保持航行体几何形状相似,式(2)可简化为

1.2 航行体过载a

a也是上述5 组控制参数的函数,并且若保持材料不变以及航行体几何形状相似,函数有以下形式

1.3 航行体头部应力σ

σ也是上述5 组控制参数的函数.同样地,若材料不变,并且保持航行体几何形状相似,函数有以下形式

为了验证上述推导的量纲分析结果,采用LSDYNA 对两种不同参数条件下的航行体出水破冰工况分别进行了数值计算.两种工况中保持材料相同以及航行体长径比10∶1 不变,第一种工况航行体直径dp=10 mm、冰层厚度hi=3 cm、冰层为无限大冰层即li趋于无穷,航行体与冰的初始距离s=5 cm,初始速度v0=40 m/s;第二种工况航行体的直径相对第一种增大1/3,即dp=13.33 mm,其余量也相应增大或不变,保持和第一种工况的无量纲参数相同.采用的计算模型及方法具体在下文2.1 介绍.图2 分别为计算的两种工况的无量纲航行体速度v*、过载a*以及头部应力σ*随无量纲时间t*的变化,由图2 可以看出,两种工况的无量纲结果吻合得较好,验证了上述推导的无量纲关系的正确性.

图2 两种工况数值计算结果对比(上角标*表示各参数对应的无量纲量)Fig.2 Comparison of the numerical results for the two cases(superscript* denotes the dimensionless quantity)

由量纲分析可以看出,除了时间t外,影响航行体速度v、过载a和头部应力σ的因素主要为无量纲的冰厚度、冰特征尺寸、航行体与冰的初始距离s*以及航行体初速度.考虑到关于航行体水下运动的研究已经很多并且相关机制较为清楚[22],本文不着重于航行体的水下运动过程,在数值模拟研究中保持航行体与冰的初始距离s*不变,重点关注其余3 个因素的影响规律.

2 流固耦合计算方法

2.1 计算模型

考虑到问题的轴对称特性,建立1/4 模型进行计算.图3 为计算模型以及局部的网格划分,模型包括航行体、水域、空气域以及冰层,均采用Truegrid进行六面体网格建模.为了提高计算效率,对称轴附近的水、冰、空气采用精细网格,而外围部分采用较大的渐变网格.在对称轴附近主要的流固耦合作用区域,水、冰的网格为0.5 mm;空气网格横向为0.5 mm,轴向为0.75 mm;航行体的网格为水、冰网格的1.6～ 1.7 倍,约0.83 mm.2.2 节给出了具体的网格无关性验证结果.对称面上采用对称边界条件,水和空气域的外围以及上下面采用无反射边界条件.对于冰层,无限大工况下在外围采用无反射边界条件,而碎冰情况下采用自由边界条件.计算单位制采用cm-g-μs.

图3 计算模型及局部网格划分(不同工况内冰层尺寸不同)Fig.3 Full view of simulation components and part of meshed model(ice sheet size is different in different cases)

各介质的计算域、计算模型以及参数分别列于下表1～ 表4 中.固体介质航行体和冰采用拉格朗日(Lagrangian)算法,并且两者之间定义侵蚀面−面接触;流体介质水和空气采用任意拉格朗日−欧拉(ALE)算法.通过*CONSTRAINED_LAGRANGE_IN_SOLID 定义流固耦合关系.ALE 算法的控制方程为[23]

表1 冰计算模型及参数[24]Table 1 Ice details[24]

表4 空气计算模型及参数[25]Table 4 Air box details[25]

式中,ρ为密度,vi为拉格朗日坐标下的速度,σij为应力张量,fi为单位质量体力,ci为对流速度,qi为热通量,e为内能.

表2 航行体计算模型及参数Table 2 Vehicle details

表3 水计算模型及参数[25]Table 3 Water box details[25]

由于冰体在高应变率载荷下呈现脆性,因此参考文献[24-25],冰采用弹脆性模型以及失效算法[24],航行体材料采用DT300 钢,水、空气分别采用标准大气压下的参数[25].

本文中水采用Gruneisen 状态方程[23],其表达式如式(9)所示,具体参数列于表5 中.空气采用线性多项式状态方程[23],表达式如式(10)所示,具体参数列于表6 中.

表5 水的Gruneisen 状态方程参数Table 5 Gruneisen EOS parameters for water

表6 空气线性多项式状态方程参数Table 6 Polynomial EOS parameters for air

式中,μ=ρ/ρ0− 1,ρ和ρ0分别为水的密度和初始密度.C为vs(vp) 曲线的截距,S1～S3是vs(vp)曲线斜率的无量纲系数,γ0为常数,E为比内能.

式中,μ=ρ/ρ0− 1,C0～C6是与气体性质有关的常数,E为比内能,ρ和ρ0分别为空气的密度和初始密度.

2.2 网格无关性验证

为了验证计算网格的无关性,将2.1 节表格所列的网格均缩小为原来的0.75 倍,分别采用这两套网格,对初速度40 m/s 的航行体出水破冰(3 cm 厚无限冰层)进行了数值计算.

图4 为计算的航行体速度以及加速度的变化历史,由图4 可以看出,无论是水中运动还是破冰过程,两套网格计算的结果基本一致,表明可以采用表1～表4 所列的以0.5 mm 为主的一套网格进行数值计算.

图4 不同网格计算的航行体速度和加速度变化Fig.4 Simulated vehicle speed history and acceleration history by different grid sizes

3 计算结果与分析

基于量纲分析,本文针对表2 所示航行体,并且保持航行体与冰的初始距离s=5 cm 不变,通过改变航行体初始发射速度v0、冰的厚度hi以及冰的尺寸li(具体值列于表7 中),对18 种不同工况的航行体出水破冰进行了数值模拟.

表7 不同工况计算条件Table 7 Calculation conditions for different cases

图5 分别为计算的不同冰层条件下初速度40 m/s 的航行体的速度变化和加速度变化历史.由图5 可以看出,当航行体在水中的运动流场稳定后,其加速度约200g～ 300g.约1.15 ms 时刻,航行体撞击冰层,此时航行体的速度为35.6 m/s.对于无限大的冰层,1 cm 厚的冰并没有对航行体造成较明显的载荷,和在水中运动时的速度衰减趋势基本相同,航行体穿冰后的剩余速度约34.8 m/s;而对于3 cm 和5 cm厚的冰层,它们对航行体造成的平均加速度约220g,其中的细小差距是由于冰层的破碎形态不同造成的,航行体穿越3 cm 冰层后的剩余速度约33.7 m/s,而穿越5 cm 冰层后的剩余速度约32.8 m/s.

图5 不同冰层时初速度40 m/s 的航行体的速度和加速度变化Fig.5 Speed and acceleration history for 40 m/s vehicle impacting with different ice targets

对于尺寸为6 cm×6 cm 的碎冰,当冰厚度为1 cm或5 cm 时,其对航行体运动的影响和无限大冰层基本相同,只有瞬时的加速度会因为冰的破碎不同有一定差别;而对于3 cm 厚的碎冰,其对航行体速度和加速度的影响不及无限冰层,平均加速度约180g,航行体剩余速度约34.3 m/s.这是因为,对于较厚的冰,无论是无限大冰层还是碎冰,在水的冲击作用下损伤较小,因此航行体撞击前其完整性均较好,两者对航行体运动特性的影响差别不大;而对于较薄的冰,无论是无限大冰层还是碎冰,在水的冲击作用下几乎完全断裂破碎,因此对航行体运动特性的影响都很小,冰的尺寸效应也不明显.只有对于3 cm 厚的冰,水对无限大冰层或碎冰的冲击破坏程度区别较大,导致两者后续与航行体的撞击作用出现明显差别.

图6 为计算的初速度40 m/s 的航行体撞击不同冰层时头部的应力变化历史.对于无限大冰层,航行体撞击1 cm 厚的冰产生的最大头部应力仅有2.2 MPa,这是薄冰在撞击前损伤较大的缘故.撞击3 cm 和5 cm 厚的无限大冰层时,产生的最大头部应力基本相同,约12.3 MPa,且在穿越5 cm 厚的冰时,出现约1 ms 的稳定侵彻阶段,以上现象表明:当无限大冰层超过3 cm 厚时,航行体头部应力只和其速度有关,而和冰的厚度无关.对于大小为6 cm×6 cm 的碎冰,当厚度为1 cm 或5 cm 时,其在航行体头部产生的最大应力和无限大冰层基本相同,而厚度为3 cm 时,相对无限大冰层,碎冰产生的最大头部应力较小,仅有7.5 MPa.

图6 不同冰层时初速度40 m/s 的航行体的头部应力历史Fig.6 Stress history at the 40 m/s vehicle head with different ice targets

图7(a),图7(b)分别为计算的初速度30 m/s 的航行体的速度变化和加速度变化历史.当航行体在水中的运动流场稳定后,其加速度约90g～ 180g.约1.65 ms 时刻,航行体撞击冰层,此时航行体的速度为26.6 m/s.由图7 中计算结果可以看出,对于3 cm和5 cm 厚的冰,无论是无限冰层还是6 cm×6 cm碎冰,它们对航行体速度和加速度的影响基本相同,加速度约130g,呈现稳定侵彻现象,航行体穿越5 cm冰层后的剩余速度约24 m/s,穿越3 cm 冰层后的剩余速度约24.9 m/s.而1 cm 厚的无限大冰层或碎冰,对航行体运动特性的影响均较小,加速度仅有100g,航行体穿冰后的剩余速度约25.7 m/s.

图7 不同冰层时初速度30 m/s 的航行体的速度和加速度变化Fig.7 Speed and acceleration history for 30 m/s vehicle impacting with different ice targets

图8 为计算的不同冰层条件下初速度30 m/s的航行体的头部应力变化历史.由计算结果可以看出,对于5 cm 厚的无限大冰层、碎冰以及3 cm 厚的无限大冰层,它们在航行体头部造成的平均应力基本相同,约7.7 MPa,明显低于初速度40 m/s 时的最大应力12.3 MPa,并且3 cm 厚的碎冰可以造成相当的最大头部应力,只是其应力脉冲略小;而对于1 cm 厚的无限大冰层或碎冰,其在航行体头部产生的最大应力为3.8 MPa,高于初速度40 m/s 时的最大应力2.2 MPa,这是因为对于1 cm 厚的薄冰,航行体速度越大,其通过水的运动冲击越容易使冰层中间开裂上鼓,所以造成航行体的头部应力较小.

图8 不同冰层时初速度30 m/s 的航行体的头部应力历史Fig.8 Stress history at the 30 m/s vehicle head with different ice targets

图9(a),图9(b)分别为计算的初速度20 m/s 的航行体的速度变化和加速度变化历史.航行体在水中的运动流场稳定后,其加速度约50g～ 90g.约2.45 ms 时刻,航行体撞击冰层,此时航行体的速度为17.3 m/s.由计算结果可以看出,不同厚度冰层对航行体加速度的影响差别较小,平均加速度约60g～80g.无论对于碎冰还是无限冰层,航行体穿越1 cm，3 cm 和5 cm 冰层后的剩余速度分别为16.7 m/s，16.3 m/s 和15.8 m/s.

图9 不同冰层时初速度20 m/s 的航行体的速度和加速度变化Fig.9 Speed and acceleration history for 20 m/s vehicle impacting with different ice targets

图10 为计算的初速度20 m/s 的航行体撞击不同冰层时头部的应力变化历史.由计算结果可以看出,对于低速撞击,不同厚度冰层对航行体头部造成的最大应力基本相同,约4.5 MPa,仅有1 cm 厚碎冰造成的头部应力较小,只有3.1 MPa,表明航行体在该速度下的水下运动对1 cm 厚的无限冰层和碎冰造成的前期破坏不同.值得注意的是,1 cm 厚的无限冰层造成的航行体头部最大应力高于初速度40 m/s和30 m/s 时的结果,这同样是由于对于薄冰,航行体初速度越小,其通过水的运动对冰层中间的前期破坏越小,致使其头部应力较大.

图10 不同冰层时初速度20 m/s 的航行体的头部应力历史Fig.10 Stress history at the 20 m/s vehicle head with different ice targets

图11 为初速度30 m/s 的航行体出水穿越3 cm厚无限冰层时不同时刻的应力分布,为了清晰可见,隐去了水和空气流体介质单元,仅显示航行体和冰层固体单元.由图11 可以看出,1 ms 时刻,由于航行体高速运动引起周围的水介质向上运动,导致冰内部出现应力,并且由于局部单元失效出现裂纹;2.2 ms 时刻,冰层断裂为碎冰块,此刻航行体头部已进入冰层,出现稳定的侵彻阶段;3.2 ms 时刻,航行体头部已完全穿出冰层,此时冰内应力已经很小,仅有航行体内部应力波的传播和衰减.

图11 初速度30 m/s 的航行体撞击3 cm 厚冰层不同时刻的应力分布Fig.11 Stress distribution of vehicle with initial velocity 30 m/s and 3 cm thick ice target at different time

图12 分别为不同工况下航行体的无量纲最大过载和最大头部应力随无量纲初始发射速度的变化情况.由计算结果可以看出,航行体速度越大,不同冰层对其过载和头部应力的影响差别越大,这主要是因为在航行体速度较大时,不同冰层受到的前期破坏程度不同.从碎冰和无限冰层尺寸效应来看,无量纲厚度为3 的两种冰对高速航行体的过载和头部应力的影响差别较大,而无量纲厚度为1 和5 的碎冰与同厚度无限冰层的影响差别均较小.当冰层较厚、较大时,航行体速度越大,其过载和头部应力越大;而对于较薄的的冰层,由于其前期破坏较严重,所以航行体过载随速度增大而增大的趋势相对较小,尤其对航行体头部应力的影响规律性较差.

图12 不同冰层时航行体最大过载与最大头部应力随速度的变化Fig.12 Variation of maximum overload and head stress of vehicle with speed for different ice targets

4 结论

本文针对航行体出水破冰的应用背景,通过量纲分析,获得了影响航行体动载荷及其头部应力的主控参数和相似律.利用LS-DYNA,对不同初速度的航行体出水穿越不同冰层进行了数值模拟研究,获得了其动力学载荷特性及作用机理.计算结果表明,航行体速度越大,通过水介质对不同冰层的前期破坏差别越大,导致冰−水−结构耦合作用的影响效果差异明显.对于较厚的冰层,航行体穿冰后期呈现稳定的侵彻现象,航行体的载荷特征仅与其速度和冰的动态力学性能相关.对于较薄的冰层,需要考虑其撞击前的破坏程度.对于中等厚度的冰,在较高冲击速度下才表现出径向尺寸效应相关的载荷作用机制,低速条件下,径向尺寸大于6 倍航行体直径的碎冰和无限大冰层的影响效果基本相同.本文研究结果可为冰区环境下航行体的出水发射提供理论参考.