前抛体对弹体入水载荷影响数值模拟研究*

2023-04-18魏继锋吉耿杰冯宇浪

爆炸与冲击 2023年4期

肖睿，魏继锋，吉耿杰，冯宇浪

（北京理工大学爆炸科学与技术国家重点实验室，北京 100081）

结构体入水过程涉及固-液-气三相耦合作用，具有高度非定常性和非线性，在救生艇下水、航天器水上回收、水上飞机降落、空投鱼雷入水等诸多领域都有着广泛应用。当结构体穿过气-液界面时，由于环境介质的密度突变，其头部会受到剧烈的冲击，结构极易损坏失效[1]。因此，研究结构体入水载荷特性及其降载技术具有重要的学术意义和工程应用价值。

最早开展的入水冲击载荷实验研究的是Bottomley[2]进行的V 型浮舟模型入水实验，而最早的入水冲击载荷理论研究则是von Karman[3]提出的楔形体和平板入水冲击理论。Wagner[4]在von Karman 入水冲击模型[3]的基础上引入水面变形因子，提出了入水冲击问题的自相似解法。Eroshin 等[5]提出钝体结构入水时，若忽略液体可压性将导致入水冲击载荷被高估，因此应考虑液体的可压性。随着对入水载荷问题认识日益深刻，入水载荷描述愈加复杂，理论解析愈发困难。数值模拟程序的研发和数值计算能力的提升，促进了数值模拟技术在入水冲击载荷研究中的广泛应用。

在入水冲击载荷降载方面，常用的降载方法有调整弹体结构、增加缓冲罩壳、加装空化器等。改善弹体头部形状可有效减小入水冲击载荷，斜锥形头部弹体受到的载荷比半球形、锥形、和球台形弹体受到的载荷更小[6]，椭球头型受到的最大冲击力与速度和头部长短轴比有关[7]，头部半球角越小入水冲击载荷越小[8]，独特的海鸟头型可有效减小入水冲击载荷[9]。缓冲罩壳在入水时的碎裂和缓冲组件的变形吸能可降低入水冲击载荷，最早由Howard[10]提出，之后的学者又研究了不同材料、结构的罩壳及其缓冲组件的降载效果[11-18]。弹体头部加装喷气空化器在入水前使水面产生空泡来降低冲击载荷，空化器通气量越大、弹体速度越高，降载效果越好[19-21]。Rabbi 等[22]提出串联入水降载思路，在进行的两小球前后入水试验中发现，前抛球可使后球的入水冲击载荷下降78%；Lyu 等[23]进行了速度为15～40 m/s 的两个相同小球串联入水试验，发现后入水小球的降载率随时间间隔的增大而减小。

目前，对串联入水体降载问题的研究刚刚起步，仅仅分析了两个小球垂直入水的情形，对于其他结构体、不同入水速度、不同入水姿态的入水特性和降载规律尚未见报道。本文据此开展平头圆柱弹体串联入水冲击数值模拟研究，采用S-ALE (structured arbitrary Lagrange-Euler)算法和罚函数流固耦合方法，深入探讨前抛体入水角度、时间间隔、入水速度等参量对入水空泡和入水冲击载荷的影响，探索串联降载机制和降载影响规律，以期为弹体降载新技术的运用提供重要的理论和数据支撑。

1 含前抛体弹体入水过程及理论分析

根据von Karman 的平板入水冲击模型[3]，平头圆柱形弹体入水时受到的冲击载荷为

式中：p为平板入水所受压力；ρw为水的密度；vp为弹体与水自由面间的相对速度；cw为水中声速。

弹体受到的加速度为

式中：a为弹体加速度；F为弹体受到的冲击力；A0为平头圆柱与水的初始接触面积；m为弹体质量。由式(2)知，决定弹体冲击载荷的因素有水的密度与声速、弹体与水自由面间的相对速度和弹体入水接触面积，因此可通过减小弹体与水面之间相对速度、降低水体密度以及减小初始接触面积的方式来降低弹体入水冲击载荷。

含前抛体的弹体入水时，前抛体改变水体的状态，使原本静止的水体具有了一定的速度，减小了弹体与水体之间的相对速度，降低了主弹体撞击水体的砰击载荷。水体在前抛体的作用下会经历空泡形成、空泡闭合（包括面闭合、浅闭合、深闭合）、射流形成等阶段。与之相对应的，主弹体可能在这之中的任意阶段入水。当然，前抛体和主弹体顺序入水时可能出现相撞的现象。如图1 所示，将含前抛体弹体入水过程中会出现的情况归结为以下四种：Ⅰ，撞击前抛体；Ⅱ, 撞击空泡壁面；Ⅲ，撞击飞溅水体；Ⅳ，撞击射流。

图1 含前抛体弹体入水情况示意图Fig.1 Projectile with front projectile entering the water

2 数值模拟算法及设置

2.1 S-ALE 与流固耦合算法

S-ALE (structured arbitraty Lagrange-Euler)算法于2015 年提出，在LS-DYNA 中得到运用[24]，与耦合了Lagrange 算法和Euler 算法的ALE 方法接近[25]：计算时先执行Lagrange 时间步，单元网格随材料产生变形；然后保持变形后的边界条件，对内部单元网格重新划分，网格拓扑关系保持不变；最后将变形网格中的单元变量和节点速度矢量输送到重分后的新网格中[26]。S-ALE 与ALE 方法的区别在于，前者的流场采用正交的结构化网格，只需在设置中输入i、j、k三个方向的节点数量及坐标即可由计算机自动生成网格，便于设置且更容易收敛，并行计算效率更高。

S-ALE 算法的流场控制方程为

式中：ρ 为流体密度；v为流体速度；w为流体与网格间的相对速度，w=v-u，u为网格移动速度；σi为应力张量，σij=-pfδij+µ(vij+vji)，其中pf为流场中的压力，µ为动力黏性系数，δij为Kronecker 函数，当i=j时为1，否则为0；b为单位体积力；e为能量；下标i、j表示不同方向的欧拉坐标。

网格控制方程组为

式中：Xi为拉格朗日坐标。

计算入水问题时气-液界面捕捉是研究的重要内容。S-ALE 采用的界面捕捉方法是根据网格单元中不同物质的体积分数梯度方向和大小来确定切割面或分界面的位置，其计算方式如图2 所示，图中η 为物质的体积分数，上标“-”表示节点的平均体积分数；V为单元体积。方向n由公式(7)确定。

图2 S-ALE 的界面捕捉法Fig.2 Surface capturing method of S-ALE

流固耦合中固体对流体的作用通过添加在流场控制方程中的“源项”来实现，流场对固体结构的作用通过结构附近的流体质点速度的插值实现[27]。其控制方程组为

式中：F(x,t)为边界传递给流场的体积力；x为流固耦合边界的位移；f为流固耦合边界产生的单位力；X(s,t)为流体与固体耦合边界的位移；s为流体与固体的边界坐标；δ 为Dirac delta 函数，该函数除0 坐标以外函数值为0，在其整个定义域上的积分为1。S-ALE 的流固耦合算法可以有效解决ALE“罚函数”耦合方法产生的“渗漏”问题，获得更准确的计算结果。

2.2 计算模型

主弹体和前抛体均为平头圆柱，主弹体尺寸为 ∅ 6 mm×24 mm，前抛体尺寸为 ∅ 2 mm×2 mm。流体计算域宽200 mm，水深300 mm，对弹体经过区域的网格进行加密。弹体采用Lagrange 网格划分，流体域用ALE 网格剖分。由于计算模型的对称性，采用1/2 模型进行计算，以主弹体入水点为原点O，水平向右方向为x轴正向，竖直向下为z轴正向，将xOz平面设为对称平面。弹体模型与流体域网格如图3 所示。进行计算时，首先设置主弹体速度为0，前抛体以预定速度入水，计算至所需时间间隔时设置主弹体速度，通过LS-DYNA 的重启动计算主弹体入水载荷。

图3 弹体模型与流体域网格Fig.3 Projectile models and fluid domain mesh

2.3 材料模型

计算模型包括流体（水和空气）、主弹体和前抛体四种材料。

水采用NULL 材料模型和Grüneisen 状态方程：

式中：pw为水的压力；µw为水的体积变化率，µw= ρw/ρw0-1；γ0为Grüneisen 系数；a1为γ0的一阶修正；S1～S3为us-up曲线斜率的无量纲系数，us为冲击波速度，up为水的质点速度；Ew为水的体积内能。表1 为水的Grüneisen 状态方程参数值。

表1 水的Grüneisen 状态方程参数Table 1 Grüneisen equation of state parameters of water

空气采用NULL 材料模型和线性多项式状态方程：

式中：pa为空气压力；µa为空气的体积变化率；C0～C6为状态方程系数；Ea为空气的体积内能。用Gamma 定律描述气体状态方程时

式中：γ 为绝热常数，空气为1.4。表2 为空气的多项式状态方程参数值，其中µa0为µa的初始值。

表2 空气的多项式状态方程参数Table 2 Polynomial equation of state parameters of air

主弹体和前抛体由钢材料制成，采用刚体材料模型，密度为7 800 kg/m3，杨氏模量为2.1×1011Pa 和泊松比为0.3。

3 无前抛体时主弹体的入水载荷

当主弹体以106.8 m/s 的速度垂直入水时，将弹体头部接触水面的时刻记为时间零点，图4 和图5 给出了主弹体入水4 ms 内的空泡形态及入水深度，与文献[28]的试验结果吻合。虽然文献[28]中并未测得弹体入水时的冲击载荷，但是通过弹体的运动状态可以反映其动力学特征,由此可以判断数值模拟结果可信。图6 为主弹体无前抛体入水过程中的受载加速度数值模拟曲线，主弹体入水瞬间的砰击载荷高达806.3 km/s2，与式(1)和式(2)的理论计算结果844.3 km/s2的误差为4.507%。

图4 弹体入水运动状态和空泡形态Fig.4 Projectile entering water motion state and cavitation shape

图5 弹体入水深度(z)Fig.5 Projectile water entry depth (z)

入水过程中，弹体头部受到水面冲击的同时，附近的水介质也受到冲击并产生一道压缩波，以声速在水中传播。当压缩波到达水面时在自由表面发生卸载并产生稀疏波，使主弹体受到的冲击载荷卸载，因此从图6 中可以看到在极高的初始冲击峰值之后主弹体的冲击加速度快速降低至零。在极短的一段时间后弹体头部压力恢复为水中动压并趋于稳定，弹体在水中开始稳定航行，弹体速度在流体阻力的作用下逐渐衰减。图6 中加速度值的波动为流固耦合时耦合节点变化产生的数值波动。

图6 主弹体的加速度(a)Fig.6 The accelerationof the main projectile (a)

4 有前抛体时主弹体的入水载荷

4.1 前抛体角度变化对主弹体载荷的影响

4.1.1 前抛体角度变化时的空泡演化特点

将前抛体的初始入水角定义为圆柱轴线与水平面夹角，如图7 所示，记为θ；主弹体的载荷变化用降载率α 表征，降载率定义为无前抛体时主弹体最大加速度与含前抛体入水时主弹体最大加速度之差与无前抛体时主弹体最大加速度的比值，即

图7 前抛体初始入水角示意Fig.7 Schematic of the initial water entry angle of front body

图8 展示了θ 为70°～90°时前抛体以100 m/s 入水产生的空泡在0～20 ms 内的演化过程。前抛体入水后与周围水体产生了剧烈的动量交换，水体被排开形成逐渐扩大的空泡，水面升高并形成了飞溅水体。随着空泡的扩大，为补充空泡内的气体，大气内的空气从空泡开口处流入空泡，开口处因为Bernoulli 效应而形成低压区，飞溅水体由于压力差的作用向中心聚拢闭合，并形成了上下两个方向相反的射流。其中，前抛体以θ=70°入水时未产生明显的射流，而是形成了一片较大的水幕，水幕形态不规则且散乱。

图8 前抛体以不同角度入水形成的空泡演化过程Fig.8 The evolution process of the cavitation formed by the front body entering the water at different angles

此时引入一个入水时间间隔的无量纲参数τ，其表达式为

式中：Δt为主弹体与前抛体入水的时间间隔，tp为空泡的闭合时间。由图8 可得当θ=90°时tp=13 ms；当θ=80°时tp=13 ms；当θ=70°时tp=10 ms。

4.1.2 前抛体90°入水时的主弹体载荷

图9 为前抛体和主弹体以100 m/s 垂直入水时的加速度载荷峰值随时间间隔的变化曲线。图10 给出了三个不同τ 时的主弹体载荷曲线，可见在不同τ 值范围内主弹体载荷具有不同的特性。当τ＜0.77时，主弹体直接进入至空泡内部，主弹体不与空泡壁面发生任何接触，与前方受到水阻力而减速的前抛体相撞，出现情况Ⅰ。此时，主弹体受到的冲击加速度远远超过无前抛体时的冲击加速度。图11 为τ=0.08 时主弹体与前抛体的速度曲线，可以看到主弹体与前抛体相撞之后主弹体减速，前抛体加速，二者再次分离。前抛体受到水阻力快速减速的同时排开周围的水体，使主弹体受到的阻力减小，主弹体的加速度较小，二者的速度差使它们再次碰撞。以上过程多次循环后二者的速度差逐渐减小，最终接触后不再分离，由主弹体推动前抛体一起运动。

图9 前抛体垂直入水时的主弹体载荷峰值曲线Fig.9 Load peak curves of the main projectile when the front body enters water vertically

图10 不同τ 时前抛体垂直入水时的主弹体载荷Fig.10 Load of the main projectile when the front body vertically enters the water while τ is different

图11 当τ=0.08 时主弹体与前抛体的速度曲线Fig.11 Velocity curves of main projectile and front body when τ=0.08

当0.77≤τ＜1 时，主弹体冲击加速度曲线中出现两组峰值。主弹体在进入空泡前先撞击水面以上的飞溅水体，主弹体受到第一次冲击，出现情况Ⅲ；之后主弹体与空泡内的前抛体相撞，出现情况Ⅰ。对比两次载荷峰值，可见情况Ⅰ产生的峰值远远大于情况Ⅲ。

当τ≥1 时，在主弹体和前抛体的撞击载荷峰值前存在一次较小的载荷，这个载荷由空泡壁面产生。主弹体入水过程中先撞击射流，然后撞击空泡壁面，最后撞击前抛体，此过程中依次出现情况Ⅳ、情况Ⅱ、情况Ⅰ。对比三次载荷，峰值由大到小依次为情况Ⅰ、情况Ⅳ、情况Ⅱ。

前抛体90°入水时，虽然τ 不同时主弹体受到的载荷特性不同，但是主弹体和前抛体的碰撞都是主弹体所受冲击载荷的最主要原因。在目前的计算时间内主弹体的入水载荷随时间间隔的增大而减小。

4.1.3 前抛体80°入水时的主弹体载荷

图12 为前抛体80°、主弹体垂直，两弹体以100 m/s 入水时主弹体的加速度载荷峰值随时间间隔的变化曲线。图13 给出了三个不同τ 时的主弹体载荷曲线，在不同τ 值范围内主弹体载荷具有不同的特性。图14 给出了τ=1 时主弹体入水过程中与水体的砰击过程。当τ＜0.7 时，主弹体直接进入空泡撞击空泡壁面，出现情况Ⅱ；当0.7≤τ＜1 时，冲击加速度曲线出现两组峰值，分别由飞溅水体和空泡壁面产生，出现情况Ⅲ、Ⅱ；当τ≥1 时，两组峰值分别由射流和空泡壁面产生，出现情况Ⅳ、Ⅱ。

图12 前抛体80°入水时的主弹体载荷峰值曲线Fig.12 Load peak curves of the main projectile when the front body enters the water at 80°

图13 不同τ 时前抛体80°入水时的主弹体载荷Fig.13 Load of main projectile when the front body enters the water at 80° while τ is different

图14 前抛体80°入水，τ=1 时的两次砰击过程Fig.14 Two slamming processes when the front body enters the water at 80 ° and τ=1

前抛体80°入水时，对主弹体的载荷由空泡壁面和飞溅水体或射流产生，其中射流产生的载荷最大，其次为空泡壁面，飞溅水体产生的载荷最小，即产生的载荷峰值由大到小依次为情况Ⅳ、情况Ⅱ、情况Ⅲ。在射流生成前主弹体入水载荷随时间间隔的增大而减小。

4.1.4 前抛体70°入水时的主弹体载荷

图15 为前抛体以70°、主弹体垂直以100 m/s 入水时主弹体冲击加速度载荷随时间间隔的变化曲线。图16 给出了两个不同τ 时的主弹体载荷曲线，在不同τ 值范围内主弹体载荷具有不同的特性。图17给出了τ=0.8 时主弹体入水过程中与水体的砰击过程。当τ＜0.6 时，主弹体不与飞溅水体相撞，只产生一次载荷，为情况Ⅱ；当τ≥0.6 时，加速度曲线中出现两组峰值，分别由飞溅水体和空泡壁面产生，出现情况Ⅲ、Ⅱ。

图15 前抛体70°入水时的主弹体载荷峰值曲线Fig.15 The load peak curves of the main projectile when the front body enters the water at 70°

图16 不同τ 时前抛体70°入水时的主弹体载荷Fig.16 The load of the main projectile when the front body enters the water at 70° while τ is different

图17 前抛体70°入水，τ=0.8 时的两次砰击过程Fig.17 Two slamming processes when the front body enters the water at 70° and τ=0.8

前抛体以70°入水时形成的空泡和飞溅水体的形态都不规律，因此主弹体的载荷随时间的变化也不规律。情况Ⅱ产生的载荷峰值大于情况Ⅲ。

4.1.5 总结

如图18 所示，三种不同的前抛体入水角度分别产生了不同的入水情况，每种入水情况下主弹体的入水载荷特性各不相同，最终的降载效果也不同。

图18 前抛体入水角度变化时主弹体载荷峰值随τ 的变化曲线Fig.18 Curves of the main projectile load peak with τ when the water-entry angle of the front body changes

前抛体90°入水产生了情况Ⅰ～Ⅳ。入水情况Ⅰ会产生较大的主弹体载荷，而且该载荷有可能远超无前抛体时的载荷，应尽量避免。载荷随τ 的增大而减小。当τ=1.46 时主弹体入水载荷最小，为575.8 km/s2，降载率为28.588%，降载效果并不理想。

前抛体80°入水产生了情况Ⅱ～Ⅳ。在形成射流前主弹体入水载荷由撞击空泡壁面（情况Ⅱ）产生，载荷随时间间隔的增大而减小。射流形成后主弹体载荷主要由撞击射流（情况Ⅲ）产生，由于射流与主弹体速度相反，它们之间的相对速度较大，产生了较大的冲击，但是射流的密度较小，冲击产生的主弹体载荷依然小于无前抛体时的载荷。射流形成前主弹体的降载率在80%～90% 之间，射流形成后在50%～60%之间，其中τ =0.85 时主弹体入水载荷最小，为80.9 km/s2，降载率为89.967%。

前抛体70°入水时产生情况Ⅱ、Ⅲ，这一条件下的空泡、飞溅水体的形态都不规律，主弹体入水载荷以空泡壁面载荷为主，随τ 的变化不规律。总体降载率在70%～85%之间，其中τ=0.9 时主弹体入水载荷最小，为131.8 km/s2，降载率为83.670%。

从以上分析可知，前抛体以80°入水、τ 处于0.8～0.9 时降载效果最好。此时主弹体先经过结构松散的飞溅水体经历第一次减速，然后再砰击扩展中的空泡壁面，主弹体和水体自由表面间的相对速度降低，主弹体受到的冲击减小。

4.2 前抛体尺寸变化对主弹体载荷的影响

根据前述分析，选取前抛体以80°、100 m/s 入水，τ=0.8 作为入水条件，研究直径为2～6 mm，长径比为1 的前抛体对主弹体入水载荷的影响。主弹体依然以100 m/s 垂直入水。图19 给出了不同前抛体尺寸下的主弹体载荷曲线，可见随着前抛体尺寸增大，其入水动能增大，产生的空泡壁面扩展速度更快，主弹体与空泡壁面间的相对速度减小，降低了主弹体的入水载荷。当前抛体直径为6 mm 时，主弹体入水冲击的降载率为93.905%。此外，由于空泡扩展速度更快，开口处的Bernoulli 效应更加明显，空泡闭合时间tp也随之减小。

图19 τ=0.8 时前抛体不同尺寸下的主弹体入水载荷曲线Fig.19 Water-entry load curve of the main projectile under different sizes of the front body while τ=0.8

4.3 速度变化对主弹体载荷的影响

为研究在不同入水速度对降载效果的影响，记主弹体和前抛体的初速度为v0，设置v0=40～100 m/s，τ=0.8。前抛体直径为6 mm，以80°入水；主弹体垂直入水。前抛体入水速度越大，空泡壁面的扩展速度越快，主弹体和空泡壁面的相对速度越小，如图20 所示主弹体载荷随入水初速度的增加而增大。如图21所示，随入水初速度的增加，主弹体的降载率增大，所有计算工况的降载率均在90%以上。