趙秉新

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

對流擴(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ù)雜的對流擴(kuò)散反應(yīng)現(xiàn)象。描述這些問題的典型方程為流體動力學(xué)方程，而對流擴(kuò)散反應(yīng)方程(CDR，convection-diffusion-reaction)是其中的基本方程。由于無法求得對流擴(kuò)散反應(yīng)方程的解析解，因此，對此類方程的數(shù)值求解具有十分重要的理論和實(shí)際應(yīng)用價值。

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

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

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

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

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

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

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

1.1 空間離散

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

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

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

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

1.2 時間方向離散

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

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

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

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

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

2 數(shù)值驗證

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

2.1 算例1

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

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

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

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

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

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

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

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

2.2 算例2

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

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

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

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

2.3 算例3

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

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

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

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

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

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

4 結(jié)束語

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

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

［6］張昆，楊茉.求解對流擴(kuò)散方程的緊致修正方法［J］.工程熱物理學(xué)報，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.