對流占優(yōu)擴散反應(yīng)方程的定制有限點法

楊紅紅, 秦新強(西安理工大學(xué) 理學(xué)院, 陜西 西安 710054)

對流擴散方程在工程、物理與應(yīng)用科學(xué)等分支中具有重要的地位。用于描述多孔介質(zhì)中多層流體的流動、油藏模擬問題、沿海鹽度擴散及溫度擴散等實際問題。對流占優(yōu)時,方程的解在邊界附近呈現(xiàn)出薄的過渡層[1],采用經(jīng)典的有限差分法和有限元法求解會產(chǎn)生數(shù)值震蕩。

近年來,對于奇異擾動的對流擴散問題[2],Kadalbajoo等[3]提出了在時間方向采用隱式Euler離散,空間方向采用B樣條配方法離散的數(shù)值方法。Clavero等[4]提出了基于隱式Euler對時間離散和Shishkin網(wǎng)格上有限差分進(jìn)行空間離散的方法。Rashidinia等[5]提出時間方向采用差分離散,空間方向采用Sinc-Galerkin離散的方法。此外,最常見的數(shù)值方法還有特征差分法[6]和特征有限元等。上述方法在一定程度上減弱了方程對流占優(yōu)時引起的數(shù)值震蕩,但算法精度不夠高。

本文采用定制有限點方法(TFPM)數(shù)值求解對流占優(yōu)問題,該方法能有效消除數(shù)值偽振蕩,且算法精度高。

TFPM基于局部基函數(shù)“量身”構(gòu)建,最早由Han等[7]提出,用于數(shù)值求解Hemker問題。隨后,Huang等[8]將TFPM用于求解拋物型問題,Tsai等[9]采用TFPM數(shù)值求解Burger方程。

本文第一部分將TFPM與有限差分有效結(jié)合,構(gòu)造了一種新的高精度算法,用于求解對流占優(yōu)擴散方程。第二部分將對流占優(yōu)擴散反應(yīng)方程進(jìn)行指數(shù)變換,對變換后的方程整體采用三角基函數(shù)構(gòu)造TFPM算法格式。第三部分給出了不同離散格式的穩(wěn)定性條件。最后,用數(shù)值算例驗證了本文算法的高效性。

1 對流占優(yōu)擴散方程的TFPM

考慮如下對流占優(yōu)擴散反應(yīng)方程:

(1)

其中Ω=I×(0,T]是有界區(qū)域,ε是擴散系數(shù),0<ε?1,q=0,p>0是區(qū)域Ω上給定的常數(shù)。

取空間步長h>0,時間步長τ>0,節(jié)點記為:

xj=jh,0≤j≤M;tn=nτ,0≤n≤N

(2)

則網(wǎng)格節(jié)點為:

(3)

在圖1模板Ij上構(gòu)造形如:

uxx|x=xj=αj-1uj-1+αjuj+αj+1uj+1

(4)

的TFPM離散格式,其中αj-1,αj,αj+1滿足某種待定關(guān)系,詳見下文。

圖1 二階導(dǎo)數(shù)構(gòu)建模板IjFig.1 Stencil Ij for construct second derivative

假設(shè)u(x)在Ij上可由以下基函數(shù)表出:

V=span{ex/ε,e-x/ε}

(5)

則在Ij上函數(shù)可表示為:

u(x)|Ij=c1je(x-xj)/ε+c2je-(x-xj)/ε

(6)

將其代入αj-1uj-1+αjuj+αj+1uj+1=0,可得:

αj-1(c1je-h/ε+c2jeh/ε)+αj(c1j+c2j)+

αj+1(c1jeh/ε+c2je-h/ε)=0

(7)

即:

(8)

方程組(8)有無窮多個解,若αj已知或給定,則求解方程組得:

(9)

則模板Ij上二階微商的TFPM離散格式:

(10)

對式(1)分別在單元Ω1,Ω2構(gòu)造TFPM顯式格式與隱式格式(見圖2～3)。

圖2 TFPM在Ω1上的節(jié)點單元Fig.2 The reference mesh point at a local cell Ω1 for TFPM

圖3 TFPM在Ω2上的節(jié)點單元Fig.3 The reference mesh point at a local cell Ω2 for TFPM

1.1 顯式TFPM

對式(1)構(gòu)造TFPM顯格式(見圖2),采用TFPM離散擴散項,對流項采用一階迎風(fēng)格式,時間方向采用Euler向前離散,則方程的顯式離散式為:

(11)

初始條件和邊界條件為u(xj,0)=u0(xj),0≤j≤M,u(x,tn)=g0(tn),x∈?I.其中αj-1,αj,αj+1由式(9)給出。

1.2 隱式TFPM

構(gòu)造TFPM隱格式(參見圖3),對式(1)采用TFPM離散擴散項,對流項采用一階迎風(fēng)格式,時間方向采用Crank-Nicolson離散,因:

(12)

故可得TFPM隱式離散格式:

(13)

將(13)寫為矩陣形式:

(14)

系數(shù)矩陣A,B分別為:

(15)

2 對流占優(yōu)擴散反應(yīng)方程的TFPM

考慮如下對流占優(yōu)擴散反應(yīng)方程:

(16)

其中0<ε?1,p≠0,q≠0是區(qū)域Ω上給定的常數(shù)。對方程進(jìn)行變換后整體采用TFPM數(shù)值求解。令u(x,t)=v(x,t)epx/2ε-(p2/4ε+q)t,則式(16)變換為:

(17)

(18)

對擴散方程(17)分別在單元Ω1,Ω2構(gòu)造TFPM顯格式(見圖2)與隱格式(見圖3)。

2.1 顯式TFPM

采用TFPM顯格式(見圖2),令:

(19)

選取基函數(shù)空間:

(20)

則:

(21)

解上述方程組得:

(22)

則可得擴散方程(17)的離散顯格式:

(23)

初始條件和邊界條件為:

(24)

2.2 隱式TFPM

采用TFPM隱格式(見圖3),令:

(25)

(26)

(27)

解方程組(26)與(27)可得:

(28)

(29)

則擴散方程(17)的離散隱格式為:

(30)

將式(30)寫為矩陣形式:

(31)

其中系數(shù)矩陣A,B分別為

(32)

3 穩(wěn)定性分析

定理1TFPM顯格式(11)條件穩(wěn)定。

證明:為了證明TFPM顯格式(11)滿足離散極大值原理,須驗證其系數(shù)矩陣對角占優(yōu),即:

(33)

結(jié)合式(9)可得

(34)

故αj滿足上式時,顯格式(11)穩(wěn)定。

定理2TFPM隱格式(13)無條件穩(wěn)定。

證明:根據(jù)穩(wěn)定性分析的矩陣方法,由隱格式(13)的矩陣格式(14),只需證明C=A-1B為正規(guī)矩陣且其Euclid范數(shù)小于1即可。

CC*=C*C顯然成立,矩陣C的特征值為

j=1,2,…,M

(35)

顯然對所有的j上式的絕對值均小于1,故隱格式(13)無條件穩(wěn)定。

定理3TFPM顯格式(23)條件穩(wěn)定。

(36)

(37)

將上述誤差項代入(23)式,得增長因子

G=eζτ=βj-1e-iph+βj+βj+1eiph=βj+(βj-1+

βj+1)cos(ph)+i(βj+1-βj-1)sin(ph)

(38)

由(22)式與(38)式可得

|G|=|βj+2βj+1cos(ph)|=|1-4βj+1sin2(ph)|

令τ1=ln(cos(h))/ε,當(dāng)τ<τ1時增長因子|G|<1,顯格式(23)穩(wěn)定。

定理4TFPM隱格式(30)無條件穩(wěn)定。

證明:根據(jù)穩(wěn)定性分析的矩陣方法,由隱格式

(30)的矩陣格式(31)得G=A-1B為對稱矩陣且其特征值為

(39)

4 數(shù)值算例

例1考慮如下對流占優(yōu)擴散方程

(40)

取t∈(0,1),τ=0.01,h=1/M,其中0<ε?1是擴散系數(shù),且f,u0,μ1和μ2由精確解確定。方程精確解為u(x,t)=eεt[x+(ex/ε-1)/(e1/ε-1)]。

分別采用4點顯式TFPM格式(11)和6點隱式格式(13)以及文獻(xiàn)[6]中的特征有限差分法(CFDM)數(shù)值求解該方程。表1給出了采用顯格式(11)與CFDM之間最大節(jié)點誤差之間的比較。表2給出了隱格式與文獻(xiàn)[8]所采用TFPM格式之間最大節(jié)點誤差之間的比較。用ETFPM表示顯式TFPM離散格式,ITFPM表示隱式TFPM離散格式,FDM為文獻(xiàn)[10]中的經(jīng)典有限差分法。

表1 ETFPM與CFDM最大節(jié)點誤差比較Tab.1 The comparison of ETFPM and CFDM maximum nodal error

表2 ITFPM與TFPM最大節(jié)點誤差比較Tab.2 The comparison of ITFPM and TFPM maximum nodal error

圖4 本文算法與差分格式(FDM)的比較Fig.4 The comparison of TFPM and FDM

圖5 特征有限差分格式(CFDM)Fig.5 Characteristic finite difference method(CFDM)

從圖4、圖5可以看到,當(dāng)對流占優(yōu)時FDM產(chǎn)生數(shù)值震蕩,而本文TFPM和CFDM可以很好地消除數(shù)值振蕩。又由表1可以看到:本文算法比CFDM的誤差精度高。由表2可以看到:隱式TFPM格式比文獻(xiàn)[8]所采用TFPM格式的誤差精度更高,更高效。由此得到結(jié)論:本文的TFPM與傳統(tǒng)的特征有限差分以及文獻(xiàn)中TFPM[8]相比,在計算誤差方面具有明顯的優(yōu)勢。

例2考慮如下對流占優(yōu)擴散反應(yīng)方程:

(41)

采用4點TFPM顯格式(23)數(shù)值求解該方程。此外,為比較本文算法的優(yōu)劣,與文獻(xiàn)[10]中的經(jīng)典有限差分法(FDM)進(jìn)行數(shù)值比較。表3給出了不同參數(shù)時兩種格式之間的最大節(jié)點誤差,從實驗數(shù)據(jù)可以看到本文所構(gòu)造方法比傳統(tǒng)的有限差分法更有效,誤差精度更高。

表3 TFPM與FDM最大節(jié)點誤差比較Tab.3 The comparison of TFPM and FDM maximum nodal error