張涵信,查 俊

(中國(guó)空氣動(dòng)力研究與發(fā)展中心,四川綿陽,621000；國(guó)家計(jì)算流體力學(xué)實(shí)驗(yàn)室,北京,100191)

0 引 言

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

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

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

而Roache[1]把誤差定義為計(jì)算值或試驗(yàn)值與真實(shí)值的差別。當(dāng)真實(shí)值不確定或者不可知時(shí),計(jì)算值或試驗(yàn)值的誤差就不能確定,這時(shí)不確定度就是誤差的估計(jì)。

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

CFD計(jì)算結(jié)果誤差的來源一般分為四類[2]:物理模型誤差、離散誤差、計(jì)算機(jī)舍入誤差、程序設(shè)計(jì)誤差。

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

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

(3)舍入誤差:源于計(jì)算機(jī)數(shù)據(jù)存儲(chǔ)字長(zhǎng)的限制。

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

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

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

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

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

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

這里又可分為:

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

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

(3)模糊邏輯(Fuzzy Logic)方法:模糊邏輯方法是對(duì)不精確、不完善的計(jì)算結(jié)果中輸入?yún)⒘康牟淮_定范圍,利用模糊邏輯和模糊規(guī)則進(jìn)行推理分析,從而確定結(jié)果的不確定性。在CFD領(lǐng)域現(xiàn)還很難使用。II.概率方法(Probabilistic Methods)

(1)常用的是Monte-Carlo方法[4,18]。該法首先假設(shè)概率密度函數(shù)計(jì)算輸入?yún)⒘康恼`差,形成樣本,然后求計(jì)算樣本值的確定性結(jié)果,再對(duì)確定輸出結(jié)果進(jìn)行統(tǒng)計(jì)分布(如均值、方差)給出不確定度。

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

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

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

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

它是在CFD確定性方程中加入輸入量的隨機(jī)變量來計(jì)算CFD模擬過程的不確定度。隨機(jī)微分方程方法已應(yīng)用于結(jié)構(gòu)力學(xué)問題中,但最近才開始應(yīng)用于CFD不確定度研究中。

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

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

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

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

1.1 CFD不確定度表達(dá)式

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

稱為不確定度。把

稱為相對(duì)不確定度。

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

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

顯然Δ應(yīng)滿足** ΔC也可表述為×10m-n+1,此時(shí)所作的結(jié)論小1倍。:

式中,1drag count=10-4。

即絕對(duì)不確定度為:

相對(duì)不確定度為:

或近似可寫成

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

對(duì)運(yùn)輸機(jī),如取n=3,即前三位真值準(zhǔn)確,CL,Cm,CD的不確定度是:

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

這里am表示CL,Cm及CD的第1位出現(xiàn)的值。大家知道,式(8)-(9)正是現(xiàn)在實(shí)驗(yàn)?zāi)苓_(dá)到的絕對(duì)和相對(duì)不確定度。

1.2 真值的估算方法

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

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

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

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

表1 n與2ΔC的關(guān)系Table 1 Relations between n and 2ΔC

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

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

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

這表明,當(dāng)數(shù)據(jù)樣本m足夠大時(shí),加權(quán)平均值

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

(1)預(yù)測(cè):

先預(yù)計(jì)一個(gè)權(quán)分布。例如在CFD的求解中,網(wǎng)格數(shù)越大,模型越好,計(jì)算方法精度越高,邊界處理越好,一般求解越準(zhǔn)確,Ni就越大。可以把 Ni先作為權(quán),于是加權(quán)平均值可表達(dá)為:

這里xi為計(jì)算值或?qū)嶒?yàn)值。

令|xi-|=Δi,因 Δi越大,偏差越大。因此可用qi=作為進(jìn)一步的權(quán)值,此時(shí)可得:

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

(2)修正

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

2 應(yīng)用舉例

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

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

求解分別用了三種方法:NND格式、三階緊致(CC3)、五階緊致(CC5),每個(gè)方法中分別用4套網(wǎng)格,由于每套網(wǎng)格中壁面附近最小網(wǎng)格Δhmin又有兩種不同,所以可以說每種算法中有八套網(wǎng)格。表2～表4分別給出計(jì)算得到的駐點(diǎn)壓力、摩阻系數(shù)及駐點(diǎn)點(diǎn)熱流的結(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)的誤差,還列出了實(shí)驗(yàn)結(jié)果。可見,計(jì)算給出的真值,與經(jīng)認(rèn)真檢驗(yàn)的實(shí)驗(yàn)值相當(dāng)接近。這說明,本文的真值估算方法是正確有效的。

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

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

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

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

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

圖1還畫出了兩位真值準(zhǔn)確數(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 二位真值準(zhǔn)確的數(shù)據(jù)帶Fig.1 Thezoneof approximating to first twodigil number of truth value

3 結(jié) 論

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

[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計(jì)算結(jié)果的不確定度問題[J].空氣動(dòng)力學(xué)學(xué)報(bào),2008,26(1):47-49.

[22]惲起麟.風(fēng)洞實(shí)驗(yàn)[M].近代空氣動(dòng)力學(xué)叢書,國(guó)防工業(yè)出版社,北京,2000.

[23]程厚梅等.風(fēng)洞實(shí)驗(yàn)干擾與修正[M].近代空氣動(dòng)力學(xué)叢書,國(guó)防工業(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.