3. 平衡凝固的核心方程与数值实现

平衡凝固是理解合金凝固行为的基准模型。它的核心假设是：固相与液相在每一温度下都达到完全扩散平衡，固相成分和液相成分始终由杠杆定律控制。本节从方程出发，逐步展开温度—固相分数关系、平衡凝固路径的逐步计算流程，以及数值实现中的收敛判据。

一、杠杆定律：平衡凝固的成分分配

对于二元合金，设合金总成分为 C₀，在温度 T 下液相成分为 C_L，固相成分为 C_S，固相分数为 f_S，液相分数为 f_L = 1 − f_S。质量守恒要求：

C₀ = f_S · C_S + (1 − f_S) · C_L

整理得杠杆定律：

f_S = (C₀ − C_L) / (C_S − C_L)
f_L = (C_S − C₀) / (C_S − C_L)

这里 C_L 和 C_S 由相图给出。在平衡条件下，两者都是温度 T 的函数：C_L = C_L(T)，C_S = C_S(T)。因此 f_S 也随温度变化，这就是温度—固相分数关系的来源。

需要特别说明参照约定：本文所有成分均以溶质质量分数表示，C_S 与 C_L 的差值方向由相图决定。对于 k₀ < 1 的合金，C_S < C_L，固相中溶质含量低于液相；对于 k₀ > 1 的合金，方向相反。全文及图示均采用这一约定。

二、温度—固相分数关系

以最简单的线性相图为例。设液相线斜率为 m_L，固相线斜率为 m_S，纯溶剂熔点为 T_m，则：

C_L(T) = (T − T_m) / m_L
C_S(T) = (T − T_m) / m_S

代入杠杆定律：

f_S(T) = [C₀ − (T − T_m)/m_L] / [(T − T_m)/m_S − (T − T_m)/m_L]

化简后可得 f_S 与 T 的显式关系。当 T = T_L（液相线温度）时，f_S = 0；当 T = T_S（固相线温度）时，f_S = 1。中间温度下，f_S 随温度降低而单调增加。

对于实际合金体系，C_L(T) 和 C_S(T) 往往不是线性关系，需要借助热力学数据库或多项式拟合。此时 f_S(T) 只能通过数值方法求解。

三、平衡凝固路径的逐步计算流程

数值实现平衡凝固路径，通常采用温度步进法。以下是可执行的逐步流程：

1. 输入参数：合金成分 C₀、相图数据（C_L(T) 和 C_S(T) 的表达式或插值表）、起始温度 T_start（通常为液相线温度以上）、终止温度 T_end（通常为固相线温度以下）、温度步长 ΔT。
2. 初始化：令 T = T_start，f_S = 0，记录初始状态。
3. 温度递减：T = T − ΔT。
4. 计算平衡成分：根据当前 T，由相图数据得到 C_L(T) 和 C_S(T)。
5. 计算固相分数：f_S = (C₀ − C_L) / (C_S − C_L)。
6. 收敛判据检查：若 |f_S − f_S_prev| < ε_f 且 |T − T_prev| < ε_T，则认为该步收敛；否则减小 ΔT 并重新计算。
7. 记录与更新：保存 (T, f_S, C_L, C_S)，更新 f_S_prev = f_S，T_prev = T。
8. 终止判断：若 f_S ≥ 1 或 T ≤ T_end，停止；否则返回步骤 3。

其中 ε_f 和 ε_T 是收敛容差，典型取值分别为 1×10⁻⁴ 和 0.01 K。对于非线性相图，步骤 4 可能需要内层迭代（如牛顿法）来求解 C_L 和 C_S 的平衡关系。

四、数值实现中的关键细节

步长选择：ΔT 过大会导致 f_S 跳跃，错过固相线附近的快速变化；ΔT 过小则计算量增加。建议在液相线附近用较大步长（如 1 K），在固相线附近自动减小到 0.1 K 或更小。

收敛判据的物理意义：ε_f 控制固相分数的变化量，ε_T 控制温度变化量。两者同时满足才认为该温度步稳定。若仅满足其一，说明该步内 f_S 对 T 过于敏感，需要细分。

边界处理：当 f_S 接近 0 或 1 时，分母 C_S − C_L 可能趋近于零，导致数值不稳定。此时应改用解析极限或插值外推，而不是直接相除。

以下是一个简化的 Python 伪代码示例，展示核心循环：

def equilibrium_path(C0, CL_func, CS_func, T_start, T_end, dT, eps_f, eps_T):
    T = T_start
    fS_prev = 0.0
    T_prev = T_start
    path = []
    while T > T_end:
        T -= dT
        CL = CL_func(T)
        CS = CS_func(T)
        if abs(CS - CL) < 1e-12:
            break
        fS = (C0 - CL) / (CS - CL)
        if abs(fS - fS_prev) < eps_f and abs(T - T_prev) < eps_T:
            path.append((T, fS, CL, CS))
            fS_prev = fS
            T_prev = T
        else:
            dT *= 0.5  # 细分步长
            T += dT   # 回退
    return path

该伪代码体现了“计算—判据—细分”的闭环逻辑。实际工程中，还需处理相图数据的插值、多组元体系的耦合求解，以及固相分数超过 1 后的截断。

五、结果判读与常见问题

计算完成后，应检查 f_S 是否从 0 单调增加到 1，C_L 和 C_S 是否随温度降低而分别沿液相线和固相线移动。若 f_S 出现非单调或超过 1 的情况，通常说明步长过大或相图数据在固相线附近不连续。

另一个常见问题是收敛判据设置过松，导致 f_S 在固相线附近跳跃。此时应减小 ε_f 或启用自适应步长。对于多组元合金，还需注意各组元杠杆定律的耦合——每个组元独立满足质量守恒，但固相分数必须一致，这要求联立求解。

掌握平衡凝固的数值实现，是后续学习 Scheil 计算和实际凝固路径分析的基础。核心在于：杠杆定律提供成分分配关系，温度步进法提供数值求解框架，收敛判据保证计算稳定性。