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

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

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

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.

1. 引言

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

本文研究粘性Cahn-Hilliard方程

$\{\begin{array}{l}{u}_t=\Delta \omega ,\\ \omega =-{\epsilon }^2\Delta u+F^{\prime }\left(u\right)+\beta u_t,\\ {\partial }_{n}u={\partial }_n\omega =0,\\ u\left(x,0\right)=u_0\left(x\right),\end{array}$ (1)

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

$E\left(u\right)={\displaystyle {\int }_{\Omega }\left(\frac{{\epsilon }^2}{2}{\left|\nabla u\right|}^2+F\left(u\right)\right)\text{d}x}.$ (2)

当 $\beta =0$ 时，方程(1.1)就变为经典的Cahn-Hilliard方程。这类四阶非线性扩散方程在理论和数值方面都已经被广泛地研究。当 $\beta >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. 预备知识

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

$\left(u,v\right)={\displaystyle {\int }_{\Omega }uv\text{d}\Omega },\Vert u\Vert =\sqrt{\left(u,u\right)},\forall u,v\in L^2\left(\Omega \right).$

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

${\Vert u\Vert }_{H^k}={\left({\displaystyle \underset{0\le \alpha \le k}{\sum }{\Vert D^{\alpha }u\Vert }^2}\right)}^{\frac{1}{2}},\forall u\in H^k\left(\Omega \right).$

2.2. 半离散数值格式

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

$\begin{array}{l}{u}_t=\Delta \omega ,\\ \omega =uq-{\epsilon }^2\Delta u+\beta u_t,\\ \frac{1}{2}{q}_t=u{u}_t.\end{array}$ (3)

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

$E\left(u\right)={\displaystyle {\int }_{\Omega }\left(\frac{{\epsilon }^2}{2}{\left|\nabla u\right|}^2+\frac{1}{4}{\left|q\right|}^2\right)\text{d}x}.$ (4)

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

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

$\left(u_t,v\right)+\left(\nabla \omega ,\nabla v\right)=0,$ (5)

$\left(\omega ,\psi \right)-{\epsilon }^2\left(\nabla u,\nabla \psi \right)-\left(uq,\psi \right)-\beta \left(u_t,\psi \right)=0,$ (6)

$\left(\frac{1}{2}{q}_t,p\right)-\left(u{u}_t,p\right)=0.$ (7)

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

$\left(\frac{u^n{}^{+1}-u^n{}^{-1}}{2\tau },v\right)+\left(\nabla {\omega }^n,\nabla v\right)=0,$ (8)

$\left({\omega }^n,\psi \right)-{\epsilon }^2\left(\nabla \frac{u^n{}^{+1}+u^n{}^{-1}}{2},\nabla \psi \right)-\left(u^n\frac{q^n{}^{+1}+q^n{}^{-1}}{2},\psi \right)-\beta \left(\frac{u^n{}^{+1}-u^n{}^{-1}}{2\tau },\psi \right)=0,$ (9)

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

2.3. 全离散数值格式

在区域 $\Omega$ 上作拟一致剖分，记作 $T_h=K$， $h_i$ 为网格大小， $h=\underset{0\le i\le N}{\max }{h}_i$， $S_h$ 是分片连续的有限元空间，定义如下

$S_h=\left\{v_h\in C\left(\Omega \right)|{v_h|}_K\in P_k\left(x,y\right),K\in T_h\right\}\subset H^1\left(\Omega \right),$

其中 $P_k\left(x,y\right)$ 是 $x,y$ 的次数不超过 $k\in Z^{+}$ 的线性多项式的集合。定义 $L_0^2=\left\{u\in L^2\left(\Omega \right)|\left(u,1\right)=0\right\}$， ${\stackrel{\dot{}}{S}}_h=S_h\cap L_0^2\left(\Omega \right)$。接下来构造模型(5)~(7)的全离散格式，给定 $\left(u_h^{n-1},u_h^n\right)$，当 $n\ge 1$ 时，求 $\left(u_h^{n+1},q_h^{n+1}\right)$

$\left(\frac{u_h^{n+1}-u_h^{n-1}}{2\tau },v\right)+\left(\nabla {\omega }_h^n,\nabla v\right)=0,$ (11)

$\left({\omega }_h^n,\psi \right)-{\epsilon }^2\left(\nabla \frac{u_h^{n+1}+u_h^{n-1}}{2},\nabla \psi \right)-\left(u_h^n\frac{q_h^{n+1}+q_h^{n-1}}{2},\psi \right)-\beta \left(\frac{u_h^{n+1}-u_h^{n-1}}{2\tau },\psi \right)=0,$ (12)

$\frac{1}{2}\left(\frac{q_h^{n+1}-q_h^{n-1}}{2\tau },p\right)-\left(u_h^n\frac{u_h^{n+1}-u_h^{n-1}}{2\tau },p\right)=0$.(13)

3. 稳定性分析

定理3.1令 $\left(u_h^{n+1},q_h^{n+1}\right)$ 是(11)~(13)的解，定义

$E^{n+1,n}=\frac{{\epsilon }^2}{2}\left({\Vert \nabla u_h^{n+1}\Vert }^2+{\Vert \nabla u_h^n\Vert }^2\right)+\frac{1}{4}\left({\Vert q_h^{n+1}\Vert }^2+{\Vert q_h^n\Vert }^2\right).$

对任意的 $\tau ,h,\epsilon >0,n\ge 1$，有下述不等式成立

$E^{n+1,n}+\frac{\beta }{2\tau }{\Vert u_h^{n+1}-u_h^{n-1}\Vert }^2+2\tau {\Vert \nabla {\omega }_h^n\Vert }^2=E^{n,n-1}.$ (14)

证明：令 $v=2\tau {\omega }_h^n$，得

$\left(u_h^{n+1}-u_h^n,{\omega }_h^n\right)+2\tau {\Vert \nabla {\omega }_h^n\Vert }^2=0.$ (15)

当 $\psi =-\left(u_h^{n+1}-u_h^{n-1}\right)$，得

$-\left({\omega }_h^n,u_h^{n+1}-u_h^{n-1}\right)+\frac{{\epsilon }^2}{2}\left({\Vert \nabla u_h^{n+1}\Vert }^2-{\Vert \nabla u_h^{n-1}\Vert }^2\right)+\left(u_h^n\frac{q_h^{n+1}+q_h^{n-1}}{2},u_h^{n+1}-u_h^{n-1}\right)+\frac{\beta }{2\tau }{\Vert u_h^{n+1}-u_h^{n-1}\Vert }^2=0.$ (16)

取 $p=\tau \left(q_h^{n+1}+q_h^{n-1}\right)$，可得

$\frac{1}{4}\left({\Vert q_h^{n+1}\Vert }^2-{\Vert q_h^{n-1}\Vert }^2\right)-\left(u_h^n\frac{u_h^{n+1}-u_h^{n-1}}{2},q_h^{n+1}+q_h^{n-1}\right)=0.$ (17)

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

$2\tau {\Vert \nabla {\omega }_h^n\Vert }^2+\frac{{\epsilon }^2}{2}\left({\Vert \nabla u_h^{n+1}\Vert }^2-{\Vert \nabla u_h^{n-1}\Vert }^2\right)+\frac{\beta }{2\tau }{\Vert u_h^{n+1}-u_h^{n-1}\Vert }^2+\frac{1}{4}\left({\Vert q_h^{n+1}\Vert }^2-{\Vert q_h^{n-1}\Vert }^2\right)=0.$ (18)

定理1证毕。

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

$\begin{array}{l}\underset{0<i\le n}{\max }\left({\Vert \nabla u_h^{i+1}\Vert }^2+{\Vert u_h^{i+1}\Vert }^2+{\Vert u_h^{i+1}-u_h^i\Vert }^2\right)\le C,\\ \underset{0<i\le n}{\max }\left({\Vert \nabla q_h^{i+1}\Vert }^2+{\Vert q_h^{i+1}\Vert }^2\right)\le C,\\ {\displaystyle \underset{i=1}{\overset{n}{\sum }}\tau {\Vert \nabla {\omega }_h^i\Vert }^2}\le C.\end{array}$ (19)

4. 误差分析

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

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

其中 $r\ge 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}}_{\omega }^{n+1}=R_h{\omega }^{n+1}-{\omega }_h^{n+1},{\tilde{e}}_{\omega }^{n+1}={\omega }^{n+1}-R_h{\omega }^{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}$ 范数 ${\Vert \cdot \Vert }_{-1,h}$ 定义如下

${\Vert \xi \Vert }_{-1,h}=\sqrt{{\left(\xi ,\xi \right)}_{-1,h}}=\underset{0\ne \chi \in {\stackrel{\dot{}}{S}}_h}{\sup }\frac{\left(\xi ,\chi \right)}{\Vert \nabla \chi \Vert },\forall \xi \in {\stackrel{\dot{}}{S}}_h.$ (21)

定义4.2 ${\Delta }_h$ 是离散的Laplacian算子， $S_h\to {\stackrel{\dot{}}{S}}_h$ 是可逆的， $\forall v_h\in S_h$， ${\Delta }_h{v}_h\in {\stackrel{\dot{}}{S}}_h$，有

$\left(\nabla \left(-{\Delta }_h^{-1}\right)v_h,\nabla \chi \right)=\left(v_h,\chi \right),\forall \chi \in {\stackrel{\dot{}}{S}}_h.$ (22)

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

$\Vert u-R_{h}u\Vert +h{\Vert u-R_{h}u\Vert }_{H^1\left(\Omega \right)}\le C{h}^{r+1}{\Vert u\Vert }_{H^{r+1}\left(\Omega \right)},\forall u\in H^{r+1}\left(\Omega \right).$ (23)

引理4.2 Ritz投影算子 $R_h:S_h\to {\stackrel{\dot{}}{S}}_h$ 满足

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

引理4.3 [15] 定义如下变分问题：给出 $\chi \in {\stackrel{\dot{}}{S}}_h\left(\Omega \right)$，求 $T_h\left(\chi \right)\in {\stackrel{\dot{}}{S}}_h$ 使得

$\left(\nabla T_h\left(\chi \right),\nabla \nu \right)=\left(\chi ,\nu \right),$ (25)

其中 $T_h:{\stackrel{\dot{}}{S}}_h\left(\Omega \right)\to {\stackrel{\dot{}}{S}}_h\left(\Omega \right)$ 是一个可逆线性算子。令 $\chi ,\psi \in {\stackrel{\dot{}}{S}}_h\left(\Omega \right)$ 且有

${\left(\chi ,\psi \right)}_{-}{}_{1,h}=\left(\nabla T_h\left(\chi \right),\nabla T_h\left(\psi \right)\right)=\left(\chi ,T_h\left(\psi \right)\right)=\left(T_h\left(\chi \right),\psi \right),$ (26)

其中 ${\left(\cdot ,\cdot \right)}_{-}{}_{1,h}$ 为定义在 ${\stackrel{\dot{}}{S}}_h\left(\Omega \right)$ 得内积。因此，对于 $\forall \chi \in {\stackrel{\dot{}}{S}}_h\left(\Omega \right)$， $\forall g\in S_h\left(\Omega \right)$。

$\left|\left(\chi ,g\right)\right|\le {\Vert \chi \Vert }_{-1,h}\Vert \nabla g\Vert .$ (27)

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

$\tau {\displaystyle \underset{i=1}{\overset{n}{\sum }}{\Vert \nabla {\hat{e}}_{\omega }^i\Vert }^2}+\frac{{\epsilon }^2}{4\tau }{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\Vert \nabla {\hat{e}}_u^{i+1}\Vert }^2}+\frac{1}{16}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\Vert \nabla {\hat{e}}_q^{i+1}-\nabla {\hat{e}}_q^{i-1}\Vert }^2}+\beta {\displaystyle \underset{i=1}{\overset{n}{\sum }}{\Vert {\delta }_{\tau }{\hat{e}}_u^{i+1}\Vert }^2}\le C{\tau }^4+C{h}^{2r}.$ (28)

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

$\left(u_t^n-\frac{u_h^{n+1}-u_h^{n-1}}{2\tau },v\right)+\left(\nabla {\hat{e}}_{\omega }^n+\nabla {\tilde{e}}_{\omega }^n,\nabla v\right)=0,$ (29)

$\begin{array}{l}\left({\hat{e}}_{\omega }^n+{\tilde{e}}_{\omega }^n,\psi \right)-{\epsilon }^2\left(\nabla u^n-\nabla \frac{u_h^{n+1}+u_h^{n-1}}{2},\nabla \psi \right)-\left(u^n{q}^n-u_h^n\frac{q_h^{n+1}+q_h^{n-1}}{2},\psi \right)\\ -\beta \left(u_t^n-\frac{u_h^{n+1}-u_h^{n-1}}{2\tau },\psi \right)=0,\end{array}$ (30)

$\frac{1}{2}\left(q_t^n-\frac{q_h^{n+1}-q_h^{n-1}}{2\tau },p\right)-\left(u^n{u}_t^n-u_h^n\frac{u_h^{n+1}-u_h^{n-1}}{2\tau },p\right)=0,$ (31)

令 $v=\tau {\hat{e}}_{\omega }^n$ 可得

$\left(u_t^n-\frac{u_h^{n+1}-u_h^{n-1}}{2\tau },\tau {\hat{e}}_{\omega }^n\right)+\tau {\Vert \nabla {\hat{e}}_{\omega }^n\Vert }^2=0.$ (32)

当 $\psi ={\delta }_{\tau }{\hat{e}}_u^{n+1}\left({\delta }_{\tau }{\hat{e}}_u^{n+1}=\frac{{\hat{e}}_u^{n+1}-{\hat{e}}_u^{n-1}}{2\tau }\right)$ 时，得

$\begin{array}{l}\left({\hat{e}}_{\omega }^n+{\tilde{e}}_{\omega }^n,{\delta }_{\tau }{\hat{e}}_u^{n+1}\right)-{\epsilon }^2\left(\nabla u^n-\nabla \frac{u_h^{n+1}+u_h^{n-1}}{2},\nabla {\delta }_{\tau }{\hat{e}}_u^{n+1}\right)\\ -\left(u^n{q}^n-u_h^n\frac{q_h^{n+1}+q_h^{n-1}}{2},{\delta }_{\tau }{\hat{e}}_u^{n+1}\right)-\beta \left(u_t^n-\frac{u_h^{n+1}-u_h^{n-1}}{2\tau },{\delta }_{\tau }{\hat{e}}_u^{n+1}\right)=0.\end{array}$ (33)

其中

$\begin{array}{l}-{\epsilon }^2\left(\nabla u^n-\nabla \frac{u_h^{n+1}+u_h^{n-1}}{2},\nabla {\delta }_{\tau }{\hat{e}}_u^{n+1}\right)\\ =\frac{{\epsilon }^2}{2}\left(\nabla u^{n+1}-2\nabla u^n+\nabla u^{n-1},\nabla {\delta }_{\tau }{\hat{e}}_u^{n+1}\right)-\frac{{\epsilon }^2}{4\tau }\left({\Vert \nabla {\hat{e}}_u^{n+1}\Vert }^2-{\Vert \nabla {\hat{e}}_u^{n-1}\Vert }^2\right).\end{array}$

取 $p=-\tau {\Delta }_h\left({\hat{e}}_q^{n+1}-{\hat{e}}_q^{n-1}\right)$，得

$\frac{1}{2}\left(q_t^n-\frac{q_h^{n+1}-q_h^{n-1}}{2\tau },-\tau {\Delta }_h\left({\hat{e}}_q^{n+1}-{\hat{e}}_q^{n-1}\right)\right)+\left(u^n{u}_t^n-u_h^n\frac{u_h^{n+1}-u_h^{n-1}}{2\tau },\tau {\Delta }_h\left({\hat{e}}_q^{n+1}-{\hat{e}}_q^{n-1}\right)\right)=0.$ (34)

于是有

$\begin{array}{l}\frac{1}{2}{\Vert \nabla {\hat{e}}_q^{n+1}-\nabla {\hat{e}}_q^{n-1}\Vert }^2-\frac{1}{2}\left(q_t^n-\frac{q^{n+1}-q^{n-1}}{2\tau }+{\delta }_{\tau }{\tilde{e}}_q^{n+1},{\Delta }_h\left({\hat{e}}_q^{n+1}-{\hat{e}}_q^{n-1}\right)\right)\\ +\left(u^n{u}_t^n-u_h^n\frac{u_h^{n+1}-u_h^{n-1}}{2\tau },\tau {\Delta }_h\left({\hat{e}}_q^{n+1}-{\hat{e}}_q^{n-1}\right)\right)=0.\end{array}$ (35)

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

$\begin{array}{l}\tau {\Vert \nabla {\hat{e}}_{\omega }^n\Vert }^2+\frac{{\epsilon }^2}{4\tau }\left({\Vert \nabla {\hat{e}}_u^{n+1}\Vert }^2-{\Vert \nabla {\hat{e}}_u^{n-1}\Vert }^2\right)+\frac{1}{2}{\Vert \nabla {\hat{e}}_q^{n+1}-\nabla {\hat{e}}_q^{n-1}\Vert }^2+\beta {\Vert {\delta }_{\tau }{\hat{e}}_u^{n+1}\Vert }^2\\ =\left(u_t^n-\frac{u_h^{n+1}-u_h^{n-1}}{2\tau },\tau {\hat{e}}_{\omega }^n\right)+\left({\hat{e}}_{\omega }^n+{\tilde{e}}_{\omega }^n,{\delta }_{\tau }{\hat{e}}_u^{n+1}\right)+\frac{{\epsilon }^2}{2}\left(\nabla u^{n+1}-2\nabla u^n+\nabla u^{n-1},\nabla {\delta }_{\tau }{\hat{e}}_u^{n+1}\right)\\ -\left(u^n{q}^n-u_h^n\frac{q_h^{n+1}+q_h^{n-1}}{2},{\delta }_{\tau }{\hat{e}}_u^{n+1}\right)-\beta \left(u_t^n-\frac{u^{n+1}-u^{n-1}}{2\tau }+{\delta }_{\tau }{\tilde{e}}_u^{n+1},{\delta }_{\tau }{\hat{e}}_u^{n+1}\right)\\ +\frac{1}{2}\left(q_t^n-\frac{q^{n+1}-q^{n-1}}{2\tau }+{\delta }_{\tau }{\tilde{e}}_q^{n+1},\tau {\Delta }_h\left({\hat{e}}_q^{n+1}-{\hat{e}}_q^{n-1}\right)\right)\\ -\left(u^n{u}_t^n-u_h^n\frac{u_h^{n+1}-u_h^{n-1}}{2\tau },\tau {\Delta }_h\left({\hat{e}}_q^{n+1}-{\hat{e}}_q^{n-1}\right)\right)\\ ={\displaystyle \underset{i=1}{\overset{7}{\sum }}{M}_i}.\end{array}$ (36)

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

$\begin{array}{l}{M}_1=\left(u_t^n-\frac{u_h^{n+1}-u_h^{n-1}}{2\tau },\tau {\hat{e}}_{\omega }^n\right)\\ =\left(u_t^n-\frac{u^{n+1}-u^{n-1}}{2\tau }+\frac{u^{n+1}-u^{n-1}}{2\tau }-\frac{u_h^{n+1}-u_h^{n-1}}{2\tau },\tau {\hat{e}}_{\omega }^n\right)\\ \le \tau \Vert u_t^n-\frac{u^{n+1}-u^{n-1}}{2\tau }\Vert \Vert {\hat{e}}_{\omega }^n\Vert +\tau \Vert {\delta }_{\tau }{\hat{e}}_u^{n+1}\Vert \Vert {\hat{e}}_{\omega }^n\Vert +\tau \Vert {\delta }_{\tau }{\tilde{e}}_u^{n+1}\Vert \Vert {\hat{e}}_{\omega }^n\Vert \\ \le C{\tau }^4+C{h}^{2r}+C{\Vert \nabla {\hat{e}}_{\omega }^n\Vert }^2+\frac{1}{8}{\Vert {\delta }_{\tau }{\hat{e}}_u^{n+1}\Vert }_{-1,h}^2.\end{array}$ (37)

$M_2\le \Vert {\hat{e}}_{\omega }^n+{\tilde{e}}_{\omega }^n\Vert \Vert {\delta }_{\tau }{\hat{e}}_u^{n+1}\Vert \le C{h}^{2r}+C\Vert \nabla {\hat{e}}_{\omega }^n\Vert +\frac{1}{8}{\Vert {\delta }_{\tau }{\hat{e}}_u^{n+1}\Vert }_{-1,h}^2.$ (38)

$M_3\le \frac{{\epsilon }^2}{2}\Vert \Delta \left(u^{n+1}-2u^n+u^{n-1}\right)\Vert \Vert {\delta }_{\tau }{\hat{e}}_u^{n+1}\Vert \le C{\tau }^4+\frac{1}{16}{\Vert {\delta }_{\tau }{\hat{e}}_u^{n+1}\Vert }^2.$ (39)

化简 $M_4$，可得

$\begin{array}{l}{M}_4=-\left(u^n{q}^n-u_h^n\frac{q_h^{n+1}+q_h^{n-1}}{2},{\delta }_{\tau }{\hat{e}}_u^{n+1}\right)\\ =\left(u^n\left({\left(u^n\right)}^2-\frac{{\left(u_h^{n+1}\right)}^2+{\left(u_h^{n-1}\right)}^2}{2}\right),{\delta }_{\tau }{\hat{e}}_u^{n+1}\right)\\ +\left(\frac{{\left(u_h^{n+1}\right)}^2+{\left(u_h^{n-1}\right)}^2}{2}\left(u^n-u_h^n\right),{\delta }_{\tau }{\hat{e}}_u^{n+1}\right)-\left(u^n-u_h^n,{\delta }_{\tau }{\hat{e}}_u^{n+1}\right)\\ ={\displaystyle \underset{i=41}{\overset{43}{\sum }}{J}_i}.\end{array}$ (40)

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

$\begin{array}{l}{J}_{41}\le \left|\left(u^n\left({\left(u^n\right)}^2-\frac{{\left(u_h^{n+1}\right)}^2+{\left(u_h^{n-1}\right)}^2}{2}\right),{\delta }_{\tau }{\hat{e}}_u^{n+1}\right)\right|\\ \le \Vert \nabla \left(u^n{\left(u^n\right)}^2-\frac{{\left(u_h^{n+1}\right)}^2+{\left(u_h^{n-1}\right)}^2}{2}\right)\Vert {\Vert {\delta }_{\tau }{\hat{e}}_u^{n+1}\Vert }_{-1,h}\end{array}$

$\begin{array}{l}\le {\Vert \nabla u^n\Vert }_{L^3}{\Vert \left({\left(u^n\right)}^2-\frac{{\left(u^{n+1}\right)}^2+{\left(u^{n-1}\right)}^2}{2}\right)\Vert }_{L^6}{\Vert {\delta }_{\tau }{\hat{e}}_u^{n+1}\Vert }_{-1,h}\\ +{\Vert \nabla u^n\Vert }_{L^3}{\Vert \left(\frac{{\left(u^{n+1}\right)}^2+{\left(u^{n-1}\right)}^2}{2}-\frac{{\left(u_h^{n+1}\right)}^2+{\left(u_h^{n-1}\right)}^2}{2}\right)\Vert }_{L^6}{\Vert {\delta }_{\tau }{\hat{e}}_u^{n+1}\Vert }_{-1,h}\\ +{\Vert u^n\Vert }_{L^{\infty }}\Vert \nabla \left({\left(u^n\right)}^2-\frac{{\left(u^{n+1}\right)}^2+{\left(u^{n-1}\right)}^2}{2}\right)\Vert {\Vert {\delta }_{\tau }{\hat{e}}_u^{n+1}\Vert }_{-1,h}\\ +{\Vert u^n\Vert }_{L^{\infty }}\Vert \nabla \left(\frac{{\left(u^{n+1}\right)}^2+{\left(u^{n-1}\right)}^2}{2}-\frac{{\left(u_h^{n+1}\right)}^2+{\left(u_h^{n-1}\right)}^2}{2}\right)\Vert {\Vert {\delta }_{\tau }{\hat{e}}_u^{n+1}\Vert }_{-1,h}\\ \le C{\tau }^4+C{h}^{2r}+\frac{1}{8}{\Vert \nabla {\hat{e}}_u^{n+1}\Vert }^2+C{\Vert {\hat{e}}_u^{n-1}\Vert }^2+4{\Vert {\delta }_{\tau }{\hat{e}}_u^{n+1}\Vert }_{-1,h}^2.\end{array}$ (41)

$J_{42}\le \Vert \nabla \left({\hat{e}}_u^n+{\tilde{e}}_u^n\right)\Vert {\Vert {\delta }_{\tau }{\hat{e}}_u^{n+1}\Vert }_{-1,h}\le C{h}^{2r}+C{\Vert \nabla {\hat{e}}_u^n\Vert }^2+\frac{1}{8}{\Vert {\delta }_{\tau }{\hat{e}}_u^{n+1}\Vert }_{-1,h}^2.$ (42)

$J_{43}\le \Vert \nabla \left({\hat{e}}_u^n+{\tilde{e}}_u^n\right)\Vert {\Vert {\delta }_{\tau }{\hat{e}}_u^{n+1}\Vert }_{-1,h}\le C{h}^{2r}+C{\Vert \nabla {\hat{e}}_u^n\Vert }^2+\frac{1}{8}{\Vert {\delta }_{\tau }{\hat{e}}_u^{n+1}\Vert }_{-1,h}^2.$ (43)

相加可得

$M_4\le C{\tau }^4+C{h}^{2r}+\frac{1}{16}{\Vert \nabla {\hat{e}}_u^{n+1}\Vert }^2+C{\Vert \nabla {\hat{e}}_u^n\Vert }^2+C{\Vert {\hat{e}}_u^{n-1}\Vert }^2+\frac{17}{4}{\Vert {\delta }_{\tau }{\hat{e}}_u^{n+1}\Vert }_{-1,h}^2.$ (44)

对于 $M_5$，根据Young不等式可得

$\begin{array}{l}{M}_5=\beta \left(u_t^n-\frac{u^{n+1}-u^{n+1}}{2\tau }+{\delta }_{\tau }{\tilde{e}}_u^{n+1},{\delta }_{\tau }{\hat{e}}_u^{n+1}\right)\\ \le \beta \Vert u_t^n-\frac{u^{n+1}-u^{n+1}}{2\tau }\Vert \Vert {\delta }_{\tau }{\hat{e}}_u^{n+1}\Vert +\beta \Vert {\delta }_{\tau }{\tilde{e}}_u^{n+1}\Vert \Vert {\delta }_{\tau }{\hat{e}}_u^{n+1}\Vert \\ \le C{\tau }^4+C{h}^{2r}+\frac{1}{16{\tau }^2}{\Vert {\delta }_{\tau }{\hat{e}}_u^{n+1}\Vert }^2+\frac{1}{4}{\Vert {\delta }_{\tau }{\hat{e}}_u^{n+1}\Vert }_{-1,h}^2.\end{array}$ (45)

对于 $M_6$，根据Young不等式可得

$\begin{array}{l}{M}_6\le \frac{1}{2}\tau \Vert \nabla \left(q_t^n-\frac{q^{n+1}-q^{n-1}}{2\tau }\right)\Vert \Vert \nabla {\hat{e}}_q^{n+1}-\nabla {\hat{e}}_q^{n-1}\Vert +\frac{1}{2}\Vert \nabla \left({\tilde{e}}_q^{n+1}-{\tilde{e}}_q^{n-1}\right)\Vert \Vert \nabla {\hat{e}}_q^{n+1}-\nabla {\hat{e}}_q^{n-1}\Vert \\ \le C{\tau }^4+C{h}^{2r}+\frac{1}{16}{\Vert \nabla {\hat{e}}_q^{n+1}-\nabla {\hat{e}}_q^{n-1}\Vert }^2.\end{array}$ (46)

化简 $M_7$，可得

$\begin{array}{l}{M}_7=-\left(u^n{u}_t^n-u_h^n\frac{u_h^{n+1}-u_h^{n-1}}{2\tau },\tau {\Delta }_h\left({\hat{e}}_q^{n+1}-{\hat{e}}_q^{n-1}\right)\right)\\ =-\left(u^n\left(u_t^n-\frac{u^{n+1}-u^{n-1}}{2\tau }\right),\tau {\Delta }_h\left({\hat{e}}_q^{n+1}-{\hat{e}}_q^{n-1}\right)\right)-\left(\frac{u^{n+1}-u^{n-1}}{2\tau }\left(u^n-u_h^n\right),\tau {\Delta }_h\left({\hat{e}}_q^{n+1}-{\hat{e}}_q^{n-1}\right)\right)\\ -\left(u_h^n\left(\frac{u^{n+1}-u^{n-1}}{2\tau }-\frac{u_h^{n+1}-u_h^{n-1}}{2\tau }\right),\tau {\Delta }_h\left({\hat{e}}_q^{n+1}-{\hat{e}}_q^{n-1}\right)\right)\\ ={\displaystyle \underset{i=71}{\overset{73}{\sum }}{J}_i},\end{array}$ (47)

其中

$J_{71}\le C{\tau }^4+\frac{1}{16}{\Vert \nabla {\hat{e}}_q^{n+1}-\nabla {\hat{e}}_q^{n-1}\Vert }^2.$ (48)

$J_{72}\le C{h}^{2r}+C{\Vert \nabla {\hat{e}}_u^n\Vert }^2+\frac{1}{16}{\Vert \nabla {\hat{e}}_q^{n+1}-\nabla {\hat{e}}_q^{n-1}\Vert }^2.$ (49)

$\begin{array}{l}{J}_{73}\le \frac{1}{2}\Vert \nabla {\hat{e}}_u^{n+1}+\nabla {\tilde{e}}_u^{n+1}+\nabla {\hat{e}}_u^{n-1}+\nabla {\tilde{e}}_u^{n-1}\Vert \Vert \nabla {\hat{e}}_q^{n+1}-\nabla {\hat{e}}_q^{n-1}\Vert \\ \le C{h}^{2r}+\frac{1}{4}{\Vert \nabla {\hat{e}}_u^{n+1}\Vert }^2+C{\Vert \nabla {\hat{e}}_u^{n-1}\Vert }^2+\frac{1}{4}{\Vert \nabla {\hat{e}}_q^{n+1}-\nabla {\hat{e}}_q^{n-1}\Vert }^2.\end{array}$ (50)

相加可得

$M_7\le C{\tau }^4+C{h}^{2r}+\frac{1}{4}{\Vert \nabla {\hat{e}}_u^{n+1}\Vert }^2+C{\Vert \nabla {\hat{e}}_u^n\Vert }^2+C{\Vert \nabla {\hat{e}}_u^{n-1}\Vert }^2+\frac{3}{8}{\Vert \nabla {\hat{e}}_q^{n+1}-\nabla {\hat{e}}_q^{n-1}\Vert }^2.$ (51)

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

$\begin{array}{l}\tau {\Vert \nabla {\hat{e}}_{\omega }^n\Vert }^2+\frac{{\epsilon }^2}{4\tau }\left({\Vert \nabla {\hat{e}}_u^{n+1}\Vert }^2-{\Vert \nabla {\hat{e}}_u^{n-1}\Vert }^2\right)+\frac{1}{2}{\Vert \nabla {\hat{e}}_q^{n+1}-\nabla {\hat{e}}_q^{n-1}\Vert }^2+\beta {\Vert {\delta }_{\tau }{\hat{e}}_u^{n+1}\Vert }^2\\ \le C{\tau }^4+C{h}^{2r}+\frac{3}{8}{\Vert \nabla {\hat{e}}_u^{n+1}\Vert }^2+\left(\frac{1}{16}+\frac{1}{16{\tau }^2}\right){\Vert {\delta }_{\tau }{\hat{e}}_u^{n+1}\Vert }^2+\frac{9}{2}{\Vert {\delta }_{\tau }{\hat{e}}_u^{n+1}\Vert }_{-1,h}^2\\ +\frac{7}{16}{\Vert \nabla {\hat{e}}_q^{n+1}-\nabla {\hat{e}}_q^{n-1}\Vert }^2+C{\Vert \nabla {\hat{e}}_u^n\Vert }^2+C{\Vert \nabla {\hat{e}}_u^{n-1}\Vert }^2+C{\Vert \nabla {\hat{e}}_{\omega }^n\Vert }^2.\end{array}$ (52)

接下来估计

$\begin{array}{l}{\Vert {\delta }_{\tau }{\hat{e}}_u^{n+1}\Vert }_{-1,h}^2=\left(u_t^n-\frac{u_h^{n+1}-u_h^{n-1}}{2\tau },T_h\left({\delta }_{\tau }{\hat{e}}_u^{n+1}\right)\right)-\left(\nabla {\hat{e}}_{\omega }^n+\nabla {\tilde{e}}_{\omega }^n,\nabla T_h\left({\delta }_{\tau }{\hat{e}}_u^{n+1}\right)\right)\\ \le \Vert u_t^n-\frac{u^{n+1}-u^{n-1}}{2\tau }\Vert \Vert T_h\left({\delta }_{\tau }{\hat{e}}_u^{n+1}\right)\Vert +\Vert {\delta }_{\tau }{\tilde{e}}_u^{n+1}\Vert \Vert T_h\left({\delta }_{\tau }{\hat{e}}_u^{n+1}\right)\Vert \\ +\Vert \nabla {\hat{e}}_{\omega }^n+\nabla {\tilde{e}}_{\omega }^n\Vert \Vert \nabla T_h\left({\delta }_{\tau }{\hat{e}}_u^{n+1}\right)\Vert \\ \le C{\tau }^4+C{h}^{2r}+C{\Vert \nabla {\hat{e}}_{\omega }^n\Vert }^2+\frac{3}{4}{\Vert {\delta }_{\tau }{\hat{e}}_u^{n+1}\Vert }_{-1,h}^2.\end{array}$ (53)

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

$\tau {\displaystyle \underset{i=1}{\overset{n}{\sum }}{\Vert \nabla {\hat{e}}_{\omega }^i\Vert }^2}+\frac{{\epsilon }^2}{4\tau }{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\Vert \nabla {\hat{e}}_u^{i+1}\Vert }^2}+\frac{1}{16}{\displaystyle \underset{i=1}{\overset{n}{\sum }}{\Vert \nabla {\hat{e}}_q^{i+1}-\nabla {\hat{e}}_q^{i-1}\Vert }^2}+\beta {\displaystyle \underset{i=1}{\overset{n}{\sum }}{\Vert {\delta }_{\tau }{\hat{e}}_u^{i+1}\Vert }^2}\le C{\tau }^4+C{h}^{2r}.$ (54)

5. 数值模拟

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

5.1. 收敛性

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

$u_0=0.25\sin \left(2x\right)\cos \left(2y\right)+\mathrm{0.48.}$

5.1.1. 时间收敛阶

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

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

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

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

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

5.1.2. 空间收敛阶

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

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

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

5.2. 无条件能量稳定

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

$E\left(u_h^{n+1}\right)={\displaystyle \underset{\Omega }{\int }\left(\frac{{\epsilon }^2}{2}{\left|\nabla u_h^{n+1}\right|}^2+\frac{1}{4}{\left|q_h^{n+1}\right|}^2\right)\text{d}x}$

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

$E_h^{n+1,n}=\frac{{\epsilon }^2}{2}\left({\Vert \nabla u_h^{n+1}\Vert }^2+{\Vert u_h^n\Vert }^2\right)+\frac{1}{4}\left({\Vert q_h^{n+1}\Vert }^2+{\Vert q_h^n\Vert }^2\right).$

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

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

图1. 离散的能量演化

5.3. 相分离

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

$u_0=2rand\left(\right)-1,$

其中 $rand\left(\right)$ 是 $[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方程的能量稳定的二阶数值格式，然后证明数值格式的稳定性和误差分析。最后通过不同的数值实验验证了前面理论的有效性，时间收敛阶与空间收敛阶表明 ${\Vert {\hat{e}}_u\Vert }_{H^1}$ 的收敛阶接近于2，粗化图显示了u随着时间增长的演化，能量耗散的数值结果表明能量随着时间的增加而减少，最后趋于稳定。

致谢

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

基金项目

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

NOTES

^*通讯作者。