2. 一维Cahn-Hilliard方程显式有限差分解

一维 Cahn-Hilliard 方程的显式有限差分格式

一维 Cahn-Hilliard 方程写作

∂c/∂t = M ∇²μ， μ = f′(c) − κ ∇²c

其中 c 为保守序参量，M 为迁移率，κ 为梯度能系数，f(c) 取双阱势 f(c)= (c²−1)²/4，于是 f′(c)=c³−c。展开后得到四阶方程：

∂c/∂t = M ∇²(c³ − c) − Mκ ∇⁴c

在一维周期域 [0,L] 上，用均匀网格 x_i = iΔx（i=0,…,N−1），显式欧拉时间推进，二阶中心差分近似空间导数。记 c_i^n 为第 n 步第 i 点的值，则离散格式为：

c_i^{n+1} = c_i^n + dt * M * ( lap_mu_i^n )
lap_mu_i^n = ( mu_{i+1}^n - 2 mu_i^n + mu_{i-1}^n ) / dx^2
mu_i^n = (c_i^n)^3 - c_i^n - kappa * ( c_{i+1}^n - 2 c_i^n + c_{i-1}^n ) / dx^2

周期边界通过索引取模实现：i+1 与 i−1 分别对 N 取模。显式格式的稳定性受四阶项限制，时间步需满足 dt ≲ dx⁴/(16 M κ) 量级，实际计算中应取更保守的值。

线性稳定性与生长率

在 c≈0 附近线性化，令 c = c₀ + δc，δc ∝ exp(ikx + σt)，代入线性化方程得到色散关系：

σ(k) = −M k² ( f″(c₀) + κ k² )

当 f″(c₀)=−1（旋节线内），σ(k) = M k² (1 − κ k²)。可见：

  • k < k_c = 1/√κ 时 σ(k) > 0，扰动增长，发生旋节线分解；
  • k = k_max = 1/√(2κ) 时 σ 最大，对应最快增长波数；
  • k > k_c 时 σ(k) < 0，扰动衰减。

本案例用显式有限差分模拟早期浓度波增长，提取不同波数下的生长率，并与上述理论曲线对比。

一维周期网格与差分模板 i−1 i i+1 … N−2 N−1 0 周期边界：i+1 与 i−1 对 N 取模 二阶中心差分（拉普拉斯） ∇²c ≈ (c_{i+1} − 2c_i + c_{i−1}) / Δx² 四阶项（双拉普拉斯） ∇⁴c ≈ (c_{i+2} − 4c_{i+1} + 6c_i − 4c_{i−1} + c_{i−2}) / Δx⁴ 色散关系 σ(k) = M k²(1 − κk²) k σ k_max = 1/√(2κ) k_c = 1/√κ σ > 0 增长区 σ < 0 衰减区 显式时间推进流程 初始条件 c_i⁰ = c₀ + ε·cos(kx_i) 计算化学势 μ μ = c³ − c − κ∇²c 计算 ∇²μ 中心差分 + 周期边界 更新 c c^{n+1} = c^n + dt·M·∇²μ 重复至目标时间，记录振幅 A(k,t) = 2|c_k(t)| 生长率提取：对早期线性阶段拟合 ln A(k,t) 的斜率 → σ_num(k)，与 σ_theory(k) = M k²(1 − κk²) 对比

显式差分实现要点

下面给出可直接运行的一维 Cahn-Hilliard 显式有限差分核心代码。参数取 M=1.0，κ=1.0，L=100，N=256，dx=L/N，dt=1e-4，总步数 20000。初始条件为 c₀=0 加小振幅单波扰动 c_i = 0.01·cos(2π k x_i / L)，k 取不同整数波数。

import numpy as np

def ch1d_explicit(N=256, L=100.0, M=1.0, kappa=1.0,
                  dt=1e-4, nsteps=20000, k_mode=3, amp=0.01):
    dx = L / N
    x = np.arange(N) * dx
    c = amp * np.cos(2 * np.pi * k_mode * x / L)

    for n in range(nsteps):
        # 周期边界索引
        ip1 = (np.arange(N) + 1) % N
        im1 = (np.arange(N) - 1) % N

        # 化学势 mu = c^3 - c - kappa * lap(c)
        lap_c = (c[ip1] - 2*c + c[im1]) / dx**2
        mu = c**3 - c - kappa * lap_c

        # 拉普拉斯 mu
        lap_mu = (mu[ip1] - 2*mu + mu[im1]) / dx**2

        # 显式更新
        c = c + dt * M * lap_mu

    return x, c

# 提取生长率：记录振幅随时间变化
def growth_rate(k_mode, N=256, L=100.0, M=1.0, kappa=1.0,
                dt=1e-4, nsteps=5000, amp=0.01):
    dx = L / N
    x = np.arange(N) * dx
    c = amp * np.cos(2 * np.pi * k_mode * x / L)
    amps = []
    times = []
    for n in range(nsteps):
        ip1 = (np.arange(N) + 1) % N
        im1 = (np.arange(N) - 1) % N
        lap_c = (c[ip1] - 2*c + c[im1]) / dx**2
        mu = c**3 - c - kappa * lap_c
        lap_mu = (mu[ip1] - 2*mu + mu[im1]) / dx**2
        c = c + dt * M * lap_mu
        if n % 100 == 0:
            # 提取该波数傅里叶分量振幅
            ck = np.fft.fft(c)[k_mode] / N
            amps.append(2 * abs(ck))
            times.append(n * dt)
    return np.array(times), np.array(amps)

生长率提取与理论对比

对每个波数 k_mode，取早期线性阶段（振幅仍远小于 1，例如 A < 0.05）的数据，对 ln A(t) 做线性拟合，斜率即为数值生长率 σ_num。理论生长率由 σ_theory = M k²(1 − κk²) 给出，其中 k = 2π k_mode / L。

下表给出一组示例参数下的预期对比（M=1，κ=1，L=100，N=256，dt=1e-4）：

k_modek = 2πk_mode/Lσ_theoryσ_num（示例拟合值）相对误差
10.06280.003930.003910.5%
20.12570.01560.01550.6%
30.18850.03460.03440.6%
40.25130.05950.05910.7%
50.31420.08770.08700.8%
60.37700.1160.1150.9%
70.43980.1410.1391.4%
80.50270.1570.1541.9%
90.56550.1580.1542.5%
100.62830.1410.1363.5%

注意：上表 σ_num 为示例拟合值，用于说明对比方法，实际运行结果会因 dt、N、拟合区间不同而有差异。当 k 接近 k_c = 1/√κ = 1 时，σ 趋近于 0，数值误差相对增大；当 k > k_c 时 σ < 0，扰动衰减，应观察到振幅单调下降。

显式格式的稳定性条件：dt ≤ dx⁴ / (16 M κ) 量级。本例 dx = 100/256 ≈ 0.39，dx⁴ ≈ 0.023，16Mκ = 16，故 dt ≤ 1.4e-3。取 dt=1e-4 有充足裕量。若增大 κ 或减小 dx，需同步减小 dt，否则出现数值振荡。
提取生长率时，必须确保处于线性增长阶段。若振幅已接近 ±1，非线性项 c³ 主导，生长率不再符合线性理论。建议拟合区间取 A ∈ [0.01, 0.05]，并检查 ln A(t) 是否保持线性。

结果判读与验证

  • 对 k < k_c 的扰动，振幅随时间指数增长，ln A(t) 呈直线，斜率 σ_num > 0；
  • 对 k > k_c 的扰动，振幅衰减，σ_num < 0；
  • σ_num(k) 随 k 先增后减，峰值在 k_max = 1/√(2κ) 附近；
  • 将 σ_num 与 σ_theory 画在同一张图上，若两者在 k < 0.8k_c 范围内相对误差小于 5%，说明显式差分格式和生长率提取流程正确。

通过本案例，可以掌握保守序参量四阶方程的显式离散方法、周期边界处理、线性稳定性理论对比以及生长率提取的完整流程，为后续半隐式谱方法和二维旋节线分解打下基础。

← 上一章 📖 返回目录 下一章 →