夏侯唐凡，陳江濤，邵志棟，吳曉軍，劉宇, 3,*

1. 電子科技大學 機械與電氣工程學院，成都 611731 2. 中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000 3. 電子科技大學 系統(tǒng)可靠性與安全性研究中心，成都 611731

計算流體力學 (Computational Fluid Dynamics, CFD) 是用數(shù)值方法對流體流動進行模擬和分析的學科，是計算力學的一個重要分支。近年來，隨著計算機技術(shù)迅猛發(fā)展，計算流體力學數(shù)值模擬技術(shù)已經(jīng)大規(guī)模地應用于航空航天、船舶、風力、水利等領(lǐng)域的復雜流體分析、設(shè)備性能參數(shù)評估、新產(chǎn)品氣動設(shè)計與優(yōu)化中。由于CFD技術(shù)成本低、周期短、能提供真實試驗無法模擬的條件等優(yōu)勢，為復雜產(chǎn)品如：飛行器、風力機、透平機械等的氣動設(shè)計帶來一場新的革命。例如，美國Boeing公司在波音787客機研制中廣泛地采用了CFD數(shù)值模擬技術(shù)，包括：高速機翼設(shè)計、發(fā)動機短艙設(shè)計、翼尖小翼設(shè)計等，其所用風洞試驗時間相比波音777和波音767分別減少了30%和55%，大幅減少了風洞試驗成本。

理論上，CFD模擬結(jié)果應與試驗結(jié)果進行對比，對模擬結(jié)果加以確認，即：模型確認 (Model Validation)，才能增加設(shè)計人員對CFD模擬結(jié)果的認可度。然而，一方面，真實試驗由于試驗條件、幾何尺寸、邊界條件等不同，導致試驗結(jié)果存在一定波動性與隨機性；另一方面，CFD模擬中也存在著大量的不確定性，包括：模型參數(shù)、數(shù)值離散、模型形式和模型預測偏差等不確定性，這些不確定性嚴重影響了CFD模擬結(jié)果的可信度。傳統(tǒng)的CFD模擬技術(shù)都是在確定性的數(shù)學模型、物理參數(shù)、邊界條件下進行數(shù)值模擬的。顯然，用確定性模擬結(jié)果與試驗結(jié)果進行對比是不合理的。采用確定性CFD模擬結(jié)果對的產(chǎn)品開展研制將造成潛在的失效風險。例如，NASA在超高速飛行器X43.A的試驗失敗后進行了氣動方面的分析，得出的結(jié)論是由于氣動設(shè)計對不確定性因素模擬不足使得設(shè)計的魯棒性不足，導致控制系統(tǒng)過高估計了設(shè)計冗余，進而引發(fā)了試驗失敗。因此，CFD模擬技術(shù)的模型確認中必須全面考慮模擬過程所面臨的各種不確定性以及這些不確定性對系統(tǒng)輸出響應的影響。

一般而言，根據(jù)不確定性的屬性不同，可以將不確定性分為隨機不確定性 (也稱固有不確定性或客觀不確定性) (Aleatory Uncertainty)和認知不確定性 (Epistemic Uncertainty)。前者表示物理現(xiàn)象中存在的固有隨機性，它無法通過增加樣本控制或減小這類隨機性；后者則是由于人們主觀認識不足、知識與數(shù)據(jù)的缺乏而導致無法精確地構(gòu)建物理模型，或是不能對某些因素/參數(shù)的不確定度用精確的概率分布進行準確描述。通常情況下，隨機不確定性可用概率分布或隨機場表征，而認知不確定性可用概率方法 (如：貝葉斯理論) 或非概率方法 (如：區(qū)間理論、D-S證據(jù)理論、概率盒(Probability-Box, P-Box) 理論、凸集模型、模糊理論等) 表征。當隨機和認知這2類不確定性共存時，被稱為混合不確定性。CFD模擬技術(shù)中，既存在模型參數(shù)或邊界條件的固有不確定性，也存在由于真實流體與物理模型不同、模型參數(shù)估計不準確、模型數(shù)值離散精度不夠等帶來的認知不確定性。

不確定性下的CFD模擬除了量化和分析所面臨的各種不確定性外，還要闡明各種不確定性如何影響系統(tǒng)的輸出響應，即不確定性傳播問題。以往的不確定性下的CFD模擬都是在隨機不確定性下通過隨機采樣或微擾動技術(shù)等方法與確定性的CFD模擬相結(jié)合實現(xiàn)的。近年來，基于譜分析的CFD模擬大量應用在隨機不確定性下的CFD模擬中。其中，基于多項式混沌方法(Polynomial Chaos Expansion)的CFD模擬被廣泛研究。Xiu和Karniadakis提出廣義多項式混沌方法，并將此方法應用于隨機不可壓槽道流動與方柱繞流問題。Knio和Maitre用嵌入式多項式混沌法(Intrusive Polynomial Chaos Expansion)研究了翼型的不可壓縮繞流問題。王曉東和康順介紹了多項式混沌展開式與流體力學方程的耦合過程，并采用此方法研究頂蓋驅(qū)動的層流方腔流動。然而，認知不確定性下的CFD模擬卻鮮有研究，其主要難點有：① 認知不確定性和隨機不確定性不同，無法用蒙特卡洛采樣的方式處理；② 認知不確定性傳播方式與隨機不確定性傳播方式有著本質(zhì)區(qū)別，一般需要采用優(yōu)化方式得到系統(tǒng)響應輸出的邊界。

模型確認度量(Model Validation Metric)是模型確認過程中對模型可信度定量評價的方法。模型確認度量方法在CFD模型確認過程中的作用如圖1所示。利用試驗數(shù)據(jù)修正后的CFD模型的模擬結(jié)果與試驗結(jié)果進行對比并用模型確認度量方法進行評價。當不滿足應用要求時，對CFD模型進行改進或重選，直到模型滿足要求為止。

圖1 CFD模型確認框架Fig.1 Model validation framework of CFD

根據(jù)系統(tǒng)響應輸出的分散性特征，模型確認度量分為：確定性下模型確認度量方法和不確定性下的模型確認度量方法。傳統(tǒng)的確定性模型確認度量方法忽略了數(shù)值模型和真實試驗數(shù)據(jù)中潛在的各種不確定性，如：模型參數(shù)不確定性、試驗測量不確定性、觀測數(shù)據(jù)不足以及模型輸入?yún)?shù)與試驗值不一一對應的現(xiàn)象。不確定性下的模型確認度量方法則是需要全面考慮模型和真實試驗數(shù)據(jù)中可能存在的不確定性，進而對比模型結(jié)果的統(tǒng)計分布與真實試驗數(shù)據(jù)結(jié)果分布的一致程度。

不確定性下的模型確認度量方法主要包括隨機不確定性框架下的模型確認度量方法和認知不確定性框架下的模型確認度量方法。本文以飛行器翼型繞流的CFD數(shù)值模擬為研究背景，主要介紹CFD模擬中面臨的各種不確定性因素，探討傳統(tǒng)的隨機不確定性下模型確認度量方法在CFD技術(shù)中的應用條件和適用范圍。著重介紹在認知不確定性下的模型確認度量方法，包括幾類重要的區(qū)間、概率盒下模型確認度量方法。最后，通過NACA0012翼型繞流問題說明CFD技術(shù)中考隨機和認知不確定性下的模型確認度量方法的有效性。

1 CFD數(shù)值模擬中的不確定性

在飛行器翼型繞流問題中，CFD數(shù)值模擬的可信度很大程度上取決于數(shù)值仿真模型能否準確地反映真實的流體運動過程。然而，由于真實物理過程的復雜性以及人們的認知偏差，CFD數(shù)值仿真模型不可避免地存在隨機和認知不確定性，主要有以下3方面來源：

1) 模型形式不確定性：這類不確定性屬于認知不確定性，產(chǎn)生于認識不充分或知識不完備情況下構(gòu)建的不準確的物理模型，主要集中在針對湍流流場的求解上。湍流是一種高度復雜的三維非穩(wěn)態(tài)、有旋的不規(guī)則運動，其包含的物理參數(shù)隨時間與空間發(fā)生隨機變化。CFD求解湍流的方法之一就是引入湍流模型，然而，現(xiàn)實中存在十余種不同形式的湍流模型，多個湍流模型的選擇就不可避免地引入模型形式的不確定性。一般而言，湍流模型包括了基于Boussinesq假設(shè)的渦黏性模型和基于雷諾應力輸運方程的二階矩封閉模型。前者包括Spalart-Allmaras、-、-、SST等湍流模型，但這些模型均忽略了湍流場應變的歷史效應，且模型性能與具體翼型相關(guān)。后者將Navier-Stokes (N-S)方程右端項(生成項、擴散項、耗散項)用平均流動的物理量和湍流的特征尺度表示，人為引入了各向同性等假設(shè)。CFD使用的湍流模型都引入了假設(shè)項，因此沒有任何一個模型能完全精確地模擬真實流體的運動。湍流模型各種假設(shè)項也從另一方面產(chǎn)生了模型形式的不確定性。然而，研究模型形式的不確定性計算量較大、處理方式復雜，相關(guān)的研究進展還非常緩慢?，F(xiàn)有處理模型形式不確定性的方法均考慮不同雷諾數(shù)來流，利用各種湍流模型展開研究，從而進行模型選擇。

2) 模型參數(shù)的不確定性：這類不確定性既包含邊界條件、幾何尺寸等參數(shù)的隨機不確定性也包含數(shù)值模型，如：湍流模型、參數(shù)的認知不確定性。首先，作為模型輸入?yún)?shù)的邊界條件，如：壓力、速度、流量、體積分數(shù)、雷諾數(shù)、攻角等物理量時常因為條件限制而無法測量或測量誤差而造成隨機不確定性。另外，由于制造誤差、工藝波動、翼型結(jié)冰等因素，無法用確定幾何描述真實翼型尺寸，從而出現(xiàn)幾何尺寸參數(shù)的隨機不確定性。另一方面，在數(shù)值模型中，模型參數(shù)存在著由于經(jīng)驗不足造成的認知不確定性。以湍流模型為例，Spalart-Allmaras模型中卡門常數(shù)的設(shè)定、-和SST模型中壁面普朗特數(shù)的設(shè)定等均存在主觀性。不同的流動特性影響著湍流模型參數(shù)的選取，不同的專家在設(shè)置這些參數(shù)時也是根據(jù)自身經(jīng)驗設(shè)定，因此很難給出一個確切的參數(shù)取值。已有文獻集中于研究模型參數(shù)不確定性如何影響CFD模擬輸出的結(jié)果。例如，陳江濤等同時考慮了Spalart-Allmaras模型中的9個模型參數(shù)，通過Sobol指標量化每個參數(shù)的不確定度對輸出不確定度的影響。王曉東和康順研究了流體黏性系數(shù)的不確定性對流場速度分布的影響，表明流體黏性系數(shù)對不確定性傳播具有阻滯作用。錢煒祺等研究了Spalart-Allmaras模型參數(shù)在不同工況下對某戰(zhàn)斗機升力系數(shù)和阻力系數(shù)的影響，結(jié)果表明不同工況下，湍流模型參數(shù)對CFD計算結(jié)果的影響規(guī)律不同。Erb和Hosder討論了3種常用的湍流模型，通過靈敏度分析研究了存在激波的情況下各個模型參數(shù)的不確定性對升力系數(shù)等響應量的影響，并將結(jié)果與低速和跨聲速情況下的結(jié)果進行了比較。

3) 模型數(shù)值離散不確定性：這類不確定性屬于認知不確定性，產(chǎn)生于對控制方程及邊界條件的離散化。CFD數(shù)值求解分為2種方式：瞬態(tài)求解和穩(wěn)態(tài)求解。對于瞬態(tài)求解，時間離散的有限精度將導致數(shù)值解與仿真精確解之間的誤差。為了能夠準確地捕捉流動的各個尺度，采用的網(wǎng)格尺度必須與最小的Kolmogorov耗散尺度相當，且計算域應大于最大的含能渦尺度。事實上，通過該方法得到的三維總網(wǎng)格數(shù)太大、不現(xiàn)實，而當網(wǎng)格密度不滿足條件時，就不可避免地帶來數(shù)值解與精確解之間的誤差，從而造成數(shù)值離散的不確定性。數(shù)值離散不確定性進一步可以分為：截斷誤差和迭代誤差。所謂迭代誤差是指在同一套網(wǎng)格上迭代終止時數(shù)值計算所得出的解與離散方程解的偏差。截斷誤差則是指不同精細的網(wǎng)格所計算的數(shù)值解與真實解之間的誤差?，F(xiàn)有的處理和分析模型數(shù)值離散不確定性方法主要是Richardson外推法。Richardson外推法旨在利用數(shù)值解關(guān)于精確解的泰勒展開冪級數(shù)，得到網(wǎng)格收斂指數(shù)(Grid Convergence Index, GCI)，并利用GCI估計網(wǎng)格離散誤差的不確定度。趙訓友等使用基于Richardson外推的GCI方法，研究了低雷諾數(shù)二維非對稱圓管繞流模型的升力系數(shù)和阻力系數(shù)，并給出它們的數(shù)值離散誤差。Xing和Stern比較了不同的GCI方法，提出了新的安全系數(shù)法，并將其應用于船舶水動力學。Eca和Hoekstra基于最小二乘估計改進了原有的GCI方法，并用4個不同領(lǐng)域的工程案例驗證了所提方法的優(yōu)良性能。

通常情況下，邊界條件、幾何尺寸等參數(shù)的隨機不確定性可用概率分布表征，而模型參數(shù)的認知不確定性可用主觀概率方法 (如：貝葉斯理論) 或非概率方法 (如：區(qū)間理論、D-S證據(jù)理論、概率盒理論、模糊理論等) 表征。通過不確定性傳播方法，可以計算CFD系統(tǒng)輸出響應的不確定性。由于CFD數(shù)值模擬的系統(tǒng)輸出帶有隨機和認知不確定性，必須擴展現(xiàn)有的模型確認度量方法以適用于各種不確定性下的CFD模擬模型確認問題。對于翼型氣動CFD仿真問題，由于試驗環(huán)境的變化，來流參數(shù)迎角、馬赫數(shù)等會發(fā)生隨機波動而具有隨機不確定性。因此，使用隨機不確定性下的模型確認度量方法進行模型確認更為合理。若進一步需要考慮由于測量儀器的精度問題帶來的來流參數(shù)的認知不確定性，此時應使用混合不確定性下的模型確認度量方法進行模型確認。

2 隨機不確定性下的模型確認度量方法

在CFD數(shù)值模擬中，模型確認度量是模型選擇和模型可信度評估的重要依據(jù)，能為數(shù)值仿真的模擬結(jié)果與真實試驗數(shù)據(jù)的一致性提供定量評判指標。僅當模型確認度量結(jié)果滿足應用要求的CFD模型才能用于分析和設(shè)計；不滿足要求的模型則需要進一步修正，甚至要對建模方法做較大改進或重新選擇其他模型。模型確認度量方法分為確定性模型確認度量方法與不確定性模型確認度量方法2類。傳統(tǒng)的確定性模型確認度量方法使用和均方根誤差 (Root Mean Squared Error, RMSE)，忽略了數(shù)值仿真模型和真實觀測數(shù)據(jù)中潛在的各種不確定性，如：模型參數(shù)不確定性、試驗測量不確定性、觀測數(shù)據(jù)不足等)以及模型輸入?yún)?shù)與試驗值不一一對應的現(xiàn)象。不確定性下的模型確認度量方法則是需要全面考慮模型和真實試驗數(shù)據(jù)中可能存在的不確定性，進而對比模型結(jié)果的統(tǒng)計分布與真實試驗數(shù)據(jù)結(jié)果分布的一致程度。以下將首先介紹不確定性下模型確認度量方法應具備的6個基本性質(zhì)。其次，介紹4種主流的不確定性模型確認度量方法及其應用條件和適用范圍。

模型確認度量是對模型輸出響應與試驗觀測數(shù)據(jù)之間一致性的定量度量。Liu等提出模型確認度量方法應該具備以下6個基本性質(zhì)：

1) 模型確認度量方法應該是定量的方法且具有客觀性。對于給定的模型輸出響應與試驗觀測數(shù)據(jù)集，不同決策者應該給出相同的模型確認結(jié)果，與決策者主觀偏好無關(guān)。

2) 確定是否接受某模型的標準應該與模型確認度量結(jié)果相獨立。

3) 模型確認度量方法需要盡可能考慮各種不確定性源。

4) 模型確認度量方法能夠提供與試驗數(shù)據(jù)量相關(guān)的置信度。

5) 模型確認度量方法應能區(qū)分不確定性大小不同的模型；換句話說，如果模型中引入過多的不確定性，模型確認度量結(jié)果不應認為該模型更好。

6) 確認度量不僅能夠?qū)δＰ晚憫驮囼炗^測進行“單點”比較，也能將“多點”試驗數(shù)據(jù)集成起來對模型的全局預測能力進行評估。

這6個性質(zhì)已成為模型確認度量方法優(yōu)劣判斷的基本準則。一個模型確認度量方法滿足上述性質(zhì)越多則越優(yōu)。下面將主要介紹4類主流的模型確認度量方法。

在隨機不確定性框架下，模型確認度量方法主要分為4種：① 經(jīng)典的假設(shè)檢驗(Classical Hypothesis Testing)；② 貝葉斯因子(Bayes Factor)；③ 頻率度量(Frequentist’s Metric)；④ 面積度量(Area Metric)。前兩者屬于假設(shè)檢驗的方法，能給出模型與真實試驗數(shù)據(jù)是否一致的結(jié)論以及置信度，但不可避免第一類錯誤(即“棄真”錯誤)和第二類錯誤(即“納偽”錯誤)；后兩者是基于距離的方法，能在統(tǒng)計意義上給出仿真模型與真實試驗數(shù)據(jù)的偏離程度，而并非直接判斷仿真模型與真實試驗數(shù)據(jù)是否一致。

經(jīng)典假設(shè)檢驗方法利用試驗數(shù)據(jù)構(gòu)造檢驗統(tǒng)計量，能夠以一定的置信度拒絕模型。它的特點是隨著試驗數(shù)據(jù)量的增加，對于預測與試驗結(jié)果偏離較大的模型，拒絕率會增加；而對于預測與試驗結(jié)果吻合的模型，拒絕率幾乎不變，因此經(jīng)典的假設(shè)檢驗有更高的機率拒絕偏差較大的模型，這對于試驗數(shù)據(jù)量較多的工程問題(如：管道的應力和應變試驗等)比較適用。然而，當可獲得的試驗數(shù)據(jù)僅有一個(如：高超聲速飛行器試驗等)，該方法不能分辨模型的優(yōu)劣，因為單個試驗數(shù)據(jù)無法計算樣本方差，此時基于單個試驗數(shù)據(jù)的模型確認結(jié)果可能不會拒絕預測結(jié)果偏差較大的模型。

與經(jīng)典的假設(shè)檢驗相反，貝葉斯因子方法給出接受模型的置信度，其計算公式為

(1)

式中：(,)是備擇假設(shè)下均值和方差的先驗概率密度函數(shù)，(,|data)是給定物理觀測值data的情況下均值和方差的后驗概率密度函數(shù)。使用貝葉斯因子方法需要知道關(guān)于試驗數(shù)據(jù)均值和方差的先驗信息，因此該方法非常適用于含有大量歷史試驗數(shù)據(jù)的工程應用；若無先驗信息可用，則人為設(shè)置無信息先驗。貝葉斯因子方法的取值范圍為∈[0,+∞)。當>1時，可以接受模型(反之則拒絕)。隨著試驗數(shù)據(jù)量的增加，對于預測與試驗結(jié)果吻合的模型，接受率會增加；對于預測與試驗結(jié)果偏離較大的模型，接受率會降低。即使只有一個試驗數(shù)據(jù)，貝葉斯因子方法也能分辨模型的優(yōu)劣，但當試驗數(shù)據(jù)量很少時，貝葉斯因子方法對備擇假設(shè)的先驗分布十分敏感，因此該方法對試驗數(shù)據(jù)量較多的工程問題也是非常適用。

頻率度量方法通過測量模擬中響應的平均值與試驗的平均值之間的距離，量化模擬與試驗數(shù)據(jù)的一致性。頻率度量距離的范圍是[0,+∞)，距離越小，表明數(shù)值模擬與真實物理過程的一致程度越高。當試驗數(shù)據(jù)缺乏時，該距離具有不確定性，可以通過置信區(qū)間來量化，置信區(qū)間可以表示為

(2)

面積度量使用響應的累積分布函數(shù)(CDF)與試驗數(shù)據(jù)的經(jīng)驗分布函數(shù)之間的面積來量化模擬與真實物理試驗的一致性，如圖2所示。2個分布之間的面積可以表示為

(3)

式中：()為響應的累積分布函數(shù)，()為試驗數(shù)據(jù)的經(jīng)驗分布函數(shù)。面積度量取值范圍是[0,+∞)，面積度量值越小，表明數(shù)值模擬

圖2 面積度量方法示意圖Fig.2 Illustration of the area metric

與真實物理過程的一致程度越高。隨著試驗數(shù)據(jù)量的增加，面積度量的結(jié)果趨向于真實的面積值。當模型參數(shù)不確定性增大時，面積度量的結(jié)果也會相應增大。因此，面積度量能夠區(qū)分含有不同不確定度的模型。面積度量在模型確認中的一個顯著優(yōu)勢是它能夠在多個確認點整合稀疏的試驗數(shù)據(jù)，以評估模型在特定目標區(qū)域的整體預測能力，因此，該方法適用于多個邊界條件下存在試驗數(shù)據(jù)的工程問題(如：渦輪機葉片翼型的風洞試驗等)。然而在試驗數(shù)據(jù)缺乏時，試驗數(shù)據(jù)的經(jīng)驗分布函數(shù)與其真實分布存在差異，將導致面積度量的結(jié)果不準確。

4類主流模型確認度量方法的優(yōu)劣對比如表1所示。從表中可以看出，面積度量方法較其他3種方法滿足更多的基本特征，且由于其能區(qū)分不同不確定度的模型和全局評估能力，近年來受到了廣泛關(guān)注?，F(xiàn)有的認知不確定性下模型確認度量方法也都是在面積度量方法上做出擴展。

表1 4類主流模型確認度量方法的優(yōu)劣對比

傳統(tǒng)的面積度量方法僅適用于一維或多維獨立的模型。工程中的模型往往是多維相關(guān)的，此時上述面積度量方法就不再適用了。所以，在面積度量方法的基礎(chǔ)上，Li等提出了基于概率積分轉(zhuǎn)換(Probability Integral Transformation, PIT)的指標法用于解決相關(guān)多輸出情況下的模型確認度量問題。趙亮和楊戰(zhàn)平通過構(gòu)造與試驗觀測數(shù)據(jù)相關(guān)的協(xié)方差矩陣，提出一種基于面積度量的多響應確認度量方法，能夠量化各響應數(shù)據(jù)間的相關(guān)性，并用熱傳導模型確認挑戰(zhàn)問題驗證了所提方法的有效性。趙錄峰等基于隨機變量的數(shù)字特征，構(gòu)建了由多輸出數(shù)學期望列陣與協(xié)方差矩陣組成的多輸出模型確認混合矩指標。PIT方法的優(yōu)點是在模型確認時考慮了多個輸出的相關(guān)性，但PIT指標法需要求解模型輸出的聯(lián)合累積分布函數(shù)，這在輸出維度很高時是很難準確求得的。為了避免求解高維聯(lián)合累積分布函數(shù)的問題，胡嘉蕊和呂震宙將該主成分分析與面積度量相結(jié)合，構(gòu)造了一種易于計算、穩(wěn)定性高的模型確認指標，克服了傳統(tǒng)多輸出模型確認度量方法中求解聯(lián)合累積分布函數(shù)的困難。張保強等引入馬氏距離對多響應模型確認問題進行降維，將多響應量的聯(lián)合累積分布函數(shù)轉(zhuǎn)化為單維的馬氏距離累積分布函數(shù)，大大降低了計算成本。Mahadevan和Rebba基于假設(shè)檢驗對多響應模型進行確認，提出了一種Box-Cox轉(zhuǎn)換方法簡化了多維概率密度函數(shù)的計算與構(gòu)造。

3 認知不確定性下的模型確認度量方法

本節(jié)將介紹考慮數(shù)值模擬或試驗數(shù)據(jù)存在認知不確定性時，傳統(tǒng)面積度量方法的擴展，包括區(qū)間和概率盒下的模型確認度量方法。

3.1 區(qū)間變量下的面積度量方法

實際工程中，由于試驗數(shù)據(jù)稀疏或不精確，往往難以獲得某個模型參數(shù)完整的統(tǒng)計信息，用區(qū)間變量描述模型參數(shù)的不確定性更為合理。肖釗等將區(qū)間變量轉(zhuǎn)化為均勻分布，利用均勻分布的累積分布函數(shù)量化了試驗區(qū)間響應與模擬區(qū)間響應，并用2條累積分布函數(shù)之間的面積作為模型確認度量，如圖3所示。模擬區(qū)間響應可以通過求解如下優(yōu)化問題獲得：

(4)

式中：表示模型不確定性參數(shù)矢量，區(qū)間變量

圖3 基于區(qū)間變量的面積度量Fig.3 Interval-based area metric

的不確定傳播問題成為在輸入變量空間為情況下，求模擬響應的區(qū)間值[,]的問題。試驗區(qū)間響應可以通過以下公式計算：

(5)

其中：min()和max()分別表示試驗數(shù)據(jù)得最小值和最大值，為試驗數(shù)據(jù)個數(shù)，試驗區(qū)間響應表示為[,]。

由此，文獻[46]中的面積度量計算公式為

(6)

式中：()和()分別為模擬區(qū)間響應、試驗區(qū)間響應的均勻分布累積分布函數(shù)。面積度量的取值范圍是[0,+∞]，越小表明數(shù)值模擬的結(jié)果越精確，如果數(shù)值模擬能夠完全反映真實的物理過程，則理論上為0。

3.2 考慮試驗數(shù)據(jù)量的區(qū)間面積度量方法

(7)

(8)

圖4 3種情況下的面積度量Fig.4 Three typical situations of area metric

(9)

3.3 概率盒的面積度量方法1

當同時存在隨機和認知不確定性時，不確定性量通常可以表示為概率盒，文獻[24]提出了一種用于比較概率盒的面積度量方法，如圖5深色陰影部分所示。該方法計算相同概率水平下試驗樣本數(shù)據(jù)概率盒的不確定性區(qū)間與數(shù)值模擬結(jié)果概率盒的不確定性區(qū)間的差值的絕對值的最小值，并在概率空間上對該最小值進行積分，即可得到數(shù)值模擬結(jié)果與真實試驗數(shù)據(jù)的偏離程度。當面積度量=0時，不能得出模型形式不確定性為零的結(jié)論，只能說明沒有證據(jù)證明試驗數(shù)據(jù)與模擬響應的不一致。因此，單一的面積度量嚴重低估了模型形式不確定性。文獻[24]克服了這一缺點，將概率盒的面積度量由單一值擴展到區(qū)間值。面積度量區(qū)間值的上下界定義為

(10)

圖5 概率盒的面積度量1Fig.5 First type of area metric under P-Box

(11)

式(11)表示當為區(qū)間中任意值、為區(qū)間中任意值時，|-|的最小值和最大值。

3.4 概率盒的面積度量方法2

趙錄峰等提出了一種新的概率盒面積度量方法，如圖6所示。該方法首先計算模擬響應與試驗數(shù)據(jù)概率盒的上邊界之間的面積，再計算模擬響應與試驗數(shù)據(jù)概率盒的下邊界之間的面積，2個面積之和為概率盒的面積度量。面積度量的計算公式為

(12)

由圖6和式(12)可以推斷當數(shù)值模型和真實物理模型完全一致時，模擬響應概率盒的上、下邊界曲線將與試驗數(shù)據(jù)概率盒的上、下邊界曲線重合，即面積度量值趨于零。如果數(shù)值模型和真實物理模型不一致且差異程度越大，那么模擬響應概率盒的上、下邊界曲線將與試驗數(shù)據(jù)概率盒的上、下邊界曲線差異越大，面積度量值也越大，反之亦然。因此，該度量方法可以客觀地評估數(shù)值模型描述真實物理模型的準確程度。

圖6 概率盒的面積度量2Fig.6 Second type of area metric under the P-Box

3.5 概率盒的面積度量方法3

McKeand等提出了一種平均曲線法，極大簡化了概率盒的面積度量計算，如圖7所示。該方法將數(shù)值模擬響應的概率盒轉(zhuǎn)換為單個累積分布函數(shù)(即平均曲線)，然后比較數(shù)值模擬響應的平均曲線與試驗數(shù)據(jù)的經(jīng)驗分布函數(shù)。面積度量定義為

(13)

求取概率盒的平均曲線的方法有很多，這里文獻[54]提出的2層嵌套式傳播算法，首先獲得概率盒的上、下邊界數(shù)據(jù)和，然后平均曲線的數(shù)據(jù)可以由式(14)獲得：

=05×(+)

(14)

圖7 概率盒的面積度量3Fig.7 Third type of area metric under P-Box

4 工程算例

4.1 NACA0012翼型繞流的模型確認問題

將上述的模型確認度量方法應用于典型機翼空氣動力學問題，可有效地評估這些方法的性能。本文研究的是二維可壓縮流動下的NACA0012翼型繞流問題(如圖8所示)，模型邊界條件是迎角和馬赫數(shù)，模型輸出是升力系數(shù)。NACA0012翼型繞流問題的數(shù)值仿真模型使用ANSYS/FLUENT19.0進行構(gòu)建。由于求解一次數(shù)值解的時間長，不利于進行不確定性的量化與傳播，故使用高斯過程(Gaussian Process, GP)模型代替FLUENT構(gòu)建的數(shù)值模型。

圖8 NACA0012翼型繞流網(wǎng)格劃分示意圖Fig.8 Illustration of flow meshing of NACA0012 airfoil

4.2 確定性模型確認度量方法的結(jié)果

為研究確定性模型確認度量方法的模型確認錯誤率，假設(shè)NACA0012翼型升力系數(shù)的試驗數(shù)據(jù)由式(15)獲得：

(15)

表2 確定性條件下預測模型的數(shù)學公式

在確定性度量的背景下，CFD數(shù)值模擬中不存在任何不確定性，因此模型1和模型2的迎角設(shè)為定值5°，模型1和模型2相應的升力系數(shù)分別為0.620 9和0.571 1。使用確定性度量指標RMSE進行模型確認。由于RMSE屬于確定性模型確認度量，未考慮CFD仿真中的各種不確定性，因此模型確認結(jié)果可能出錯(即模型確認結(jié)果表明模型2優(yōu)于模型1)。該算例使用式(15)生成1 000個試驗數(shù)據(jù)，使用每個試驗數(shù)據(jù)對模型1和模型2進行模型確認，計算得模型確認的錯誤率為35%。該算例表明使用確定性模型確認度量方法對CFD模型進行模型確認將有35%的概率選擇錯誤的模型。

4.3 隨機不確定性下模型確認度量方法的結(jié)果

為了分析隨機不確定下模型確定度量的性質(zhì)，假設(shè)NACA0012翼型的試驗數(shù)據(jù)通過式(15)獲得。其中，～(5,05)表示迎角服從正態(tài)分布，=05表示馬赫數(shù)真實值為0.5。此外，設(shè)計了2個預測模型，如表3所示。模型1與試驗模型完全一致，可認為是正確的預測模型；模型2的馬赫數(shù)為0.3，代表錯誤的預測模型。

表3 隨機不確定性下預測模型的數(shù)學公式

為了分析經(jīng)典假設(shè)檢驗方法的性質(zhì)，我們在置信度為95%進行1 000次假設(shè)檢驗，以獲得模型的拒絕率。表4結(jié)果表明隨著試驗數(shù)目的增加，模型1的拒絕率在7%上下波動，而模型2的拒絕率顯著提高，這表明增加試驗數(shù)據(jù)，有更大的幾率拒絕錯誤的模型。當試驗數(shù)目為20時，模型2的拒絕率高達99.7%，表明接受錯誤預測模型的概率僅為0.3%，且與確定性模型確認度量RMSE的結(jié)果相比，錯誤率顯著降低。此外，在相同的試驗數(shù)目下，模型2比模型1的拒絕率更高，這表明經(jīng)典的假設(shè)檢驗可以辨別模型的優(yōu)劣。

表4 模型拒絕率隨試驗數(shù)目的變化

為了研究貝葉斯因子的性質(zhì)，假設(shè)升力系數(shù)均值的先驗分布服從區(qū)間為[-,+]的均勻分布，標準差的先驗分布服從區(qū)間為[05,2]的均勻分布，和分別表示預測模型升力系數(shù)的均值和標準差。在相同的試驗數(shù)目下計算1 000次貝葉斯因子的值，統(tǒng)計>1的次數(shù)，/1 000為模型的接受率。表5展示了隨著試驗數(shù)目的增加，模型1的接受率增加，模型2的接受率降低，并且相同的試驗數(shù)目下，模型1比模型2的接受率更高。這表明增加試驗數(shù)據(jù)，有更大的幾率接受正確的模型，并且貝葉斯因子能夠辨別模型的優(yōu)劣。

表5 模型接受率隨試驗數(shù)目的變化

為了研究頻率度量的性質(zhì)，在相同的試驗數(shù)目下計算1 000次預測模型均值與試驗模型均值之間估計誤差的95%置信區(qū)間，區(qū)間上下界為1 000次 的平均值。并研究隨著試驗數(shù)目的增加，估計誤差的置信區(qū)間如何逼近真實誤差。模型1的模型形式與試驗模型完全一致，所以它們均值之間的真實誤差為0；模型2與試驗模型的馬赫數(shù)不同，它們均值之間的真實誤差為0.06。圖9給出了估計誤差的95%置信區(qū)間隨著試驗數(shù)目的變化趨勢，可以看出，隨著試驗數(shù)目的增加，模型1和模型2的區(qū)間寬度變窄，并始終包絡(luò)各自的真實誤差。此外，模型2的估計誤差的區(qū)間始終在模型1之上，這表明頻率度量也能區(qū)分模型的好壞。

圖9 估計誤差的95%置信區(qū)間隨試驗數(shù)目的變化Fig.9 95% confidence intervals of estimation error with the number of experimental data

為了研究面積度量的性質(zhì)，在相同試驗數(shù)目下計算了1 000次面積度量值，并繪制了面積度量值的概率密度函數(shù)，如圖10所示。模型1和模型2的真實的面積度量值分別是0和0.06。隨著試驗數(shù)目的增加，模型1的面積度量值的分布逐漸靠近其真實值0，但始終比0大(高估了模型偏差)。這是因為經(jīng)驗分布函數(shù)與累積分布函數(shù)之間的面積始終大于0。隨著試驗數(shù)目的增加，模型2的面積度量值的分布變窄并包絡(luò)其真實值0.06。此外，比較相同試驗數(shù)目下模型1與模型2的面積度量值的分布，模型2面積度量值的分布位于模型1的右側(cè)，這表明面積度量能夠區(qū)分模型的好壞。

圖10 面積度量值的分布Fig.10 Distributions of area metrics

4.4 基于區(qū)間變量的面積度量結(jié)果

為了研究基于區(qū)間變量的面積度量的有效性，假設(shè)NACA0012翼型升力系數(shù)的試驗數(shù)據(jù)由式(15)獲得。迎角∈[3°,6°]為一區(qū)間變量，=05表示馬赫數(shù)的真實值為0.5。此外，設(shè)計了與該物理試驗對應的2個預測模型(即生成模擬響應的模型)，如表6所示。

由表6可知，模型1與試驗模型完全一致，代表正確的預測模型，因此模型1確認結(jié)果的真實值為0；模型2的馬赫數(shù)為0.3，與試驗模型的馬赫數(shù)不一致，代表不準確的預測模型，模型2確認結(jié)果的真實值為0.044 4。模型1的確認結(jié)果應比模型2的確認結(jié)果好。

表6 認知不確定性下預測模型的數(shù)學公式

表7 基于區(qū)間變量的面積度量結(jié)果均值Table 7 Average area metric results by interval variables

4.5 考慮試驗數(shù)據(jù)量的區(qū)間面積度量結(jié)果

假設(shè)NACA0012翼型升力系數(shù)的試驗數(shù)據(jù)由式(15)獲得，其中迎角～(5,05)具有隨機不確定性，=05表示馬赫數(shù)的真實值為0.5。生成2個預測模型，如表3所示。模型1與試驗模型完全一致，代表正確的預測模型，模型1確認結(jié)果的真實值為0；模型2與試驗模型的馬赫數(shù)不一致，模型2確認結(jié)果的真實值為0.06。

如表8所示，相同試驗數(shù)目下，模型1的確認結(jié)果比模型2的確認結(jié)果好，這表明區(qū)間面積度量能夠判斷模型的優(yōu)劣。如圖12和圖13所示，大多數(shù)傳統(tǒng)面積度量值落在相應的區(qū)間面積度量所形成區(qū)間內(nèi)，這表明區(qū)間面積度量是對傳統(tǒng)面積度量的區(qū)間估計。對于模型1，區(qū)間面積度量的下界隨著試驗數(shù)目的增加趨近與真實值0，但始終大于0，這是因為經(jīng)驗分布函數(shù)與累積分布函數(shù)之間的面積始終大于0，這一結(jié)論與Liu等的結(jié)論一致。因此，區(qū)間面積度量在確認正確的模型時會高估真實的模型偏差。對于模型2，真實的面積度量值0.06都落在相應的區(qū)間面積度量所形成區(qū)間內(nèi)，這一結(jié)論能夠很好地應用于工程中。因為在工程實際中，仿真模型總存在模型不足，用試驗數(shù)據(jù)來確認仿真模型的情形類似于這里介紹的確認模型2(錯誤的預測模型)的情形。因此，在工程中進行隨機不確定性下的模型確認時，建議采用區(qū)間面積度量方法，它能給出真實模型偏差的區(qū)間估計。值得說明的是，區(qū)間面積度量方法是對傳統(tǒng)面積度量的改進，其認識不確定性是由于對比的試驗數(shù)據(jù)不足所造成的，而并非來源于仿真模型的參數(shù)。

表8 區(qū)間面積度量的模型確認結(jié)果均值

圖11 Ne= 5時模型1的模型確認Fig.11 Model validation of model 1 when Ne= 5

圖12 Ne= 5時區(qū)間面積度量與傳統(tǒng)面積度量的比較Fig.12 Comparison of interval-valued area metric and traditional area metric when Ne= 5

圖13 Ne= 20時區(qū)間面積度量與傳統(tǒng)面積度量的比較Fig.13 Comparison of interval-valued area metric and traditional area metric when Ne= 20

4.6 概率盒下面積度量結(jié)果

假設(shè)NACA0012翼型的試驗數(shù)據(jù)由式(15)獲得，其中=05表示馬赫數(shù)的真實值為0.5， 迎角～([3,6],05)表示迎角服從標準差為0.5的正態(tài)分布，其均值未知，但在區(qū)間[3°,6°]內(nèi)，這是一個典型的概率盒變量。升力系數(shù)的概率盒可通過外層隨機內(nèi)層認知的2層嵌套式傳播算法獲得。同時設(shè)置2個與試驗模型對照的預測模型，如表9所示。模型1與試驗模型完全一致，代表正確的預測模型，模型2與試驗模型馬赫數(shù)不同，代表錯誤的預測模型。

表9 概率盒框架下預測模型的數(shù)學公式

概率盒的面積度量方法1給出試驗概率盒與預測模型概率盒的區(qū)間面積度量，如果試驗數(shù)目足夠多，如圖14所示，則模型1與試驗模型的概率盒完全重合，根據(jù)式(10)，其真實的區(qū)間為[0,0.36]；如圖15所示，模型2與試驗模型的概率盒部分重合，根據(jù)式(10)，其真實區(qū)間為[0,0.40]。 概率盒的面積度量方法2計算的是試驗與預測概率盒上界之間的面積與下界之間的面積之和，由圖14和圖15所示，模型1確認結(jié)果的真實值為0，模型2確認結(jié)果的真實值為0.99。概率盒的面積度量方法3以試驗與預測概率盒平均曲線之間的面積作為度量，模型1與模型2確認結(jié)果的真實值分別為0和0.05。在相同的試驗數(shù)目下使用相同的概率盒的面積度量方法進行1 000次 模型確認，取其均值記錄于表10。

圖14 模型1與試驗模型的概率盒(試驗數(shù)目足夠多)Fig.14 P-Boxes of model 1 and experiential model with enough experiential data

圖15 模型2與試驗模型的概率盒(試驗數(shù)目足夠多)Fig.15 P-Boxes of model 2 and experiential model with enough experiential data

表10展示了3種概率盒的面積度量隨試驗數(shù)目增加的變化結(jié)果?？梢钥闯?，對于同一方法、同種模型，隨著試驗數(shù)目增加，概率盒的面積度量趨于真實值；對于同一方法、同一試驗數(shù)，模型1的確認結(jié)果均小于模型2的確認結(jié)果，這表明3種概率盒的面積度量方法均能夠判斷模型的優(yōu)劣。再看同一模型、同一試驗數(shù)下的3種不同的度量方法，可以看出方法3的度量值偏小，方法2的度量值偏大，而方法1的區(qū)間面積度量包含方法2與方法3的度量值，這表明概率盒的面積度量方法1能給出模型偏差的保守估計。因此，在工程中進行混合不確定性下的模型確認時，建議使用概率盒的面積度量方法1，因為它能給出模型偏差的保守估計。

表10 3種概率盒的面積度量結(jié)果Table 10 Validation results of three different P-Box area metrics

5 結(jié) 論

1) 分析了CFD數(shù)值模擬技術(shù)中存在的模型形式、模型參數(shù)和數(shù)值離散不確定性以及各種不確定性的屬性與表征方式。

2) 總結(jié)了單系統(tǒng)響應和多相關(guān)系統(tǒng)響應下主流的模型確認度量方法，給出了隨機不確定性框架下模型確認度量方法的應用條件和適用范圍。

3) 總結(jié)了認知不確定性框架下模型確認度量方法，主要包括區(qū)間和概率盒下的模型確認度量方法，并以NACA0012翼型繞流問題闡明了認知不確定性框架下模型確認度量方法的有效性。