位错动力学在极端环境力学中的发展及应用

崔一南，柳占立，胡剑桥，刘凤仙，庄 茁

（1. 清华大学航天航空学院应用力学实验室，北京 100086；2. 中国科学院力学研究所非线性力学国家重点实验室，北京 100190）

随着现代科技的发展，越来越多的材料和结构需要在极端环境中服役。金属材料在高应变率、高温、高压等极端条件下的力学行为与常态下有着显著的区别，变形失效机理非常复杂，给材料在极端条件下服役时的性能评价和安全评估带来了挑战。在极端环境下，相关实验研究往往需要巨大的经费投入，并且仍然显著受限于时空分辨率，难以解耦复杂影响因素。近些年来，多尺度数值计算方法的迅猛发展为理解极端环境下材料与结构的变形失效机理提供了前所未有的机遇。

如图1 所示，材料变形和破坏的内在多尺度物理过程可以借助不同时空尺度的计算方法加以研究。通常来说，计算的尺度越小，捕捉的物理过程越完备，但是也意味着更显著的复杂性，其所能研究的时空尺度也就越小。目前广泛采用的纳观尺度计算方法主要是分子动力学(Molecular dynamics，MD)模拟方法。该方法将原子看作质点，通过求解牛顿运动方程，获得系统内所有原子的运动轨迹，基于对大量离散原子运动性质的统计平均得到宏观意义上的物理量。这一方法能够从最基本的原子尺度刻画晶格迁移、缺陷演化等信息，通过大规模原子模拟，探寻材料变形的原子级机理[1-3]。该方法对于揭示缺陷形核及基本相互作用机制等提供了很多有价值的信息。例如，Zhu 等[4]通过分子动力学模拟系统地研究了纳米单晶材料在不同应变率加载下的变形过程，研究结果表明位错的表面形核机制对纳米材料的屈服起着主导作用，继而从原子热激活角度出发提出了在给定温度和应力状态下的位错形核率；Liu 等[5]揭示了单个位错线与离散的辐照缺陷之间的相互作用机理。但是，由于受计算能力的限制，该方法能够处理的问题在时间和空间尺度上仍然具有很大局限性。经典的MD 模拟目前能够处理的问题尚在皮秒(10−12s)量级。此外，为了缩短计算时间，通常需要增加载荷速率，模拟过程中的应变率远高于实验值，而高应变下的微结构演化特征可能失真。

图1 材料变形失效行为研究在不同时空尺度下通常采用的计算方法[6]Fig. 1 Simulation methods at different temporal and spatial scales to investigate material deformation and failure mechanisms[6]

为了拓展研究问题的时空尺度，大量采用连续介质力学有限元计算研究材料的力学响应，针对离散的微观结构进行统计平均，得到宏观意义上的物理量，通过在连续介质框架中引入与微结构演化相关的内禀参数来描述微结构演化。对材料微塑性行为研究通常基于两种理论框架。一种是基于J2 塑性理论，研究者通过在屈服应力中引入考虑几何必需位错的Taylor 硬化模型，建立了应变梯度塑性理论[7-8]。该方法采用各向同性硬化，不包含晶体中滑移系的信息，因而无法反映晶体取向的影响。另一种是晶体塑性理论，研究者以位错密度为演化量建立每个滑移系的位错密度演化方程或者直接建立各个滑移系上塑性量的演化方程来描述材料的塑性变形行为[9-12]。该方法较好地预测了单轴压缩、纳米压痕、变形形态学、动态冲击等问题，能模拟的时间和空间尺度相比MD 方法有较大的提升，然而模拟过程通常是引入以位错密度为代表的内变量，对位错微结构等细节的刻画具有一定的局限性。

此外，MD 模拟难以为晶体塑性理论所考虑的大量微结构累积行为描述提供直接信息。为解决这一问题，必须借助微介观尺度的离散位错动力学(Discrete dislocation dynamics，DDD)方法向下联系MD，向上联系晶体塑性理论。DDD 方法以位错线为研究对象，可以从离散层次处理比原子尺度更大的塑性问题，其特征时间和空间尺度在亚微米、微米和秒量级，与原位投射电镜（TEM）实验中位错运动的时空尺度相一致。因此，离散位错动力学被认为是一种“计算实验”[13]，是揭示微尺度塑性形变失效机制非常有效的手段[14-16]。

本文中首先简要介绍位错动力学方法，然后分别介绍近些年该方法在高应变率、高温及辐照极端环境下的发展和应用。

1 位错动力学方法简介

DDD 方法主要是基于位错的弹性理论发展起来的，通过求解弹性体中位错线的运动方程来计算塑性应变。如图2[17]所示，其基本思路为首先将位错线离散为具有一定长度的位错段，然后计算位错段所受到的力，继而计算在一定的温度等条件下的位错运动速度，据此更新位错的位置。当位错之间的距离很近时，根据能量法则判断其是否发生相互反应。此外，对于螺性位错还需要判断其是否发生交滑移机制，继而更新外载条件，并且计算位错产生的塑性应变。此过程不断循环，直到加载完成。

早在1968 年，Foreman[18]就基于二维位错动力学研究了位错和刚性障碍物的相互作用。20 世纪80 年代以来，随着计算条件的不断改善，大量的二维离散位错动力学模拟(2D-DDD)相继涌现。它们通常只考虑无限长的直刃位错组，从垂直位错线的平面上看，每根位错线都对应一个带有滑移系特征的离散的点。Deshpande 等[19-20]发展了有限变形框架下的2D-DDD，模拟了微柱拉伸试验，结果表明当限制拉伸轴的旋转时，会引起几何必需位错密度的累积，弱化尺寸效应而强化包辛格效应。我国华中科技大学李振环课题组利用2D-DDD 模拟，系统地研究了微柱压缩[21]和纳米压痕[22]中的尺寸效应。此外，他们在2D-DDD 中引入了位错攀移和晶界穿透机制，研究了多晶镍基合金的疲劳裂尖特征等[23]。但是，通常2D-DDD 无法考虑三维位错连接的形成和破坏，所捕捉到的应变硬化行为多是由GND 的存储引起的[24]。为解决这一问题，Benzerga[25]发展了2.5D-DDD，引入了一些表征三维位错反应的法则，捕捉到了微柱压缩过程中的位错匮乏硬化现象。但是，此类模型涉及到诸多简化，难以考虑复杂的位错反应，也难以重现真实的位错构型。

图2 位错动力学基本流程图[17]Fig. 2 Flowchart of dislocation dynamics method[17]

三维离散位错动力学(3D-DDD)则可以更全面的考虑不同类型的位错反应，并且自然捕捉真实的三维位错形貌。如图3 所示，在现有的3D-DDD 中，位错线有两种离散方式。一种是依赖于潜在晶格的描述，将位错线离散为刃型-螺型位错段或刃型-混合型-螺型位错段。力的计算在位错段的中点上进行，即图3(a)中的空心点。它起源于Kubin 等[26]的模型，便于保证位错段的运动始终沿着晶体学取向进行，计算效率高。另一种是连续化描述，将位错段离散为通过节点连接的直线位错段[27]或曲线位错段[28]，位错段的节点位置和取向可以连续变化，力的计算在位错段的节点上进行，即图3(b)中的实心点。该方法易于处理复杂形状的位错线[29]。

图3 三维位错线的不同离散方式：（a）基于晶格的描述，（b）连续化描述Fig. 3 Different discretization of three-dimensional dislocation lines:(a) lattice based description，(b) continuum description

单纯的3D-DDD 模拟往往存在一定的局限性：(1)通常采用周期性边界条件，位错之间的相互作用基于位错在无限大体积中的弹性应力场得到，难以有效计算复杂形状表面或界面的镜像力；(2) 3DDDD 模型多采用固壁假设，适用于模拟小变形情况，无法考虑变形过程中试样形貌的变化；(3)位错核附近的强烈非线性效应难以精确描述；(4)计算成本随着位错段数量的增加而急剧增加，模拟所能达到的尺寸范围和应变范围仍然比较有限。为了解决上述问题，前人开展了大量工作将DDD 与其他尺度的计算方法(如分子动力学、连续介质力学)结合，实现多尺度模拟。到目前为止，DDD 与有限元方法[12,29-30]或边界元方法[28,31]的耦合主要有两种思路：一种是叠加法(Superposition method，SPM)，另一种是离散连续晶体塑性计算方法(Discrete-continuous model，DCM)。

叠加法最初是由Giessen 等[30]提出的。如图4 所示，有限晶体试样中的总应力场σ包括两个部分

式中：σ∞为位错在无限介质中产生的应力场，通常可借助解析解得到；为附加应力场，用于抵消该面力 的作用并修正真实的位移边界条件。σ∞会在真实试样的表面产生面力

式中：ns为表面的法向。

叠加法可以很好地捕捉位错间的短程相互作用[32]，但是每一根位错的解析应力场都需要在每个计算步重新计算，占用较长的计算时间。同时，该方法基于两个弹性场的叠加，难以直接体现塑性应变的概念。

图4 叠加法示意图（a）和DCM 方法变量传递示意图（b）[33]Fig. 4 Schematic of superposition method (a) and the variable transfer process in discrete continuous method (b)[33]

DCM 基于细观力学中的特征应变思想，直接根据DDD 计算的塑性应变，在统一的有限元框架下求解边值问题的解[33-35]。如图4(b)所示，主要包括4 个变量传递过程[33]。(1)利用DDD 计算位错的滑移量，进而计算塑性应变εp，并将其传递给有限元模型(FEM)，以此来代替传统唯象的本构关系

式中：M是质量矩阵，u为位移场，N是形函数向量，ρ是材料密度，Ω代表体积，Be=grad[N]，fint是内部节点力向量，fext是外加面力T0和初始应力场σ0引起的外部节点力。这里，σ0是由试样中初始存在的静态位错引起的。随后，将FEM 计算的应力场σ传递给DDD，并作为外力驱动位错运动。由于FEM 无法准确计算位错核附近的应力场，位错间的近程相互作用力需要通过解析解加以修正。(3)FEM 计算得到的位移场u传递给DDD，以保证位错动力学模型和有限元模型同步变形。(4)FEM 计算得到的变形梯度场传递给DDD，以考虑滑移系的转动效应。

如何将DDD 中计算得到的离散塑性应变场εp合理分配至连续介质力学框架中的材料点上，是DCM 的核心问题。研究者们对此开展了大量研究，并提出了多种不同的塑性局部化分配方案[12,32,35]，其主要差别在于塑性局部化时所选的特征体积以及向有限元积分点分配塑性应变量的权函数不同。Cui 等[33]通过定量对比各种局部化分配方案的优缺点，整合各方法的优势，提出了一种新的塑性局部化方案，该方法考虑了局部化参数对滑移系信息的依赖性，表现出良好的计算精度。同时，还提出了有效的应力插值方案，弥补了DCM 相对于SPM 方法在计算镜像力方面精度不足的缺陷。

2 位错动力学在高应变率效应研究中的进展与应用

与常态下相比，高速冲击导致的极端高温、高压、高应变率条件使得金属材料的物理和力学性质均会发生显著的改变，这给获得材料在冲击波加载下的状态方程、冲击强度等动力学参数带来了挑战。在冲击波加载下，材料的变形从受力处向材料的内部快速传播，冲击波的波形则能够反映材料的动态响应特征，包括材料物态和强度等的变化。目前材料动态本构模型的建立及参数确定主要依赖于宏观屈服强度的拟合，相关模型的适用性受到极大的限制。因此，大量研究者试图深入研究极端条件下材料的塑性变形机理，构建基于微介观机理的动态本构模型。

随着计算机数值模拟技术的发展，各种数值模拟方法已经成为获取材料在冲击载荷等极端条件下内部微结构演化的重要手段。目前以MD 模拟和DDD 为核心的多尺度计算方法在研究金属材料冲击问题中得到了广泛的应用。本节将主要介绍DDD 方法在高应变率条件下的发展及应用。

DDD 方法在研究高应变率乃至冲击加载下材料的动力学行为时面临新的挑战，主要包括：(1)高应变率下位错的运动规律发生变化；(2)高应变率下位错的形核萌生过程起着重要作用；(3)高应变率往往伴随着高温和高压效应，其中冲击波压缩会引起高温。此外，位错运动导致大量的塑性热产生，使得材料温升进而引起热力耦合过程；(4)超高应变率下，位错动力学的准静态应力场解不再适用，亟需发展高速位错的弹性动力学理论，揭示其相互作用规律由于时间依赖性及历史依赖性带来的新机理。对上述问题的研究与认识常常需要借助大量MD 的研究、冲击实验以及深入的理论分析，下面将分别展开进行具体介绍。

2.1 高应变率下的位错运动规律

传统的DDD 模拟常常研究材料在准静态加载下的变形过程，位错的运动通常被认为处于过阻尼状态，即忽略位错的惯性效应[37]。Hirth 等[38]从热力学角度出发推导了位错在不同速度下的有效质量；Gurrutxaga-Lerma 等[39]进一步总结了各种高速位错的运动方程，并详细探讨了不同运动方程对材料动力学响应的影响。Wang 等[40]通过DDD 模拟对比了考虑包含位错惯性效应和不包含位错惯性效应的运动方程对材料力学行为的影响。研究结果表明，当加载应变率超过103s−1时，位错惯性效应对材料力学行为的影响开始凸显，从而不可被忽略；同时Wang 等也详细探讨了质量项的相对论效应(即当位错速度接近剪切波速时的质量修正)，发现即便当应变率高达106s−1时，大多数位错的运动速度仍小于50%的声速，因此可以在一定程度上忽略速度变化对位错等效质量造成的影响。因此，在考虑高应变率效应时，常采用包含位错惯性效应的运动方程

式中：me为单位长度位错段的等效质量；Bz为黏性阻尼系数，f为位错段的受力。

2.2 高应变率下的位错形核

随着加载应变率的增加，主导变形过程的位错机制会发生变化，合理引入初始的位错缺陷(包括均匀形核过程和非均匀形核过程)对模拟动态变形行为的准确性十分重要。Shehadeh 等[41]通过基于DDD 的多尺度计算模型研究了冲击载荷下位错的均匀形核和增殖过程，结果表明：位错的均匀形核是冲击载荷下主导材料动力学行为的位错机制，而且随着变形的持续进行，材料最终会达到一种静水压的载荷状态。Gurrutxaga-Lerma 等[42]探讨了不同应变率加载下(10～1010s−1)Frank-Read 位错源的开动过程，结果表明当加载应变率高于5 × 107s−1时，位错的均匀形核过程将会先于Frank-Read 源开动；Gurrutxaga-Lerma 等[43]通过对两种形核过程进行进一步统计，发现均匀形核对应的位错数目高于Frank-Read 位错源两个数量级，并详细探讨了几种典型金属材料(铝、铁、镍与钼)在冲击载荷下的动态响应，建立了位错微结构演化规律与动态力学响应之间的关联。McDowell 等[44-45]通过MD 方法系统地研究了单晶铜中的位错均匀形核过程，并结合参数化的位错密度演化模型[46-48]对动态冲击问题展开了大量研究。尽管诸多研究者都认识到在冲击载荷下位错形核的重要作用，但是在模拟过程中，常常都是通过人为预设位错源数目的方式引入位错均匀形核过程的，这一过程仍有待更加深入的探讨。

通过采用考虑惯性效应的位错运动方程以及在模型中引入基于MD 模拟得到的位错表面形核判据[4]，Hu 等[49]通过位错动力学方法系统研究了单晶微柱在不同应变率下的力学响应以及对应的微结构演化规律，并建立了预测单晶微柱中流动应力的解析模型，且结合该模型系统地分析了单晶微柱中位错源数目的影响，给出了亚微米单晶柱变形的超曲面相图。

2.3 高应变率伴随的高温效应

冲击波压缩引起的温度升高可以从热力学定律出发加以推导[50]，但是在现有的位错动力学方法中还鲜有考虑。Shehadeh 等[41]和Hu 等[51]考虑了动态冲击过程中剧烈塑性变形引起的热效应。在冲击载荷作用下，波阵面上的材料先经历弹性变形后经历塑性变形过程。在塑性变形过程中，位错会在极短时间内剧烈运动产生大量塑性功，这些塑性功会转化成材料的内能，从而导致局部温度升高。因此，变形过程中的能量转化

式中：ρ为材料密度，Tw为温度，c和λt分别为材料比热容和热传导系数，γ代表塑性功转化为内能的比例因子，σ为应力张量，为由位错运动产生的塑性应变率。

通过计算材料域内的温度场可得材料的热弹性本构关系

式中：E、ν和α分别为材料的杨氏模量、泊松比和线膨胀系数。最后，在统一的连续介质框架下，利用FEM 求解动量方程

式中：M、C、K分别为质量矩阵、阻尼矩阵和刚度矩阵，U为节点位移，F为力矢量。

Hu 等[51]在图4(b)的计算框架基础上，利用FEM 求解了温度场的演化，传递给DDD 实现了对动态冲击过程的有效计算。此外，通过系统的MD 模拟得到了位错均匀形核的相关参数，并通过粗粒化方法将这一位错形核过程引入到耦合模型的计算中。通过该耦合模型研究了单晶铜柱在不同冲击载荷下的动态力学行为，如图5 所示，模拟结果能够有效捕捉冲击载荷下一些典型的位错微结构，如亚稳态位错环和位错增殖形成的剪切带等。

图5 冲击载荷下的典型位错微结构：（a）DDD-FEM 模拟结果[51]；（b）实验中观测到的剪切带[52]；（c）位错线缺陷与位错环结构[53]Fig. 5 Typical dislocation microstructure under shock loading: （a） DDD-FEM simulation results[51]（a）；（b）shear bands observed in experiments[52]；（c）dislocation lines and loops microstructure[53]

在冲击波加载过程中，冲击波波阵面上的位错形核与增殖会使应力水平有所衰减。Gurrutxaga等[54]通过基于位错动力学的多尺度耦合模型研究了不同应变率加载下冲击波波阵面上的应力衰减与位错密度变化过程，如图6 所示，研究表明弹性前驱的衰减与冲击波前沿的位错形核密切相关。模拟结果定量地讨论这一相关性并与实验结果对比良好。值得一提的是，在位错动力学模拟中，考虑到位错的高速运动，Gurrutxaga-Lerma 等考虑了位错的动态应力场。随后，Gurrutxaga-Lerma 等[55]通过该方法进一步研究了在两种典型的面心立方(FCC)与体心立方(BCC)材料中，温度对弹性前驱衰减的影响。需要指出，上述面向高应变率问题的研究过程中，往往只考虑了温度效应对于位错滑移过程的影响，在第3 节中将进一步展示高温下位错也将借助攀移产生塑性变形。

2.4 高应变率伴随的高压效应

金属材料在高压下也表现出不同的力学行为，但是目前的研究仍然比较有限。Yang 等[56]从原子尺度出发，建立了研究极端高压条件下单晶塑性变形的多尺度模型，并通过ParaDis(一套开源DDD 程序)研究了Ta、Mo 等金属材料在不同压强、温度下的流动应力，其中单晶Mo 材料的力学响应如图7所示。一方面，单晶Mo 的屈服强度随着温度的升高而降低；另一方面，静水压可提升流动应力和模量。Hu 等[57]通过位错动力学研究了静水压加载下位错密度的演化和位错结构特点，结果表明：加载率对位错密度的演化及最终稳定位错密度存在一定的影响，且静水压下位错的形态与单轴压下存在明显的区别。

图6 （a）冲击方向应力衰减，（b）不同加载率下加载时间为75 ps 时的位错微结构演化，（c）位错密度随时间的演化，（d）应变率为1010 s−1 时不可动位错和湮灭位错[54]Fig. 6 (a) stress attenuation in the direction of shock loading, (b) evolution of dislocation microstructure at 75 ps under different loading rates , (c) evolution of dislocation density, (d) immobile and annihilated dislocations under strain rate of 1010 s−1[54]

图7 单晶Mo 在不同温度、压强下的应力-应变曲线[56]Fig. 7 Stress-strain curves of single crystal Mo under different temperatures and pressures[56]

2.5 超高应变率下位错的弹性动力学效应

在极高应变率加载下，位错速度接近声速(亚声速位错)，与位错自身应力场的传播速度相当。因此，基于位错的准静态弹性解来描述位错的长程和短程相互作用已经不再适用[58-61]。单纯引入惯性项已经无法描述超高应变率下的亚声速位错行为。就像接近光速的粒子会呈现相对论效应，亚声速位错也会产生类似的相对论效应，与低速位错的物理机制截然不同。此时，必须考虑位错的弹性动力学应力场。最近的研究表明：弹性波和位错的弹性动力场之间的相互作用可以解释弹性波波前的衰变[54]，入射平面波激活的位错解的叠加可用于确定波振幅的衰减[62]。

针对位错的弹性动力学场的研究早在1949 年就已经展开，但是增加时间维度带来了显著的理论和计算复杂性，因此现有研究多局限于二维问题[56,58]。Xiong 等[60]通过粗粒化二维的MD 模拟结果，建立了原子-连续介质方法来研究弹性动力学效应的影响。Gurrutxaga-Lerma 等[61]建立了考虑弹性动力学框架的二维位错动力学方法，揭示了一对纯刃型位错的相互作用机理。然而，二维计算方法无法揭示真实的三维位错结构的演化规律，及其长程和短程相互作用机理。

时至今日，真实的弹性动力学场下的三维缺陷相互作用机制这一基本科学问题尚无答案。为解决这一问题，Cui 等[63]基于迟滞位错思想来处理耦合的时空积分，推导了任意三维位错环在弹性动力学场下的位移、应变、应力等的理论解，继而首次建立了三维位错弹性动力学方法。如图8 所示，该方法有效捕捉了三维位错环的弹性动力学应力场分布的历史依赖性。在弹性动力学场下，经典的位错理论已经不再适用。比如，两个异号的位错环并不一定永远互相吸引，其相互作用力的大小和符号都呈现显著的历史依赖性[63]。三维位错弹性动力学理论的进一步发展有望带来颠覆传统认知的新的位错理论和动态变形理论。

图8 亚音速位错环应力场σ11 的等值面（V 为位错环的扩张速度，cT 为横向声速[63]）Fig. 8 Isosurface of stress component σ11 of subsonic dislocation(V is dislocation velocity, cT is transverse wave speed[63])

3 位错动力学在高温位错攀移问题上的研究进展

研究金属材料的高温力学行为对材料的设计与加工、航天航空装备、军事工程设计以及安全防护等诸多领域均具有重要意义。强动载荷下，由于存在大量的塑性热，材料局域也会进入高温高压状态。高温热激活作用使得金属的物理和力学性能相比常温情况下发生显著改变，这给获得材料在高温下的微结构演化规律带来了极大挑战。研究早期，高温下晶体塑性形变理论模型的建立主要依赖于唯象的方法[64]，相关模型并不能体现热激活条件下微结构演化速率的温度和时间依赖性。对此，Amodeo 等[65]在唯象假设基础上建立了热激活的攀移运动与外加应力和温度的函数关系，研究了位错运动的动态响应。

DDD 方法为理解高温变形机理提供了新的机遇。在热激活作用下，由于晶格中点缺陷的形成能和迁移能大大降低，点缺陷浓度和扩散速率显著提高。位错线可以通过吸收或释放这些点缺陷来实现垂直于晶体滑移面方向的非保守运动——位错攀移[66]，如图9 所示。这种垂直滑移面方向的攀移运动对高温[67-68]甚至常温[69-71]下的晶体塑性行为有非常重要的影响。在过去的几十年间，位错的攀移过程引起了研究者们的广泛关注[66,72]。

作为一种热激活运动，位错攀移明显区别于位错滑移。首先，位错攀移发生在垂直于滑移面方向，相比于滑移需要更高的激活能；其次，由于攀移过程伴随着物质的迁移，是点缺陷与线缺陷耦合作用的结果，因此攀移过程受点缺陷扩散速率主导[73]；此外，由于扩散过程通常比较缓慢，并且受限于温度、点缺陷浓度等因素，而位错滑移速率接近于材料的剪切波速，导致攀移与滑移在运动速率上存在高达105～1010的巨大差异[74]。因此，位错攀移的出现，使得高温热激活下晶体材料的力学行为具有明显的多物理场、多时间尺度特性：在空间上，涉及点缺陷扩散场与应力场的相互作用；在时间上，同时跨越多个时间尺度。这种由攀移带来的多物理场、多时间尺度特性也给热激活下晶体塑性形变规律的研究带来了极大的挑战。在DDD 中引入位错攀移机制主要需要解决两个问题：(1)如何解决位错滑移、攀移两种运动形式在时间尺度上数量级差异的难题，在统一框架下进行计算；(2)如何建立准确的基于物理机制的攀移模型，将点缺陷扩散与位错攀移运动耦合起来[74-76]。

图9 多余半原子面发生攀移运动的示意图[78]Fig. 9 Schematic showing the climb of edge dislocation[78]

对于位错攀移和滑移速率在时间尺度上的巨大差异问题，研究者主要提出以下几种解决方案：在最早的唯象法中，研究者们将攀移的黏滞力系数作为参数，讨论攀移速率大小对塑性行为的影响，但是这种方法并不物理，不能真正描述攀移的影响[77]。此后，采用顺序耦合的方式来计算，即攀移一个时间增量步然后滑移一个时间增量步，其中攀移时间增量步较大，滑移较小[73]。或者为了方便计算，直接将攀移速率放大至与滑移速率相当[79]，然后二者叠加来进行计算。这两种方法都一定程度上放大了攀移的作用，难以准确描述攀移的实际影响。此外，有学者提出自适应时间增量步的方法[74,76,80-81]：让位错线先进行滑移，此时采用较小的时间增量；当滑移受到塞积，位错开始攀移，此时采用较大的时间增量；当位错攀移到新的滑移面，再开始滑移，如此循环。该方法对于处理局部区域的塞积等问题，计算效率较高，但是在较大规模计算中，位错线并不一定同时发生塞积，需要较长的等待时间，因此一定程度上弱化了攀移作用。

对于在DDD 中引入位错攀移机制的问题，Mordehai 等[79,82]首次将空位扩散机制引入DDD 中，将位错线看作空位的源或汇，通过求解无限长直位错线周围空位的柱状扩散方程获得位错攀移速率，然后基于DDD 建立了离散位错攀移模型。Bakó等[75]通过计算攀移对模型中平均空位浓度的影响，对上述模型进行了修正。Gao 等[83-84]又进一步将这种方法扩展到空位管扩散主导的位错攀移过程，并在后续工作中建立了可同时考虑管扩散和体扩散的攀移模型。Clouet[85]通过对比上述攀移模型与原子模拟得到的攀移速率，证明在一定位错密度下，上述攀移模型在一定范围内能够准确高效地模拟位错的攀移运动。

近年来，该方法被广泛应用于定性解释一系列攀移相关的力学行为。Huang 等[86-87]利用此方法研究了高温下位错攀移对蠕变速率以及多晶体晶粒尺寸效应的影响；Davoudi 等[77]则利用此方法研究了高温下位错攀移对Au 薄膜中流动应力的影响。然而，上述模型均利用了位错线周围点缺陷柱状扩散的解析解，而在得到这一解析解的过程，通常需要假设：(1)单根直位错线，(2)位错线无限长，(3)位错线周围点缺陷满足稳态柱状扩散。因此，此类模型对于任意形状位错线并不适用，而且模型中仅考虑扩散的短程作用，对于位错线周围空位扩散对邻近位错线运动的影响并未考虑，是一种仅适用于位错密度较低情况的“局部”模型。对此，Gu 等[81]考虑了扩散的长程作用，建立了基于格林函数公式的位错攀移的全局模型，适用于任意密度、任意形状的位错线，并证明了局部模型在位错密度较高情况下的局限性。Keralavarma 等[74]建立了2D-DDD 与FEM 结合的模型，利用FEM 求解空位扩散方程，2D-DDD 求解位错微结构的演化，同时考虑位错攀移与滑移，是在以上计算模型上的一大突破，但是在求解扩散方程时采用空位源均匀分布假设，不足以反映位错运动对空位扩散的影响，而且二维模型对于真实情况的模拟有很大局限性。此外，上述基于位错动力学的离散位错塑性描述需要精确到每一条位错线的演化，计算量巨大且难以进行定量表征。对此，Geers 等[73]扩展了传统的晶体塑性理论，在位错运动方程中加入垂直滑移面方向的攀移，在本构方程中加入攀移引起的塑性应变率，同时考虑位错攀移和滑移对塑性形变的贡献。此外，采用局部-全局模型来描述位错运动后晶体中空位浓度的重分配，一定程度上阐述了亚微米尺度晶体塑性变形理论，但是该模型仍仅适用于二维情况，对于三维实际情况的模拟具有很大局限性。Ghoniem 等[88]通过开发三维参数化位错动力学(PDD)，在Gibbs 自由能中引入能够反映点缺陷吸收与扩散的化学势，来处理热激活下的位错构型演化。

图10 离散-连续位错攀移模型中参数传递示意图[89]Fig. 10 Schematically showing the variable transferring in discrete-continuous climb model[89]

Liu 等[78,89]建立了攀移-滑移全耦合的三维离散-连续位错理论模型和数值方法。在图4(b)所示模型的基础上，增加了两个变量传递过程，如图10所示：(1)利用DDD 求解位错段的受力，进而求解其对应的平衡空位浓度。然后将位错段周围平衡空位浓度局部化分配至FEM 中单元节点上，作为空位扩散的初始值与边界条件；(2)根据相应内部与外部边界条件，在连续介质力学框架中，利用FEM 求解空位扩散方程，得到位错核心周围的空位体积通量Ji，并将其传递至DDD 中位错段上，用以求解位错攀移速率，进而在DDD 中模拟位错攀移过程。Liu 等利用该方法诠释了Au 单晶柱应变软化和高温退火硬化的反常规力学机制，揭示了位错攀移对高温退火硬化的主导作用。

4 位错动力学在强辐照环境下的进展与应用

强辐照环境在核电站、核聚变能、空间推进技术等先进能源与动力装置中广泛存在。当材料或结构受到高能中子或离子等的强辐照作用后，其内部产生大量的辐照缺陷。如图11 所示，在面心立方(FCC)晶体中，典型的辐照缺陷是层错四边形(SFT)和间隙环；在体心立方(BCC)晶体中，主要是间隙环和纳米孔洞。这些辐照缺陷的存在带来了服役构件的辐照损伤、脆化、蠕变、肿胀、表面失稳等失效模式，极大限制了装备的服役寿命。MD 只能研究非常少的辐照缺陷和材料中其他微结构的相互作用，而连续介质力学往往只能在已有认知的基础上建立连续化描述，难以揭示新的微观机理，位错动力学方法的发展对理解强辐照环境下累积微结构相互作用主导的材料变形失效行为提供了有效的手段[17,90]。

图11 位错与辐照缺陷相互作用示例：（a） 在受辐照的FCC 晶体中位错被层错四边形（SFT）截获，（b） 在受辐照的BCC 晶体中位错与间隙位错环反应生成位错锁（不同颜色代表不同的柏氏矢量）[16, 91]Fig. 11 Interaction between dislocation and irradiation defects： (a) trappiest of dislocation by SFT in irradiated FCC crystal ,(b) formation of dislocation lock due to the reaction between dislocation and interstitial loop in irradiated BCC crystal (Different colors represent different burgers vector)[16, 91]

在位错动力学框架下考虑辐照缺陷主要有3 种思路[17]。(1)直接对辐照缺陷进行建模，考虑离散的辐照缺陷与位错线之间的相互作用。这是DDD 中引入辐照缺陷效应的最直接也最常用的处理方法，为理解位错与辐照缺陷的相互作用提供了很多有价值的信息[92-94]。研究发现，当一个SFT 被滑动位错扫过时，SFT 可能被剪切，被部分吸收形成位错锁或完全回复[91,95]。当一个间隙位错环被滑动位错扫过时，如果它的柏氏矢量与扫过位错的柏氏矢量共面，这个间隙位错环会被拖拽[91]。否则，这个间隙位错环可能与位错线发生反应，继而形成位错锁，此时间隙位错环可能被位错线完全吸收，或留下一个具有不同柏氏矢量的位错环[96]。该方法已经可以捕捉近原子级的位错-辐照缺陷相互作用机制[91]。(2)通过考虑单个辐照缺陷的应力场，考虑其与位错的相互作用。比如基于无限小环近似，Ghoniem 等给出了棱柱形位错环的应力场公式，基于此可推得间隙位错环产生的力和力矩[17]。该方法相对第一种方法的计算效率更高，但是捕捉微观机理的能力也有所下降。(3)前两种方法在计算高剂量辐照材料时都需要非常大的计算成本，因为此时辐照缺陷的数量密度达到1021～1022m−3。为解决这一问题，Cui 等[97]提出了用连续化辐照缺陷场来描述辐照缺陷效应，通过FEM 求解辐照缺陷场的演化(式(13))，然后传递给DDD，考虑其对位错运动和变形行为的影响(式(12))。为此，他们发展了两个模型来耦合离散位错和连续辐照缺陷场：第1 个模型是辐照硬化模型

式中：τ是辐照缺陷引起的滑移阻力；µ是剪切模量；b是柏氏矢量；N是辐照缺陷数量密度；a是辐照缺陷尺寸；k是无量纲系数，依赖于缺陷强度和分布，借助大量离散位错和离散辐照缺陷计算结果的统计分析得到[97-99]。第2 个是辐照缺陷密度和尺寸演化的动力学模型。该模型可以考虑位错与辐照缺陷的相互作用。研究表明，当辐照缺陷与滑动位错的距离小于y/2 时，可能会被吸收或拖拽[100-102]。因此，滑动位错引起的辐照缺陷破坏率可以表示为[97]

式中：ΔN是辐照缺陷数量变化值，ΔNs是辐照缺陷产生率，ΔA是位错滑移面积增量，V是特征体积，λ是辐照缺陷破坏系数。λy的值可以通过拟合应力-应变曲线得到[103]，也可以通过DDD 与离散辐照缺陷耦合的计算得到[97]。

这一杂交化方法为研究高剂量辐照材料的塑性失稳问题提供了新的机遇。比如，辐照脆化问题中非常经典的“位错通道”形成(如图12 所示，有时也称为“无辐照缺陷通道”)问题得以系统化的研究。基于此，Cui 等[104]提出了定量表征任意复杂三维变形的局部化程度及其空间关联性的方法，建立了预测辐照材料中典型的位错通道形成及其宽化的理论模型，有效预测了辐照条件及关键微结构参数影响下的塑性失稳相图[105]。相关研究对于设计抗辐照材料及核安全应用具有指导意义。

图12 位错动力学与辐照缺陷场耦合计算高剂量辐照材料中的位错通道形成[97]Fig. 12 Formation of dislocation channel in high-dose irradiated material calculated by DDD coupled with irradiation defect fields[97]

5 结 论

极端环境下材料与结构的变形破坏问题是一个涉及多时空尺度耦合、多物理场耦合的复杂问题。迄今为止，研究者在纳观尺度和宏观尺度开展的计算研究较多，可以基于较为成熟的商业软件开展研究。然而，在微介观尺度，尚无成熟的商业软件，研究成熟度也最低。位错动力学方法作为联系纳观尺度和宏观尺度的桥梁，是揭示位错主导的材料塑性变形及失效行为的强有力手段。通过与连续介质有限元方法的耦合，位错动力学方法的应用领域不断拓展，从单纯的力学行为研究，发展到多物理场耦合行为的研究。

本文分别针对高应变率、高温、强辐照几种极端环境，系统介绍了位错动力学方法的发展思路、研究进展及其对理解微观机理、发展理论模型上的贡献。极端环境下材料与结构的变形失效问题仍有诸多尚未解决的难题，迫切需要进一步发展DDD 及其与其他方法的耦合算法，为理解极端环境下的变形破坏机理建立连续化理论预测模型提供强有力的工具。