【研究意义】 断裂和高温力学性质变化是金属材料在使用过程中遇到的重要实际问题，其中研究材料断裂的宏观规律与微观结构之间的关系是金属物理学的重要内容之一^[1]。材料断裂不仅使构件完全失效，还往往产生灾难性的后果。因此，前人很早就发现断裂与材料微观结构有直接关系，并从不同的角度广泛研究断裂问题。【前人研究进展】目前已经有许多断裂模拟研究，其中既包括实验观测也有计算模拟，如Hosford等^[2]对铜的滑移线进行了观测，张跃等^[3]进行脆断微裂纹形核的原位观察，刘晓波等^[4]对铝裂纹扩展行为进行分子动力学模拟，Song等^[5]用相场法研究铁电体材料的裂纹尖端，Adelzadeh等^[6]研究裂纹与孔洞的相互作用和韧—脆断裂转变等。相对于实验模拟，计算模拟不仅能还原试验模拟现象，还便于观察实际实验中难以定格观察的瞬时的实验现象，所以说计算模拟的发展对于断裂模拟研究具有深远的实际意义。在这些模拟研究方法中，包含了单尺度分子动力学模拟^[7]，还有较为复杂的多尺度技术模拟。分子动力学模拟有一定局限性，即空间尺度上它通常只能模拟纳米尺度范围的原子、分子运动，而宏观过程一般发生在毫米、厘米甚至更大空间尺度上；且时间尺度上只能模拟飞秒(1 fs=10^-15 s)至皮秒(1 ps=10^-12 s)尺度上的原子、分子运动，而与宏观过程存在10⁶量级以上的差异^[8-9]。而相对于其它方法，Elder等^[10-12]基于密度泛函理论，结合分子动力学和相场方法的优势提出了晶体相场(PFC)方法模型。之后在此模型基础上，高英俊等^[13]研究韧性材料的微裂纹扩展与分叉，Hirouchi等^[14]在拉伸变形模拟中开发新的数值技术等。【本研究切入点】晶体相场(PFC)方法^[15-16]作为目前最先进的计算模拟方法，相对于分子动力学和传统相场方法而言，在结合材料微裂纹扩展和宏观情况的研究中更适合用于微纳米尺度上的裂纹扩展研究^[17-19]。【拟解决的关键问题】应用PFC方法研究15°初始晶向倾角条件下^[12]，不同方向单轴拉伸应变作用的位错发射与裂纹扩展规律，揭示纳米级韧性裂纹扩展机理和对材料断裂的影响。

PFC模型是通过对原子密度场变量取极小值得到的自由能密度函数。模型所构建的自由能函数产生周期对称的原子密度状态，能自洽地耦合弹性能、各向异性以及晶体学对称特性，还能反映原子空间尺度和扩散时间尺度上的现象，揭示原子效应。对于固态金属材料，原子的排列是有序的。因此，原子的周期性排列特征需要相场变量反映出来；同时，液相中的原子位置因为随时间随机变化，所以取时间的平均值看作常量。为了满足这两个要求的相场变量定义，我们用保守的原子密度场函数作为相场变量，表达式^[7]为

其中，等号右边第1项反映的是均匀相(如液相)的原子密度均匀分布，等号右边第2项和第3项反映晶格原子的周期排列结构特征。此时系统无量纲的自由能函数^[7]可以写为

式中，ρ为原子密度，∇²为拉普拉斯算符，γ表征系统的温度参数。对二维体系的的极小自由能的二维解作单模近似，可以得到公式(2)的一个稳定特解为

式中，$\mathit{q}=\frac{2\pi }{a}$，是在能量函数取极小值时所得；$A=\frac{4}{5}\left( {{\rho }_{0}}+\frac{1}{3}\sqrt{-15r-36{{\rho }_{0}}} \right)$，反映晶相原子密度周期结构的振幅；ρ₀为平均密度。

联立液相、条状相和三角相的原子密度函数和自由能公式，求解各相自由能的极小值，即可得到3种对应的平衡相。在二维系统中，通过单模近似可以获得平衡时各相的自由能密度，再利用自由能公切线法确定相图(图 1^[11])。

图 1 (Fig. 1)

S, T, L分别表示条状相，三角相，液相；(b)为(a)中虚线框的放大图 S, T and L represent stripe phase, triangular phase and liquid phase, respectively; (b) is enlarged map of the dashed block in (a) 图 1 PFC模型的二维相图[11] Fig.1 2D phase diagram of PFC model[11]

采用保守场Cahn-Hilliard动力学方程^[19-20]描述原子密度随时间的演化，该方程如下：

其中，ζ为高斯随机噪声，具有零平均值。本研究不需考虑ζ的作用。对公式(4)采用傅里叶转换方法，将实空间的数据转换到傅里叶空间中求解，再转换到实空间，最后通过MATLAB软件显示结果图像。详细求解步骤，见文献[20]。

在受到拉伸作用后，样品中所有格点的坐标位置都会发生相应的变化。因此，可以通过格点的坐标位置变换来反映应变的作用所引起的样品变形。设变形前的样品中格点位置的坐标用(x, y)表示，受拉伸作用引起变形后，格点位置的坐标用(x₁, y₁)表示。样品分别受y轴和x轴单方向拉伸作用引起形变后，格点位置的坐标分别用(x, y₁), (x₁, y)表示。样品在y方向和x方向分别受到应变ε_y和ε_x的拉伸作用下，变形前后坐标变换关系为y₁=y/(1+ε_x)和x₁=x/(1+ε_x)，其中应变$\varepsilon =\dot{\varepsilon }n\mathit{\Delta }t$，${\dot{\varepsilon }}$为无量纲的应变率，n为施加拉力应变的时间步数，Δt为时间步长。详细步骤见文献[21]。模拟中不涉及具体材料的物性参数，所用参数均已无量纲化处理^[11]，并将连续空间离散为四方格子，采用周期性边界条件。计算区域尺寸为1 024*512 gp(gp代表格子点数)，空间步长设为dx=dy=π/3，时间步长Δt=0.05，其它参数在下节中介绍。

根据图 1给出的相图，选取体系自由能密度函数的温度参数γ=-1.0，原子密度参数为ρ₀=0.49(图 1中的A点)，位于液相与三角格子固相共存区。在样品的中心位置，挖出一个缺口作为初始缺陷的位置，缺口为20*10 gp的长方形。缺口处的参数设为γ=-1.0，ρ₀=0.79(图 1中的B点)，其正好处在液相与三角格子固相共存区附近的液相区。选取该点参数有利于三角晶相结构向发射位错结构转变。本研究模拟所用单晶样品的原子排列方向与y轴夹角θ为15°。为观察裂纹的萌生和扩展情况，对设置的含有裂口的样品施加拉应变作用。首先经过50 000步弛豫，得到含缺口的稳态样品，为原子排列方向与y轴夹角θ为0°的缺口尺寸为20*10 gp的长方形(图 2，其放大区域给出实验的坐标体系)。

图 2 (Fig. 2)

图 2 初始样品 Fig.2 Initial sample

经过50 000步弛豫后，初始样品自由能达到稳定，即系统达到平衡(图 3a)。图 4为样品在y方向单轴拉伸作用下经过500 000步演化的位错发射情况，其中y方向单轴拉伸经历的时间步数t分别为350 000步、375 000步、425 000步、500 000步，对应的应变值为ε_y=10.5%、11.25%、12.75%、15%，对应步长的自由能曲线如图 3b所示。

图 3 (Fig. 3)

(a)自由能随时间变化曲线; (b)应力-时间变化曲线; (c)局部的应力-时间曲线 (a)Free energy change curve over time; (b)the curves of stress and time; (c)the local stress and time curve 图 3 自由能-时间曲线与应力-时间曲线 Fig.3 Free energy-time curve over time and stress-time curve

图 4 (Fig. 4)

(a)t=350 000;(b)t=375 000;(c)t=425 000;(d)t=500 000 图 4 y方向单轴拉伸作用下样品的裂纹扩展演化 Fig.4 y direction of the crack extension evolution of the sample under uniaxial tension

开始阶段，在y方向的拉伸作用下，在时间步数达到350 000步前样品所受应力不断增大，但在裂口处由于原子密度不均匀，裂口处阻力大，导致裂口应力集中，所以在350 000步前系统应力呈递减式增长(图 3b~c)；当时间步数达到350 000步，即应变ε_y达到10.5%，系统应力达到临界值时，裂口向两边发射滑移位错，这是典型的局部塑性变形表现。由图 3c可知，应力曲线呈振荡变化，曲线每两个相邻的峰谷差值反映了位错每滑移一个原子距离的驱动力，结合图 3c与图 4d可知，曲线峰谷间隔约为32峰谷值，与裂口左右两边位错几乎同时滑移约32个原子距离相对应。由图 4a可见，距裂口约2个原子距离的左右两边各出现1个位错，左右两边位错矢量方向相反，其柏氏矢量与水平方向夹角为15°，且此时没有出现与裂口相接的裂纹。随着y方向拉伸时间不断增大，裂口两边的位错沿着$[2\bar{1}\bar{1}]$和$[\bar{2}11]$方向的原子排列方向继续滑移；位错在滑移过程中沿滑移线诱发生成一系列孤立的空位，空位发展成空洞，然后空洞长大连通，形成裂纹，即裂纹随着位错运动方向而扩展运动，且此时裂纹边缘呈不光滑形状(图 4c~d)。在空位前端，左边的刃型位错滑移方向为向下左偏沿$[\bar{2}11]$方向，其柏氏矢量与水平方向夹角为15°；右边的刃型位错滑移方向为向上右偏沿$[2\bar{1}\bar{1}]$方向，其柏氏矢量与水平方向夹角为15°(图 4b~d)。

图 5为样品在x方向单轴拉伸作用下经过500 000步演化的位错发射情况，其中x方向单轴拉伸经历的时间步数分别为t=400 000步、420 000步、450 000步、500 000步，对应的应变值为ε_x=12%、12.6%、13.5%、15%，对应步数的自由能曲线如图 3a所示。与y方向单轴拉伸情况类似(图 4)，在400 000步前系统应力呈递减式增长；当时间步长达到400 000步，即应变ε_x达到12%，系统应力达到临界值时，裂口发生应力集中，裂口向上下两边发射滑移位错。距裂口约2个原子距离的上下两边各出现1个位错，上下两边位错矢量方向相反，柏氏矢量与y方向夹角为15°，且此时没有出现与裂口相接的裂纹(图 5a)。随着x方向拉伸时间不断增大，裂口上下两边的位错继续滑移；位错在滑移过程中沿滑移线诱发生成一系列孤立的空位，空位发展成空洞，然后空洞长大连通，形成裂纹，即裂纹随着位错运动方向而扩展运动，且此时裂纹边缘呈不光滑形状(图 5c~d)。在空位前端，裂口下边的刃型位错移运动方向为向下与y轴反方向、$[\bar{2}11]$方向夹角都为30°，其柏氏矢量与y轴反方向夹角为15°；裂口上边的刃型位错运动方向为向上与y轴正方向、$[2\bar{1}\bar{1}]$方向夹角都为30°，其柏氏矢量与y轴正方向夹角为15°(图 5b~d)。与y方向单轴拉伸情况类似同理，结合图 3c与图 5d可知，x方向单轴拉伸下应力曲线峰谷差值约为19，与裂口左右两边位错几乎同时滑移约19个原子距离相对应，应力曲线每两个相邻的峰谷差值反映了位错每滑移一个原子距离的驱动力。

图 5 (Fig. 5)

(a)t=400 000;(b)t=420 000;(c)t=450 000;(d)t=500 000 图 5 x方向单轴拉伸作用下样品的裂纹扩展演化 Fig.5 x direction of the crack extension evolution of the sample under uniaxial tension

晶向角为15°的样品，在x、y方向单轴拉应变作用下裂口起裂方向不同，并且都沿着拉应力方向作韧性裂纹扩展；拉应力单轴拉伸方向不影响裂纹韧脆扩展模式，但却影响韧性裂纹扩展方向。

裂纹作韧性裂纹扩展时，在裂尖前端发生位错滑移运动，并且位错在沿着拉应力单轴拉伸的原子排列方向的滑移轨道上诱发出一系列空位，空位发展成空洞，然后空洞生长扩大并连通，形成裂纹，裂纹随着位错运动方向扩展运动，此时裂纹边缘呈不光滑形状，此为典型的韧性裂纹扩展模式。