基于Richardson插值法的CFD驗證和確認方法的研究

陳堅強， 張益榮

（中國空氣動力研究與發展中心，四川 綿陽621000）

0 引 言

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

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

本文通過采取多套不同尺度但拓撲一致的網格對T2－97模型的計算，利用Richardson插值法，得到“精確解”與離散誤差，主要考察不完全迭代收斂誤差和空間離散誤差，通過網格收斂性研究完成了CFD的驗證，并繼而與實驗數據進行確認。

1 驗證和確認的方法

1.1 驗證方法

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

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

其中，εi，j＝fi－fj，p指收斂精度，fk指k網格上的值，“精確解”exact指網格解的趨近值，gi指不隨網格變化的第i階誤差項系數，hk指k網格網格尺度，為不失其一般性，可將細網格尺度定為h1＝1，網格細化比rk，k＋1＝hk＋1／hk。為方便向非結構網格推廣應用，有將網格單元總數作為網格尺度衡量標準的［5－6］，即h＝，NPTS是網格點總數。

而對于連續變化的三套網格，GRE求解過程與兩套網格相似，只是式（1）中k＝1，2，3，且收斂精度p需經迭代計算，當r12＝r23＝r時，De Vahl Davis［7］給出求解表達式：

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

當r＝2時，有：

觀察收斂精度有：

最后，定義收斂率［8］Rk＝εk21／εk32，只有當0＜Rk＜1時，解才是單調收斂的，Richardson插值法才可用于估計誤差和不確定度。

Christopher［4］指出，在使用混合階精度數值格式時，如二階精度格式在激波附近降為一階，則當一階誤差項與二階誤差項出現反號而相互抵消時，做網格精細化研究時就可能會出現解的非單調收斂情況。

1.2 確認方法

1.2.1 實驗數據處理方法

一般地，實驗數據包含偏差和隨機誤差，本文中假定實驗測量偏差為零，且認為隨機誤差為標準分布［9］，即高斯分布，則真值Y（xi）可用N次測量平均值（xi）代替：

定義標準差：

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

1.2.2 計算數據的確認方法

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

其中，D為實驗真值估計值，通常即平均值，S為計算值，δD為實驗誤差，δSM為數值模型誤差，δSN為數值計算誤差，所有誤差均是與真值相比的。而數值計算誤差可表示為：

最后是確認過程，如果｜E｜＜UV，則可認為所有的D與S的誤差都在UV范圍內，確認可以達到UV水平；如果｜E｜?UV，則E≈δSM，就應該對數值計算模型進行改進。

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

為了考慮全局性質，Oberkampf［9］提出了另一種確認不確定度尺度：

其中，I指需確認點的總數；y（xi）、Y（xi）分別代表xi位置處的計算值和實驗值。V越接近1，則實驗值與計算值的一致性越好。

2 算例模型與計算方法

2.1 算例模型

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

2.2 計算方法

計算分析由成熟的高超聲速計算軟件平臺完成。其中層流粘性系數由Sutherland公式確定，理想氣體，γ＝1.4。數值方法采用 MUSCL格式，無粘通量采取AUSMPW通量分裂格式，使用Minmod限制器，粘性項采取中心差分格式，時間迭代方式采取標準LU－SGS方法，CFL＝10。具體描述見文獻［11］。

表1 實驗條件Table 1 The conditions of experiments

圖1 實驗模型T2－97尺寸Fig.1 The size of the T2－97model for experiments

3 網格生成

利用基礎結構網格H1生成其他網格：H2～H16網格從H1抽取，下一層網格點數量N2與上一層N1關系如式（16），因軸對稱流動，周向網格四分之一圓周內取19個網格點不變。表2給出了網格規模及壁面網格雷諾數，圖2～圖4給出了典型網格分布。

表2 網格分布Table 2 Grids

4 計算結果

4.1 不完全迭代收斂誤差

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

表3 各網格Stanton數迭代收斂情況Table 3 Iteration condition of different grids'Stanton numbers

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

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

實驗數據由ONERA Chalais－Meudon中心的R5CH風洞的多次運行結果提供。分離點與再附著點位置及分離區大小數據見表4。不難發現，密網格的計算值與實驗值符合令人滿意。圖6給出了壓力、Stanton數及摩擦阻力沿壁面的分布，計算結果也表明與實驗符合較好。圖7為流場結構圖，可以看出，由于激波與邊界層之間的干擾，導致流動分離，從流線圖也可以看到明顯的分離區。

表4 分離點、再附著點位置及分離區大小Table 4 The position of separation and reattachment points and the size of the separation region

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

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

5 驗證過程

5.1 “精確值”的獲得

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

表5 特征點選取Table 5 The characteristic points

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

5.2 離散誤差

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

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

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

此外，還有一種對于網格收斂特性的評價方法，即 Measure－of－Merit（MOM）［9］：

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

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

圖9 收斂性評價Fig.9 Measure of merit

5.3 不同形式離散誤差比較

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

圖10 各特征點離散誤差估計Fig.10 The estimated spatial error of different characteristic points

在粗網格上表達式：

不同的空間離散誤差的求取方式或者說是表達形式的結果在量值有差別，尤其在密網格上，所以在討論離散誤差的時候需要注意的是離散誤差的定義方式。圖11為誤差變化。雖然在量值上有區別，但變化趨勢基本一致，當網格變密時，數值解趨于收斂值，表明數值解得到確認。

6 確認過程

選取特征點計算，分別有不同位置的熱流Stanton數，壓力系數以及分離點與再附點位置，分別標記如表6。基于驗證過程的討論，表7中計算結果處理均采用MOE方法。

按照 Wilson［8］的確認方式，表中數據 C1、C2、C3、C4、C5、C7有｜E｜＜UV，可認為達到了UV水平的確認，而C6沒有，可認為模型誤差較大，需改進數值計算模型，比如計算格式，邊界條件等等，當然前提是風洞測量無偏差。

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

表6 特征點選取Table 6 The characteristic points

表7 不確定度估計Table 7 The estimatied uncertainties

同時，我們采用Oberkampf提出的有關不確定度的確認尺度，選取I＝22個有實驗數據的點，進行Stanton數的實驗數據與計算數據的確認。由式（15）得到不確定度的確認尺度V＝0.8634，V越接近于1，說明計算數據與實驗數據的一致性越好。

7 結 論

采用多套連續變化的網格對強粘性干擾占優的T2－97模型進行了研究，通過對“精確值”、誤差、不確定度的分析及與實驗的對比，完成了CFD的驗證與確認。得到以下結論：

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

（2）從確認過程看，除個別特征點外，模型的計算數據基本可以認為是得到了確認。

CFD的驗證和確認工作是一項長期而艱巨的研究工作，影響其結論的因素很多，本文只在這方面進行了嘗試，下一步擬在模型、不確定度準則及驗證與確認的規范等方面開展研究。

［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］劉化勇.超聲速引射器的數值模擬方法及其引射特性研究［D］.［博士學位論文］.中國空氣動力研究與發展中心，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.