装满炸药的圆柱套筒是破片型武器最常用的结构。在一端中心起爆后，圆柱套筒会发生膨胀变形并最终破裂成破片，当爆轰产物作用在破片上的压力与破片受到的空气阻力平衡时，破片的速度趋于稳定并达到最大值，此时破片的速度称为破片初速。破片初速是衡量破片杀伤威力及评估防护结构性能的重要指标，因此大量的研究学者采用理论分析、数值模拟及试验等方式对炸药驱动下的圆柱套筒产生的破片初速进行了相关研究及预测。在预测破片初速诸多公式中，Gurney^[1]基于能量守恒定律提出了最为典型的破片初速计算公式（Gurney公式），其破片初速与装药质量和套筒质量的比值及炸药类型相关。Gurney公式可表示为：

式中：$\sqrt {2E}$为格尼能，${v_0}$是破片初速，β是装药质量与套筒质量比值。虽然Gurney公式假定装药瞬间起爆并、炸药的内能全部转换为圆柱套筒与爆轰产物气体的动能且所有的破片沿圆柱套筒轴向分布的初速相同，但一些研究(如文献[2-6])表明，在较大长径比（L/D）下的圆柱套筒中，Gurney公式的计算结果与试验得到的最大破片初速吻合良好，相对误差在10%之内。然而，真实中圆柱套筒并非无限长，端部的稀疏波会降低两端的破片初速，Gurney公式只能预测最大破片初速，并不能预测两端的破片初速。因此，许多研究人员研究了两端稀疏波对破片初速的影响，并基于Gurney公式提出了计算精度更高的修正公式。

Zulkouski^[2]通过大量实验，在Gurney公式中添加了一个指数形式的修正项，以修正起爆端和非起爆端轴向稀疏波对破片初速的影响，修正后的公式可以表示为：

式中：x是圆柱套筒上距离起爆端的轴向距离，D是装药的直径，L是装药的长度。修正公式表明了轴向稀疏波对破片初速的影响与装药直径、距起爆点的距离及装药长度有关。

Charron等^[3]提出了一个破片初速轴向分布的修正公式，假设在圆柱装药的两端分别挖去两个锥体，其中在起爆端锥体的长度等于两倍的装药半径，非起爆端锥体的长度等于一倍装药半径，以此作为圆柱套筒起爆的等效模型，得到速度分布结果：

式中：F(x)为修正系数，可以表示为$F(x) = 1 - {\text{min}}\left\{ {x/2R,1,\left( {L - x} \right)/R} \right\}$，其中R是装药半径。

式(2)计算出的起爆端的速度为0，式(3)计算出的起爆端和非起爆端的速度为0，这显然与试验结果不符。为修正轴向稀疏波对破片初速分布的影响，解决端部为0的问题，Huang等^[4]利用X射线实验清楚的拍摄到了不同时刻下破片轴向速度分布情况，并根据试验结果，在Gurney公式的基础上，添加了指数形式的修正函数：

式中：${v_0}$为Gurney公式计算的破片初速；起爆端修正项${F_1}(x/D) = 1 - {A_1}{{\text{e}}^{ - {B_1}x/D}}$，非起爆端修正项${F_2}((L - x)/D) = 1 - {C_1}{{\text{e}}^{ - {D_1}(L - x)/D}}$，A₁、B₁为起爆端的修正参数，C₁、D₁为非起爆端的修正参数。Huang等^[4]根据L/D=3.27圆柱套筒的试验数据，在不考虑端部稀疏波相互作用条件下，采用最小二乘法拟合得到参数A₁、B₁、C₁、D₁的固定值，进而得到很高的计算精度的破片初速修正公式：

Gao等^[5]基于对有效装药质量的考虑，采用与文献[4]类的方法，通过缩减有效装药质量来修正轴向稀疏波对破片初速分布的影响，得到如下公式：

式(5)和式(6)虽然解决了起爆端和非起爆端为0的问题，但都是从L/D＞3试样的试验数据推导出来并进行验证的，对L/D较小的套筒轴向初速分布预测并不准确，适用性较窄。因此，Liu等^[6]采用正交闪光X射线摄影技术，对不同L/D下圆柱套筒的轴向破片初速分布进行了测量和分析，在式(5)形式上添加了考虑L/D影响的修正项：

式(7)虽然考虑了L/D对破片轴向初速分布的影响，但当试样参数改变时，其破片初速预测误差较大，且L/D适用范围较窄。因此，应该提出一个更准确的用于描述不同L/D下圆柱套筒破片初速预测公式。本文首先通过试验验证的数值模型研究L/D对破片初速的影响，在此基础上构造与L/D相关的受轴向稀疏波影响的修正函数，并构建适用于L/D≥1圆柱套筒的轴向速度初速分布计算模型，最后通过试验和数值模拟验证所提模型的准确性。

1.   数值模拟

1.1   数值模拟方法

采用AUTODYN软件^[7-9]的ALE（arbitraty Lagrangian-Eulerain）算法来数值模拟L/D对破片初速的影响，模型见图1。套筒采用采用AISI1045钢，炸药选择B炸药^[4-5,10]。套筒的Lagrange网格尺寸设置为0.3 mm^[11]。由于结构和载荷的对称性，为了节省计算时间，建立四分之一的有限元模型。考虑边界对套筒动态响应的影响^[12]，在空气域边界设置非反射边界条件，避免压力在边界上的反射。为获得不同轴向位置的破片初速分布，沿轴向不同位置等距设计若干个高斯点。为了研究L/D对破片初速分布的影响，设置了8个不同L/D的数值模拟试样，试样的相关参数见表1，δ为圆柱套筒的厚度。

图  1  数值计算模型

Figure  1.  Numerical calculation model

下载: 全尺寸图片 幻灯片

表  1  不同L/D的数值模拟试样

Table  1.  Simulation samples with different L/D

模型	D/mm	L/mm	L/D	δ/mm	模型	D/mm	L/mm	L/D	δ/mm
1	23.6	23.6	1.0	3.04	5	23.6	59.0	2.5	3.04
2	23.6	29.5	1.25	3.04	6	23.6	70.8	3.0	3.04
3	23.6	35.4	1.5	3.04	7	23.6	94.4	4.0	3.04
4	23.6	47.2	2.0	3.04	8	23.6	118	5.0	3.04

下载: 导出CSV 

| 显示表格

1.2   材料模型

Johnson-Cook模型^[13]常用于描述金属材料在高应变、高应力以及高温情况下的弹塑性变形过程，因此采用Johnson-Cook模型来描述ASI1045钢的动态变形过程：

式中：$\sigma$表示材料屈服应力；$\varepsilon _{\text{p}}^n$为等效塑性应变；${\varepsilon ^ * }$为规范化等效塑性应变；${T^ * } = \left( {T - {T_{\text{r}}}} \right)/\left( {{T_{\text{m}}} - {T_{\text{r}}}} \right)$表示约化温度（无量纲），${T_{\text{m}}}$和${T_{\text{r}}}$分别为材料的融化温度和室温；A_JC、B_JC、C_JC、n和m均为材料相关常数。

Johnson-Cook失效模型常用于描述金属的失效断裂过程：

式中：${\varepsilon _{\text{f}}}$为失效时刻应变值，${\sigma ^ * }$表示压力与有效应力的比值，${D_1}$～${D_5}$则为失效参数，由试验测得。AISI1045材料参数见表2。

表  2  AISI1045的J-C本构模型参数^[15-16]

Table  2.  J–C constitutive model parameters of AISI1045^[15-16]

密度/(g·cm−3)	AJC/MPa	BJC/MPa	n	CJC	m	D1	D2	D3	D4	D5
7.83	507	320	0.28	0.064	1.06	0.15	0.72	1.66	0.005	-0.84

下载: 导出CSV 

| 显示表格

采用JWL状态方程^[14]描述B炸药的动态力学性能：

式中：p为压力，A、B、R₁、R₂及ω为常数，可通过炸药圆筒试验得到，V为爆轰产物的相对比容，E为初始体积内能。B炸药的材料参数见表3。

表  3  B炸药的材料参数^[4]

Table  3.  Material parameters of Comp B^[4]

爆速/(m·s−1)	爆轰压力/GPa	E/(kJ·m−3)	AJWL/GPa	BJWL/GPa	R1	R2	ω
7980	29	8.5×106	542	7.68	4.2	1.1	0.24

下载: 导出CSV 

| 显示表格

空气使用AUTODYN软件的理想气体状态方程进行描述：

式中：$e$为空气的比内能，空气密度为${\rho _{{\text{air}}}}$=1.225 kg/m³，$\gamma$=1.4为空气的绝热指数^[17]。

1.3   准确性验证

本文主要研究不同L/D下破片初速的分布，而破片初速的分布主要受轴向（x方向）稀疏波的影响，因此只要验证轴向稀疏波对圆柱套筒动态响应的影响，就能验证模型的准确性。Huang等^[4]利用X光试验机清晰地拍摄到破片沿轴向位置的分布情况。本文根据Huang等^[4]的试验工况建立数值计算模型，并与Huang等X光试验机拍摄到的破片沿轴向位置的分布情况进行对比，验证数值模型的准确性，试验样本参数见表4。

表  4  黄广炎的试验样本参数^[4]

Table  4.  Parameters of HUANG G Y's test sample ^[4]

工况	套筒材料	装药材料	D/mm	L/mm	δ/mm
1	ASI1045	B炸药	23.60	77.30	3.04
2	ASI1045	B炸药	23.56	77.15	6.69

下载: 导出CSV 

| 显示表格

从图2(a)中可以看出，破片初速随着Euler网格尺寸的减小而增加，而网格尺寸为0.4和0.5 mm的破片初速差异不大。此外，数值模拟的初速分布（Euler网格尺寸为0.5 mm）与两组试验得到的初速分布具有较好的一致性。在图2(b)中可以看出，由于起爆端和非起爆端轴向稀疏波的影响，起爆端附近的套筒发生较少的轴向断裂，非起爆端附近的套筒发生较多的轴向断裂，而套筒中部未发生轴向断裂，数值模拟（Euler网格尺寸为0.5 mm）的套筒断裂结果与文献[4]中第一组试验结果吻合性较好。因此，尺寸为0.5 mm的Euler网格足以提供准确的结果，且该数值模型能够模拟轴向稀疏波对破片初速的影响。

图  2  试验结果与数值模拟结果对比

Figure  2.  Comparison of experimental results and numerical results

下载: 全尺寸图片 幻灯片

2.   数值计算结果及分析

2.1   长径比对破片初速的影响

为研究L/D对圆柱套筒初速的影响，统计了不同L/D下不同轴向位置对应的破片初速，破片初速分布如图3所示。图3(a)中各个L/D下的装药质量与套筒质量比值相同，按照Gurney公式^[1]进行计算，不同L/D下每个轴向位置的破片速度均是相同的，但显然数值模拟结果和Gurney公式^[1]计算结果不一致，且不同L/D下的同一位置数值模拟结果也不一致，因此不能忽略L/D对破片初速沿轴向分布的影响。

图  3  不同L/D下破片初速分布

Figure  3.  Fragment initial velocity distribution with different L/D

下载: 全尺寸图片 幻灯片

从图3中可以明显看出，不同L/D的套筒初速都受到轴向稀疏波的影响，破片初速分别均呈现出相似的中间高，两端低的变化趋势。炸药起爆后，起爆端轴向稀疏波会立即传播到起爆端，从而削弱了冲击波与爆轰产物对起爆端套筒的加速作用，而非起爆端的轴向稀疏波只有在爆轰波传播到非起爆端时才会产生。爆轰波传播到非起爆端时冲击波已经加速其他位置套筒一段时间，且起爆端轴向稀疏波在爆轰产物中随传播过程迅速减弱。因此，非起爆端的套筒初速受轴向稀疏波的影响较小，导致不同L/D下起爆端的破片初速均低于非起爆端。

L/D越小，破片初速越低，如图3(b)所示。为探究L/D对破片初速的影响机理，分析了L/D对装药中稀疏波传播过程及套筒加速过程。图4为轴向稀疏波的传播过程，T₁为起爆端轴向稀疏波传播到非起爆端的时间，T₂为非起爆端轴向稀疏波传播到起爆端的时间，基于爆轰产物膨胀理论^[18]可知，起爆端轴向稀疏波和非起爆端轴向稀疏波的传播速度分别为装药爆速的1倍和0.5倍，因此，T₂=2T₁，与数值模拟结果一致，这再次证明了数值模拟方法的准确性。

图  4  稀疏波传播过程

Figure  4.  Process of rarefaction wave propagation

下载: 全尺寸图片 幻灯片

圆柱套筒从速度零加速到破片初速的时间为T₃，当L/D较小时，T₁+T₂较小，则T₃＞T₁+T₂，说明非端轴向稀疏波传播到起爆端时，破片加速未完成；当L/D较大时，T₁+T₂较大，则T₃＜T₁+T₂，说明部分破片只受到起爆端轴向稀疏波的影响。这是因为，当L/D逐渐增大时，T₁+T₂逐渐增大，只受到起爆端轴向稀疏波影响的破片数量增多，破片初速出现增加的趋势。由于起爆端轴向稀疏波在传播过程中是逐渐减弱的^[4]。因此，当L/D足够大时，起爆端轴向稀疏波和非起爆端轴向稀疏波中间存在基本不受稀疏波影响的区域，稀疏波未影响区域内的破片初速非常接近Gurney公式^[1]计算结果，如图5所示：最大破片初速随着L/D的增加呈现指数型增长；L/D越大，最大破片初速增量越少；当L/D达到5时，模拟结果与Gurney公式的相对误差仅为0.15%。

图  5  不同L/D下破片的最大初速

Figure  5.  Maximum initial velocity of fragments for different L/D

下载: 全尺寸图片 幻灯片

2.2   已有的破片初速计算公式适用性验证

如上文所述，L/D通过影响轴向稀疏波传播来影响破片加速历程，进而影响破片初速分布。因此，有必要验证已有破片初速计算公式在不同L/D下的适用性。式(2)和式(3)在端部的速度计算结果为0，与实际不符合，本节不再继续验证其适用性。本节主要通过数值结果来验证式(5)、式(6)、式(7)的适用性。不同L/D下破片初速分布及与数值结果的相对误差分布如图6～10所示。可以看出，式(5)、式(6)能够准确预测L/D=3, 4, 5的试样的破片初速分布，且L/D越大，预测值越接近数值计算值；但不能准确预测L/D=1、2试样的轴向初速分布，L/D越小，预测值越偏离数值计算值，当L/D=1时，式(5)、式(6)的最大相对误差分别为48.7%、46.7%，平均相对误差分别为31.2%、26.5%，这是因为式(5)、式(6)在推导过程中未考虑L/D对破片初速的影响。式(7)虽然添加了考虑L/D影响的修正项，在预测L/D=1的试样时具有较高的准确性，但最大相对误差仍高达28.1%。文献给出了式(7)适用的L/D范围（0.5～3），但当L/D=2和3时，式(7)预测值的精度远远不如式(5)、式(6)，最大相对误差分别为29.4%、39.3%。因此，需要提出一个更准确的用于描述不同L/D下圆柱套筒破片初速预测公式。

图  6  L/D=1时下破片速度分布和相对误差分布对比情况

Figure  6.  Comparison of fragment velocity distribution and relative error distribution while L/D=1

下载: 全尺寸图片 幻灯片

图  7  L/D=2时下破片速度分布和相对误差分布对比情况

Figure  7.  Comparison of fragment velocity distribution and relative error distribution while L/D=2

下载: 全尺寸图片 幻灯片

图  8  L/D=3时下破片速度分布和相对误差分布对比情况

Figure  8.  Comparison of fragment velocity distribution and relative error distribution while L/D=3

下载: 全尺寸图片 幻灯片

图  9  L/D=4时下破片速度分布和相对误差分布对比情况

Figure  9.  Comparison of fragment velocity distribution and relative error distribution while L/D=4

下载: 全尺寸图片 幻灯片

图  10  L/D=5时下破片速度分布和相对误差分布对比情况

Figure  10.  Comparison of fragment velocity distribution and relative error distribution while L/D=5

下载: 全尺寸图片 幻灯片

3.   轴向初速分布模型

通过上文分析可知，L/D会导致破片速度分布存在差异，而Huang等^[4]的模型中，A₁、B₁、C₁、D₁是与L/D不相关的常数，因此式(5)存在较大的误差。从图3中可以看出，L/D虽然影响破片的初速分布，但初速分布的形式基本不发生变化。因此，本文继续采用式(4)指数函数的修正项修正不同L/D下破片初速分布，但此时A₁、B₁、C₁、D₁是随L/D变化的变化量，即：

式中:考虑L/D影响的起爆端修正项${F_1}(x/D) = 1 - {A_1}(L/D){{\text{e}}^{ - {B_1}(L/D) \cdot x/D}}$,考虑L/D影响的非起爆端修正项${F_2}((L - x)/D) = 1 - {C_1}{\text{(}}L/D{\text{)}}{{\text{e}}^{ - {D_1}(L/D) \cdot (L - x)/D}}$，A₁、B₁、C₁、D₁在不同L/D条件下的值可以通过最小二乘法拟合数值结果得到，然后再分析其与L/D的关系。

为得到Gurney公式^[1]计算的破片初速，需计算炸药的格尼能$\sqrt {2E}$。Hossein和Keshavarz^[19]在拟合实验数据的基础上建立了只与炸药化学式${{\text{C}}_a}{{\text{H}}_b}{{\text{N}}_c}{{\text{O}}_d}$中原子有关的计算公式：

式中：${\rho _0}$为炸药的密度。式(13)无需对爆轰产物成分及爆炸热量进行假设，更能反映炸药的本质，计算结果具有很高的准确性。

当L/D较小时，两端轴向稀疏波相互影响，因此不能通过单独分析起爆端和非起爆端得到系数A₁、B₁、C₁、D₁的值^[20]，只能通过整体拟合获得各个系数的值。由此得到不同L/D下如表5所示的系数，系数与L/D的关系如图7所示，显然，在不同的L/D条件下，系数发生了变化。

表  5  7个数值模型得到的A₁、B₁、C₁、D₁系数

Table  5.  Coefficients A₁, B₁, C₁ and D₁ obtained from seven numerical models

L/D	A1	B1	C1	D1	L/D	A1	B1	C1	D1
1.0	0.465	2.282	0.425	0.823	2.5	0.465	1.293	0.296	1.970
1.25	0.465	1.983	0.326	1.121	3.0	0.465	1.267	0.303	2.030
1.5	0.465	1.607	0.299	1.399	4.0	0.466	1.175	0.301	2.072
2.0	0.466	1.314	0.291	1.83	5.0	0.466	1.137	0.297	2.086

下载: 导出CSV 

| 显示表格

在图11中能看到：当L/D较大时，系数都会趋于平缓，这是因为L/D较大时，两端轴向稀疏波相互不影响且破片初速基本不发生变化，L/D的变化并不引起系数的显著变化，这进一步说明修正公式的合理性。根据图11的曲线，参数的函数采用如下形式：

图  11  系数A₁、B₁、C₁、D₁与L/D的关系

Figure  11.  Relationship between coefficients A₁, B₁, C₁, D₁ and L/D

下载: 全尺寸图片 幻灯片

式中：q₁、q₂、q₃、q₄、m₁、m₂、m₃、m₄、n₁、n₂、n₃、n₄为待定参数，可以通过图11中A₁、B₁、C₁、D₁和L/D的关系拟合得到。最终，可获得计算不同L/D下圆柱套筒的破片轴向初速分布公式：

式(15)虽然是L/D在1～5范围内推导出来的，但L/D越大，轴向稀疏波对式(15)中的系数影响越小，当L/D大于4时，式(15)中的系数变化量很小，基本可以忽略L/D对其的影响，这就意味着式(15)可以适用于更大L/D下圆柱套筒的初速分布预测，因此式(15)适用于L/D≥1的圆柱套筒。

4.   新公式普适性验证

为验证提出的破片初速分布计算模型的准确性，采用文献[4]和文献[21]的试验结果进行验证。文献[4]的两组试样参数见表4；文献[21]中$\beta$=0.429、L/D=2、$\sqrt {2E}$=2840 m/s。图12为文献结果和式(15)的计算结果对比情况。在图12在能明显看出公式计算结果和试验结果吻合性较好，平均相对误差不超过6%，表明所提公式具有足够的精度。

图  12  试验结果和式(15)计算结果对比

Figure  12.  Comparison of experimental results and calculation results of Eq. (15)

下载: 全尺寸图片 幻灯片

为进一步验证式(15)的可靠性，对另外9个带有不同L/D、装药类型和$\beta$的圆柱套筒进行了数值计算，用于数值计算和公式计算的相关参数见表6，HMX炸药的数值计算参数见表7。三组工况（V1、V2、V3）下破片初速分布的计算结果和数值结果的对比分别如图13、图14、图15所示。可以看出，式(15)计算结果与数值计算结果吻合较好，最大相对误差不超过10%，平均相对误差不超过4%。以上表明，式(15)能够准确预测圆柱套筒在不同L/D下的破片初速分布情况。

表  6  9个用于验证公式的试样参数

Table  6.  Nine sample parameters used to validate the formula

工况	装药	D/mm	L/mm	D/L	δ/mm	工况	装药	D/mm	L/mm	D/L	δ/mm	工况	装药	D/mm	L/mm	D/L	δ/mm
V1-1	B炸药	20	20	1.0	3.0	V2-1	B炸药	30	30	1.0	4.0	V3-1	HMX	20	20	1.0	3.0
V1-2	B炸药	20	40	1.2	3.0	V2-2	B炸药	30	60	1.2	4.0	V3-2	HMX	20	40	1.2	3.0
V1-3	B炸药	20	60	3.0	3.0	V2-3	B炸药	30	90	3.0	4.0	V3-3	HMX	20	60	3.0	3.0
注：所有工况中套筒均材料均为ASI1045.

下载: 导出CSV 

| 显示表格

表  7  HMX的材料参数^[22]

Table  7.  The material parameters of HMX^[22]

爆速/(m·s−1)	爆轰压力/GPa	E/(kJ·m−3)	AJWL/GPa	BJWL/GPa	R1	R2	ω
9110	42.0	10.5×106	778.28	7.07	4.2	1.0	0.30

下载: 导出CSV 

| 显示表格

图  13  仿真结果与式(15)计算结果对比（Cases V1）

Figure  13.  Comparison of numerical results and calculation results of Eq. (18) (Cases V1)

下载: 全尺寸图片 幻灯片

图  14  仿真结果与式(15)计算结果对比（Cases V2）

Figure  14.  Comparison of numerical results and calculation results of Eq. (18) (Cases V2)

下载: 全尺寸图片 幻灯片

图  15  仿真结果与式(15)计算结果对比（Cases V3）

Figure  15.  Comparison of numerical results and calculation results of Eq. (18) (Cases V3)

下载: 全尺寸图片 幻灯片

5.   结　论

本文通过试验验证的数值模拟研究了不同长径比L/D下圆柱套筒的破片初速分布。

不同L/D下的破片轴向速度分布均呈现出相似的中间高，两端低的变化趋势，且起爆端的破片初速均低于非起爆端。

L/D直接影响着轴向稀疏波传播及破片加速历程，L/D越小，破片初速衰减量越大，导致破片初速分布不能用已有公式进行准确计算。

本文在之前的公式基础上，提出了与L/D相关的受轴向稀疏波影响的修正函数，通过拟合数值计算结果得到修正函数中的系数。

此外，通过试验和数值模拟对所提出的公式进行了进一步验证，验证结果表明，平均误差不超过6%，本文共识的计算结果与试验结果和数值计算结果吻合较好。因此，本文提出的公式能够准确计算不同L/D下的破片初速分布。