梯度结构是由一种形态逐渐向另一种形态渐变的结构。在很多生物结构中都能够观察到梯度结构的存在，比如骨骼、竹子、贝壳等^[1]。生物材料中的梯度结构是经过自然界长时间选择和进化的结果，具有一些特定的优异性能，如轻质、高强。受生物梯度结构材料的启发，梯度纳米结构金属材料的设计应运而生。金属材料的梯度结构主要包括：晶粒尺寸梯度、位错密度梯度、孪晶层片厚度梯度、固溶浓度梯度、相梯度以及混合梯度等结构^[2,3]。近年来，对梯度晶粒结构材料的制备、力学性能、疲劳性能、摩擦磨损行为、变形机理和固体本构理论等的研究比较丰富^[4,5,6,7]，梯度纳米孪晶金属中的额外强化与加工硬化研究也取得了重要突破^[8,9]。但是，关于位错密度梯度结构的研究却鲜有报道。

Sakai等^[10,11]在纯Ni热加工过程的动态再结晶(dynamically recrystallized，DRX)晶粒中发现了位错密度梯度。2017年，He等^[12]报道了纯Ni在扭转变形时存在一定的位错密度梯度，并且发现位错密度梯度结构对材料的力学性能有影响。根据描述材料流动应力(τ_s)和位错密度(ρ)间关系的Taylor强化机制(τ_s$\propto \sqrt[]{\rho }$)^[6]可知，位错密度越大，材料发生塑性流动的应力越大。但是Taylor关系只适用于揭示均匀位错密度材料的微结构与宏观力学性能的关联，缺乏对位错密度梯度的描述。

离散位错动力学(discrete dislocation dynamics，DDD)模拟在研究位错的动态演化和位错间的相互作用方面极具优势^{[13,14,15,16,17,18]}。近20年来，大量学者基于离散位错动力学框架研究了位错微结构演化与小尺寸晶体材料力学行为之间的联系，揭示材料的塑性变形机制。Weinberger和Cai^[19]研究了表面位错形核对微柱压缩尺寸效应的影响。Senger等^[20]模拟了微柱压缩的塑性响应与试样的高宽比和初始位错微结构的关系。El-awady^[21]基于离散位错动力学模拟结果，揭示了不同位错密度会导致不同的硬化机制。Li等^[22]和Huang等^[23]系统地研究了位错与内界面的交互作用对材料力学行为的影响，以及金属多层膜的高温力学行为^[24]和镍基高温合金的力学性能^[25]。Liu等^[26]和Cui等^[27]通过离散位错动力学模拟揭示了微柱压缩应变陡增现象和尺寸效应的机理。Fan等^[28]研究了Mg中孪晶的尺寸和取向对材料屈服强度和硬化率的影响。郭祥如等^[29]基于三维离散位错动力学研究了Ni单晶微柱在不同应变率下压缩的微观变形机制及力学响应。但是，关于位错密度梯度结构材料微柱压缩行为的研究还未见报道。本工作基于三维离散位错动力学模拟平台ParaDiS^[13]，研究位错密度梯度结构Cu单晶微柱的压缩行为，并进一步揭示加载方向对压缩响应的影响。

离散位错动力学模拟框架主要包括4部分：位错线的离散化、计算位错节点受力、求解位错动力学方程和处理位错的短程相互作用并进行拓扑改变^[13]。美国劳伦斯国家实验室开发的位错动力学模拟平台ParaDiS具有高效的并行计算能力，可开展较大规模的计算，并大大地减少了计算时间^[30]。本工作针对ParaDiS中节点力和塑性应变的计算进行简单介绍，详细描述可以参考文献[13]。

在离散位错动力学模拟中，位错会被离散成很多小的位错直线段，如图1所示。位错节点主要分为2类：一类是真实存在的物理节点，包括位错连接点(红色)和表面节点(绿色)。位错连接点一般是3个或3个以上位错交汇的节点，或者是2个不同滑移面的位错交汇的节点；另一类是虚拟节点，即同一滑移面的2个位错线段的交点(蓝色)，这些虚拟的节点可以通过一定的规则进行重划分。

每个节点i的受力(F_i)可以由位错系统整体储能(φ)的空间偏导数得到，即^[13]:

式中，X_i为节点i所在的空间坐标，φ是关于所有节点位置(X_j)、与节点j相连的位错段的Burgers矢量(b_jk)和外部施加表面牵引力(T^s)的函数，分为2部分：一部分是与位错芯的局部原子组态相关的能量(φ^c)；另一部分是与长程弹性变形相关的能量(φ^el)，即^[13]：

所以F_i的计算也可以分为2部分，即^[13]：

式中，F_i^c、F_i^el分别为节点i受到的位错芯能产生的力和长程弹性变形产生的力。

式(2)中的φ^c可以表示为^[13]：

式中，t_ij=l_ij /_{$‖{\mathit{l}}_{ij}‖$为位错段lij的单位向量，i和j表示位错上各个节点的编号，ω^c(bij, tij)表示与Burgers矢量和位错线方向相关的位错芯能。}

通常情况下，F_i^c可以认为是与节点i相连接的所有位错段l_ij由于晶格变形产生对节点i的作用力(f_ij^c)之和，结合式(1)和式(4)可得到^[13]：

式中，I₂为二阶单位张量。等式(6)右边的第一部分表示使位错段收缩的线张力，沿位错线方向；第二部分是使位错芯能最小而作用的旋转力，垂直于位错线方向。

由于长程弹性变形而作用在位错节点i的力(F_i^el)为^[31]：

式中，f_ij^el表示位错段l_ij受到的弹性作用力^[13]：

式中，N表示位错的形函数，l为位错段长度，f_ij^pk表示位错段l_ij在应力场中受到的Peach-Koehler力，包括外部载荷产生的应力场(σ_ext)、位错段自身产生的应力场(σ_s)以及其它位错段产生的应力场($\stackrel{n-\mathrm{1}}{\sum _{k=\mathrm{1}}}\stackrel{n}{\sum _{l=k+\mathrm{1}}}{\mathit{\sigma }}_{ij}^{kl}$)的叠加^[13]：

式中，[k, l]≠[i, j] or [j, i]。

在离散位错动力学中，由于设置的加载步长足够使位错达到稳定的状态，故而可以忽略位错段的惯性力的作用，位错段的运动方程通常可以简化为^[29]：

式中，B为阻尼系数，v_i为位错段节点i的速度。通过此运动方程，可追踪和描述每个位错段的时空演化。

当一个Burgers矢量为b的位错段在法向量为n的滑移面上滑移时，通过计算位错段滑移的面积可得位错滑移对塑性应变的贡献。单个位错段m在计算时间步dt内的滑移对整个计算模型的塑性应变增量(Δε_p^m)的贡献为^[32]:

式中，V为模型的体积，dS是位错段滑移时滑过的面积。计算时间步dt内位错滑移产生的总塑性应变增量(Δε_p)为：

本工作模拟的梯度位错密度结构材料的模型尺寸为1 μm×1 μm×3 μm，设置位错密度沿Z轴呈梯度变化。Senger等^[33]发现，当试样底边长度大于1 μm时，高宽比引起的硬化效应可以忽略不计。因此，进一步将试样沿Z轴方向分为3层，每层为边长为1 μm的正方体。设置单个位错源长度约为400 nm，在模型的底层、中层和高层随机生成100、200和300个Frank-Read位错源，进而形成一定的位错密度梯度结构。第a层的位错源密度(ρ_a)为：

式中，n_a为第a层的位错源数量，l_FRS为单个位错源长度。位错平均间距(d_a)为：

Zhang等^[34]的研究揭示了位错源密度对材料的变形机制有重要的影响，研究表明：当Frank-Read位错源长度远小于位错平均间距时，材料的塑性响应是位错源控制机制；当Frank-Read位错源长度远大于位错间距时(通常为2.5倍以上)，强化机制主要为林位错硬化。表1给出了各层的位错源数量、初始位错密度和位错间距。结合位错源长度为约400 nm可知，位错平均间距的最大值约为162 nm，因此，本工作中各层的塑性变形机制均为林位错硬化控制。

林位错硬化控制材料塑性变形过程时，塑性变形的临界分切应力(τ_cr)与位错密度(ρ)满足Taylor关系^[34]：

式中，α为Taylor强化因子，通常取值为0.2~0.5；μ为材料的剪切模量；b为Burgers矢量模。目前的离散位错动力学模拟框架采用的是各向同性假设，即忽略了材料弹性性质的各向异性。各向同性弹性常数(μ和Poisson比v)与fcc结构单晶材料各向异性弹性常数(C₁₁、C₁₂、C₄₄)之间满足Voigt-Reuss-Hill等效^[35]：

对于单晶Cu，C₁₁=168.4 GPa、C₁₂=121.4 GPa、C₄₄=75.4 GPa，可得Cu的剪切模量为54.6 GPa、Poisson比为0.324。fcc结构Cu单晶微柱压缩模拟的其它材料参数见表2。

Cu单晶微柱的X、Y、Z轴分别对应晶向[1$\overline{\mathrm{1}}$0]、[111]、[$\overline{\mathrm{1}}$$\overline{\mathrm{1}}$2]，位错密度沿Z轴呈梯度变化。本工作通过对位错密度梯度结构Cu单晶微柱分别施加垂直于梯度方向和平行于梯度方向的压缩载荷，研究加载方向对位错密度梯度结构材料力学性能的影响。图2为位错密度梯度结构模型在不同加载方向下的初始位错构型，图2a的加载方向沿着X轴方向(垂直于位错密度梯度方向)，图2b的加载方向沿着Z轴方向(平行于位错密度梯度方向)。研究^[29]表明，当应变率范围为1000~10000 s^-1时，材料力学行为的离散位错动力学模拟不会有明显的应变率效应。本研究通过应变率控制加载，应变率恒定为5000 s^-1，模型的所有表面都设置为当位错到达时即可逃逸的自由表面。

图3给出了位错密度梯度结构Cu单晶微柱在不同加载方向下的应力-应变曲线及位错密度演化曲线。沿X轴加载(加载方向垂直于位错密度梯度方向)时，Cu单晶微柱会在更低的应力水平时进入塑性屈服，而沿Z轴加载(加载方向平行于位错密度梯度方向)时，Cu单晶微柱会在较高的应力水平屈服。但是，在后续的塑性变形过程中，2个微柱模型的流动应力会达到同一个相对稳定的水平。

外部加载应力(σ)与作用在位错线上的分切应力(τ)之间的关系可以通过Schmid定律来表达^[34]：

式中，M是Schmid因子，ϕ为加载轴与滑移面法向的夹角，θ为加载轴与滑移方向的夹角。加载方向沿X ([1$\overline{\mathrm{1}}$0])和Z ([$\overline{\mathrm{1}}$$\overline{\mathrm{1}}$2])轴时，各个滑移系的Schmid因子见表3。图4进一步通过Thompson四面体描述不同加载方向时滑移系的开动情况。图中，红色的滑移面表示主滑移面，Thompson四面体的红色棱边表示主滑移方向。当沿X轴压缩时，有4个主滑移系，分别是：ACD滑移面(法向为($\overline{\mathrm{1}}$11)的滑移面)上的AD方向(滑移方向为[101])和AC方向(滑移方向为[0$\overline{\mathrm{1}}$1])，及BCD滑移面(法向为(1$\overline{\mathrm{1}}$1)的滑移面)上的BD方向(滑移方向为[011])和BC方向(滑移方向为[$\overline{\mathrm{1}}$01])，它们的Schmid因子均为$\mathrm{1}/\sqrt[]{\mathrm{6}}\approx$0.408。而沿Z轴压缩时，有2个主滑移系，为ACD滑移面(法向为($\overline{\mathrm{1}}$11)的滑移面)上的AC方向(滑移方向为[0$\overline{\mathrm{1}}$1])和BCD滑移面(法向为(1$\overline{\mathrm{1}}$1)的滑移面)上的BC方向(滑移方向为[$\overline{\mathrm{1}}$01])，Schmid因子同样为0.408。几何分析表明：尽管X轴加载和Z轴加载的最大Schmid因子均为0.408，但是塑性变形的临界屈服应力显著不同。根据单臂位错源模型，在相同位错结构下，主滑移系越多，位错源越容易开动，此时材料的屈服应力愈小。因此，沿垂直于位错密度梯度的X轴方向压缩时，4个主滑移系同时开动，材料整体更早进入塑性屈服。为了对这种现象的产生有一个更清晰的认识，接下来进一步对塑性变形和位错密度的时空演化进行分析，以揭示梯度位错结构材料压缩响应的各向异性。

图5显示的是当加载方向为X轴方向时，加载应变为0.05%、0.10%、0.15%、0.20%以及0.25%时的位错结构。可以看出，当应变较低时，材料主要发生弹性变形，位错在外力作用下产生可逆的弓曲变形。而当应变加载到约0.15%时，模型开始进入塑性阶段，位错源开动并增殖，从而导致了位错密度的增加，根据前文中的Taylor硬化关系，位错密度增加使得位错更难开动，从而产生了硬化。通过不同应变时刻的位错微结构可以发现，Burgers矢量为[110]和[1$\overline{\mathrm{1}}$0] (即图4中Thompson四面体的AB和CD方向)的位错在整个变形过程中没有开动，而其它滑移方向的位错段均有运动并滑移出表面，这与前文对各个滑移系位错的Schmid因子的分析结果相吻合。

图6显示了X轴加载时不同应变时刻的塑性应变云图。值得注意的是，顶部的塑性应变集中区是一种偶然现象，这是由于位错间的反应，生成了很长的容易开动的位错源。这种一端固定，另一端沿着2个滑移面交线运动，进而能够逐渐变长的位错源，在很小的应力下就可以开动并且增殖，而这一部分位错的滑移对其它位错的影响可以忽略。从图6可以看出，当模型加载到达塑性阶段时，位错的激活和滑移首先发生在位错密度最小的底层；随应变不断增加，位错密度较高的中间层和位错密度最高的顶层相继发生位错源的激活与开动，塑性变形以多个滑移系开动的形式进行。图7进一步给出了各层塑性应变和位错密度随应变的演化过程。从图7a可以发现，位错密度最小的底层在应变加载达到0.10%之后有一个较大的塑性应变增量，而中间层在0.15%的应变之后有一个较明显的塑性应变的增加，而最顶层的塑性应变的增加发生在应变0.2%之后。这与图6分析得出的塑性变形过程吻合。而从图7b可以发现，各层的位错密度随着应变加载的增加并没有很明显的变化。通过分析发现，位错密度在塑性变形过程中几乎不变的主要原因是，位错源激活导致位错密度增加和位错滑出试样表面导致的位错密度减少2个过程达到了动态的平衡。

图8给出了沿Z轴加载时，不同应变时刻的位错微结构。当应变较小时，与图5类似，位错段发生可逆的弓曲变形，材料整体处于弹性变形阶段。而当加载应变大于0.15%后，有大量Burgers矢量为[10$\overline{\mathrm{1}}$]和[0$\overline{\mathrm{1}}$1]的位错开动并逃逸出表面，而Burgers矢量为[1$\overline{\mathrm{1}}$0]的位错没有变化，其它Burgers矢量的位错段也有一定的激活并逃逸出表面的情况。这种现象也是由不同滑移系位错的Schmid因子不同而引起的，与表3给出的Schmid因子大小符合。图9进一步给出了不同应变下各个区域的塑性应变分布。从图中可以观察到，当应变达到0.15%时有位错滑移出表面，然后随着应变的继续增加，有更多的位错滑移出表面导致了塑性应变的增加，并且当加载方向为Z轴时，位错的滑移区域主要集中在(1$\overline{\mathrm{1}}$1)滑移面。

图10为位错密度梯度结构模型沿Z轴加载时，各层塑性应变和位错密度随应变的演化过程。对比图7a和图10a可以发现，与X轴加载时位错源逐层开动的情况不同，Z轴加载时，各层的塑性应变在应变较小时差异较小，当应变达到0.20%之后，各层的塑性应变的增长产生明显的差异，中间层的塑性应变增长率最大，其次是位错密度最低的底层，增长率最低的是位错密度最高的顶层；与此同时，Z轴加载时，位错密度也是中间层有较明显的增加，而其它2层变化很小。

通过对2种加载方向下的结果进行对比，可以得到不同加载方向对位错密度梯度结构材料的力学响应的差异：X轴加载时，加载方向垂直于位错密度梯度方向，各层在加载方向的受力相对比较独立，而底层位错密度较低，临界分切应力低，所以首先发生塑性变形。当底层的位错滑移到达中间层时，由于中间层的位错密度更大而需要更大的临界剪切应力才能开动，而此时底层的其它滑移系的位错也由于达到临界剪切应力而开始开动，而这也是产生多个滑移带的原因。随着加载应变的增加，位错受到的应力也不断增加，当应力水平达到中层的临界剪切应力时，中间层的位错也随之开动。最后，顶层的位错才会被激活并增殖。因此，微柱是逐层屈服的。Z轴加载时，加载方向平行于位错密度梯度方向，底层位错向着Z轴正方向滑移，随后会遇到更高位错密度层的阻碍，这也就导致了加载方向为Z轴时其临界剪切应力较大的现象。随着应力的增加，中间层位错的受力达到了临界剪切应力而发生了激活并滑移，此时由于滑移面的宽度限制，位错从中间层滑移到侧表面逃离，导致塑性变形首先在中间层发生，而由于加载方向为Z轴时一个滑移面只有一个最易开动的滑移系，这也就导致了单滑移的发生，进而整个模型只有一个滑移带。

从上述对比分析可知，加载方向对位错密度梯度结构的各个位错密度层的影响差异较大，当加载方向垂直于位错密度梯度方向时，位错的开动顺序是由位错密度最低的层逐层发展到位错密度更高层，并且伴随有多个滑移带产生；加载方向平行于位错密度梯度方向时，位错的开动主要发生在中间层，并逐渐向相邻层扩散，并且只有一个较集中的滑移带产生。

(1) 建立了位错密度梯度结构Cu单晶微柱的压缩变形的3D-DDD模型，该模型能展示位错在外载作用下的演化过程，从而更好地理解材料的变形机制。

(2) 加载方向对位错密度梯度结构Cu单晶微柱压缩响应有较大的影响，当加载方向垂直于位错密度梯度方向(X轴)时，会产生多个滑移带，整体的塑性变形更加均匀。而当加载方向平行于位错密度梯度方向(Z轴)时，位错的滑移主要发生在一个滑移面上。但是，当加载应变达到一定值后，2种不同加载方向下的流动应力的差异消失。

(3) 加载方向会对位错密度梯度结构材料不同位错密度层的塑性行为产生影响。当加载方向垂直于位错密度梯度方向时，位错开动顺序是由位错密度最低的层逐层发展到位错密度更高的层，并且在较低的应力水平就进入应变硬化阶段；当加载方向平行于位错密度梯度方向时，位错的开动主要发生在位错密度较高的中间层，并逐渐向相邻层扩散，进而使其需要更高的临界应力才能进入应变硬化阶段。