对流换热仿真-主题073_相场法模拟相变-相场法模拟相变

kkchenjj

49人浏览 · 2026-04-20 21:24:48

kkchenjj · 2026-04-20 21:24:48 发布

主题073：相场法模拟相变

一、引言

相变是自然界和工程领域中广泛存在的物理现象。从水的结冰到金属的凝固，从聚合物的结晶到合金的沉淀，相变过程伴随着潜热的释放或吸收、界面的移动以及复杂的微观结构演化。准确模拟相变过程对于材料科学、能源工程、气候科学等众多领域具有重要意义。

传统的相变模拟方法，如Stefan问题的前追踪法（Front Tracking），需要显式地追踪相界面位置，这在处理复杂界面形貌（如枝晶生长、多晶交汇）时会遇到巨大困难。相场法（Phase Field Method）作为一种基于连续介质理论的界面追踪方法，通过引入一个连续的序参量（相场变量）来描述不同相之间的过渡，避免了显式追踪界面的复杂性，成为模拟相变过程的有力工具。

本主题将系统介绍相场法的基本理论、数学模型和数值方法，重点探讨Allen-Cahn方程和Cahn-Hilliard方程在相变模拟中的应用，并通过Python实例展示如何实现凝固、熔化、晶粒生长等过程的数值模拟。
在这里插入图片描述

二、相场法基础理论

2.1 相场法的基本概念

2.1.1 什么是相场法

相场法是一种基于热力学原理的界面追踪方法。其核心思想是引入一个连续的相场变量 $ϕ(x,t)\phi(\mathbf{x}, t)$ 来区分不同的相：

$ϕ=1\phi = 1$ 表示一相（如固相）
$ϕ=−1\phi = -1$ 表示另一相（如液相）
$\phi < 1$ 表示界面过渡区域

与传统的尖锐界面模型不同，相场法将界面视为具有一定厚度的扩散区域。这种"扩散界面"的近似使得相场法具有以下优势：

无需显式追踪界面：界面位置由相场变量的等值面隐式定义
自动处理拓扑变化：界面的合并、断裂等复杂行为可以自然发生
易于处理复杂几何：可以模拟枝晶、多晶等复杂微观结构
热力学一致性：基于自由能泛函变分推导，保证物理守恒

2.1.2 自由能泛函

相场法的理论基础是Ginzburg-Landau理论。系统的总自由能泛函可以写为：

$F[ϕ,T]=∫Ω[f(ϕ,T)+ϵ22∣∇ϕ∣2]dΩF[\phi, T] = \int_\Omega \left[ f(\phi, T) + \frac{\epsilon^2}{2} |\nabla \phi|^2 \right] d\Omega$

其中：

$f(ϕ,T)f(\phi, T)$ 是体自由能密度，描述两相的热力学性质
$ϵ22∣∇ϕ∣2\frac{\epsilon^2}{2} |\nabla \phi|^2$ 是梯度能项，惩罚相场变量的快速变化，控制界面厚度
$ϵ\epsilon$ 是界面厚度参数

体自由能密度通常采用双阱势函数形式：

$f(ϕ,T)=14(ϕ2−1)2+λg(ϕ)ΔTf(\phi, T) = \frac{1}{4}(\phi^2 - 1)^2 + \lambda g(\phi) \Delta T$

其中：

$14(ϕ2−1)2\frac{1}{4}(\phi^2 - 1)^2$ 是双阱势，在 $ϕ=±1\phi = \pm 1$ 处有两个稳定极小值
$g(ϕ)g(\phi)$ 是插值函数，耦合温度场
$λ\lambda$ 是耦合系数
$ΔT=T−Tm\Delta T = T - T_m$ 是过冷度， $T_m$ 是熔点

2.1.3 界面厚度与表面张力

相场模型中的界面厚度 $W$ 与参数 $ϵ\epsilon$ 的关系为：

$2\sqrt{2}\epsilon$

界面能（表面张力） $σ\sigma$ 与相场参数的关系为：

$σ=223ϵ\sigma = \frac{2\sqrt{2}}{3} \epsilon$

在实际计算中，需要选择适当的界面厚度：

界面太厚：无法捕捉界面细节
界面太薄：需要极细网格，计算成本高

通常选择界面厚度远小于特征长度但足以分辨界面结构。

2.2 Allen-Cahn方程

2.2.1 方程推导

Allen-Cahn方程描述非守恒序参量的演化，适用于凝固、熔化等相变过程。它可以从自由能泛函的变分导出：

$∂ϕ∂t=−MϕδFδϕ\frac{\partial \phi}{\partial t} = -M_\phi \frac{\delta F}{\delta \phi}$

其中 $MϕM_\phi$ 是迁移率（mobility），控制相场演化的速率。

计算变分：

$δFδϕ=∂f∂ϕ−ϵ2∇2ϕ\frac{\delta F}{\delta \phi} = \frac{\partial f}{\partial \phi} - \epsilon^2 \nabla^2 \phi$

对于双阱势 $\frac{1}{4}(\phi^2 - 1)^2$ ：

$∂f∂ϕ=ϕ(ϕ2−1)=ϕ3−ϕ\frac{\partial f}{\partial \phi} = \phi(\phi^2 - 1) = \phi^3 - \phi$

因此，Allen-Cahn方程为：

$∂ϕ∂t=Mϕ(ϵ2∇2ϕ−ϕ3+ϕ)\frac{\partial \phi}{\partial t} = M_\phi \left( \epsilon^2 \nabla^2 \phi - \phi^3 + \phi \right)$

2.2.2 物理意义

Allen-Cahn方程可以重写为：

$∂ϕ∂t=Mϕ(ϵ2∇2ϕ−f′(ϕ))\frac{\partial \phi}{\partial t} = M_\phi \left( \epsilon^2 \nabla^2 \phi - f'(\phi) \right)$

其中：

$ϵ2∇2ϕ\epsilon^2 \nabla^2 \phi$ 是扩散项，使界面平滑
$−f′(ϕ)=ϕ−ϕ3-f'(\phi) = \phi - \phi^3$ 是反应项，驱使系统趋向双阱势的极小值

方程的稳态解（一维情况）为：

$ϕ(x)=tanh⁡(x2ϵ)\phi(x) = \tanh\left(\frac{x}{\sqrt{2}\epsilon}\right)$

这是一个光滑的界面轮廓，从 $ϕ=−1\phi = -1$ 过渡到 $ϕ=1\phi = 1$ ，界面中心位于 $x = 0$ 。

2.2.3 与温度的耦合

在凝固问题中，需要将相场方程与温度场耦合。温度场方程为：

$∂T∂t=α∇2T+Lcp∂ϕ∂t\frac{\partial T}{\partial t} = \alpha \nabla^2 T + \frac{L}{c_p} \frac{\partial \phi}{\partial t}$

其中：

$α\alpha$ 是热扩散系数
$L$ 是潜热
$c_p$ 是比热容
$Lcp∂ϕ∂t\frac{L}{c_p} \frac{\partial \phi}{\partial t}$ 是相变热源项

相场方程也需要包含温度驱动项：

$∂ϕ∂t=Mϕ(ϵ2∇2ϕ−ϕ3+ϕ)+Mϕλ(T−Tm)g′(ϕ)\frac{\partial \phi}{\partial t} = M_\phi \left( \epsilon^2 \nabla^2 \phi - \phi^3 + \phi \right) + M_\phi \lambda (T - T_m) g'(\phi)$

其中 $g(ϕ)=ϕ−ϕ33g(\phi) = \phi - \frac{\phi^3}{3}$ 是插值函数。

2.3 Cahn-Hilliard方程

2.3.1 守恒型相场方程

Cahn-Hilliard方程描述守恒序参量的演化，适用于成分分离、晶粒粗化等过程。它保证序参量的总量守恒：

$∂ϕ∂t=∇⋅(M∇δFδϕ)\frac{\partial \phi}{\partial t} = \nabla \cdot \left( M \nabla \frac{\delta F}{\delta \phi} \right)$

展开后得到：

$∂ϕ∂t=∇⋅(M∇(f′(ϕ)−ϵ2∇2ϕ))\frac{\partial \phi}{\partial t} = \nabla \cdot \left( M \nabla (f'(\phi) - \epsilon^2 \nabla^2 \phi) \right)$

或写成：

$∂ϕ∂t=M(f′′(ϕ)∇2ϕ−ϵ2∇4ϕ)\frac{\partial \phi}{\partial t} = M \left( f''(\phi) \nabla^2 \phi - \epsilon^2 \nabla^4 \phi \right)$

2.3.2 与Allen-Cahn方程的区别

特征	Allen-Cahn方程	Cahn-Hilliard方程
守恒性	非守恒	守恒
方程阶数	二阶	四阶
主要应用	凝固/熔化	相分离/晶粒粗化
界面速度	与曲率无关	与曲率相关
数值稳定性	较易处理	需要特殊处理

2.3.3 应用实例

旋节分解（Spinodal Decomposition）：
当二元合金从高温淬火到混溶隙以下时，会发生成分分离。Cahn-Hilliard方程可以描述这种自发的相分离过程。

晶粒粗化（Grain Coarsening）：
多晶材料中，大晶粒会吞噬小晶粒以降低总界面能。使用多序参量相场模型可以模拟这一过程。

2.4 各向异性与枝晶生长

2.4.1 界面能各向异性

实际晶体材料的界面能通常具有各向异性。在相场模型中，可以通过使界面厚度参数 $ϵ\epsilon$ 或迁移率 $M$ 依赖于界面法向方向来引入各向异性。

界面法向量为：

$n=−∇ϕ∣∇ϕ∣\mathbf{n} = -\frac{\nabla \phi}{|\nabla \phi|}$

对于二维四重对称的枝晶生长，界面能可以表示为：

$ϵ(θ)=ϵ0(1+δ4cos⁡(4θ))\epsilon(\theta) = \epsilon_0 (1 + \delta_4 \cos(4\theta))$

其中 $θ\theta$ 是界面法向与晶体学方向的夹角， $δ4\delta_4$ 是各向异性强度。

2.4.2 枝晶生长模型

枝晶生长是凝固过程中常见的现象。相场法可以自然地模拟枝晶的形成和演化：

噪声引入：在相场方程中添加热噪声项，触发界面失稳
各向异性驱动：界面能各向异性导致优先生长方向
热扩散限制：潜热释放改变局部温度场，反馈影响生长

枝晶生长的相场方程：

$τ∂ϕ∂t=W2∇2ϕ+ϕ(1−ϕ)(ϕ−12+m(T))+各向异性项+噪声项\tau \frac{\partial \phi}{\partial t} = W^2 \nabla^2 \phi + \phi(1-\phi)\left(\phi - \frac{1}{2} + m(T)\right) + \text{各向异性项} + \text{噪声项}$

其中 $m (T)$ 是过冷度函数。

三、相场法数值方法

3.1 空间离散方法

3.1.1 有限差分法

有限差分法是相场法最常用的空间离散方法。对于二维问题，使用均匀网格 $(i, j)$ ，空间步长为 $Δx\Delta x$ 。

拉普拉斯算子（五点格式）：

$∇2ϕi,j=ϕi+1,j+ϕi−1,j+ϕi,j+1+ϕi,j−1−4ϕi,jΔx2\nabla^2 \phi_{i,j} = \frac{\phi_{i+1,j} + \phi_{i-1,j} + \phi_{i,j+1} + \phi_{i,j-1} - 4\phi_{i,j}}{\Delta x^2}$

双调和算子（Cahn-Hilliard方程需要）：

$∇4ϕi,j=∇2ϕi+1,j+∇2ϕi−1,j+∇2ϕi,j+1+∇2ϕi,j−1−4∇2ϕi,jΔx2\nabla^4 \phi_{i,j} = \frac{\nabla^2 \phi_{i+1,j} + \nabla^2 \phi_{i-1,j} + \nabla^2 \phi_{i,j+1} + \nabla^2 \phi_{i,j-1} - 4\nabla^2 \phi_{i,j}}{\Delta x^2}$

3.1.2 谱方法

对于周期性边界条件的问题，谱方法（特别是傅里叶谱方法）具有高精度优势。

在傅里叶空间，微分算子变为代数运算：

$F(∇2ϕ)=−∣k∣2ϕ^\mathcal{F}(\nabla^2 \phi) = -|\mathbf{k}|^2 \hat{\phi}$

其中 $k\mathbf{k}$ 是波矢量， $ϕ^\hat{\phi}$ 是傅里叶变换。

Cahn-Hilliard方程在傅里叶空间的半隐式格式：

$ϕ^n+1−ϕ^nΔt=−M∣k∣2(f′(ϕ^n)+ϵ2∣k∣2ϕ^n+1)\frac{\hat{\phi}^{n+1} - \hat{\phi}^n}{\Delta t} = -M |\mathbf{k}|^2 \left( f'(\hat{\phi}^n) + \epsilon^2 |\mathbf{k}|^2 \hat{\phi}^{n+1} \right)$

3.2 时间积分方法

3.2.1 显式欧拉法

最简单的显式格式：

$ϕn+1=ϕn+Δt⋅RHS(ϕn)\phi^{n+1} = \phi^n + \Delta t \cdot \text{RHS}(\phi^n)$

稳定性限制：

Allen-Cahn方程： $Δt<CΔx2ϵ2\Delta t < C \frac{\Delta x^2}{\epsilon^2}$
Cahn-Hilliard方程： $Δt<CΔx4ϵ2\Delta t < C \frac{\Delta x^4}{\epsilon^2}$

Cahn-Hilliard方程的时间步长限制非常严格，通常需要隐式或半隐式方法。

3.2.2 半隐式格式

为了提高稳定性，可以将线性项隐式处理，非线性项显式处理。

Allen-Cahn方程的半隐式格式：

$ϕn+1−ϕnΔt=Mϕ(ϵ2∇2ϕn+1−(ϕn)3+ϕn)\frac{\phi^{n+1} - \phi^n}{\Delta t} = M_\phi \left( \epsilon^2 \nabla^2 \phi^{n+1} - (\phi^n)^3 + \phi^n \right)$

这可以重写为：

$M_\phi \epsilon^2 \Delta t \nabla^2) \phi^{n+1} = \phi^n + M_\phi \Delta t \left( - (\phi^n)^3 + \phi^n \right)$

左边形成Helmholtz方程，可以用快速算法求解。

3.2.3 自适应时间步长

相变过程中，界面运动速度会随时间变化。自适应时间步长可以提高计算效率：

$Δtn+1=min⁡(Δtmax⁡,ϵmax⁡∣∂ϕ/∂t∣)\Delta t_{n+1} = \min\left( \Delta t_{\max}, \frac{\epsilon}{\max|\partial \phi / \partial t|} \right)$

3.3 边界条件处理

3.3.1 绝热边界

对于绝热边界，使用Neumann边界条件：

$∂ϕ∂n=0,∂T∂n=0\frac{\partial \phi}{\partial n} = 0, \quad \frac{\partial T}{\partial n} = 0$

在有限差分中，使用镜像法：

$ϕ−1,j=ϕ1,j,ϕN+1,j=ϕN−1,j\phi_{-1,j} = \phi_{1,j}, \quad \phi_{N+1,j} = \phi_{N-1,j}$

3.3.2 周期性边界

周期性边界条件：

$ϕ(0,y)=ϕ(L,y),ϕ(x,0)=ϕ(x,L)\phi(0, y) = \phi(L, y), \quad \phi(x, 0) = \phi(x, L)$

适合模拟无限大系统或周期性结构。

3.3.3 固定温度边界

Dirichlet边界条件：

$T = T_{wall}$

直接设置边界格点的温度值。

四、Python仿真实例

4.1 案例1：一维相界面演化

本案例模拟一维情况下相界面的平衡轮廓和弛豫过程。

"""
案例1：一维相界面演化
模拟相场变量的平衡轮廓和弛豫过程
"""
import numpy as np
import matplotlib.pyplot as plt

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

# 参数设置
Nx = 200  # 网格数
L = 10.0  # 计算域长度
dx = L / Nx  # 空间步长
epsilon = 0.5  # 界面厚度参数
M_phi = 1.0  # 迁移率
dt = 0.001  # 时间步长
Nt = 5000  # 时间步数

# 创建网格
x = np.linspace(-L/2, L/2, Nx)

# 初始条件：阶梯函数
phi = np.ones(Nx)
phi[x < 0] = -1

# 添加小扰动
phi += 0.1 * np.random.randn(Nx) * (np.abs(x) < 2)

# 存储结果
phi_history = [phi.copy()]
t_history = [0]

# 时间推进
for n in range(Nt):
    # 计算拉普拉斯算子（中心差分）
    laplacian = np.zeros_like(phi)
    laplacian[1:-1] = (phi[2:] - 2*phi[1:-1] + phi[:-2]) / dx**2
    laplacian[0] = laplacian[1]  # Neumann边界
    laplacian[-1] = laplacian[-2]
    
    # Allen-Cahn方程
    dphi_dt = M_phi * (epsilon**2 * laplacian - phi**3 + phi)
    phi = phi + dt * dphi_dt
    
    # 保存结果
    if n % 500 == 0:
        phi_history.append(phi.copy())
        t_history.append((n+1) * dt)

# 理论平衡解
phi_theory = np.tanh(x / (np.sqrt(2) * epsilon))

# 绘制结果
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 图1：相场演化
ax = axes[0]
for i, (phi_snap, t) in enumerate(zip(phi_history[::2], t_history[::2])):
    ax.plot(x, phi_snap, label=f't={t:.1f}')
ax.plot(x, phi_theory, 'k--', linewidth=2, label='理论平衡解')
ax.set_xlabel('位置 x')
ax.set_ylabel('相场变量 φ')
ax.set_title('相界面弛豫过程')
ax.legend()
ax.grid(True)
ax.set_ylim(-1.2, 1.2)

# 图2：自由能演化
def free_energy(phi, dx, epsilon):
    """计算自由能"""
    # 梯度能
    grad_phi = np.gradient(phi, dx)
    E_grad = 0.5 * epsilon**2 * np.sum(grad_phi**2) * dx
    
    # 体自由能
    E_bulk = np.sum(0.25 * (phi**2 - 1)**2) * dx
    
    return E_grad + E_bulk

E_history = [free_energy(phi, dx, epsilon) for phi in phi_history]
ax = axes[1]
ax.plot(t_history, E_history, 'b-', linewidth=2)
ax.set_xlabel('时间 t')
ax.set_ylabel('自由能 F')
ax.set_title('自由能随时间演化')
ax.grid(True)

plt.tight_layout()
plt.savefig('d:\\文档\\500仿真领域\\工程仿真\\对流换热仿真\\主题073_相场法模拟相变\\仿真结果\\case1_1d_interface.png', dpi=150)
plt.close()

print("案例1完成！")

4.2 案例2：二维凝固过程

本案例模拟二维空间中的凝固过程，展示平面界面的失稳和胞状结构的形成。

"""
案例2：二维凝固过程
模拟过冷熔体中的凝固过程
"""
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation

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

# 参数设置
Nx, Ny = 150, 150  # 网格数
Lx, Ly = 30.0, 30.0  # 计算域尺寸
dx = Lx / Nx  # 空间步长

# 相场参数
epsilon = 0.5  # 界面厚度
M_phi = 1.0  # 相场迁移率
tau = 1.0  # 弛豫时间

# 热力学参数
alpha = 1.0  # 热扩散系数
Tm = 1.0  # 熔点
L_latent = 1.0  # 潜热
cp = 1.0  # 比热容

# 过冷度
T_undercool = 0.55  # 初始过冷度

# 创建网格
x = np.linspace(0, Lx, Nx)
y = np.linspace(0, Ly, Ny)
X, Y = np.meshgrid(x, y)

# 初始化
# 底部为固相，顶部为液相
phi = -np.tanh((Y - 5) / (np.sqrt(2) * epsilon))
T = np.ones((Ny, Nx)) * (Tm - T_undercool)  # 初始均匀过冷

# 在界面附近添加噪声触发失稳
noise_region = np.abs(Y - 5) < 3 * epsilon
phi[noise_region] += 0.1 * np.random.randn(np.sum(noise_region))

# 时间参数
dt = 0.005
Nt = 3000
save_interval = 300

# 存储结果
phi_history = [phi.copy()]
T_history = [T.copy()]

print("开始计算二维凝固过程...")

# 时间推进
for n in range(Nt):
    # 计算相场拉普拉斯
    laplacian_phi = np.zeros_like(phi)
    laplacian_phi[1:-1, 1:-1] = (
        phi[2:, 1:-1] + phi[:-2, 1:-1] + phi[1:-1, 2:] + phi[1:-1, :-2] - 4*phi[1:-1, 1:-1]
    ) / dx**2
    
    # 计算温度拉普拉斯
    laplacian_T = np.zeros_like(T)
    laplacian_T[1:-1, 1:-1] = (
        T[2:, 1:-1] + T[:-2, 1:-1] + T[1:-1, 2:] + T[1:-1, :-2] - 4*T[1:-1, 1:-1]
    ) / dx**2
    
    # 相场方程（简化模型）
    m_T = (Tm - T) / T_undercool  # 过冷度驱动项
    dphi_dt = (laplacian_phi + phi * (1 - phi) * (phi - 0.5 + m_T)) / tau
    
    # 温度方程
    dT_dt = alpha * laplacian_T + L_latent / cp * dphi_dt
    
    # 更新
    phi = phi + dt * dphi_dt
    T = T + dt * dT_dt
    
    # 边界条件
    phi[0, :] = 1  # 底部固相
    phi[-1, :] = -1  # 顶部液相
    T[0, :] = Tm  # 底部恒温
    
    # 保存结果
    if n % save_interval == 0:
        phi_history.append(phi.copy())
        T_history.append(T.copy())
        print(f"时间步 {n}/{Nt}, 固相分数: {np.sum(phi > 0) / (Nx*Ny):.3f}")

# 绘制结果
fig, axes = plt.subplots(2, 2, figsize=(14, 12))

# 图1：最终相场分布
ax = axes[0, 0]
im1 = ax.contourf(X, Y, phi_history[-1], levels=20, cmap='RdYlBu_r', vmin=-1, vmax=1)
plt.colorbar(im1, ax=ax, label='相场变量 φ')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('最终相场分布')
ax.set_aspect('equal')

# 图2：最终温度分布
ax = axes[0, 1]
im2 = ax.contourf(X, Y, T_history[-1], levels=20, cmap='hot')
plt.colorbar(im2, ax=ax, label='温度 T')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('最终温度分布')
ax.set_aspect('equal')

# 图3：相场演化过程
ax = axes[1, 0]
for i, phi_snap in enumerate(phi_history[::2]):
    ax.contour(X, Y, phi_snap, levels=[0], colors=[plt.cm.viridis(i/len(phi_history[::2]))], linewidths=2)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('界面演化过程')
ax.set_aspect('equal')

# 图4：固相分数演化
ax = axes[1, 1]
solid_fraction = [np.sum(phi > 0) / (Nx*Ny) for phi in phi_history]
times = np.arange(len(solid_fraction)) * save_interval * dt
ax.plot(times, solid_fraction, 'b-', linewidth=2)
ax.set_xlabel('时间 t')
ax.set_ylabel('固相分数')
ax.set_title('凝固进程')
ax.grid(True)

plt.tight_layout()
plt.savefig('d:\\文档\\500仿真领域\\工程仿真\\对流换热仿真\\主题073_相场法模拟相变\\仿真结果\\case2_2d_solidification.png', dpi=150)
plt.close()

print("案例2完成！")

4.3 案例3：枝晶生长模拟

本案例模拟各向异性条件下的枝晶生长，展示典型的枝晶形貌。

"""
案例3：枝晶生长模拟
模拟各向异性条件下的枝晶生长
"""
import numpy as np
import matplotlib.pyplot as plt

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

# 参数设置
Nx, Ny = 200, 200  # 网格数
L = 40.0  # 计算域尺寸
dx = L / Nx

# 相场参数
epsilon_0 = 0.3  # 基准界面厚度
delta_4 = 0.05  # 四重各向异性强度
tau = 0.3  # 弛豫时间

# 热力学参数
alpha = 2.0  # 热扩散系数
Tm = 1.0  # 熔点
L_latent = 2.0  # 潜热
cp = 1.0  # 比热容

# 过冷度
T_undercool = 0.55

# 创建网格
x = np.linspace(-L/2, L/2, Nx)
y = np.linspace(-L/2, L/2, Ny)
X, Y = np.meshgrid(x, y)

# 初始化：中心为固相晶核
r = np.sqrt(X**2 + Y**2)
phi = -np.tanh((r - 2) / (np.sqrt(2) * epsilon_0))
T = np.ones((Ny, Nx)) * (Tm - T_undercool)

# 时间参数
dt = 0.01
Nt = 4000
save_interval = 400

# 存储结果
phi_history = [phi.copy()]

print("开始计算枝晶生长...")

# 时间推进
for n in range(Nt):
    # 计算梯度
    grad_phi_x = np.zeros_like(phi)
    grad_phi_y = np.zeros_like(phi)
    grad_phi_x[1:-1, :] = (phi[2:, :] - phi[:-2, :]) / (2*dx)
    grad_phi_y[:, 1:-1] = (phi[:, 2:] - phi[:, :-2]) / (2*dx)
    
    # 计算界面法向和角度
    grad_mag = np.sqrt(grad_phi_x**2 + grad_phi_y**2)
    theta = np.arctan2(grad_phi_y, grad_phi_x)
    
    # 各向异性界面厚度
    epsilon = epsilon_0 * (1 + delta_4 * np.cos(4 * theta))
    
    # 计算拉普拉斯
    laplacian_phi = np.zeros_like(phi)
    laplacian_phi[1:-1, 1:-1] = (
        phi[2:, 1:-1] + phi[:-2, 1:-1] + phi[1:-1, 2:] + phi[1:-1, :-2] - 4*phi[1:-1, 1:-1]
    ) / dx**2
    
    laplacian_T = np.zeros_like(T)
    laplacian_T[1:-1, 1:-1] = (
        T[2:, 1:-1] + T[:-2, 1:-1] + T[1:-1, 2:] + T[1:-1, :-2] - 4*T[1:-1, 1:-1]
    ) / dx**2
    
    # 相场方程（含各向异性）
    m_T = (Tm - T) / T_undercool
    noise = 0.05 * np.random.randn(Ny, Nx) * (np.abs(phi) < 0.5)
    dphi_dt = (epsilon**2 * laplacian_phi + phi * (1 - phi) * (phi - 0.5 + m_T)) / tau + noise
    
    # 温度方程
    dT_dt = alpha * laplacian_T + L_latent / cp * dphi_dt
    
    # 更新
    phi = phi + dt * dphi_dt
    T = T + dt * dT_dt
    
    # 边界条件
    phi[0, :] = phi[1, :]
    phi[-1, :] = phi[-2, :]
    phi[:, 0] = phi[:, 1]
    phi[:, -1] = phi[:, -2]
    
    # 保存结果
    if n % save_interval == 0:
        phi_history.append(phi.copy())
        tip_position = np.max(np.sqrt(X[phi > 0.5]**2 + Y[phi > 0.5]**2)) if np.sum(phi > 0.5) > 0 else 0
        print(f"时间步 {n}/{Nt}, 枝晶尖端位置: {tip_position:.2f}")

# 绘制结果
fig, axes = plt.subplots(2, 2, figsize=(14, 12))

# 图1：最终相场分布
ax = axes[0, 0]
im1 = ax.contourf(X, Y, phi_history[-1], levels=20, cmap='RdYlBu_r', vmin=-1, vmax=1)
plt.colorbar(im1, ax=ax, label='相场变量 φ')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('枝晶形貌')
ax.set_aspect('equal')

# 图2：温度场
ax = axes[0, 1]
im2 = ax.contourf(X, Y, T, levels=20, cmap='hot')
plt.colorbar(im2, ax=ax, label='温度 T')
ax.contour(X, Y, phi_history[-1], levels=[0], colors='white', linewidths=1)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('温度场与界面')
ax.set_aspect('equal')

# 图3：界面演化
ax = axes[1, 0]
colors = plt.cm.viridis(np.linspace(0, 1, len(phi_history)))
for i, phi_snap in enumerate(phi_history):
    ax.contour(X, Y, phi_snap, levels=[0], colors=[colors[i]], linewidths=1.5)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('枝晶生长过程')
ax.set_aspect('equal')

# 图4：各向异性函数
ax = axes[1, 1]
theta_plot = np.linspace(0, 2*np.pi, 100)
epsilon_plot = epsilon_0 * (1 + delta_4 * np.cos(4 * theta_plot))
r_plot = epsilon_plot / epsilon_0
ax.plot(theta_plot * 180 / np.pi, r_plot, 'b-', linewidth=2)
ax.axhline(y=1, color='r', linestyle='--', label='各向同性')
ax.set_xlabel('角度 θ (度)')
ax.set_ylabel('ε(θ)/ε₀')
ax.set_title('界面能各向异性')
ax.legend()
ax.grid(True)

plt.tight_layout()
plt.savefig('d:\\文档\\500仿真领域\工程仿真\对流换热仿真\主题073_相场法模拟相变\仿真结果\case3_dendrite_growth.png', dpi=150)
plt.close()

print("案例3完成！")

4.4 案例4：Cahn-Hilliard相分离

本案例使用Cahn-Hilliard方程模拟二元合金的旋节分解过程。

"""
案例4：Cahn-Hilliard相分离
模拟二元合金的旋节分解
"""
import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft2, ifft2, fftfreq

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

# 参数设置
Nx, Ny = 128, 128  # 网格数
L = 2 * np.pi  # 计算域尺寸
dx = L / Nx

# Cahn-Hilliard参数
epsilon = 0.1  # 界面厚度参数
M = 1.0  # 迁移率

# 创建网格
x = np.linspace(0, L, Nx)
y = np.linspace(0, L, Ny)
X, Y = np.meshgrid(x, y)

# 波矢量（用于谱方法）
kx = 2 * np.pi * fftfreq(Nx, dx)
ky = 2 * np.pi * fftfreq(Ny, dx)
KX, KY = np.meshgrid(kx, ky)
k_squared = KX**2 + KY**2
k_fourth = k_squared**2

# 初始条件：均匀混合物添加小扰动
c0 = 0.0  # 平均浓度
phi = c0 + 0.1 * np.random.randn(Ny, Nx)

# 时间参数
dt = 0.001
Nt = 10000
save_interval = 1000

# 存储结果
phi_history = [phi.copy()]

print("开始计算Cahn-Hilliard相分离...")

# 时间推进（谱方法半隐式格式）
for n in range(Nt):
    # 傅里叶变换
    phi_hat = fft2(phi)
    
    # 计算化学势的导数 f'(phi) = phi^3 - phi
    f_prime = phi**3 - phi
    f_prime_hat = fft2(f_prime)
    
    # 半隐式更新（傅里叶空间）
    # phi^{n+1} = (phi^n - dt * M * k^2 * f'^n) / (1 + dt * M * epsilon^2 * k^4)
    phi_hat_new = (phi_hat - dt * M * k_squared * f_prime_hat) / (1 + dt * M * epsilon**2 * k_fourth)
    
    # 反变换回物理空间
    phi = np.real(ifft2(phi_hat_new))
    
    # 保存结果
    if n % save_interval == 0:
        phi_history.append(phi.copy())
        print(f"时间步 {n}/{Nt}, 浓度范围: [{np.min(phi):.3f}, {np.max(phi):.3f}]")

# 绘制结果
fig, axes = plt.subplots(2, 2, figsize=(14, 12))

# 图1：初始状态
ax = axes[0, 0]
im1 = ax.contourf(X, Y, phi_history[0], levels=20, cmap='RdYlBu_r')
plt.colorbar(im1, ax=ax, label='浓度 c')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('初始状态')
ax.set_aspect('equal')

# 图2：中间状态
ax = axes[0, 1]
im2 = ax.contourf(X, Y, phi_history[len(phi_history)//2], levels=20, cmap='RdYlBu_r')
plt.colorbar(im2, ax=ax, label='浓度 c')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('中间状态')
ax.set_aspect('equal')

# 图3：最终状态
ax = axes[1, 0]
im3 = ax.contourf(X, Y, phi_history[-1], levels=20, cmap='RdYlBu_r')
plt.colorbar(im3, ax=ax, label='浓度 c')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title('最终状态')
ax.set_aspect('equal')

# 图4：浓度分布直方图演化
ax = axes[1, 1]
for i, phi_snap in enumerate([phi_history[0], phi_history[len(phi_history)//2], phi_history[-1]]):
    ax.hist(phi_snap.flatten(), bins=30, alpha=0.5, label=f't={i*save_interval*dt*len(phi_history)//2:.1f}')
ax.set_xlabel('浓度 c')
ax.set_ylabel('频率')
ax.set_title('浓度分布演化')
ax.legend()

plt.tight_layout()
plt.savefig('d:\\文档\\500仿真领域\\工程仿真\\对流换热仿真\\主题073_相场法模拟相变\\仿真结果\\case4_cahn_hilliard.png', dpi=150)
plt.close()

print("案例4完成！")

4.5 案例5：熔化过程GIF动画

本案例创建一个熔化过程的GIF动画，展示固相在加热条件下的熔化过程。

"""
案例5：熔化过程GIF动画
模拟固相在加热条件下的熔化过程
"""
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
import io

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

# 参数设置
Nx, Ny = 120, 120
L = 20.0
dx = L / Nx

# 相场参数
epsilon = 0.4
M_phi = 1.0
tau = 0.5

# 热力学参数
alpha = 1.5
Tm = 1.0
L_latent = 1.5
cp = 1.0

# 创建网格
x = np.linspace(0, L, Nx)
y = np.linspace(0, L, Ny)
X, Y = np.meshgrid(x, y)

# 初始化：中心为方形固相
phi = -np.ones((Ny, Nx))
solid_mask = (np.abs(X - L/2) < 4) & (np.abs(Y - L/2) < 4)
phi[solid_mask] = 1
phi = np.tanh(phi / (np.sqrt(2) * epsilon))

# 初始温度（略高于熔点，开始熔化）
T = np.ones((Ny, Nx)) * (Tm + 0.1)

# 边界加热
T[0, :] = Tm + 0.5  # 底部加热
T[-1, :] = Tm + 0.5  # 顶部加热
T[:, 0] = Tm + 0.5  # 左侧加热
T[:, -1] = Tm + 0.5  # 右侧加热

# 时间参数
dt = 0.005
Nt = 2000
save_interval = 100

# 存储帧
frames = []

print("生成熔化过程动画...")

# 时间推进
for n in range(Nt):
    # 计算拉普拉斯
    laplacian_phi = np.zeros_like(phi)
    laplacian_phi[1:-1, 1:-1] = (
        phi[2:, 1:-1] + phi[:-2, 1:-1] + phi[1:-1, 2:] + phi[1:-1, :-2] - 4*phi[1:-1, 1:-1]
    ) / dx**2
    
    laplacian_T = np.zeros_like(T)
    laplacian_T[1:-1, 1:-1] = (
        T[2:, 1:-1] + T[:-2, 1:-1] + T[1:-1, 2:] + T[1:-1, :-2] - 4*T[1:-1, 1:-1]
    ) / dx**2
    
    # 相场方程
    m_T = (Tm - T) / 0.5  # 过热驱动熔化
    dphi_dt = (epsilon**2 * laplacian_phi + phi * (1 - phi) * (phi - 0.5 + m_T)) / tau
    
    # 温度方程
    dT_dt = alpha * laplacian_T + L_latent / cp * dphi_dt
    
    # 更新
    phi = phi + dt * dphi_dt
    T = T + dt * dT_dt
    
    # 保持边界温度
    T[0, :] = Tm + 0.5
    T[-1, :] = Tm + 0.5
    T[:, 0] = Tm + 0.5
    T[:, -1] = Tm + 0.5
    
    # 保存帧
    if n % save_interval == 0:
        fig, axes = plt.subplots(1, 2, figsize=(12, 5))
        
        # 相场
        ax = axes[0]
        im1 = ax.contourf(X, Y, phi, levels=20, cmap='RdYlBu_r', vmin=-1, vmax=1)
        plt.colorbar(im1, ax=ax, label='相场变量 φ')
        ax.contour(X, Y, phi, levels=[0], colors='black', linewidths=1)
        ax.set_xlabel('x')
        ax.set_ylabel('y')
        ax.set_title(f'熔化过程 (t={n*dt:.1f})')
        ax.set_aspect('equal')
        
        # 温度
        ax = axes[1]
        im2 = ax.contourf(X, Y, T, levels=20, cmap='hot')
        plt.colorbar(im2, ax=ax, label='温度 T')
        ax.contour(X, Y, phi, levels=[0], colors='white', linewidths=1)
        ax.set_xlabel('x')
        ax.set_ylabel('y')
        ax.set_title('温度场')
        ax.set_aspect('equal')
        
        plt.tight_layout()
        
        # 转换为图像
        buf = io.BytesIO()
        fig.savefig(buf, format='png', dpi=100, bbox_inches='tight')
        buf.seek(0)
        frames.append(Image.open(buf))
        plt.close()
        
        solid_fraction = np.sum(phi > 0) / (Nx*Ny)
        print(f"帧 {n//save_interval + 1}/{(Nt//save_interval)}, 固相分数: {solid_fraction:.3f}")

# 保存GIF
print("保存GIF动画...")
frames[0].save('d:\\文档\\500仿真领域\\工程仿真\\对流换热仿真\\主题073_相场法模拟相变\\仿真结果\\case5_melting.gif',
               save_all=True, append_images=frames[1:], duration=100, loop=0)

print("案例5完成！")

AtomGit开源社区

AtomGit 是由开放原子开源基金会联合 CSDN 等生态伙伴共同推出的新一代开源与人工智能协作平台。平台坚持“开放、中立、公益”的理念，把代码托管、模型共享、数据集托管、智能体开发体验和算力服务整合在一起，为开发者提供从开发、训练到部署的一站式体验。

更多推荐

为搜索引擎 Agent 设计 Harness 查询改写与扩展

本文将带你深入探讨如何为搜索引擎 Agent 设计一套强大的Harness（框架），专门用于查询改写（Query Rewriting）与查询扩展（Query Expansion）。我们将从理论基础出发，一步步构建一个实用的系统，包括同义词扩展、意图识别、错误纠正、上下文理解等核心功能。架构设计：我们设计了一个从预处理到候选生成再到排序的完整流程查询预处理：实现了文本清洗、分词、去停用词和词形还原查