郭 媛, 王旦霞, 張建文

(太原理工大學(xué) 數(shù)學(xué)學(xué)院, 太原 030024)

0 引 言

經(jīng)典Cahn-Hilliard方程用于描述非均勻體系中的相分離和粗化現(xiàn)象[1-3]. 黏性Cahn-Hilliard方程[4]是對經(jīng)典Cahn-Hilliard方程的推廣. 目前, 關(guān)于黏性Cahn-Hilliard方程數(shù)值解法的研究已得到廣泛關(guān)注. 文獻[5]基于標量輔助變量方法構(gòu)造了黏性Cahn-Hilliard方程的一階和二階數(shù)值格式; 文獻[6]給出了時間雙層網(wǎng)格的有限元數(shù)值方法; 文獻[7]對帶有非恒定梯度能量系數(shù)的黏性Cahn-Hilliard方程建立了有限元數(shù)值格式; 文獻[8]使用凸分裂方法提出了有限差分格式, 并證明了所提格式是無條件能量穩(wěn)定的.

求解黏性Cahn-Hilliard方程的關(guān)鍵是如何在保持能量穩(wěn)定性的條件下, 對非線性項進行線性離散. 本文采用文獻[9]的Lagrange乘子方法, 在黏性Cahn-Hilliard方程中構(gòu)造線性數(shù)值格式. 引入Lagrange乘子黏性Cahn-Hilliard方程如下:

(1)

其邊界條件和初值條件分別為

模型(1)的能量函數(shù)定義[10]為

滿足能量耗散定律

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

本文首先給出模型(1)的半離散格式和全離散格式; 其次給出能量穩(wěn)定性分析及所提格式的二階收斂估計; 最后給出一些數(shù)值算例證明所提格式的精確性和有效性.

1 離散格式

1.1 半離散格式

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

把時間區(qū)間[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的次數(shù)不超過k∈+的多項式集合.

2 穩(wěn)定性分析

對方程組(8)-(10)求和得

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

證畢.

證畢.

(12)

證明: 將式(11)從1～n求和即可得式(12).

3 誤差分析

為簡單, 引入下列符號:

對于(u,r), 做如下正則性假設(shè):

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]假設(shè)u是方程(1)的解, 則有如下估計:

成立, 其中CT,ε表示常數(shù)C與T和ε有關(guān).

證明: 當t=n+1時, 方程組(2)-(4)減去方程組(5)-(7), 得

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

其中

下面依次估計Mi.根據(jù)Young不等式[13]、 Cauchy-Schwarz不等式和引理1, 得

其中

當n=0時, 有

根據(jù)Cauchy-Schwarz不等式、 引理1及Young不等式, 得

當n=0時, 對M8估計如下:

把式(17)～(27)代入式(16), 并將兩邊同乘4τ: 當n≥1時, 有

其中

當n=0時, 有

證畢.

4 數(shù)值分析

下面通過數(shù)值算例[14-15]對理論誤差估計和能量穩(wěn)定性進行驗證, 其中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