求解對流擴(kuò)散反應(yīng)方程的一種高精度緊致差分方法

楊苗苗， 葛永斌

（寧夏大學(xué) 數(shù)學(xué)統(tǒng)計學(xué)院，寧夏 銀川750021）

對流擴(kuò)散反應(yīng)方程作為一種基本的運(yùn)動方程，在很多科學(xué)和工程技術(shù)領(lǐng)域都有著非常廣泛地應(yīng)用，可用于描述流體的流動與傳熱、燃燒與爆炸、核工業(yè)中反應(yīng)堆的冷卻及工業(yè)生產(chǎn)中化學(xué)堆的反應(yīng)等現(xiàn)象.因此對于該類方程精確數(shù)值解的研究具有很高的理論價值和實際作用.

求解該類方程的方法有很多種，最常用的是有限差分法.目前已經(jīng)發(fā)展了很多具有高精度并且緊致的差分格式，如針對對流擴(kuò)散方程：Sun等［1］利用理查森外推法與算子插值法改進(jìn)了四階緊致差分公式，使格式具有六階精度；田振夫［2］在迎風(fēng)變換的基礎(chǔ)上，將對流擴(kuò)散方程轉(zhuǎn)化成含源純擴(kuò)散方程，并利用Hermite法推導(dǎo)出指數(shù)型的差分格式，具有四階精度；梁昌弘等［3］利用3個點上的未知函數(shù)值及四階Padé格式將一階導(dǎo)數(shù)代換的方法，借助顯式緊致格式和隱式緊致格式的思想，得到了一種性能更優(yōu)的四階混合型格式；Tian等［4］建立了一種可用于求解定常對流擴(kuò)散方程的高精度指數(shù)型有限差分格式，適合大梯度問題或邊界層問題的求解；Chen等［5］在差分格式的構(gòu)造中將攝動方法與一階迎風(fēng)格式相結(jié)合，得到了求解二維定常對流擴(kuò)散方程的一種高精度格式；Gupta等［6］采用泰勒級數(shù)展開法得到具有四階精度的有限差分格式；Spotz［7］則利用截斷誤差余項修正的方法給出了一種多項式型的四階緊致差分格式.盡管文獻(xiàn)［6-7］格式的推導(dǎo)方法不一樣，但本質(zhì)上是同一格式；Zhao等［8］將該方法推廣到變系數(shù)的對流擴(kuò)散反應(yīng)方程中，還給出了所得格式的收斂性分析；劉明會［9］基于原方程代入和四階緊致差分公式法，將二階導(dǎo)數(shù)代替后得到了一種新的格式.文獻(xiàn)［8-9］格式盡管推導(dǎo)過程及方法不一樣，但本質(zhì)上仍是同一格式.田振夫［10］用截斷誤差余項修正法，將余項中的三階導(dǎo)數(shù)和四階導(dǎo)數(shù)用一階和二階導(dǎo)數(shù)替換，代入原方程整理便得到四階精度的差分格式，將該方法推廣到二維后，則得文獻(xiàn)［11］中的格式；其他針對對流擴(kuò)散反應(yīng)方程的數(shù)值求解方法還有，魏劍英［12］首先由常系數(shù)變易法得到3個點上的指數(shù)型差分格式，再利用源項在中心點做二階泰勒級數(shù)展開，得到求解常系數(shù)的四階差分格式；田芳等［13］借助常系數(shù)指數(shù)型高精度緊致差分格式，采用殘量修正法得到變系數(shù)對流擴(kuò)散反應(yīng)方程的緊致差分格式；田芳等［14］還基于泰勒級數(shù)展開，給出了一、二階導(dǎo)數(shù)的表達(dá)式，構(gòu)造了非均勻網(wǎng)格上的對流擴(kuò)散反應(yīng)方程的多項式型的差分格式；蘭斌等［15］采用坐標(biāo)變換法將原方程由物理空間的非均勻網(wǎng)格變換為計算空間的均勻網(wǎng)格，再結(jié)合中心差分格式得到了求解對流擴(kuò)散反應(yīng)方程的緊致差分格式.

本文針對一維、二維定常對流擴(kuò)散反應(yīng)方程，基于截斷誤差余項修正的方法推導(dǎo)建立了一種新的緊致差分格式，具有四階精度.盡管本文方法與文獻(xiàn)［7-11］一樣采用到了截斷誤差余項修正法，但本文僅用到了一階導(dǎo)數(shù)的表達(dá)式來修正未知函數(shù)中的三階和四階導(dǎo)數(shù)項，而沒有用到二階導(dǎo)數(shù)項，這樣做的目的是使格式具有更小的耗散誤差.

1 一維問題差分格式的建立

首先，考慮如下一維對流擴(kuò)散反應(yīng)方程：

其中a＞0為擴(kuò)散項系數(shù)，一般是常數(shù)；u（x）是未知函數(shù)，p（x）是對流項系數(shù)，s（x）是反應(yīng)項系數(shù)，f（x）為源項，并且要求u（x）、p（x）、s（x）、f（x）具有充分的光滑性.當(dāng)s（x）≡0時，模型方程（1）為對流擴(kuò)散方程.

將區(qū)間［c，d］分為N等份：x i＝c＋ih，i＝0，1，…，N，定義網(wǎng)格步長h＝（d－c）／N；中心差分算子：

由泰勒級數(shù)展開可得：

將（3）和（4）式代入（1）式得

由（1）式得

對（6）式兩邊關(guān)于x求導(dǎo)可得

則（7）式為文獻(xiàn)［7-10］中處理三階導(dǎo)數(shù)項的方法，即三階導(dǎo)數(shù)由一階和二階導(dǎo)數(shù)來表示.本文為了減小耗散誤差，將（6）式代入（7）式并整理可得

（8）式為本文處理三階導(dǎo)數(shù)項的方法，即u xxx i僅由一階導(dǎo)數(shù)來表示.這是本文格式區(qū)別于文獻(xiàn)［7-10］的地方.同理，對四階導(dǎo)數(shù)的處理如下：

即u xxxx i也僅用到一階導(dǎo)數(shù)u x i，沒有如文獻(xiàn)［7-10］一樣用到二階導(dǎo)數(shù)u xx i.

將（8）和（9）式代入（5）式并舍去 Ο（h4）得

對上式中的對流項和反應(yīng)項系數(shù)的導(dǎo)數(shù)利用中心差分公式（2）離散可得由該過程可知，本文格式的截斷誤差為Ο（h4），即格式（10）在空間上可達(dá)四階精度.對于f（x）的導(dǎo)數(shù)項的計算，可直接代入精確解，也可用中心差分離散，并不會影響格式的整體精度.另外，注意到本文一維格式只用到3個網(wǎng)格點，形成三對角型的代數(shù)方程組，可采用追趕法進(jìn)行求解.

2 誤差分析

下面分別對格式（10）的耗散誤差和色散誤差進(jìn)行分析.

為了簡化計算，將本文格式（10）化為常系數(shù)

則有

同理

將（12）～（14）式代入（11）式得

化簡得

因此，本文差分格式的特征函數(shù)可表示為

其中

另外，為了進(jìn)行對比分析，通過計算可得文獻(xiàn)［3］中MHOC格式的特征函數(shù)為

其中

文獻(xiàn)［8］中FOC格式的特征函數(shù)為

其中

圖1～4分別給出了當(dāng)η＝1，Pe取不同數(shù)值時，MHOC格式［3］、FOC格式［8］和本文格式的耗散性及色散性的誤差分析.不難發(fā)現(xiàn)，對耗散誤差而言，文本格式比MHOC格式和FOC格式的耗散誤差均??；但就色散性而言，本文格式比FOC格式要好，但不及MHOC格式.

圖1 當(dāng)Pe＝10時耗散誤差Fig.1 The dissipation error when Pe＝10

圖2 當(dāng)Pe＝1 000時耗散誤差Fig.2 The dissipation error when Pe＝1 000

圖3 當(dāng)Pe＝10時色散誤差Fig.3 The dispersion error when Pe＝10

圖4 當(dāng)Pe＝1 000時色散誤差Fig.4 The dispersion error when Pe＝1 000

3 二維問題差分格式的建立

考慮如下二維對流擴(kuò)散反應(yīng)方程

其中，a，b＞0為擴(kuò)散項系數(shù)，一般是常數(shù)；u（x，y）是未知函數(shù)，p（x，y）和q（x，y）是對流項系數(shù)，s（x，y）是反應(yīng)項系數(shù)，f（x，y）為源項，并且要求u（x，y）、p（x，y）、q（x，y）、s（x，y）、f（x，y）具有充分的光滑性.

將區(qū)間［c，d］分為N等份：x i＝c＋ih，y j＝c＋jh，i，j＝0，1，…，N，定義網(wǎng)格步長h＝（d－c）／N.

將（15）式可轉(zhuǎn)化成如下2個一維形式的方程：

顯然方程（16）和（17）從形式上看可認(rèn)為是一維的.因此，可直接利用一維對流擴(kuò)散反應(yīng)方程的四階精度差分方法對其進(jìn)行離散，即可得在x方向有

其中

對（16）式中f1（x，y）兩邊關(guān)于x求偏導(dǎo)，可得：

將（20）和（21）式代入（19）式整理得

同理可得在y方向有

其中

于是有

注意（15）式為一種輔助關(guān)系，由（18）和（23）式的離散整理可得

其中

定義中心差分算子：

將（26）式代入（25）式進(jìn)行離散：

整理可得（15）式的離散格式，即為本文格式：

其中

其中，A0，A1，…，A8中的函數(shù)p、q、s及其導(dǎo)數(shù)均在（i，j）點處取值.關(guān)于其一階和二階導(dǎo)數(shù)項的計算，均采用中心差分公式進(jìn)行計算.

由推導(dǎo)過程可知，該格式的截斷誤差為O（h4），即格式（27）在空間上可達(dá)四階精度.另外，注意到本文二維格式只用到9個網(wǎng)格點，故所得格式為緊致格式.

4 數(shù)值算例

為了驗證本文一維和二維格式的精確性和有效性，現(xiàn)將以下有精確解的數(shù)值算例分別與MHOC格式［3］、Chen格式［5］、HOC格式［7］和FOC格式［8］的最大絕對誤差的數(shù)值結(jié)果進(jìn)行比較.其中，將一維和二維問題的最大絕對誤差和收斂階定義如下：

式中的uc和ue分別代表數(shù)值解和精確解，h1和h2代表不同空間步長，Error1和Error2分別對應(yīng)步長為h1和h2時的最大絕對誤差.

問題1［3］

該問題的精確解為

問題2［8］

該問題的精確解為

表2 問題2的最大絕對誤差和收斂階Tab.2 The maximum absolute error and convergence rate for Problem 2

問題3［7］

該問題的精確解為

表3 問題3的最大絕對誤差和收斂階Tab.3 The maximum absolute error and convergence rate for Problem 3

問題4

該問題的精確解為

表4 問題4的最大絕對誤差和收斂階Tab.4 The maximum absolute error and convergence rate for Problem 4

表1～4列出了針對一維問題和二維問題1～4的4種格式在不同h下的最大絕對誤差和收斂階.不難發(fā)現(xiàn)，盡管以上4種格式的精度基本上都達(dá)到了四階，但是本文計算所得的最大絕對誤差要更小些，即本文格式的計算結(jié)果更接近于精確解，也具有更高的計算精度.需要說明的是，對于MHOC格式［3］、Chen格式［5］、HOC格式［7］和FOC格式［8］的計算結(jié)果，采用了Phoebe Solver軟件進(jìn)行計算得到.Phoebe Solver軟件為Web版，網(wǎng)址為www.phoebesolver.com.對于問題4沒有采用Chen格式［5］和HOC格式［7］計算是因為該問題是對流擴(kuò)散反應(yīng)方程，而這2篇文獻(xiàn)中的方程模型為對流擴(kuò)散方程，不包括反應(yīng)項.因此，這2篇文獻(xiàn)所提格式不適用于該問題的求解.

表1 問題1的最大絕對誤差和收斂階Tab.1 The maximum absolute error and convergence rate for Problem 1

5 結(jié)論

本文首先針對方程（1）提出了一種新的緊致差分格式，具有四階精度，即在空間上采取泰勒級數(shù)展開法和截斷誤差余項修正法，先利用一階導(dǎo)數(shù)項來表示由原方程得到的未知函數(shù)的二階導(dǎo)數(shù)項，再表示出誤差余項中的三階和四階導(dǎo)數(shù)項，最終得到了一個三點四階的緊致差分格式，格式對應(yīng)的三對角線型方程組可采用追趕法進(jìn)行求解.接下來，對本文格式的耗散誤差和色散誤差進(jìn)行了分析，與文獻(xiàn)中格式的對比結(jié)果表明本文格式具有更小的耗散誤差.然后將該方法推廣到二維，得到了求解二維對流擴(kuò)散反應(yīng)方程的一種具有四階精度的緊致差分格式.最后給出了數(shù)值實驗，通過與文獻(xiàn)中格式數(shù)值結(jié)果的比較發(fā)現(xiàn)本文格式具有更小的計算誤差和更高的計算精度.

致謝寧夏自治區(qū)重點研發(fā)項目（2018BEE03007）和寧夏大學(xué)研究生創(chuàng)新項目（GIP2019010）對本文給予了資助，謹(jǐn)致謝意.