基于變精度模型的變外形飛行器彈道優化

王健磊, 陳曉宇, 洪厚全, 李春娜, 龔春林, 付俊興

(1.西北工業大學 航天學院, 陜西 西安 710072; 2.江西洪都航空工業集團有限責任公司, 江西 南昌 330024)

在未來戰爭中,越發復雜的戰場環境和飛行任務對飛行器提出了更嚴格的要求。但傳統飛行器外形固定不變的模式,無法保證飛行器在大速度包絡任務中具有很好的表現。因此,變外形飛行器的概念應運而生。

在變外形飛行器變形規律研究方面,Oliviu等[1]給出了上表面可變形翼的數值模擬和風洞試驗結果,對變形規律的研究有一定的參考價值。Koreanschi等[2-3]通過對翼型的優化設計,研究了翼型變形對減阻的影響。這些對于變外形飛行器氣動性能的探索,為后來的飛行器變形機構和變形規律研究奠定了基礎。呂吉嬋等[4]將遺傳算法與CAD和CFD軟件相結合,得到了變后掠機翼在中、低空下的最優變后掠規律。劉璐等[5]將CFD數值模擬與Matlab軟件相結合,研究了不同高度下變后掠翼飛行器的最佳變形規律。但如果結合彈道來研究變形規律,直接采用試驗或CFD數值模擬提供氣動性能會導致計算負擔過重或計算耗時過長。近年來蓬勃發展的變精度代理模型技術可用于變外形飛行器的氣動性能建模[6-7],能夠在彈道優化求解過程中提供氣動輸入,有效節約計算資源,減少計算耗時。

在彈道優化方面,直接打靶法、粒子群算法和偽譜法等都是常用的最優控制方法。王方鵬等[8]將攻角作為控制變量,利用打靶法得到了最優攻角變化規律。但直接打靶法僅適用于較簡單、精度要求不高的最優控制問題[9],李偉明和郭杰等[10-11]利用粒子群算法,研究了變后掠飛行器的增程能力與發射條件的關系。但粒子群算法收斂精度較低、收斂速度慢,會進一步增加變外形飛行器彈道優化問題的計算負擔。偽譜法在求解最優控制問題時融合了間接法和直接法的優點,在飛行器彈道優化方面逐漸成為研究熱門。其中,hp自適應偽譜法對初始猜測值不敏感,適應性強,能夠有效節省計算時間,適用于處理控制變量較多的問題,在高超聲速飛行器最優彈道和再入彈道優化,以及變外形飛行器變形規律研究等方面都有較好的效果。龔春林等[12]基于hp自適應偽譜方法,對翼型分布式變體問題進行了彈道優化,獲得了更省燃料質量的飛行方案,但該飛行器的變形部位是翼型,變形量較小。王娜等[13]針對后掠角和攻角2個控制變量,采用hp自適應偽譜法,優化得到了飛行器的變形規律,但優化時并沒有考慮飛行器的發動機性能。

本文基于較多的低精度CFD仿真結果和少量的高精度CFD仿真結果,建立變精度代理模型,以此模型預測變外形飛行器在不同變形狀態和工況下的氣動性能,運用hp自適應偽譜法對變外形飛行器開展了以燃料消耗最少為優化目標的彈道優化,并與固定翼飛行器進行了對比分析。

1 問題描述

本文提出了一種能夠進行低亞聲速和跨聲速巡航的變外形飛行器,其彈道剖面如圖1所示。

圖1 彈道剖面示意圖

飛行器首先以較低的初始速度從地面起飛,然后爬升到高速巡航高度;完成高速巡航任務后,飛行器減速并下降到低速巡航高度;完成低速巡航任務后,飛行器再次下降并加速對目標進行壓制打擊。根據飛行任務要求,彈道參數要求為:①高速巡航時高度大于14 km,巡航速度不小于0.9Ma;②低速巡航時高度不大于3 km,巡航速度不大于0.35Ma;③壓制打擊高度為0.1 km,打擊速度為0.8Ma;④總航程不小于600 km;⑤總巡航時間不低于40 min。因此,需要進行優化的飛行彈道可以分為3個主要階段:加速爬升階段、減速下降階段和加速下降階段。

為完成飛行任務,本文提出了一種后掠角和軸向位置可變的變形機翼,機翼后掠角的變化可以使飛行器兼顧亞/跨聲速性能,而機翼軸向位置的變化可以保證飛行器的縱向靜穩定性。考慮到機翼結構強度和指標限制要求,機翼后掠角變化范圍為0°～42°,軸向位置變化范圍為0～370 mm。機翼通過剪切變形來改變形狀,因此翼型始終與來流保持平行。為了使飛行器具有合理的縱向靜穩定度,當后掠角增大時,機翼向前移動,反之亦然。該飛行器的外形示意圖如圖2所示。機翼變形的極限位置俯視圖如圖3所示,分別為后掠角為0°和42°情況。另外,除飛行器外形參數外,在進行彈道優化時,還需考慮飛行器發動機控制參數和彈道上的飛行攻角。

圖2 變外形飛行器外形示意圖圖3 機翼極限位置俯視圖

2 變外形飛行器最優彈道研究方法

2.1 求解流程

對于變外形飛行器,要研究其最優彈道,首先要獲得各個外形在不同飛行狀態下的氣動性能,然后進行彈道優化。本文中研究的變外形飛行器,機翼的后掠角和軸向位置在飛行過程中都可以改變,且飛行狀態和幾何外形的變化范圍較大,會導致直接基于CFD仿真的彈道優化的計算量過大。為了提高彈道優化效率,本文提出一種基于變精度Kriging模型的流程,采用大量低精度CFD仿真數據和少量高精度CFD仿真數據建立氣動性能模型,以此來預測各個外形在不同飛行狀態下飛行器的氣動性能。

基于變精度模型的變外形飛行器彈道優化,是利用上述建立的變精度模型提供氣動輸入,利用hp自適應偽譜法控制飛行器攻角、發動機參數和變形參數,獲得最優變形規律。研究的具體流程如圖4所示。

圖4 彈道優化研究流程

1) 根據確定任務剖面,選擇典型飛行狀態,包括飛行速度和飛行攻角;

2) 在選擇的典型飛行狀態下,根據變形翼的2個變形參數(后掠角和機翼軸向位置)選擇用于高、低精度仿真分析的變形構型;

3) 分別采用Euler和N-S求解器對低精度和高精度分析的構型進行CFD仿真,得到高、低氣動力數據;

4) 在不同飛行狀態下,建立變外形飛行器的變精度Kriging模型;

5) 分析飛行狀態和變形參數對飛行器氣動性能的影響,確定控制變量與狀態變量;

6) 將氣動性能模型作為輸入,利用偽譜法對飛行器的彈道進行優化。

2.2 變精度Kriging模型

對于實際工程應用,即使采用代理模型直接對高精度數據建模,有時也是難以承受的,因為高精度分析的耗時非常大。變精度模型提供了一種更為高效的建模方法[14-19],采用大量代價少的低精度數據預測模型趨勢,采用少量昂貴的高精度數據修正模型的局部精度。變精度Kriging代理模型(multi-fidelity Kriging,MFK)采用了具有一階自動遞歸的MFK模型來建模[20-21],與多層Kriging(hierarchical Kriging,HK)模型類似[22-24]。在本文中,采用無黏歐拉方程求解的氣動力作為低精度數據,采用N-S方程求解的氣動力作為高精度數據,只有2級精度。因此建立的MFK模型可以簡化成如下形式

yhigh(x)=ρ(x)·ylow(x)+δ(x)

(1)

式中:ρ(x)為高、低精度數據間的尺度/相關因子,可以是常數、線性函數或二次函數;δ(x)為差異函數;ylow(x)為低精度數據建立的代理模型;yhigh(x)為高精度數據建立的代理模型。

在(1)式中,代理模型采用了Kriging模型,用于高、低精度數據的建模。Kriging模型對非線性問題具有很好的預測能力,而且能夠在建模的同時獲得模型的預測誤差。Kriging模型的具體建模過程參考文獻[25-26]。δ(x)是與低精度模型ylow(x)無關的高斯過程。ρ(x)代表了高、低精度模型之間的相關程度,體現了高、低精度模型尺度的關系,可以通過(2)式進行計算

(2)

式中:cov[yhigh(x),ylow(x)]為高、低精度模型之間的協方差矩陣;var[ylow(x)]為低精度模型的方程矩陣。

采用Kriging模型替代高、低精度模型有利于協方差和方程計算,適用于該變精度建模方法。

2.3 基于hp自適應偽譜法的最優彈道問題求解

考慮到變外形飛行器彈道優化問題設計空間復雜,選擇直接法中計算性能較高的hp自適應偽譜法[27-28]進行彈道優化計算。hp自適應偽譜法結合了偽譜法和hp型有限元法的優點,能夠最大限度地保證近似精度,同時降低計算量。

常規偽譜法將狀態變量和控制變量在若干配點上離散,然后以各離散配點為節點,利用Lagrange插值基函數來全局近似狀態變量和控制變量。這樣,連續的最優控制問題可以轉化為離散的非線性規劃(nonlinear programming,NLP)問題。為了避免常規偽譜法設計變量數目過多的問題,自適應偽譜法將整個連續問題分成多個片段,在每個片段內應用常規偽譜法,其迭代示意圖如圖5所示,步驟簡要介紹如下:

1) 初始化求解問題,分割初始網格時間片段,對配置點數和容許偏差等參數進行設置;

2) 采用常規偽譜法對各個時間片段求解NLP問題;

3) 檢驗每個時間片段中得到的變量能否滿足容許偏差,若滿足則計算停止，否則執行步驟4);

4) 對于不滿足容許偏差的時間片段,若殘差向量為一致型則增加配點數量，否則在該時刻增加時間片段;

5) 更新所有時間片段,返回步驟2)。

圖5 hp自適應偽譜法迭代流程

3 結果分析

3.1 氣動性能建模結果

在建立MFK模型時,在飛行包絡內,選取的馬赫數為{0.35,0.5,0.65,0.8,0.95},攻角為{0°,2°,4°,6°,8°},共25種飛行狀態。針對每種狀態,考慮到飛行器的氣動性能以及后掠角和機翼軸向位置的限制,選擇不同外形進行CFD計算。經過分析共選擇了240個點進行無黏氣動性能計算,110個點來進行有黏氣動性能計算。本節將這些無黏和有黏氣動性能計算的結果分別作為低精度和高精度數據來進行氣動性能建模。

對于0.35Ma情況時,機翼全展開為最優外形;對于0.95Ma情況時,機翼達到最大后掠角為最優外形。另外考慮到壓力中心的限制,在0.35Ma和0.95Ma情況時,機翼應分別位于可移動范圍的最后端和最前端,即2 600 mm處和2 230 mm處。由此確定了這2個馬赫數下氣動外形,無需建立這2種工況下的變精度Kriging模型。對于其他飛行狀態,在每個狀態下分別建立升力系數、阻力系數和壓力中心3個變精度Kriging模型,共建立45個模型。另外,升阻比模型不是直接建模得到的,而是用升力系數和阻力系數模型計算得到。

在建立變精度Kriging模型時,使用了開放源碼工具箱(surrogate model toolbox,SMT)[29]。以0.65Ma,0°攻角情況為例,建立的氣動性能代理模型如圖6所示。從圖中可以看出,變外形飛行器的升阻力系數隨后掠角的增大而減小,機翼軸向位置對升阻力系數的影響很小,對壓力中心影響明顯。

圖6 0.65Ma時0°攻角的MFK模型

MFK模型建立后,采用留一交叉驗證方法對所有的飛行狀態下升力系數、阻力系數和壓力中心的模型精度進行驗證。不同飛行狀態的評價指標RMSE和RRMSE如表1～6所示。

從表中可以看出,多數RMSE和RRMSE較小,其中75%以上的RRMSE不超過5%,且最大的RRMSE為8.06%,在工程上可以接受。因此,變精度Kriging模型具有足夠的精度,可以用于后續的優化和氣動輸入。

表1 不同飛行狀態下Cl的RMSE

表2 不同飛行狀態下Cd的RMSE

表3 不同飛行狀態下xp的RMSE

表4 不同飛行狀態下Cl的RRMSE %

表5 不同飛行狀態下Cd的RRMSE %

表6 不同飛行狀態下xp的RRMSE %

3.2 彈道優化條件設置

本文進行彈道優化時,基于瞬時平衡假設,只研究飛行器質心在縱向平面內的運動。分析3.1中飛行器的氣動性能可知,影響變外形飛行器氣動性能的因素主要是后掠角,機翼軸向位置對氣動性能的影響很小。在本節中,僅研究后掠角對飛行器彈道的影響。為保證靜穩定性,將飛行器機翼位置定義為后掠角的線性函數,后掠角增大時機翼由最末端移向最前端。

飛行器彈道優化的一個重要指標是燃料消耗,故以最小化燃料消耗為目標函數,即

J=maxms(α,n,χ)

(3)

約束條件為:

1) 初始約束條件:飛行器某段彈道投放初始時刻t0時狀態變量的初值,包括飛行器初始點的位置(x0,y0)、初始速度v0和初始彈道傾角θ0,根據飛行任務給定約束,即

(4)

2) 終端約束條件:飛行器彈道末端時刻狀態變量的值,包括飛行器末點的位置(xf,yf)、末點速度vf和末點彈道傾角θf,根據飛行任務給定約束,即

(5)

3) 過程約束條件:飛行器實際飛行過程中狀態變量的值,在飛行過程中過大的飛行速度會導致過載過大,破壞飛行器機體結構;同時,由于希望飛行器平穩上升和下降,要對飛行器的彈道傾角進行約束,故過程約束包括飛行速度vmax和彈道傾角θ,即

(6)

4) 控制變量約束條件:適當限制控制變量的變化范圍可以保證彈道的平滑和良好的飛行性能,除飛行器的攻角α和發動機轉速n外,還包括飛行器的變形量后掠角χ,機翼軸向位置xa由后掠角χ得到,根據飛行任務給定約束,即

(7)

式中,ξ指各個變量到參考點的偏差值。對于3段需要進行優化的彈道來說,每段彈道約束條件分別設置如下。

1) 加速爬升段:初始條件m0=1 000 kg,x0=0,y0=0,v0=0,θ0=20°;末端約束設置為xf≤100 km,yf≥14 km,vf≥0.9Ma,|θf|≤20°;過程約束條件設置為vf≤1.2Ma,|θf|≤60°;控制變量約束條件設置為8°≤α≤14°,0°≤χ≤42°,28 000 r/min≤n≤36 000 r/min。

2) 減速下降段:初始條件m0=888 kg,x0=0,y0=14 km,v0=0.95Ma,θ0=-10°;末端約束設置為xf≤100 km,yf≤3 km,vf≤0.35Ma,|θf|≤20°;過程約束條件設置為vf≤1.2Ma,|θf|≤60°;控制變量約束條件設置為4°≤α≤14°, 0°≤χ≤42°, 24 000 r/min≤n≤36 000 r/min。

3) 加速下降段:初始條件m0=800 kg,x0=0,y0=3 km,v0=0.35Ma,θ0=40°;末端約束設置為xf≤100 km,yf≤0.1 km,vf≥0.8Ma,|θf|≤50°;過程約束條件置為vf≤1.2Ma,|θf|≤60°;控制變量約束條件設置為-4°≤α≤0°,0°≤χ≤42°,24 000 r/min≤n≤36 000 r/min。

3.3 彈道優化結果

下文對固定翼飛行器和變外形飛行器的3段優化彈道及彈道參數進行對比分析。其中,固定翼飛行器是變外形飛行器機翼后掠角為42°、軸向位置在最前端時(見圖3b))的外形。

1) 加速爬升段

固定翼飛行器和變外形飛行器加速爬升段的彈道及參數變化如圖7～13所示。

圖7 彈道對比 圖8 速度對比 圖9 質量對比

圖10 攻角對比 圖11 發動機轉速對比 圖12 后掠角變化

圖13 軸向位置變化

從圖中可以看出,變外形飛行器能夠在較短時間和射程內迅速加速爬升到14 km高度,且能達到0.9Ma的巡航速度要求。與固定翼飛行器相比, 節省了40%的爬升時間。固定翼飛行器和變外形飛行器加速爬升段末尾質量分別為930和950 kg,變外形飛行器能夠節省28.57%的燃料。2種飛行器攻角和發動機轉速雖然最后都穩定于同一定值,但在彈道的前半段固定翼飛行器的攻角和發動機轉速整體略高于變外形飛行器。這是由于飛行器發射速度較低,固定翼飛行器42°的后掠角使其升力較小,需要更大的攻角和推力來完成加速爬升任務。

2) 減速下降段

固定翼飛行器和變外形飛行器減速下降段的彈道及參數變化如圖14～20所示。

圖14 彈道對比 圖15 速度對比 圖16 質量對比

圖17 攻角對比 圖18 發動機轉速對比 圖19 后掠角變化

圖20 軸向位置變化

變外形飛行器能夠在較短時間和射程內迅速減速下降到飛行任務要求的高度,也能達到巡航速度要求,節省了25%的下降時間。固定翼飛行器和變外形飛行器加速爬升段末尾的質量分別為860和882 kg,分別消耗了10和6 kg的燃料,故在減速下降段變外形飛行器能夠節省40%的燃料。變外形飛行器的攻角整體高于固定翼飛行器,波動更少。這是由于變外形飛行器的燃料消耗較少,質量高于固定翼飛行器,因此需要更大的升力。

3) 加速下降段

固定翼飛行器和變外形飛行器加速下降段的彈道及參數變化如圖21～27所示。

圖21 彈道對比 圖22 速度對比 圖23 質量對比

圖24 攻角對比 圖25 發動機轉速對比 圖26 后掠角變化

圖27 軸向位置變化

從圖中可以看出,變外形飛行器能夠在較短時間內下降到飛行任務要求的高度,也能達到巡航速度要求,節省了20.83%的時間,能夠更迅速地實施突防打擊。固定翼飛行器和變外形飛行器加速下降段末尾的質量分別為767和790 kg,分別消耗了11和10 kg的燃料,故在加速下降段變外形飛行器能夠節省9.09%的燃料。變外形飛行器加速性能更好,在彈道開始前期就已經加速到了較大的馬赫數,因此整體射程更長。變外形飛行器的攻角整體低于固定翼飛行器,且波動更少。這是由于變外形飛行器機翼處于全展開狀態,整體升力更大,若攻角過大飛行器無法滿足下降要求。

綜合以上對比結果可知,變外形飛行器能夠通過靈活改變自身外形來改變氣動特性,從而能夠花費較短的時間滿足不同飛行任務的要求。同時,從完整的飛行彈道來看,相同飛行目標下變外形飛行器比固定翼飛行器節省9.87%燃料。故而變外形飛行器在寬速域、多任務的方案飛行彈道中具有更大的優勢,通過兼顧高速和低速飛行,能夠執行更復雜的戰場任務。

4 結 論

為避免浪費過多的計算資源,提高變外形飛行器彈道優化效率,本文提出了一套基于MFK建模的最優彈道求解流程,通過留一交叉驗證,證明了MFK模型對預測氣動性能的準確性,并采用了hp自適應偽譜法對固定翼飛行器和變外形飛行器進行了彈道優化,研究結果表明:

1) 與直接采用CFD仿真相比,采用MFK模型的計算成本大大降低,計算時間減少86%。

2) 較傳統固定翼飛行器,變外形飛行器通過兼顧變形參數、攻角和發動機控制參數等變量的控制,能夠保持在多任務條件下最優氣動特性飛行。

3) 與固定翼飛行器的最優彈道相比,變外形飛行器能節省9.87%的燃料,且飛行器的上升下降效率更高。

4) 該流程為研究變外形飛行器沿彈道的最優變形規律提供了一種有效途徑。