相场理论基础-Foundation of Phase Field Modeling

Foundation of Phase Field Modeling

你好！ 这是你第一次使用 Markdown编辑器 所展示的欢迎页。如果你想学习如何使用Markdown编辑器, 可以仔细阅读这篇文章，了解一下Markdown的基本语法知识。

基本原理

对于特定物理体系的状态或过程，平衡状态可以视为有约束（或无约束）条件下的全局能量最小化的结果。 基于广义力或通量的运动学模型是用来表示微结构演化的经典模型。

相场参量（序参量）

在相场方法中，根据相场参量的取值而区分不同的相。 对于特定的物理问题，需先选择一个变量作为相场参量，它可以是模型中的具体物理量，也可以是一组变量的函数。通过定义好的相场参量，可以反应体系的物理性质，比如结晶度（degree of crystallinity）、原子排列（atomic ordering）、晶体对称性（crystal symmetry）、状态（state）、磁化强度（magnetization）、电极化量（electric polarization）、浓度（concentration）、组成（composition）、损伤程度（damage level）等。

因为相场法研究的更多是微观/介观尺度，而非原子/电子尺度，所以相场参量是相的一种特定热力学性质的平均值。相场参量也被称为序参量。

扩散界面模型

在传统相变模型中，常因为追踪不同相之间界面的变化而导致大量繁琐的计算。相场法假定相之间的界面是扩散性（diffusive）的，介于相邻的相之间的扩散界面由用 $\varphi$ 表示的相场参量的连续插值描述。

在扩散界面模型下，微观结构类似于相场变量的一个分布函数 $\varphi (x)$，其中 $x$是位置向量。这种模型提供了用公式表示体积内能量密度的可能。在传统的尖锐界面下，界面对能量的贡献与界面能密度和界面拓扑结构有关，而界面形态（interface morphology）轨迹，特别是三维情况下，可能会十分复杂。而在扩散界面模型下，界面能够通过相场参量的梯度进行识别。因此界面形态和界面能量都能够很自然地基于梯度项定义。

自由能泛函

每个相的自由能密度都依赖于相的热力学状态，因此每个相的体自由能密度（bulk free energy density）都能够由相场参量$\varphi$ 用公式描述。界面上的自由能密度是不同于相内部的，所以界面位置能够用相场参量的梯度来描述。系统 A 的自由能 $F$ 能够由以下相场自由能泛函给定：
$$F\left(\varphi ,\nabla \varphi \right)={\int }_A\left[f\left(\varphi \right)+\frac{1}{2}\kappa {\left|\nabla \varphi \right|}^2\right]dV$$其中，$f\left(\varphi \right)$ 是体自由能密度， $\frac{1}{2}\kappa {\left|\nabla \varphi \right|}^2$是包含参量$\kappa$的梯度能密度。

体积自由能密度$f\left(\varphi \right)$ 依赖于局域的相场参量。例如，一个多相系统的体积自由能密度图像通常呈多势阱的形状。梯度项$\frac{1}{2}\kappa {\left|\nabla \varphi \right|}^2$是非局域的，并且会在界面处引起能量损耗。需要注意的是，仅当$\kappa \to 0$ 时，尖锐界面才是可能的。

另外，在某些问题中，还需要考虑梯度的高阶项，例如$|\nabla \varphi {|}^4$。由于界面能的空间对称性，更高阶项应该是偶数阶的。

相场法的优势

在相场理论中，系统的状态和动力学由能量最小原理决定，系统处于平衡状态时，序参量${\varphi }_{eq}\left(x\right)$应当使得自由能 $F$ 取得最小值，即自由能的一阶变分导数为零：$$\frac{\delta F\left(\varphi ,\nabla \varphi \right)}{\delta \varphi }=0$$

对于不受约束的热力学量，可以假定微结构演化过程是以整体能量或势最快降低的方向进行的。在相场理论中，这种关于微结构变的整体能量变化即由这个一阶变分导数给定。因此相场参量下的动力学模型能够直接由该导数用公式明确表示。注意到该一阶变分导数能够由下述变分得到：

$$\begin{array}{rl}\delta F& =\delta {\int }_A\left[f\left(\varphi \right)+\frac{1}{2}\kappa {\left|\nabla \varphi \right|}^2\right]\mathrm{d}V\\ & ={\int }_A\left[\delta f\left(\varphi \right)+\frac{1}{2}\kappa \delta {\left|\nabla \varphi \right|}^2\right]\mathrm{d}V\\ & ={\int }_A\left[\frac{\partial f\left(\varphi \right)}{\partial \varphi }\delta \varphi +\kappa {\partial }_i\varphi {\partial }_i\delta \varphi \right]\mathrm{d}V\\ & ={\int }_A\left\{\frac{\partial f\left(\varphi \right)}{\partial \varphi }\delta \varphi +\kappa \left[{\partial }_i\left({\partial }_i\varphi \delta \varphi \right)-{\partial }_i{\partial }_i\varphi \delta \varphi \right]\right\}\mathrm{d}V\\ & ={\int }_A\left[\frac{\partial f\left(\varphi \right)}{\partial \varphi }-\kappa {\partial }_{ii}\varphi \right]\delta \varphi \mathrm{d}V+{\int }_A\kappa {\partial }_i\left({\partial }_i\varphi \delta \varphi \right)\mathrm{d}V\\ & ={\int }_A\left[\frac{\partial f\left(\varphi \right)}{\partial \varphi }-\kappa {\partial }_{ii}\varphi \right]\delta \varphi \mathrm{d}V+{\int }_{\partial A}\kappa {\partial }_i\varphi n_i\delta \varphi \mathrm{d}S\\ & ={\int }_A\left[\frac{\partial f\left(\varphi \right)}{\partial \varphi }-\kappa \Delta \varphi \right]\delta \varphi \mathrm{d}V+{\int }_{\partial A}\kappa \nabla \varphi \cdot \mathbf{n}\delta \varphi \mathrm{d}S\end{array}$$

给定$\kappa \nabla \varphi \cdot \mathbf{n}=0$，则由上式可得：
$$\frac{\delta F}{\delta \varphi }=\frac{\partial f\left(\varphi \right)}{\partial \varphi }-\kappa \Delta \varphi$$ 其中，$\Delta$是拉普拉斯算子。

体自由能密度

体自由能密度$f\left(\varphi \right)$依赖于相场参量的意义。下面举例一些常用的体积自由能。

1. 二元混合物中的相分离

在这个问题中，我们将其中一种组分的密度作为相场参量$\varphi$ ，基于二元混合物的平均场自由能理论，体积自由能采用以下形式：
$$f\left(\varphi \right)=\frac{\epsilon \nu }{2}\varphi \left(1-\varphi \right)+k_{B}T\left[\varphi \ln \varphi +\left(1-\varphi \right)\ln \left(1-\varphi \right)\right]+b\varphi -\frac{\nu {\epsilon }_{AA}}{2}$$ 其中， $k_B$ 是玻尔兹曼常量，$T$是温度， $\epsilon ,{\epsilon }_{AA},\nu$ 是材料参数。

如上图左可见，随着温度 $T$ 升高并超过临界温度 $T_C=\frac{\epsilon \nu }{4k_B}$ ，体积自由能也从双势阱变为了单势阱。低于临界温度时，存在两种平衡状态，分别对应两个势阱。当 $b=0$时，这两个平衡状态实际上是局部最小值状态。而对于$b\ne 0$，这两种平衡状态能够被余切定理控制，如上图右所示。当$T\to T_C$时，这两个状态连续地汇聚为一个。

2. 相变的朗道自由能泛函

在朗道理论中，假定自由能在相 $\varphi =0$ 处的级数展开：
$$f\left(\varphi \right)=f\left(T,\varphi =0\right)+\sum _{n=2}^M\frac{a_n\left(T\right)}{n}{\varphi }^n$$
对于对称相变，可以使用
$$f\left(\varphi \right)=a\left(T\right)+\frac{a_2\left(T\right)}{2}{\varphi }^2+\frac{a_4\left(T\right)}{4}{\varphi }^4+O\left(\varphi \right)$$
其中，系数$a_n\left(T\right)$也能够被级数展开为线性项。例如
$$\begin{array}{rl}& a_2\left(T\right)=a_2^0\left(T-T_C\right)\\ & a_4\left(T\right)=a_4^0+b_4^0\left(T-T_C\right)\end{array}$$ 高于临界温度时，$a_2\left(T\right),a_4\left(T\right)>0$。

当温度高于临界温度时，只存在一个稳态 $\varphi =0$，而当温度低于临界温度时，存在两个对称的稳态 $$\varphi \approx \pm \sqrt{\frac{a_2^0}{a_4^0}\left(T_C-T\right)}$$
当$T\to T_C$ 时，两个稳态向 $\varphi =0$ 汇聚。上述代表了一个典型的二阶相变。

通过增加奇数阶项，可以获得非对称相图的相变体积自由能。此外，通过调节$a_n\left(T\right)$ 的符号也能获得一阶相变的多项式。

3. 凝固

对纯材料单晶的凝固过程，体自由能密度可以设为以下形式：
$$f\left(\varphi ,T\right)=f_0-S_L\left(T-T_m\right)+H{\varphi }^2{\left(1-\varphi \right)}^2+\left[S_L-\frac{L}{T_m}\left(3-2\varphi \right){\varphi }^2\right]{\varphi }^2\left(T-T_m\right)$$ 其中，$S_L$ 是液相的体积熵密度； $T_m$是熔点；$H$ 是双势阱高度； $L$ 是熔化潜热。
右式第三项形成了一个双势阱，而最后一项破坏了关于熔点的对称性，如下图所示：

实际上，使用多项式来形成双势阱体自由能泛函在相场模型中是非常常见的。

平衡

平衡微结构的系统全局能量或势最小化过程依赖于约束条件。一般来说，在平衡状态下，相场参量的驱动力或通量会变为 。而在相场模型下，驱动力或通量可以表示为自由能泛函关于相场参量的一阶变分导数，所以平衡状态微结构应该满足以下关系式：
$$\frac{\delta F}{\delta \varphi }=\frac{\partial f\left(\varphi \right)}{\partial \varphi }-\kappa \Delta \varphi =0$$
以一个简单的体积自由能的多项式双势阱函数$f\left(\varphi \right)=4H{\varphi }^2{\left(1-\varphi \right)}^2$ 为例。它有两个稳定的相$\varphi =0$，$\varphi =1$，并且能量势垒为$H$。一个边界条件为${\varphi |}_{x=-\infty }=0$，${\varphi |}_{x=+\infty }=1$的一维问题，对上一个等式的解为：
$${\varphi }_{eq}\left(x\right)=\frac{1}{2}\left[1+\tanh \left(\frac{x}{2l}\right)\right]$$
式中，$$l=\sqrt{\frac{\kappa }{H}}$$

体自由能和序参量的图像如下
该模型中包括两个序参量：$H$代表稳态之间的势垒，$\kappa $为梯度能项的系数。另外，模型中包含两个界面相关的物理参量：$l$为界面厚度，$G$为界面能密度。

对于给定的双势阱函数$f\left(\varphi \right)=4H{\varphi }^2{\left(1-\varphi \right)}^2$，有$l=\sqrt{\frac{\kappa }{H}}$。

另外，因为两个稳态处的体自由能密度为$0$，所以体自由能密度$f\left({\varphi }_{eq}\right)+\frac{1}{2}\kappa {\left|\nabla {\varphi }_{eq}\left(x\right)\right|}^2$ 在整个域上的积分即为界面上的界面能：$$G={\int }_{-\infty }^{+\infty }\left[f\left({\varphi }_{eq}\left(x\right)\right)+\frac{1}{2}\kappa {\left|\nabla {\varphi }_{eq}\left(x\right)\right|}^2\right]\mathrm{d}x=\frac{1}{3}\sqrt{\kappa H}$$
对于一般形式的体积自由能，界面参量依赖于模型参量：
$$l\sqrt{\frac{\kappa }{H}}，l\sqrt{\kappa H}$$

用界面参量$G$，$l$ 来表达的相场模型有更加明晰的物理解释，比如裂纹扩展的相场模型经常用如下形式表示：$$F={\int }_A\frac{G}{l}\left(g\left(\varphi \right)+Gl{\left|\nabla {\varphi }_{eq}\left(x\right)\right|}^2\right)\mathrm{d}V$$

如果需要同时考虑超过一种物理特性，可以直接在体积自由能中增加相应的能量项，如果需要的话也可以添加耦合项，所以，相场模型能够被扩展应用到解决多物理特性问题上。

相场演化动力学

微观结构的演化能够通过定义了与时间相关的相场变化率的运动学模型进行模拟。根据序参量所受到的约束不同，对应的运动学模型也不同。

Allen-Cahn方程

根据 Allen-Cahn 模型（也称 Langevin-type 模型），对于无约束的序参量，微结构$\varphi \left(x\right)$ 以最陡的泛函梯度方向使全局自由能最小化（或全局系统熵最大化），即
$$\frac{\partial \varphi }{\partial t}=-M\frac{\partial F}{\partial \varphi }=-M\left(\frac{\delta f\left(\varphi \right)}{\delta \varphi }-\kappa \Delta \varphi \right)$$ 其中，$M$ 是迁移率参数。上式是一个二阶偏微分方程。

值得注意的是，该模型只适用于纯弛豫过程的序参量。例如，非保守序参量的磁化强度、损伤变量和凝固。

Cahn-Hilliard 方程

如果序参量遵循以下形式的通量守恒法则：
$$\frac{\partial \varphi }{\partial t}=-\nabla \cdot \mathbf{j}，\mathbf{j}=-D\nabla \mu$$
则称序参量是保守的。

其中，$\mathbf{j}$是通量向量，它依赖于化学势 $\mu$ 的梯度。在这种情况下，假定化学势采取以下一阶变分导数形式：
$$\mu =\frac{\delta F}{\delta \varphi }$$
因此
$$\frac{\partial \varphi }{\partial t}=\nabla \cdot M\nabla \frac{\partial F}{\partial \varphi }$$ 这是一个四阶微分方程。

值得说明的是，浓度是典型的保守序参量。

Landau-Lifshitz-Gilbert 方程

在一些弛豫过程中，相场向量$M$ 的模恒定。特别地，对于温度远低于居里点$T_C$ 的磁化过程，磁化强度几乎不变。因此该磁化仅仅在微结构演化中旋转。对应的演化方程由 L-L-G 方程给定：
$$\frac{\partial M}{\partial t}=-{\gamma }_{1}M\times \frac{\delta F}{\delta M}-{\gamma }_{2}M\times \left(M\times \frac{\delta F}{\delta M}\right)$$
叉乘保证了磁化向量的驱动力与该向量保持垂直。

随机涨落

如果演化动力学方程中包含了热扰动或化学干扰，则需要考虑序参量的涨落，此时可以引入一个随机噪声项 。例如，对于 A-C 型方程：
$$\frac{\partial \varphi }{\partial t}=-M\frac{\partial F}{\partial \varphi }+\eta \left(x,t\right)$$
根据扰动产生的原因，噪声项应该遵循特定的时空随机分布。

界面特性

界面的各向异性和迁移率对演化的影响十分重要。在特定的情况下，梯度项系数 和迁移率参数 可以由晶格取向，甚至由相场参量梯度方向表示的界面取向决定。例如：$$M\to M\left(\theta \right)，\kappa \to \kappa \left(\theta \right)，\theta =\arctan \left(\nabla \varphi \right)$$

相场模型的计算

有限差分算法

有限差分算法为偏微分方程的数值解提供了比任何其他方法都更直接的方法。有限差分算法编码简单、计算经济、易于并行，适用于分布式计算环境。然而，它们在精度和施加复杂边界条件方面也有缺点。

有限差分算法的核心思想是将每个导数替换为一个差商。对于一维网格上连续且相对光滑的函数$u$，网格点间距为$h$，有限差分算法采用如下公式：

后向差分：$${\left(\Delta u\right)}_i^{-}=\frac{u_i-u_{i-1}}{h}$$
前向差分：$${\left(\Delta u\right)}_i^{+}=\frac{u_{i+1}-u_i}{h}$$
中心差分：$${\left(\Delta u\right)}_i^{\pm }=\frac{u_{i+1}-u_{i-1}}{2h}$$
二阶中心差分：$${\left({\Delta }^{2}u\right)}_i=\frac{u_{i+1}-2u_i+u_{i-1}}{h^2}$$

经过泰勒展开后，各差分的误差项可表示为如下形式：
$$\begin{array}{rl}& {\left(\Delta u\right)}_i^{-}-u^{\prime }\left(x_i\right)=O\left(h\right)\\ & {\left(\Delta u\right)}_i^{+}-u^{\prime }\left(x_i\right)=O\left(h\right)\\ & {\left(\Delta u\right)}_i^{\pm }-u^{\prime }\left(x_i\right)=O\left(h^2\right)\\ & {\left({\Delta }^{2}u\right)}_i-u^{\prime \prime }\left(x_i\right)=O\left(h^2\right)\end{array}$$

这些差分公式，特别是近似拉普拉斯算子的中心二阶差分，将用于相场模型的数值解。
对于二维网格，在网格内点$\left(x_i,y_j\right)$处，拉普拉斯算子为：
$$\begin{array}{rl}{\left({\nabla }^{2}u\right)}_{i,j}& ={\left({\Delta }_{xx}^{2}u\right)}_{i,j}+{\left({\Delta }_{yy}^{2}u\right)}_{i,j}\\ & =\frac{u_{i+1,j}-2u_{i,j}+u_{i-1,j}}{h_x^2}+\frac{u_{i,j+1}-2u_{i,j}+u_{i,j-1}}{h_y^2}\end{array}$$

这种近似称为五点法，只涉及五个格点，这可以从方程和图1中推断出来。

九点法可以离散化误差降低到$O\left(h^4\right)$，因此可以使用更大的$h$，而误差小于五点法。然而，这种情况会导致计算成本的增加。因此，在相场模型求解时，通常选用五点法来进行空间离散化。

如果在空间中 $x$ 和 $y$ 方向上的格点间距彼此相等，$h=h_x=h_y$，则前式可以继续化简为：
$${\left({\nabla }^{2}u\right)}_{i,j}=\frac{u_{i+1,j}+u_{i-1,j}+u_{i,j+1}+u_{i,j-1}-4u_{i,j}}{h^2}$$

对于图1所示等间距的9个格点，在不考虑任何边界条件的情况下，可以使用矩阵向量积的形式将上式重写为：

$${\left({\nabla }^{2}u\right)}_{i,j}=\frac{1}{h^2}\left(\begin{array}{ccccccccc}-4& 1& 0& 1& 0& 0& 0& 0& 0\\ 1& -4& 1& 0& 1& 0& 0& 0& 0\\ 0& 1& -4& 0& 0& 1& 0& 0& 0\\ 1& 0& 0& -4& 1& 0& 1& 0& 0\\ 0& 1& 0& 1& -4& 1& 0& 1& 0\\ 0& 0& 1& 0& 1& -4& 0& 0& 1\\ 0& 0& 0& 1& 0& 0& -4& 1& 0\\ 0& 0& 0& 0& 1& 0& 1& -4& 1\\ 0& 0& 0& 0& 0& 1& 0& 1& -4\end{array}\right)\left(\begin{array}{c}{u}_{1,1}\\ u_{2,1}\\ u_{3,1}\\ u_{1,2}\\ u_{2,2}\\ u_{3,2}\\ u_{1,3}\\ u_{2,3}\\ u_{3,3}\end{array}\right)\text{\hspace{1em}(1)}$$

通常，该公式可以重写为：
$${\left({\nabla }^{2}u\right)}_{i,j}=\frac{1}{h^2}M{U}_{u_{i,j}}$$
式中，
$$M=\left(\begin{array}{cccccc}A& I& O& & O& O\\ I& A& I& \cdots & O& O\\ O& I& A& & O& O\\ & ⋮& & \ddots & ⋮& \\ O& O& O& \cdots & I& A\end{array}\right)$$ 其中，$I$是一个$3\times 3$ 的单位阵，$O$是一个$3\times 3$ 的零矩阵，而A为：
$$A=\left(\begin{array}{ccc}-4& 1& 0\\ 1& -4& 1\\ 0& 1& -4\end{array}\right)$$
值得注意的是，$M$具有一些独特的性质：

• $M$是一个稀疏矩阵，每行最多只有5个非零元素；
• 它是一个块三对角矩阵（block tri-diagonal），有三对角线和对角线块；
• 对称；
• 它是一个对角占优矩阵（diagonally dominant matrix）：矩阵中每个主对角元素的模都大于与它同行的其他元素的模的总和；
• $M$的对角元素为负，而其他元素为正；
• 它是一个负定矩阵：其特征值均为负。

应用

一维瞬态热传导

$$\rho c_{\rho }\frac{\partial T}{\partial t}=\frac{\partial }{\partial x}\left(\lambda \frac{\partial T}{\partial x}\right)$$ $$T_i^{n+1}=T_i^n+\mu \Delta t\left(\frac{T_{i+1}^n-2T_i^n+T_{i-1}^n}{\Delta x^2}\right)$$

二元合金的斯皮诺达分解（Spinodal Decomposition）

使用有限差分算法在相场模型中求解保守Cahn-Hilliard方程。

晶粒生长

使用有限差分算法在相场模型中求解多分量非守恒Allen-Cahn方程。

固相烧结

使用有限差分法求解守恒Cahn-Hillard方程和多个非守恒Allen-Cahn方程的耦合解。
该模型基于两类场变量对微观结构的表示。一个是守恒密度场$\rho$，它在固相处为1，在孔隙处为0，在固-孔界面处迅速变化。另一个是非守恒阶参数${\eta }_i$，用于区分微观结构中的不同颗粒。

枝晶凝固

处理界面能和晶界迁移率的各向异性。

多细胞系统的结构演化

通过修改自由能泛函，使相场模型中用非守恒序参数描述的每个单元在保留其原始体积的同时，能够发生大变形和形状变化。公式仍遵循总能量最小化的原则，但是当细胞聚集在一起时，它们会发生集体变形，而不是收缩、生长或相互吸收。

傅里叶谱法

应用

多元合金的析出行为

在固态微观结构演变过程中，多元合金系统会同时发生相变和相分离。相变采用非守恒的Allen - Cahn公式，相分离过程中合金成分守恒采用Cahn - hilliard公式。

二元合金相分离行为

在微观组织演化的过程中，演化相的力学性能和施加的应力场的不匹配将导致新的动态力学平衡出现。

合金中的晶格缺陷

在考虑弹性不均匀性的基础上，还将考虑晶格缺陷：位错和晶界。

有限元方法