黏性Cahn-Hilliard方程的二阶BDF数值格式

郭 媛, 王旦霞, 张建文

(太原理工大学 数学学院, 太原 030024)

0 引 言

经典Cahn-Hilliard方程用于描述非均匀体系中的相分离和粗化现象[1-3]. 黏性Cahn-Hilliard方程[4]是对经典Cahn-Hilliard方程的推广. 目前, 关于黏性Cahn-Hilliard方程数值解法的研究已得到广泛关注. 文献[5]基于标量辅助变量方法构造了黏性Cahn-Hilliard方程的一阶和二阶数值格式; 文献[6]给出了时间双层网格的有限元数值方法; 文献[7]对带有非恒定梯度能量系数的黏性Cahn-Hilliard方程建立了有限元数值格式; 文献[8]使用凸分裂方法提出了有限差分格式, 并证明了所提格式是无条件能量稳定的.

求解黏性Cahn-Hilliard方程的关键是如何在保持能量稳定性的条件下, 对非线性项进行线性离散. 本文采用文献[9]的Lagrange乘子方法, 在黏性Cahn-Hilliard方程中构造线性数值格式. 引入Lagrange乘子黏性Cahn-Hilliard方程如下:

(1)

其边界条件和初值条件分别为

模型(1)的能量函数定义[10]为

满足能量耗散定律

并且是质量守恒的, 即(u(·,t),1)=(u0,1).

本文首先给出模型(1)的半离散格式和全离散格式; 其次给出能量稳定性分析及所提格式的二阶收敛估计; 最后给出一些数值算例证明所提格式的精确性和有效性.

1 离散格式

1.1 半离散格式

模型(1)的混合弱形式为

把时间区间[0,T]做一致划分0=t0

考虑模型(1)的半离散格式, 即给定un-1,un, 求un+1满足

其中

1.2 全离散格式

Sh={vh∈C(Ω)||vh|K∈Pk(x,y),K∈T}⊂H1(Ω),

这里Pk(x,y)是x,y的次数不超过k∈+的多项式集合.

2 稳定性分析

对方程组(8)-(10)求和得

根据2a·(3a-4b+c)=a2-b2+(2a-b)2-(2b-c)2+(a-2b+c)2, 得

证毕.

证毕.

(12)

证明: 将式(11)从1～n求和即可得式(12).

3 误差分析

为简单, 引入下列符号:

对于(u,r), 做如下正则性假设:

u∈w3,∞(0,T;L2(Ω))∩w1,∞(0,T;Hq+1(Ω)),r∈w3,∞(0,T;L2(Ω))∩w1,∞(0,T;Hq+1(Ω)).

定义1[11]Ritz算子Rh:H1(Ω)→Sh满足

((u-Rhu),u)=0, ∀v∈Sh, (Rhu-u,1)=0,

并且Ritz投影算子满足以下估计:

‖u-Rhu‖+h‖u-Rhu‖H1(Ω)≤Chq+1‖u‖Hq+1(Ω).

引理1[12]假设u是方程(1)的解, 则有如下估计:

成立, 其中CT,ε表示常数C与T和ε有关.

证明: 当t=n+1时, 方程组(2)-(4)减去方程组(5)-(7), 得

当n=0时, 由2a·(a-b)=a2-b2+(a-b)2得,

其中

下面依次估计Mi.根据Young不等式[13]、 Cauchy-Schwarz不等式和引理1, 得

其中

当n=0时, 有

根据Cauchy-Schwarz不等式、 引理1及Young不等式, 得

当n=0时, 对M8估计如下:

把式(17)～(27)代入式(16), 并将两边同乘4τ: 当n≥1时, 有

其中

当n=0时, 有

证毕.

4 数值分析

下面通过数值算例[14-15]对理论误差估计和能量稳定性进行验证, 其中u,w,r取P2元[16]有限元空间.

4.1 空间收敛阶

表1 当ε=0.1时的空间收敛阶Table 1 Spatial convergence order of when ε=0.1

u0=0.5+0.17cos(πx)cos(2πy)+0.2cos(3πx)cos(πy).

(30)

4.2 时间收敛阶

表2 当β=0.04时的时间收敛阶Table 2 Time convergence order of when β=0.04

4.3 能量耗散

图1 能量随时间的演化曲线Fig.1 Evolution curves of energy with time

4.4 相分离

(A) T=0.001 s; (B) T=0.1 s; (C) T=0.5 s; (D) T=1.5 s; (E) T=3 s; (F) T=5 s.图3 当β=0.1时的相分离过程Fig.3 Phase separation process when β=0.1