高温下奥克托今单晶冲击响应数值计算

丁凯，王昕捷，黄亨建，吴艳青，黄风雷

(1.北京理工大学 爆炸科学与技术国家重点实验室，北京 100081；2.中国工程物理研究院 化工材料研究所，四川 绵阳 621999)

0 引言

高聚物粘结炸药(PBX)的点火起爆始终是评价其安全性和应用其良好爆轰性能的核心问题。奥克托今(HMX)单晶作为PBX的主要成分，其在高温下的冲击响应以及细观变形机制，极大地影响PBX在复杂环境下的冲击响应及点火起爆特性。

炸药单晶的点火起爆与其动态响应和力学变形机制相关。Dick等[1-5]对太恩(PETN)及HMX等单晶进行一系列平板撞击实验，研究表明炸药单晶的冲击感度呈现明显的各向异性，并指出基于位错滑移的塑性是单晶变形的主要机制。Winey等[6]发展了考虑非线性弹性的位错模型，可描述PETN单晶的弹塑性冲击响应。Barton等[7]发展了基于位错的模型表征HMX单晶的冲击压缩响应，模型同时考虑了热激活以及声子拖曳两种位错运动机制，通过与平板撞击实验中粒子速度历史曲线对比标定了模型参数。Wu等[8]发展了动高压晶体热力学本构模型，引入状态方程来描述高压条件下炸药单晶非线性变形，同时考虑弹性模量的压力和温度相关性，分析不同取向HMX单晶的滑移塑性应变分布以及位错密度的变化。Wang等[9-10]对不同取向的HMX和黑索今(RDX)单晶进行常温下的平板撞击实验，并构建各向异性晶体位错塑性模型，可较好地模拟常温下单晶冲击热力学响应。以上发展的本构模型大多忽略温度对弹性及位错塑性的影响，很难用于预测高温下单晶各向异性冲击热力学响应，同时模型包含大量参数，不适应后续开展工程计算的需要。

利用位错塑性模型研究高温下单晶的冲击响应最早开始于非含能单晶。Krasnikov等[11]建立了考虑热激活和声子拖曳的位错运动模型，通过数值计算再现了受冲击单晶铝Hugoniot弹性极限(HEL)的热硬化现象，并认为声子拖曳位错机制为主导原因。Gurrutxaga-Lerma等[12]通过动态离散位错动力学模型模拟发现，冲击波作用下弹性先驱波的演化主要受波阵面上位错形核产生的辐射阻尼影响，辐射阻尼与位错运动速度呈负相关。金属铝波阵面上位错运动速度趋近于剪切波速，而剪切波速受剪切模量控制并且随温度升高而减小，因此辐射阻尼效应增强，导致HEL随温度升高而增大。Austin[13]通过建立的连续位错动力学模型研究冲击加载条件下金属铝的热硬化效应，研究表明弹性先驱波的演化受声子散射与辐射阻尼效应控制，二者随着温度升高而增强。目前，关于高应变率条件下金属晶体屈服强度的温度相关性研究较多[14-16]，而炸药单晶在热力耦合作用下本构模型的发展由于相关实验开展的难度较大而发展缓慢。

因此，本文发展基于热激活和声子拖曳位错滑移机制的非线性热弹黏塑性模型，并基于前期开展的高温(373 K、423 K)下HMX单晶平板撞击实验，数值研究相变温度以下HMX单晶冲击热力学响应及相应细观变形机制。

1 热弹黏塑性模型

针对高温、高应变率的加载条件，基于Austin等[17]的位错模型，发展考虑高压下晶体的非线性热弹性，及基于热激活和声子拖曳位错机制的热弹黏塑性模型。

1.1 热弹性

根据广义胡克定律，热弹性本构模型[18]可表示为

(1)

(2)

式中：E0为T0=300 K时的弹性模量；θ为体积压缩比，θ=V/Vi，V和Vi分别为比容和初始比容。另外，将应力中球量部分以Grüneisen状态方程pe替代:

(3)

式中：K为体积模量；s为应力偏量部分。pe由Hugoniot压力的Taylor展开[20]得到：

(4)

式中：ρ0为初始密度；c0为体积声速；s为冲击波波速D和波后质点速度u之间线性关系的斜率；η为体积相对变化，η=1-θ.热弹性模型参数如表1所示。

表1 热弹性模型参数[9，21-22]Tab.1 Parameters for thermoelasticity model[9，21-22]

1.2 位错塑性

基于位错滑移的塑性本构模型，假设位错滑移满足Schmid定律，剪应力τ[23]可表示为

(5)

(6)

(7)

(8)

式中：k为玻尔兹曼常数；υe为有效越过障碍的频率；T为当前温度；ΔG为热激活自由能，可表示为

ΔG=gtμb3[1-(τ/τ0)P]Q，

(9)

gt为正则化的总自由能系数，P和Q描述障碍分布[26]，τ0为声子拖曳机制控制位错滑移的临界剪应力。

(10)

(11)

式中：BT和Tm为归一化因子；ak为拟合系数。

根据Taylor关系[28]，τ0可表示为

(12)

式中：τp0为0 K时的晶格阻力；α表征位错纠缠的长程相互作用系数；ρt为总位错密度。

τμ可表示为

(13)

式中：β为比例常数；μ0为0 K时材料的剪切模量。

根据位错理论，位错可分为可移动和不可移动位错两部分，可移动位错密度ρm占总位错密度ρt的比例分数f[29]可表示为

f=exp(-Hγp/τ)，

(14)

式中：H为位错硬化系数；γp为滑移剪应变。

(15)

(16)

图2 计算流程图Fig.2 Computational flowchart

表2 位错塑性模型参数Tab.2 Parameters for dislocation plasticity model

图1 模型参数敏感性分析曲线Fig.1 Analytical curves of model parameter sensitivity

2 数值模型

2.1 平板撞击实验

表3 不同初温下HMX单晶的实验结果[9，34]Tab.3 Experimental results of HMX at different initial temperatures[9,34]

2.2 有限元模型

利用ABAQUS有限元软件建立HMX单晶平板撞击实验的二维有限元模型如图3所示，并使用VUMAT定义HMX单晶本构行为。飞片、砧板和窗口均使用Mie-Grüneisen状态方程描述，材料参数如表4所示。2024铝飞片尺寸为8 mm×20 mm，单元尺寸为0.01 mm并划分为800个网格。石英尺寸为4 mm×20 mm，单元尺寸为0.008 mm并划分为500个网格。HMX单晶尺寸为1 mm×20 mm，单元尺寸为0.005 mm并划分为200个网格。LiF窗口尺寸为4 mm×20 mm，单元尺寸为0.01 mm并划分为400个网格。通过前期网格敏感性分析表明，在此网格密度下能得到精确的计算结果。单元类型为四节点双线性减缩积分单元(CPE4R)。在分析步中设置线性体积黏性参数为0.15，从而避免数值振荡。对2024铝飞片整体施加初始速度，单晶的环境温度分别设置为实验初始温度300 K、373 K和423 K.

表4 飞片及靶板系统材料参数Tab.4 Parameters for flyer and target materials

图3 二维有限元模型Fig.3 Two-dimensional finite element model

3 结果和讨论

为更好地对比不同初温下实验和计算得到的HMX/LiF界面粒子速度历史，将初温373 K和423 K的实验和计算结果分别向右平移0.1 μs和0.4 μs，如图4所示为平移后不同初温下界面粒子速度历史的实验与计算结果对比。不同初温下界面粒子速度曲线第1个波峰对应的粒子速度为HEL对应的界面粒子速度。从图4中可以看出，不同初温下计算得到的弹塑性双波结构与实验结果吻合较好。对于300 K，模型能够很好地描述弹性先驱波到达界面以及随后塑性波的到达，但计算出的HEL略低于实验结果，这可能是由于较快的位错速度或较高的位错密度，导致相应塑性松弛速率偏高。而对于373 K，计算得到的HEL和塑性波的到达均与实验结果一致。对于423 K，计算出的HEL与实验结果相同，但随后的应力松弛更为明显，这是由于位错的扩展和形核过程与实验存在一定差异。另外，该模型能较好地描述初始温度较高时HEL呈现的热硬化效应。

图4 不同初温下HMX/LiF界面粒子速度历史的实验和计算结果对比Fig.4 Comparison of experimental and calculated particle velocities at interface of HMX/LiF at different initial temperatures

3.1 热硬化影响机制

图5 不同初温下剪应力τ与塑性剪应变率的关系Fig.5 Relationship between resolved shear stress τ and plastic shear strain rate at different initial temperatures

图6表示不同初温下界面处可移动位错密度ρm以及热激活临界剪应力τμ历史曲线。在模型中考虑了位错产生的两种机制，即已有位错环的扩展以及新位错的形核，因此当剪应力达到热激活临界剪应力τμ，会引起位错的增殖和运动。在图6(a)中300 K、373 K和423 K对应HEL的可移动位错密度依次为3.40×1010m-2、4.52×1010m-2和4.80×1010m-2，由(6)式可知随着可移动位错密度增大，塑性剪应变率增大，HEL会减小，但HEL表现出热硬化，表明位错演化的温度相关性并不对热硬化效应起决定性作用。由(13)式可知，可移动位错密度增大会导致τμ增大，但由图6(b)可知τμ整体随初温升高而减小，这是因为随初温升高、剪切模量减小的程度大于可移动位错密度增大的程度。同时根据(12)式可知声子拖曳临界剪应力τ0也减小。由于临界剪应力τμ和τ0的热软化效应，位错更容易启动进入热激活状态，以及更容易从热激活状态转变为声子拖曳状态。因此，位错演化并不直接导致单晶HEL的热硬化效应，剪切模量的温度相关性决定临界剪应力的大小，从而影响位错滑移机制的演变。

图6 不同初温下界面处可移动位错密度ρm以及热激活临界剪应力τμ历史曲线Fig.6 Mobile dislocation density ρm and thermal threshold shear stress τμ at different initial temperatures

图7 不同初温下界面处平均位错速度以及剪应力τ历史曲线Fig.7 Mean dislocation velocity and resolved shear stress τ at different initial temperatures

声子拖曳系数B随初温升高而增大与声子散射系数B0以及剪切模量μ的温度效应有关，因此，HEL热硬化效应的实质来源于由声子散射和辐射阻尼控制的位错运动。为量化声子散射和辐射阻尼对热硬化效应的影响，计算获得不同初温下考虑声子散射、辐射阻尼以及二者共同作用时的HEL相对值如图8所示。图8中HEL相对值表示相对300 K时HEL的大小。基准线是在B0(T0)和∂μ/∂T=0条件下计算得到。当只考虑辐射阻尼对热硬化效应的影响时，声子散射系数B0无温度相关性；当只考虑声子散射对热硬化效应的影响时，剪切模量μ无温度相关性；当二者共同作用时，需同时考虑B0和μ的温度相关性。从图8可以看出，声子散射和辐射阻尼效应随着初温升高而增强，但由于初始温度对剪切模量的影响较小，导致辐射阻尼效应对HEL热硬化效应的贡献小于声子散射效应。

图8 不同初温下考虑声子散射、辐射阻尼以及二者共同作用时HEL相对值Fig.8 Relative value of HEL provided separately by phonon scattering,radiation damping and their combined effect at different initial temperatures

3.2 位错滑移机制的演变

图9 不同初温下与τ/τ0的关系Fig.9 Relationship between and τ/τ0 at different initial temperatures

3.3 单晶的热力学响应

图10(a)为不同初温下界面轴向应力-s11历史曲线，从中可以看出，300 K、373 K和423 K的动态屈服极限依次为0.904 GPa、1.529 GPa及1.722 GPa，随着初温升高动态屈服极限增大。平台阶段轴向应力差别不大，无明显温度相关性。图10(b)为423 K时单晶不同位置的轴向应力-s11历史曲线，从中可以看出，由于波的传播，单晶内部形成弹性波及塑性波，随着离左端面的距离增大，弹性波在传播过程中不断衰减，使动态屈服极限依次减小，但当弹性波到达单晶右端面时在界面处发生反射使动态屈服极限突然增大。

图10 不同初温下界面轴向应力-s11和423 K时单晶不同位置的轴向应力-s11历史曲线Fig.10 Longitudinal stresses -s11 at different initial temperatures and different distances at 423 K

HMX单晶受到冲击时的温升主要来源于冲击压缩体积功以及位错滑移塑性功，当应力水平较高时，会直接引起较高的温升。图11为不同初温下界面处温升历史曲线，从中可以看出，随着初温升高，HEL对应温升逐渐增大，而平台阶段由于应力差别不大，对应温升无明显温度相关性。

图11 不同初温下界面处温升历史曲线Fig.11 Temperature rise curves at interface at different initial temperatures

图12为初温373 K下不同初始可移动位错密度对位错演化的影响。从图12可以看出，塑性波到达界面之前，可移动位错密度为初始可移动位错密度，随着塑性波传播，可移动位错密度增大，从图7(b)看出最后剪应力小于图6(b)对应热激活临界剪应力，位错最后停止增殖达到饱和状态。当初始可移动位错密度较高时，位错在较短时间内达到饱和位错密度，导致图1(g)中应力松弛时间和塑性波上升时间缩短。同时，从图1(g)可以看出，较低的初始可移动位错密度对应HEL较高，这是由于较低初始可移动位错密度产生较低塑性剪应变率，表现出较高HEL.

图12 初温373 K下不同初始可移动位错密度对位错演化的影响Fig.12 Effect of different initial mobile dislocation densities on dislocation evolution at 373 K

4 结论

本文发展了考虑热激活和声子拖曳位错滑移机制的非线性热弹黏塑性模型，基于前期开展的高温(373 K、423 K)下HMX单晶平板撞击实验，可数值再现HMX单晶HEL的热硬化效应。通过定量分析声子散射和辐射阻尼对热硬化效应的影响并研究位错滑移机制演变和热力学响应。得到主要结论如下：

1) 声子散射和辐射阻尼随着初温升高而增大导致声子拖曳系数增大，使可移动位错黏性摩擦增强。此时，平均位错速度由2 237 m/s减小至1 537 m/s，进而产生较低的塑性剪应变率和较高的流动应力，引起HMX单晶HEL的热硬化效应。

2) 在高温(373 K、423 K)下，仅考虑辐射阻尼时对应HEL相对值分别为1.068和1.112，由于剪切模量随着初始温度升高变化较小，导致辐射阻尼对应HEL相对值均小于声子散射对应HEL相对值(1.39、1.63)。因此，辐射阻尼对HEL热硬化效应的贡献小于声子散射。

参考文献(References)

[1] DICK J J.Effect of crystal orientation on shock initiation sensitivity of pentaerythritol tetranitrate explosive[J].Applied Physics Letters,1984,44(9):859-861.

[2] DICK J J,MULFORD R N,SPENCER W J,et al.Shock response of pentaerythritol tetranitrate single crystals[J].Journal of Applied Physics,1991,70(7):3572-3587.

[3] DICK J,RITCHIE J.Molecular mechanics modeling of shear and the crystal orientation dependence of the elastic precursor shock strength in pentaerythritol tetranitrate[J].Journal of Applied Physics,1994,76(5):2726-2737.

[4] DICK J J,MARTINEZ A R.Elastic precursor decay in HMX explosive crystals[J].AIP Conference Proceedings,2002,620:817-820.

[5] DICK J J,HOOKS D E,MENIKOFF R,et al.Elastic-plastic wave profiles in cyclotetramethylene tetranitramine crystals[J].Journal of Applied Physics,2004,96(1):374-379.

[6] WINEY J M,GUPTA Y M.Anisotropic material model and wave propagation simulations for shocked pentaerythritol tetranitrate single crystals[J].Journal of Applied Physics,2010,107(10):103505.

[7] BARTON N R,WINTER N W,REAUGH J E.Defect evolution and pore collapse in crystalline energetic materials[J].Modelling and Simulation in Materials Science and Engineering,2009,17(3):035003.

[8] WU Y Q,HUANG F L.Thermal mechanical anisotropic constitutive model and numerical simulations for shocked β-HMX single crystals[J].European Journal of Mechanics A/Solids,2012,36:66-82.

[9] WANG X J,WU Y Q,HUANG F L,et al.Dynamic anisotropic response ofβ-HMX andα-RDX single crystals using plate impact experiments at ～1 GPa[J].Propellants,Explosives,Pyrotechnics,2018,43(8):759-770.

[10] WANG X J,WU Y Q,HUANG F L,et al.An anisotropic elastoviscoplasticity model of thermomechanical responses of shockedβ-HMX andα-RDX single crystals[J].Propellants,Explosives,Pyrotechnics,2019,44(7):870-888.

[11] KRASNIKOV V S,MAYER A E,YALOVETS A P.Dislocation based high-rate plasticity model and its application to plate-impact and ultra short electron irradiation simulations[J].International Journal of Plasticity,2011,27(8):1294-1308.

[12] GURRUTXAGA-LERMA B,SHEHADEH A M,BALINT D S,et al.The effect of temperature on the elastic precursor decay in shock loaded FCC aluminum and BCC iron[J].International Journal of Plasticity,2017,96:135-155.

[13] AUSTIN R.Elastic precursor wave decay in shock-compressed aluminum over a wide range of temperature[J].Journal of Applied Physics,2018,123(3):035103-1-035103-14.

[14] YAO S L,PEI X Y,LIU Z L,et al.Numerical investigation of the temperature dependence of dynamic yield stress of typical BCC metals under shock loading with a dislocation-based constitutive model[J].Mechanics of Materials,2019,140:103211.

[15] KOSITSKI R,MORDEHAI D.On the origin of the stress spike decay in the elastic precursor in shocked metals[J].Journal of Applied Physics,2019,126(8):085901.

[16] KANEL G I,SAVINYKH A S,GARKUSHIN G V,et al.Effects of temperature on the flow stress of aluminum in shock waves and rarefaction waves[J].Journal of Applied Physics,2020,127(3):035901.

[17] AUSTIN R A,MCDOWELL D L.A dislocation-based constitutive model for viscoplastic deformation of fcc metals at very high strain rates[J].International Journal of Plasticity,2011,27(1):1-24.

[18] LUSCHER D J,BUECHLER M A,WALTERS D J,et al.On computing the evolution of temperature for materials under dynamic loading[J].International Journal of Plasticity,2018,111:188-210.

[19] STEINBERG D J,COCHRAN S G,GUINAN M W,et al.A constitutive model for metals applicable at high-strain rate[J].Journal of Applied Physics,1980,51(3):1498-1504.

[20] LUKYANOV A A.Constitutive behaviour of anisotropic materials under shock loading[J].International Journal of Plasticity,2008,24(1):140-167.

[21] PALMER S J P,FIELDJ E.The deformation and fracture ofβ-HMX[J].Proceedings of the Royal Society of London .Series A,Mathematical and Physical Sciences,1982,383(1785):399-407.

[22] CUI H L,JI G F,CHEN X R,et al.Phase transitions and mechanical properties of octahydro-1，3，5，7-tetranitro-1，3，5，7-tetrazocine in different crystal phases by molecular dynamics simulation[J].Journal of Chemical &Engineering Data,2010,55(9):3121-3129.

[23] BISHOP J F W,HILL R.XLVI.A theory of the plastic distortion of a polycrystalline aggregate under combined stresses[J].The London,Edinburgh and Dublin Philosophical Magazine and Journal of Science,1951,42(327):414-427.

[24] GILMAN J J,JOHNSTON W G.Dislocations in lithium fluoride crystals[J].Solid State Physics,1962,13:147-222.

[25] FOLLANSBEE P S,KOCKS U F.A constitutive description of the deformation of copper based on the use of the mechanical threshold stress as an internal state variable[J].Acta Metallurgica,1988,36(1):81-93.

[26] KOCKS U F,MECKING H.Physics and phenomenology of strain hardening:the FCC case[J].Progress in Materials Science,2003,48(3):171-273.

[27] KOCKS U F,ARGON A S,ASHBY M F.Thermodynamics and kinetics of slip[M]∥Chalmers B,Christian J W,Massalski T B.Progress in Materials Science.Oxford,UK:Pergamon Press,1975:141-145.

[28] TAYLOR G I.The mechanism of plastic deformation of crystals.Part Ⅱ.Comparison with observations[J].Proceedings of the Royal Society.Series A,Containing Papers of a Mathematical and Physical Character,1934,145(855):388-404.

[29] WINEY J M,GUPTA Y M.Nonlinear anisotropic description for the thermomechanical response of shocked single crystals:inelastic deformation[J].Journal of Applied Physics,2006,99(2):023510.

[30] GUPTA Y M,DUVALL G E,FOWLES G R． Dislocation me-chanisms for stress relaxation in shocked LiF[J].Journal of Applied Physics,1975,46(2):532-546.

[31] HULL D,BACON D J.Introduction to dislocations[J].Materials Today,2011,19(12):91-92.

[32] CLIFTON R J.On the analysis of elastic/visco-plastic waves of finite uniaxial strain[C]∥Proceedings of Sagamore Army Materials Research Conference.Syracuse,NY,US:Syracuse University Press,1971:73-116.

[33] ARMSTRONG R W,ELBAN W L.Materials science and technology aspects of energetic (explosive) materials[J].Materials Science and Technology,2013,22(4):381-395.

[34] 丁凯,王昕捷,吴艳青,等.高温下奥克托今单晶的冲击响应特性[J].兵工学报,2020,41(9):1783-1791.

DING K,WANG X J,WU Y Q,et al.Shock response of HMX single crystal at elevated temperatures[J].Acta Armamentarii,2020,41(9):1783-1791.(in Chinese)

[35] MARSH S P.LASL shock Hugoniot data[M].Berkeley,CA,US:University of California Press,1980:165-166.

[36] BARRIOS M A,BOEHLY T R,HICKS D G,et al.Precision equation-of-state measurements on National Ignition Facility ablator materials from 1 to 12 Mbar using laser-driven shock waves[J].Journal of Applied Physics,2012,111(9):093515.

[37] JONES S C,GUPTA Y M.Refractive index and elastic properties of z-cut quartz shocked to 60 kbar[J].Journal of Applied Physics,2000,88(10):5671-5679.

[38] AO T,KNUDSON M D,ASAY J R,et al.Strength of lithium fluoride under shockless compression to 114 GPa[J].Journal of Applied Physics,2009,106(10):103507

[39] LALONE B M,FATYANOV O V,ASAY J R,et al.Velocity correction and refractive index changes for [100] lithium fluoride optical windows under shock compression,recompression,and unloading[J].Journal of Applied Physics,2008,103(9):093505.