求解一維對流擴散反應方程的高階緊致格式

趙秉新

(寧夏大學數學計算機學院，銀川 750021)

對流擴散及反應現象廣泛存在于自然界和工程應用中，如在大氣、河流污染中污染物質的擴散，受溫鹽擴散影響的大洋環流，室內空氣的調節，核工業中核反應堆的冷卻及工業生產中的化學氣相沉積中均存在復雜的對流擴散反應現象。描述這些問題的典型方程為流體動力學方程，而對流擴散反應方程(CDR，convection-diffusion-reaction)是其中的基本方程。由于無法求得對流擴散反應方程的解析解，因此，對此類方程的數值求解具有十分重要的理論和實際應用價值。

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

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

1 差分格式的構造

本文考慮一維非定常對流擴散反應方程:

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

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

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

1.1 空間離散

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

其中Si為二階偏導數φxx在i點的逼近值。該格式的離散矩陣為對角占優三對角矩陣，可直接采用追趕法求解。

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

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

1.2 時間方向離散

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

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

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

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

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

2 數值驗證

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

2.1 算例1

在x∈［0，1］上考慮如下對流擴散反應問題:

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

本文對以下3種工況進行計算:

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

2)a=10-4，c=1，r=102，反應占優情形;

3)a=10-4，c=102，r=102，對流及反應占優情形。

對于上述3種工況，數值計算的結果均與精確解吻合得很好。圖1(a)、(b)分別給出在工況1)的對流占優情形下，步長h=1/20和1/40時，問題(11)分別在t=1.5和3.0時刻數值結果與解析解的對比。可以看出，在粗細2種網格下，數值解均與精確解均吻合得很好。

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

2.2 算例2

考慮如下的非線性含源對流擴散反應方程:

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

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

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

2.3 算例3

考慮對流擴散方程:

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

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

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

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

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

4 結束語

本文給出了一種求解對流擴散反應問題的四階緊致迎風差分格式，格式構造過程簡單易懂。數值算例的驗證表明，本文格式能有效控制格式的數值振蕩，適用于對流占優問題的模擬。另外，該格式具有構造簡單，易于編程等特點，且可直接推廣到高維問題。對于定常問題，將?φ/?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］王彩華.一維對流擴散方程的一類新型高精度緊致差分格式［J］.水動力學研究與進展:A輯，2004，19(5):655-663.

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