對(duì)流占優(yōu)擴(kuò)散方程的幾種差分解法比較

馮立偉

摘要：對(duì)流占優(yōu)擴(kuò)散方程是描述對(duì)流和擴(kuò)散過(guò)程的模型方程。對(duì)于求解對(duì)流占優(yōu)擴(kuò)散方程的幾種差分格式進(jìn)行分析，并編制了MATLAB程序。使用算例進(jìn)行了數(shù)值試驗(yàn)，差分格式消除了偽數(shù)值振蕩，實(shí)現(xiàn)了對(duì)流占優(yōu)擴(kuò)散方程的求解，并比較了幾種解法的優(yōu)劣。

關(guān)鍵詞： 對(duì)流占優(yōu)擴(kuò)散方程； 有限差分法； 迎風(fēng)； 薩馬爾斯基；指數(shù)

中圖分類(lèi)號(hào)：TP18 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：1009-3044（2015）20-0155-03

Comparing with Numerical Solution of Several Difference Methods for Convection-dominated Diffusion Equation

FENG Li-wei

（Shenyang University of Chemical Technology ， Shenyang 110142， China）

Abstract： Convection-dominated diffusion equation describes convection and diffusion process. It analyzed finite difference methods solving Convection-dominated diffusion equation， and writes MATLAB programs. It makes numerical experiment using example， these experiments reveals the finite difference methods eliminating false numerical oscillation， And it compares advantages and disadvantages of the several finite difference methods.

Key words： convection-dominated diffusion equation； finite difference method； upwind； samarskii； exponential

對(duì)流擴(kuò)散方程描述了在自然界中大量出現(xiàn)的對(duì)流和擴(kuò)散現(xiàn)象，在流體力學(xué)、環(huán)境科學(xué)以及能源開(kāi)發(fā)等諸多領(lǐng)域有著廣泛的應(yīng)用。因此研究對(duì)流擴(kuò)散問(wèn)題的數(shù)值計(jì)算方法就尤為重要。當(dāng)對(duì)流占優(yōu)即小ε值的情形，傳統(tǒng)中心差分格式會(huì)產(chǎn)生數(shù)值振蕩，對(duì)流項(xiàng)的處理上以及對(duì)小ε值的情形處理，成為求解對(duì)流擴(kuò)散方程的關(guān)鍵問(wèn)題。在迎風(fēng)差分格式[1]中對(duì)流項(xiàng)采用一階偏上游的迎風(fēng)差商離散，修正中心差分格式， 薩馬爾斯基格式，指數(shù)差分格式[2，3，4]通過(guò)添加攝動(dòng)項(xiàng)來(lái)消除數(shù)值振蕩。

1 對(duì)流擴(kuò)散方程差分格式

一維對(duì)流占優(yōu)擴(kuò)散方程

[?u?t+a?u?x=ε?2u?x2+fx，t，x，t∈0，L×0，T]

其中[ux，t]為待求量，[a]是對(duì)流系數(shù)，[ε]是擴(kuò)散系數(shù)，[f（x，t）]為源匯。

1.1 差分格式的建立

首先對(duì)[x-t]平面進(jìn)行網(wǎng)格剖分。分別取[h，τ]為空間步長(zhǎng)與時(shí)間步長(zhǎng)，用兩族平行直線[xj=jh ]，[tk=kτ ]將矩形區(qū)域[0， T×0， L]分割成矩形網(wǎng)格。采用中心差商對(duì)導(dǎo)數(shù)進(jìn)行離散可得到中心差分格式，但在[ε]比較小時(shí)擴(kuò)散效應(yīng)比較小主要表現(xiàn)為對(duì)流效應(yīng)，該格式會(huì)出現(xiàn)偽數(shù)值振蕩。消除偽數(shù)值振蕩的一種方法是添加攝動(dòng)項(xiàng)，不同的攝動(dòng)項(xiàng)導(dǎo)致了不同的差分格式，另一種消除偽數(shù)值振蕩，對(duì)流項(xiàng)采用偏心差商的離散，從而得到迎風(fēng)格式。

修正中心差分格式

[Uk+1j-Ukjτ+aUkj+1-Ukj-12h=ε+τa22Ukj+1-2Ukj+Ukj-1h2+fkj]

迎風(fēng)差分格式

[Uk+1j-Ukjτ+aUkj-Ukj-1h=εUkj+1-2Ukj+Ukj-1h2+fkj]

薩馬爾斯基格式

[Uk+1j-Ukjτ+aUkj-Ukj-12h=11+ha/2εεUkj+1-2Ukj+Ukj-1h2+fkj]

指數(shù)差分格式

[ Uk+1j-Ukjτ+aUkj+1-Ukj-12h=-ah21+eah/ε1-eah/εUkj+1-2Ukj+Ukj-1h2+fkj] 其中[k=1，2，…，N-1， j=1，2，…，M-1]。

由Taylor公式容易得出：它們都與一維對(duì)流擴(kuò)散方程相容，其截?cái)嗾`差分別為[O（τ+h2）]，[O（τ+h）]，[O（τ+h2）]和[O（τ+h2）]。

1.2 差分格式的穩(wěn)定性分析

對(duì)流擴(kuò)散方程的中心差分格式的穩(wěn)定性條件是[ετh2≤12，τ≤2εa2]【1】。對(duì)于修正中心差分格式，用[ε+τ2a2]代替中心差分格式穩(wěn)定性條件中的[ε]得到穩(wěn)定性條件[ε+τ2a2τh2≤12，]，[τ≤2a2ε+τ2a2]而后式恒成立，所以只需[ε+τ2a2τh2≤12]【1】。迎風(fēng)差分格式可改寫(xiě)為

[Uk+1j-Ukjτ+aUkj+1-Ukj-12h=ε+ah2Ukj+1-2Ukj+Ukj-1h2+fkj]

用[ε+ah2]代替中心差分格式穩(wěn)定性條件中的[ε]得到穩(wěn)定性條件[ε+ah2τh2≤12，][τ≤2a2ε+ah2] 而前式可推出后式，所以穩(wěn)定性條件為[ε+ah2τh2≤12]。

對(duì)于薩馬爾斯基格式用[ε1+ha/2ε]代替迎風(fēng)差分格式穩(wěn)定性條件中的[ε]得到穩(wěn)定性為[ah2+2ε22ε+ahτh2≤12]。在指數(shù)差分格式中 [1+eah/ε1-eah/ε=σε]，令[σ=ah2εcothah2ε]為擬合因子，指數(shù)格式可改寫(xiě)為

[Uk+1j-Ukjτ+aUkj+1-Ukj-12h=σεUkj+1-2Ukj+Ukj-1h2+fkj]，用[σε]代替中心差分格式穩(wěn)定性條件中的[ε]得到穩(wěn)定性條件[εστh2≤12，τ≤2εσa2]。

1.3差分格式的求解

令[λ=aτh]，[μ=ετh2]，將上述各個(gè)差分格式改寫(xiě)成便于計(jì)算的形式修正中心格式

[Uk+1j=Ukj-λ2Ukj+1-Ukj-1+μ+12λ2Ukj+1-2Ukj+Ukj-1+τfkj] 迎風(fēng)差分格式

[Uk+1j=Ukj-λUkj-Ukj-1+μUkj+1-2Ukj+Ukj-1+τfkj]

薩馬爾斯基格式

[Uk+1j=Ukj-λUkj-Ukj-1+11+ha/2εμUkj+1-2Ukj+Ukj-1+τfkj]

指數(shù)差分格式

[Uk+1j=Ukj-aτ2hUkj+1-Ukj-1-aτ2h21-eah/ε-1Ukj+1-2Ukj+Ukj-1+τfkj]

其中[k=1，2，…，N-1， j=1，2，…，M-1]。因?yàn)楦鱾€(gè)格式均為顯式格式，可以逐層計(jì)算。將差分格式與離散化的初邊值條件結(jié)合，由前一時(shí)間層上的結(jié)果可以計(jì)算下一時(shí)間層上的結(jié)果，從而得到所有的網(wǎng)格點(diǎn)[xj，tk]上的近似解[Ukj]。

1.4 數(shù)值實(shí)驗(yàn)

對(duì)前面的幾種不同差分格式分別編制Matlab程序進(jìn)行求解初邊值問(wèn)題。定義網(wǎng)格點(diǎn)[xj，tk]上的誤差[ekj=u（xj，tk）-Ukj]，使用指標(biāo)[ek1=1Mj=1Mekj]， [ek2=1Mj=1Mekj2]和[ekc=maxjekj]來(lái)度量第k個(gè)時(shí)間層上的總體誤差。

算例1：

[?u?t+?u?x=ε?2u?x2]，相應(yīng)的初邊值條件為

[u（0，t）=exp1-1εx，u（1，t）=exp1ε+1-1εt]，[u（x，0）=expxε]，該問(wèn)題的解析解為

[u（x，t）=exp1εx+1-1εt]

取[ε=0.1]，空間網(wǎng)格步長(zhǎng)為[h=0.1]，時(shí)間網(wǎng)格步長(zhǎng)為[τ=0.01]，計(jì)算結(jié)果繪制成圖像如下：

圖1 t=3時(shí)刻的解

圖2 2-范數(shù)誤差

圖3 1-范數(shù)誤差

圖4 c-范數(shù)誤差

取[ε=0.05]，空間網(wǎng)格步長(zhǎng)為[h=0.1]，時(shí)間網(wǎng)格步長(zhǎng)為[τ=0.01]，計(jì)算結(jié)果繪制成圖像如下：

圖5 t=3時(shí)刻的解

圖6 2-范數(shù)誤差

算例2：

[?u?t+?u?x=ε?2u?x2+πεe-π2tεcos（πxε）]，初邊值條件為

[u（x，0）=sin（πx）ε，u（0，t）=0， u（1，t）=0]，該問(wèn)題的解析解為

[u（x，t）=1εe-π2εtsin（πx）]。

取[ε=0.1]，空間網(wǎng)格步長(zhǎng)[h =0.1]，時(shí)間網(wǎng)格步長(zhǎng)為[τ=0.01]，計(jì)算到3秒時(shí)結(jié)果如圖5-圖8所示。

圖7 t=3時(shí)刻的解

圖8 2-范數(shù)誤差

圖9 1-范數(shù)誤差

圖10 c-范數(shù)誤差

取[ε=0.01]，空間網(wǎng)格步長(zhǎng)為[h=0.1]，時(shí)間網(wǎng)格步長(zhǎng)為[τ=0.01]，計(jì)算結(jié)果繪制成圖像如下

圖11 t=3時(shí)刻的解

圖12 2-范數(shù)誤差

取[ε=0.001]，空間網(wǎng)格步長(zhǎng)為[h=0.1]，時(shí)間網(wǎng)格步長(zhǎng)為[τ=0.01]，計(jì)算結(jié)果繪制成圖像如下。

（下轉(zhuǎn)第215頁(yè)）

（上接第157頁(yè)）

圖13 t=3時(shí)刻的解

圖14 2-范數(shù)誤差

2 結(jié)論

通過(guò)上述算例實(shí)驗(yàn)，可以看到四種差分格式都可抑制數(shù)值振蕩，實(shí)現(xiàn)對(duì)流占優(yōu)擴(kuò)散方程的求解。在文中定義三種范數(shù)度量下，指數(shù)差分格式和修正中心差分格式誤差最小，其次是薩馬爾斯基格式，迎風(fēng)格式效果最差。

參考文獻(xiàn)：

[1] 陸金甫，關(guān)治. 偏微分方程數(shù)值解法[M].北京： 清華大學(xué)出版社，2003：168-170.

[2] 胡健偉，湯懷民. 微分方程數(shù)值方法[M].2版.北京：科學(xué)出版社，2007：131-138.

[3]Morton K W， Mayers D F， Numerical Solution of Partial Differential Equations[M].2th ed. Cambridge UK： Cambridge University Press，2005：19-26.

[4] 劉國(guó)慶， 傅冬生. 非線性奇異攝動(dòng)問(wèn)題一致收斂指數(shù)型擬合差分格式[J]. 高等學(xué)校計(jì)算數(shù)學(xué)學(xué)報(bào)，1998，10（2）： 173-184.

[5] 陳昊. 對(duì)流擴(kuò)散方程差分格式[D]. 成都： 電子科技大學(xué)數(shù)學(xué)科學(xué)學(xué)院，2003：12-15.

[6] 王同科. 一維對(duì)流擴(kuò)散方程CRANK-NICOLSON特征差分格式[J]. 應(yīng)用數(shù)學(xué)，2001，14（4）： 55-60.

[7] 鄭阿奇. MATLAB實(shí)用教程[M].北京： 電子工業(yè)出版社，2005：175-182.

[8] 彭芳麟. 數(shù)學(xué)物理方程的MATLAB解法與可視化[M].北京： 清華大學(xué)出版社，2004：156-160.