📚 课程目录第4章相场与CALPHAD的耦合接口设计

4. 相场与CALPHAD的耦合接口设计

相场方程描述的是“界面如何演化”，CALPHAD描述的是“给定成分与温度下体系有多稳定”。两者要真正耦合，关键不在求解器，而在接口：相场模型中的每一个热力学量，都必须能从CALPHAD数据库中以一致、可微、可插值的方式取出，并映射到相场变量与成分场上。接口设计一旦含糊，就会出现驱动力符号错误、界面异常展宽、溶质截留或守恒性破坏等问题。

一、耦合的物理量链条

在典型的二元或多元合金相场模型中，总自由能泛函通常写成：

F = ∫ [ f_chem(c, η, T) + f_grad(∇η) + f_elastic(...) ] dV

其中 f_chem 是化学自由能密度，η 是相场变量（如 η=0 表示基体相，η=1 表示析出相），c 是成分场。CALPHAD提供的是各相摩尔Gibbs自由能 G_m^φ(T, x)，因此接口的第一层任务，是把“摩尔量”转换为“体积量”，把“相成分 x”转换为“相场加权成分 c”。

设摩尔体积为 V_m，则体积自由能密度为：

f_chem = (1/V_m) · Σ_φ h(η_φ) · G_m^φ(T, x_φ)

这里 h(η_φ) 是插值函数，x_φ 是各相内部成分。注意：G_m^φ 必须对 x_φ 可微，否则化学势无法解析求导，只能数值差分，会显著降低稳定性。

二、化学势与扩散势的接口

相场演化中真正驱动扩散的不是浓度梯度，而是化学势梯度。对每一相，化学势定义为：

μ_i^φ = ∂G_m^φ / ∂x_i^φ

在耦合框架中，通常使用“扩散势”形式：

M_ij = Σ_φ h(η_φ) · M_ij^φ

其中 M_ij^φ 是原子迁移率，可由示踪扩散系数或互扩散系数转换得到。扩散通量写为：

J_i = - Σ_j M_ij ∇μ_j

接口设计的关键是：μ_j 必须来自同一套CALPHAD自由能模型，而不是用理想溶液近似替代。否则在界面处会出现虚假的化学势跳跃，导致溶质在界面异常富集或贫化。

工程上常把 μ 写成“理想项 + 过剩项”：μ_i = RT ln x_i + μ_i^E。理想项解析可导，过剩项由CALPHAD的Redlich-Kister多项式或CALPHAD模型给出。这样既保留物理意义，又便于求导。

三、插值函数的选择

插值函数 h(η) 决定相场变量如何“混合”各相自由能。常用形式有：

h(η) = η²(3-2η)，即平滑阶跃函数，满足 h(0)=0、h(1)=1、h'(0)=h'(1)=0；
h(η) = η³(10-15η+6η²)，更高阶平滑，适合多相场；
h(η) = η，线性插值，简单但界面处导数不连续，易产生伪力。

插值函数必须满足两个条件：在纯相处只保留该相自由能；在界面处自由能单调过渡。若 h(η) 选择不当，界面能会随界面宽度漂移，甚至出现“双井”被抹平的现象。

四、相场变量与成分场的映射

多相场中，每个相有自己内部的成分 x_φ，而实际可观测的是总成分 c。二者关系为：

c = Σ_φ h(η_φ) · x_φ / Σ_φ h(η_φ)

当 Σ h(η_φ)=1 时，简化为 c = Σ h(η_φ) x_φ。这个约束必须在每个时间步满足，否则质量不守恒。常见做法是：先由扩散方程更新 c，再通过“相成分回代”求解各 x_φ，使总成分与各相成分自洽。

下面这张图给出从CALPHAD数据库到相场方程右端项的完整接口链路，重点标出自由能密度、化学势、扩散势、插值函数与相场变量映射之间的数据流向。

五、实际实现中的关键细节

在代码层面，耦合接口通常封装成一个“热力学求值器”，输入为 T、c、η 数组，输出为 f_chem、μ、∂μ/∂c 以及迁移率。以二元系为例，一个可用的接口函数签名可以写成：

def evaluate_thermo(T, c, eta):
    # 1. 由 eta 计算各相分数 h_phi
    # 2. 由 c 和 h_phi 反解各相成分 x_phi
    # 3. 调用 CALPHAD 模型计算 G_m^phi(T, x_phi)
    # 4. 计算 f_chem = (1/V_m) * sum(h_phi * G_m_phi)
    # 5. 计算 mu_i = dG_m/dx_i，并组装扩散势
    # 6. 返回 f_chem, mu, dmu_dc, M
    return f_chem, mu, dmu_dc, M

其中第2步“反解各相成分”是数值上最敏感的部分。常用做法是：对每个相，用牛顿迭代求解 ∂G_m^φ/∂x_φ = μ_ref，其中 μ_ref 由总成分约束确定。若迭代不收敛，可退化为理想溶液近似作为初值，再逐步修正。

注意：CALPHAD自由能通常以摩尔为基准，而相场方程以体积为基准。若 V_m 随成分变化，必须把 V_m(x) 也纳入求导链，否则化学势会缺少一项 ∂V_m/∂x 的贡献，导致界面处驱动力偏差。

六、一个可复现的示例参数

假设某二元合金在 T=800 K 下，基体相 α 与析出相 β 的摩尔Gibbs自由能分别为：

G_m^α = x_A G_A^α + x_B G_B^α + RT(x_A ln x_A + x_B ln x_B) + L^α x_A x_B
G_m^β = x_A G_A^β + x_B G_B^β + RT(x_A ln x_A + x_B ln x_B) + L^β x_A x_B

取 G_A^α=0、G_B^α=-8000 J/mol、G_A^β=-2000 J/mol、G_B^β=-4000 J/mol、L^α=5000 J/mol、L^β=8000 J/mol、V_m=7×10⁻⁶ m³/mol。用上述接口计算 f_chem 与 μ_B，可以观察到：在 η 从0过渡到1的界面区域内，μ_B 应保持连续；若出现明显跳跃，说明插值函数或相成分反解存在问题。

判断接口是否正确的三个实用检查：

纯相处：η=0 时 f_chem 应等于 G_m^α/V_m，η=1 时应等于 G_m^β/V_m；
界面处：μ_B 沿界面法向应单调或平滑变化，不应出现振荡；
守恒性：每步更新后，∫c dV 的变化应仅由边界通量决定，内部总质量守恒。

接口设计的目标不是“能算”，而是“算得稳、算得准、算得可复现”。把自由能密度、化学势、扩散势、插值函数与相场变量映射这五个环节逐一显式化，并在代码中保持一致的参照系与单位，才能让相场与CALPHAD的耦合真正服务于工程预测。