錐導乘波體構型的氣動特性不確定度分析

宋賦強，閻超，馬寶峰，鞠勝軍

北京航空航天大學 航空科學與工程學院，北京 100083

先進氣動布局是高超聲速飛行器發展的關鍵技術之一，在傳統的高超聲速氣動布局優化設計中，只考慮提高確定狀態的氣動性能[1-2]，沒有研究不確定性因素帶來的影響。而在實際的高超聲速飛行中，不可避免會遇到一些不確定性因素的干擾，如來流參數變化、氣彈載荷變化等，這些不確定性因素對飛行器的氣動特性會產生較大影響，從而降低優化結果的可靠性。近年來，將穩健性優化設計(又稱魯棒優化設計)的概念應用于飛行器的氣動設計中，逐漸引起研究人員的重視[3-5]。對于穩健性優化設計而言，其前提條件是不確定性量化分析,在NASA出版的白皮書中明確指出要在飛行器優化設計中考慮不確定性因素的影響[6]。

然而，單純的不確定性分析僅能獲得不確定性因素對輸出量的影響，在實際問題中還應該關注哪些是影響顯著的變量，這就需要對不確定性因素進行靈敏度分析。靈敏度分析主要分成兩類，包括局部靈敏度和全局靈敏度分析，局部靈敏度分析只檢驗某個選定參數的變化對輸出量的影響程度，而全局靈敏度分析不僅可以考察單個參數對輸出量的影響，還能分析每一個參數以及參數之間的耦合作用[7]。近年來，全局靈敏度分析方法也得到了廣泛的研究和應用，不但可以深入研究每個隨機輸入變量對輸出的貢獻程度，還容易求出不同變量之間的交互影響[8-9]。

傳統方法在進行飛行器氣動特性的不確定度分析時，存在計算量巨大、計算效率低、過程復雜等問題，而采用非嵌入式混沌多項式(NIPC)不確定度分析方法結合CFD工具則是一種有效的研究手段[10]。近年來，非嵌入式混沌多項式方法在流體力學領域的不確定性分析中得到廣泛的應用。Westiv等[11]采用非嵌入式的多項式混沌方法，對返回艙后體進行不確定度及全局靈敏度分析，降低隨機問題維度的同時提高了預測的精度。Weaver等[12]研究了高超聲速數值模擬的流場不確定度分析方法，得到了飛行器表面熱流的不確定度。劉全等[13]研究了混沌多項式法的數學理論基礎，給出了采用混沌多項式方法進行不確定度量化的主要流程。鄔曉敬等[14]對翼型跨聲速流場進行了不確定度分析，找出了氣動特性波動的主要因素，并對非嵌入式混沌多項式方法和蒙特卡羅方法進行比較，發現非嵌入式混沌多項式方法計算效率更高。

乘波體是一種典型的高超聲速氣動構型，其最顯著的特征是在設計狀態激波完全附著于整個前緣，從而產生很高的升阻比。自乘波體的概念提出以來，國內外都已開展了比較成熟的研究，廣泛應用于吸氣式高超聲速巡航飛行器、再入飛行器和吸氣式高超聲速導彈中[15]。然而，乘波體是在特定的條件下設計的，在偏離設計條件下，氣動性能可能會出現一定程度的下降。所以，對乘波體進行氣動特性的不確定度分析，再通過穩健性設計方法尋求寬速域內性能較好的乘波體構型具有重要的實際意義。目前，以寬速域乘波體飛行器為目標的設計研究才剛剛起步，國內有王發民[16]、李世斌[17]、段焰輝[18]等進行了相關研究。而且，針對乘波體這種構型的氣動特性的來流影響分析，對單個因素如迎角、馬赫數等的研究比較多，而多個因素的交叉影響則鮮有涉及，因此在多變量綜合影響下對乘波體進行氣動特性不確定性分析是很有必要的。

本文首先進行錐導乘波體的參數化建模，然后在可變的來流條件下，采用拉丁超立方試驗設計構建樣本，通過計算流體力學方法對每個樣本點進行計算。最后，基于混合的混沌多項式方法，對服從不同分布的多個流動變量進行乘波體氣動特性的不確定度量化和全局靈敏度分析，研究了不同變量之間的耦合情況，得到影響最顯著的變量，并具體研究來流不確定度對流場的影響，獲得多個變量綜合作用下的流場變化情況，分析了乘波體氣動特性變化的流動機理。

1 錐導乘波體的設計

1.1 錐導乘波體的生成原理

圖1 錐導乘波體示意圖[19]Fig.1 Schematic diagram of cone-derived waverider[19]

錐導乘波體是根據圓錐繞流生成的乘波體，其原理見圖1[19]。首先選定設計馬赫數以及一個基準圓錐，然后求解該圓錐繞流的解。其次，選取一個柱面作為流動捕獲管(FCT)，與圓錐激波面相交形成乘波體前緣。然后在前緣曲線上選定若干個點，求出每個點在過激波后流場的流線，組成流面，從而構成乘波體的下表面，而上表面的確定則是通過自由流面來形成的。這樣組合起來，就構成了乘波體的外形。因而，不同的圓錐角、流動捕獲管、流動捕獲管與圓錐激波面交線的軸向位置，能夠獲得不同的乘波體外形。

圓錐流場的流動物理量可以通過Runge-Kutta方法積分求解Taylor-Maccoll方程得到[20]，即

(1)

式中：vr為圓錐射線的速度分量；φ為圓錐射線與中心軸的夾角；γ為比熱比。

1.2 乘波體的外形參數化

根據錐導乘波體的生成原理，確定來流馬赫數Ma，半圓錐角δ、無黏流場的長度L以及乘波體上表面的底部型線后，即可以獲得唯一的乘波體外形。這里選取適應性比較好的四次曲線y=az4+bz2+c作為上表面的底部型線，乘波體外形的參數化見圖2，β為激波角，參數k=c/y0，ψ為OA軸與z軸夾角，θ為底部型線在A點處的切線和水平方向的夾角，Lw為乘波體長度，由底部型線向前延伸得到。因此，決定乘波體外形生成的參數一共有6個，分別是Ma、L、δ、ψ、θ和k。

在控制外形的各個參數確定后，可以通過Visual Basic編程進行CATIA二次開發來生成乘波體外形。先求解錐型流場得到前緣，然后進行流線追蹤得到下底面流線，再通過曲面擬合來生成乘波體的下表面，上表面可以直接由前緣沿自由來流拉伸到底部位置產生。生成錐導乘波體的具體參數見表1，外形三視圖見圖3。

圖2 乘波體外形參數化示意圖Fig.2 Parametric diagram of waverider configuration

表1 乘波體設計參數Table 1 Design parameters of waverider

參數MaL/mδ/(°)ψ/(°)θ/(°)數值61025510

圖3 乘波體外形三視圖Fig.3 Three views of waverider configuration

2 混沌多項式理論

混沌多項式的基本思想是用包含獨立隨機變量的正交多項式之和來模擬隨機過程，將輸入函數分解成隨機和確定兩部分，這是一個無窮項，需要進行截斷，得到[21]

(2)

式中：α*為確定獨立變量X和n維隨機變量ξ=[ξ1ξ2…ξn]的函數，ξ服從設定的概率分布；αi(X)和Ψi(ξ)分別為第i階模態的確定和隨機部分；P為所截斷的有限階模態階數，即

(3)

式中：n為隨機變量的維數；s為混沌多項式的階數?；煦缍囗検椒椒ǖ年P鍵是求出多項式系數，對于j項，式(2)做如下變換：

(4)

式中：〈·〉表示內積，函數f(ξ)和g(ξ)在定義域上的內積可表示為

(5)

式中：w(ξ)為基函數。不同的概率分布類型對應的正交多項式基函數是不一樣的。表2給出了2種分布類型和正交多項式之間的對應關系。

對于滿足正態分布的隨機變量，基函數是Hermite正交多項式,則有

(6)

進而可以推導出

(7)

考慮n維的Hermite多項式：

H(ξ1,ξ2,…,ξn)=

(8)

基函數：

(9)

系統總方差可以表示為[11]

(10)

表2 標準分布類型與正交多項式的對應關系

而根據混沌多項式理論，局部方差為

1≤i1<…

(11)

然后，靈敏度Sobol指數可以定義為

(12)

求解方程多項式系數有嵌入式與非嵌入式兩種方法。嵌入式方法是將響應量映射到基函數Ψ(ξ)上，需要更改計算程序，處理過程復雜且計算量龐大，非常不方便。而非嵌入式方法把系統作為黑箱處理，通過回歸分析來求解系數α(X)，從而建立近似的數學模型，不需要對控制方程進行修改，這樣能夠充分利用已經編好的數值模擬程序，大幅減小不確定性分析所需要的工作量[22]。

3 流場求解及分析

3.1 求解方法及網格無關性驗證

本文采用文獻[23]提出的MI-CFD軟件平臺對流場進行求解?？臻g格式為Roe格式；時間格式選擇LU-SGS格式推進；湍流模型采用Menter-SST模型，使用minmod限制器。參考長度為5.085 m,參考面積為1.109 m2，力矩取矩點坐標為坐標原點，壁面條件為絕熱壁。

為確保計算的可靠性，采用粗、中、密3套網格進行無關性驗證，網格量分別是52萬、112萬和309萬，第1層網格高度的y+及計算結果見表3。其中，計算狀態為設計狀態，即馬赫數Ma=6.0，迎角α=0°，飛行高度h=26 km?？梢缘贸觯ο禂礐L、 阻力系數CD和俯仰力矩系數Cm的相對誤差ε都在3%以內，而且中網格和細網格的結果更為接近(1%以內)。因此，綜合考慮計算精度和計算效率等問題，采用中網格進行計算，計算網格見圖4。

表3 網格無關性驗證結果Table 3 Results of grid-independent validation

圖4 計算網格Fig.4 Calculation grid

3.2 計算結果分析

本節主要研究4個獨立變量(來流速度v、迎角α、來流密度ρ和來流溫度T)的隨機擾動對乘波體氣動特性的影響，并進行全局非線性靈敏度分析，找出對氣動特性影響較大的因素。假設來流速度v服從N(1 795, 1502)的正態分布，迎角α服從N(0, 22)的正態分布,來流密度ρ服從U(0.018, 0.05)的均勻分布, 來流溫度T服從U(217,228)的均勻分布，對應的飛行高度約為20～30 km。采用拉丁超立方試驗設計方法產生了30個樣本，樣本設計表及計算結果見表4，響應量分別為升力系數CL、阻力系數CD、俯仰力矩系數Cm和升阻比K。然后采用混沌多項式方法,進行線性最小二次回歸分析，得到4個參數對響應量的Sobol指數,計算見式(12)。Sobol指數的大小表征了該參數對于乘波體氣動特性不確定度的靈敏程度，其歸一化結果見圖 5。可以直觀看出，在4個參數中，迎角對于乘波體的氣動特性的靈敏程度是最高的，占主導作用,其次是來流速度，來流密度和來流溫度對升力系數和俯仰力矩系數的影響較小，但這兩項對于升阻比的影響有所增大。此外，升力系數和俯仰力矩系數的Sobol指數大小很接近，這跟取矩點為原點關系密切。

表4 樣本設計及結果Table 4 Sample design and results

由于迎角和來流速度對乘波體的氣動特性不確定度影響較大，下面重點分析這兩個參數的影響。表5 給出了氣動力系數的不確定性及靈敏度統計分析結果，S是Sobol指數的簡寫， 如Sv -α表示來流速度和迎角耦合項對于各響應量的Sobol指數。從表中可以看出，升力系數、阻力系數和俯仰力矩系數的不確定度均較大，在21.5%以上，而升阻比的不確定度較小，只有8.06%。由此可見，雖然升阻力系數波動較大，但乘波體的高升阻比特性仍保持得比較良好。靈敏度分析結果表明：① 在統計的分量中，迎角一階分量的靈敏度指標是最高的，而且對于阻力系數和升阻比，迎角二階分量的靈敏度指標也比速度一階分量的高，這說明了迎角對乘波體氣動特性變化貢獻是最大的；② 速度一階分量對于響應量的貢獻遠大于來流速度二階分量，整體而言，來流速度對響應量的波動貢獻不算大；③ 速度和迎角的耦合項的靈敏度指標最大只有0.036 24，說明來流速度和迎角的耦合作用對氣動特性波動的貢獻很小。

圖5 來流參數的Sobol指數歸一化比較Fig.5 Comparison of Sobol indices normalization of flow parameters

表5 來流速度與迎角對氣動力系數的不確定度及靈敏度結果

Table5Uncertaintyandsensitivityanalysisresultsofaerodynamiccoefficientforflowvelocityandangleofattack

項目CLCDCmK平均值0.513700.100830.834695.0381標準差0.1360.02190.2240.406不確定度/%26.521.726.98.06Sα0.837110.653990.836860.35423Sv0.087590.085780.084520.02625S2α0.023290.189200.027130.43070S2v0.006930.004790.006250.00948Sv-α0.036240.018200.034960.003091

圖6 流場馬赫數云圖的均值和標準差分布Fig.6 Mean value distribution and standard deviation distribution of Ma contour in flow field

目前，大多數氣動特性不確定度的研究只關注不確定性因素對氣動力系數的影響，這不利于從流場細節分析不確定性參數對整個流場分布的作用機理。因此，本文分析了來流參數不確定度對整個流場馬赫數和壓強的分布影響。為了便于直觀理解，將壓強p無量綱化為p′=p/p∞，p∞為無窮遠處的來流壓強。圖6和圖7分別是流場的馬赫數以及壓強的均值和標準差分布圖，其中，標準差數值大的地方意味著流場變化劇烈。從這兩個圖可以看出，乘波體上表面被自由來流包裹著，除了一道較弱的頭部激波外，流場變化相對平緩。在乘波體的下表面，激波依附在前緣上，構成了典型的“乘波”效應。對于整個流場，來流條件的不確定性主要影響的是乘波體下表面的附著激波，具體表現為削弱前緣對下表面高壓氣體往上表面泄露的抵抗能力，導致乘波體氣動性能的下降。從圖 7的壓力云圖可以發現，在乘波體的底部，壓力泄露已經造成上表面局部壓強增高，對乘波體壁面壓強分布產生了明顯的污染。

圖7 流場壓強云圖的均值和標準差分布Fig.7 Mean value distribution and standard deviation distribution of pressure contour in flow field

圖8給出了乘波體不同截面的壓強標準差分布，結果對比表明：不同截面的下表面壓強標準差近似相等，其中曲線峰值對應的是前緣線處的壓強波動，可以看出，隨著x增大，前緣線的壓強波動程度加劇。在底面x=10.00 m處，上表面靠近前緣處的壓強已經超過了下表面，這正是下表面高壓氣體泄露所造成的影響。

圖8 不同截面壓強標準差分布Fig.8 Distribution of pressure standard deviation in different cross-sections

3.3 靈敏度分析結果驗證

根據靈敏度分析結果可以知道，迎角對乘波體氣動特性的不確定度影響最大，來流密度和溫度則影響很小。為進一步驗證流場分析的結果，在設計馬赫數Ma=6狀態下，分別研究了不同迎角和來流密度下乘波體氣動特性的變化規律。由于升力系數和俯仰力矩系數的靈敏度結果接近，為簡化分析，這里只給出升力系數、阻力系數和升阻比的結果，不再分析俯仰力矩系數的變化趨勢。

圖9為不同迎角下乘波體的氣動特性曲線，計算迎角α=-10°～10°，間隔2°，其他來流參數不變?？梢钥吹?，乘波體的升力系數隨著迎角呈近似的線性變化，阻力系數則先減小后增大。在負迎角的情況下，升阻比先減后增，正迎角的時候則是先增后減，而且，負迎角時的升阻比變化范圍比較大，可見迎角對乘波體氣動特性影響顯著。在迎角α=0°～2°之間，達到最大的升阻比，且在正迎角范圍內，升阻比均大于3。總體而言，升力系數、阻力系數和升阻比的變化幅度都很大，驗證了迎角對于乘波體的氣動特性變化非常靈敏。

圖9 不同迎角下乘波體的氣動特性曲線Fig.9 Aerodynamic characteristic curves of waverider at different angles of attack

圖10為不同來流密度下乘波體的氣動特性曲線，計算來流密度ρ=0.015～0.06 kg/m3，間隔0.005 kg/m3，對應的飛行高度約為18～35 km,其他來流參數不變?？梢钥闯觯谟嬎銇砹髅芏确秶鷥?，升力系數和阻力系數的變化很小，升阻比隨著密度增加而增加，但變化幅度也只有5.6%。由此可以驗證，來流密度對于乘波體氣動特性的敏感度很小。

圖10 不同來流密度下乘波體的氣動特性曲線Fig.10 Aerodynamic characteristic curves of waverider at different flow densities

4 結 論

1) 對乘波體偏離設計狀態的氣動特性進行不確定度分析,在迎角和來流參數(來流速度、來流密度和來流溫度)可變的情況下，相比于氣動力系數而言，升阻比的不確定度較小，乘波體高升阻比特性受來流參數影響相對不大。

2) 靈敏度分析結果表明，在流動參數的變化范圍內，迎角對乘波體氣動特性的不確定度的影響最大，起主導作用，其次是來流速度，來流密度和來流溫度影響較小，但迎角和來流速度之間的交叉耦合作用不明顯，在一定程度上可以簡化乘波體來流不確定度的分析。

3) 通過對流場不確定性的分析，發現來流狀態的不確定性加劇了前緣處的壓力脈動，雖然下表面高壓維持較好，但是靠近底面的高壓氣體的泄露污染了上表面壓力分布，造成乘波體氣動性能的改變。

以上對乘波體氣動特性來流不確定度的研究，能夠提供更為可靠的數值預測，提高優化設計的魯棒適應性，從而獲得寬速域內氣動性能更為穩定的乘波體構型，這也是本文下一步的研究目標。

[1] AHN J, KIM H J, LEE D H, et al. Response surface method for airfoil design in transonic flow[J]. Journal of Aircraft, 2001, 38(2): 231-238.

[2] VAVALLE A, QIN N. Iterative response surface based optimization scheme for transonic airfoil design[J]. Journal of Aircraft, 2007, 44(2): 365-376.

[3] PADULO M, MAGINOT J, GUENOV M , et al. Airfoil design under uncertainty with robust geometric parameterization: AIAA-2009-2270[R]. Reston, VA: AIAA, 2009.

[4] PADULO M, CAMPOBASSO M S, GUENOV M D. Novel uncertainty propagation method for robust aerodynamic design[J]. AIAA Journal, 2011, 49(3): 530-543.

[5] CROICU A M, HUSSAINI M Y, JAMESON A, et al. Robust airfoil optimization using maximum expected value and expected maximum value approaches[J]. AIAA Journal, 2012, 50(9): 1905-1919.

[6] ZANG T A, HEMSCH M J, HILBURGER M W, et al. Needs and opportunities for uncertainty-based multidisciplinary design methods for aerospace vehicles: NASA/TM-2002-211462[R]. Washington, D.C.: NASA, 2002.

[7] 徐崇剛，胡遠滿，常禹, 等. 生態模型的靈敏度分析[J]. 應用生態學報, 2004,15(6): 1056-1062.

XU C G, HU Y M, CHANG Y, et al. Sensitivity analysis in ecological modeling[J]. Chinese Journal of Applied Ecology, 2004, 15(6): 1056-1062 (in Chinese).

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

[9] SALTELLI A, RATTO M, ANDRES T, et al. Global sensitivity analysis: The primer[M]. New York: John Wiley & Sons, 2008: 138-153.

[10] ZHANG Y. Efficient uncertainty quantification in aerospace analysis and design [D]. Missouri: Missouri University of Science and Technology, 2013: 13-40.

[11] WESTIV T K, JOHNSTON C O, HOSDER S. Uncertainty and sensitivity analysis of afterbody radiative heating predictions for earth entry[J]. Journal of Thermophysics and Heat Transfer, 2017, 31(2): 294-306.

[12] ALEXEENKO A, WEAVER A B, GREENDYKE R B, et al. Flowfield uncertainty analysis for hypersonic CFD simulations: AIAA-2010-1180 [R]. Reston, VA: AIAA, 2010.

[13] 劉全, 王瑞利, 林忠. 非嵌入式多項式混沌方法在拉氏計算中的應用[J]. 固體力學學報, 2013, 10(33): 224-233.

LIU Q, WANG R L, LIN Z. Uncertainty quantification for Lagrangian cmputation using non-instrusive polynomial chaos [J]. Chinese Journal of Solid Mechanics, 2013, 10(33): 224-233 (in Chinese).

[14] 鄔曉敬, 張偉偉, 宋述芳, 等. 翼型跨聲速氣動特性的不確定性及全局靈敏度分析 [J].力學學報, 2015,47(4): 588-595.

WU X J, ZHANG W W, SONG S F, et al. Uncertainty quantification and global sensitivity analysis of transonic aerodynamics about airfoil [J]. Chinese Journal of Theoretical and Applied Mechanics, 2015, 47(4): 588-595 (in Chinese).

[15] 趙志, 宋文艷, 肖隱利. 高超聲速錐導乘波體非設計點性能研究 [J].飛行力學, 2009, 27(1): 47-50.

ZHAO Z, SONG W Y, XIAO Y L. Numerical simulation on off-design performance of hypersonic cone-derived waverider[J]. Flight Dynamics, 2009, 27(1): 47-50 (in Chinese).

[16] 王發民，丁海河，雷麥芳． 乘波布局飛行器寬速域氣動特性與研究[J].中國科學： 技術科學，2009，39 (11) : 1828-1835

WANG F M, DING H H, LEI M F. Wide-velocity aerodynamic characteristics and study of waverider configuration aircrafts[J]. Science China: Technical Science, 2009, 39(11): 1828-1835 (in Chinese).

[17] 李世斌，羅世彬，黃偉, 等. 新型寬速域高超聲速飛行器氣動特性研究[J]. 固體火箭技術, 2012, 35(5): 588-592.

LI S B，LUO S B，HUANG W，et al. Investigation on aerodynamic performance for a novel wide-ranged hypersonic vehicle[J]. Journal of Solid Rocket Technology, 2012, 35(5): 588-592 (in Chinese).

[18] 段焰輝，范召林，吳文華. 定后掠角密切錐乘波體的生成和設計方法[J]. 航空學報, 2016, 37(10): 3023-3034.

DUAN Y H, FAN Z L, WU W H. Generation and design methods of osculating cone waverider with constant angle of sweepback[J] Acta Aeronautica et Astronautica Sinica, 2016, 37(10): 3023-3034 (in Chinese).

[19] 苗萌. 高超聲速飛行器氣動布局優化設計[D]. 北京：北京航空航天大學，2012: 43-45.

MIAO M. Optimization design of aerodynamic configuration of hypersonic aircraft[D]. Beijing: Beihang University, 2012: 43-45 (in Chinese).

[20] 耿永兵, 劉宏, 姚文秀, 等.錐形流乘波體優化設計研究[J]. 航空學報, 2006, 27(1): 23-28.

GENG Y B，LIU H，YAO W X, et al. Viscous optimized design of waverider derived from cone flow[J]. Acta Aeronautica et Astronautica Sinica, 2006, 27(1): 23-28 (in Chinese).

[21] 熊芬芬, 楊樹興, 劉宇, 等. 工程概率不確定性分析方法[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).

[22] XIU D, KARNIADAKIS G E. Modeling uncertainty in flow simulations via generalized polynomial chaos[J]. Journal of Computational Physics, 2003, 187(1): 137-167.

[23] 甘文彪, 閻超. 乘波飛行器一體化構型設計[J]. 空氣動力學學報, 2012, 30(1): 68-73.

GAN W B, YAN C. Waverider design of integrated configuration[J]. Acta Aerodynamica Sinica, 2012, 30(1): 68-73 (in Chinese).