來流導致的高超聲速氣動熱不確定度量化分析

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

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

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

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

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

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

1 計算方法

1.1 控制方程

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

(1)

(2)

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

(3)

1.2 離散方法

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

1.3 物理模型

真實氣體效應計算過程中單一組分的黏性系數、平動溫度熱傳導系數和振動溫度熱傳導系數分別由Blotter曲線擬合和Eucken關系式計算。混合氣體的黏性系數、熱傳導系數由Wilke半經驗公式計算。

1.4 化學反應模型

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

(4)

式中：Xi為單位體積摩爾數的化學組分或催化體，Xi=ρi/Mi，Mi為氣體組分的摩爾質量，αr和βr分別為化學反應式中反應物和生成物的計量系數。

則ωi具體形式為

(5)

整理后源項表達式為

(6)

1.5 不確定度分析方法

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

(7)

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

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

(8)

根據隨機輸入參數的分布類型，需要選取不同的多項式基函數。當輸入參數服從正態分布時選取Hermite正交多項式基函數；當輸入參數服從均勻分布時選取Legendre正交多項式基函數。

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

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

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

(9)

可以簡寫為

Aα=G

定義：

(10)

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

(11)

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

(12)

總方差可以分解為

(13)

部分方差表示為

(14)

Sobol指數定義為

(15)

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

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

(16)

2 算例驗證與結果分析

2.1 算例驗證

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

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

圖2 ELECTRE標模幾何與網格示意圖Fig.2 Schematic of geometry and computational mesh of ELECTRE vehicle

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

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

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

2.2 結果分析

2.2.1 網格無關性驗證與基準狀態分析

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

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

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

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

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

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

圖7 基準狀態平動溫度與壓力分布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個變量作為隨機輸入來研究輸出響應的不確定度傳播和敏感性分析。4個變量的不確定度大小如表2所示，取值范圍參考了文獻[10]，來流速度變化范圍為±120 m/s(±2%)，來流溫度、壁面溫度和來流密度變化范圍為±10%，并假設均服從均勻分布。

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

表2 輸入變量不確定度

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

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

駐點和肩部熱流密度的平均值、標準差、95%置信不確定區間和不確定度如表3所示。

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

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

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

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

2.2.3 敏感性分析

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

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

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

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

表4 駐點和肩部各輸入變量Sobol指數

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

3 結 論

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

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

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

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

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

參考文獻 (References)

[1] 閻超.計算流體力學方法及應用[M].北京:北京航空航天大學出版社: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] 周禹.高超聲速熱化學非平衡流場數值模擬研究[D].北京：北京航空航天大學，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].北京:科學出版社,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.