气孔与晶界相互作用的相场模拟

图1 相场模型示意图

Fig.1 Schematic of phase field model (C—concentration field variable; C_α—concentration in grain; C_β—concentration in pore; η $_{1}^{α}$ , η $_{2}^{α}$ —orientation field variable of grains 1, 2, respectively; η^β—orientation field variable of pore )

相场变量η $_{i}^{α}$ 、η^β、C在不同取向的晶粒以及晶粒与气孔之间是连续变化的，即相场变量在晶界和相界上具有一个扩散型边界^[19]。例如：在晶粒1和晶粒2的晶界处，η $_{1}^{α}$ 的取值从1连续变为0，η $_{2}^{α}$ 的取值从0连续变为1；在晶粒1和气孔的相界处，η $_{1}^{α}$ 的取值从1连续变为0，η^β的取值从0连续变化为1，C的取值从C_α连续变化为C_β。

1.2 自由能密度函数的构造

根据Ginzburg-Landau理论^[26]，当相场变量在空间中变化时，体系的总自由能密度F包含与相场变量有关的体自由能密度和与相场变量梯度有关的扩散界面自由能密度。因此，F的形式如下：

F = \int_{v} [f (η_{i}^{α}, η^{β}, C) + \frac{κ_{α}}{2}

\sum_{i = 1}^{p} {|\nabla η_{i}^{α}|}^{2} +

\frac{κ_{β}}{2} {|\nabla η^{β}|}^{2} + \frac{κ_{c}}{2} {|\nabla C|}^{2}] d^{3} r

(1)

式中，等号右边积分号内的第一项f为体自由能密度函数，后面三项为梯度自由能密度；κ_α、κ_β和κ_c为梯度参数；r为位置坐标。

f应保证满足以下的极值条件：当相场变量(C, η $_{1}^{α}$ , η $_{2}^{α}$ , …, η $_{p}^{α}$ , η^β)=(C_α, ±1, 0, …, 0, 0)，(C_α, 0, ±1, …, 0, 0)，…，(C_α, 0, 0, …, ±1, 0)和(C_β, 0, 0, …, 0, ±1)时，f取到2(p+1)个的极小值。换句话说，p个不同取向的晶粒以及气孔是f的势阱。为满足这个极值条件，本工作在Chen等^[27]提出的经典相场模型的体自由能密度函数基础上，构造了新的体自由能密度函数f_new，其形式如下：

f_{n e w} (η_{i}^{α}, η^{β}, C) = A (C - C_{α})^{2} {(C - C_{β})}^{2} + \frac{D_{α}}{4} {(C - C_{β})}^{4} +

\sum_{i = 1}^{p} (- (\frac{γ_{α}}{2}) {(C - C_{β})}^{2} (η_{i}^{α})^{2} + (\frac{δ_{α}}{4}) (η_{i}^{α})^{4}) +

\frac{D_{β}}{4} {(C - C_{a})}^{4} - (\frac{γ_{β}}{2}) {(C - C_{α})}^{2} (η^{β})^{2} +

(\frac{δ_{β}}{4}) (η^{β})^{4} + \frac{ε}{2} (\sum_{i = 1}^{p} (η_{i}^{α})^{2} (η^{β})^{2} + \sum_{i = 1}^{p} \sum_{j > i}^{p} (η_{i}^{α})^{2} (η_{j}^{α})^{2})

(2)

式中，A、D_α、γ_α、δ_α、D_β、γ_β、δ_β、ε均为唯象参数。模拟时需选择一组合适的唯象参数以保证f_new满足极值条件。下文1.4节将对极值条件进行详细讨论。

1.3 相场方程的演化

η $_{i}^{α}$ 、η^β和C的演化方程分别为Allen-Cahn方程和Cahn-Hilliard扩散方程^[28]：

\frac{\partial η_{i}^{α}}{\partial t} = - L^{α} \frac{δ F}{δ η_{i}^{α}} = - L^{α} [\frac{\partial f_{n e w}}{\partial η_{i}^{α}} - κ_{α} \nabla^{2} η_{i}^{α}]

(i = 1, 2, \dots, p)

(3)

\frac{\partial η^{β}}{\partial t} = - L^{β} \frac{δ F}{δ η^{β}} = - L^{β} [\frac{\partial f_{n e w}}{\partial η^{β}} - κ_{β} \nabla^{2} η^{β}]

(4)

\frac{\partial C}{\partial t} = \nabla \cdot (M \nabla \frac{δ F}{δ C})

(5)

式中，t为时间；L^α、L^β和M为动力学参数，M与张量形式的表面扩散系数D^s有关：

M = \frac{ν_{m} D^{s}}{R T}

(6)

式中，ν_m为摩尔体积；R为气体常数；T为温度。

本工作中采用的张量形式的D^s如下：

\{\begin{array}{l} D^{s} = D^{s} C^{2} (1 - C^{2}) T^{s} \\ T^{s} = I - n_{s} \otimes n_{s} \\ n_{s} = \frac{\nabla C}{|\nabla C|} \end{array}

(7)

式中，D^s是标量形式的表面扩散系数；T^s是表面投影张量，它的作用是保证表面扩散仅发生在晶粒与气孔接触表面的切线方向；I为单位张量；n_s为表面单位法向量。

本工作基于有限差分方法编写程序对上述模型进行模拟，Allen-Cahn方程和Cahn-Hilliard方程的求解采用显式Euler算法，Allen-Cahn方程中Laplace项的求解采用了五点差分法。

本工作引入一个可视化变量φ^[18]，结合扩散界面的特征可知：在晶粒内部φ=1；在扩散界面处0<φ<1；在气孔内部φ=0。通过将φ的不同数值转化为不同的颜色，即可直观地展示晶粒的生长演化过程。

φ = \sum_{i = 1}^{p} {(η_{i}^{α})}^{2}

(8)

1.4 模型中唯象参数的选择

为了保证f_new在不同取向晶粒内部和气孔内部取得相同极小值，需要选择一组合适的唯象参数。为了便于讨论，需要将式(2)拆分为3部分：f₁、f₂、f₃，令f_new=f₁+f₂+f₃，即：

\{\begin{array}{l} f_{1} = A (C - C_{α})^{2} {(C - C_{β})}^{2} \\ f_{2} = \frac{D_{α}}{4} {(C - C_{β})}^{4} + \sum_{i = 1}^{p} (- (\frac{γ_{α}}{2}) {(C - C_{β})}^{2} (η_{i}^{α})^{2} + \\ (\frac{δ_{α}}{4}) {(η_{i}^{α})}^{4}) + \frac{D_{β}}{4} {(C - C_{a})}^{4} - \\ (\frac{γ_{β}}{2}) {(C - C_{α})}^{2} (η^{β})^{2} + (\frac{δ_{β}}{4}) (η^{β})^{4} \\ f_{3} = \frac{ε}{2} (\sum_{i = 1}^{p} (η_{i}^{α})^{2} (η^{β})^{2} + \sum_{i = 1}^{p} \sum_{j > i}^{p} (η_{i}^{α})^{2} (η_{j}^{α})^{2}) \end{array}

(9)

再令f₂=g^α+g^β，g^α、g^β分别为：

g^{α} (C, η_{i}^{α}) = \frac{D_{α}}{4} {(C - C_{β})}^{4} + \sum_{i = 1}^{p} (- (\frac{γ_{α}}{2}) {(C - C_{β})}^{2} \cdot

(η_{i}^{α})^{2} + (\frac{δ_{α}}{4}) (η_{i}^{α})^{4})

(10)

g^{β} (C, η^{β}) = \frac{D_{β}}{4} {(C - C_{a})}^{4} - (\frac{γ_{β}}{2}) {(C - C_{α})}^{2} (η^{β})^{2} + (\frac{δ_{β}}{4}) (η^{β})^{4}

(11)

假设p=1，此时函数g^α可简化为g^α(C, η^α)=D_α(C-C_β)⁴/4-γ_α(C-C_β)²(η^α)²/2+δ_α(η^α)⁴/4。二元函数g^α(C, η^α)在C=C_α、|η^α|=1取得极小值的必要条件是2个偏导数均为0，即：

{\frac{\partial g}{\partial C}|}_{\begin{array}{l} |η^{α}| = 1 \\ C = C_{α} \end{array}} = D_{α} {(C - C_{β})}^{3} - γ_{α} (C - C_{β}) (η^{α})^{2} = 0

(12)

{\frac{\partial g}{\partial η^{α}}|}_{\begin{array}{l} |η^{α}| = 1 \\ C = C_{α} \end{array}} = - γ_{α} {(C - C_{β})}^{2} η^{α} + δ_{α} (η^{α})^{3} = 0

(13)

求解上述2个方程可得唯象参数D_α、γ_α、δ_α满足的关系如下：

D_{α} {(C_{α} - C_{β})}^{2} = γ_{α}, γ_{α} {(C_{α} - C_{β})}^{2} = δ_{α}

(14)

根据对称性可知，唯象参数D_β、γ_β、δ_β满足的关系如下：

D_{β} {(C_{β} - C_{α})}^{2} = γ_{β}, γ_{β} {(C_{β} - C_{α})}^{2} = δ_{β}

(15)

将上述关系代入函数g^α、g^β，得到：

g^{α} (C, η^{α}) = {[{(\frac{C - C_{β}}{C_{α} - C_{β}})}^{2} - {(η^{α})}^{2}]}^{2}

(16)

g^{β} (C, η^{β}) = {[{(\frac{C - C_{α}}{C_{β} - C_{α}})}^{2} - {(η^{β})}^{2}]}^{2}

(17)

此时，函数g^α、g^β分别在(η^α)²=[(C-C_β)/(C_α-C_β)]²，(η^β)²=[(C-C_α)/(C_β-C_α)]²条件下取到极小值0，但此条件并不满足f_new的极值条件。根据极值条件，f_new应该只在C=C_α且|η^α|=1或C=C_β且|η^β|=1处取到极小值。因此需要对f₂=g^α+g^β的极值条件进行限制。

为了限制f₂=g^α+g^β的极值条件，引入f₁=A(C-C_α)²·(C-C_β)²，A是一个大于0的常数。可以证明，f₁是一个双势阱函数，分别在C=C_α和C=C_β处取得极小值0。因此f₁+f₂的极值范围缩小到C=C_α且|η^α|=1或C=C_β且|η^β|=1。

上述对f₁+f₂的讨论都是基于p=1的假设，但是在通常的情况下p>1。此时不难证明f₁+f₂的值在C=C_α，|η_i^α|=1 (i=1, 2, 3, …, p)时最小。例如：当p=2时，g(C=C_α, η $_{1}^{α}$ =1, η $_{2}^{α}$ =1)=-δ_α/4<g(C=C_α, η $_{1}^{α}$ =1, η $_{2}^{α}$ =0)=g(C=C_α, η $_{1}^{α}$ =0, η $_{2}^{α}$ =1)=0，这显然不符合f_new的极值条件。为解决此问题，f₃引入了取向场变量的交叉项，交叉项前面的参数ε应满足以下条件^[18]：

ε > δ_{α} / 2

(18)

因此，f_new=f₁+f₂+f₃满足极值条件的充分必要条件如下：

\{\begin{array}{l} D_{α} {(C_{α} - C_{β})}^{2} = γ_{α}, D_{β} {(C_{β} - C_{α})}^{2} = γ_{β} \\ γ_{α} {(C_{α} - C_{β})}^{2} = δ_{α}, γ_{β} {(C_{β} - C_{α})}^{2} = δ_{β} \\ A > 0, ε > δ_{α} / 2 \end{array}

(19)

2 模型分析

在本工作研究的体系中共有2种界面：一种是2个晶粒接触位置附近的晶界，另一种是晶粒相与气孔相接触位置附近的相界。根据Ginzburg-Landau理论^[26]，体系中扩散界面的界面宽度和能量密度与自由能密度函数中梯度项参数相关。晶界的能量密度σ_gb在二维形式下的计算表达式^[27]为：

σ_{g b} = \int_{- \infty}^{\infty} [Δ f_{n e w} (η_{i}^{α}, η_{j}^{α}, C) + \frac{κ_{α}}{2} {(\frac{d η_{i}^{α}}{d x})}^{2} +

\frac{κ_{α}}{2} {(\frac{d η_{j}^{α}}{d x})}^{2} + \frac{κ_{c}}{2} {(\frac{d C}{d x})}^{2}] d x

(20)

Δ f_{n e w} (η_{i}^{α}, η_{j}^{α}, C) = f_{n e w} (η_{i}^{α}, η_{j}^{α}, C) - f_{n e w} (η_{i, e}^{α}, η_{j, e}^{α}, C_{α}) -

(C - C_{α}) {(\frac{\partial f_{n e w}}{\partial C})}_{η_{i, e}^{α}, η_{j, e}^{α}, C_{α}}

(21)

式中，η $_{i, e}^{α}$ 、η $_{j, e}^{α}$ 分别为η $_{i}^{α}$ 、η $_{j}^{α}$ 在平衡位置的取值；f_new (η $_{i, e}^{α}$ , η $_{j, e}^{α}$ , C_α)为f_new在平衡位置的极小值。

相界的能量密度σ_int在二维形式下的计算表达式^[27]为：

σ_{i n t} = \int_{- \infty}^{\infty} [Δ f_{n e w} (η_{i}^{α}, η^{β}, C) + \frac{κ_{α}}{2} {(\frac{d η_{i}^{α}}{d x})}^{2} +

\frac{κ_{β}}{2} {(\frac{d η^{β}}{d x})}^{2} + \frac{κ_{c}}{2} {(\frac{d C}{d x})}^{2}] d x

(22)

Δ f_{n e w} (η_{i}^{α}, η^{β}, C) = f_{n e w} (η_{i}^{α}, η^{β}, C) - f_{n e w} (η_{i, e}^{α}, η_{e}^{β}, C_{α}) -

(C - C_{α}) {(\frac{\partial f_{n e w}}{\partial C})}_{η_{i, e}^{α}, η_{e}^{β}, C_{α}}

(23)

利用有限差分数值方法研究晶界参数、相界参数和梯度参数的关系。模拟区域的大小为40×40格点，模拟区域的边界条件为：左右为非周期性边界，上下为周期性边界。模拟选取的无量纲参数如表1所示，模拟采用“固定变量法”，即模拟某一个梯度参数的影响时，其余的梯度参数保持初始值不变。

表1 模拟选取的无量纲参数

Table 1 Dimensionless parameters in simulation

Parameter	Symbol	Value
Phenomenological parameter	A	2.5
	D_α, D_β	1.08
	γ_α, γ_β	1.04
	ε	3
Kinetic parameter	L^α, L^β	1
	M	5
Gradient parameter	κ_α, κ_β, κ_c	2
Space step	Δx	1
Time step	Δt	0.001

Note:M corresponds to the scalar form of M

新窗口打开| 下载CSV

2.1 晶界宽度及晶界能量

为研究晶界参数与梯度参数的关系，将模拟区域等分为左右2份，分别表示2个不同取向的晶粒。左右2侧相场变量的初始值分别为：左侧η $_{1}^{α}$ =1，η $_{2}^{α}$ =0，C=C_α，右侧η $_{1}^{α}$ =0，η $_{2}^{α}$ =1，C=C_α。待体系演化到稳态时，做出晶界附近相场变量的变化函数，如图2所示。其中，图2a为η $_{i}^{α}$ 、η $_{j}^{α}$ 的变化函数，图2b为C的变化函数。可以看出，κ_α的取值对η $_{1}^{α}$ 、η $_{2}^{α}$ 影响显著，κ_α越大η $_{1}^{α}$ 、η $_{2}^{α}$ 变化越平缓，晶界宽度越大，κ_α对C影响很小。根据对称性，κ_β对η^β具有相同的影响。

图2

图2 不同梯度参数κ_α下取向场变量η $_{1}^{α}$ 、η $_{2}^{α}$ 和浓度场变量C的平衡取值

Fig.2 Equilibrium profiles for phase field variables at different κ_α(a) orientation field variables η $_{1}^{α}$ , η $_{2}^{α}$ (b) concentration field variable C

图3为σ_gb和晶界宽度l与κ_α平方根的线性拟合。可以看出，σ_gb和l与κ_α平方根成正比。Moelans等^[29]利用数值方法研究了单相相场模型(体系中不含第二相气孔或者颗粒)晶界宽度和能量密度与梯度项参数的关系，模拟结果表明晶界宽度和能量密度与梯度项参数的平方根成正比，这与本工作的模拟结果基本一致。然而，本工作中的3类梯度项参数的平方根对l和σ_gb影响的贡献并不完全相同。

图3

图3 晶界能量密度σ_gb和晶界宽度l与κ $_{α}^{1 / 2}$ 的线性拟合

Fig.3 Linear fitting curves of energy density and width of grain boundary vsκ $_{α}^{1 / 2}$ (l—width of grain boundary, σ_gb—energy density of grain boundary)

2.2 相界宽度及相界能量

为研究相界参数与梯度参数的关系，将模拟区域等分为左右2份，分别表示晶粒相和气孔相。左右2侧相场变量的初始值分别为：左侧η $_{1}^{α}$ =1，η^β=0，C=C_α，右侧η $_{1}^{α}$ =0，η^β=1，C=C_β。

图4为在不同κ_c下体系达到平衡状态时相界附近相场变量的变化函数。其中，图4a为η $_{1}^{α}$ 、η^β的变化函数，图4b为C的变化函数。图5表示相界能量密度与κ_α、κ_β、κ_c 3者平方根的线性拟合。

图4

图4 不同梯度参数κ_c下相场变量的平衡取值

Fig.4 Equilibrium profiles for phase field variables at different κ_c(a) orientation field variables η $_{1}^{α}$ , η^β(b) concentration field variable C

图5

图5 相界能量密度σ_int与κ $_{α}^{1 / 2}$ (κ $_{β}^{1 / 2}$ )、κ $_{c}^{1 / 2}$ 的线性拟合

Fig.5 Linear fitting curves of energy density of phase boundary σ_intvsκ $_{α}^{1 / 2}$ , κ $_{β}^{1 / 2}$ or κ $_{c}^{1 / 2}$

对比图2a和图4a可以看出，κ $_{α}^{1 / 2}$ 对η $_{i}^{α}$ 在扩散界面处的取值有较大影响，根据方程(3)、(4)的对称性可知，κ $_{β}^{1 / 2}$ 对η^β在扩散界面处的取值有较大影响；结合图2b和图4b分析可得，κ_c对C在扩散界面处的取值有较大影响。进一步分析图4a可知：κ_α和κ_β相同的情况下改变κ_c，相界界面宽度不变且与晶界界面宽度相同。对比图3和图5可以看出，晶界能量密度和晶界宽度与κ $_{α}^{1 / 2}$ 成正比。相界能量密度与κ $_{α}^{1 / 2}$ 、κ $_{β}^{1 / 2}$ 、κ $_{c}^{1 / 2}$ 均成正比，其与κ $_{α}^{1 / 2}$ 、κ $_{β}^{1 / 2}$ 的比例系数相同且略小于与κ_c^1/2的比例系数。

3 气孔与晶界的相互作用

3.1 晶界曲率对其移动速率的影响

为研究自由晶界(不含气孔)的移动规律，本工作构建一个大小为75×75格点的模拟区域，以模拟区域中心为圆心，做一个半径为20格点的圆将模拟区域分为2部分，圆内和圆外分别表示2个不同取向的晶粒。圆内相场变量的初始值为η $_{i}^{α}$ =1，η $_{j}^{α}$ =0，C=C_α，圆外相场变量的初始值为η $_{i}^{α}$ =0，η $_{j}^{α}$ =1，C=C_α。本部分模拟中选取无量纲参数M=25，其余参数见表1。

图6是这个双晶粒体系演化的相场模拟。其中黄色部分表示晶粒，蓝绿色圆圈表示晶界。如图所示，随着时间的延长，圆形晶界不断朝向圆心移动，直径逐渐减小直到消失，2个晶粒最终完全“融合”成一个晶粒。图7是晶界移动速率和中心晶粒尺寸随时间演化的动力学曲线。可以看出，随着时间的延长，圆形晶粒的面积减小，晶界的曲率增加，晶界移动速率逐渐增大，这一结果与Hsueh等^[30]的理论研究基本一致。

图6

图6 双晶粒体系演化的相场模拟

Color online

(a) 1×10⁴ step (b) 3×10⁴ step (c) 4×10⁴ step (d) 5×10⁴ step

Fig.6 Phase-field simulation results of the microstructure evolution of two grain system

图7

图7 晶界移动速率和中心晶粒尺寸随时间演化的动力学曲线

Fig.7 Dynamic curves of grain boundary velocity and central grain size vs dimensionlesstimescale t_d

3.2 气孔与晶界的相互作用

本工作通过向模拟区域中添加气孔以研究气孔与晶界的相互作用，每个气孔的半径为5个格点。采用与3.1节相同的模拟区域和参数。在气孔相内部，η^β=1，C=C_β。

3.2.1 移动晶界与气孔的分离

图8是在晶界附近引入单个气孔后体系演化的相场模拟。其中黄色区域表示晶粒，绿色区域表示晶界，蓝色区域表示气孔。可以看出，气孔附近的晶界与气孔结合在一起，无法正常移动。气孔的引入会阻碍晶界的移动并导致体系演化速率减慢。图9展示了双气孔与晶界体系演化的相场模拟。可以看出，在演化过程中气孔和晶界发生了分离现象。

图8

图8 单气孔体系演化的相场模拟

Color online

(a) 1×10⁴ step (b) 3×10⁴ step (c) 4×10⁴ step (d) 5×10⁴ step

Fig.8 Phase-field simulation results of the microstructure evolution of one pore system

图9

图9 双气孔体系演化的相场模拟

Color online

(a) 1×10⁴ step (b) 5×10⁴ step (c) 7.5×10⁴ step (d) 10×10⁴ step

Fig.9 Phase-field simulation results of the microstructure evolution of two pore system

Ahmed等^[3]在对颗粒与晶界相互作用的模拟研究中，也发现了存在颗粒与晶界结合在一起运动和颗粒与晶界发生分离2种情况。进一步证明了相场模型中自由晶界(无气孔)的移动速率v_b和气孔移动速率v_p分别为^[31]:

v_{b} = - σ_{g b} M_{b} κ_{b}

(24)

v_{p} = \frac{σ_{i n t} D_{s} ν_{m} l_{s}}{R T} \nabla_{s}^{2} κ_{s}

(25)

式中，M_b为晶界迁移系数；κ_b为晶界曲率；l_s为气孔与晶粒接触表面的宽度； $\nabla_{s}^{2}$ 为气孔与晶粒接触表面的Laplace算子；κ_s为相界曲率。式(24)中等号右边的负号表示v_b的符号与κ_b的符号相反。在这里正负号的定义为：正号表示朝向界面的外(凸)侧，负号表示朝向界面的内(凹)侧。

从式(24)和(25)可以看出，当体系的物理参数σ_gb和M_b不变时，v_b与κ_b成正比；当气孔的形状(∇_s²κ_s)保持不变时，v_p与D_s成正比。一般情况下，v_b远大于v_p，因此气孔如同第二相钉扎颗粒^[3]一样，会严重阻碍晶界的运动。

结合这一理论对模拟结果进行分析：晶界运动是受曲率驱动的，κ_b为晶界运动提供了一个动力F_D，这个动力与κ_b成正比，即F_D∝κ_b。而气孔作为阻碍晶界运动的因素，对晶界施加了阻力，假设这个阻力的最大值为F_Zmax，当F_Zmax≥F_D，晶界与气孔钉扎在一起运动，当F_Zmax<F_D时，晶界与气孔便会分离^[9]。这便是气孔与晶界的分离条件。

利用分离条件可以解释图9中观察的现象：在演化初期，晶界半径较大，曲率较小，此时F_D<F_Zmax，晶粒和气孔结合在一起。随着时间的演化，气孔上下的晶界逐渐朝小晶粒的中心位置运动，这导致气孔附近的κ_b增大，F_D也随之增大，最终超过了F_Zmax=F_D这个临界，满足F_Zmax<F_D的条件，晶界气孔发生分离。虽然双气孔体系中晶粒和气孔发生了分离现象，但是相同模拟条件下该体系演化的速率仍然小于单气孔体系，这表明气孔数目的增加会降低晶粒的生长速率。

3.2.2 气孔对移动晶界的钉扎

图10是四气孔体系演化的相场模拟。可以看出，当体系演化时间步长为10×10⁴ step时，远离气孔的晶界曲率已经变为0，而靠近气孔的晶界曲率又不足以克服气孔的阻力，因此晶界被气孔牢牢钉扎，只能随着气孔一起缓慢移动(气孔并非不动，而是以非常缓慢的速率移动)，在这种情况下，体系的演化速率主要取决于气孔的移动速率，因此演化速率会显著减慢。

图10

图10 四气孔体系演化的相场模拟

Color online

(a) 1×10⁴ step (b) 2×10⁴ step (c) 10×10⁴ step (d) 25×10⁴ step

Fig.10 Phase-field simulation results of the microstructure evolution of four pore system

在气孔和晶界的相互作用过程中，气孔既可能与晶界一起运动，也可能与晶界分离，这取决于气孔附近晶界的曲率大小，若体系中靠近气孔的晶界曲率较大，足以克服气孔的最大阻力时，体系的演化将为晶界主导，发生气孔与晶界分离的现象；若体系中靠近气孔的晶界曲率很小，不足以克服气孔最大阻力时，体系的演化将由晶界主导变为气孔主导，而气孔的移动速率远远小于晶界的移动速率，因此体系演化速率显著减慢，不会发生气孔与晶界相互分离的现象。但上述2种情况中，气孔都会阻碍其附近晶界的运动。

3.3 含气孔多晶UO₂的晶粒生长

对2000 K下含气孔多晶UO₂体系的晶粒生长过程进行了模拟。模拟区域为256×256格点，区域内随机产生128个不同取向的晶粒。模拟中时间步长Δt，空间步长Δx和唯象参数A、D_α、D_β、γ_α、γ_β、δ_α、δ_β、ε见表1，动力学参数L=10，M=25，由于含气孔UO₂多晶体系2000 K下σ_int=0.6 J/m²，σ_gb=0.3 J/m^{2 [32]}，为了保证σ_int/σ_gb=2.0的条件，选取κ_α=κ_β=2，κ_c=8。

图11为气孔率f_p=2%时含气孔多晶UO₂体系演化的相场模拟。可以看出，由于含气孔附近的晶界受到气孔的阻碍，移动非常缓慢，而未与气孔发生作用的晶界移动相对较快。由于晶界的移动，小尺寸晶粒逐渐被大尺寸晶粒吞并，晶粒的总数目下降，晶粒平均面积增加。进一步分析发现，该模拟条件下未发现晶界与气孔相分离的现象，晶内的气孔对自由晶界的移动有很强的阻碍作用(如椭圆形区域内的气孔)。随着模拟时间的延长，很多原来处于晶界线周围的气孔移动到三叉晶界处，这主要是由于气孔的移动速率远远小于晶界的移动速率。

图11

图11 气孔率f_p=2%时含气孔多晶UO₂演化的相场模拟

Color online

(a) 0.1×10⁴ step (b) 0.5×10⁴ step (c) 1×10⁴ step (d) 2×10⁴ step

Fig.11 Phase-field simulation results of microstructure evolution of multi-grains UO₂ at porosity f_p=2% (The pores in the red elliptical areas hinder the movement of grain boundaries)

图12是不同气孔率下晶粒平均直径的变化规律。可以看出，相同演化时间下，晶粒的平均直径随气孔率的增大而下降，这表明气孔率的增加会减慢晶粒生长速率。这一结果与Kundin等^[12]的模拟结果类似。

图12

图12 不同气孔率下的晶粒生长曲线

Fig.12 Grain growth curves at different f_p

根据Rahaman^[1]提出的晶粒生长理论，多晶体系晶粒生长满足幂函数的规律：

D^{n} (t) - D_{0}^{n} = k t

(26)

式中，D为晶粒平均直径；D₀为初始晶粒平均直径；n为生长指数；k为比例系数。

表2是利用式(26)对图12中各曲线拟合得到的结果。可以看出，当气孔率由0增加到8%，n由2.05增加到4.19。Kundin等^[12]研究表明，随着气孔率的增加，n将由2变为4。这与本工作的研究结果基本一致。

表2 不同气孔率下的晶粒生长指数

Table 2 Grain growth exponents at different f_p

f_p	n	R²
0	2.05	0.990
2%	2.44	0.997
4%	3.15	0.992
8%	4.19	0.995

Note：n—grain growth exponent, R²—correlation coefficient

新窗口打开| 下载CSV

4 结论

(1) 本工作建立的相场模型能够成功模拟气孔与晶界相互作用的过程。通过对模型的分析发现：相场自由能密度中取向场梯度项前的参数对晶界和相界的界面宽度的影响较大，而浓度场梯度项前的参数对其影响很小；晶粒取向场梯度项前的参数对晶界界面能的影响较大，而浓度场梯度项前的参数对相界界面能的影响较大。

(2) 气孔和晶界相互作用的相场模拟结果表明：晶界的曲率是晶界移动的动力，晶界移动速率与晶界的曲率成正比；气孔会对移动的晶界施加一定的阻力，当气孔施加的最大阻力大于等于晶界移动的动力时，气孔会随晶界一起运动；而当气孔施加的最大阻力小于晶界移动的动力时，气孔与晶界分离。若气孔与晶界未发生分离，体系的演化将由晶界主导转变为气孔主导，演化速率显著下降。

(3) 含气孔UO₂多晶体系的晶粒生长的相场模拟结果表明：UO₂平均晶粒直径与时间成幂函数关系，幂指数随气孔率的增大而增大。气孔率越大，晶粒生长速率越慢。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

Rahaman

M N

. Ceramic Processing and Sintering [M]. New York: Marcel Dekker, 1995: 540

[本文引用: 3]

[2]

Turnbull

J A

The effect of grain size on the swelling and gas release properties of UO₂ during irradiation

[J]. J. Nucl. Mater., 1974, 50: 62

[3]

Ahmed

, Tonks

, Zhang

Y F

, et al.

Particle-grain boundary interactions: A phase field study

[J]. Comput. Mater. Sci., 2017, 134: 25

DOI URL [本文引用: 6]

[4]

Zhao

, Zhang

H Y

, Wei

, et al.

A phase field study for scaling rules of grain coarsening in polycrystalline system containing second-phase particles

[J]. Acta Metall. Sin., 2013, 49: 981

Ra and the time t from the non-linear relationship t=ARam+B. The kinetic exponent m also increased with the increasing volume fraction of second-phase particles. No matter whether the second-phase particles existed or not in the system studied, the scaling rule had been satisfied at the late stage of grain coarsening. The increase in the volume fraction of the second-phase particles would cause the decrease in the peak value of structure factor profile. When the value of the wave vector k increased to a certain value, the structure factor curve of the studied system was essentially coincident. With the increase in the volume fraction of second-phase particles, The peak values of scaling function decreased and the peak width became wider. According to structure factor and scaling function, it was known that with the increase in the volume fraction of second-phase particles, the interaction among grains weakens and the grain size would become more uniform during the grain coarsening.]]>

(赵彦, 张洪宇, 韦华等.

相场法研究含第二相颗粒多晶体系的晶粒粗化标度律

[J]. 金属学报, 2013, 49: 981)

. 结果表明:随着第二相颗粒体积分数增加, 晶粒生长阻力增大,晶粒平均半径Ra随时间t的演化关系偏离t=ARam+B非线性关系,且动力学指数m随第二相颗粒体积分数增加而增大.在晶粒粗化阶段,无论第二相颗粒是否存在, 及第二相颗粒体积分数的大小, 系统均满足标度律.随着第二相颗粒体积分数增加, 系统结构因子曲线的峰值逐渐降低,且当波矢k值增大到某一值后, 系统结构因子曲线基本重合;随着第二相颗粒体积分数增加, 系统标度函数峰值降低, 峰宽变宽. 在较大的k值情况下,同一k值对应的标度函数值随第二相颗粒体积分数的增加而增大.由结构因子和标度函数可知, 第二相颗粒体积分数增加, 晶粒间的相互作用减弱,晶粒粗化过程中的尺寸均匀性将更好. ]]>

[5]

Zhou

G Z

, Wang

Y X

, Chen

Phase-field method simulation of the effect of hard particles with different shapes on two-phase grain growth

[J]. Acta Metall. Sin., 2012, 48: 227

The effects of hard particles with different shapes, volume fractions and sizes on two-phase grain growth have been systematically investigated by phase-field method. The results showed that most of the spherical hard particles located at the intersection of tricrystal boundary, while flaky hard particles distributed along the grain boundary. Particles of different shapes have not obvious effect on the alpha phase grain growth, and the effect of hard particles with different shapes on the beta phase grain growth depends on the number of particles. The flaky particles have stronger pinning effect on the beta phase grain growth than the spherical particles when hard particles reach enough number. The pinning effect of the hard particles is enhanced when the volume fraction increased or the size of hard particles reduced. The greater the volume fraction or the smaller the size of hard particles is, the smaller the grains'size is.

(周广钊, 王永欣, 陈铮.

相场法模拟不同形状的硬质颗粒对两相晶粒长大的影响

[J]. 金属学报, 2012, 48: 227)

采用相场方法研究了不同颗粒体积分数及尺寸条件下不同形状的硬质颗粒对两相系统晶粒长大的影响, 结果表明: 球形颗粒大多处于三角晶界处, 片状颗粒处于晶界处且沿晶界分布. 不同形状的硬质颗粒对体积占优的α相晶粒长大无明显影响, 对体积分数较小的β相晶粒长大的影响主要取决于颗粒数目. 颗粒数目较少时, 不同形状的硬质颗粒对β相晶粒长大无明显影响; 颗粒数目较多时, 片状颗粒比球形颗粒对β相晶粒长大的阻碍作用强烈. 颗粒体积分数越大,颗粒对晶界的钉扎作用越强, 稳态时晶粒的半径越小; 颗粒尺寸越大,单个颗粒对晶界的钉扎作用越强, 但总的钉扎作用越弱, 稳态时晶粒的半径越大.

[6]

Gao

Y J

, Zhang

H L

, Jin

, et al.

Phase-field simulation of two-phase grain growth with hard particles

[J]. Acta Metall. Sin., 2009, 45: 1190

URL [本文引用: 1]

Grain growth, due to its importance in controlling the physical properties of a wide variety of materials, has been extensively investigated. Second–phase particles have the capacity to "pin" grain boundaries and therefore affect the grain growth behavior of polycrystalline materials profoundly. They reduce the mobility of grain boundaries and eventually, when a critical grain size is reached, arrest grain growth. Based on a diffuse–interface description, a computer simulation model for studying the microstructural evolution in two–phase solid has been developed. For a grain system with hard particles, the kinetics of two–phase grain growth with the third hard particles was investigated by phase field model with a continuum diffuse–interface field. A polycrystalline microstructure of temporal and spatial evolution of the three–phase–solid system was obtained by solving three kinetics equations. It is found that the pinning effect is enhanced with the increase of the size and the volume fraction of third–phase particles. The greater the volume fraction and size of third–phase particles are, the smaller the limited sizes of grain growth are. If the volume fraction of third–phase particle maintains a constant and the size of third–phase particles is smaller, then the pinning effect of third–phase particles is stronger. When third particles with two different sizes under the same volume fraction are introduced in the system of grain growth, the pinning effect of the particles is the best. The power growth law, grain morphology, critical grain size, grain growth dynamics and topology structure of two–phase polycrystalline materials simuated by phase–fielmodel are in well accordnce with the experimental results and theoretical results of other simulations.

(高英俊, 张海林, 金星等.

相场方法研究硬质颗粒钉扎的两相晶粒长大过程

[J]. 金属学报, 2009, 45: 1190)

URL [本文引用: 1]

采用相场方法模拟第三相颗粒钉扎的两相耦合的晶粒长大过程, 系统地研究了第三相颗粒体积分数和尺寸大小对两相晶粒长大过程的影响. 模拟结果表明, 第三相颗粒体积分数越大, 对晶界的钉扎作用越强, 且极限晶粒尺寸越小. 单个第三相颗粒尺寸越大, 对晶界钉扎作用越强. 但当体积分数一定时第三相颗粒尺寸越小时, 颗粒数目会越多, 此时总的钉扎效果会越好, 晶粒极限尺寸也越小. 若晶粒长大系统同时引入两种不同大小的第三相钉扎颗粒, 且两种颗粒所占比例相同时, 钉扎效果最好. 相场方法模拟所得到的二相多晶材料晶粒组织演化规律和晶粒生长指数、晶粒形态、生长动力学和拓扑结构特征与已有实验和理论结果相符合

[7]

Ahmed

, Yablinsky

C A

, Schulte

, et al.

Phase field modeling of the effect of porosity on grain growth kinetics in polycrystalline ceramics

[J]. Model. Simul. Mater. Sci. Eng., 2013, 21: 065005

[8]

Moelans

, Blanpain

, Wollants

Pinning effect of second-phase particles on grain growth in polycrystalline films studied by 3-D phase field simulations

[J]. Acta Mater., 2007, 55: 2173

DOI URL [本文引用: 2]

AbstractThree-dimensional simulations of grain growth in thin films containing finely dispersed second-phase particles were performed using a phase field model. The simulations show that although the growth behavior of the columnar grain structures in thin films is essentially two-dimensional, the interaction between the particles and the grain boundaries is three-dimensional. Grain boundaries can therefore more easily break free from the particles than in purely two-dimensional systems, resulting in fewer grain boundary–particle intersections and a larger final grain size. For a given volume fraction fV and size of the particles r, the final grain size increases with film thickness. Moreover, it was found that particles located in the middle of the film are most efficient in pinning grain boundaries. A classical Zener type relation cannot describe these effects.]]>

[9]

Moelans

, Blanpain

, Wollants

A phase field model for the simulation of grain growth in materials containing finely dispersed incoherent second-phase particles

[J]. Acta Mater., 2005, 53: 1771

DOI URL [本文引用: 2]

AbstractA phase field model is presented for simulating grain growth in materials containing small incoherent second-phase particles that are constant in time. The model of [L.-Q. Chen, W. Yang, Phys. Rev. B 50 (1994) 15752] for grain growth in single phase materials, that represents a polycrystalline microstructure with a set of phase field variables ηi(r,t), has been modified to incorporate the presence of second-phase particles by adding an extra term to the free energy expression, with Φ a spatially dependent parameter that equals one inside the particles, and zero elsewhere and ? a positive constant. The interaction between a particle and a diffuse interphase is analyzed from an energetic point of view and the effect of pinning particles on the shape and the evolution of an initially spherical grain is studied for two-dimensional and three-dimensional systems.]]>

[10]

Chang

, Feng

W M

, Chen

L Q

Effect of second-phase particle morphology on grain growth kinetics

[J]. Acta Mater., 2009, 57: 5229

AbstractSecond-phase particles are often employed to inhibit grain growth in polycrystalline metals and ceramics. In this work, we studied the effect of second-phase particle morphology on the effectiveness of inhibiting grain boundary migration using the phase-field method. We employed a multi-order parameter phase-field model in combination with an efficient memory allocation strategy which allows large-scale and coalescence-free grain growth simulations. We analyzed the dependence of pinning forces on the particle size and shape, and performed computer simulations of grain growth in the presence of second-phase particles with different sizes and varying aspect ratios. We also discuss the relationship between the pinned grain size and size distributions.]]>

[11]

Hötzer

, Rehn

, Rheinheimer

, et al.

Phase-field study of pore-grain boundary interaction

[J]. J. Ceram. Soc. Jpn., 2016, 124: 329

[12]

Kundin

, Sohaib

, Schiedung

, et al.

Phase-field modelling of pores and precipitates in polycrystalline systems

[J]. Model. Simul. Mater. Sci. Eng., 2018, 26: 065003

DOI URL [本文引用: 4]

[13]

Nichols

F A

Theory of grain growth in porous compacts

[J]. J. Appl. Phys., 1966, 37: 4599

[14]

Brook

R J

Pore-grain boundary interactions and grain growth

[J]. J. Am. Ceram. Soc., 1969, 52: 56

[15]

Riedel

, Svoboda

A theoretical study of grain growth in porous solids during sintering

[J]. Acta Metall. Mater., 1993, 41: 1929

[16]

Petrishcheva

, Renner

Two-dimensional analysis of pore drag and drop

[J]. Acta Mater., 2005, 53: 2793

AbstractWe consider stationary coupled migration (drag) and conditions for separation (drop) of pores situated on grain boundaries or triple junctions in two-dimensional (2D) geometry. Pore mobility is realized by surface diffusion and boundaries migrate owing to surface tension. A small velocity approximation yields pore mobility and estimates for the velocity at which separation occurs. The estimate is refined by numerical solution. In contrast to previous expectations, separation occurs in 2D and the critical velocity for pores with circular cross section is non-zero. Our numerics are further confirmed by bifurcation analysis. The critical velocities for 2D pores at triple junctions are considerably smaller than for 2D pores on grain boundaries, which in turn are smaller by a factor of two to three compared to critical velocities of three-dimensional lenticular pores on grain boundaries. We also determine the coupled pore-boundary mobility and demonstrate that the boundary mobility is practically always reduced by pores.]]>

[17]

Klinger

, Rabkin

, Shvindlerman

L S

, et al.

Grain growth in porous two-dimensional nanocrystalline materials

[J]. J. Mater. Sci., 2008, 43: 5068

[18]

Fan

, Chen

L Q

Computer simulation of grain growth using a continuum field model

[J]. Acta Mater., 1997, 45: 611

DOI URL [本文引用: 3]

[19]

Moelans

, Blanpain

, Wollants

An introduction to phase-field modeling of microstructure evolution

[J]. Calphad, 2008, 32: 268

AbstractThe phase-field method has become an important and extremely versatile technique for simulating microstructure evolution at the mesoscale. Thanks to the diffuse-interface approach, it allows us to study the evolution of arbitrary complex grain morphologies without any presumption on their shape or mutual distribution. It is also straightforward to account for different thermodynamic driving forces for microstructure evolution, such as bulk and interfacial energy, elastic energy and electric or magnetic energy, and the effect of different transport processes, such as mass diffusion, heat conduction and convection. The purpose of the paper is to give an introduction to the phase-field modeling technique. The concept of diffuse interfaces, the phase-field variables, the thermodynamic driving force for microstructure evolution and the kinetic phase-field equations are introduced. Furthermore, common techniques for parameter determination and numerical solution of the equations are discussed. To show the variety in phase-field models, different model formulations are exploited, depending on which is most common or most illustrative.]]>

[20]

Krivilyov

M D

, Mesarovic

S D

, Sekulic

D P

Phase-field model of interface migration and powder consolidation in additive manufacturing of metals

[J]. J. Mater. Sci., 2017, 52: 4155

[21]

Lin

M Z

, Zhang

R J

, Fang

, et al.

Phase field simulation of sintering process in biphasic porous material

[J]. Acta Metall. Sin., 2012, 48: 1207

Sintering is a process of bonding between solid particles which typically occurs under high temperature. Currently, simulation of sintering process is mainly concentrated on the single-phase polycrystalline materials. As there are a lot of materials which are biphasic porous system, it is of practical significance to simulate the microstructural evolution of biphasic porous system during sintering process. In this work, a new phase field model is established to simulate sintering process in biphasic porous system. The evolution of the component is governed by Cahn-Hilliard equation, while the orientation field by the time-dependent Allen-Calm equation. A. function is established to describe the relationship between atomic diffusion coefficient and grain boundary diffusion, surface diffusion and volume diffusion. A group of phenomenological coefficients are obtained by analyzing the characteristic of the phase-field model. The simulation results show that the new phase-field model can effectively simulate the sintering process in biphasic porous system. The formation and growth of sintering neck, the seal spheroidization and disappearance of pores as well as the mergence and growth of grains are observed during simulation. The sintering necks between the parent phase and the second phase grow very fast at the early stage of simulation, while at the late stage, because of the pinning effect, the growth rate of the sintering neck slows down obviously, pores become isolated by the grains, and its shape change from concave to convex, the relative small pores are eliminated, which leads to densification. As the sintering proceeds, the grain size of the second phase gradually decreases and the parent-phase grains are wrapped by the second-phase grains Because of the pinning effect of the second phase, the migration rate of the grain boundary of the parent phase is restrained. The evolution course of pores depends largely on the interaction between the second phase and the pores. The evolution rate of pores is quantitatively compared between the biphasic porous system and the single-phase system. In the case of biphasic porous system, the evolution rate of pores is slower than that in single-phase system. The simulating growth exponents of the parent phase are calculated with different volume fractions of the second phase. As the volume fractions of the second phase increase from 15% to 25%, the grain growth exponent changes from 2.9 to 3.4.

(刘明治, 张瑞杰, 方伟等.

相场法模拟两相多孔组织烧结

[J]. 金属学报, 2012, 48: 1207)

建立了新的模拟两相多孔材料烧结过程的相场模型, 采用Cahn-Hillard方程和Allen-Cahn方程来控制相对密度场和长程取向场的变化, 通过分析相场方程的特点, 对模型进行数学处理得到一组相场模型的唯象系数, 建立了原子扩散系数与晶界扩散、表面扩散和体积扩散的函数关系式. 模拟结果表明: 该模型能够有效地模拟两相多孔材料的烧结过程, 通过分析模拟图像可以很好地观察到两相多孔材料的微观组织演化过程.

[22]

Chen

, Kang

X H

, Xiao

N M

, et al.

Phase field modelling of grain growth in polycrystalline material

[J]. Acta Phys. Sin., 2009, 58: 124

(陈云, 康秀红, 肖纳敏等.

多晶材料晶粒生长粗化过程的相场方法模拟

[J]. 物理学报, 2009, 58: 124)

ψ和晶体学取向θ来表示多晶粒结构的相场模型,利用自适应有限元方法模拟了多晶材料等温过程中的晶粒粗化现象.模拟结果显示,在曲率作用下,通过晶界迁移弯曲晶界逐渐平直化,小晶粒逐渐被大晶粒吞并,当晶界之间的取向差较小时,满足一定能量和几何条件的两晶粒在界面能作用下会发生转动,合并为单个晶粒.模拟结果与实验结果符合较好.因此,该相场模型可以很好地用来模拟固态相变中多晶材料的生长粗化等现象.]]>

[23]

Jing

X N

, Ni

, He

L H

, et al.

2-D phase-field simulation of pore evolution in sintering ceramics

[J]. J. Inorg. Mater., 2002, 17: 1078

(景晓宁, 倪勇, 何陵辉等.

陶瓷烧结过程孔隙演化的二维相场模拟

[J]. 无机材料学报, 2002, 17: 1078)

URL [本文引用: 1]

采用二维相场模型模拟陶瓷烧结过程中颗粒间孔隙的演化过程.选取四方堆积颗粒间气孔作为对象.通过连续的密度场和长程取向场（LRO）描述烧结体的微结构,密度场的演化由Cahn-Hillard（CH）方程控制,而颗粒的取向场演化由时间相关的Ginzburg-Laudau（TDGL）方程控制.上述非线性演化方程利用半隐傅立叶频域法求解.模拟结果反映了颗粒间接触,烧结颈生长和气孔球化的微观过程.量化计算烧结颈生长率以及在不同晶界和表面迁移率比值时的烧结率,较好地符合理论分析的趋势.

[24]

Deng

A phase field model of sintering with direction-dependent diffusion

[J]. Mater. Trans., 2012, 53: 385

[25]

Sun

Z Y

, Yang

, Liu

W B

Phase-field simulations of the sintering process of uranium dioxide

[J]. Acta Metall. Sin., 2020, 56: 1295

(孙正阳, 杨超, 柳文波.

UO₂烧结过程的相场模拟

[J]. 金属学报, 2020, 56: 1295)

[26]

, Jiang

, Xu

Z Z

. Solid State Physics [M]. 2nd Ed., Shanghai: Shanghai Scientific and Technical Publishers, 2010: 364

[本文引用: 2]

(陆栋, 蒋平, 徐至中. 固体物理学 [M]. 第2版. 上海: 上海科学技术出版社, 2010: 364)

[本文引用: 2]

[27]

Chen

L Q

, Fan

D N

Computer simulation model for coupled grain growth and Ostwald ripening—Application to Al₂O₃-ZrO₂ two-phase systems

[J]. J. Am. Ceram. Soc., 1996, 79: 1163

DOI URL [本文引用: 3]

[28]

Cahn

J W

On spinodal decomposition

[J]. Acta Metall., 1961, 9: 795

[29]

Moelans

, Blanpain

, Wollants

Quantitative analysis of grain boundary properties in a generalized phase field model for grain growth in anisotropic systems

[J]. Phys. Rev., 2008, 78B: 024113

[30]

Hsueh

C H

, Evans

A G

, Coble

R L

Microstructure development during final/intermediate stage sintering—I. Pore/grain boundary separation

[J]. Acta Metall., 1982, 30: 1269

[31]

Ahmed

, Allen

, El-Azab

Phase field modeling for grain growth in porous solids

[J]. J. Mater. Sci., 2016, 51: 1261

[32]

Ahmed

, Pakarinen

, Allen

, et al.

Phase field simulation of grain growth in porous uranium dioxide

[J]. J. Nucl. Mater., 2014, 446: 90