虽然1943年就发现晶体材料的塑性变形是由各个滑移系上位错滑移引起的^[47–50]，但力学有限元模拟方法仍长期只被用于各向同性材料模型。直到1982年，Peirce等^[51]第一次利用CPFE研究单晶的拉伸行为，且模型只包含2个对称滑移系。随后Harren等^[52–53]将其拓展到包含2～3个滑移系的二维模型。1991年Becker针对三维的柱状多晶^[54]和单晶^[55]，首次考虑了面心立方（Face-Centered Cubic，FCC）晶体的全部12个滑移系。随着计算能力的提升，越来越多的CPFE模型被用于研究复杂二维和三维晶体结构在晶粒或亚晶粒尺度上的变形问题^{[31, 56–62]}。在宏观尺度，均匀化处理方式可使CPFE考虑金属成形过程中的织构演化^[63–67]。

然而，上述CPFE大多数利用唯象型的本构方程，且将位错滑移视作唯一的变形机制。当遇到小尺度变形或者有界面效应时，这种框架下的CPFE将失效。另外，对于应用很广的孪晶塑性（Twinning-Induced Plasticity，TWIP）钢和相变塑性（Transformation-Induced Plasticity，TRIP）钢，除了位错滑移外，还必须考虑增加孪晶和相变所致的变形和强化机制。尺寸效应最初被引入CPFE是通过唯象型的应变梯度理论实现的^[68–70]。由于应变梯度可以通过几何必须位错(Geometrically Necessary Dislocations，GNDs)来描述，故新的内变量本构方程用位错密度代替应变^{[6, 13–14, 19–20, 71]}，使模型具有物理意义。最近的CPFE都能够引进新的变形机制，如晶界机制^{[6, 13–14, 19–20, 71]}、破坏机制^[72]、TWIP钢和TRIP钢里面的孪晶和相变^{[6, 13–14,19–20, 27, 71, 73–74]}。

本节将简单介绍CPFE的基础原理和方法：首先通过连续介质力学联系变形与微观剪切，进而将本构方程分为两类进行描述，最后是其他物理机制的引入。具体而详尽的阐述参考Roters等^[33]的综述。本文主要关注材料动态变形中的塑性机制，对有限元的数值处理方法不作详细介绍。

在连续力学中，利用变形梯度${{F}}$联系当前构型${{x}}$和参照构型${{X}}$，即${\rm d}{{x}} = {{F}}{\rm d}{{X}}$。可以引入变形中间态对整个变形进行分解，常用的分解方式为

式中：${{{F}}_{\rm{e}}}$为弹性变形梯度，${{{F}}_{\rm{p}}}$为塑性变形梯度。为了对整个变形的动力学进行描述，必须引进一个描述变形梯度关于时间的变化率的物理量。变形速度可以表示为 ${{v}} = {\dot { u}}$，其中 ${{u}} = {{x}} - {{X}}$。定义速度的梯度张量为

式中：${{{L}}_{\rm{e}}} = {{\dot { F}}_{\rm{e}}}{{F}}_{\rm{e}}^{ - 1}$为弹性速度梯度，${{{L}}_{\rm{p}}} = {{\dot { F}}_{\rm{p}}}{{F}}_{\rm{p}}^{ - 1}$为塑性速度梯度。

如果塑性变形基于位错滑移，则塑性速度梯度可以表示为

式中：${{{m}}^{\alpha} }$ 和 ${{{n}}^{\alpha} }$ 分别为 $ {\alpha}$ 滑移系的滑移方向和滑移面法向的单位矢量，${{\dot \gamma} ^{\alpha} }$为此滑移系的剪切率，n为激活滑移系的个数。这样就把宏观变形和微观晶向剪切联系起来。注意到此时已经将滑移面上的位错滑移均匀化，才能利用一个剪切率来描述。

本构方程主要寻求外加应力状态与剪切率的关系。现有本构方程可以分为两类——唯象型和物理型，具体分类标准基于模型采用的假设条件，例如可以粗略地以是否采用位错密度作为内变量为判断基准。

塑性模型一般利用临界分切应力 $\tau _{\rm{c}}^{\alpha} $ 表示滑移面 $\alpha $ 的状态，即此面上的滑移阻力。剪切滑移率通常表示为施加在滑移上的分切应力（Schmid应力）${\tau ^{\alpha} }$与临界分切应力的函数

式中：Schmid应力表示为 ${\tau ^{\alpha} } = {{\sigma }}:{{{m}}^{\alpha} } \otimes {{{n}}^{\alpha} }$，${{\sigma }}$ 为Cauchy应力张量。此方程描述流变规则。表征材料状态的 $\tau _{\rm{c}}^{\alpha} $是剪切 $\gamma $ 和剪切率 ${\dot \gamma} $ 的函数

此方程即为硬化规则。必须注意到（5）式中剪切和剪切率没有加上标以代表特定的面，它们是总体的剪切和剪切率。

应用最广的流变准则是Rice^[11]、Hutchinson^[75]以及Peirce等^{[51, 76]}针对FCC金属体系提出的幂律形式

式中：${{\dot \gamma} _0}$和 $m$ 是材料参数，分别为参考剪切率和滑移的率敏感系数；sgn为符号函数；$m \to 0$和$m \to 1$为两个特殊状态，对应率无关和黏弹性。最常用的硬化模型为

其中 ${h_{\alpha \beta }}$ 为硬化矩阵

式中：${h_0}$、 $a$ 和 ${\tau _{\rm{s}}}$ 均为滑移硬化系数，是由所有滑移系所确定的值，基于滑移系内位错反应特征；${q_{\alpha \beta }}$ 是对潜在硬化系数的估计值，对于共面的滑移系，其取值为1.0，其他情况则取值为1.4。此关系能够经验地描述不同滑移系之间的相互作用。

在大多数文献中，都将（6）式和（7）式作为模型的流变规则和硬化关系。但也有文献利用sinh函数替代（6）式中的幂律函数^[54]，有些则利用更一般的Voce方程代替（7）式^{[39, 59, 77]}。在这些唯象型动力学方程组中，只采用了唯一的临界分切应力来描述材料的状态，而非晶格缺陷。因此，基于物理机制的本构方程在需求的推动下开始发展，以切实考虑微结构对变形的影响。

不同于唯象型本构方程，基于物理的模型依赖于内变量。由于位错是承载塑性变形的主要方式，因此位错密度是最重要的微结构变量。已经有大量研究者以位错密度演化来计算流变应力^{[6, 14, 17,20–21, 71]}。下面介绍Ma等^[18–20]建立的应用最广泛的一种基于位错密度的模型。

可动位错密度 $\rho _{\rm{m}}^{\alpha} $ 被用于代表滑移面 $ {\alpha}$ 上承担的部分外载塑性变形。为了可动，这些位错必须克服另两种位错的应力场：平行位错引起的穿越应力以及林位错的交错阻碍。所以需要引进另两种位错密度：与滑移面平行的平行位错密度 $\rho _{\rm{P}}^{\alpha} $ 和与滑移面垂直的林位错密度 $\rho _{\rm{F}}^{\alpha} $。基于晶体中总的不可动位错密度 $\rho _{\rm{SSD}}^{\alpha} $，可以得到

式中：${\chi ^{\alpha \beta }}$为不同滑移面间相互作用强度，包括自作用强度、共面作用强度、交滑移强度、可动结强度、Hirth锁强度以及Lomer-Cottrell锁强度。在这些方程中，只有刃位错被视为无面外速度，其线方向为${{{t}}^\beta }$。在基于位错密度的模型中，Orowan方程被用来替代（6）式，因为它联系了剪切率与可动位错，能够将连续力学和位错物理性质耦合在一起。其方程为

式中：$b$ 为位错Burgers矢量的大小，$v$ 为可动位错的平均速度。可动位错密度可利用简单的标度率给出

其中T为绝对温度，${B^*}$由下式给出

式中：${k_{\rm{B}}}$为Boltzmann常数，$G$为剪切模量，${c_1}$～${c_3}$为位错密度相关的常数。假设林位错交错是率相关的过程，则可动位错速度可以定义为

式中：${\lambda ^{\alpha} }$为与林位错间距成比例的浮动距离；${Q_{\rm{slip}}}$为位错滑移的有效激活能；${V^{\alpha} } = {c_3}{\lambda ^{\alpha} }{b^2}$为激活体积；${\nu _{\rm{att}}}$为接触频率；$\tau _{\rm{eff}}^{\alpha} $是有效剪切应力，可以根据分剪切应力和穿越应力计算

唯象型模型中的硬化规则可以利用位错密度的演化代替。为此，需要建立以单独位错反应为基础的率方程。在此应用广泛的模型中，位错密度的增加进程通过位错锁或位错对的生成而实现，恢复过程则通过非热或热的激活湮灭。

在此简写一些率方程，如位错锁形成和非热激活湮灭分别为

式中：c₄和c₅为拟合常数。

上述基于局部位错密度的模型在处理多晶织构演化等方面非常有效^[78]。然而，对于小尺度的实验现象如纳米压痕等，局域化模型很难处理其中的尺度效应。流变应力对晶粒尺寸的依赖关系最初由Hall-Petch经验方程描述^[79–80]。许多研究都给出了细晶强化原理，基本都是由于晶界周围形成大量塑性变形，可动位错堆积在晶界处引起滑移阻力的上升和界面处的应变梯度。应变梯度被认为是通过增加位错密度而增加滑移阻力^[15]。如前所述，应变梯度与GNDs有关，唯象型模型无法处理GNDs，而基于位错的物理模型则可以轻易将GNDs当作本构框架的一部分。但引入GNDs有一个难点，即应变梯度为非局域化理论，加入GNDs将会增加方程的自由度和额外的边界条件，这极大地增加了特别是复合加载模式下的模型复杂性。Ma等^[19]给出了一种可用且统一的方案。Nye位错张量^[81]可以将应变密度转化为GNDs

注意梯度符号是作用在参考构型上的。据此位错张量，GNDs密度的变化率为

本构方程中，简单直接地将GNDs加入林位错和平行位错之中。为了更加方便，GNDs被分为3类：切向量平行于滑移方向${{{d}}^{\alpha} }$的螺位错、切向量分别平行于滑移面法向${{{n}}^{\alpha} }$和线方向${{{t}}^{\alpha} }$的刃位错。其表达式分别为

它们满足${\left( {{\dot \rho} _{\rm{GND}}^{\alpha} } \right)^2} = {\left( {{\dot \rho} _{\rm{GNDs}}^{\alpha} } \right)^2} + {\left( {{\dot \rho} _{\rm{GNDet}}^{\alpha} } \right)^2} + {\left( {{\dot \rho} _{\rm{GNDen}}^{\alpha} } \right)^2}$。将其加入林位错和平行位错之中

这样一来，引入GNDs描述应变梯度引起尺寸效应的模型已经建立。

晶界能够作为位错运动的有效障碍，可动位错会在界面前端累积并造成应力集中。这些位错-界面相互作用机制很难被唯象型的模型引入；但在处理较大尺度的塑性应变时，可以利用均匀化处理手段将其引入晶体塑性理论。将晶界机制引入基于位错的本构CPFE方法主要有两种：一种将晶界视为位错不可穿越的界面^[16]，一种将其视为部分“透明”界面^[20]。位错不可穿越晶界的假设其实即为添加晶界处的边界条件，虽然能够增强硬化，但并不能预测屈服应力的上升，即不能捕捉Hall-Petch现象。为了应对这个问题，Evers等^[16]引入了晶界位错(Grain Boundary Dislocations，GBDs)表示晶界处初始的GNDs，其密度与晶界两边的晶向以及滑移系有关。另一种处理方式是将位错穿越界面视作一个热激活的概率过程，此热激活能作为可动位错滑移激活能的附加部分^[20]。每个穿越事件将发生在最小可能能量假设的基础上。然而，在任意穿越事件中，进入界面的位错未必等于射出的位错，从而导致界面失配位错。这些失配位错会成为额外的障碍，必须再加入额外的穿越障碍能。一般处理方式并非严格地一对一校准此能量，更多时候利用统计的方法将失配位错视为刚化效应，具体的例子为

（25）式与（14）式唯一的不同是用有效激活能代替了滑移激活能。此时有效激活能为$Q_{\rm{eff}}^{\alpha} = Q_{\rm{slip}}^{} + Q_{\rm{GB}}^{\alpha} $，而晶界激活能则为寻求晶界射出滑移系中的最小值。对于FCC中的[111]和[110]扭转晶界的激活能，可以参见文献[33]中的图表。

以上介绍的都是基于位错机制的CPFE模型。然而，材料的塑性变形并不仅仅由位错滑移所致，还有位移型相变，例如前述的TRIP钢、TWIP钢等。这些相变与位错类似，会产生剪切，一般情况下主要讨论两种：马氏体相变^[22–23]和孪生变形^[24–27]。马氏体相变是剪切所致晶体结构发生改变，从而引起体积的变化；孪生变形则是剪切导致与基体晶向成镜面关系的孪晶。

1965年Greenwood和Johnson^[82]开始研究钢中的TRIP效应，发现体积膨胀引起的额外塑性变形是由于奥氏体向马氏体发生转变所致。同年，Patel和Cohen^[83]发现这些相变倾向于发生在相变驱动力最大的方向上。此后，大量关于马氏体相变的本构关系被提出和发展，但大多数模型都是基于小变形框架下，这些强制假设条件会产生较大的误差，特别是单晶中的各向异性无法体现。基于晶向热力学的模型，可以考虑多相TRIP钢中的变形行为。Suiker和Turteltaub^{[74, 84–85]}利用这种相变模型处理FCC奥氏体向体心四方马氏体的转变。

在相变机制的CPFE中，总的变形梯度不再只有弹性和塑性部分，还有相变部分${{{F}}_{\rm{tr}}}$，即

相变部分的变形梯度为

式中：${{{b}}^i}$为相变系统形状应变向量，${{{d}}^i}$为其惯习面的法向单位矢量，${\xi ^i}$为其相对于参考构型的体积分数。体积分数的演化与相变驱动力${f^i}$有关

式中：$\dot \xi _0^i$ 和 $\nu $ 为材料相关的参数，决定相变发生的频率；$f_{\rm{cr}}^i$ 为相变的临界驱动力；总驱动力${f^i}$包括弹塑性应变、热、缺陷以及表面能引起的驱动力。这些驱动力与滑移状态及外力有关，可以类似位错机制的本构关系通过一定的硬化规则表达出来。

孪晶变形是另一种承载塑性变形的有效方式。一个孪晶可以看作是母相中被剪切部分，该部分晶格取向正好与母相成镜像对称。孪晶也可看作惯习面上的非方向性剪切。在此框架下，FCC里面的孪晶面主要为$ \left\{ 111\right\}\left\langle {112} \right\rangle $类型，体心立方（Body-Centered Cubic，BCC）晶体为$\left\{112\right\}\left\langle {111} \right\rangle $，密排六方（Hexagonal Close-Packed，HCP）晶体则为$\{ {10\bar 12}\}\langle {10\bar 11}\rangle $。孪晶变形的许多物理机制仍不清楚，但其影响因素已被很多研究给出，如温度、应变率、晶粒尺寸和层错能等。由于孪晶变形在一些材料中处于重要地位，故将孪晶机制加入CPFE框架是非常必要的。最早将唯象型孪晶机制加入CPFE的是Doquet^[86]，尔后有Schlogl等^[87]和Mecking等^[88]。

一般可以利用转动矩阵 ${{Q}}$表示孪晶从初始构型变为新构型

式中：${{n}}$为孪晶面法向。对于总的速度梯度，还包含孪晶产生的部分

式中：${f^\beta }$为 $\beta $孪晶相的体积分数，${n_{\rm{t}}}$为孪晶系的个数，${\gamma _{\rm{t}}}$为本征孪晶剪切应变。然而（30）式并不包含孪晶相中的滑移，为此Kalidindi^[24]对其进行了改写

由于现在孪晶的物理机制尚不清晰，故只有唯象型的模型描述孪晶体积分数的演化

孪晶演化符合幂律规则，类似于滑移剪切的幂律，参数的定义也使用同样的方式。硬化规则为

注意到这里的${\tilde \alpha} $表示只有不与滑移面共面的孪晶才会对位错形成阻碍。

连续的损伤断裂力学已经能够处理一些损伤的演化，但这些方法基于已经知晓损伤的形核位置，例如预置裂纹和空洞。当损伤点在材料中随机分布时，各向同性的有限元方法能够预测剪切局域化导致损伤的排布。虽然连续力学方法能够描述特定的损伤点，但一般在不模拟损伤的情况下引入材料的本构模型，比如常用的模型即为损伤降低材料的弹性模量。在微结构层次上，已有少量CPFE从晶粒晶界的晶向性出发来模拟裂纹的扩展和孔洞的粗化^[89–90]。一般情况下，CPFE能够直接模拟实验观察到的损伤出现的特征微结构，并评估损伤形核和初始生长理论。具体微结构所致损伤的模拟主要包含：非均质塑性变形、界面、内聚区边界、晶界穿越、实验的断裂始发准则、损伤的应变能准则。需要考虑损伤断裂的典型材料为多孔晶体，图1所示为引入损伤机制的变形梯度张量分解示意图。

CPFE模拟材料在外加力载荷作用下的响应行为，无论是准静态、低应变率还是爆炸冲击，其基础还是单晶塑性模型的力学本构行为。为此现在以本构方程是唯象型还是物理型进行分类，再对其在材料动态响应方面的应用加以介绍。由于这些本构方程只针对位错滑移，故最后再介绍存在其他机制的模型在材料动态行为模拟方面的进展。

描述唯象型流变规则的方程(6)式本来就体现了剪切率(应变率)的效应关系。当率敏感因子为零时，即为率无关(Rate-independent)的本构关系，也称刚性塑性关系；当率敏感因子非零时，即为率相关(Rate-dependent)本构关系，也称黏塑性关系。与此流变规则方程类似的经典宏观本构模型已经能够考虑应力-应变关系的应变率效应，但直到有限元方法的应用才涉及晶向相关各向异性的微结构特征。Peirce等^[51]最早利用有限元对率相关和率无关的单晶本构方程进行分析，发现剪切局域化受率敏感系数、应变硬化和潜在硬化因素的影响，然而这些结果都是在准静态加载条件下得到的。Zikry和Nemat-Nasser^[92]利用率相关方程研究了拉伸应变率100～2000 s^–1加载条件下单晶的剪切局域化特征，数值处理方面用到了特殊的五阶Runge-Kutta方法，其结果也显示率敏感性和应变率硬化对应变局域化有重要的影响。

动态条件下，温度是影响材料响应的重要因素。Schoenfeld^[93]利用单晶模型研究了多晶钽（Ta）的应变率和温度效应，其应变率在${10^3}\;{{\rm{s}}^{ - 1}}$量级，为早期应用CPFE研究高应变率响应的代表性工作。此研究模型虽然仍采用(6)式的幂律流变规则，但滑移阻力除了(7)式硬化率提供的部分外，还包括Peierls障碍，这些阻力项用一个热激活相关的系数联系起来

式中：$\tau _{\rm{a}}^{\alpha} $为远程分阻力；$\tau _{\rm{p}}^{\alpha} $为Peierls力；$\dot \tau _{\rm{\gamma}} ^{\alpha} = \displaystyle\sum\limits_{\beta = 1}^n {{h_{\alpha \beta }}\left| {{{{\dot \gamma} }^\beta }} \right|} $为硬化率相关的阻力；${S_{\rm{p}}}$和${S_{\rm{\gamma}} }$是关于温度T和剪切率${{\dot \gamma} ^{\alpha} }$的函数，为热激活相关系数。引入Peierls力的本构方程能够很好地模拟Ta在不同温度下的动态响应规律，原因是Ta为BCC结构。与FCC晶体中的平面位错结构不同(Peierls力非常小)，BCC晶体中的很多位错呈现非平面结构，导致其位错的启动阻力很大，因而必须加以考虑这部分滑移阻力。此后，BCC结构材料的温度和应变率相关CPFE模型中都考虑了这种非平面结构导致的滑移阻力^[94–96]，即non-Schmid效应。

由于具有简单实用的优势，唯象型本构方程仍不断被改进和完善，用以描述不同材料的响应特性。对于高应变率的加载情况，有些模型由低应变率下的方程直接外推得到，有些根据实验数据拟合得到，有些则是从物理机制上进行部分改进。近年来，这些应用于材料动态响应的CPFE模型主要关注极端条件下材料的非线性弹性行为以及高速位错的拖曳机制。由于高应变率下变形严重局部化，所以发展不同参照系框架下的模型也非常重要。下面主要从这3个方面加以叙述。

鉴于以前的CPFE模型均应用在准静态或者应变率低于 ${10^4}\;{{\rm{s}}^{ - 1}}$的情况，Becker^[97]将其应用到爆炸冲击条件下的高应力高应变率情况，其中应力高达10～100 GPa，应变率超过${10^7}\;{{\rm{s}}^{ - 1}}$。在此晶体塑性模型中，采用了率无关的流变规则，即流变应力由屈服准则决定

硬化规则未采用（7）式，而是根据幂律硬化格式给出了应变率敏感系数相关的滑移阻力表达式

式中：$\mu /{\mu _0}$为相对剪切模量，是由于压力引起的模量变化；$\tau _0^{}$为参考剪切障碍力；$B'{\text{、}}\!\!k{\text{、}}\!\!l$为材料参数；${\bar \gamma} $为时间平均的剪切；${{\dot \gamma} _{\rm{off}}}$为零应变时的偏移量；$\exp \left( { - \lambda \Delta E} \right)$为变形热效应。此模型真正的重点在于弹性变形部分，如（36）式中的剪切模量，有效地考虑了随应力增加而增加的弹性常数。弹性常数的变化采用Simmons和Wang^[98]的结果，即弹性常数与压力成正比，${\rm{d}}{C_{{ij}}}/{\rm{d}}p = c_{{ij}}^*$。图2给出了此模型计算的压力脉冲沿垂直于表面的方向传播的情况，结果显示波前的应力梯度有利于应变的重新分布。

连续尺度的单晶模型都视弹性变形为线性的小量，然而在极端的高压力环境下，材料的可逆变形并不小。Becker^[97]重点考察了弹性模量在高压力下增大的情况，其本质就是考虑了非线性弹性变形。Luscher等^[99]更进一步将完整的p-V-T状态方程融合进晶体塑性模型，利用经典的内状态变量理论，能够有效应对热-弹性耦合（$\Delta {T_{\rm{r}}}$）、冲击热耗散（$\Delta {T_{\rm{d}}}$）、塑性热耗散（$\Delta {T_{\rm{p}}}$）等导致的温度改变。在塑性本构关系上，必须将热激活因素引入(6)式以改进流变准则

式中：E为热障碍能，$\tau _{\rm{c}}^{\alpha} $由硬化方程(7)式决定，$\tau _{\rm{cl}}^{\alpha} $ 为$\alpha $ 面上的本征晶格阻碍力。结果发现，冲击下的非线性热弹性响应与塑性变形同等重要，如图3所示，一维冲击的波前温度${T_1}$及各个温差随压力变化巨大。将模型用于描述3～110 GPa下单晶铜的变形响应，与实验结果吻合得很好。

在准静态或低应变率情况下，位错速度较低，位错运动主要由热激活主导；在高应变率下，高速位错运动会受到黏性拖曳，如声子拖曳。由于切实关系到位错的运动机制，所以应该考虑用物理型的本构关系来引入拖曳效应，这些将在3.2节加以叙述。然而，物理型本构方程涉及的物理参数和物理过程过于繁多，在工程实验中应用并不方便。鉴于此，De等^[100]在研究含能材料RDX的冲击响应时仍采用唯象型模型。首先在弹性变形部分，考虑了非线性的热弹性效应；在塑性变形部分，采用了率无关的屈服准则作为流变规则，对于滑移障碍阻力则在(36)式的基础上真正考虑了声子拖曳部分的贡献

式中：$f_{\alpha}^{}= f_{\alpha} ^{\rm h}f_{\alpha} ^{\rm s} + f_{\alpha} ^{\rm p}$，$f_{\alpha} ^{\rm h}$、$f_{\alpha} ^{\rm s}$和$f_{\alpha} ^{\rm p}$分别为位错相互作用和生成、热软化和湮灭以及声子拖曳对于可动位错的影响，前两种影响因素的具体表达式为

(39)式和(40)式其实就是(36)式的中间部分。(38)式与(36)式唯一的不同在于额外的影响因素$f_{\alpha} ^{\rm p}$，即高速位错所遇到的声子拖曳。在该模型中，为了吻合实验，采用基于经验的声子拖曳力

式中：${\dot \gamma} _{\rm{p0}}^{\alpha} = \rho _{\rm{m}}^{\alpha} {b^2}\tau _{\rm{yT}}^{\alpha} /B$，$\tau _{\rm{yT}}^{\alpha} $为高温下的屈服应力，B为拖曳系数。将此模型应用到RDX时，能够捕捉到实验所观察到的弹-塑性波特征。图4所示为RDX在[111]方向受到1.25 GPa加载时的粒子速度剖面，可以看到明显分开的弹性波和塑性波、弹性波后的应力弛豫、弹性波衰退以及随着进程逐渐加大的弹-塑性波间距。

大多数CPFE都是在Lagrangean框架下建立的，在处理高应变率等极端条件下的材料响应时需要计算所有格点的变形，而此时一些结构材料的变形经常会严重局部化。鉴于此，Cazacu和Ionescu^[102]构造了基于Eulerian框架的动态晶体塑性模型，在面对高应变率时有很大的优势，如图5所示等径角板牙的变形情况。此模型并非完全采用（6）式的幂律规则，而是引进了黏性系数来描述流变应变的演化情况

式中：${\eta _{\rm{s}}}$即黏性系数。对于剪切阻力和剪切变形的演化，均在Eulerian框架下

当 ${\eta _{\rm{s}}}$ 取特定值时，(42)式会退化成(6)式。将此模型运用到特殊结构的等径角板牙时，能够有效预测其变形情况，如图5所示。研究结果表明：即使在较低应变率下，变形局部化仍主要由黏性性质决定；而在高应变率下，塑性性质和塑性各向异性并不重要，此时处于主导地位的是惯性效应。

Clayton^[103]在描述非线性热弹性时利用对数弹性理论方法，而非Eulerian或Lagrangean理论。相对于Eulerian或Lagrangean框架，此对数框架下的非线性弹性对于蓝宝石、金刚石及石英具有更高的准确度，例如需求更高阶的弹性常数。

基于唯象本构方程的模型在忽略位错物理机制的基础上，提供了简单而有效的方式模拟预测工程中遇到的材料宏观效应，其参数很多是根据经验或者实验数据拟合得到的。然而，正如3.1节所述，极端条件下材料的响应特征主要是由于位错滑移物理机制的转变，所以必须发展基于物理本质的CPFE才能更有效地捕捉材料动态响应中微观结构的变形及其对整体性能的影响。本节仍然从非线性弹性、位错拖曳机制以及参照系框架方面加以叙述。

Winey和Gupta^[104]发展了各向异性的连续框架来描述受冲击单晶的非线性热弹性响应。而后，他们将此非线性热弹性方程应用到晶体塑性模型中，模拟受任意方向冲击的单晶中的弹性波和塑性波^[105]。其方程建立在一维冲击条件下，利用Orowan方程描述基于位错滑移的剪切率，但本构关系却利用了非常强的假设。首先，各个滑移系之间的本构关系相互独立(即(11)式可以完全略去上标$\alpha $)。其次，可动位错密度只取决于当前的剪切应变和切应力

式中：${\rho _{\rm{m0}}}$为初始位错密度，M为累积因子，H为硬化常数。此位错密度累积公式只适用于中小塑性变形$\left( {{\gamma ^{\rm p}} < 30{\text{%}}} \right)$，不过一般的冲击变形都符合此条件。最后，位错速度基于当前切应力与滑移阻力

式中：${v_0}$为滑移面剪切波声速，B为黏性阻尼系数，${\tau _0}$为位错启动的临界力。此塑性模型并未考虑热激活效应，滑移阻力也为常值。图6为LiF受[100]方向冲击时的纵波历史，实线为实验测量值，虚线为模拟值。图中所示的弹性波衰减是率相关、应力弛豫材料的典型响应特征。

以热激活位错滑移为基础的内变量模型能够很好地描述应变率在 ${10^4}\;{{\rm{s}}^{ - 1}}$以下的黏塑性变形。然而，在高压高应变率情况下，黏塑性流变受许多机制的控制，特别是可动位错的动态生成和阻滞效应必须加以考虑。一般认为位错运动的转变机制在 ${10^5}\;{{\rm{s}}^{ - 1}}$左右，低于此应变率时位错启动及穿越障碍是在热扰动的帮助下进行的，而高于此应变率时则由拖曳阻滞和相对论效应控制着位错的生成和连续滑移。对基于位错机制的本构模型做出代表性改进的是Austin和McDowell^[106]，他们将包含位错亚结构特征的本构方程应用到弱冲击(根据材料而不同，但均低于25 GPa)、应变率为10⁴～10⁸ s⁻¹的情况中。在其模型的本构关系中，用(11)式的Orowan方程联系剪切率与位错速度及密度。对于位错速度，采用了Clifton^[107]模型

式中：${\bar L}$为相继障碍之间的平均滑移距离，${t_{\rm{w}}}$为位错处于障碍处等待热激活的时间，${t_{\rm{r}}}$为障碍之间位错的“奔跑”速度。对于热激活跨越障碍有

式中：${\gamma_{\rm{G}}}$为位错等待障碍的频率。类似于(46)式，对于某一滑移系（略去上标$\alpha $），忽略惯性效应，“奔跑”位错速度有

(49)式中的有效应力与(15)式有所不同。黏滞系数B的表达式基于相对论效应

式中：${B_0}$为黏滞系数的本征值，${c_{\rm{s}}}$为剪切波波速。(50)式使位错速度无法超越波速。结合(48)式和(49)式，可以得到

至此，位错速度的完整表达式可以给出

其中位错滑移阻力可以由Taylor关系的(15)式给出。对于位错密度的演化，考虑了位错的形核、累积、湮灭、捕获和恢复。利用此模型预测的稳定塑性波、粒子速度、波前剪切力和塑性速率、Hugoniot材料强度都与实验吻合。图7给出了此模型预测的Hugoniot塑性变形的临界剪切应力，以及与实验数据的对比。

Austin-McDowell模型有效考虑了高应变率下位错拖曳机制的问题，此后被很多研究者引用、改进和完善。Lloyd等^[108]将Austin-McDowell的晶体塑性模型拓展并与高阶热弹性方程结合，获得能够体现晶体学各向异性的热动力学平面波方程。除了弹性部分，他们对塑性本构关系中不可动位错密度的描述也做了改进，即将不可动位错细分为林位错和平行位错，能够更有效地描述各向异性性质。图8所示是沿不同晶向冲击下单晶铝的粒子速度剖面，很明显地看到冲击方向对于材料高应变率响应的影响巨大。

Luscher等^[109]在Austin- McDowell模型和Lloyd模型的基础上，考虑了非局域化的CPFEM模拟材料的冲击响应。非局域化模型主要针对非局域的位错-位错、位错-界面之间的相互作用，其中很重要的一个问题是在一定的初始及边界条件下位错密度的守恒准则。此模型在处理位错密度演化时耦合了连续位错输运方程，单晶铜的平板冲击结果说明了位错输运和位错堆积的重要性。此后该模型被用于模拟RDX的各项异性冲击响应，模拟结果与实验吻合，证明了RDX中位错主导塑性机制的重要性^[110]。

高、低应变率下的位错机制有所不同，基于不同机制的CPFE有其局限性，因此寻求一个能适用于广泛应变率下较统一的模型显得非常重要。Hansen等^[111]发展了一种适用于多种应变率的基于位错的CPFE方程，将位错按其速度分为3类：自由可动位错、应力才能驱动的堆积阻塞的位错、热激活才能驱动的位错碎片。高速位错不再考虑热激活因素，而考虑拖曳效应；低速位错则只考虑热激活效应。模型所用参数不是从实验中拟合得到，而是具有实际物理意义，且可以从原子模拟获取。此模型被用来研究不同应变率下铜的流变应力，并与大量实验和模拟做了对比，如图9所示，从中可以明显看到存在一个临界的转变应变率。

Shahba和Ghosh^[112]基于Austin- McDowell模型，将低应变率的位错热激活机制与高应变率的位错拖曳机制相结合，力求建立一个统一且适用范围广的晶体塑性模型。在其模型中，位错速度被分解开来考虑，对 ${t_{\rm{w}}}$ 和 ${t_{\rm{r}}}$ 进行比较，较小的被忽略掉。另外，位错密度演化中考虑了GNDs密度的贡献，能够处理尺寸效应。此模型被应用于研究Ti-7Al合金，能够捕捉到位错运动的机制转变在10⁴～${10^5}\;{{\rm{s}}^{ - 1}}$。

上述Lloyd等^[108]发展的高应变率热弹性-塑性模型是在Lagrangian框架下获得的。而实际上在高应变率大变形下，复合的Lagrangian-Eulerian或全Eulerian热弹性理论能够更好地描述微结构的演化特征。之后，Lloyd等^[113]采用Eulerian应变构建了热弹性方程，具有更快收敛的特性。此模型的模拟结果体现了初始详细的微观结构对于材料的高应变率变形具有重要影响。

位错是承载塑性变形的主要方式，但却不是唯一的方式，特别是在极端条件下。在高应变率条件下，塑性变形随时间快速变化，位错萌生、运动以及湮灭等活动受位错本身缺陷性质的控制，从而导致其对塑性变形的贡献不足。此时，其他变形机制如相变机制、孪晶机制就会发挥巨大的作用。然而在更极端的情况下，材料以损伤破坏的形式对快速变形作出响应。本节将介绍引入这些额外变形机制的CPFE在材料动态响应中的进展。

如2.4.1节所述，塑性模型已经能够有效引入相变机制，相变部分的变形梯度是基于马氏体相变理论得到的，这些模型已经广泛用于准静态变形模拟。对于高应变率情况下的变形模拟，Barton等^[114]在处理铁冲击响应时考虑了$\alpha \to \varepsilon $相变。在此模型中：塑性部分采用唯象型的本构方程；对于相变部分，其相变驱动力包括机械力、基于热的化学自由能的改变、各组分之间的相互作用3部分。模拟结果显示变形机制包含相变和弹塑性变形，与实验结果相吻合。

最近，Addessio等^[115]发展了高应变率下单晶RDX的塑性模型。此模型包含了非线性热弹性、相变以及塑性滑移，能考虑大变形和各向异性性质。塑性滑移部分采用引入热激活现象的唯象本构方程；相变部分基于马氏体相变理论发展而来，主要考虑了$\alpha $多形体向$\gamma $多形体的转变。当处于相变临界压力之上时，模拟结果能够有效描述实验观察的现象。

不同于FCC和BCC，HCP金属中经常见到孪晶。孪晶机制已被引入晶体塑性模型以研究HCP金属的变形机制。然而，应变率对HCP孪晶变形影响机制的相关研究最近才发展开来。

Knezevic等^[116]利用实验和晶体塑性模型研究应变率及温度对HCP锆中主要和次要孪晶系统的影响，温度为76～673 K，应变率为0.001～4500 s^–1。在塑性变形部分采用了唯象型的幂律定则作为流变规则。通过模拟发现，主要和次要孪晶系统的选择基于加载方向、温度和应变率。但其中的原因并非孪晶变形对应变率敏感，而是由于塑性滑移的应变率敏感性。

由于镁合金的广泛应用，利用引入孪晶机制的CPFE研究镁合金的动态响应更受关注。Ardeljan等^[117]建立了多尺度的晶体塑性模型：在介观尺度将速度梯度分解为塑性部分和孪晶部分，宏观尺度的变形则利用Taylor形式的均匀化处理；塑性部分和孪晶变形部分的流变规则都采用唯象型的幂律定则。他们利用此模型研究了77～423 K、10⁻⁴～3000 s⁻¹情况下AZ31镁合金的微结构演化和力学响应，发现应力-应变曲线、孪晶体积分数和织构演化都与实验现象吻合，且各个变形系统对应变率和温度并不十分敏感。Meredith等^[118]利用率相关的本构方程描述非基面滑移，率无关的本构方程描述孪晶系统和基面滑移，将模型用于研究AMX602镁合金的准静态和动态响应，并与实验作比较，结果也显示孪晶变形对应变率不敏感。García-Grajales等^[119]采用新的应变率相关模型研究了AZ31镁合金的准静态和动态响应，模型中考虑了所有滑移系和孪晶系统的率相关性，准静态条件下的结果与率无关模型吻合，而高应变率下能够准确描述变形系统的率敏感性，收敛性亦被证实。

最近潘昊等^[120]研究了金属铍受冲击时孪晶对其加卸载动力学性能的影响，与上述模型不同的是，此CPFE模型考虑了非线性的热弹性变形，结果发现，非线弹性变形过程是导致材料准弹性卸载的主要原因。图10给出了加卸载过程中的孪晶份额变化，发现卸载过程也能导致孪晶份额的增长。

损伤破坏也是材料动态响应的重要形式。现阶段利用CPFE模型研究材料在中高应变率下的破坏时，主要考虑的是界面破坏和孔洞的演化。下面围绕这两点叙述最近的相关进展。

Clayton^[121]将界面破坏机制引入CPFE来研究钨合金在高应变率下的拉伸变形。该模型中：热弹性和唯象型塑性本构关系能够捕捉到有限变形、应变率效应、热膨胀、热软化、热传导和弹性产热等典型动态现象，内聚区模型描述钨-钨晶界以及钨-基体界面的破坏，而晶粒内部的破坏则被忽略。结果发现，晶粒取向对于破坏机制和等效刚度有显著的影响。Vogler和Clayton^[122]在此模型的基础上结合实验研究了不同冲击速度下钨合金的动态响应。虽然实验与数值模拟很难在同一尺度上实现，但是他们力求在微结构特征上进行两者的对比。通过比较实验和模拟结果，证实了晶间断裂是材料破碎的重要机制，图11即为晶界破坏的模拟结果。

Nguyen等^[91]为单晶动态韧性断裂发展了第一个基于位错的CPFE模型，并将其应用到商业的有限元软件包中。其晶体塑性部分基于Lloyd等^[113]建立的单晶位错黏塑性模型，能够考虑动态变形中运动位错的声子拖曳效应和相对论效应。损伤演化基于孔洞生长，且受微惯性和位错动力学、位错亚结构演化的控制；受益于孔洞动力学的半解析解，近似方程极大地增强了计算的效率。晶体塑性和动态孔洞生长的耦合通过位错结构和孔隙率演化实现。此耦合模型能够描述位错林交错和高速位错对孔洞生长的本征阻碍作用。

相较于宏观动态本构模型，CPFE能够从微结构演化出发研究材料的动态响应。很多材料的动态各向异性就是源于晶体结构之间的晶向关系及其导致的织构特征，而CPFE能够从各个滑移系的滑移剪切与宏观变形的关系出发处理微观至宏观的变形特征，是研究材料动态响应的强有力工具。CPFE具有极大的灵活性，能够有效引入其他变形机制，如孪晶、相变和损伤失效机制，而这些其他机制在动态变形中起到非常重要的作用，特别是在高速冲击等极端情况下。唯象型的CPFE模型虽然缺乏物理机制，但简单易解，在工程应用方面非常有效。物理型的CPFE模型则能够揭示材料应变率效应的物理本质，解释材料动态响应的一般现象，有利于材料的设计和制造。

基于位错剪切滑移的CPFE对于材料的动态响应主要集中在非线性热弹性、位错拖曳机制和非Lagrangean模型框架这3方面。高压高应变率使得弹性变形不再为小量，且热效应明显，这对于材料的变形影响很大。高应变率下位错速度极高，位错滑移由热激活机制控制转变为位错拖曳机制，需要对基于位错的物理型CPFE模型进行改进和完善。对于具有特殊结构的材料，高应变率下变形局域化严重，基于Eulerian框架的CPFE模型比Lagrangean框架下的模型更具优势。另一方面，高应变率下位错承载塑性变形的能力有限，此时其他变形机制被激活，将相变、孪晶和损伤破坏等物理过程耦合到位错滑移的CPFE模型中是必需的。

虽然近年来CPFE模型在研究材料的动态响应方面取得了较大进展，但仍然存在很多难题和挑战，主要可归纳为两方面：本构关系背后物理机制相关研究，与实验或工程结合方面的难题。

传统唯象型模型缺乏对物理机制的足够认识，而且高、低应变率下的物理过程不尽相同，因此必须发展物理型模型。基于位错的模型虽然能够认识到高、低速位错阻力的转变，但目前仍缺乏对高速位错的动力学机制的全面认识，例如晶界处位错的形核^[123–125]、高速位错的交错作用机制及交滑移机制的改变、非惯习面的激活等。而且材料动态变形经常包含其他变形机制，而现有的孪晶等机制均基于现象型描述，缺乏物理基础，对这些孪晶和位错微结构相互作用的研究更加困难。动态损伤失效是极端条件下材料的重要问题，但其失效破坏准则十分复杂，这也是限制CPFE模型在动态响应方面发展的重要因素。解决这些背后物理机制问题的有效方式需要更低尺度的模拟为物理型CPFE模型提供支撑，这就涉及到多尺度模拟问题，以及低尺度向高尺度机制和参数传递的均匀化问题。

与为了探索背后的物理机制而向低尺度模拟结合相反，为了实现实验比较且能够有效服务于工程实践，则需要与宏观尺度模型相耦合。CPFE方法本来就是横跨微结构尺度和宏观尺度的多尺度模拟工具，但也涉及到均匀化的问题。唯象型CPFE模型由于具有模型简单且数值化计算效率高的特点，容易集成到各种商业软件中，因此应该受到同样的重视和发展，然而合理处理唯象型与动态物理型CPFE模型的结合关系也是未来应该关注的难点之一。