張涵信,查 俊
(中國空氣動力研究與發(fā)展中心,四川綿陽,621000;國家計算流體力學(xué)實驗室,北京,100191)
隨著計算機技術(shù)、計算格式及網(wǎng)格技術(shù)等的發(fā)展,計算流體力學(xué)(CFD)取得了長足的進步,在基礎(chǔ)研究及工程的應(yīng)用方面日趨廣泛。然而CFD方法的的可信度(不確定度)或可靠性一直是關(guān)心的問題。
AIAA在1998年發(fā)布的《Guide for the Verification and Validation of CFD Simulations》中對誤差(error)和不確定度(uncertainty)給出了如下見解:
誤差是建模和模擬過程中可認知的缺陷,不是由于知識缺乏導(dǎo)致的。(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]把誤差定義為計算值或試驗值與真實值的差別。當真實值不確定或者不可知時,計算值或試驗值的誤差就不能確定,這時不確定度就是誤差的估計。
盡管國外的工作已有很多關(guān)于這方面的研究[1-15],但是對于CFD不確定度方面的一些基本概念還是缺乏明確和好用的定義。
CFD計算結(jié)果誤差的來源一般分為四類[2]:物理模型誤差、離散誤差、計算機舍入誤差、程序設(shè)計誤差。
(1)物理模型誤差:其主要源于不精確的物理模型,也就是說,控制方程和邊界條件不能充分地描述要模型化的物理現(xiàn)象。例如湍流模型、流動由層流到湍流的轉(zhuǎn)捩模式的誤差、氣體狀態(tài)方程與真實情況之間的誤差以及邊界條件表述的誤差等。
(2)離散誤差:其主要來源與各種數(shù)值方法對控制方程及邊界條件的離散化,因空間離散和時間離散的有限精度以及有限分辨率導(dǎo)致數(shù)值解與所求解方程的精確解之間存在誤差;空間網(wǎng)格及表面網(wǎng)格不夠密和不夠光滑所帶來的誤差[3,16]。
(3)舍入誤差:源于計算機數(shù)據(jù)存儲字長的限制。
(4)程序設(shè)計誤差:這是簡單的失誤,可以歸為上述未提及的誤差。并且一般可在使用某些方法或者在程序驗證過程中發(fā)現(xiàn)這些錯誤。
此外,還有迭代如何判定收斂的誤差。
CFD計算模擬的不確定度可以分為以下兩類[17]:
(1)模型形式不確定度:它是指數(shù)學(xué)模型描述實際物理系統(tǒng)的真實行為時的不確定度,也稱為結(jié)構(gòu)不確定度或者非參量化不確定度。這種類型的不確定度很難用概率密度函數(shù)來描述。在CFD中,湍流模擬即屬此類。
(2)參量不確定度:它是CFD中某些參量(包括網(wǎng)格、算法中的系數(shù))其精確的結(jié)果無法得到而產(chǎn)生。
在文獻中,不確定度分析方法有以下幾種[17]:I.非概率方法(Non-probabilistic Methods)
這里又可分為:
(1)區(qū)間分析(Interval Analysis)方法:給出計算值的上限和下限,從而確定不確定度。但要求可能的計算值不能遺漏。
(2)誤差傳播的敏感性導(dǎo)數(shù)(Sensitivity Derivatives)方法:其基本思想是量化輸出結(jié)果對輸入?yún)?shù)微小變化的靈敏度,從而分辨出各種輸入量對于輸出結(jié)果不確定度的影響。
(3)模糊邏輯(Fuzzy Logic)方法:模糊邏輯方法是對不精確、不完善的計算結(jié)果中輸入?yún)⒘康牟淮_定范圍,利用模糊邏輯和模糊規(guī)則進行推理分析,從而確定結(jié)果的不確定性。在CFD領(lǐng)域現(xiàn)還很難使用。II.概率方法(Probabilistic Methods)
(1)常用的是Monte-Carlo方法[4,18]。該法首先假設(shè)概率密度函數(shù)計算輸入?yún)⒘康恼`差,形成樣本,然后求計算樣本值的確定性結(jié)果,再對確定輸出結(jié)果進行統(tǒng)計分布(如均值、方差)給出不確定度。
Monte-Carlo方法的計算成本相當高昂,因此發(fā)展了許多修正方法。但由于輸入量的概率密度函數(shù)的確定是否合宜,Monte-Carlo方法的結(jié)果仍有問題。
(2)矩量法(Moment Method)[5,19,20]
該法給定可用的一階(均值)、二階(方差)、三階(偏度)和四階矩(峰度)來表示概率分析的特征,然后來確定概率分布,從而確定不確定度。
III.隨機微分方程方法(Stochastic Differential Equantion Method)[6]
它是在CFD確定性方程中加入輸入量的隨機變量來計算CFD模擬過程的不確定度。隨機微分方程方法已應(yīng)用于結(jié)構(gòu)力學(xué)問題中,但最近才開始應(yīng)用于CFD不確定度研究中。
以上情況表明,誤差的定義、來源是明確的。由于真值不能確定,不確定度就定義為對誤差的估計,但什么叫誤差的估計也是不明確的。關(guān)于不確定度的理論、方法,雖已提出不少,但能夠具體應(yīng)用的不多。這就給飛行器CFD的實驗和確認造成了困難。
飛行器的CFD驗證、確認,一般是這樣進行的:選定飛行器外形(根據(jù)需要確定一個或多個),指定來流條件和要驗證的量,分別委托各實驗部門和CFD部門,各進行實驗和計算。一般各實驗部門的實驗數(shù)據(jù)各不相同,各計算部門提供的計算數(shù)據(jù)也不一致。實驗部門必須給出實驗結(jié)果的真值和不確定度(或可信度)。計算部門也必須給出計算結(jié)果的真值和不確定度,然后兩者進行比較,完成驗證與確認。但實際上,由于真值和不確定度現(xiàn)在還沒有理論確定,因此,驗證、確認的結(jié)果,只是給軟件提供者給出他的結(jié)果與實驗及其他計算結(jié)果的對比情況。對其計算軟件和真值差多遠,得不到確定的結(jié)論。這是現(xiàn)有驗證、確認存在的不足。
本文的目的是在不確定度、誤差等術(shù)語[21]和試驗中使用的定義[22,23]相同的情況下,給不確定度一個新的解釋方法,從而可給不確定度一個新的表達式,以此為基礎(chǔ),可給出真值的估計方法以及真值的近似表達式。這樣就容易在CFD驗證和確認中作出應(yīng)用和判斷。
同實驗所用的概念一樣,我們把數(shù)值解 xC與真值C之差的絕對值的最大值,即:
稱為不確定度。把
稱為相對不確定度。
CFD的不確定度我們可以提出另一種說法來表達。大家知道,工程制造上是按設(shè)計數(shù)據(jù)加工的,常有一種說法,加工能達到形狀,可準確到設(shè)計數(shù)據(jù)前幾位,這也是表示準確度的一種說法,我們不妨也采用這種方法來表述計算的準確度。即采用計算值可達到真值的前n位來表示計算結(jié)果的不確定度。設(shè)氣動系數(shù)C可表達為:
如果要求前n位真值準確,它可寫成:
顯然Δ應(yīng)滿足** ΔC也可表述為×10m-n+1,此時所作的結(jié)論小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)在實驗?zāi)苓_到的絕對和相對不確定度。
設(shè)Ct為計算值,C為真值,則Ct-ΔC≤C≤Ct+ΔC。這表明,計算結(jié)果滿足前n位真值準確的數(shù)據(jù)帶的帶寬應(yīng)該為:
如果一個量有多種計算方法給出計算結(jié)果,將它們的數(shù)據(jù)畫出,就可給出一個數(shù)據(jù)帶。若這個數(shù)據(jù)帶滿足式(10),這里面的數(shù)據(jù)就滿足前n位真值準確。這樣我們就確定了真值的前n位。
做為例子,對運輸機,表1給出了運輸機的數(shù)值帶寬2ΔCL,2ΔCm及2ΔCD與n的關(guān)系。當數(shù)據(jù)帶落在n的范圍內(nèi)時,真值前n位就可確定。
用這種方法,我們可以進一步估算真值。事實上要更準確的估算真值(計算值或者實驗值),需作大量的計算或?qū)嶒?。即需要大的?shù)據(jù)樣本,此時可引用大數(shù)定律和統(tǒng)計理論。
表1 n與2ΔC的關(guān)系Table 1 Relations between n and 2ΔC
設(shè)多個計算或多個實驗值給出的離散數(shù)據(jù)ξ為xi,i=1,2,…,且 P(ξ=xi)=Pi為其出現(xiàn)的概率,則加權(quán)平均值x
是ξ的數(shù)學(xué)期望,對任意的ε>0,其出現(xiàn)的概率滿足:
若稱Δi=|xi-x|為方差,其
這表明,當數(shù)據(jù)樣本m足夠大時,加權(quán)平均值
是樣本中最可能出現(xiàn)的。當已知數(shù)值滿足前n位真值準確后,以此可給出真值的估算方法,建議分兩步進行:
(1)預(yù)測:
先預(yù)計一個權(quán)分布。例如在CFD的求解中,網(wǎng)格數(shù)越大,模型越好,計算方法精度越高,邊界處理越好,一般求解越準確,Ni就越大。可以把 Ni先作為權(quán),于是加權(quán)平均值可表達為:
這里xi為計算值或?qū)嶒炛怠?/p>
令|xi-|=Δi,因 Δi越大,偏差越大。因此可用qi=作為進一步的權(quán)值,此時可得:
式(13)即可作為預(yù)測的真值,從而確定數(shù)據(jù)(x1,x2…xm)接近真值幾位,可以證明,它的誤差是:
(2)修正
在決定各計算值或?qū)嶒炛到咏嬷档念A(yù)測的位數(shù)后,例如n=2,在這種情況下,我們把預(yù)測中已經(jīng)準確到2位的數(shù)據(jù)集中起來,略去不到2位真值準確的數(shù)據(jù)。然后利用預(yù)測步的第二步公式進行重新計算,這樣得到的結(jié)果可能是2位真值準確的最優(yōu)結(jié)果。我們稱之為最優(yōu)解或最優(yōu)值。用這個解,再回頭評價各個CFD軟件或?qū)嶒灲Y(jié)果,從而分別給出對它們精度的評價。
為簡單,這里僅討論高超聲速圓柱繞流。來流條件為:M∞=8.03,T∞=124.94K,Tw=294.44K,Re=1.835×105,壁面采用等溫壁條件。這個例子有實驗結(jié)果[24]。
求解分別用了三種方法:NND格式、三階緊致(CC3)、五階緊致(CC5),每個方法中分別用4套網(wǎng)格,由于每套網(wǎng)格中壁面附近最小網(wǎng)格Δhmin又有兩種不同,所以可以說每種算法中有八套網(wǎng)格。表2~表4分別給出計算得到的駐點壓力、摩阻系數(shù)及駐點點熱流的結(jié)果。
表2 NND的結(jié)果Table 2 Computed results using NND schemes
表 3 CC3的結(jié)果Table 3 Computed results using CC3 method
表 4 CC5的結(jié)果Table 4 Computed results using CC5 method
表5是利用上節(jié)理論給出的最優(yōu)解及其相應(yīng)的誤差,還列出了實驗結(jié)果??梢?計算給出的真值,與經(jīng)認真檢驗的實驗值相當接近。這說明,本文的真值估算方法是正確有效的。
表6是根據(jù)最優(yōu)解給出各個計算結(jié)果的比較。
這里僅引用DLR-F6無發(fā)動機艙的各家阻力計算數(shù)據(jù)[25,26]。
表7是網(wǎng)格數(shù)與相應(yīng)阻力系數(shù)的數(shù)據(jù)表(它是根據(jù)圖讀出來的),計算條件是:M∞=0.75,Re=3×106,CL=0.5。
表8給出了根據(jù)上節(jié)理論給出的最優(yōu)值,它與實驗值相當接近。這再次說明,本文真值估算方法是滿意的。
圖1還畫出了兩位真值準確數(shù)據(jù)帶。
表5 最優(yōu)解的結(jié)果Table 5 Optimum computed results
表6 各軟件結(jié)果的比較Table6 Comparison of computational results given by different codes
表7 網(wǎng)格與阻力系數(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)解的結(jié)果Table 8 Optimum computed results
圖1 二位真值準確的數(shù)據(jù)帶Fig.1 Thezoneof approximating to first twodigil number of truth value
本文回顧了CFD驗證、確認中不確定度的概念和研究方法,CFD的不確定度尚無表達式可以使用。本文也討論了現(xiàn)在正在進行的實驗驗證,對各個參加驗證的軟件,如何作出定量的精度評價也缺乏原則。針對這些情況,我們在不改變不確定度定義的前提下,對不確定度作了新的解讀,即不確定度可解讀為計算值或?qū)嶒炛蹬c真值準確到前n位,從而可給出不確定度的表達式和真值估算的原則。并根據(jù)大樣本數(shù)據(jù)的統(tǒng)計理論,對真值認為接近數(shù)學(xué)期望,從而給出準確到n位真值的計算方法。這個方法,可用于計算結(jié)果的檢驗,例如當模型一定時,可用此法尋求計算方法的真值,對算法進行檢驗;如算法一定模型改變時,也可用于檢驗?zāi)P偷目煽啃?。利用這種方法,在沒有實驗結(jié)果的情況下,也可評價各計算軟件的質(zhì)量。文中方法當然也可以運用處理實驗數(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]張涵信.關(guān)于CFD計算結(jié)果的不確定度問題[J].空氣動力學(xué)學(xué)報,2008,26(1):47-49.
[22]惲起麟.風(fēng)洞實驗[M].近代空氣動力學(xué)叢書,國防工業(yè)出版社,北京,2000.
[23]程厚梅等.風(fēng)洞實驗干擾與修正[M].近代空氣動力學(xué)叢書,國防工業(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.