趙 飛，蔡志權(quán)，葛永斌

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

0 引言

在處理流體流動和傳熱、溫度擴(kuò)散、海水鹽度及地下水污染物質(zhì)的擴(kuò)散等實際問題中，很多都可歸結(jié)為求解對流擴(kuò)散方程.由于實際問題往往比較復(fù)雜，精確解很難求出，因此研究其精確、穩(wěn)定和高效的數(shù)值方法具有十分重要的意義，有限差分方法是最常用的一種數(shù)值方法［1-13］.針對1維非定常對流擴(kuò)散方程，文獻(xiàn)［1］構(gòu)造了1種對角元嚴(yán)格占優(yōu)的Crank-Nicolson差分格式，文獻(xiàn)［2］提出了 1種Crank-Nicolson特征差分格式，文獻(xiàn)［3］提出了1種攝動有限差分方法，文獻(xiàn)［4］建立了1類指數(shù)型交替顯式方法.遺憾的是這些方法在空間上的整體精度均不超過2階.近年來，已經(jīng)發(fā)展了許多高階緊致差分格式.比如，文獻(xiàn)［5］通過引入1個綜合變換，將對流擴(kuò)散方程變?yōu)榧償U(kuò)散方程，然后利用求解擴(kuò)散方程已有的4階格式，最終得到了1種求解對流擴(kuò)散方程的2層高階緊致隱格式，文獻(xiàn)［6］利用4階精度的3次樣條公式，提出了1種求解非線性對流擴(kuò)散方程的2層緊致隱格式，文獻(xiàn)［7］由定常對流擴(kuò)散方程的4階緊致格式出發(fā)，將時間導(dǎo)數(shù)項作為非齊次項處理，建立了求解1維非定常對流擴(kuò)散方程的1種2層高階緊致隱式差分格式，文獻(xiàn)［8］給出了求解1維含源項非定常對流擴(kuò)散方程的1種3層全隱格式，即可適用于線性又可適用于非線性對流擴(kuò)散問題的求解，文獻(xiàn)［9］基于Hermite插值多項式的構(gòu)造思路，推導(dǎo)出了1維非定常對流擴(kuò)散方程無條件穩(wěn)定的指數(shù)型高精度緊致差分格式.這些格式在空間上均具有4階精度，時間上均具有2階精度.

本文針對1維非定常對流擴(kuò)散方程，結(jié)合已發(fā)展的求解非定常對流擴(kuò)散方程的高階緊致差分格式的基本思想，建立1種無條件穩(wěn)定的有理型高階緊致差分格式，給出格式的穩(wěn)定性分析和證明，并通過數(shù)值算例驗證本文方法的精確性和可靠性.

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

考慮如下1維非定常對流擴(kuò)散方程:

初始條件為u(x，0)= φ(x);邊界條件為u(b1，t)=g0(t)，u(b2，t)=g1(t).(1)式中，a為擴(kuò)散項系數(shù)，p為對流項系數(shù).假設(shè)函數(shù)u(x，t)在給定計算區(qū)域內(nèi)足夠光滑.

將空間求解區(qū)域［b1，b2］進(jìn)行網(wǎng)格剖分為N個子區(qū)間:b1=x0＜x1＜x2＜…＜xN=b2，空間步長定義為h=(b2-b1)/N，用τ表示時間步長.

為了方便格式推導(dǎo)，首先考慮如下1維定常對流擴(kuò)散方程:

a和p為常數(shù)，u，s均為x的函數(shù)，將點ui+1和ui-1在點xi處利用泰勒級數(shù)展開，可得

將(3)式和(4)式相減，整理后可得

將(3)式和(4)式相加，整理后可得

由此定義空間1階導(dǎo)數(shù)和2階導(dǎo)數(shù)的中心差分算子為

將(5)式與(6)式代入(2)式，并利用關(guān)于u的1階導(dǎo)數(shù)和2階導(dǎo)數(shù)的定義，可得

為了得到高精度的差分格式，將(9)式中的3階導(dǎo)數(shù)項和4階導(dǎo)數(shù)項進(jìn)行處理，為此對原方程(2)關(guān)于x求導(dǎo)數(shù)，整理可得

將(10)式代入(11)式，整理可得

將(10)式和(12)式代入(9)式，并利用(7)式和(8)式對u的1階導(dǎo)數(shù)和2階導(dǎo)數(shù)進(jìn)行離散，整理可得

事實上，(13)式即是求解1維定常對流擴(kuò)散方程的多項式型高階緊致差分格式［10-12］.然而，文獻(xiàn)［13］指出了此格式不能用于求解高雷諾數(shù)問題和純對流問題.為此，把(2)式代入(12)式，整理可得

為了得到求解1維定常對流擴(kuò)散方程的有理型高階緊致差分格式，對(9)式做如下處理.首先，將(10)式和(14)式代入(9)式，整理可得

接著，將(5)式和(6)式代入(15)式，并利用(7)式和(8)式對u的1階導(dǎo)數(shù)和2階導(dǎo)數(shù)進(jìn)行離散，整理可得

再將(10)式和(12)式代入(16)式，并利用(5)式和(6)式代替u的1階導(dǎo)數(shù)項和2階導(dǎo)數(shù)項，用(7)式和(8)式對u的1階導(dǎo)數(shù)和2階導(dǎo)數(shù)進(jìn)行離散，整理可得

最后，將(18)式中s的1階導(dǎo)數(shù)項和2階導(dǎo)數(shù)項用中心差分近似，略去4階導(dǎo)數(shù)項，化簡整理可得(2)式的高階緊致差分格式為

用(- ?u?)ti代替(19)式中的si，并且考慮方程在n+1/2時刻的值，則(19)式可以寫為

從推導(dǎo)過程可知，此格式的計算精度為O(h4)+O(τ2)，即此格式在空間上具有4階精度，時間上具有2階精度.

2 穩(wěn)定性分析

下面采用von Neumann方法分析上述格式的穩(wěn)定性.將數(shù)值解表示為Fourier級數(shù)的展開形式，其特征項為

將(23)式代入(22)式得放大因子為

其中，ξ1=1-(4α2h2)sin2(θ 2)，ξ2=(2ατ h2)·sin2(θ 2)，ξ3=(α1h)sin θ，ξ4= ［pτ(2h)］sin θ，sin2(θ 2)∈［0，1］.所以

1)當(dāng)p=0，a ＞0時，由(20)式可知，α=a，α1=0，α2=h212，所以

同理可證，2)當(dāng)p≠0，a＞0時，3)當(dāng)p≠0，a=0時，都有，即針對非定常對流擴(kuò)散方程和非定常純對流方程，本文格式均是無條件穩(wěn)定的.

3 數(shù)值算例

為了驗證本文格式的精確性和可靠性，對3個有精確解的非定常純擴(kuò)散問題、非定常對流擴(kuò)散問題及非定常純對流問題進(jìn)行數(shù)值實驗，并與Crank-Nicolson(C-N)格式和指數(shù)型高階緊致(EHOC)格式［9］的計算結(jié)果進(jìn)行比較.針對非定常純擴(kuò)散問題，本文格式與EHOC格式是同一格式，因此對于該問題，只給出了本文格式的計算結(jié)果.另外，需要說明的是，指數(shù)型高階緊致(EHOC)格式不能被用于求解非定常純對流問題.

例1 非定常純擴(kuò)散問題:

該問題的精確解為 u(x，t)=e-π2tsin(πx).

對于例1，表1給出了當(dāng)τ=h2，T=1時，C-N格式與本文格式在不同空間步長下的最大絕對誤差、2范數(shù)誤差和收斂階.從計算結(jié)果可以看出，本文格式在空間上具有4階精度，而C-N格式只有2階精度，并且隨著空間步長的減小，本文格式的計算誤差要比C-N格式減小得更快.表2給出當(dāng)T=0.2時，針對不同的網(wǎng)格比r(τ=rh2)，采用C-N格式與本文格式計算的時間推進(jìn)步數(shù)、最大絕對誤差和2范數(shù)誤差.可以看出，2種格式均是無條件穩(wěn)定的.當(dāng)取相同的τ和h時，本文格式的計算結(jié)果明顯優(yōu)于C-N格式;當(dāng)取相同的h時，本文格式的計算誤差隨著τ的減小而減小，而C-N格式的計算誤差隨著τ的減小而增大.

表1 當(dāng)τ=h2，T=1時，C-N格式和本文格式的最大絕對誤差、2范數(shù)誤差及收斂階

表2 當(dāng)T=0.2時，C-N格式和本文格式的推進(jìn)步數(shù)、最大絕對誤差及2范數(shù)誤差

例2 非定常對流擴(kuò)散問題:

對于例2，表3給出了當(dāng)a=0.01時，對不同的p、T、τ和 h，C-N 格式、EHOC 格式［9］及本文格式數(shù)值結(jié)果的最大絕對誤差和2范數(shù)誤差.可以看出，對于給定的p、T、τ和h，本文格式的計算誤差要比C-N格式或EHOC格式更小.這充分表明了，針對1維非定常對流擴(kuò)散問題，本文格式比C-N格式和EHOC格式具有更好的計算效果.

表3 C-N格式、EHOC格式及本文格式的最大絕對誤差和2范數(shù)誤差

表3 (續(xù))

例3 非定常純對流問題:

該問題的精確解為 u(x，t)=sin［π(x-t)］.

對于例3，表4給出了當(dāng)τ=h2，T=1時，C-N格式與本文格式在不同空間步長下的最大絕對誤差、2范數(shù)誤差和收斂階.計算結(jié)果表明，本文格式在空間上具有4階精度，而C-N格式只有2階精度，并且在相同空間步長下，本文格式的計算效果明顯優(yōu)于C-N格式.表5給出了當(dāng)T=0.2時，針對不同的網(wǎng)格比r(τ=rh2)，計算的推進(jìn)步數(shù)、最大絕對誤差和2范數(shù)誤差比較.可以看出，2種格式均是無條件穩(wěn)定的.

表4 當(dāng)τ=h2，T=1時，C-N格式與本文格式的最大絕對誤差、2范數(shù)誤差及收斂階

表5 當(dāng)τ=rh2，T=0.2時，C-N格式與本文格式的推進(jìn)步數(shù)、最大絕對誤差及2范數(shù)誤差

4 結(jié)論

本文建立了1種求解1維非定常對流擴(kuò)散方程的有理型高階緊致差分格式，其空間具有4階精度、時間具有2階精度，并且采用von Neumann分析方法證明了該格式針對非定常純擴(kuò)散方程、非定常對流擴(kuò)散方程及非定常純對流方程均是無條件穩(wěn)定的.最后通過數(shù)值實驗，并與C-N格式和EHOC格式的計算結(jié)果進(jìn)行比較，充分驗證了本文方法的精確性和可靠性.

［1］劉揚.對流擴(kuò)散方程的新型Crank-Nicholson差分格式［J］.數(shù)學(xué)雜志，2005，25(4):463-467.

［2］王同科.一維對流擴(kuò)散方程Crank-Nicolson特征差分格式［J］.應(yīng)用數(shù)學(xué)，2001，14(4):55-60.

［3］Gao Zhi.An infinite-order accurate upwind compact difference scheme for the convective-diffusion equati-on［C］∥Proceeding ofAsia Workshop on computational Fluid Dynamics.Chengdu:SichuanUniversity Press，1994:18-24.

［4］蔡國洋，田敏，馮秀芳.一類求解變系數(shù)對流擴(kuò)散方程的指數(shù)型顯式方法［J］.寧夏大學(xué)學(xué)報:自然科學(xué)版，2013，34(1):1-3.

［5］楊志峰，陳國謙.含源項非定常對流擴(kuò)散問題緊致四階差分格式 ［J］.科學(xué)通報，1993，38(2):113-116.

［6］林建國，許維德，陶堯森.含源項非定常非線性對流擴(kuò)散方程的三次樣條四階差分格式［J］.水動力學(xué)研究與進(jìn)展:A 輯，1994，9(2):599-602.

［7］Spotz W F，Carey G F.Extension of high-order compact schemes to time-dependent problems ［J］.Numerical Methods for Partial Differential Equations，2001，17(2):657-672.

［8］葛永斌，田振夫，吳文權(quán).含源項非定常對流擴(kuò)散方程的高精度緊致隱式差分方法［J］.水動力學(xué)研究與進(jìn)展:A 輯，2006，21(5):619-625.

［9］田振夫.一維對流擴(kuò)散方程的四階精度有限差分方法［J］.寧夏大學(xué)學(xué)報:自然科學(xué)版，1995，16(1):32-35.

［10］Karaa S，Zhang Jun.High orderADImethod for solving unsteady convection-diffusion problems［J］.Journal of Computational Physics，2004，198(1):1-9.

［11］MacKinnon R J，Johnson R M.Differential equation based representation of truncation errors for accurate numerical simulation ［J］.International Journal for Numerical Methods in Fluids，1991，13(6):739-757.

［12］Spotz W F，Carey G F.High-order compact scheme for the steady stream-function vorticity equations［J］.International Journal for Numerical Methods in Engineering，1995，38(3):3497-3512.

［13］You D.A High order PadéADImethod for unsteady convection-diffusion equations［J］.Journal of Computational Physics，2006，214(1):1-11.