王秋菊，任玉新

（清華大學(xué) 航天航空學(xué)院，北京 100084）

關(guān)于非均勻網(wǎng)格高精度格式的誤差分析方法

王秋菊，任玉新

（清華大學(xué) 航天航空學(xué)院，北京 100084）

以三次樣條重構(gòu)有限體積方法為例，研究非均勻網(wǎng)格上截斷誤差的分析方法。通過推導(dǎo)得到了分析非均勻網(wǎng)格上截斷誤差的基本準(zhǔn)則，即在非均勻網(wǎng)格的截斷誤差分析中，要保證不顯式或者隱含地改變數(shù)值方法對應(yīng)的模板點——當(dāng)不滿足這一準(zhǔn)則時，誤差分析會得到不自洽的結(jié)果；而滿足這一準(zhǔn)則時，可保證分析結(jié)果的正確性。利用正確的誤差分析結(jié)果，可發(fā)展進(jìn)一步提高計算精度的措施。據(jù)此發(fā)展了擴散項在非均勻網(wǎng)格上達(dá)到三階精度計算方法，從而可以使對流項和擴散項的計算達(dá)到一致三階精度。

三次樣條；誤差分析；有限體積方法

0 引 言

樣條插值方法是Schoenberg［1］在20世紀(jì)40年代提出的。三次樣條插值函數(shù)具有非常良好的特性，可用于求解偏微分方程［2-3］。Rubin［4］總結(jié)了樣條函數(shù)在數(shù)值求解偏微分方程方面所具有的優(yōu)點，包括只需求解三對角方程，計算效率高；對函數(shù)值、一階導(dǎo)數(shù)和二階導(dǎo)數(shù)都具有良好的逼近效果；邊界條件處理非常方便等。許多學(xué)者將樣條插值方法應(yīng)用到非線性雙曲守恒方程的求解，如，McCartin和Jameson［5］引入指數(shù)樣條來處理非光滑的數(shù)值解，并采用文獻(xiàn)［6-9］提出的方法計算張力因子；王建立等［10］研究了基于樣條逼近有限體積法的通量分裂算法。

緊致格式具有較好的譜特性能分辨流場中的小尺度結(jié)構(gòu)［11］；Hermite插值方法在非均勻網(wǎng)格上容易達(dá)到高階精度［12-15］。樣條插值可以看做基于Hermite插值的某種緊致格式。因此，基于樣條插值的算法有望在非均勻網(wǎng)格上實現(xiàn)高精度、高分辨率。三次樣條插值方法可以直接應(yīng)用于非均勻網(wǎng)格，無需將非均勻網(wǎng)格通過坐標(biāo)變換至均勻網(wǎng)格的計算空間。在非均勻網(wǎng)格上，三次樣條插值對對流-擴散方程的對流項和擴散項的離散可分別達(dá)到三階精度和二階精度。在發(fā)展求解任意曲線坐標(biāo)非均勻網(wǎng)格可壓縮流動的高精度樣條逼近有限體積方法［16］的過程中發(fā)現(xiàn)，三次樣條插值的截斷誤差具有局部特性，不恰當(dāng)?shù)耐茖?dǎo)有可能導(dǎo)致錯誤結(jié)果。因此，提出了在非均勻網(wǎng)格上，分析截斷誤差的基本準(zhǔn)則，即不顯式或者隱含地改變數(shù)值方法對應(yīng)的模板點。如果不滿足這一準(zhǔn)則，將得到不自洽的截斷誤差分析結(jié)果。在滿足這一準(zhǔn)則的前提下，得出了非均勻網(wǎng)格上用三次樣條重構(gòu)擴散項所產(chǎn)生的截斷誤差，并利用該截斷誤差分析結(jié)果，發(fā)展了擴散項三階精度離散的計算方法。從而使三次樣條重構(gòu)在非均勻網(wǎng)格上對對流項和擴散項的離散達(dá)到一致三階精度。本文的精度分析所遇到的問題及其解決方案，在非均勻網(wǎng)格高精度數(shù)值方法的構(gòu)造中，具有普遍性，因此對此進(jìn)行研究是有意義的。

本文將在第2節(jié)中給出三次樣條重構(gòu)有限體積方法離散對流擴散方程的具體過程；在第3節(jié)中，首先比較截斷誤差的兩種分析方法，通過具體推導(dǎo)過程分析對比得出非均勻網(wǎng)格上分析截斷誤差的基本準(zhǔn)則；在第4節(jié)中，將利用截斷誤差分析，給出擴散項三階精度的離散方法。本文第5節(jié)將給出結(jié)論。

1 三次樣條重構(gòu)方法

采用一維標(biāo)量線性對流擴散方程來說明有限體積框架下樣條重構(gòu)的構(gòu)造方法，方程為：

其中，ūj是單元平均值，數(shù)值通量為：

其中和分別是因變量及其一階導(dǎo)數(shù)在控制體界面的近似值，通常由某種重構(gòu)方法得到。

本文的有限體積重構(gòu)方法，采用基于原函數(shù)的三次樣條重構(gòu)。已知某時刻u（x）的平均值，定義

u（x）的原函數(shù)定義為：

單元界面xj＋1／2處的原函數(shù)的值為：

因此單元界面處原函數(shù)的值可以根據(jù)單元平均值計算，三次樣條函數(shù)S（x）來插值原函數(shù)W（x），插值條件為：

對于三次樣條多項式，有4N（N為網(wǎng)格節(jié)點數(shù)）個待定系數(shù)，根據(jù)式（8）～式（10）有3（N-1）個關(guān)系式，再加上式（7）的插值條件，共有4N-2個關(guān)系，因此另外需要兩個條件就可以確定待定系數(shù)。通常在求解域的兩個端點［x1/2，xN－1/2］上各加一個邊界條件，邊界條件根據(jù)實際問題確定。

在實際實施中，有兩種等價的方法可以進(jìn)行三次樣條重構(gòu)。第一種是求解所謂“三轉(zhuǎn)角”方程，即：

其中，N是網(wǎng)格節(jié)點數(shù)，并且

其中mj＋1／2是原函數(shù)一階導(dǎo)數(shù)S′（xj＋1／2）的近似，也可認(rèn)為是因變量在控制體邊界處的函數(shù)值，即：

因變量的一階導(dǎo)數(shù)采用下面公式計算，即：

第二種方法求解“三彎矩”方程，即：

Mj＋1／2是因變量一階導(dǎo)數(shù)的近似，即式（13）。因變量本身在控制體界面處的值可用式（12）計算，其中：

把式（12）、式（13）代入式（3）即可得到數(shù)值通量。這種方法可以推廣到求解Euler和Navier-Stokes方程；也可以應(yīng)用于多維曲線坐標(biāo)非均勻網(wǎng)格的情況。關(guān)于格式在上述情況的具體形式，以及穩(wěn)定性、激波捕捉等問題的分析和處理，請參見文獻(xiàn)［16］。本文主要研究在非均勻網(wǎng)格下誤差分析的問題。

2 誤差分析方法

根據(jù)第2節(jié)，在進(jìn)行重構(gòu)時，我們既可以通過求解三轉(zhuǎn)角方程計算，也可以根據(jù)三彎矩方程計算。這兩種方法在邊界條件相容時，是完全等價的。所以，下面只分析求解三轉(zhuǎn)角方程的情形。首先研究m的精度。將式（11）等式兩邊在j＋1／2處進(jìn)行泰勒級數(shù)展開。式（11）等式左邊mj－1／2，mj＋3／2的泰勒級數(shù)展開式分別為：

式（11）等式右邊，原函數(shù)的泰勒級數(shù)展開式分別為：

將式（17）～（20）代入到式（11）中，則得到等式左邊為：

等式右邊為：

在求解三轉(zhuǎn)角方程的條件下，二階導(dǎo)數(shù)M可根據(jù)式（14）或式（15）計算。為了便于區(qū)分，式（14）及式（15）分別寫為：

從數(shù)值上，上述兩個公式得到的計算結(jié)果完全相同，因此，這兩種方法的誤差也應(yīng)該是相同的。根據(jù)式（24）和式（25）對M的誤差分析有兩種方法，下面將給出這兩種方法，并對其分析過程進(jìn)行對比。

2.1 分析方法1

根據(jù)式（23）一階導(dǎo)數(shù)mj＋1／2的表達(dá)式，可推知mj－1／2和mj＋3／2的表達(dá)式分別為：

和

把上面兩個公式分別在j＋1／2處做泰勒級數(shù)展開：

同時，原函數(shù)Wj－1／2和Wj＋3／2在j＋1／2處的泰勒級數(shù)展開式為式（19）～式（20）。將式（23）、式（28）～式（29）和（19）～式（20）代入到式（24）、式（25），有

2.2 分析方法2

將式（11）左側(cè)在j＋1／2處做泰勒級數(shù)展開，可得式（23）。同理，我們也可以把式（11）左側(cè)在j－1／2和j＋3／2處展開。當(dāng)式（11）左側(cè)在j－1／2處展開時，等式左邊mj＋1／2、mj＋3／2的泰勒級數(shù)展開式分別為：

式（11）右邊，Wj－1／2和Wj＋3／2在j＋1／2處泰勒級數(shù)展開式為式（19）～式（20）。將式（19）～式（20）和式（32）～式（33）代入到式（11）中，則得到等式左邊為：

同理，將式（11）等式左邊在j＋3／2處進(jìn)行泰勒級數(shù)展開，等式右邊在j＋1／2處進(jìn)行泰勒級數(shù)展開。式（11）等式左邊mj－1／2，mj＋1／2的泰勒級數(shù)展開式分別為：

式（11）右邊，Wj－1／2和Wj＋3／2在j＋1／2處泰勒級數(shù)展開式為式（19）～式（20）。將式（19）～式（20）和式（36）～式（37）代入到式（11）中，則得到等式左邊為：

將式（19）、式（23）和式（35）代入到式（24），有：

將式（20），式（23）和式（39）代入到式（25），則有：

對比式（40）、式（41）中的和可以發(fā)現(xiàn)，二者余項完全相等。這個結(jié)果與和在數(shù)值上恒等是自洽的，因而是正確的。

2.3 兩種方法的比較分析

方法1和方法2的不同之處在于：方法1中mj－1／2和mj＋3／2在j＋1／2處的泰勒級數(shù)展開式是在mj＋1／2（式（23））的基礎(chǔ)上分別向左和向右移動一個單元指標(biāo)。方法2中mj－1／2和mj＋3／2在j＋1／2處的泰勒級數(shù)展開式是在m關(guān)系式（式（11））的基礎(chǔ)上，通過循環(huán)消去法得到。

為什么方法1不正確而方法2正確呢？這是因為，非均勻網(wǎng)格上的數(shù)值方法與均勻網(wǎng)格數(shù)值方法不同，其誤差估計具有局部特性。例如，式（40）或者式（41）表明Mj＋1／2的大小與局部空間步長hj、hj＋1有關(guān)。如果改變了模板點，則會引入其他模板的長度尺度；而又由于非均勻網(wǎng)格不同單元（或者控制體）的尺度互不相關(guān)，因而采用不同模板，會得到不同的結(jié)果；從而造成分析的不自洽。因此，我們可以得到非均勻網(wǎng)格誤差分析的一個重要原則，即在非均勻網(wǎng)格的誤差分析中，只有保證不顯式或者隱含地改變數(shù)值方法對應(yīng)的模板點，才能得到正確結(jié)果。

為了進(jìn)一步說明這個問題，下面回顧一下三次樣條重構(gòu)中三轉(zhuǎn)角方程的推導(dǎo)。我們知道在控制體Ij＝［xj－1／2，xj＋1／2］上，三次樣條函數(shù)S（x）可以表示為

其中，αj±1／2（x）和βj±1／2（x）是三次Hermite插值基函數(shù)，分別為：

和

根據(jù)式（42），可得式（24）；同理根據(jù)Ij＋1＝［xj＋1／2，xj＋3／2］上的三次樣條函數(shù)S（x）的表達(dá)式，可得式（25）。樣條函數(shù)要求＝，由式（24）、式（25）即可得到三轉(zhuǎn)角方程式（11）。由此可知，式（11）的模板為單元Ij，Ij＋1，因此在進(jìn)行精度分析時，不能采用在這兩個單元以外定義的關(guān)系。

方法1在推導(dǎo)過程中，使用了式（23）、式（26）和式（27），且式（26）和式（27）由式（23）平移得到。由于式（23）的模板為單元Ij，Ij＋1，則式（26）和式（27）對應(yīng)的模板分別為Ij－1，Ij和Ij＋1，Ij＋2。也就是說，在推導(dǎo)出式（30）和式（31）的過程中，二者的模板是不同的，因而余項的結(jié)果也不同。在方法2中，所有的推導(dǎo)基于式（11），保證了模板為單元Ij，Ij＋1。這與三轉(zhuǎn)角方程的推導(dǎo)相一致，從而可以得到正確的結(jié)果。另外，我們注意到，在均勻網(wǎng)格上，兩種方法得到的結(jié)果是相同的。

3 擴散項的三階精度離散方法

首先，計算Mj－1／2和Mj＋3／2在j＋1／2處的泰勒級數(shù)展開式，

根據(jù)式（45）和式（46）可以計算得到：

將式（48）～式（50）代入到式（47）中，得到：

將式（51）代入到式（40）得到：

令：

則有：

4 結(jié) 論

本文研究了有限體積方法框架下，在非均勻網(wǎng)格上采用三次樣條重構(gòu)方法離散方程對流項和擴散項所產(chǎn)生的截斷誤差的分析方法。通過具體推導(dǎo)過程的比較分析，提出了非均勻網(wǎng)格上截斷誤差分析的基本準(zhǔn)則，即保證不顯式或隱含地改變數(shù)值方法對應(yīng)的模板點。若不滿足該準(zhǔn)則，將得到不自洽的誤差分析結(jié)果。本文在該準(zhǔn)則的基礎(chǔ)上得出了對擴散項離散的截斷誤差，提出了提高擴散項離散精度的方法，發(fā)展了在非均勻網(wǎng)格上擴散項三階精度離散的計算方法。從而使三次樣條重構(gòu)在非均勻網(wǎng)格上對方程對流項和擴散項的離散達(dá)到一致三階精度。

［1］SCHOENBERG I J．Cardinal spline interpolation［M］．In：CBMSNSE region conference series in applied mathematics 12．Philadelphia，PA：Society for industrial and applied mathematics，1973．

［2］AHLBERG J，NILSON E，WALSH J．The theory of splines and their applications，academic press［M］．New York：1967．

［3］RUBIN S G，GRAVES R A．Cubic spline approximation for problems in fluid mechanics［R］．NASA TR R-436，1975．

［4］RUBIN S G，KHOSLA P K．Higher-order numerical solutions using cubic splines［J］．AIAA Journal，1976，14（7）：

［5］MCCARTIN B J，JAMESON A．Numerical solution of nonlinear hyperbolic conservation laws using exponential splines［J］．Computational Mechanics，1990，6；77-91．

［6］MCCARTIN B J．Theory，computation，and application of exponential splines［R］．DOE／ER／03077-171（1981）．

［7］MCCARTIN B J．Applications of exponential splines in computational fluiddynamics［J］．AIAA J．，1983，21：1059-1065．

［8］MCCARTIN B J．Theory of exponential splines［J］．J．Approx．Theory，1991，66（1）：1-23．

［9］MCCARTIN B J．Computation of exponential splines［J］．SIAM J．Sci．Stat．Comp．，1990，11（2）：242-262．

［10］王建立，任玉新，沈孟育．基于樣條逼近有限體積法的通量分裂算法［J］．計算物理，2000，17（4）：421-425．

［11］LELE S K．Compact finite-difference schemes with spectral-like resolution［J］．Journal of Computational Physics，1992，103（1）：16-42．

［12］ZHONG X，TATINENI M．High-order non-uniform grid schemes for numerical simulation of hypersonic boundary-layer stability and transition［J］．J．Comput．Phys．，2003，190：419-458．

［13］MAHESH K．A family of high order finitedifference schemes with good spectral resolution［J］．J．Comput．Phys．，1998，145：332-358．

［14］CHU P C，EAN C．A three-point combined compactdifference scheme［J］．J．Comput．Phys．，1998，140：370．

［15］CHU P C，EAN C．A three-point sixth-order nonuniform combined compactdifference［J］．J．Comput．Phys．，1999，148：663．

［16］WANG Q J，REN Y X，An accurate and robust finite volume scheme based on the spline interpolation for solving the Euler and Navier-Stokes equations on non-uniform curvilinear grids［J］．Journal of Computational Physics，（submitted）．

Error analysis of high order schemes on non-uniform grids

WANG Qiuju，REN Yuxin
（School of Aerospace，Tsinghua University，Beijing 100084，China）

The derivation of truncation error on non-uniform grid for cubic spline reconstruction in finite volume framework is analyzed in this paper．The basic rule for thederivation of truncation error on non-uniform grid is proposed which is that allderivations in the truncation error analysis should base on the same stencil．If this rule is not satisfied，the truncation error will not be self-consistent．Eurthermore，the thirdorderdiscretization of thediffusion terms is obtained based on the correct truncation error analysis which makes both the convection anddiffusion terms to achieve consistent third-order accuracy．

cubic spline；truncation error analysis；finite volume method

V211.3

Adoi：10.7638／kqdlxxb-2014.0108

0258-1825（2014）06-0741-07

2014-09-16；

2014-10-11

國家自然科學(xué)基金（11172153）

任玉新，教授，研究方向：非定常流動的數(shù)值模擬．E-mail:ryx＠tsinghua.edu.cn

王秋菊，任玉新．關(guān)于非均勻網(wǎng)格高精度格式的誤差分析方法［J］．空氣動力學(xué)學(xué)報，2014，32（6）：741-747．

10．7638／kqdlxxb-2014．0108 WANG Q J，REN Y X．Error analysis of high order schemes on non-uniform grids［J］．ACTA Aerodynamica Sinica，2014，32（6）：741-747．