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

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

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

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

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

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

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

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

由于非參數(shù)化概率盒模型耦合了輸入變量的隨機(jī)和認(rèn)知不確定性，非參數(shù)化概率下的靈敏度分析挑戰(zhàn)較大。本文將輸入非參數(shù)化概率盒中的隨機(jī)不確定性和認(rèn)知不確定性分離，將非參數(shù)化概率盒下靈敏度分析轉(zhuǎn)化為單獨(dú)研究輸入非參數(shù)化概率盒的隨機(jī)不確定性（或認(rèn)知不確定性）對(duì)系統(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)知不確定性對(duì)響應(yīng)隨機(jī)/認(rèn)知不確定性的影響。針對(duì)挑戰(zhàn)①，本文分別提出了基于格點(diǎn)法和期望值法分離非參數(shù)化概率盒的隨機(jī)與認(rèn)知不確定性。針對(duì)挑戰(zhàn)②，本文將構(gòu)建概率盒最大方差和面積度量指標(biāo)分別衡量系統(tǒng)輸入隨機(jī)、認(rèn)知不確定性對(duì)輸出隨機(jī)、認(rèn)知不確定性的影響。最后，針對(duì)NACA0012翼型繞流問(wèn)題，分析了來(lái)流參數(shù)和湍流模型設(shè)計(jì)參數(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可以定義為一個(gè)五元組和為概率盒的上、下邊界，概率盒包絡(luò)的分布函數(shù)FX(x)隸屬于概率盒空間FX，且滿(mǎn)足關(guān)系和VFX分別為概率盒包絡(luò)分布的均值和方差所屬區(qū)間［19］：

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

式中：u(x)為階躍函數(shù)，x≥0時(shí)為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的概率盒可通過(guò)雙層嵌套蒙特卡洛方法得到［20］。如圖3所示，雙層嵌套方法包含2層循環(huán)：外層通過(guò)蒙特卡洛采樣計(jì)算累積概率的笛卡爾積集p：

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

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

為了避免分布擬合誤差［22］，本文采用經(jīng)驗(yàn)分布函數(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組 成 一 個(gè) 區(qū) 間為非 參數(shù)化概率盒分離認(rèn)知不確定性的結(jié)果。顯然，和即是概率盒上、下邊界分布的期望。如圖4所示，以式（3）定義的非參數(shù)化概率盒為例，分離認(rèn)知不確定性后，概率盒退化為區(qū)間數(shù)［?1.11，0.84］。

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

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

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

和

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

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

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

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

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

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

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

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

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

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

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

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

考慮輸入變量的認(rèn)知不確定性對(duì)響應(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ù)化概率盒變量，利用格點(diǎn)法分離Xs的隨機(jī)不確定性，使Xs退化為概率盒包絡(luò)的一條任意分布形式的分布函數(shù)，并通過(guò)雙層嵌套不確定性傳播方法構(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)知不確定性對(duì)輸出響應(yīng)認(rèn)知不確定性的影響，指標(biāo)越大表明Xs的認(rèn)知不確定性對(duì)輸出響應(yīng)的認(rèn)知不確定性的影響越大。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 算例分析

4.1 數(shù)值算例

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

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

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

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

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

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

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

從圖 13中可以發(fā) 現(xiàn)，當(dāng)nep=600時(shí)，S的 均值收斂且置信區(qū)間寬度較小。因此，本算例將隨機(jī)選擇nep=600條分布函數(shù)曲線(xiàn)開(kāi)展后續(xù)靈敏度分析。單獨(dú)分離表 3給出的7個(gè)參數(shù)的隨機(jī)不確定性，得到輸出升阻比概率盒如圖 14所示。構(gòu)建式（18）中的基于面積的分離式靈敏度指標(biāo)，以nep=600時(shí)S的 均 值S?eps-ep作 為 靈 敏 度 指 標(biāo) 的 估計(jì)，并給出均值的95%置信區(qū)間寬度WCISeps?ep，如表4所示。從表4中可以發(fā)現(xiàn)，自由來(lái)流馬赫數(shù)Ma∞和迎角α是對(duì)輸出響應(yīng)升阻比的認(rèn)知不確定性貢獻(xiàn)最大的2組模型輸入?yún)?shù)，表明采用陣風(fēng)減緩等措施以減小自由來(lái)流馬赫數(shù)的認(rèn)知區(qū)間，同時(shí)利用顫振主動(dòng)抑制減小迎角的浮動(dòng)范圍對(duì)降低升阻比預(yù)測(cè)結(jié)果的認(rèn)知不確定性有顯著影響。針對(duì)湍流模型參數(shù)，κ，cv1，cb1，cw2的影響均比較顯著，但相對(duì)來(lái)流參數(shù)對(duì)升阻比的認(rèn)知不確定性影響較弱，有必要重點(diǎn)考慮卡門(mén)常數(shù)κ和cv1的認(rèn)知不確定性對(duì)升阻比的認(rèn)知不確定性的影響。同時(shí)，分別單獨(dú)分離表 3給出的7個(gè)模型參數(shù)的認(rèn)知不確定性，得到輸出升阻比概率盒如圖15所示，建立式（30）中的基于最大方差的分離式靈敏度指標(biāo)Sals-al（s=1，2，…，7）。表5說(shuō)明自由來(lái)流馬赫數(shù)Ma∞和迎角α依舊是最主要的2個(gè)模型參數(shù)，它們的隨機(jī)不確定性對(duì)升阻比的隨機(jī)不確定性影響最大。湍流模型參數(shù)中κ，cw2，cv1，cb1的影響均比較顯著，需要重點(diǎn)考慮卡門(mén)常數(shù)κ的隨機(jī)不確定性對(duì)升阻比的隨機(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é)果，可以得出來(lái)流參數(shù)的不確定性是影響升阻比預(yù)測(cè)結(jié)果不確定性的關(guān)鍵因素，湍流模型參數(shù)的不確定性對(duì)升阻比的影響相比來(lái)流參數(shù)較小但不可忽略，這是因?yàn)镾palart-Allmaras模型參數(shù)綜合影響了渦黏性系數(shù)的生成、擴(kuò)散和耗散，而渦黏性系數(shù)是基于雷諾平均方程一類(lèi)湍流模型中?；牧髁康年P(guān)鍵參數(shù)，將影響通過(guò)數(shù)值模擬預(yù)測(cè)升阻比的結(jié)果。在湍流模型參數(shù)當(dāng)中，有必要重點(diǎn)關(guān)注卡門(mén)常數(shù)κ和參數(shù)cv1，cw2的不確定性對(duì)升阻比預(yù)測(cè)結(jié)果的不確定性的影響。本文所提方法指出了來(lái)流參數(shù)和一部分湍流模型參數(shù)的不確定性對(duì)升阻比預(yù)測(cè)結(jié)果的不確定性存在一定影響，對(duì)于靈敏度排序靠后的湍流模型參數(shù)，可以通過(guò)賦常數(shù)值等處理方法，降低數(shù)值模擬的計(jì)算量。例如，按表 4和表 5的靈敏度排序，從低到高將重要性排名較低的湍流模型參 數(shù)σ，cw2，cb1，cv1賦為其試驗(yàn)標(biāo)定值 （如表 3所示），計(jì)算輸出升阻比概率盒的面積和最大方差，并計(jì)算和未賦常數(shù)時(shí)升阻比概率盒的面積和最大方差的相對(duì)誤差，結(jié)果如表6所示?？梢园l(fā)現(xiàn)，分別將湍流模型參數(shù)σ和cw2以及同時(shí)將σ和cw2賦為試驗(yàn)標(biāo)定值時(shí)，升阻比概率盒的面積和最大方差相較未賦試驗(yàn)標(biāo)定值的情況相對(duì)誤差都不超過(guò)3%。表明σ和cw2的認(rèn)知、隨機(jī)不確定性對(duì)升阻比的影響非常小，可以賦為常數(shù)。而其他參數(shù)賦為試驗(yàn)標(biāo)定值時(shí)，輸出的相對(duì)誤差都超過(guò)3%，不能簡(jiǎn)單賦為試驗(yàn)標(biāo)定值處理。因此，本文所提出方法可以降低復(fù)雜系統(tǒng)的計(jì)算負(fù)擔(dān)。

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

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

表7 輸入?yún)?shù)的隨機(jī)和認(rèn)知不確定性成分對(duì)升阻比預(yù)測(cè)結(jié)果的影響對(duì)比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)知不確定性，通過(guò)巴氏距離的變化量化輸入?yún)?shù)對(duì)輸出概率盒的靈敏度結(jié)果。將本文所提出的2個(gè)分離式靈敏度指標(biāo)結(jié)果與巴氏距離結(jié)果進(jìn)行對(duì)比，如表 8所示?？梢园l(fā)現(xiàn)，文獻(xiàn)［20］的靈敏度排序與分離認(rèn)知不確定性的靈敏度排序相同，而與分離隨機(jī)不確定性的靈敏度排序并不相同。這是因?yàn)榘褪暇嚯x度量了概率盒上下邊界分布間距，側(cè)重于量化輸出概率盒的認(rèn)知不確定性。通過(guò)對(duì)比可以說(shuō)明，分離輸入變量混合不確定性中的隨機(jī)與認(rèn)知不確定性成分更能區(qū)分輸入?yún)?shù)的隨機(jī)、認(rèn)知不確定性成分對(duì)輸出的貢獻(xiàn)差異。

表8 本文方法與文獻(xiàn)［20］中的靈敏度方法對(duì)比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)知不確定性對(duì)輸出隨機(jī)和認(rèn)知不確定性的影響。

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

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