張麗劍，羅躍生，張文平

（1.哈爾濱工程大學動力與能源工程學院，黑龍江哈爾濱150001；2.哈爾濱工程大學理學院，黑龍江哈爾濱150001）

變限積分的有限體積法解決對流擴散方程

張麗劍1，羅躍生2，張文平1

（1.哈爾濱工程大學動力與能源工程學院，黑龍江哈爾濱150001；2.哈爾濱工程大學理學院，黑龍江哈爾濱150001）

針對一維對流擴散方程，基于變限積分的有限體積法，提出一種高精度全離散方法。該方法在控制體內對方程進行變限積分，引入變限因子，然后分別對上下限再進行積分，從而將微分方程轉化為積分方程，最后運用插值的方法對目標函數進行近似代替。該方法提高了計算精度，結果得到一維的收斂精度。采用Fourier分析法證明了格式為條件穩(wěn)定。最后給出了非穩(wěn)態(tài)和穩(wěn)態(tài)2種情況下的數值算例，驗證了所提出的格式具有高精度和易于編程計算的優(yōu)點。

變限積分；高精度；對流-擴散方程；Fourier分析法；非穩(wěn)態(tài)；一維

對流擴散方程在許多物理系統，特別是流體力學與傳熱學中有廣泛的應用［1?3］。它描述的是諸如質量、熱量、能量、渦度等的對流擴散的物理量的運動機理［4?11］。對于初邊值條件一維對流擴散方程解析解的研究，幾乎沒有取得任何進展［1，4］。因此，很多的研究都致力于尋找對流擴散方程穩(wěn)定的和相對精確的數值解［12?13，23?26］。近年，高階緊致差分格式被用來求解方程的數值解［14?16，27］。格式最高具有四階空間和時間精度［17?19，27］。本文提出了一種高精度全離散格式，其主要思想是將微分方程轉化為積分方程來求解。

1 格式的構造

一維非穩(wěn)態(tài)對流擴散方程：

式中：（0，T］為時間域，φ、g1和g2均為已知的充分光滑函數，k和r分別為相速度和粘性系數，均為已知的常系數。

對于問題（1）～（3），本文首先將區(qū)域劃分成M×N的網格，其中網格xi＝ih，i＝0，1，2，…，M，h＝l／M；tn＝nτ，n＝0，1，2，…，N，τ＝T／N。記u( ih，kτ)＝uni表示函數u在網格(ih，kτ)上的值。

在控制區(qū)域Ω＝[y1，y2]×[tn，tn+1]內，首先運用有限體積法思想對方程（1）兩邊作變限積分，其中y1，y2∈(0，l)，最后得到化簡式：

然后在上下限的控制區(qū)域Ω′＝[xi-ε2，xi]× [xi，xi+ε1]內，分別對方程（4）兩邊進行積分：

式中：ε1、ε2分別為控制上下限積分區(qū)域大小的變限因子。式（5）中，沿x方向上的積分區(qū)域如圖1所示。

圖1 沿x方向上的積分區(qū)域Fig.1 The domain of integration of the proposed for?mat about x direction

通過上述積分，將微分方程（1）轉化為積分方程（5），于是用二元六點拉格朗日插值法，對u x，t()進行插值。整理并舍去小項，可得方程（1）的全離散格式：

則式（6）、（7）即為問題（1）～（3）的全離散格式。

2 精度估計

由于在對方程（1）進行離散時，一共對x進行了3次積分，對t進行了1次積分，則等式兩邊應除以h3τ，從而可得離散格式（6）的截斷誤差為O h4／τ+h3+hτ+τ2／h+τ3／h2

( )，顯然，該截斷誤差大于二階。

3 穩(wěn)定性分析

假設邊界條件精確滿足，采用Fourier分析法對全離散格式進行穩(wěn)定性分析，用μ表示計算u時產生的誤差，代入離散格式（6），則滿足齊次方程：

顯然可得G′1

2≤G′22，即G2≤1，由此當ε1≥ε2時可得全離散格式穩(wěn)定。

4 算法分析

基于變限積分的有限體積法是在有限體積法的基礎上提出的新算法，該方法結合了有限元和有限差分法的思想，同時又繼承了有限體積法的優(yōu)點。

該算法首先對整體計算域進行控制區(qū)域的劃分，然后對目標方程在局部控制區(qū)域內進行變限積分。將微分方程完全轉化為積分方程，使其消除偏導數影響。最后，選擇關于網格點的插值多項式對積分方程中的被積函數進行插值，從而得到一組含有變限因子ε的全離散格式。

該方法相較于有限差分法，明顯提高了網格剖分的靈活性，使其能夠適用于非規(guī)則網格。能夠通過選取高精度的插值多項式來提高離散格式的精度。通過選取不同的ε，使得在達到同樣精度的條件下，得到的全離散格式滿足某種較好的性質。本文算法的優(yōu)勢主要體現在如下3個方面：

1）關于積分區(qū)域的選取。在對目標方程進行變限積分時，能夠以非規(guī)則網格區(qū)域為積分區(qū)域進行積分。這是網格剖分靈活性的體現。

2）關于變限因子ε的選取。選取不同的ε可得到不同的離散格式，從而能夠有目的地挑選其中滿足某種性質的離散格式，例如穩(wěn)定性。通過選取不同的ε，達到的計算效果也會有所不同。

3）關于插值多項式的選取。選取不同精度的插值多項式進行插值能夠直接影響離散格式的精度，因此，使用該方法能在理論上達到任意精度。甚至還能夠使用2種以上的插值多項式進行加權混合使用。

5 數值算例

對流擴散問題包含非穩(wěn)態(tài)和穩(wěn)態(tài)情況，下面運用本文離散格式分別對非穩(wěn)態(tài)和穩(wěn)態(tài)對流擴散方程進行數值實驗。

5.1 非穩(wěn)態(tài)情況

考慮一維非穩(wěn)態(tài)對流擴散方程，利用格式（6），解決如下初邊值問題：

其定態(tài)解析解為：u x，t()＝e5x-4t，方程的初邊值條件由解析解中得到。

于是，令L＝1，T＝1，用本文所提出的新算法對于不同的步長h和τ進行試算，可得結果如表1。從表1中能夠得到，隨著步長h和τ的減少，L2模誤差逐漸減小。變限因子ε1和ε2對實際計算結果有很大的影響，因此在表2中，在取定h和τ時，選取不同的變限因子進行誤差比較。rate為收斂階。

表1 ε1＝ε2＝h／2 000時的誤差比較Table1 The error comparison when ε1＝ε2＝h／2 000

表2 h＝τ＝0.01時的誤差比較Table2 The error comparison when h＝τ＝0.01

圖2 h＝τ＝0.01，不同ε1、ε2時的解析解和數值解Fig.2 The analytical and numerical solu?tions of different ε1，ε2，h＝τ＝0.01

圖3 h＝τ＝0.005，T＝5，ε1＝ε2＝h／240時的解析解和數值解Fig.3 The analytical and numerical solutions，h＝τ＝0.005，T＝5，ε1＝ε2＝h／240

圖2 為取時間步長h和空間步長τ均為0.01，取T＝1，取不同變限因子ε1、ε2，將解析解和數值解進行比較。圖3為取時間步長h和空間步長τ均為0.005，取T＝5，取變限因子ε1＝ε2＝h／240，將解析解和數值解進行比較。

圖4為取時間步長h和空間步長τ均為0.005，取變限因子ε1＝ε2＝h／240時，絕對誤差隨時間變化趨勢。從圖中可知，誤差主要集中在邊界。

圖4 h＝τ＝0.005，ε1＝ε2＝h／240，絕對誤差隨時間變化Fig.4 The absolute value of error changing with time，h＝τ＝0.005，ε1＝ε2＝h／240

計算結果表明，式（6）滿足理論精度。隨著時間的推移，數值解和解析解之間的絕對誤差逐漸減小，符合之前的證明分析。

5.2 穩(wěn)態(tài)情況

考慮一維穩(wěn)態(tài)對流擴散方程，利用格式（6），解決如下邊值問題：

5.2.1 流速u＝0.1 m／s，h＝0.2 m

由u＝0.1 m／s，h＝0.2 m，ρ＝1.0 kg／m3，Γ＝0.1 kg／（m·s），取ε1＝ε2時，得到離散方程組：

此時：

對離散方程組進行求解，結果列于表3中，傳統格式是指二階中心差分格式，由解可以看出，所求得的結果在物理上是真實的。

表3 u＝0.1 m／s，h＝0.2 m時，解析解、數值解與相對誤差Table3 Analytical solutions，numerical solutions and rela?tive error，u＝0.1 m／s，h＝0.2 m

5.2.2 流速u＝2 m／s，h＝0.2 m

求解過程5.2.1中的其他條件不變，流速由0.1 m／s提高到2 m／s。此時：

得到離散方程組：

表4 u＝2 m／s，h＝0.2 m時，解析解、數值解與相對誤差Table4 Analytical solutions，numerical solutions and rela?tive error，u＝2 m／s，h＝0.2 m

對離散方程組進行求解，結果列于表4中，傳統格式是指二階中心差分格式，由解可以看出，所求得的結果在物理上是不真實的，出現了振蕩。求解過程5.2.2節(jié)與求解過程5.2.1節(jié)相比較，只提高了流速，其他任何條件都未改變。

5.2.3 流速u＝2 m／s，h＝0.05 m

求解過程5.2.2節(jié)中的其他條件不變，加密網格。此時：

得到離散方程組：

對離散方程組進行求解，結果列于表5中，可見，求解過程5.2.2節(jié)中出現解的振蕩；求解過程5.2.3節(jié)中只是加密了網格，又得到了物理上真實的解。

表5 u＝2 m／s，h＝0.05 m時，解析解、數值解與相對誤差Table5 Analytical solutions，numerical solutions and rela?tive error，u＝2 m／s，h＝0.05 m

5.2.4 結果分析

在求解過程5.2.1中F／D＝0.1／0.5＝0.2；求解過程5.2.2中F／D＝2／0.5＝4；在求解過程5.2.3中F／D＝2／2＝1，見表5。傳統格式是指二階中心差分格式，可見在F／D較大時可能出現解的不真實，本文格式比傳統格式精度提高了許多。

6 結束語

本文針對一維對流擴散方程，并以一維非穩(wěn)態(tài)方程為例，基于變限積分的有限體積法，建立了一維的新的高精度數值離散格式，其精度能夠達到O h4／τ+h3+hτ+τ2／h+τ3／h2

( )，利用Fourier分析法對全離散格式進行了穩(wěn)定性分析，發(fā)現所構造的格式能在變限因子ε1≥ε2條件下能夠達到穩(wěn)定，通過數值算例證實了上述結論。ε的選取目前只是通過實驗獲得的經驗，還需要進一步總結和證明理論上的規(guī)則。本文在理論上給出了離散格式構造的方法以及精度估計和穩(wěn)定性分析的證明，通過算例驗證對工程實際有借鑒意義。進一步研究本課題將對高精度離散格式的構造具有重大意義。

［1］DEHGHAN M.Weighted finite difference techniques for the one?dimensional advection diffusion equation［J］.Appl Math Comput，2004，147：307?319.

［2］ROACHE P J.Computational fluid dynamics［M］.Hermosa Publishers，1972.

［3］SPALDING D B.A novel finite difference formulation for dif?ferential expressions involving both first and second deriva?tives［J］.International Journal for Numerical Methods in En?gineering，1972，4（4）：551?559.

［4］DEHGHAN M.Numerical solution of the three?dimensional advection?diffusion equation［J］.Applied Mathematics and Computation，2004，150（1）：5?19.

［5］ISENBERG J，GUTFINGER C.Heat transfer to a draining film［J］.International Journal of Heat and Mass Transfer，1973，16（2）：505?512.

［6］PARLANGE J Y.Water transport in soils［J］.Annual Re?view of Fluid Mechanics，1980，12（1）：77?102.

［7］FATTAH Q N，HOOPES J A.Dispersion in anisotropic，homogeneous，porous media［J］.Journal of Hydraulic Engi?neering，1985，111（5）：810?827.

［8］CHATWIN P C，ALLEN C M.Mathematical models of dis?persion in rivers and estuaries［J］.Annual Review of Fluid Mechanics，1985，17（1）：119?149.

［9］HOLLY J R F M，USSEGLIO?POLATERA J M.Dispersion simulation in two?dimensional tidal flow［J］.Journal of Hy?draulic Engineering，1984，110（7）：905?926.

［10］KUMAR N.Unsteady flow against dispersion in finite por?ous media［J］.Journal of Hydrology，1983，63（3）：345?358.

［11］GUVANASEN V，VOLKER R E.Numerical solutions for solute transport in unconfined aquifers［J］.International Journal for Numerical Methods in Fluids，1983，3（2）：103?123.

［12］KARAHAN H.Unconditional stable explicit finite differ?ence technique for the advection?diffusion equation using spreadsheets［J］.Advances in Engineering Software，2007，38（2）：80?86.

［13］KARAHAN H.Implicit finite difference techniques for the advection?diffusion equation using spreadsheets［J］.Ad?vances in Engineering Software，2006，37（9）：601?608.

［14］SUN H，ZHANG J.A high?order compact boundary value method for solving one?dimensional heat equations［J］.Nu?merical Methods for Partial Differential Equations，2003，19（6）：846?857.

［15］TIAN Z F.A rational high?order compact ADI method for unsteady convection?diffusion equations［J］.Computer Physics Communications，2011，182（3）：649?662.

［16］HIRSH R S.Higher order accurate difference solutions of fluid mechanics problems by a compact differencing tech?nique［J］.Journal of Computational Physics，1975，19（1）：90?109.

［17］CIMENT M，LEVENTHAL S H，WEINBERG B C.The op?erator compact implicit method for parabolic equations［J］.Journal of Computational Physics，1978，28（2）：135?166.［18］BAZAN F S V.Chebyshev pseudospectral method for com? puting numerical solution of convection?diffusion equation［J］.Applied Mathematics and Computation，2008，200（2）：537?546.

［19］MOHEBBI A，DEHGHAN M.High?order compact solution of the one?dimensional heat and advection?diffusion equa?tions［J］.Applied Mathematical Modelling，2010，34（10）：3071?3084.

［20］DEHGHAN M.Implicit collocation technique for heat equa?tion with non?classic initial condition［J］.International Journal of Nonlinear Sciences and Numerical Simulation，2006，7（4）：461?466.

［21］KARAA S，ZHANG J.High order ADI method for solving unsteady convection?diffusion problems［J］.Journal of Computational Physics，2004，198（1）：1?9.

［22］TIAN Z F，GE Y B.A fourth?order compact ADI method for solving two?dimensional unsteady convection?diffusion problems［J］.Journal of Computational and Applied Mathe?matics，2007，198（1）：268?286.

［23］王紅梅.對流擴散方程的特征有限元方法［D］.濟南：山東大學，2012.WANG Hongmei.Characteristics finite element methods for convection?diffusion problems［D］.Jinan：Shandong Uni?versity，2012.

［24］丁曉燕，馮秀芳.求解二維非定常對流擴散方程的高精度指數型差分方法［J］.寧夏大學學報：自然科學版，2014，35（1）：6?13.DING Xiaoyan，FENG Xiufang.A high?order exponential method for solving unsteady convection?diffusion equation［J］.Journal of Ningxia University：Natural Science Edi?tion，2014，35（1）：6?13.

［25］秦新強，王志剛，王全九，等.對流占優(yōu)擴散方程的楔形基無網格法［J］.工程數學學報，2013，30（5）：721?730.QIN Xinqiang，WANG Zhigang，WANG Quanjiu，et al.Meshless method with ridge basis functions for convection?dominated diffusion equations［J］.Chinese Journal of Engi?neering Mathematics，2013，30（5）：721?730.

［26］馬亮亮.變系數空間分數階對流?擴散方程的有限差分解法［J］.沈陽大學學報：自然科學版，2013，25（4）：341?344.MA Liangliang.Finite difference methods for space fraction?al convection?diffusion equation with variable coefficients［J］.Journal of Shenyang University：Natural Science，2013，25（4）：341?344.

［27］楊錄峰，李春光.一種求解對流擴散反應方程的高階緊致差分格式［J］.寧夏大學學報：自然科學版，2013，34（2）：101?104.YANG Lufeng，LI Chunguang.A high?order compact finite difference scheme for solving the convection diffusion reac?tion equations［J］.Journal of Ningxia University：Natural Science Edition，2013，34（2）：101?104.

Solving convection?diffusion equation by using the finite volume method with variable limit integral

ZHANG Lijian1，LUO Yuesheng2，ZHANG Wenping1
（1.College of Power and Energy Engineering，Harbin Engineering University，Harbin 150001，China；2.College of Science，Harbin Engineering University，Harbin 150001，China）

In this paper，the one?dimensional convection?diffusion equation using one?dimensional unsteady equations as an example figures out a new high precision fully discrete format that is based on the finite volume method with variable limit integral.This method first integrates the equation with variable limit integral in the control body and introduces a variable limit integral factor.After that it integrates the upper and lower limits to transform the differential equation into the integral equation and then uses the interpolation method to approximately replace the objective function.In this way，the accuracy of the diffusion term is improved，deriving the one?dimensional convergence precision.The Fourier analysis method is used to prove the conditional stability of formation.Finally，the numerical examples of non?stable and stable conditions are given to verify the advantages of high accuracy and the easiness of programming calculation.

variable limit integral；high precision；convection diffusion model；Fourier analysis；unsteady state；one?dimensional

10.3969／j.issn.1006?7043.201410049

http：／／www.cnki.net／kcms／detail／23.1390.U.20150309.1453.002.html

O551

A

1006?7043（2015）03?0427?05

2014?10?20.網絡出版時間：2015?03?09.

國家自然科學基金資助項目（51206031，51479038）.

張麗劍（1979?），女，講師，博士研究生；羅躍生（1960?），男，教授，博士生導師；張文平（1956?），男，教授，博士生導師.

張麗劍，E?mail：756778282＠qq.com.