基于力热耦合材料模型的Steven试验数值模拟方法

楼建锋; 张延耿; 周婷婷; 洪滔

doi:10.11883/1001-1455(2017)05-0807-06

基于力热耦合材料模型的Steven试验数值模拟方法

doi: 10.11883/1001-1455(2017)05-0807-06

北京应用物理与计算数学研究所，北京 100094

基金项目:

国家自然科学基金项目 11302031

国家自然科学基金项目 11402031

中国工程物理研究院科学技术发展基金项目 2014B0101014

中国工程物理研究院安全弹药研发中心开放基金项目 RMC2014B02

详细信息

作者简介:
楼建锋(1980—)，男，博士，副研究员

通讯作者:
张延耿，zhang_yangeng@iapcm.ac.cn

中图分类号: O381
计量
- 文章访问数: 4667
- HTML全文浏览量: 1322
- PDF下载量: 250
- 被引次数: 6
出版历程
- 收稿日期: 2016-01-29
- 修回日期: 2016-06-12
- 刊出日期: 2017-09-25

Numerical method for simulating Steven test basedon thermo-mechanical coupled material model

Institute of Applied Physics and Computational Mathematics, Beijing 100094, China

摘要

摘要: 在力热耦合材料模型中，增加炸药自热放能模型，建立了Steven试验的力-热-化耦合的数值模拟方法。数值计算模型中，应力应变关系采用双线性硬化弹塑性模型，炸药受力后的热作用采用各向同性热材料模型，炸药的化学反应采用Arrhenius反应率函数，同时还考虑了升温和熔化对材料力学、热学性能的影响。针对标准Steven试验，通过数值分析得到了靶板的变形情况和炸药点火的速度阈值，将计算结果与实验数据进行了比较，两者符合较好。表明该方法可以较好地模拟Steven试验，而且与以往的分析模型和方法相比，本文的方法不需要增加经验性的点火准则和判据，具有更广泛的适用性，可以为研究低速撞击条件下炸药的力热响应和局域化点火问题提供参考。
- Steven试验 /
- 力热耦合 /
- 低速撞击 /
- 非冲击点火
Abstract: In this paper, we figured out a numerical simulation method involving the mechanical, thermal and chemical properties of the Steven test based on the thermo-mechanical coupled material model to simulate the Steven test of the plastic bonded explosive 9501. In this model, the stress-strain relationship is described by the dynamic plasticity model, the impact-induced thermal effect depicted by the isotropic thermal material model, the chemical reaction is described by the Arrhenius reaction rate law, with the effects of heating and melting on mechanical properties and thermal properties of materials also taken into account. Specific to the standard Steven test, the numerical model was validated by comparing the obtained deformation of the target and the ignition threshold of the PBX 9501 with the experimental data in the references. The calculated results are in good agreement with the experimental data, suggesting that this method is capable of simulating the Steven test. Compared to the previous models, this model does not need to incorporate experiential ignition criterion and therefore can be used more widely in the study of thermo-mechanical responses and local ignition of explosives subjected to low velocity impact.
- Steven test /
- thermo-mechanical coupled model /
- impact /
- non-shock initiation

HTML全文

混凝土、岩石类脆性材料在土木工程中大量使用，对于一些重要的防护结构，如大坝、核电站安全壳等，设计时必须考虑强冲击载荷作用。而这些材料往往含有大量的原生缺陷，在冲击载荷作用下的动态响应与静载荷下显著不同，因此研究其动态力学性能，特别是与裂纹扩展相关的动态断裂韧性十分重要。

由于理论和实验技术的困难，高应变率下脆性材料动态断裂行为的研究一直进展缓慢。M.Nakano等^[1]和Y.Yamauchi等^[2]在分离式Hopkinsion压杆(SHPB)上采用含中心直裂纹巴西圆盘试样(CCBD)进行了复合型动态断裂韧性实验。由于实验技术的不成熟，得到的实验数据较离散，但仍然发现结果与静态实验的显著不同，不仅受实验应变率的影响，同时与试样大小和裂纹长度密切相关。为了更好地测量脆性材料的断裂韧性，有学者采用含裂纹或缺口的半圆弯曲(SCB)试样代替传统的CCBD试样，并结合有限元模拟得到静态断裂韧性的计算公式^[3-5]。R.Chen等^[6]和F.Dai等^[7-8]将SCB试样引入脆性材料动态断裂韧性的测量，并采用先进的激光测量系统测量试样裂纹或缺口的张开位移，假设实验中试样达到应力平衡，忽略惯性效应的影响，沿用静态分析结果计算动态断裂韧性。然而在SHPB实验中，大多数脆性材料很难达到应力平衡，惯性效应的影响始终存在，上述研究中没有明确提出动态加载下静态分析的适用范围和可靠性。

本文中采用有限元方法，首先对静态含裂纹SCB实验进行数值模拟，结合量纲分析推导适用于复合型断裂的Ⅰ型应力强度因子(K_Ⅰ)拟合公式；对基于SHPB技术的动态复合型SCB实验进行数值模拟，计算试样裂尖的动态应力强度因子(K_Ⅰd)，分析试样半径、支座间距、相对裂纹长度及裂纹倾角对断裂韧性的影响，讨论静态公式的适用范围。

1. 静态复合型加载SCB实验模拟及Ⅰ型应力强度因子拟合公式

1.1 含裂纹SCB实验

采用有限元软件ANSYS程序对静态复合型断裂SCB实验进行了数值模拟，实验示意图如图 1所示，其中：a为裂纹长度，R为试样半径，F试样顶部集中载荷，β为裂纹倾角(裂纹与试样中心轴间的夹角)，s为半支座间距。

图 1 静态复合型加载SCB实验简图

Figure 1. Diagram of static SCB test under mixed mode loading

下载: 全尺寸图片幻灯片

试样采用线弹性本构，弹性模量为30GPa，密度2500kg/m³，泊松比0.2。试样半径为50mm，厚25mm，支座间距为60mm，裂纹相对长度a/R从0.01变化至0.8，裂纹倾角β从0°变化至45°。

依据图 1建立有限元模型，单元类型为solid 186六面体实体单元。为了更好地反映裂纹应力奇异的影响，围绕裂尖一圈采用命令KSCON实现退化奇异等参元。应力强度因子采用相互作用积分法^[9]计算，该法是基于J积分的一种新型应力强度因子计算方法，其精度与数值外插法相当。对于静态问题，ANSYS提供了支持该法的命令流，适于批量操作，对于而动态问题则没有，因此本文中采用不同的方法计算静态和动态应力强度因子。

1.2 Ⅰ型应力强度因子拟合公式

根据量纲分析并结合文献[3-8]进行计算，选择Ⅰ型应力强度因子K_Ⅰ拟合公式形式为：

$K_{\mathrm{I}}=\cos ^{2} \beta \frac{F \sqrt{a}}{B R}\left[g\left(\frac{a}{R}, \beta\right) \frac{s}{R}+h\left(\frac{a}{R}, \beta\right)\right]$

(1)

式中：B为试样厚度，g、h为拟合函数。式(1)的前半部满足应力强度因子的量纲形式，且满足当a趋于0时，公式退化为弯曲梁单边斜裂纹的K_I公式。后半部为量纲一形状因子，反映了裂纹长度a、支座间距s和裂纹倾角β的影响，拟合函数g、h采用多项式形式，表达式为：

$\left\{ {\begin{array}{*{20}{l}} {g\left( {\frac{a}{R},\beta } \right) = {d_1}(\beta ) + {d_2}(\beta ){{\left( {\frac{a}{R}} \right)}^{0.5}} + {d_3}(\beta ){{\left( {\frac{a}{R}} \right)}^{1.5}} + {d_4}(\beta ){{\left( {\frac{a}{R}} \right)}^{2.5}} + {d_5}(\beta ){{\left( {\frac{a}{R}} \right)}^{6.5}}\;\;\;(2)}\\ {h\left( {\frac{a}{R},\beta } \right) = {e_1}(\beta )\left( {\frac{a}{R}} \right) + {e_2}(\beta ){{\left( {\frac{a}{R}} \right)}^2} + {e_3}(\beta ){{\left( {\frac{a}{R}} \right)}^3}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;(3)} \end{array}} \right.$

式中：d_n和e_n(n=1, 2, …, 5)为裂纹倾角β的拟合系数，均采用直线方程拟合：

$\left\{\begin{array}{ll}{d_{1}(\beta)=0.129 \beta+2.24, } & {d_{2}(\beta)=-2.09 \beta+77.60,} \\ {d_{3}(\beta)=-65.59 \beta+2649.61, } & {d_{4}(\beta)=-149.11 \beta+5928.51,} \\ {d_{5}(\beta)=-5.05 \beta+185.89, } & {e_{1}(\beta)=9.96 \beta-400.21,} \\ {e_{2}(\beta)=82.67 \beta-3327.49, } & {e_{3}(\beta)=38.87 \beta-1520.91}\end{array}\right.$

(4)

采用本文公式模拟不同倾角β下的应力强度因子，结果如图 2所示。最大拟合误差小于10%。从图 2可以看出K_Ⅰ值随裂纹倾角β增加而降低，当倾角大于40°时，K_Ⅰ值趋向0，反映出断裂由纯Ⅰ型向复合型直至纯Ⅱ型转变的过程。

图 2 静态复合型加载SCB实验应力强度因子拟合结果

Figure 2. Fitted effect of stress intensity factor in static SCB test under mixed mode loading

下载: 全尺寸图片幻灯片

2. 基于SHPB装置的动态SCB实验的模拟

2.1 有限元模型和动态应力强度因子计算

模拟采用LS-DYNA程序，其中试样的材料常数与静态完全相同，SHPB装置的入射杆和透射杆采用线弹性材料：弹性模量200GPa，密度7850kg/m³，泊松比0.3。直径100mm，长度2000mm。所用单元类型均采用Solid 164六面体实体单元。有限元模型如图 3所示。载荷为加在左端入射杆端面的三角形速度脉冲，如图 4所示。

图 3 基于SHPB装置的动态SCB实验有限元模型

Figure 3. Finite element model of dynamic SCB test based on SHPB device

下载: 全尺寸图片幻灯片

图 4 三角形速度脉冲

Figure 4. Triangular velocity pulse

下载: 全尺寸图片幻灯片

由于动态模拟无法实现退化奇异单元，在裂尖附近采用稠密单元，减小单元网格尺寸影响。采用2种方法计算裂尖的动态应力强度因子K_Ⅰd：(1)由裂纹面张开位移通过数值外插法^[10]得到，记为K_Ⅰd1；(2)根据试样左右端面的接触力-时间曲线平均后带入式(1)计算得到，记为K_Ⅰd2。K_Ⅰd1为裂尖实际动态应力强度因子，K_Ⅰd2为采用静态公式得到的动态应力强度因子。

2.2 动态纯Ⅰ型加载SCB实验的模拟及分析

首先对裂纹倾角β=0°的纯Ⅰ型加载SCB实验进行了模拟，研究相对裂纹长度a/R、试样半径R、以及支座间距s对裂尖实际动态应力强度因子K_Ⅰd1的影响。图 5所示为固定试样半径R=50mm、相对支座间距s/R=0.6、裂纹相对长度a/R从0.1变化至0.8时的K_Ⅰd1。图 6所示为a/R=0.4，s/R=1.2、试样半径R从40mm变化至70mm时的K_Ⅰd1。图 7所示为试样半径R=50mm、a/R=0.4、相对支座间距s/R从0.4变化至0.8时的K_Ⅰd1。

图 5 不同相对裂纹长度下的应力强度因子变化规律

Figure 5. Variation of stress intensity factor for different relative crack lengths

下载: 全尺寸图片幻灯片

图 6 不同试样半径下应力强度因子的变化规律

Figure 6. Variation of stress intensity factor for different sample radiuses

下载: 全尺寸图片幻灯片

图 7 不同支座间距下应力强度因子的变化规律

Figure 7. Variation of stress intensity factor for different distances between two supports

下载: 全尺寸图片幻灯片

由图 5~7可以看出：

(1) K_Ⅰd1随裂纹相对长度a/R、支座间距s的增大而增大，随试样半径R的增大而减小。

(2) K_Ⅰd1-t曲线在加载过程中均单调增加，而达到峰值后出现震荡；随a/R、试样半径R增加，震荡加剧，较大和较小的支座间距也会导致震荡加剧。这种震荡反映了试样应力不平衡导致的惯性效应的影响。当a/R <0.4、试样半径小于60mm、s/R≈1.2时，试样应力不平衡影响较小。

为了进一步分析静态公式的适用范围，对a/R在0.1~0.6范围内裂尖实际K_Ⅰd1和静态公式得到的K_Ⅰd2进行了比较，结果如图 8所示，从图中可以看出，两者在a/R=0.4时吻合最好。按下式得到不同a/R情况下两者的相对误差δ并绘图，如图 9所示，其中δ的表达式为

$\delta=\frac{K_{\mathrm{IdI}}-K_{\mathrm{Id} 2}}{K_{\mathrm{Idl}}} \times 100 \%$

(5)

图 8 K_Ⅰd1和K_Ⅰd2的比较(a/R=0.4)

Figure 8. Comparison of K_Ⅰd1 and K_Ⅰd2 (a/R=0.4)

下载: 全尺寸图片幻灯片

从图 9可以看出，当a/R & #60;0.4时，δ趋近于0，两者相对误差在10%左右，即此时惯性效应的影响较小，静态公式具有足够的精度。

图 9 静态公式的相对误差随a/R的变化

Figure 9. Relative error of static formula for different a/R

下载: 全尺寸图片幻灯片

2.3 动态复合型加载SCB实验模拟及分析

为了研究复合型加载SCB实验中K_Ⅰd的变化规律，在纯I型研究的基础上，固定试样半径R=50mm、支座间距s=1.2R，裂纹倾角β从10°变化至40°，每个倾角下相对裂纹长度a/R从0.2变化至0.8建立有限元模型进行研究。为了便于分析，取裂尖K_Ⅰd1-t曲线加载段的峰值K_Ⅰd1max进行比较。图 10所示为不同a/R情况下K_Ⅰd1max随裂纹倾角β的变化规律，图 11所示为不同裂纹倾角β情况下K_Ⅰd1max随a/R的变化规律。

图 10 K_Ⅰd1max随裂纹倾角β的变化规律

Figure 10. Variation of K_Ⅰd1max with crack angle β

下载: 全尺寸图片幻灯片

图 11 K_Ⅰd1max随裂纹相对长度a/R的变化规律

Figure 11. Variation of K_Ⅰd1max with relative crack length a/R

下载: 全尺寸图片幻灯片

由图 10可以发现，对应相同的裂纹相对长度，K_Ⅰd1max随裂纹倾角β的增大而减小，且a/R越大，下降速率越大。由图 11可以发现，对应相同裂纹倾角，K_Ⅰd1max随裂纹相对长度a/R的增大而增大，而当裂纹倾角β达到40°时，K_Ⅰd1max随裂纹相对长度a/R的增大而减小，变化规律与静态模拟结果类似，如图 2所示。

为了进一步分析复合加载条件下静态拟合公式的适用性，按式(5)计算了裂尖实际应力强度因子峰值K_Ⅰd1max与静态公式得到的应力强度因子峰值K_Ⅰd2max之间的相对误差δ_m，如表 1所列。可以看出，当裂纹倾角β在10°~40°之间时，裂纹相对长度a/R在0.2~0.4之间时，静态公式的相对误差小于10%，精度较高，而超出此范围，相对误差增大至20%~50%，即静态公式失效。

表 1 复合型加载下静态公式计算K_Ⅰdmax的相对误差δ_m

Table 1. Relative errorδ_m of K_Ⅰdmax in the static formula under mixed-mode loading

β/(°)	δ_m/%
β/(°)	a/R=0.2	a/R=0.4	a/R=0.6	a/R=0.7	a/R=0.8
10	9.6	6.6	-0.7	3.7	26.8
20	6.4	9.8	19.8	23.2	40.3
30	5.4	5.6	24.1	39.2	49.1
40	9.3	0.9	17.2	34.8	138.0

下载: 导出CSV

| 显示表格

3. 结论

本文中通过有限元分析方法，对静态和动态复合型加载SCB实验进行了数值模拟研究，得到如下结论：

(1) 根据静态实验的模拟结果，结合量纲分析得到了适合于裂纹倾角在0°~45°，裂纹相对长度在0.1~0.8范围内的应力强度因子拟合公式，计算结果最大误差小于10%。

(2) 对于含纯Ⅰ型中心裂纹的动态SCB实验，试样裂尖的动态断裂韧性随试样半径增大而减小；随支座间距和裂纹相对长度的增大而增大。

(3) 对于含复合型裂纹的动态SCB实验，当裂纹相对长度在0.2~0.7范围内时，裂尖的动态断裂韧性随着裂纹倾角的增大而减小；当裂纹倾角在10°~30°范围内时，裂尖的动态断裂韧性随着裂纹相对长度的增大而增大。

(4) 通过采用位移外插法和静态拟合公式2种方法计算裂尖的动态断裂韧性，并计算相对误差发现，对于含纯Ⅰ型中心裂纹的动态SCB实验，当裂纹相对长度在0.1~0.4范围内时，可以忽略惯性效应，直接采用静态KI公式计算裂尖的最大断裂韧性，相对误差约10%。对于含复合型斜裂纹的动态SCB实验，当裂纹相对长度在0.2~0.4，裂纹倾角在10°~40°范围内时，也可以直接采用静态公式计算裂尖的最大断裂韧性，相对误差小于10%。

图 1 Steven试验装置结构示意图

Figure 1. Configuration for Steven test

下载: 全尺寸图片幻灯片

图 2 计算模型的初始图像

Figure 2. Initial configuration of simulation model

下载: 全尺寸图片幻灯片

图 3 被弹丸撞击后靶板的计算图像

Figure 3. Deformed state of target impacted by projectile

下载: 全尺寸图片幻灯片

图 4 样品盒背板中心应变的对比

Figure 4. Comparison of strain at center of holder's back plate

下载: 全尺寸图片幻灯片

图 5 典型时刻温度云图对比

Figure 5. Comparison of temperature profile at typical time

下载: 全尺寸图片幻灯片

图 6 典型时刻压力云图对比

Figure 6. Comparison of pressure profile at typical time

下载: 全尺寸图片幻灯片

图 7 不同撞击速度下点火位置的温升

Figure 7. Temperature variation of ignition point at different impact velocities

下载: 全尺寸图片幻灯片

表 1 炸药点火反应的速度阈值与实验数据比较

Table 1. Predicted threshold of impact velosity inducing explosive ignisioncompared experimental data

文献	v/(m·s^-1)
文献	实验^[4]	计算^[4]	本文计算
[3]	54~56	44~45	45~48
[2]	43~53	45~50	45~49

下载: 导出CSV

参考文献(19)

[1]	Chidester S K, Green L G, Lee C G. A frictional work predictive method for the initiation of solid high explosives from low-pressure impacts[C]//Proceeding of 10th International Detonation Symposium. 1993: 785-792.
[2]	Chidester S K, Tarver C M, Garza R. Low amplitude impact testing and analysis of pristine and aged solid high explosives[C]//Proceeding of 11th International Detonation Symposium. 1998: 93-100.
[3]	Idar D J, Lucht R A, Straight J W, et al. Low amplitude insult project: PBX 9501 high explosive violent reaction experiments[C]//Proceeding of 11th International Detonation Symposium. 1998: 102-110.
[4]	Scammon R J, Browning R V, Middleditch J, et al. Low amplitude insult project: Structural analysis and prediction of low order reaction[C]//Proceeding of 11th International Detonation Symposium. 1998: 111-118.
[5]	Vandersall K S, Chidester S K, Forbes J W, et al. Experimental and modeling studies of crush, puncture, and perforation scenarios in the steven impact test[C]//Proceeding of 12th International Detonation Symposium. 2002: 131-139.
[6]	Wortley S, Jones A, Cartwright M, et al. Low speed impact of pristine and aged solid high explosive[C]//Proceeding of 12th International Detonation Symposium. 2002: 399-408.
[7]	Switzer L L, Vandersall K S, Chidester S K, et al. Threshold studies of heated HMX-based energetic material targets using the Steven impact test[C]//Proceeding of Shock Compression of Condensed Matter. 2003: 1045-1048.
[8]	Lee E L, Tarver C M. Phenomenological model of shock initiation in heterogeneous explosives[J]. Physics of Fluids, 1980, 23(12):2362-2372. doi: 10.1063/1.862940
[9]	Murphy M J, Lee E L, Weston A M, et al. Modeling shock initiation in composition B[C]//Proceeding of the 10th International Detonation Symposium. 1993: 786-792.
[10]	Vandersall K S, Tarver C M, Garcia F, et al. Shock initiation experiments on PBX9501 explosive at 150℃ for ignition and growth modeling[C]//Proceeding of Shock Compression of Condensed Matter. 2005: 1127-1130.
[11]	Tarver C M, Lefrancois A S, Lee R S, et al. Shock initiation of the PETN-based explosive LX-16[C]//Proceeding of 13th International Detonation Symposium. 2006: 139.
[12]	Vandersall K S, Tarver C M, Garcia F, et al. Low amplitude single and multiple shock initiation experiments and modeling of LX-04[C]//Proceeding of 13th International Detonation Symposium. 2006: 145.
[13]	Urtiew P A, Vandersall K S, Tarver C M, et al. Shock initiation experiments and modeling of composition B and C-4[C]//Proceeding of 13th International Detonation Symposium. 2006: 147.
[14]	Tarver C M, Chidester S K. Ignition and growth modeling of detonating TATB cones and arcs[C]//Proceeding of Shock Compression of Condensed Matter. 2007: 429-432.
[15]	Chidester S K, Garcia F, Vandersall K S, et al. Shock initiation experiments plus ignition and growth modeling of damaged LX-04 charges[C]//Proceeding of Shock Compression of Condensed Matter. 2009: 271-274.
[16]	May C M, Tarver C M. Modeling short shock pulse duration initiation of LX-16 and LX-10 charges[C]//Proceeding of Shock Compression of Condensed Matter. 2009: 275-278.
[17]	Vandersall K S, Tarver C M, Garcia F, et al. On the low pressure shock initiation of octahydro-1, 3, 5, 7-tetranitro-1, 3, 5, 7-tetrazocine based plastic bonded explosives[J]. Journal of Applied Physics, 2010, 107(9):094906. doi: 10.1063/1.3407570
[18]	Hallquist J O. LS-DYNA Theoretical manual[Z]. Livemore Software Technology Corporation, 1998.
[19]	Livemore Software Technology Corporation. LS-DYNA Keyword user's manual: Version 970[Z]. 2003.

施引文献

期刊类型引用(5)

1.	张晓婧，郁红陶. 含细观缺陷固体推进剂非冲击点火性能研究进展. 兵器装备工程学报. 2024(07): 217-223 . 百度学术
2.	周静静，李子涵，黄蒙，祝艳龙，丁黎，常海. 炸药热-力耦合响应机制研究进展. 固体火箭技术. 2024(04): 446-455 . 百度学术
3.	Hai-jiao Xue，Yan-qing Wu，Kun Yang，Yi Wu. Microcrack-and microvoid-related impact damage and ignition responses for HMX-based polymer-bonded explosives at high temperature. Defence Technology. 2022(09): 1602-1621 . 必应学术
4.	张琪林，段卓平，孟凡星，南海，王昕捷，黄风雷. 浇注炸药PBX-1侵彻安定性试验与数值模拟. 含能材料. 2021(02): 107-113 . 百度学术
5.	黄彬彬，傅华，喻寅，刘仓理. 基于有限元-离散元结合方法的Steven实验三维数值模拟. 含能材料. 2020(10): 995-1002 . 百度学术