一體化飛行器面向控制的建模與彈道規劃*

黃 俊，劉知貴,3，劉志勤，王慶鳳，張 勇

(1 中國工程物理研究院電子工程研究所，四川綿陽 621900； 2 西南科技大學計算機科學與技術學院，四川綿陽 621010；3 西南科技大學信息工程學院，四川綿陽 621010； 4 中國空氣動力研究與發展中心高超中心，四川綿陽 621000)

0 引言

乘波構形、氣動推進一體化設計[1]在升阻特性和進氣道流量捕獲特性上的優勢，是吸氣式高超聲速飛行器氣動布局的重要選形，該選形飛行器簡稱一體化飛行器。美國的X-51A[2]的一項重要技術驗證即是乘波構型一體化設計且被命名為乘波者。美澳聯合的HIFiRE[3]、法德聯合的JAPHAR[4]等項目開展一體化飛行器的基礎技術研究。我國賀旭照[5]、衛峰[6]等對乘波進氣道的設計方法開展深入研究并取得較好成果。一體化飛行器在設計方法不斷發展的現狀下，對其氣動與控制特性的研究需求日益顯著。文獻[7]對采用少量CFD數據結構工程計算方法對HIFiRE6的飛行動力學進行分析。文獻[8]給出了X-51A的特定狀態的靜穩定性分析。文獻[9]在參數化方法生成一體化飛行器基礎上采用工程計算和多項式擬合方法面向控制建立氣動力模型。相關研究在逐步拓展到控制領域，因此適時開展一體化飛行器的可控制性、可飛性、飛行可達性的度量是必要的。文中建立一體化飛行器的“計算、建模、分析”全程快速方法，作為一體化飛行器設計初期的快速性能評估工具。

近年X-51A相似飛行器是一體化飛行器氣動特性的計算主要研究對象。周正[10]借助X-51A飛行器影像資料構建了相似飛行器并對巡航條件下的氣動特性分析；文獻[11-12]對X-51A相似飛行器的巡航氣動特性進行分析。這些研究結果給出飛行器在巡航條件下的升阻特性，但單一條件的氣動數據，無法支撐控制相關研究工作開展。氣動推進一體化高超聲速飛行器的氣動、推進、控制相互耦合，氣動數據的獲取是一個巨大的難題。例如，美國X-51A[8]項目，在NAART、PSWT等風洞中試驗，風洞時間超1 720 h，運行達3 200次。采用CART3D-Euler等解算器進行數值模擬，總共使用78套網格，運行1 982次。飛行器軌跡、制導和控制等問題研究人員難以獨立完成。因此，面向控制相關問題，有必要開展快速計算型號數據集與建模的方法研究。

意大利航空航天中心(CIRA)[13]利用增量法(build-up approach)對PRPRA-USV1高超聲速飛行器實現風洞試驗與數值計算數據的融合建模，并闡述了方法的有效性。增量法是控制建模經典方法，如美國X34[14]航天飛機中應用增量法建立數據集。文中借鑒上述具體型號研制過程的數據建模思路，通過數值方法模擬X-51A相似外形飛行器的數據集構建并建模，為控制方法研究提供數據基礎。在X-51A研制單位美國空軍研究實驗室(AFRL)的公開報告[8]中對X-51A進行了靜穩定性分析，文中在氣動力模型輸出數據基礎上分析靜穩定性數據，并與之對比，獲得相同結論，表明文中氣動力模型具備該類飛行器氣動特性。

彈道規劃是總體設計中對飛行器可控性、可飛性、可達性分析的重要途徑。文中提出基于偽譜法計算X-51A相似飛行器的最大航程滑翔彈道，同時驗證建立模型的完備性，仿真結果表明建模方法與彈道求解的可行性與有效性。

1 計算模型與方法

以一種與X-51A氣動外形相似的氣動推進一體化飛行器模型為研究對象，如圖1所示。模型參數見文獻[10]。根據后文氣動模型建立的最小數據集需求，參考文獻[8]中彈道特性分析飛行走廊，按飛行走廊設置計算條件，計算條件設置如表1所示。

圖1 一體化飛行器數模

表1 計算條件表

Ma∞H/kmα/(°)δp/(°)3.054.5185.0215.5246.028-10～10間隔20(通氣)-15～15(不通氣)

為了獲得不同高度的氣動數據，對每一個馬赫數補充計算兩組不同高度數據，如表2所示。

表2 計算條件表

文中以Fluent軟件作為計算工具，采用SST方程k-ω模型進行數值求解N-S方程。

舵是最重要的氣動控制機構，為了獲得不同舵偏角下的氣動特性，制作了不同舵偏角模型網格。為了提高計算速度，采用三維半模網格，如圖2所示；同時在計算舵偏數據時采用不通氣半模網格；舵偏角每5°生成一套網格，其中零舵偏有通氣與不通氣半模網格，因此共8套網格，如圖3所示。

圖2 δp=0° 局部體網格

圖3 總共生成8套數模(網格)

2 增量建模方法

首先，定義氣動模型無量綱輸出分別為升力系數、阻力系數、俯仰力矩系數，如式(1)。其中Sref=1 m2,bref=1 m，力矩參考點取(0.55L,0,0)。

(1)

式中：FL為升力，FD為阻力，M為俯仰力矩。

按照前一章節計算，則可獲得數據集：

{Ma，Re,α,δp,CL,CD,Cm}

(2)

式中：Ma表示馬赫數；Re表示雷諾數；α表示攻角；δp表示俯仰通道數值舵偏角。文中模型如圖1所示，采用了全動平尾和雙立尾，物理舵面的控制方式為：

(3)

按照增量法思想，縱向全局氣動力系數可表示為：

(4)

式中：考慮一體化飛行器的耦合特性，右端各項分別代表了基準、俯仰通道的貢獻，兩者的表達式為：

(5)

借鑒X-51A地面試驗經驗，攻角和舵偏角作為基本組合，具有較多的維度值，因此將各個貢獻表示為以下回歸形式：

(6)

針對氣動布局的特點，舵偏角的貢獻進一步簡化，文中忽略雷諾數的影響，忽略側滑角對偏航通道控制的影響，忽略升降舵位置對差動效率的影響，由此得到：

(7)

在基準貢獻中，雷諾數外推是非常關鍵的環節。基本的雷諾數修正數據可通過數值模擬方法得到，那么，變雷諾數氣動力數據表示為以下回歸形式：

Ci(Ma,Re,

(8)

于是，通過平移變換，基準貢獻可以表示為：

Ci(Ma,

(9)

對于不同高度的氣動數據，采用雷諾數修正獲得。飛行高度、大氣參數、雷諾數建立關聯關系：

(10)

式中:ρ、v、L、μ、T、c分別表示密度、速度、飛行器長度、粘性系數、大氣溫度、當地聲速；atmos為參照1976版美國標準大氣表編寫的插值函數。

可見，增量法分別得到各個貢獻的氣動力“數據立方體”。同時，增量法提出了氣動力建模的最小數據集需求，這就為地面試驗、數值模擬的試驗設計提供了依據。

3 計算結果與模型建立

圖4 基準數據集和

圖5 舵效數據集

(11)

按照前述增量建模理論，設計如圖6所示流程構建增量模型。由于氣動數據集變量較多、樣本少，且變量間存在相互耦合的相關關系，建模過程中的擬合采用一種基于偏最小二乘法的多維擬合法來實現。

圖6 增量模型構建流程圖

對于基準數據，考慮雷諾數外推，不可避免帶來誤差。如圖7，在馬赫數為4.5時，源數據是18 km，通過雷諾數外推獲得17 km和19 km的氣動阻力系數，其誤差在10-4級別，符合控制問題研究需求。對于舵效數據已計算全維度的笛卡爾積數據集，直接查表內插值即可獲得任意維度的舵效增益數據，圖8給出模型在馬赫數為4.5時舵效增益。

圖7 基于基準數據模型誤差

圖8 舵效增益模型輸出(Ma=4.5)

4 模型縱向靜穩定度分析

文中采用“瞬時平衡”假定飛行軌跡是定態飛行的。為了使飛行器在某一飛行攻角下縱向平衡，俯仰舵必須偏轉相應角度。如圖9(a)中給出了各個速度和對應高度下的配平曲線δp=f(α)，曲線給出了攻角在-10°～10°范圍的配平舵偏，表明了在此姿態角范圍內，舵偏角能夠使飛行器縱向平衡。

圖9 穩定性分析

5 滑翔彈道規劃

彈道規劃為導彈性能指標和總體設計的可行性提供重要參考依據。前文建立的氣動力模型為彈道計算的重要基礎。下文以最大滑翔距離優化目標的彈道規劃為建立模型的應用算例。滑翔初始高度為27 km、速度為6Ma；末端高度為1 km、速度為3Ma；根據穩態配平結果，全程攻角約束為-8°～8°。建立最優控制問題如下：

(12)

(13)

式中：[h,v,γ]為軌跡狀態量的離散向量，離散點(節點)數為n；FD、FL為離散的阻力和升力向量，由前述增量模型查表插值獲得；wi為積分算子，D為微分矩陣。設置節點數n=32。

Matlab仿真結果如圖10所示，滑翔時間351.2 s，滑翔水平距離484.8 km；δp是對應飛行姿態下保持縱向平衡的俯仰舵偏；系統偽譜法離散正解與仿真反解的對比證明了獲得解的可行性，同時也驗證了文中所建立的氣動力模型足夠支撐該類飛行器的控制特性研究。

如圖11所示，文中的偽譜法改進了控制量插值，采用分段3次埃爾米特插值(PCHIP)，解決了控制量的Lagrange多項插值因Gipps現象使控制量超出控制約束邊界的問題。

6 結論

面向飛行器控制特性研究需求，系統完成了對一種一體化飛行器的氣動數據計算、建模與分析應用的集成研究，獲得結論如下：1)文中的增量法能夠指導算例設計以較小的代價獲得足夠構建面向控制建模的數據集，并建立面向控制的氣動力模型；2)靜穩定度分析結果表明，文中模型具有與X-51A相同的氣動特性，是靜不穩定的，為了維持飛行的姿態平衡，需要及時控制舵偏加以配平；3)提出改進偽譜法計算出可配平攻角約束下，以最大航程為目標優化的可行彈道，估計出文中模型的滑翔部分彈道滑翔水平距離可達484.8 km，驗證了建模方法與彈道求解的可行性與有效性。

圖10 滑翔彈道偽譜法的解及仿真驗證

圖11 PCHIP克服Gipps現象引起控制量越界問題

綜上，提出了一套面向控制的氣動模型建立與彈道規劃方法，解決了多因素耦合的一體化飛行器控制相關研究中數據不足的問題。進一步的工作是開展熱態數據計算工作，建立熱態控制模型，對飛行器的帶動力彈道及控制進行相關分析。