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

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

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

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

式中：M为弹体质量，${v_{\rm{b}}}$为弹体在$\tau$时刻瞬时速度，${\rho _{\rm{w}}}$为液体密度，${A_0}$为弹体在运动方向上的投影面积，${C_{\rm{d}}}$为阻力因数。

对式（1）进行变换可得：

可得阻力因数${C_{\rm{d}}}$与瞬时速度${v_{\rm{b}}}$的关系：

式中：$v_{\rm{b}}^{(i)}$为i时刻弹丸的瞬时速度，$v_{\rm{b}}^{(i + 1)}$为i+1时刻弹丸的瞬时速度，$x_{\rm{b}}^{(i)}$为i时刻弹丸的瞬时距离，$x_{\rm{b}}^{(i + 1)}$为i+1时刻弹丸的瞬时距离。

目前对于阻力因数${C_{\rm{d}}}$，学者们提出了不同的取值方式。

沈晓乐等^[4]认为当弹体形状为立方体时${C_{\rm{d}}} = 0.15$。

Garrett等^[9]认为阻力因数与空化数$\sigma$满足如下关系：

式中：${C_0}$的取值为0.82～0.83。空化数：

式中：${p_\infty }$为水压力；${p_{\rm{c}}}$为空泡压力，一般取水的饱和蒸气压（$2.3\;{\rm{ kPa}}$）。

孔祥韶等^[10]根据Lecysyn等^[11]给出的有关雷诺数Re的计算公式：

得出，当弹体的运动速度在200～2 000 m/s之间时，${C_{\rm{d}}}$取值在0.490 8～0.491 7之间，并提出在简化计算时${C_{\rm{d}}}$可取定值$0.491$。

Lee等^[2]对于球形弹体入水问题，采用分段的阻力因数：

式中：$Ma = {v_0}/{c_{\rm{w}}}$，${v_0}$为弹体入水初速度，${c_{\rm{w}}}$为水中声速。

Zhao等^[6]在考虑雷诺数的基础上对球形弹入水阻力因数进行拟合：

李营等^[3]在Lee等^[2]工作的基础上针对圆柱形平头弹开展了研究，并考虑长径比参数的影响，对阻力因数进行了修正。

Zhang等^[7]在Lee等^[2]工作的基础上针对长径比在0.1～0.5之间的饼形弹，考虑入水角度的变化，提出了弹体侵彻水介质过程中速度衰减的修正公式。

但是在以上的研究中，都是将阻力因数${C_{\rm{d}}}$作为只与初始入水速度${v_0}$有关的常量来处理，弹体的受力面是按最大截面来考虑，没有考虑到弹体头部形状的差异，与实际物理过程不符，本文中针对这一问题开展研究。

2.   数值计算

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

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

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

图  1  截锥形弹体示意图

Figure  1.  Schematic diagram of a truncated cone-shaped projectile

下载: 全尺寸图片 幻灯片

图  2  不同头形因数截锥形弹体有限元模型

Figure  2.  Finite element models for truncated cone-shaped projectiles with different head type coefficients

下载: 全尺寸图片 幻灯片

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

式中：${\sigma _{\rm{y}}}$为材料的动态屈服应力；$A$为材料的静态屈服应力；$B$为材料的硬化参数；${\varepsilon _{\rm{p}}}$为等效塑性应变；$n$为硬化指数；$C$为应变率参数；$\dot \varepsilon$为等效塑性应变率；${\dot \varepsilon _0}$为参考应变率，取${\dot \varepsilon _0} = 1\;{{\rm{s}}^{ - 1}}$；${T^ * } = (T - {T_{\rm{r}}})/({T_{\rm{m}}} - {T_{\rm{r}}})$，$T$为温度，${T_{\rm{r}}}$为室温，${T_{\rm{m}}}$为材料熔化温度，$m$为温度指数。本文中弹体材料选用45钢^[13]，$A$为506 MPa，$B$为320 MPa，n为0.28，C为0.064，m为1.06，${T_{\rm{m}}}$为1 750 K。

液体水介质采用多项式状态方程，在水介质受压时（$\mu {\text{＞}} 0$）：

在水介质受拉时（$\mu {\text{＜}}0$）：

式中：${A_1}$、${A_2}$、${A_3}$、${B_0}$、${B_1}$、${T_1}$、${T_2}$为材料常数；$\mu = \rho /{\rho _0} - 1$，$\rho$为水的瞬时密度，${\rho _0}$为水的参考密度；$E = (\rho gh + {p_0})\rho {B_0}$为质量内能，g为重力加速度，h为弹体在水中的深度，${p_0}$为大气压力。具体参数如表1所示^[14]。

表  1  水的状态方程各项参数^[14]

Table  1.  Parameters of equation of state for water^[14]

${A_1}{\rm{/GPa}}$	${A_{\rm{2}}}{\rm{/GPa}}$	${A_3}{\rm{/GPa}}$	${B_0}$	${B_1}$	${T_1}{\rm{/GPa}}$	${T_2}{\rm{/GPa}}$	${\rho }_{0}\rm{/(kg}·{\rm{m}}^{\rm{-3}}\rm{)}$
2.2	9.54	14.57	0.28	0.28	2.2	0	1 000

下载: 导出CSV 

| 显示表格

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

表  2  数值计算工况

Table  2.  Numerical calculation conditions

工况	$\psi$	v0/(m·s−1)		工况	$\psi$	v0/(m·s−1)		工况	$\psi$	v0/(m·s−1)		工况	$\psi$	v0/(m·s−1)
1-1	0	400		2-1	1/3	400		3-1	2/3	400		4-1	1	400
1-2	0	650		2-2	1/3	650		3-2	2/3	650		4-2	1	650
1-3	0	900		2-3	1/3	900		3-3	2/3	900		4-3	1	900
1-4	0	1 150		2-4	1/3	1 150		3-4	2/3	1 150		4-4	1	1 150
1-5	0	1 400		2-5	1/3	1 400		3-5	2/3	1 400		4-5	1	1 400

下载: 导出CSV 

| 显示表格

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时，试验与数值计算的空泡尺寸

Figure  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的弹体在两种入水速度工况下位移和速度的变化

Figure  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的弹体在两种入水速度工况下位移和速度的变化

Figure  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，通过对比分析图中数据可以发现：

图  6  不同头形因数的截锥形弹体在不同入水速度工况下速度随入水距离的衰减

Figure  6.  Velocity attenuation of truncated cone-shaped projectiles with different head coefficients with water-entry distance at different initial velocities of water entry

下载: 全尺寸图片 幻灯片

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

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

3.   速度衰减分析模型

3.1   模型建立

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

图  7  不同入水速度下，$\psi = 0$的截锥形弹体阻力因数与${v_{\rm{b}}}/{v_0}$的关系

Figure  7.  Resistance factor varying with ${v_{\rm{b}}}/{v_0}$ for the truncated cone-shaped projectile with $\psi = 0$ at different water-entry velocities

下载: 全尺寸图片 幻灯片

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

式中：${a_0}$和${a_1}$为未知参数。

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

表  3  不同工况下的参数${a_0}$和${a_1}$的数值

Table  3.  Values of parameters ${a_0}$ and ${a_1}$ under different working conditions

$\psi$	拟合参数	v0/(m·s−1)
		400	650	900	1150	1400
0	a0	0.617	0.612	0.633	0.587	0.367
	a1	−0.217	−0.212	−0.223	−0.153	0.128
1/3	a0	0.329	0.459	0.477	0.440	0.288
	a1	−0.019	−0.158	−0.168	−0.113	0.077
2/3	a0	0.768	0.707	0.784	0.714	0.553
	a1	−0.387	−0.297	−0.394	−0.289	−0.082
1	a0	1.216	1.204	1.527	1.382	1.416
	a1	−0.344	−0.345	−0.823	−0.580	−0.640

下载: 导出CSV 

| 显示表格

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

式中：a₀₁、a₀₂、a₀₃、a₀₄、a₀₅、a₁₁、a₁₂、a₁₃、a₁₄、a₁₅为待定参数，结果如图8所示。

图  8  对于头形因数不同的截锥形弹体，参数a₀和a₁与马赫数$Ma$的关系

Figure  8.  Changes of parameters a₀ and a₁ with Mach number $Ma$ for the truncated cone-shaped projectiles with different head shape factors

下载: 全尺寸图片 幻灯片

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

表  4  参数${a_0}$、${a_1}$的拟合结果

Table  4.  Fitting results of parameters ${a_0}$ and ${a_1}$

$\psi$	a01	a02	a03	a04	a05	a11	a12	a13	a14	a15
0	0.847	−1.459	2.575	−0.707	−0.887	−0.376	0.918	−1.096	−1.059	1.743
1/3	−0.611	6.078	−13.025	12.901	5.055	1.060	−7.042	15.285	−15.196	5.971
2/3	2.810	−15.067	37.573	−38.342	13.579	−3.246	21.137	−52.932	54.463	−19.504
1	8.033	−53.017	141.081	−153.418	58.737	−10.547	79.685	−213.159	232.929	−89.548

下载: 导出CSV 

| 显示表格

基于上述结果，a系列参数仅与头形因数ψ有关。采用二阶方程对a系列参数进行最小二乘法拟合：

式中：${a}_{ij0}{\text{、}}{a}_{ij1}{\text{、}}{a}_{ij2}{\text{、}}{a}_{ij3}$为待定参数，结果如图9所示。

图  9  $a$系列参数与头形因数$\psi$的关系

Figure  9.  Series parameters of a varyied with head shape factor $\psi$

下载: 全尺寸图片 幻灯片

整体来看拟合曲线同样具有较好的对称关系。由于参数${a_1}$及${a_2}$的变化趋势相反，所以拟合出的a系列参数具有高度对称的特性，这是由二次函数基本特性所决定的，拟合结果见表5。

表  5  ${a_{ij}}$系列参数的拟合结果

Table  5.  Fitting results of ${a_{ij}}$ series parameters

${a_{ij}}$	${a_{ij}}_0$	${a_{ij}}_1$	${a_{ij}}_2$
${a_0}_1$	0.693	−7.536	15.031
${a_{02}}$	−0.864	49.585	−102.332
${a_{03}}$	1.908	−128.117	267.957
${a_{04}}$	−0.654	136.681	−289.498
${a_{05}}$	0.814	−32.000	88.221
${a_1}_1$	−0.239	9.212	−19.657
${a_1}_2$	0.628	−70.273	149.620
${a_1}_3$	−1.463	185.983	−397.312
${a_1}_4$	0.187	−201.803	433.298
${a_1}_5$	1.001	77.280	−167.088

下载: 导出CSV 

| 显示表格

综上所述，当头形因数$\psi$不同的截锥形弹体垂直侵彻入水时可以通过如下公式计算实时阻力因数${C_{\rm{d}}}$：

式中：${a_{ij}} = {a_{ij}}_0 + {a_{ij}}_1\psi + {a_{ij}}_2{\psi ^2}$，具体参数列于表5中。

将式（17）与式（1）联立，可得截锥形弹体速度衰减计算模型：

3.2   算例分析及验证

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

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

图  10  入水速度为400 m/s，头形因数不同的弹丸，速度随入水距离衰减的模拟结果与公式计算结果的比较

Figure  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

下载: 全尺寸图片 幻灯片

图  11  头形因数为1/6，入水速度不同的弹丸，速度随入水距离衰减的模拟结果与公式计算结果的比较

Figure  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

下载: 全尺寸图片 幻灯片

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

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

4.   结　论

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

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

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

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