冰柱冲击问题的数值仿真分析

王 超， 杨 波， 张 媛， 郭春雨， 叶礼裕

(哈尔滨工程大学 船舶工程学院，哈尔滨 150001)

冰的撞击过程和随之引起的动态载荷是关系极地科学领域和冰区科学问题的重要研究内容之一，例如极地船舶的设计和建造、极区科考探索以及资源开采、航空和星际空间任务过程遭遇环境.因此研究冰在撞击过程的破坏过程以及动态载荷变化，一方面能够为上述实际问题提供安全作业依据，另一方面有助于为各种有与冰发生碰撞可能的结构物提供规范设计指导.

目前研究冰冲击破坏过程有试验研究和数值研究方法.试验方法能够揭示冰真实撞击过程的破坏机理，然而冰的冲击试验并没有标准的规范章程，因此冰模型的制备具有不统一性且很难定量把控试验过程的对比工况参数.可靠的数值模型计算方法可以弥补试验中的缺陷，研究定量参数对冰冲击过程的影响，基于商用LS-DYNA软件中网格方法的冰冲击过程预报是一种最为常见的数值方法[1-4].然而软件库中的本构模型并不能完全描述海冰复杂的力学性能，后期学者们在用户自定义模块改进了冰的本构模型[5-7]，增加温度、应变率及孔隙率等因素的影响.有学者基于黏聚单元模型构建了海冰的均化弹塑性本构，采用LS-DYNA计算了海洋结构物与平整冰相互作用的冰载荷[8-9]，然而传统的网格方法在处理断裂问题方面仍旧具有一定的局限性.后来，无网格粒子法逐渐被应用于冰的冲击破坏问题，例如光滑流体动力学(SPH)[10]和近场动力学 (PD)[11]方法.

近场动力学方法是近二十年新兴的一种非局部、无网格粒子、积分形式的计算方法，在断裂问题域和非断裂问题域上都可以连续求解[12-14].该方法最早由Silling提出并建立推导了其控制方程，给出了基本的数值实现程序案例[15].PD方法主要有两种表现形式，经推导和比较，两种表达方式大同小异且具有相类似的计算精度[16].PD发展至今主要有键型、常规状态型和非常规状态型3种形式，目前已成功应用于材料的弹性、塑性、热力学等多种物理问题计算分析领域[17-18].然而PD应用于冰计算领域刚刚起步，多数文献采用基于键型的弹性模型来模拟冰的破坏特性，忽略了泊松比的真实状态[19-20].在键型PD中，物体中两点之间的力被假定为位移和这些点的初始位置的函数.由于忽略了其他键变形对材料点对相互作用的影响，该公式存在如下缺点：泊松比的限制(对于二维平面应力情况，将泊松比固定为1/3，对于平面应变和三维情况，将泊松比固定为1/4)以及仅经历体积应变的多孔材料的塑性建模能力[21].

本文的主要目的是建立基于状态型PD方法的冰的弹性本构模型，消除键型对泊松比的限制，并研究和讨论冰本身物理参数、数值预报方法参数的敏感性特征.首先通过冰柱冲击试验验证了数值模型的准确性，分析了泊松比对冰体裂纹扩展的影响，然后进行了数值计算工况的变参数预报分析，即参数敏感性分析.

1 理论基础

1.1 常规状态型近场动力学

图1 物质点x与其临近物质点x′相互作用Fig.1 Interaction of a material point x with its neighboring points x′

考虑到力矢量状态，基于OSB-PD的积分微分方程描述为[23]

t′[x,t]〈x-x′〉]dVx′+b(x,t)

(1)

(2)

(3)

(4)

(5)

式中：V′、V分别为材料点x和x′的体积.

OSB-PD的最终表达式可以表示为

(6)

式中：a、b、d分别为PD专属参数；x′-x为变形前材料点相对位置；y′-y为变形后材料点相对位置；Λ为PD附加参数；θ和θ′分别为材料点x和x′处的膨胀变形；s为材料点之间变形伸长量(简称键).

1.2 刚性碰撞载荷

冰物质点和碰撞体有相对位移时，海冰物质点将与刚性物体表面发生接触.接触以后，物质点会渗入到刚体内，如图2所示.为反映真实的物理情况，渗入到刚体内的物质点需要进行位置再分配，将该物质点分配到临近处的刚体表面，如图2(c)所示.

图2 物质点的重新分配Fig.2 Relocation of material particle

新分配物质点k的坐标可以通过下式来计算：

(7)

(8)

(9)

(10)

(11)

1.3 失效准则

当键处于拉伸状态，s为正值.由于s不依赖于变形方向，故这种材料是各向同性的.为在本构模型中引入材料破坏的概念，使用一个比较简单判定条件，即某一时刻键伸长率超过一定值时，认为物质点之间的键永久发生断裂.从此刻起，此对物质点之间力密度为0.该条件称为近场动力学的破坏准则.

设μ是历史变形判断标准，当键伸长率达到破坏准则，键发生断裂，μ=0；反之，μ=1，可定义为

(12)

式中：s0为键破坏时键伸长率的极限值，可称为极限伸长率.对于典型的标准弹脆性材料(PMB)，s0可由下式计算得到[24]：

s0=

(13)

(14)

式中：G0为裂缝扩展时的能量释放率，可通过断裂韧度得到[25]；R为切变模量；K为体积模量；KI为冰的I型破坏断裂韧度；E为冰的弹性模量.

尽管弹脆性材料在初始条件下是各向同性的，但是有些键断裂之后会引起材料的各向异性，为此引入了键的破坏水平参数

(15)

2 数值执行过程

本文数值过程通过Fortran语言实现，OSB-PD的程序参考文献[26]中的键型PD的基本程序，在此基础上进行改进和扩展.

(16)

(17)

式中：j为领域粒子;s(k)(j)和Λ(k)(j)分别为粒子k和粒子j之间的伸长量和附加参量；V为当前材料点占据的体积.新增膨胀项通过式(18)的离散形式进行计算(即第k个粒子看i子的膨胀θk).

应变能密度的修正项Sm(k)为

(18)

整个数值实现过程如图3所示.

图3 数值模拟程序框图Fig.3 Flow chart for numerical simulation

3 冰柱冲击

3.1 计算模型与工况输入

本文选用Carney等[27]进行的冰柱冲击刚性板试验作为对照，建立了柱状冰冲击过程的近场动力学模型，如图4所示.模拟过程的参数设定如表1所示.

图4 柱状冰的近场动力学模型Fig.4 Peridynamic model of cylindrical ice

表1 数值模拟中的参数设定Tab.1 Parameter setting in numerical simulation

试验中，柱状冰直径D=18 mm，厚度H=42 mm.PD模型中柱状冰近似为各向同性匀质材料.数值过程中，冰柱离散为均匀分布的材料点.本文定义冰的力学参数为弹性模量E=1.8 GPa，泊松比υ=0.33，密度ρ=900 kg/m3，失效准则s0=0.003，对应断裂韧度K1=150 (kPa·m)0.5(s0与K1的关系将在后文详细说明).模型不考虑刚性界面的变形，并且同时忽略了重力、摩擦阻力以及气动载荷带来的影响.

3.2 收敛性分析与数值方法验证

3.2.1时间步长收敛分析 时间步长是影响数值模拟收敛性的一个重要因素，为验证数值方法收敛性和时间步长的关系，分别设定了dt=3.682×10-7s、dt=1.841×10-7s、dt=1.38×10-7s、dt=9.20×10-8s 4组工况进行计算，冲击速度选择61.4 m/s.计算结果如图5所示.从图5中可以看出，随着时间步长的改变，只有dt=3.682×10-7s的曲线与其他曲线有一定的差别，其他曲线重叠程度较高，表明时间步长在1.841×10-7s以下时，冲击力的计算值基本趋于收敛.

图5 不同时间步长下的冲击力计算结果Fig.5 Impact force calculation results at different time steps

图6所示为不同时间步长下的模拟现象，不同色彩为公式(16)计算得到的破坏程度.选取t=0.27 ms时刻.从图中可以看出，对不同的时间步长，同一时刻模拟的现象基本一致，除dt=1.841×10-7s的破坏程度略大一些外，其余无明显差别，可以认为当时间步长在1.841×10-7s以下时，模拟现象结果趋于收敛.因此，考虑到计算的效率，本文将时间步长选定为dt=1.841×10-7s.

图6 不同时间步长下的模拟现象Fig.6 Simulation phenomena at different time steps

3.2.2粒子间距收敛性分析 粒子间距是影响数值模拟收敛性的另一重要因素，粒子间距的选择会影响到结果的准确性与计算效率.本节选择冲击速度为61.4 m/s为计算工况，将柱形冰离散为D/10、D/15、D/20、D/25、D/30 5种粒子间距进行计算分析.图7所示为不同粒子间距下的冲击力计算值的比较曲线，由图可知，粒子间距对冲击力的计算结果具有较大的影响，冲击力峰值出现的时间和曲线变化过程甚至都会不同.从曲线中也可以看出，当粒子间距在D/25和D/30时，冲击力曲线基本一致，因此可以认为，当粒子间距设定为D/25～D/30时，冲击力计算结果趋于收敛.

图7 不同粒子间距下的冲击力计算结果Fig.7 Calculation results of impact force at different particle spacings

为了讨论不同粒子间距下柱状冰的破坏情况，图8给出了t=0.27 ms时不同粒子间距下柱状冰的破坏图.图8中的颜色变化说明了柱状冰的破损程度，蓝色表示冰柱未破坏，红色表示冰柱完全破坏.从图中可以看出，不同粒子间距条件下柱状冰破坏程度有所不同，粒子间距较大时破坏程度较低，随着粒子间距逐渐变小，当粒子间距在D/25和D/30时，破坏程度基本趋于一致，因此结合上文所述，可以认为粒子间距在D/25以上时，计算结果与模拟现象都趋于收敛，考虑到计算效率，本文的粒子间距选定为D/25.

图8 不同粒子间距下的模拟现象Fig.8 Simulation phenomenon at different particle spacings

3.2.3数值方法验证 本节的方法验证手段是将数值计算结果和现象与文献[27]中试验结果进行比较分析，冰柱的冲击速度为vice=91.4 m/s，海冰参数与文献[27]中测量结果一致.图9所示为冰柱冲击过程数值计算和试验的冲击载荷.从图9可以看出，冰柱冲击过程中，冲击载荷在极短的时间内快速上升，达到其最大破坏载荷后，冲击载荷逐渐下降为0.试验中，冲击载荷在0.12 ms时达到最大值，为3.4 kN，数值计算中，冲击载荷在0.1 ms时达到最大值3.54 kN.说明本文数值模型计算得到的冲击载荷与试验结果基本一致，故证明了本文采用状态型近场动力学方法模拟海冰的冲击破坏过程的准确性.卸载阶段数值计算结果与试验结果有一定的误差，主要是由于试验中冰材料整体发生了结构性破坏，因此卸载过程载荷变化并不单调.

图9 vice=91.4 m/s时柱状冰冲击载荷时间历程曲线Fig.9 Impact force of cylindrical ice at a time history of vice=91.4 m/s

图10所示为本文数值模拟的现象与试验现象[27]的对比，图中不同色彩表示公式(16)计算得到的破坏程度.在数值模拟的现象中，红色表示冰柱已完全破坏，随着颜色的变浅，表示冰柱破坏程度越低.冰柱自接触刚性板之后，自接触面向上逐渐发生粉碎性破坏，无明显裂纹产生，破坏后的部分呈现向四周发散的特点，该现象与试验现象具有极高的一致性，且能够模拟冰柱冲击过程的破坏程度变化与破坏传递.因此本文的数值方法可以高效地模拟出冰柱的粉碎性破坏过程.

图10 冰柱冲击过程数值模拟现象与试验现象对比Fig.10 Comparison of numerical simulation phenomenon and experimental phenomenon

3.3 计算工况参数敏感分析

3.3.1冲击速度 为研究不同冲击速度对冰柱冲击破坏过程的影响，本文在前文选定的计算参数基础上，改变了冲击速度进行了数值计算，冲击力计算结果如图11所示.

图11 不同冲击速度下的冲击力曲线Fig.11 Impact force at different impact speeds

可以看出，随着冲击速度的增加，冲击力的峰值有着明显的增加，从曲线形态上看，随着冲击速度的增加，曲线的波动情况也更加明显，而vice=30 m/s、vice=45 m/s的曲线是平滑的曲线，没有波动情况.同时，随着冲击速度的增加，冲击力卸载速度也越快，可能是因为速度越大冲击力越大，冰柱破坏程度越高，因此卸载速度也越快.冲击速度改变会导致冲击力峰值和卸载速度改变.

图12所示为不同冲击速度下的模拟现象，图中全部选取t=0.27 ms的时刻.可以看出，随着冲击速度的增大，冰柱的破坏程度更加明显，但是当冲击速度大于105 m/s时，破坏现象差别不大，主要是由于状态型近场动力学在模拟裂纹扩展方面具有优势，冰柱破坏主要是粉碎性破坏，没有明显裂纹，现象对比不明显.综上，在冰柱冲击问题中，状态型近场动力学在冲击速度45～100 m/s时，对速度这一参数有较高的敏感性.

图12 不同速度下的模拟现象Fig.12 Simulated phenomena at different speeds

3.3.2冰柱尺寸 冰柱尺寸是影响冰冲击破坏过程的另一个重要计算参数，本节改变冰柱形状，分析了不同尺寸冰柱的冲击破坏过程.其中冲击力计算结果如图13所示.

图13 不同大小的冰柱冲击力计算结果Fig.13 Calculation results of impact force of icicles of different sizes

从图13可以看出，随着冰柱尺寸增加，冲击力峰值明显增加，但曲线的形态相差不大，卸载速度也基本一致，因此冰柱形状只会导致冲击力峰值变化.

图14所示为不同大小的冰柱冲击现象，图中全部选取t=0.27 ms的时刻.可以看出，冰柱的尺寸变化后，其破坏程度、现象均无明显区别，因此可以认为状态型近场动力学方法对冰柱的直径的尺寸这一计算参数并不敏感.

图14 不同大小的冰柱冲击现象Fig.14 Impact of icicles of different sizes

3.4 冰材料参数敏感分析

3.4.1弹性模量 弹性模量是冰材料的重要属性，弹性模量与海冰的温度、盐度、卤水体积分数等有关系.Weeks等[28]等给出了海冰常见的弹性模量分布范围，其数值在1.7～9.1 GPa之间，本节改变了海冰的弹性模量，分别选取E=1.8，3，5，7，9 GPa进行了计算，其冲击力计算结果如图15所示.

图15 不同弹性模量的冲击力计算结果Fig.15 Calculation results of impact force at different elastic modulus values

可以看出，弹性模量会影响到冲击力的峰值与其卸载速度，尤其是卸载速度会随着弹性模量的增加明显增加，从曲线中可以明显看出，前半段在上方的曲线，在t=0.2 ms之后，都到了下方，分析原因是材料内部的应力传递速度只与材料的密度与弹性模量有关，弹性模量越大，应力传递速度越大，因此卸载的过程越快.

图16所示为不同弹性模量下的冲击效果，时间选为t=0.18 ms.从图中可以看出，随着弹性模量的增加，冰柱在同一时间的破坏程度逐渐增加，由此可以发现，弹性模量对于冰的冲击载荷峰值有极大的影响，对于破坏过程的传递速度影响较小.

图16 不同弹性模量下的冰柱冲击现象Fig.16 Icicle impact phenomenon at different elastic modulus values

3.4.2断裂韧度 断裂韧度是冰的重要力学属性，在本文采用的计算方法中，其决定了冰的失效判定标准.目前关于冰的断裂韧度的研究并没有给出统一的结论，尤其对于动态冲击问题，其断裂韧度的测量与研究更少，特别的在动态作用问题中，冰的断裂韧度并不是固定值，文献[29]指出，对于当年冰，其常见断裂韧度在115～250 (kPa·m)0.5.对此，本文改变了文中直接与断裂韧度相关的参数极限伸长率s0，来探究断裂韧度的影响，本文设定的s0对应的断裂韧度分别为108、150、198、252 (kPa·m)0.5.图17所示为不同断裂韧度条件下的冲击力计算结果.

图17 不同断裂韧度条件下的冲击力计算结果Fig.17 Calculation results of impact force under different fracture toughness conditions

从图中可以看出，改变断裂韧度对于冲击力的峰值没有影响，对于达到冲击载荷的时间影响也较小.断裂韧度越大，冲击力卸载越慢，分析原因可能是由于断裂韧度较大，冰柱破坏的程度相对较低，因此冲击力卸载速度较慢.此外，s0=0.001 5 的曲线和s0=0.003的曲线相差很小，说明当s0<0.003时，状态型近场动力学对断裂韧度这个参数不敏感.

图18所示为不同断裂韧度下冰柱撞击的现象，图中选取t=0.27 ms时刻.从图中可以看出，随着s0的逐渐增加，冰柱破坏程度逐渐降低.图18(a)和图18(b)现象差别不大，这也证明了冲击力分析的结果，状态型近场动力学在s0>0.003时，对s0的敏感度较高，更为适用.

图18 不同断裂韧度条件下的冰柱撞击现象Fig.18 Icicle impact phenomenon at different fracture toughness values

3.4.3泊松比 对于冰材料来说，其泊松比与盐度、温度、孔隙度等因素有关，不是固定值.在3-D键型近场动力学方法中，冰材料的泊松比只能设定为0.25，但是在状态型近场动力学方法中，对于材料的泊松比没有限制，因此可以对不同泊松比的冰柱冲击问题进行数值模拟研究.Timco等[29]在研究中指出，对于纵向冲击问题，冰材料常用的泊松比是0.33和0.42，因此本文设定了0.25、0.33、0.42这3种泊松比，同时对不同冲击速度下泊松比对计算结果的影响进行了计算.图19所示为3种冲击速度

图19 3种冲击速度下的不同泊松比计算结果Fig.19 Calculation results of different Poisson’s ratios at three impact speeds

下的计算结果.

从图19中可以看出，当冲击速度较低时，改变泊松比对计算结果的影响并不大，vice=30 m/s时的3条曲线基本贴合，没有发散现象.当冲击速度较高时，冲击力计算值的峰值随泊松比的增加而降低，同时卸载过程中的冲击力卸载速度也明显降低，以冲击速度vice=150 m/s的3条曲线为例，曲线前半段，泊松比υ=0.42的曲线在3条曲线中的最下方，但是在曲线后半段，υ=0.42的曲线处于3条曲线的最上方，说明泊松比增大，曲线下降的速度降低，即冲击力卸载速度越慢.此现象与材料中应力传递的横向弥散效应有关，横向弥散效应是由于材料泊松比不同导致的应力横向传播现象，泊松比越大，弥散效应越明显，应力传递速度越慢，这与计算结果一致.

图20所示为3种冲击速度下的不同泊松比的破坏现象，图中给出的现象均选取自t=0.27 ms时刻.从图20中可以看出，在冲击速度vice=30 m/s下的冰柱破坏现象上，破坏程度与模式受泊松比影响较小，破坏现象基本一致；在vice=60，90 m/s条件下，冰柱破坏程度随泊松比增加逐渐减小，对应图20曲线中其冲击力峰值逐渐减小，主要原因是随着泊松比的增加，物体在收到冲击力时其横向发生的位移会更大，导致其横截面积增大，冲击力沿横向传播更多，由于同一速度下冲击产生的能量是一定的，故纵向冲击力峰值降低，破坏程度也就越低，表现为随着泊松比的增加，破坏程度有一定的降低.

图20 3种冲击速度下的不同泊松比破坏效果Fig.20 Different Poisson’s ratio damage effects at three impact speeds

综上，泊松比对冰柱冲击破坏过程有一定的影响，这种影响表现为冲击速度一定时，泊松比越大，冰柱破坏程度越低，冲击力卸载速度越慢.这种影响在冲击速度较低时(vice=30 m/s)几乎不存在，随着冲击速度的增加，这种影响的表现逐渐明显，因此，状态型近场动力学方法在模拟高速的冰柱冲击问题时，较键型近场动力学方法更加符合真实情况.

4 结论

本文采用状态型近场动力学方法，对冰柱高速冲击问题进行了数值模拟，开展了数值方法验证与收敛性分析，而后分别针对计算工况参数和冰参数进行了敏感性分析，得出结论如下：

(1) 与试验结果相比，本文的数值模型能够较准确地模拟冰柱的破坏过程和冲击载荷，可以用于冰冲击破坏的研究.

(2) 针对本文的计算模型，当计算时间步长选取dt=1.841×10-7，粒子间距选取D/25及以下时，计算结果趋于收敛.

(3) 工况参数方面，冲击速度在45～100 m/s的区间内冰柱冲击破坏过程变化明显，可以明显观测冲击力与现象的不同.然而冲击破坏过程与冰柱的直径敏感相关性不高，冲击力变化趋势与冲击现象基本一致.

(4) 在冰参数方面，弹性模量对冰冲击破坏过程的破坏传递影响不大，即不同弹性模量的柱状冰可在同一时间内达到冲击载荷峰值，峰值的大小有明显的区别，弹性模量越大，峰值越高.海冰的断裂韧度对冲击载荷的峰值和冲击过程的现象基本没有影响.海冰的泊松比在高速(大于90 m/s)冲击时，对模拟计算结果与现象有影响，在低速冲击时影响不大.