有限元Crank-Nicolson格式高階求解非穩(wěn)態(tài)擴(kuò)散方程

王曉瑩, 劉 雪, 江 山

(南通大學(xué)理學(xué)院, 江蘇 南通 226019)

非穩(wěn)態(tài)擴(kuò)散問題廣泛存在于流體力學(xué)、化工制造、環(huán)境科學(xué)和能源開發(fā)等領(lǐng)域．在實(shí)際應(yīng)用中, 主要采用有限差分法[1-2]、有限體積法[3-4]、譜方法[5-6]或有限元法[7-10]等數(shù)值方法進(jìn)行求解．上述方法大多需要一定的限制條件, 如Gowrisankar等[1]在解決非穩(wěn)態(tài)拋物型初邊值問題時(shí), 須運(yùn)用迎風(fēng)有限差分格式在時(shí)間層的等分布原理; Quenjel[3]提出的全隱格式有限體積法, 需要在多邊形網(wǎng)格上才能更好地保證數(shù)值解的非負(fù)性．然而, 有限元法的優(yōu)勢(shì)在于可以針對(duì)不同的微分方程建立對(duì)應(yīng)的變分形式以適應(yīng)不同的問題, 在有限單元離散剖分背景下參考單元與實(shí)際單元之間形成必要的仿射變換, 并通過獨(dú)立構(gòu)造的多項(xiàng)式基函數(shù)有效耦合局部與全局間的空間尺度信息, 得到精度更好、收斂更快的模擬結(jié)果．如Cheng等[7]利用局部間斷有限元法處理奇異攝動(dòng)非穩(wěn)態(tài)問題, 給出全離散時(shí)間在三階龍格-庫(kù)塔法的誤差分析; Lin等[8]提出一種對(duì)流擴(kuò)散反應(yīng)方程的流線迎風(fēng)Petrov-Galerkin(streamline-upwind Petrov-Galerkin, SUPG)穩(wěn)定化時(shí)空有限元法, 得到了時(shí)空數(shù)值格式的誤差估計(jì); Varma等[10]針對(duì)非穩(wěn)態(tài)對(duì)流擴(kuò)散反應(yīng)方程, 分析了自適應(yīng)有限元法的后驗(yàn)誤差．近年來, 利用多種數(shù)值方法相結(jié)合進(jìn)行優(yōu)勢(shì)互補(bǔ)的措施備受關(guān)注．如Das等[11]針對(duì)奇異攝動(dòng)延遲拋物方程, 在Shishkin網(wǎng)格應(yīng)用有限元與有限差分的混合計(jì)算格式得到一致收斂; Li[12], Mekonnen[13]等分別利用迎風(fēng)格式的本征正交分解與混合有限差分格式有效模擬了雙參數(shù)非穩(wěn)態(tài)模型．此外, 通過構(gòu)造特殊函數(shù)與邊界元技術(shù), 在二維情形緊湊型的高階計(jì)算格式也取得相關(guān)進(jìn)展[14-15]．截至目前, 尚未出現(xiàn)文獻(xiàn)報(bào)道在空間上應(yīng)用高次有限元法與時(shí)間上應(yīng)用θ=0.5的Crank-Nicolson有限差分隱格式相結(jié)合求解時(shí)空間非穩(wěn)態(tài)的對(duì)流擴(kuò)散反應(yīng)方程．本文擬基于高次有限元法結(jié)合Crank-Nicolson有限差分格式求解非穩(wěn)態(tài)對(duì)流擴(kuò)散反應(yīng)方程,以期獲得精確化的高階一致收斂模擬結(jié)果．

1 問題提出

考慮非穩(wěn)態(tài)對(duì)流擴(kuò)散反應(yīng)方程

(1)

其中u為原時(shí)空問題的真解,x為空間變量,t為時(shí)間變量,α,β,χ為方程系數(shù),f為方程右端項(xiàng)．設(shè)定初邊值條件:

u(x,0)=u0,x∈[a,b];

(2)

u(a,t)=ua,u(b,t)=ub, 0

(3)

本文研究目的是得到真解u的精確化高階數(shù)值逼近結(jié)果uh．

2 空間尺度的有限元法

2.1 有限元變分形式

根據(jù)虛功原理, 在問題(1)左右兩邊同時(shí)乘以檢驗(yàn)函數(shù)v, 在區(qū)間I=[a,b]進(jìn)行積分,得

(4)

限定檢驗(yàn)函數(shù)滿足v(a)=v(b)=0, 在此基礎(chǔ)上利用分部積分化簡(jiǎn)式(4)得

(5)

(6)

2.2 半離散Galerkin格式

Galerkin方法通過有限維空間近似無限維空間．記有限元空間Uh為無限維空間H1(I)的子空間, 即Uh?H1(I), 再記uh為空間的有限元解函數(shù), 其中h為空間步長(zhǎng), 則問題(1)的半離散Galerkin格式為: 尋找uh∈H1(0,t;Uh), 滿足

(7)

在半離散Galerkin格式(7)中, 有限維空間Uh是在有限元格式下由分片多項(xiàng)式構(gòu)造的基函數(shù)張成．本文中, 有限單元上的基函數(shù)通過“參考元→局部元→全局元”的模式構(gòu)造, 現(xiàn)以二次基函數(shù)為例描述其構(gòu)造過程．

(8)

(9)

(10)

(11)

(12)

由此可得局部有限元空間Sh(Ihn)=span{ψn1,ψn2,ψn3}．

4) 在剖分形成的每一個(gè)全局節(jié)點(diǎn)xj(j=1,2,…,Nb)處定義全局基函數(shù), 使得φj|Ihn∈Sh(Ihn)且滿足條件

(13)

其中φj(j=1,2,…,Nb)為所張成有限維空間的基函數(shù)數(shù)目．類似地, 可構(gòu)造線性有限元基函數(shù)及其張成的有限維空間,且線性有限元基函數(shù)張成有限維空間時(shí), 多項(xiàng)式次數(shù)和節(jié)點(diǎn)數(shù)據(jù)對(duì)應(yīng)關(guān)系較二次有限元更簡(jiǎn)單．

(14)

其待定系數(shù)為u1(t),u2(t),…,uNb(t)．

將方程組(14)進(jìn)一步以向量和矩陣形式表示, 則可得到對(duì)應(yīng)的代數(shù)系統(tǒng)為

MX′(t)+A(t)X(t)=b(t),

(15)

3 時(shí)間尺度的全離散格式

采用有限差分法對(duì)系統(tǒng)(15)的時(shí)間尺度進(jìn)行全離散．將問題(1)中時(shí)間(0,T]一致剖分成Nm個(gè)步長(zhǎng)為ht的時(shí)間層單元, 時(shí)間節(jié)點(diǎn)為tm=mht,m=0,1,…,Nm．于是, 經(jīng)過時(shí)間尺度離散, 可得系統(tǒng)(15)的全離散格式

(16)

其中θ為具體格式的參數(shù), 其取值決定著時(shí)間層的差分格式．當(dāng)θ=0時(shí)為向前差分格式, 這種顯格式雖可直接迭代且無須解方程組, 但其穩(wěn)定性差,會(huì)導(dǎo)致誤差傳播不可控; 而當(dāng)θ=1、θ=0.5時(shí)分別為向后差分格式和Crank-Nicolson六點(diǎn)對(duì)稱格式, 均為隱格式,雖然都需要求解方程組,但是穩(wěn)定性和數(shù)值精度更佳．故本文分別選擇向后差分格式和Crank-Nicolson對(duì)稱格式進(jìn)行時(shí)間離散．

令θ≠0, 定義Xm+θ=θXm+1+(1-θ)Xm, 則有

(17)

將式(17)代入全離散Galerkin格式(16), 化簡(jiǎn)得

(18)

(19)

求解系統(tǒng)(19), 可得未知向量組Xm+θ．

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

(20)

為了清晰有效地對(duì)比驗(yàn)證數(shù)值解的真實(shí)效果,在空間處理的有限元格式分別采用線性有限元法和二次有限元法, 由此形成從局部單元到全局單元的構(gòu)造與組裝．在時(shí)間處理的有限差分格式分別取θ=1及θ=0.5的2種隱格式差分, 由此對(duì)全離散系統(tǒng)(19)進(jìn)行求解, 最終得到數(shù)值解uh．為了度量uh對(duì)方程(20)的真解u的逼近效果, 在t時(shí)刻分別計(jì)算二者誤差的L∞范數(shù)、L2范數(shù)及H1半范數(shù):

‖u-uh‖∞=supx∈I|u-uh|,

(21)

(22)

(23)

再利用上下層空間剖分形成的粗網(wǎng)格和細(xì)網(wǎng)格對(duì)應(yīng)的誤差范數(shù)值EN,E2N, 計(jì)算3種范數(shù)各自的收斂階數(shù)

(24)

進(jìn)而根據(jù)收斂階數(shù)驗(yàn)證數(shù)值方法的穩(wěn)定性和收斂速度．

將空間步長(zhǎng)h與時(shí)間步長(zhǎng)ht的關(guān)系固定為ht=h2, 采用線性有限元法分別結(jié)合θ=1、θ=0.5兩種時(shí)間全離散隱格式, 計(jì)算問題(20)在空間剖分?jǐn)?shù)N不斷偶數(shù)倍增大時(shí)的數(shù)值解uh,并根據(jù)式(21)～(24)分別計(jì)算兩種差分格式下數(shù)值解與真解之間的3種誤差范數(shù)及相應(yīng)收斂階, 結(jié)果如表1～2所示．由表1～2可見:兩種差分格式下, 3種誤差范數(shù)值隨空間剖分?jǐn)?shù)N偶數(shù)倍增大而不斷遞減,表明數(shù)值解與真解的逼近越來越好;L∞范數(shù)和L2范數(shù)均有清晰的二階收斂,H1半范數(shù)有清晰的一階收斂,保證了線性有限元法結(jié)合時(shí)間層計(jì)算的可靠性和穩(wěn)定性;θ=0.5的Crank-Nicolson差分格式下數(shù)值解的逼近效果稍優(yōu)于θ=1的向后差分格式．

表1 θ=1時(shí)線性有限元的誤差范數(shù)和收斂階Tab.1 Error norm and convergence order of linear finite elements, θ=1

表2 θ=0.5時(shí)線性有限元的誤差范數(shù)和收斂階Tab.2 Error norm and convergence order of linear finite elements, θ=0.5

類似地, 取空間步長(zhǎng)h與時(shí)間步長(zhǎng)ht的關(guān)系為ht=h2, 當(dāng)空間剖分?jǐn)?shù)N偶數(shù)倍增大時(shí), 采用二次有限元法分別結(jié)合θ=1、θ=0.5兩種時(shí)間全離散隱格式求解問題(20), 并計(jì)算相應(yīng)的誤差范數(shù)和收斂階, 結(jié)果如表3～4所示．由表3～4可見:θ=1的向后差分格式下,L∞范數(shù)、L2范數(shù)及H1半范數(shù)均有二階收斂;θ=0.5的Crank-Nicolson對(duì)稱差分格式下,L∞范數(shù)、L2范數(shù)達(dá)到三階收斂, 而H1半范數(shù)也能達(dá)二階收斂; 應(yīng)用二次有限元法結(jié)合兩種不同的時(shí)間全離散隱格式求解時(shí),θ=0.5的收斂性結(jié)果明顯優(yōu)于θ=1時(shí)的, 二次有限元法在θ=0.5的Crank-Nicolson隱格式下精度更高,收斂更快．

表3 θ=1時(shí)二次有限元的誤差范數(shù)和收斂階Tab.3 Error norm and convergence order of quadratic finite elements, θ=1

表4 θ=0.5時(shí)二次有限元的誤差范數(shù)和收斂階Tab.4 Error norm and convergence order of quadratic finite elements, θ=0.5

分別選定ht=h2,θ=1和ht=h2,θ=0.5, 對(duì)比分析線性有限元法和二次有限元法在空間尺度求解問題(20)時(shí)的誤差精度及收斂階．由表1和表3可知, 二次有限元法的誤差范數(shù)值明顯優(yōu)于線性有限元法的, 其結(jié)果更為精確, 誤差精度更高．由表2和表4可知, 二次有限元法的數(shù)值精度較線性有限元法高103倍或以上, 且3種范數(shù)的收斂階全部提升了一階,L∞范數(shù)、L2范數(shù)達(dá)三階,H1半范數(shù)達(dá)二階．

圖1給出了t=1 s時(shí)分別采用線性有限法和二次有限元法結(jié)合Crank-Nicolson對(duì)稱差分格式求解問題(20)得到的數(shù)值解uh與真解u．由圖1可見,無論是采用線性有限元法還是二次有限元法,在最簡(jiǎn)便且較稀疏的一致網(wǎng)格剖分?jǐn)?shù)N=32上計(jì)算皆可實(shí)現(xiàn)數(shù)值解uh對(duì)真解u的完美逼近．其中,由于二次有限元法在對(duì)空間尺度進(jìn)行剖分時(shí)自由度加入了單元中點(diǎn), 所以使得數(shù)值解uh在圖1(b)顯示更密集, 這表明本文方法求解一維非穩(wěn)態(tài)對(duì)流擴(kuò)散反應(yīng)問題的有效性與穩(wěn)定性較高．

圖1 當(dāng)θ=0.5、t=1 s時(shí),線性有限元(a)與二次有限元(b)在N=32下的真解和數(shù)值解Fig.1 When θ=0.5, t=1 s, the exact solution and the numerical solution of linear finite element (a) and quadratic finite element (b) on N=32

圖2給出了真解u與數(shù)值解uh之間的絕對(duì)誤差|u-uh|．由圖2可見, 在t=1 s時(shí)采用線性有限元法和二次有限元法模擬結(jié)果的精確度均較高．由于二次有限元法是基于線性有限元法, 在空間尺度進(jìn)行單元剖分時(shí)再取同一單元步長(zhǎng)的中點(diǎn)作為單元新節(jié)點(diǎn),所以依次連接所有離散節(jié)點(diǎn)的誤差值時(shí)會(huì)出現(xiàn)圖2(b)所示的上下振蕩情況,且二次有限元法的誤差在較粗剖分?jǐn)?shù)N=32上更是達(dá)到了10-8數(shù)量級(jí), 再次驗(yàn)證二次有限元求解非穩(wěn)態(tài)問題的優(yōu)越性．

圖2 當(dāng)θ=0.5、t=1 s時(shí), 線性有限元(a)與二次有限元(b)在N=32下的絕對(duì)誤差Fig.2 When θ=0.5, t=1 s, the absolute error of linear finite element (a) and quadratic finite element (b) on N=32

綜上分析,本文利用空間上的二次有限元法,結(jié)合時(shí)間上的Crank-Nicolson有限差分隱格式, 精確有效地模擬了一維非穩(wěn)態(tài)對(duì)流擴(kuò)散反應(yīng)問題．?dāng)?shù)值算例表明,在一致網(wǎng)格進(jìn)行等距的時(shí)空間步長(zhǎng)離散后可以獲得理想結(jié)果,二次有限元結(jié)合Crank-Nicolson有限差分隱格式求解時(shí)的計(jì)算精度和收斂階數(shù)更高,應(yīng)用優(yōu)勢(shì)更明顯．