袁冬芳，曹富軍，葛永斌

(1.內(nèi)蒙古科技大學(xué) 理學(xué)院，內(nèi)蒙古 包頭 014010;2.寧夏大學(xué) 數(shù)學(xué)統(tǒng)計(jì)學(xué)院，寧夏 銀川 750021)

0 引言

對流擴(kuò)散方程廣泛應(yīng)用于許多實(shí)際物理問題的建模，如核廢料污染、滲流驅(qū)動、海水入侵等，所以研究其精確、穩(wěn)定和高效的數(shù)值方法具有十分重要的意義.當(dāng)物理解足夠光滑時(shí)一致網(wǎng)格上的高階格式就可以得到滿意的數(shù)值解，然而當(dāng)物理解存在間斷或大梯度時(shí)，采用均勻網(wǎng)格計(jì)算足夠的精度就需要比較多的網(wǎng)格點(diǎn),從而增加了存儲空間、計(jì)算量和計(jì)算時(shí)間.合理的做法是在邊界層和大梯度附近分布比較多的網(wǎng)格點(diǎn)，在物理解變化比較平緩的區(qū)域分布較少網(wǎng)格點(diǎn).針對帶邊界層和大梯度問題的對流擴(kuò)散方程高階格式的研究已經(jīng)引起了國內(nèi)外研究者的廣泛關(guān)注.文獻(xiàn)[1-2]基于非均勻網(wǎng)格上函數(shù)的泰勒級數(shù)展開,結(jié)合殘量修正法,推導(dǎo)了非均勻網(wǎng)格上對流擴(kuò)散方程的高階指數(shù)型緊致差分格式,結(jié)果表明該格式兼有高精度和高分辨率的優(yōu)點(diǎn),能夠很好地適用于大梯度變化.文獻(xiàn)[3]基于泰勒級數(shù)展開，構(gòu)造了一種非均勻網(wǎng)格三點(diǎn)四階精度的緊致差分格式，并對Burgers方程和對流方程進(jìn)行求解.文獻(xiàn)[4]基于泰勒級數(shù)展開法提出了求解一維定常對流擴(kuò)散方程非均勻網(wǎng)格上的具有3～4階精度的緊致差分格式.曹廣滿等[5]針對一維定常對流擴(kuò)散方程提出了一種二階非等距網(wǎng)格上的差分格式.獻(xiàn)文[6-9]提出了非均勻網(wǎng)格上求解一維非定常對流擴(kuò)散方程的緊致差分格式.文獻(xiàn)[10-11]提出了求解二維定常對流擴(kuò)散方程非均勻網(wǎng)格上的高精度緊致差分格式.然而，以上文獻(xiàn)采用的非均勻網(wǎng)格均由網(wǎng)格分布函數(shù)來確定，且對于不同的算例需要指定不同的網(wǎng)格分布函數(shù).同時(shí)需要不斷調(diào)整網(wǎng)格伸縮系數(shù)從而對網(wǎng)格的疏密程度進(jìn)行控制，網(wǎng)格分布函數(shù)及伸縮參數(shù)的選擇直接影響數(shù)值計(jì)算的精度，很大程度上限制了該方法的適用性.

自適應(yīng)網(wǎng)格算法根據(jù)實(shí)際問題需要在物理解變動較大的區(qū)域網(wǎng)格自動密集，而在物理解變化平緩區(qū)域網(wǎng)格相對稀疏，對于合理分布網(wǎng)格及高效和精確計(jì)算起著極其重要的作用[12-15].據(jù)我們所知，基于自適應(yīng)網(wǎng)格的高階緊致格式求解邊界層對流擴(kuò)散方程的研究鮮有報(bào)道.文中針對帶邊界層對流擴(kuò)散問題提出自適應(yīng)網(wǎng)格方法，根據(jù)物理解的特征對網(wǎng)格進(jìn)行自適應(yīng)加密，然后結(jié)合高精度緊致差分格式對一維邊界層對流擴(kuò)散方程進(jìn)行求解.

1 高精度緊致格式

考慮定常對流擴(kuò)散方程

將區(qū)域[a,b]非均勻剖分為N個(gè)子區(qū)間:a≤x0

將(2)和(3)式相減，整理可得

將(2)和(3)式分別乘以xf和xb后相加，整理得

非均勻網(wǎng)格上的一階和二階導(dǎo)數(shù)的中心差分算子可表示為

將(6)和(7)式代入(1)式可得

(8)

其中，pi=p(xi),fi=f(xi)，τi為截?cái)嗾`差,且

這里

為了得到高階精度，對截?cái)嗾`差項(xiàng)(9)中的三階和四階導(dǎo)數(shù)項(xiàng)進(jìn)行處理，并利用方程(1)可得

將(10)和(11)式代入(9)可得

將(12)式代入(8)式并整理可得非均勻網(wǎng)格上的高階緊致差分格式

βui-1+αui+γui+1=Fi,

(13)

其中

由泰勒級數(shù)展開與截?cái)嗾`差分析可知，格式(13)具有3～4階精度.當(dāng)xf≠xb時(shí)，格式(13)具有3階精度；當(dāng)xf=xb時(shí)退化為均勻網(wǎng)格，格式具有4階精度.

類似地，在y方向定義hy=(d-c)/Ny,hfy=yj+1-yj=θfyhy,hby=yj-yj-1=θbyhy.同時(shí)令αy=θfyθby,βy=θfy+θby,γy=θfy-θby,則二維對流擴(kuò)散方程在非均應(yīng)網(wǎng)格上的緊致格式為[10-11]:

(14)

其中Aij,Bij,Cij,Dij,Gij,Hij,Kij,Lij分別為

顯見，格式(14)在非均勻網(wǎng)格(θx≠θy)上至少具有三階精度,且網(wǎng)格為均勻網(wǎng)格(θx=θy=1)時(shí)退化為標(biāo)準(zhǔn)四階緊致格式[10-11].

2 自適應(yīng)網(wǎng)格

2.1 一維問題

(ω-1ξx)x=0,

其等價(jià)于

(ωxξ)ξ=0,

(15)

(16)式將根據(jù)物理解的特點(diǎn)對網(wǎng)格進(jìn)行自適應(yīng)加密，迭代求解(16)式至收斂將得到新的網(wǎng)格分布.

以上網(wǎng)格自適應(yīng)算法過程描述如下：

2.2 二維問題

考慮矩形計(jì)算區(qū)域Ω=[a,b]×[c,d]，將區(qū)域進(jìn)行等距離散并將網(wǎng)格表示為{coord[0]}，即

記均勻網(wǎng)格上的初始數(shù)值解為{u[0]}.令sgr(j)為第j行上每個(gè)區(qū)間上梯度的和，即

類似地，第i列上每個(gè)區(qū)間上梯度的和可表示為{sgc(i),i=0,1,…,Nx}.令jmax和imax分別表示為行和列上梯度變化最大的行和列號，即

記rmax和cmax分別表示jmax行和imax列上梯度最大的網(wǎng)格，即

記ylay=rmax/sgr(jmax),xlay=cmax/sgc(imax),ylay和xlay用于監(jiān)測邊界層的存在且根據(jù)值的不同，可以分為以下4種情形:

以情形 4為例，在x,y方向都存在邊界層，則需要分別在x和y方向上交替實(shí)施一維網(wǎng)格自適應(yīng)算法.記

當(dāng)網(wǎng)格更新為{coord[1]}時(shí)，將初始值{u[0]}更新為{u(1)}，即

重復(fù)以上過程，直到

||coord[v+1]-coord[v]||<ε2,

則可以得到新的網(wǎng)格.

對于情形2和情形3可以采用類似的方式分別在x和y方向上實(shí)施一維自適應(yīng)算法.

3 數(shù)值算例

為了驗(yàn)證文中方法的精確性和可靠性，在定義域[0,1]內(nèi)針對以下兩個(gè)有精確解的對流擴(kuò)散問題進(jìn)行數(shù)值實(shí)驗(yàn)，并比較均勻網(wǎng)格和自適應(yīng)網(wǎng)格上計(jì)算結(jié)果的精度和收斂階.

算例1考慮線性常系數(shù)對流擴(kuò)散問題

該問題的精確解為

該問題在x=1附近有一個(gè)邊界層，因此采用如下變換函數(shù)生成非均勻網(wǎng)格:

其中,N是區(qū)域剖分后子區(qū)間的個(gè)數(shù)，λ是伸縮變化系數(shù)，用來調(diào)節(jié)網(wǎng)格點(diǎn)在某一區(qū)域的密集程度[4-8].

圖1給出了ε=10-4時(shí)均勻網(wǎng)格和自適應(yīng)網(wǎng)格上的計(jì)算結(jié)果，看出均勻網(wǎng)格的計(jì)算結(jié)果在邊界層附近產(chǎn)生很大的誤差，自適應(yīng)網(wǎng)格上的計(jì)算結(jié)果與精確解吻合得很好.表1給出了ε=10-2時(shí)均勻網(wǎng)格、非均勻網(wǎng)格和自適應(yīng)網(wǎng)格上的最大誤差和收斂階.從表1可以看出，在相同網(wǎng)格數(shù)下非均勻網(wǎng)格和自適應(yīng)網(wǎng)格上的誤差明顯比均勻網(wǎng)格上的誤差小2個(gè)量級，如均勻網(wǎng)格上網(wǎng)格數(shù)N=128時(shí)的最大誤差為1.984×10-4，然而自適應(yīng)網(wǎng)格上網(wǎng)格數(shù)N=32時(shí)的誤差即可達(dá)到相同的量級，因此自適應(yīng)網(wǎng)格能夠大大減少計(jì)算網(wǎng)格，充分說明了其優(yōu)越性.從表1的收斂階可見，自適應(yīng)網(wǎng)格與非均勻網(wǎng)格上的收斂階基本保持4階精度，明顯高于均勻網(wǎng)格上的精度.

(a)均勻網(wǎng)格 (b)自適應(yīng)網(wǎng)格

圖1 均勻網(wǎng)格和自適應(yīng)網(wǎng)格上的數(shù)值解(N=65,ε=10-4)

Fig 1 The exact solution and numerical solutions under uniform mesh and adaptive mesh (N=65,ε=10-4)

表1 均勻網(wǎng)格與自適應(yīng)網(wǎng)格上的最大誤差與收斂階,ε=10-2

算例2考慮二維對流擴(kuò)散問題

該問題的精確解為

u(x,y)=ey-x+21/ε(1+y)1+1/ε.

圖2給出了網(wǎng)格Nx=Ny=33,ε=10-3時(shí)精確解以及均應(yīng)網(wǎng)格、非均勻網(wǎng)格和自適應(yīng)網(wǎng)格上的數(shù)值解.圖3比較了x=0,Nx=33,ε=10-3時(shí)均應(yīng)網(wǎng)格和自適應(yīng)網(wǎng)格上精確解和數(shù)值解的比較，從圖中可見，自適應(yīng)網(wǎng)格在邊界層附近分布更多的網(wǎng)格點(diǎn)，且數(shù)值解與精確解一致.

圖2 精確解以及不同網(wǎng)格上的數(shù)值解

圖3 均勻網(wǎng)格與自適應(yīng)網(wǎng)格上的精確解與數(shù)值解

表2比較了ε=10-2和ε=10-3時(shí)不同網(wǎng)格下的最大誤差、CPU時(shí)間和收斂階，其中非均勻網(wǎng)格的伸縮系數(shù)分別為λx=0.0,λy=0.55和λx=0.0,λy=0.9.從表2可見,當(dāng)ε=10-2時(shí)，3種網(wǎng)格上的最大誤差和收斂階都接近于四階精度;當(dāng)ε=10-3時(shí)，均勻網(wǎng)格喪失了精度，非均勻和自適應(yīng)網(wǎng)格依然能夠得到比較合理的誤差和收斂階.需要指出的是,非均勻網(wǎng)格需要在已知邊界層位置的情況下根據(jù)不同的網(wǎng)格和參數(shù)指定合理的網(wǎng)格伸縮系數(shù)，自適應(yīng)網(wǎng)格在不知邊界層位置的情況下可以根據(jù)初始值對網(wǎng)格進(jìn)行迭代從而使其重新分布.

表2 不同網(wǎng)格上的最大誤差、CPU時(shí)間及收斂階

4 結(jié)束語

提出了二維計(jì)算區(qū)域上正交網(wǎng)格的自適應(yīng)方法，并結(jié)合高精度緊致差分格對一維和二維對流擴(kuò)散方程大梯度和邊界層問題進(jìn)行求解.該方法可以在現(xiàn)有代碼中實(shí)現(xiàn)，也可以應(yīng)用于邊界層與坐標(biāo)軸平行的問題.其根據(jù)初值采用迭代的方式進(jìn)行網(wǎng)格自適應(yīng)調(diào)整，有效改進(jìn)了現(xiàn)有方法中事先指定邊界層位置，并采用網(wǎng)格分布函數(shù)調(diào)整網(wǎng)格控制參數(shù)進(jìn)行非均勻網(wǎng)格生成的缺點(diǎn).數(shù)值算例表明：自適應(yīng)網(wǎng)格方法能夠合理地調(diào)節(jié)網(wǎng)格疏密，同時(shí)結(jié)合自適應(yīng)網(wǎng)格方法與高階格式對邊界層問題進(jìn)行求解可以提高數(shù)值解的精度，減少計(jì)算網(wǎng)格，降低計(jì)算量.

參考文獻(xiàn):

[1] 田芳.一維對流擴(kuò)散方程非均勻網(wǎng)格上的指數(shù)型高階緊致差分格式[J].數(shù)學(xué)的實(shí)踐與認(rèn)識,2015(4):268.

[2] 田芳,田振夫.非均勻網(wǎng)格上求解對流擴(kuò)散問題的高階緊致差分方法[J].寧夏大學(xué)學(xué)報(bào)(自然科學(xué)版),2009,30(3):209.

[3] 孫建安，賈偉，吳廣智.一種非均勻網(wǎng)格上的高精度緊致差分格式[J].西北師范大學(xué)學(xué)報(bào)(自然科學(xué)版),2014,51(4):31.

[4] 薛文強(qiáng),蘭斌,葛永斌.一維定常對流擴(kuò)散方程非均勻網(wǎng)格上的高階緊致差分格式[J].西北師范大學(xué)學(xué)報(bào)(自然科學(xué)版),2013,50(4):16.

[5] 曹廣滿，王彩華，齊海濤.對流擴(kuò)散方程的非一致網(wǎng)格有限差分方法[J].天津師范大學(xué)學(xué)報(bào)(自然科學(xué)版)，2010,30(1):7.

[6] 黃雪芳,郭 銳,葛永斌.一維非定常對流擴(kuò)散方程非均勻網(wǎng)格上的高精度緊致差分格 式[J].工程數(shù)學(xué)學(xué)報(bào),2014(3):371.

[7] 趙飛，陳建華，葛永斌.一維非定常對流擴(kuò)散方程非均勻網(wǎng)格上的高階緊致差分格式[J].西安理工大學(xué)學(xué)報(bào)，2013，29(4):475.

[8] MOHEBBI A,DEHGHAN M.High-order compact solution of the one-dimensional heat and advection-diffusion equations[J].AppliedMathematicalModelling,2010,34(10):3071.

[9] TIAN Z F,YU P X.A high-order exponential scheme for solving 1D unsteady convection-diffusion equations[J].JournalofComputationalandAppliedMathematics,2011,235(8):2477.

[10] KALITA C,DASS A K,DALAL D C.A transformational-free HOC scheme for steady convection-diffusion on non-uniform grids[J].InternationalJournalforNumericalMethodsinFluids,2004,44(1):33.

[11] GE Y B,CAO F J.Multigrid method based on the transformation-free HOC scheme on nonuniform grids for 2D convection diffusion problems[J].JournalofComputationalPhysics,2011,230(10):4051.

[12] TANG H Z,TANG T.Adapative mesh methods for one-and two-dimensional hyperbolic conservation laws[J].SIAMNumericalAnalysis,2003,41(2):487.

[13] YANG X,HUANG W,QIU J.A moving mesh WENO method for one-dimensional conservation laws[J].SIAMJournalofScientificComputing,2012,34:2317.

[14] HU F X,WANG R,CHEN X.Y,et al.An adaptive mesh method for 1D hyperbolic conservation laws[J].AppliedNumericalMathematics,2015,91:11.

[15] Van DAM A,ZEGELING P A.A robust moving mesh finite volume method applied to 1D hyperbolic conservation laws from magneto-hydrodynamics[J].JournalofComputationalPhysics,2006,216(2):526.