高精度迎风偏斜格式的比较与分析

第31卷第2期2007年3月

大　气　科　学

ChineseJournalofAtmosphericSciencesVol131　No12

Mar.2007

高精度迎风偏斜格式的比较与分析

冯涛1,2　李建平1

1　1000292中国科学院研究生院,北京　100039

摘　要　,并利用Fourier分析法。结果表明,偶阶精度格式的数值相速度快于实际相速度,而奇。并且,偶阶精度格式的耗散误差与频散误差低于相邻的奇阶精度格式。,在一维问题上进行了应用。首先,考虑恒定风场条件下的一维平流试验。主要选择两种不同的初始条件来评价数值格式的精度,这两种试验问题分别是高斯函数、方波函数。试验结果表明,随着数值格式精度的提高,数值格式的误差逐渐减小。而对于高于六阶精度的格式来说,改进的程度并不是很大。其次,应用各阶格式到具有两种不同初始条件的无粘Burgers方程。数值结果表明,随着数值格式阶数的增加,数值结果也得到了明显改进。而对于高于六阶精度的格式来说,进一步的变化并不明显。总之,在兼顾效率与精度条件下六阶迎风偏斜格式是最好的。关键词　显式有限差分公式　迎风偏斜格式　平流方程

文章编号　10069895(2007)02024509　　　中图分类号　P435　　　文献标识码　A

AComparisonandAnalysisofHighOrderUpwind2BiasedSchemes

FENGTao1,2andLIJian2Ping1

1StateKeyLaboratoryofNumericalModelingforAtmosphericSciencesandGeophysicalFluidDynamics,InstituteofAtmos2

phericPhysics,ChineseAcademyofSciences,Beijing　100029

2GraduateUniversityoftheChineseAcademyofSciences,Beijing　100049

Abstract　Highorderupwind2biasedschemesbasedonthegeneralexplicitdifferenceformulaswitharbitraryorderaccuracyfromLi(2005)aredeveloped1Thedissipationanddispersionerrorsoftheseupwind2biasedschemesareas2sessedbymeansoftheFourieranalysis1Theresultsshowthatthenumericalphasespeedsforevenorderschemesarefasterthantherealphasespeed1Butthenumericalphasespeedsforoddorderschemesareslowerthantherealphasespeed1Moreover,thedissipationanddispersionerrorsoftheevenorderschemearesmallerthanthenextlower2orderoddscheme1Fortheupwind2biasedschemesderived,thereisanincreaseofthedissipationerrorinthehighwavenumberrange1However,theincreaseofnumericaldissipationinthelowwavenumberrangeisverysmall1Somenumericaldissipationisneededinthehighwavenumberrangeforsomeflowproblems1Thedecisionthatanumericalmodelershouldalwaysmakeishowmuchnumericaldissipationisneed1Theupwind2biasedschemesprovidesomeadditionalalternativesforadjustingnumericaldissipation1Theseschemesareappliedtoaone2dimensionalmodelinordertotesttheircomputationalperformance1

Firstly,considerone2dimensionaladvectioninaconstantvelocityfield1ThetwotypesoftestproblemsthatareusedtoevaluatetheaccuracyofanumericalschemearetheGaussianfunctionandthesquarewavefunction1Thefirstexampleistheadvectionofthesquarewavefunc2收稿日期　20050711,20060205收修定稿

资助项目　国家自然科学基金资助项目40325015,中国科学院创新团队国际合作伙伴计划项目“气候系统模式研发及应用研究”作者简介　冯涛,男,1977年出生,硕士,主要从事计算地球流体力学与数值天气预报等方面的研究。E2mail:[email protected]

246

大　气　科　学

ChineseJournalofAtmosphericSciences　31卷

Vol1

tion1ThisfunctionrevealsanumericalmethodπscapabilitytohandleGibbπsoscillationsthatariseinthevicinityofdiscontinuities1ThesecondexampleistheadvectionoftheGaussianfunction1Thisfunctionprovidesademonstra2tionofeachschemeπsabilitytotransportwellresolved,smoothlyvaryingfunctionsoverlargerdistances1Testingresultsshowthattheerrorsofnumericalschemesgraduallydecreaseastheorderoftheschemeincreasesandfurtherchangesaremoderateforhigherthansixthorderscheme1Ahighorderscheme,althoughmaintainingapositivedefinitefield,willreducetheproblemofnegativequantities1inviscidBurgersπequation1TheinviscidBurgersπequationisusedasthetestisthethatallowsscalecollapse(shockformation),andbecauseithasanalytic1importantphenomenaintheat2mospherethatareassociatedwithsharpgradients1methodforthescalecollapsephenomenawouldseemtobeanabilityverysharpgradientinthecomputationaldomain1Inaddition,anynumericaltheofashockshouldnotcontaminatethesmoothsolutionawayfromtheshowthatthereisadrasticimprovementwhengoingfromthefirsttosixthorder,whilefurtheraremoderate1Theresultsindicatetheusefulnessoftheupwind2biasedschemeinhan2dlingtheadvectionneardiscontinuitiesandthedevelopmentofscalecollapseregionsinatmosphericnumericalmod2els1Inconclusion,thesixth2orderupwind2biasedschemeappearstobethebestbalancebetweenefficiencyandaccu2racy1Itappearsthatthesixth2orderupwind2biasedschemeiswellsuitedformanyatmosphericmodelingapplicationswhereadvectionplaysasignificantrole1

Keywords　explicitfinitedifferenceformulas,upwind2biasedscheme,advectionequation

1　引言

有限差分格点所能分辨的最短波长是2倍格距波。高频波的存在不但会引入计算误差,而且会导致数值不稳定。迎风有限差分法可以通过引入数值耗散来有效衰减高波数部分。然而,数值耗散有可能超越物理扩散。这样,迎风格式需要合理地设计来产生适宜的数值耗散。为了减小数值耗散,已经发展了二阶Crowley格式[1]、基于三次迎风插值的三阶格式[2]、单调保持的TVD(TotalVariationDiminishing)格式[3]、应用JFNK(Jacobin2freeNewton2Krylov)方法求解的二阶全隐式差分格

点上的值通过Lagrange多项式插值得到。最近,Li[7]给出了具有任意阶精度的一般显式有限差分公

式,此公式是以泰勒展开式为基础,利用范德蒙行

列式性质而得到的。利用此公式可以直接得到基于等距与非等距网格的一阶或更高阶导数的数值近似值,还可以构造出微分方程和偏微分方程的差分格式。本文利用此一般显式有限差分公式,构造出基于等距网格的一阶至十阶精度迎风偏斜格式,同时对这些迎风偏斜格式进行比较分析,并得到一个最优的格式。

2　高精度有限差分方法

考虑如下一维线性平流方程及对应的半离散化方程:

+c=0,tx

+c=0,

Δxt

(1)(2)

式[4]等等。由于需要控制数值耗散与改进数值解精度,高精度迎风格式的发展显得尤其重要。高阶迎风有限差分格式因其较小的数值耗散已经被应用于湍流数值模拟。Rai和Moin[5]利用具有五阶精度的迎风偏斜差分格式,直接模拟不可压槽道湍流问题,其五阶迎风偏斜逼近式是由非等距Lagrange插值多项式而得。由于Lagrange插值多项式的系数求解非常复杂,因此用它来构造高阶精度有限差分逼近格式是困难的[5]。Tremback等[6]提出了一般时间向前差迎风平流算子,比较了一阶到十阶迎风平流格式,在同时考虑效率与精度时得出六阶格式是一种最好的格式。此迎风平流格式是利用网格

式中,u=u(x,t),c取为常数,本文c选为大于

零,Fj/Δx为5u/5x的差分逼近式。

一阶导数的一般差分逼近式可以直接通过在网点上函数值的线性组合来表示。(2)式中的导数逼近式Fj/Δx可直接通过离散函数uj的线性组合来表示。Li[7]给出如下具有任意阶精度的一般显式有限差分公式,

2期　No12冯涛等:高精度迎风偏斜格式的比较与分析

FENGTaoetal.AComparisonandAnalysisofHighOrderUpwind2BiasedSchemes

　　　247

(m)

(xi)=

j=0

∑

(m)(m)n-m+1

),(3)dn+1,i,jf(xj)+On,j(h

(xj)=f′

(

fj-2-4fj-1+3fj),2h

(11)

其中,xi=x0+ih(i=0,1,…,n,h≠0),f是所选定点的函数。(1)

d2,0,1=1,　　n=1,

(4)(1)

d2,1,0=-1,　　n=1,

(m-1)

(m-j)(m)

dn+1,i,j=,　j≠i,n>1,(5)

j!(n-j)!

(m)

n+1,i,j

三阶迎风偏斜格式:

(xj)=(ff′

j-2

-6f

j-1

3fj+2fj+1),(12)

(j)f′

(

-j-3+fj-2-18fj-1+10fj+h

(13)

j1,

(dn+1,i,i=-(

)

(m)

j=0,j≠i

∑d

,　n≥1,

(xj)=f′

m-1

an-1,i,j=am-1(-i,ki,i),

　　　　　　　　(k,…n,　k≠i,j)(0)

an=a0(x1,x2,…,xn)=x1x2…xn,

(1)n

(7)

(

-2fj-3+15fj-2-60fj-1+60h

(14)

　　20fj+30fj+1-3fj+2),六阶迎风偏斜格式:

(xj)=f′

=a1(x1,x2,…,xn)=x1x2…xn-1+=a2(x1,x2,…,xn)=x1x2…xn-2+…+

　　x1x2…xn-2xn+…+x1x3…xn+x2x3…xn,

(2)

(

fj-4-8fj-3+30fj-2-80fj-1+60h

(15)

　　35fj+24fj+1-2fj+2),七阶迎风偏斜格式:

(xj)=f′

　　x3x4…xn,…

(n-2)

=an-2(x1,x2,…,xn)=x1x2+

=an-1(x1,x2,…,xn)=x1+x2+…+xn,=an(x1,x2,…,xn)=1,

(8)

(

3fj-4-28fj-3+126fj-2-420h

(16)

　　x2x3+…+xn-1xn,

anan

(n-1)(n)

　　420fj-1+105fj+252fj+1-42fj+2+4fj+3),八阶迎风偏斜格式:

(xj)=f′

余项

(m)

n,j

(

-3fj-5+30fj-4-140fj-3+

840h

(17)

(n-m+1)

n-mn-m+1

)=・

(n+1)!

(n+1)

(m-1)n+1

(-1)jj!(n-j)!

　　420fj-2-1050fj-1+378fj+420fj+1-　　60fj+2+5fj+3),

(9)

j=0,j≠i

∑

九阶迎风偏斜格式:

(xj)=f′

　　可根据精度的需要,选择适当个数的网格基架

点数。构成二阶精度的逼近式需要三个网格基架点,三个网格基架点可以构造三种二阶精度的差分逼近式。构造N阶精度的格式需要N+1个网格基架点,N+1个网格基架点可以构造N+1种N阶精度的差分逼近式。在本文中,导数的逼近式中更多地利用了上游点上的函数值(指相对于波的传播方向而言),这类格式称为迎风偏斜格式。这一点从物理上考虑是很容易理解的,因为波是从上游传下来的。下面给出具体的逼近一阶导数的一阶至十阶迎风偏斜格式

:一阶迎风偏斜格式:

(xj)=(-fj-1+fj),(10)f′

(

-4fj-5+45fj-4-240fj-3+

2520h

(18)

　　840fj-2-2520fj-1+504fj+1680fj+1-　　360fj+2+60fj+3-5fj+4),十阶迎风偏斜格式:

(xj)=f′

(

2fj-6-24fj-5+135fj-4-2520h

(19)

　　480fj-3+1260fj-2-3024fj-1+924fj+　　1440fj+1-270fj+2+40fj+3-3fj+4)1

3　数值格式的Fourier分析

对线性问题,当求解区域为[-∞,+∞]或在有限区域[-L,+L]的边界上满足周期条件,则问题的解可通过Fourier级数来表示。Fourier分析法[8]为取初值为单个Fourier分量,用以研究数值解的精度与行为。通过Fourier分析法可了解不同

二阶迎风偏斜格式:

248

大　气　科　学

ChineseJournalofAtmosphericSciences　31卷

Vol131

尺度量离散导数的逼近精度。在单纯进行精度分析时并不要求被离散的函数是线性的。现在考虑线性问题,取初值为u(x,0)=exp(ikx),方程(1)的准确解为

　u(x,t)=exp[

ik(x-ct)],准确解为　u(xj,t)=exp

-ck

texpa

ikxjta

,(20)

中高波分量有较好的模拟能力。此外,偶阶精度迎风偏斜格式的频散性明显小于相邻的奇阶精度迎风偏斜格式。从图2可以看出,各阶格式在低波段有较小的数值耗散,。取u(xj,0)=exp(ikxj)为初值,则差分方程(2)的

(21)

4　线性平流试验

为了进一步比较一阶至十阶迎风偏斜格式的性能,我们以方波和高斯波的初值问题为例进行计算。

时间导数的差分逼近的方法有:(1)采用多层结构,这就要求知道前时间层的值,如三层结构;(2)Lax2Wendroff型的多步格式。该方法对于常系数的情况,可构造简单且计算量小的高阶精度差分格式,但在一般非常系数的情况下,方法极为复杂;(3)R2K型的多步法。该方法发展较为成熟,逻辑简单,且可占用较少的内存。

时间导数的差分逼近采用三阶精度的R2K三步方法[9]。设L(u)为5u/5x的差分逼近式,则对于方程(1),三阶精度的三步方法有以下形式

uuu

(1)

Δx,h,p与数值波的传其中a=k

播速度有关h与p的分

析表达式来。对照式(20)和式(21)可知,对取定波数k,数值解逼近准确解,要求

Δx),　h/a→0,p/a→1　(a=k

(22)

它们的逼近程度反映了格式的精度。下面给出一阶至十阶迎风偏斜格式h与p随a的变化情况。

在图1和图2中,给出了一阶至十阶迎风偏斜格式的频散误差和耗散误差,其中曲线1、2、3、4、5、6、7、8、9、10依次对应为一阶格式、二阶格式、

三阶格式、四阶格式、五阶格式、六阶格式、七阶格式、八阶格式、九阶格式、十阶格式,曲线0为p

Δx曲线,它对应于微分方程的准确解。从=a=k

图1可以看出,对低波分量,各种格式都能较好地逼近准确解。并且,偶阶精度格式的数值相速度快于实际相速度,而奇阶精度格式的数值相速度慢于实际相速度。与低阶精度格式相比,高精度格式对

=u==

(0)

ΔtL(u(0)),+c

(23)

)

(0)(1)(1)

u+u+cΔtL(u

),444(0)(2)(2)

u+u+cΔtL(u),333

(3)

图1　一阶至十阶迎风偏斜格式的频散误差

Fig11　Thedispersionerrorsofupwind2biasedschemesfororders1through10

图2　一阶至十阶迎风偏斜格式的耗散误差

Fig12　Thedissipationerrorsofupwind2biasedschemesfororders1through10

2期　No12冯涛等:高精度迎风偏斜格式的比较与分析

FENGTaoetal.AComparisonandAnalysisofHighOrderUpwind2BiasedSchemes

　　　249

空间导数的差分逼近选择空间一阶至十阶迎风偏斜格式。

给定方波的初始条件:

0,x

u(x)|

t=0

=1,0,

3/32≤x≤9/32,

x>9/32,

(24)

Δt/Δx=计算时选取常数c为1,柯朗数为Cu=c

011,空间步长为Δx=1/128,计算区域[0,1]用周期边界条件,计算1280个周期,状态保持一致。。

u(x)|t=0=exp[-400(x-015)]1(25)

Δt/Δx=计算时选取常数c为1,柯朗数为Cu=c

011,空间步长为Δx=1/128,计算区域[0,1]选用周期边界条件,计算1280个时间步长刚好是一个周期,也就是说在此理想情况下最终状态与初始状态保持一致。计算结果如图4所示。

在图3和图4中,给出了一阶至六阶迎风偏斜

格式在两个不同初始条件下的数值结果。从图3与图4可以看出,一阶格式具有很强的耗散性,而二阶格式的耗散性减小,但是存在一定程度的频散性。对于方波分布,高于二阶精度的高精度迎风偏,相对于频散性来说,,而偶阶精度迎,六阶精度迎风偏斜格式的数值解。而对于高斯分布,六阶迎风偏斜格式的数值解与解析解更为接近。

Takacs[10]给出三种指标来衡量差分格式的优劣,这三种指标分别是总误差,耗散误差,频散误差。总误差由均方差来定义,

(qM∑

-qD),

(26)

这里,qT为准确解,qD为数值解,M为空间网格点的总数。(26)式可以写为

222

σ(qT)σ(qD)+( E=σ(qT)+σ(qD)-2qT- qD),

(27)

图3　方波平流迎风偏斜格式的数值结果:(a)一阶格式;(b)二阶格式;(c)三阶格式;(d)四阶格式;(e)五阶格式;(f)六阶格式

Fig13　Numericalresultsofupwind2biasedschemesfortheadvectionofasquarewave:(a)first2orderscheme;(b)second2orderscheme;(c)third2orderscheme;(d)fourth2orderscheme;(e)fifth2orderscheme;(f)sixth2orderscheme

250

大　气　科　学

ChineseJournalofAtmosphericSciences　31卷

Vol1

图4　高斯平流迎风偏斜格式的数值结果,其余同图3

Fig14　SameasFig13,butfornumericalresultsofupwind2biasedschemesfortheadvectionofaGaussianwave

表1　方波问题的数值误差比较

Table1　Comparisonofnumericalerrorsforasquarewaveprofile

指标

Index

ESP

迎风偏斜格式的阶数Ordernumberoftheupwind2biasedschemes

4163×10-22115×10-21116×10-28153×10-39178×10-38122×10-39169×10-38166×10-39190×10-32153×10-12104×10-24186×10-42170×10-41101×10-47120×10-54168×10-53160×10-52188×10-52130×10-58135×10-22158×10-22110×10-21113×10-28143×10-39170×10-38117×10-39165×10-38163×10-39187×10-31169×10-1

表2　高斯波问题的数值误差比较

Table2　ComparisonofnumericalerrorsforaGaussianwaveprofile

指标

Index

ESP

迎风偏斜格式的阶数Ordernumberoftheupwind2biasedschemes

2143×10-25136×10-31163×10-31150×10-31150×10-31150×10-31150×10-31150×10-31150×10-34120×10-31119×10-21105×10-41179×10-58125×10-81140×10-81150×10-104150×10-118190×10-127140×10-123175×10-41124×10-25126×10-31161×10-31150×10-31150×10-31150×10-31150×10-31150×10-31150×10-33182×10-3

这里,σ为方差,ρ为qT与qD的相关系数,σ为标准差, qT与 qD分别为qT和qD的平均值。(27)式又可写为

(qT)-σ(qD)]2+( E=[σqT- qD)+

)σ(qT)σ(qD),(28)2(1-ρ

这样耗散误差与频散误差分别定义为

(qT)-σ(qD)]2+(qS=[σ ,T-qD)

)σ(qT)σ(qD),P=2(1-ρ

(29)(30)

下面给出以上方波问题与高斯波问题相应的总误

差,耗散误差与频散误差。

2期　No12冯涛等:高精度迎风偏斜格式的比较与分析

FENGTaoetal.AComparisonandAnalysisofHighOrderUpwind2BiasedSchemes

　　　251

在表1和表2中,给出了一阶至十阶迎风偏斜格式在两个不同初始条件下的数值误差。从表1和表2可以看出,高于二阶精度的高精度迎风偏斜格式具有较小的总误差。偶阶精度迎风偏斜格式总误差小于相邻的奇阶精度迎风偏斜格式。相对于频散误差来说,各阶格式的耗散误差是较小的,频散误差的变化趋势与总误差保持一致。并且,随着数值格式精度的提高,数值格式的误差逐渐减小于高于六阶精度的格式来说,大。也就是说,误差的阶数决定,,虽然截断误差逐渐减小,,这样,必然存在着一个数值格式使得总误差达到最小。因此,在同时考虑效率与精度的情况下可以得到六阶精度迎风偏斜格式是最优的格式。

迎风偏斜格式对该方程进行求解。

考虑无粘Burgers方程

(31)+u=0,t在有限区域(-1,1),(I):　u(x,=f--1(x-x

0),(32)(33)u(x-u(x,t)t),t-1

x=x0+ ut

→-∞　(当t→1时)1(34)

　　该方程在(x0- ut,1)出现奇异点。计算中所选各参数为 u=110,x0=-110,Δt=010078125,格点数N=32。计算结果如图5所示。

Bermejo等[12]给出另外一个平滑的初始条件(II):

u(x,0)=

5　无粘Burgers方程试验

Kuo等[11]给出了用内插方法求解的一维无粘Burgers方程的解。该方程在(x0- ut,1)附近出

+sinπx,42

(35)

现较大梯度。为了便于比较,我们采用一阶至十阶

此时t为0166。计算中所选各参数为 u=110,x0=

-110,Δt=01015,格点数N=81。计算结果如图6所示

图5　无粘Burgers方程在t=110的解,其余同图3

Fig15　SameasFig13,butforsolutionsoftheinviscidBurgersπequationatt=110

252

大　气　科　学

ChineseJournalofAtmosphericSciences　31卷

Vol1

图6　无粘Burgers方程在t=0166的解,其余同图3

Fig16　SameasFig13,butforsolutionsoftheinviscidBurgersπequationatt=0166

表3　初始条件(I)的数值误差比较

Table3　Comparisonofnumericalerrorsfortheinitialcondition(I)

指标

Index

ESP

迎风偏斜格式的阶数Ordernumberoftheupwind2biasedschemes

9102×10-32193×10-31141×10-37189×10-48182×10-45117×10-46179×10-44105×10-45173×10-43143×10-49196×10-42146×10-51122×10-62136×10-63144×10-63128×10-64117×10-63178×10-64135×10-64104×10-68102×10-32191×10-31141×10-37187×10-48178×10-45114×10-46175×10-44101×10-45169×10-43139×10-4

表4　初始条件(II)的数值误差比较

Table4　Comparisonofnumericalerrorsfortheinitialcondition(II)

指标

Index

ESP

迎风偏斜格式的阶数Ordernumberoftheupwind2biasedschemes

8123×10-46112×10-55134×10-54107×10-58138×10-51116×10-51167×10-41113×10-42197×10-41187×10-51130×10-43105×10-67103×10-77168×10-72120×10-63112×10-75103×10-74152×10-61115×10-51157×10-76193×10-45182×10-55127×10-54100×10-58116×10-51113×10-51167×10-41109×10-42185×10-41185×10-5

　　下面给出以上两个不同初始条件相对应的总误

差,耗散误差与频散误差。

在图5和图6中,给出了一阶至六阶迎风偏斜格式在两个不同初始条件下的数值结果,在表3和表4中,给出了一阶至十阶迎风偏斜格式在两个不同初始条件下的数值误差。从图5和图6可以看出,各阶格式在突变点误差都比其他点大;而六阶迎风偏斜格在突变点误差最小。这表明六阶格式能

2期　No12冯涛等:高精度迎风偏斜格式的比较与分析

FENGTaoetal.AComparisonandAnalysisofHighOrderUpwind2BiasedSchemes

　　　253

够比较准确地描述具有较大梯度的问题[13]。从表3和表4可以得到,随着数值格式阶数的增加,数值结果也得到了明显改进。而对于高于六阶精度的格式来说,改进的程度并不是很大,因此,在同时考虑精度与效率的情况下,六阶迎风偏斜格式是最优的格式。

参考文献(References)

[1]　CrowleyWP.Numericaladvectionexperiments.Mon.Wea.

Rev.,1968,96:1～11

[2]　LeonardBP.Astableandaccurateconvectivemodelingpro2

cedurebasedoninterpolation.Comput.

MetAp.,19:59～98

[A.forhyperbolicconserva2

ut.Phys.,1983,49(3):357～393

[,纪立人,陈长胜,等.JFNK方法概述及其在大气全

6　结论

程,通过Fourier对低波分量,偶阶精度格式的耗散误差小于相邻的奇阶精度格式,偶阶精度格式与奇阶精度格式都具有较小的频散误差。并且,偶阶精度格式的数值相速度快于实际相速度,而奇阶精度格式的数值相速度慢于实际相速度。对于高波分量,偶阶精度格式比相邻的奇阶精度格式具有较强的耗散效应,偶阶精度格式的频散性小于相邻的奇阶精度格式。奇、偶阶精度格式的行为差异的原因在于用差分替导数时会带来截断误差,截断误差中的偶阶导数产生数值耗散(数值粘性),奇阶导数产生数值频散(色散)。为了检验这些格式的计算性能,在一维问题上进行了应用。首先,考虑恒定风场条件下的一维线性平流试验。主要选择两种不同的初始条件来评价数值格式的精度,这两种试验问题分别是高斯函数、方波函数。对于方波分布,高于二阶精度的高精度迎风偏斜格式的耗散误差明显减小,同时存在着一定程度的数值频散。而对于高斯分布,高精度迎风偏斜格式的耗散性与频散性得到明显改善。数值结果表明,随着数值格式精度的提高,数值格式的误差逐渐减小。而对于高于六阶精度的格式来说,改进的程度并不是很大。其次,应用各阶格式到具有两种不同初始条件的无粘Burgers方程。试验结果表明,随着数值格式阶数的增加,数值结果也得到了明显改进。而对于高于六阶精度的格式来说,进一步的变化并不明显。也就是说,数值格式的性能并不能完全由截断误差的阶数决定,随着数值格式精度的提高,虽然截断误差逐渐减小,但是舍入误差却逐渐增加,这样,必然存在着一个数值格式使得总误差达到最小。总之,在兼顾效率与精度条件下六阶迎风偏斜格式是最好的。

隐式非静力模式中的应用方案.大气科学,2006,30(5):

821～833

ChenJiabin,JiLiren,ChenChangsheng,etal.AreviewofJFNKmethodsanditsapplicationsinatmosphericnon2hydro2staticalmodel.ChineseJournalofAtmosphericSciences(inChinese),2006,30(5):821～833

[5]　RaiMM,MoinP.Directsimulationofturbulentflowusingfinite

differenceschemes.J.Comput.Phys.,1991,96(1):15～53

[6]　TrembackCJ,PowellJ,CottonWR.Theforward2in2time

upstreamadvectionscheme:Extensiontohigherorders.

Mon.Wea.Rev.,1986,115:540～555

[7]　LiJP.Generalexplicitdifferenceformulasfornumericaldif2

ferentiation.J.Comput.Appl.Math.,2005,183(1):29

～52,doi:10.1016/j.cam.2004.12.026

[8]　马延文,傅德薰.高精度有限差分法与复杂流动的数值模

拟.自然科学进展,2002,12(8):785～793

MaYanwen,FuDexun.

Highorderaccuratedifference

schemeandnumericalsimulationofcomplexflows.Progress

inNaturalScience(inChinese),2002,12(8):785～793

[9]　傅德薰,马延文.物理问题的数值模拟及高精度差分格式.

计算物理,1992,9(4):501～505

FuDexun,MaYanwen.Numericalsimulationofphysicalquestionandhighorderaccuratedifferencescheme.Chinese

JournalofComputationalPhysics(inChinese),1992,9

(4):501～505

[10]　TakacsLL.Atwo2stepschemefortheadvectionequation

withminimizeddissipationanddispersionerrors.Mon.Wea.

Rev.,1985,113:1050～1065

[11]　KuoH2C,WilliansRT.Semi2Lagrangiansolutiontotheinviscid

Burgersπequation.Mon.Wea.Rev.,1990,118:1278～1288

[12]　BermejoR,StaniforthA.Theconversionofsemi2Lagrangian

advectionschemestoquasi2monotoneschemes.Mon.Wea.

Rev.,1992,120:2622～2632

[13]　王军,陈嘉滨.完全非内插半拉格朗日格式在一维Burgers

方程及二维浅水方程上的应用.大气科学,2000,24(4):

493～508

WangJun,ChenJiabin.Anapplicationofnoninterpolatingsemi2Lagrangianschemetoone2dimensionalBurgersequationandtwo2dimensionalshallowwaterequations.ChineseJournalofAtmos2

phericSciences(inChinese),2000,24(4):493～508

第31卷第2期2007年3月

大　气　科　学

ChineseJournalofAtmosphericSciencesVol131　No12

Mar.2007

高精度迎风偏斜格式的比较与分析

冯涛1,2　李建平1

1　1000292中国科学院研究生院,北京　100039

文章编号　10069895(2007)02024509　　　中图分类号　P435　　　文献标识码　A

AComparisonandAnalysisofHighOrderUpwind2BiasedSchemes

FENGTao1,2andLIJian2Ping1

1StateKeyLaboratoryofNumericalModelingforAtmosphericSciencesandGeophysicalFluidDynamics,InstituteofAtmos2

phericPhysics,ChineseAcademyofSciences,Beijing　100029

2GraduateUniversityoftheChineseAcademyofSciences,Beijing　100049

246

大　气　科　学

ChineseJournalofAtmosphericSciences　31卷

Vol1

Keywords　explicitfinitedifferenceformulas,upwind2biasedscheme,advectionequation

1　引言

点上的值通过Lagrange多项式插值得到。最近,Li[7]给出了具有任意阶精度的一般显式有限差分公

式,此公式是以泰勒展开式为基础,利用范德蒙行

2　高精度有限差分方法

考虑如下一维线性平流方程及对应的半离散化方程:

+c=0,tx

+c=0,

Δxt

(1)(2)

式中,u=u(x,t),c取为常数,本文c选为大于

零,Fj/Δx为5u/5x的差分逼近式。

2期　No12冯涛等:高精度迎风偏斜格式的比较与分析

FENGTaoetal.AComparisonandAnalysisofHighOrderUpwind2BiasedSchemes

　　　247

(m)

(xi)=

j=0

∑

(m)(m)n-m+1

),(3)dn+1,i,jf(xj)+On,j(h

(xj)=f′

(

fj-2-4fj-1+3fj),2h

(11)

其中,xi=x0+ih(i=0,1,…,n,h≠0),f是所选定点的函数。(1)

d2,0,1=1,　　n=1,

(4)(1)

d2,1,0=-1,　　n=1,

(m-1)

(m-j)(m)

dn+1,i,j=,　j≠i,n>1,(5)

j!(n-j)!

(m)

n+1,i,j

三阶迎风偏斜格式:

(xj)=(ff′

j-2

-6f

j-1

3fj+2fj+1),(12)

(j)f′

(

-j-3+fj-2-18fj-1+10fj+h

(13)

j1,

(dn+1,i,i=-(

)

(m)

j=0,j≠i

∑d

,　n≥1,

(xj)=f′

m-1

an-1,i,j=am-1(-i,ki,i),

　　　　　　　　(k,…n,　k≠i,j)(0)

an=a0(x1,x2,…,xn)=x1x2…xn,

(1)n

(7)

(

-2fj-3+15fj-2-60fj-1+60h

(14)

　　20fj+30fj+1-3fj+2),六阶迎风偏斜格式:

(xj)=f′

=a1(x1,x2,…,xn)=x1x2…xn-1+=a2(x1,x2,…,xn)=x1x2…xn-2+…+

　　x1x2…xn-2xn+…+x1x3…xn+x2x3…xn,

(2)

(

fj-4-8fj-3+30fj-2-80fj-1+60h

(15)

　　35fj+24fj+1-2fj+2),七阶迎风偏斜格式:

(xj)=f′

　　x3x4…xn,…

(n-2)

=an-2(x1,x2,…,xn)=x1x2+

=an-1(x1,x2,…,xn)=x1+x2+…+xn,=an(x1,x2,…,xn)=1,

(8)

(

3fj-4-28fj-3+126fj-2-420h

(16)

　　x2x3+…+xn-1xn,

anan

(n-1)(n)

　　420fj-1+105fj+252fj+1-42fj+2+4fj+3),八阶迎风偏斜格式:

(xj)=f′

余项

(m)

n,j

(

-3fj-5+30fj-4-140fj-3+

840h

(17)

(n-m+1)

n-mn-m+1

)=・

(n+1)!

(n+1)

(m-1)n+1

(-1)jj!(n-j)!

　　420fj-2-1050fj-1+378fj+420fj+1-　　60fj+2+5fj+3),

(9)

j=0,j≠i

∑

九阶迎风偏斜格式:

(xj)=f′

　　可根据精度的需要,选择适当个数的网格基架

:一阶迎风偏斜格式:

(xj)=(-fj-1+fj),(10)f′

(

-4fj-5+45fj-4-240fj-3+

2520h

(18)

　　840fj-2-2520fj-1+504fj+1680fj+1-　　360fj+2+60fj+3-5fj+4),十阶迎风偏斜格式:

(xj)=f′

(

2fj-6-24fj-5+135fj-4-2520h

(19)

　　480fj-3+1260fj-2-3024fj-1+924fj+　　1440fj+1-270fj+2+40fj+3-3fj+4)1

3　数值格式的Fourier分析

二阶迎风偏斜格式:

248

大　气　科　学

ChineseJournalofAtmosphericSciences　31卷

Vol131

尺度量离散导数的逼近精度。在单纯进行精度分析时并不要求被离散的函数是线性的。现在考虑线性问题,取初值为u(x,0)=exp(ikx),方程(1)的准确解为

　u(x,t)=exp[

ik(x-ct)],准确解为　u(xj,t)=exp

-ck

texpa

ikxjta

,(20)

(21)

4　线性平流试验

为了进一步比较一阶至十阶迎风偏斜格式的性能,我们以方波和高斯波的初值问题为例进行计算。

时间导数的差分逼近采用三阶精度的R2K三步方法[9]。设L(u)为5u/5x的差分逼近式,则对于方程(1),三阶精度的三步方法有以下形式

uuu

(1)

Δx,h,p与数值波的传其中a=k

播速度有关h与p的分

析表达式来。对照式(20)和式(21)可知,对取定波数k,数值解逼近准确解,要求

Δx),　h/a→0,p/a→1　(a=k

(22)

它们的逼近程度反映了格式的精度。下面给出一阶至十阶迎风偏斜格式h与p随a的变化情况。

在图1和图2中,给出了一阶至十阶迎风偏斜格式的频散误差和耗散误差,其中曲线1、2、3、4、5、6、7、8、9、10依次对应为一阶格式、二阶格式、

三阶格式、四阶格式、五阶格式、六阶格式、七阶格式、八阶格式、九阶格式、十阶格式,曲线0为p

Δx曲线,它对应于微分方程的准确解。从=a=k

=u==

(0)

ΔtL(u(0)),+c

(23)

)

(0)(1)(1)

u+u+cΔtL(u

),444(0)(2)(2)

u+u+cΔtL(u),333

(3)

图1　一阶至十阶迎风偏斜格式的频散误差

Fig11　Thedispersionerrorsofupwind2biasedschemesfororders1through10

图2　一阶至十阶迎风偏斜格式的耗散误差

Fig12　Thedissipationerrorsofupwind2biasedschemesfororders1through10

2期　No12冯涛等:高精度迎风偏斜格式的比较与分析

FENGTaoetal.AComparisonandAnalysisofHighOrderUpwind2BiasedSchemes

　　　249

空间导数的差分逼近选择空间一阶至十阶迎风偏斜格式。

给定方波的初始条件:

0,x

u(x)|

t=0

=1,0,

3/32≤x≤9/32,

x>9/32,

(24)

Δt/Δx=计算时选取常数c为1,柯朗数为Cu=c

011,空间步长为Δx=1/128,计算区域[0,1]用周期边界条件,计算1280个周期,状态保持一致。。

u(x)|t=0=exp[-400(x-015)]1(25)

Δt/Δx=计算时选取常数c为1,柯朗数为Cu=c

在图3和图4中,给出了一阶至六阶迎风偏斜

Takacs[10]给出三种指标来衡量差分格式的优劣,这三种指标分别是总误差,耗散误差,频散误差。总误差由均方差来定义,

(qM∑

-qD),

(26)

这里,qT为准确解,qD为数值解,M为空间网格点的总数。(26)式可以写为

222

σ(qT)σ(qD)+( E=σ(qT)+σ(qD)-2qT- qD),

(27)

图3　方波平流迎风偏斜格式的数值结果:(a)一阶格式;(b)二阶格式;(c)三阶格式;(d)四阶格式;(e)五阶格式;(f)六阶格式

250

大　气　科　学

ChineseJournalofAtmosphericSciences　31卷

Vol1

图4　高斯平流迎风偏斜格式的数值结果,其余同图3

Fig14　SameasFig13,butfornumericalresultsofupwind2biasedschemesfortheadvectionofaGaussianwave

表1　方波问题的数值误差比较

Table1　Comparisonofnumericalerrorsforasquarewaveprofile

指标

Index

ESP

迎风偏斜格式的阶数Ordernumberoftheupwind2biasedschemes

表2　高斯波问题的数值误差比较

Table2　ComparisonofnumericalerrorsforaGaussianwaveprofile

指标

Index

ESP

迎风偏斜格式的阶数Ordernumberoftheupwind2biasedschemes

这里,σ为方差,ρ为qT与qD的相关系数,σ为标准差, qT与 qD分别为qT和qD的平均值。(27)式又可写为

(qT)-σ(qD)]2+( E=[σqT- qD)+

)σ(qT)σ(qD),(28)2(1-ρ

这样耗散误差与频散误差分别定义为

(qT)-σ(qD)]2+(qS=[σ ,T-qD)

)σ(qT)σ(qD),P=2(1-ρ

(29)(30)

下面给出以上方波问题与高斯波问题相应的总误

差,耗散误差与频散误差。

2期　No12冯涛等:高精度迎风偏斜格式的比较与分析

FENGTaoetal.AComparisonandAnalysisofHighOrderUpwind2BiasedSchemes

　　　251

迎风偏斜格式对该方程进行求解。

考虑无粘Burgers方程

(31)+u=0,t在有限区域(-1,1),(I):　u(x,=f--1(x-x

0),(32)(33)u(x-u(x,t)t),t-1

x=x0+ ut

→-∞　(当t→1时)1(34)

　　该方程在(x0- ut,1)出现奇异点。计算中所选各参数为 u=110,x0=-110,Δt=010078125,格点数N=32。计算结果如图5所示。

Bermejo等[12]给出另外一个平滑的初始条件(II):

u(x,0)=

5　无粘Burgers方程试验

Kuo等[11]给出了用内插方法求解的一维无粘Burgers方程的解。该方程在(x0- ut,1)附近出

+sinπx,42

(35)

现较大梯度。为了便于比较,我们采用一阶至十阶

此时t为0166。计算中所选各参数为 u=110,x0=

-110,Δt=01015,格点数N=81。计算结果如图6所示

图5　无粘Burgers方程在t=110的解,其余同图3

Fig15　SameasFig13,butforsolutionsoftheinviscidBurgersπequationatt=110

252

大　气　科　学

ChineseJournalofAtmosphericSciences　31卷

Vol1

图6　无粘Burgers方程在t=0166的解,其余同图3

Fig16　SameasFig13,butforsolutionsoftheinviscidBurgersπequationatt=0166

表3　初始条件(I)的数值误差比较

Table3　Comparisonofnumericalerrorsfortheinitialcondition(I)

指标

Index

ESP

迎风偏斜格式的阶数Ordernumberoftheupwind2biasedschemes

表4　初始条件(II)的数值误差比较

Table4　Comparisonofnumericalerrorsfortheinitialcondition(II)

指标

Index

ESP

迎风偏斜格式的阶数Ordernumberoftheupwind2biasedschemes

　　下面给出以上两个不同初始条件相对应的总误

差,耗散误差与频散误差。

2期　No12冯涛等:高精度迎风偏斜格式的比较与分析

FENGTaoetal.AComparisonandAnalysisofHighOrderUpwind2BiasedSchemes

　　　253

参考文献(References)

[1]　CrowleyWP.Numericaladvectionexperiments.Mon.Wea.

Rev.,1968,96:1～11

[2]　LeonardBP.Astableandaccurateconvectivemodelingpro2

cedurebasedoninterpolation.Comput.

MetAp.,19:59～98

[A.forhyperbolicconserva2

ut.Phys.,1983,49(3):357～393

[,纪立人,陈长胜,等.JFNK方法概述及其在大气全

6　结论

隐式非静力模式中的应用方案.大气科学,2006,30(5):

821～833

ChenJiabin,JiLiren,ChenChangsheng,etal.AreviewofJFNKmethodsanditsapplicationsinatmosphericnon2hydro2staticalmodel.ChineseJournalofAtmosphericSciences(inChinese),2006,30(5):821～833

[5]　RaiMM,MoinP.Directsimulationofturbulentflowusingfinite

differenceschemes.J.Comput.Phys.,1991,96(1):15～53

[6]　TrembackCJ,PowellJ,CottonWR.Theforward2in2time

upstreamadvectionscheme:Extensiontohigherorders.

Mon.Wea.Rev.,1986,115:540～555

[7]　LiJP.Generalexplicitdifferenceformulasfornumericaldif2

ferentiation.J.Comput.Appl.Math.,2005,183(1):29

～52,doi:10.1016/j.cam.2004.12.026

[8]　马延文,傅德薰.高精度有限差分法与复杂流动的数值模

拟.自然科学进展,2002,12(8):785～793

MaYanwen,FuDexun.

Highorderaccuratedifference

schemeandnumericalsimulationofcomplexflows.Progress

inNaturalScience(inChinese),2002,12(8):785～793

[9]　傅德薰,马延文.物理问题的数值模拟及高精度差分格式.

计算物理,1992,9(4):501～505

FuDexun,MaYanwen.Numericalsimulationofphysicalquestionandhighorderaccuratedifferencescheme.Chinese

JournalofComputationalPhysics(inChinese),1992,9

(4):501～505

[10]　TakacsLL.Atwo2stepschemefortheadvectionequation

withminimizeddissipationanddispersionerrors.Mon.Wea.

Rev.,1985,113:1050～1065

[11]　KuoH2C,WilliansRT.Semi2Lagrangiansolutiontotheinviscid

Burgersπequation.Mon.Wea.Rev.,1990,118:1278～1288

[12]　BermejoR,StaniforthA.Theconversionofsemi2Lagrangian

advectionschemestoquasi2monotoneschemes.Mon.Wea.

Rev.,1992,120:2622～2632

[13]　王军,陈嘉滨.完全非内插半拉格朗日格式在一维Burgers

方程及二维浅水方程上的应用.大气科学,2000,24(4):

493～508

WangJun,ChenJiabin.Anapplicationofnoninterpolatingsemi2Lagrangianschemetoone2dimensionalBurgersequationandtwo2dimensionalshallowwaterequations.ChineseJournalofAtmos2

phericSciences(inChinese),2000,24(4):493～508

高精度迎风偏斜格式的比较与分析

相关文章