非參數(shù)化概率盒下隨機(jī)與認(rèn)知不確定性的分離式靈敏度分析

吳沐宸，陳江濤，夏侯唐凡，趙煒，劉宇，3，*

1．電子科技大學(xué) 機(jī)械與電氣工程學(xué)院，成都 611731

2．中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000

3．電子科技大學(xué) 系統(tǒng)可靠性與安全性研究中心，成都 611731

復(fù)雜工程系統(tǒng)及其數(shù)值計算模型，如：飛行器及其計算流體力學(xué)（Computational Fluid Dy?namics， CFD）數(shù)值模擬模型，其性能響應(yīng)受眾多輸入變量的影響［1］。工程設(shè)計人員需要從大量的輸入變量中篩選出對系統(tǒng)響應(yīng)有影響的變量，從而簡化計算模型并降低計算復(fù)雜度。靈敏度分析是一種量化輸入變量對輸出響應(yīng)不確定性貢獻(xiàn)度的方法，可實現(xiàn)對輸入變量的重要性排序。一方面，通過靈敏度分析可以將不重要變量固定為一個確定的常數(shù)，從而減小輸入變量維數(shù)；另一方面，通過降低重要輸入變量的不確定性可以大幅降低輸出響應(yīng)的不確定性，從而為復(fù)雜工程系統(tǒng)的穩(wěn)健設(shè)計提供重要依據(jù)。

靈敏度分析方法主要包括以下3類［2］：局部靈敏度分析（Local Sensitivity Analysis， LSA）、區(qū)域靈敏度分析（Regional Sensitivity Analysis，RSA）和全局靈敏度分析（Global Sensitivity Analysis， GSA）。其中，全局靈敏度分析研究輸入變量在整個取值域的不確定性對輸出響應(yīng)不確定性的影響程度［3-4］，相較于局部靈敏度分析的結(jié)果更為準(zhǔn)確，且適用于相對復(fù)雜的非線性系統(tǒng)。常用的全局靈敏度分析方法有掃描法［5］、微分法［6］、基于方差的靈敏度分析［7］、矩獨立靈敏度分析［8］、基于信息量的靈敏度分析［9］等。

值得注意的是，現(xiàn)有的全局靈敏度分析方法及其拓展方法都是建立在輸入變量只存在隨機(jī)不 確 定 性（Aleatory Uncertainty）［10］或 認(rèn) 知 不 確定性（Epistemic Uncertainty）［11］的條件上，對 于只含單一類型不確定性的輸入變量可以得到較為準(zhǔn)確的靈敏度結(jié)果。然而，由于對復(fù)雜工程系統(tǒng)的物理規(guī)律認(rèn)識不清、待估計參數(shù)樣本貧乏以及輸入?yún)?shù)固有的隨機(jī)性，隨機(jī)和認(rèn)知不確定性常常耦合在一起，稱之為混合不確定性［12］。概率盒 模型（如圖 1所示，圖中p為概率；X1~X3為參數(shù)變量）是一種典型的混合不確定性量化方法，包括：參數(shù)化概率盒（Parametric P-box）和非參數(shù)化概率盒（Non-Parametric P-box）。參數(shù)化概率盒由一族同類型的累積分布函數(shù)（Cumu?lative Distribution Function， CDF）組成，而非參數(shù)化概率盒由邊界分布包絡(luò)的不同類型CDF曲線組成，其包絡(luò)的分布函數(shù)的形式任意，并不局限于經(jīng)典的分布函數(shù)［12］，因此，在工程中應(yīng)用更加廣泛。本文主要針對非參數(shù)化概率盒，發(fā)展混合不確定性下的全局靈敏度指標(biāo)。

要解決非參數(shù)化概率盒下全局靈敏度分析問題，首先需要解決非參數(shù)化概率盒下混合不確定性傳播問題。非參數(shù)化概率盒下混合不確定性傳播與概率框架下不確定性傳播方法［13-14］不同，因為非參數(shù)化概率盒變量僅邊界分布已知，內(nèi)部包絡(luò)的CDF形式均未知。近年來，學(xué)者們提出了高效的非參數(shù)化概率盒混合不確定性傳播方法。Zhang等［15］提出了一種區(qū)間蒙特卡洛算法，也稱作雙層嵌套蒙特卡洛算法。外層對累積概率分布采樣得到輸入變量區(qū)間，內(nèi)層通過外層采樣區(qū)間內(nèi)尋優(yōu)獲得響應(yīng)概率盒的邊界分布。Li等［12］提出了一種基于累積分布函數(shù)離散化的不確定性傳播方法，計算響應(yīng)概率盒的前4階統(tǒng)計矩取值范圍，并利用Johnson分布擬合［16］和基于百分?jǐn)?shù)優(yōu)化的方式構(gòu)建響應(yīng)概率盒的邊界分布。本文通過雙層嵌套蒙特卡洛算法得到輸出響應(yīng)的區(qū)間值序列，并利用經(jīng)驗累積概率分布［17］（Empirical CDF）構(gòu)建響應(yīng)的邊界分布。

由于非參數(shù)化概率盒模型耦合了輸入變量的隨機(jī)和認(rèn)知不確定性，非參數(shù)化概率下的靈敏度分析挑戰(zhàn)較大。本文將輸入非參數(shù)化概率盒中的隨機(jī)不確定性和認(rèn)知不確定性分離，將非參數(shù)化概率盒下靈敏度分析轉(zhuǎn)化為單獨研究輸入非參數(shù)化概率盒的隨機(jī)不確定性（或認(rèn)知不確定性）對系統(tǒng)響應(yīng)的隨機(jī)不確定性（或認(rèn)知不確定性）的影響，從而指導(dǎo)后續(xù)系統(tǒng)參數(shù)的隨機(jī)或認(rèn)知不確定性降低甚至消除。然而，非參數(shù)化概率盒下構(gòu)建靈敏度指標(biāo)面臨2大挑戰(zhàn)：①如何分離非參數(shù)化概率盒變量的隨機(jī)與認(rèn)知不確定性；②如何構(gòu)建指標(biāo)分別衡量輸入隨機(jī)/認(rèn)知不確定性對響應(yīng)隨機(jī)/認(rèn)知不確定性的影響。針對挑戰(zhàn)①，本文分別提出了基于格點法和期望值法分離非參數(shù)化概率盒的隨機(jī)與認(rèn)知不確定性。針對挑戰(zhàn)②，本文將構(gòu)建概率盒最大方差和面積度量指標(biāo)分別衡量系統(tǒng)輸入隨機(jī)、認(rèn)知不確定性對輸出隨機(jī)、認(rèn)知不確定性的影響。最后，針對NACA0012翼型繞流問題，分析了來流參數(shù)和湍流模型設(shè)計參數(shù)中的混合不確定性，給出了相應(yīng)的靈敏度分析結(jié)果和工程指導(dǎo)意義。

1 非參數(shù)化概率盒及其不確定性傳播

1.1 非參數(shù)化概率盒

如 圖2所 示，概 率 盒（Probability-Box， Pbox）是概率分布與區(qū)間的廣義表現(xiàn)形式［18］。非參數(shù)化概率盒包絡(luò)的分布函數(shù)形式任意，相較參數(shù)化概率盒更具有一般性。非參數(shù)化概率盒FX可以定義為一個五元組和為概率盒的上、下邊界，概率盒包絡(luò)的分布函數(shù)FX(x)隸屬于概率盒空間FX，且滿足關(guān)系和VFX分別為概率盒包絡(luò)分布的均值和方差所屬區(qū)間［19］：

式中：ΩX為樣本域。一般來說，非參數(shù)化概率盒的邊界分布已知，邊界分布內(nèi)包絡(luò)無數(shù)條非參數(shù)化的分布函數(shù)［12］。如圖2所示，非參數(shù)化概率盒的 邊 界 分 布 為 分 段 高 斯 分 布 函 數(shù)FX；μ，σ2(x)=N(μ，σ2)：

式中：u(x)為階躍函數(shù)，x≥0時為1；否則為0。

圖2 非參數(shù)化概率盒示例Fig. 2 Illustration of non-parametric P-box

1.2 混合不確定性傳播

假 定 系 統(tǒng) 響 應(yīng) 函 數(shù)Y=g(X1，X2，…，Xn)，輸入X=[ ]X1，X2，…，Xn為非參數(shù)化概率盒向量，則系統(tǒng)輸出響應(yīng)Y的概率盒可通過雙層嵌套蒙特卡洛方法得到［20］。如圖3所示，雙層嵌套方法包含2層循環(huán)：外層通過蒙特卡洛采樣計算累積概率的笛卡爾積集p：

式 中：pi，j表 示 對Xi采 樣 得 到 的 第j個 概 率 樣 本。系統(tǒng)輸入的區(qū)間笛卡爾積集［X］可通過p進(jìn)行如下計算得到：

圖3 非參數(shù)化概率盒下雙層嵌套不確定性傳播流程圖Fig. 3 Double loop procedure with non-parametric Pboxes

為了避免分布擬合誤差［22］，本文采用經(jīng)驗分布函數(shù)構(gòu)建響應(yīng)概率盒的邊界分布［17］：

2 非參數(shù)化概率盒下隨機(jī)與認(rèn)知不確定性分離方法

2.1 基于期望值法的認(rèn)知不確定性分離

非參數(shù)化概率盒FX由邊界分布包絡(luò)的任意分布形式的分布函數(shù)構(gòu)成，分離FX的認(rèn)知不確定性指概率盒包絡(luò)的任意分布FX(x)退化為其期望值所 有EFX組 成 一 個 區(qū) 間為非 參數(shù)化概率盒分離認(rèn)知不確定性的結(jié)果。顯然，和即是概率盒上、下邊界分布的期望。如圖4所示，以式（3）定義的非參數(shù)化概率盒為例，分離認(rèn)知不確定性后，概率盒退化為區(qū)間數(shù)［?1.11，0.84］。

2.2 基于格點法的隨機(jī)不確定性分離

非參數(shù)化概率盒的隨機(jī)不確定性分離等價于將概率盒退化為內(nèi)部的一條分布函數(shù)［23］。根據(jù)分布函數(shù)性質(zhì)，分離隨機(jī)不確定性后得到的分布函數(shù)應(yīng)滿足：

圖4 非參數(shù)化概率盒認(rèn)知不確定性分離示意圖Fig. 4 Illustration for separating epistemic uncertainty of non-parametric P-box

和

式中：Ng為樣本數(shù)量；xi為樣本取值。非參數(shù)化概率盒包絡(luò)的任意分布應(yīng)同時滿足式（10）和式（11）。然而，由于式（10）線性不等式約束的系數(shù)矩陣為大型稀疏矩陣，直接求解式（10）和式（11）構(gòu)建的可行域不可行。因此，為了高效求解非參數(shù)化概率盒分離隨機(jī)不確定性后的分布函數(shù)，本文提出了基于格點法的非參數(shù)化概率盒隨機(jī)不確定性分離方法。

格 點(xi，F(xiàn)X(xi)) (i=0，1，…，Ng+1)由 樣本點xi與對應(yīng)累積概率FX(xi)組成，分別稱為x值和p值。如圖5所示，格點法包括4個主要步驟：

步驟1 從 區(qū) 間和區(qū)間中均勻隨機(jī)產(chǎn)生首末2個x值x0和xNg+1，通 過 輪 盤 賭 方 法［24］保 證0.5的概率得到和兩個區(qū)間邊界值，其余x值由等間距插值得到：

圖5 基于格點法的非參數(shù)化概率盒隨機(jī)不確定性分離Fig. 5 Grid point method for separating aleatory uncer?tainty of non-parametric P-boxes

此外，隨機(jī)產(chǎn)生Ng個在1～Ng范圍內(nèi)的整數(shù)，并用Ψ=?存放已經(jīng)產(chǎn)生的格點標(biāo)簽。

步驟2 按J中整數(shù)序列在區(qū)間FˉX(xj1)]中 均 勻 隨 機(jī) 產(chǎn) 生 第1個p值FX(xj1)，并將標(biāo)簽j1放入集合Ψ中。

步驟3 產(chǎn)生第i個p值時，在Ψ中搜索最鄰近ji且大于ji的標(biāo)簽jl和小于ji的標(biāo)簽js，則第i個p值FX(xji)的取值范圍：

式中：當(dāng)js不存在時FX(xjs)=0；當(dāng)jl不存在時FX(xjl)=1。執(zhí)行本步驟Ng次得到所有格點。

步驟4 將格點按標(biāo)簽1～Ng的順序順次排列，并 加入FX(x0)=0和FX(xNg+1)=1這2個 格點，得到格點集{(xi，F(xiàn)X(xi))|i=0，1，…，Ng+1}。格 點 法 可 在 區(qū) 間［0，1］均 勻 隨 機(jī) 產(chǎn) 生Ns個FX(x)，通過線性插值得到FX(x)對應(yīng)的樣本

步驟5 線性插值后得到的分布函數(shù)FX(x)可能超出概率盒邊界分布包絡(luò)的范圍。利用概率盒邊界分布對FX(x)進(jìn)行修正，若超出其取值范圍，則將移至距離其最近的邊界，保證了產(chǎn)生的分布函數(shù)FX(x)位于非參數(shù)化概率盒包絡(luò)的范圍之內(nèi)。圖5為格點法的示例，基于Ng=103個格點產(chǎn)生了3條分布函數(shù)。需要注意，分離隨機(jī)不確定性后非參數(shù)化概率盒退化為非參數(shù)化概率盒內(nèi)任意一條分布函數(shù)，因此結(jié)果不唯一。

3 隨機(jī)與認(rèn)知不確定性分離式靈敏度

3.1 基于面積的分離式靈敏度指標(biāo)

考慮輸入變量的認(rèn)知不確定性對響應(yīng)認(rèn)知不確定性的影響，本文采用面積指標(biāo)度量輸出響應(yīng)認(rèn)知不確定性，構(gòu)建基于面積的認(rèn)知不確定性分離式靈敏度指標(biāo)。圖6為基于面積的輸出概率盒認(rèn)知不確定性量化示意圖。

圖6 輸出響應(yīng)的認(rèn)知不確定性量化示意圖Fig. 6 Quantification of output epistemic uncertainty

不失一般性，令Xs為輸入非參數(shù)化概率盒變量，利用格點法分離Xs的隨機(jī)不確定性，使Xs退化為概率盒包絡(luò)的一條任意分布形式的分布函數(shù)，并通過雙層嵌套不確定性傳播方法構(gòu)建下面的優(yōu)化模型得到輸出響應(yīng)區(qū)間[Yj]：

式中：s≠i；s，i=1，2，…，n；j=1，2，…，Ns。利用式（7）和式（8）可構(gòu)建輸出響應(yīng)概率盒，進(jìn)而得到分離Xs的隨機(jī)不確定性后輸出響應(yīng)概率盒的面積：

式中：Seps-ep∈[0，1]反映了輸入變量Xs的認(rèn)知不確定性對輸出響應(yīng)認(rèn)知不確定性的影響，指標(biāo)越大表明Xs的認(rèn)知不確定性對輸出響應(yīng)的認(rèn)知不確定性的影響越大。

由于分離輸入非參數(shù)化概率盒變量中的隨機(jī)不確定性結(jié)果為概率盒包絡(luò)的任意一條分布函數(shù)，Seps-ep的取值不穩(wěn)定。假設(shè)通過格點法產(chǎn)生了nep條分布函數(shù)，記nep個靈敏度指標(biāo)結(jié)果為本文進(jìn)一步研究指標(biāo)隨nep的變化規(guī)律，計算S的均值和均值的95%置信區(qū)間，選取使S的均值收斂、置信區(qū)間寬度較小的nep值，并將對應(yīng)的均值作為指標(biāo)Seps-ep的估計。S的均值和95%置信區(qū)間（CI）為

式中：σS為nep個靈敏度指標(biāo)結(jié)果S的標(biāo)準(zhǔn)差：

3.2 基于最大方差的分離式靈敏度指標(biāo)

為了量化輸入非參數(shù)化概率盒的隨機(jī)不確定性對輸出隨機(jī)不確定性的影響，本文采用最大方差指標(biāo)度量輸出隨機(jī)不確定性，構(gòu)建基于最大方差的分離式靈敏度指標(biāo)。

圖7 輸出響應(yīng)概率盒的最大方差和最小方差示例Fig. 7 Maximum and minimum variances of output response

方差指標(biāo)已經(jīng)被廣泛應(yīng)用在靈敏度分析中，用于衡量輸出隨機(jī)不確定性的大小。對于輸出為概率盒而言，其方差為一個包含最大值與最小值的區(qū)間數(shù)［19］。本文將概率盒方差的上界定義為概率盒最大方差，將概率盒方差的下界定義為概率盒最小方差，如圖 7所示。然而，并非概率盒最大方差和最小方差都適用于衡量輸入非參數(shù)化概率盒的隨機(jī)不確定性對輸出隨機(jī)不確定性的影響。原因在于用于衡量輸入概率盒的隨機(jī)不確定性對輸出隨機(jī)不確定性的指標(biāo)必須滿足單調(diào)性，即：降低任意輸入概率盒的隨機(jī)不確定性，量化輸出概率盒的隨機(jī)不確定性指標(biāo)必須滿足單調(diào)不增。顯然概率盒最小方差指標(biāo)不滿足該單調(diào)性。如圖7所示，分離輸入概率盒的認(rèn)知不確定性后，輸出的概率盒更加“立陡”（圖 7中紅色虛線包含的概率盒），其隨機(jī)不確定性必然降低，但是通過計算發(fā)現(xiàn)，最小方差在分離輸入認(rèn)知不確定性前后始終為0，靈敏度指標(biāo)不變。因此，概率盒最小方差不能用作衡量輸入非參數(shù)化概率盒的隨機(jī)不確定性對輸出隨機(jī)不確定性的影響。而最大方差在分離輸入概率盒的認(rèn)知不確定性后會降低，能夠有效評判輸入的隨機(jī)不確定性對輸出隨機(jī)不確定性的影響。因此，本文提出建立基于概率盒最大方差的分離式靈敏度指標(biāo)。

求解概率盒的最大方差等價于求解概率盒的一條分布函數(shù)，并使其方差最大，本文稱該分布函數(shù)為“最大方差分布函數(shù)”。式（10）表明難以通過遍歷概率盒的任意分布求解其最大方差。因此本文提出了一種高效的最大方差求解方法。分布的方差可通過生成［0，1］間均勻分布的概率樣本，經(jīng)逆分布函數(shù)采樣獲得。令p(j)(j=0，1，…，Ns)為均勻分布概率樣本由小到大排列的順序統(tǒng)計量，則響應(yīng)的方差為

式中：y(j)為p(j)通過逆分布函數(shù)采樣得到的輸出樣本和為 樣 本集的函數(shù)：

式（22）表明，輸出響應(yīng)的方差是樣本y(m)的二次函數(shù)，而因此當(dāng)且僅當(dāng)全部樣本y(j)取值為其邊界時，才能得到概率盒的最大方差。考慮到最大方差分布函數(shù)需滿足單調(diào)遞增性以及離散程度盡可能大2個要求，如圖 8所示，最大方差分布應(yīng)由3段分布函數(shù)組成：（?。┊?dāng)y∈[y(0)，y(k)]時，輸出響應(yīng)概率盒的上邊界時，累積概率位于 區(qū) 間均 勻 分 布 函 數(shù)；（ⅲ）當(dāng)y∈[y(k+1)，y(Ns)]時，輸出響應(yīng)概率盒的下邊界至此，遍歷式（10）的可行域求解最大方差的問題簡化為求解跳躍點y(k)的問題。

圖8 輸出響應(yīng)概率盒最大方差示意圖Fig. 8 Maximum variance of output response

下面將闡述概率盒最大方差的求解過程：首先由式（6）產(chǎn)生了Ns+1個輸出響應(yīng)區(qū)間第ⅱ段的累積概率可通過端點累積概率p(k)和p(k+1)進(jìn)行線性插值得到，則各段的累積概率為

式中：Nii表示第ⅱ段含有的累積概率值個數(shù)，取決于第ⅱ段的累積概率p(k+1)?p(k)：

進(jìn)一步地，以式（27）為目標(biāo)函數(shù)構(gòu)建優(yōu)化模型，尋找方差最大時的跳躍點y(k)以及對應(yīng)的最大方 差

分離輸入非參數(shù)化概率盒變量Xs的認(rèn)知不確定性使概率盒FXs退化為概率盒邊界分布期望組成的區(qū)間數(shù)[EFˉXs(x)，E-F Xs(x)]，并利用雙層嵌套不確定性傳播計算輸出響應(yīng)區(qū)間[Yj]：

式 中：s≠i；s，i=1，2，…，n；j=1，2，…，Ns，利用式（7）和式（8）構(gòu)建輸出響應(yīng)概率盒并利用式（28）計算分離輸入非參數(shù)化概率盒變量Xs的認(rèn)知不確定性后輸出概率盒的最大方差。令表示輸入變量Xs的認(rèn)知不確定性分離前輸出響應(yīng)概率盒的最大方差，則基于最大方差的分離式靈敏度指標(biāo)定義為

其中，Sals-al∈[0，1]反映了輸入變量Xs的隨機(jī)不確定性對輸出響應(yīng)隨機(jī)不確定性的影響，指標(biāo)越大表明Xs的隨機(jī)不確定性對輸出響應(yīng)的隨機(jī)不確定性的影響越大。

4 算例分析

4.1 數(shù)值算例

本算例采用簡單的多項式方程驗證本文所提出分離式靈敏度分析方法的正確性。具體的，有如下響應(yīng)函數(shù)：

式中：X1~X3均為非參數(shù)化概率盒變量，如表1所示。表1中列舉了Case1~Case5共5種類型概率盒邊界分布，分別是高斯分布N(μ，σ2)、對數(shù)正態(tài)分布LN(μ，σ2)、伽瑪 分布Γ(a，b)、多峰高斯分布［14］以 及 包 含 前3種 分 布 形式的混合邊界分布。設(shè)定2層嵌套不確定性傳播方法外層采樣量為Ns=105，格點法產(chǎn)生Ng=103個格點，隨機(jī)選擇nep=600條分布函數(shù)曲線，分別計算5種邊界分布形式下的分離式靈敏度指標(biāo)的均值和Sals-al（s=1， 2， 3），結(jié)果如圖 9所示。從圖中可以看出，分離輸入隨機(jī)不確定性和分離認(rèn)知不確定性的靈敏度排序結(jié)果分離均為X3>X2>X1。輸入X3的隨機(jī)不確定性和認(rèn)知不確定性分別對輸出Y的隨機(jī)不確定性和認(rèn)知不確定性的影響最大，而輸入X1的隨機(jī)不確定性和認(rèn)知不確定性分別對輸出Y的影響都是最小的。

表1 輸入非參數(shù)化概率盒變量的邊界分布Table 1 Boundary distributions of input non-parametric P-box variables

圖9 5種邊界分布形式下分離式靈敏度指標(biāo)Fig. 9 Separating sensitivity analysis indices for five boundary distribution types

4.2 NACA0012翼型繞流模擬

圖10 NACA0012翼型繞流模擬網(wǎng)格劃分Fig. 10 Grid setup of flow over NACA0012 airfoil simulation

圖11 無量綱壁面距離求解Fig. 11 Dimensionless wall distance resolution

表2 CFD升阻系數(shù)預(yù)測與風(fēng)洞試驗對比Table 2 Comparison between CFD results and wind tunnel test data

本文提出的分離式靈敏度分析方法將應(yīng)用于NACA0012翼型關(guān)鍵氣動參數(shù)升阻比的預(yù)測，分析來流參數(shù)和湍流模型參數(shù)的混合不確定性對升阻比預(yù)測結(jié)果不確定性的影響。圖10為NACA0012二維繞流模擬網(wǎng)格劃分。利用ANSYS/ICEM CFD19.0軟件構(gòu)建總格點數(shù)為203×216的C型網(wǎng)格，其中翼型表面共有208個格點，近壁面第1層網(wǎng)格高度為2×10?6。翼型的弦長c=1 m，前緣和后緣分別距離上游流場、下游流場的距離為10c和20c。首先，給定自由來流馬赫數(shù)Ma∞=0.7，迎角α=1.55°，雷諾數(shù)Rec=9×106，來 流 靜 溫T=288.15 K，采 用Spalart-Allmaras一 方 程 湍 流 模 型［25］，利 用ANSYS/FLUENT19.0計算無量綱壁面距離Δy+。Δy+的計算結(jié)果如圖11所示，可以發(fā)現(xiàn)最大無量綱壁面距離為4.63，滿足Δy+<5的黏性底層的條件［26］，并預(yù)測升阻系數(shù)與風(fēng)洞試驗結(jié)果［27］對比如表 2所示。通過表 2發(fā)現(xiàn)，本文構(gòu)建的CFD模型預(yù)測結(jié)果與風(fēng)洞試驗結(jié)果誤差較小，該CFD模型可以用于計算翼型氣動參數(shù)。由于翼型繞流包含激波-邊界層分離、轉(zhuǎn)捩等復(fù)雜的空氣動力學(xué)現(xiàn)象，采用確定性的湍流模型參數(shù)往往會帶來較大的預(yù)測誤差［28］。同時，真實飛行過程中翼面繞流機(jī)制較為復(fù)雜，存在大氣紊流和陣風(fēng)而對自由來流馬赫數(shù)和迎角等來流參數(shù)產(chǎn)生一定干擾，因此來流參數(shù)也存在不確定性［29］?？紤]到來流參數(shù)受大氣擾流的影響，同時經(jīng)典湍流模型參數(shù)的標(biāo)定來源于試驗、假設(shè)或?qū)＜医?jīng)驗，因此，實際工程中，來流參數(shù)和湍流模型參數(shù)用混合不確定性變量描述更加貼切。本文將NACA0012繞流模擬的邊界條件設(shè)定為來流參數(shù)（自由來流馬赫數(shù)Ma∞和迎角α）和Spalart-Allmaras一方程湍流模型參數(shù)［25］（卡門常數(shù)κ和cv1等5個參數(shù)），輸出為升阻比CLCD。NACA0012翼型繞流模擬模型由ANSYS/FLUENT19.0構(gòu)建，采用基于壓力的耦合算法，梯度計算采用格林高斯節(jié)點方法，使用二階迎風(fēng)Roe格式計算無黏通量［30］。由于求解一次數(shù)值解的時間過長，不利于不確定性量化與傳播，因此本文采用Kriging模型［31］代替FLUENT求解器構(gòu)建數(shù)值模型，并利用風(fēng)洞試驗［27］測得的394組升阻比數(shù)據(jù)近似代替CFD仿真模型驗證代理模型的預(yù)測精度，計算平均絕對百分比誤差為7.73%。表 3為來流參數(shù)、湍流模型參數(shù)非參數(shù)化概率盒的邊界分布，其中湍流模型參數(shù)的試驗標(biāo)定和取值范圍由文獻(xiàn)［32］給定。本文采用的2層嵌套不確定性傳播方法外層采樣量［33］為Ns=105，得到輸出升阻比概率盒如圖12所示，并利用式（15）和式（28）計算輸出升阻比概率盒的面積和最大方差，分別為10.738 0和86.973 9。

表3 NACA0012繞流模擬邊界條件設(shè)置Table 3 Boundary condition setup for flow over NACA0012 airfoil simulation

下面構(gòu)建基于面積的分離式靈敏度指標(biāo)。格點法產(chǎn)生Ng=103個格點，利用式（14）進(jìn)行線性插值得到分離概率盒隨機(jī)不確定性后的分布函數(shù)曲線。如圖 13所示，假定nep=200～1 000，按式（19）~式（21）求解表 3中7個參數(shù)的nep個靈敏度指標(biāo)結(jié)果S的均值和95%置信區(qū)間。

圖12 輸出升阻比的概率盒Fig. 12 Output lift-to-drag ratio P-box

從圖 13中可以發(fā) 現(xiàn)，當(dāng)nep=600時，S的 均值收斂且置信區(qū)間寬度較小。因此，本算例將隨機(jī)選擇nep=600條分布函數(shù)曲線開展后續(xù)靈敏度分析。單獨分離表 3給出的7個參數(shù)的隨機(jī)不確定性，得到輸出升阻比概率盒如圖 14所示。構(gòu)建式（18）中的基于面積的分離式靈敏度指標(biāo)，以nep=600時S的 均 值S?eps-ep作 為 靈 敏 度 指 標(biāo) 的 估計，并給出均值的95%置信區(qū)間寬度WCISeps?ep，如表4所示。從表4中可以發(fā)現(xiàn)，自由來流馬赫數(shù)Ma∞和迎角α是對輸出響應(yīng)升阻比的認(rèn)知不確定性貢獻(xiàn)最大的2組模型輸入?yún)?shù)，表明采用陣風(fēng)減緩等措施以減小自由來流馬赫數(shù)的認(rèn)知區(qū)間，同時利用顫振主動抑制減小迎角的浮動范圍對降低升阻比預(yù)測結(jié)果的認(rèn)知不確定性有顯著影響。針對湍流模型參數(shù)，κ，cv1，cb1，cw2的影響均比較顯著，但相對來流參數(shù)對升阻比的認(rèn)知不確定性影響較弱，有必要重點考慮卡門常數(shù)κ和cv1的認(rèn)知不確定性對升阻比的認(rèn)知不確定性的影響。同時，分別單獨分離表 3給出的7個模型參數(shù)的認(rèn)知不確定性，得到輸出升阻比概率盒如圖15所示，建立式（30）中的基于最大方差的分離式靈敏度指標(biāo)Sals-al（s=1，2，…，7）。表5說明自由來流馬赫數(shù)Ma∞和迎角α依舊是最主要的2個模型參數(shù)，它們的隨機(jī)不確定性對升阻比的隨機(jī)不確定性影響最大。湍流模型參數(shù)中κ，cw2，cv1，cb1的影響均比較顯著，需要重點考慮卡門常數(shù)κ的隨機(jī)不確定性對升阻比的隨機(jī)不確定性的影響。

圖13 靈敏度指標(biāo)Seps-ep的均值和95%置信區(qū)間Fig. 13 Mean value and 95% confidence interval of sensitivity analysis indexSeps-ep

圖14 分離輸入?yún)?shù)的隨機(jī)不確定性前后輸出升阻比的概率盒Fig. 14 Lift-to-drag ratio before and after separating aleatory uncertainty of input parameters

表4 基于面積的分離式靈敏度分析結(jié)果Table 4 Separating sensitivity analysis results based on area metric

圖15 分離輸入?yún)?shù)的認(rèn)知不確定性前后輸出升阻比概率盒Fig. 15 Lift-to-drag ratio before and after separating epistemic uncertainty of input parameters

表5 基于最大方差的分離式靈敏度分析結(jié)果Table 5 Separating sensitivity analysis results based on maximum variance metric

綜合表 4和表 5的結(jié)果，可以得出來流參數(shù)的不確定性是影響升阻比預(yù)測結(jié)果不確定性的關(guān)鍵因素，湍流模型參數(shù)的不確定性對升阻比的影響相比來流參數(shù)較小但不可忽略，這是因為Spalart-Allmaras模型參數(shù)綜合影響了渦黏性系數(shù)的生成、擴(kuò)散和耗散，而渦黏性系數(shù)是基于雷諾平均方程一類湍流模型中?；牧髁康年P(guān)鍵參數(shù)，將影響通過數(shù)值模擬預(yù)測升阻比的結(jié)果。在湍流模型參數(shù)當(dāng)中，有必要重點關(guān)注卡門常數(shù)κ和參數(shù)cv1，cw2的不確定性對升阻比預(yù)測結(jié)果的不確定性的影響。本文所提方法指出了來流參數(shù)和一部分湍流模型參數(shù)的不確定性對升阻比預(yù)測結(jié)果的不確定性存在一定影響，對于靈敏度排序靠后的湍流模型參數(shù)，可以通過賦常數(shù)值等處理方法，降低數(shù)值模擬的計算量。例如，按表 4和表 5的靈敏度排序，從低到高將重要性排名較低的湍流模型參 數(shù)σ，cw2，cb1，cv1賦為其試驗標(biāo)定值 （如表 3所示），計算輸出升阻比概率盒的面積和最大方差，并計算和未賦常數(shù)時升阻比概率盒的面積和最大方差的相對誤差，結(jié)果如表6所示?？梢园l(fā)現(xiàn)，分別將湍流模型參數(shù)σ和cw2以及同時將σ和cw2賦為試驗標(biāo)定值時，升阻比概率盒的面積和最大方差相較未賦試驗標(biāo)定值的情況相對誤差都不超過3%。表明σ和cw2的認(rèn)知、隨機(jī)不確定性對升阻比的影響非常小，可以賦為常數(shù)。而其他參數(shù)賦為試驗標(biāo)定值時，輸出的相對誤差都超過3%，不能簡單賦為試驗標(biāo)定值處理。因此，本文所提出方法可以降低復(fù)雜系統(tǒng)的計算負(fù)擔(dān)。

表6 輸入?yún)?shù)賦為試驗標(biāo)定值時輸出概率盒的面積和最大方差以及相對誤差Table 6 Comparison of lift-to-drag P-box’s area and maximum variance after and before setting in?puts as constants

同時，對比同一組輸入?yún)?shù)的隨機(jī)和認(rèn)知不確定性對升阻比預(yù)測的影響結(jié)果如表 7所示。表 7說明同一湍流模型參數(shù)的隨機(jī)不確定性成分與認(rèn)知不確定性成分對系統(tǒng)的不確定性貢獻(xiàn)并不完全一致。因此，工程中需要將混合不確定性參數(shù)的隨機(jī)和認(rèn)知不確定性分離，分別研究其對輸出結(jié)果的影響，才能為后續(xù)混合不確定性降低提供決策依據(jù)。

表7 輸入?yún)?shù)的隨機(jī)和認(rèn)知不確定性成分對升阻比預(yù)測結(jié)果的影響對比Table 7 Comparison of effects of aleatory and epis?temic uncertainties on lift-to-drag ratio

此外，文獻(xiàn)［20］提出了一種不分離隨機(jī)和認(rèn)知不確定性下靈敏度分析方法。該方法利用巴氏距離統(tǒng)一量化輸出的隨機(jī)和認(rèn)知不確定性，通過巴氏距離的變化量化輸入?yún)?shù)對輸出概率盒的靈敏度結(jié)果。將本文所提出的2個分離式靈敏度指標(biāo)結(jié)果與巴氏距離結(jié)果進(jìn)行對比，如表 8所示?？梢园l(fā)現(xiàn)，文獻(xiàn)［20］的靈敏度排序與分離認(rèn)知不確定性的靈敏度排序相同，而與分離隨機(jī)不確定性的靈敏度排序并不相同。這是因為巴氏距離度量了概率盒上下邊界分布間距，側(cè)重于量化輸出概率盒的認(rèn)知不確定性。通過對比可以說明，分離輸入變量混合不確定性中的隨機(jī)與認(rèn)知不確定性成分更能區(qū)分輸入?yún)?shù)的隨機(jī)、認(rèn)知不確定性成分對輸出的貢獻(xiàn)差異。

表8 本文方法與文獻(xiàn)［20］中的靈敏度方法對比Table 8 Comparison of proposed method and reported method in Ref.［20］

5 結(jié) 論

1）提出了非參數(shù)化概率盒下隨機(jī)與認(rèn)知不確定性分離式靈敏度分析方法。分離輸入非參數(shù)化概率盒變量的隨機(jī)和認(rèn)知不確定性成分，分別研究輸入隨機(jī)和認(rèn)知不確定性對輸出隨機(jī)和認(rèn)知不確定性的影響。

2）構(gòu)建了非參數(shù)化概率盒下基于格點法的隨機(jī)不確定性分離方法和基于期望值的認(rèn)知不確定性分離方法，可有效地分離輸入?yún)?shù)混合不確定性中隨機(jī)和認(rèn)知不確定性成分。

3）提出了基于面積和概率盒最大方差的靈敏度指標(biāo)，說明了利用概率盒最大方差構(gòu)建靈敏度指標(biāo)的合理性，并通過與傳統(tǒng)方法對比說明構(gòu)建分離式靈敏度指標(biāo)的優(yōu)勢。