陳江濤，章超，吳曉軍，趙煒

1. 中國(guó)空氣動(dòng)力研究與發(fā)展中心 空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，綿陽(yáng) 621000

2. 中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所，綿陽(yáng) 621000

近些年來(lái)，隨著流場(chǎng)解算器、網(wǎng)格生成、前后置處理和高性能計(jì)算機(jī)等學(xué)科的快速發(fā)展，計(jì)算流體力學(xué)(Computational Fluid Dynamics，CFD)在工業(yè)領(lǐng)域得到了廣泛應(yīng)用，也發(fā)揮了日益重要的作用。但目前中國(guó)在CFD軟件體系化、工程化及大規(guī)模推廣應(yīng)用等方面還與國(guó)外先進(jìn)水平存在較大差距。因此中國(guó)2018年底啟動(dòng)了國(guó)家數(shù)值風(fēng)洞(National Numerical Windtunnel，NNW)工程。該工程由中國(guó)空氣動(dòng)力研究與發(fā)展中心牽頭建設(shè)，旨在建成擁有自主知識(shí)產(chǎn)權(quán)、中國(guó)開(kāi)放共享的空氣動(dòng)力數(shù)值模擬平臺(tái)，滿足航空航天等領(lǐng)域?qū)ψ灾鰿FD軟件產(chǎn)品的迫切需求。經(jīng)過(guò)兩年多的努力，NNW工程在模型算法、軟件研發(fā)等一系列關(guān)鍵技術(shù)上取得了長(zhǎng)足的進(jìn)展[1]。

在復(fù)雜工程問(wèn)題的數(shù)值模擬過(guò)程中，存在著多種來(lái)源的不確定性因素，包括網(wǎng)格、物理模型、數(shù)值離散方法等。在CFD中，數(shù)值、模型選擇和模型預(yù)測(cè)偏差是3類重要的不確定因素。數(shù)值離散引入數(shù)值誤差，但對(duì)于工程問(wèn)題數(shù)學(xué)模型的真實(shí)解通常未知，因此不能準(zhǔn)確估計(jì)數(shù)值離散誤差，只能采用數(shù)值不確定度進(jìn)行表征。CFD中模型選擇引入的不確定性非常普遍。對(duì)同一種自然現(xiàn)象，通常存在多種可能的物理模型模化，如各種湍流模型、亞格子模型和化學(xué)反應(yīng)模型等，不存在唯一的模型。人們?cè)噲D在可能的備選模型中選擇最好的模型，但由于工程問(wèn)題的復(fù)雜性以及使用者認(rèn)知的局限性，這種選擇本身也有一定的不確定性。對(duì)于某一個(gè)特定的模型，即使能夠準(zhǔn)確求解，其預(yù)測(cè)結(jié)果也可能與真實(shí)現(xiàn)象存在偏差。很多情況下這種偏差不能準(zhǔn)確地量化(因缺少試驗(yàn)數(shù)據(jù)或者試驗(yàn)存在顯著的不確定度)，轉(zhuǎn)化為模型預(yù)測(cè)偏差的不確定度。

不確定因素使CFD模擬結(jié)果也存在著不可忽視的不確定度。傳統(tǒng)CFD應(yīng)用過(guò)程中，從業(yè)者基于工業(yè)部門提供的數(shù)模，從經(jīng)驗(yàn)或理論分析出發(fā)生成分布直覺(jué)上合理的網(wǎng)格，使用常用的數(shù)學(xué)模型和數(shù)值格式運(yùn)行程序得到模擬結(jié)果。在整個(gè)過(guò)程中，沒(méi)有考慮到不確定性因素導(dǎo)致的結(jié)果的不確定度。這也給工業(yè)產(chǎn)品設(shè)計(jì)優(yōu)化、性能評(píng)估和模型確認(rèn)等活動(dòng)帶來(lái)潛在的風(fēng)險(xiǎn)?？茖W(xué)、定量地描述這些因素對(duì)模擬結(jié)果的影響已成為科學(xué)計(jì)算領(lǐng)域的重點(diǎn)方向。NASA[2]認(rèn)為，誤差和不確定度管理是未來(lái)CFD需要具備的六大基本能力之一。由于CFD中存在顯著的不確定因素，無(wú)法準(zhǔn)確地給出關(guān)注的氣動(dòng)力結(jié)果(給出的是某種模型、參數(shù)和離散方法下的數(shù)值結(jié)果)，但通過(guò)系統(tǒng)的不確定度量化工作可以給出氣動(dòng)力的置信區(qū)間或概率分布。這也可能改變CFD給工業(yè)部門提供數(shù)據(jù)的傳統(tǒng)方式。

針對(duì)模型選擇引入的不確定度，通常可用兩種方法進(jìn)行量化：調(diào)整參數(shù)法[3-5]和貝葉斯模型平均方法[6-8]。調(diào)整參數(shù)法在缺少試驗(yàn)數(shù)據(jù)的情況下，通過(guò)專家經(jīng)驗(yàn)給備選模型中的每一個(gè)模型分配相應(yīng)的概率，通過(guò)平均得到目標(biāo)量的均值和標(biāo)準(zhǔn)差，以此量化模型選擇的不確定度；該方法雖然取得了部分成功應(yīng)用，但是過(guò)分依賴于專家的經(jīng)驗(yàn)，不同專家可能給出不同的概率分配，這進(jìn)一步引入了不確定性。貝葉斯模型平均方法在獲得試驗(yàn)結(jié)果的情況下，通過(guò)構(gòu)建模型似然函數(shù)評(píng)價(jià)模型預(yù)測(cè)和試驗(yàn)結(jié)果的吻合程度，使用貝葉斯定理更新模型的概率；該方法能夠比較客觀地分配每個(gè)模型的概率，得到了廣泛的應(yīng)用。

貝葉斯模型平均方法依賴于模型預(yù)測(cè)結(jié)果分配模型后驗(yàn)概率，但在CFD領(lǐng)域通常無(wú)法獲得模型的真實(shí)解，CFD得到的是該模型在某種特定數(shù)值離散(格式、網(wǎng)格)下的數(shù)值解。使用數(shù)值解代替模型真實(shí)解會(huì)給模型后驗(yàn)概率計(jì)算帶來(lái)未知的不確定性，但現(xiàn)有應(yīng)用算例均未考慮數(shù)值離散對(duì)模型后驗(yàn)概率的影響。為考查這一影響，本文發(fā)展考慮數(shù)值離散誤差的貝葉斯模型平均方法，使用概率的上、下限表征數(shù)值離散和模型選擇兩者的總不確定度。首先對(duì)貝葉斯模型平均方法和數(shù)值離散誤差估計(jì)方法進(jìn)行介紹，然后提出考慮數(shù)值離散誤差的模型選擇不確定度量化方法，最后通過(guò)NACA0012翼型和CHN-T1模擬兩個(gè)算例分析湍流模型選擇中的不確定度。

1 貝葉斯模型平均方法

在存在多種可能模型的情況下引入模型概率的概念，定義為在備選模型集合中某一模型是最好模型的可信度。給定試驗(yàn)數(shù)據(jù)集D和包含n個(gè)備選模型的集合，使用貝葉斯定理將模型的先驗(yàn)概率更新為后驗(yàn)概率：

(1)

式中：P(Mk|D)為給定試驗(yàn)數(shù)據(jù)集D的情況下模型Mk的后驗(yàn)概率；P(D|Mk)為在給定試驗(yàn)數(shù)據(jù)D的情況下模型Mk的似然函數(shù)，即該模型預(yù)測(cè)結(jié)果為試驗(yàn)結(jié)果的概率，這項(xiàng)是貝葉斯模型平均方法的關(guān)鍵，需要建立計(jì)算數(shù)據(jù)和試驗(yàn)數(shù)據(jù)的概率關(guān)系，計(jì)為L(zhǎng)(Mk|D)；P(Mk)為模型Mk的先驗(yàn)概率，即在獲得試驗(yàn)數(shù)據(jù)之前根據(jù)專家意見(jiàn)或問(wèn)題分析分配的概率，當(dāng)對(duì)模型沒(méi)有明顯偏好時(shí)，可認(rèn)為每個(gè)模型的先驗(yàn)概率相等，即P(Mk)=1/n。通常做法是假定模型預(yù)測(cè)值和試驗(yàn)值之間的差異為均值為0的高斯隨機(jī)變量，即

(2)

(3)

(4)

此時(shí)

(5)

可以看到似然函數(shù)與模型預(yù)測(cè)值和試驗(yàn)的偏差成反比，與直覺(jué)想象一致。

考慮到試驗(yàn)數(shù)據(jù)有m個(gè)狀態(tài)點(diǎn)的情況，假定每個(gè)狀態(tài)都是獨(dú)立無(wú)關(guān)的，有

(6)

現(xiàn)可以通過(guò)式(1)給出模型的后驗(yàn)概率。綜合各個(gè)模型的后驗(yàn)概率得到一個(gè)新的融合模型，其預(yù)測(cè)值的概率描述為

(7)

由于每個(gè)模型輸出的概率表示P(y|Mk,D)均為獨(dú)立的高斯隨機(jī)變量，因此總的輸出也是高斯分布，其均值和方差為

(8)

在實(shí)際工程中，比較關(guān)注的是預(yù)測(cè)值在一定置信水平下的置信區(qū)間。由于輸出變量滿足高斯分布，根據(jù)其均值和方差得到一定置信水平下的置信區(qū)間非常方便，比如95%的置信區(qū)間為

(9)

值得注意的是，通過(guò)簡(jiǎn)單運(yùn)算可知，輸出變量的統(tǒng)計(jì)特性只與試驗(yàn)值和備選模型集合中每個(gè)模型的預(yù)測(cè)值有關(guān)。這也為第3節(jié)構(gòu)建目標(biāo)量的概率邊界提供了方便。

2 數(shù)值離散誤差估計(jì)方法

數(shù)值離散包括離散格式和計(jì)算域離散，必然引入數(shù)值誤差。所有針對(duì)模型的驗(yàn)證與確認(rèn)活動(dòng)(如模型確認(rèn)、模型平均等)都必須考慮數(shù)值誤差的影響，將模型結(jié)果和模擬結(jié)果區(qū)分開(kāi)。數(shù)值誤差估計(jì)就是通過(guò)一系列的數(shù)值計(jì)算估計(jì)模型方程真實(shí)解的情況。但對(duì)于工程問(wèn)題，準(zhǔn)確的數(shù)值誤差估計(jì)非常困難，此時(shí)數(shù)值誤差轉(zhuǎn)化成數(shù)值不確定度問(wèn)題。數(shù)值不確定度屬于認(rèn)知不確定度，因?yàn)閿?shù)值解和模型真實(shí)解的誤差是定值，但是由于認(rèn)知水平局限不知道模型真實(shí)解，所以也無(wú)法準(zhǔn)確估計(jì)離散誤差。通過(guò)一系列逐漸加密的自相似網(wǎng)格進(jìn)行Richardson插值分析是數(shù)值離散誤差估計(jì)的主要手段[9-11]。

Richardson插值是基于Taylor展開(kāi)的，假定離散解滿足

(10)

式中：fi為在該套網(wǎng)格上關(guān)注輸出量的離散解；f0為關(guān)注輸出量的真實(shí)解；α為未知常數(shù)；hi為網(wǎng)格的特征尺度，在二維計(jì)算中，通常使用N-1/2近似(N為單元總量)，在三維計(jì)算中，使用N-1/3近似；p為網(wǎng)格加密時(shí)數(shù)值誤差的收斂速度。該展開(kāi)有3個(gè)未知數(shù)(f0、α和p)，理論上只需要在3套滿足自相似和全計(jì)算域同一比例加密的網(wǎng)格上進(jìn)行計(jì)算即可。但是要求3套網(wǎng)格的離散解都位于真實(shí)解的漸進(jìn)收斂區(qū)域內(nèi)。假定粗、中、細(xì)3套網(wǎng)格計(jì)算得到的目標(biāo)變量分別為f3、f2、f1，3套網(wǎng)格的網(wǎng)格尺度分別為h3、h2、h1。通過(guò)簡(jiǎn)單的證明可知，只有當(dāng)0<(f1-f2)/(f2-f3)0，該展開(kāi)才是合理的。

但即使實(shí)際收斂速度p>0，也要考慮與理論收斂速度的差異。當(dāng)實(shí)際收斂速度與理論收斂速度差異較大時(shí)，會(huì)對(duì)誤差估計(jì)產(chǎn)生顯著的影響。CFD軟件使用二階的空間離散格式。密網(wǎng)格離散解的相對(duì)誤差定義為

(11)

通過(guò)簡(jiǎn)單的運(yùn)算可以得到相對(duì)誤差，也可以寫為

(12)

當(dāng)實(shí)際收斂速度0.5≤p≤2.0時(shí)，通過(guò)式(10) 進(jìn)行離散誤差估計(jì)是合理的；當(dāng)p>2.0時(shí)，數(shù)值誤差估計(jì)過(guò)小；當(dāng)0

因使用軟件的理論精度為二階，考慮使用基于一階、二階和一階/二階混合展開(kāi)進(jìn)行誤差估計(jì)，分別如式(13)～式(15)所示。

一階展開(kāi)為

fi=f0+αhi

(13)

二階展開(kāi)為

(14)

一階/二階混合展開(kāi)為

(15)

式中：α′、β和β′均為展開(kāi)常數(shù)。

式(10)、式(13)～式(15)給出了4種離散誤差展開(kāi)形式，具體使用哪種展開(kāi)需要根據(jù)離散解收斂性質(zhì)和展開(kāi)擬合誤差決定。

首先定義4種展開(kāi)的擬合誤差。當(dāng)使用ng(>3)套網(wǎng)格時(shí)，擬合誤差標(biāo)準(zhǔn)差的無(wú)偏估計(jì)分別為

(16)

式中:下標(biāo)1、2、12和p分別表示一階、二階、一階/二階混合展開(kāi)及使用式(10)的p階展開(kāi)。

當(dāng)只有3套網(wǎng)格時(shí)，σp和σ12均為0。

已知至少3套自相似和全計(jì)算域同一比例加密的網(wǎng)格上的離散解fi和相應(yīng)網(wǎng)格的全局網(wǎng)格尺度hi，對(duì)空間離散理論精度為二階的CFD軟件進(jìn)行數(shù)值離散誤差估計(jì)的流程如下：

1) 根據(jù)式(10)擬合得到實(shí)際收斂階p。

2) 確定擬合關(guān)系式：

① 當(dāng)0.5≤p≤2.0時(shí)，確定式為最終的展開(kāi)擬合關(guān)系式。

② 當(dāng)p>2.0時(shí),在一階展開(kāi)和二階展開(kāi)中，選擇擬合誤差最小的展開(kāi)為最終的擬合關(guān)系式。

③ 當(dāng)0

④ 當(dāng)p≤0時(shí),在一階展開(kāi)、二階展開(kāi)和一階/二階混合展開(kāi)中，選擇擬合誤差最小的展開(kāi)為最終的擬合關(guān)系式。

3) 根據(jù)選定的展開(kāi)關(guān)系式進(jìn)行最小二乘擬合，得到f0、擬合值fi,fit和擬合誤差標(biāo)準(zhǔn)差

求解的最小二乘問(wèn)題分別為

(17)

4) 根據(jù)擬合結(jié)果和實(shí)際收斂階p判斷誤差估計(jì)過(guò)程是否可靠，并給出安全因子Fs取值。

5) 給出考慮數(shù)值不確定度的模型真解的95%置信區(qū)間：

[fi-Δ,fi+Δ]

(18)

當(dāng)σ<Δ時(shí)，擬合誤差可接受，此時(shí)定義

Δ=Fs|fi-f0|+σ+|fi-fi,fit|

(19)

當(dāng)σ>Δ時(shí)，擬合誤差較大，此時(shí)定義

(20)

雖然步驟5)可以針對(duì)網(wǎng)格序列中的所有網(wǎng)格離散解，但一般在細(xì)網(wǎng)格上進(jìn)行，即模型真解的置信區(qū)間為

[f1-Δ,f1+Δ]

(21)

在實(shí)際應(yīng)用中，直覺(jué)上應(yīng)該是密網(wǎng)格的計(jì)算結(jié)果比粗網(wǎng)格更可信。因此在求解的最小二乘問(wèn)題中考慮網(wǎng)格尺度的權(quán)重wi，修改為

(22)

3 考慮數(shù)值離散誤差的模型選擇不確定度量化方法

貝葉斯模型平均方法給出了每個(gè)模型輸出為定值情況下的模型后驗(yàn)概率；數(shù)值離散誤差估計(jì)給出了模型真實(shí)解的置信區(qū)間。模型真實(shí)解是不確定的，在這個(gè)區(qū)間內(nèi)都有可能。此時(shí)模型的后驗(yàn)概率和輸出的統(tǒng)計(jì)特性都是不確定的。

為考慮數(shù)值不確定度對(duì)模型概率更新的影響，借鑒認(rèn)知隨機(jī)混合不確定性因素量化的二階概率框架[12]。該方法將認(rèn)知和隨機(jī)不確定性因素分離開(kāi)，通過(guò)嵌套的迭代方法構(gòu)建概率盒。其內(nèi)層循環(huán)是隨機(jī)不確定性因素，外層循環(huán)是認(rèn)知不確定性因素。當(dāng)固定認(rèn)知不確定性因素為某一確定取值時(shí)，只有隨機(jī)不確定性因素發(fā)揮作用，可以通過(guò)經(jīng)典的概率方法量化關(guān)注輸出量的隨機(jī)響應(yīng)，獲得其概率密度函數(shù)和累計(jì)分布函數(shù)。當(dāng)認(rèn)知不確定性因素在給定區(qū)間內(nèi)變化時(shí)，可以得到若干條可能的累積分布曲線。取該曲線簇的上、下限即為目標(biāo)量的概率盒。

提出考慮數(shù)值離散誤差的模型選擇不確定度量化方法，具體流程如下：

1) 針對(duì)備選模型集合中的每個(gè)模型估計(jì)其數(shù)值離散誤差，得到各自的置信區(qū)間。

2) 當(dāng)每個(gè)模型在其置信區(qū)間內(nèi)循環(huán)取值時(shí)，根據(jù)貝葉斯模型平均方法得到累積分布曲線簇，如圖1中樣本1～樣本5所示。

3) 通過(guò)簡(jiǎn)單的代數(shù)運(yùn)算得到目標(biāo)變量累積分布函數(shù)的上、下限，如圖1所示，即概率分布邊界，將其定義為目標(biāo)變量的概率盒。

圖1 NACA0012翼型累積分布曲線簇和分布邊界Fig.1 Clusters and bounds of cumulative distribution functions for NACA0012 airfoil

4) 通過(guò)概率盒可以分析一定置信水平下目標(biāo)變量的置信區(qū)間，此時(shí)區(qū)間的端點(diǎn)值不是固定值，而是一個(gè)區(qū)間值。

4 計(jì)算結(jié)果與分析

在CFD工程應(yīng)用中，通常使用湍流模型?；字Z應(yīng)力對(duì)平均流動(dòng)的影響。學(xué)者們從渦黏性假設(shè)或雷諾應(yīng)力輸運(yùn)方程出發(fā)，發(fā)展了數(shù)十種的湍流模型；如果加上各種湍流模型的修正形式，總共有上百種湍流模型。但現(xiàn)階段沒(méi)有對(duì)所有問(wèn)題普適的湍流模型。許多CFD使用者還是基于以往的經(jīng)驗(yàn)選擇某種湍流模型完成計(jì)算。

針對(duì)湍流模型選擇問(wèn)題，通過(guò)NACA0012翼型低速繞流和CHN-T1翼身組合體跨聲速繞流兩個(gè)算例證明第3節(jié)提到的模型選擇不確定度量化框架。需要指出的是，量化模型選擇不確定度時(shí)，需要計(jì)算所有可能的備選模型。出于計(jì)算可行性考慮，選用兩種使用最為廣泛的湍流模型，即SA(Spalart-Allmaras)[13]和Menterk-ωSST(Shear Stress Transfer)[14]模型。

計(jì)算使用的程序是課題組自行發(fā)展的MFlow[15-17]。該程序采用基于格心的有限體積法，能夠處理多種網(wǎng)格單元類型(六面體、四面體、三棱柱、金字塔以及使用幾何多重過(guò)程中融合產(chǎn)生的多面體)。單元內(nèi)使用線性重構(gòu)達(dá)到二階精度，梯度計(jì)算使用基于格點(diǎn)的格林高斯方法[18]，使用Venkatakrishnan限制器[19]抑制間斷附近的振蕩，使用Roe格式計(jì)算無(wú)黏通量，使用一階歐拉后差推進(jìn)到收斂狀態(tài)，使用局部時(shí)間步長(zhǎng)加速收斂過(guò)程，使用一階迎風(fēng)格式得到Jacobian通量。

4.1 NACA0012翼型低速繞流

第1個(gè)算例是NACA0012二維翼型繞流問(wèn)題，計(jì)算狀態(tài)為Ma=0.15、臨界雷諾數(shù)Rec=6.0×106，攻角為-3.99°～17.13°，試驗(yàn)結(jié)果可參考文獻(xiàn)[20]。計(jì)算使用了4套網(wǎng)格，均為四邊形單元。最密的一套網(wǎng)格使用了917 504個(gè)單元，C型網(wǎng)格拓?fù)?，流向和徑向單元?shù)分別是1 792 和512，標(biāo)示為Grid0。在密網(wǎng)格基礎(chǔ)上通過(guò)相鄰4個(gè)單元聚合形成下一級(jí)粗網(wǎng)格，網(wǎng)格尺度加倍。最粗的網(wǎng)格標(biāo)識(shí)為Grid3。

4套網(wǎng)格計(jì)算的升、阻力系數(shù)如圖2所示，直觀上感覺(jué)最粗網(wǎng)格計(jì)算結(jié)果與其他網(wǎng)格計(jì)算結(jié)果差異較大，特別是在中到大攻角下。在進(jìn)行數(shù)值離散誤差估計(jì)時(shí)，由于要同時(shí)協(xié)調(diào)所有網(wǎng)格計(jì)算結(jié)果完成數(shù)據(jù)擬合，最粗網(wǎng)格的結(jié)果導(dǎo)致產(chǎn)生非常大的置信區(qū)間，失去了研究問(wèn)題的意義。因此該算例使用較密的3套網(wǎng)格進(jìn)行后續(xù)分析。此時(shí)通過(guò)第2節(jié)的流程方法可以得到升、阻力系數(shù)的置信區(qū)間，如圖3所示。在小攻角時(shí)，置信區(qū)間很窄，說(shuō)明此時(shí)數(shù)值不確定度較?。辉诖蠊ソ菚r(shí)，置信區(qū)間變寬，說(shuō)明此時(shí)數(shù)值離散引入的不確定度較大，離散格式或網(wǎng)格會(huì)對(duì)模擬結(jié)果產(chǎn)生顯著影響。

圖2 NACA0012翼型不同網(wǎng)格下SA模型升力和阻力預(yù)測(cè)結(jié)果Fig.2 Lift and drag prediction results with SA model using different grids for NACA0012 airfoil

圖3 NACA0012翼型SA和SST模型升力和阻力預(yù)測(cè)的置信區(qū)間Fig.3 Confidence interval of lift and drag prediction with SA and SST models for NACA0012 airfoil

兩種湍流模型計(jì)算產(chǎn)生比較顯著的差異，如圖4所示。對(duì)于升力系數(shù)預(yù)測(cè)，直觀上無(wú)法判斷哪種模型更接近試驗(yàn)結(jié)果；對(duì)于阻力系數(shù)預(yù)測(cè)，SA模型結(jié)果更接近試驗(yàn)值，特別是在中到大攻角下。

圖4 NACA0012翼型密網(wǎng)格下SA和SST模型升力和阻力預(yù)測(cè)結(jié)果Fig.4 Lift and drag prediction results with SA and SST models using finest grid for NACA0012 airfoil

使用貝葉斯模型平均方法可以定量地給出兩種湍流模型的后驗(yàn)概率，如圖5所示。對(duì)于升力系數(shù)預(yù)測(cè)，大多數(shù)攻角下都是SST概率更大；對(duì)于阻力系數(shù)，大多數(shù)攻角下SA概率更大。

值得注意的是，兩個(gè)模型的后驗(yàn)概率在負(fù)攻角和對(duì)應(yīng)的正攻角區(qū)域不盡一致。這是因?yàn)樨惾~斯模型平均方法與模型預(yù)測(cè)結(jié)果和試驗(yàn)結(jié)果密切相關(guān)，數(shù)據(jù)的偏差會(huì)顯著影響后驗(yàn)概率。表1給出了在攻角-1.98°和2.00°的情況下使用最密網(wǎng)格時(shí)升力系數(shù)的計(jì)算結(jié)果和試驗(yàn)結(jié)果。此時(shí)試驗(yàn)結(jié)果或計(jì)算結(jié)果的微小誤差都會(huì)影響模型概率，不過(guò)兩種模型的計(jì)算結(jié)果差別不大，概率的偏差也不會(huì)過(guò)多影響輸出的統(tǒng)計(jì)特性。

表1 計(jì)算和試驗(yàn)升力預(yù)測(cè)對(duì)比

從圖5可以看出，網(wǎng)格對(duì)模型后驗(yàn)概率計(jì)算還是存在一定的影響。后驗(yàn)概率的差異必然會(huì)導(dǎo)致輸出統(tǒng)計(jì)性質(zhì)的不同。圖6給出了不同網(wǎng)格下，考慮模型選擇和模型預(yù)測(cè)偏差不確定度給出的關(guān)注輸出量的累積分布函數(shù)，可以估算出阻力系數(shù)一定置信水平下的置信區(qū)間，如在攻角為11.13°時(shí)，根據(jù)密網(wǎng)格結(jié)果估計(jì)阻力系數(shù)的95%置信區(qū)間為[0.011 04, 0.016 61]，根據(jù)粗網(wǎng)格結(jié)果估計(jì)的置信區(qū)間為[0.010 52, 0.018 51]?？煽闯鼍W(wǎng)格因素對(duì)置信區(qū)間估計(jì)影響顯著。

圖5 NACA0012翼型不同網(wǎng)格下SA模型升力和阻力預(yù)測(cè)的后驗(yàn)概率Fig.5 Posterior probabilities of lift and drag prediction using SA model with different grids for NACA0012 airfoil

圖6 NACA0012翼型不同網(wǎng)格得到的累積分布曲線Fig.6 Cumulative distribution functions using different grids for NACA0012 airfoil

如要考查模型本身的后驗(yàn)概率，拋掉數(shù)值離散的影響，需使用第3節(jié)介紹的方法。此時(shí)，SA和SST模型的真實(shí)解無(wú)法確定，只能假設(shè)其在置信區(qū)間內(nèi)都有可能，因此模型的后驗(yàn)概率也在一定的區(qū)間內(nèi)都有可能，如圖7中的灰色區(qū)域所示。

圖7 NACA0012翼型考慮數(shù)值離散誤差的SA模型升力和阻力預(yù)測(cè)的后驗(yàn)概率Fig.7 Posterior probabilities of lift and drag prediction using SA model considering discretization errors for NACA0012 airfoil

當(dāng)只能確定模型后驗(yàn)概率的區(qū)間時(shí)，輸出的統(tǒng)計(jì)特性也呈現(xiàn)一定的不確定度，可以通過(guò)概率的上、下限(如圖8所示)給出置信區(qū)間端值的區(qū)間。如在攻角為13.31°時(shí)，升力系數(shù)的95%置信區(qū)間的左端值為區(qū)間[1.335 5, 1.344 0]，右端值為區(qū)間[1.407 0, 1.439 5]；在攻角為11.13°時(shí)，阻力系數(shù)的95%置信區(qū)間的左端值為區(qū)間[0.010 92, 0.011 25]，右端值為區(qū)間[0.015 86, 0.017 35]。

圖8 NACA0012翼型累積分布函數(shù)邊界Fig.8 Bounds of cumulative distribution function for NACA0012 airfoil

4.2 CHN-T1跨聲速繞流

第2個(gè)算例是CHN-T1翼身組合體跨聲速繞流[21-22]。該外形是中國(guó)空氣動(dòng)力研究與發(fā)展中心自行設(shè)計(jì)的單通道運(yùn)輸機(jī)標(biāo)模，采用機(jī)翼/機(jī)身/平尾/立尾構(gòu)型。計(jì)算狀態(tài)為Ma=0.78，Rec=3.3×106，攻角為2.0°～4.0°，間隔0.5°。計(jì)算使用官方提供的非結(jié)構(gòu)混合網(wǎng)格，單元數(shù)目從6 518 485(標(biāo)示為Grid3)提升至162 253 613(標(biāo)示為Grid0)。具體網(wǎng)格信息可以參考文獻(xiàn)[23]。

通過(guò)第2節(jié)給出的離散誤差估計(jì)方法可得兩種湍流模型各自的95%置信區(qū)間，如圖9所示。置信區(qū)間給出了模型真實(shí)解的估計(jì)情況，與離散數(shù)值解有不同的含義。

圖9 CHN-T1外形SA和SST模型升力和阻力預(yù)測(cè)的置信區(qū)間Fig.9 Confidence interval of lift and drag prediction with SA and SST models for CHN-T1 configuration

通過(guò)貝葉斯模型平均方法可計(jì)算兩種模型的后驗(yàn)概率，如圖10所示。考慮數(shù)值離散的不確定度之后，無(wú)法直接判斷哪個(gè)模型的后驗(yàn)概率更大。

圖10 CHN-T1外形考慮數(shù)值離散誤差的SA模型升力和阻力預(yù)測(cè)的后驗(yàn)概率Fig.10 Posterior probabilities of lift and drag prediction using SA model considering discretization errors for CHN-T1 configuration

考慮數(shù)值離散、模型選擇、模型預(yù)測(cè)偏差3種不確定性因素后，得到升、阻力系數(shù)累積分布函數(shù)的上、下限，如圖11所示。此時(shí)可以估計(jì)其置信區(qū)間的范圍。比如對(duì)于升力系數(shù)，在2°攻角下，其95%的置信區(qū)間的左端值為區(qū)間[0.370 90, 0.372 75]，右端值為區(qū)間[0.433 0, 0.461 4]。

圖11 CHN-T1外形累積分布曲線簇和邊界Fig.11 Clusters and bounds of cumulative distribution curves for CHN-T1 configuration

5 結(jié) 論

將數(shù)值離散、模型選擇、模型預(yù)測(cè)偏差3種不確定性因素統(tǒng)一考慮，發(fā)展了考慮數(shù)值離散誤差的貝葉斯模型平均方法，給出關(guān)注目標(biāo)量累積分布函數(shù)的上、下限。該方法拋掉了傳統(tǒng)貝葉斯模型平均方法使用的數(shù)值離散解，通過(guò)嵌套方法構(gòu)建概率盒。

針對(duì)湍流模型選擇問(wèn)題，通過(guò)NACA0012翼型低速繞流和CHN-T1翼身組合體跨聲速繞流兩個(gè)算例給出了該方法實(shí)施的具體示例。

關(guān)于“已知試驗(yàn)結(jié)果，通過(guò)計(jì)算結(jié)果估算其置信區(qū)間的意義”問(wèn)題，需要指出的是模型確認(rèn)是確定模型在預(yù)期用途內(nèi)是否能表征真實(shí)流體系統(tǒng)的過(guò)程。當(dāng)存在高可信度試驗(yàn)結(jié)果時(shí)，可以比較計(jì)算與試驗(yàn)結(jié)果對(duì)模型進(jìn)行評(píng)價(jià)；當(dāng)不存在試驗(yàn)結(jié)果時(shí)，需要根據(jù)已有的模型確認(rèn)結(jié)論進(jìn)行推論。本文發(fā)展的方法可以在有試驗(yàn)結(jié)果的情況下，計(jì)算關(guān)注目標(biāo)量的置信區(qū)間。研究結(jié)論可以推廣到?jīng)]有試驗(yàn)結(jié)果的物理過(guò)程或狀態(tài)，通過(guò)有限的數(shù)值計(jì)算估算關(guān)注量的置信區(qū)間。這給工業(yè)部門合理使用CFD模擬結(jié)果提供了新的思路。