Cahn-Hilliard方程迭代低正则性方法的稳定性与收敛性

《Journal of Differential Equations》：Stability and convergence of an iterative low-regularity method for the Cahn-Hilliard equation

　　徐新宇|罗照南|王胜 中国上海复旦大学智能复杂系统研究所

摘要
在偏微分方程研究中，低正则性解是非常重要的研究对象。1949年，诺贝尔奖得主、物理学家和化学家昂萨格提出猜想[38]：正则性的丧失可能导致“异常能量耗散”或“能量级联”，进而引发流体力学中的湍流现象。因此，计算这类低正则性解具有极其重要的意义。本文我们采用低正则性积分器方法（LGI）来研究Cahn-Hilliard方程，该方法无需假设高正则性，同时仍能保持能量稳定性。现有文献中已有部分研究针对不同偏微分方程模型探讨了多种LGI方法，并通过要求非线性源项满足强利普希茨条件或数值解满足一定的L∞范数界限（最大原理）来证明能量耗散现象。然而对于诸如（非局部）Cahn-Hilliard方程这样的模型，最大原理已不再适用。因此，长期以来，如何利用LGI类方法求解此类模型一直是个挑战。本文旨在通过构建稳定的LGI类方法体系，解决那些不存在最大原理的模型问题，具体方法是通过对带稳定项的Cahn-Hilliard方程进行求解并证明其能量耗散特性。我们所提出的低正则性积分器方法无需利普希茨条件或L∞范数限制。此外，我们还通过一种新开发的迭代方法，在仅要求u属于H6阶正则性的极弱条件下，证明了解的L2收敛性，同时还能保持良好的稳定性。

引言
许多描述物理和生物现象的数学模型都可以用偏微分方程来表述，而这些方程在材料科学、生物学以及动力学理论等领域都具有至关重要的作用。在本研究中，我们重点关注两种经典的偏微分方程模型：Allen–Cahn方程和Cahn–Hilliard方程。Allen–Cahn模型最初由[2]提出，用于研究二元合金中晶体晶粒取向之间的竞争现象；而Cahn–Hilliard方程则是在[10]中为模拟相分离与相变过程而发展起来的。它们的精确形式如下：

{?tu=νΔu?f(u), (x,t)∈Ω×(0,∞), u(x,0)=u0},
{?tu=Δ(?νΔu+f(u)), (x,t)∈Ω×(0,∞), u(x,0)=u0}。

其中，实值函数u(x,t)表示两相的混合状态，u=?1和u=+1分别对应各相的纯态。参数ν通常取较小值，我们有时也用ε=ν来表示过渡层宽度。空间域Ω通常取为二维环面T=(R/2πZ)2，不过我们的分析也可以扩展到具有狄利克雷或诺伊曼边界条件的有界区域。非线性势函数通常选取为f(u)=F′(u)=u3?u，F(u)=1/4(u2?1)2，由此可得到能量泛函E(u)=∫Ω(1/2ν|?u|2+F(u))dx。需要说明的是，为简化问题，本文的所有结果均基于此类f(u)形式，但实际上也可推广到其他双井势函数情形。

这两种方程都属于梯度流，会随着时间推移而耗散能量（1.1）。这一性质为全局适定性提供了H1范数上的先验界限，同时也是数值稳定性的重要依据。当ε→0时，Allen–Cahn方程会收敛为平均曲率流，而Cahn–Hilliard方程则会收敛为Mullins–Sekerka流。关于它们的详细渐近行为和严格分析，可参考[30]、[39]、[1]、[26]。正因如此，Allen-Cahn方程和Cahn-Hilliard方程在微分几何领域（如极小曲面等）也发挥着重要作用，相关例子可参阅[22]、[23]、[24]。此外，Cahn–Hilliard方程在流体力学中也具有重要意义，尤其是与纳维-斯托克斯方程结合使用时，能够以热力学上一致的方式模拟多相流动。通过扩散界面法，纳维-斯托克斯-Cahn-Hilliard系统（以及纳维-斯托克斯-Q张量系统等相关模型）无需显式追踪界面，即可捕捉到复杂的界面现象，如相分离、液滴聚合与破碎以及动态接触线等。另一方面，研究低正则性解是流体力学偏微分方程理论中的基础性问题。1949年，诺贝尔奖得主昂萨格提出猜想[38]：正则性的丧失可能导致异常的能量耗散或能量级联，进而引发流体力学中的湍流现象。因此，准确计算这类“粗糙”的解具有非常重要的意义。本研究旨在在合适的低正则性框架下，探讨Cahn-Hilliard方程的极限行为。

目前已有许多数值方法被用于模拟Cahn–Hilliard方程及相关相场模型。这些方法可根据其对时间离散化和空间近似方式的不同进行分类。在时间步进策略方面，既有完全显式的方法（如前向欧拉法[4]），也有半隐式（隐式-显式）方案[15]、[29]、[36]、[34]、[35]，它们都已被广泛研究。为了在保持结构特性的同时实现更大的时间步长，人们还引入了凸分裂法[25]以及稳定的半隐式方案[29]、[42]、[48]、[49]、[35]、[36]、[34]、[16]、[17]。在空间离散化方面，常见的方法包括谱方法[8]、[14]、[35]、[46]、[15]、[29]、[37]，有限元方法[9]、[27]、[31]，以及有限差分方法[6]、[7]。此外，这里还列举了一些其他相关梯度流模型的研究结果[11]、[20]、[28]、[21]、[40]、[44]、[45]等。这些数值近似方法能够成功捕捉到解的定量和定性特征，其中尤为重要的是能够保持能量耗散的特性。

在这类梯度流的数值分析中，能量稳定性是一个核心问题。例如，冯和普罗尔[27]为Cahn–Hilliard方程的有限元离散化给出了误差估计，但这些估计结果严重依赖于界面参数ν。而完全显式方法通常存在较大的时间步长限制，且往往无法保证能量衰减；相比之下，半隐式方法——即线性部分以隐式方式处理，非线性部分以显式方式处理——则提供了一种更为有效的替代方案。不过，标准的半隐式方案在较大时间步长下仍可能变得不稳定。为了解决这一问题，人们开发了带有辅助稳定项的稳定化版本，从而在更宽松的条件下确保能量衰减，参见[29]、[48]、[49]、[43]。然而，这些方法中的很多仍然依赖于强利普希茨条件或对数值解的先验L∞范数限制。最近的进展[34]、[35]、[16]则借助在[5]、[33]中发展出的更精妙的分析工具，放宽了这些要求。

受半隐式方案进展的启发，所谓低正则性积分器（LGI）作为一种有效方法应运而生，尤其适用于那些正则性有限的方程[32]、[3]、[12]、[18]、[41]、[47]。这类方法旨在即使在低正则性条件下，也能准确捕捉能量耗散及其他结构特性。简而言之，LGI方法是对非线性偏微分方程解的半群表示形式进行离散化。以Allen-Cahn方程为例，众所周知，该方程的平滑解可通过半群的温和形式（即杜阿梅尔公式）给出：

u(t)=e^νtΔu(t?)+∫_t?^t e^ν(t?σ)Δf(u(σ))dσ，

其中etΔ表示标准的热半群（热核）。经典积分器方法的思路在于，在每个时间区间(tn,tn+1)内，通过将积分函数近似为常数来对半群表示形式进行离散化（参见[32]）：

u(tn+1)=e^ντΔu(tn)+∫_t_n^{tn+1} e^ν(tn+1?σ)Δf(u(σ))dσ≈e^ντΔu(tn)+∫_t_n^{tn+1} e^ν(tn+1?σ)Δf(u(tn))dσ=e^ντΔu(tn)+1?e^ν(tn+1?tn)Δ?νΔf(u(tn))，

或者等价地，将u(tn)≈un，u(tn+1)≈un+1，且τ=tn+1?tn，这样就可以得到Allen-Cahn方程的LGI方法形式：

un+1=e^ντΔun+e^ντΔ?1/νΔf(un)。

LGI方法(1.4)的主要优势之一在于其对精确解的正则性要求更低。为此，我们首先回顾一下常用的半隐式时间步进方法：

un+1?un/τ=νΔun?f(un)。

需要指出的是，半隐式时间步进方法(1.5)的截断误差属于O(τ??ttu)阶；因此，为保证收敛性，最低的正则性要求是?ttu≈Δ2u∈L2，即u∈H?。接下来我们简要分析LGI方法(1.4)的截断误差。从(1.3)式可以明显看出：

u(tn+1)=e^ντΔu(tn)+e^ν(tn+1?tn)Δ?1/νΔf(u(tn))+Rn，

其中Rn=∫_t_n^{tn+1} e^ν(tn+1?σ)Δ(f(u(σ))?f(u(tn)))dσ。根据Allen-Cahn方程的最大原理（参见[16]），我们可以估计Rn的范数：

‖Rn‖2?τmax_σ∈[tn,tn+1]‖u(σ)?u(tn)‖2，

这里的‖?‖2指的是第2节中将定义的L2范数。根据连续性方程，我们有：

u(σ)?u(tn)=ν∫_tn^σ Δu(σ?)dσ??∫_tn^σ f(u(σ?))dσ?。

因此，误差的范数可表示为：

‖Rn‖2?τ2max_σ∈[0,T]‖Δu(σ)‖2。

显然，此时的正则性要求为Δu∈L2或u∈H2。与常规的半隐式时间步进方法相比，这种LGI方法所需的正则性要低得多。不难看出，在Cahn-Hilliard方程中，半隐式方法要求u∈H?，而LGI方法仅需u∈H?即可。

不过，上述思路要么需要对非线性源项施加强利普希茨条件，要么需要对数值解设定一定的L∞范数限制。Allen-Cahn方程由于存在标准的最大原理，满足这些要求（相关示例可参阅[34]、[16]）。但对于其他模型，比如（非局部）Cahn-Hilliard方程，最大原理已不再成立。因此，如何在不再适用强利普希茨条件或最大原理的情况下，利用LGI方法求解偏微分方程，长期以来一直是一个未解之题。本文的一个目标，就是为这类模型应用LGI类方法提供系统的解决方案。另一个主要难点在于，在尽可能低的正则性假设下证明解的L2收敛性。因为如果直接应用常规的离散格伦沃尔不等式，可能需要更高的正则性假设才能得到所需的O(τ)阶误差。因此，我们的目标是证明在最低正则性假设下依然能够实现L2收敛，同时还能保持良好的稳定性。

更具体地说，我们对Cahn-Hilliard方程采用以下稳定的LGI方法：

{un+1=e^?νΔ2τun+1?e^?νΔ2τνΔΠ_N[f(un)]?SτΔ2(un+1?un), u0=Π_Nu0},

其中τ为时间步长，S>0为O(τ)阶正则化项的系数，Π_N则表示将L2函数的所有傅里叶模式截断到|k|∞≤N的范围，此处k=(k1,k2)∈Z2，|k|∞=max{|k1|,|k2|}。实际上，上述公式(1.10)可以重新整理为如下形式：

{un+1?un/τ+SΔ2(un+1?un)=e^?νΔ2τ?1/τun+1?e^?νΔ2ττνΔΠ_N[f(un)], u0=Π_Nu0}。

**注1.1** 值得一提的是，从公式(1.11)可以看出，我们的积分器方法与显式前向欧拉方法的特性非常相似，因为：

e^?νΔ2τ?1/τ=?νΔ2+1/2ν2Δ?τ+O(τ2)，
1?e^?νΔ2ττνΔ=Δ?1/2νΔ3τ+O(τ2)。

因此，非正式地讲，公式(1.11)可以被理解为在原有显式前向欧拉方法的基础上加入了稳定项的形式：

un+1?un/τ+SΔ2(un+1?un)=?νΔ2un+ΔΠ_N[f(un)]+O(ντ,τ2)。

**注1.2** 稳定项SΔ2(un+1?un)并非唯一。实际上，许多其他的稳定项也同样有效，比如Se^?Δ(un+1?un)等等。另一种特殊的选择是ν2Δ2(un+1?un)?SΔ(un+1?un)。这种选择有助于分析能量耗散现象。更准确地说，我们需要ν2Δ2(un+1?un)来抵消“显式前向”扩散所带来的负面影响；实际上真正的稳定项应该是?SΔ(un+1?un)。我们选择SΔ2(un+1?un)只是为了简化表达。不过由于这种选择，所需的正则性会更高，更多细节可见注1.11。上述讨论的稳定项都属于O(τ)阶，当然也可以引入O(τ2)阶的稳定项，比如τSΔ2(un+1?un)；不过这样一来，能量稳定性会降低，因为时间步长τ需要更小（小于某个与ν相关的量）。尽管这类稳定项并非唯一，但它们却是不可或缺的。正如第5节将通过数值实验所证明的，如果不使用稳定项，解的能量将会不断增长。

**注1.3** 理论上的显式前向欧拉方法似乎很难实现能量稳定性，因为它很难得到统一的索伯列夫Hs阶范数估计。据我们所知，这个问题源于显式前向欧拉方法的本质——由于“显式前向”扩散项?Δ2un的存在，步骤n处的正则性要高于步骤n+1处的正则性。这似乎正是我们的积分器方法(1.11)与没有稳定项SΔ2(un+1?un)以及更高阶O(ντ,τ2)修正项的纯粹显式前向欧拉方法(1.12)之间的关键技术差异。

**注1.4** 如前文注1.1所述，我们的方法(1.11)可以被看作是显式前向欧拉方法的一种形式。众所周知，当在尺寸为h的网格上进行空间离散化后，任何显式方法的时间步长τ都受到限制，对于四阶方程来说，τ需满足τ0，且假设u0∈H2(T2)。那么存在一个仅与初始能量E0=E(u0)相关的常数β0，使得只要满足：

S≥β?(‖u0‖H22+ν?1|log ν|+ν?1log N+ν)，

就能保证β≥β0。那么对于所有的n≥0，以及任意选定的时间步长τ，都有E(un+1)≤E(un)，其中能量E的定义见(1.1)。注1.7：此处所选的S并非最优值；实际上，我们建议读者参考第5节，那里有使用更小S值的数值实验结果。确定稳定器的最优界限可能是一个非常有趣的问题，但这并非本文的重点，我们将这个问题留给读者去探讨。定理1.8：固定ν>0，设u0∈Hs，且s≥6。设0<τ≤d，其中d>0。设u(t)为具有初始条件u0的2D Cahn-Hilliard方程(CH)的连续解。设um（1≤m≤M，M为固定的迭代次数）根据(1.11)式定义，初始条件也为u0。假设S满足定理1.6中的相同条件。定义t0=0，且对于m≥1，有tm=mτ。那么对于所有满足1≤m≤M的m，都有‖um?u(tm)‖2≤(1+S)?C2?e^C1tm(N?s+τ+τs^8+τ?N?s+4)，其中C1>0是依赖于ν和u0的常数；C2>0是依赖于s、ν和u0的常数。注1.9：值得一提的是，定理1.6中要求S≥βν?1log?N这一条件在数值方案的收敛性方面会带来麻烦。然而如定理1.8所述，随着N→∞，S?N?s→0，因此这样的S选择不会破坏收敛性。另一方面，为了使数值方案能够收敛，我们确实需要(τ+τs^8)?S→0当τ→0时，这就要求τ<(log?N)?1。正如注1.4中所讨论的，这仍然是一个相当宽松的限制，因为如果采用空间离散化方法，为保证数值稳定性，对时间步长τ的大小限制通常要小得多（为1N的多项式级）。注1.10：这里我们需要确保τ不能任意大。但实际上这不是什么大问题，因为我们总是使用较小的（或至少不是任意大的）时间步长。注1.11：此处选择s≥6主要是为了获得高阶稳定性的效果；实际上，如果允许较低的稳定性（即更小的时间步长τ），那么通过截断误差分析，我们甚至可以在比所给条件更弱的假设下得到L2误差估计。另一方面，收敛速率τs^8是一个技术性指标。实际上，如果我们像通常的半隐式时间步进那样允许s≥8，那么就可以实现所需的一阶收敛性。另外还需要指出的是，已有诸如[32]之类的研究为处理初始数据规律性较差的系统，例如不可压缩Navier-Stokes系统，开发出了相应的策略。事实上，由于Navier-Stokes方程的不可压缩性，能量耗散是自然而然发生的，不再构成挑战。注1.12：证明定理1.8的主要困难在于最低要求的H6阶规律性假设。实际上，如果直接应用常规的离散Gr?nwall不等式，我们可能需要H10阶的假设才能得到所需的O(τ)阶误差。出现这种情况主要是因为无法控制L2Gn项。为了解决这一难题，我们开发了一种新的迭代离散Gr?nwall框架，成功地将要求降低到了H6阶。更多相关讨论请参阅第5节。本文的结构如下：第2节列出符号说明及预备知识，包括一些有用的引理。第3节展示Cahn-Hilliard方程的长期渐近行为及其平滑效应。第4节将阐述2D Cahn-Hilliard方程的LGI方法的能量稳定性，而误差估计则在第5节给出。附录A将呈现数值实验结果。第6节为进一步讨论和总结性意见。