一維非定常對流擴散反應方程的高精度緊致差分格式

楊曉佳，田芳

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

?

一維非定常對流擴散反應方程的高精度緊致差分格式

楊曉佳，田芳

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

針對一維非定常對流擴散反應方程，首先推導了一種新的2層高精度緊致差分隱格式，其截斷誤差為O(τ2+τh2+h4)，即當τ=O(h2)時，格式空間具有四階精度；然后采用Fourier分析方法分析了格式的穩(wěn)定性；最后通過數(shù)值算例驗證了本文格式的精確性和可靠性.

對流擴散反應方程；非定常；緊致差分格式；隱式格式；高精度

MSC 2010：O 241.82

對流擴散反應方程是描述黏性流體流動的非線性方程的線性化模型方程，該方程可用于描述許多自然現(xiàn)象，例如水和大氣中污染物質濃度的擴散，海水鹽度﹑溫度的擴散，流體流動和流體傳熱等過程，因此尋求此類方程穩(wěn)定、實用的數(shù)值方法有著重要的理論與實際意義.目前大多數(shù)的文獻報道都針對的是對流擴散方程模型，如葛永斌[1]等人通過時間導數(shù)項采用二階的向后歐拉差分公式，空間二階導數(shù)采用四階緊致差分公式，最終得到了一種求解含源項非定常變系數(shù)對流擴散方程的3層四階緊致隱格式；Ding和Zhang[2]針對一維對流擴散方程，采用半離散的方法和Padé逼近法，推導了一種截斷誤差為O(τ5+h4)的無條件穩(wěn)定的 隱格式.Tian和Yu[3]采用四階緊致指數(shù)差分公式和(2,2)Padé公式分別對空間變量和時間變量進行離散，推導了一種求解一維非定常對流擴散方程的四階精度的無條件穩(wěn)定的緊致差分隱格式；趙飛[4]等人針對一維非定常對流擴散方程，考慮方程在n+1/2時刻的值，時間上利用中心差分格式，空間上采用截斷誤差余項修正的方法及加權平均的思想，最終推導出了非均勻網(wǎng)格上的一種2層高精度全隱式緊致差分格式；Piao等[5]針對一維非定常對流擴散方程，采用3點緊致差分格式和快速單對角隱式Runge-Kutta方法分別對空間和時間進行離散，得到了截斷誤差為O(τ4+h4)的無條件穩(wěn)定的隱格式.也有一些科研工作者針對對流擴散反應方程模型展開研究，如Wo等[6]針對一維非定常對流擴散反應方程，采用Godunov方法，推導出了一種高精度的顯式差分格式.魏劍英[7]提出了一種求解非定常對流擴散反應方程的無條件穩(wěn)定的隱式差分格式，此格式的時間和空間方向均具有二階精度.Liao[8]先將一維非定常對流擴散方程轉化為擴散方程，然后采用四階緊致差分格式，最終推導出了一種截斷誤差為O(τ4+h4)的無條件穩(wěn)定的隱格式，該論文利用此格式計算了對流擴散反應問題.趙秉新[9]針對一維對流擴散反應方程，通過指數(shù)變換將原方程變換為對流擴散方程，對變換后的方程中的對流項和擴散項分別采用高階迎風緊致差分格式和對稱緊致格式進行離散，在時間上采用四階Runge-Kutta方法進行推進，從而得到一種截斷誤差為O(τ4+h4)的緊致差分格式.楊錄峰等[10]針對對流擴散反應方程，采用消除對流項的處理技巧，結合四階Padé公式推導出了一種無條件穩(wěn)定的高精度差分格式.

本文針對一維非定常變系數(shù)對流擴散反應方程，首先采用泰勒級數(shù)展開法和余項修正法，得到一種高精度緊致隱格式.然后采用Fourier分析法分析了格式的穩(wěn)定性.最后通過數(shù)值實驗對本文方法的精確性和可靠性進行驗證.

1 高階緊致差分格式的建立

考慮如下一維變系數(shù)非定常對流擴散反應方程

(1)

其中，a為正常數(shù)，p(x,t)為對流項系數(shù)，q(x,t)為反應項系數(shù)，f(x,t)為源項.以τ表示時間步長，空間取等間距網(wǎng)格，步長用h表示.網(wǎng)格點為(xi,tn)，xi=ih,tn=nτ，i=0，1，…,N,n≥0.

考慮方程(1)在n第時刻的值并對時間和空間導數(shù)項分別采用向前和中心差分離散可得

(2)

).

(3)

為了得到時間二階和空間四階精度的緊致差分格式，需對式(3)中的時間二階導數(shù)項、空間三階和四階導數(shù)項進行處理，為此對原方程(1)兩邊關于和分別求導并整理得

(4)

(5)

(6)

在式(4)、(5)和(6)中，全部變量關于空間的一階和二階導數(shù)項采用中心差分公式離散，未知量關于時間的一階導數(shù)項采用向前差分公式離散，各項系數(shù)p,q及右端項f關于時間的一階導數(shù)項采用二階向后歐拉差分公式離散，并整理得

(7)

(8)

(9)

(10)

將式(7)、(8)和(9)代入式(3)中并整理，可得到方程(1)的高階緊致差分格式

(11)

其中，

(12)

由推導過程可知，差分格式(11)的截斷誤差為O(τ2+τh2+h4)，即當τ=h2時，格式在空間具有四階精度.注意到差分格式(11)對未知量u只涉及到了2層，但是對各項系數(shù)p,q和右端源項f涉及到了3層，由于這些函數(shù)都是已知的，因此可直接計算，不需要另外構造格式啟動.

2 穩(wěn)定性分析

下面對格式(11)的穩(wěn)定性進行分析.

考慮齊次的非定常對流擴散反應方程，即令方程(1)式右端項f=0，變形整理可得

(13)

取p(x,t)/a的上確界為Q，q(x,t)/a的上確界為C，令K=1/a，則(13)式變?yōu)?/p>

(14)

根據(jù)格式式(11)的推導過程，可得式(14)的離散格式為

(15)

其中，

(16)

(17)

則

(18)

其中，

A=(wi-1+wi+1)cosr+wi,B=(wi-1-wi+1)sinr,

F1=(gi-1+gi+1)cosr+gi,F2=(gi-1-gi+1)sinr.

在r-v平面內作出|G|的等值線，圖1-4中分別給出了Pe=1、10、100、1 000，均取R=1時，|G|的等值線，|G|<1的區(qū)域表示格式穩(wěn)定，|G|>1的區(qū)域表示格式不穩(wěn)定.從圖中可以看出，差分格式式(11)對小波數(shù)以及大部分中小波數(shù)是穩(wěn)定的，圖3和圖4表明，隨著Pe數(shù)的增加，差分格式式(11)的穩(wěn)定區(qū)域不再發(fā)生變化.

3 數(shù)值算例

為了驗證本文方法的精確性和可靠性，下面選取以下4個算例進行數(shù)值實驗，算例中的右端項f(x,t)及初邊值條件均由精確解給出， 并與已有文獻計算結果進行比較.計算是用Fortran 77語言進行編程，且在PC機上雙精度制下進行的，離散后在每一個時間層上所得到的代數(shù)方程組可采用追趕法進行求解.

圖1 當R=1,Pe=1時，在r-v平面內|G| 的等值線

圖2 當R=1,Pe=10時，在r-v平面內|G| 的等值線

圖3 當R=1,Pe=100時，在r-v平面內|G| 的等值線

圖4 當R=1,Pe=1 000時，在r-v平面內|G| 的等值線 Fig.4 Contours of the amplification factor |G| in the r-v plane for R=1,Pe=100

其中Ui表示點xi處的精確解，ui表示點xi處的數(shù)值解，Error(N1)和Error(N2)分別表示網(wǎng)格數(shù)為N1和N2時對應的最大絕對誤差.

算例1[1]

其精確解為u(x,t)=e-tsin(x).

表1 問題1當τ=h2,t=0.5時的最大絕對誤差及收斂階

算例2[2,9]

其精確解為u(x,t)=e-tsin(x).

表2 問題2當τ=0.001,t=20時的L2誤差及收斂階

算例3

其精確解為u(x,t)=e5x+t(0.25+0.1π2)sinπx.

表3 問題3當τ=h2,t=0.5時的最大絕對誤差及收斂階

表4 問題3當N=32,t=1時對不同的網(wǎng)格比的最大絕對誤差

算例4

表5 問題4當τ=h2,t=1時的最大絕對誤差及收斂階

算例1是對流擴散方程，表1給出了算例1當τ=h2,t=0.5時的最大絕對誤差和收斂階的比較.從表中的計算結果可以看出，古典隱格式和C-N格式空間只具有二階精度，本文格式(11)空間具有四階精度，這與格式的理論分析相吻合.文獻[1]對此問題空間也具有四階精度，但本文格式的計算結果比文獻[1]的計算結果更加精確，并且文獻[1]中的格式為3層格式，需要另外構造格式來啟動，而本文格式為2層格式，在計算過程中不需要另外構造格式來啟動，這使得本文格式使用時更加簡便.算例2是對流擴散方程，表2給出了算例2當τ=0.001,t=20時的L2誤差和收斂階的比較.從表中的計算結果可以看出，本文格式與文獻[1,2,9]中的格式空間均具有四階精度，但文獻[1]中的格式是一個3層格式，需要另外構造格式來啟動，文獻[2]中的格式在求解過程中需要進行矩陣求逆和乘積等運算，當網(wǎng)格數(shù)很大時，這2種格式的求解效率會明顯下降，本文格式的L2誤差比文獻[1,2,9]中格式的L2誤差小，這說明本文格式具有更高的計算精度.算例3是一個對流擴散反應方程，由于文獻[1,2]中未給出關于此類方程的計算格式，所以在此處未與文獻[1,2]做比較而僅與古典隱格式和C-N格式進行了比較.表3和表4分別給出了算例3的計算結果，其中表3給出了當τ=h2,t=0.5時，算例3的最大絕對誤差和收斂階的比較.從表中的計算結果可以看出，本文格式(11)對此問題空間仍具有四階精度，而古典隱格式和C-N格式空間只具有二階精度，并且本文格式的最大絕對誤差比古典隱格式小1～4個數(shù)量級，比C-N格式小1～3個數(shù)量級.表4給出了當t=1時，網(wǎng)格數(shù)N取32，網(wǎng)格比λ分別取0.4,0.8,1.6,3.2,6.4時，算例3的最大絕對誤差，表中的計算結果與穩(wěn)定性的理論分析是一致的.算例4也是一個對流擴散反應問題，表5給出了當τ=h2,t=1時，算例4的最大絕對誤差和收斂階的比較.從表中的計算結果可以看出，本文格式對此問題空間具有四階精度，而古典隱格式和C-N格式對此問題的空間只具有二階精度，并且本文格式的最大絕對誤差比古典隱格式小2～4個數(shù)量級，比C-N格式小1～4個數(shù)量級，當網(wǎng)格數(shù)較少時，古典隱格式和C-N格式的最大絕對誤差較大.

4 結論

首先建立了一種數(shù)值求解一維非定常對流擴散反應方程的2層全隱式緊致差分格式，格式在每個時間層上只用到了3個網(wǎng)格點，所以差分離散得到的代數(shù)方程組是三對角方程組，故可采用追趕法進行求解.然后，通過Fourier分析方法分析了格式的穩(wěn)定性.最后選取了4個典型數(shù)值算例進行了數(shù)值實驗.數(shù)值實驗結果表明，當τ=O(h2)時，本文格式在空間上可以達到4階精度，與文獻中的計算結果的比較顯示，本文方法總體上具有更高的計算精度.

此外，本文雖然只給出了一維對流擴散反應方程的高精度緊致隱格式，但此方法可以推廣到高維問題的求解.對于高維方程的隱格式的求解，所得到的方程組不再是三對角線型的，所以不能直接采用追趕法，一般需要迭代求解.然而傳統(tǒng)的迭代法收斂速度很慢，為此可采用多重網(wǎng)格方法來加速收斂.文獻[11,12]給出了高維對流擴散方程的高精度隱式差分格式及其多重網(wǎng)格算法，取得了較好的計算效果.從理論上講，將本文方法推廣到高維的對流擴散反應方程后，同樣可以采用多重網(wǎng)格方法進行加速求解.

[1] 葛永斌,田振夫,吳文權.含源項非定常對流擴散方程的高精度緊致隱式差分方法[J].A輯,水動力研究與進展,2006,21(5):619-625. GE Y B,TIAN Z F,WU W Q.A high-order compact implic it difference method for the unsteady convection-diffusion equation with source term[J].Journal of Hydrodynamics,2006,21 (5):619-625.DOI:10.3969/j.issn.1000-4874.2006.05.010.

[2] DING H,ZHANG Y.A new difference scheme with high accuracy and absolute stability for solving convection-diffusion equatons[J].Comput Appl Math,2009,230:600-606.DOI.org/10.1016/j.cam.2008.12.015.

[3] TIAN Z F,YU P X.A high-order exponential scheme for solving 1D unsteady convection-diffusion Equations[J].J Comput Appl Math,2011,235:2477-2491.DOI.org/10.1016/j.cam.2010.11.001.

[4] 趙飛,陳建華,葛永斌.一維非定常對流擴散方程非均勻網(wǎng)格上的高階緊致差分格式[J].西安理工大學學報,2013,29 (4):475-480. ZHAO F,CHEN J H,GE Y B.A High order compact difference scheme on non-uniform grids for the 1D unsteady convection diffusion equation[J].Journal of Xi`an University of Technology,2013,29 (4):475-480.

[5] PIAO X,CHOI H J,KIM S D,et al.A fast singly diagonally implicit Runge-Kutta method for solving 1D unsteady convection-diffusion equations[J].Numer Methods Partial Differential Eq,2013,30:788-812.DOI.org/10.1002/num.21832.

[6] WO S,CHEN B M,WANG J.A high-order Godunov method for one-dimensional convection-diffusion-reaction problems[J].Numer Methods Partial Differential Eq,2000,16:495-512.DOI.org/10.1002/1098-2426(200011)16:6<495:aid-num1>3.0.co;2-s.

[7] 魏劍英.求解對流擴散反應方程的一種隱式差分格式[J].四川理工學院學報(自然科學版),2011,24(5):580-582. WEI J Y.An implicit scheme of the 1D convection-diffusion-reaction equation[J].Journal of Sichuan University of Science & Engineering(Natural Science Edition),2011,24(5):580-582.

[8] LIAO W.A compact high-order finite difference method for unsteady convection-diffusion equation[J].Intern J Comput Meth Eng Sci Mech,2012,13:135-145.DOI.org/10.1080/15502287.2012.660227.

[9] 趙秉新.一種求解一維對流擴散反應方程的高階緊致格式[J].重慶理工大學學報(自然科學版),2012,26(7),100-104. ZHAO B X.A High-order compact difference scheme for solving 1D convection-diffusion-reaction equation[J].Journal of Chongqin University of Technology(Natural Science),2012,26(7),100-104.

[10] 楊錄峰,李春光.一種求解對流擴散反應方程的高精度緊致差分格式[J].寧夏大學學報(自然科學版),2013,34(2):101-109. YANG L F,LI C G.A High-order compact finite difference scheme for solving the convection diffusion reaction equations[J].Journal of Ningxia University(Natural Science Edition),2013,34(2):101-109

[11] 葛永斌,吳文權,田振夫.二維對流擴散方程的高精度全隱式多重網(wǎng)格方法[ J ].計算力學學報,2004,21 (6):678-682. GE Y B,WU W Q,TIAN Z F.Multigrid method based on the high accuracy full implicit scheme of the convection diffusion equation[J].Chinese Journal of Computational Mechanics,2004,21 (6):678-682.DOI:10.3969/j.issn.1007-4708.2004.06.007.

[12] 葛永斌,田振夫,吳文權.三維對流擴散方程的高精度全隱式多重網(wǎng)格方法[ J ].計算力學學報,2007,24 (2):181-186. GE Y B,TIAN Z F,WU W Q.Multigrid method based on high accuracy full implicit scheme of 3D convection diffusion equation[J].Chinese Journal of Computational Mechanics,2007,24 (2):181-186.DOI:10.3969/j.issn.1007-4708.2007.02.010.

(責任編輯：梁俊紅)

High order compact difference scheme for the one-dimensional unsteady convection diffusion reaction equation

YANG Xiaojia,TIAN Fang

(School of Mathematics and Statistis Science,Ningxia University,Yinchuan 750021,China)

A two-level high order compact finite difference implicit scheme is proposed to solve the one-dimensional unsteady convection diffusion reaction equation.The local truncation error of the scheme is O(τ2+τh2+h4),i.e.the scheme is the fourth order accuracy for space when τ=O(h2).Then,Fourier analysis method is used to prove the stability of the scheme.Finally,numerical experiments are conducted to verify the accuracy and the reliability of the present scheme.

convection diffusion reaction equation;unsteady;compact difference scheme;implicit scheme;high accuracy

10.3969/j.issn.1000-1565.2017.01.002

2016-03-28

國家自然科學基金資助項目(11361045;11161036);寧夏大學自然科學基金項目資助(ZR15014);寧夏大學研究生創(chuàng)新項目(GIP201620)

楊曉佳(1988—)，女，寧夏吳忠人，寧夏大學在讀碩士研究生，E-mail：yang_xiaoj@sina.com

田芳(1979—)，女，寧夏中寧人，寧夏大學副教授，主要從事偏微分方程數(shù)值解法的研究. E-mail：tianfsunny@126.com

O241.82

A

1000-1565(2017)01-0005-08