粘性Cahn-Hilliard方程的二阶稳定Crank-Nicolson-Leapfrog格式 A Stable Second-Order Crank-Nicolson-Leapfrog Scheme for the Viscous Cahn-Hilliard Equation

doi:10.12677/AAM.2023.121037

Advances in Applied Mathematics
Vol. 12 No. 01 ( 2023 ), Article ID: 60967 , 13 pages
10.12677/AAM.2023.121037

粘性Cahn-Hilliard方程的二阶稳定 Crank-Nicolson-Leapfrog格式

刘静，王旦霞^*，王李靖

●How to Cite this Article

太原理工大学数学系，山西晋中

收稿日期：2022年12月28日；录用日期：2023年1月24日；发布日期：2023年1月31日

摘要

本文研究了粘性Cahn-Hilliard方程的一个二阶逼近。对于双阱势函数，本文通过引入拉格朗日乘子得到一个等价形式。其次，使用Crank-Nicolson-Leapfrog格式进行时间离散，使用有限元方法进行空间离散，从而得出了一个二阶线性无条件稳定的数值格式。然后，证明了数值格式的无条件能量稳定性和误差分析。最后，给出了几个数值模拟，验证了格式的数值精度。

关键词

粘性Cahn-Hilliard，Crank-Nicolson-Leapfrog，无条件稳定，误差估计

A Stable Second-Order Crank-Nicolson-Leapfrog Scheme for the Viscous Cahn-Hilliard Equation

Jing Liu, Danxia Wang^*, Lijing Wang

Department of Mathematics, Taiyuan University of Technology, Jinzhong Shanxi

Received: Dec. 28^th, 2022; accepted: Jan. 24^th, 2023; published: Jan. 31^st, 2023

ABSTRACT

In this paper, we present a second-order approximation of the viscous Cahn-Hilliard equation. Firstly, an equivalent form of the system has been obtained by introducing a Lagrange multiplier for the double-well potential function. Secondly, a second-order linear unconditionally stable numerical scheme is proposed by using the Crank-Nicolson-Leapfrog scheme and mixed finite element method, respectively, to discrete the time and space. Furthermore, we prove that the scheme is second-order convergent. Finally, numerical examples are performed to show that the numerical accuracy of the proposed scheme is accurate and effective.

Keywords:Viscous Cahn-Hilliard, Crank-Nicolson-Leapfrog, Unconditionally Stable, Error Estimates

This work is licensed under the Creative Commons Attribution International License (CC BY 4.0).

http://creativecommons.org/licenses/by/4.0/

1. 引言

Cahn-Hilliard方程是1958年由Cahn和Hilliard [1] 提出的，是用来描述热力学中两种物质(如合金，聚合物等等)之间相互扩散现象的数学模型。目前，Cahn-Hilliard方程的数值算法及模拟已成为科学计算领域中的一个国际热点问题，国内外诸多学者已经提出了许多Cahn-Hilliard方程的数值求解方法，其研究成果颇多 [2] [3] [4] [5]。

本文研究粘性Cahn-Hilliard方程

${\begin{cases} u_{t} = Δ ω, \\ ω = - ε^{2} Δ u + F^{'} (u) + β u_{t}, \\ \partial_{n} u = \partial_{n} ω = 0, \\ u (x, 0) = u_{0} (x), \end{cases}$ (1)

其中 $F (u) = \frac{1}{4} {(u^{2} - 1)}^{2}$ ，u是混合物中某种物质的浓度， $ε$ 是描述界面厚度的参数， $β \geq 0$ 表示粘性参数。其能量泛函定义为：

$E (u) = \int_{Ω} (\frac{ε^{2}}{2} {| \nabla u |}^{2} + F (u)) d x .$ (2)

当 $β = 0$ 时，方程(1.1)就变为经典的Cahn-Hilliard方程。这类四阶非线性扩散方程在理论和数值方面都已经被广泛地研究。当 $β > 0$ 时，为粘性Cahn-Hilliard方程，它源于动力学，是玻璃和聚合体系统中相位差现象的闭联集模型。Novick-Cohen在文献 [6] 中提出了该模型，给出了推导过程和在物理上的诸多应用。目前，在数值研究方面，有限差分法，谱方法和有限元方法有了不少研究工作。在文献 [7] 中，作者提出了带有浓度迁移率的粘性Cahn-Hilliard方程数值格式。在文献 [8] 中，作者讨论了粘性Cahn-Hilliard方程的一个二阶数值格式。文献 [9] 给出了一个线性化差分格式并验证了格式满足唯一可解性。文献 [10] 提出了一个二阶对数势格式。

研究此类方程的主要困难是Cahn-Hilliard方程的非线性项，针对此困难，已经有了很多研究。比如，在 [11] 中使用了凸分裂方法，在 [12] 中使用了SAV方法，在 [13] 中使用了IEQ方法。本文基于拉格朗日乘子的方法，提出了一个新的二阶格式，并对此格式进行了严格的分析。

本文的其余部分安排如下。在第2节，通过引入拉格朗日乘子将势函数线性化，得到其等价形式，并分别给出了方程的半离散格式和全离散格式；在第3节和第4节中，证明了所提格式满足无条件能量稳定性，并证明了误差估计。第5节，通过几个数值模拟，验证了上述数值格式的准确性。第6节给出了全文的总结。

2. 模型及数值格式

2.1. 预备知识

给定有界开基 $Ω \in R^{d} (d = 2, 3)$ ，其边界 $\partial Ω$ 光滑，C是一个不依赖于h和 $τ$ 的一般正常数，设 $L^{2} (Ω)$ 表示在 $Ω$ 上Lebesgue平方可测的函数空间。内积和范数分别为：

$(u, v) = \int_{Ω} u v d Ω, ‖ u ‖ = \sqrt{(u, u)}, \forall u, v \in L^{2} (Ω) .$

对于 $k \in N$ ， $H^{k} (Ω)$ 是通常的Sobolev空间。 $H^{k} (Ω)$ 的范数定义如下：

${‖ u ‖}_{H^{k}} = {(\sum_{0 \leq α \leq k} {‖ D^{α} u ‖}^{2})}^{\frac{1}{2}}, \forall u \in H^{k} (Ω) .$

2.2. 半离散数值格式

本小节对于双阱势函数 $F^{'} (u) = u (u^{2} - 1)$ ，我们引入一个辅助变量q。令 $q = u^{2} - 1$ ，使得 $F^{'} (u) = u q$ 。进一步我们将q对t求导得 $q_{t} = 2 u u_{t}$ ，则模型可重写为如下形式

$\begin{array}{l} u_{t} = Δ ω, \\ ω = u q - ε^{2} Δ u + β u_{t}, \\ \frac{1}{2} q_{t} = u u_{t} . \end{array}$ (3)

相应的能量泛函 $E (u)$ 为

$E (u) = \int_{Ω} (\frac{ε^{2}}{2} {| \nabla u |}^{2} + \frac{1}{4} {| q |}^{2}) d x .$ (4)

可以发现(4)式仍然保持能量耗散。

构造系统(3)的弱解形式

$(u_{t}, v) + (\nabla ω, \nabla v) = 0,$ (5)

$(ω, ψ) - ε^{2} (\nabla u, \nabla ψ) - (u q, ψ) - β (u_{t}, ψ) = 0,$ (6)

$(\frac{1}{2} q_{t}, p) - (u u_{t}, p) = 0.$ (7)

设 $N \in Z^{+}$ ， $0 = t_{0} < t_{1} < \dots < t_{N} = T$ 为 $[0, T]$ 上的剖分， $[0, T]$ 为时间区间，任给函数 $u (t)$ ，定义 $u^{n}$ 为 $u (n τ)$ 的近似，其中 $τ : = t_{n + 1} - t_{n}, n = 0, 1, \dots, N$ 。构造系统(5)~(7)在时间上离散的Crank-Nicolson-Leapfrog格式，给定 $(u^{n - 1}, u^{n})$ ，当 $n \geq 1$ 时，求 $(u^{n + 1}, q^{n + 1})$

$(\frac{u^{n}^{+ 1} - u^{n}^{- 1}}{2 τ}, v) + (\nabla ω^{n}, \nabla v) = 0,$ (8)

$(ω^{n}, ψ) - ε^{2} (\nabla \frac{u^{n}^{+ 1} + u^{n}^{- 1}}{2}, \nabla ψ) - (u^{n} \frac{q^{n}^{+ 1} + q^{n}^{- 1}}{2}, ψ) - β (\frac{u^{n}^{+ 1} - u^{n}^{- 1}}{2 τ}, ψ) = 0,$ (9)

$\frac{1}{2} (\frac{q^{n}^{+ 1} - q^{n}^{- 1}}{2 τ}, p) - (u^{n} \frac{u^{n}^{+ 1} - u^{n}^{- 1}}{2 τ}, p) = 0.$ (10)

2.3. 全离散数值格式

在区域 $Ω$ 上作拟一致剖分，记作 $T_{h} = K$ ， $h_{i}$ 为网格大小， $h = \max_{0 \leq i \leq N} h_{i}$ ， $S_{h}$ 是分片连续的有限元空间，定义如下

$S_{h} = {v_{h} \in C (Ω) | {v_{h} |}_{K} \in P_{k} (x, y), K \in T_{h}} \subset H^{1} (Ω),$

其中 $P_{k} (x, y)$ 是 $x, y$ 的次数不超过 $k \in Z^{+}$ 的线性多项式的集合。定义 $L_{0}^{2} = {u \in L^{2} (Ω) | (u, 1) = 0}$ ， ${\dot{S}}_{h} = S_{h} \cap L_{0}^{2} (Ω)$ 。接下来构造模型(5)~(7)的全离散格式，给定 $(u_{h}^{n - 1}, u_{h}^{n})$ ，当 $n \geq 1$ 时，求 $(u_{h}^{n + 1}, q_{h}^{n + 1})$

$(\frac{u_{h}^{n + 1} - u_{h}^{n - 1}}{2 τ}, v) + (\nabla ω_{h}^{n}, \nabla v) = 0,$ (11)

$(ω_{h}^{n}, ψ) - ε^{2} (\nabla \frac{u_{h}^{n + 1} + u_{h}^{n - 1}}{2}, \nabla ψ) - (u_{h}^{n} \frac{q_{h}^{n + 1} + q_{h}^{n - 1}}{2}, ψ) - β (\frac{u_{h}^{n + 1} - u_{h}^{n - 1}}{2 τ}, ψ) = 0,$ (12)

$\frac{1}{2} (\frac{q_{h}^{n + 1} - q_{h}^{n - 1}}{2 τ}, p) - (u_{h}^{n} \frac{u_{h}^{n + 1} - u_{h}^{n - 1}}{2 τ}, p) = 0$ .(13)

3. 稳定性分析

定理3.1令 $(u_{h}^{n + 1}, q_{h}^{n + 1})$ 是(11)~(13)的解，定义

$E^{n + 1, n} = \frac{ε^{2}}{2} ({‖ \nabla u_{h}^{n + 1} ‖}^{2} + {‖ \nabla u_{h}^{n} ‖}^{2}) + \frac{1}{4} ({‖ q_{h}^{n + 1} ‖}^{2} + {‖ q_{h}^{n} ‖}^{2}) .$

对任意的 $τ, h, ε > 0, n \geq 1$ ，有下述不等式成立

$E^{n + 1, n} + \frac{β}{2 τ} {‖ u_{h}^{n + 1} - u_{h}^{n - 1} ‖}^{2} + 2 τ {‖ \nabla ω_{h}^{n} ‖}^{2} = E^{n, n - 1} .$ (14)

证明：令 $v = 2 τ ω_{h}^{n}$ ，得

$(u_{h}^{n + 1} - u_{h}^{n}, ω_{h}^{n}) + 2 τ {‖ \nabla ω_{h}^{n} ‖}^{2} = 0.$ (15)

当 $ψ = - (u_{h}^{n + 1} - u_{h}^{n - 1})$ ，得

$- (ω_{h}^{n}, u_{h}^{n + 1} - u_{h}^{n - 1}) + \frac{ε^{2}}{2} ({‖ \nabla u_{h}^{n + 1} ‖}^{2} - {‖ \nabla u_{h}^{n - 1} ‖}^{2}) + (u_{h}^{n} \frac{q_{h}^{n + 1} + q_{h}^{n - 1}}{2}, u_{h}^{n + 1} - u_{h}^{n - 1}) + \frac{β}{2 τ} {‖ u_{h}^{n + 1} - u_{h}^{n - 1} ‖}^{2} = 0.$ (16)

取 $p = τ (q_{h}^{n + 1} + q_{h}^{n - 1})$ ，可得

$\frac{1}{4} ({‖ q_{h}^{n + 1} ‖}^{2} - {‖ q_{h}^{n - 1} ‖}^{2}) - (u_{h}^{n} \frac{u_{h}^{n + 1} - u_{h}^{n - 1}}{2}, q_{h}^{n + 1} + q_{h}^{n - 1}) = 0.$ (17)

将(15)，(16)和(17)结合，可以得到

$2 τ {‖ \nabla ω_{h}^{n} ‖}^{2} + \frac{ε^{2}}{2} ({‖ \nabla u_{h}^{n + 1} ‖}^{2} - {‖ \nabla u_{h}^{n - 1} ‖}^{2}) + \frac{β}{2 τ} {‖ u_{h}^{n + 1} - u_{h}^{n - 1} ‖}^{2} + \frac{1}{4} ({‖ q_{h}^{n + 1} ‖}^{2} - {‖ q_{h}^{n - 1} ‖}^{2}) = 0.$ (18)

定理1证毕。

推论1存在常数 $C > 0$ ，设 $E^{1, 0} \leq C$ ，对任意的 $τ, h > 0$ ，有下列估计式成立

$\begin{array}{l} \max_{0 < i \leq n} ({‖ \nabla u_{h}^{i + 1} ‖}^{2} + {‖ u_{h}^{i + 1} ‖}^{2} + {‖ u_{h}^{i + 1} - u_{h}^{i} ‖}^{2}) \leq C, \\ \max_{0 < i \leq n} ({‖ \nabla q_{h}^{i + 1} ‖}^{2} + {‖ q_{h}^{i + 1} ‖}^{2}) \leq C, \\ \sum_{i = 1}^{n} τ {‖ \nabla ω_{h}^{i} ‖}^{2} \leq C . \end{array}$ (19)

4. 误差分析

在本节中，将详细给出误差估计。假设弱解具有以下正则性

$\begin{array}{l} u \in L^{\infty} (0, T; W^{1, \infty} (Ω)) \cap H^{1} (0, T; H^{r}^{+ 1} (Ω)), \\ ω \in L^{\infty} (0, T; W^{1, 6} (Ω)), \\ q \in L^{\infty} (0, T; W^{1, 6} (Ω)) . \end{array}$ (20)

其中 $r \geq 1$ 。

为了书写简便，给出如下符号：

$\begin{array}{l} {\hat{e}}_{u}^{n + 1} = R_{h} u^{n + 1} - u_{h}^{n + 1}, {\tilde{e}}_{u}^{n + 1} = u^{n + 1} - R_{h} u^{n + 1}, \\ {\hat{e}}_{ω}^{n + 1} = R_{h} ω^{n + 1} - ω_{h}^{n + 1}, {\tilde{e}}_{ω}^{n + 1} = ω^{n + 1} - R_{h} ω^{n + 1}, \\ {\hat{e}}_{q}^{n + 1} = R_{h} q^{n + 1} - q_{h}^{n + 1}, {\tilde{e}}_{q}^{n + 1} = q^{n + 1} - R_{h} q^{n + 1} . \end{array}$

定义4.1 [14] $H^{- 1}$ 范数 ${‖ \cdot ‖}_{- 1, h}$ 定义如下

${‖ ξ ‖}_{- 1, h} = \sqrt{{(ξ, ξ)}_{- 1, h}} = \sup_{0 \neq χ \in {\dot{S}}_{h}} \frac{(ξ, χ)}{‖ \nabla χ ‖}, \forall ξ \in {\dot{S}}_{h} .$ (21)

定义4.2 $Δ_{h}$ 是离散的Laplacian算子， $S_{h} \to {\dot{S}}_{h}$ 是可逆的， $\forall v_{h} \in S_{h}$ ， $Δ_{h} v_{h} \in {\dot{S}}_{h}$ ，有

$(\nabla (- Δ_{h}^{- 1}) v_{h}, \nabla χ) = (v_{h}, χ), \forall χ \in {\dot{S}}_{h} .$ (22)

引理4.1 Ritz投影算子满足下面估计

$‖ u - R_{h} u ‖ + h {‖ u - R_{h} u ‖}_{H^{1} (Ω)} \leq C h^{r + 1} {‖ u ‖}_{H^{r + 1} (Ω)}, \forall u \in H^{r + 1} (Ω) .$ (23)

引理4.2 Ritz投影算子 $R_{h} : S_{h} \to {\dot{S}}_{h}$ 满足

$(\nabla (u - R_{h} u), \nabla v) = 0, (R_{h} u - u, 1) = 0.$ (24)

引理4.3 [15] 定义如下变分问题：给出 $χ \in {\dot{S}}_{h} (Ω)$ ，求 $T_{h} (χ) \in {\dot{S}}_{h}$ 使得

$(\nabla T_{h} (χ), \nabla ν) = (χ, ν),$ (25)

其中 $T_{h} : {\dot{S}}_{h} (Ω) \to {\dot{S}}_{h} (Ω)$ 是一个可逆线性算子。令 $χ, ψ \in {\dot{S}}_{h} (Ω)$ 且有

${(χ, ψ)}_{-}_{1, h} = (\nabla T_{h} (χ), \nabla T_{h} (ψ)) = (χ, T_{h} (ψ)) = (T_{h} (χ), ψ),$ (26)

其中 ${(\cdot, \cdot)}_{-}_{1, h}$ 为定义在 ${\dot{S}}_{h} (Ω)$ 得内积。因此，对于 $\forall χ \in {\dot{S}}_{h} (Ω)$ ， $\forall g \in S_{h} (Ω)$ 。

$| (χ, g) | \leq {‖ χ ‖}_{- 1, h} ‖ \nabla g ‖ .$ (27)

定理4.1假设 $(u, q)$ ， $(u_{h}^{n + 1}, q_{h}^{n + 1})$ 分别为(5)~(7)和(11)~(13)的解，则对于任意的 $τ, h > 0$ ，存在常数 $C > 0$ 使得

$τ \sum_{i = 1}^{n} {‖ \nabla {\hat{e}}_{ω}^{i} ‖}^{2} + \frac{ε^{2}}{4 τ} \sum_{i = 1}^{n} {‖ \nabla {\hat{e}}_{u}^{i + 1} ‖}^{2} + \frac{1}{16} \sum_{i = 1}^{n} {‖ \nabla {\hat{e}}_{q}^{i + 1} - \nabla {\hat{e}}_{q}^{i - 1} ‖}^{2} + β \sum_{i = 1}^{n} {‖ δ_{τ} {\hat{e}}_{u}^{i + 1} ‖}^{2} \leq C τ^{4} + C h^{2 r} .$ (28)

证明：当 $t = n$ 时，方程(5)~(7)减去方程(11)~(13)，

$(u_{t}^{n} - \frac{u_{h}^{n + 1} - u_{h}^{n - 1}}{2 τ}, v) + (\nabla {\hat{e}}_{ω}^{n} + \nabla {\tilde{e}}_{ω}^{n}, \nabla v) = 0,$ (29)

$\begin{array}{l} ({\hat{e}}_{ω}^{n} + {\tilde{e}}_{ω}^{n}, ψ) - ε^{2} (\nabla u^{n} - \nabla \frac{u_{h}^{n + 1} + u_{h}^{n - 1}}{2}, \nabla ψ) - (u^{n} q^{n} - u_{h}^{n} \frac{q_{h}^{n + 1} + q_{h}^{n - 1}}{2}, ψ) \\ - β (u_{t}^{n} - \frac{u_{h}^{n + 1} - u_{h}^{n - 1}}{2 τ}, ψ) = 0, \end{array}$ (30)

$\frac{1}{2} (q_{t}^{n} - \frac{q_{h}^{n + 1} - q_{h}^{n - 1}}{2 τ}, p) - (u^{n} u_{t}^{n} - u_{h}^{n} \frac{u_{h}^{n + 1} - u_{h}^{n - 1}}{2 τ}, p) = 0,$ (31)

令 $v = τ {\hat{e}}_{ω}^{n}$ 可得

$(u_{t}^{n} - \frac{u_{h}^{n + 1} - u_{h}^{n - 1}}{2 τ}, τ {\hat{e}}_{ω}^{n}) + τ {‖ \nabla {\hat{e}}_{ω}^{n} ‖}^{2} = 0.$ (32)

当 $ψ = δ_{τ} {\hat{e}}_{u}^{n + 1} (δ_{τ} {\hat{e}}_{u}^{n + 1} = \frac{{\hat{e}}_{u}^{n + 1} - {\hat{e}}_{u}^{n - 1}}{2 τ})$ 时，得

$\begin{array}{l} ({\hat{e}}_{ω}^{n} + {\tilde{e}}_{ω}^{n}, δ_{τ} {\hat{e}}_{u}^{n + 1}) - ε^{2} (\nabla u^{n} - \nabla \frac{u_{h}^{n + 1} + u_{h}^{n - 1}}{2}, \nabla δ_{τ} {\hat{e}}_{u}^{n + 1}) \\ - (u^{n} q^{n} - u_{h}^{n} \frac{q_{h}^{n + 1} + q_{h}^{n - 1}}{2}, δ_{τ} {\hat{e}}_{u}^{n + 1}) - β (u_{t}^{n} - \frac{u_{h}^{n + 1} - u_{h}^{n - 1}}{2 τ}, δ_{τ} {\hat{e}}_{u}^{n + 1}) = 0. \end{array}$ (33)

其中

$\begin{array}{l} - ε^{2} (\nabla u^{n} - \nabla \frac{u_{h}^{n + 1} + u_{h}^{n - 1}}{2}, \nabla δ_{τ} {\hat{e}}_{u}^{n + 1}) \\ = \frac{ε^{2}}{2} (\nabla u^{n + 1} - 2 \nabla u^{n} + \nabla u^{n - 1}, \nabla δ_{τ} {\hat{e}}_{u}^{n + 1}) - \frac{ε^{2}}{4 τ} ({‖ \nabla {\hat{e}}_{u}^{n + 1} ‖}^{2} - {‖ \nabla {\hat{e}}_{u}^{n - 1} ‖}^{2}) . \end{array}$

取 $p = - τ Δ_{h} ({\hat{e}}_{q}^{n + 1} - {\hat{e}}_{q}^{n - 1})$ ，得

$\frac{1}{2} (q_{t}^{n} - \frac{q_{h}^{n + 1} - q_{h}^{n - 1}}{2 τ}, - τ Δ_{h} ({\hat{e}}_{q}^{n + 1} - {\hat{e}}_{q}^{n - 1})) + (u^{n} u_{t}^{n} - u_{h}^{n} \frac{u_{h}^{n + 1} - u_{h}^{n - 1}}{2 τ}, τ Δ_{h} ({\hat{e}}_{q}^{n + 1} - {\hat{e}}_{q}^{n - 1})) = 0.$ (34)

于是有

$\begin{array}{l} \frac{1}{2} {‖ \nabla {\hat{e}}_{q}^{n + 1} - \nabla {\hat{e}}_{q}^{n - 1} ‖}^{2} - \frac{1}{2} (q_{t}^{n} - \frac{q^{n + 1} - q^{n - 1}}{2 τ} + δ_{τ} {\tilde{e}}_{q}^{n + 1}, Δ_{h} ({\hat{e}}_{q}^{n + 1} - {\hat{e}}_{q}^{n - 1})) \\ + (u^{n} u_{t}^{n} - u_{h}^{n} \frac{u_{h}^{n + 1} - u_{h}^{n - 1}}{2 τ}, τ Δ_{h} ({\hat{e}}_{q}^{n + 1} - {\hat{e}}_{q}^{n - 1})) = 0. \end{array}$ (35)

将(32)，(33)，(35)相加可得

$\begin{array}{l} τ {‖ \nabla {\hat{e}}_{ω}^{n} ‖}^{2} + \frac{ε^{2}}{4 τ} ({‖ \nabla {\hat{e}}_{u}^{n + 1} ‖}^{2} - {‖ \nabla {\hat{e}}_{u}^{n - 1} ‖}^{2}) + \frac{1}{2} {‖ \nabla {\hat{e}}_{q}^{n + 1} - \nabla {\hat{e}}_{q}^{n - 1} ‖}^{2} + β {‖ δ_{τ} {\hat{e}}_{u}^{n + 1} ‖}^{2} \\ = (u_{t}^{n} - \frac{u_{h}^{n + 1} - u_{h}^{n - 1}}{2 τ}, τ {\hat{e}}_{ω}^{n}) + ({\hat{e}}_{ω}^{n} + {\tilde{e}}_{ω}^{n}, δ_{τ} {\hat{e}}_{u}^{n + 1}) + \frac{ε^{2}}{2} (\nabla u^{n + 1} - 2 \nabla u^{n} + \nabla u^{n - 1}, \nabla δ_{τ} {\hat{e}}_{u}^{n + 1}) \\ - (u^{n} q^{n} - u_{h}^{n} \frac{q_{h}^{n + 1} + q_{h}^{n - 1}}{2}, δ_{τ} {\hat{e}}_{u}^{n + 1}) - β (u_{t}^{n} - \frac{u^{n + 1} - u^{n - 1}}{2 τ} + δ_{τ} {\tilde{e}}_{u}^{n + 1}, δ_{τ} {\hat{e}}_{u}^{n + 1}) \\ + \frac{1}{2} (q_{t}^{n} - \frac{q^{n + 1} - q^{n - 1}}{2 τ} + δ_{τ} {\tilde{e}}_{q}^{n + 1}, τ Δ_{h} ({\hat{e}}_{q}^{n + 1} - {\hat{e}}_{q}^{n - 1})) \\ - (u^{n} u_{t}^{n} - u_{h}^{n} \frac{u_{h}^{n + 1} - u_{h}^{n - 1}}{2 τ}, τ Δ_{h} ({\hat{e}}_{q}^{n + 1} - {\hat{e}}_{q}^{n - 1})) \\ = \sum_{i = 1}^{7} M_{i} . \end{array}$ (36)

接下来依次估计 $M_{i} (i = 1, 2, 3, 4, 5, 6, 7)$ 。根据引理4.1~4.3，Cauchy-Schwarz和Young不等式有

$\begin{array}{l} M_{1} = (u_{t}^{n} - \frac{u_{h}^{n + 1} - u_{h}^{n - 1}}{2 τ}, τ {\hat{e}}_{ω}^{n}) \\ = (u_{t}^{n} - \frac{u^{n + 1} - u^{n - 1}}{2 τ} + \frac{u^{n + 1} - u^{n - 1}}{2 τ} - \frac{u_{h}^{n + 1} - u_{h}^{n - 1}}{2 τ}, τ {\hat{e}}_{ω}^{n}) \\ \leq τ ‖ u_{t}^{n} - \frac{u^{n + 1} - u^{n - 1}}{2 τ} ‖ ‖ {\hat{e}}_{ω}^{n} ‖ + τ ‖ δ_{τ} {\hat{e}}_{u}^{n + 1} ‖ ‖ {\hat{e}}_{ω}^{n} ‖ + τ ‖ δ_{τ} {\tilde{e}}_{u}^{n + 1} ‖ ‖ {\hat{e}}_{ω}^{n} ‖ \\ \leq C τ^{4} + C h^{2 r} + C {‖ \nabla {\hat{e}}_{ω}^{n} ‖}^{2} + \frac{1}{8} {‖ δ_{τ} {\hat{e}}_{u}^{n + 1} ‖}_{- 1, h}^{2} . \end{array}$ (37)

$M_{2} \leq ‖ {\hat{e}}_{ω}^{n} + {\tilde{e}}_{ω}^{n} ‖ ‖ δ_{τ} {\hat{e}}_{u}^{n + 1} ‖ \leq C h^{2 r} + C ‖ \nabla {\hat{e}}_{ω}^{n} ‖ + \frac{1}{8} {‖ δ_{τ} {\hat{e}}_{u}^{n + 1} ‖}_{- 1, h}^{2} .$ (38)

$M_{3} \leq \frac{ε^{2}}{2} ‖ Δ (u^{n + 1} - 2 u^{n} + u^{n - 1}) ‖ ‖ δ_{τ} {\hat{e}}_{u}^{n + 1} ‖ \leq C τ^{4} + \frac{1}{16} {‖ δ_{τ} {\hat{e}}_{u}^{n + 1} ‖}^{2} .$ (39)

化简 $M_{4}$ ，可得

$\begin{array}{l} M_{4} = - (u^{n} q^{n} - u_{h}^{n} \frac{q_{h}^{n + 1} + q_{h}^{n - 1}}{2}, δ_{τ} {\hat{e}}_{u}^{n + 1}) \\ = (u^{n} ({(u^{n})}^{2} - \frac{{(u_{h}^{n + 1})}^{2} + {(u_{h}^{n - 1})}^{2}}{2}), δ_{τ} {\hat{e}}_{u}^{n + 1}) \\ + (\frac{{(u_{h}^{n + 1})}^{2} + {(u_{h}^{n - 1})}^{2}}{2} (u^{n} - u_{h}^{n}), δ_{τ} {\hat{e}}_{u}^{n + 1}) - (u^{n} - u_{h}^{n}, δ_{τ} {\hat{e}}_{u}^{n + 1}) \\ = \sum_{i = 41}^{43} J_{i} . \end{array}$ (40)

根据(20)，引理4.1以及Young不等式可得

$\begin{array}{l} J_{41} \leq | (u^{n} ({(u^{n})}^{2} - \frac{{(u_{h}^{n + 1})}^{2} + {(u_{h}^{n - 1})}^{2}}{2}), δ_{τ} {\hat{e}}_{u}^{n + 1}) | \\ \leq ‖ \nabla (u^{n} {(u^{n})}^{2} - \frac{{(u_{h}^{n + 1})}^{2} + {(u_{h}^{n - 1})}^{2}}{2}) ‖ {‖ δ_{τ} {\hat{e}}_{u}^{n + 1} ‖}_{- 1, h} \end{array}$

$\begin{array}{l} \leq {‖ \nabla u^{n} ‖}_{L^{3}} {‖ ({(u^{n})}^{2} - \frac{{(u^{n + 1})}^{2} + {(u^{n - 1})}^{2}}{2}) ‖}_{L^{6}} {‖ δ_{τ} {\hat{e}}_{u}^{n + 1} ‖}_{- 1, h} \\ + {‖ \nabla u^{n} ‖}_{L^{3}} {‖ (\frac{{(u^{n + 1})}^{2} + {(u^{n - 1})}^{2}}{2} - \frac{{(u_{h}^{n + 1})}^{2} + {(u_{h}^{n - 1})}^{2}}{2}) ‖}_{L^{6}} {‖ δ_{τ} {\hat{e}}_{u}^{n + 1} ‖}_{- 1, h} \\ + {‖ u^{n} ‖}_{L^{\infty}} ‖ \nabla ({(u^{n})}^{2} - \frac{{(u^{n + 1})}^{2} + {(u^{n - 1})}^{2}}{2}) ‖ {‖ δ_{τ} {\hat{e}}_{u}^{n + 1} ‖}_{- 1, h} \\ + {‖ u^{n} ‖}_{L^{\infty}} ‖ \nabla (\frac{{(u^{n + 1})}^{2} + {(u^{n - 1})}^{2}}{2} - \frac{{(u_{h}^{n + 1})}^{2} + {(u_{h}^{n - 1})}^{2}}{2}) ‖ {‖ δ_{τ} {\hat{e}}_{u}^{n + 1} ‖}_{- 1, h} \\ \leq C τ^{4} + C h^{2 r} + \frac{1}{8} {‖ \nabla {\hat{e}}_{u}^{n + 1} ‖}^{2} + C {‖ {\hat{e}}_{u}^{n - 1} ‖}^{2} + 4 {‖ δ_{τ} {\hat{e}}_{u}^{n + 1} ‖}_{- 1, h}^{2} . \end{array}$ (41)

$J_{42} \leq ‖ \nabla ({\hat{e}}_{u}^{n} + {\tilde{e}}_{u}^{n}) ‖ {‖ δ_{τ} {\hat{e}}_{u}^{n + 1} ‖}_{- 1, h} \leq C h^{2 r} + C {‖ \nabla {\hat{e}}_{u}^{n} ‖}^{2} + \frac{1}{8} {‖ δ_{τ} {\hat{e}}_{u}^{n + 1} ‖}_{- 1, h}^{2} .$ (42)

$J_{43} \leq ‖ \nabla ({\hat{e}}_{u}^{n} + {\tilde{e}}_{u}^{n}) ‖ {‖ δ_{τ} {\hat{e}}_{u}^{n + 1} ‖}_{- 1, h} \leq C h^{2 r} + C {‖ \nabla {\hat{e}}_{u}^{n} ‖}^{2} + \frac{1}{8} {‖ δ_{τ} {\hat{e}}_{u}^{n + 1} ‖}_{- 1, h}^{2} .$ (43)

相加可得

$M_{4} \leq C τ^{4} + C h^{2 r} + \frac{1}{16} {‖ \nabla {\hat{e}}_{u}^{n + 1} ‖}^{2} + C {‖ \nabla {\hat{e}}_{u}^{n} ‖}^{2} + C {‖ {\hat{e}}_{u}^{n - 1} ‖}^{2} + \frac{17}{4} {‖ δ_{τ} {\hat{e}}_{u}^{n + 1} ‖}_{- 1, h}^{2} .$ (44)

对于 $M_{5}$ ，根据Young不等式可得

$\begin{array}{l} M_{5} = β (u_{t}^{n} - \frac{u^{n + 1} - u^{n + 1}}{2 τ} + δ_{τ} {\tilde{e}}_{u}^{n + 1}, δ_{τ} {\hat{e}}_{u}^{n + 1}) \\ \leq β ‖ u_{t}^{n} - \frac{u^{n + 1} - u^{n + 1}}{2 τ} ‖ ‖ δ_{τ} {\hat{e}}_{u}^{n + 1} ‖ + β ‖ δ_{τ} {\tilde{e}}_{u}^{n + 1} ‖ ‖ δ_{τ} {\hat{e}}_{u}^{n + 1} ‖ \\ \leq C τ^{4} + C h^{2 r} + \frac{1}{16 τ^{2}} {‖ δ_{τ} {\hat{e}}_{u}^{n + 1} ‖}^{2} + \frac{1}{4} {‖ δ_{τ} {\hat{e}}_{u}^{n + 1} ‖}_{- 1, h}^{2} . \end{array}$ (45)

对于 $M_{6}$ ，根据Young不等式可得

$\begin{array}{l} M_{6} \leq \frac{1}{2} τ ‖ \nabla (q_{t}^{n} - \frac{q^{n + 1} - q^{n - 1}}{2 τ}) ‖ ‖ \nabla {\hat{e}}_{q}^{n + 1} - \nabla {\hat{e}}_{q}^{n - 1} ‖ + \frac{1}{2} ‖ \nabla ({\tilde{e}}_{q}^{n + 1} - {\tilde{e}}_{q}^{n - 1}) ‖ ‖ \nabla {\hat{e}}_{q}^{n + 1} - \nabla {\hat{e}}_{q}^{n - 1} ‖ \\ \leq C τ^{4} + C h^{2 r} + \frac{1}{16} {‖ \nabla {\hat{e}}_{q}^{n + 1} - \nabla {\hat{e}}_{q}^{n - 1} ‖}^{2} . \end{array}$ (46)

化简 $M_{7}$ ，可得

$\begin{array}{l} M_{7} = - (u^{n} u_{t}^{n} - u_{h}^{n} \frac{u_{h}^{n + 1} - u_{h}^{n - 1}}{2 τ}, τ Δ_{h} ({\hat{e}}_{q}^{n + 1} - {\hat{e}}_{q}^{n - 1})) \\ = - (u^{n} (u_{t}^{n} - \frac{u^{n + 1} - u^{n - 1}}{2 τ}), τ Δ_{h} ({\hat{e}}_{q}^{n + 1} - {\hat{e}}_{q}^{n - 1})) - (\frac{u^{n + 1} - u^{n - 1}}{2 τ} (u^{n} - u_{h}^{n}), τ Δ_{h} ({\hat{e}}_{q}^{n + 1} - {\hat{e}}_{q}^{n - 1})) \\ - (u_{h}^{n} (\frac{u^{n + 1} - u^{n - 1}}{2 τ} - \frac{u_{h}^{n + 1} - u_{h}^{n - 1}}{2 τ}), τ Δ_{h} ({\hat{e}}_{q}^{n + 1} - {\hat{e}}_{q}^{n - 1})) \\ = \sum_{i = 71}^{73} J_{i}, \end{array}$ (47)

其中

$J_{71} \leq C τ^{4} + \frac{1}{16} {‖ \nabla {\hat{e}}_{q}^{n + 1} - \nabla {\hat{e}}_{q}^{n - 1} ‖}^{2} .$ (48)

$J_{72} \leq C h^{2 r} + C {‖ \nabla {\hat{e}}_{u}^{n} ‖}^{2} + \frac{1}{16} {‖ \nabla {\hat{e}}_{q}^{n + 1} - \nabla {\hat{e}}_{q}^{n - 1} ‖}^{2} .$ (49)

$\begin{array}{l} J_{73} \leq \frac{1}{2} ‖ \nabla {\hat{e}}_{u}^{n + 1} + \nabla {\tilde{e}}_{u}^{n + 1} + \nabla {\hat{e}}_{u}^{n - 1} + \nabla {\tilde{e}}_{u}^{n - 1} ‖ ‖ \nabla {\hat{e}}_{q}^{n + 1} - \nabla {\hat{e}}_{q}^{n - 1} ‖ \\ \leq C h^{2 r} + \frac{1}{4} {‖ \nabla {\hat{e}}_{u}^{n + 1} ‖}^{2} + C {‖ \nabla {\hat{e}}_{u}^{n - 1} ‖}^{2} + \frac{1}{4} {‖ \nabla {\hat{e}}_{q}^{n + 1} - \nabla {\hat{e}}_{q}^{n - 1} ‖}^{2} . \end{array}$ (50)

相加可得

$M_{7} \leq C τ^{4} + C h^{2 r} + \frac{1}{4} {‖ \nabla {\hat{e}}_{u}^{n + 1} ‖}^{2} + C {‖ \nabla {\hat{e}}_{u}^{n} ‖}^{2} + C {‖ \nabla {\hat{e}}_{u}^{n - 1} ‖}^{2} + \frac{3}{8} {‖ \nabla {\hat{e}}_{q}^{n + 1} - \nabla {\hat{e}}_{q}^{n - 1} ‖}^{2} .$ (51)

将上述不等式相加，我们得到

接下来估计

$\begin{array}{l} {‖ δ_{τ} {\hat{e}}_{u}^{n + 1} ‖}_{- 1, h}^{2} = (u_{t}^{n} - \frac{u_{h}^{n + 1} - u_{h}^{n - 1}}{2 τ}, T_{h} (δ_{τ} {\hat{e}}_{u}^{n + 1})) - (\nabla {\hat{e}}_{ω}^{n} + \nabla {\tilde{e}}_{ω}^{n}, \nabla T_{h} (δ_{τ} {\hat{e}}_{u}^{n + 1})) \\ \leq ‖ u_{t}^{n} - \frac{u^{n + 1} - u^{n - 1}}{2 τ} ‖ ‖ T_{h} (δ_{τ} {\hat{e}}_{u}^{n + 1}) ‖ + ‖ δ_{τ} {\tilde{e}}_{u}^{n + 1} ‖ ‖ T_{h} (δ_{τ} {\hat{e}}_{u}^{n + 1}) ‖ \\ + ‖ \nabla {\hat{e}}_{ω}^{n} + \nabla {\tilde{e}}_{ω}^{n} ‖ ‖ \nabla T_{h} (δ_{τ} {\hat{e}}_{u}^{n + 1}) ‖ \\ \leq C τ^{4} + C h^{2 r} + C {‖ \nabla {\hat{e}}_{ω}^{n} ‖}^{2} + \frac{3}{4} {‖ δ_{τ} {\hat{e}}_{u}^{n + 1} ‖}_{- 1, h}^{2} . \end{array}$ (53)

将上式结果代回(52)式，并从1到n求和。根据Gronwall不等式可得，

5. 数值模拟

在本节，采用数值算例验证前面理论分析的结果。具体安排如下：在第一部分，我们得到了格式(11)~(13)在不同参数下的收敛结果。第二部分给出了无条件能量耗散图，验证了系统的稳定性。最后，粘性Cahn-Hilliard模型的相分离图被给出。所有数值算例均由FreeFem++编程实现。

5.1. 收敛性

首先，给出格式(11)~(13)时间收敛阶和空间收敛阶的结果。计算区域为 $[- 1, 1] \times [- 1, 1]$ ，且初值设为

$u_{0} = 0.25 \sin (2 x) \cos (2 y) + 0.48.$

5.1.1. 时间收敛阶

在本小节中，对浓度u的时间收敛阶进行研究。表1，表2中，选择时间步长 $τ = \frac{1}{64}, \frac{1}{128}, \frac{1}{256}$ 。误差 ${‖ {\hat{e}}_{u}^{n} ‖}_{H^{1}}$ 较小，其时间收敛阶总是趋近于2，与误差分析得到的时间收敛阶一致，当变化 $ε$ 和 $β$ 时，得到的结果变化较小。

Table 1. Error and the temporal convergence ( ε = 0.3 and β = 2 , 4 )

表1. 误差和时间收敛阶( $ε = 0.3$ and $β = 2, 4$ )

Table 2. Error and the temporal convergence ( ε = 0.4 and β = 4 , 6 )

表2. 误差和时间收敛阶( $ε = 0.4$ and $β = 4, 6$ )

5.1.2. 空间收敛阶

在本小节中，对浓度u的时间收敛阶进行研究。在表3中，固定参数 $τ = 0.01$ ，网格步长为 $h = \frac{1}{8}, \frac{1}{16}, \frac{1}{64}$ 。误差 ${‖ {\hat{e}}_{u} ‖}_{H^{1}}$ 较小，其空间收敛阶总是趋近于2，与误差分析得到的时间收敛阶一致，当变化 $ε$ 和 $β$ 时，得到的结果变化较小。

Table 3. Error and the special convergence ( β = 0.01 and ε = 0.2 , 0.3 )

表3. 误差和空间收敛阶( $β = 0.01$ and $ε = 0.2, 0.3$ )

5.2. 无条件能量稳定

接着，验证了上述所提格式的稳定性。粘性Cahn-Hilliard模型(3)的能量泛函(4)被离散为

$E (u_{h}^{n + 1}) = \int_{Ω} (\frac{ε^{2}}{2} {| \nabla u_{h}^{n + 1} |}^{2} + \frac{1}{4} {| q_{h}^{n + 1} |}^{2}) d x$

且格式(11)~(13)修正的能量为

$E_{h}^{n + 1, n} = \frac{ε^{2}}{2} ({‖ \nabla u_{h}^{n + 1} ‖}^{2} + {‖ u_{h}^{n} ‖}^{2}) + \frac{1}{4} ({‖ q_{h}^{n + 1} ‖}^{2} + {‖ q_{h}^{n} ‖}^{2}) .$

选择参数 $T = 1, ε = 0.1$ 。当 $β = 0.1$ 时，图1表明能量由最初递减到逐渐趋于稳定。当 $β = 0.01$ ，从图1中发现，随着时间的推移，能量仍然是由最初递减到逐渐趋于稳定。还可以看到不同的 $β$ 对结果的影响微乎其微，从而验证了格式(11)~(13)的能量是耗散的，通过该数值算例验证了定理3.1的有效性。

Figure 1. Diagram of the practical teaching system of automation major

图1. 离散的能量演化

5.3. 相分离

在本小节中，为了更好的观察数值解的变化过程，模拟粘性Cahn-Hilliard模型的相分离。计算区域为 $[0, 1] \times [0, 1]$ ，模拟参数为 $ε = 0.1$ 和 $β = 1$ ，初始条件为

$u_{0} = 2 r a n d () - 1,$

其中 $r a n d ()$ 是 $[0, 1]$ 的一个随机数。从图2给出了不同时刻下相分离过程的快照。图中表明随着时间的增加，相分离现象越来越明显。

t = 0.01 t = 0.05 t = 0.1 t = 0.5 t = 1 t = 2

Figure 2. Snapshot of coarsening

图2. 粗化快照

6. 结论

本文主要提出了一种基于拉格朗日乘子法的粘性Cahn-Hilliard方程的能量稳定的二阶数值格式，然后证明数值格式的稳定性和误差分析。最后通过不同的数值实验验证了前面理论的有效性，时间收敛阶与空间收敛阶表明 ${‖ {\hat{e}}_{u} ‖}_{H^{1}}$ 的收敛阶接近于2，粗化图显示了u随着时间增长的演化，能量耗散的数值结果表明能量随着时间的增加而减少，最后趋于稳定。

致谢

在此感谢山西省回国留学人员科研资助项目(2021-029)，山西省科技合作交流专项项目(202104041101019)，山西省自然科学基金面上项目(202203021211129)的支持。

基金项目

山西省回国留学人员科研资助项目(2021-029)，山西省科技合作交流专项项目(202104041101019)，山西省自然科学基金面上项目(202203021211129)。

文章引用

刘静,王旦霞,王李靖. 粘性Cahn-Hilliard方程的二阶稳定Crank-Nicolson-Leapfrog格式
A Stable Second-Order Crank-Nicolson-Leapfrog Scheme for the Viscous Cahn-Hilliard Equation[J]. 应用数学进展, 2023, 12(01): 339-351. https://doi.org/10.12677/AAM.2023.121037

参考文献

1. Cahn, J.W. and Hilliard, J.E. (1958) Free Energy of a Nonuniform System. I. Interfacial Free Energy. Journal of Chemi-cal Physics, 28, 258-267. https://doi.org/10.1063/1.1744102

2. Puri, S. and Binder, K. (1991) Phenomenological Theory for the Formation of Interfaces via the Interdiffusion of Layers. Physical Review B, 44, 9735-9738. https://doi.org/10.1103/PhysRevB.44.9735

3. Jabbari, E. and Peppas, N.A. (1995) A Model for Interdiffusion at In-terfaces of Polymers with Dissimilar Physical Properties. Polymer, 36, 575-586. https://doi.org/10.1016/0032-3861(95)91567-Q

4. 叶兴德, 程晓良. Cahn-Hilliard方程的Legendre谱逼近[J]. 计算数学, 2003, 25(2): 157-170.

5. Novick-Cohen, A. (1988) On the Viscous Cahn-Hilliard Equation. In: Ball, J.M., Ed., Material Instabilities in Continuum and Related Mathematical Problem, Oxford University Press, Oxford.

6. 李亚楠, 王旦霞, 任永华. 具有浓度迁移率和对数势能的粘性Cahn-Hilliard方程的有限元算法[J]. 数学的实践与认识, 2021, 51(2): 241-250.

7. Yang, X.F., Zhao, J. and He, X.M. (2018) Linear, Second Order and Unconditionally Energy Stable Schemes for the Viscous Cahn-Hilliard Equation with Hyperbolic Relaxation Using the Invariant Energy Quadratization Method. Journal of Computational and Applied Mathematics, 343, 80-97. https://doi.org/10.1016/j.cam.2018.04.027

8. 李娟. 粘性Cahn-Hilliard方程的高精度线性化差分方法[J]. 西南大学学报(自然科学版), 2020, 42(1): 51-58.

9. 王志丽, 王旦霞, 贾宏恩. 具有对数势能的粘性Cahn-Hilliard方程的有限元算法[J]. 重庆师范大学学报(自然科学版), 2021, 38(5): 81-89.

10. Elliott, C.M. and Stuart, A.M. (1993) The Global Dynamics of Discrete Semi-Linear Parabolic Equations. SIAM Journal on Numerical Analysis, 30, 1622-1663. https://doi.org/10.1137/0730084

11. Shen, J., Xu, J. and Yang, J. (2018) The Scalar Auxiliary Varia-ble (SAV) Approach for Gradient Flows. Journal of Computational Physics, 353, 407-416. https://doi.org/10.1016/j.jcp.2017.10.021

12. Yang, X.F. (2016) Linear, First and Second-Order, Unconditionally Energy Stable Numerical Schemes for the Phase Field Model of Homopolymer Blends. Journal of Computational Phys-ics, 327, 294-316. https://doi.org/10.1016/j.jcp.2016.09.029

13. Guo, Y.Y., Jia, H.E., Li, J.C. and Li, M. (2020) Numerical Analysis for the Cahn-Hilliard-Hele-Shaw System with Variable Mobility and Logarithmic Flory-Huggins Potential. Applied Nu-merical Mathematics, 150, 206-221. https://doi.org/10.1016/j.apnum.2019.09.014

14. Wang, S.Q. and Shi, Q. (1993) Interdiffusion in Binary Polymer Mixtures. Macromolecules, 26, 1091-1096. https://doi.org/10.1021/ma00057a033

15. Liu, Y., Chen, W.B., Wang, C. and Wise, S.M. (2017) Error Analysis of a Mixed Finite Element Method for a Cahn-Hilliard-Hele-Shaw System. Numerische Mathematik, 135, 679-709. https://doi.org/10.1007/s00211-016-0813-2

NOTES

^*通讯作者。

期刊菜单