粘性Cahn-Hilliard方程Crank-Nicolson/Adams-Bashforth格式的混合有限元方法

Crank-Nicolson/Adams-Bashforth Scheme of Mixed Finite Element Method for the Viscous Cahn-Hilliard Equation

摘要: 本文基于线性的Crank-Nicolson/Adams-Bashforth格式，建立了有效稳定的数值算法求解粘性Cahn-Hilliard方程。在该算法中，非线性体积力被显示处理，导致求解具有常系数的线性系统，从而提高算法效率。此外，该算法的稳定性以及误差估计被详细证明。最终，通过数值实验证明了该格式的稳定性及收敛阶。

Abstract: In this paper, an efficient stabilized algorithm based on the linear Crank-Nicolson/Adams-Bashforth scheme is proposed to solve the viscous Cahn-Hilliard equation. In this algorithm, the nonlinear bulk force is treated explicitly with linear stabilization term. This treatment leads to solve linear systems with constant coefficients that can improve algorithm efficiency. Further, the stability analysis and priori error estimates on proposed method are provided in detail. Finally, a series of numerical experiments are implemented to illustrate the theoretical analysis.

1. 引言

由Cahn和Hilliard在1958年提出的Cahn-Hilliard方程，常用于描述二元合金在某种不稳定状态时相的分离和粗化现象 [1]。粘性Cahn-Hilliard方程是由Novick Cohen等人 [2] [3] 在研究带粘性的二物质相互扩时提出，主要描述在冷却两种溶液如合金、玻璃及聚合物的混合体时出现的粘性相变。

关于Cahn-Hilliard方程的建模，算法，推导和计算方法已经被广泛研究。值得注意的是，设计高效、准确的数值方案是必须的。其中，显式格式通常会导致严重的时间步长限制，而且通常不满足能量衰减特性；隐显格式的优点是是求解常系数代数方程，能使得求解更容易实现 [4] [5] [6]。然而，半隐式格式通常具有较大的截断误差，因此需要较小的时间步骤来保证精度和能量稳定性。另一方面，满足能量稳定性和较小截断误差的隐式格式，主要缺点是需要在每个时间步长上求解非线性系统。

针对粘性Cahn-Hilliard方程的四阶导数和非线性特点，基于文 [7] [8] 提出的稳定的二阶隐显格式，我们研究了混合有限元的Crank-Nicolson/Adams-Bashforth (CN/AB)格式。本文要研究的粘性Cahn-Hilliard方程如下：

$\{\begin{array}{l}{u}_t=\Delta w,\left(x,t\right)\in \Omega \times \left(0,T\right)\\ w=-{\epsilon }^2\Delta u+f\left(u\right)+\beta u_t,\left(x,t\right)\in \Omega \times \left(0,T\right)\\ {\partial }_{n}u={\partial }_{n}w=0,\left(x,t\right)\in \partial \Omega \times \left(0,T\right)\\ u\left(x,0\right)=u_0\left(x\right),x\in \Omega ,\end{array}$ (1.1)

其中 $u_t=\frac{\partial u}{\partial t},\Omega \subset R^d\left(d\le 3\right)$ 是有界区域且具有光滑边界 $\partial \Omega$，n表示边界上 $\partial \Omega$ 的单位外法向量，u是两

相物质之一的浓度，参数 $\epsilon$ 是一个与界面厚度有关的参数， $\beta \ge 0$ 是粘性系数。其中函数 $f\left(u\right)=F^{\prime }\left(u\right)=u^3-u$ 满足如下条件 [9]：即存在一个常数L，有

$\underset{u\in \Re }{\max }\left|f^{\prime }\left(u\right)\right|\le L$ (1.2)

其中， $F\left(u\right)=\frac{1}{4}{\left(u^2-1\right)}^2$ 是双势阱函数。该方程的自由能被定义为

$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}$ (1.3)

2. 预备知识

首先我们引入一些记号。 $H^m$ 表示标准的Sobolev空间并且范数定义为 ${\Vert u\Vert }_m={\left({{\displaystyle \underset{\alpha \le m}{\sum }\Vert D^{\alpha }u\Vert }}^2\right)}^{\frac{1}{2}}$。 $\left(\cdot ,\cdot \right)$ 和 $\Vert \cdot \Vert$ 分别表示 $L^2=H^0$ 上的内积和范数。定义 $L^2\left(0,T;H^m\left(\Omega \right)\right)=\left\{u\left(x,t\right)\in H^m;{\displaystyle {\int }_0^T{\Vert u\Vert }_m^2<\infty }\right\}$。

式(1.1)的混合变分形式为：

$\{\begin{array}{l}\left(u_t,v\right)+\left(\nabla w,\nabla v\right)=0,\forall v\in H^1\left(\Omega \right)\\ \left(w,\varphi \right)={\epsilon }^2\left(\nabla u,\nabla \varphi \right)+\left(f\left(u\right),\varphi \right)+\beta \left(u_t,\varphi \right),\forall \varphi \in H^1\left(\Omega \right)\end{array}$ (2.1)

令 $T_h$ 为 $\Omega$ 的拟一致三角剖分，h为空间网格大小， $S_h$ 为定义在 $T_h$ 上的分段多项式空间：

$S_h=\left\{v\in C^0\left(\Omega \right)|{v|}_K\in P_r\left(K\right),\forall K\in T_h\right\}\subset H^1\left(\Omega \right)$，其中 $P_r\left(k\right)$ 是次数不超过r的多项式空间。

我们引入Ritz投影算子 $R_h:H^1\left(\Omega \right)\to S_h$，满足：

$\begin{array}{l}\left(\nabla \left(R_h\varphi -\varphi \right),\nabla \chi \right)=0,\forall \chi \in S_h,\\ \left(R_h\varphi -\varphi ,1\right)=0,\end{array}$ (2.2)

并且算子满足如下不等式 [10]：

${\Vert \varphi -R_h\varphi \Vert }_p+h{\Vert \varphi -R_h\varphi \Vert }_{1,p}\le C{h}^{q+1}{\Vert \varphi \Vert }_{q+1,p},\forall 1\le p\le \infty .$ (2.3)

定义截断误差 $R^{n+\frac{1}{2}}=R_1^{n+\frac{1}{2}}+R_2^{n+\frac{1}{2}}$，其中

$R_1^{n+\frac{1}{2}}:=\left(\frac{u\left(t_{n+1}\right)-u\left(t_n\right)}{\tau }\right)-u_t\left(t_{n+\frac{1}{2}}\right),$ (2.4)

$R_2^{n+\frac{1}{2}}:=-\Delta \left(\frac{u\left(t_{n+1}\right)+u\left(t_n\right)}{2}-u\left(t_{n+\frac{1}{2}}\right)\right),$ (2.5)

并且满足如下不等式 [11]：

${\Vert R_1^{n+\frac{1}{2}}\Vert }_s^2\le {\tau }^3{\displaystyle {\int }_{t_n}^{t_{n+1}}{\Vert u_{ttt}\left(t\right)\Vert }_s^2\text{d}t,s=-1,0}$ (2.6)

${\Vert R_2^{n+\frac{1}{2}}\Vert }_s^2\le {\tau }^3{\displaystyle {\int }_{t_n}^{t_{n+1}}{\Vert u_{ttt}\left(t\right)\Vert }_{s+2}^2\text{d}t,s=-1,0}$ (2.7)

3. 二阶Crank-Nicolson/Adams-Bashforth格式

将时间区间 $\left[0,T\right]$ 剖为： $0=t_0<t_1<\cdots <t_n=T$， $t_{n+1}-t_n=\tau =T/N$，N是一个正整数。全离散格式逼近问题(2.1)可以表述为下述方程，给定 $u_h^{n-1},u_h^n\in S_h$，求 $u_h^{n+1},w_h^{n+\frac{1}{2}}\in S_h$，有

$\{\begin{array}{l}\left({\delta }_t{u}_h^{n+1},v_h\right)+\left(\nabla w_h^{n+\frac{1}{2}},\nabla v_h\right)=0\text{,}\forall v_h\in S_h\text{(3}\text{.1)}\\ \left(w_h^{n+\frac{1}{2}},{\varphi }_h\right)={\epsilon }^2\left(\nabla u_h^{n+\frac{1}{2}},\nabla {\varphi }_h\right)+S\left(u_h^{n+1}-2u_h^n+u_h^{n-1},{\varphi }_h\right)\\ +\left(\frac{3}{2}f\left(u_h^n\right)-\frac{1}{2}f\left(u_h^{n-1}\right),{\varphi }_h\right)+\beta \left({\delta }_t{u}_h^{n+1},{\varphi }_h\right),\forall {\varphi }_h\in S_h\text{(3 .2)}\end{array}$

其中 ${\delta }_t{u}_h^{n+1}:=\frac{u_h^{n+1}-u_h^n}{\tau }$， $u_h^{n+\frac{1}{2}}:=\frac{u_h^{n+1}+u_h^n}{2}$， $\tau$ 是时间步长。

定理3.1 设 $\left(u_h^{n+1},w_h^{n+\frac{1}{2}}\right)$ 是格式(3.1)-(3.2)的解，在 $\tau \le \beta /L\left(\beta ,L\ge 0\right)$ 条件下，该格式是稳定的并且满足下面的能量法则：

$E\left(u_h^{n+1}\right)+\left(\frac{\beta }{\tau }-\frac{3L}{4}\right){\Vert u_h^{n+1}-u_h^n\Vert }^2\le E\left(u_h^n\right)+\frac{L}{4}{\Vert u_h^n-u_h^{n-1}\Vert }^2.$

证明 在(3.1)中取 $v_h=w_h^{n+\frac{1}{2}}$，有

$\left(u_h^{n+1}-u_h^n,w_h^{n+\frac{1}{2}}\right)+\tau {\Vert \nabla w_h^{n+\frac{1}{2}}\Vert }^2=0$ (3.3)

在(3.2)中取 ${\varphi }_h=u_h^{n+1}-u_h^n$，有

$\begin{array}{l}\left(w_h^{n+\frac{1}{2}},u_h^{n+1}-u_h^n\right)\\ ={\epsilon }^2\left(\nabla u_h^{n+\frac{1}{2}},\nabla \left(u_h^{n+1}-u_h^n\right)\right)+S\left(u_h^{n+1}-2u_h^n+u_h^{n-1},u_h^{n+1}-u_h^n\right)\\ +\left(\frac{3}{2}f\left(u_h^n\right)-\frac{1}{2}f\left(u_h^{n-1}\right),u_h^{n+1}-u_h^n\right)+\beta \left({\delta }_t{u}_h^{n+1},u_h^{n+1}-u_h^n\right)\\ :={\prod }_1+{\prod }_2+{\prod }_3+{\prod }_4\end{array}$

接下来我们对 ${\prod }_1~{\prod }_4$ 进行估计，

${\prod }_1=\frac{{\epsilon }^2}{2}\left(\nabla \left(u_h^{n+1}+u_h^n\right),\nabla \left(u_h^{n+1}-u_h^n\right)\right)=\frac{{\epsilon }^2}{2}{\Vert \nabla u_h^{n+1}\Vert }^2+\frac{{\epsilon }^2}{2}{\Vert \nabla u_h^n\Vert }^2$

利用公式 $2a\left(a-b\right)=a^2-b^2+{\left(a-b\right)}^2$，则

$\begin{array}{c}{\prod }_2=\left(\left(u_h^{n+1}-u_h^n\right)-\left(u_h^n-u_h^{n-1}\right),u_h^{n+1}-u_h^n\right)\\ =\frac{1}{2}\left({\Vert u_h^{n+1}-u_h^n\Vert }^2-{\Vert u_h^n-u_h^{n-1}\Vert }^2+{\Vert u_h^{n+1}-2u_h^{n+1}+u_h^{n-1}\Vert }^2\right)\end{array}$

利用泰勒公式，则

$\begin{array}{c}{\prod }_{\text{3}}=\frac{1}{2}\left(f\left(u_h^n\right)-f\left(u_h^{n-1}\right),u_h^{n+1}-u_h^n\right)+\left(f\left(u_h^n\right),u_h^{n+1}-u_h^n\right)\\ +\left(\frac{1}{2}{f}^{\prime }\left(\xi \right){\left(u_h^{n+1}-u_h^n\right)}^2,1\right)-\left(\frac{1}{2}{f}^{\prime }\left(\xi \right){\left(u_h^{n+1}-u_h^n\right)}^2,1\right)\\ =\left(\frac{1}{2}{f}^{\prime }\left(\xi \right)\left(u_h^n-u_h^{n-1}\right),u_h^{n+1}-u_h^n\right)+\left(F\left(u_h^{n+1}\right)-F\left(u_h^n\right),1\right)\\ -\left(\frac{1}{2}{f}^{\prime }\left(\xi \right){\left(u_h^{n+1}-u_h^n\right)}^2,1\right)\end{array}$

$\begin{array}{l}{\prod }_4=\frac{\beta }{\tau }\left(u_h^{n+1}-u_h^n,u_h^{n+1}-u_h^n\right)\\ =\frac{\beta }{\tau }{\Vert u_h^{n+1}\Vert }^2-\frac{\beta }{\tau }{\Vert u_h^n\Vert }^2\end{array}$

结合(3.3)和以及对的估计，有

$\begin{array}{l}\frac{{\epsilon }^2}{2}\left({\Vert \nabla u_h^{n+1}\Vert }^2-{\Vert \nabla u_h^n\Vert }^2\right)+\left(F\left(u_h^{n+1}\right)-F\left(u_h^n\right),1\right)\\ +\frac{S}{2}\left({\Vert u_h^{n+1}-u_h^n\Vert }^2-{\Vert u_h^n-u_h^{n-1}\Vert }^2+{\Vert u_h^{n+1}-2u_h^n+u_h^{n-1}\Vert }^2\right)\\ +\left(F\left(u_h^{n+1}\right)-F\left(u_h^n\right),1\right)-\left(w_h^{n+\frac{1}{2}},u_h^{n+1}-u_h^n\right)+\frac{\beta }{\tau }{\Vert u_h^{n+1}-u_h^n\Vert }^2\\ =-\left(\frac{1}{2}{f}^{\prime }\left(\xi \right)\left(u_h^n-u_h^{n-1}\right),u_h^{n+1}-u_h^n\right)+\left(\frac{1}{2}{f}^{\prime }\left(\xi \right){\left(u_h^{n+1}-u_h^n\right)}^2,1\right)\end{array}$ (3.4)

利用Cauchy-Schwarz不等式和(1.2)，即等式(3.4)的右边项有：

$\begin{array}{l}\left(\frac{f^{\prime }\left(\xi \right)}{2}{\left(u_h^{n+1}-u_h^n\right)}^2,1\right)-\left(\frac{f^{\prime }\left(\xi \right)}{2}\left(u_h^n-u_h^{n-1}\right),\left(u_h^{n+1}-u_h^n\right)\right)\\ \le \frac{L}{2}{\Vert u_h^{n+1}-u_h^n\Vert }^2+\frac{L}{2}\Vert u_h^n-u_h^{n-1}\Vert \Vert u_h^{n+1}-u_h^n\Vert \\ \le \frac{3L}{4}{\Vert u_h^{n+1}-u_h^n\Vert }^2+\frac{L}{4}{\Vert u_h^n-u_h^{n-1}\Vert }^2\end{array}$ (3.5)

根据能量函数(1.3)的定义，结合(3.4)和(3.5)，该定理得证。

定理3.2 设 $\left(u\left(t_n\right),w\left(t_n\right)\right)$ 和 $\left(u_h^n,w_h^{n+\frac{1}{2}}\right)$ 分别是方程(2.1)和方程(3.1)，(3.2)的解。且假设

$u,w\in C\left(0,T;H^1\left(\Omega \right)\right)$， $u_t\in L^2\left(0,T;H^1\left(\Omega \right)\cap L^2\left(\Omega \right)\right)$， $u_{tt}\in L^2\left(0,T;H^2\left(\Omega \right)\right)$， $u_{ttt}\in L^2\left(0,T;H^{-1}\left(\Omega \right)\cap L^2\left(\Omega \right)\right)$，则满足下述估计式

$\Vert u\left(t_n\right)-u_h^n\Vert \le C\left({\tau }^2+h^2\right).$

证明：首先给一些符号

$\begin{array}{l}{\bar{e}}^n=u\left(t_n\right)-R_{h}u\left(t_n\right),e^n=u_h^n-R_{h}u\left(t_n\right)\\ {\bar{\eta }}^n=w\left(t_n\right)-R_{h}w\left(t_n\right),{\eta }^n=w_h^n-R_{h}w\left(t_n\right)\\ {\bar{\eta }}^{n+\frac{1}{2}}=w\left(t_{n+\frac{1}{2}}\right)-R_{h}w\left(t_{n+\frac{1}{2}}\right),{\eta }^n=w_h^{n+\frac{1}{2}}-R_{h}w\left(t_{n+\frac{1}{2}}\right)\end{array}$

然后全离散格式利用上面的关系可以写成下面的等价式：

$\left(\frac{e^{n+1}-e^n}{\tau },v_h\right)+\left(\nabla {\eta }^{n+\frac{1}{2}},\nabla v_h\right)=\left(\frac{{\bar{e}}^{n+1}-{\bar{e}}^n}{\tau },v_h\right)-\left(R_1^{n+\frac{1}{2}},v_h\right)+\left(\nabla {\bar{\eta }}^{n+\frac{1}{2}},\nabla v_h\right)$ (3.6)

$\begin{array}{l}{\epsilon }^2\left(\nabla e^{n+\frac{1}{2}},\nabla {\varphi }_h\right)+{\epsilon }^2\left(R_2^{n+\frac{1}{2}},{\varphi }_h\right)+S\left(e^{n+1}-2e^n+e^{n-1},{\varphi }_h\right)\\ -S\left(u\left(t_{n+1}\right)-2u\left(t_n\right)+u\left(t_{n-1}\right),{\varphi }_h\right)+S\left(\left(I-R_h\right)\left(u\left(t_{n+1}\right)-2u\left(t_n\right)+u\left(t_{n-1}\right)\right),{\varphi }_h\right)\\ +\beta \left(e^{n+1}-e^n-\left({\bar{e}}^{n+1}-{\bar{e}}^n\right),{\varphi }_h\right)+\beta \left(R_1^{n+\frac{1}{2}},{\varphi }_h\right)\\ =\left(f\left(u\left(t_{n+\frac{1}{2}}\right)\right)-\frac{3}{2}f\left(u_h^n\right)+\frac{1}{2}f\left(u_h^{n-1}\right),{\varphi }_h\right)+\left({\eta }^{n+\frac{1}{2}}-{\bar{\eta }}^{n+\frac{1}{2}},{\varphi }_h\right)\end{array}$ (3.7)

在(3.6)中分别取 $v_h=2\tau e^{n+\frac{1}{2}},\frac{2\beta \tau }{{\epsilon }^2}{\eta }^{n+\frac{1}{2}}$ 以及在(3.7)中取 ${\varphi }_h=-\frac{2\tau {\eta }^{n+\frac{1}{2}}}{{\epsilon }^2}$ 得到三个等式，对其求和有，

$\begin{array}{l}{\Vert e^{n+1}\Vert }^2-{\Vert e^n\Vert }^2+\frac{2\tau }{{\epsilon }^2}{\Vert {\eta }^{n+\frac{1}{2}}\Vert }^2+\frac{2\beta \tau }{{\epsilon }^2}{\Vert \nabla {\eta }^{n+\frac{1}{2}}\Vert }^2\\ =\left(\frac{{\bar{e}}^{n+1}-{\bar{e}}^n}{\tau },2\tau e^{n+\frac{1}{2}}\right)+\left({\bar{\eta }}^{n+\frac{1}{2}},\frac{2\tau }{{\epsilon }^2}{\eta }^{n+\frac{1}{2}}\right)+S\left(e^{n+1}-2e^n+e^{n-1},\frac{2\tau }{{\epsilon }^2}{\eta }^{n+\frac{1}{2}}\right)\\ -\left(f\left(u\left(t_{n+\frac{1}{2}}\right)\right)-\frac{3}{2}f\left(u_h^n\right)+\frac{1}{2}f\left(u_h^{n-1}\right),\frac{2\tau }{{\epsilon }^2}{\eta }^{n+\frac{1}{2}}\right)-\left(R_1^{n+\frac{1}{2}},\tau e^{n+\frac{1}{2}}\right)+\left(R_2^{n+\frac{1}{2}},2\tau {\eta }^{n+\frac{1}{2}}\right)\\ -S\left(u\left(t_{n+1}\right)-2u\left(t_n\right)+u\left(t_{n-1}\right),\frac{2\tau }{{\epsilon }^2}{\eta }^{n+\frac{1}{2}}\right)+S\left(\left(I-R_h\right)\left(u\left(t_{n+1}\right)-2u\left(t_n\right)+u\left(t_{n-1}\right),\frac{2\tau }{{\epsilon }^2}{\eta }^{n+\frac{1}{2}}\right)\right)\\ :={\displaystyle \underset{i}{\overset{8}{\sum }}{\prod }_i}\end{array}$ (3.8)

接下来，我们对 ${\prod }_1~{\prod }_9$ 进行估计，其中 ${\prod }_1~{\prod }_4$ 利用Cauchy-Schwarz不等式，Young不等式有，

${\prod }_1\le \left|\left(\frac{{\bar{e}}^{n+1}-{\bar{e}}^n}{\tau },2\tau e^{n+\frac{1}{2}}\right)\right|\le 2\tau \Vert {\delta }_t{\bar{e}}^{n+\frac{1}{2}}\Vert \Vert e^{n+\frac{1}{2}}\Vert \le 2\tau {\Vert {\delta }_t{\bar{e}}^{n+\frac{1}{2}}\Vert }^2+\frac{\tau }{2}{\Vert e^{n+\frac{1}{2}}\Vert }^2$

${\prod }_2\le \left|\left({\bar{\eta }}^{n+\frac{1}{2}},\frac{2\tau }{{\epsilon }^2}{\eta }^{n+\frac{1}{2}}\right)\right|\le \frac{2\tau }{{\epsilon }^2}\Vert {\bar{\eta }}^{n+\frac{1}{2}}\Vert \Vert {\eta }^{n+\frac{1}{2}}\Vert \le \frac{8\tau }{{\epsilon }^2}{\Vert {\bar{\eta }}^{n+\frac{1}{2}}\Vert }^2+\frac{\tau }{8{\epsilon }^2}{\Vert {\eta }^{n+\frac{1}{2}}\Vert }^2$

$\begin{array}{c}{\prod }_3\le \left|S\left(e^{n+1}-2e^n+e^{n-1},\frac{2\tau }{{\epsilon }^2}{\eta }^{n+\frac{1}{2}}\right)\right|\\ \le \left|S\left(e^{n+1}-e^n,\frac{2\tau }{{\epsilon }^2}{\eta }^{n+\frac{1}{2}}\right)\right|+\left|S\left(e^{n-1}-e^n,\frac{2\tau }{{\epsilon }^2}{\eta }^{n+\frac{1}{2}}\right)\right|\\ \le \frac{2S\tau }{{\epsilon }^2}\Vert e^{n+1}-e^n\Vert \Vert {\eta }^{n+\frac{1}{2}}\Vert +\frac{2S\tau }{{\epsilon }^2}\Vert e^{n-1}-e^n\Vert \Vert {\eta }^{n+\frac{1}{2}}\Vert \\ \le \frac{8S^2\tau }{{\epsilon }^2}{\Vert e^{n+1}-e^n\Vert }^2+\frac{8S^2\tau }{{\epsilon }^2}{\Vert e^{n-1}-e^n\Vert }^2+\frac{2\tau }{8{\epsilon }^2}{\Vert {\eta }^{n+\frac{1}{2}}\Vert }^2\end{array}$

$\begin{array}{c}{\prod }_4\le \left|\left(f\left(u\left(t_{n+\frac{1}{2}}\right)\right)-\frac{3}{2}f\left(u_h^n\right)+\frac{1}{2}f\left(u_h^{n-1}\right),\frac{2\tau }{{\epsilon }^2}{\eta }^{n+\frac{1}{2}}\right)\right|\\ \le \frac{8\tau }{{\epsilon }^2}{\Vert f\left(u\left(t_{n+\frac{1}{2}}\right)\right)-\frac{3}{2}f\left(u_h^n\right)+\frac{1}{2}f\left(u_h^{n-1}\right)\Vert }^2+\frac{\tau }{8{\epsilon }^2}{\Vert {\eta }^{n+\frac{1}{2}}\Vert }^2\\ =\frac{8\tau }{{\epsilon }^2}{\prod }_9+\frac{\tau }{8{\epsilon }^2}{\Vert {\eta }^{n+\frac{1}{2}}\Vert }^2\end{array}$

${\prod }_5,{\prod }_6$ 利用截断误差(2.6)和(2.7)，Cauchy-Schwarz不等式，Young不等式有，

$\begin{array}{c}{\prod }_5\le \left|-\left(R_1^{n+\frac{1}{2}},\tau e^{n+\frac{1}{2}}\right)\right|\le \tau \Vert R_1^{n+\frac{1}{2}}\Vert \Vert e^{n+\frac{1}{2}}\Vert \le \frac{\tau }{2}{\Vert R_1^{n+\frac{1}{2}}\Vert }^2+\frac{\tau }{2}{\Vert e^{n+\frac{1}{2}}\Vert }^2\\ \le \frac{{\tau }^4}{2}{\displaystyle {\int }_{t_n}^{t_{n+1}}{\Vert u_{ttt}\left(s\right)\Vert }^2\text{d}s}+\frac{\tau }{2}{\Vert e^{n+\frac{1}{2}}\Vert }^2\end{array}$

$\begin{array}{c}{\prod }_6\le \left|\left(R_2^{n+\frac{1}{2}},2\tau {\eta }^{n+\frac{1}{2}}\right)\right|\le 2\tau \Vert R_2^{n+\frac{1}{2}}\Vert \Vert {\eta }^{n+\frac{1}{2}}\Vert \le 8{\epsilon }^2\tau {\Vert R_2^{n+\frac{1}{2}}\Vert }^2+\frac{\tau }{8{\epsilon }^2}{\Vert {\eta }^{n+\frac{1}{2}}\Vert }^2\\ \le 8{\epsilon }^2{\tau }^4{\displaystyle {\int }_{t_n}^{t_{n+1}}{\Vert u_{tt}\left(s\right)\Vert }_2^2\text{d}s}+\frac{\tau }{8{\epsilon }^2}{\Vert {\eta }^{n+\frac{1}{2}}\Vert }^2\end{array}$

对于 ${\prod }_7,{\prod }_8,{\prod }_9$ 的估计，利用泰勒展示，Young不等式有，Cauchy-Schwarz不等式，有

$\begin{array}{c}{\prod }_{\text{7}}\le \left|-S\left(u\left(t_{n+1}\right)-2u\left(t_n\right)+u\left(t_{n-1}\right),\frac{2\tau }{{\epsilon }^2}{\eta }^{n+\frac{1}{2}}\right)\right|\\ \le \frac{C{S}^2\tau }{{\epsilon }^2}{\Vert {\displaystyle {\int }_{t_n}^{t_{n+1}}{u}_{tt}\left(s\right)\left(t_{n+1}-s\right)\text{d}s}\Vert }^2+\frac{C{S}^2\tau }{{\epsilon }^2}{\Vert {\displaystyle {\int }_{t_{n-1}}^{t_n}{u}_{tt}\left(s\right)\left(s-t_n\right)\text{d}s}\Vert }^2+\frac{2\tau }{8{\epsilon }^2}{\Vert {\eta }^{n+\frac{1}{2}}\Vert }^2\\ \le \frac{C{S}^2{\tau }^3}{{\epsilon }^2}{\Vert {\displaystyle {\int }_{t_{n-1}}^{t_n}{u}_{tt}\left(s\right)\text{d}s}\Vert }^2+\frac{C{S}^2{\tau }^3}{{\epsilon }^2}{\Vert {\displaystyle {\int }_{t_{n-1}}^{t_n}{u}_{tt}\left(s\right)\text{d}s}\Vert }^2+\frac{2\tau }{8{\epsilon }^2}{\Vert {\eta }^{n+\frac{1}{2}}\Vert }^2\\ \le \frac{C{S}^2{\tau }^4}{{\epsilon }^2}{\displaystyle {\int }_{t_{n-1}}^{t_{n+1}}{\Vert u_{tt}\left(s\right)\Vert }^2\text{d}s}+\frac{2\tau }{8{\epsilon }^2}{\Vert {\eta }^{n+\frac{1}{2}}\Vert }^2\end{array}$

${\prod }_8\le \left|S\left(\left(I-R_h\right)\left(u\left(t_{n+1}\right)-2u\left(t_n\right)+u\left(t_{n-1}\right)\right),\frac{2\tau }{{\epsilon }^2}{\eta }^{n+\frac{1}{2}}\right)\right|\le \frac{C{S}^2{\tau }^4}{{\epsilon }^2}{\displaystyle {\int }_{t_{n-1}}^{t_{n+1}}{\Vert \left(I-R_h\right)u_{tt}\left(s\right)\Vert }^2\text{d}s}+\frac{\tau }{8{\epsilon }^2}{\Vert {\eta }^{n+\frac{1}{2}}\Vert }^2$

$\begin{array}{c}{\prod }_9\le 2({\Vert f\left(u\left(t_{n+\frac{1}{2}}\right)\right)-\frac{3}{2}f\left(u\left(t_n\right)\right)+\frac{1}{2}f\left(u\left(t_{n-1}\right)\right)\Vert }^2\\ +{\Vert \frac{3}{2}f\left(u\left(t_n\right)\right)-\frac{3}{2}f\left(u_h^n\right)-\frac{1}{2}f\left(u\left(t_{n-1}\right)\right)+\frac{1}{2}f\left(u_h^{n-1}\right)\Vert }^2)\\ \le C{\Vert f\left(u\left(t_{n+\frac{1}{2}}\right)\right)-\frac{1}{2}\left(f\left(u\left(t_{n+1}\right)\right)+f\left(u\left(t_n\right)\right)\right)\Vert }^2\\ +C{\Vert \frac{1}{2}\left(f\left(u\left(t_{n+1}\right)\right)+f\left(u\left(t_n\right)\right)\right)-\frac{3}{2}f\left(u\left(t_n\right)\right)+\frac{1}{2}f\left(u\left(t_{n-1}\right)\right)\Vert }^2\\ +C{L}^2\left({\Vert u\left(t_n\right)-u_h^n\Vert }^2+{\Vert u\left(t_{n-1}\right)-u_h^{n-1}\Vert }^2\right)\\ \le CL{\tau }^3{\displaystyle {\int }_{t_{n-1}}^{t_{n+1}}{\Vert u_{tt}\left(s\right)\Vert }^2}\text{d}s+C{L}^2\left({\Vert e^n\Vert }^2+{\Vert {\bar{e}}^n\Vert }^2+{\Vert e^{n-1}\Vert }^2+{\Vert {\bar{e}}^{n-1}\Vert }^2\right)\end{array}$

综上，将 ${\prod }_1\text{~}{\prod }_9$ 不等式带入(3.8)可以得到，

$\begin{array}{l}{\Vert e^{n+1}\Vert }^2-{\Vert e^n\Vert }^2\\ \le 2\tau {\Vert {\delta }_t{\bar{e}}^{n+1}\Vert }^2+\tau {\Vert e^{n+\frac{1}{2}}\Vert }^2+\frac{{\tau }^4}{2}{\displaystyle {\int }_{t_n}^{t_{n+1}}{\Vert u_{ttt}\left(s\right)\Vert }^2}\text{d}s+\frac{8{\tau }^4}{{\epsilon }^2}{\displaystyle {\int }_{t_n}^{t_{n+1}}{\Vert u_{tt}\left(s\right)\Vert }_2^2\text{d}s}+\frac{8S^2\tau }{{\epsilon }^2}{\Vert e^{n+1}-e^n\Vert }^2\\ +\frac{16S^2{\tau }^4}{{\epsilon }^2}{\displaystyle {\int }_{t_{n-1}}^{t_{n+1}}{\Vert u_{tt}\left(s\right)\Vert }^2}\text{d}s+\frac{8S^2\tau }{{\epsilon }^2}{\Vert e^{n-1}-e^n\Vert }^2+\frac{16S^2{\tau }^4}{{\epsilon }^2}{\displaystyle {\int }_{t_{n-1}}^{t_{n+1}}{\Vert \left(I-R_h\right)u_{tt}\left(s\right)\Vert }^2}\text{d}s\\ +\frac{8{\tau }^4}{{\epsilon }^2}\Vert {\bar{\eta }}^{n+\frac{1}{2}}\Vert +\frac{CL{\tau }^4}{{\epsilon }^2}{\displaystyle {\int }_{t_{n-1}}^{t_{n+1}}{\Vert u_{tt}\left(s\right)\Vert }^2}\text{d}s+\frac{C{L}^2\tau }{{\epsilon }^2}\left({\Vert e^n\Vert }^2+{\Vert {\bar{e}}^n\Vert }^2+{\Vert e^{n-1}\Vert }^2+{\Vert {\bar{e}}^{n-1}\Vert }^2\right)\end{array}$ (3.9)

将(3.9)式两端n从0到 $n-1$ 求和，可得

$\begin{array}{l}{\Vert e^n\Vert }^2-{\Vert e^0\Vert }^2\\ \le \tau {\displaystyle \underset{n=0}{\overset{n-1}{\sum }}\Vert e^n\Vert }+C\tau {\displaystyle \underset{n=0}{\overset{n-1}{\sum }}\left({\Vert {\delta }_t{\bar{e}}^{n+1}\Vert }^2+{\Vert {\bar{\eta }}^{n+\frac{1}{2}}\Vert }^2\right)}\\ +C{\tau }^4({\Vert u_{ttt}\left(s\right)\Vert }_{L^2\left(0,T;L^2\right)}^2+{\Vert u_{tt}\left(s\right)\Vert }_{L^2\left(0,T;H^2\right)}^2\\ +{\Vert u_{tt}\left(s\right)\Vert }_{L^2\left(0,T;L^2\right)}^2+{\Vert \left(I-R_h\right)u_{tt}\left(s\right)\Vert }_{L^2\left(0,T;L^2\right)}^2)\\ +\frac{\tau }{2}{\Vert e^n\Vert }^2+C{L}^2\tau {\displaystyle \underset{n=0}{\overset{n-1}{\sum }}{\Vert {\bar{e}}^n\Vert }^2}\end{array}$

限制 $\tau \le \min \left\{1,\frac{L}{\beta }\right\}$，利用离散Gronwall引理和 $e^0$ 定义，可得

$\begin{array}{l}\Vert e^n\Vert \le C\tau {\displaystyle \underset{n=0}{\overset{n-1}{\sum }}\left({\Vert {\delta }_t{\bar{e}}^{n+1}\Vert }^2+{\Vert {\bar{\eta }}^{n+\frac{1}{2}}\Vert }^2+{\Vert {\bar{e}}^n\Vert }^2\right)}\\ +C{\tau }^4({\Vert u_{ttt}\left(s\right)\Vert }_{L^2\left(0,T;L^2\right)}^2+{\Vert u_{tt}\left(s\right)\Vert }_{L^2\left(0,T;H^2\right)}^2\\ +{\Vert u_{tt}\left(s\right)\Vert }_{L^2\left(0,T;L^2\right)}^2+{\Vert \left(I-R_h\right)u_{tt}\left(s\right)\Vert }_{L^2\left(0,T;L^2\right)}^2)\end{array}$

利用三角不等式和Ritze投影算子(2.3)，有

$\Vert u\left(t_n\right)-R_{h}u\left(t_n\right)\Vert \le C{h}^2,$

$\Vert w\left(t_n\right)-R_{h}w\left(t_n\right)\Vert \le C{h}^2,$

$\Vert {\delta }_t\left(u\left(t_n\right)-R_{h}u\left(t_n\right)\right)\Vert \le C{h}^2,$

最终，则可得 ${\Vert u\left(t_n\right)-u_h^n\Vert }^2\le C\left({\tau }^4+h^4\right)$，即 $\Vert u\left(t_n\right)-u_h^n\Vert \le C\left({\tau }^2+h^2\right)$，定理得证。

4. 数值实验

在本部分，采用数值实例验证理论分析的正确性和有效性。在下面的算例中，用 $P_1$ 元构建有限元空间，选择参数 $\epsilon =0.01$， $\beta =0.01$， $T=0.01$， $\tau =0.0001$，计算区域为 $\left[-1,1\right]\times \left[-1,1\right]$，初始条件为：

$u_0\left(x,y\right)=0.1\left(\sin \left(3x\right)\sin \left(2y\right)+\sin \left(5x\right)\sin \left(5y\right)\right).$

表1~3列出了全离散格式在范数 $L^2$ 和范数 $H^1$ 下的误差估计和收敛阶，可以观察到空间收敛阶在 $L^2$ 范数下趋向于2.0，在 $H^1$ 范数下趋向于1.0，与理论结果一致。其次，为了证明格式的稳定性，我们假设最

终时间 $T=0.1$，时间步长 $\tau =0.0001$， $S=0.5,1,2$， $h=\frac{1}{64}$，图1观察到能量是随时间递减的。最终，我们分别刻画了图2~4在时间 $t=0.001,0.05,0.1$ 时数值解的等值线 $S=0.5$， $h=\frac{1}{64}$，观察数值解的变化过程，随着时间的推移，最终等值线会达到一个稳定。

Table 1. Numerial results for the spatial discretization for time step τ = 0.0001 , ( P 1 , P 1 ) mixed finite

表1. 固定时间步长 $\tau =0.0001$, $\left(P_1,P_1\right)$ 混合有限元空间收敛阶

Table 2. Numerial results for the spatial discretization for time step τ = 0.0001 , ( P 1 , P 1 ) mixed finite

表2. 固定时间步长 $\tau =0.0001$, $\left(P_1,P_1\right)$ 混合有限元空间收敛阶

Table 3. Numerial results for the spatial discretization for time step τ = 0.0001 , ( P 1 , P 1 ) mixed finite

表3. 固定时间步长 $\tau =0.0001$, $\left(P_1,P_1\right)$ 混合有限元空间收敛阶

Figure 1. The evolution of the discrete energy for the viscous Cahn-Hilliard equation for $\tau =0.0001$, s = 0.5, 1, 2

图1. $\tau =0.0001$, s = 0.5, 1, 2时，粘性Cahn-Hilliard方程的能量演变

Figure 2. The contour lines at $t=0.001$

图2. 当 $t=0.001$ 时的等值线

Figure 3. The contour lines at $t=0.05$

图3. 当 $t=0.05$ 时的等值线

Figure 4. The contour lines at $t=0.1$

图4. 当 $t=0.1$ 时的等值线

5. 结论

本文研究了粘性的Cahn-Hilliard方程线性CN/AB格式混合有限元方法，达到了二阶收敛阶。最后，通过数值算例对理论分析进行验证，算例结果与理论分析的结果相一致。

基金项目

国家自然科学基金项目(NO. 11172194)；山西省自然科学基金(NO. 201801D121016)。