修正的Cahn-Hilliard 方程的大時(shí)間步長(zhǎng)方法

胡歡歡, 李 楊, 賈宏恩

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

1 引言

在二元合金中，為了模擬現(xiàn)象學(xué)中的失穩(wěn)分解[1,2]，Cahn 和Hilliard 于1950 年提出Cahn-Hilliard 方程.為了抑制粗化現(xiàn)象，Aristotelous 等人[3]提出了修正的Chan-Hilliard 方程.修正的Cahn-Hilliard 方程具有如下形式

并且Ω ∈Rd, d=2,3.u 是指混合物中兩種物質(zhì)之一的濃度，稱為相變量.

當(dāng)θ =0 時(shí)，方程(1)是經(jīng)典的Cahn-Hilliard 方程[4].許多學(xué)者對(duì)經(jīng)典Cahn-Hilliard方程的數(shù)值解進(jìn)行了研究，例如，Zhang 和Wang[5]提出結(jié)合凸分裂方法的全離散格式，此格式滿足質(zhì)量守恒及原始問(wèn)題的能量耗散；Guill′en-Gonz′alez 和Tierra[6]用不同方法去逼近雙勢(shì)阱項(xiàng)，并分析了關(guān)于時(shí)間分別是一階和二階的線性格式；Elliott 等人[7,8]利用非協(xié)調(diào)有限元Morely 元，得到了最佳L2誤差估計(jì)；Liu 等人[9]用傳統(tǒng)有限元和混合有限元兩重網(wǎng)格方法來(lái)解Cahn-Hilliard 方程；Du 和Nicolaides[10]提出一種有限元格式來(lái)解帶有Dirichlet 邊界條件的Cahn-Hilliard 方程，并且證明了這種格式是穩(wěn)定的.當(dāng)θ = 1 時(shí)，方程(1)為特殊的修正的Cahn-Hilliard 方程[11].此時(shí)方程(1)與經(jīng)典的Cahn-Hilliard 方程有很大的不同，修正的Cahn-Hilliard 方程仍用來(lái)描述相分離和粗化現(xiàn)象的模型[12-16].Lee 等人[17]利用隱式方法由二維截面圖像重構(gòu)三維實(shí)體模型；Gillette[18]應(yīng)用凸分裂和譜方法對(duì)修正的Cahn-Hilliard 方程進(jìn)行了研究；Choi 等人[19]利用譜方法研究了修正的Cahn-Hilliard 方程.

本文的主要工作如下：首先給出修正的Cahn-Hilliard 方程的半離散數(shù)值格式，并證明此格式的穩(wěn)定性；其次，給出全離散格式及其誤差估計(jì)；最后，通過(guò)數(shù)值算例來(lái)驗(yàn)證理論部分的正確性與有效性.

2 半離散格式及穩(wěn)定性

2.1 基本理論及符號(hào)

設(shè)L2(Ω)表示平方可積函數(shù)，其內(nèi)積和范數(shù)分別為

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

記

接下來(lái)定義H-1(Ω)，用〈·,·〉表示H-1(Ω)和H1(Ω)上的對(duì)偶內(nèi)積，記

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

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

則

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

2.2 半離散格式

對(duì)于劃分[0,T] : 0 = t0＜t1＜··· ＜tM= T, tn+1- tn= Δt = T/M，這里M ＞0 為整數(shù)，則

其中

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

其中Δt 是時(shí)間步長(zhǎng)，tn= nΔt，并且un是u(x,tn)的近似值.在數(shù)值模擬時(shí)，當(dāng)參數(shù)ν 較小時(shí)，格式(4)無(wú)法在較大的時(shí)間步長(zhǎng)上計(jì)算.為了解決這一問(wèn)題，加O(Δtut)到格式(4)中

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

2.3 穩(wěn)定性分析

在證明(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，因此有

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

利用引理1，有

定理證明完成.

3 全離散格式及誤差估計(jì)

3.1 全離散格式

且存在不依賴h 的常數(shù)c ＞0，滿足逆不等式

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

3.2 誤差分析

從不等式

可以得到

根據(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，當(dāng)h 足夠小時(shí)，則存在不依賴h, Δt, n 的C =C(u)滿足

證明 在證明之前，先給出一個(gè)先驗(yàn)假設(shè)[23]：若0 ＜h ＜h0，則存在h0滿足

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

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

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

結(jié)合方程(3),(16)和(20)，有下面的估計(jì)

其中用到下面的不等式

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

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

則有

根據(jù)離散的Gronwall 引理，則有‖μn+1‖ ≤C(h4+Δt)，因此，結(jié)合(15)和三角不等式，所以有

‖un+1-Un+1‖≤‖μn+1‖+‖ηn+1‖≤C(h4+Δt).

4 數(shù)值實(shí)驗(yàn)

本小節(jié)，利用數(shù)值算例來(lái)驗(yàn)證理論分析的準(zhǔn)確性和有效性.

4.1 誤差和收斂性

考慮二維修正的Cahn-Hilliard 方程，計(jì)算區(qū)域?yàn)棣?[0,2π]2，初值為u0=0.2 sin(x)sin(y).由于方程(1)的精確解未知，我們選取Δt = 0.0001 和N = 128 時(shí)，所計(jì)算的數(shù)值解作為精確解.

表1 給出了當(dāng)ν = 0.03 時(shí)，A, θ, Δt 取不同值時(shí)的L2誤差.通過(guò)觀察表1 的數(shù)據(jù)，可以看出當(dāng)A 取固定值時(shí)，全離散格式是穩(wěn)定的，并且關(guān)于時(shí)間是一階收斂，和理論分析一致.

表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

當(dāng)ν = 0.03 時(shí)，表2、表3、表4 分別呈現(xiàn)了A 取0,1,2 時(shí)，θ 分別取0,2,5 的空間收斂階.通過(guò)觀察表2、表3、表4，可以發(fā)現(xiàn)空間收斂階和理論分析相符.

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

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

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

4.2 能量穩(wěn)定

圖1 能量曲線

4.3 等值線

圖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 會(huì)抑制相分離.通過(guò)觀察圖2 和圖3，可知：“A-項(xiàng)”確實(shí)起到了增加時(shí)間步長(zhǎng)的作用，但不恰當(dāng)?shù)腁 值會(huì)導(dǎo)致數(shù)值解的發(fā)散，即當(dāng)計(jì)算時(shí)間步長(zhǎng)比較大時(shí)，大時(shí)間步長(zhǎng)方法得到的解可能不是真實(shí)解.

5 結(jié)論

本文中，我們對(duì)于修正的Cahn-Hilliard 方程的大時(shí)間步長(zhǎng)方法進(jìn)行了研究.為解決由于非線性和小參數(shù)帶來(lái)的影響，提出了穩(wěn)定的離散格式，從理論上證明了方法的穩(wěn)定性，并給出了誤差估計(jì).最后，通過(guò)數(shù)值實(shí)驗(yàn)證明了方法的有效性.