黏性烧结过程的相场模拟

图1 半交错网格示意图

Fig.1 Schematic illustration of staggered mesh (u is the velocity component in the x direction, and v is the velocity component in the y direction. The red points are used for the discrete format of u, and the blue points are used for the discrete format of v. The black points are used for the discrete format of scalars, including viscosity μ, density ρ and phase field c)

{\tilde{u}}_{(I, J)} = \frac{{\tilde{u}}_{(i + 1, J)} + {\tilde{u}}_{(i, J)}}{2}

(19)

一号坐标与三号坐标系之间的映射同理。二号坐标系与三号坐标系的映射公式如下所示：

{\tilde{v}}_{(i, J)} = \frac{{\tilde{v}}_{(I, j - 1)} + {\tilde{v}}_{(I + 1, j - 1)} + {\tilde{v}}_{(I, j)} + {\tilde{v}}_{(I + 1, j)}}{4}

(20)

不同坐标系下的映射完成后对方程进行离散处理，下文介绍在半交错网格下对Cahn-Hilliard方程和Navier-Stokes方程的离散方式。

2.3.1 Navier-Stokes方程差分格式

首先对Navier-Stokes方程中的对流项进行差分，对流项为等式(14)左边第二项 $\tilde{ρ} \tilde{V} \tilde{\nabla} \cdot \tilde{V}$ 。由于该项为非线性项，对微分方程的数学性质影响较大，因此它的离散格式主要决定数值求解的精确性和稳定性。该项差分格式采用迎风格式，演化中由于每点的速度方向不确定，根据运动速度方向的不同选择相应的差分表达，将其水平竖直方向进行分解，x方向表示为 $(\tilde{u} {\tilde{u}}_{x}^{'} + \tilde{v} {\tilde{u}}_{y}^{'})$ ，y方向表示为 $(\tilde{u} {\tilde{v}}_{x}^{'} + \tilde{v} {\tilde{v}}_{y}^{'})$ ，代入坐标系之间的映射关系有：

(\tilde{u} {\tilde{u}}_{x}^{'} + \tilde{v} {\tilde{u}}_{y}^{'})_{(i, J)} = {\tilde{u}}_{(i, J)} {\tilde{u}}_{x (i, J)}^{'} + {\tilde{v}}_{(i, J)} {\tilde{u}}_{y (i, J)}^{'} = {\tilde{u}}_{(i, J)} {\tilde{u}}_{x (i, J)}^{'} + \frac{{\tilde{v}}_{(I, j - 1)} + {\tilde{v}}_{(I + 1, j - 1)} + {\tilde{v}}_{(I, j)} + {\tilde{v}}_{(I + 1, j)}}{4} {\tilde{u}}_{y (i, J)}^{'}

(21)

(\tilde{u} {\tilde{v}}_{x}^{'} + \tilde{v} {\tilde{v}}_{y}^{'})_{(I, j)} = {\tilde{v}}_{(I, j)} {\tilde{v}}_{y (I, j)}^{'} + {\tilde{u}}_{(I, j)} {\tilde{v}}_{x (I, j)}^{'} = {\tilde{v}}_{(I, j)} {\tilde{v}}_{y (I, j)}^{'} + \frac{{\tilde{u}}_{(i - 1, J)} + {\tilde{u}}_{(i - 1, J + 1)} + {\tilde{u}}_{(i, J)} + {\tilde{u}}_{(i, J + 1)}}{4} {\tilde{v}}_{x (I, j)}^{'}

(22)

式中， ${\tilde{u}}_{x}^{'}$ 、 ${\tilde{u}}_{y}^{'}$ 分别表示u在x、y方向的导数， ${\tilde{v}}_{x}^{'}$ 、 ${\tilde{v}}_{y}^{'}$ 分别表示v在x、y方向的导数，根据u₍_i_,_J₎、v₍_i_,_J₎、u₍_I_,_j₎、v₍_I_,_j₎的正负选用相应方向的迎风格式进行差分化。对于其中偏导数，根据速度方向进行迎风方向的一阶差分即：

{\tilde{u}}_{x (i, J)}^{'} = \{\begin{matrix} \frac{{\tilde{u}}_{(i, J)} - {\tilde{u}}_{(i - 1, J)}}{2} ({\tilde{u}}_{(i, J)} > 0) \\ \frac{{\tilde{u}}_{(i + 1, J)} - {\tilde{u}}_{(i, J)}}{2} (o t h e r w i s e) \end{matrix}

(23)

{\tilde{u}}_{y (i, J)}^{'} = \{\begin{matrix} \frac{{\tilde{u}}_{(i, J)} - {\tilde{u}}_{(i, J - 1)}}{2} \\ ({\tilde{v}}_{(I, j - 1)} + {\tilde{v}}_{(I + 1, j - 1)} + {\tilde{v}}_{(I, j)} + {\tilde{v}}_{(I + 1, j)} > 0) \\ \frac{{\tilde{u}}_{(i, J + 1)} - {\tilde{u}}_{(i, J)}}{2} (o t h e r w i s e) \end{matrix}

(24)

将Navier-Stokes方程中的扩散项差分化，速度Laplace算子使用五点差分法求解。半交错网格中，速度点上的动力黏性系数 $\tilde{μ} (c)$ 用中心格式进行求解。

2.3.2 Cahn-Hilliard方程差分格式

Cahn-Hilliard方程中自由能变分项 $\tilde{\nabla} (\tilde{M} \tilde{\nabla} \frac{δ \tilde{F}}{δ c})$ 在空间的离散类似于Navier-Stokes方程中的黏性项增加方程数学上的稳定性，可以直接使用中心差分法进行计算。Cahn-Hilliard方程中引入速度场增加对流项 $\tilde{\nabla} \cdot (\tilde{V} c)$ 对数值模拟稳定性影响较大。Cahn-Hilliard方程数学形式为椭圆形方程。因此对流项也应使用迎风差分格式进行离散，即：

\tilde{\nabla} \cdot (\tilde{V} c) = \frac{\partial (\frac{{\tilde{u}}_{(i + 1, J)} + {\tilde{u}}_{(i, J)}}{2} c_{(I, J)})}{\partial \tilde{x}} + \frac{\partial (\frac{{\tilde{v}}_{(I, j + 1)} + {\tilde{v}}_{(I, j)}}{2} c_{(I, J)})}{\partial \tilde{y}}

(25)

式(25)中， $若 \frac{{\tilde{u}}_{(i + 1, J)} + {\tilde{u}}_{(i, J)}}{2} \geq 0$ ，

\frac{\partial (\frac{{\tilde{u}}_{(i + 1, J)} + {\tilde{u}}_{(i, J)}}{2} c_{(I, J)})}{\partial \tilde{x}} = \frac{{\tilde{u}}_{(i + 1, J)} + {\tilde{u}}_{(i, J)}}{2} c_{(I, J)} - \frac{{\tilde{u}}_{(i, J)} + {\tilde{u}}_{(i - 1, J)}}{2} c_{(I - 1, J)}

(26)

若 \frac{{\tilde{u}}_{(i + 1, J)} + {\tilde{u}}_{(i, J)}}{2} < 0

，

\frac{\partial (\frac{{\tilde{u}}_{(i + 1, J)} + {\tilde{u}}_{(i, J)}}{2} c_{(I, J)})}{\partial x} = \frac{{\tilde{u}}_{(i + 2, J)} + {\tilde{u}}_{(i + 1, J)}}{2} c_{(I + 1, J)} - \frac{{\tilde{u}}_{(i + 1, J)} + {\tilde{u}}_{(i, J)}}{2} c_{(I, J)}

(27)

若 \frac{{\tilde{v}}_{(I, j + 1)} + {\tilde{v}}_{(I, j)}}{2} \geq 0

，

\frac{\partial (\frac{{\tilde{v}}_{(I, j + 1)} + {\tilde{v}}_{(I, j)}}{2} c_{(I, J)})}{\partial \tilde{y}} = \frac{{\tilde{v}}_{(I, j + 1)} + {\tilde{v}}_{(I, j)}}{2} c_{(I, J)} - \frac{{\tilde{v}}_{(I, j)} + {\tilde{v}}_{(I, j - 1)}}{2} c_{(I, J - 1)}

(28)

若 \frac{{\tilde{v}}_{(I, j + 1)} + {\tilde{v}}_{(I, j)}}{2} < 0

，

\frac{\partial (\frac{{\tilde{v}}_{(I, j + 1)} + {\tilde{v}}_{(I, j)}}{2} c_{(I, J)})}{\partial \tilde{y}} = \frac{{\tilde{v}}_{(I, j + 2)} + {\tilde{v}}_{(I, j + 1)}}{2} c_{(I, J + 1)} - \frac{{\tilde{v}}_{(I, j + 1)} + {\tilde{v}}_{(I, j)}}{2} c_{(I, J)}

(29)

其余项中Laplace算子采用五点差分法，梯度和散度使用中心差分法进行计算。

2.3.3 投影法求解压强

根据时间离散与空间离散求解出速度与相场变量预测步的值并非精确解，上面的求解并未考虑式(13)对速度的限制，求解出预测步的速度无法满足式(13)散度为0的条件，因此需要考虑压强对速度进行修正，求解出满足式(13)限制条件的速度场。

不可压缩流体的流场求解压强并没有直接的计算表达形式，只能通过动量守恒方程和质量守恒方程联立求得压强场的形式。本工作使用投影法，将u和v分别进行求解，暂不考虑压强项的影响得出预测步的x方向速度 $u^{*}$ ：

\tilde{ρ} (\frac{{\tilde{u}}^{*} - {\tilde{u}}^{n}}{d \tilde{t}} + (\tilde{u} {\tilde{u}}_{x} + \tilde{v} {\tilde{u}}_{y})^{n}) = \tilde{\nabla} \cdot (\tilde{μ} (c) (\tilde{\nabla} \tilde{u})) + {\tilde{s}}_{x}^{n}^{}

(30)

式中， ${\tilde{s}}_{x}^{}$ 为无量纲冲量源项 $\tilde{s}$ 的x方向分量。将离散的动量方程代入连续性方程得到压力方程，二维有限差分格式下经过计算整理得出压强的Poisson方程：

- \tilde{\nabla} \cdot (\frac{\tilde{\nabla} {\tilde{p}}^{n + 1}}{{\tilde{ρ}}^{n}}) = - \frac{1}{Δ \tilde{t}} \tilde{\nabla} \cdot {\tilde{u}}^{*}

(31)

压强求解后，对预测步x方向的u进行修正：

{\tilde{ρ}}^{n} (\frac{1}{Δ \tilde{t}} \tilde{\nabla} \cdot {\tilde{u}}^{n + 1} - \frac{1}{Δ \tilde{t}} \tilde{\nabla} \cdot {\tilde{u}}^{*}) = - \tilde{\nabla} \cdot {\tilde{p}}^{n + 1}

(32)

对压强与速度的值进行迭代计算，直到满足散度为0的条件即可输出，对y方向同理可以得到v的值。

3 结果与讨论

3.1 烧结过程的形貌演化

利用构造的黏性烧结模型对等大双颗粒的烧结过程进行模拟。模拟区域差分网格尺寸为256 × 256，采用周期性边界条件，模拟颗粒的形貌、速度以及压强在整个区域的演化过程。等大双圆形颗粒烧结过程的形貌演变如图2所示。可以看出，随着模拟时间的延长，由最初等大两圆形颗粒逐渐演化为表面能更低的单个圆形颗粒。黏性烧结演化过程的驱动力主要为曲率梯度驱动。曲率梯度驱动是由烧结颈处曲率和液体颗粒曲率的差值引起的传质。最初阶段为液体颗粒黏弹性接触，此时界面间附着力使其进一步演化。接下来在液体颗粒接触点形成凹凸结合的马鞍状烧结颈，在晶体颗粒曲率梯度的作用下液体原子从曲率较小的液体边界处移动到曲率较大的烧结颈处，原子的沉积和迁移使得烧结颈半径增大。随着烧结的进行，烧结颈的曲率逐渐减小，曲率梯度逐渐减小，驱动力减小导致物质通量下降，反应过程变得缓慢。该相场模拟结果与文献[25]黏弹性接触结果类似。

图2

图2 等大双圆形颗粒烧结过程的形貌演化

(a) 0 step (b) 2 × 10⁴ step (c) 2 × 10⁵ step (d) 1.2 × 10⁶ step (e) 6 × 10⁶ step (f) 4 × 10⁷ step

Fig.2 Microstructure evolutions of the two equal circle grains during sintering

3.2 颈球比与收缩率的变化

图3和4分别为不同迁移率下的颈球比和收缩率随演化时间步的变化关系。其中，T表示时间步；颈球比为X / D (其中，X为t时刻烧结颈的长度，D为t时刻的粒径)；收缩率为(2D- L) / D (其中，L为双颗粒演化过程中的长度)。可以看出，颈球比和收缩率变化趋势一致，烧结开始阶段变化较快，随着时间的延长变得缓慢。由于开始时不平衡，驱动力较大，在2 × 10⁷时间步后，变化速度明显变慢。随着迁移率的增大，整体演化速率加快，整体演化趋势不变，且最终稳定的值并不发生改变。进一步分析发现，模拟后期的颈球比为1.411 (接近 $\sqrt[]{2}$ )，而收缩率为0.588 (接近2 - $\sqrt[]{2}$ ）。因此本工作提出的模型满足质量守恒以及相场变量守恒 $\frac{d}{d t} \int c d Ω = 0$ 。Frenkel基于原子扩散和黏性流动的类比，提出的烧结模型的模拟结果颗粒半径不发生变化^[7]；且颈球比随时间的变化为 ${(\frac{X}{D})}^{2} \propto t$ 。本工作中，颗粒演化过程的模拟结果，接近黏弹性接触(viscoelastic contact)，与文献[26]中颈球比随时间变化函数相似。

图3

图3 不同迁移率下颈球比随时间的变化关系

Fig.3 Normalized contact radius growth curves under different mobilities ( $\tilde{M}$ —mobility, T—time step)

图4

图4 不同迁移率下收缩率随时间的变化关系

Fig.4 Normalized shrinkage vs time step curves under different mobilities

3.3 烧结过程速度场与压强场的变化

图5为烧结过程中沿水平方向(即x方向)的速度分量u的分布。可以看出，u非零的部分主要分布在液体内部和气/液边界处。在黏性烧结模拟过程中，液相颗粒内部的速度形成主要是由相的迁移引起的。烧结过程中烧结颈处，曲率梯度最大，导致该处质量迁移最大，因此表现为图5中颜色最深处，水平方向速率最大。速度演化过程中，速度开始在气液接触界面产生，逐渐传递到液体颗粒内部。双颗粒中液体分子水平速度方向均沿着颗粒合并的方向。

图5

图5 颗粒形貌演变与对应的水平方向的速度场分布

Fig.5 Microstructure evolutions (a-c) and the corresponding horizontal velocity distributions (d-f) of the two equal circle grains (Blue represents negative values indicating the direction of velocity to the left; yellow is a positive value representing the right direction of velocity, ranging from -3.3 × 10^-3 m/s to 3.3 × 10^-3 m/s in the velocity cloud map) (a, d) 6 × 10⁵ step (b, e) 2 × 10⁶ step (c, f) 8 × 10⁶ step

图6为沿水平方向对称轴上的相场变量和轴向速度场的变化曲线，其中水平方向速度为正，则速度方向水平向右；反之则水平向左。可以看出，在液体颗粒接触界面上存在一个滞止点(stagnation point)，该点水平方向的速率为0且在该点附近水平方向速率与距滞止点的距离成线性关系，线性变化率与轴向应变率有关且与接触面的距离成正比，该区域称为速度的接触应变区。接触应变区随着烧结颈的增大而增大。远离接触界面的区域称为刚体运动区域，距离接触面较远，因此演化过程中烧结颈变化引起的应力变化对该部分影响较小，其速度变化为单纯的黏性流体流动所引发的，速度变化较慢。沿远离对称轴的方向，液相区转变为气相区，此时轴向速度很小且迅速减为0。轴向速度模拟结果与Wakai等^[27]的理论分析结果一致。

图6

图6 沿水平方向对称轴上的相场变量和轴向速度场的变化曲线

Fig.6 Distributions of phase-field variable and velocity field along horizontal symmetry axis of the two circle grains

图7为烧结过程中压强在空间的分布，可以看出由压强差造成液体颗粒的合并与演化过程。在液相颗粒的两边压强较大，中间压强较小。在液体颗粒接触面上压强分布并不均匀，而是在接触面上根据烧结颈曲率的大小进行变化。烧结颈为内凹的马鞍状，因此其曲率大于0，导致压强小于周围环境压强。伴随着时间步长的推进，烧结颈处曲率逐渐减小，烧结颈部压强与外界压差也逐渐减小，最终演化为一个液体球形颗粒时曲率变为负，此时液体颗粒内部也不存在负压强点。演化结束后，形成单个液相颗粒，内部的压强均匀分布，与文献[27]压强模拟结果类似。

图7

图7 等大双圆形颗粒烧结过程的压强分布

Fig.7 Pressure (P) distributions of two equal circle grains

(a) 2 × 10⁴ step (b) 1 × 10⁵ step (c) 2 × 10⁵ step (d) 1.2 × 10⁶ step (e) 1 × 10⁷ step (f) 4 × 10⁷ step

3.4 多颗粒烧结过程模拟

设定三颗粒夹角为60°，其形貌演化过程如图8所示。由于颗粒接触的限制，三颗粒合并过程中，并没有产生过大的位移，在表面能减少驱动作用下，三颗粒最终形成单个表面能最小的圆形颗粒。可以看到，由于三颗粒之间孔隙较小，很快发生球化，液相在演化过程中填充孔隙，液相颗粒致密化。在烧结过程中，由于小孔隙的表面能较高，因此优先被填充。该处孔隙存在浮力作用，孔隙会缓慢向上迁移到达烧结顶部，导致气孔消失。模拟结果与文献[26]中提出颗粒黏性烧结4个阶段相同：自由流动、接触、液桥增长、气孔消失。

图8

图8 等大三颗粒烧结过程的形貌演化

Fig.8 Microstructure evolutions of three particles with the initial contact angles 60°

(a) 0 step (b) 6 × 10⁴ step (c) 2 × 10⁵ step (d) 1 × 10⁶ step (e) 6 × 10⁶ step (f) 2 × 10⁷ step

多颗粒烧结过程的形貌演变图如图9所示。在模拟的开始阶段，选择16个尺寸和分布均随机圆形颗粒作为初始组织，颗粒之间的孔隙形状各不相同。演化过程中，在黏性烧结的初始阶段烧结颈形成，接触点扩大。中间阶段，图中蓝色圆圈中的孔隙，在该阶段先进行球化，之后被液相填充。中间阶段的颈部生长较快，通过相邻颈部的融合使得孔隙进行变化，较小孔隙迅速被液相填充，较大孔隙则发生球化现象。中间阶段结束时，孔隙闭合圆状，总孔隙率下降，此时较大孔隙并未形成完美的圆形，即红圈中的孔隙在T = 1 × 10⁵ step时，边界变得圆润但并未形成完美的圆形。最终阶段，孔隙闭合致密度增加，导致大孔隙被挤压成封闭的圆形，此时孔隙数量减少、平均尺寸增加。

图9

图9 多颗粒烧结过程的形貌演化

Fig.9 Microstructure evolutions of multi-particles (The red circle shows the region of pore spheroidization, and the blue circle shows the region of pore disappearance)

(a) 0 step (b) 2 × 10⁴ step (c) 1 × 10⁵ step (d) 2 × 10⁶ step (e) 2.4 × 10⁶ step (f) 3 × 10⁶ step

4 结论

(1) 建立了黏性烧结相场模型，通过有限差分法和投影法求解了Cahn-Hilliard方程与Navier-Stokes方程，并对等大双圆形颗粒的烧结过程进行模拟。形貌演化结果为黏弹性接触满足质量守恒；颈球比与收缩率在模拟的开始阶段变化迅速，随着演化时间的延长驱动力下降，缓慢趋近于稳定值；增大M可以加快演化的形成，但对稳定状态的形貌影响不大。

(2) 对速度场和压强场的模拟表明，颗粒内部速度分为接触应变区和刚体运动区，接触应变区速度变化为纯应变运动，轴向速率与接触面距离成正比，刚体运动区为单纯黏性流动，速度变化较为缓慢；颗粒界面附近压强与外界压强差与该处曲率成正比，颗粒内部压强均匀分布。

(3) 多颗粒黏性烧结经历了接触点增大、烧结颈形成、相邻颈部融合等微观过程；孔隙尺寸对球化的影响很大，尺寸较小的孔隙快速发生球化并被液体填充，而尺寸较大的孔隙发生球化的速度较慢。

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

German

R M

, translated by Jia

C C

, Chu

, Liu

B W

, et al. Sintering: From Empirical Observations to Scientific Principles [M]. Beijing: Chemical Industry Press, 2021: 5

兰德尔·杰曼著, 贾成厂, 褚克, 刘博文

等译. 烧结实践与科学原理 [M]. 北京: 化学工业出版社, 2021: 5

[2]

German

R M

, Suri

, Park

S J

Review: Liquid phase sintering

[J]. J. Mater. Sci., 2009, 44: 1

[3]

H L

, Shi

X L

, Xu

X X

Effect of SiC-ZrC coating prepared by SiZr liquid phase sintering on the oxidation resistance of C/SiC composites

[J]. Acta Metall. Sin., 2021, 57: 791

DOI [本文引用: 1]

Owing to their low density and good mechanical properties at high temperatures, C/SiC composites are increasingly used in the aerospace industry. They are also being proposed as thermal-structural materials in the hypersonic field; however, C/SiC composites are easily oxidized in high-temperature air environments. In this study, a C/SiC composite was coated with a SiC-ZrC oxidation-resistant layer by a two-step sintering method using Si-25%Zr (mass fraction) alloy, and the phase evolution of the coating was studied during the sintering. The oxidation resistance of the material was then tested at 1400oC in an air environment. The microstructural changes of the coating before and after oxidation and the effect of oxidation on the bending properties of C/SiC were analyzed. After the reaction with carbon, Si and ZrSi2 disappear in the coating, leaving only pure ZrC and SiC. The ZrC phase refined the structure of the reactive SiC layer. The grain size of the sintered SiC was 2 μm, versus 5-20 μm for SiC sintered from pure Si. The refined grains created a dense and continuous SiO2 film during the oxidation process. As the oxidation time was increased at 1400oC, the C/SiC composite with the SiC-ZrC coating began losing weight at 200 s, but began gaining weight at 500 s as a dense SiO2 film was formed. After 1000 s of oxidation, the flexural strength of the C/SiC composites was 335 MPa, only 5% lower than that of the initial C/SiC composite. According to this result, the sintered SiC-ZrC oxidation-resistant film effectively protected the mechanical properties of the C/SiC composite during the oxidation process.

皮慧龙, 石小磊, 徐兴祥.

SiZr液相烧结法制备SiC-ZrC涂层对C/SiC复合材料抗氧化性能的影响

[J]. 金属学报, 2021, 57: 791

DOI [本文引用: 1]

使用Si-25%Zr (质量分数)合金通过液相烧结法在C/SiC复合材料表面制备了SiC-ZrC抗氧化涂层，研究了烧结过程中涂层的相结构演化，并测试了1400℃及空气气氛下材料的抗氧化性能，分析了涂层氧化前后显微结构的变化，以及氧化对C/SiC复合材料弯曲性能的影响。结果表明，Zr元素在涂层中以ZrC相的形式存在，ZrC颗粒的引入细化了反应SiC层的组织，在氧化过程中形成致密连续的SiO2薄膜，并在氧化500 s后试样出现增重，制备了SiC-ZrC抗氧化涂层的C/SiC样品在1000 s的氧化实验后弯曲强度下降低于5%。

[4]

Kraft

, Riedel

Numerical simulation of solid state sintering; model and application

[J]. J. Eur. Ceram. Soc., 2004, 24: 345

[5]

Van de Vorst

G A L

, Mattheij

R M M

Numerical analysis of a 2-D viscous sintering problem with non-smooth boundaries

[J]. Computing, 1992, 49: 239

[6]

Kuiken

H K

Viscous sintering: The surface-tension-driven flow of a liquid form under the influence of curvature gradients at its surface

[J]. J. Fluid Mech., 1990, 214: 503

[7]

Frenkel

J J

Viscous flow of crystalline bodies under the action of surface tension

[J]. J. Phys., 1945, 9: 385

[本文引用: 3]

[8]

Martínez-Herrera

J I

, Derby

J J

Viscous sintering of spherical particles via finite element analysis

[J]. J. Am. Ceram. Soc., 1995, 78: 645

[9]

Hiram

, Nir

A simulation of surface tension driven coalescence

[J]. J. Colloid Interface Sci., 1983, 95: 462

[10]

Gross

, Steinbach

, Raabe

, et al.

Viscous coalescence of droplets: A lattice Boltzmann study

[J]. Phys. Fluids, 2013, 25: 052101

[11]

Braginsky

, Tikare

, Olevsky

Numerical simulation of solid state sintering

[J]. Int. J. Solids Struct., 2005, 42: 621

[12]

Yang

Q C

, Kirshtein

, Ji

Y Z

, et al.

A thermodynamically consistent phase-field model for viscous sintering

[J]. J. Am. Ceram. Soc., 2019, 102: 674

DOI [本文引用: 10]

A thermodynamically consistent phase-field model for viscous sintering is proposed. It is based on an energetic variational formulation that allows the governing equations to be analytically derived from a defined energy law. The conservation of mass is satisfied through the incompressibility assumption and the assumption that mass density is uniform initially within the particle compact while the balance of linear momentum is formulated from an energy dissipation law. The morphological changes of particles are described by the temporal and spatial evolution of a phase-field variable governed by a modified Cahn-Hilliard equation, and the motion of viscous mass flow is controlled by the Stokes equation incorporating the surface tension effect. The application of the phase-field model is illustrated by examining the effect of particle shape, initial contact angle and rearrangement effects on viscous sintering.

[13]

Landau

L D

, Khalatikow

I M

. The Selected Works of L.D. Landau (English Translation) [M]. Oxford: Pergamon, 1963: 1

[14]

Sun

Z Y

, Yang

, Liu

W B

Phase field simulations of the sintering process of UO₂

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

孙正阳, 杨超, 柳文波.

UO₂烧结过程的相场模拟

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

DOI [本文引用: 1]

利用相场模型对UO2陶瓷粉末的烧结过程进行了模拟。在修正的相场模型中，不仅考虑了表面扩散、晶界扩散和晶格扩散3种各向异性的扩散机制对烧结组织形貌和烧结动力学的影响，而且考虑了不同陶瓷颗粒之间的界面能对烧结形貌的影响。基于实验条件和热力学物性参数，对UO2陶瓷粉末在2000 K的烧结过程进行了模拟。模拟结果显示：初始形貌为圆形的陶瓷粉末有利于烧结过程的进行；烧结过程中存在大晶粒吞噬小晶粒的现象；晶界扩散机制是UO2烧结过程中的主导机制；晶界能的改变导致晶界与相界之间的平衡二面角发生改变。在此基础上，模拟了多晶UO2陶瓷粉末的烧结过程，模拟结果与实验结果吻合较好。

[15]

Wang

Y U

Computer modeling and simulation of solid-state sintering: A phase field approach

[J]. Acta Mater., 2006, 54: 953

[16]

Chockalingam

, Kouznetsova

V G

, Van der Sluis

, et al.

2D phase field modeling of sintering of silver nanoparticles

[J]. Comput. Meth. Appl. Mech. Eng., 2016, 312: 492

[17]

Zhang

, Liao

Y L

A phase-field model for solid-state selective laser sintering of metallic materials

[J]. Powder Technol., 2018, 339: 677

[18]

Kumar

, Fang

Z Z

, Fife

P C

Phase field simulations of grain growth during sintering of two unequal-sized particles

[J]. Mater. Sci. Eng., 2010, A528: 254

[19]

Tegze

, Pusztai

, Gránásy

Phase field simulation of liquid phase separation with fluid flow

[J]. Mater. Sci. Eng., 2005, A413-414: 418

[20]

Soligo

, Roccon

, Soldati

Mass-conservation-improved phase field methods for turbulent multiphase flow simulation

[J]. Acta Mech., 2019, 230: 683

[21]

Gao

Y J

, Lu

Y J

, Kong

L Y

, et al.

Phase field crystal model and its application for microstructure evolution of materials

[J]. Acta Metall. Sin., 2018, 54: 278

DOI [本文引用: 2]

With the rapid development of computer technology, the roles of computer numerical simulation technology in materials are more and more prominent. Computer numerical simulation technology, real experimental observation and theoretical model analysis are the same important and are known as three great scientific research methods since the 20th century. In this paper, several important computational numerical simulation methods are briefly compared, firstly, in the spatial characteristic resolution scale and the characteristic time scale, for example, for molecular dynamics (MD), traditional phase field (TPF), and phase field crystal (PFC) method. For simulation of microstructure evolution in nano-scale, the PFC method is of the advantage on the characteristic time scale. Then, the PFC model, and its physical and mathematical basises for establishment, as well as the special feature of the method, are introduced. Next, the development of the PFC models are presented, including the PFC model of binary and multi-element alloys, of gas-liquid-solid three systems, of two-mode and multimode systems, as well as the key technology and the main procedure of the numerical calculation of the dynamic equation solution. After that, combining with the research works of the authors' group in the microstructure evolution of materials, several examples of important aspects of application of the PFC model are presented, including the nanostructure of defects of materials, dendritic growth and heterogenous epitxial growth, premelting under deformation at high temperature and dynamic recovery, extension and bifurcation of cracks on nanoscale, matalllic glass transition, defect structures of graphene, voids formation of electromigration in metal interconnects, microstructure in multiferroic composite matrials, and the formation of the structure of the metal foams. Finally, a summary is given and the development direction and future emphasis application and new fields of the PFC model are pointed out.

高英俊, 卢昱江, 孔令一

等.

晶体相场模型及其在材料微结构演化中的应用

[J]. 金属学报, 2018, 54: 278

DOI [本文引用: 2]

随着计算机技术的快速发展,计算机模拟实验在材料科学中的作用越来越突出。计算机数值模拟技术已经和实验观测、理论模型分析并称为20世纪以来的三大科学研究方法。本文首先简要地从空间特征分辨尺度和时间特征尺度比较了几种重要的计算模拟方法——分子动力学(MD)、传统相场方法(TPF)和晶体相场(PFC)方法的各自适应的特征尺度范围和特点。在模拟纳观尺度的材料微结构演化,PFC在特征时间尺度上更具优势。其次,介绍了PFC模型,及其建立的物理基础和数学基础,以及该方法的特色优势。同时,介绍该PFC模型的拓展与推广,包括二元和多元体系、气-液-固三相体系、双模和多模体系的PFC模型,以及求解PFC模型的动力学方程数值计算的关键技术与主要步骤。再次,结合作者在材料微结构演化方面的研究,着重介绍PFC模型的几个重要方面的应用例子,包括材料纳观缺陷结构演化、凝固的枝晶生长和晶体外延生长、高温预熔化变形和动态回复、纳观尺度的裂纹扩展与分叉、无序-有序金属玻璃转变、石墨烯的缺陷结构、金属互联线电迁移空洞、多铁复合材料的畴结构、金属泡沫结构的生成等。最后,总结并指出PFC模型的拓展方向与今后应用的重点方面和新领域。

[22]

Zografos

, Afonso

A M

, Poole

R J

, et al.

A viscoelastic two-phase solver using a phase-field approach

[J]. J. Non-Newton. Fluid Mech., 2020, 284: 104364

[23]

Balemans

, Hulsen

M A

, Anderson

P D

Sintering of two viscoelastic particles: A computational approach

[J]. Appl. Sci., 2017, 7: 516

[本文引用: 9]

[24]

Kim

Phase-field models for multi-component fluid flows

[J]. Commun. Comput. Phys., 2012, 12: 613

[25]

Lin

Y Y

, Hui

C Y

, Jagota

The role of viscoelastic adhesive contact in the sintering of polymeric particles

[J]. J. Colloid Interface Sci., 2001, 237: 267

[26]

Kamyabi

, Sotudeh-Gharebagh

, Zarghami

, et al.

Principles of viscous sintering in amorphous powders: A critical review

[J]. Chem. Eng. Res. Des., 2017, 125: 328

[27]

Wakai

, Katsura

, Kanchika

, et al.

Sintering force behind the viscous sintering of two particles

[J]. Acta Mater., 2016, 109: 292