4. 二维旋节线分解的随机初值与结构因子演化

问题设定与数值方案

本节在二维周期域 Ω=[0,L]×[0,L] 上求解 Cahn-Hilliard 方程：

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

取 L=64、N=128（Δx=0.5）、M=1、κ=1，则界面宽度 ε=√(2κ)=√2。初值 c₀(x,y)=c̄+η(x,y)，其中 c̄=0（临界淬火），η 为零均值高斯白噪声，幅度 A=0.05。采用半隐式傅里叶谱方法：非线性项 f′(c)=c³−c 在实空间计算后变换到谱空间，线性项 −κk⁴ 隐式处理，时间步 Δt=0.05。

更新格式（谱空间）：ĉ^{n+1} = [ĉ^n − Δt·M·k²·N̂(f′(c^n))] / [1 + Δt·M·κ·k⁴]

其中 N̂ 表示非线性项的傅里叶变换。该格式对刚性线性项无条件稳定，允许较大时间步。

结构因子与特征波数

结构因子定义为 S(k,t)=|ĉ(k,t)|²/N²，对同一 |k| 的模做环向平均得到 S(k,t)。旋节线分解早期，S(k,t) 在特征波数 k_m 处出现峰值，k_m 由线性稳定性理论给出：k_m=√(−f″(c̄)/(2κ))=√(1/2)≈0.707（对应波长 λ_m=2π/k_m≈8.89）。

随时间演化，峰值位置 k_m 基本不变，但峰值强度指数增长后饱和；后期结构因子呈现标度行为：S(k,t)~t^α·F(k/k_m(t))，其中 k_m(t) 随时间向低波数移动（粗化），α 为标度指数。

可运行实现（Python / NumPy）

以下代码在二维周期域上以高斯随机场为初值，用半隐式傅里叶谱方法推进 Cahn-Hilliard 方程，并计算环向平均结构因子。

import numpy as np

# ---------- 参数 ----------
N = 128          # 网格数
L = 64.0         # 域尺寸
dx = L / N       # 网格间距
dt = 0.05        # 时间步
M = 1.0          # 迁移率
kappa = 1.0      # 梯度系数
c_bar = 0.0      # 平均浓度（临界淬火）
A = 0.05         # 噪声幅度
n_steps = 2000   # 总步数
save_every = 200 # 每隔多少步保存结构因子

# ---------- 波数 ----------
kx = np.fft.fftfreq(N, d=dx) * 2 * np.pi
ky = np.fft.fftfreq(N, d=dx) * 2 * np.pi
KX, KY = np.meshgrid(kx, ky, indexing='ij')
K2 = KX**2 + KY**2
K4 = K2**2

# ---------- 初值：高斯随机场 ----------
rng = np.random.default_rng(42)
c = c_bar + A * rng.standard_normal((N, N))

# ---------- 半隐式谱方法 ----------
# 线性项隐式因子
denom = 1.0 + dt * M * kappa * K4

def compute_Sk(c_hat):
    """计算环向平均结构因子 S(k)"""
    S2d = np.abs(c_hat)**2 / N**2
    k_mag = np.sqrt(K2)
    k_bins = np.arange(0, k_mag.max() + 0.5, 0.5)
    S_k = np.zeros(len(k_bins) - 1)
    for i in range(len(k_bins) - 1):
        mask = (k_mag >= k_bins[i]) & (k_mag < k_bins[i+1])
        if mask.sum() > 0:
            S_k[i] = S2d[mask].mean()
    k_centers = 0.5 * (k_bins[:-1] + k_bins[1:])
    return k_centers, S_k

# 记录
S_history = []
time_history = []

for step in range(n_steps + 1):
    if step % save_every == 0:
        c_hat = np.fft.fft2(c)
        kc, Sk = compute_Sk(c_hat)
        S_history.append(Sk)
        time_history.append(step * dt)

    if step == n_steps:
        break

    # 实空间计算非线性项 f'(c) = c^3 - c
    dfdc = c**3 - c
    dfdc_hat = np.fft.fft2(dfdc)

    # 谱空间更新
    c_hat = np.fft.fft2(c)
    c_hat_new = (c_hat - dt * M * K2 * dfdc_hat) / denom
    c = np.real(np.fft.ifft2(c_hat_new))

# ---------- 输出特征波数 ----------
k_centers, Sk_final = compute_Sk(np.fft.fft2(c))
k_peak = k_centers[np.argmax(Sk_final)]
print(f"最终特征波数 k_m ≈ {k_peak:.3f}")
print(f"理论预测 k_m = {np.sqrt(0.5):.3f}")

结果判读与验证

运行上述代码后，可观察以下特征：

• 早期峰值位置：结构因子在 k≈0.7 附近出现峰值，与线性稳定性理论预测的 k_m=√(1/2)≈0.707 一致。
• 峰值增长：S(k_m,t) 在早期呈指数增长，增长速率 R(k_m)=−Mk_m²(f″(c̄)+κk_m²)。代入 f″(0)=−1、κ=1、M=1，得 R(k_m)=−k_m²(−1+k_m²)=k_m²(1−k_m²)。在 k_m=0.707 处 R≈0.25。
• 后期粗化：随着时间推进，峰值向低波数移动，k_m(t) 逐渐减小，结构因子曲线在低 k 区域抬升，呈现标度行为。
• 标度检验：将不同时刻的 S(k,t) 按 k/k_m(t) 和 S(k,t)/S(k_m,t) 归一化后，曲线应趋于重合，验证动态标度假设。