王曉瑩, 劉 雪, 江 山

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

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

1 問題提出

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

(1)

其中u為原時空問題的真解,x為空間變量,t為時間變量,α,β,χ為方程系數(shù),f為方程右端項．設(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ù)v, 在區(qū)間I=[a,b]進行積分,得

(4)

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

(5)

(6)

2.2 半離散Galerkin格式

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

(7)

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

(8)

(9)

(10)

(11)

(12)

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

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

(13)

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

(14)

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

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

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

(15)

3 時間尺度的全離散格式

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

(16)

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

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

(17)

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

(18)

(19)

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

4 數(shù)值實驗

(20)

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

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

(21)

(22)

(23)

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

(24)

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

將空間步長h與時間步長ht的關(guān)系固定為ht=h2, 采用線性有限元法分別結(jié)合θ=1、θ=0.5兩種時間全離散隱格式, 計算問題(20)在空間剖分?jǐn)?shù)N不斷偶數(shù)倍增大時的數(shù)值解uh,并根據(jù)式(21)～(24)分別計算兩種差分格式下數(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é)合時間層計算的可靠性和穩(wěn)定性;θ=0.5的Crank-Nicolson差分格式下數(shù)值解的逼近效果稍優(yōu)于θ=1的向后差分格式．

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

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

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

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

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

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

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

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

圖2 當(dāng)θ=0.5、t=1 s時, 線性有限元(a)與二次有限元(b)在N=32下的絕對誤差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é)合時間上的Crank-Nicolson有限差分隱格式, 精確有效地模擬了一維非穩(wěn)態(tài)對流擴散反應(yīng)問題．?dāng)?shù)值算例表明,在一致網(wǎng)格進行等距的時空間步長離散后可以獲得理想結(jié)果,二次有限元結(jié)合Crank-Nicolson有限差分隱格式求解時的計算精度和收斂階數(shù)更高,應(yīng)用優(yōu)勢更明顯．