混合嵌入式间断伽辽金法求解相场晶体方程

1. 混合嵌入式间断伽辽金法在相场晶体方程中的应用解析

相场晶体（Phase Field Crystal, PFC）方程作为连接原子尺度与连续介质尺度的重要桥梁，已成为模拟材料微观结构演化的核心工具。该方程通过引入原子密度场变量，能够自然刻画晶体生长、晶界迁移和位错动力学等复杂现象。然而其六阶非线性偏微分方程的特性，使得传统数值方法面临两大核心挑战：高阶导数离散带来的计算复杂度激增，以及非线性项导致的数值稳定性问题。

针对这些挑战，我们发展了一套基于混合嵌入式间断伽辽金（Hybridizable/Embedded Discontinuous Galerkin, HDG/EDG）的创新数值框架。该方法通过将原始六阶方程重构为四个矢量场和三个标量场组成的一阶方程组，从根本上规避了高阶导数的直接计算。特别设计的连续面函数与静态凝聚技术，使全局耦合自由度数量降低至传统局部间断伽辽金（LDG）方法的12%，在保持高阶精度的同时显著提升计算效率。

关键创新：通过引入辅助变量r=-∇ϕ、ψ=-∇·r、q=-∇ψ等，将原始六阶方程分解为仅含一阶导数的耦合系统。这种重构不仅简化了离散过程，更使得HDG/EDG方法特有的面函数设计成为可能。

2. 方法设计与数学框架

2.1 相场晶体方程的重构

原始PFC方程描述为：

∂tϕ = ∇·(M(ϕ)∇μ) μ = ϕ³ + (1-ε)ϕ + 2Δϕ + Δ²ϕ

通过引入辅助变量，我们将其重构为以下一阶系统：

r = -∇ϕ ψ = -∇·r q = -∇ψ μ = f(ϕ) + 2ψ - ∇·q p = -∇μ s = M(ϕ)p ∂tϕ = -∇·s

这种重构具有三个显著优势：

1. 计算复杂度从O(h⁻⁶)降至O(h⁻²)
2. 各变量物理意义明确（r对应梯度，ψ对应拉普拉斯算子）
3. 自然适应混合有限元离散格式

2.2 HDG/EDG离散策略

在单元K∈𝒯_h上，我们定义离散空间：

• 标量场：W_h = {w∈L²(𝒯_h): w|_K∈P_k(K)}
• 矢量场：V_h = {v∈[L²(𝒯_h)]^d: v|_K∈[P_k(K)]^d}
• 面空间：M_h = {ŵ∈L²(ℰ_h)∩C(ℰ_h): ŵ|_e∈P_k(e)}

离散格式的核心在于数值通量的设计：

r̂·n = r·n + τ₁(ϕ-ϕ̂) q̂·n = q·n + τ₂(ψ-ψ̂) - τ₄(ϕ-ϕ̂) ŝ·n = s·n + τ₃(μ-μ̂)

其中τ₁-τ₄为稳定化参数，通过理论分析我们确定其满足τ₂=τ₁、τ₄=2τ₁时，格式具备无条件稳定性。

3. 无条件能量稳定性证明

3.1 能量函数构造

系统总能量定义为：

E = ∫[1/4ϕ⁴ + (1-ε)/2ϕ² - |r|² + 1/2ψ²]dx

通过能量变分可得连续层面的能量耗散律：

dE/dt + ∫M(ϕ)|∇μ|²dx = 0

3.2 半离散能量稳定性

定理1：当稳定化参数满足τ₂=τ₁、τ₄=2τ₁时，半离散格式满足：

dE_h/dt + ||M(ϕ_h)^{1/2}p_h||² ≤ 0

证明要点：

1. 对重构方程分别取合适的测试函数
2. 利用数值通量的特殊设计消去交叉项
3. 通过格林公式恢复能量导数形式

3.3 全离散凸分裂格式

采用一阶凸分裂时间离散：

(ϕ_h^n - ϕ_h^{n-1})/Δt = ∇·(M(ϕ_h^{n-1})p_h^n) μ_h^n = (ϕ_h^n)³ + (1-ε)ϕ_h^n + 2ψ_h^{n-1} - ∇·q_h^n

定理2：全离散格式满足能量递减：

E_h^n + Δt||M(ϕ_h^{n-1})^{1/2}p_h^n||² ≤ E_h^{n-1}

4. 解的存在唯一性分析

4.1 不动点定理应用

通过构造映射F: W_h→W_h，将离散方程转化为不动点问题。在满足小数据条件：

7-4ε||ϕ_h^{n-1}||² + 11-3ε||ψ_h^{n-1}||² ≤ Ch^d

时，应用Brouwer不动点定理证明解的存在性。

4.2 唯一性证明

利用能量估计和多项式逆不等式，证明当τ_i>0时，两个解的差满足：

||ψ_{1,2}||² + (1-ε)||ϕ_{1,2}||² ≤ 0

从而ϕ_{1,2}=0，解唯一。

5. 数值验证与性能分析

5.1 收敛性测试

在周期边界条件下，采用制造解析解方法验证收敛阶。当选用k次多项式时，观察到：

• L²误差：O(h^{k+1})
• H¹误差：O(h^k)

5.2 自由度对比

EDG方法通过强制面函数连续，使自由度降至接近连续伽辽金方法，同时保持DG方法的局部守恒性。

5.3 复杂微观结构模拟

模拟二维六方晶体生长时，该方法成功捕捉到：

• 晶界缺陷演化
• 枝晶生长各向异性
• 位错攀移动力学

计算中时间步长可超越CFL限制达100倍，验证了无条件稳定的优势。

6. 关键实现细节

6.1 静态凝聚技术

通过以下步骤大幅降低求解规模：

1. 在单元内部消去(r_h,q_h,p_h,s_h)
2. 全局系统仅求解(ϕ̂,ψ̂,μ̂)
3. 后处理恢复内部变量

6.2 稳定化参数选择

推荐取值：

τ₁ = τ₂ = 1/h, τ₃ = h, τ₄ = 2/h

这种缩放关系保证：

• 足够大的τ₁,τ₂维持稳定性
• 较小的τ₃避免过扩散
• τ₄=2τ₁满足能量条件

6.3 非线性迭代策略

采用牛顿迭代处理立方非线性项：

1. 预测步：用ϕ^{n-1}计算M(ϕ)
2. 校正步：更新雅可比矩阵
3. 收敛标准：||δϕ||<10⁻⁶h²

7. 扩展应用与优化方向

本方法可自然推广到：

• 多组分PFC系统
• 电化学耦合模型
• 三维多晶生长模拟

未来优化包括：

• 采用hp自适应网格
• 开发GPU并行算法
• 结合机器学习加速非线性求解

这种混合嵌入式框架为高阶相场模型的高效模拟提供了通用范式，其核心思想也可应用于其他高阶偏微分方程的数值求解。