主题027：相变与凝固模拟

目录

1. 引言
2. 核心理论
3. 实例一：Stefan问题基础
4. 实例二：相变潜热处理
5. 实例三：移动边界追踪
6. 总结与习题

---

引言

1.1 工程背景与应用

相变与凝固现象广泛存在于自然界和工程实践中。从冰川的形成到金属铸件的凝固，从食品冷冻保存到相变储能材料的应用，相变过程对材料性能、产品质量和能源利用效率都有着决定性影响。

典型应用领域：

• 铸造工业：金属凝固过程直接决定铸件的微观组织和力学性能。铸件中的缩孔、偏析等缺陷都与凝固过程密切相关。
• 焊接工艺：熔池的凝固行为影响焊缝的晶粒结构和力学性能。
• 材料制备：单晶硅生长、定向凝固等先进材料制备技术都依赖于对凝固过程的精确控制。
• 能源存储：相变储能材料利用固液相变潜热实现高效热能存储。
• 环境科学：冰川融化、冻土退化等气候变化研究需要准确的相变模型。
• 食品工业：冷冻食品的冻结速率影响产品品质和保质期。

1.2 学习目标

通过本主题的学习，读者将能够：

1. 理解相变问题的数学描述和物理本质
2. 掌握Stefan问题的相似解方法
3. 学会使用焓方法处理相变潜热
4. 了解移动边界追踪的数值方法
5. 能够编写Python代码模拟一维凝固过程

---

核心理论

2.1 相变问题的数学描述

2.1.1 控制方程

相变问题的核心是热传导方程与相变条件的耦合。对于一维凝固问题，固相区的温度场满足：

$$\frac{\partial T}{\partial t}=\alpha \frac{{\partial }^{2}T}{\partial x^2},0<x<s(t)$$

其中：

• $T$ 是温度
• $t$ 是时间
• $x$ 是空间坐标
• $s(t)$ 是固液界面位置（随时间变化）
• $\alpha =k/(\rho c)$ 是热扩散系数

2.1.2 边界条件

在固定边界 $x=0$ 处，通常给定温度或热流边界条件：

$$T(0,t)=T_0<T_m$$

在移动界面 $x=s(t)$ 处，温度为熔点：

$$T(s(t),t)=T_m$$

2.1.3 Stefan条件（能量守恒）

界面处的能量守恒给出了界面移动速度的控制方程：

$$\rho L\frac{ds}{dt}=-k{\frac{\partial T}{\partial x}|}_{x=s(t)}$$

这就是著名的Stefan条件，它表明：

• 界面释放的潜热等于传导到固相的热量
• 温度梯度越大，界面移动越快
• 潜热越大，界面移动越慢

2.2 Stefan数

Stefan数是相变问题中最重要的无量纲参数：

$$St=\frac{c(T_m-T_0)}{L}$$

物理意义：

• 表示显热与潜热的比值
• $St\ll 1$：潜热主导，界面移动缓慢
• $St\gg 1$：显热主导，界面移动迅速

对于水/冰系统，$St\approx 0.1-0.2$，属于潜热主导的情况。

2.3 相似解方法

对于半无限大区域的凝固问题，存在相似解。引入相似变量：

$$\eta =\frac{x}{2\sqrt{\alpha t}}$$

假设温度分布为 $T(x,t)=f(\eta )$，则热方程转化为常微分方程：

$$f^{\prime \prime }+2\eta f^{\prime }=0$$

其解为误差函数形式：

$$T(x,t)=T_0+(T_m-T_0)\frac{\text{erf}(\eta )}{\text{erf}(\lambda )}$$

界面位置满足：

$$s(t)=2\lambda \sqrt{\alpha t}$$

其中 $\lambda$ 由超越方程确定：

$$\lambda e^{{\lambda }^2}\text{erf}(\lambda )=\frac{St}{\sqrt{\pi }}$$

2.4 数值方法概述

2.4.1 焓方法（Enthalpy Method）

焓方法通过引入焓变量 $H$，将移动边界问题转化为固定区域问题：

$$H=\int \rho cdT+\rho L\cdot f_l$$

其中 $f_l$ 是液相分数。焓方法的优点：

• 无需显式追踪界面
• 自然处理糊状区
• 易于扩展到多维问题

2.4.2 坐标变换法

通过坐标变换 $\xi =x/s(t)$，将物理域映射到固定计算域：

$$\frac{\partial T}{\partial t}=\frac{\alpha }{s^2}\frac{{\partial }^{2}T}{\partial {\xi }^2}+\frac{\xi }{s}\frac{ds}{dt}\frac{\partial T}{\partial \xi }$$

第二项是坐标变换引入的对流项。

2.4.3 水平集法

水平集法通过隐式函数 $\varphi (x,t)$ 描述界面：

$$\frac{\partial \varphi }{\partial t}+v|\nabla \varphi |=0$$

界面由 $\varphi (x,t)=0$ 定义。优点：

• 自然处理拓扑变化
• 易于扩展到多维
• 几何属性计算方便

---

实例一：Stefan问题基础

3.1 问题描述

本实例求解一维半无限大区域的凝固问题，使用相似解方法获得解析解，并与数值结果对比。

3.2 完整代码实现

"""
主题027：相变与凝固模拟
实例一：Stefan问题基础

求解一维凝固问题的相似解
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.special import erf
from scipy.optimize import fsolve
import matplotlib
matplotlib.use('Agg')
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False


class StefanProblem:
    """
    Stefan问题求解器 - 精确解（相似解）
    """
    def __init__(self, T_m, T_0, k_s, rho_s, c_s, L):
        """
        参数:
        T_m: 熔点温度 [K]
        T_0: 边界温度 [K] (T_0 < T_m 表示凝固)
        k_s: 固相导热系数 [W/(m·K)]
        rho_s: 固相密度 [kg/m³]
        c_s: 固相比热容 [J/(kg·K)]
        L: 相变潜热 [J/kg]
        """
        self.T_m = T_m
        self.T_0 = T_0
        self.k_s = k_s
        self.rho_s = rho_s
        self.c_s = c_s
        self.L = L
        
        # 固相热扩散系数
        self.alpha_s = k_s / (rho_s * c_s)
        
        # Stefan数 (无量纲参数)
        self.St = c_s * (T_m - T_0) / L

    def solve_lambda(self):
        """求解相似变量 λ (界面位置参数)"""
        from scipy.optimize import fsolve
        
        St = self.St
        
        def equation(lam):
            """Stefan条件方程"""
            if lam <= 0:
                return float('inf')
            return lam * np.exp(lam**2) * erf(lam) - St / np.sqrt(np.pi)
        
        # 初始猜测
        lam_initial = np.sqrt(St / 2)
        
        # 求解
        lam = fsolve(equation, lam_initial)[0]
        return lam

    def temperature(self, x, t, lam):
        """计算温度分布（精确解）"""
        if t <= 0:
            return np.full_like(x, self.T_m)
        
        # 相似变量 η = x / (2*sqrt(α_s*t))
        eta = x / (2 * np.sqrt(self.alpha_s * t))
        
        # 温度分布 (固相区)
        T = self.T_0 + (self.T_m - self.T_0) * erf(eta) / erf(lam)
        
        # 液相区温度为熔点
        s = 2 * lam * np.sqrt(self.alpha_s * t)  # 界面位置
        T = np.where(x <= s, T, self.T_m)
        
        return T

    def interface_position(self, t, lam):
        """计算固液界面位置 s(t) = 2*λ*sqrt(α_s*t)"""
        return 2 * lam * np.sqrt(self.alpha_s * np.maximum(t, 1e-10))

    def interface_velocity(self, t, lam):
        """计算界面移动速度 ds/dt = λ*sqrt(α_s/t)"""
        return lam * np.sqrt(self.alpha_s / np.maximum(t, 1e-10))

    def heat_flux(self, x, t, lam):
        """计算热流密度 q = -k * dT/dx"""
        # 处理标量或数组输入
        t = np.asarray(t)
        x = np.asarray(x)
        
        # 对于t<=0的情况返回0
        t_scalar = np.isscalar(t) or (isinstance(t, np.ndarray) and t.size == 1)
        if t_scalar and t <= 0:
            return np.zeros_like(x, dtype=float)
        
        # dT/dx = (T_m - T_0) / erf(λ) * d[erf(η)]/dη * dη/dx
        # d[erf(η)]/dη = 2/sqrt(π) * exp(-η²)
        # dη/dx = 1 / (2*sqrt(α_s*t))
        
        eta = x / (2 * np.sqrt(self.alpha_s * t))
        
        dT_dx = (self.T_m - self.T_0) / erf(lam) * \
                (2 / np.sqrt(np.pi)) * np.exp(-eta**2) / (2 * np.sqrt(self.alpha_s * t))
        
        q = -self.k_s * dT_dx
        
        # 液相区热流为零
        s = self.interface_position(t, lam)
        if x.size == s.size:
            q = np.where(x <= s, q, 0)
        
        return q

3.3 代码深度解析

3.3.1 Stefan数的计算

self.St = c_s * (T_m - T_0) / L

这行代码计算Stefan数，它是相变问题中最重要的无量纲参数。Stefan数反映了显热与潜热的相对重要性：

• 当 $St\ll 1$ 时，潜热效应占主导
• 当 $St\gg 1$ 时，显热效应占主导

3.3.2 相似变量λ的求解

def equation(lam):
    """Stefan条件方程"""
    if lam <= 0:
        return float('inf')
    return lam * np.exp(lam**2) * erf(lam) - St / np.sqrt(np.pi)

lam = fsolve(equation, lam_initial)[0]

这段代码求解超越方程得到相似变量λ。方程来源于Stefan条件与相似解的结合，需要通过数值方法求解。

3.3.3 温度分布计算

eta = x / (2 * np.sqrt(self.alpha_s * t))
T = self.T_0 + (self.T_m - self.T_0) * erf(eta) / erf(lam)

这里使用误差函数表示温度分布，这是Stefan问题相似解的标准形式。

3.4 运行结果预期

运行代码后将生成以下可视化结果：

1. 温度分布演化图：展示不同时刻固相区的温度分布，呈现典型的误差函数形状
2. 界面位置曲线：界面位置 $s(t)$ 随时间呈 $\sqrt{t}$ 增长
3. 边界热流衰减：热流密度随时间按 $t^{-1/2}$ 衰减
4. Stefan数影响：不同Stefan数下界面移动速度的对比

---

实例二：相变潜热处理

4.1 焓方法原理

焓方法是处理相变问题最常用的数值方法之一。其核心思想是引入焓变量：

$$H=\int \rho cdT+\rho L\cdot f_l$$

其中 $f_l$ 是液相分数。通过焓-温度关系，可以将移动边界问题转化为固定区域问题。

焓-温度本构关系：

$$T=\{\begin{array}{ll}{T}_m+\frac{H-H_s}{\rho c_s}& H<H_s(\text{固相})\\ T_m& H_s\le H\le H_l(\text{糊状区})\\ T_m+\frac{H-H_l}{\rho c_l}& H>H_l(\text{液相})\end{array}$$

4.2 完整代码实现

"""
主题027：相变与凝固模拟
实例二：相变潜热处理 - 焓方法
"""

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.use('Agg')
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False


class EnthalpyMethod:
    """
    焓方法求解器 - 一维相变问题
    """
    def __init__(self, L, T_m, c_s, c_l, rho, k_s, k_l, H_s, H_l):
        """
        参数:
        L: 相变潜热 [J/kg]
        T_m: 熔点温度 [K]
        c_s: 固相比热 [J/(kg·K)]
        c_l: 液相比热 [J/(kg·K)]
        rho: 密度 [kg/m³]
        k_s: 固相导热系数 [W/(m·K)]
        k_l: 液相导热系数 [W/(m·K)]
        H_s: 固相焓值 [J/kg]
        H_l: 液相焓值 [J/kg]
        """
        self.L = L
        self.T_m = T_m
        self.c_s = c_s
        self.c_l = c_l
        self.rho = rho
        self.k_s = k_s
        self.k_l = k_l
        self.H_s = H_s
        self.H_l = H_l
        
        # 计算相变区间
        self.delta_H = H_l - H_s
    
    def temperature_from_enthalpy(self, H):
        """
        从焓计算温度（本构关系）
        """
        T = np.zeros_like(H)
        
        # 固相区
        solid_mask = H < self.H_s
        T[solid_mask] = self.T_m + (H[solid_mask] - self.H_s) / (self.rho * self.c_s)
        
        # 糊状区 (相变区)
        mushy_mask = (H >= self.H_s) & (H <= self.H_l)
        T[mushy_mask] = self.T_m
        
        # 液相区
        liquid_mask = H > self.H_l
        T[liquid_mask] = self.T_m + (H[liquid_mask] - self.H_l) / (self.rho * self.c_l)
        
        return T
    
    def liquid_fraction(self, H):
        """
        计算液相分数
        
        f_l = 0          (H < H_s, 固相)
        f_l = (H-H_s)/ΔH (H_s ≤ H ≤ H_l, 糊状区)
        f_l = 1          (H > H_l, 液相)
        """
        f_l = np.zeros_like(H)
        
        mushy_mask = (H >= self.H_s) & (H <= self.H_l)
        f_l[mushy_mask] = (H[mushy_mask] - self.H_s) / self.delta_H
        
        liquid_mask = H > self.H_l
        f_l[liquid_mask] = 1.0
        
        return f_l
    
    def solve_1d_solidification(self, Nx, L_domain, T_initial, T_boundary, 
                                 dt, t_end):
        """
        求解一维凝固问题
        """
        # 空间离散
        dx = L_domain / (Nx - 1)
        x = np.linspace(0, L_domain, Nx)
        
        # 初始化焓场
        H = np.ones(Nx) * self.rho * self.c_l * (T_initial - self.T_m) + self.H_l
        
        # 边界条件（左端固定温度）
        H[0] = self.rho * self.c_s * (T_boundary - self.T_m) + self.H_s
        
        # 存储结果
        times = [0]
        H_history = [H.copy()]
        T_history = [self.temperature_from_enthalpy(H)]
        f_l_history = [self.liquid_fraction(H)]
        
        # 时间推进
        n_steps = int(t_end / dt)
        
        for n in range(n_steps):
            H_new = self._explicit_step(H, dx, dt, Nx)
            
            # 应用边界条件
            H_new[0] = self.rho * self.c_s * (T_boundary - self.T_m) + self.H_s
            H_new[-1] = H[-1]
            
            H = H_new
            
            # 存储结果
            if (n + 1) % max(1, n_steps // 20) == 0:
                times.append((n + 1) * dt)
                H_history.append(H.copy())
                T_history.append(self.temperature_from_enthalpy(H))
                f_l_history.append(self.liquid_fraction(H))
        
        return {
            'x': x,
            'times': np.array(times),
            'H_history': np.array(H_history),
            'T_history': np.array(T_history),
            'f_l_history': np.array(f_l_history)
        }
    
    def _explicit_step(self, H, dx, dt, Nx):
        """显式时间推进"""
        H_new = H.copy()
        
        for i in range(1, Nx - 1):
            T = self.temperature_from_enthalpy(H)
            k_eff = self.effective_thermal_conductivity(H)
            
            # 热扩散系数
            alpha_eff = k_eff[i] / self.rho
            
            # 显式格式
            H_new[i] = H[i] + dt * alpha_eff * (T[i+1] - 2*T[i] + T[i-1]) / dx**2
        
        return H_new

4.3 糊状区（Mushy Zone）

对于合金材料，凝固过程存在一个温度区间（固相线到液相线之间），在此区间内固液共存，称为糊状区。

糊状区的工程意义：

• 影响溶质再分配和宏观偏析
• 决定铸件的致密性和力学性能
• 与热裂敏感性密切相关

4.4 运行结果预期

运行代码后将生成：

1. 温度分布演化：展示凝固过程中温度场的变化
2. 液相分数分布：清晰显示糊状区的形成和演化
3. 界面位置曲线：固液界面的移动轨迹
4. 糊状区宽度：糊状区随时间的变化规律

---

实例三：移动边界追踪

5.1 坐标变换法

坐标变换法通过引入新的空间坐标，将移动区域映射到固定计算域。

坐标变换：

$$\xi =\frac{x}{s(t)},\xi \in [0,1]$$

变换后的控制方程：

$$\frac{\partial T}{\partial t}=\frac{\alpha }{s^2}\frac{{\partial }^{2}T}{\partial {\xi }^2}+\frac{\xi }{s}\frac{ds}{dt}\frac{\partial T}{\partial \xi }$$

第二项是由于坐标变换引入的"伪对流"项。

5.2 水平集法

水平集法通过隐式函数描述界面位置：

$$\varphi (x,t)=0\text{定义界面}$$

水平集方程：

$$\frac{\partial \varphi }{\partial t}+v|\nabla \varphi |=0$$

其中 $v$ 是界面法向速度。

5.3 完整代码实现

"""
主题027：相变与凝固模拟
实例三：移动边界追踪 - 坐标变换法与水平集法
"""

import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.use('Agg')
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False


class CoordinateTransformationMethod:
    """
    坐标变换法求解器
    """
    def __init__(self, T_m, T_0, alpha, k, rho, L):
        self.T_m = T_m
        self.T_0 = T_0
        self.alpha = alpha
        self.k = k
        self.rho = rho
        self.L = L
        self.St = (T_m - T_0) * 2100 / L
    
    def solve(self, N, t_end, dt):
        """求解移动边界问题"""
        # 计算域 ξ ∈ [0, 1]
        xi = np.linspace(0, 1, N)
        dxi = 1 / (N - 1)
        
        # 初始条件
        s = 1e-6
        T = self.T_m - (self.T_m - self.T_0) * xi
        
        # 存储结果
        times = [0]
        s_history = [s]
        T_history = [T.copy()]
        
        # 时间推进
        n_steps = int(t_end / dt)
        
        for n in range(n_steps):
            # 计算界面速度
            dT_dxi_at_interface = (T[-1] - T[-2]) / dxi
            dT_dx_at_interface = dT_dxi_at_interface / s
            dsdt = -self.k * dT_dx_at_interface / (self.rho * self.L)
            
            # 更新界面位置
            s_new = s + dt * dsdt
            if s_new < 1e-6:
                s_new = 1e-6
            
            # 更新温度场
            T_new = self._update_temperature(T, xi, dxi, s, dsdt, dt)
            
            # 应用边界条件
            T_new[0] = self.T_0
            T_new[-1] = self.T_m
            
            s = s_new
            T = T_new
            
            # 存储结果
            if (n + 1) % max(1, n_steps // 50) == 0:
                times.append((n + 1) * dt)
                s_history.append(s)
                T_history.append(T.copy())
        
        return {
            'xi': xi,
            'times': np.array(times),
            's_history': np.array(s_history),
            'T_history': np.array(T_history)
        }
    
    def _update_temperature(self, T, xi, dxi, s, dsdt, dt):
        """更新温度场"""
        N = len(T)
        T_new = T.copy()
        
        for i in range(1, N - 1):
            # 扩散项
            diffusion = self.alpha / s**2 * (T[i+1] - 2*T[i] + T[i-1]) / dxi**2
            
            # 对流项（坐标变换引入）
            convection = xi[i] * dsdt / s * (T[i+1] - T[i-1]) / (2 * dxi)
            
            T_new[i] = T[i] + dt * (diffusion - convection)
        
        return T_new


class LevelSetMethod:
    """
    水平集法求解器
    """
    def __init__(self, T_m, alpha, k, rho, L, v_normal):
        self.T_m = T_m
        self.alpha = alpha
        self.k = k
        self.rho = rho
        self.L = L
        self.v_normal = v_normal
    
    def initialize_level_set(self, x, interface_pos):
        """初始化水平集函数"""
        return x - interface_pos
    
    def solve(self, N, L_domain, t_end, dt, interface_pos_init):
        """求解水平集方程"""
        dx = L_domain / (N - 1)
        x = np.linspace(0, L_domain, N)
        
        # 初始化
        phi = self.initialize_level_set(x, interface_pos_init)
        T = np.ones(N) * self.T_m
        T[0] = self.T_m - 20
        
        # 存储结果
        times = [0]
        phi_history = [phi.copy()]
        T_history = [T.copy()]
        interface_history = [interface_pos_init]
        
        # 时间推进
        n_steps = int(t_end / dt)
        
        for n in range(n_steps):
            # 计算界面速度
            interface_idx = np.argmin(np.abs(phi))
            if interface_idx > 0 and interface_idx < N - 1:
                dT_dx = (T[interface_idx + 1] - T[interface_idx - 1]) / (2 * dx)
                v = -self.k * dT_dx / (self.rho * self.L)
            else:
                v = 0
            
            # 更新水平集函数
            phi_new = self._update_level_set(phi, dx, dt, v)
            
            # 更新温度场
            T_new = self._update_temperature(T, phi, dx, dt)
            
            # 边界条件
            T_new[0] = self.T_m - 20
            T_new[-1] = self.T_m
            
            phi = phi_new
            T = T_new
            
            # 找到界面位置
            interface_idx = np.argmin(np.abs(phi))
            interface_pos = x[interface_idx]
            
            # 存储结果
            if (n + 1) % max(1, n_steps // 50) == 0:
                times.append((n + 1) * dt)
                phi_history.append(phi.copy())
                T_history.append(T.copy())
                interface_history.append(interface_pos)
        
        return {
            'x': x,
            'times': np.array(times),
            'phi_history': np.array(phi_history),
            'T_history': np.array(T_history),
            'interface_history': np.array(interface_history)
        }
    
    def _update_level_set(self, phi, dx, dt, v):
        """更新水平集函数"""
        N = len(phi)
        phi_new = phi.copy()
        
        for i in range(1, N - 1):
            if v > 0:
                dphi_dx = (phi[i] - phi[i-1]) / dx
            else:
                dphi_dx = (phi[i+1] - phi[i]) / dx
            
            phi_new[i] = phi[i] - dt * v * np.abs(dphi_dx)
        
        return phi_new
    
    def _update_temperature(self, T, phi, dx, dt):
        """更新温度场"""
        N = len(T)
        T_new = T.copy()
        
        for i in range(1, N - 1):
            T_new[i] = T[i] + dt * self.alpha * (T[i+1] - 2*T[i] + T[i-1]) / dx**2
        
        return T_new

5.4 方法对比

特性	坐标变换法	水平集法
计算域	固定（变换后）	固定
界面表示	显式	隐式
拓扑变化	困难	自然处理
多维扩展	复杂	直接
数值稳定性	好	需重新初始化
计算效率	高	中等

5.5 运行结果预期

运行代码后将生成：

1. 坐标变换法结果：计算域和物理域中的温度分布、界面位置、界面速度
2. 水平集法结果：水平集函数演化、温度分布、界面位置
3. 方法对比：两种方法的界面位置对比和特性比较

---

总结与习题

6.1 知识点回顾

本主题涵盖了相变与凝固模拟的核心内容：

1. Stefan问题：移动边界热传导问题的经典模型，具有相似解
2. Stefan数：无量纲参数，表征显热与潜热的相对重要性
3. 焓方法：将移动边界问题转化为固定区域问题的有效方法
4. 糊状区：合金凝固中固液共存的区域，对材料性能有重要影响
5. 坐标变换法：通过坐标变换处理移动边界
6. 水平集法：通过隐式函数描述界面运动

6.2 关键公式总结

Stefan条件：

$$\rho L\frac{ds}{dt}=-k{\frac{\partial T}{\partial x}|}_{x=s(t)}$$

相似解：

$$s(t)=2\lambda \sqrt{\alpha t},\lambda e^{{\lambda }^2}\text{erf}(\lambda )=\frac{St}{\sqrt{\pi }}$$

焓定义：

$$H=\int \rho cdT+\rho L\cdot f_l$$

6.3 常见报错与解决方案

报错信息	原因	解决方案
ValueError: The truth value of an array...	数组与标量比较	使用np.any()或np.all()
RuntimeWarning: overflow encountered	数值溢出	减小时间步长或检查参数单位
ValueError: x and y must have same first dimension	数组维度不匹配	检查输入数组形状
fsolve did not converge	非线性方程求解失败	调整初始猜测值或检查方程

6.4 进阶挑战

1. 多维扩展：将一维焓方法扩展到二维凝固问题
2. 合金凝固：考虑溶质扩散的凝固模型（Scheil方程）
3. 对流效应：在液相区考虑自然对流的影响
4. 各向异性：考虑晶体生长各向异性的相场模型

附录：环境准备

依赖库安装

pip install numpy scipy matplotlib