析出过程的相场模拟模拟合金中第二相颗粒（沉淀物）嵌入在基体中的情况。由于两相晶格常数不同，沉淀物在形成时会发生晶格畸变，

第一部分：模型的核心思想与总自由能泛函

1.1 物理问题背景

模拟合金中第二相颗粒（沉淀物）嵌入在基体中的情况。由于两相晶格常数不同，沉淀物在形成时会发生晶格畸变，产生内应力。这种弹性应变能会显著影响沉淀物的平衡形状和演化动力学。相场法的优势在于能自然地模拟复杂形状的变化，如分裂、合并。

1.2 序参量 (Order Parameter) `η`

定义：引入一个连续的标量场 η(x, y, z, t) 来描述系统的相。
物理意义：
- η = 0：代表基体相 (Matrix Phase)。
- η = 1：代表沉淀物相 (Precipitate Phase)。
- 0 < η < 1：代表两相之间的扩散界面 (Diffuse Interface)。这与现实中的原子级别过渡层相对应，是相场法区别于 sharp-interface 方法的核心。
为什么需要它：它使我们能用一组在空间上光滑、且易于进行数值计算的场变量来描述复杂的微观结构。

1.3 总自由能泛函 (Total Free Energy Functional) `ℱ`

系统的总能量是系统中所有能量贡献的积分。教程中给出的公式是基石：

$\mathcal{F} = \int_{V} \left[ \underbrace{f_{\text{bulk}}(\eta)}_{\text{体化学能}} + \underbrace{\frac{\kappa}{2}|\nabla \eta|^{2}}_{\text{梯度界面能}} + \underbrace{f_{\text{el}}(\eta)}_{\text{弹性应变能}} \right] dV$

物理意义：系统总是趋向于使总自由能 ℱ 最小化的状态。这个泛函决定了 η 场和位移场 uᵢ 将如何演化。

f_bulk(η)：体化学自由能密度。它代表了不考虑界面效应时，材料本身是基体相还是沉淀相所具有的化学能。
(κ/2)|∇η|²：梯度能量密度。它惩罚 η 场在空间上的剧烈变化，从而定义了界面。|∇η| 越大（即界面越“陡峭”），这项能量越高。系统为了降低能量，会倾向于使界面变得平滑，从而产生一个具有一定宽度的扩散界面。κ 与界面能 γ 和界面宽度 l 直接相关。
f_el(η)：弹性应变能密度。这是本问题的核心，代表了由于晶格错配而产生的弹性畸变所储存的能量。

第二部分：各项能量的详细推导与解释

2.1 体化学能项 `f_bulk(η)` 的详细解释

数学形式： $f_{\text{bulk}} = w \sum_{j=0}^{10} a_j \eta^j$

w：势垒高度系数。控制双势阱的能垒高度，单位是 aJ/nm³。它影响了相变的驱动力。
a_j：多项式系数。这些系数是经过精心设计的，以满足以下关键边界条件：
1. f_bulk(0) = 0：基体相的参考化学能为零。
2. f_bulk(1) = 0：沉淀物相的参考化学能为零。（这意味着两相在无应变、平界面时处于平衡）。
3. f'_{bulk}(0) = 0 和 f'_{bulk}(1) = 0：在平衡点处化学势为零。这是平衡态的热力学要求。
4. 在 η=0 和 η=1 处是能量极小值。
5. 在两个极小值之间是一个很高的能垒，这确保了在大多数区域内，η 的值会非常稳定地保持在 0 或 1 附近，只有在界面处才发生快速变化。高阶多项式（10阶）可以构造出非常“深”且“窄”的势阱。

物理意义：f_bulk 就像一个“双坑”地形。小球（系统）自然倾向于待在坑底（η=0 或 η=1）。这个项本身不决定平衡形状，但它提供了相变的倾向性。

2.2 梯度能项 `(κ/2)|∇η|²` 的详细解释

数学形式：简单明了，即梯度的平方。

κ：梯度能量系数。单位是 aJ/nm。它是一个材料参数，与界面能 γ 和界面宽度 l 有关。大致关系为 γ ∝ √(κ w)。
|∇η|：序参量梯度的模。在界面处最大，在相内部为零。

物理意义：这项能量代表了创建界面所需要的能量。系统为了降低总能量，会倾向于减少界面面积（这是导致沉淀物粗化的驱动力），但同时受限于界面宽度。这项与 f_bulk 共同决定了平衡界面的轮廓和能量。

2.3 弹性应变能项 `f_el(η)` 的详细推导

这是最复杂的部分，我们一步步拆解。

第1步：弹性应变能密度的一般形式 在弹性力学中，单位体积内储存的应变能（弹性势能）为： $f_{\text{el}} = \frac{1}{2} \sigma_{ij} \epsilon_{ij}^{\text{el}}$ 其中采用了爱因斯坦求和约定。对于线性弹性材料，这等价于 (1/2) Cᵢⱼₖₗ εᵢⱼᵉˡ εₖₗᵉˡ。

第2步：本构关系 (Constitutive Relation) - 胡克定律 应力 σ 和弹性应变 εᵉˡ 通过刚度张量 C 联系起来： $\sigma_{ij} = C_{ijkl} \epsilon_{kl}^{\text{el}}$ 关键点：在相场模型中，刚度张量 C 是序参量 η 的函数！因为界面处的材料性能是两相的混合。

插值方法： $C_{ijkl}(\eta) = C_{ijkl}^{\text{matrix}} \left[1 - h(\eta)\right] + C_{ijkl}^{\text{precip}} h(\eta)$

C_{ijkl}^{matrix}：纯基体相的刚度张量。
C_{ijkl}^{precip}：纯沉淀物相的刚度张量。
h(η)：插值函数。例如 h(η) = η³(6η² - 15η + 10)。

为什么是这个奇怪的函数？ 因为它满足以下完美的性质：

h(0) = 0, h(1) = 1。
h'(0) = 0, h'(1) = 0。这意味着在纯基体 (η=0) 和纯沉淀物 (η=1) 内部，刚度对 η 的变化率是零。这可以防止在相内部产生虚假的“微观结构”，确保材料性能的变化只集中在界面区域。

第3步：弹性应变、总应变与本征应变 (Eigenstrain) 这是弹性理论的核心。弹性应变 εᵉˡ 并不等于总应变 εᵗᵒᵗᵃˡ。

$\epsilon_{ij}^{\text{el}} = \epsilon_{ij}^{\text{total}} - \epsilon_{ij}^{0}(\eta)$

总应变 εᵗᵒᵗᵃˡ：由位移场 u 导出的几何上的应变。 $\epsilon_{ij}^{\text{total}} = \frac{1}{2} \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right)$ 这是小变形下的柯西应变张量。
本征应变 ε⁰(η)：应力自由应变。想象一下，如果从材料中切下一个微小单元，没有任何外部约束，它“想要”自发产生的应变。对于相变问题，这通常就是晶格错配应变 εᵀ。

$\epsilon_{ij}^{0}(\eta) = \epsilon_{ij}^{T} \, h(\eta)$

εᵀ 是常数张量，表示纯沉淀物 (η=1) 相对于基体 (η=0) 的错配。比如ε₁₁ᵀ = ε₂₂ᵀ = 0.005 (0.5%), ε₁₂ᵀ = 0。
同样使用 h(η) 进行插值，确保在相内部没有不必要的变化。

物理图像：在沉淀物内部 (η=1)，材料“想要”应变达到 εᵀ。但因为它被周围的基体牢牢粘住，无法自由膨胀，实际产生的总应变 εᵗᵒᵗᵃˡ 小于 εᵀ，这就产生了压缩应力 (εᵉˡ = εᵗᵒᵗᵃˡ - εᵀ < 0)。相反，基体会被沉淀物拉伸，产生拉伸应力。整个系统处于一种相干应变 (Coherent Strain) 状态。

第三部分：演化方程与求解策略

系统通过演化来降低总能量 ℱ。

3.1 序参量演化：Cahn-Hilliard 方程

$\frac{\partial \eta}{\partial t} = \nabla \cdot \left[ M \nabla \left( \frac{\delta \mathcal{F}}{\delta \eta} \right) \right]$

M：迁移率 (Mobility)。是一个动力学系数，控制界面移动的快慢。
δℱ/δη：变分导数/化学势 μ。是系统总自由能对序参量场变化的“敏感度”。

化学势 μ 的详细计算： $\mu \equiv \frac{\delta \mathcal{F}}{\delta \eta} = \frac{\partial f_{\text{bulk}}}{\partial \eta} + \frac{\partial f_{\text{el}}}{\partial \eta} - \kappa \nabla^{2} \eta$

∂f_bulk/∂η：体化学能的驱动力。试图将 η 拉向 0 或 1。
-κ∇²η：界面能的贡献。试图平滑化 η 场，减少曲率。
∂f_el/∂η：弹性驱动力。这是最复杂的一项。 $\frac{\partial f_{\text{el}}}{\partial \eta} = \frac{\partial}{\partial \eta} \left( \frac{1}{2} C_{ijkl}(\eta) \epsilon_{ij}^{\text{el}} \epsilon_{kl}^{\text{el}} \right)$ 利用链式法则，这项包含两部分：
- 一部分来自于刚度 C(η) 随 η 的变化。
- 另一部分来自于本征应变 ε⁰(η) 随 η 的变化（因为它影响了 εᵉˡ）。 这项解释了为什么弹性应力会改变平衡形状，例如，它可能使沉淀物沿着某些晶向拉长以降低总弹性能。

Cahn-Hilliard 方程是守恒的，∫η dV 不随时间变化，这模拟了沉淀物中原子总量的守恒。

3.2 机械平衡：应力平衡方程

弹性弛豫的速度远快于扩散过程（声速 vs. 原子扩散速率）。因此，在每个时间步，我们都假设机械平衡是瞬间建立的。这意味着对于当前的 η 场，位移场 u 总是调整到使弹性能最小。

这由以下方程描述： $\nabla \cdot \sigma_{ij} = 0 \quad \text{或} \quad \frac{\partial \sigma_{ij}}{\partial x_j} = 0$ 这就是弹性力学中的平衡方程，在无体力的条件下表示内力（应力）的散度为零，即每个体积元都处于静力平衡状态。

求解策略 (Operator Splitting)：

给定一个初始的 η 场。
求解力学平衡方程 ∇·σ=0，得到位移场 u 和应力场 σ。
计算当前的化学势 μ（它依赖于 σ 和 η）。
根据 Cahn-Hilliard 方程演化一个时间步 Δt，得到新的 η 场。
回到第2步，循环直至系统总能量不再显著下降（达到平衡）。

第四部分：参数与数据的物理意义详解

能量单位 aJ (attojoule)：1 aJ = 10⁻¹⁸ J。在纳米尺度模拟中，这个单位很方便。
刚度单位 aJ/nm³：1 aJ/nm³ = 10⁻¹⁸ J / 10⁻²⁷ m³ = 10⁹ J/m³ = 1 GPa。这与常用的弹性模量单位一致。
界面能 50 aJ/nm²：50 aJ/nm² = 50 × 10⁻¹⁸ J / 10⁻¹⁸ m² = 50 mJ/m²，这是一个合理的表面能数值。
错配应变 0.5%：εᵀ = 0.005，是一个典型的数值，足以产生显著的弹性效应，又不会大到导致失配位错。
初始值 η₀^matrix = 0.0065：这是一个数值技巧。因为弹性能的引入改变了化学势，如果严格将基体设为 0，沉淀物可能会完全溶解以使基体化学势达到新的平衡值。稍微提高基体的 η 初始值，相当于预先设定了一个“背景浓度”，可以防止这种非物理的溶解，同时锁定沉淀物的体积大致不变。

第五部分：总结与物理图像

整个模型可以这样理解：系统有两个主要的能量竞争方：

界面能：倾向于让沉淀物变成球形（因为相同体积下，球形的表面积最小）。
弹性应变能：倾向于让沉淀物变成某种非球形（例如，正方形或椭圆形），以更好地适应晶格错配，降低整体的应变能。

最终的平衡形状是这两种能量（以及可能的各向异性）竞争和妥协的结果。您的相场模拟，就是通过求解上述耦合方程组，来揭示这一最小能量状态及其演化路径。

希望这份极其详尽的解释能帮助您彻底理解这个模型。如果您对某个特定部分的代码实现有疑问，我们可以继续深入探讨。

参考文献

根据 PFHub: The Phase Field Community Hub ，经过deepseek生成后，笔者手动校验并发布。

析出过程的相场模拟