【研究意义】纳米晶体材料是近年材料学方面的研究热点，它具有良好的形变特性，因此其形变机理引起了人们的广泛重视^[1]。要深入研究材料的变形及其特性与机理，需要把宏观分析与微纳观分析结合起来，在更深的层次上找到其变形机制^[2]。现在，对材料变形的研究早已深入到微纳观层次。纳米级微裂纹的形核与扩展在金属材料微观缺陷中普遍存在，并严重降低了金属材料的使用寿命^[3]，因此，研究纳米级裂纹的萌生和扩展对预防材料的断裂，提高材料的使用寿命具有重要意义^[4]。【前人研究进展】近几年，基于密度泛函理论建立的晶体相场方法(PFC)^[5]，能很好地用于描述晶界和位错在扩散时间尺度下的运动特征^[6]，可以用于模拟纳米级的微观结构和演化过程^[6-8]，并用于研究晶体弹性及粘弹性行为。【本研究切入点】在当前的实验测量条件下，对材料的纳米级行为很难原位观测^[9]，因此，发挥计算模拟实验的优势，应用其研究微纳米尺度的结构就显得极为迫切和重要^[10-12]。【拟解决的关键问题】应用PFC方法模拟材料弹塑性变形，揭示该过程的原子位移运动特征，研究材料的微结构和弹性及粘弹性行为特征。

PFC模型系统无量纲的自由能函数可以写成^[8]

式中，ρ为局域原子密度；r为与温度有关的唯象参数；∇²为拉普拉斯算子。在单模近似下，可以求得式子(1) 的一个稳定特解

式中，ρ₀为平均原子密度， $A=\frac{4}{5}{{\rho }_{0}}+\frac{4}{15}\sqrt{-15r-36{{\rho }_{0}}^{2}}$ 反映固相原子密度周期结构的振幅，由能量函数取极小值，得到ρ。

采用保守场Cahn-Hilliard动力学方程^[13]描述原子密度随时间的演化。该方程具体表示如下：

式中，ρ为局域原子密度，r为与温度有关的唯象参数，表征体系的过冷度。∇²为拉普拉斯算子，t为时间变量。α是一个运动的参量，β是一个可调的阻尼参数。

对无量纲动力学演化方程(3) 采用二次半隐式傅里叶伪谱方法^[14]求解，其离散形式为

式中，${\hat{\rho }}$_k(t)，${\hat{\rho }}$_k(t+Δt)，${\hat{\rho }}$_k(t-Δt)为傅里叶空间下的原子密度，k为傅里叶空间的波矢^[13]。

模拟设置一个长方形[0,L=360]×[0,H=100]，使它满足四周边界都是零通量条件。固相基体用三角相点阵表示，选取晶体平均原子密度为ρ₀=0.31，温度参数r=-0.4。计算模拟区域为360Δx×100Δy，单位格子长度Δx=Δy=π/8。对模拟所用参数进行无量纲化处理，并将连续空间离散为正方格子分布。计算时四周边界采用周期性边界条件，采用弛豫方法使系统能量达到稳定状态，样品弛豫时间步长为50 000步，最终得到的初始样品如图 1所示。

图 1 (Fig. 1)

图 1 初始样品的二维图 Fig.1 Two dimensional diagram of initial sample

在变形过程中，如图 2所示，x方向空间步长保持不变，y方向空间步长随着应变速率在每一时间步长下都有一增量d=${\dot{\varepsilon }}$nΔt，其中应变速率${\dot{\varepsilon }}$=6.5×10^-4,n为时间步长数，Δt为时间步长，Δt=0.001。采用补偿函数方法^[14]施加应力作用于样品。假设原子的垂直线最初位于x=L的位置，并进一步假设初始原子密度轮廓沿着已知的原子线n_clamp(y)=n_OM(x=L,y)。为实现水平速度ω=ω(t)e_x，通过调整F=F₀+V(t)，施加拉应力作用于样品。表 1中,A样品表示进行弹性模拟,B样品表示粘弹性模拟。

图 2 (Fig. 2)

图 2 沿x轴方向施加拉应变的示意图 Fig.2 Schematic diagram of tensile strain applied along x axis

表 1(Table 1)

表 1 不同的样品模拟参数 Table 1 Simulation parameters of different samples 样品Sample r ρ0 Δt α β A -0.4 0.31 0.001 15 0.9 B -0.4 0.31 0.001 15 9
	表 1 不同的样品模拟参数 Table 1 Simulation parameters of different samples

由图 3可见：1) 当t=120 000时，对应图 3a，大约2/3处位移为0，从右往左，位移慢慢增加，应变随之增加，位移从0变化到4，最大位移为4，对应图 4中曲线1。2) 当t=140 000，对应图 3b，位移从0增到5，最大位移为5，对应图 4中曲线2。3) 当t=160 000时，对应图 3c，位移从0变化到6，最大位移为6，对应图 4中曲线3。由图 4可知：位移随原子位置增加而减小。综上所述，当β=0.9时，随着力F的施加，位置变化小，位移变化也小，位置变化大，位移也相应变大，位置越远，应变越大，位移和位置成正比关系，形成弹性响应关系，晶体响应是有弹性的。

图 3 (Fig. 3)

图 3 弹性模拟样品不同时刻原子位置与相对应的位移 Fig.3 Atomic position and relative displacement at different moments of elastic sample

图 4 (Fig. 4)

图 4 不同时刻的原子位置和位移关系曲线 Fig.4 The atomic positions and displacement curves at different times

由图 5可见：1) 当t=120 000时，从右到左，位移缓慢增加，位移从0变化到3，最大位移为3，对应图 6中曲线1。2) 当t=140 000，对应图 5b，位移从0增到4，最大位移为4，对应图 6中曲线2。3) 当t=160 000时，对应图 5c，位移从0变化到7，最大位移为7，对应图 6中曲线3。由图 6可知，位移随原子位置增加而减小，但在位置3处，出现转折变化。综上所述，当β=9时，随着力F的施加，刚开始时，位置变化小，位移变化也小且几乎不变，当该变量达到一定程度，位移有一个跃迁的过程，位移变化相应也大，曲线类似于折线。β从0.9增加到9导致晶体出现粘弹性行为，形成粘弹性响应关系。

图 5 (Fig. 5)

图 5 粘弹性模拟样品不同时刻原子位置与相对应的位移 Fig.5 Atomic position and relative displacement at different moments of viscoelastic sample

图 6 (Fig. 6)

图 6 不同时刻的原子位置和位移关系曲线 Fig.6 The atomic positions and displacement curves at different times

应用PFC模型模拟研究材料弹塑性变形，得到如下结论：1) 当β=0.9时，在力F的作用下，此时晶体响应是有弹性的。位置距离越远，其对应的应变越大，位移变化与原子对应位置大致成正比，具有弹性关系；2) 当β=9时，在力F的作用下，位移变化与原子对应位置关系不服从弹性关系而是服从粘弹性响应关系。表明改变阻尼参数β后，可以采用晶体相场模型模拟晶体的弹性及粘弹性行为。