一種求解三維非穩(wěn)態(tài)對流擴(kuò)散反應(yīng)方程的高精度有限差分格式*

魏劍英， 葛永斌

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

引 言

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

有限差分方法作為經(jīng)典實(shí)用的數(shù)值方法之一，在計(jì)算流體力學(xué)中發(fā)揮著重要作用.目前關(guān)于非穩(wěn)態(tài)對流擴(kuò)散方程已經(jīng)有非常多的研究報(bào)道，如文獻(xiàn)[1-11].而關(guān)于非穩(wěn)態(tài)對流擴(kuò)散反應(yīng)方程的報(bào)道相對還較少，并且大多是針對一維和二維問題的.如崔翔鵬等[12]針對一維問題，在單調(diào)迭代的差分方法上進(jìn)行改進(jìn)，構(gòu)造了一種預(yù)估-校正差分方法，該方法的截?cái)嗾`差為O(τ2+h4)；Karaa[13]提出了求解二維問題的高精度局部一維格式，其截?cái)嗾`差為O(τ2+h4)，且是無條件穩(wěn)定的；趙秉新[14]將一維方程利用指數(shù)變換消去反應(yīng)項(xiàng)，采用高階迎風(fēng)緊致差分公式對對流項(xiàng)進(jìn)行離散，采用高階對稱差分公式對擴(kuò)散項(xiàng)進(jìn)行離散，時(shí)間方向采用四階Runge-Kutta 方法計(jì)算，所提格式的截?cái)嗾`差為O(τ4+h4)，且該格式適用于對流占優(yōu)問題的計(jì)算；楊錄峰等[15]對一維方程消去對流項(xiàng)后，利用四階Padé公式，構(gòu)造了一種無條件穩(wěn)定的差分格式，其截?cái)嗾`差為O(τ2+h4)；楊曉佳等[16]針對一維非穩(wěn)態(tài)方程，采用Taylor 級數(shù)展開式結(jié)合余項(xiàng)修正法，推導(dǎo)出一種截?cái)嗾`差為O(τ2+τh2+h4)的隱式差分格式，但該格式是條件穩(wěn)定的.已有的文獻(xiàn)中關(guān)于三維非穩(wěn)態(tài)對流擴(kuò)散反應(yīng)方程的研究則更加少見.Kaya[17]利用添加一些特殊網(wǎng)格點(diǎn)的處理技巧，結(jié)合算子分裂技術(shù)，在時(shí)間方向用Crank-Nicolson(C-N)格式或向后差分公式離散，構(gòu)造了求解三維對流擴(kuò)散反應(yīng)方程的二階精度的隱格式，且該格式適合于對流和反應(yīng)占優(yōu)問題的計(jì)算；張亞剛[18]針對三維非穩(wěn)態(tài)方程，空間利用四階緊致差分公式，時(shí)間利用四階向后差分公式，提出了一種截?cái)嗾`差為O(τ4+h4)的無條件穩(wěn)定的緊致差分格式，但該格式為五層格式，所以除初始條件外，還需分別構(gòu)造第一、第二、第三時(shí)間層上的格式，才能使時(shí)間向下推進(jìn).綜上所述，關(guān)于三維非穩(wěn)態(tài)對流擴(kuò)散反應(yīng)方程的高階緊致有限差分格式的研究還很少見，有待于進(jìn)一步深入研究.為此，本文將針對含有變對流項(xiàng)和反應(yīng)項(xiàng)系數(shù)的三 維非穩(wěn)態(tài)對流擴(kuò)散反應(yīng)方程，建立一種高精度緊致的有限差分格式.

1 差分格式的建立

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

初始條件為

邊界條件為

然后，利用向前差分離散上式右端的時(shí)間導(dǎo)數(shù)項(xiàng)，利用式(5)計(jì)算空間方向的一階導(dǎo)數(shù)項(xiàng)，利用式(6)離散二階導(dǎo)數(shù)項(xiàng)，可以得到

再將式(10)代入式(8)，得到

最后把式(11)代入式(7)，舍去高階項(xiàng)，得到

類似地，可以給出uy和uz在邊界上的計(jì)算格式.

2 穩(wěn)定性分析

下面利用Fourier 穩(wěn)定性分析法對格式(12)的穩(wěn)定性進(jìn)行分析.在方程(1)中，假設(shè)對流項(xiàng)和反應(yīng)項(xiàng)系數(shù)均為常數(shù)，分別為且f項(xiàng)精確成立，則格式(12)變?yōu)?/p>

在網(wǎng)格點(diǎn)(xi,yj,zk,tn)處，令

將式(16)、(17)分別代入式(15)，進(jìn)行化簡，得到

根據(jù)式(5)有

3 算 法

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

本文將利用以下三個(gè)問題進(jìn)行數(shù)值實(shí)驗(yàn)，驗(yàn)證上述高精度緊致差分格式(12)的穩(wěn)定性和有效性，文中涉及到的誤差及收斂階定義如下:

右端項(xiàng)f(x,y,z,t)和 初邊值條件均由精確解u(x,y,z,t)=etsin(πx)sin(πy)sin(πz)給出.

問題1 是含有非齊次項(xiàng)的三維非穩(wěn)態(tài)問題，本文將對對流項(xiàng)系數(shù)取兩種不同函數(shù)的情況進(jìn)行計(jì)算：①p=tx，q=ty，r=tz； ②p=t(x+y+z)，q=t2xyz，r=sin[π(x+y+z+t)].對于上述兩種情況，本文格式的計(jì)算結(jié)果均符合理論分析.表1 給出了當(dāng) τ=h2，t=0．25時(shí) 的L∞范 數(shù)誤差、L2范數(shù)誤差以及收斂階的比較.可以看出，C-N 格式和BTCS 格式在空間上達(dá)到的精度是二階，而本文格式在空間上達(dá)到的精度是四階.表2 給出了當(dāng)h=0．031 25,t=0．5時(shí) ，不同時(shí)間步長τ 下 的L∞范 數(shù)誤差、L2范數(shù)誤差、收斂階的比較及計(jì)算時(shí)間.從表中可以看出，本文方法的兩種誤差在時(shí)間上均達(dá)到了二階精度，當(dāng)時(shí)間步長較小時(shí)計(jì)算收斂慢，所用的計(jì)算時(shí)間相對較長.表3 給出了當(dāng)t=1， 網(wǎng)格數(shù)取 3 2×32×32， 網(wǎng)格比 λ 分別取0.8，1.6，3.2，6.4 時(shí)，該問題的L∞范 數(shù)誤差和L2范數(shù)誤差.計(jì)算結(jié)果表明，對于給定的所有λ 值，本文格式同C-N 格式和BTCS 格式一樣，仍保持收斂，即本文格式是無條件穩(wěn)定的，這與理論分析的結(jié)果一致.

表1 問題1 當(dāng)τ =h2,t=0．25 時(shí) 的 L∞誤 差和 L2誤差及收斂階Table 1 L∞ errors, L2 errors and convergence rates for τ =h2, t=0．25 in problem 1

表2 問題1 當(dāng)h =0．031 25,t=0．5時(shí) ，本文格式的 L∞誤 差、L 2誤差、收斂階及CPU 時(shí)間Table 2 L∞ errors, L2 errors, convergence rates and CPU time for h =0．031 25, t=0．5 in problem 1

表3 問題1 當(dāng)h =1/32,t=1時(shí) ，對不同網(wǎng)格比λ =τ/h2的 L∞誤 差和 L2誤差Table 3 L∞ errors and L2 errors for h =1/32,t=1 for different mesh ratio λ =τ/h2 values in problem 1

問題2考慮如下Gauss 脈沖問題：

問題2 是常系數(shù)問題，表4 給出了當(dāng) h=0．025時(shí) 對不同的 t,p,q,r,τ，本文格式與幾種交替方向隱格式(ADI)的 L∞范 數(shù)誤差和 L2范數(shù)誤差的比較.其中 Douglas-Gunn ADI格式[1]為二階精度無條件穩(wěn)定的格式，HOC-ADI格 式[3]和 PHOC-ADI格式[4]均為時(shí)間二階、空間四階精度格式，HOC-ADI 格式是條件穩(wěn)定的，而PHOC-ADI 格式是無條件穩(wěn)定的，EHOC-ADI 格式[5]是時(shí)間二階、空間四階精度的無條件穩(wěn)定的指數(shù)型格式，RHOC-ADI 格式[6]是時(shí)間二階、空間四階精度的無條件穩(wěn)定的有理型格式.從表中可以看出，當(dāng)p=q=r=0．8 時(shí) ，本文格式的 L∞范 數(shù)誤差和 L2范數(shù)誤差比 D ouglas-Gunn ADI格式、EHOC-ADI 格式和RHOCADI 格式的都小，計(jì)算結(jié)果更精確；隨著 p,q,r的增大，本文格式的計(jì)算誤差比前四種ADI 格式的計(jì)算誤差都小，但比RHOC-ADI 格式的計(jì)算誤差稍大一些；當(dāng) p,q,r增大到800 和8000 時(shí)，HOC-ADI 格式的計(jì)算結(jié)果發(fā)散.表5 給出了當(dāng)t =0．25,τ=h2,p=q=r=0．8時(shí) ，本文格式的 L∞范 數(shù)誤差和 L2范數(shù)誤差及收斂階，表中結(jié)果表明本文格式在空間上達(dá)到了四階精度.圖1和圖2 給出了問題2 在計(jì)算區(qū)域1 ．2≤x,y≤1．8內(nèi) ，截面 z =1．5上對不同的 t ,p,q,r,τ，不同格式的計(jì)算解和解析解的等值線，其中虛線為解析解，實(shí)線為計(jì)算解.從圖1 可以看出，當(dāng)h=0．025,t=1．25, p=q=r=0．8,τ=6．25×10?3時(shí) , P HOC-ADI 格 式、 E HOC-ADI格 式、 R HOC-ADI格式和本文格式的數(shù)值解均能與解析解吻合得很好，能準(zhǔn)確捕獲以(1.5,1.5)為中心的移動脈沖.從圖2 可以看出，當(dāng)h=0．025,t=1．25×10?4, p=q=r=8 000,τ=6．25×10?7時(shí)，RHOC-ADI 格式的數(shù)值解曲線基本與解析解曲線一致，本文格式的數(shù)值解與解析解之間的吻合度只產(chǎn)生了微弱的偏差，EHOC-ADI 格式的數(shù)值解與解析解之間產(chǎn)生的偏差稍微大一些，而PHOC-ADI 格式的數(shù)值解與解析解之間產(chǎn)生了相當(dāng)大的偏差.圖中結(jié)果表明，同為高精度緊致差分格式，當(dāng) h,α 固 定不變時(shí)，隨著 p,q,r值的不斷增大，相較于EHOC-ADI 格式和PHOC-ADI 格式，本文格式和RHOC-ADI 格式更適合求解此類波的傳播問題.

圖2 問題2 當(dāng)h =0．025, t=1．25×10?4, p=q=r=8 000, τ=6．25×10?7 時(shí) ，在區(qū)域1 ．2≤x,y≤1．8內(nèi)，截面z=1.5 上的等值線圖：(a) PHOC-ADI 格式和精確解；(b) EHOC-ADI 格式和精確解；(c) RHOC-ADI 格式和精確解；(d) 本文格式和精確解Fig.2 Numerical contour lines compared with exact solutions at z=1．5 in region 1．2≤x,y≤1．8 in problem 2 for h=0．025, t=1．25×10?4,p=q=r=8 000, τ=6．25×10?7: (a) PHOC-ADI scheme and exact solutions; (b) EHOC-ADI scheme and exact solutions; (c) RHOC-ADI scheme and exact solutions;(d) the present scheme and exact solutions

表4 問題2 當(dāng)h =0．025 時(shí)，在不同參數(shù)下的 L∞誤 差和 L2誤差Table 4 L ∞ errors and L 2 errors for h =0．025 in problem 2

表5 問題2 當(dāng)t =0．25,τ=h2,p=q=r=0．8時(shí) ，本文格式的 L∞誤 差和L 2誤差及收斂階Table 5 L∞ errors, L2 errors and convergence rates for t =0．25,τ=h2,p=q=r=0．8 in problem 2

圖1 問題2 當(dāng)h =0．025,t=1．25,p=q=r=0．8,τ=6．25×10?3 時(shí)，在區(qū)域 1 ．2≤x,y≤1．8內(nèi) ，截面 z =1．5上的等值線圖：(a) PHOC-ADI 格式和精確解；(b) EHOC-ADI 格式和精確解；(c) RHOC-ADI 格式和精確解;(d) 本文格式和精確解Fig.1 Numerical contour lines compared with exact solutions at z=1．5 in region 1．2≤x,y≤1．8 in problem 2 for h=0．025, t=1．25,p=q=r=0．8,τ=6．25×10?3: (a) PHOC-ADI scheme and exact solutions; (b) EHOC-ADI scheme and exact solutions; (c) RHOC-ADI scheme and exact solutions; (d) the present scheme and exact solutions

問題3考慮如下非線性問題：

其中

精確解為u(x,y,z,t)=?2α[a2e?γαtnxπcos(nxπx)sin(nyπy)sin(nzπz)]d0， 初邊值條件由其給出，其中，d0=1/[a1+a2e?γαtsin(nxπx)sin(nyπy)sin(nzπz)].

問題3 是三維Burgers 方程，在計(jì)算中取nx=ny=nz=3,a1=1,a2=0．1.表6 給出了該問題在不同空間網(wǎng)格數(shù)下，當(dāng)τ =h2,α=0．1,t=0．5時(shí) 的L∞范數(shù)誤差及收斂階的比較.DHOC 格式[20]是條件穩(wěn)定的，且時(shí)間上的精度為二階、空間上的精度為四階.從表中可以看出，DHOC 格式和本文格式在空間上的計(jì)算精度均達(dá)到了四階.表7 給出了該問題當(dāng) τ =0．001,t=0．1時(shí) ，對不同的 α本 文格式與DHOC 格式的L∞范數(shù)誤差的比較，從表中可以看出，本文格式的計(jì)算誤差與DHOC 格式的計(jì)算誤差具有相同的數(shù)量級，但本文格式的計(jì)算誤差更小一些.表8 給出了當(dāng)h=0．062 5,α=0．1,t=2， 網(wǎng)格比 λ分別取1，2，4，6，8 時(shí)，針對該問題本文格式與DHOC 格式的L∞范 數(shù)誤差和L2范數(shù)誤差的比較.計(jì)算結(jié)果表明，隨著λ 不斷增大，DHOC 格式的計(jì)算結(jié)果發(fā)散，即格式是條件穩(wěn)定的，而本文格式的計(jì)算結(jié)果仍保持收斂，即格式是無條件穩(wěn)定的，這與兩種格式各自的理論分析完全一致.充分說明了同為高精度緊致格式，相較于DHOC 格式，本文方法更適合求解此類非線性問題.

表6 問題3 當(dāng)τ =h2, α=0．1, t=0．5時(shí) 的 L∞ 誤 差和 L2誤差及收斂階Table 6 L∞ errors, L2 errors and convergence rates for τ =h2, α=0．1,t=0．5 in problem 3

表7 問題3 當(dāng)τ =0．001, t=0．1時(shí) ，對不同α 的 L∞誤差Table 7 L∞ errors for τ =0．001, t=0．1 with different α values in problem 3

表8 問題3 當(dāng)h =0．062 5, α=0．1, t=2 時(shí) ，對不同網(wǎng)格比λ =τ/h2 的 L∞誤 差和 L2誤差Table 8 L ∞ errors and L 2 errors for h =0．062 5, α=0．1, t=2 with different mesh ratio λ =τ/h2 values in problem 3

5 結(jié) 論

本文針對含有變對流項(xiàng)和反應(yīng)項(xiàng)系數(shù)的三維非穩(wěn)態(tài)對流擴(kuò)散反應(yīng)方程，空間一階導(dǎo)數(shù)采用四階Padé公式計(jì)算，二階導(dǎo)數(shù)利用一階導(dǎo)數(shù)來逼近，時(shí)間采用Taylor 級數(shù)展開和余項(xiàng)修正，構(gòu)造了一種高精度緊致有限差分格式，該格式被證明是無條件穩(wěn)定的，其截?cái)嗾`差為即格式具有空間四階精度，時(shí)間二階精度.由于格式是兩層的，因此可直接利用初始條件啟動，不需要如文獻(xiàn)[18]一樣另外構(gòu)造啟動層的計(jì)算格式，從而使得計(jì)算過程更為簡便.最后進(jìn)行的數(shù)值實(shí)驗(yàn)驗(yàn)證了本文的理論結(jié)果.