胡歡歡, 李 楊, 賈宏恩
(太原理工大學數學學院,太原 030024)
在二元合金中,為了模擬現象學中的失穩分解[1,2],Cahn 和Hilliard 于1950 年提出Cahn-Hilliard 方程.為了抑制粗化現象,Aristotelous 等人[3]提出了修正的Chan-Hilliard 方程.修正的Cahn-Hilliard 方程具有如下形式


并且Ω ∈Rd, d=2,3.u 是指混合物中兩種物質之一的濃度,稱為相變量.
當θ =0 時,方程(1)是經典的Cahn-Hilliard 方程[4].許多學者對經典Cahn-Hilliard方程的數值解進行了研究,例如,Zhang 和Wang[5]提出結合凸分裂方法的全離散格式,此格式滿足質量守恒及原始問題的能量耗散;Guill′en-Gonz′alez 和Tierra[6]用不同方法去逼近雙勢阱項,并分析了關于時間分別是一階和二階的線性格式;Elliott 等人[7,8]利用非協調有限元Morely 元,得到了最佳L2誤差估計;Liu 等人[9]用傳統有限元和混合有限元兩重網格方法來解Cahn-Hilliard 方程;Du 和Nicolaides[10]提出一種有限元格式來解帶有Dirichlet 邊界條件的Cahn-Hilliard 方程,并且證明了這種格式是穩定的.當θ = 1 時,方程(1)為特殊的修正的Cahn-Hilliard 方程[11].此時方程(1)與經典的Cahn-Hilliard 方程有很大的不同,修正的Cahn-Hilliard 方程仍用來描述相分離和粗化現象的模型[12-16].Lee 等人[17]利用隱式方法由二維截面圖像重構三維實體模型;Gillette[18]應用凸分裂和譜方法對修正的Cahn-Hilliard 方程進行了研究;Choi 等人[19]利用譜方法研究了修正的Cahn-Hilliard 方程.
本文的主要工作如下:首先給出修正的Cahn-Hilliard 方程的半離散數值格式,并證明此格式的穩定性;其次,給出全離散格式及其誤差估計;最后,通過數值算例來驗證理論部分的正確性與有效性.
設L2(Ω)表示平方可積函數,其內積和范數分別為

空間L∞(Ω)和Hm(Ω)的范數分別為

記

接下來定義H-1(Ω),用〈·,·〉表示H-1(Ω)和H1(Ω)上的對偶內積,記


(ζ,ξ)H-1:=(?T(ζ),?T(ξ))=(ζ,T(ξ))=(T(ζ),ξ),

引理1(離散的Gronwall 引理)[20]設C0, Δt 是正數,并ak, bk, ck, dk是滿足下面條件的非負序列

則

方程(1)保持能量耗散,若定義能量泛函

對于劃分[0,T] : 0 = t0<t1<··· <tM= T, tn+1- tn= Δt = T/M,這里M >0 為整數,則

其中

首先考慮修正的Cahn-Hilliard 方程的一階半隱格式

其中Δt 是時間步長,tn= nΔt,并且un是u(x,tn)的近似值.在數值模擬時,當參數ν 較小時,格式(4)無法在較大的時間步長上計算.為了解決這一問題,加O(Δtut)到格式(4)中

其中A 是正常數,格式(5)的弱形式為

在證明(6)的穩定性之前,我們限制Φ′(u)滿足下面的條件[21]:存在常數L 使得

定理1 若A >0,則半隱格式(6)是穩定的,即滿足

證明 在方程(6)中,令v =un+1,則

使用等式2a(a-b)=a2-b2+(a-b)2,因此有


化簡上述不等式并兩邊乘以2Δt,對n 從0 到k(0 ≤k ≤M -1)進行求和,得到

利用引理1,有

定理證明完成.


且存在不依賴h 的常數c >0,滿足逆不等式



類似地,定義方程(1)的半離散格式:求uh(t):(0,T]→S3h,使得



從不等式

可以得到

根據雙調和方程的有限元分析[22],得到



在方程(16)中,令t=tn+1并減去方程(13),則有

在方程(17)中

故有

定理2 記u(t)和Un+1分別是方程(11)和(13)的解,如果u(0)∈H4(Ω),滿足‖u(0)-U0‖≤Ch4‖u(0)‖4,網格比Δt/h2≤c,當h 足夠小時,則存在不依賴h, Δt, n 的C =C(u)滿足

證明 在證明之前,先給出一個先驗假設[23]:若0 <h <h0,則存在h0滿足

在方程(18)中令vh=μn+1,則有

使用Cauchy 不等式和Young 不等式,所以


利用δtμn+1的定義和Cauchy 不等式,可以得到

結合方程(3),(16)和(20),有下面的估計

其中用到下面的不等式

把上述不等式帶入方程(21),得到

上式從n=1 加到M,注意到

則有


根據離散的Gronwall 引理,則有‖μn+1‖ ≤C(h4+Δt),因此,結合(15)和三角不等式,所以有
‖un+1-Un+1‖≤‖μn+1‖+‖ηn+1‖≤C(h4+Δt).
本小節,利用數值算例來驗證理論分析的準確性和有效性.
考慮二維修正的Cahn-Hilliard 方程,計算區域為Ω=[0,2π]2,初值為u0=0.2 sin(x)sin(y).由于方程(1)的精確解未知,我們選取Δt = 0.0001 和N = 128 時,所計算的數值解作為精確解.
表1 給出了當ν = 0.03 時,A, θ, Δt 取不同值時的L2誤差.通過觀察表1 的數據,可以看出當A 取固定值時,全離散格式是穩定的,并且關于時間是一階收斂,和理論分析一致.
表1 L2 誤差:ν =0.03, T =1, h=

表1 L2 誤差:ν =0.03, T =1, h=
A Δt θ =0 θ =2 θ =5 0 0.01 0.00382036 0.0152565 0.167332 0.005 0.00212644 0.00702886 0.0877675 0.0025 0.00108075 0.00339221 0.0441148 0.00125 0.00052603 0.00168505 0.0213931 1 0.01 0.0107831 0.0530077 0.311978 0.005 0.00623478 0.0219408 0.178414 0.0025 0.00340485 0.00992049 0.0943016 0.00125 0.00173112 0.00473141 0.047645
續表1 L2 誤差:ν =0.03, T =1, h=

續表1 L2 誤差:ν =0.03, T =1, h=
A Δt θ =0 θ =2 θ =5 2 0.01 0.0189934 0.102308 0.414267 0.005 0.00985323 0.0402944 0.2551 0.0025 0.00554085 0.0172588 0.140868 0.00125 0.00291495 0.0079702 0.0730369
當ν = 0.03 時,表2、表3、表4 分別呈現了A 取0,1,2 時,θ 分別取0,2,5 的空間收斂階.通過觀察表2、表3、表4,可以發現空間收斂階和理論分析相符.

表2 收斂階:ν =0.03, θ =0

表3 收斂階:ν =0.03, θ =2

表4 收斂階:ν =0.03, θ =5

圖1 能量曲線


圖2 等值線:ν =0.03, θ =2,左:A=0, Δt=0.001,中:A=1, Δt=0.005,右:A=2, Δt=0.01

圖3 等值線:ν =0.03, A=1,左:θ =0, Δt=0.001,中:θ =2, Δt=0.005,右L:θ =5, Δt=0.01
在圖2 中,第1 列,θ = 2, A = 0, Δt = 0.001;第2 列,θ = 2, A = 1, Δt =0.005;第3 列,θ = 2, A = 2, Δt = 0.01.同樣地,在圖3 中,第1 列,θ = 0, A =1, Δt = 0.001;第2 列,θ = 2, A = 1, Δt = 0.005;第3 列,θ = 5, A = 1, Δt =0.01.觀察圖3,不難發現,θ >0 會抑制相分離.通過觀察圖2 和圖3,可知:“A-項”確實起到了增加時間步長的作用,但不恰當的A 值會導致數值解的發散,即當計算時間步長比較大時,大時間步長方法得到的解可能不是真實解.
本文中,我們對于修正的Cahn-Hilliard 方程的大時間步長方法進行了研究.為解決由于非線性和小參數帶來的影響,提出了穩定的離散格式,從理論上證明了方法的穩定性,并給出了誤差估計.最后,通過數值實驗證明了方法的有效性.