張偉， 王小永， 于劍， 閻超

(北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院, 北京 100083)

隨著計算機技術(shù)和數(shù)值計算方法的發(fā)展，計算流體力學(xué)[1]在工程中起到了越來越重要的作用。目前大部分的數(shù)值模擬都是確定性的，即在計算中需要給出確定的來流條件、幾何尺寸、邊界條件以及精確的物理模型。然而高度復(fù)雜的航天飛行器的分析和設(shè)計中不確定性因素是普遍存在的，例如幾何模型在載荷作用下的變形和大氣來流條件變化等。不確定性因素可分為偶然不確定度(aleatory uncertainty)和認(rèn)知不確定度(epistemic uncertainty)兩大類[2]。偶然不確定度是由于事物固有的隨機性導(dǎo)致的不確定度；而認(rèn)知不確定度是由于缺乏理論知識導(dǎo)致的不確定度。不確定度量化(uncertainty quantification)分析對復(fù)雜航天系統(tǒng)的魯棒設(shè)計和可靠性分析是非常重要的，并對飛船返回艙的熱防護系統(tǒng)(Thermal Protection System，TPS)設(shè)計具有重要參考意義[3]。

飛船返回艙在再入大氣層過程中飛行高度逐漸降低，大氣來流等條件也在不斷發(fā)生變化。Champion[4]的研究表明高度為48～69 km的大氣密度波動約為5%。因此為得到可靠的氣動熱數(shù)值預(yù)測，必須考慮來流參數(shù)的不確定性和由此帶來的氣動熱預(yù)測的不確定度。目前已經(jīng)發(fā)展出很多不確定度量化方法，如已經(jīng)在工程中應(yīng)用較多的蒙特卡羅(Monte Carlo)方法[5-6]，這種方法的主要缺陷是需要大量樣本，對計算資源消耗較大。近年來，采用多項式混沌(Polynomial Chaos,PC)方法在CFD不確定分析得到了廣泛應(yīng)用。多項式混沌方法又可分為嵌入式和非嵌入式[7]。本文采用的是非嵌入式多項式混沌(Non-Intrusive Polynomial Chaos,NIPC)方法，即選擇一組正交多項式作為空間的無限基函數(shù)，將變量分解為確定和隨機兩部分從而構(gòu)造出多項式混沌展開式。該方法相比于嵌入式多項式混沌方法只需把CFD求解器當(dāng)做黑盒子，基于確定性的解來近似多項式混沌展開式的系數(shù)。在NIPC方法的基礎(chǔ)上，可采用Sobol指數(shù)來衡量輸入變量的敏感性。

目前，NIPC方法被廣泛應(yīng)用到CFD不確定性分析中。Hosder等[7-8]將NIPC方法應(yīng)用于求解楔角變化的隨機斜激波問題和來流馬赫數(shù)與迎角變化的三維NACA65A004機翼氣動力問題，研究了單個隨機輸入變量和多維不確定性輸入變量的情況，結(jié)果表明NIPC方法具有較高的效率和準(zhǔn)確性。Loeven和Bijl[9]開展了NACA5412翼型數(shù)值模擬的幾何不確定度量化分析，研究了翼型最大彎度、最大彎度弦向位置以及厚度3個不確定性因素對升力系數(shù)的影響。Weaver和Alexeenko[10]對FIRE2飛行器進行了高超聲速數(shù)值模擬，研究了來流速度、來流密度、碰撞系數(shù)和來流溫度對氣動熱和流場的影響。Bettis和Hosder[2]對高超聲速再入流動的不確定度問題進行了分析，研究混合不確定性輸入變量對氣動熱的影響。相比于國外，國內(nèi)類似的工作還未開展，相關(guān)的研究還未見公開文獻報道。

本文以Apollo飛船返回艙為研究對象開展氣動熱數(shù)值預(yù)測，選取了來流速度、來流溫度、壁面溫度和來流密度4個不確定性輸入變量，其中來流速度變化范圍為±120 m/s(±2%)，來流溫度、壁面溫度和來流密度變化范圍為±10%，利用NIPC方法對氣動熱進行不確定度量化和敏感性分析，為飛船返回艙氣動熱防護設(shè)計提供參考依據(jù)，進一步用于高超聲速飛行器的魯棒優(yōu)化設(shè)計和可靠性分析。

1 計算方法

1.1 控制方程

采用有限體積法求解層流Navier-Stokes方程組，曲線坐標(biāo)系(τ,ξ,η,ζ)下，三維無量綱守恒形式的熱化學(xué)非平衡控制方程為

(1)

(2)

S=1/J(ωi,0,0,0,0,Sv)

(3)

1.2 離散方法

時間離散采用隱式LUSGS方法；無黏通量空間離散采用Roe格式，二階MUSCL重構(gòu)，采用minmod限制器；黏性通量離散采用二階中心格式。

1.3 物理模型

真實氣體效應(yīng)計算過程中單一組分的黏性系數(shù)、平動溫度熱傳導(dǎo)系數(shù)和振動溫度熱傳導(dǎo)系數(shù)分別由Blotter曲線擬合和Eucken關(guān)系式計算?；旌蠚怏w的黏性系數(shù)、熱傳導(dǎo)系數(shù)由Wilke半經(jīng)驗公式計算。

1.4 化學(xué)反應(yīng)模型

化學(xué)非平衡源項主要反應(yīng)了組分間化學(xué)反應(yīng)的影響。根據(jù)組分選取的不同，空氣的化學(xué)反應(yīng)模型一般分為2、5、7和11組分化學(xué)反應(yīng)模型。對于ns個組分的混合氣體，設(shè)化學(xué)反應(yīng)總數(shù)為nr，其化學(xué)反應(yīng)方程為如下通式[12]：

(4)

式中：Xi為單位體積摩爾數(shù)的化學(xué)組分或催化體，Xi=ρi/Mi，Mi為氣體組分的摩爾質(zhì)量，αr和βr分別為化學(xué)反應(yīng)式中反應(yīng)物和生成物的計量系數(shù)。

則ωi具體形式為

(5)

整理后源項表達式為

(6)

1.5 不確定度分析方法

NIPC方法主要基于混沌多項式展開(Polynomial Chaos Expansions，PCE)理論，具有堅實的數(shù)學(xué)基礎(chǔ)，是一種非常有效的基于隨機展開的不確定度分析方法。PCE方法的一個重要特性是可以將任一隨機變量分解為確定和隨機2個獨立部分。在不確定性的流體力學(xué)問題中，任意隨機變量α*(如壓力、溫度和壁面熱流等)可以表示為

(7)

式中：αj(x)為確定部分的耦合系數(shù)；Ψj(λ)為以n維隨機變量λ=(λ1,λ2,…,λn)為自變量的j階隨機基函數(shù)。

在式(7)中，Ψj(λ)的項數(shù)是無窮的，出于計算量的考慮，PCE模型通常在某階p被截斷。相應(yīng)地，p階PCE模型的項數(shù)可以表示為階數(shù)p和隨機變量的維數(shù)n的函數(shù)[13]為

(8)

根據(jù)隨機輸入?yún)?shù)的分布類型，需要選取不同的多項式基函數(shù)。當(dāng)輸入?yún)?shù)服從正態(tài)分布時選取Hermite正交多項式基函數(shù)；當(dāng)輸入?yún)?shù)服從均勻分布時選取Legendre正交多項式基函數(shù)。

目前，主要有2種形式的基于非嵌入式混沌多項式的不確定度分析方法：隨機響應(yīng)面法(Stochastic Response Surface Method，SRSM)和基于Galerkin投影法。本文采用隨機響應(yīng)面法來求解PCE系數(shù)，基本流程如圖1所示。

為了計算PCE模型中的未知PCE系數(shù)αj(j=0,1,…,P)，首先需要選取有效的樣本點，這里將回歸中用到的樣本點數(shù)表示為N。通常采用過采樣(over sampling)策略，文獻[8]中推薦用兩倍于未知PCE系數(shù)個數(shù)的樣本(即N=2Nt=2(P+1))可以得到比較滿意的結(jié)果。采用最小二次回歸來求解PCE系數(shù)，將樣本λ=[λ1,λ2,…,λj,λj+1,…,λN]T和相應(yīng)的函數(shù)響應(yīng)值G=[g(X1),g(X2),…,g(XN)]T分別代入到PCE模型右端和左端得

圖1 不確定度分析方法流程Fig.1 Flowchart of uncertainty analysis method

(9)

可以簡寫為

Aα=G

定義：

(10)

利用線性最小二次回歸，使K(αj)最小，可以求得PCE系數(shù)為

(11)

敏感性是表征隨機物理問題中每個輸入變量對輸出變量不確定度的相對貢獻大小。基于方差的Sobol分解方法可以用于敏感性分析，PCE系數(shù)確定之后，Sobol指數(shù)可以通過以下推導(dǎo)獲得，輸出變量總方差可以表示為[14]

(12)

總方差可以分解為

(13)

部分方差表示為

(14)

Sobol指數(shù)定義為

(15)

輸入變量i的Sobol指數(shù)(STi)定義為包含變量i的所有部分Sobol指數(shù)的和，即

Li={(i1,i2,…,is):?k,1≤k≤s,ik=i}

(16)

2 算例驗證與結(jié)果分析

2.1 算例驗證

本文首先選取ELECTRE標(biāo)模[15]對本文的熱化學(xué)非平衡程序進行驗證。ELECTRE球錐標(biāo)模幾何外型如圖2(a)所示。來流速度u∞=4 230 m/s (Ma=12.9)，來流密度ρ∞=6.944×10-4kg/m3，來流溫度T∞=265 K，空氣質(zhì)量分?jǐn)?shù)組成為YN=0.77，YO=0.23，YN和YO分別為N2和O2的質(zhì)量分?jǐn)?shù)，采用雙溫度熱力學(xué)模型和Park的11組分21反應(yīng)化學(xué)反應(yīng)動力學(xué)模型，壁面為等溫壁(Tw=343 K)，并且設(shè)為完全催化壁。網(wǎng)格示意圖如圖2(b)所示，網(wǎng)格數(shù)為80×90，第一層網(wǎng)格法向間距取為1×10-5m。

圖3為本文程序的熱流計算結(jié)果與飛行試驗的對比。計算結(jié)果與飛行試驗數(shù)據(jù)相吻合。在頭部駐點區(qū)域，本文程序計算的熱流與Fay-Riddell熱流公式的結(jié)果吻合，當(dāng)略高于飛行試驗結(jié)果。產(chǎn)生差異的原因是在實際飛行條件下ELECTRE標(biāo)模的熱防護設(shè)備減小了壁面催化效果[16]。

圖2 ELECTRE標(biāo)模幾何與網(wǎng)格示意圖Fig.2 Schematic of geometry and computational mesh of ELECTRE vehicle

圖3 ELECTRE標(biāo)模壁面熱流分布Fig.3 Heat flux distribution on wall surface of ELECTRE vehicle

圖4對比了本文與文獻[15]中ELECTRE標(biāo)模的駐點線溫度分布。由圖4可知，本文程序的計算結(jié)果與文獻[15]的計算結(jié)果相吻合，圖中T為平動溫度，Tv為振動溫度。

圖4 ELECTRE標(biāo)模駐點線溫度分布Fig.4 Temperature distribution along stagnation line of ELECTRE vehicle

2.2 結(jié)果分析

2.2.1 網(wǎng)格無關(guān)性驗證與基準(zhǔn)狀態(tài)分析

計算模型為鈍頭體的Apollo飛船返回艙前體[15]，Apollo飛船返回艙幾何外型如圖5(a)所示，基準(zhǔn)狀態(tài)流動條件如表1所示。計算中采用雙溫度熱力學(xué)模型，選用Park的11組分21反應(yīng)化學(xué)反應(yīng)動力學(xué)模型。壁面溫度取為800 K，并且設(shè)為非催化壁面條件。網(wǎng)格如圖5(b)所示，第一層網(wǎng)格高度為2×10-4m。

為了保證計算在一定網(wǎng)格密度下不受網(wǎng)格數(shù)量的影響，本文選取了3個水平的網(wǎng)格數(shù)量來開展網(wǎng)格無關(guān)性研究：A(100×140)，B(75×110)，C(50×80)。圖6給出了3個網(wǎng)格水平下壁面熱流結(jié)果對比?？梢钥闯鼍W(wǎng)格A與網(wǎng)格B結(jié)果非常接近，而網(wǎng)格C與其有一定區(qū)別，故網(wǎng)格B已經(jīng)可以滿足網(wǎng)格無關(guān)性要求，本文采用網(wǎng)格B(75×110)計算。

圖7為基準(zhǔn)狀態(tài)下平動溫度和壓力的等值線圖，圖8為平動溫度及振動溫度沿駐點線變化。激波過后，平動溫度首先被激發(fā)而升高。隨后平動能和振動能之間產(chǎn)生能量松弛過程，振動溫度也隨之增加。一段距離后，平動溫度和振動溫度達到平衡。

圖5 Apollo返回艙幾何與網(wǎng)格示意圖Fig.5 Schematic of geometry and computational mesh of Apollo returning capsule

參 數(shù)高度/kmu∞/(m·s-1)T∞/Kρ∞/(kg·m-3)YN/YO數(shù) 值67.360402259.8×10-50.77/0.23

圖6 不同網(wǎng)格水平下壁面熱流對比Fig.6 Wall heat flux comparison under different grid levels

圖7 基準(zhǔn)狀態(tài)平動溫度與壓力分布Fig.7 Translational temperature and pressure contours for baseline case

圖8 Apollo返回艙駐點線溫度分布Fig.8 Temperature distribution along stagnation line of Apollo returning capsule

2.2.2 不確定度量化

選取大氣來流條件(來流速度、來流溫度與來流密度)和邊界條件(壁面溫度)共4個變量作為隨機輸入來研究輸出響應(yīng)的不確定度傳播和敏感性分析。4個變量的不確定度大小如表2所示，取值范圍參考了文獻[10]，來流速度變化范圍為±120 m/s(±2%)，來流溫度、壁面溫度和來流密度變化范圍為±10%，并假設(shè)均服從均勻分布。

PCE模型在二階截斷，由式(8)計算得到樣本數(shù)為30，采用拉丁超立方(Latin Hypercube Design，LHD)[12]方法進行抽樣。95%置信度條件下，定義輸出響應(yīng)壁面熱流的不確定度為UQ%=100×1.96σR/μR。

表2 輸入變量不確定度

圖9為基準(zhǔn)狀態(tài)下的壁面熱流和NIPC方法計算的壁面均值熱流。壁面熱流沿徑向距離總體呈先減小后增大再減小的趨勢，在駐點和肩部附近出現(xiàn)峰值，在肩部附近，曲率突然減小且膨脹波氣流加速換熱效應(yīng)起主導(dǎo)作用，故其壁面熱流較高?；鶞?zhǔn)狀態(tài)下壁面熱流和NIPC方法計算的壁面均值熱流在駐點處差異為3.2%，在肩部峰值處差異約為0.78%。

圖10為相應(yīng)的壁面熱流不確定度和輸入變量敏感性指數(shù)分布。壁面熱流的不確定度沿徑向距離先減小后增大再減小，在駐點和肩部存在峰值分別約為19.8%和17.3%，相應(yīng)的熱流不確定度區(qū)間大小分別約為0.087 MW/m2和0.076 MW/m2。

駐點和肩部熱流密度的平均值、標(biāo)準(zhǔn)差、95%置信不確定區(qū)間和不確定度如表3所示。

不確定度量化結(jié)果表明，不確定性輸入變量會引起壁面熱流明顯的變化，其不確定度范圍約為15.9%～19.8%，并且在駐點(19.8%)和肩部(17.3%)達到峰值。

圖9 壁面熱流對比Fig.9 Wall heat flux comparison

圖10 壁面熱流不確定度和輸入變量Sobol指數(shù)分布Fig.10 Uncertainty of wall heat flux and Sobol index distribution of input variable

量化結(jié)果駐點SR=0肩部SR=1817mm平均值μR/(MW·m-2)0.22010.2230標(biāo)準(zhǔn)差σR/(MW·m-2)0.02230.01952-σ區(qū)間(0.1764,0.2638)(0.1848,0.2612)不確定度/%19.817.3

2.2.3 敏感性分析

各輸入變量的Sobol指數(shù)沿徑向距離變化如圖9所示，在來流速度、來流溫度、壁面溫度和來流密度4個不確定性變量中，各變量的Sobol指數(shù)在駐點處稍大,沿徑向距離變化較為平緩，故各變量對不同位置的壁面熱流影響規(guī)律較為一致。來流密度的不確定度對壁面熱流變化貢獻最大，其Sobol指數(shù)約為0.72。其次是來流速度，其Sobol指數(shù)約為0.27。相比而言壁面熱流幾乎不受來流溫度和壁面溫度影響。

駐點和肩部各輸入變量的Sobol指數(shù)如表4所示。

各輸入不確定度對輸出不確定度相對貢獻如圖11所示。在駐點處，來流密度和來流速度對壁面熱流不確定度貢獻分別為62.29%和23.56%，相應(yīng)的不確定度和不確定度區(qū)間大小分別為12.3%、4.8%和0.054、0.02 MW/m2。在肩部，來流密度和來流速度對壁面熱流不確定度貢獻分別為72.76%和26.74%，相應(yīng)的不確定度和不確定度區(qū)間大小分別為12.6%、4.7%和0.056、0.02 MW/m2。來流溫度、壁面溫度及其他交互影響項幾乎不對壁面熱流變化起作用。

敏感性分析表明，在高超聲速飛行器氣動熱數(shù)值模擬和試驗中，需要重點關(guān)注來流密度和來流速度的變化。來流密度變化(10%)和來流速度變化(2%)對壁面熱流的不確定度貢獻最大分別約為70%和25%，引起的不確定度分別約為12%和5%，故來流密度和來流速度的變化會對氣動熱數(shù)值預(yù)測產(chǎn)生不可忽略甚至明顯的影響，而來流溫度和壁面溫度幾乎不產(chǎn)生影響。

表4 駐點和肩部各輸入變量Sobol指數(shù)

圖11 駐點和肩部各輸入變量對總體不確定度的相對貢獻Fig.11 Relative importance of each input variable on overall uncertainty at stagnation and shoulder

3 結(jié) 論

本文對高超聲速地球再入飛行器進行了熱化學(xué)非平衡氣動熱數(shù)值預(yù)測，選取了來流速度、來流溫度、壁面溫度和來流密度4個不確定性輸入變量，其中來流速度變化范圍為±120 m/s(±2%)，來流溫度、壁面溫度和來流密度變化范圍為±10%，采用非嵌入式多項式混沌(NIPC)方法對氣動熱進行了不確定度量化分析和敏感性分析，得到的結(jié)論如下：

1) ELECTRE標(biāo)模熱流計算結(jié)果與飛行試驗值相吻合，表明熱化學(xué)非平衡計算程序具有較高的氣動熱預(yù)測精度。Apollo飛船返回艙再入時由于高溫高馬赫數(shù)，流場表現(xiàn)出典型的熱化學(xué)非平衡特性，壁面熱流在駐點和肩部附近存在峰值。

2) 采用基于NIPC方法的不確定性分析方法，得到了壁面各點熱流的均值、標(biāo)準(zhǔn)差和不確定度。來流條件的變化會引起壁面熱流明顯變化，壁面熱流的不確定度沿壁面先減小后增大再減小，在駐點和肩部存在2個峰值，其不確定度和不確定度區(qū)間大小分別約為19.8%、17.3%和0.087、0.076 MW/m2。

3) 采用基于Sobol指數(shù)的敏感性方法，對計算結(jié)果進行了敏感性分析。各變量對不同位置的壁面熱流影響規(guī)律一致，來流密度變化和來流速度變化對熱流不確定度貢獻最大，而來流溫度變化和壁面溫度變化幾乎不產(chǎn)生影響，在高超聲速氣動熱預(yù)測中應(yīng)重點關(guān)注來流密度和來流速度的變化。

4) 不確定度量化和敏感性分析對數(shù)值模擬可信度評估有重要意義，為高超聲速氣動熱數(shù)值預(yù)測提供了可靠的參考，可進一步用于高超聲速飛行器的魯棒性優(yōu)化設(shè)計和可靠性分析。后續(xù)工作將對熱化學(xué)非平衡數(shù)值模擬中的化學(xué)反應(yīng)模型進行不確定度分析。

參考文獻 (References)

[1] 閻超.計算流體力學(xué)方法及應(yīng)用[M].北京:北京航空航天大學(xué)出版社:2006:2-6.

YAN C.Computational fluid dynamic’s methods and applications[M].Beijing:Beihang University Press,2006:2-6(in Chinese).

[2] BETTIS B,HOSDER S.Uncertainty quantification in hypersonic reentry flows due to aleatory and epistemic uncertainties[C]∥49th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition.Reston:AIAA,2011.

[3] BRUNE A J,WEST T K,HOSDER S,et al.A review of uncertainty analysis for hypersonic inflatable aerodynamic decelerator design[C]∥21st AIAA International Space Planes and Hypersonics Technologies Conference.Reston:AIAA,2017:2373.

[4] CHAMPION K S W.Middle atmosphere density data and comparison with models[J].Advances in Space Research,1990,10(6):17-26.

[5] ECKHARDT R.Stan Ulam,John Von Neumann,and the Monte Carlo method[J].Los Alamos Science Special Issue,1987(15Special):131-137.

[6] ZAREMBA S K.The mathematical basis of Monte Carlo and quasi-Monte Carlo methods[J].SIAM Review,1968,10(3):303-314.

[7] HOSDER S,WALTERS R,PEREZ R.A non-intrusive polynomial chaos method for uncertainty propagation in CFD simulations[C]∥44th AIAA Aerospace Sciences Meeting and Exhibit.Reston:AIAA,2006：891.

[8] HOSDER S,WALTERS R,BALCH M.Efficient sampling for non-intrusive polynomial chaos applications with multiple uncertain input variables[C]∥Non-Deterministic Approaches Conference.Reston:AIAA,2007：1939.

[9] LOEVEN A,BIJL H.Airfoil analysis with uncertain geometry using the probabilistic collocation method[C]∥49th AIAA/ASME/ASCE/AHS/ASC Structures,Structural Dynamics,and Materials Conference.Reston:AIAA,2008：2070.

[10] WEAVER A B,ALEXEENKO A.Flowfield uncertainty analysis for hypersonic CFD simulations[C]∥48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition.Reston:AIAA,2010：1180.

[11] 周禹.高超聲速熱化學(xué)非平衡流場數(shù)值模擬研究[D].北京：北京航空航天大學(xué)，2009：13-19.

ZHOU Y.Numerical simulation of hypersonic thermal and chemical nonequilibrium flows[D].Beijing:Beihang University,2009:13-19(in Chinese).

[12] WANG X,YAN C,ZHENG W L,et al.Laminar and turbulent heating predictions for mars entry vehicles[J].Acta Astronautica,2016,128:217-228.

[13] 熊芬芬,楊樹興,劉宇,等.工程概率不確定性分析方法[M].北京:科學(xué)出版社,2015:115-117.

XIONG F F,YANG S X,LIU Y,et al.Engineering probability uncertainty analysis method[M].Beijing:Science Press,2015:115-117(in Chinese).

[14] SOBOL I M.Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[J].Mathematics & Computers in Simulation,2001,55(1-3):271-280.

[15] MUYLAERT J,WALPOT L,HAUSER J.Standard model testing in the European high enthalpy facility F4 and extrapolation to flight[C]∥Joint Propulsion Conference and Exhibit.Reston:AIAA,1992:S396.

[16] HAO J A,WANG J Y,LI C H.Numerical study of hypersonic flows over reentry configurations with different chemical nonequilibrium models[J].Acta Astronautica,2016,126:1-10.