3. 二元体系规则溶液模型参数优化与相图计算

规则溶液模型与参数优化思路

二元规则溶液模型把液相和固相（或两个固溶体相）的摩尔Gibbs能写成理想混合项与一个成分对称的过剩项之和。对A-B二元体系，以1 mol原子为基准：

G_m^φ = x_A·G_A^φ + x_B·G_B^φ
        + R·T·(x_A·ln x_A + x_B·ln x_B)
        + x_A·x_B·L^φ

其中 L^φ = a^φ + b^φ·T   （规则溶液：L 与成分无关）

纯组元项 G_A^φ、G_B^φ 来自前一步的纯组元Gibbs能函数，本例视为已知固定值；待优化量是各相的相互作用参数 a^φ 与 b^φ。相界位置由两相Gibbs能公切线条件决定，因此相图对 L 的取值非常敏感：L 越负，液相越稳定、液相线越低；L 越正，越容易在低温出现互溶间隙。

用实验相界数据反推参数

设某A-B体系液相为规则溶液，固相（端际固溶体）近似为理想溶液（L^s≈0）。已知实验测得的液相线温度点如下（示例数据，用于演示优化流程）：

优化目标函数取液相线温度残差平方和：

F(a_L, b_L) = Σ_i [ T_calc(x_i; a_L, b_L) − T_exp,i ]²

T_calc 由公切线条件数值求解：
  μ_A^L(x_L) = μ_A^S(x_S)
  μ_B^L(x_L) = μ_B^S(x_S)
  联立解出 x_L, x_S, T

由于只有两个待定参数，可用网格搜索加局部细化：先以 a_L∈[−40000, 0] J/mol、b_L∈[−20, 20] J/(mol·K) 粗扫，再在最优邻域用单纯形法细化。若残差仍大，说明规则溶液模型不足，需要引入亚规则项（下一例处理）。

参数—相界关系与优化回路

可执行的优化脚本（Python）

下面给出一个最小可运行示例，使用 scipy 的 minimize 对 a_L、b_L 进行优化。纯组元Gibbs能函数用简化线性形式代替，实际使用时替换为前一步拟合结果。

import numpy as np
from scipy.optimize import minimize

R = 8.314  # J/(mol·K)

# 示例：纯组元熔点与熔化熵（简化，仅用于演示）
TmA, TmB = 1700.0, 1200.0
SmA, SmB = 10.0, 8.0   # J/(mol·K)

def G_pure(T, Tm, Sm, phase):
    """简化纯组元Gibbs能：液相为0参考，固相为熔化驱动"""
    if phase == 'L':
        return 0.0
    return Sm * (Tm - T)  # 固相相对液相的Gibbs能差

def G_mix(xB, T, L):
    """规则溶液摩尔Gibbs能（相对纯液相参考）"""
    xA = 1.0 - xB
    # 理想混合项
    G_ideal = R * T * (xA * np.log(xA + 1e-30) + xB * np.log(xB + 1e-30))
    # 过剩项
    G_excess = xA * xB * L
    return G_ideal + G_excess

def G_phase(xB, T, L, phase):
    """某相的摩尔Gibbs能（相对纯液相参考）"""
    xA = 1.0 - xB
    G = xA * G_pure(T, TmA, SmA, phase) + xB * G_pure(T, TmB, SmB, phase)
    G += G_mix(xB, T, L)
    return G

def chemical_potentials(xB, T, L, phase):
    """数值求化学势：μ_A = G - xB*dG/dxB, μ_B = G + xA*dG/dxB"""
    dx = 1e-6
    G0 = G_phase(xB, T, L, phase)
    G1 = G_phase(xB + dx, T, L, phase)
    dGdx = (G1 - G0) / dx
    muA = G0 - xB * dGdx
    muB = G0 + (1.0 - xB) * dGdx
    return muA, muB

def solve_tie_line(T, L_L, L_S):
    """给定温度，求液/固两相平衡成分"""
    from scipy.optimize import fsolve
    def eqs(vars):
        xL, xS = vars
        muA_L, muB_L = chemical_potentials(xL, T, L_L, 'L')
        muA_S, muB_S = chemical_potentials(xS, T, L_S, 'S')
        return [muA_L - muA_S, muB_L - muB_S]
    try:
        sol = fsolve(eqs, [0.5, 0.5], full_output=False)
        if 0 < sol[0] < 1 and 0 < sol[1] < 1:
            return sol
    except Exception:
        pass
    return None

def liquidus_T(xB_exp, a_L, b_L):
    """对给定成分，求液相线温度（固相L_S=0）"""
    from scipy.optimize import brentq
    def f(T):
        L_L = a_L + b_L * T
        res = solve_tie_line(T, L_L, 0.0)
        if res is None:
            return 1e6
        xL, xS = res
        return xL - xB_exp
    try:
        return brentq(f, 800, 2000)
    except Exception:
        return np.nan

# 实验数据
x_exp = np.array([0.20, 0.40, 0.60, 0.80])
T_exp = np.array([1580.0, 1480.0, 1400.0, 1340.0])

def objective(params):
    a_L, b_L = params
    err = 0.0
    for x, T in zip(x_exp, T_exp):
        Tc = liquidus_T(x, a_L, b_L)
        if np.isnan(Tc):
            return 1e9
        err += (Tc - T) ** 2
    return err

# 粗扫 + 局部优化
best = None
for a in np.linspace(-40000, 0, 9):
    for b in np.linspace(-20, 20, 9):
        val = objective([a, b])
        if best is None or val < best[0]:
            best = (val, a, b)

res = minimize(objective, [best[1], best[2]], method='Nelder-Mead',
               options={'xatol': 1e-3, 'fatol': 1e-3, 'maxiter': 2000})
print("最优 a_L = %.1f J/mol, b_L = %.3f J/(mol·K)" % (res.x[0], res.x[1]))
print("残差平方和 F = %.1f" % res.fun)

结果判读与验证

• 残差量级：若 F 对应的均方根偏差小于实验温度不确定度（通常 5–15 K），说明规则溶液模型在该体系可接受。
• 参数符号：a_L 为负表示液相互吸引、液相线被压低；a_L 为正且较大时，低温可能出现液相互溶间隙，需检查是否与实验一致。
• 相界形状：规则溶液模型的液相线与固相线在成分-温度图上呈对称拱形。若实验相界明显不对称，说明需要亚规则项（L 随成分变化）。
• 外推检查：用优化后的参数计算 x_B=0.1 和 0.9 处的相界，与实验或文献对比，确认参数没有过拟合中间成分点。

交付物：参数优化记录表