基于氣動-彈道一體化模型的飛行器外形優化設計*

張海瑞，秦 夢，周國峰，王 浩

(中國運載火箭技術研究院，北京 100076)

新型高超聲速飛行器具有飛行速度快、突防能力強、機動范圍大等優勢，是世界主要軍事大國的研究熱點之一。外形優化設計是影響這一類新型飛行器總體性能的關鍵技術之一，旨在尋找一定參數約束條件下飛行器總體性能最優的氣動外形，由于其本身的基礎性和全局性，在總體設計中占有重要地位[1]。結合新型飛行器大空域、寬速域特點，有必要考慮飛行剖面的影響，開展飛行器氣動-彈道一體化外形優化設計。

傳統飛行器外形優化設計方法通常針對典型工況開展氣動外形設計，給出該工況條件下的外形設計方案，在工程中有著廣泛的應用。馬洋[2]基于類型函數/形狀函數轉換(Class/Shape function Transformation,CST)技術，以升阻比、容積率和壓心位置為目標函數，開展了氣動外形的多目標優化設計；Deng等[3]基于二次曲線參數化模型，以升阻比為目標函數，并以其他總體性能參數為約束條件，開展了總體氣動布局優化設計。然而，隨著新型飛行器總體性能的不斷提升，傳統方法難以綜合考慮大空域、寬速域飛行剖面特征對氣動外形適應性的影響，嚴重制約了飛行器綜合性能的有效提升。針對此問題，在傳統氣動外形優化模型基礎上，進一步耦合了彈道的影響。

本文以新型飛行器為研究對象，以外形參數為設計變量，以最大縱程為優化目標，綜合利用CST參數化方法、氣動工程估算方法和Radau偽譜法，建立了飛行器氣動-彈道一體化模型。由于耦合了彈道的影響，極大增加了一體化模型的復雜程度。為進一步提升優化求解效率，引入了Kriging代理模型，提出了一種基于氣動-彈道一體化模型的飛行器外形高效全局優化設計方法。

1 氣動-彈道一體化模型

1.1 氣動-彈道學科分析模型

在新型飛行器的方案設計階段，外形設計參數決定飛行器的氣動特性，并最終影響彈道特性。為了描述這些設計參數對彈道的影響，需要建立一個既能反映系統特點，又能快速分析和優化的學科分析模型。本文的學科分析模型主要包括：外形參數化模型、氣動分析模型和彈道分析模型。

1)外形參數化模型通過一組參數表達和約束飛行器幾何外形[2]。采用CST方法[4-6]描述飛行器底部截面曲線η(ψ)，通過類型函數C(ψ) 定義底部截面曲線的種類，利用形狀函數S(ψ) 修正類型函數定義的幾何外形，采用冪函數[7]控制飛行器俯視面輪廓曲線，進而生成參數化幾何外形，涉及的曲線控制方程如下：

(1)

(2)

其中，(ξ,ψ,η)為飛行器物理坐標(x,y,z)歸一化后的參數坐標。式(1)表示升力體底部截面輪廓控制曲線，Nc為類型函數指數；式(2)表示飛行器俯視面輪廓控制曲線，xn為俯視面輪廓控制指數。

2)氣動分析模型根據飛行器外形以及飛行狀態參數等計算氣動力系數。目前常用的氣動分析方法主要有風洞試驗、數值計算方法和快速工程估算方法。前兩者能夠更為精確地獲取氣動參數，缺點在于成本高、計算周期長。為兼顧精度和快速性的計算要求，本文采用快速工程估算方法求解氣動特性，將飛行器表面劃分面元網格，計算各面元壓力系數，對所有面元上的氣動力進行求和，給出飛行器氣動特性。其中迎風面的計算使用修正牛頓理論，背風面的計算使用Prandtl-Meyer法[8]。

3)彈道分析模型是在氣動分析模型的基礎上，根據飛行器動力學方程實現飛行過程仿真。為引入全空、速域變化對氣動外形優化設計的影響，本文將彈道學科模型處理為軌跡優化問題，在一定的約束條件下，以飛行器縱向射程最大為優化目標，利用Radau偽譜法實現軌跡優化。采用簡化的縱向運動方程作為彈道計算模型，飛行器在此階段只受氣動力和重力影響，不考慮地球的自轉與非球形攝動項等因素帶來的影響。用飛行速度V、速度傾角θ、地心距r和射程角φ描述速度系下的飛行器運動方程：

(3)

式中，CL、CD分別為升力系數和阻力系數，m為飛行器質量，q為飛行動壓，SM為飛行器參考面積。針對縱向動力學方程組，采用Radau偽譜法開展快速軌跡優化。Radau偽譜法[9-11]的基本原理是先將狀態變量和控制變量在(-1,1]上的N個LGR (Legendre-Gauss-Radau)配點上進行離散，然后以N個配點和初始點τ0=-1為節點構造Lagrange插值多項式來逼近狀態和控制變量，把最優控制問題轉化成非線性規劃問題，再使用某種非線性規劃算法求得使性能指標最優的參數，最后代回原問題進行求解。

1.2 一體化優化模型

1)目標函數。外形設計參數對飛行器容積率ηv和氣動性能具有顯著影響[2]，為了結果的可比較性，將飛行器的容積率固定，以最大縱程L為優化目標開展氣動-彈道一體化優化設計。

2)設計變量。CST參數化外形的設計參數中，對升阻比影響最大的5個參數依次為底部截面高度xHl、xHu，俯視面輪廓控制指數xn，底部截面類型函數指數xNcu、xNcl[2]。文獻[10]指出，射程最大的飛行控制方式是最大升阻比，因此設計變量選為X=(xHu,xHl,xNcu,xNcl,xn)。考慮飛行器總體規模及外形尺寸限制，設飛行器最大寬度和最大長度為定值。

(4)

式中：Kn為飛行器相關經驗參數，與飛行器頭部駐點曲率半徑、飛行器外形及材料相關；ρ為大氣密度。

綜上，優化模型可表示為：

(5)

式中，Re為地球平均半徑。

2 高效全局優化方法

由于引入了Radau偽譜法軌跡優化模塊，極大地增加了一體化模型的復雜度，增大了單次迭代的計算時間。針對這一問題，引入Kriging代理模型代替設計空間內的響應值。Kriging代理模型[12-13]是一種基于隨機過程的無偏估計模型，其預測的響應值服從正態分布，不僅能給出對未知函數的預估值，還能給出預估值的誤差估計，在工程中有著廣泛的應用。

采用高效全局優化(Efficient Global Optimization,EGO)方法[14-15]進行求解與優化。EGO方法作為一種全局優化算法，具有很高的效率和精度[16]，其基本原理是通過期望改善函數(Expected Improvement Function,EIF)描述代理模型預測值優于真實值的可能性，并在EIF最大處加入樣本點以更新代理模型。

圖1 高效全局優化方法流程Fig.1 Flow chart of efficient global optimization

然后，以初始樣本點為輸入、對應的最大縱程L為輸出，構建Kriging代理模型為：

(6)

(7)

式中的均值和方差都可以由式(8)[14]計算：

(8)

(9)

EIF(X)=E(P(X))

(10)

(11)

3 算例

以升力體構型飛行器開展優化設計，飛行器質量為m=900 kg，參考面積SM考慮飛行器底部最大截面積，容積率ηv不小于0.292。飛行器最大寬度和長度分別為4 m和2.5 m，外形設計變量取值范圍如表1所示。

表1 外形設計變量及變化范圍Tab.1 Design variables and ranges of variation

表2 偽譜法狀態變量初值和終值設置Tab.2 Initial and final value settings of state variables in pseudospectral method

利用高效全局優化方法進行氣動-彈道一體化優化設計。在初始階段采用全因子采樣方案對初始參數設計空間進行等密度采樣，構建初始Kriging代理模型。為對比驗證，利用高效全局優化方法，開展Ma=10、α=8°工況下的氣動外形優化設計，構建初始代理模型的方法同上，進而將優化結果代入彈道模型，采用Radau偽譜法以最大射程為目標開展快速軌跡優化。優化結果對比如表3所示，外形對比如圖2所示。

表3 氣動優化與氣動-彈道一體化優化結果對比Tab.3 Comparison of results of aerodynamic-only optimization with aerodynamic and trajectory integrated optimization

由表3和圖2可以看出，兩種方法給出的俯視面輪廓控制指數xn相同，俯視面輪廓外形完全一致；一體化優化方法給出的底部截面高度相對較低，通過優化底部截面類型函數，在保證容積率相同條件下，獲取較優氣動外形。

為了進一步說明氣動-彈道一體化模型的優勢，給出AeroOpt-8和JointOpt升力體的升阻比隨攻角變化曲線，如圖3所示。可以看到，當攻角α=7.71°時，AeroOpt-8升力體的升阻比最大，為4.604 8；當攻角α=7.39°時，JointOpt升力體的升阻比最大，為4.602 3。AeroOpt-8和JointOpt兩者的最大升阻比相同，然而JointOpt升力體的攻角-升阻比曲線更平緩，這是因為AeroOpt-8只考慮了攻角α=8°時的優化外形，而氣動-彈道一體化模型耦合了彈道，考慮了飛行過程攻角擾動對外形設計的影響，得到全空、速域下的優化外形。

(a) 氣動優化外形(a) Aerodynamic optimization shape

圖3 升阻比隨攻角變化曲線對比Fig.3 Comparison of lift-to-drag ratio curves

分別求解AeroOpt-8和JointOpt升力體的氣動特性，最大射程對比如表4所示，偽譜法控制飛行攻角的變化曲線對比如圖4所示。JointOpt的最大升阻比略小于AeroOpt-8的最大升阻比，但前者的最大縱向射程比后者提升了28.04 km；JointOpt升力體的飛行攻角整體小于AeroOpt-8，二者的飛行攻角分別在各自最大升阻比對應攻角附近小幅度振蕩。結合圖表可以說明，當飛行攻角在最大升阻比對應攻角附近發生擾動時，JointOpt升阻比的降低量比AeroOpt-8更小，前者受到攻角擾動的影響比AeroOpt-8更小，說明采用高效全局優化算法針對氣動-彈道一體化模型開展優化設計可以提升飛行器對工況變化的適應能力，避免因攻角擾動而影響升阻比最終導致射程的降低，進一步驗證了氣動-彈道一體化模型能夠綜合考慮全空、速域下優化外形的結論。

圖4 飛行攻角變化曲線對比Fig.4 Comparison of the angle of attack

表4 最大升阻比及最大縱向射程對比Tab.4 Comparison of the maximum lift-to-drag ratio and maximum trajectory

4 結論

1)以新型飛行器為研究對象，為考慮大空域、寬速域飛行剖面特征對氣動外形適應性的影響，采用CST參數化建模方法結合快速工程估算方法獲取氣動特性，通過偽譜法給出彈道特性，建立飛行器氣動-彈道一體化模型，提升了外形優化模型精度。

2)針對這一模型的特點，以最大縱向射程為優化目標，以外形參數為設計變量，并明確約束條件，給出了基于Kriging的氣動外形優化方法，實現在多參數約束條件下的飛行器氣動外形高效全局優化，提升了復雜外形優化模型的求解效率。

3)以升力體構型飛行器為例，開展氣動外形優化設計，結果表明，氣動-彈道一體化模型能夠較好描述大空域、寬速域的氣動和彈道特征，有效提升新型飛行器氣動外形設計精度和水平，具有一定的工程應用價值。