張麗劍，羅躍生，張文平

（1.哈爾濱工程大學(xué)動(dòng)力與能源工程學(xué)院，黑龍江哈爾濱150001；2.哈爾濱工程大學(xué)理學(xué)院，黑龍江哈爾濱150001）

變限積分的有限體積法解決對(duì)流擴(kuò)散方程

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

（1.哈爾濱工程大學(xué)動(dòng)力與能源工程學(xué)院，黑龍江哈爾濱150001；2.哈爾濱工程大學(xué)理學(xué)院，黑龍江哈爾濱150001）

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

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

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

1 格式的構(gòu)造

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

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

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

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

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

式中：ε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）轉(zhuǎn)化為積分方程（5），于是用二元六點(diǎn)拉格朗日插值法，對(duì)u x，t()進(jìn)行插值。整理并舍去小項(xiàng)，可得方程（1）的全離散格式：

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

2 精度估計(jì)

由于在對(duì)方程（1）進(jìn)行離散時(shí)，一共對(duì)x進(jìn)行了3次積分，對(duì)t進(jìn)行了1次積分，則等式兩邊應(yīng)除以h3τ，從而可得離散格式（6）的截?cái)嗾`差為O h4／τ+h3+hτ+τ2／h+τ3／h2

( )，顯然，該截?cái)嗾`差大于二階。

3 穩(wěn)定性分析

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

顯然可得G′1

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

4 算法分析

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

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

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

1）關(guān)于積分區(qū)域的選取。在對(duì)目標(biāo)方程進(jìn)行變限積分時(shí)，能夠以非規(guī)則網(wǎng)格區(qū)域?yàn)榉e分區(qū)域進(jìn)行積分。這是網(wǎng)格剖分靈活性的體現(xiàn)。

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

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

5 數(shù)值算例

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

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

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

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

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

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

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

圖2 h＝τ＝0.01，不同ε1、ε2時(shí)的解析解和數(shù)值解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時(shí)的解析解和數(shù)值解Fig.3 The analytical and numerical solutions，h＝τ＝0.005，T＝5，ε1＝ε2＝h／240

圖2 為取時(shí)間步長(zhǎng)h和空間步長(zhǎng)τ均為0.01，取T＝1，取不同變限因子ε1、ε2，將解析解和數(shù)值解進(jìn)行比較。圖3為取時(shí)間步長(zhǎng)h和空間步長(zhǎng)τ均為0.005，取T＝5，取變限因子ε1＝ε2＝h／240，將解析解和數(shù)值解進(jìn)行比較。

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

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

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

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

考慮一維穩(wěn)態(tài)對(duì)流擴(kuò)散方程，利用格式（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時(shí)，得到離散方程組：

此時(shí)：

對(duì)離散方程組進(jìn)行求解，結(jié)果列于表3中，傳統(tǒng)格式是指二階中心差分格式，由解可以看出，所求得的結(jié)果在物理上是真實(shí)的。

表3 u＝0.1 m／s，h＝0.2 m時(shí)，解析解、數(shù)值解與相對(duì)誤差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。此時(shí)：

得到離散方程組：

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

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

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

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

得到離散方程組：

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

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

5.2.4 結(jié)果分析

在求解過程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。傳統(tǒng)格式是指二階中心差分格式，可見在F／D較大時(shí)可能出現(xiàn)解的不真實(shí)，本文格式比傳統(tǒng)格式精度提高了許多。

6 結(jié)束語(yǔ)

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

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

［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］王紅梅.對(duì)流擴(kuò)散方程的特征有限元方法［D］.濟(jì)南：山東大學(xué)，2012.WANG Hongmei.Characteristics finite element methods for convection?diffusion problems［D］.Jinan：Shandong Uni?versity，2012.

［24］丁曉燕，馮秀芳.求解二維非定常對(duì)流擴(kuò)散方程的高精度指數(shù)型差分方法［J］.寧夏大學(xué)學(xué)報(bào)：自然科學(xué)版，2014，35（1）：6?13.DING Xiaoyan，F(xiàn)ENG 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］秦新強(qiáng)，王志剛，王全九，等.對(duì)流占優(yōu)擴(kuò)散方程的楔形基無網(wǎng)格法［J］.工程數(shù)學(xué)學(xué)報(bào)，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］馬亮亮.變系數(shù)空間分?jǐn)?shù)階對(duì)流?擴(kuò)散方程的有限差分解法［J］.沈陽(yáng)大學(xué)學(xué)報(bào)：自然科學(xué)版，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］楊錄峰，李春光.一種求解對(duì)流擴(kuò)散反應(yīng)方程的高階緊致差分格式［J］.寧夏大學(xué)學(xué)報(bào)：自然科學(xué)版，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.網(wǎng)絡(luò)出版時(shí)間：2015?03?09.

國(guó)家自然科學(xué)基金資助項(xiàng)目（51206031，51479038）.

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

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