如何在FiPy的双温模型模拟中添加潜热项

方案稳定性对比 #

• 等效热容法（C_ph加高斯峰）：稳定性更优，该方案将潜热折算为热容的温度依赖项，属于瞬态项的原生系数，不需要额外引入动态开关的源项，避免了源项突变带来的迭代收敛问题。仅需注意合理设置相变温度区间宽度，区间过窄易出现数值振荡，过宽会偏离物理实际。
• 源项汇法：逻辑更贴合物理过程，直接在熔点位置扣除潜热，但需要引入温度阈值判断，源项的突变很容易导致迭代残差飙升甚至发散，需要配合更小的时间步长和更强的欠松弛，稳定性弱于等效热容法。

具体实现代码修改 #

1. 等效热容法（推荐） #

首先定义潜热相关参数，修改原C_ph的定义即可，其余代码无需调整：

# 新增潜热相关参数
T_melt = 2700  # 熔点，单位K
L = 360000     # 潜热，单位J
dT = 10        # 相变温度区间宽度，可根据精度需求在5~20K范围内调整

# 替换原C_ph定义
C_ph_base = 4.10446 - 3.886 * numerix.exp(-T_ph / 373.8)
# 高斯峰形式的等效热容增量，积分结果等于潜热L
gauss_peak = L / (dT * numerix.sqrt(2 * numerix.pi)) * numerix.exp(-(T_ph - T_melt)**2 / (2 * dT**2))
C_ph = C_ph_base + gauss_peak

如果需要更接近阶跃相变，可适当缩小dT，同时将时间步长的增长系数从1.01下调到1.005左右，避免迭代发散。

2. 源项汇法（可选） #

如果需要更严格贴合相变物理过程，可采用该方案，需修改方程定义和时间步循环逻辑：

# 新增潜热相关参数
T_melt = 2700
L = 360000
# 新增变量存储上一步声子温度，用于判断是否处于升温阶段
T_ph_old = CellVariable(mesh=mesh, value=300)

# 替换原eq1定义，添加潜热汇项
latent_source = - L * numerix.where((T_ph > T_melt - 1) & (T_ph < T_melt + 1) & (T_ph > T_ph_old), 1.0 / dt, 0)
eq1 = (TransientTerm(var=T_ph, coeff=C_ph) == DiffusionTerm(var=T_ph, coeff=k_ph) + ImplicitSourceTerm(var=T_e, coeff=G) - ImplicitSourceTerm(coeff=G, var=T_ph) + latent_source)

# 时间步循环新增T_ph_old更新逻辑，同时调小欠松弛系数提升稳定性
for step in range(steps):
    T_e.updateOld()
    T_ph.updateOld()
    T_ph_old.setValue(T_ph.value) # 保存上一步温度
    vi.plot()
    res = 1e100
    dt *= 1.005 # 下调时间步长增速
    count = 0
    while res > 1:
        res = eq.sweep(dt=dt, underRelaxation=0.3) # 欠松弛系数从0.5下调到0.3
        print(t, res)
    t.setValue(t + dt)

内容的提问来源于stack exchange，提问作者Salah