求解對流擴散反應方程的一種高精度緊致差分方法

楊苗苗， 葛永斌

（寧夏大學 數學統計學院，寧夏 銀川750021）

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

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

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

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

首先，考慮如下一維對流擴散反應方程：

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

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

由泰勒級數展開可得：

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

由（1）式得

對（6）式兩邊關于x求導可得

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

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

即u xxxx i也僅用到一階導數u x i，沒有如文獻［7-10］一樣用到二階導數u xx i.

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

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

2 誤差分析

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

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

則有

同理

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

化簡得

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

其中

另外，為了進行對比分析，通過計算可得文獻［3］中MHOC格式的特征函數為

其中

文獻［8］中FOC格式的特征函數為

其中

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

圖1 當Pe＝10時耗散誤差Fig.1 The dissipation error when Pe＝10

圖2 當Pe＝1 000時耗散誤差Fig.2 The dissipation error when Pe＝1 000

圖3 當Pe＝10時色散誤差Fig.3 The dispersion error when Pe＝10

圖4 當Pe＝1 000時色散誤差Fig.4 The dispersion error when Pe＝1 000

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

考慮如下二維對流擴散反應方程

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

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

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

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

其中

對（16）式中f1（x，y）兩邊關于x求偏導，可得：

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

同理可得在y方向有

其中

于是有

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

其中

定義中心差分算子：

將（26）式代入（25）式進行離散：

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

其中

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

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

4 數值算例

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

式中的uc和ue分別代表數值解和精確解，h1和h2代表不同空間步長，Error1和Error2分別對應步長為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下的最大絕對誤差和收斂階.不難發現，盡管以上4種格式的精度基本上都達到了四階，但是本文計算所得的最大絕對誤差要更小些，即本文格式的計算結果更接近于精確解，也具有更高的計算精度.需要說明的是，對于MHOC格式［3］、Chen格式［5］、HOC格式［7］和FOC格式［8］的計算結果，采用了Phoebe Solver軟件進行計算得到.Phoebe Solver軟件為Web版，網址為www.phoebesolver.com.對于問題4沒有采用Chen格式［5］和HOC格式［7］計算是因為該問題是對流擴散反應方程，而這2篇文獻中的方程模型為對流擴散方程，不包括反應項.因此，這2篇文獻所提格式不適用于該問題的求解.

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

5 結論

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

致謝寧夏自治區重點研發項目（2018BEE03007）和寧夏大學研究生創新項目（GIP2019010）對本文給予了資助，謹致謝意.