裂纹等微结构对材料性质有重要的影响，如裂纹扩展会引起材料的断裂等。为防止裂纹的扩展以及提高材料的使用寿命，人们开始重视孔洞裂纹愈合研究。例如，对陶瓷^[1-3]、聚合物^{[4, 5]}和金属材料^[6]进行加热可使其裂纹愈合；对TiB₂-SiC陶瓷^[7]可以通过氧化反应形成氧化膜来防止其表面氧化，而对聚合物^[8]则可通过氧化还原反应刺激溶胶-凝胶相变来控制样品愈合；Kumar等^[9]使用脉冲电流修复钢疲劳裂纹孔洞，Chen等^[10]采用高温空气热处理方法研究Al₂O₃/TiC/TiB₂陶瓷材料的裂纹愈合效果，等等。但在严苛的实验条件如高温、高压下，宏观实验过程中不能实时观察材料微观结构的变化情况，而计算机模拟方法具有安全性高且可以补充实验细节等优点，因此已成为研究材料微观结构演化的有效方法之一。例如，Li等^[11]通过分子动力学模拟实现了剪切力作用下铜板的裂纹愈合，Li等^[12]通过分子动力学模拟研究了压应力作用下金属铜内椭球孔洞的愈合，Xu等^[13]通过向错诱导实现纳米晶材料中的孔洞裂纹愈合。分子动力学模拟^{[14, 15]}是研究材料微观结构的重要方法之一，但该方法适用的时间尺度主要是原子振动的时间尺度(10^-15-10^-12 s)，难以拓展到原子扩散的时间尺度(10^-6 s)，导致模拟的加载应变速率达到10⁷-10⁹量级^[16]，远远偏离实际工业生产的应变速率(10^-2 s^-1)。

本研究使用晶体相场模型(Phase Field Crystal model, PFC模型)对材料纳米孔洞愈合特征进行分析。PFC模型是Elder等^{[17, 18]}基于密度泛函理论提出的，该模型已有广泛应用，如材料的外延生长^{[19, 20]}、相结构转变^[21-23]、裂纹扩展^[24-26]等现象的研究，但是目前有关PFC模型模拟孔洞愈合方面的文献报道较少。本研究利用PFC模型模拟在温度不变的情况下对样品施加压应变的纳米孔洞修复实验，对孔洞愈合的微观细节进行观察，旨在揭示受压应变下的材料纳米孔洞愈合的微观机理。

1 模拟方法 1.1 PFC模型

PFC模型不仅能够揭示晶体点阵周期性结构特性，而且还能模拟原子尺度下的纳微结构演化^[27-29]。对于固态晶体材料，可引入周期性原子密度函数ρ(r)作为场变量。原子密度函数的极大值对应于原子点阵的中心位置，其原子的位置呈周期性排列，则周期性原子密度函数ρ(r)的表达式为

$ \rho(r)=\sum\limits_\nolimits {n \cdot m} a_{n, m} e^{i G \cdot r}+\rho_0, $

(1)

式中，等号右边第1项反映了晶格原子的周期性排列结构特征，其中a_{n, m}为Fourier系数，i为虚数，G为倒格子矢量，r为空间坐标矢量；第2项ρ₀为均匀相(如液相)，表明原子密度均匀分布。系统无量纲的自由能函数可以写成：

$ F=\int\left\{\frac{\rho}{2}\left[\gamma+\left(1+\nabla^2\right)^2\right] \rho+\frac{\rho^4}{4}\right\} \mathrm{d} r^3, $

(2)

式中，γ为与温度有关的参数，$\nabla$²为拉普拉斯算符。在单模近似下，将六角晶相的自由能取极小值，可以求得方程(2)的一个稳定特解，即六角晶相的原子密度函数：

$ \begin{aligned} & \rho(r)=A\left[\cos (q x) \cos \left(\frac{q y}{\sqrt{3}}\right)-\frac{1}{2} \cos \left(\frac{2 q y}{\sqrt{3}}\right)\right]+ \\ & \rho_0, \end{aligned} $

(3)

式中，$A=\frac{4}{5} \rho_0+\frac{4}{15} \sqrt{-15 \gamma-36 \rho_0^2}$，反映了六角晶相原子密度的振幅；x、y表示密度波函数的空间坐标；$q=\frac{\sqrt{3}}{2}$，ρ₀为平均密度。本研究采用保守场动力学方程^[30]描述原子密度随时间的演化，该方程具体如下：

$ \frac{\partial \rho}{\partial t}=\nabla^2 \frac{\delta F}{\delta \rho}, $

(4)

式中，t为时间变量。采用半隐式傅里叶谱方法^[31]求解该动力学演化方程。

1.2 模拟实验样品的制备

模拟实验样品区域的边界条件设置成周期性边界条件，模拟区域为1 024×2 048的均匀网格，空间步长$\Delta x=\Delta y=\frac{\pi}{4}$，时间步长Δt=0.05，孔洞半径r=24。其中样品的温度参数γ=-0.5，原子平均密度ρ₀=0.36，孔洞处密度ρ₁=0.6；样品的弛豫时间步数为n=50 000 ts；晶粒取向角为0°。沿y轴对样品施加压应变，应变速率为$ \dot{\varepsilon}=6 \times 10^6$，应变增量$ \Delta \varepsilon=\dot{\varepsilon} \Delta t$，总施加应变步数n=50 000 ts。经过n时间步长的压缩，可以得到变形后的x轴和y轴方向的空间步长，分别如公式(5)和公式(6)所示：

$ \Delta x^{\prime}=\Delta x, $

(5)

$ \Delta y^{\prime}=\frac{\Delta y}{1+\varepsilon}=\frac{\Delta y}{1+n \dot{\varepsilon} \Delta t} 。$

(6)

2 结果与分析

如图 1所示，蓝色虚线表示总体系的自由能变化；黑色实线是自由能斜率曲线，表示样品在施加压应变过程中体系内部平均应力的变化趋势。曲线上的点(A₀, A, B, C……)对应后续图中的样品孔洞愈合演化时刻。自由能斜率曲线A₀点到A点阶段表示随着压应变的持续施加，曲线线性上升，该阶段变形服从虎克定律，表现为弹性变形。B点到K点阶段表示位错形核并发射，孔洞收缩，该过程表现为塑性变形。K点到L点阶段表示孔洞修复完成。后续按3个阶段分析材料纳米孔洞愈合的模拟结果。

2.1 纳米孔洞的形貌变化

图 2(a)为未加力前初始孔洞样品情况，对应于图 1中A₀点。在样品中心处预设一缺口，插图为缺口放大示意图，黄线区域即是缺口边界，该缺口是边长为4个格点的正六边形(l=4a)。由该插图放大图可知，缺口竖直方向厚度为7层原子层，用蓝色箭头标出。从时刻A₀点开始，沿y轴方向在缺口处加压应变。在A₀点→A点的过程中，随着压应变的增加，系统的自由能呈单调上升趋势。图 2(b)对应于图 1中A点时的样品情况。当加力时间为171 000步时，在缺口竖直方向上的第3层和第5层原子面开始凸出，第4层原子面凹陷，缺口呈阶梯状分布。如图 2(b)的插图(A2)所示，样品中心预设的缺口两侧出现两组位错对，用粉色的“⊥”刃型位错表示，并在附近用蓝色辅助线显示出位错的位置。

2.2 位错发射与孔洞收缩

图 3对应图 1中B点时的样品情况。如图 3(B1)所示，在181 000步时缺口的原第3层和原第5层原子面向第4层原子面挤压，原第4层原子面在x方向上两侧扩展并形成第一组位错团，原第4层原子面开始消失，整个缺口面积减小。图 3(B2)为将要发射的两个位错团，该位错团分别由3个位错组成，用粉色的“⊥”刃型位错表示。

图 4(a)对应图 1中C点时的样品情况。如图 4(a)所示，在压应变挤压作用下，缺口的第3层和第4层原子面开始发生扩展，长出凸口。B点时形成的位错团以位错发射形式从凸口两侧沿x轴方向水平滑移，构成滑开型位错。如图 4(a)的插图(C2)所示，a₁、a₂分别表示发射的2个位错团。位错发射后，两侧的a₁、a₂位错团与缺口之间出现狭缝连线区域，该区域形成4个双位错对。位错团a₁、a₂的滑移方向与压应变方向垂直，如黑色箭头所示。图 4(b)为C点时所对应样品的ε_yy应变分布情况，拉应变区域为黄色至红色，对应值为0-0.2；压应变区域为淡蓝色至深紫色，对应值为-0.2-0；绿色为无应变区域，数值为0，后续图中的应变场以该色标为准。位错团a₁、a₂的多余原子面区域颜色为蓝色，代表此区域原子受到压应变作用；位错团头部区域颜色为红色，则该区域受到拉应变作用。图 4(b)插图(C1)是ε_yy应变分布图叠加原子密度的放大图，黄色线围成的孔洞区域在ε_yy应变分布图中本应该是绿色无应变区域，但由于初始设置孔洞的密度并不为零，所以该区域存在应变。压应变在缺口附近的位错团处形成应变集中，提供位错团滑移所需要的能量。作用在位错团上的总应力为斥力，当斥力大于滑移阻力时位错团从凸口处发射，位错团运动方向与被修复原子面方向保持一致，位错团做加速运动远离凸口。

如图 5(a)所示，随着压应变的持续作用，位错团a₁、a₂与缺口之间的狭缝连线区域愈合，表面能下降，第一层原子面修复过程结束。孔洞的第3层和第4层原子面持续扩展并形成梯形结构。如图 5(a)的插图(D2)所示，在凸口处狭缝区域愈合过程中，位错团a₁、a₂转变成由两个位错组成的位错对。在裂纹扩展尖端形成位错对b₁、b₂，开启第二层原子面修复过程。如图 5(b)所示，已发射位错对仍处于应力集中状态，受压应变作用。

图 6(a)表示孔洞的第二原子面收缩完成。由图 6(a)的插图(E1)可见，样品中心缺口厚度降低为5层原子面，形成边长为3个格点的正六边形。在图 6(a)的插图(E2)中，共发射4个位错对a₁、a₂、b₁和b₂，位错对沿x轴方向滑移。受压应变作用，b₁和b₂位错对与缺口之间的表面缝隙提前愈合。图 6(b)中缺口两侧的位错对a₁、a₂、b₁和b₂之间的距离逐渐增加。如图 6(b)的插图(G6)所示，缺口的第2层和第3层原子面开始扩展，此时缺口扩展的长度比之前扩展的长度长。凸口扩展的尖端形成2个位错对c₁、c₂，在凸口扩展的尖端和初始端之间还存在2个位错，凸口形成的位错对会对凸口产生一个回复力，因此扩展的尖端宽度大于初始端宽度。随着压应变的增加，作用在位错上的总应力增大，当其大于滑移阻力时，缺口扩展的尖端处会发射位错对c₁、c₂，该过程与修复第一层原子面阶段相似。凸口发射位错后，在压应变作用下，裂纹张开位移减小，凸口初始端与尖端之间的间隙提前愈合，导致凸口扩展的驱动力下降，凸口扩展区域被逐渐挤压而消失，而凸口的闭合导致自由能增长趋势的降低。E点到G点阶段，在凸口附近的区域内，由于位错受吸引力影响很难发射，因此位错发射需要克服一定的晶格阻力。

图 7(a)对应图 1中H点的位置情况。如图 7(a)的插图(H5)所示，G→H阶段持续进行凸口扩展区域的愈合过程，样品表面能减少，发射位错对间的距离与被修复原子面克服滑移阻力所需做功增加，生成新的位错对d₁和d₂。新位错对发射前，凸口的扩展区域便已经愈合。图 7(a)的插图(H6)显示出样品中心缺口厚度减少为3层原子面。图 7(b)中有8处应力集中处，当压应变达到发射位错的临界应力强度时，开始发射位错。

图 8(a)对应图 1中I点的位置情况。由图 8(a)中插图可知，由于样品的边界性周期条件，发射的位错对a₁和a₂与边界接触后运动方向发生变化，导致两者按原来相反的滑移方向运动。如图 8(a)的插图(I7)和(I8)中所示，孔洞已经收缩为两层原子面，上底为4个格点，下底为3个格点，发射新的位错对e₁、e₂，在缺口底部存在两个位错。第5层原子面修复过程非常短，没有产生大幅度的裂纹扩展现象。图 8(b)中有10处应力集中区域，这些区域间的距离随应变施加而增加。

图 9对应图 1中J点的位置情况。如图 9所示，中心缺口形状变为上底为3个格点、下底为5个格点的双层阶梯。I→J点过程孔洞的第一层原子面收缩，第二层原子面发生扩展，并在端点处生成2个位错对。

2.3 孔洞修复完成

由图 10中插图(K7)和(K8)可见，缺口附近发射f₁和f₂位错对，修复了孔洞的第6层原子面；而孔洞最后一层原子面的修复不依靠位错发射完成，而是由其他平面的原子逐渐向缺口处迁移实现。发射位错对f₁和f₂所在区域受压应变作用，该原子面提前闭合，扩展驱动力下降，最后孔洞消失，裂纹修复完成。图 11表示到444 000步时，厚度为7层原子面的孔洞修复完成。由图 11(a)中插图(L5)可见，原孔洞所在虚线圆区域已愈合成完整的单晶。图 11(b)中有8处应力集中区域，与图 11(a)已经发射的位错对位置一致。

3 结论

裂纹孔洞愈合过程分为3个阶段：缺口两侧局部凸起发生孔洞扩展与位错的产生、位错发射导致缺口厚度的台阶数减少、缺口区域的消失。位错发射机制是孔洞愈合的主要表现形式，当样品缺口厚度大于2层原子面时，缺口两侧会同时发射位错以实现孔洞缩口，最后1层原子面的修复会以其他平面原子迁移填充的方式进行。当压应变达到临界应力条件时，凸口会发射位错，并出现几个格点大小的应力集中区域，该区域的位置会随着位错运动而移动。