高溫下奧克托今單晶沖擊響應數值計算

丁凱，王昕捷，黃亨建，吳艷青，黃風雷

(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.