邊界層對流擴散方程在自適應(yīng)網(wǎng)格上的高精度緊致格式

袁冬芳，曹富軍，葛永斌

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

0 引言

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

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

1 高精度緊致格式

考慮定常對流擴散方程

將區(qū)域[a,b]非均勻剖分為N個子區(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為截斷誤差,且

這里

為了得到高階精度，對截斷誤差項(9)中的三階和四階導(dǎo)數(shù)項進行處理，并利用方程(1)可得

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

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

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

(13)

其中

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

類似地，在y方向定義hy=(d-c)/Ny,hfy=yj+1-yj=θfyhy,hby=yj-yj-1=θbyhy.同時令αy=θfyθby,βy=θfy+θby,γy=θfy-θby,則二維對流擴散方程在非均應(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)時退化為標(biāo)準(zhǔn)四階緊致格式[10-11].

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

2.1 一維問題

(ω-1ξx)x=0,

其等價于

(ωxξ)ξ=0,

(15)

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

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

2.2 二維問題

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

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

類似地，第i列上每個區(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方向上交替實施一維網(wǎng)格自適應(yīng)算法.記

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

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

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

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

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

3 數(shù)值算例

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

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

該問題的精確解為

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

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

圖1給出了ε=10-4時均勻網(wǎng)格和自適應(yīng)網(wǎng)格上的計算結(jié)果，看出均勻網(wǎng)格的計算結(jié)果在邊界層附近產(chǎn)生很大的誤差，自適應(yīng)網(wǎng)格上的計算結(jié)果與精確解吻合得很好.表1給出了ε=10-2時均勻網(wǎng)格、非均勻網(wǎng)格和自適應(yīng)網(wǎng)格上的最大誤差和收斂階.從表1可以看出，在相同網(wǎng)格數(shù)下非均勻網(wǎng)格和自適應(yīng)網(wǎng)格上的誤差明顯比均勻網(wǎng)格上的誤差小2個量級，如均勻網(wǎng)格上網(wǎng)格數(shù)N=128時的最大誤差為1.984×10-4，然而自適應(yīng)網(wǎng)格上網(wǎng)格數(shù)N=32時的誤差即可達到相同的量級，因此自適應(yīng)網(wǎng)格能夠大大減少計算網(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考慮二維對流擴散問題

該問題的精確解為

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

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

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

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

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

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

4 結(jié)束語

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

參考文獻:

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

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

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

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

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

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

[7] 趙飛，陳建華，葛永斌.一維非定常對流擴散方程非均勻網(wǎng)格上的高階緊致差分格式[J].西安理工大學(xué)學(xué)報，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.