關于CFD驗證確認中的不確定度和真值估算

張涵信,查 俊

(中國空氣動力研究與發(fā)展中心,四川綿陽,621000；國家計算流體力學實驗室,北京,100191)

0 引 言

隨著計算機技術、計算格式及網格技術等的發(fā)展,計算流體力學(CFD)取得了長足的進步,在基礎研究及工程的應用方面日趨廣泛。然而CFD方法的的可信度(不確定度)或可靠性一直是關心的問題。

AIAA在1998年發(fā)布的《Guide for the Verification and Validation of CFD Simulations》中對誤差(error)和不確定度(uncertainty)給出了如下見解:

誤差是建模和模擬過程中可認知的缺陷,不是由于知識缺乏導致的。(A recognizable deficiency in any phase or activity of modeling and simulation that is not due to lack of knowledge.)而不確定度是由于知識的缺乏,在建模和模擬過程中潛在的缺陷。(A potential deficiency in any phase or activity of the modeling process that is due to lack of knowledge.)

而Roache[1]把誤差定義為計算值或試驗值與真實值的差別。當真實值不確定或者不可知時,計算值或試驗值的誤差就不能確定,這時不確定度就是誤差的估計。

盡管國外的工作已有很多關于這方面的研究[1-15],但是對于CFD不確定度方面的一些基本概念還是缺乏明確和好用的定義。

CFD計算結果誤差的來源一般分為四類[2]:物理模型誤差、離散誤差、計算機舍入誤差、程序設計誤差。

(1)物理模型誤差:其主要源于不精確的物理模型,也就是說,控制方程和邊界條件不能充分地描述要模型化的物理現(xiàn)象。例如湍流模型、流動由層流到湍流的轉捩模式的誤差、氣體狀態(tài)方程與真實情況之間的誤差以及邊界條件表述的誤差等。

(2)離散誤差:其主要來源與各種數(shù)值方法對控制方程及邊界條件的離散化,因空間離散和時間離散的有限精度以及有限分辨率導致數(shù)值解與所求解方程的精確解之間存在誤差；空間網格及表面網格不夠密和不夠光滑所帶來的誤差[3,16]。

(3)舍入誤差:源于計算機數(shù)據(jù)存儲字長的限制。

(4)程序設計誤差:這是簡單的失誤,可以歸為上述未提及的誤差。并且一般可在使用某些方法或者在程序驗證過程中發(fā)現(xiàn)這些錯誤。

此外,還有迭代如何判定收斂的誤差。

CFD計算模擬的不確定度可以分為以下兩類[17]:

(1)模型形式不確定度:它是指數(shù)學模型描述實際物理系統(tǒng)的真實行為時的不確定度,也稱為結構不確定度或者非參量化不確定度。這種類型的不確定度很難用概率密度函數(shù)來描述。在CFD中,湍流模擬即屬此類。

(2)參量不確定度:它是CFD中某些參量(包括網格、算法中的系數(shù))其精確的結果無法得到而產生。

在文獻中,不確定度分析方法有以下幾種[17]:I.非概率方法(Non-probabilistic Methods)

這里又可分為:

(1)區(qū)間分析(Interval Analysis)方法:給出計算值的上限和下限,從而確定不確定度。但要求可能的計算值不能遺漏。

(2)誤差傳播的敏感性導數(shù)(Sensitivity Derivatives)方法:其基本思想是量化輸出結果對輸入參數(shù)微小變化的靈敏度,從而分辨出各種輸入量對于輸出結果不確定度的影響。

(3)模糊邏輯(Fuzzy Logic)方法:模糊邏輯方法是對不精確、不完善的計算結果中輸入參量的不確定范圍,利用模糊邏輯和模糊規(guī)則進行推理分析,從而確定結果的不確定性。在CFD領域現(xiàn)還很難使用。II.概率方法(Probabilistic Methods)

(1)常用的是Monte-Carlo方法[4,18]。該法首先假設概率密度函數(shù)計算輸入參量的誤差,形成樣本,然后求計算樣本值的確定性結果,再對確定輸出結果進行統(tǒng)計分布(如均值、方差)給出不確定度。

Monte-Carlo方法的計算成本相當高昂,因此發(fā)展了許多修正方法。但由于輸入量的概率密度函數(shù)的確定是否合宜,Monte-Carlo方法的結果仍有問題。

(2)矩量法(Moment Method)[5,19,20]

該法給定可用的一階(均值)、二階(方差)、三階(偏度)和四階矩(峰度)來表示概率分析的特征,然后來確定概率分布,從而確定不確定度。

III.隨機微分方程方法(Stochastic Differential Equantion Method)[6]

它是在CFD確定性方程中加入輸入量的隨機變量來計算CFD模擬過程的不確定度。隨機微分方程方法已應用于結構力學問題中,但最近才開始應用于CFD不確定度研究中。

以上情況表明,誤差的定義、來源是明確的。由于真值不能確定,不確定度就定義為對誤差的估計,但什么叫誤差的估計也是不明確的。關于不確定度的理論、方法,雖已提出不少,但能夠具體應用的不多。這就給飛行器CFD的實驗和確認造成了困難。

飛行器的CFD驗證、確認,一般是這樣進行的:選定飛行器外形(根據(jù)需要確定一個或多個),指定來流條件和要驗證的量,分別委托各實驗部門和CFD部門,各進行實驗和計算。一般各實驗部門的實驗數(shù)據(jù)各不相同,各計算部門提供的計算數(shù)據(jù)也不一致。實驗部門必須給出實驗結果的真值和不確定度(或可信度)。計算部門也必須給出計算結果的真值和不確定度,然后兩者進行比較,完成驗證與確認。但實際上,由于真值和不確定度現(xiàn)在還沒有理論確定,因此,驗證、確認的結果,只是給軟件提供者給出他的結果與實驗及其他計算結果的對比情況。對其計算軟件和真值差多遠,得不到確定的結論。這是現(xiàn)有驗證、確認存在的不足。

本文的目的是在不確定度、誤差等術語[21]和試驗中使用的定義[22,23]相同的情況下,給不確定度一個新的解釋方法,從而可給不確定度一個新的表達式,以此為基礎,可給出真值的估計方法以及真值的近似表達式。這樣就容易在CFD驗證和確認中作出應用和判斷。

1 CFD不確定度和真值的估算方法

1.1 CFD不確定度表達式

同實驗所用的概念一樣,我們把數(shù)值解 xC與真值C之差的絕對值的最大值,即:

稱為不確定度。把

稱為相對不確定度。

CFD的不確定度我們可以提出另一種說法來表達。大家知道,工程制造上是按設計數(shù)據(jù)加工的,常有一種說法,加工能達到形狀,可準確到設計數(shù)據(jù)前幾位,這也是表示準確度的一種說法,我們不妨也采用這種方法來表述計算的準確度。即采用計算值可達到真值的前n位來表示計算結果的不確定度。設氣動系數(shù)C可表達為:

如果要求前n位真值準確,它可寫成:

顯然Δ應滿足** ΔC也可表述為×10m-n+1,此時所作的結論小1倍。:

式中,1drag count=10-4。

即絕對不確定度為:

相對不確定度為:

或近似可寫成

可見,如果n=2,即前兩位真值準確,相對不確定度為[10/(am+am-1 10-1)]%；如果n=3,相對不確定度為[10/(am+am-1 10-1)]‰；如果n=4,為萬分之10/(am+am-110-1)。這種表示方法對實驗測量值和計算值均可運用。

對運輸機,如取n=3,即前三位真值準確,CL,Cm,CD的不確定度是:

若給出 δCL,δCm,δCD,它們分別是 :

這里am表示CL,Cm及CD的第1位出現(xiàn)的值。大家知道,式(8)-(9)正是現(xiàn)在實驗能達到的絕對和相對不確定度。

1.2 真值的估算方法

設Ct為計算值,C為真值,則Ct-ΔC≤C≤Ct+ΔC。這表明,計算結果滿足前n位真值準確的數(shù)據(jù)帶的帶寬應該為:

如果一個量有多種計算方法給出計算結果,將它們的數(shù)據(jù)畫出,就可給出一個數(shù)據(jù)帶。若這個數(shù)據(jù)帶滿足式(10),這里面的數(shù)據(jù)就滿足前n位真值準確。這樣我們就確定了真值的前n位。

做為例子,對運輸機,表1給出了運輸機的數(shù)值帶寬2ΔCL,2ΔCm及2ΔCD與n的關系。當數(shù)據(jù)帶落在n的范圍內時,真值前n位就可確定。

用這種方法,我們可以進一步估算真值。事實上要更準確的估算真值(計算值或者實驗值),需作大量的計算或實驗。即需要大的數(shù)據(jù)樣本,此時可引用大數(shù)定律和統(tǒng)計理論。

表1 n與2ΔC的關系Table 1 Relations between n and 2ΔC

設多個計算或多個實驗值給出的離散數(shù)據(jù)ξ為xi,i=1,2,…,且 P(ξ=xi)=Pi為其出現(xiàn)的概率,則加權平均值x

是ξ的數(shù)學期望,對任意的ε＞0,其出現(xiàn)的概率滿足:

若稱Δi=|xi-x|為方差,其

這表明,當數(shù)據(jù)樣本m足夠大時,加權平均值

是樣本中最可能出現(xiàn)的。當已知數(shù)值滿足前n位真值準確后,以此可給出真值的估算方法,建議分兩步進行:

(1)預測:

先預計一個權分布。例如在CFD的求解中,網格數(shù)越大,模型越好,計算方法精度越高,邊界處理越好,一般求解越準確,Ni就越大?？梢园?Ni先作為權,于是加權平均值可表達為:

這里xi為計算值或實驗值。

令|xi-|=Δi,因 Δi越大,偏差越大。因此可用qi=作為進一步的權值,此時可得:

式(13)即可作為預測的真值,從而確定數(shù)據(jù)(x1,x2…xm)接近真值幾位,可以證明,它的誤差是:

(2)修正

在決定各計算值或實驗值接近真值的預測的位數(shù)后,例如n=2,在這種情況下,我們把預測中已經準確到2位的數(shù)據(jù)集中起來,略去不到2位真值準確的數(shù)據(jù)。然后利用預測步的第二步公式進行重新計算,這樣得到的結果可能是2位真值準確的最優(yōu)結果。我們稱之為最優(yōu)解或最優(yōu)值。用這個解,再回頭評價各個CFD軟件或實驗結果,從而分別給出對它們精度的評價。

2 應用舉例

2.1 模型方程確定,用CFD求解NS方程給出駐點壓力、熱流和摩阻的驗證

為簡單,這里僅討論高超聲速圓柱繞流。來流條件為:M∞=8.03,T∞=124.94K,Tw=294.44K,Re=1.835×105,壁面采用等溫壁條件。這個例子有實驗結果[24]。

求解分別用了三種方法:NND格式、三階緊致(CC3)、五階緊致(CC5),每個方法中分別用4套網格,由于每套網格中壁面附近最小網格Δhmin又有兩種不同,所以可以說每種算法中有八套網格。表2～表4分別給出計算得到的駐點壓力、摩阻系數(shù)及駐點點熱流的結果。

表2 NND的結果Table 2 Computed results using NND schemes

表 3 CC3的結果Table 3 Computed results using CC3 method

表 4 CC5的結果Table 4 Computed results using CC5 method

表5是利用上節(jié)理論給出的最優(yōu)解及其相應的誤差,還列出了實驗結果?？梢?計算給出的真值,與經認真檢驗的實驗值相當接近。這說明,本文的真值估算方法是正確有效的。

表6是根據(jù)最優(yōu)解給出各個計算結果的比較。

2.2 DPW第二次驗證會議阻力數(shù)據(jù)的分析驗證

這里僅引用DLR-F6無發(fā)動機艙的各家阻力計算數(shù)據(jù)[25,26]。

表7是網格數(shù)與相應阻力系數(shù)的數(shù)據(jù)表(它是根據(jù)圖讀出來的),計算條件是:M∞=0.75,Re=3×106,CL=0.5。

表8給出了根據(jù)上節(jié)理論給出的最優(yōu)值,它與實驗值相當接近。這再次說明,本文真值估算方法是滿意的。

圖1還畫出了兩位真值準確數(shù)據(jù)帶。

表5 最優(yōu)解的結果Table 5 Optimum computed results

表6 各軟件結果的比較Table6 Comparison of computational results given by different codes

表7 網格與阻力系數(shù)的數(shù)據(jù)表Table7 Drag coeff icients and grids

5.252E+06 0.029298 -0.36 2 6.270E+06 0.029865 1.53 2 4.114E+06 0.027384 -7.38 1 2.836E+06 0.025825 -13.8 1 6.271E+06 0.029865 1.53 2 6.510E+06 0.029014 -1.35 2 8.188E+06 0.029440 0.11 3 8.667E+06 0.028660 -2.60 2 8.827E+06 0.028589 -2.85 2 9.086E+06 0.028235 -4.14 1 1.010E+07 0.029794 1.30 2 9.945E+06 0.029298 -0.36 2 9.546E+06 0.029156 -0.85 2 9.945E+06 0.028660 -2.60 2 1.130E+07 0.027810 -5.73 1 1.322E+07 0.028306 -3.88 1 1.258E+07 0.029227 -0.61 2 2.312E+07 0.029014 -1.35 2 2.256E+07 0.028306 -3.88 1 2.256E+07 0.027951 -5.20 1

表8 最優(yōu)解的結果Table 8 Optimum computed results

圖1 二位真值準確的數(shù)據(jù)帶Fig.1 Thezoneof approximating to first twodigil number of truth value

3 結 論

本文回顧了CFD驗證、確認中不確定度的概念和研究方法,CFD的不確定度尚無表達式可以使用。本文也討論了現(xiàn)在正在進行的實驗驗證,對各個參加驗證的軟件,如何作出定量的精度評價也缺乏原則。針對這些情況,我們在不改變不確定度定義的前提下,對不確定度作了新的解讀,即不確定度可解讀為計算值或實驗值與真值準確到前n位,從而可給出不確定度的表達式和真值估算的原則。并根據(jù)大樣本數(shù)據(jù)的統(tǒng)計理論,對真值認為接近數(shù)學期望,從而給出準確到n位真值的計算方法。這個方法,可用于計算結果的檢驗,例如當模型一定時,可用此法尋求計算方法的真值,對算法進行檢驗；如算法一定模型改變時,也可用于檢驗模型的可靠性。利用這種方法,在沒有實驗結果的情況下,也可評價各計算軟件的質量。文中方法當然也可以運用處理實驗數(shù)據(jù)。因為CFD中計算模型是人為建立的,雖然可以檢驗它的解是否正確,但與物理情況是否一致,并未得到回答。因此,開展實驗驗證及確認是必需的。

[1]ROACHE P J.Verification of codes and calculations[J].AIAA Journal,1998,36(5):696-702.

[2]OBERKAMPF W L,BLOTTNER F G.Issues in computational fluid dynamics code verification and validation[J].AIAA Journal,1998,36:687-695.

[3]ROACHE P J.Quantification of uncertainty in computational fluid dynamics[J].Annual Review of Fluid Mechanics,1997,29:123-160.

[4]WALTERSR W,HUYSE L.Uncertainty analysis for fluid mechanics with applications[R].NASA/CR-2002-211449,ICASE Report No.2002-1,2002.

[5]PUTKO M M,NEWMAN P A,TAYLOR A C,GREEN L L.Approach for uncertainty propagation and robust design in CFD using sensitivity derivatives[R].AIAA Paper,2001:2001-2558.

[6]MATHELIN L,HUSSAINI M Y,et al.Uncertainty propagation for turbulent,compressiblef low in a quasi-1D nozzle using stochastic methods[A].16thAIAA Computational Fluid Dynamics Conference[C].Orlando,Florida,AIAA,2003:2003-4240.

[7]LUCOR D,XIU D,et al.Predictability and uncertainty in CFD[J].Int.J.Numer.Meth.Fluids,2003,43(5):483-505.

[8]OBERKAM PF W L,TRUCANO T G.Verification and validation in computational fluid dynamics[J].Prog.Aero.Sci.,2002,38:209-272.

[9]LUCKRING J M,HEMSCH M J,MORRISON J H.Uncertainty in computational aerodynamics[R].AIAA-2003-0409,2003.

[10]RAYMOND R,et al.Theimportanceof uncertainty estimation in computational f luid dynamics[R].AIAA-2003-0406,2003.

[11]FREITAS C J,GHIA U,CELIK I,ROACHE P,RAAD P.AMSE'S quest to quantify numerical uncertainty[R].AIAA-2003-627,2003.

[12]ROACHE J.Need for control of numerical accuracy[J].J.Spacecraf t and Rockets,1990,27(2):98-102.

[13]COLEMAN H W,STERN F.Uncertainties and CFD code validation[J].Journal of Fluids Engineering,1997,119(4):795-803.

[14]Quantifying uncertainty in CFD[J].Journal of Fluids Engineering,2002,124(1):2-3.

[15]B.DE VOLDER,GLIMM J,GROVE JW,KANG Y,LEEY,PAO K,SHARP D H,YE K.Uncertainty quantification for multiscalesimulations[J].Journal of Fluids Engineering,2002,124(1):29-40.

[16]ROACHE P J.Quantification of uncertainty in computational fluid dynamics[J].Annual Review of Fluid Mechanics,1997,29:123-160.

[17]FA RAGHER.Probabilistic methods for the quantification of uncertainty and error in computational fluid dynamics simulations[R].DSTO-TR-1633,2004.

[18]HAMMERSLEY J M,HANDSCOMB D C.Monte Carlo methods,methuen's monographs on applied probability and statistics[M].Flether&Son Ltd.,Norwich,1964.

[19]HUYSE L.Free-form airfoil shape optimization under uncertainty using maximum expected value and secondorder second-moment strategies[R].Tech.Report,ICASE Report 2001-18/NASA CR 2001-211020,2001.

[20]HUYSE L,LEWIS RM.Aerodynamic shapeoptimization of two-dimensional airfoils under certain conditions[R].Tech.Report,ICASE Report 2001-1/NASA CR 2001-210648,2001.

[21]張涵信.關于CFD計算結果的不確定度問題[J].空氣動力學學報,2008,26(1):47-49.

[22]惲起麟.風洞實驗[M].近代空氣動力學叢書,國防工業(yè)出版社,北京,2000.

[23]程厚梅等.風洞實驗干擾與修正[M].近代空氣動力學叢書,國防工業(yè)出版社,北京,2003.

[24]WIETING A R.Experimental study of shock wave interference heating on a cylindrical leading edge[R].NASA TM-100484,1987.

[25]LAFLIN R,et al.Summary of data from the second AIAA CFD drag prediction workshop(Invited)[R].AIAA-2004-0555,2004.

[26]HEMSCH M J,M ORRISON J H.Statistical analysis of CFD solutions from 2nd drag prediction workshop[R].AIAA 2004-556,2004.