胡歡歡, 李 楊, 賈宏恩
(太原理工大學數(shù)學學院,太原 030024)
在二元合金中,為了模擬現(xiàn)象學中的失穩(wěn)分解[1,2],Cahn 和Hilliard 于1950 年提出Cahn-Hilliard 方程.為了抑制粗化現(xiàn)象,Aristotelous 等人[3]提出了修正的Chan-Hilliard 方程.修正的Cahn-Hilliard 方程具有如下形式
并且Ω ∈Rd, d=2,3.u 是指混合物中兩種物質(zhì)之一的濃度,稱為相變量.
當θ =0 時,方程(1)是經(jīng)典的Cahn-Hilliard 方程[4].許多學者對經(jīng)典Cahn-Hilliard方程的數(shù)值解進行了研究,例如,Zhang 和Wang[5]提出結(jié)合凸分裂方法的全離散格式,此格式滿足質(zhì)量守恒及原始問題的能量耗散;Guill′en-Gonz′alez 和Tierra[6]用不同方法去逼近雙勢阱項,并分析了關于時間分別是一階和二階的線性格式;Elliott 等人[7,8]利用非協(xié)調(diào)有限元Morely 元,得到了最佳L2誤差估計;Liu 等人[9]用傳統(tǒng)有限元和混合有限元兩重網(wǎng)格方法來解Cahn-Hilliard 方程;Du 和Nicolaides[10]提出一種有限元格式來解帶有Dirichlet 邊界條件的Cahn-Hilliard 方程,并且證明了這種格式是穩(wěn)定的.當θ = 1 時,方程(1)為特殊的修正的Cahn-Hilliard 方程[11].此時方程(1)與經(jīng)典的Cahn-Hilliard 方程有很大的不同,修正的Cahn-Hilliard 方程仍用來描述相分離和粗化現(xiàn)象的模型[12-16].Lee 等人[17]利用隱式方法由二維截面圖像重構三維實體模型;Gillette[18]應用凸分裂和譜方法對修正的Cahn-Hilliard 方程進行了研究;Choi 等人[19]利用譜方法研究了修正的Cahn-Hilliard 方程.
本文的主要工作如下:首先給出修正的Cahn-Hilliard 方程的半離散數(shù)值格式,并證明此格式的穩(wěn)定性;其次,給出全離散格式及其誤差估計;最后,通過數(shù)值算例來驗證理論部分的正確性與有效性.
設L2(Ω)表示平方可積函數(shù),其內(nèi)積和范數(shù)分別為
空間L∞(Ω)和Hm(Ω)的范數(shù)分別為
記
接下來定義H-1(Ω),用〈·,·〉表示H-1(Ω)和H1(Ω)上的對偶內(nèi)積,記
(ζ,ξ)H-1:=(?T(ζ),?T(ξ))=(ζ,T(ξ))=(T(ζ),ξ),
引理1(離散的Gronwall 引理)[20]設C0, Δt 是正數(shù),并ak, bk, ck, dk是滿足下面條件的非負序列
則
方程(1)保持能量耗散,若定義能量泛函
對于劃分[0,T] : 0 = t0<t1<··· <tM= T, tn+1- tn= Δt = T/M,這里M >0 為整數(shù),則
其中
首先考慮修正的Cahn-Hilliard 方程的一階半隱格式
其中Δt 是時間步長,tn= nΔt,并且un是u(x,tn)的近似值.在數(shù)值模擬時,當參數(shù)ν 較小時,格式(4)無法在較大的時間步長上計算.為了解決這一問題,加O(Δtut)到格式(4)中
其中A 是正常數(shù),格式(5)的弱形式為
在證明(6)的穩(wěn)定性之前,我們限制Φ′(u)滿足下面的條件[21]:存在常數(shù)L 使得
定理1 若A >0,則半隱格式(6)是穩(wěn)定的,即滿足
證明 在方程(6)中,令v =un+1,則
使用等式2a(a-b)=a2-b2+(a-b)2,因此有
化簡上述不等式并兩邊乘以2Δt,對n 從0 到k(0 ≤k ≤M -1)進行求和,得到
利用引理1,有
定理證明完成.
且存在不依賴h 的常數(shù)c >0,滿足逆不等式
類似地,定義方程(1)的半離散格式:求uh(t):(0,T]→S3h,使得
從不等式
可以得到
根據(jù)雙調(diào)和方程的有限元分析[22],得到
在方程(16)中,令t=tn+1并減去方程(13),則有
在方程(17)中
故有
定理2 記u(t)和Un+1分別是方程(11)和(13)的解,如果u(0)∈H4(Ω),滿足‖u(0)-U0‖≤Ch4‖u(0)‖4,網(wǎng)格比Δ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 不等式,可以得到
結(jié)合方程(3),(16)和(20),有下面的估計
其中用到下面的不等式
把上述不等式帶入方程(21),得到
上式從n=1 加到M,注意到
則有
根據(jù)離散的Gronwall 引理,則有‖μn+1‖ ≤C(h4+Δt),因此,結(jié)合(15)和三角不等式,所以有
‖un+1-Un+1‖≤‖μn+1‖+‖ηn+1‖≤C(h4+Δt).
本小節(jié),利用數(shù)值算例來驗證理論分析的準確性和有效性.
考慮二維修正的Cahn-Hilliard 方程,計算區(qū)域為Ω=[0,2π]2,初值為u0=0.2 sin(x)sin(y).由于方程(1)的精確解未知,我們選取Δt = 0.0001 和N = 128 時,所計算的數(shù)值解作為精確解.
表1 給出了當ν = 0.03 時,A, θ, Δt 取不同值時的L2誤差.通過觀察表1 的數(shù)據(jù),可以看出當A 取固定值時,全離散格式是穩(wěn)定的,并且關于時間是一階收斂,和理論分析一致.
表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
續(xù)表1 L2 誤差:ν =0.03, T =1, h=
續(xù)表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 分別呈現(xiàn)了A 取0,1,2 時,θ 分別取0,2,5 的空間收斂階.通過觀察表2、表3、表4,可以發(fā)現(xiàn)空間收斂階和理論分析相符.
表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,不難發(fā)現(xiàn),θ >0 會抑制相分離.通過觀察圖2 和圖3,可知:“A-項”確實起到了增加時間步長的作用,但不恰當?shù)腁 值會導致數(shù)值解的發(fā)散,即當計算時間步長比較大時,大時間步長方法得到的解可能不是真實解.
本文中,我們對于修正的Cahn-Hilliard 方程的大時間步長方法進行了研究.為解決由于非線性和小參數(shù)帶來的影響,提出了穩(wěn)定的離散格式,從理論上證明了方法的穩(wěn)定性,并給出了誤差估計.最后,通過數(shù)值實驗證明了方法的有效性.