季廷煒，查旭，謝芳芳，吳雨思，張鑫帥，蔣逸陽，杜昌平，鄭耀

(浙江大學(xué) 航空航天學(xué)院，浙江 杭州 310027)

氣動(dòng)性能評估是空天飛行器[1-6]設(shè)計(jì)階段的重要環(huán)節(jié)，其中氣動(dòng)力數(shù)據(jù)的獲取是空天飛行器氣動(dòng)布局、結(jié)構(gòu)設(shè)計(jì)、軌跡仿真等相關(guān)工作的基礎(chǔ).傳統(tǒng)上，空天飛行器氣動(dòng)力數(shù)據(jù)可以通過地面風(fēng)洞實(shí)驗(yàn)、飛行實(shí)驗(yàn)、計(jì)算流體動(dòng)力學(xué)（computational fluid dynamics，CFD）數(shù)值模擬、工程估算等多種途徑獲得.其中，工程估算計(jì)算結(jié)果置信度不高，依賴于風(fēng)洞試驗(yàn)結(jié)果修正，對不同外形的適應(yīng)性差，一般只用于簡單外形或飛行器概念設(shè)計(jì)的氣動(dòng)估算.地面試驗(yàn)存在靜溫低、雷諾數(shù)低、有洞壁干擾、支架干擾等問題.飛行試驗(yàn)產(chǎn)生的數(shù)據(jù)是客觀的，但數(shù)據(jù)量通常偏少，且來流條件存在擾動(dòng)以及傳感器精度問題也會(huì)導(dǎo)致辨識的數(shù)據(jù)存在偏差.雖然準(zhǔn)確可靠的氣動(dòng)力數(shù)據(jù)可以通過建造更先進(jìn)的風(fēng)洞、采用更高性能的數(shù)值模擬系統(tǒng)以及進(jìn)行更多更精細(xì)的飛行試驗(yàn)來獲取，但在短時(shí)間內(nèi)采用單一方式提高氣動(dòng)力數(shù)據(jù)精準(zhǔn)度的程度有限，不能完全滿足未來飛行器研制的需要，而且將付出高昂的成本[7].研究者希望可以將這些復(fù)雜來源的氣動(dòng)力數(shù)據(jù)進(jìn)行融合、建模，降低氣動(dòng)數(shù)據(jù)獲取成本，最大限度的提升氣動(dòng)模型建模精度.

對于氣動(dòng)力多源數(shù)據(jù)融合方法的研究，目前主要圍繞以下2 類問題展開：1）氣動(dòng)力數(shù)據(jù)源之間存在明顯精度差異的數(shù)據(jù)融合問題（有相對真解的數(shù)據(jù)融合問題）.受到計(jì)算成本、變量多樣、魯棒性等條件限制，高精度的氣動(dòng)數(shù)據(jù)（比如飛行試驗(yàn)）難以充分獲得，而低精度氣動(dòng)數(shù)據(jù)的獲取較便利.考慮到低精度氣動(dòng)數(shù)據(jù)與高精度氣動(dòng)數(shù)據(jù)的匹配特性，低精度和高精度模型具有相近的或者有關(guān)聯(lián)的內(nèi)在特征，因此可以通過混合高精度與低精度數(shù)據(jù)（或模型）建立數(shù)據(jù)融合模型以逼近更高的數(shù)據(jù)精度[8-13].2）氣動(dòng)數(shù)據(jù)精度區(qū)分不明確、狀態(tài)不匹配的數(shù)據(jù)融合問題.在采用數(shù)值計(jì)算氣動(dòng)數(shù)據(jù)、風(fēng)洞試驗(yàn)數(shù)據(jù)與飛行試驗(yàn)數(shù)據(jù)對飛行器氣動(dòng)力或模型進(jìn)行確定時(shí)，須考慮不同狀態(tài)、不同數(shù)據(jù)來源的影響.不同氣動(dòng)數(shù)據(jù)間的狀態(tài)與參數(shù)很難達(dá)到統(tǒng)一，這時(shí)就須針對不同狀態(tài)、不同來源的數(shù)據(jù)進(jìn)行融合，提升多源氣動(dòng)力數(shù)據(jù)庫的一致性和精度.近年來，機(jī)器學(xué)習(xí)在多源數(shù)據(jù)融合工作中逐步得到了嘗試和應(yīng)用[14-18].它不僅能通過融合不同精度的數(shù)據(jù)實(shí)現(xiàn)高精確度預(yù)測，而且可以對模型的不確定性進(jìn)行量化估計(jì).因此，將傳統(tǒng)方法與機(jī)器學(xué)習(xí)技術(shù)緊密結(jié)合，為空天飛行器氣動(dòng)問題求解，甚至流體力學(xué)的發(fā)展都提供了嶄新的思路.

傳統(tǒng)的工程估算方法計(jì)算速度快，但是無法考慮飛行中的諸多特殊復(fù)雜的物理現(xiàn)象，預(yù)測精度低；CFD 數(shù)值模擬可以高分辨率解析飛行中的物理現(xiàn)象，但是計(jì)算成本大、周期長.本研究基于工程估算和CFD 數(shù)值模擬方法，引入高斯過程回歸理論，提出面向空天飛行器的多精度氣動(dòng)建模方法，以求在較低計(jì)算成本下實(shí)現(xiàn)高精度氣動(dòng)特性預(yù)測.隨后，將該方法應(yīng)用于FTB 空天飛行器的氣動(dòng)建模問題中，構(gòu)建了該飛行器的多精度氣動(dòng)模型，并將氣動(dòng)模型與再入軌跡仿真模型相耦合.通過對比分析仿真結(jié)果，凸顯多精度氣動(dòng)建模方法的工程價(jià)值.

1 空天飛行器氣動(dòng)分析方法

1.1 CFD 數(shù)值模擬

基于物理學(xué)三大守恒定律可以對流體運(yùn)動(dòng)進(jìn)行嚴(yán)格的數(shù)學(xué)描述，從而推及由質(zhì)量守恒方程、動(dòng)量守恒方程、能量守恒方程構(gòu)成的流體力學(xué)基本方程形式即N-S 方程組.由于空天飛行器的高速運(yùn)動(dòng)特性，忽略其中黏性項(xiàng)，進(jìn)而得到三維非定?？蓧嚎s流體的Euler 方程組，其微分形式如下：

式中：W為守恒變量，F(xiàn)為二階對流通量張量.

采用由美國斯坦福大學(xué)開發(fā)的開源CFD 求解器SU2，它是一款基于C++編程語言和非結(jié)構(gòu)網(wǎng)格的高精度偏微分方程求解器，能夠解決非結(jié)構(gòu)網(wǎng)格拓?fù)湎碌亩辔锢韴龇治龊蛢?yōu)化問題.SU2 求解器的流動(dòng)分析能力在RAM-C II 高超音速飛行器的試驗(yàn)案例中[19]已經(jīng)過嚴(yán)格驗(yàn)證和確認(rèn)，可確保對復(fù)雜外型的氣動(dòng)特性實(shí)現(xiàn)準(zhǔn)確評估.

1.2 基于工程估算的快速預(yù)測方法

針對空天飛行器的復(fù)雜外形，在應(yīng)用物面壓力估算方法進(jìn)行氣動(dòng)力計(jì)算前，須將其表面劃分為很多小的網(wǎng)格面元.一般劃分為三角形或四邊形網(wǎng)格面元，本研究采用典型的三角形非結(jié)構(gòu)網(wǎng)格面元.

在高超聲速條件下，基于工程估算理論求解飛行器的氣動(dòng)力系數(shù)的基本流程如圖1 所示.該流程主要包括4 個(gè)部分，通過幾何生成軟件定義氣動(dòng)外形，再利用網(wǎng)格生成軟件作非結(jié)構(gòu)網(wǎng)格面元的劃分.將飛行器分解為機(jī)身、機(jī)翼區(qū)域，再結(jié)合迎風(fēng)/背風(fēng)情況，在不同區(qū)域的網(wǎng)格面元上采用適合的方法計(jì)算壓力系數(shù)并求和，最后進(jìn)行數(shù)據(jù)轉(zhuǎn)化，將氣動(dòng)力系數(shù)由機(jī)體系映射為慣性系，轉(zhuǎn)換完畢后存儲(chǔ)或輸出.其中，機(jī)身迎風(fēng)面采用修正牛頓法：

圖1 快速預(yù)測方法的基本流程Fig.1 Basic process of fast predicting method

式中：C p 為壓力系數(shù)，δ 為物面傾角，C pmax為正激波后滯止點(diǎn)的壓力系數(shù).機(jī)身背風(fēng)面采用Prandtl-Meyer 法：

式中：Ma∞為來流馬赫數(shù)，γ0為比熱比.

機(jī)翼背風(fēng)面采用激波膨脹波法，迎風(fēng)面采用切楔/切錐法.在存在攻角 α 的情況下，為了修正物面近似后得到的錐體的物面傾角 δ，采用近似后等價(jià)錐的半頂角 δTC取代物面傾角δ：

式中：φ為等價(jià)錐的徑向角.

采用自主開發(fā)的高超聲速飛行器氣動(dòng)力快速預(yù)測平臺(tái)，它遵循工程估算中氣動(dòng)力的一般求解流程.該平臺(tái)是基于非結(jié)構(gòu)網(wǎng)格文件進(jìn)行計(jì)算求解的，它主要由輸入模塊、幾何處理模塊、求解模塊、數(shù)據(jù)轉(zhuǎn)換模塊、輸出模塊等模塊構(gòu)成.

2 多精度高斯過程回歸

2.1 高斯過程回歸

高斯過程回歸（Gaussian process regression，GPR）是通過對樣本數(shù)據(jù)構(gòu)建高斯過程，再基于貝葉斯理論求解后驗(yàn)概率分布實(shí)現(xiàn)的，本質(zhì)上是利用高斯過程先驗(yàn)知識對樣本數(shù)據(jù)進(jìn)行回歸分析的非參數(shù)回歸模型[20].訓(xùn)練集為觀測值輸入為x，輸出為y.輸入與輸出之間存在一個(gè)未知的映射關(guān)系：

假設(shè)該映射關(guān)系服從均值為0 的高斯分布，則

式中：x'為核函數(shù)中心，k為由超參數(shù) θ 確定的核函數(shù).

由核函數(shù)k確定其協(xié)方差矩陣：

式中：k(xi,xj)為樣本xi、xj之間的協(xié)方差.

協(xié)方差矩陣由核函數(shù)k確定，一般選用徑向基函數(shù)作為高斯過程的核函數(shù)，定義如下：

式中：l0和 σ0為該核函數(shù)的超參數(shù).該超參數(shù)可以通過最大化邊緣對數(shù)似然（marginal log-likelihood）來找到最優(yōu)值.該對數(shù)似然函數(shù)表達(dá)式如下：

式中：n為樣本數(shù)量，p(y|x,θ) 為條件x、θ 下出現(xiàn)y的概率.

本研究采用L-BFGS 算法來求解該優(yōu)化問題的最優(yōu)參數(shù).

2.2 基于高斯過程的多精度回歸實(shí)現(xiàn)

針對不同信息源或不同模型提供的數(shù)據(jù)，根據(jù)模型本身在計(jì)算精度和計(jì)算成本方面的特征，可以將數(shù)據(jù)劃分為多個(gè)精度的數(shù)據(jù)集.其中，多精度高斯過程回歸方法最適合處理由具備如圖2所示特征的模型生成的數(shù)據(jù)，它能夠利用不同精度數(shù)據(jù)實(shí)現(xiàn)多精度建模，獲得的代理模型精度不僅高于低精度模型，而且求解速度遠(yuǎn)快于高精度計(jì)算模型，這是多精度高斯過程回歸模型的優(yōu)勢所在.

圖2 高低精度模型的計(jì)算精度與成本分布Fig.2 Calculation fidelity versus cost distribution for high fidelity and low fidelity model

多精度高斯過程回歸是在高斯過程回歸的基礎(chǔ)上，基于不同精度的數(shù)據(jù)集實(shí)現(xiàn)的.假設(shè)不同精度數(shù)據(jù)集為Dt={(xt,yt)},t=1,···,s.其中，ys表示精度最好的數(shù)據(jù)，y1表示精度最低的數(shù)據(jù).采用自回歸的形式來融合不同精度的數(shù)據(jù)，表達(dá)式如下：

式中：ft(x) 和ft-1(x) 分別表示精度t和精度t-1 的高斯過程回歸模型；ρ 為參數(shù)，將不同精度數(shù)據(jù)集間聯(lián)系起來；δt為隨機(jī)選取的高斯分布，與高斯分布ft-1相互獨(dú)立.假設(shè) δt服從高斯分布，即

式中：μ 表示高斯分布對應(yīng)的均值.基于以上假設(shè)，此時(shí)，精度t的ft(x) 只與下一級精度t-1相關(guān)，而與其他精度的模型無關(guān).因此，該多精度模型在新的輸入點(diǎn)x*的后驗(yàn)分布表達(dá)式如下：

式中：X為觀測值輸入集合.

根據(jù)貝葉斯定理，該后驗(yàn)分布在精度t中的均值和方差表達(dá)式如下：

式中：nt表示精度t的訓(xùn)練數(shù)據(jù)點(diǎn)個(gè)數(shù)，k**=k(x*,x*),.同樣的，該模型的超參數(shù)也可以由最大化式(9)的邊緣對數(shù)似然函數(shù)得到.

3 空天飛行器再入段運(yùn)動(dòng)模型

3.1 運(yùn)動(dòng)模型

在再入飛行過程中，飛行器可以視作質(zhì)量不變的剛體，對其運(yùn)動(dòng)模型作簡化處理，僅考慮它的平移運(yùn)動(dòng).所采用的算例初始再入高度為70 km，飛行器發(fā)動(dòng)機(jī)推力降低為0，由于科氏加速度較小忽略不計(jì)，此時(shí)地球自轉(zhuǎn)影響可忽略不計(jì)，飛行器受到的外力僅為氣動(dòng)力和重力.在再入過程中假設(shè)側(cè)滑角為0°，速度矢量始終保持在飛行器的縱向平面內(nèi).因此，可以建立飛行器在航跡坐標(biāo)系下的三自由度運(yùn)動(dòng)模型如下：

式中：m為飛行器質(zhì)量，V為速度，r為飛行器與地心之間的距離，L為升力，D為阻力，λ 為經(jīng)度，Φ為緯度，σ 為傾斜角，μ 為地球引力常量，γ 為航跡角，ψ 為航向角.其中，升力和阻力表達(dá)式如下：

式中：Cl 為飛行器的升力系數(shù)，Cd 為阻力系數(shù)，Sref為參考面積.

3.2 大氣模型

為了便于軌跡仿真中大氣參數(shù)的獲取，采用基于1976 年標(biāo)準(zhǔn)大氣表發(fā)展出的表達(dá)公式[21]，該數(shù)學(xué)模型以海平面為基準(zhǔn)，將91 km 下的大氣層分為8 層，在每層以高度h為自變量定義了大氣參數(shù)如密度、溫度、氣壓、聲速等.其中，部分高度內(nèi)的大氣密度模型如圖3 所示.隨著高度的下降，大氣密度呈指數(shù)增長，其顯著變化發(fā)生在高度約為30 km 處.

圖3 0～70 km 高度的大氣密度模型Fig.3 Density model of atmosphere between 0 and 70 km

4 仿真分析

4.1 FTB 模型描述

FTB 是由意大利開發(fā)的可重復(fù)使用運(yùn)載飛行器（reusable launch vehicle，RLV），如圖4 所示為FTB 模型的三視圖.機(jī)體總長為7.15 m，翼展為3.6 m，機(jī)翼面積為5.18 m2，鼻錐半徑為0.12 m.機(jī)翼前緣后掠角65°且與機(jī)身腹部呈5°二面角，該設(shè)計(jì)可以增強(qiáng)機(jī)體的橫向穩(wěn)定性.機(jī)身的上下側(cè)具有簡單的錐形、球體幾何形狀和光滑的流線型表面，采用三角形機(jī)翼和V 形尾翼，機(jī)身和機(jī)翼均有尖銳的前緣[21-24].

圖4 FTB 模型三視圖Fig.4 Three views of FTB model

4.2 氣動(dòng)分析方法比較

4.2.1 計(jì)算精度 以FTB 模型為例，比較2 種氣動(dòng)分析方法的求解精度，在馬赫數(shù)等于16 的條件下CFD 數(shù)值模擬、工程估算結(jié)果與文獻(xiàn)數(shù)據(jù)[22]的對比如圖5 所示.2 種方法的結(jié)果隨攻角的變化趨勢基本與文獻(xiàn)數(shù)據(jù)吻合，其中CFD 數(shù)值模擬結(jié)果與文獻(xiàn)數(shù)據(jù)更吻合，工程估算結(jié)果則與文獻(xiàn)數(shù)據(jù)之間存在較大偏差.

圖5 升、阻力系數(shù)的計(jì)算結(jié)果對比Fig.5 Comparison of calculation results of lift and drag coefficients

式中：nt為測試集點(diǎn)數(shù)；為第i個(gè)測試點(diǎn)通過CFD 數(shù)值模擬或工程估算獲得的氣動(dòng)系數(shù)，即升、阻力系數(shù)；為第i個(gè)測試點(diǎn)對應(yīng)的測試集中的數(shù)據(jù)；分別為CFD 數(shù)值模擬和工程估算的平均誤差.2 種方法求解FTB 模型氣動(dòng)系數(shù)的平均誤差如表1 所示.可以看出，升、阻力系數(shù)的計(jì)算中CFD 數(shù)值模擬的精度均比工程估算的高.在升力系數(shù)的計(jì)算中兩者差距更小，而在阻力系數(shù)方面CFD 數(shù)值模擬相對于工程估算的精度提升比率為9.56%.

表1 2 種方法求解氣動(dòng)系數(shù)的平均誤差對比Tab.1 Comparison of average errors of solving aerodynamic coefficients by two methods

4.2.2 計(jì)算成本 2 種方法對于計(jì)算成本的要求如表2 所示.表中，Ns為網(wǎng)格數(shù)量，NC為CPU 使用數(shù)，t為計(jì)算時(shí)間.CFD 數(shù)值模擬中的網(wǎng)格數(shù)量更多，且在時(shí)間、空間維度上迭代求解，故在計(jì)算設(shè)備配置和計(jì)算時(shí)間方面的要求均遠(yuǎn)高于工程估算的，尤其在計(jì)算時(shí)間方面，CFD 數(shù)值模擬所需時(shí)間是工程估算的104數(shù)量級倍.綜合計(jì)算方法的平均誤差和計(jì)算成本的對比結(jié)果，CFD 數(shù)值模擬方法相較于基于工程估算的快速預(yù)測方法可以作為高精度計(jì)算方法使用.

表2 2 種方法求解氣動(dòng)系數(shù)的計(jì)算成本對比Tab.2 Comparison of calculation cost of solving aerodynamic coefficients by two method

4.3 多精度氣動(dòng)建模

根據(jù)前面的結(jié)論，將通過CFD 數(shù)值模擬和基于工程估算的快速預(yù)測方法獲得的數(shù)據(jù)分別用作高精度數(shù)據(jù)和低精度數(shù)據(jù).多精度氣動(dòng)建模如圖6所示，將2 種信息源數(shù)據(jù)通過GPR 完成多精度建模.由496 個(gè)低精度樣本點(diǎn)和42 個(gè)高精度樣本點(diǎn)共同組成多精度氣動(dòng)建模的數(shù)據(jù)集，從數(shù)據(jù)集中隨機(jī)選取42 個(gè)高精度樣本構(gòu)建高精度訓(xùn)練集，選取256 個(gè)低精度樣本構(gòu)建低精度訓(xùn)練集.在-5°≤α ≤25°,5 ≤Ma≤20的參數(shù)空間內(nèi)構(gòu)建二維氣動(dòng)力系數(shù)的多精度預(yù)測模型.

4.3.1 擬合效果 選取高精度訓(xùn)練集中的數(shù)據(jù)為參考值，通過繪制高精度樣本與多精度預(yù)測值間的相關(guān)性圖展現(xiàn)多精度高斯過程回歸對數(shù)據(jù)的擬合效果，如圖7 所示.圖中，C lHF表示升力系數(shù)的高精度樣本，ClLF、ClMF分別為升力系數(shù)的低精度樣本、多精度預(yù)測值，實(shí)線表示完全相關(guān)曲線y=x，散點(diǎn)表征高精度樣本與多精度預(yù)測值的聯(lián)合位置.當(dāng)散點(diǎn)遠(yuǎn)離完全相關(guān)曲線時(shí)，表明在該處的擬合精度較低，反之，擬合精度較高.其中，高精度樣本數(shù)N=0 時(shí)實(shí)際上為僅對256 個(gè)低精度樣本的單精度回歸結(jié)果.在多精度GPR 的結(jié)果圖中，散點(diǎn)與完全相關(guān)曲線幾乎完全重合，表明相較于僅用低精度數(shù)據(jù)的單精度GPR 結(jié)果，多精度GPR 對升力系數(shù)擬合效果更好.

將預(yù)測結(jié)果與訓(xùn)練集中高精度樣本作差得到升力系數(shù)的擬合誤差eCl，并以云圖形式呈現(xiàn).如圖8所示為2 種高精度樣本數(shù)下，預(yù)測模型給出的升力系數(shù)擬合誤差.當(dāng)N=0 時(shí)，低精度代理模型在參數(shù)空間內(nèi)的最大擬合誤差穩(wěn)定在約0.12；當(dāng)N=42 時(shí)，多精度代理模型在參數(shù)空間內(nèi)的最大誤差僅為0.002 5.多精度升力系數(shù)模型的擬合精度得到了極大的提升，將升力系數(shù)的擬合誤差降低至原來的1/48，相對誤差小于1%.

4.3.2 預(yù)測效果 在參數(shù)空間內(nèi)隨機(jī)選取12 個(gè)有別于訓(xùn)練集的位置點(diǎn)，基于CFD 數(shù)值模擬構(gòu)建升力系數(shù)多精度預(yù)測模型的測試集，將該測試集作為參照用以說明多精度高斯過程回歸對參數(shù)空間的預(yù)測能力.如圖9 所示為N=0，42 時(shí)升力系數(shù)測試集數(shù)據(jù)與多精度預(yù)測值之間的相關(guān)性結(jié)果.當(dāng)N=42 時(shí)升力系數(shù)預(yù)測結(jié)果幾乎與完全相關(guān)曲線重合，多精度代理模型對升力系數(shù)表現(xiàn)出良好的預(yù)測能力.

圖9 不同數(shù)量高精度樣本下升力系數(shù)的測試集數(shù)據(jù)與多精度預(yù)測值的相關(guān)性Fig.9 Correlation between data of test set and multi-fidelity predicted values for lift coefficient under different numbers of high-fidelity samples

進(jìn)一步地，以測試集中高精度樣本數(shù)據(jù)為參照，對多精度氣動(dòng)模型的誤差進(jìn)行定量的分析.其中，多精度氣動(dòng)模型的誤差（平方誤差的平均值）定義如下：

基于測試集數(shù)據(jù)計(jì)算得到的多精度氣動(dòng)模型的平方誤差的平均值及平方誤差的方差隨訓(xùn)練集使用的高精度樣本數(shù)變化，如表3、4 所示.同時(shí)，僅使用訓(xùn)練集中高精度樣本通過單精度GPR 構(gòu)建代理模型，再基于測試集數(shù)據(jù)計(jì)算該模型的與多精度結(jié)果形成對比.根據(jù)對比結(jié)果，在FTB 氣動(dòng)外形的多精度升力系數(shù)建模過程中，多精度高斯過程回歸方法基本上遵循隨著使用的高精度樣本數(shù)增加，模型精度也得到提升的規(guī)律，更重要的是它比僅用高精度數(shù)據(jù)構(gòu)建的單精度模型更具精確性和穩(wěn)定性.當(dāng)訓(xùn)練集高精度樣本數(shù)為42、低精度樣本數(shù)為256 時(shí)，升力系數(shù)代理模型平方誤差的平均值降到10-5數(shù)量級水平.此時(shí)，多精度氣動(dòng)模型的預(yù)測誤差波動(dòng)很小，預(yù)測穩(wěn)定性高.

表3 多精度升力系數(shù)模型的平方誤差的平均值Tab.3 Mean of squared error for multi-fidelity model of lift coefficient

表4 多精度升力系數(shù)模型的平方誤差的方差Tab.4 Variance of squared error for multi-fidelity model of lift coefficient

4.3.3 氣動(dòng)建模結(jié)果 類似地，對FTB 模型的阻力系數(shù)也進(jìn)行多精度建模，當(dāng)訓(xùn)練集高精度樣本數(shù)為42、低精度樣本數(shù)為256 時(shí)，升、阻力系數(shù)的多精度建模結(jié)果如圖10 所示.此時(shí)阻力系數(shù)與升力系數(shù)模型預(yù)測效果相仿，在高精度樣本數(shù)小于等于42 范圍內(nèi)達(dá)到最優(yōu).

圖10 升、阻力系數(shù)的多精度建模結(jié)果云圖Fig.10 Cloud figure of multi-fidelity modeling results for lift and drag coefficients

4.4 基于多精度氣動(dòng)模型的軌跡仿真

針對空天飛行器的再入問題，將已經(jīng)構(gòu)建的多精度氣動(dòng)模型應(yīng)用在軌跡仿真中，模擬出FTB 模型在控制攻角條件下的飛行狀態(tài)變化.根據(jù)飛行動(dòng)力學(xué)方程求解流程，構(gòu)建軌跡仿真模型，再入軌跡仿真流程如圖11 所示.其中，初始參數(shù)包括飛行器質(zhì)量m、初始高度h、飛行速度V、參考面積S、航跡角γ、航向角ψ、經(jīng)度λ、緯度Φ、傾斜角σ 等.

圖11 再入軌跡仿真流程Fig.11 Process of reentry trajectory simulation

在再入過程中，存在嚴(yán)重的氣動(dòng)加熱現(xiàn)象、巨大的動(dòng)壓以及過載情況，超過一定限制會(huì)對機(jī)體本身產(chǎn)生破壞，一般氣動(dòng)加熱取駐點(diǎn)處熱流率密度作為參考，熱流率密度、動(dòng)壓 q、過載n定義如下：式中：Rn為駐點(diǎn)曲率半徑，g0為重力加速度，K為常數(shù).

4.4.1 問題描述 在再入過程中，通常會(huì)預(yù)先定義飛行器制導(dǎo)策略，而定義制導(dǎo)策略常見的方式就是使用預(yù)先定義的攻角曲線.這里給出攻角控制方案，假設(shè)它是隨速度變化的函數(shù)，表達(dá)式如下：

在再入過程初期飛行速度和初始高度都較大，初始仿真條件如下：機(jī)體駐點(diǎn)曲率半徑為0.12 m，質(zhì)量為2 000 kg，參考面積為5.18 m2，初始高度為70 km，速度為6 km/s，航跡角為-2°，航向角為0°，起始經(jīng)緯度均為0°，傾斜角σ 為0°.根據(jù)Detra-Hidalgo 模型[21]，計(jì)算駐點(diǎn)熱流率密度的常數(shù)取K=5.16×10-5.由于多精度氣動(dòng)模型基于參數(shù)空間-5°≤α ≤25°，5 ≤Ma≤20 中的高低精度數(shù)據(jù)構(gòu)建，為了保證氣動(dòng)力系數(shù)預(yù)測的有效性，仿真結(jié)束條件設(shè)置為飛行馬赫數(shù)小于5.

4.4.2 軌跡仿真結(jié)果 將僅用42 個(gè)高精度樣本、僅用256 個(gè)低精度樣本以及使用256 個(gè)低精度樣本和42 個(gè)高精度樣本構(gòu)建的3 種氣動(dòng)模型分別作為再入軌跡仿真的氣動(dòng)數(shù)據(jù)輸入，并根據(jù)仿真條件和初始參數(shù)設(shè)置再入軌跡模型.在多精度氣動(dòng)建模中，本研究已經(jīng)從定性和定量的角度證明了多精度相對于單精度模型具有更好的預(yù)測效果，因此這里以基于多精度模型的軌跡仿真結(jié)果為參照，將單精度的仿真結(jié)果分別與多精度的仿真結(jié)果作差獲得再入特征變量的偏差，仿真結(jié)果及偏差對比如圖12 所示.

圖12 基于單、多精度氣動(dòng)模型的再入特征曲線及偏差比較Fig.12 Curve of reentry characteristic and comparison of errors based on single-fidelity and multi-fidelity aerodynamic models

根據(jù)仿真結(jié)果，高度隨時(shí)間的變化呈震蕩曲線下降，這是空天飛行器再入過程中的典型路徑特征，它是下降過程中升阻力系數(shù)、速度、空氣密度等綜合影響的結(jié)果.根據(jù)熱流率密度曲線、過載曲線、動(dòng)壓曲線，熱流率密度和過載為再入初期影響飛行的主要因素，隨著飛行速度下降和高度下降，空氣密度急速上升，動(dòng)壓對飛行器狀態(tài)影響逐漸顯著.偏差對比結(jié)果顯示，通過多精度GPR、高精度GPR 構(gòu)建的氣動(dòng)模型分別作為輸入獲得的再入曲線相差較小，通過低精度GPR 構(gòu)建的氣動(dòng)模型作為輸入獲得的再入曲線與前面兩者相差較大，采用低精度GPR 的最大仿真偏差約為采用高精度GPR 的10～100 倍.具體地，初始高度、速度較大，熱流率密度、動(dòng)壓為速度的函數(shù)，因此低精度GPR 構(gòu)建的氣動(dòng)模型對再入高度、速度以及熱流率密度、動(dòng)壓偏差均帶來顯著影響.其中，熱流率密度偏差絕對值最大至105數(shù)量級，多精度氣動(dòng)建模的工程意義也因此得到體現(xiàn).由于不能獲取軌跡仿真所需要的全部高精度氣動(dòng)數(shù)據(jù)，無法對多精度氣動(dòng)模型獲得的再入仿真結(jié)果給出定量的誤差分析.

5 結(jié)論

(1)以FTB 氣動(dòng)外形為研究對象，分析得出相較于工程估算，CFD 數(shù)值模擬方法求解阻力系數(shù)的精度提升9.56%，而升力系數(shù)的結(jié)果差距略小.但是，CFD 數(shù)值模擬計(jì)算成本要遠(yuǎn)高于工程估算的.

(2)采用多精度高斯過程回歸方法，能以低計(jì)算成本完成該空天飛行器的多精度氣動(dòng)建模，且模型精度較高.當(dāng)訓(xùn)練集使用42 個(gè)高精度樣本、256 個(gè)低精度樣本時(shí)，所提出的多精度氣動(dòng)模型均方誤差低至10-5數(shù)量級，且比僅用相同數(shù)量高精度數(shù)據(jù)構(gòu)建的單精度模型精確度和穩(wěn)定性更高.

(3)將多精度氣動(dòng)模型與再入軌跡仿真耦合，模擬再入過程中的重要狀態(tài)參數(shù)變化，對比分析單、多精度氣動(dòng)模型的仿真結(jié)果，發(fā)現(xiàn)氣動(dòng)模型的微小差異對再入高度、速度、熱流密度、動(dòng)壓影響顯著，進(jìn)一步驗(yàn)證了所提出的多精度氣動(dòng)建模方法對于空天飛行器再入過程軌跡高精度仿真的工程意義.

(4)本研究在多精度氣動(dòng)建模中，并未對訓(xùn)練集和測試集的選用策略作深入考量，最終采用模型是否為全局最優(yōu)模型也未分析論證.在未來研究中可將優(yōu)化過程和多精度氣動(dòng)建模方法相結(jié)合，選用不同采樣策略，通過迭代尋優(yōu)給出目標(biāo)誤差下的最優(yōu)氣動(dòng)模型.