1 引言

在高温高压下, 含能材料会发生相变, 从而影响物理和化学性质。奥克托今(HMX)是一种比较常用的固态含能材料, 有四种相态, 记为α, β, 水合物γ, δ, 四种相态的稳定性是β＞α＞γ＞δ^[1], 这同时是四种相态的密度减小的顺序。常温下最稳定的相态是单斜晶β相^[2], 高温下最稳定的是六方晶δ相。HMX化学反应包括四个反应步, 其中β→δ相变是第一步^[3]。

关于HMX相变的实验和分子动力学研究有很多。研究HMX相变的主要实验方法有:拉曼光谱, 热重分析, 傅里叶法(FTIR), 二次谐波(SHG), 差式扫描量热法(DSC), X射线衍射法(XRD), 金刚石压腔法(DAC)。Henson等^[4-5]使用SHG技术研究含能材料的相变动力学, 根据SHG信号的强度变化对两相进行了定量分析, 并基于HMX的β→δ相变的全过程, 提出了不可逆的二级反应的假设, 建立了研究相变动力学的模型。Smilowitz等^[6]通过实验发现, 当将烤燃后的HMX样品温度降低到δ相稳定温度下时, δ-HMX能完全变回β-HMX, 并且晶体表面出现裂纹。烤燃过程中, 温度梯度将产生压力梯度, 从而导致相变和化学分解, 并影响点火时间和位置^[7]。陈朗等^[8]对含HMX的炸药进行了不同升温速率下的多点测温烤燃实验, 并使用多步反应动力学模型模拟化学反应产物(β-HMX, δ-HMX, 最终气体产物)的演化过程。他们发现HMX的相变对炸药温度有很大影响, 从而影响炸药的点火。龙和陈^[9]通过建立声子-电子自由能模型研究HMX的相变和热动力学性质, 在较大的温度和压力范围内, 得到了不同相HMX的体积模量, 热膨胀系数, Hugoniot曲线和相变曲线。周婷婷等^[10]通过计算HMX不同晶向的分解剪应力, 得到了三个冲击面的感度为(0 1 0)＞(1 1 0)＞(0 1 1), 因此对于敏感的(0 1 0)平面, 化学反应进行得更快, 与较不敏感的(1 1 0)和(0 1 1)平面相比, 在更短的时间内可形成更多的产物。目前相变主要是通过实验^[4-8]研究相变机理, 从而得到不同温度和压力条件下的HMX相变历史, 以及粒度和添加剂等对HMX相变的影响。大量分子动力学^[9-10]计算可得到HMX晶体不同相的热力学参数。对于HMX相变在力学变形中的影响, 国内外的研究较少。

为此, 本研究在晶体塑性本构的基础上, 发展了考虑相变的细观本构模型, 并研究了HMX单晶在烤燃作用下的各向异性力学响应及相变。通过对HMX晶体的烤燃模拟, 描述了HMX晶体的可逆相变, 及不同升温速率下β→δ相变对炸药晶体变形和感度的影响, 为更好地预测HMX点火奠定基础。

2 模型部分

本构模型是基于吴和黄^[11-12]的模型发展的, 用于描述Dick等^[13]的HMX平板撞击实验和高聚物粘结炸药(PBX)的热力学撞击响应。模型考虑了晶体塑性^[14]和相变^[5], 用于模拟烤燃下HMX的相变。

2.1 基于位错理论的晶体塑性本构

基于乘法分解和极分解, 变形梯度F可表示为:

$\begin{eqnarray} \boldsymbol{F}=\boldsymbol{F}^{\rm{e}}\boldsymbol{F}^{\rm{p}}=\boldsymbol{R}^{\rm{e}}\boldsymbol{U}^{\rm{e}}\boldsymbol{F}^{\rm{p}} \end{eqnarray}$

(1)

式中, F^e表示弹性变形, F^p表示塑性变形, F^e可分解为R^e和U^e, 分别表示热弹性晶格右拉伸张量和晶格旋转张量。

通过(1), 速度梯度L可表示为:

$\begin{eqnarray} \boldsymbol{L}=\dot{\boldsymbol{F}}\boldsymbol{F}^{-1}=\mathit{\boldsymbol{Ω}}+\boldsymbol{R}^{\rm{e}}\dot{\boldsymbol{L}}\boldsymbol{R}^{\rm{e}^{-1}}; \mathit{\boldsymbol{Ω}}=\dot{\boldsymbol{R}}^{\rm{e}}\boldsymbol{R}^{\rm{e}^{-}1} \end{eqnarray}$

(2)

$\begin{eqnarray} \bar{\boldsymbol{L}}&=&\bar{\boldsymbol{L}}^{\rm{e}}+\bar{\boldsymbol{L}}^{\rm{p}}\\ &=&\dot{\boldsymbol{U}}^{\rm{e}}\boldsymbol{U}^{\rm{e}^{-1}}+\boldsymbol{U}^{\rm{e}}\dot{\boldsymbol{F}}^{p}\boldsymbol{F}^{\rm{p}^{-1}}\boldsymbol{U}^{\rm{e}^{-1}}\\ &=&\dot{\boldsymbol{U}}^{\rm{e}}\boldsymbol{U}^{\rm{e}^{-1}}+∑\limits^{n}_{α=1}(\bar{\boldsymbol{s}}^{(α)}⊗\bar{\boldsymbol{m}}^{(α)})\dot{γ}^{(α)} \end{eqnarray}$

(3)

式中, L表示无转动的速度梯度。本文中所有描述符号上方有“—”均与L描述相同。滑移系α的滑移方向为s(α), 滑移面法向为m(α), $\dot{γ}$^(α)表示滑移系α的剪应变速率。

通过滑移系描述晶体塑性本构方程, 基于Schmid定律的分解剪应力可表示为:

$\begin{eqnarray} &&\boldsymbol{τ}^{(α)}=\bar{\boldsymbol{P}}^{(α)}:\bf{σ}^{L};\bar{\boldsymbol{P}}^{(α)}=\frac{1}{2}\left[\bar{\boldsymbol{s}}^{(α)}\bar{\boldsymbol{m}}^{(α)}+\bar{\boldsymbol{m}}^{(α)}\bar{\boldsymbol{s}}^{(α)}\right];\\ &&\bf{σ}^{L}=\boldsymbol{R}^{\rm{e}^{T}}σ\boldsymbol{R}^{\rm{e}} \end{eqnarray}$

(4)

式中, P(α)为Schimid因子, σ^L表示柯西应力。

Orowan方程^[37]可以描述率相关晶体塑性的剪应变:

$\begin{eqnarray} \dot{γ}^{(α)}=ρ^{\left(α\right)}_{m}b^{(α)}\bar{v}^{(α)} \end{eqnarray}$

(5)

式中, b^(α)为柏氏向量, $ρ^{(α)}_{m}$为位错密度, v^(α)为滑移系α的平均位错速度。位错速度被认为与分解剪应力有关^[15]:

$\begin{eqnarray} \bar{v}^{(α)}=v^{(α)}_{0}\rm{exp}\left[\frac{-s_{d}}{(τ^{(α)}-τ_{0})}\right] \end{eqnarray}$

(6)

式中, $v^{(α)}_{0}$为剪切波速, s_d为牵引力, τ₀为位错运动的临界剪应力。位错速度小于相应剪切滑移平面处的剪切波速度。

为描述位错密度的演化规律, 考虑位错产生的两种模型:已有位错环的扩展和新位错的形核^[16]。应变率较低时, 位错密度的增大主要是由于多交叉位错滑动引起已有位错环的增值^[17]; 应变率较高时, 较高的剪应力集中会引起位错形核^[18]。位错密度速率为:

$\begin{eqnarray} &&\dot{ρ}^{(α)}_{m}=M\dot{γ}^{(α)}\rm{exp}(\frac{-H\bar{γ}}{τ^{(α)}})+A(τ-τ_{c})\dot{τ};\\ &&\bar{γ}=∑\limits^{n}_{α=1}∫^{t}_{0}\left|\dot{γ}^{(α)}\right|{\rm{d}}t \end{eqnarray}$

(7)

式中, γ表示每个滑移系的累积剪应变; τ_c是位错形核的临界剪应力, M是位错增殖系数, H为位错硬化系数, A为位错形核系数。临界剪应力τ₀和τ_c受温度影响, 当温度升高时, 临界剪应力减小。

2.2 非线性热弹性

非线性热弹性模型包括各向异性弹性和热效应, 考虑有限体积变化, 用对数应变测量法描述弹性^[19]:

$\begin{eqnarray} \bar{ε}^{\rm{e}}=\rm{ln}\boldsymbol{U}^{\rm{e}};~~\dot{\bar{\boldsymbol{ε}}}^{\rm{e}}=\boldsymbol{D}^{\rm{e}} \end{eqnarray}$

(8)

弹性应变和熵可确定热动力学状态, 单位质量的内能随着弹性应变张量和熵增大^[20]。等熵条件下, 可以使用高阶弹性常数来近似描述比内能增量;

$\begin{eqnarray} ρ_{0}e(S,\bar{ε}^{\rm{e}})=ρ_{0}e(S,0)+\frac{1}{2}\bar{ε}^{\rm{e}}:\bar{\boldsymbol{C}}:\bar{\boldsymbol{ε}}^{\rm{e}} \end{eqnarray}$

(9)

式中, C(p, T)是压力和温度相关的二阶等熵弹性张量。无刚体转动下的柯西应力率$\dot{\bf{σ}}$^L可通过弹性应变率D^e得到:

$\begin{eqnarray} \bar{\dot{\bf{σ}}}^{L}=\bar{\boldsymbol{C}}:\bar{\boldsymbol{D}}^{\rm{e}}+\frac{1}{2}\bar{ε}^{\rm{e}}:\frac{∂\bar{\boldsymbol{C}}}{∂\bar{ε}^{\rm{e}}}:\bar{\boldsymbol{D}}^{\rm{e}}-ρT\mathit{Γ}\dot{S}; \mathit{Γ}=-\frac{1}{ρT}\frac{∂\bf{σ}^{L}}{∂S} \end{eqnarray}$

(10)

式中，Г为各向异性Grüneisen张量, $\dot{S}$为熵增加速率。

因此, 完整的率相关本构模型方程可表示为:

$\begin{eqnarray} &&\dot{\bf{σ}}^{L}=\bar{\boldsymbol{C}}:\bar{\boldsymbol{D}}^{\rm{e}}+\frac{1}{2}\bar{ε}^{\rm{e}}:\frac{∂\bar{\boldsymbol{C}}}{∂\bar{ε}^{\rm{e}}}:\bar{\boldsymbol{D}}^{\rm{e}}-ρT\mathit{Γ}\dot{S}-\\ &&∑\limits_{α}(\boldsymbol{W}(α)\bf{σ}^{L}-\bf{σ}^{L}\boldsymbol{W}(α)) \end{eqnarray}$

(11)

式中，W为反对称自旋张量。式(13)中需要单晶β-HMX的高阶弹性常数, 而目前获得的数据较少。假设β-HMX单晶的非线性弹性响应与温度和压力有关, $\frac{∂\bar{\boldsymbol{C}}}{∂\bar{ε}^{\rm{e}}}=\frac{∂\bar{\boldsymbol{C}}}{∂p}\frac{∂p}{∂\bar{ε}^{\rm{e}}}+\frac{∂\bar{\boldsymbol{C}}}{∂T}\frac{∂T}{∂\bar{ε}^{\rm{e}}}=\frac{∂\bar{\boldsymbol{C}}}{∂p}K+\frac{∂\bar{\boldsymbol{C}}}{∂T}α$, 其中, K为体积模量, α为热膨胀系数矩阵。体积模量K与温度和压力相关, 为了得到$\frac{∂\bar{\boldsymbol{C}}}{∂p}$和$\frac{∂\bar{\boldsymbol{C}}}{∂T}$, 崔等^[21]通过分子动力学对β-HMX在COMPASS力场下模拟, 得到了0 GPa, 5~555 K和298 K, 0~40 GPa的体积模量和弹性系数, 由此可得压力和温度相关的弹性系数。

相变过程中, 应力与β-HMX和δ-HMX百分含量有关:

$\begin{eqnarray} σ=λσ_{δ-HMX}+(1-λ)σ_{β-HMX} \end{eqnarray}$

(12)

为完整描述热动力学性质, 温度升高速率为:

$\begin{eqnarray} ρc_{V}\dot{T}=k_{c}(\frac{∂^{2}T}{∂x^{2}}+\frac{∂^{2}T}{∂y^{2}}+\frac{∂^{2}T}{∂z^{2}})+ρc_{V}T\mathit{Γ}:\dot{ε}{e}+\bar{σ}:\dot{ε}{p}+\dot{Q}_{β-δ} \end{eqnarray}$

(13)

式中, c_V为等容比热容, 与温度相关c_V=667.7+0.88T(J·kg^-1·K^-1) ^[22]。式(13)中第一项表示热传导过程产生的热量, 第二项表示由体积应变引起的体积功, 第三项表示塑性功。

模型中所需的参数如表 1和表 2所示。

2.3 HMX晶体的β→δ相变模型

奥克托今HMX单晶相变会经历两个过程:成核和增长。Henson等^[5]基于实验中HMX的β→δ相变全过程, 建立二阶相变反应动力学模型, 一阶描述成核过程, 二阶描述增长过程。该模型的公式与二次谐波SHG试验[11]所获得的数据能得到较好的吻合, 因此能较准确地预测分析HMX相变过程。HMX相变反应过程如下所示:

$\begin{eqnarray} &β-HMX\mathop{\rightleftharpoons}\limits^{k_{1}}_{k_{-1}}δ-HMX\\ &β-HMX+δ-HMX\mathop{\to}^{k_{2}}δ-HMX\\ &β-HMX+δ-HMX\mathop{\to}^{k_{-2}}β-HMX \end{eqnarray}$

(14)

式中, k₁, k₂, k_-1, k_-2是四个反应速率常数, 均与温度和压力相关, 表示如下:

$\begin{eqnarray} k_{i}=\frac{k_{B}T}{h}Q_{i}\left[\frac{TS_{i}-(U_{i}+pV_{i})}{RT}\right] \end{eqnarray}$

(15)

式中, k_B是玻尔兹曼常熟, h是普朗克常数, U_i, S_i, V_i分别表示内能、熵和体积, R是气体常数, Q_i是平衡常数, 与活化过渡态的分子浓度有关。因此, 相变速度可由式(14)得到:

$\begin{eqnarray} v_{ph}=k_{1}+\left[β_{0}(k_{2}-k_{-2})-(k_{1}+k_{-1})\right]x+β_{0}(k_{-2}-k_{2})x^{2} \end{eqnarray}$

(16)

式中, x=[δ]/β₀表示δ-HMX的摩尔分数, β₀为β-HMX摩尔密度, β₀=0.0064 mol·cm^-3。

HMX的β→δ相变过程是吸热过程, 潜热为ΔE=9.382 kJ·mol^-1, 式(13)中最后一项为相变吸热。

相变模型中的参数如表 3所示。

3 结果与讨论 3.1 模型建立

采用ABAQUS有限元软件建立1 mm×1 mm×1 mm HMX烤燃试验的三维模型如图 1所示, 并使用子程序VUMAT定义材料本构行为, 用ABAQUS/explicit计算模拟烤燃和静水压下HMX相变。HMX单晶的网格尺寸为50 μm, 总共有8000个单元, 网格类型为C3D8RT(三维六面体热耦合线性插值减少积分单元)。

3.2 相变的可逆性

对HMX晶体六个表面施加如图 2所示温度边界, 样品从300 K开始以2 K·min^-1升温, 当温度达到500 K时, 对样品以2 K·min^-1速率开始降温至300 K。

晶体的相变百分比, 即δ-HMX的摩尔分数随时间变化如图 3所示。当晶体发生β→δ相变时, 晶体密度由β相的1.90 g·cm^-3变为δ相的1.76 g·cm^-3, 会出现晶格体积膨胀, 从而产生裂纹, 而当晶体回到β相时, 裂纹不会消失。因此, 相变虽然可逆, 但路径不同。除此之外, 当晶体由δ相变回β相时, 会有残余应变, 如图 4所示。当晶体变回β相时, 三个主应变均不为零, 与初始时不同。这是由于HMX单晶在高温下的热膨胀是各向异性的, 并且δ相的热膨胀系数与δ相的热膨胀系数不同。

四个不同时刻的相变云图和从晶体边界到中心的温度梯度如图 5和图 6所示, 这四个时刻对应图 2中的A、B、C、D四点, 且图 6的四条曲线分别对应图 5中的四幅云图。由图 5可看出, 晶体的β→δ相变和δ→β相变均从边界开始, 逐渐向晶体内扩散。从图 6中可以看出, 在70 min时, 晶体最高温度为438 K, 在晶体边界处, 这时晶体边界开始发生β→δ相变。由于热传导, 晶体中心的温度会逐渐升高至超过相变温度, 故相变是由边界开始, 不断向中心扩展, 最终晶体完全转变为δ相。100 min后, 晶体外表面温度开始降低, 当温度低于相稳定温度时, 晶体开始发生逆相变, 随着热传导, 中心温度不断降低, 最终晶体重新回到β相。随着晶体表面温度的改变, 由于热传导, 晶体内会出现温度梯度, 温度梯度对相变有很大的影响。当温度梯度较大时, 晶体内部温度不均匀, 而相变依赖于温度, 晶体内的相变也不均匀。相变在晶体内的传播与温度的传播一致。

当升温速率不够高时, 相变不能完全进行。升温和降温速率为1.5 K·min^-1时, 相变百分比的演化如图 7所示。从图 7中可以看出, 当升温速率不够高时, 如果达到峰值温度后立刻降温, β相不能完全转化为δ相。当温度维持在峰值温度一段时间时, 由于热传导, 晶体内温度会逐渐均匀化, 相变百分比会增大, 而当维持时间足够长时, 相变能完全发生。由此可见, 温度均匀化有利于相变的发生。

3.3 不同升温速率下的β→δ相变

烤燃下的β→δ相变对HMX晶体变形有很大的影响, 这部分主要研究不同升温速率下, 相变对晶体变形的影响。对晶体的六个表面施加如图 8所示温度边界, 晶体分别以0.1, 0.2, 0.5, 1 K·min^-1的升温速率从300 K加热至500 K, 当温度达到500 K时保持恒温不变。

图 9a和图 9b分别描述了不同升温速率下相变百分比随时间和温度边界的演化历史。由图 9a可看出, 升温速率越快, 相变速度越快。从3.2节可知, 温度均匀化更利于相变的发生, 升温速率越快, 由于热传导, 晶体内温度梯度越大, 相变越不均匀, 而当晶体内部温度达到相变温度时, 晶体表面的温度会更高。因此, 升温速率越快, 相变完全时边界温度越高, 与图 9b的结果相符。

不同升温速率下HMX晶体的三个主应变的演化历史如图 10所示, 由于HMX单晶是各向异性的, 三个主应变的演化历史是不同的。晶体的总应变是由弹性应变加上塑性应变和热应变组成的, 烤燃过程中, 热应变起主要作用。升温速率越高, 晶体的应变越大。从图 10可看出, 当升温速率分别为0.1, 0.2, 0.5, 1 K·min^-1时, 应变曲线分别在1300, 650, 260, 120 min处出现明显拐点, 这是由于β相的热膨胀系数和δ相的不同, 而这些拐点对应了不同升温速率下的相变结束点。其中, ε₁₁和ε₃₃曲线的拐点比e₂₂更为明显, 且ε₁₁和ε₃₃曲线在拐点后是增大, 而ε₂₂曲线拐点后是略微减小, 这是由于β相的热膨胀系数在第一和第三主方向上小于δ相, 而第二主方向上略大于δ相。在晶体的第三主方向上, 晶体先处于压缩状态, 而相变后变为拉伸状态, 晶体在该方向上易产生裂纹。相变后晶体在第三方向长时间处于较大的拉应力状态下, 会发生损伤, 最终导致产生裂纹。由此可知, HMX的β→e相变会导致晶体更易产生裂纹。从图 10中还可看出, 最大主应变为ε₂₂, 最小主应变为ε₃₃, 而δ相的热膨胀系数α₂₂与α₃₃相同, 因此相变后的ε₂₂-ε₃₃应该保持恒定不变, 而相变前则是随着时间增大。不同升温速率下最大主应变和最小主应变之差如图 11所示, 从图 11中可看出, 最大主应变和最小主应变之差保持不变一段时间后会增大, 这是晶体温度升高导致的。而晶体温度升高主要是由于环境温度和晶体的体积功以及塑性功造成的, 图 11中曲线后半部分的升高主要是由于晶体内的塑性滑移系启动, 从而产生塑性功导致晶体温度升高。

不同升温速率下晶体的体积应变如图 12所示, 随着温度升高, 晶体呈现明显的体积膨胀。体积膨胀主要是由弹性膨胀和热膨胀导致, 以热膨胀为主。从图 12可看到明显的拐点, 这些拐点代表相变点, 相变后, 晶体体积应变增长更快, 这证明了相变会导致明显的体积膨胀。当晶体内的剪应力足够大时, 晶体的滑移系会启动, 产生塑性应变。不同升温速率下晶体的累积塑性剪应变的演化历史如图 13所示。由图 13可知, 升温速率越高, 晶体内部应力越大, 塑性滑移系启动越早, 且累积剪应变越大。

4 结论

考虑非线性弹性、晶体塑性和温度相关的固固相变的晶体塑性本构模型, 研究了烤燃下的HMX相变及各向异性相应。

(1) HMX相变是可逆的, 但是路径不同。因为当晶体重新回到β相时, 会有残余应变和残余应力。此外, β→δ伴随着明显的体积膨胀, 会产生裂纹等损伤, 而当晶体回到β相时, 这些损伤不会消失。同时, 当高温持续时间较长时, 相变更完全, 温度均匀化有利于相变的发生。

(2) 当HMX晶体的温度高于临界温度时, HMX开始β→δ相变, 相变会导致体积膨胀, 并且相变后δ相的由于热膨胀会产生更大的热应变, 从而产生晶体应变增大。相变后晶体内应力也会增加, 当晶体内分解剪应力足够大时, 晶体的滑移系启动。升温速率越高, 相变越快, 同时晶体内应力和应变更大。相变后晶格膨胀会导致晶体更易生成裂纹, 因此升温速率越高, HMX晶体越容易相变, 并且HMX晶体越容易点火。