4. 二元体系亚规则溶液模型对液相互溶间隙的再现

为什么规则溶液模型画不出液相互溶间隙

规则溶液模型中，液相摩尔Gibbs能写作：

G^L_m = x_AG_A^0,L + x_BG_B^0,L + RT(x_Alnx_A + x_Blnx_B) + x_Ax_BL₀

其中相互作用参数 L₀ 为常数。当 L₀ 足够大且为正时，G_m–x 曲线出现双谷，但对称的规则溶液模型只能给出对称的互溶间隙，且临界点固定在 x_B = 0.5。实际体系（如Al–Zn、Cu–Ni、Fe–Cr等）的液相互溶间隙往往不对称，临界成分偏离0.5，必须引入成分依赖的亚规则项。

亚规则溶液模型与参数设置

亚规则溶液模型将相互作用参数展开为成分的线性函数：

L(x_B) = L₀ + L₁(x_A − x_B) = L₀ + L₁(1 − 2x_B)

此时过剩Gibbs能 G^E = x_Ax_B[L₀ + L₁(x_A − x_B)]。化学势差（即公切线条件）为：

Δμ_B = G_B^0,L − G_A^0,L + RT ln(x_B/x_A) + G^E + (1 − x_B)∂G^E/∂x_B

互溶间隙边界由 Δμ_B(x_B^α) = Δμ_B(x_B^β) 且 Δμ_A(x_B^α) = Δμ_A(x_B^β) 确定。临界点满足 ∂²G_m/∂x_B² = 0 且 ∂³G_m/∂x_B³ = 0。

互溶间隙与临界点的计算流程

以下Python代码使用亚规则模型，在给定温度下求解公切线方程，得到互溶间隙两端的成分，并扫描温度找到临界点。

import numpy as np
from scipy.optimize import fsolve, brentq

R = 8.314
L0 = 18000.0   # J/mol
L1 = -6000.0   # J/mol

def dG_dx(xB, T):
    """亚规则模型 dGm/dxB (除以RT归一化前的原始值)"""
    xA = 1.0 - xB
    # 理想混合项导数: RT*ln(xB/xA)
    ideal = R * T * np.log(xB / xA)
    # 过剩项 GE = xA*xB*(L0 + L1*(xA-xB))
    # dGE/dxB = (1-2xB)*(L0+L1*(1-2xB)) + xA*xB*(-2*L1)
    s = 1.0 - 2.0 * xB
    dGE = s * (L0 + L1 * s) + xA * xB * (-2.0 * L1)
    return ideal + dGE

def common_tangent(x, T):
    """公切线条件: dG/dxB 在 x[0] 和 x[1] 处相等, 且截距相等"""
    x1, x2 = x
    # 条件1: 两处斜率相等
    eq1 = dG_dx(x1, T) - dG_dx(x2, T)
    # 条件2: 两处切线截距相等 (用 dG/dxB 和 G 的关系)
    # G(x) = RT[xA ln xA + xB ln xB] + xA*xB*(L0+L1*(xA-xB))
    def Gm(xB):
        xA = 1.0 - xB
        ideal = R * T * (xA * np.log(xA) + xB * np.log(xB))
        s = xA - xB
        ge = xA * xB * (L0 + L1 * s)
        return ideal + ge
    # 截距: G(x) - x * dG/dx 在两点相等
    b1 = Gm(x1) - x1 * dG_dx(x1, T)
    b2 = Gm(x2) - x2 * dG_dx(x2, T)
    eq2 = b1 - b2
    return [eq1, eq2]

# 在 T = 700 K 求解互溶间隙
T_test = 700.0
x_guess = [0.15, 0.85]
sol = fsolve(common_tangent, x_guess, args=(T_test,))
x_alpha, x_beta = sorted(sol)
print(f"T = {T_test} K: 互溶间隙 = {x_alpha:.4f} ~ {x_beta:.4f}")

# 扫描温度找临界点: 两相成分合并
def gap_width(T):
    try:
        s = fsolve(common_tangent, [0.15, 0.85], args=(T,), full_output=True)
        xa, xb = sorted(s[0])
        if s[2] == 1 and 0 < xa < xb < 1:
            return xb - xa
    except Exception:
        pass
    return None

T_scan = np.linspace(600, 900, 61)
widths = []
for T in T_scan:
    w = gap_width(T)
    widths.append(w if w is not None else np.nan)

# 找宽度趋近0的温度作为临界温度估计
valid = [(T, w) for T, w in zip(T_scan, widths) if w is not None and w > 1e-4]
if valid:
    T_c_est = valid[-1][0]
    print(f"临界温度估计: T_c ≈ {T_c_est:.1f} K")
    # 在临界温度附近细化
    for T in np.linspace(T_c_est - 5, T_c_est + 5, 21):
        w = gap_width(T)
        if w is not None:
            print(f"  T = {T:.1f} K, 间隙宽度 = {w:.5f}")

参数对临界点位置的影响

通过改变 L₀ 和 L₁，可以观察临界温度和临界成分的变化规律：

• 增大 L₀：互溶间隙扩大，临界温度升高。规则溶液模型下 T_c = L₀/(2R)，亚规则模型下临界温度近似为 T_c ≈ (L₀ + L₁/3)/(2R)（小 L₁ 近似）。
• 引入 L₁ ≠ 0：临界成分偏离0.5。当 L₁ < 0 时，临界成分向 x_B < 0.5 方向移动；当 L₁ > 0 时向 x_B > 0.5 方向移动。
• 临界成分的解析近似：在临界点附近，x_c ≈ 0.5 − L₁/(6L₀ + 2L₁)，可用于快速估算。

结果判读与验证

计算完成后，需要检查以下要点：

1. 间隙宽度随温度单调减小：从低温到临界温度，两相成分应逐渐靠近，宽度趋近于零。
2. 临界点处二阶和三阶导数同时为零：可用数值微分验证 ∂²G_m/∂x_B² 和 ∂³G_m/∂x_B³ 在临界点附近变号。
3. 公切线条件残差：将求得的 x_α 和 x_β 代回公切线方程，残差应小于 1×10⁻⁶。
4. 与实验相图对比：若已知实验互溶间隙数据，检查计算边界是否落在实验误差范围内。

实操检查清单

检查项	通过条件	常见问题
公切线残差	< 1×10−6	初值选取不当导致收敛到平凡解
间隙宽度随T变化	单调递减至0	参数过大导致数值溢出
临界点导数条件	∂²G/∂x² 和 ∂³G/∂x³ 同时变号	仅检查二阶导数，遗漏三阶条件
临界成分偏移方向	与L1符号一致	符号约定错误
与实验对比	边界落在实验误差带内	未考虑固相或其它相的竞争