對(duì)流占優(yōu)擴(kuò)散方程的一種新C—N 緊致差分格式

馮立偉+席偉

摘要：給出了對(duì)流擴(kuò)散方程一種新的C-N緊致差分格式，其截?cái)嗾`差為時(shí)間二階空間四階，且為無(wú)條件穩(wěn)定的。編制了MATLAB程序，數(shù)值試驗(yàn)表明了格式的有效性。

關(guān)鍵詞：對(duì)流擴(kuò)散方程；緊致格式；C-N格式

中圖分類號(hào)：TP311 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：1009-3044（2015）30-0173-04

A New C-N Compact Difference Scheme for Solving Convection-Diffusion Equation

FENG Li-wei1，XI Wei2

（1. ShenYang University of Chemical Technology ， Shenyang 110142， China； 2. ShenYang University of Chemical Technology ， Shenyang 110142， China）

Abstract： It gives a new C-N compact finite difference scheme for solving Convection-diffusion equation. Its error is two order in time and four order in space， and is unconditionally stable， and writes MATLAB programs. Finally， a numerical experiment shows the effectiveness of the scheme.

Key words：convection-diffusion equation； compact scheme；C-N scheme；

對(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ì)算方法就尤為重要[1-3]。田振夫、葛永斌等[4-6]使用Hennite 插值思路給出了求解對(duì)流擴(kuò)散方程的空間四階差分格式。楊志峰和陳國(guó)謙等[7-9]使用綜合變換建立了求解對(duì)流擴(kuò)散方程的一種兩層隱式四階緊致差分格式。肖建英等[10]通過(guò)引入指數(shù)變換造了一種高精度的緊致隱式差分格式，本文采用指數(shù)變換將含源項(xiàng)的對(duì)流擴(kuò)散方程變?yōu)榧償U(kuò)散方程，擴(kuò)散項(xiàng)使用pade逼近的緊致差分離散，時(shí)間層上采用C-N格式，得到了一種不同于[10]文的新緊致差分格式。

1 差分格式的建立

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

[?u?t+a?u?x=ε?2u?x2+fx，t] （ 1 ）

對(duì)[x-t]平面進(jìn)行矩形網(wǎng)格剖分，分別取[h，τ]為空間步長(zhǎng)與時(shí)間步長(zhǎng)， [xj=jh ]，[tk=kτ ]

作指數(shù)變換[u=vea2εx-a24εt]，對(duì)流擴(kuò)散方程變?yōu)閿U(kuò)散方程

[?v?t=ε?2v?x2+F] （ 2 ）

其中[F=ea24εt-a2εxf]

下面對(duì)（ 2 ）式推導(dǎo)差分格式

把[δ2x1+h212δ2xvn+1/2j=?2v?x2n+1/2j+Oh4]代入式（2）得到

[δ2x1+h212δ2xvn+1/2j=1ε?v?t-Fn+1/2j+Oh4]

即

[δ2xvn+1/2j=1ε?v?t-Fn+1/2j+h212εδ2x?v?t-Fn+1/2j+Oh4] （ 3 ）

將右端偏導(dǎo)項(xiàng)離散

[?v?tn+1/2j=δtvn+1/2j+Oτ2]，

[δ2x?v?tn+1/2j=?3v?t?x2n+1/2j+Oh2=?3v?x2?tn+1/2j+Oh2 =δt?2v?x2n+1/2j+Oτ2+h2=δtδ2xvn+1/2j+Oτ2+h2]

將上兩式代入（ 3 ）式得

[δ2xvn+1/2j=1εδtvn+1/2j+h212εδtδ2xvn+1/2j-1εFn+1/2j-h212εδ2xFn+1/2j+Oτ2+h2]

即

[ε2δ2xvn+1j+δ2xvnj=δtvn+1/2j+h212δtδ2xvn+1/2j-12Fn+1j+Fnj-h224δ2xFn+1j+δ2xFnj+Oτ2+h2]令[ετh2=r]，略去高階無(wú)窮小項(xiàng)后得差分格式

[Aj-1vn+1j-1+Ajvn+1j+Aj+1vn+1j+1=Bj-1vnj-1+Bjvnj+Bj+1vnj+1+Cj] （4）

其中

[Aj-1=1-6r，Aj=10+12r，Aj+1=1-6r]

[Bj-1=1+6r，Bj=10-12r，Bj+1=1+6r]

[Cj=τ2Fn+1j+1+10Fn+1j+Fn+1j-1+Fnj+1+10Fnj+Fnj-1]

利用反變換[v=uea24εt-a2εx]得到對(duì)流擴(kuò)散方程的C-N緊致差分格式

[Aj-1Un+1j-1+AjUn+1j+Aj+1Un+1j+1=Bj-1Unj-1+BjUnj+Bj+1Unj+1+Cj] （5）

其中

[Aj-1=1-6rea2τ4ε+ah2ε，Aj=10+12rea2τ4ε，Aj+1=1-6rea2τ4ε-ah2ε]

[Bj-1=1+6reah2ε，Bj=10-12r，Bj+1=1+6re-ah2ε]

[Cj=τ2ea2τ4ε-ah2εfn+1j+1+10ea2τ4εfn+1j+ea2τ4ε+ah2εfn+1j-1+e-ah2εfnj+1+10fnj+eah2εfnj-1]

其中[k=1，2，…，N-1， j=1，2，…，M-1]。

由上面的推導(dǎo)得到

定理1 對(duì)流擴(kuò)散方程（1）的C-N緊致差分格式（5）的截?cái)嗾`差為[O（τ2+h4）]

2 數(shù)值特性

定理2 對(duì)流擴(kuò)散方程（1）的C-N緊致差分格式（5）是無(wú)條件穩(wěn)定的

證明：差分格式（4）對(duì)應(yīng)的齊次差分格式為

[Aj-1vn+1j-1+Ajvn+1j+Aj+1vn+1j+1=Bj-1vnj-1+Bjvnj+Bj+1vnj+1]

以[vnj=Vnexpikxj]代入上式消去公因子[expikxj]得

[Vn+11-6reikh+e-ikh+10+12r=1+6reikh+e-ikh+10-12rVn]

得到增長(zhǎng)因子為

[G=10-12r+1+6r2coskh10+12r+1-6r2coskh=10+2coskh-12r1-coskh10+2coskh+12r1-coskh]

[G≤1]，差分格式（4）是無(wú)條件穩(wěn)定的，所以差分格式（5）也是無(wú)條件穩(wěn)定的

3 差分格式的計(jì)算

將差分格式（5）改寫成便于計(jì)算的形式，并代入邊界條件后，得到

[A0A1A-1A0A1???A-1A0A1A-1A0 Un+11Un+12?Un+1M-2Un+1M-1=B0B1B-1B0B1???B-1B0B1B-1B0 Un1Un2?UnM-2UnM-1+Cn1+B-1Un0-A-1Un+10Cn2?CnM-2CnM-1+B1UnM-A1Un+1M][Cn1Cn2?CnM-2CnM-1=τ2C0C1C-1C0C1???C-1C0C1C-1C0 fn+11fn+12?fn+1M-2fn+1M-1+D0D1D-1D0D1???D-1D0D1D-1D0 fn1fn2?fnM-2fnM-1+D-1fn0+C-1fn+100?0D1fnM+C1fn+1M]因?yàn)椴罘指袷较禂?shù)矩陣為對(duì)角占優(yōu)陣，所以差分方程的解存在且唯一。

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

使用前面的差分格式（5）求解下面的對(duì)流擴(kuò)散方程的初邊值問(wèn)題。定義網(wǎng)格點(diǎn)[xj，tn]上的誤差[enj=u（xj，tn）-Unj]，使用指標(biāo) [en2=1Mj=1Menj2]來(lái)度量第[n]個(gè)時(shí)間層上的總體誤差。

算例1：

[?u?t+?u?x=ε?2u?x2，x，t∈-1，1×0，1]，

該問(wèn)題的準(zhǔn)確解為[u（x，t）=exp1εx+1-1εt]

相應(yīng)的初邊值條件為

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

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

例2，考慮以下一維常系數(shù)模型問(wèn)題

[?u?t+a?u?x=ε?2u?x2+fx，t∈0，1×0，Tu（x，0）=sinπx 0≤x≤1u（0，t）=0，u（1，t）=0 0≤t≤1]

解析解設(shè)為[u（x，t）=e-tsinπx]，右端函數(shù)由解析解確定

取[a=1，ε=0.1]，空間和時(shí)間網(wǎng)格步長(zhǎng)[h =0.1，τ=0.01]，計(jì)算到[t=2]。

5 結(jié)論

通過(guò)上述兩個(gè)數(shù)值實(shí)驗(yàn)，可看出本格式可實(shí)現(xiàn)對(duì)流擴(kuò)散方程的求解。

參考文獻(xiàn)：

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

[2] 王文洽. 變系數(shù)對(duì)流擴(kuò)散方程的交替分段Crank-Nicolson方法[J]. 應(yīng)用數(shù)學(xué)和力學(xué)，2004，24（1）：29-36.

[3] 由同順. 對(duì)流擴(kuò)散方程的本質(zhì)非振蕩特征差分方法[J]. 應(yīng)用數(shù)學(xué)， 2000，13（4）：89-94.

[4] 田振夫.一維對(duì)流擴(kuò)散方程的四階精度有限差分法[J].寧夏大學(xué)學(xué)報(bào)：自然科學(xué)版，1995，16 （1）： 30-35.

[5] 葛永斌，田振夫，詹詠，等.求解擴(kuò)散方程的一種高精度緊致隱式差分方法[J].上海理工大學(xué)學(xué)報(bào)，2005， 27（2）：107-110.

[6] 葛永斌，田振夫，吳文權(quán).含源項(xiàng)非定常對(duì)流擴(kuò)散方程的高精度緊致隱式差分方法[J].水動(dòng)力學(xué)研究與進(jìn)展（輯），2006， 121（5）：619-625.

[7] 楊志峰，陳國(guó)謙.含源擴(kuò)散方程的一種高精度差分方法[J].北京師范大學(xué)學(xué)報(bào)：自然科學(xué)版，1992，28（3）： 315-316

[8] 楊志峰，陳國(guó)謙.含源匯非定常對(duì)流擴(kuò)散問(wèn)題緊致四階差分格式[J].科學(xué)通報(bào)，1993，38 （2）：113-116.

[9] 王煊，楊志峰. 基于非均勻網(wǎng)格求解非線性對(duì)流擴(kuò)散問(wèn)題的一種高精度差分格式[J].北京師范大學(xué)學(xué)報(bào)，自然科學(xué)版，2003，39（1）：131-137.

[10] 肖建英 劉小華 李永濤.非定常對(duì)流擴(kuò)散方程得高階差分格式[J]. 西南石油大學(xué)：自然科學(xué)版，2012，34（3）：145-149.