截锥形弹体在液体介质中运动速度衰减规律分析*

孔祥韶，石 干，王旭阳,3，周红昌，吴卫国

（1. 武汉理工大学绿色智能江海直达船舶与邮轮游艇研究中心，湖北 武汉 430063；2. 武汉理工大学交通学院，湖北 武汉 430063；3. 西安航天动力研究所液体火箭发动机技术重点实验室，陕西 西安 710100）

大型水面舰艇通常采用液舱来防护反舰武器爆炸破片对内部结构和人员的威胁[1]。弹体在侵彻液舱过程中与液体的作用机理十分复杂，弹体在液体介质中的贯彻距离、速度衰减规律以及高速弹体引起的水锤效应，是广大学者长期关注的问题，对此已开展了大量的研究。Lee 等[2]在经典冲击压力模型的基础上，将弹体简化为点源，对高速弹体在液体中引起的空腔效应展开了研究，推导建立了空腔动力学方程。李营等[3]在考虑长径比的影响下，将阻力因数作为与初速度有关的常量来处理，拟合得到了阻力因数计算经验公式。沈晓乐等[4]考虑弹体的墩粗变形进行了实验研究，对方形弹体入水的阻力因数进行了修正。郭子涛[5]分析了弹头形状对弹体水弹道稳定性以及速度衰减规律的影响，提出了一种简化的弹体头型阻力因数估算公式。Zhao 等[6]在考虑雷诺数的基础上对球形弹入水阻力因数进行拟合，提出了球形弹入水阻力因数经验公式。Zhang 等[7]在Lee 等[2]研究工作的基础上针对长径比在0.1～0.5 之间的饼形弹进行研究，并对饼形弹入水阻力因数计算公式进行了修正。这些研究普遍简化了阻力因数的影响因素，与实际物理过程有较大的区别。

本文中，拟采用理论分析和数值模拟的方法，针对以上存在的问题，开展截锥形弹体垂直穿透液体介质时的速度衰减特性研究，重点分析不同弹体头形因数对其速度衰减规律的影响，以期建立阻力因数与弹体的瞬时速度之间的关系，并引入弹体形状参数的影响，建立截锥形弹体在液体中的速度衰减分析模型。

1 弹体侵彻液体介质速度衰减规律

当弹体在液体介质中高速运动时，弹体损失的动能主要转化为液体的动能及空泡的压力势能，而空泡的存在极大地减小了弹体与液体的接触面积，从而弹体运动所受的摩擦阻力减小，因此压差阻力是影响弹体阻力的主要因素。则根据牛顿第二定律建立弹体运动方程如下[8]：

2 数值计算

2.1 数值计算工况设计及模型

本节设计4 种不同头形因数的截锥形弹体，采用显式动力学非线性有限元程序Autodyn 建立弹体的仿真模型并开展系列工况的数值计算，分析头形因数对弹体垂直入水后速度衰减的影响规律。

图1 为截锥形弹体示意图，定义弹体头部圆台直径D2与主体最大直径D1的比值D2/D1为截锥形弹体头形因数 ψ 。各工况中弹体主体长度L1=30 mm ，锥体长度L2=10 mm ，主体直径D1=12 mm 保持不变，4 种弹体头部圆台直径分别为0、4、8、12 mm，对应头形因数 ψ1～ψ4分别为0、1/3、2/3、1。为了控制变量的数量，不同头形因数截锥形弹体的质量均调整为32.5 g，有限元模型如图2 所示。本次数值模拟采用二维轴对称建模方法，水域网格尺寸为400 mm×400 mm，采用边长为1 mm 的四边形欧拉单元进行离散，四周设置透射边界条件，模拟无限水域，弹体用拉格朗日网格离散。

数值模拟中，采用S h o c k 状态方程和Johnson-Cook 本构关系描述弹体材料的动态力学行为[12]：

图 1 截锥形弹体示意图Fig. 1 Schematic diagram of a truncated cone-shaped projectile

式中： σy为材料的动态屈服应力；A为材料的静态屈服应力；B为材料的硬化参数； εp为等效塑性应变；n为硬化指数；C为应变率参数；ε ˙ 为等效塑性应变率；ε ˙0为参考应变率，取 ε ˙0=1 s-1；T*=(T-Tr)/(Tm-Tr) ，T为温度，Tr为室温，Tm为材料熔化温度，m为温度指数。本文中弹体材料选用45 钢[13]，A为506 MPa，B为320 MPa，n为0.28，C为0.064，m为1.06，Tm为1 750 K。

图 2 不同头形因数截锥形弹体有限元模型Fig. 2 Finite element models for truncated cone-shaped projectiles with different head type coefficients

表 1 水的状态方程各项参数[14]Table 1 Parameters of equation of state for water[14]

本文中对每种不同头形因数的截锥形弹体各设置了5 种不同的入水初速度开展数值模拟，共包含20 种工况条件，如表2 所示。

表 2 数值计算工况Table 2 Numerical calculation conditions

2.2 有限元方法验证

采用本文的数值模拟方法对试验过程[5]开展计算，通过对比分析数值计算结果和试验数据来验证数值方法对高速弹体入水过程分析的可靠性和准确性。文献[5]中使用轻气炮发射圆柱形平头弹体，该弹体直径D=12.65 mm，长度L=25.4 mm，入水初速度为603 m/s，试验中使用高速摄影仪记录了弹体运动轨迹及空泡变化，测量了不同时刻的空泡尺寸以及弹体速度等数据。将数值模拟得到的高速弹体入水过程中与试验[5]中相应的数据进行对比，弹体入水后空泡尺寸对比如图3 所示，可以发现在t=0.222 ms和t=0.440 ms 时刻相同位置处的空泡直径数值计算结果与试验测量数据吻合较好

图 3 针对弹体直径为12.65 mm，长度为25.4 mm，入水初速度为603 m/s 时，试验与数值计算的空泡尺寸Fig. 3 Comparison of cavitation sizes obtained experimentally and numerically for the projectile with the diameter of 12.65 cm and the length of 25.4 cm, water entering at 603 m/s

此外，进一步开展尺寸D=12.65 mm，L=25.4 mm 的平头圆柱形弹体入水速度603 m/s 及397m/s 工况；D=12.65 mm，L=38.1 mm 的平头圆柱形弹体入水速度为498 m/s 及414 m/s 工况的数值计算，将弹体位移随时间的变化以及弹体速度衰减的数值计算结果与相应的试验测量数据进行对比，如图4～5 所示。可以看出，本文数值模拟方法在计算弹体的运动特性方面具有较好的精度，可为本文的进一步的研究工作开展提供了准确可靠的分析手段。

图 4 长度为25.4 mm 的弹体在两种入水速度工况下位移和速度的变化Fig. 4 Changes of displacement and velocity of the projectile with the length of 25.4 mm at two initial velocities of water entry

图 5 长度为38.1 mm 的弹体在两种入水速度工况下位移和速度的变化Fig. 5 Changes of displacement and velocity of the projectile with the length of 38.1 mm at two initial velocities of water entry

2.3 数值模拟结果及分析

采用验证后的数值模拟方法开展表2 中20 种工况的数值计算，将头形因数不同的截锥形弹体在不同入水速度下的速度衰减与入水距离的关系绘制于图6，通过对比分析图中数据可以发现：

（1）当截锥形弹体的头形因数不变时，弹体入水后速度变化曲线的斜率随着入水初始速度的升高而逐渐增大，说明入水速度越高，速度衰减越快；

（2）弹体入水初速度不变、头形因数增大时，弹体剩余速度并非单纯随头形因数增大而降低，这一特性在v0=1 400 m/s 时表现最明显。4 种头形因数的弹体均以1 400 m/s 的初速度垂直入水时，运动400 mm后弹体剩余速度与初速度的比值vb/v0分别为0.699 0、0.771 0、0.709 9、0.512 7。说明入水初速度不变时，随着弹体头形因数的增大，vb/v0呈现先增大后减小的趋势，在 ψ =1/3 时达到最大，为0.771。

图 6 不同头形因数的截锥形弹体在不同入水速度工况下速度随入水距离的衰减Fig. 6 Velocity attenuation of truncated cone-shaped projectiles with different head coefficients with water-entry distance at different initial velocities of water entry

3 速度衰减分析模型

3.1 模型建立

根据式（4）可以计算得出截锥形弹体垂直入水后的瞬时阻力因数Cd，并将其作为纵坐标，将弹体的瞬时速度vb与入水初速度v0的无量纲比值作为横坐标，绘制Cd- (vb/v0) 关系图。由于篇幅有限，本文中仅绘出ψ=0 的弹体Cd- (vb/v0) 关系图，如图7 所示。从图7 可以发现，阻力因数Cd随着弹体速度的衰减呈振荡趋势，且当截锥形弹体速度较低时，随着侵彻过程的进行，速度衰减率出现较小幅度的上升。究其原因是因为随着速度的降低，弹体诱导空泡的尺寸也在不断减弱，弹液间的接触面积增大，摩擦阻力增大，导致弹体的速度衰减率和阻力因数Cd小幅增大。

图 7 不同入水速度下， ψ =0 的截锥形弹体阻力因数与 v b/v0 的关系Fig. 7 Resistance factor varying with v b/v0 for the truncated cone-shaped projectile with ψ=0 at different water-entry velocities

从完整的侵彻过程来看，阻力因数Cd变化过程较复杂，阻力因数随着弹体速度的降低，呈现如上述图像所示的震荡变化趋势。在以往的研究方法中将其作为常量来处理，不具备广泛的适用性。从简化计算的角度出发，本文中采用线性关系进行拟合：

式中：a0和a1为未知参数。

针对图7 中的数据，采用最小二乘法进行拟合，拟合的结果如图7 中红线所示。同时，对头形因数ψ=1/3 ， ψ =2/3 ， ψ =1 的 截锥形弹体数据按照同样方式进行拟合，确定的参数a0和a1列于表3。

表 3 不同工况下的参数 a0 和 a1 的数值Table 3 Values of parameters a0 and a1 under different working conditions

水中声速取1 400 m/s，用最小二乘法对参数a0和a1进行拟合：

式中：a01、a02、a03、a04、a05、a11、a12、a13、a14、a15为待定参数，结果如图8 所示。

由拟合结果来看，四阶方程可以较好的拟合a0、a1关于无量纲数M的关系，可以发现同一头形因数下参数a0及a1的高度相关，但变化趋势相反，近似关于某一水平线y=n1对称。随着头形因数的增大，参数a0及a1的变化逐渐趋于复杂，这是因为当头形因数增大时，弹体的垂直迎流面积增大，侵彻时镦粗变形效应变得明显，导致阻力因数的变化趋于复杂。参数a0、a1的拟合结果列于表4。

图 8 对于头形因数不同的截锥形弹体，参数a0 和a1 与马赫数 Ma 的关系Fig. 8 Changes of parameters a0 and a1 with Mach number Ma for the truncated cone-shaped projectiles with different head shape factors

表 4 参数 a0 、 a1 的拟合结果Table 4 Fitting results of parameters a0 and a1

图 9 a 系列参数与头形因数 ψ 的关系Fig. 9 Series parameters of a varyied with head shape factor ψ

表 5 ai j 系列参数的拟合结果Table 5 Fitting results of ai j series parameters

综上所述，当头形因数 ψ 不同的截锥形弹体垂直侵彻入水时可以通过如下公式计算实时阻力因数Cd：

3.2 算例分析及验证

通过2.2 节中已验证的计算方法开展系列模拟计算，通过对公式（18）的编程求解，得到弹体的速度衰减。将数值计算结果与上述速度衰减模型计算结果进行对比，验证其可靠性。

验证计算中，设置v0=400 m/s 速度下1/6、1/2、5/6 等3 种头形因数工况以及头形因数为 ψ =1/6 的弹体在100～600 m/s 区间内6 种入水速度工况进行验证计算，结果如图10～11 所示。

图 10 入水速度为400 m/s，头形因数不同的弹丸，速度随入水距离衰减的模拟结果与公式计算结果的比较Fig. 10 Comparison of velocity attenuation with distance between numerical simulation and formula calculation for the projectiles with different head shape factors at the initial water-entry velocity of 400 m/s

由计算可以发现，对于弹径相同、入水初速相同而头形因数不同的弹体，其速度衰减规律有着显著差异。本文的计算分析模型中考虑了因头形因数不同而引起的阻力因数变化对弹速衰减的影响，与数值计算结果吻合较好。当弹体入水初速为400 m/s，头形因数为1/2 时，数值计算数据与本文公式计算结果之间的误差约为6.8%，其他两种头形因数工况下误差更小；在弹体头形因数为1/6、不同速度工况中，公式计算与数值模拟结果之间的最大误差约为8.71%。

图 11 头形因数为1/6，入水速度不同的弹丸，速度随入水距离衰减的模拟结果与公式计算结果的比较Fig. 11 Comparison of velocity attenuation with distance between numerical simulation and formula calculation for the projectile with the head shape factor of 1/6 at different initial water-entry velocities

从上述对比结果可以发现，本文的速度衰减模型可对不同头形因数的截锥形弹体垂直侵彻液体介质过程中的速度衰减进行快速预报，弥补了以往研究中将阻力因数Cd当作常量计算的不足，更符合弹体垂直侵彻液体介质的实际物理过程。

4 结 论

本文中提出了与弹体形状和无量纲速度相关的阻力因数的计算方法，通过建立运动方程得到了截锥形弹体垂直侵彻液体介质速度衰减的计算分析模型。并通过算例分析，验证了本文计算方法的可靠性。得到以下结论：

（1）高速弹体在液体介质中运动方程的阻力因数与弹体形状和无量纲速度有关；该因数对弹体高速运动过程中与液体的接触面积进行了修正，将弹体头形因数和速度因素考虑在内，更接近实际的物理过程。

（2）将阻力因数引入运动方程后得到的计算结果与数值模拟数据吻合较好, 说明本文提出的分析模型可对不同头形因数的截锥形弹体垂直侵彻液体介质过程中的速度衰减进行快速预报。

（3）由于弹体侵彻液体时的作用机理非常复杂，且弹体高速侵彻时会发生弹体镦粗和变形，同时流体具有可压缩性，对弹速衰减影响较大。本文的分析中未考虑其影响，在后续工作中将进一步的深入研究。