基于Richardson插值法的CFD驗(yàn)證和確認(rèn)方法的研究

陳堅(jiān)強(qiáng)， 張益榮

（中國(guó)空氣動(dòng)力研究與發(fā)展中心，四川 綿陽(yáng)621000）

0 引 言

隨著計(jì)算流體力學(xué)（CFD）在航空航天等領(lǐng)域的廣泛應(yīng)用，飛行器設(shè)計(jì)對(duì)CFD所提供結(jié)果的可信度要求越來越高。但由于絕大多數(shù)情況下問題的“真值”未知，而且數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果往往存在差異，這就要求確定數(shù)值計(jì)算結(jié)果的不確定度或者是誤差，而驗(yàn)證和確認(rèn)方法正是評(píng)價(jià)數(shù)值解精度和可信度的主要手段。

Boehm［1］和 Blottner［2］給出了驗(yàn)證和確認(rèn)的簡(jiǎn)單定義：驗(yàn)證過程即是否正確地解方程，確認(rèn)過程就是是否求解正確的方程。驗(yàn)證一般采用精確解比較方法、制造解比較方法以及網(wǎng)格收斂性研究等方法進(jìn)行。其中前兩種只對(duì)簡(jiǎn)單問題成立，對(duì)于較為復(fù)雜的問題，網(wǎng)格收斂性研究是數(shù)值計(jì)算驗(yàn)證的最有效手段。確認(rèn)的基本內(nèi)容是指出和量化在概念模型和計(jì)算模型中的誤差和不確定度，量化數(shù)值解的數(shù)值誤差，評(píng)估實(shí)驗(yàn)的不確定度，最后進(jìn)行計(jì)算結(jié)果和實(shí)驗(yàn)數(shù)據(jù)的比較。

本文通過采取多套不同尺度但拓?fù)湟恢碌木W(wǎng)格對(duì)T2－97模型的計(jì)算，利用Richardson插值法，得到“精確解”與離散誤差，主要考察不完全迭代收斂誤差和空間離散誤差，通過網(wǎng)格收斂性研究完成了CFD的驗(yàn)證，并繼而與實(shí)驗(yàn)數(shù)據(jù)進(jìn)行確認(rèn)。

1 驗(yàn)證和確認(rèn)的方法

1.1 驗(yàn)證方法

驗(yàn)證過程主要考慮不完全迭代收斂誤差和空間離散誤差。

對(duì)于網(wǎng)格收斂性的分析方法很多［3］，其中基于誤差行列展開式的Richardson插值法最為常用。利用兩套網(wǎng)格的一般Richardson插值法（Generalized Richardson Extrapolation）［4］（GRE）求解過程為（帶～為近似值）：

其中，εi，j＝fi－fj，p指收斂精度，fk指k網(wǎng)格上的值，“精確解”exact指網(wǎng)格解的趨近值，gi指不隨網(wǎng)格變化的第i階誤差項(xiàng)系數(shù)，hk指k網(wǎng)格網(wǎng)格尺度，為不失其一般性，可將細(xì)網(wǎng)格尺度定為h1＝1，網(wǎng)格細(xì)化比rk，k＋1＝hk＋1／hk。為方便向非結(jié)構(gòu)網(wǎng)格推廣應(yīng)用，有將網(wǎng)格單元總數(shù)作為網(wǎng)格尺度衡量標(biāo)準(zhǔn)的［5－6］，即h＝，NPTS是網(wǎng)格點(diǎn)總數(shù)。

而對(duì)于連續(xù)變化的三套網(wǎng)格，GRE求解過程與兩套網(wǎng)格相似，只是式（1）中k＝1，2，3，且收斂精度p需經(jīng)迭代計(jì)算，當(dāng)r12＝r23＝r時(shí)，De Vahl Davis［7］給出求解表達(dá)式：

另外一種常用的方法即混合一階＋二階精度插值 法 （Mixed 1st＋ 2ndOrder Extrapolation）［4］（MOE），其表達(dá)式：

當(dāng)r＝2時(shí)，有：

觀察收斂精度有：

最后，定義收斂率［8］Rk＝εk21／εk32，只有當(dāng)0＜Rk＜1時(shí)，解才是單調(diào)收斂的，Richardson插值法才可用于估計(jì)誤差和不確定度。

Christopher［4］指出，在使用混合階精度數(shù)值格式時(shí)，如二階精度格式在激波附近降為一階，則當(dāng)一階誤差項(xiàng)與二階誤差項(xiàng)出現(xiàn)反號(hào)而相互抵消時(shí)，做網(wǎng)格精細(xì)化研究時(shí)就可能會(huì)出現(xiàn)解的非單調(diào)收斂情況。

1.2 確認(rèn)方法

1.2.1 實(shí)驗(yàn)數(shù)據(jù)處理方法

一般地，實(shí)驗(yàn)數(shù)據(jù)包含偏差和隨機(jī)誤差，本文中假定實(shí)驗(yàn)測(cè)量偏差為零，且認(rèn)為隨機(jī)誤差為標(biāo)準(zhǔn)分布［9］，即高斯分布，則真值Y（xi）可用N次測(cè)量平均值（xi）代替：

定義標(biāo)準(zhǔn)差：

σ值越小，則（xi）與Y（xi）偏差越小，精度越高，即平均值可信賴程度越高。Stern［10］定義不確定度是對(duì)誤差的估計(jì)，表示包含95%的誤差。高斯分布當(dāng)誤差｜α｜≤2σ時(shí)的概率正好為95%，則文中認(rèn)為實(shí)驗(yàn)不確定度為UD＝2σ。

1.2.2 計(jì)算數(shù)據(jù)的確認(rèn)方法

Wilson［8］給出一種確認(rèn)的方法，首先定義比較誤差（The Comparison Error）E：

其中，D為實(shí)驗(yàn)真值估計(jì)值，通常即平均值，S為計(jì)算值，δD為實(shí)驗(yàn)誤差，δSM為數(shù)值模型誤差，δSN為數(shù)值計(jì)算誤差，所有誤差均是與真值相比的。而數(shù)值計(jì)算誤差可表示為：

最后是確認(rèn)過程，如果｜E｜＜UV，則可認(rèn)為所有的D與S的誤差都在UV范圍內(nèi)，確認(rèn)可以達(dá)到UV水平；如果｜E｜?UV，則E≈δSM，就應(yīng)該對(duì)數(shù)值計(jì)算模型進(jìn)行改進(jìn)。

對(duì)于數(shù)值計(jì)算不確定度USN，Wilson建議采用修正因子［8］，而非Roache［11］采取的安全因子，其表達(dá)式為：

為了考慮全局性質(zhì)，Oberkampf［9］提出了另一種確認(rèn)不確定度尺度：

其中，I指需確認(rèn)點(diǎn)的總數(shù)；y（xi）、Y（xi）分別代表xi位置處的計(jì)算值和實(shí)驗(yàn)值。V越接近1，則實(shí)驗(yàn)值與計(jì)算值的一致性越好。

2 算例模型與計(jì)算方法

2.1 算例模型

在本研究中，選擇了具有較多實(shí)驗(yàn)、計(jì)算數(shù)據(jù)及詳細(xì)來流條件的空心圓柱裙問題（記為T2－97）為研究對(duì)象。該模型具有詳細(xì)的幾何尺寸及流動(dòng)參數(shù)條件（見表1及圖1），所反映的流動(dòng)特征主要是激波－激波干擾，激波－邊界層干擾及激波／分離區(qū)相互干擾。研究這一類問題，對(duì)準(zhǔn)確評(píng)估空間飛行器表面熱流／摩阻、控制面在大偏折角下的效率和熱環(huán)境有重要意義。

2.2 計(jì)算方法

計(jì)算分析由成熟的高超聲速計(jì)算軟件平臺(tái)完成。其中層流粘性系數(shù)由Sutherland公式確定，理想氣體，γ＝1.4。數(shù)值方法采用 MUSCL格式，無(wú)粘通量采取AUSMPW通量分裂格式，使用Minmod限制器，粘性項(xiàng)采取中心差分格式，時(shí)間迭代方式采取標(biāo)準(zhǔn)LU－SGS方法，CFL＝10。具體描述見文獻(xiàn)［11］。

表1 實(shí)驗(yàn)條件Table 1 The conditions of experiments

圖1 實(shí)驗(yàn)?zāi)Ｐ蚑2－97尺寸Fig.1 The size of the T2－97model for experiments

3 網(wǎng)格生成

利用基礎(chǔ)結(jié)構(gòu)網(wǎng)格H1生成其他網(wǎng)格：H2～H16網(wǎng)格從H1抽取，下一層網(wǎng)格點(diǎn)數(shù)量N2與上一層N1關(guān)系如式（16），因軸對(duì)稱流動(dòng)，周向網(wǎng)格四分之一圓周內(nèi)取19個(gè)網(wǎng)格點(diǎn)不變。表2給出了網(wǎng)格規(guī)模及壁面網(wǎng)格雷諾數(shù)，圖2～圖4給出了典型網(wǎng)格分布。

表2 網(wǎng)格分布Table 2 Grids

4 計(jì)算結(jié)果

4.1 不完全迭代收斂誤差

圖5給出了不同網(wǎng)格下計(jì)算的殘差收斂情況，表3是五套網(wǎng)格的計(jì)算結(jié)果在接近迭代收斂過程中結(jié)果的振蕩變化幅度。Stmax－min取壁面幅值變化最大位置的振幅，所有網(wǎng)格Stmax－min位置都在x／L＝0附近，而頭部之后的量級(jí)均在E－7左右，這與St在壁面沿流向分布的最小值E－3相比是在0.01%量級(jí)，則不完全迭代收斂誤差可以忽略。

表3 各網(wǎng)格Stanton數(shù)迭代收斂情況Table 3 Iteration condition of different grids'Stanton numbers

圖5 不同網(wǎng)格的殘差收斂情況Fig.5 Convergence history of different grids

4.2 壁面壓力、熱流及摩阻分布

實(shí)驗(yàn)數(shù)據(jù)由ONERA Chalais－Meudon中心的R5CH風(fēng)洞的多次運(yùn)行結(jié)果提供。分離點(diǎn)與再附著點(diǎn)位置及分離區(qū)大小數(shù)據(jù)見表4。不難發(fā)現(xiàn)，密網(wǎng)格的計(jì)算值與實(shí)驗(yàn)值符合令人滿意。圖6給出了壓力、Stanton數(shù)及摩擦阻力沿壁面的分布，計(jì)算結(jié)果也表明與實(shí)驗(yàn)符合較好。圖7為流場(chǎng)結(jié)構(gòu)圖，可以看出，由于激波與邊界層之間的干擾，導(dǎo)致流動(dòng)分離，從流線圖也可以看到明顯的分離區(qū)。

表4 分離點(diǎn)、再附著點(diǎn)位置及分離區(qū)大小Table 4 The position of separation and reattachment points and the size of the separation region

圖6 （左）壓力系數(shù)分布；（中）Stanton數(shù)分布；（右）摩擦阻力系數(shù)分布Fig.6 The distribution of pressure coefficient（left）；Stanton number（middle）；and the friction coefficient（right）

圖7 （左）馬赫數(shù)等值線；（中）壓力等值線；（右）流線Fig.7 Mach contours（left）；Pressure contours（middle）；Streamlines pattern（right）

5 驗(yàn)證過程

5.1 “精確值”的獲得

為了分析，分別選取有熱流Stanton數(shù)，壓力系數(shù)，分離點(diǎn)與重附著點(diǎn)位置等特征點(diǎn)，分別標(biāo)記如表5。圖8為各網(wǎng)格下數(shù)值解和"精確解"（由GRE或MOE插值得到）分布。網(wǎng)格H1、H2、H4的解有明顯收斂趨勢(shì)，H8、H16解可認(rèn)為已出漸進(jìn)收斂域，Richardson插值已不適用，下文計(jì)算便將其舍棄。兩種插值方式從V2、V4可以看出似乎是MOE要比GRE更接近于網(wǎng)格收斂的解，Christopher［4］也更推薦前者。

表5 特征點(diǎn)選取Table 5 The characteristic points

圖8 各特征點(diǎn)在不同網(wǎng)格下的解與精確解Fig.8 The solutions of different grids at different characteristic points and the exact solutions

5.2 離散誤差

定義不帶安全因子的離散誤差估計(jì)［14－15］：

exact可用GRE或MOE獲得。另外，Roache提出網(wǎng)格收斂指數(shù)（GCI）［11］：

其中安全因子Fs一般可保守取3，在三套或更多網(wǎng)格時(shí)取1.25便已足夠［15］。

此外，還有一種對(duì)于網(wǎng)格收斂特性的評(píng)價(jià)方法，即 Measure－of－Merit（MOM）［9］：

其中，x＝NPTS－，y表示函數(shù)值，其下標(biāo)F、M、C分別表示細(xì)網(wǎng)格、中網(wǎng)格以及粗網(wǎng)格。MOM值越小，則網(wǎng)格解的收斂性越好，如圖9。圖10為各特征點(diǎn)的誤差估計(jì)，其中空間離散誤差的exact由GRE得到，一階、二階以及一階＋二階誤差項(xiàng)，由 MOE得到：

從圖可看到，雖然數(shù)據(jù)點(diǎn)只有三個(gè)，但是基本上還是可以看出在各個(gè)網(wǎng)格上何種誤差起主要作用。如發(fā)生最大壓力系數(shù)位置的捕捉為二階精度，而其值大小的捕捉卻最多只有一階精度。模擬再附點(diǎn)的精度（二階）比模擬分離點(diǎn)的精度（一階）要高。這些結(jié)果與圖8、圖9所反映的是一致的。

圖9 收斂性評(píng)價(jià)Fig.9 Measure of merit

5.3 不同形式離散誤差比較

采用三套網(wǎng)格的GRE、MOE和采用兩套網(wǎng)格的一階GRE和二階GRE的離散誤差表達(dá)式仍為（17）式。另外，分別定義不帶安全因子的一階GCI和二階GCI，在細(xì)網(wǎng)格上表達(dá)式：

圖10 各特征點(diǎn)離散誤差估計(jì)Fig.10 The estimated spatial error of different characteristic points

在粗網(wǎng)格上表達(dá)式：

不同的空間離散誤差的求取方式或者說是表達(dá)形式的結(jié)果在量值有差別，尤其在密網(wǎng)格上，所以在討論離散誤差的時(shí)候需要注意的是離散誤差的定義方式。圖11為誤差變化。雖然在量值上有區(qū)別，但變化趨勢(shì)基本一致，當(dāng)網(wǎng)格變密時(shí)，數(shù)值解趨于收斂值，表明數(shù)值解得到確認(rèn)。

6 確認(rèn)過程

選取特征點(diǎn)計(jì)算，分別有不同位置的熱流Stanton數(shù)，壓力系數(shù)以及分離點(diǎn)與再附點(diǎn)位置，分別標(biāo)記如表6?；隍?yàn)證過程的討論，表7中計(jì)算結(jié)果處理均采用MOE方法。

按照 Wilson［8］的確認(rèn)方式，表中數(shù)據(jù) C1、C2、C3、C4、C5、C7有｜E｜＜UV，可認(rèn)為達(dá)到了UV水平的確認(rèn)，而C6沒有，可認(rèn)為模型誤差較大，需改進(jìn)數(shù)值計(jì)算模型，比如計(jì)算格式，邊界條件等等，當(dāng)然前提是風(fēng)洞測(cè)量無(wú)偏差。

圖11 各特征點(diǎn)離散誤差估計(jì)比較（T2－97－Axisymmetric）Fig.11 The comparison of estimated spatial error at different characteristic points（T2－97－Axisymmetric）

表6 特征點(diǎn)選取Table 6 The characteristic points

表7 不確定度估計(jì)Table 7 The estimatied uncertainties

同時(shí)，我們采用Oberkampf提出的有關(guān)不確定度的確認(rèn)尺度，選取I＝22個(gè)有實(shí)驗(yàn)數(shù)據(jù)的點(diǎn)，進(jìn)行Stanton數(shù)的實(shí)驗(yàn)數(shù)據(jù)與計(jì)算數(shù)據(jù)的確認(rèn)。由式（15）得到不確定度的確認(rèn)尺度V＝0.8634，V越接近于1，說明計(jì)算數(shù)據(jù)與實(shí)驗(yàn)數(shù)據(jù)的一致性越好。

7 結(jié) 論

采用多套連續(xù)變化的網(wǎng)格對(duì)強(qiáng)粘性干擾占優(yōu)的T2－97模型進(jìn)行了研究，通過對(duì)“精確值”、誤差、不確定度的分析及與實(shí)驗(yàn)的對(duì)比，完成了CFD的驗(yàn)證與確認(rèn)。得到以下結(jié)論：

（1）從驗(yàn)證過程可認(rèn)為H1網(wǎng)格在大多特征點(diǎn)上的解已很好地收斂。由于激波，激波／邊界層干擾的存在，計(jì)算格式的精度往往達(dá)不到二階，所以誤差項(xiàng)表現(xiàn)出來既有二階項(xiàng)占主導(dǎo)，也有一階項(xiàng)占主導(dǎo)。離散誤差形式很多，但相互間大多差別不大。

（2）從確認(rèn)過程看，除個(gè)別特征點(diǎn)外，模型的計(jì)算數(shù)據(jù)基本可以認(rèn)為是得到了確認(rèn)。

CFD的驗(yàn)證和確認(rèn)工作是一項(xiàng)長(zhǎng)期而艱巨的研究工作，影響其結(jié)論的因素很多，本文只在這方面進(jìn)行了嘗試，下一步擬在模型、不確定度準(zhǔn)則及驗(yàn)證與確認(rèn)的規(guī)范等方面開展研究。

［1］BOEHM B W.Software engineering economics［R］.Prentice－Hall，Englewood Cliffs，NJ，1981.

［2］BLOTTNER F G.Accurate Navier－Stokes results for the hypersonic flow over a spherical nosetip［J］.JournalofSpacecraftandRockets，1990，27（2）：113－122.

［3］WILSON R，STERN F.Verification and validation for RANS simulation of a naval surface combatant［R］.AIAA Paper 2002－0904，2002.

［4］CHRISTOPHER J ROY.Grid convergence error analysis for mixed－order Numerical schemes［R］.AIAA Paper 2001－2606，2001.

［5］MORRISON J H，HEMSCH M J.Statistical analysis of CFD solutions from the third AIAA Drag Prediction Workshop（invited）［R］.AIAA 2007－254，2007.

［6］VASSBERG J C，TINOCO E N，MANI M，et al.Summary of the Third AIAA CFD Drag Prediction Workshop［R］.AIAA 2007－260，2007.

［7］De VAHL DAVIS G.Natural covection of air in a square cavity：a bench mark numerical solution［J］.InternationalJournalforNumericalMethodsinFluids，1983，3（3）：249－264.

［8］WILSON R，STERN F.Verification and validation for RANS simulation of a naval surface combatant［R］.AIAA 2002－0904，2002.

［9］OBERKAMPF W L，TRUCANO T G.Validation methodology in computational fluid dynamics（invited）［R］.AIAA 2002－2549，2002.

［10］STERN F，WILSON R V，COLEMAN H，et al.Comprehensive approach to verification and validation of CFD simulation－part 1：methodology and procedures［J］.ASMEJ.FluidsEng.，2001.

［11］劉化勇.超聲速引射器的數(shù)值模擬方法及其引射特性研究［D］.［博士學(xué)位論文］.中國(guó)空氣動(dòng)力研究與發(fā)展中心，2009.

［12］ROACHE P J.Perspective：a method for uniform reporting of grid refinement studies［J］.ASMEJournal ofFluidsEngineering，1994，116：405－413.

［13］CANDLER G V，NOMPELIS I，DRUGUET M C，etal.CFD validation for hypersonic flight：hypersonic double－cone flow simulations［R］.AIAA 2002－0581，2002.

［14］MATTHEW MACLEAN，MICHAEL HOLDEN.Validation and comparison of WIND and DPLP results for hypersonic，laminar problems［R］.AIAA Paper 2004－0529，2004.

［15］JOHNSON R W，HUGHES E D.Quantification of uncertainty in computational fluids dynamics－1995［A］.Joint JSME－ASME Fluid Mechanics Meeting ［C］.ASME FED，Vol.213，American Society of Mechanical Engineers，1995.