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

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

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

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

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

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

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

1 錐導(dǎo)乘波體的設(shè)計

1.1 錐導(dǎo)乘波體的生成原理

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

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

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

(1)

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

1.2 乘波體的外形參數(shù)化

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

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

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

表1 乘波體設(shè)計參數(shù)Table 1 Design parameters of waverider

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

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

2 混沌多項式理論

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

(2)

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

(3)

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

(4)

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

(5)

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

對于滿足正態(tài)分布的隨機變量，基函數(shù)是Hermite正交多項式,則有

(6)

進而可以推導(dǎo)出

(7)

考慮n維的Hermite多項式：

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

(8)

基函數(shù)：

(9)

系統(tǒng)總方差可以表示為[11]

(10)

表2 標(biāo)準(zhǔn)分布類型與正交多項式的對應(yīng)關(guān)系

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

1≤i1<…

(11)

然后，靈敏度Sobol指數(shù)可以定義為

(12)

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

3 流場求解及分析

3.1 求解方法及網(wǎng)格無關(guān)性驗證

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

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

表3 網(wǎng)格無關(guān)性驗證結(jié)果Table 3 Results of grid-independent validation

圖4 計算網(wǎng)格Fig.4 Calculation grid

3.2 計算結(jié)果分析

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

表4 樣本設(shè)計及結(jié)果Table 4 Sample design and results

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

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

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

Table5Uncertaintyandsensitivityanalysisresultsofaerodynamiccoefficientforflowvelocityandangleofattack

項目CLCDCmK平均值0.513700.100830.834695.0381標(biāo)準(zhǔn)差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 流場馬赫數(shù)云圖的均值和標(biāo)準(zhǔn)差分布Fig.6 Mean value distribution and standard deviation distribution of Ma contour in flow field

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

圖7 流場壓強云圖的均值和標(biāo)準(zhǔn)差分布Fig.7 Mean value distribution and standard deviation distribution of pressure contour in flow field

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

圖8 不同截面壓強標(biāo)準(zhǔn)差分布Fig.8 Distribution of pressure standard deviation in different cross-sections

3.3 靈敏度分析結(jié)果驗證

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

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

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

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

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

4 結(jié) 論

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

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

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

以上對乘波體氣動特性來流不確定度的研究，能夠提供更為可靠的數(shù)值預(yù)測，提高優(yōu)化設(shè)計的魯棒適應(yīng)性，從而獲得寬速域內(nèi)氣動性能更為穩(wěn)定的乘波體構(gòu)型，這也是本文下一步的研究目標(biāo)。

[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] 徐崇剛，胡遠滿，常禹, 等. 生態(tài)模型的靈敏度分析[J]. 應(yīng)用生態(tài)學(xué)報, 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] 劉全, 王瑞利, 林忠. 非嵌入式多項式混沌方法在拉氏計算中的應(yīng)用[J]. 固體力學(xué)學(xué)報, 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].力學(xué)學(xué)報, 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] 趙志, 宋文艷, 肖隱利. 高超聲速錐導(dǎo)乘波體非設(shè)計點性能研究 [J].飛行力學(xué), 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] 王發(fā)民，丁海河，雷麥芳． 乘波布局飛行器寬速域氣動特性與研究[J].中國科學(xué)： 技術(shù)科學(xué)，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]. 固體火箭技術(shù), 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] 段焰輝，范召林，吳文華. 定后掠角密切錐乘波體的生成和設(shè)計方法[J]. 航空學(xué)報, 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] 苗萌. 高超聲速飛行器氣動布局優(yōu)化設(shè)計[D]. 北京：北京航空航天大學(xué)，2012: 43-45.

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

[20] 耿永兵, 劉宏, 姚文秀, 等.錐形流乘波體優(yōu)化設(shè)計研究[J]. 航空學(xué)報, 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]. 北京: 科學(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).

[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] 甘文彪, 閻超. 乘波飛行器一體化構(gòu)型設(shè)計[J]. 空氣動力學(xué)學(xué)報, 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).