求解一維對(duì)流擴(kuò)散反應(yīng)方程的高階緊致格式

趙秉新

(寧夏大學(xué)數(shù)學(xué)計(jì)算機(jī)學(xué)院，銀川 750021)

對(duì)流擴(kuò)散及反應(yīng)現(xiàn)象廣泛存在于自然界和工程應(yīng)用中，如在大氣、河流污染中污染物質(zhì)的擴(kuò)散，受溫鹽擴(kuò)散影響的大洋環(huán)流，室內(nèi)空氣的調(diào)節(jié)，核工業(yè)中核反應(yīng)堆的冷卻及工業(yè)生產(chǎn)中的化學(xué)氣相沉積中均存在復(fù)雜的對(duì)流擴(kuò)散反應(yīng)現(xiàn)象。描述這些問題的典型方程為流體動(dòng)力學(xué)方程，而對(duì)流擴(kuò)散反應(yīng)方程(CDR，convection-diffusion-reaction)是其中的基本方程。由于無(wú)法求得對(duì)流擴(kuò)散反應(yīng)方程的解析解，因此，對(duì)此類方程的數(shù)值求解具有十分重要的理論和實(shí)際應(yīng)用價(jià)值。

目前，求解CDR方程的常用數(shù)值方法為有限差分法、有限元法及有限體積法等，其中有限差分法作為一種發(fā)展成熟且易于應(yīng)用的方法受到了廣泛的關(guān)注。數(shù)值求解對(duì)流擴(kuò)散反應(yīng)型方程不同程度地存在著產(chǎn)生數(shù)值擴(kuò)散和振蕩解的可能，如經(jīng)典二階中心格式會(huì)產(chǎn)生非物理振蕩，導(dǎo)致計(jì)算結(jié)果不可靠［1-2］等。減小數(shù)值振蕩的有效方法是構(gòu)造穩(wěn)定的、且能有效反映對(duì)流項(xiàng)的迎風(fēng)效應(yīng)的數(shù)值格式，要求差分格式的設(shè)計(jì)要保持流體流動(dòng)的傳輸特征。常用的一階迎風(fēng)格式是一種無(wú)條件穩(wěn)定的格式，但卻只有一階精度，且存在無(wú)法有效捕捉微小物理量等缺陷。因此，要想對(duì)問題進(jìn)行較成功的數(shù)值模擬，同時(shí)滿足實(shí)際問題中高精度的要求，就必須對(duì)計(jì)算格式的求解效率及精度提出更高要求。高階緊致迎風(fēng)差分格式通過使用較少的網(wǎng)格基架點(diǎn)便能達(dá)到提高精度的要求，保持了迎風(fēng)格式的優(yōu)勢(shì)，具有邊界無(wú)需特殊處理以及能達(dá)到與譜方法相近的分辨率等優(yōu)點(diǎn)，是數(shù)值計(jì)算方法研究的一個(gè)重要方向。Noye等［3］利用加權(quán)修正技術(shù)，提出了一種三階半隱格式，但其是條件穩(wěn)定的。陳國(guó)謙等［4］通過對(duì)對(duì)流系數(shù)和源項(xiàng)作2階修正，消除截?cái)嗾`差中的低階項(xiàng)，給出了對(duì)流擴(kuò)散方程的4階指數(shù)型格式。王彩華［5］通過在低精度格式基礎(chǔ)上進(jìn)行簡(jiǎn)單修正得到了針對(duì)一維無(wú)源對(duì)流擴(kuò)散方程的4種高精度差分格式。通過在低階格式中的源項(xiàng)中引入緊致修正項(xiàng)，楊茉等［6］提出了一種緊致修正方法。Ding等［7］提出了一種O(h4+τ4)階無(wú)條件穩(wěn)定的格式，但適用于對(duì)流占優(yōu)問題的求解［8］，且存在求解效率不高的問題。

本文通過將原方程作簡(jiǎn)單轉(zhuǎn)換，對(duì)得到的對(duì)流擴(kuò)散方程中的對(duì)流項(xiàng)采用四階緊致迎風(fēng)格式離散，擴(kuò)散項(xiàng)采用四階Padé格式離散;之后，對(duì)于空間半離散格式，時(shí)間方向采用四階龍格庫(kù)塔方法計(jì)算，從而整體達(dá)到了O(h4+τ4)階精度。經(jīng)數(shù)值驗(yàn)證，該格式具有良好性能。該格式構(gòu)造方法簡(jiǎn)單，容易理解，且易于推廣到高維情形。

1 差分格式的構(gòu)造

本文考慮一維非定常對(duì)流擴(kuò)散反應(yīng)方程:

其中:Ω =［0，l］×［0，T］;φ(x)為充分光滑的函數(shù);常擴(kuò)散系數(shù)a＞0;c(x，t)表示相速度;r為反應(yīng)速度;g1，g2為常數(shù)。

以τ=Δt為時(shí)間步長(zhǎng)，空間取均勻等間距網(wǎng)格，步長(zhǎng)為 h= Δx=1/m，網(wǎng)格點(diǎn)為(xi，tn)，其中xi=ih，i=0，1，…，m，tn=nτ，n≥0。對(duì)方程(1)做如下變換:令u=exp(-rt)φ，代入方程(1)中消去反應(yīng)項(xiàng)，對(duì)定解條件(2)和(3)作相應(yīng)變換，得

方程(4)為未知量φ的對(duì)流擴(kuò)散方程，記Q(x，t)=e-rtf(x，t)為源項(xiàng)。以下給出式(4)～(6)的差分格式，回代求解方程(1)～(3)。

1.1 空間離散

將方程(4)中的擴(kuò)散項(xiàng)φxx采用四階精度Padé格式離散:

其中Si為二階偏導(dǎo)數(shù)φxx在i點(diǎn)的逼近值。該格式的離散矩陣為對(duì)角占優(yōu)三對(duì)角矩陣，可直接采用追趕法求解。

對(duì)對(duì)流項(xiàng)進(jìn)行離散時(shí)，考慮迎風(fēng)效應(yīng)而采用四階緊致迎風(fēng)格式［13］離散。例如對(duì)流項(xiàng)c(x，t)φx可表示成

其中Si利用式(7)計(jì)算。該格式基于3個(gè)網(wǎng)格基架點(diǎn)，具有四階精度。

1.2 時(shí)間方向離散

Runge-Kutta方法最初用于常微分方程的數(shù)值求解。當(dāng)把時(shí)間處理成獨(dú)立變量時(shí)，可將式(4)進(jìn)行空間離散后的半離散方程視為只依賴于時(shí)間的常微分方程，此時(shí)，可采用Runge-Kutta方法進(jìn)行時(shí)間推進(jìn)求解。本研究將式(4)的半離散格式寫為

其中Lh為空間半離散算子，且

Qi為源項(xiàng)Q(x，t)在i點(diǎn)的離散值。

對(duì)半離散格式(8)，采用四步四階Runge-Kutta方法求解，具體為:

其中Lh(φn)為φ在n時(shí)間層時(shí)式(9)的值。

2 數(shù)值驗(yàn)證

為驗(yàn)證格式的性能，考察3個(gè)有精確解的算例，并與已有格式的結(jié)果進(jìn)行對(duì)比。

2.1 算例1

在x∈［0，1］上考慮如下對(duì)流擴(kuò)散反應(yīng)問題:

其中源項(xiàng) f(x，t)=e-t[(r-1)x2+2(cx-a)]，其解析解為 u(x，t)=x2e-t。

本文對(duì)以下3種工況進(jìn)行計(jì)算:

1)a=10-4，c=102，r=1，對(duì)流占優(yōu)情形;

2)a=10-4，c=1，r=102，反應(yīng)占優(yōu)情形;

3)a=10-4，c=102，r=102，對(duì)流及反應(yīng)占優(yōu)情形。

對(duì)于上述3種工況，數(shù)值計(jì)算的結(jié)果均與精確解吻合得很好。圖1(a)、(b)分別給出在工況1)的對(duì)流占優(yōu)情形下，步長(zhǎng)h=1/20和1/40時(shí)，問題(11)分別在t=1.5和3.0時(shí)刻數(shù)值結(jié)果與解析解的對(duì)比?？梢钥闯?，在粗細(xì)2種網(wǎng)格下，數(shù)值解均與精確解均吻合得很好。

圖1 工況(1)下數(shù)值解與精確解的對(duì)比

2.2 算例2

考慮如下的非線性含源對(duì)流擴(kuò)散反應(yīng)方程:

精確解取為u(x，t)=tanh( Re(1-2x)/4)et，源項(xiàng)f(x，t)可由原方程確定。

圖2給出了Re=1 000，T=0.02時(shí)的精確解及CDS、SCD4和本文格式在中心點(diǎn)附近的解(步長(zhǎng)取 h=1/320，τ=10-5)?？梢?，CDS格式和SCD4格式在該網(wǎng)格尺度下均產(chǎn)生明顯的數(shù)值振蕩，而本文格式能夠獲得很好的數(shù)值解，體現(xiàn)了其求解非線性對(duì)流擴(kuò)散反應(yīng)方程時(shí)的良好性能。

圖2 算例2不同格式在間斷位置附近的解

2.3 算例3

考慮對(duì)流擴(kuò)散方程:

精確解為 u(x，0)=e5x-t(0.25+0.01π2)sin(πx)。

時(shí)間步長(zhǎng)取τ=0.001，表1給出了T=20時(shí)，本文格式(Present)與Crank-Nicolson格式及Ding格式［7］在不同空間網(wǎng)格步長(zhǎng)下的L2誤差和收斂階(Rate)的對(duì)比情況。收斂階由下式定義:

其中e1、e2分別表示網(wǎng)格步長(zhǎng)為 h1、h2時(shí)的 L2誤差。

通過對(duì)比表明:本文格式與Ding的格式［7］均具有四階空間精度，高于Crank-Nicolson格式的二階精度;但Ding格式均存在求解效率不高的問題，即在求解過程中需要進(jìn)行矩陣求逆和乘積等運(yùn)算，當(dāng)網(wǎng)格數(shù)很大時(shí)，求解效率會(huì)明顯降低，特別是在求解非線性問題時(shí)，由于每個(gè)迭代步的迭代矩陣均不相同，這樣每個(gè)迭代步都必須進(jìn)行矩陣求逆運(yùn)算，在高維或多網(wǎng)格節(jié)點(diǎn)下，計(jì)算效率會(huì)明顯降低。

表1 算例1中不同網(wǎng)格尺度下L2誤差和收斂階對(duì)比

4 結(jié)束語(yǔ)

本文給出了一種求解對(duì)流擴(kuò)散反應(yīng)問題的四階緊致迎風(fēng)差分格式，格式構(gòu)造過程簡(jiǎn)單易懂。數(shù)值算例的驗(yàn)證表明，本文格式能有效控制格式的數(shù)值振蕩，適用于對(duì)流占優(yōu)問題的模擬。另外，該格式具有構(gòu)造簡(jiǎn)單，易于編程等特點(diǎn)，且可直接推廣到高維問題。對(duì)于定常問題，將?φ/?t視為人工時(shí)間項(xiàng)，計(jì)算程序無(wú)需太多修改即可直接求解。

［1］Brandt A，Yavneh I.Inadequacy of first-order upwind difference schemes for some recirculating flows［J］.Journal of Computational Physics，1991，93(1):128-143.

［2］Zhang J.Accelerated multigrid high accuracy solution of the convection-diffusion equation with high Reynolds number［J］.Numerical Methods for Partial Differential E-quations，1997，13(1):77-92.

［3］Noye B J，Tan H H.A third-order semi-implicit finite difference method for solving the one-dimensional convection-diffusion equation［J］.International Journal for Numerical Methods in Engineering，1988，26(7):1615-1629.

［4］Chen G Q，Gao Z，Yang Z F.A perturbational h4exponential finite difference scheme for the convective diffusion equation［J］.Journal of Computational Physics，1993，104(1):129-139.

［5］王彩華.一維對(duì)流擴(kuò)散方程的一類新型高精度緊致差分格式［J］.水動(dòng)力學(xué)研究與進(jìn)展:A輯，2004，19(5):655-663.

［6］張昆，楊茉.求解對(duì)流擴(kuò)散方程的緊致修正方法［J］.工程熱物理學(xué)報(bào)，2011，32(3):459-461.

［7］Ding H，Zhang Y.A new difference scheme with high accuracy and absolute stability for solving convection-diffusion equations［J］.Journal of Computational and Applied Mathematics，2009，230(2):600-606.

［8］Tian Z F，Yu P X.A high-order exponential scheme for solving 1D unsteady convection-diffusion equations［J］.Journal of Computational and Applied Mathematics，2011，235(8):2477-2491.

［9］Tian Z F，Liang X，Yu P X.A higher order compact finite difference algorithm for solving the incompressible Navier-Stokes equations［J］.International Journal for Numerical Methods in Engineering，2011，88(6):511-532.