吳 昊，楊 濤，豐志偉，葛健全，張青斌，黃 浩

（國防科技大學(xué)空天科學(xué)學(xué)院，湖南長沙 410073）

由于速度和突防方面的優(yōu)勢，高超聲速滑翔飛行器受到了各軍事強(qiáng)國的廣泛關(guān)注。鄰近空間內(nèi)長時間超高速飛行所帶來的嚴(yán)酷氣動加熱給飛行器防熱及結(jié)構(gòu)設(shè)計帶來了諸多挑戰(zhàn)［1］。文獻(xiàn)［2］指出，飛行器長時間受氣動加熱作用后，熱變形可能會引起氣動特性的變化，從而導(dǎo)致飛行彈道和姿態(tài)的變化，影響飛行器的整體性能。目前的數(shù)值仿真以確定性仿真為主，即需要給出精確的來流條件、邊界條件和物理模型。而氣動加熱熱流受來流條件、幾何尺寸、邊界條件和物理模型的影響，這些不確定性因素是普遍存在的，當(dāng)這些因素與實(shí)際情況存在偏差時，實(shí)際熱流就會偏離預(yù)示值。不確定性因素的來源主要可分為隨機(jī)不確定性（Aleatory Uncertainty）和認(rèn)知不確定性（Epistemic Uncertainty）兩大類。前者是由物理系統(tǒng)及其環(huán)境固有的隨機(jī)性產(chǎn)生的；而后者是由于人的認(rèn)識不足或者信息缺乏所造成的。氣動熱環(huán)境的不確定性量化評估對于高超聲速滑翔飛行器飛行試驗方案的確定具有重要參考意義。

Bettis和Hosder［3］對高超聲速再入流動的不確定度問題進(jìn)行了分析，研究混合不確定性輸入變量對氣動熱的影響。Brian A.Lockwood等進(jìn)一步研究了認(rèn)知和隨機(jī)不確定性混合傳播的優(yōu)化使用［4］，并發(fā)展出一種區(qū)間統(tǒng)計方法，結(jié)合Fay-Riddell駐點(diǎn)熱流公式證明其在流體仿真計算中的適用性；Lockwood在進(jìn)一步研究中提出一種多變量輸入情況下省時高效的基于梯度信息的代理模型［5］，研究了理想氣體和不平衡氣體條件下高超聲速黏流中熱流的不確定性。相比于國外，國內(nèi)對氣動熱不確定度量化評估的研究還處于起步階段。張偉等［6］運(yùn)用基于非侵入式多項式混沌（NIPC）的方法和基于Sobol指標(biāo)的方法對Apollo飛船返回艙開展了氣動熱不確定度量化分析和敏感性分析，得到了駐點(diǎn)和肩部熱流不確定度隨徑向距離變化的規(guī)律，為飛船返回艙氣動熱防護(hù)設(shè)計提供參考依據(jù)，證實(shí)了NIPC方法與Sobol敏感性分析方法在熱不確定性問題中的適用性，對進(jìn)一步用于高超聲速飛行器的魯棒優(yōu)化設(shè)計和可靠性分析提供了思路。

本文以某滑翔飛行器為研究對象開展氣動熱數(shù)值預(yù)示，選取來流速度、來流密度、來流溫度和端頭半徑4個不確定性變量，利用多項式混沌展開方法對氣動熱進(jìn)行不確定度量化和敏感性分析，為飛行試驗方案的確定提供理論依據(jù)。

1 氣動熱分析方法

高超聲速滑翔飛行器氣動熱環(huán)境可分解為駐點(diǎn)、飛行器前緣和大面積等部位，下面分別給出具體的計算方法。

1.1 駐 點(diǎn)

球頭駐點(diǎn)熱環(huán)境預(yù)測采用壓縮性修正的Fay-Riddle公式，即

式（3）中，壓縮性系數(shù)Z與駐點(diǎn)壓力和溫度相關(guān)，采用壓縮性系數(shù)函數(shù)表插值計算。

1.2 飛行器前緣

飛行器前緣的熱環(huán)境預(yù)測采用后掠圓柱理論，即

式（4）中，qy為翼前緣駐點(diǎn)線熱流，qs為與翼前緣相同半徑的駐點(diǎn)熱流，λe為有效后掠角，λ為后掠角，α為攻角，指數(shù)n取1.2。

1.3 大面積

飛行器大面積的熱環(huán)境預(yù)測采用流線法，計算中考慮變熵效應(yīng)。本文選用Zoby公式［7］計算熱流，該公式是基于Eckert參考焓的可壓縮平板摩擦系數(shù)關(guān)系得到的，計算結(jié)果與實(shí)驗數(shù)據(jù)和飛行數(shù)據(jù)均吻合較好。下面分別給出層流、湍流和轉(zhuǎn)捩區(qū)的熱流計算公式。

（1）層流熱流，即

式（5）中，qs為采用Fay-Riddle公式得到的駐點(diǎn)熱流。

（2）湍流熱流，即

式（6）中，c1，c2，c3，c4，m是速度剖面指數(shù)N的函數(shù)，分別為

在通常選用的1/7次速度分布中，N=7，但真實(shí)情況下N是變化的，由軸對稱噴管壁面實(shí)驗數(shù)據(jù)擬合出的N曲線［8］，即

式（8）中，Reθ是基于動量厚度的雷諾數(shù)，表達(dá)式為

本文根據(jù)式（8）和（9），采用迭代方法求解N與Reθ。

（3）轉(zhuǎn)捩區(qū)熱流，主要采用間歇因子方法計算轉(zhuǎn)捩區(qū)熱流，即

式（10）中，W為間歇因子，變化范圍是0＜W＜1，它與轉(zhuǎn)捩起始位置和轉(zhuǎn)捩結(jié)束位置有關(guān)，取

定義邊界層外緣當(dāng)?shù)乩字Z數(shù)等于轉(zhuǎn)捩雷諾數(shù)Retrb的位置為轉(zhuǎn)捩起始點(diǎn)，其位置為Strb。文獻(xiàn)［9］總結(jié)碳-酚醛端頭的飛行試驗數(shù)據(jù)，給出了與邊界層外緣馬赫數(shù)相關(guān)聯(lián)的公式，即

轉(zhuǎn)捩結(jié)束點(diǎn)位置為

2 不確定性量化評估方法和敏感性分析

本文采用多項式混沌（Polynomial Chaos，PC）法作為均值隨機(jī)展開，對滑翔飛行器氣動熱環(huán)境進(jìn)行不確定性量化分析。多項式混沌法基于不確定性的譜表示，相比傳統(tǒng)的蒙特卡洛仿真具有更高的效率。

2.1 多項式混沌不確定性量化方法

不確定性譜表示的基本思想是將響應(yīng)值或隨機(jī)函數(shù)α*分解為由確定性和隨機(jī)性分量表示的形式，即

其中，αi是確定性分量，Ψi是與第i模態(tài)對應(yīng)的隨機(jī)變量基函數(shù)，α*假設(shè)為獨(dú)立變量x和n維標(biāo)準(zhǔn)隨機(jī)變量ξ的函數(shù)。式（14）是一個無窮序列，實(shí)際使用時對其進(jìn)行截斷處理并保留部分輸出模態(tài)作為輸出。對于p階多項式混沌展開（Polynomial Chaos Expansion，PCE）和n維隨機(jī)變量，可計算總階數(shù)Nt，即

PCE方法的核心是確定展開系數(shù)αi，通?？墒褂们秩牖蚍乔秩氲姆绞接嬎愕玫健Ｇ秩胧椒椒ㄔ诶碚撋细鼮橹苯?，但是對于復(fù)雜問題該過程具有耗時長、代價大等缺點(diǎn)。而非侵入式方法，則不需要修改原確定性模型，容易構(gòu)建表示復(fù)雜模型的代理模型。

對于非侵入性多項式混沌（Nonintrusive Polynomial Chaos，NIPC），目前已發(fā)展了多種方法，其中在航空航天領(lǐng)域應(yīng)用最廣泛的是配點(diǎn)NIPC法。該方法首先將隨機(jī)響應(yīng)面或隨機(jī)函數(shù)近似表示為PCE；然后在隨機(jī)空間中選擇Nt個矢量，每個矢量即為一個采樣點(diǎn)，采用確定性模型評估這些采樣點(diǎn)，得到式（14）左邊；接下來，形成式（16）給出的Nt個方程的線性方程系統(tǒng)；最后求解隨機(jī)變量的譜模態(tài)。

需要說明的是，Nt是用于確定PCE系數(shù)的確定性采樣的最小數(shù)目。如果存在更多的采樣且線性無關(guān)，則系統(tǒng)為超定方程，可使用最小二乘方法求解。超過所需最小個數(shù)的采樣使用過采樣率（Oversampling Ratio，OSR）表示，定義為實(shí)際采樣與最小需要采樣的比值。一般而言，配點(diǎn)數(shù)量可以通過式（15）乘以O(shè)SR確定，PCE的精度依賴于配點(diǎn)數(shù)量。

2.2 Sobol靈敏度分析方法

Sobol靈敏度分析方法是基于方差的全局敏感性分析方法，通過計算一個或多個輸入變量對輸出方差的貢獻(xiàn)量，以此來評價該輸入變量在不確定性模型中的重要性。Sobol指標(biāo)需要通過Sobol重積分分解（以下簡稱Sobol分解）計算［10］。

假設(shè)通過PCE方法得到多項式函數(shù)f(x)，x=(x1，x2，…，xn)，定義域為In，通過Sobol分解可將該函數(shù)分解為如下形式，即

其中，各項按照以下公式分解得到

得到式（17）后，可求得總方差，即

各階偏方差為

總方差和偏方差滿足

各階Sobol指標(biāo)Si1…is的定義為

Sobol指標(biāo)提供了敏感性的度量，即，每個輸入不確定變量Si的貢獻(xiàn)以及混合貢獻(xiàn)Si1，...，is。輸入?yún)?shù)的總組合定義為部分Sobol指標(biāo)的和，即

即將所有與變量xi有關(guān)的Sobol指標(biāo)相加，當(dāng)n=4時，變量x1對總方差的總貢獻(xiàn)為

Sobol指標(biāo)可用于給出每個輸入不確定性變量對總方差貢獻(xiàn)的相對排序，且能夠考慮輸入變量和輸出指標(biāo)的非線性相關(guān)性。

3 仿真算例

以某高超聲速滑翔飛行器端頭為對象，開展氣動熱環(huán)境的不確定性量化評估和全局敏感性分析。其中，熱環(huán)境以輻射平衡熱流表示，邊界層外緣屬性按照激波后屬性計算［11］，邊界層黏性用Sutherland公式計算［12］，且壁面壓力可假設(shè)為邊界層外緣壓力。

飛行環(huán)境及端頭半徑如表1所示。分析中取4個隨機(jī)不確定變量：來流速度U∞、來流密度ρ∞、來流溫度T∞和端頭半徑Rn，由于飛行器加工過程中制造誤差和飛行環(huán)境的不確定性，來流條件的隨機(jī)波動是存在的。這些變量假設(shè)為關(guān)于某個均值的正態(tài)分布，變量分布和方差系數(shù)見表2。

表1 飛行環(huán)境及端頭半徑數(shù)值Tab.1 Aviation environment and the radius of curvature of the body

表2 隨機(jī)不確定性隨機(jī)變量分布Tab.2 Probability distribution of aleatory variables

由于該問題的不確定變量維數(shù)較小，可以使用最小二乘解來比較。對隨機(jī)不確定性在給定采樣規(guī)模下分析響應(yīng)面提供累積分布概率（CDF）的上下邊界。當(dāng)過采樣率取2時，能夠克服傳統(tǒng)配點(diǎn)法的不穩(wěn)定性，最大程度減少各組配點(diǎn)對響應(yīng)面的影響。圖1為過采樣率為2時熱流分布的P-box圖。相應(yīng)的95%置信區(qū)間為表3所示的下邊界2.5%概率水平到上邊界97.5%概率水平。

圖1 過采樣率取2時的熱流分布P-box圖Fig.1 P-box plots for heat flow(OSR=2)

圖2 隨機(jī)變量的Sobol指標(biāo)Fig.2 Sobol indices for aleatory variables

表3 熱流密度95%置信區(qū)間Tab.3 95%Confidence interval of heat flux for selected samples

4 總 結(jié)

本文以高超聲速滑翔飛行器為對象開展氣動熱環(huán)境數(shù)值預(yù)示，選取來流速度、來流密度、來流溫度和端頭半徑等4個不確定變量，采用NIPC方法對端頭半徑氣動熱進(jìn)行不確定量化評估和敏感性分析，得到如下結(jié)論：

（1）利用非侵入式多項式混沌方法與大量蒙特卡洛采樣所得不確定度區(qū)間吻合，表明利用NIPC方法能夠在縮短不確定性量化評估時間的同時，得到可信的結(jié)果。

（2）給定變量不確定條件下，端頭熱流的均值大約為305.05W/cm2，不確定度不小于7.03%。

（3）采用基于Sobol指標(biāo)的敏感性方法，對計算結(jié)果進(jìn)行了敏感性分析，結(jié)果表明來流速度變化和來流密度變化對熱流不確定度貢獻(xiàn)最大，端頭半徑對熱流的不確定度有一定影響，而來流溫度變化幾乎不產(chǎn)生影響。在高超聲速氣動熱預(yù)測中應(yīng)重點(diǎn)關(guān)注來流速度和來流密度的變化。