在连续介质力学理论中，物质占据空间几何区域所构成的图形，称为物质构型。在变形或刚体运动过程中，物质构型将发生改变。为了对晶体材料的变形过程进行描述，需要建立晶体运动学框架。初始时刻未经变形的晶体材料占据的空间几何图形称为初始构型。采用拉格朗日坐标(X)来描述初始构型中任意物质点的位置。晶体材料在变形过程中的瞬时空间几何图形称为现时构型。对于初始构型中位于坐标X的物质点，在现时构型中采用欧拉坐标(x)描述其瞬时空间位置。显然，物质点在初始和现时构型之间具有x = x(X, t)的位置对应关系。描述材料变形程度的基本物理量通常选择变形梯度F = ∂x/∂X，其物理含义是现时构型相对于初始构型的梯度张量。物质点的运动速度通常定义在拉格朗日坐标系中v(X, t) = ∂x(X, t)/∂t，也可采用欧拉坐标表示为v(x, t) = v[X(x, t), t]，通常定义速度梯度L = ∂v/∂x描述物质点速度在现时构型中的梯度。

经典的运动学构型如图1所示。晶体材料由初始构型经变形梯度(F)变换至现时构型，如图1中的黑色箭头所示。Hill等^[21]、Asaro等^[22]建立中间构型Ⅰ，将弹性和塑性变形进行乘和分解

式中：F^e、F^p分别为弹性和塑性变形梯度。材料由初始构型经塑性变形(F^p)变换至中间构型Ⅰ，再经弹性变形(F^e)变换至现时构型，如图1中的红色箭头所示。弹性变形包括刚体旋转和晶格畸变，塑性变形可由位错滑移、相变、孪生等引起。将现时构型中速度梯度张量进行弹、塑性分解

式中：$ {\dot{{{F}}}}^{\rm{e}}$、${\dot{{{F}}}}^{\rm{p}}$分别为弹性和塑性变形梯度的变化率，L^e、L^p分别为现时构型中弹性和塑性速度梯度，L₀^p为中间构型Ⅰ中由塑性滑移引起的速度梯度。对于弹性变形部分，若采用对数应变$ {{{\varepsilon }}}^{\rm{e}}={\rm{l}}{\rm{n}}{{U}}^{\rm{e}}$作为应变度量，一般还需进一步建立中间构型Ⅱ，采用极分解将弹性变形分解为拉伸和旋转部分

式中：R^e为旋转张量，U^e为右伸长张量。材料由中间构型Ⅰ，经弹性拉伸(U^e)变换至中间构型Ⅱ，再经晶格旋转(R^e)变换至现时构型，如图1中的蓝色箭头所示。

材料在动态冲击加载下往往存在显著的温升现象，而温升引起的热膨胀会导致材料内部应力、应变的变化，因此如何考虑热膨胀效应的影响是动态晶体塑性模型的一个重要问题。

在晶体塑性模型框架中，通常有两种考虑热膨胀效应的方法。第1种方法是在弹塑性两构型的基础上额外引入热构型，将热变形与弹性变形解耦。第2种方法不额外引入热构型，而在现有框架下基于热力学自由能进行推导，在求解超弹性本构过程中考虑热膨胀效应(详见2.2.2节)。本节主要介绍第1种引入热构型的方法。Clayton^[23]、Vogler等^[24]、Hansen等^[25]、De等^[26]以及Shahba等^[27]在构型分解时加入了热构型F^θ，如图2所示，将变形梯度分解为塑性变形、热变形以及弹性变形3部分

式中：F^e、F^θ、F^p分别表示弹性变形梯度、热变形梯度和塑性变形梯度。

相应地，现时构型中的速度梯度可分解为

式中：L^θ、${{{L}}}_{0}^{\theta } $分别为现时构型和热构型中的热变形速度梯度，${ {{{L}}}_{0}^{\theta }={\dot{{{F}}}}^{\theta }({{{F}}}^\theta)^{ -1},\;{{{L}}}_{0}^{\rm{p}}={\dot{{{F}}}}^{\rm{p}}({{{F}}}^{\rm{p}})^{-1}}$。${{{L}}}_{0}^{\theta } $是温度变化率($ \dot{T}$)的函数

式中：α为热膨胀系数张量。

通常在晶体塑性模型中有超弹性和次弹性两种描述材料弹性行为的本构关系。超弹性本构是指材料存在一个与应变相关的能量函数，其应力与应变全量之间的关系由能量函数决定，具有做功独立于加载路径的特点。次弹性本构描述应力率与应变率之间的关系，由于大变形下弹性模量是随应变全量变化的，其实际应用比较复杂，通常只用于小弹性变形情况^[28](这种情况下弹性模量可视为常数)。考虑到实验中通常采用全量形式的高压固体状态方程描述材料的非线性容变关系，因此研究材料在冲击载荷下的大弹性变形行为时，往往采用超弹性本构来描述。

材料在冲击载荷下的超弹性本构关系主要包含基于热力学自由能的超弹性应力表达式以及考虑热弹性耦合效应的温升表达式两方面。其中，特别需要考虑热膨胀对应力的影响、弹性模量随压强和温度的变化，以及非线性弹性容变关系(固体状态方程)等。

超弹性本构通常根据热力学第一、第二定律，建立热力学势函数所对应的耗散不等式，得到相应的热力学应力表达式。Becker^[29]、Pi等^[30]、De等^[26]基于内能构建热力学势函数并建立相应的超弹性本构框架，而更多的研究学者，如Clayton^[23]、Vogler等^[24]、Luscher等^[31-32]、Santos等^[33]、Addessio等^[14]则基于亥姆霍兹自由能构建相应的热力学势函数。分析表明，基于内能或亥姆霍兹自由能构建热力学势函数并推导所得的热力学应力表达式是等价的。本节将以形式更简洁的亥姆霍兹自由能为例，详细介绍超弹性本构框架。

根据热力学第二定律，对于任何孤立系统(包括系统A与环境Sur)，整体的熵增恒大于或等于零

式中：无限大环境的熵增与可逆过程的热温熵一致，即${\text{Δ}} {S_{{\rm{Sur}}}} = - Q/T$。式(7)对于系统A的任意微小体积单元及其外部环境同样适用，此时系统微小体积单元的热变化Q由系统内热源生成的热与从系统外流入的热量两部分构成，熵增定律最终可转化为率形式的Clausius-Duhem耗散不等式

式中：$ {\tilde \rho }$为参考构型中的材料密度；$ \dot{s}$为系统局部的比熵增率；R表示单位质量的热生成率，产热来自于内部的热源；H表示由热传导引起的通过单位面积上的热流矢量；D_int为材料的内禀耗散项，表示材料内部与热源、热传导无关的熵增项。

将能量守恒表达式$ {\tilde \rho }\dot{U}={{\sigma}}:{\dot{{{\varepsilon }}}}+{\tilde \rho }R-{\rm{div}}{{H}}$代入式(8)，同时根据亥姆霍兹自由能与内能的关系式ψ = U – sT，耗散不等式最终可用亥姆霍兹自由能表示为

单位体积亥姆霍兹自由能(比亥姆霍兹自由能)通常可表示为弹性变形ε^e、温度T以及内状态变量q_n的函数$ \psi \left({{{\varepsilon }}}^{\rm{e}},T,{q}_{n}\right)$，其率形式可表示为

最终将比亥姆霍兹自由能变化率$ \dot{\psi }$的表达式代入耗散不等式(式(9))可得

式中：$ {\tilde \rho }{\displaystyle\sum \limits_{n=0}^{M}\dfrac{\partial \psi }{\partial {q}_{n}}{\dot{q}}_{n}}$表示材料中由塑性变形引起的与位错等缺陷相关的微观畸变能的变化率^[26]，在实际计算过程中通常取畸变能变化率为塑性功率的某一比例。如该比例取10%^{[27, 34]}时，塑性功的90%经过耗散转换成热能贡献于温升，其余10%则转换为畸变能存储在位错等微观结构中。

显然，上述耗散不等式对于任意的$ {\dot{{{\varepsilon }}}}^{\rm{e}}$与$ \dot{T}$均成立，因此可得到热力学应力与自由能的关系以及内禀耗散项的最终表达式

材料在冲击加载下的高压变形往往表现出一定的黏性，因此在热力学应力之外还需额外考虑材料黏性所产生的变形阻力。黏性反映了物体在体积发生剧烈变化时所受到的阻力与体积变化速度之间的关系。人工黏性最初是由von Neumann与Richtmyer^[35]为了解决有限差分法在模拟冲击波传播过程中的振荡问题所提出的，后来被逐渐引入冲击波的流体方程组中，变成了一个具有物理意义的附加项，即黏性耗散项^[31]，与压强的量纲相同

式中：J = V/V₀表示体积变形，$ \dot{J}$表示体积变形率，c₀是与温度、压强相关的系数，c₁、c₂为常数。

在考虑黏性耗散压的情况下，总应力由准保守的热力学应力与非保守的耗散应力两部分组成

热弹性耦合指应力、弹性应变与温度的相互影响，包含变形引起的温升以及温度对变形及应力的影响。

温升引起的热效应是冲击加载过程中不可忽视的重要因素。在弹性本构中，温升引起的热膨胀会导致应力、应变的变化；在塑性本构中，温度也是影响位错滑移、相变等塑性变形行为的重要因素。

由经典热力学关系可知，状态方程中的独立变量仅有两个。温度T可表示为比熵与弹性应变的状态函数$ T\left(s,{{{\varepsilon }}}^{\rm{e}}\right)$，通过泰勒展开可得

其率形式为

式中：$ {\left(\dfrac{\partial T}{\partial s}\right)}_{V}=\dfrac{T}{{c}_{V}}$，$ {\left(\dfrac{\partial T}{\partial {{{\varepsilon }}}^{\rm{e}}}\right)}_{s}=\dfrac{1}{{\tilde \rho }}{\left(\dfrac{\partial \sigma }{\partial s}\right)}_{V}=-T{{\varGamma}}$为Maxwell关系式，$ {{\varGamma}}=\dfrac{1}{{\tilde \rho }}{\left(\dfrac{\partial p}{\partial U}\right)}_{V}=\dfrac{\alpha {K}_{T}}{{c}_{V}{\tilde \rho }}$为Grüneisen系数，U为比内能。准静态加载下弹性变形较小，往往忽略弹性变形引起的等熵温升，而只考虑温升与熵增之间的热容关系。但是材料在冲击载荷下会发生较大的弹性变形，弹性变形引起的等熵温升不再可以忽略。

对于历时极短的冲击加载过程，由于来不及进行热交换，可将变形材料视为绝热体系，这时内禀耗散项$ {D}_{\rm{int}}={\tilde \rho }T\dot{s}-\left({\tilde \rho }R-{\rm{div}}{{H}}\right)$在绝热且无热源的情况下退化为$ {\tilde \rho }T\dot{s}$项。根据Clausius-Duhem耗散不等式(式(13))，在考虑黏性耗散作用的情况下，熵增项$ {\tilde \rho }T\dot{s}$可表示为

将熵增表达式(式(18))代入式(17)可得温度的变化率，其中温度的变化包括塑性变形、微观畸变、黏性耗散和弹性变形4部分的贡献

热膨胀是热弹性耦合的另一方面，表现为温度变化引起了应力与弹性应变的改变，即产生了额外的热应变与热应力，这里的弹性应变为弹塑性构型中的表观弹性应变，包括温升引起的热膨胀变形和应力引起的纯弹性变形两部分。在晶体塑性模型框架中，通常有两种等价的考虑温升对变形影响的方法。

第1种方法是在弹塑性两构型分解中引入热膨胀构型(详见1.2节)。在此三构型框架下，求得的弹性应变$ {\bar{{\varepsilon }} }^{\rm{e}}$为表观弹性应变与热应变之差$ \left({{{\varepsilon }}}^{\rm{e}}-{{\alpha}}{\text{Δ}}T\right)$，代表了真实的不包含热膨胀的产生应力的应变部分。第2种方法则仍然采用传统的弹塑性两构型分解方法，不额外引入热膨胀构型。由于热膨胀效应可体现在比熵与应变的偏导关系中，它对应力的影响可在求解应力表达式$ {{\sigma}}={\left(\dfrac{\partial {\tilde \rho }\psi }{\partial {{{\varepsilon }}}^{\rm{e}}}\right)}_{T}$的过程中推出。

本节以形式更为简洁的第2种方法为例，考虑热膨胀效应推导热力学应力$ {{\sigma}}={\left(\dfrac{\partial {\tilde \rho }\psi }{\partial {{{\varepsilon }}}^{\rm{e}}}\right)}_{T}$的表达式^[31]。系统的比亥姆霍兹自由能$ {\tilde \rho }\psi \left({{{\varepsilon }}}^{\rm{e}},T,{q}_{n}\right)$可表示为弹性应变、温度以及内变量q_n的函数

基于上述亥姆霍兹自由能，热力学应力可根据式(12)进行求解。需要注意的是，热力学应力在求解过程中需考虑弹性模量随压强的变化以及比熵与应变之间的耦合关系(即热膨胀效应)

式中：热弹性耦合体现为Maxwell方程中熵增与弹性应变的关系式，即$ {\tilde \rho }{\left(\dfrac{\partial s}{\partial {{{\varepsilon }}}^{\rm{e}}}\right)}_{T}=-{\left(\dfrac{\partial {{\sigma}}}{\partial T}\right)}_{V}=-\dfrac{{c}_{V}}{T}{\left(\dfrac{\partial T}{\partial {{{\varepsilon }}}^{\rm{e}}}\right)}_{s}=$${c}_{V}{{\varGamma}}={{C}}:{{\alpha}} $(考虑各向异性的Grüneisen张量)，最终热力学应力表达式为

由式(22)可知，表达式中的$ {{{\varepsilon }}}^{\rm{e}}-{{\alpha}}{\text{Δ}}T$为真实的弹性应变，是表观弹性应变与热膨胀应变之差。这一结果与引入热构型方法得到的真实弹性应变一致。需要指出的是，式(22)中第2项中的ε^e是表观弹性应变，而引入热构型的方法中表观弹性应变ε^e被真实弹性应变$ {{{\varepsilon }}}^{\rm{e}}-{{\alpha}}{\text{Δ}}T$替代，两者的差别属于高阶小量，可忽略不计。

材料的本构关系通常可分为描述体积变化的球量部分(容变关系)以及描述形状变化的偏量部分(畸变关系)。由于材料在冲击载荷下的体积变化很大，不再满足线弹性关系，因此通常采用高压状态方程描述材料的非线性弹性容变关系(压强与比容的关系)。本节将介绍高压冲击下常用的状态方程以及将状态方程引入超弹性本构的方法。

所谓状态方程，即物质系统从原有的热力学平衡状态转变为另一热力学平衡状态过程中各个状态参数之间的关系。在固体材料的高压、高温问题的研究中，常采用压力-温度-比容状态方程作为基本关系式

式中：p为压力，V为比容，T为温度。基于该方程，通过热力学关系推导，可得其他所有状态参量。

在固体状态方程研究中，一般采用经典的三项式物态方程形式。将固体材料中的能量和压力分解为与温度无关的“冷”贡献部分(即冷能和冷压)，以及依赖于温度的“热”贡献部分(即热能和热压)，其中“热”贡献包括晶格热振动贡献和自由电子热振动贡献。因此，能量和压力可表达为

式中：E_C为冷能，E_T为热能，E_TN为晶格热贡献能量，E_TE为电子热贡献能量；p_C为冷压，p_T为热压，p_TN为晶格热贡献压力，p_TE为电子热贡献压力。一般情况下，电子热贡献压力项可忽略，仅在极高压阶段(100 TPa以上)才需考虑自由电子热振动的影响。下面将具体介绍常用的冷贡献状态方程，如Murnaghan状态方程、Vinet状态方程，以及晶格热贡献状态方程，如Grüneisen状态方程、基于Debye模型的状态方程。

（1）Murnaghan状态方程

Murnaghan状态方程是描述冷贡献的一个重要方程。当加载过程中冲击压力较低时，体积模量与压力之间的依赖关系可近似地采用一阶泰勒展开式描述

式中：K(p)为依赖于压力的瞬态体积模量，K₀、K′分别为零压体积模量和体积模量关于压力的导数。通过对式(25)中第2个等式进行积分，可得到Murnaghan状态方程

式中：p_Mur为Murnaghan状态方程描述的冷贡献压力，V₀、V分别为零压比容和变形中的瞬态比容。根据弹性体积应变和比容之间的关系($ V/{V}_{0}=1+{\varepsilon }_{v}^{\rm{e}}$)以及指数函数幂级数展开式的一阶近似，可得到Murnaghan状态方程的等价形式

式中：体积应变$ {\varepsilon }_{v}^{\rm{e}}={{{\varepsilon }}}^{\rm{e}}:{{I}}$，ε^e为弹性应变张量，I为二阶单位张量。

对于通常的固体材料，Murnaghan状态方程适用的压力范围可达到太帕量级。对于更高压力范围的动态研究，需要采用高阶泰勒展开式代替一阶近似来描述体积模量与压力的关系，如三阶Birch-Murnaghan状态方程

能够对极高压力下的压力-比容关系进行更准确的数值描述。

（2）Vinet状态方程

材料的冷压除了可通过对体积模量的微分形式积分得到之外，也可通过冷能即晶格结合能对比容的偏导关系获得。例如Vinet状态方程^[36]采用亥姆霍兹自由能来描述固体材料的晶格结合能即冷能

式中：K₀和K′分别表示零压体积模量以及瞬时体积模量对压强的导数，V₀、V分别为零压比容和变形中的瞬态比容。

根据系统冷压与冷能之间的关系$ {p}_{\rm{C}}=-{\left(\dfrac{\partial {\psi }_{\rm{C}}}{\partial V}\right)}_{T}$，可得到Vinet形式的冷压表达式

（1）Grüneisen状态方程

Grüneisen状态方程是描述晶格热贡献的一个重要方程。引入热力学参量Γ用于描述等体积条件下压力关于内能的变化率，称为Grüneisen系数。考虑到等体积条件下冷压p_C(V)为常数，可将Γ表示为

假设Γ对温度不敏感，仅为比容V的函数，在比容V恒定的情况下对式(32)中定义的热压和热能关系式进行积分，即可得到Grüneisen状态方程

式中：p_G为Grüneisen状态方程描述的晶格热振动贡献压力，晶格热能$ {E}_{T}={\int }_{{T}_{0}}^{T}{c}_{V}{\rm{d}}T$，T₀为参考温度，c_V为定容比热容。

关于Grüneisen系数的计算，大致有3种方法。传统的方法是基于热力学基本公式，推导Grüneisen系数与热力学参数之间的关系

式中：α为热膨胀系数，K_T、K_S分别为等温和等熵体积模量，c_V、c_p分别为定容比热容和定压比热容。事实上Grüneisen系数并非常数，会随着材料体积的变化而变化。第2种方法基于Debye模型描述的微观热振动方法(Migault公式)，将Grüneisen系数对比容的依赖表示为

式中：a为常数，且有3种理论公式可选用。当a = 0时，式(35)为Slater公式；当a = 2/3时，为Dugdale-MacDonald（D-M）公式；当a = 4/3时，为自由体积公式。然而，特定变形条件下对a的选择令学者们感到困惑。李欣竹^[37]基于冲击实验结果建立了第3种方法：首先将实验得到的比容、冲击波速和物质点运动速度分别代入冲击波质量守恒方程和能量守恒方程，得到冲击压力和冲击能量，再由Born-Meyer势得到冷压和冷能，代入式(33)即可得到Grüneisen系数，即Γ = V(p_H–p_C)/(E_H –E_C)。

（2）基于Debye模型的状态方程

材料的晶格热压除了可直接借助压力与热能之间的关系即Grüneisen系数来获得，也可用热能与比容之间的偏导关系推导得到。Debye模型是常用的描述晶格振动与热能关系的模型，该模型把晶格振动的频率看作连续分布的函数，晶格振动能可表示为

式中：R、M_mol分别为摩尔气体常数和摩尔质量，Debye函数$ D\left(x\right) = {\left(\dfrac{1}{x}\right)}^{3}{\int }_{0}^{x}\dfrac{{z}^{3}}{{\rm{e}}^{z}-1}{\rm{d}}z$，T_D为Debye温度。

根据系统热压与热能之间的关系$ {p}_{T}=-{\left(\dfrac{\partial {\psi }_{\rm{TN}}}{\partial V}\right)}_{T}$，可得到基于Debye模型的热压表达式

式中：$ \varGamma \left(V\right) = -\dfrac{\partial {\rm{l}}{\rm{n}}{T}_{\rm{D}}}{\partial {\rm{l}}{\rm{n}}V}$表示与体积相关的Grüneisen系数。

此外，Cawkwell等^[38]认为，在分子晶体中，除了晶格声子振动外，还存在分子内振动，且两种振动模式对温度和体积的依赖性有所区别。因此对于分子晶体，晶格振动能可由两部分Debye模型组成，即Ψ_TN = Ψ_D,1 + Ψ_D,2

式中：N = R/k_BM_mol表示单位质量的原子数量，k_B为玻尔兹曼常数，Q表示大多数的低频声子数，3N – Q则为分子内振动的声子数。

为了将状态方程引入应力表达式中，需要考虑应变球量和偏量的解耦。可直接解耦的应变张量应满足在等体积变形时对角元素之和(球量)为零，这样可用小变形下相同的方法将球偏量解耦。定义在中间构型Ⅱ(见1.1节)上的对数应变ε^e = ln U^e满足此要求，可直接做球偏量分解，并引入状态方程。但是定义在中间构型Ⅰ(见1.1节)上的格林-拉格朗日应变E^e = (F^eTF^e–I)/2不满足此要求，无法直接做球偏量分解。Luscher等^[31]指出实际上并不存在与对数应变ε^e完全功共轭的应力张量，因此本节以严格功共轭的格林-拉格朗日应变E^e和PKⅡ应力为例，详细介绍对E^e球偏量分解的方法以及将状态方程引入应力表达式的过程。

Luscher等^{[31-32, 39]}提出了将与纯体积变形相关的应变部分从格林-拉格朗日应变E^e中解耦的方法。由于E^e无法直接分解，先将变形梯度F^e分解为体积变形部分和等体积变形部分

式中：$ \overline {{{{F}}^{\rm{e}}}}={J}^{1/3}{{I}}$为体积变形梯度，$ \widehat {{{{F}}^{\rm{e}}}}={J}^{-1/3}{{{F}}}^{\rm{e}}$为等体积变形梯度。

根据格林-拉格朗日应变的一般表达式$ {{E}}^{\rm{e}}=\left({{{F}}}^{\rm{eT}}{{{F}}}^{\rm{e}}-{{I}}\right)/2$，相应的等体积应变与体积应变分别表示为

显然$ {{E}}^{\rm{e}}\ne \widehat{{{E}}^{\rm{e}}}+\overline {{{{E}}^{\rm{e}}}}$，这是由于等体积应变$ \widehat{{{E}}^{\rm{e}}}$定义在体积变形构型中，而格林-拉格朗日应变E^e和体积应变$ \overline {{{{E}}^{\rm{e}}}}$均定义在中间构型Ⅰ中。将$ \widehat{{{E}}^{\rm{e}}}$后拉至中间构型Ⅰ得到

最终定义在相同构型(即中间构型Ⅰ)中的等体积应变与体积应变满足

式中：$ \delta =\dfrac{3}{2}\left({J}^{2/3}-1\right)$表示总体积应变。

为了将状态方程引入应力表达式，需将式(44)中仅描述材料体积变化的自由能和应力部分与描述形变的部分解耦出来。式(44)可由亥姆霍兹自由能式(20)，将比熵s泰勒展开得到

式中：与弹性变形无关的自由能均归纳在$ {\tilde \rho }{\psi }_{{\rm{r}}{\rm{ef}}}\left(T,{q}_{n}\right)$项中。将式(44)中的应变E^e分解为等体积应变$ \widehat{\widehat{{{E}}^{\rm{e}}}}$与体积应变$ \overline {{{{E}}^{\rm{e}}}}$，将热膨胀系数α分解为球量与偏量部分$ {{\alpha}}=\dfrac{1}{3}{\alpha }_{v}{{I}}+{{{\alpha}}'}$。显然，仅仅由体积膨胀引起的自由能变化描述的是与体积变化相关的自由能部分

式中：K为体积模量，$ K=\dfrac{1}{9}{{I}}:{{C}}:{{I}}$。需要指出的是，亥姆霍兹自由能关系式(式(44))中体积应变与等体积应变的耦合项$ \overline {{{{E}}^{\rm{e}}}}:{{C}}:\widehat{\widehat{{{E}}^{\rm{e}}}}$中虽然包含了体积应变$ \widehat{\widehat{{{E}}^{\rm{e}}}}$，但这部分自由能与描述材料仅在球量变形状态下的压强(状态方程)无关。

从这部分自由能可得到与体积变化相关的应力为

式中：右柯西-格林变形张量$ {{{C}}}^{\rm{e}}={{F}}^{\rm{eT}}{{F}}^{\rm{e}}$，体积变化J对格林应变的偏导数∂J/∂E^e = –J(C^e)⁻¹，压强对格林应变的偏导数∂p/∂E^e = –JK(C^e)^-1。

显然，这部分描述体积变化的应力与描述压力p变化规律的状态方程相关，可转化为用状态方程p_EOS表示的更为精确的表达式S_EOS。而其他描述形状变化的应力则可用应力(式(22))减去式(46)得到。最终PKⅡ应力S = S_EOS + S_dev可表示为

式中：状态方程$ {p}_{\rm{EOS}}=-\dfrac{\partial {\tilde \rho }{\psi }_{v}}{\partial V}$包含了冷压与热压两部分，$ {p}_{\rm{EOS}}={p}_{\rm{C}}+{p}_{T}$，$ \dfrac{\partial V}{\partial {{E}}^{\rm{e}}}=-J\left({{{{C}}}^{\rm{e}}}\right)^{-1}$。

晶体塑性本构框架一般包括流动准则和硬化方程两部分。准静态晶体塑性框架中，流动准则主要有率相关、率无关的流动准则两类；硬化方程通常考虑材料的应变硬化、应变率硬化和热软化效应。动态晶体塑性本构框架中，考虑到高温、高压、高应变速率的动态冲击过程中，固体材料的黏塑性即应变率效应变得尤为显著^[40]，一般需要建立率相关的流动准则，并且高速变形过程中的塑性滑移阻力与声子黏性拖曳密切相关，因此动态硬化方程中还应当考虑声子拖曳效应。

动态晶体塑性本构框架可分为唯象的动态晶体塑性框架和基于位错行为的动态晶体塑性框架两类，下面分别对这两类晶体塑性框架进行详细介绍。

唯象的动态晶体塑性框架包括流动准则和硬化本构方程两部分，它们描述了冲击变形过程中滑移速率、流动应力等材料状态变量的演变规律，其唯象描述与具体物理机制越接近，对材料动态冲击变形的模拟越准确，相应的模型预测能力越强^[27]。本节主要介绍高速冲击加载条件下，唯象的动态晶体塑性有限元研究中建立的流动准则和硬化本构方程。

率相关的流动准则描述了塑性应变速率与滑移系分切应力、塑性滑移阻力以及微观结构相关的其他内变量之间的关系。现有的率相关流动准则主要包括幂形式黏塑性模型和热激活模型等。

幂形式的黏塑性流动准则由Hutchinson^[41]在研究准静态变形问题时提出，将滑移系α上的滑移速率$ {\dot{\gamma }}^{\alpha }$表示为滑移系上分切应力和滑移阻力比值的幂函数

式中：$ {\dot{\gamma }}_{0}$为参考滑移速率，τ^α为滑移系α上的分切应力，g^α为滑移系α上的临界分切应力即滑移阻力，m为应变率敏感性因子。当应变率敏感因子m→0时，即为率无关流动准则。Li等^[42]、Khan等^[43]将幂形式黏塑性模型应用于动态单轴压缩变形模拟，Zhang等^[44]、Lim等^[45]应用于Taylor冲击变形的晶体塑性有限元模拟，Bobbili等^[46]应用于分离式霍普金森压杆(SHPB)试验的动态晶体塑性模拟研究中。

Becker^[29]、De等^[26]将幂形式黏塑性流动方程(式(49))进行转化，成为形式上的率无关流动准则。但是为了考虑材料的应变率效应，他们对滑移阻力进行了率相关的修正，从而建立形式上率无关但本质上率相关的流动方程

式中：$ {\widehat{\tau }}^{\alpha}$为最初率无关形式的流动准则中定义的滑移阻力。率相关形式的流动准则中定义的滑移阻力g^α与$ {\widehat{\tau }^{\alpha }}$之间的关系为$\widehat{\tau}^\alpha = {g^\alpha }{\left( {{{\dot \gamma }^\alpha }/{{\dot \gamma }_0}} \right)^m}{\rm{sgn}}\left( {{{\dot \gamma }^\alpha }} \right)$，这样式(50)与式(49)实际上是等价的。当滑移系上的分切应力小于临界分切应力或滑移阻力时，该滑移系不发生塑性滑移；一旦滑移系上的分切应力达到临界分切应力，则该滑移系启动，随之开始发生塑性滑移。动态冲击变形过程与准静态变形过程相比，塑性滑移速率$ {\dot{\gamma }}^{\alpha }$显著增大，参考滑移速率$ {\dot{\gamma }}_{0}$相应地增大。

Kocks等^[47]基于Arrhenius方程提出热激活流动准则，能够显式地反映塑性滑移热激活过程对温度的依赖。滑移系α上的热激活流动准则可以表示为

式中：Q^α为滑移系α上发生塑性滑移的热激活能，T为瞬态材料温度。参考塑性滑移速率为$ {\dot{\gamma }}_{0}^{\alpha }$，成功跨越势垒Q^α使得塑性滑移发生的概率为exp(–Q^α/k_BT)，两者的乘积即为滑移系α上真实的塑性滑移速率$ {\dot{\gamma }}^{\alpha }$。显然，环境温度或变形引起的温升越大，热起伏提供的能量越大，跨越势垒所需的力越小，如图3(a)所示。$ {\dot{\gamma }}_{0}^{\alpha }$通常是与加载压力和应变速率相关的常数，因此跨越势垒的概率越大，相应的塑性滑移速率越大。

塑性滑移过程的热激活能依赖于局部流动应力和滑移阻力

式中：${Q}_{0}^{\alpha} $为参考激活能；m₁、m₂是与滑移阻力相关的材料参数；g^α为滑移系α上克服短程势垒的滑移阻力，对温度和应变率敏感；${g}_{\rm{ath}}^{\alpha} $为滑移系α上克服与缺陷相关长程势垒的滑移阻力，对温度和应变率不敏感，如图3(b)所示。随动态变形过程中材料局部流动应力的增大，热激活能减小，跨越势垒的概率增大，相应的塑性滑移速率增大。Khan等^[43]、Li等^[48]分别将热激活流动准则应用于模拟动态单轴压缩和单轴拉伸变形中，Bronkhorst等^[49]、Chandra等^[50]应用于分离式霍普金森压杆(SHPB)试验的动态晶体塑性模拟研究中，Luscher等^[31]、Addessio等^[14]应用于模拟平板撞击试验中的冲击变形行为。

建立适用于高温、高压、高应变速率变形条件下材料的硬化模型，需要考虑的机制包括：应变硬化效应、应变率硬化效应、声子拖曳效应、热软化效应等。动态晶体塑性框架中，在每个晶体滑移系上分别建立热黏塑性硬化方程，其一般形式为$ {g}^{\alpha }=f\left({\gamma }^{\alpha },{\dot{\gamma }}^{\alpha },T\right)$。唯象的动态晶体塑性框架中建立的硬化方程大致可以分为两类：应变硬化方程和位错密度硬化方程。

经典的唯象硬化方程考虑材料的应变硬化效应。晶体材料在外加载荷作用下，当滑移系上的分切应力达到临界分切应力即滑移阻力时，该滑移系开始发生塑性滑移。随后的变形过程中滑移阻力不断增大，维持该滑移系继续塑性滑移所需的分切应力相应增大，称为滑移系的自硬化现象；由于滑移系之间的交互作用，该滑移系上发生塑性滑移能够影响其他滑移系硬化的能力，称为滑移系的潜硬化现象^[53]。基于此，引入滑移系硬化模量(h_αβ)来描述滑移阻力随时间的演化

进一步引入参数r和h_β建立瞬态硬化模量的表达式^[41]

可见，当β = α时，h_ββ = h_β(不求和)为滑移系自硬化模量；当β ≠ α时，h_αβ = rh_β为滑移系潜硬化模量。

自硬化模量(h_β)的形式有很多。Bittencourt^[54]、Bobbili等^[46]将滑移系自硬化模量简化为一个常数，并假设晶体中不同滑移系的硬化能力相同

式中：h₀为常数。然而，考虑真实的材料硬化特征，学者们更多地采用变化的自硬化模量。Peirce等^[55]、Li等^[42]、Khan等^[43]引入名义饱和滑移阻力(g_sn)，采用双曲函数型自硬化模量

式中：g₀为初始滑移阻力，$ \gamma ={\int }_{0}^{t}\sum\limits_{\beta }\left|{\dot{\gamma }}^{\beta }\right|{\rm{d}}t$为所有滑移系上的累积塑性滑移量。随着塑性变形的增大，一方面双曲函数型自硬化模量逐渐减小；另一方面滑移阻力逐渐增大并且趋于饱和值，甚至可以略微超过名义饱和滑移阻力g_sn。Zhang等^[44]、Chandra等^[50]、Li等^[48]引入实际饱和滑移阻力(g_s)，采用幂形式的自硬化模量

式中：幂指数a < 1。随着塑性变形的增大，滑移阻力逐渐增大，相应的自硬化模量逐渐减小；当滑移阻力达到饱和值之后不再改变，此后的自硬化模量始终为零。考虑到不同滑移系上的饱和滑移阻力未必相同，Acharya等^[56]、Bronkhorst等^[49]、Luscher等^[31]、Addessio等^[14]引入滑移系上的饱和滑移阻力($g_{\rm{s}}^{\beta} $)。滑移系β上的自硬化模量h_β = f(g^β)在g^β∈[g₀, $g_{\rm{s}}^{\beta} $]区间内有定义，并且区间端点值为${{h}}_{\beta}|_{{g}^{\beta}=g_0}={{h_0}} $，${{h}}_{\beta}|_{{g}^{\beta}=g_{\rm{s}}^{\beta}}=0$，他们采用线性插值函数

式中：g_s0为滑移系β上的参考滑移阻力，k_B为玻尔兹曼常数，Q为参考激活能。随着塑性滑移速率的增大，滑移系上的饱和滑移阻力相应地增大，这也体现了应变率硬化效应。此外，Khan等^[43]考虑晶粒尺寸效应，对经典的唯象硬化方程进行改进，将滑移系自硬化模量表示为晶粒尺寸、滑移量和滑移速率的函数

式中：d为平均晶粒尺寸，A、k、n₀为Hall-Petch系数，B、n₁、n₂为应变率硬化和应变硬化相关的材料参数。

动态冲击变形过程中有显著的温升，这是不同于准静态变形的重要特点之一，所以动态晶体塑性硬化模型中需要考虑温升引起的热软化。将热软化机制引入唯象的动态晶体塑性框架的思路是，在温度无关的硬化方程的基础上，对变形过程中的滑移阻力进行温度相关的修正。基于式(53)定义的塑性滑移阻力g^β，Clayton^{[23, 57-58]}、Vogler等^[24]建立了温度相关的塑性滑移阻力g^β(T)与变形温度之间的幂形式关系

式中：T、T_r分别为瞬时温度和参考温度，C为热软化指数。由于动态冲击过程中的剧烈温升可能使变形温度接近材料熔点，Tajalli等^[59]、Khan等^[43]、Li等^[60-61]引入材料熔点参量，对热软化机制起作用的温度范围规定了上限，将滑移系自硬化模量表示为

式中：T_m为材料熔点。Becker^[29]基于Steinberg-Guinan硬化模型和Arrhenius形式的热软化模型，将滑移系α上的滑移阻力表示为

式中：μ₀、μ分别为零压和瞬态的材料剪切模量，B、n为应变硬化参数，γ为所有滑移系上的累积滑移量，γ₀为初始塑性滑移量(通常为零)，λ为热软化修正系数，${\text{Δ}} {E_T}$为比热能变化，τ₀为参考屈服应力。他们评估了热软化效应对多晶Ta斜波加载冲击变形的影响，如图4^[29]所示。结果表明，基于依赖于温度的材料本构计算得到的塑性应变有所提高，并且该现象在晶界处更为显著。

在极高变形速率下(大于10⁴ s^–1)，塑性滑移阻力与声子黏性拖曳密切相关，针对动态冲击问题进行模拟时要求在唯象的硬化模型中引入声子拖曳机制。在Becker发展的动态晶体塑性硬化模型的基础上，De等^[26]考虑声子拖曳效应，建立滑移系α上的滑移阻力方程

式中：$ {\dot{\gamma }}_{0}$为参考滑移速率，$ {\dot{\gamma }}_{\rm{off}}$为滑移系率硬化修正参数，τ₀、τ_T分别为常温和瞬态参考屈服应力，θ为线性声子拖曳系数。他们在Becker模型的基础上，一方面对滑移阻力进行了率相关的修正，同时考虑了应变硬化效应、应变率硬化效应和热软化效应对滑移阻力的贡献，即式(63)中的第1项；另一方面引入了声子拖曳效应对滑移阻力的贡献，即式(63)中的第2项。需要指出的是，他们建立的流动方程和晶体塑性本构均难以简单地写出滑移系上的塑性应变率($ {\dot{\gamma }}^{\alpha }$)和滑移阻力(g^α)的显式方程。他们通过对α-RDX单晶的平板撞击过程进行模拟，发现未考虑声子拖曳效应的晶体塑性模型对极高应变速率下的滑移临界分切应力预测过低，如图5^[26]所示。

基于位错密度建立硬化模型，相比于传统的唯象应变硬化模型更准确，并且物理意义更明确。Khan等^[43]、Lim等^[45]将滑移阻力表达为位错密度的函数

式中：τ₀为参考屈服应力，μ为剪切模量，b为柏氏矢量，H_αβ为硬化矩阵，ρ^β为滑移系β上的位错密度。考虑位错密度的增殖和湮灭机制，将位错密度变化率表示为瞬态位错密度和塑性滑移速率的函数，建立唯象的位错密度演化方程

式中：k₁、k₂分别为位错增殖和湮灭项系数，A_αβ为位错增殖相互作用的张量分量，${{{\dot \gamma }^\alpha }}$为滑移系α上的塑性滑移速率。该模型实质上是以位错密度为中间变量的应变硬化模型，3.2.2节中将进一步考虑位错的增殖、形核、缠结与湮灭等机制，建立更加完善的位错密度演化方程。

唯象的晶体塑性本构通常存在适用范围的局限性，只能在一定应变率区间内描述材料的塑性变形规律。因此为了满足不同应变率的加载条件，并同时描述各微观结构的演化，通常采用Orowan方程作为流动法则的模型框架

式中：${\rho}_{\rm{m}}^{\alpha } $表示可动位错密度，v^α表示平均位错运动速度。

Orowan方程模型框架主要包含可动位错密度的演化以及平均位错速度与各内变量之间的关系两部分内容。可动位错密度演化部分需要考虑位错的形核、增殖、缠结、湮灭等微观过程。位错平均滑移速度部分可由位错滑移障碍物之间的平均距离除以在障碍物之间运动的平均时间得到。而位错在障碍物之间的平均运动时间可分成跨越障碍物的等待时间以及在障碍物之间的滑移时间两部分，位错的两个运动过程可分别由Arrhenius热激活机制和黏性拖曳机制描述。

位错速度模型描述了位错滑移速度与滑移系分切应力以及位错滑移阻力之间的关系，主要可分为3类：(1)简单指数模型，即基于实验中观察到的单根位错的运动规律，建立位错速度与分切应力的简单指数关系；(2)单一机制模型，即聚焦于位错在低应力区间内或高应力区间内的滑移规律，只采用单一的热激活模型或线性拖曳模型；(3)统一模型，即将位错的滑移过程分为跨越障碍物的等待过程以及在障碍物之间的滑移过程，统一了位错在不同应力区间内的滑移规律。

第1类模型是基于单根位错实验结果的简单拟合，只描述了位错运动速度与应力之间的关系，不仅忽略了位错网络对滑移的阻碍作用，也没有考虑温度对滑移的影响。第2类模型只描述位错在某一应力区间与温度范围内的单一规律，适用范围较窄，缺乏对位错完整滑移过程的描述。第3类模型描述了位错滑移的整个物理过程，适用于不同的应力区间和温度范围。本节将简单介绍前两类模型，重点介绍适用范围最广的第3类模型。

（1）简单指数模型

Johnston与Gilman^[62]利用蚀坑法测量了LiF单晶中单根螺位错与刃位错在不同应力下的位错速度，根据位错运动速度与应力关系的实验结果得到以下结论：其一，刃位错的滑移速度比螺位错更快，但两者的差距在接近声速时逐渐消失；其二，位错运动的最高速度不能超过声速；其三，存在一个位错开始滑移的最低临界应力；其四，在低温低应力下，位错速度与温度近似满足负指数关系。这些结论已被证实在其他材料中同样适用。

基于上述实验结果，Gilman^[63]提出了一种简单的位错速度与应力关系的指数模型

式中：D为特征阻力；v^*为位错运动的极限速度，通常取材料的声速(剪切波速)。

上述简单的指数模型满足实验上位错速度的边界条件：当τ→0时，v→0；当τ→∞时，v→v^*。但该模型不仅缺乏对位错运动机制的解释，也没有反映出位错在低应力区间的转变特征。

为了将模型应用于更广的应力区间，Gilman^[64]提出了在单指数模型的基础上添加一项在低应力区间内的位错速度

式中：${v}_{\rm{s}}^{*} $表示低应力下的位错参考速度，S为低应力下的阻力，${v}_{\rm{D}}^{*} $表示高应力下的位错参考速度。此改进模型更好地解释了位错运动速度在不同应力区间下的转变现象。

此外，Johnson与Barker^[65]相继提出了关于位错运动速度的指数模型，但这些模型均只是根据位错的运动趋势提出的唯象模型，缺乏相应的物理机制，且适用范围较窄。

（2）单一机制模型

单一机制模型包括只适用于低应力区间内的热激活机制以及只适用于高应力区间的线性拖曳机制。Alankar等^[66]、Zhang等^[67]、Monnet等^[68]采用热激活机制来描述位错在低应力区间内的滑移速度

式中：λ为障碍物之间的平均距离，v_d为Debye频率，${\text{Δ}} {G_0}$为参考激活能。此外，Hansen等^[25]、Nguyen等^[69]、Grilli等^[70]均采用线性拖曳模型来描述位错的滑移速度

式中：B为拖曳系数。

单一机制模型描述了位错在某一应力区间内的滑移规律，适用于模拟材料在特定加载区间内的变形过程。但实际上材料在不同时间或空间尺度上的应变率往往存在很大差距，单一机制模型难以描述这种复杂的变形过程。

（3）统一物理模型

目前使用最为广泛的模型基于Frost与Ashby^[71]首先提出的平均位错滑移速度模型框架。该框架将位错在障碍物之间的平均运动时间分成了跨越障碍物的等待时间以及在障碍物之间的滑移时间两部分，可分别采用热激活模型与声子拖曳模型来描述

式中：λ^α表示位错滑移过程中障碍物的平均距离，${t}_{\rm{w}}^{\alpha } $为热激活过程主导的位错脱钉所需要的等待时间，${t}_{\rm{r}}^{\alpha} $为声子拖曳主导的位错在障碍物之间黏性滑移所需要的时间。

对于滑移较为困难的体心立方(bcc)晶体材料，Shahba等^[27]、Lim等^[72]采用螺位错滑移的Kink-Pair机制描述上述位错模型的两个滑移过程。Kink-Pair滑移机制可分为两个阶段，如图6^[27]所示。在第1阶段，位错在障碍物(Peierls能垒)前有一段等待时间，直到成功通过热激活方式跨越障碍物。在第2阶段，位错分离成一对扭折并以黏性滑移的方式运动到下一个障碍物。对于面心立方(fcc)等晶体材料，也可根据与林位错的交互机制将位错滑移分为类似的等待过程与滑移过程两部分。

位错在障碍物前的等待时间${t}_{\rm{w}}^{\alpha } $可用Arrhenius方程形式的热激活模型^[71]描述

式中：ν_d为尝试频率，${\text{Δ}} G\left( \tau \right)$是与分切应力相关的激活能。${\text{Δ}} G\left( \tau \right)$可表述为

式中：阶跃函数〈x〉=(|x|+x)/2；参考激活能${\text{Δ}} {G_0} = {g_0}\mu {b^3}$，g₀为材料的能垒系数；滑移阻力g^α可采用简单的硬化方程式(64)表示。

热激活过程的等待时间以及激活能存在各种不同的表达形式，但这些表达式均表现为等待时间与温度和分切应力满足负相关的关系。

位错在障碍物之间的滑移时间${t}_{\rm{r}}^{\alpha} $取决于位错的黏性滑移速度(${t}_{\rm{r}}^{\alpha} $ = λ^α/${v}_{\rm{r}}^{\alpha} $)，而位错的黏性滑移过程则受声子拖曳效应的影响。声子拖曳效应指位错在高速运动时，位错与声子之间发生的散射作用^{[63, 73]}，位错滑移速率越大，声子对位错滑移的阻碍作用越强。此外，位错在高速运动下也受到其他阻碍作用，例如位错核的辐射作用^[74]。

位错的平均黏性滑移速度通常采用准线性拖曳关系^{[39, 75-76]}描述

式中：B(v,T)是与位错速度和温度相关的黏性拖曳系数，${g}_{\rm{ath}}^{\alpha} $为无热的长程阻力。

Austin等^[75]为了体现位错速度不能超过声速的物理限制，提出了耦合相对论效应的拖曳系数表达式

式中：B₀(T) = k_BT/v_s(b^α)²是与温度线性相关的初始拖曳系数，v表示位错平均速度，v_s为材料的剪切波声速。

位错的滑移速度在不同应力区间内的主导因素不同。对于低应变率加载，位错的运动由热激活过程主导，位错在障碍物之间的黏性滑移时间相对于热激活时间可忽略不计。对于高应变率加载，位错的运动由黏性拖曳过程主导，位错在障碍物前等待的热激活时间相对于黏性滑移时间可忽略不计。通过比较不同应力状态下平均位错运动速度v与热激活速度${v}_{\rm{w}}^{\alpha} $ = λ^α/${t}_{\rm{w}}^{\alpha} $、拖曳速度${v}_{\rm{r}}^{\alpha} $之间的关系，可得到两种主导机制的转变过程，如图7^[27]所示。

位错密度对塑性滑移的影响主要有两方面：一方面，位错滑移受到的平均阻力取决于总的位错密度；另一方面，塑性变形只来源于可动位错密度的贡献。所以通常将总位错密度分为可动位错和不可动位错两部分。材料在准静态下的可动位错密度相对于不可动位错来说很少，并且可近似为常数，在计算塑性滑移速率时可以忽略可动位错密度的演化。因此，在唯象的晶体塑性框架中通常不区分材料中的可动位错和不可动位错^[77]，在分析滑移速率演化时仅考虑位错滑移速度的变化，而分析材料强度变化时仅关注总位错密度的演化。但是Wang等^[78]指出高应变率下材料中的可动位错密度是总位错密度的主要部分，这是由可动位错密度在高应变率下迅速增殖导致的。

材料中可动位错密度ρ_m与不可动位错密度ρ_i的演化过程主要与位错的形核、增殖、缠结以及湮灭等局部微观机制相关

式中：${\dot{\rho }}_{\rm{nuc}}^{\alpha }$、$ {\dot{\rho }}_{\rm{mult}}^{\alpha }$、${\dot{\rho }}_{\rm{trap}}^{\alpha }$、$ {\dot{\rho }}_{\rm{ann}}^{\alpha }$分别为位错的形核速率、增殖速率、缠结速率和湮灭速率。位错的形核过程$ {\dot{\rho }}_{\rm{nuc}}^{\alpha }$包括均匀形核过程$ {\dot{\rho }}_{\rm{hom}}^{\alpha }$和非均匀性形核过程${\dot{\rho }}_{\rm{het}}^{\alpha }$。均匀形核指位错环等缺陷在基体中无择优地均匀成核，通常只在极高的应变率下并且分切应力接近晶格的理论剪切强度下才可能发生，此时高应力状态下位错形核的激活能显著地降低。Tschopp等^[79]、Cawkwell等^[8]、Shehadeh等^[80]均在冲击加载下的分子动力学(MD)模拟中观察到了位错的均匀形核现象，其中位错均匀形核速率的热激活模型为

式中：N表示单位体积内形核点的平均数量，ν_d为特征尝试频率(Debye频率)，${\text{Δ}}G$是与分切应力以及温度相关的位错环形核的自由能。

非均匀形核则指在晶界、第二相界面以及裂纹等应力集中处发生的位错形核现象，非均匀形核所需要的应力远小于晶格的理论剪切强度。Gupta最初在模拟LiF单晶的冲击响应中引入了非均匀形核模型，非均匀形核过程通常与基体中的应力集中现象有关。Austin等^[75]在Gupta模型的基础上提出了适用性更广的非均匀形核模型

式中：α_het为单位体积内位错源形核的平均长度，$\left| {{{\dot \tau }^\alpha }} \right|$为应力率的绝对值，m是用于描述位错源形核的临界应力分布的参数，τ_max和τ_min分别为发生非均匀形核的临界最大应力和最小应力。此外，也有模型将位错的非均匀形核过程考虑在位错的增殖过程$ {\dot{\rho }}_{\rm{mult}}^{\alpha }$中。

位错的增殖机制包括弗兰克-里德源增殖机制、双交滑移增殖机制以及攀移增殖机制等。其中弗兰克-里德位错源增殖机制指位错在与障碍物交互过程中通过外加应力或者高速运动下的惯性效应所形成位错环及其扩展的过程。双交滑移机制指螺位错经过双交滑移后形成了刃型割阶，并对原位错产生钉扎作用，成为一个新的弗兰克-里德源。攀移增殖机制通常指刃位错在攀移后形成了super-jog，并成为两个单臂弗兰克-里德源。这些增殖机制可统一抽象为位错与林位错之间的交互作用^[81-82]

式中：C_m为与各种增殖机制相关的增殖系数，${\rho}_{\rm{for}}^{\alpha} $为与滑移系α发生交互作用的林位错密度。

式中：交互矩阵$ {A}_{\rm{for}}^{\alpha \beta }=\left({{n}}^{\alpha } \cdot {{s}}^{\beta }+{{n}}^{\alpha } \cdot {{q}}^{\beta }\right)/2$在基于刃位错与螺位错数量相等的假设下可用滑移系的滑移方向s^α与滑移面法向n^α表示，q^β = n^β × s^β表示刃位错线的方向，s^β表示螺位错线的方向。

位错的缠结过程代表了可动位错与林位错作用被缠结并转化为不可动位错的过程

式中：C_T为缠结系数。

位错的湮灭指符号相反的位错在一定捕获距离内相遇所发生的消亡现象，决定了位错饱和密度的存在。Luscher等^[39]、Alankar等^[66]根据上述湮灭过程建立了简单的湮灭模型

式中：C_A为湮灭系数，d_a为捕获距离。

为了分析冲击加载下位错各演化机制在不同压力下的主导地位转变，Lloyd等^[83]模拟了5～20 GPa加载下位错各个演化机制的占比情况，如图8^[83]所示。在中低压(5～15 GPa)下，位错的增殖主要依赖于非均匀形核过程$ {\dot{\rho }}_{\rm{het}}^{\alpha }$与增殖过程$ {\dot{\rho }}_{\rm{mult}}^{\alpha }$；当压力大到一定程度(大于20 GPa)时，材料发生剧烈的均匀形核过程，成为位错密度演化的主导因素。

实际上除了可动位错和不可动位错之外，位错密度还有许多其他区分方式。例如，Alankar等^{[66, 84]}为了体现刃位错和螺位错滑移速率的差异，将位错密度分为螺位错ρ_s和刃位错ρ_e，并解释了基面取向的单晶Ti的硬化率转变现象。此外，Keshavarz等^[85]、Liang等^[86]为了考虑几何必须位错引起的晶界效应，将位错区分为统计存储位错ρ_ssd和几何必须位错ρ_GND，分析了在Al双晶、Ni基合金等材料中晶界对材料的硬化作用，并与实验结果较好地吻合。

在高温、高压、高应变率动态冲击条件下，位错滑移往往不再是材料唯一的塑性变形方式，例如Fe、Ti和RDX等晶体材料也可通过相变与孪晶等方式实现变形。固体相变按原子的迁移方式可分为扩散型相变和位移型(切变)相变。而在高应变率动态变形下，材料发生的往往是切变型马氏体相变，整个相变过程无原子扩散，通过原子协同进行短距离迁移，并以切变的形式完成相变。本节重点介绍如何将马氏体相变机制引入晶体塑性有限元框架，以及在动态冲击与静态加载下相变驱动力与演化准则的差异。

当然除了在晶体塑性有限元框架中直接引入马氏体相变机制之外，也可采用相场与动态晶体塑性有限元结合的方法来模拟变形过程中的相变过程。不同于晶体塑性框架中的局部相变模型，相场法在自由能的基础上考虑了相变量在空间的分布及其所引起的梯度能。相场与动态晶体塑性模型相结合的方法已被广泛应用于模拟准静态下材料的相变、孪生与再结晶等变形过程^[87-91]，该方法有望发展至模拟动态冲击下的相变、孪生等塑性变形过程。

关于马氏体相变介观模型的研究由来已久^[92-96]。Thamburaja和Anand^[97]最初将马氏体相变机制引入晶体塑性框架中，以研究Ni-Ti形状记忆合金的宏观力学响应，模型中总变形分为弹性变形和相变变形，忽略了位错滑移产生的塑性变形。在此基础上发展起来的更为广泛使用的模型考虑了塑性滑移，将变形梯度^{[13, 98-99]}分解为

式中：F^e、F^p、F^tr分别表示弹性变形、塑性变形以及相变对总变形梯度的贡献。

塑性变形包括母相与新相中的塑性滑移总和

式中：$ {v}_{\rm{p}}=1-{v}_{\rm{N}}$表示母相的体积分数，$ {v}_{\rm{N}}=\sum\limits_{t=1}^{{N}_{\rm t}}{v}_{t}$为相变后新相的总体积分数，N_t为变体的数量，N_α和N_β分别代表母相和新相的滑移系，$ {{m}}_{\rm{s}}^{\alpha }$和$ {{n}}_{\rm{s}}^{\alpha }$分别表示母相中的滑移方向和滑移面法向，$ {{m}}_{\rm{s}}^{\beta }$和$ {{n}}_{\rm{s}}^{\beta }$分别为新相中的滑移方向和滑移面法向。当相变产生的新相的位错滑移能力较弱时，往往只计及母相中的塑性滑移^[100]。

马氏体相变过程的变形梯度为

式中：b^t和m^t分别为相变过程的切变向量和相变惯习面的法向。

相变模型的关键在于如何描述相变过程的驱动力以及相变的演化规律。本节将重点介绍在准静态和动态冲击下几种常用的相变驱动力与相变演化规律的模型。

相变驱动力是与相变量功共轭的参量，取决于材料相变前后自由能的变化。在现有的研究中通常有两种方法描述相变驱动力与系统自由能的关系。第一种方法较为常见，直接将相变过程中Gibbs自由能的变化作为相变的驱动力。第二种方法以Turteltaub和Suiker^[101-103]为代表，认为相变驱动力在Gibbs自由能的变化之外还应额外加上一项相变熵增所导致的热驱动力，并指出这部分相变熵增并没有引起系统Gibbs自由能的变化。上述两种方法在相变驱动力与Gibbs自由能的关系表达式上存在区别，这一区别实际上来源于自由能定义的具体表达式和物理含义的不同，但是最终得到的相变驱动力是一致的，反映的物理规律是完全相同的。

相变驱动力与各热力学势函数之间的关系不完全相同^[98]

式中：f_t表示与相变量功共轭的相变驱动力，τ_m = F^eTF^eSF^tr–T:b^t$ \otimes $m^t表示相变过程的机械驱动力，ψ、G分别为系统的Helmholtz自由能与Gibbs自由能。

Manchiraju等^[104]为研究Ni-Ti形状记忆合金中的相变过程，给出了系统的亥姆霍兹自由能表达式

式中：θ、θ₀、θ_T分别表示瞬时温度、初始温度和相变温度，λ_T为单位体积的比热，h_tu表示不同相变变体之间的相容关系。式(88)中：第1项表示弹性应变能；第2项表示与热应变对应的自由能；第3项表示相变过程中的化学能；第4项表示马氏体变体之间的相互作用能，反映了变体之间的相互协调能力。

相变驱动力f_t可根据式(87)中与亥姆霍兹自由能的关系得到

式中：S为PKⅡ应力。

此外，Feng等^[105]、Greeff等^[106]认为在冲击高温高压作用下，材料的相变驱动力主要依赖于系统的压强和温度。Cawkwell等^[38]根据材料在高压下的状态方程建立了一个评估两相Gibbs自由能的模型，并给出了含能材料RDX的α与γ相在不同温度与压强下的Gibbs自由能之差，见图9^[38]。该模型简化了相变驱动力的组成部分，忽略了两相的错配能以及相变的界面阻力等各个因素，通常只适用于材料的高压冲击相变过程。

相变过程中新相体积分数的演化依赖于相变的驱动力f_t和相变的临界阻力f_c。通常将相变过程中各种难以量化的阻力集合在一个总的相变临界阻力f_c。因此，马氏体相变启动的临界条件为系统可提供的自由能达到马氏体相变的能垒

Turteltaub等^[101-102]、Tjahjanto等^[107]认为相变启动之后新相的生长速率$ {\dot{v}}_{{t}}$依赖于驱动力与相变临界阻力的差值(f_t – f_c)，并依此提出了一个唯象的相变演化的动力学关系

式中：k为相变演化系数。

Thamburaja等^{[97, 108-110]}、Manchiraju等^[104]则倾向于认为一旦相变驱动力f_t超过临界阻力f_c，马氏体相变过程会立即发生并释放系统的自由能，从而使相变驱动力永远满足自洽关系f_t = f_c，因此新相的生成速率$ {\dot{v}}_{{t}}$可根据相变的自洽关系求解

Barton等^[111]将上述唯象的相变演化动力学关系式(式(91))与固体状态方程相结合，模拟了Fe在高压下的冲击相变$ \left(\alpha \leftrightarrow \varepsilon\right)$，如图10^[111]所示。

此外，考虑动态冲击过程中相变驱动力(Gibbs自由能)依赖于压力和温度，Greeff等^[106]根据实验观察到的在一定压力区间内新相生长速率和压力满足指数关系，提出了适用于冲击高压下的相变演化模型

式中：${\text{Δ}} G\left( p \right)$为两相的Gibbs自由能之差，A、B是与相变演化相关的材料参数。

对于Ti^[112-113]、Zr^[10]等密排六方结构(hcp)材料，在冲击加载过程中，孪生变形与位错滑移是系统发生塑性变形和释放能量的主要途径，而相变往往需要更大的加载压力和应变速率^[105]。孪生变形通常出现在滑移受阻引起的应力集中区域，通过均匀剪切的方式形成镜面对称的两部分晶体。由于孪生的应变速率敏感性通常低于位错滑移的应变速率敏感性^[114]，即孪生阻力随应变率变化不大，材料在高应变率下比准静态下更倾向于发生孪生变形而非塑性滑移。因此，Cu^[12]、Al^{[11, 115]}等滑移系较多的fcc材料在高压冲击过程中也能观察到显著的孪生变形。

Kalidindi^[116]最先将孪生变形引入晶体塑性框架中，模拟了低层错能的多晶fcc金属α-黄铜以及hcp金属Zr中的织构演化。Kalidindi模型沿用了变形梯度的弹塑性分解

而塑性速度梯度包含材料的塑性滑移与孪生变形两部分

对于体积分数通常较小的fcc和bcc孪晶而言，往往只考虑母相中位错滑移而忽略孪晶内部的塑性滑移。但对于孪晶体积分数较大的情况如某些hcp金属，孪晶中的位错滑移则不能忽略。因此，塑性变形梯度中应同时包含母相和孪晶中的位错滑移

式中：$ {f}^{\beta }$表示孪晶的体积分数，$ {{m}}^{\alpha }$和$ {{n}}^{\alpha }$分别表示柏氏矢量与滑移面法向，$ {{Q}}^{\beta }$为孪生变形前后晶体位相的旋转矩阵

式中：$ {{\text{δ}} _{ij}} $为克罗内克符号。

孪生变形本身造成的速度梯度L_tw可由孪晶扩展方向${{m}}_{\rm{tw}}^\beta $与孪晶面法向${{n}}_{\rm{tw}}^\beta $表示

式中：γ_tw为孪生的特征剪应变，fcc和bcc晶体中${\gamma }_{\rm{tw}}=\dfrac{\sqrt{2}}{2}$，hcp晶体中$ {\gamma }_{\rm{tw}}=\dfrac{{c}^{3}/{a}^{3}-3}{\sqrt{3}c/a}$^[117]。

此外，Clayton^[118]采用将孪生变形从塑性变形中独立出来的三构型分解方法

式中：F^e表示可逆的热弹性变形梯度，F^p表示体积守恒的塑性滑移变形梯度，F^I表示孪生变形以及由空位等缺陷引起的不可逆体积变形。

材料的孪生过程与相变过程类似，孪晶的体积分数演化取决于孪生的驱动力和孪生阻力，但晶体塑性模型中孪生模型与相变模型对关键因素的处理却存在差异。相变模型细化了材料在不同情况下的各种相变驱动力，将各种难以量化的相变阻力简单地归纳为一个临界相变阻力常数。而孪生模型的驱动力即为材料内部在孪生方向上的分切应力，其孪生阻力则需要考虑孪晶自身的硬化、滑移和孪生相互作用导致的硬化以及孪晶饱和体积分数等因素。

目前缺乏简洁明确的物理方程以及相关理论模型描述孪晶的演化过程，因此在晶体塑性框架中通常采用率相关的唯象模型作为孪晶的演化法则^[19]

式中：$ {\dot{\gamma }}_{0}$为参考滑移速率，τ^β为孪晶方向上的分切应力，${s_{\rm{tw}}^\beta} $表示孪晶系的孪生阻力。

存在变形孪晶的材料，其应变硬化通常是一个多因素影响的复杂过程。Salem等^[119-120]总结了hcp材料中滑移与孪生的主要硬化效应：（1）变形孪晶对所有滑移系和孪晶系均有硬化效果；（2）孪晶界具有与晶界相似的Hall-Petch效应，对母相中的位错滑移存在硬化作用；（3）Basinski硬化机制(孪晶扩展过程中将可动位错固化缠结)会导致孪晶强度大于基体母相，并提高材料整体强度；（4）孪晶体积分数通常存在饱和值，材料在孪晶体积分数达到饱和后将采用位错滑移作为单一的塑性变形方式。

基于上述分析结果，Kalidindi^[121]提出了饱和型位错滑移和孪生变形的唯象硬化模型，将位错滑移阻力随时间的变化率表示为

式中：h_s为孪晶初始硬化率，a为硬化指数，系数C表示Basinski硬化机制对滑移系的硬化作用，s_pr描述了孪晶界引起的Hall-Pech效应，s_sa0为没有孪生时的滑移饱和阻力。

类似地，孪生阻力的演化过程可表示为

式中：h_tw、h_tw-s分别表示孪晶和位错滑移对孪生阻力的硬化系数，b、d为硬化指数。式(103)中第1项表示孪晶内部的硬化作用，第2项表示位错滑移与孪生过程相互作用引起的硬化。Salem等^[120]指出由于式(103)中第2项的快速硬化作用，孪晶体积分数随着应变增大逐渐达到饱和值。

滑移阻力与孪生阻力的硬化过程实际上取决于材料中不同微观结构的演化以及变形中的应变率和温度^[122-123]。位错在滑移过程中除了受到晶格阻力和林位错的硬化阻力之外，还需考虑晶界阻力的作用，对于晶界阻力(即Hall-Petch效应)，则需要额外考虑孪晶界引起的阻碍作用。Beyerlein等^[124]在位错模型中考虑了林位错硬化、位错碎片(Dislocation Debris)硬化以及Hall-Petch效应等对位错滑移的硬化作用，将滑移阻力分为

式中：${S}_{0}^{\alpha} $为初始滑移阻力，${S}_{\rm{for}}^{\alpha} $为林位错硬化阻力，关于位错密度的演化详见3.2.2节；ρ_deb为位错碎片密度，通常指具有一定分布规律的固定的不可动位错，由于位错碎片分布的非随机性，如式(105)所示，需要对硬化公式中平均障碍距离($ d \approx 1/\sqrt{\rho }$)进行一定的校正^[125]；${S}_{\rm{HP}}^{\alpha} $是由晶界与孪晶界共同引起的Hall-Petch效应阻力，H^α为硬化常数，b^α为柏氏矢量。Ardeljan等^[122]提出Hall-Petch效应中位错的平均自由程d_mfp与孪晶的体积分数相关

式中：f为孪晶的体积分数，d_c为孪晶间距，θ为滑移面与孪晶面之间的夹角，d_g为初始晶粒尺寸。

孪晶的扩展过程与位错滑移过程类似，均受到晶界与位错的阻碍作用。但是不同于位错在晶界处的塞积作用，Hall-Petch效应对孪晶的硬化作用目前并没有明确的物理模型。此外，孪晶扩展与位错之间的相互作用也不再是简单的林位错硬化，而是同时取决于位错密度以及滑移系与孪晶系的柏氏矢量。孪生阻力$S_{{\rm{twin}}}^\beta $可分为初始孪生阻力、Hall-Petch效应以及与位错滑移相互作用3部分的贡献

式中：$S_{0}^\beta $为初始孪生阻力，$S_{{\rm{HP}}}^{\rm t} $为Hall-Petch效应对孪晶系的阻碍作用，$S_{{\rm{slip}}}^\beta $代表位错滑移与孪晶扩展之间的相互作用，C^αβ为位错滑移对孪生的硬化矩阵。位错滑移对孪晶扩展并不只有抑制作用(C^αβ > 0)，在一些特殊情况下，位错滑移能促进孪晶的形核与扩展^[126-128](C^αβ < 0)。

高压、高应变速率冲击加载条件下的材料破坏过程与微观机理等问题，是动态变形与破坏领域的另一个重要问题。层裂和绝热剪切断裂是两种典型的材料动态破坏方式^[52]。

层裂是冲击波经样品自由表面反射产生的反向卸载稀疏波与来自加载面的正向卸载稀疏波相互作用(迎面相遇)，产生拉应力，从而引起材料发生宏观脆性断裂的现象。层裂过程包括微孔洞形核、长大并汇合形成断裂面，微孔洞所在区域存在应变局域化特征。

绝热剪切带是材料在高速动态变形过程中由于局部温升过高引起热软化，在微观上发生剪切变形局域化的现象。由于变形时间极短，变形产热来不及扩散，剪切变形局域化几乎在绝热条件下进行。绝热剪切带表现为细长的窄带，其演变包括萌生、扩展和断裂3个阶段。

材料在动态破坏过程中的微观结构动态演化缺少直接的测量结果，动态晶体塑性有限元方法对动态破坏力学响应和微观结构演化的预测及分析对于相关实验研究尤为重要。关于动态层裂破坏的模拟，主要介绍两种将损伤引入晶体塑性模型的方法；关于绝热剪切带的模拟，分别从忽略损伤和直接引入损伤两种思路对动态晶体塑性模拟方法进行具体介绍。

层裂是材料动态破坏的一种特殊形式。受高速冲击加载时，材料内部两个卸载稀疏波相互作用，一旦材料中局部应力达到材料的动态强度极限即发生断裂，如图11^[52]所示。因此，层裂是结构动态响应（波的交互作用）和材料动态响应（材料动态强度极限或破坏准则）共同作用的结果^[51]。

采用动态晶体塑性有限元对冲击层裂现象的研究主要有两种方法：（1）在动态晶体塑性模型中引入层裂破坏准则，通过显式地删除单元实现材料的损伤形核及其演化；（2）将材料视为有微孔洞的多孔晶体材料，在动态晶体塑性模型中基于孔洞体积在变形过程中的演化引入材料的损伤演化。下面具体介绍这两种动态冲击变形中层裂破坏的模拟方法。

直接引入材料损伤是晶体塑性模拟动态层裂破坏的一个重要方法。依据损伤形核的难易程度，将晶体材料划分为不易损伤的晶内材料和易出现损伤形核点的结构(如晶界、第二相粒子处等)，对两者分别建立不同的材料模型进行研究。对晶内材料建立动态晶体塑性有限元模型，模拟动态冲击压缩过程中的晶内微观塑性变形行为；对晶体材料中易出现损伤形核点的结构，建立各向同性的宏观模型进行动态变形模拟，并结合层裂破坏准则对损伤形核进行定量预测。一旦模型中某个单元满足层裂破坏准则即成为损伤形核点，人为地将相应的单元抹去，即在模型中引入微孔洞。随着动态冲击变形的继续发生，微孔洞将会长大、合并。

材料失效模型通常有两类：内聚应力准则和损伤积累准则^[129]。内聚应力准则假定材料失效是瞬态行为，当材料中的局部应力达到晶体的内聚应力即理论断裂强度时，材料立刻发生断裂。损伤积累准则假定材料失效依赖于变形历史，断裂的发生不仅与应力大小相关，还与应力持续时间相关。现有的动态层裂破坏晶体塑性有限元研究均简单地采用内聚应力准则，考虑到高应变速率变形中存在明显的断裂滞后现象，后续研究有待将损伤积累准则引入动态晶体塑性模型中。下面介绍两种典型的内聚应力准则：最大正拉应力准则和复合模式的最大应力准则。

最大正拉应力准则可表示为^[130]

当材料所受正拉应力(σ)达到临界值(σ_c)时，则发生层裂破坏。Zhang等^[131]、Lloyd等^[132]分别以晶界和第二相粒子作为损伤形核点，采用该准则判断裂纹形核，通过动态晶体塑性有限元模拟冲击变形卸载过程中的层裂破坏现象。

复合模式的最大应力准则可表示为^{[23-24, 57-58]}

式中：s、τ分别为界面上的正应力和切应力，s₀、τ₀分别为参考正应力和参考切应力，${\text{Δ}}T$为变形过程中的温升，s₁、τ₁分别为依赖于温升的临界正应力增量和临界切应力增量。当材料局部正应力或切应力超过对温度线性依赖的内禀界面强度时，均会引起材料的局部裂纹形核，随后由正应力和切应力主导的损伤演变略有不同。Clayton^{[23, 57-58]}、Vogler等^[24]将晶界视为损伤易形核点，采用该准则判断裂纹形核，模拟结果如图12(a)和图12(b)所示。样品右侧区域层裂开始形核，该区域材料承受着较大的拉伸应力，即层裂区域附近具有较大的弹性能。动态层裂破坏表现为晶间断裂，并且飞片的冲击加载速度越大，材料的层裂破坏越强烈，如图12(c)和图12(d)所示。

不同于在无缺陷材料中逐步引入材料损伤(如6.1.1节)，将材料视为有微孔洞的多孔晶体材料，以孔洞的体积变化表征材料微观损伤演变，是研究材料动态层裂破坏的另一个重要思路。

多孔晶体在动态弹塑性大变形过程中，不再满足塑性变形体积不变假设(Ling等^[133])，即det(F^p) ≠ 1，需要考虑孔洞引起的不可逆塑性体积变形。Nguyen等^[134-135]在多孔晶体变形运动学中将塑性变形分为体积变形和等体积变形两部分，将变形梯度表示为

式中：F^e、F^p、F^d分别为弹性变形梯度、塑性等体积变形梯度和塑性体积变形梯度，如图13^[134]所示。通过显式考虑动态冲击变形过程中孔洞引起的不可逆体积变化det(F^d)，可以保证塑性等体积变形依然满足det(F^p) = 1。引入孔隙度(${\varphi } $)状态变量表示孔洞在块体材料中的体积比，用于描述孔洞引起的不可逆体积变形梯度

式中：${\varphi }_0 $ = $V_0^{\rm{voids}}$/V₀、${\varphi } $ = $V_{\rm{d}}^{\rm{voids}} $/V_d分别为参考构型(B₀)和损伤构型(B_d)中的孔隙度，$V_0^{\rm{voids}} $、$V_{\rm{d}}^{\rm{voids}} $分别为参考构型和损伤构型中孔洞的体积，V₀、V_d分别为参考构型和损伤构型中块体材料的体积。同时，假设块体材料中的实心部分经历塑性变形时体积不变，即V₀ – $V_0^{\rm{voids}} $ = V_d – $V_{\rm d}^{\rm{voids}} $。初始孔隙度的下限为零，表示无孔洞的实心材料。

损伤对动态晶体塑性本构方程的影响主要包括3方面：孔隙度对弹性本构的影响、孔隙度对塑性本构的影响以及多孔晶体中孔隙度本身的演变。首先，孔隙度对弹性本构的影响主要表现为弹性模量对孔隙度的依赖^[134-135]

式中：C、C^d分别为无损伤和有损伤材料的弹性模量，E^e、S分别为格林-拉格朗日应变和PKⅡ应力。随动态冲击变形过程中孔隙度${\varphi} $的增大，多孔晶体材料的弹性模量C^d不断减小。其次，孔隙度对塑性本构的影响主要表现为塑性滑移速率对孔隙度及其变化率的依赖

式中：与孔隙度无关的塑性滑移速率$ {\dot{\varGamma }}^{\alpha }$由速度梯度的塑性分量确定($ {{{L}}}^{\rm{p}}=\sum\limits_{\alpha }{\dot{\varGamma }}^{\alpha }{{s}}_{0}^{\alpha }\otimes {{m}}_{0}^{\alpha }$)；M为泰勒因子；n_slip为滑移系数目；$ {\dot{\gamma }}_{\rm{eq}}$是与微观尺度晶格振动相关的各向同性等效塑性应变速率；孔隙度${\varphi } $是孔洞半径(a)和相邻孔洞中心距(l)之比的函数，${\varphi } $ = (a/l)³；$ \dot{\varphi }$为孔隙度随时间的变化率。式(117)实质上是以$ {\dot{\gamma }}_{\rm{eq}}$作为桥梁，建立各向同性的等效塑性应变速率$ {\dot{\gamma }}_{\rm{eq}}$与孔隙度的关系，对与孔隙度无关的塑性滑移速率$ {\dot{\varGamma }}^{\alpha }$进行修正，从而建立滑移系α上塑性滑移速率${\dot{\gamma }}^{\alpha }$与孔隙度的关系。最后，基于孔洞径向动量守恒方程，得到孔隙度随时间的演变方程

式中：ρ为基体材料密度，l为相邻孔洞之间的中心距，Σ_m为平均宏观柯西应力，R为材料的等体积变形阻力。

上述方法能够有效地对多孔晶体中孔洞的长大进行模拟，Ling等^[133]通过比较多孔晶体塑性模拟得到的横向应变和纵向应变的相对大小判定孔洞合并。他们发现：一方面某些取向的单晶材料在应力三轴度很小且应变很高时仍无孔洞合并发生，这显然与实际不符；另一方面孔洞合并临界应变依赖于晶粒取向和应力三轴度，而临界孔隙度仅依赖于晶粒取向，如图14^[133]所示。因此，在动态多孔晶体塑性模拟中引入孔洞合并准则是有必要的，相比临界应变，采用仅依赖于晶粒取向的临界孔隙度作为判据更为合适。

绝热剪切带是材料在动态高应变速率加载变形过程中常发生的剪切变形局域化现象，最终会导致材料动态断裂^[136-137]。Zener和Hollomon^[138]最早提出绝热剪切带的形成是热软化克服材料加工硬化的结果，并将微观绝热剪切带的形成与材料宏观热塑性软化(即变形后期应力-应变曲线发生下降)联系在一起，见图15^[51]中Stage 3。绝热剪切带的基本特征包括：微观上观察到的绝热剪切带、宏观应力-应变曲线上发生热塑性失稳、高速变形过程接近绝热过程^[51]。剪切局域化机理包含两方面竞争因素：促进形成剪切局域化的正面因素(即热软化和几何软化)和不利于剪切局域化的负面因素(即应变硬化和应变率硬化)^[52]。

高速冲击变形过程中绝热剪切带的晶体塑性模拟研究主要有两种方法：（1）忽略损伤，直接模拟冲击变形过程中的塑性变形和温度演化，根据模拟结果在后处理中观测和分析绝热剪切带现象；（2）引入材料微观损伤，模拟动态冲击变形过程中绝热剪切带的形核及演化。下面分别对这两种方法进行具体介绍。

高速变形下绝热剪切的动态晶体塑性有限元研究可以追溯到20世纪90年代，采用忽略损伤的动态晶体塑性有限元能够直接对高温、高压、高应变速率下绝热剪切带的形成及演变进行模拟。Hines等^[139]、Lee等^[140]基于单晶体弹塑性理论，忽略热效应，采用经典的应变率相关的唯象硬化本构方程，粗略地模拟了动态加载下的剪切带现象。Bronkhorst等^[141]、Tajalli等^[59]考虑塑性变形产热以及热软化效应，基于高速变形中的绝热假设进行热力耦合计算，采用唯象的动态晶体塑性有限元模拟分析绝热剪切局域化现象。考虑到动态变形温升显著并且对材料弹塑性本构具有重要影响^[142]，Bargmann等^[143]进一步考虑热弹性耦合效应和热传导效应，对动态变形温升进行更加精确的计算。他们同时考虑塑性产热、弹性产热以及热传导对温升的贡献，采用动态梯度晶体塑性有限元模拟10⁴～10⁷ s⁻¹应变速率范围的材料冲击变形行为。通过比较不同应变速率、相同变形程度的模拟结果，发现相比于高应变速率(10⁷ s⁻¹)的动态变形，材料在相对较低应变速率(5×10⁵ s⁻¹)下的动态变形过程中热传导效果明显，其温度梯度更为平缓，如图16所示。

采用上述动态晶体塑性有限元模型已经能够模拟绝热剪切带的形成，根据模拟得到的应力-应变曲线上热塑性失稳转折点的力学状态(如应变、应变率、剪切模量或变形能等)，进一步给出热塑性失稳等效判据，能够有效地对微观上可观察到的绝热剪切带形成的临界点进行预测和分析。当材料受到动态高速载荷时，虽然应力在屈服点附近可能存在短暂下降，但是宏观应力-应变曲线整体上首先经历一段上升，随后曲线发生持续性下降，如图15所示。通常认为曲线上峰值转折点(即dσ/dε = 0或dτ/dγ = 0)处宏观上开始发生热塑性失稳，同时意味着微观上开始形成绝热剪切带。发生宏观热塑性失稳时，除宏观应力之外，其他的变形或状态参量也会达到相应的临界值，相应地就有其他等效判据，共有4种：临界应变判据^[144-147]、临界应变率判据^[148]、临界剪切模量判据^[149-151]、临界变形能判据^[152]。Zhang等^[147]通过动态晶体塑性有限元模拟得到的等效塑性应变云图上出现了绝热剪切带，如图17所示。进一步采用临界应变判据能够预测宏观热塑性失稳临界点，即微观绝热剪切带形成的临界状态。他们根据模拟得到的宏观应力-应变曲线，对热塑性失稳转折点处的应变状态进行统计(见图18)，发现动态压缩过程中形成绝热剪切带的临界应变约4.5%，动态剪切过程中的临界应变约为3.7%。

从微观角度考虑损伤的形核和演变机制是研究动态绝热剪切破坏的另一种重要方法。将损伤作为材料状态变量引入动态晶体塑性模型时，需要考虑3个方面：引入损伤的弹性本构、损伤驱动力以及变形过程中损伤的演变。Boubaker等^[153]在超弹性本构方程中引入材料损伤，将弹性自由能视为弹性应变、温度和损伤变量的函数，将塑性自由能视为滑移速率和损伤变量的函数，将应力对损伤的依赖归因于弹性模量随损伤变量的变化，从而得到依赖于损伤的自由能和应力

式中：ψ为比自由能，ψ^e、ψ^p分别为比自由能的弹性和塑性分量，S、E^e分别为PKⅡ应力和弹性格林-拉格朗日应变张量，T为瞬态温度，D为损伤变量，$ {\dot{\gamma }}^{\alpha }$为滑移系α上的塑性滑移速率，ρ₀为初始材料密度，$ {{{C}}}^{*}\left(T,D\right)$为依赖于温度和损伤的弹性模量，${{\alpha }}{\text{Δ}} T$表示热膨胀应变。这里与式(116)相比增加了热膨胀构型，因此在应力表达式中采用了纯弹性应变，详见2.2.2节。其次，热力学损伤驱动力可以通过自由能对损伤的依赖性得到

式中：Y为损伤驱动力，g^α为滑移系α上的滑移阻力。显然，损伤驱动力和损伤变量是功共轭的。最后，损伤变量在动态变形过程中的演变可以通过损伤驱动力与损伤演化阻力比值的幂函数表示

式中：Y_r为参考损伤演化阻力，L^p为速度梯度张量的塑性分量。在初始未变形条件下，损伤驱动力和损伤变量均为零；当材料受冲击载荷开始变形，则存在损伤驱动力，损伤变量不再为零，意味着材料发生损伤形核；在随后的变形过程中，损伤变量按式(121)随时间而变化，表示材料中损伤的演变规律。通过这样的方式，将连续介质损伤模型与晶体塑性模型进行耦合，能够有效地模拟动态绝热剪切破坏现象。

动态晶体塑性有限元对于研究高温、高压、高应变速率加载条件下的材料动态冲击力学响应和微观结构演化具有重要作用，但是如何对材料在动态冲击下微观变形的物理机制进行准确的描述仍然存在许多的发展空间。

在理论模型方面，动态晶体塑性有限元能够有效地反映晶体塑性滑移的各向异性、对微观组织的依赖性以及考虑动态冲击下各种影响变形的宏观因素与微观机制。与准静态晶体塑性模型相比，动态晶体塑性模型考虑了材料在高压下的非线性弹性响应、显著温升引起的热弹性耦合效应以及高速运动位错具有的声子拖曳效应等。但是目前对于冲击下位错高速滑移过程中位错之间的相互作用与位错形核、增殖的演化规律缺乏足够的认识，在未来的研究中可基于多尺度模拟框架，从分子动力学与位错动力学出发，探究位错在高速滑移的滑移特征和演化规律。

在应用研究方面，动态晶体塑性有限元能够应用于位错滑移、相变、孪晶和动态损伤破坏等变形与破坏机制研究中。现有的动态晶体塑性模型通过引入局部相变体积分数或孪晶体积分数及其演化规律，研究相变与孪晶对材料变形的影响，但这种方法只能均匀化地将相变与孪晶的影响简化为对材料点应力的影响，无法反映真实的相变过程与孪晶演化过程。在未来研究中可考虑将动态晶体塑性模型与相场模型相结合描述动态冲击下的相变和孪生过程。除了将动态晶体塑性模型应用于动态损伤破坏的定性研究，还可分别结合宏观破坏准则和微观损伤形核及演化准则，对动态冲击变形过程中的层裂和绝热剪切带进行定量预测和分析。从微观损伤演化的角度，对动态层裂和绝热剪切破坏的动态晶体塑性研究尚处于发展阶段，目前的损伤主要以引入损伤系数进行唯象的模拟为主，后续研究中需要更多地考虑基于物理的晶体动态损伤机制。