二維非穩(wěn)態(tài)對流擴散方程的高階緊致差分格式

魏劍英，王燕，葛永斌

寧夏大學(xué) 數(shù)學(xué)統(tǒng)計學(xué)院，銀川 750021

非穩(wěn)態(tài)對流擴散方程是一類基本的發(fā)展方程，在生態(tài)環(huán)境、流體力學(xué)、生物數(shù)學(xué)等領(lǐng)域都有廣泛應(yīng)用．由于解析解通常很難求出，因此，研究各種有效實用的數(shù)值方法對此類方程進行求解就顯得極為重要。

有關(guān)非穩(wěn)態(tài)對流擴散方程的有限差分法的文獻報道已有很多，如文獻[1]利用半離散法和Padé逼近，推導(dǎo)了一維方程的截斷誤差為O(τ5+h4)的隱格式，且是無條件穩(wěn)定的； 文獻[2]利用Taylor級數(shù)展開和待定系數(shù)法構(gòu)造了一致四階緊致格式并結(jié)合三階TVD-Runge-Kutta方法求解了一維方程； 文獻[3]針對二維方程，對空間導(dǎo)數(shù)項和時間導(dǎo)數(shù)項分別采用四階離散和加權(quán)平均的方法，推導(dǎo)出一種截斷誤差為O(τ2+h4)的格式； 文獻[4]采用半離散的方法，從一維定常問題出發(fā)，利用四階緊致差分算子和Crank-Nicolson(C-N)格式，推導(dǎo)了一種求解二維常系數(shù)非穩(wěn)態(tài)對流擴散方程的差分格式，截斷誤差為O(τ2+h4)，且是無條件穩(wěn)定的； 文獻[5]針對此類方程，構(gòu)造了有理型緊致交替方向隱式(ADI)格式，截斷誤差為O(τ2+h4)，且是無條件穩(wěn)定的； 文獻[6]還提出了該方程無條件穩(wěn)定的指數(shù)型緊致ADI差分格式，其截斷誤差依然為O(τ2+h4)； 文獻[7]利用指數(shù)變換消除對流項后轉(zhuǎn)化為擴散方程，再利用緊致公式和擴展的Simpson公式構(gòu)造了一維方程的高階格式，截斷誤差為O(τ4+h4)，且是無條件穩(wěn)定的； 文獻[8]針對二維非穩(wěn)態(tài)常系數(shù)方程，空間方向直接采用六階組合緊致差分公式進行計算，時間方向用C-N格式離散，所提格式的截斷誤差為O(τ2+h6)，盡管該格式空間達到了六階精度，但是由于其時間只有二階精度，因此為了保證空間精度達到六階，其計算所需要采取的時間步長必須為O(h3)，即必須采用較小的時間步長，并且該方法僅適用于常系數(shù)問題。

本文針對二維非穩(wěn)態(tài)變系數(shù)對流擴散方程，首先對空間二階導(dǎo)數(shù)采用一階導(dǎo)數(shù)的四階逼近公式，而一階導(dǎo)數(shù)采用四階Padé逼近，時間項采用二階向后差分公式(BDF)，得到一個時間二階、空間四階精度的緊致差分格式； 然后利用Taylor級數(shù)展開，對一階和二階空間導(dǎo)數(shù)項采用六階緊致差分公式，時間項采用三階BDF得到一個時間三階、空間六階精度的緊致差分格式．由于采用了向后差分，因此兩種格式均是無條件穩(wěn)定的。

1 差分格式的建立

考慮如下二維非穩(wěn)態(tài)對流擴散方程：

ut-α(uxx+uyy)+p(x，y，t)ux+q(x，y，t)uy=f(x，y，t) (x，y)∈Ω，t≥0

(1)

初始條件：

u(x，y，0)=g(x，y) (x，y)∈Ω

(2)

邊界條件：

u(x，y，t)=s(x，y，t) (x，y)∈?Ω，t>0

(3)

其中：Ω∈[a1，a2]×[b1，b2]為R2上的矩形區(qū)域，?Ω為Ω的邊界，u(x，y，t)為待求未知量，p(x，y，t)，q(x，y，t)分別為x，y方向的對流項系數(shù)，α為擴散項系數(shù)(常數(shù))，f(x，y，t)為源項，且p(x，y，t)，q(x，y，t)，f(x，y，t)，g(x，y)，s(x，y，t)均為已知函數(shù)。

(4)

(5)

1.1 四階精度緊致格式

將方程(1)改寫為

ut=f(x，y，t)+α(uxx+uyy)-p(x，y，t)ux-q(x，y，t)uy

(6)

對二階導(dǎo)數(shù)uxx和uyy采用如下逼近：

(7)

(8)

將(7)，(8)式代入(6)式，考慮其在點(xi，yj)的值，得到

(9)

為了使(9)式空間保持四階精度，對ux和uy的計算采用如下四階Padé公式[9]

(10)

(11)

考慮(9)式在第n+1時間層的離散值，對ut采用如下二階BDF[10-11]進行離散

(12)

舍去高階項，可得

(13)

(14)

其中G為(9)式的右端項。

在文獻[12]中已經(jīng)證明當(dāng)K≤6(K代表精度)時，BDF是無條件穩(wěn)定的．另外，由式(13)可以看出，在每個時間步上，該格式為5點模板的全隱格式．為了得到u(x，y，t)的計算值，(14)式中出現(xiàn)的ux和uy除需知道其內(nèi)點值，還需知道其在邊界點處的值，為保持與內(nèi)點差分格式同樣的精度，對邊界上的ux，uy采用一致四階邊界條件[2]計算。

1.2 六階精度緊致格式

為了得到方程(1)的六階精度的格式，需由Taylor級數(shù)展開得到一、二階導(dǎo)數(shù)的六階近似公式

(15)

(16)

(17)

(18)

將(15)-(18)式代入(6)式，有

(19)

(20)

(21)

(22)

(23)

(24)

(25)

將(20)-(25)式代入(19)式，考慮其在點(xi，yj)的值，即可得到如下半離散的緊致格式

(26)

為了使(26)式具有六階精度，對ux和uy采用文獻[9]中的六階差分公式計算，uxx和uyy采用文獻[13]中的六階差分公式計算

(27)

(28)

(29)

(30)

考慮(26)式在第n+1時間層的值，對ut采用三階BDF[10]進行離散

(31)

略去高階項，可得

(32)

2 數(shù)值實驗

本文將利用以下4個問題進行數(shù)值實驗，驗證(2，4)格式和(3，6)格式的穩(wěn)定性和有效性，文中涉及到的誤差及收斂階定義如下：

1) 最大絕對誤差：

2)L2范數(shù)誤差：

3) 收斂階：

問題1考慮如下非齊次對流擴散問題

其精確解為u(x，t)=e-txy(1-x)(1-y)，初邊值條件及右端項f(x，y，t)均由精確解給出。

表1給出了問題1當(dāng)τ=h2，t=0.5時的計算結(jié)果．從表中可以看出，本文(2，4)格式和文獻[6]格式在空間上均達到了四階精度，但本文(2，4)格式的計算結(jié)果比BTCS格式和文獻[6]格式更好； 本文(3，6)格式在空間上達到了六階精度(網(wǎng)格為64時的數(shù)值已經(jīng)達到機器精度，所以收斂階受到影響不準(zhǔn)確)，且最大絕對誤差比(2，4)格式的小2～3個數(shù)量級，計算結(jié)果更精確．當(dāng)h=0.01，t=0.5時，不同時間步長τ下的L∞誤差、收斂階及CPU時間由表2給出．從表中可以看出，本文(2，4)格式時間上達到的精度為二階，本文(3，6)格式時間上達到的精度為三階，這與理論分析一致； 另外當(dāng)τ較大時計算收斂慢，所用的CPU時間相對較長，計算過程中(3，6)格式需要計算在節(jié)點處的一、二階導(dǎo)數(shù)值，因此所用的計算時間比(2，4)格式長。

表1 問題1當(dāng)τ=h2，t=0.5時的最大絕對誤差L∞及收斂階

表2 問題1當(dāng)h=0.01，t=0.5時的最大絕對誤差L∞，收斂階及CPU時間

問題2考慮如下變系數(shù)對流擴散問題

其精確解為u(x，t)=txy(1-x)(1-y)ex+y，初邊值條件及右端項f(x，y，t)均由精確解給出。

表3給出了問題2當(dāng)τ=h2，t=0.25時本文格式與C-N格式和BTCS格式的L∞誤差及收斂階的比較．從表中可以看出，C-N格式和BTCS格式在空間上達到的精度均為二階，本文(2，4)格式在空間上達到的精度為四階，而本文(3，6)格式空間上達到的精度為六階，且(3，6)格式的計算結(jié)果誤差更小，充分驗證了本文兩種格式的精確性和有效性．另外，上述兩個例子的計算結(jié)果進一步說明低精度的啟動步并不影響格式的整體精度。

表3 問題2當(dāng)τ=h2，t=0.25時的L∞誤差及收斂階

問題3考慮如下高斯脈沖問題

圖1 t=1，τ=2.5×10-3，Pe=2，1.2≤x，y≤1.8時的等值線圖

圖2 t=0.001，τ=2.5×10-6，Pe=2000，1.2≤x，y≤1.8時的等值線圖

表4 問題3當(dāng)h=0.02時在不同參數(shù)下的L∞誤差和L2誤差

問題4考慮如下非線性方程

表5給出了問題4當(dāng)t=1，τ=h2時，本文(2，4)格式和(3，6)格式在α分別取1和0.1時的最大絕對誤差與收斂階．從表中可以看出，本文(2，4)格式在空間上達到的精度是四階，(3，6)格式達到的精度是六階(網(wǎng)格為64，α=1時最大絕對誤差的值已經(jīng)達到機器精度)，這與理論推導(dǎo)是一致的，從表中還可以看出，α=1時所用的CPU時間比α=0.1時的長，(3，6)格式的CPU時間比(2，4)格式的長．表6給出了當(dāng)N=7，τ=0.01，α=1時，t=15和t=20的計算解與精確解的絕對誤差．表中數(shù)據(jù)表明，本文(2，4)格式的計算誤差比文獻[14]格式小5～6個數(shù)量級，與文獻[15]格式的計算結(jié)果大致相當(dāng)； 本文(3，6)格式的計算結(jié)果比文獻[14]格式小7～9個數(shù)量級，比文獻[15]格式小2～3個數(shù)量級．因此，對于此類非齊次邊界的Burgers方程，本文兩種高精度緊致差分格式非常有效．圖3給出了N=20，t=1，τ=0.01，α=0.1時本文兩種格式的絕對誤差和文獻[15]格式的絕對誤差比較，從圖中可以看出，本文(2，4)格式的絕對誤差與文獻[15]格式的絕對誤差大致相當(dāng)，但本文(3，6)格式的計算結(jié)果更加精確，絕對誤差更小。

表5 問題 4 當(dāng)τ=h2，t=1時的最大絕對誤差L∞，收斂階和CPU時間

表6 問題4當(dāng)N=7，τ=0.01，α=1時的絕對誤差

續(xù)表6

圖3 N=20，t=1，τ=0.01，α=0.1時的精確解與絕對誤差

3 結(jié)論