黃胤錚，王曉喆，王柳青

（1.北京航空航天大學(xué) 飛行學(xué)院；2.北京航空航天大學(xué) 無人系統(tǒng)研究院，北京 100083）

0 引言

無人機(jī)因具有體積小、成本低、對惡劣環(huán)境適應(yīng)性好等優(yōu)點(diǎn)，在民用和軍用領(lǐng)域均有巨大的應(yīng)用潛力［1］。但無人機(jī)的作戰(zhàn)能力會受到飛機(jī)結(jié)構(gòu)、氣動布局、任務(wù)載荷、巡航時間等因素影響，無人機(jī)任務(wù)完成度以及生存概率也與其飛行性能和作戰(zhàn)性能密切相關(guān)。隨著航空技術(shù)的發(fā)展，未來空域的環(huán)境復(fù)雜度都會有不同程度的上升，對無人機(jī)的整體性能提出了更高要求。

隨著有效載荷不斷增加，無人機(jī)的結(jié)構(gòu)重量系數(shù)不斷下降，結(jié)構(gòu)柔度不斷增加，氣動、結(jié)構(gòu)等學(xué)科的耦合關(guān)系進(jìn)一步增強(qiáng)，若在優(yōu)化過程中采用串列優(yōu)化方法會造成學(xué)科間的反復(fù)迭代，在降低設(shè)計(jì)效率的同時，也影響了整體性能的提升。氣動彈性優(yōu)化充分考慮了氣動與結(jié)構(gòu)之間的耦合效應(yīng)，并以此挖掘設(shè)計(jì)潛力［2］，實(shí)現(xiàn)飛行器在復(fù)雜環(huán)境下經(jīng)濟(jì)性和使用性能的提升［3］。在大型客機(jī)［4-5］和超音速客機(jī)［6］的設(shè)計(jì)中，氣動彈性優(yōu)化已經(jīng)得到了較為廣泛的應(yīng)用?；诿舳刃畔⒌膬?yōu)化算法在求解大規(guī)模的多學(xué)科優(yōu)化問題時，收斂較快［7-8］，而以遺傳算法為代表的進(jìn)化式算法需要大量樣本，不適合直接應(yīng)用于多設(shè)計(jì)變量的氣動/結(jié)構(gòu)耦合優(yōu)化問題中［9］。1988 年，Jameson［10］首次將伴隨方程法應(yīng)用于氣動目標(biāo)函數(shù)梯度的求解過程中，并對機(jī)翼進(jìn)行優(yōu)化設(shè)計(jì)。由于伴隨理論的計(jì)算消耗與設(shè)計(jì)變量的數(shù)量無關(guān)，適合處理具有大規(guī)模設(shè)計(jì)變量的多學(xué)科優(yōu)化問題，因此近年來得到了廣泛關(guān)注［11-14］。Anderson 等［15］采用多學(xué)科建模優(yōu)化方法和伴隨理論，實(shí)現(xiàn)了風(fēng)力渦輪機(jī)葉片的高保真結(jié)構(gòu)優(yōu)化和負(fù)載應(yīng)力最小化；Wang 等［16］采用基于伴隨方法的優(yōu)化設(shè)計(jì)框架對機(jī)翼受動載荷進(jìn)行了分析優(yōu)化；Batay 等［17］采用多學(xué)科設(shè)計(jì)優(yōu)化工具，結(jié)合伴隨求解器，對風(fēng)力渦輪機(jī)進(jìn)行了并行氣動設(shè)計(jì)優(yōu)化；Wilke［18］采用自主開發(fā)的多目標(biāo)優(yōu)化框架實(shí)現(xiàn)了對直升機(jī)旋翼槳葉的優(yōu)化設(shè)計(jì)；李潤澤等［19］在多目標(biāo)優(yōu)化過程中結(jié)合伴隨方法，提高了超臨界機(jī)翼氣動優(yōu)化設(shè)計(jì)效率；劉曉東等［20］采用基于伴隨理論的氣動優(yōu)化設(shè)計(jì)方法，實(shí)現(xiàn)了飛翼布局飛行器的氣動優(yōu)化。

隨著現(xiàn)代無人機(jī)對高機(jī)動以及長航時性能的需求不斷提高，采用基于伴隨理論的氣動彈性優(yōu)化方法深度挖掘多學(xué)科耦合所帶來的設(shè)計(jì)潛力是必要且迫切的。本文結(jié)合渦格法和有限元分析法，通過松耦合的方式實(shí)現(xiàn)氣動彈性分析，并基于序列二次規(guī)劃和伴隨理論開展非線性敏度優(yōu)化。分別針對長航時和高機(jī)動兩類典型無人機(jī)，以最小燃油消耗為目標(biāo)，對機(jī)翼的氣動外形和空間梁結(jié)構(gòu)等參數(shù)進(jìn)行優(yōu)化設(shè)計(jì)，并對優(yōu)化設(shè)計(jì)的結(jié)果進(jìn)行對比分析。

1 氣動彈性仿真

1.1 氣動渦格法

對于飛行器的氣動分析本質(zhì)上是對空間N-S 方程的求解。渦格法在位勢流方程的求解過程中較為常見，與空間N-S 方程的求解相比，利用基于求解位勢流的氣動計(jì)算方法能夠以更高的計(jì)算效率對啟動導(dǎo)數(shù)進(jìn)行比較準(zhǔn)確的求解［21］。

根據(jù)邊界不可穿透定理，可以在控制點(diǎn)處添加馬蹄渦與氣動表面相切的條件，如式（1）所示。其中，CC為控制點(diǎn)處誘導(dǎo)速度的系數(shù)矩陣，Γ為每個渦線渦強(qiáng)組成的列向量，V∞為控制點(diǎn)來流速度矩陣，N為控制點(diǎn)處單位法向量。

通過求解式（1）可得到每個渦的渦流強(qiáng)度和氣動力，如式（2）所示。其中，F(xiàn)n為渦段上的升力，ρ∞為來流密度，V∞為來流速度，Vn為當(dāng)前渦段中點(diǎn)誘導(dǎo)速度，Γn為當(dāng)前渦段上的渦強(qiáng)。

渦格法渦線布置方法如圖1所示。

Fig.1 Vortex lattice method vortex line layout method圖1 渦格法渦線布置方法

1.2 結(jié)構(gòu)有限元法

對于組成機(jī)翼的各部分而言，機(jī)翼上下表面的蒙皮主要承載機(jī)翼上的彎曲載荷，機(jī)翼的翼梁和副梁主要承載剪切載荷，蒙皮和梁共同承載扭轉(zhuǎn)載荷。在分析過程中，將機(jī)翼翼梁簡化為空間梁單元的有限元模型進(jìn)行分析?？臻g梁單元的有限元模型如圖2所示［22］。

Fig.2 Finite element model of space beam element圖2 空間梁單元有限元模型

通過對空間梁模型的分析，可以獲得結(jié)構(gòu)分析部分的控制方程如下：

其中，f為作用在節(jié)點(diǎn)上的力和力矩，d為節(jié)點(diǎn)位移列陣。通過求解式（3），可以得到各節(jié)點(diǎn)上的靜力學(xué)載荷分布。

1.3 仿真模型

氣動彈性的耦合求解采用松耦合方法，氣動與結(jié)構(gòu)之間的數(shù)據(jù)傳遞如圖3、式（4）和式（5）所示［23］。

Fig.3 Aerodynamic load transfer圖3 氣動載荷傳遞

其中，i=1，2，分別代表氣動網(wǎng)格對應(yīng)有限元模型中的左右兩個結(jié)構(gòu)節(jié)點(diǎn)，T為氣動網(wǎng)格上通過渦格法計(jì)算得到的合力，rsai為結(jié)構(gòu)節(jié)點(diǎn)指向氣動力作用點(diǎn)的方向向量，F(xiàn)i、Mi分別為通過該牽引方案得到的節(jié)點(diǎn)力和節(jié)點(diǎn)扭矩。

在完成氣動載荷傳遞后，通過有限元分析可以得到對應(yīng)的結(jié)構(gòu)位移，并傳遞回氣動網(wǎng)格。

其中，i=1，2，分別代表氣動網(wǎng)格對應(yīng)有限元模型中的左右兩個結(jié)構(gòu)節(jié)點(diǎn)，di為平動位移，θi為轉(zhuǎn)動位移，rsai同式（5），為結(jié)構(gòu)節(jié)點(diǎn)指向氣動力作用點(diǎn)的方向向量。在此基礎(chǔ)上，結(jié)合開源框架OpenMDAO 建立無人機(jī)的氣動、結(jié)構(gòu)以及氣動彈性耦合分析模型［24］。

2 優(yōu)化算法及框架

2.1 伴隨理論

傳統(tǒng)的導(dǎo)數(shù)計(jì)算方法以有限差分法和復(fù)變量差分法為主，在梯度優(yōu)化中，都需要對控制方程反復(fù)求解，計(jì)算量大。而伴隨理論算法可以在保持精度的前提下，一次性求解出目標(biāo)函數(shù)對所有設(shè)計(jì)變量的導(dǎo)數(shù)［3］。其計(jì)算量僅與目標(biāo)函數(shù)和約束條件的數(shù)量有關(guān)，與設(shè)計(jì)變量的數(shù)量無關(guān)，適合多設(shè)計(jì)變量的優(yōu)化設(shè)計(jì)問題。

本文的無人機(jī)氣動彈性優(yōu)化以氣動結(jié)構(gòu)耦合系統(tǒng)作為分析系統(tǒng)，以氣動參數(shù)和結(jié)構(gòu)參數(shù)作為設(shè)計(jì)變量。優(yōu)化目標(biāo)函數(shù)I和控制方程殘差R可設(shè)定為：

其中，x表示設(shè)計(jì)變量；w、u表示氣動/結(jié)構(gòu)耦合系統(tǒng)狀態(tài)變量，其中w具體表現(xiàn)為渦格法中網(wǎng)格單元的氣動負(fù)載，u具體表現(xiàn)為有限元分析中結(jié)構(gòu)有限元分析得到的結(jié)構(gòu)位移；Ra對應(yīng)氣動分析的殘差，Rs對應(yīng)結(jié)構(gòu)分析的殘差。

令優(yōu)化目標(biāo)函數(shù)I和控制方程殘差R分別對設(shè)計(jì)變量x進(jìn)行求導(dǎo)，并利用恒等變換和伴隨算子ψ得到伴隨方程。通過求解伴隨方程，可以一次性求解出目標(biāo)函數(shù)針對所有設(shè)計(jì)變量的導(dǎo)數(shù)，如式（8）所示：

2.2 序列二次規(guī)劃

常用的梯度優(yōu)化方法包括BFGS 擬牛頓算法、共軛梯度法、序列二次規(guī)劃算法（Sequential Quadratic Programming，SQP）等。其中，序列二次規(guī)劃法對大規(guī)模數(shù)據(jù)的計(jì)算效率較高，在工程中得到了廣泛應(yīng)用。具體形式如式（9）所示：

其中，f（x）為目標(biāo)函數(shù)，h（x）為等式約束函數(shù)，g（x）為不等式約束函數(shù)。

序列二次規(guī)劃法是一種用于求解非線性最優(yōu)化問題的方法，其基本思路是將原問題轉(zhuǎn)化為一系列二次規(guī)劃（Quadratic Programming，QP）問題進(jìn)行迭代求解。在每次迭代中，通過求解二次規(guī)劃問題得到下一次迭代的值，直到序列｛xk｝收斂于極值點(diǎn)。具體來說，對于某一次迭代值xk，NLP 會被近似為該點(diǎn)處的二次規(guī)劃問題。通過求解該二次規(guī)劃問題，可以得到下一次迭代的值xk+1。隨著迭代次數(shù)k的增加，序列｛xk｝將收斂于極值點(diǎn)［25］。本文采用SLSQP 算法（Sequential Least SQuares Programming），該算法是由Kraft［26］于1988 年提出的一種非線性約束優(yōu)化算法，通過將優(yōu)化問題轉(zhuǎn)化為一系列線性或二次規(guī)劃子問題，用于求解無約束或約束的非線性優(yōu)化問題。具體來說，SLSQP 算法采用牛頓法求解每個子問題，并利用輔助函數(shù)法處理約束條件。

2.3 氣動彈性優(yōu)化框架

本文所采用的氣動彈性優(yōu)化框架如圖4 所示。具體流程如下：①設(shè)定優(yōu)化目標(biāo)并設(shè)置目標(biāo)參數(shù)，明確目標(biāo)函數(shù)、約束函數(shù)、設(shè)計(jì)變量與狀態(tài)變量等相關(guān)變量及數(shù)據(jù)；②根據(jù)初始數(shù)據(jù)建立氣動和結(jié)構(gòu)的物理模型；③對目標(biāo)函數(shù)和約束函數(shù)關(guān)于所有設(shè)計(jì)變量的總梯度進(jìn)行計(jì)算；④利用序列二次規(guī)劃算法進(jìn)行尋優(yōu)；⑤通過不斷更新氣動和結(jié)構(gòu)的物理模型，直至得到收斂解，完成優(yōu)化設(shè)計(jì)。

Fig.4 Basic process of optimization design based on gradient optimization algorithm圖4 基于梯度優(yōu)化算法的優(yōu)化設(shè)計(jì)基本流程

3 無人機(jī)氣動彈性優(yōu)化

3.1 大過載無人機(jī)

在有人機(jī)上，8g是受過專業(yè)訓(xùn)練的航天員一般可以承受的最大過載［27］。而無人機(jī)沒有相關(guān)限制，可以做出更多大過載的機(jī)動動作來規(guī)避潛在危險。

3.1.1 工況計(jì)算

大過載無人機(jī)具體的飛行工況如表1 所示。選取NACA0015 翼型，機(jī)翼初始參數(shù)如表2 所示，并選取7075 鋁合金作為翼梁結(jié)構(gòu)材料。

Table 1 High overload UAV flight state parameters表1 大過載無人機(jī)飛行狀態(tài)參數(shù)

Table 2 Initial shape parameters of high overload UAV wing表2 大過載無人機(jī)機(jī)翼初始外形參數(shù)

目標(biāo)函數(shù)為巡航狀態(tài)下的燃油消耗量，對于飛機(jī)燃油消耗，采用Breguet 距離方程進(jìn)行計(jì)算。具體方程形式如下：

其中，Wf為燃油重量，W0為飛機(jī)空重，Ws為結(jié)構(gòu)重量，R為最大航程，V為巡航速度，SFC為燃油消耗率。設(shè)計(jì)變量為機(jī)翼梢根比、機(jī)翼后掠角、機(jī)翼扭轉(zhuǎn)角、機(jī)翼翼梁厚度、機(jī)動狀態(tài)下的機(jī)翼迎角。設(shè)計(jì)變量具體數(shù)量及數(shù)值范圍如表3 所示。約束條件為：機(jī)動狀態(tài)下機(jī)翼結(jié)構(gòu)應(yīng)力接近但不超過應(yīng)力極限、機(jī)動狀態(tài)下升力等于重力。

Table 3 Parameter range and number of design variables for high overload UAV表3 大過載無人機(jī)設(shè)計(jì)變量參數(shù)范圍及數(shù)量

3.1.2 結(jié)果優(yōu)化

經(jīng)過優(yōu)化設(shè)計(jì)后，可得到如表4 所示的3 種優(yōu)化設(shè)計(jì)結(jié)果。無人機(jī)機(jī)翼的外形經(jīng)過優(yōu)化設(shè)計(jì)后的變化如圖5所示。

Table 4 Optimal design results of high overload UAV表4 大過載無人機(jī)優(yōu)化設(shè)計(jì)結(jié)果

Fig.5 Optimal design results of UAV with different wing root chord lengths圖5 不同機(jī)翼翼根弦長無人機(jī)優(yōu)化設(shè)計(jì)結(jié)果

通過對不同機(jī)翼翼根弦長情況下機(jī)翼優(yōu)化設(shè)計(jì)結(jié)果的綜合分析，可以看出，在機(jī)翼經(jīng)過優(yōu)化設(shè)計(jì)后，大過載無人機(jī)的燃油消耗量相較于初始狀態(tài)均降低了20%～30%。以機(jī)翼翼根弦長5.5 m 為例，如圖6 所示，在經(jīng)過優(yōu)化設(shè)計(jì)后，可以看到馮米塞斯應(yīng)力在巡航狀態(tài)下遠(yuǎn)低于應(yīng)力極限，符合基本邏輯和設(shè)計(jì)要求。在扭轉(zhuǎn)角方面，機(jī)翼翼尖和翼根的相對扭轉(zhuǎn)角為2.89°，符合機(jī)翼外形設(shè)計(jì)邏輯，改善了升力分布，有利于提升飛機(jī)的巡航性能［28］。設(shè)計(jì)方案總體上較為合理，符合飛行器機(jī)翼外形的設(shè)計(jì)邏輯。

Fig.6 Optimal design results of wing with root chord length of 5.5m圖6 5.5m機(jī)翼翼根弦長優(yōu)化設(shè)計(jì)結(jié)果

3.2 長航時無人機(jī)

3.2.1 工況計(jì)算

對于長航時無人機(jī)，具體的飛行工況如表5 所示。選取NACA0015 翼型，采用如圖7 所示的矩形機(jī)翼，初始參數(shù)如表6所示，并選取7075鋁合金作為翼梁結(jié)構(gòu)材料。

Table 5 Flight state parameters of long-endurance UAV表5 長航時無人機(jī)飛行狀態(tài)參數(shù)

Table 6 Initial parameters of long-endurance UAV wing表6 長航時無人機(jī)機(jī)翼初始參數(shù)

目標(biāo)函數(shù)為燃油消耗量，設(shè)計(jì)變量為機(jī)翼梢根比、機(jī)翼后掠角、機(jī)翼扭轉(zhuǎn)角、機(jī)翼翼梁厚度、巡航狀態(tài)下的機(jī)翼迎角。設(shè)計(jì)變量的具體數(shù)量及數(shù)值范圍如表7 所示。約束條件為：機(jī)翼結(jié)構(gòu)應(yīng)力不超過應(yīng)力極限、巡航狀態(tài)下升力等于重力。

Table 7 Parameter range and number of high overload UAV flight status表7 高過載無人機(jī)飛行狀態(tài)參數(shù)范圍及數(shù)量

3.2.2 結(jié)果優(yōu)化

長航時無人機(jī)優(yōu)化后的機(jī)翼幾何外形如圖8所示。

Fig.8 Optimized wing geometry of long-endurance UAV圖8 長航時無人機(jī)優(yōu)化后的機(jī)翼幾何外形

初始狀態(tài)下的燃油消耗量為3 898 kg，長航時優(yōu)化后的燃油消耗量為3 004 kg，優(yōu)化后的結(jié)構(gòu)重量為1 761 kg。與初始狀態(tài)相比，燃油消耗量降低了22.93%。優(yōu)化設(shè)計(jì)后的具體參數(shù)如表8 所示。長航時情況下的優(yōu)化設(shè)計(jì)結(jié)果如圖9 所示。在扭轉(zhuǎn)角方面，機(jī)翼翼尖和翼根的相對扭轉(zhuǎn)角為5.15°，處于合理范圍內(nèi)。在馮米塞斯應(yīng)力方面，由于需要獲得最佳的優(yōu)化設(shè)計(jì)結(jié)果，因此需要在設(shè)計(jì)過程中盡可能逼近結(jié)構(gòu)極限。本文的優(yōu)化設(shè)計(jì)方案未超出應(yīng)力極限，符合飛行器外形設(shè)計(jì)邏輯。

Table 8 Optimized design results of long endurance UAV表8 長航時無人機(jī)優(yōu)化設(shè)計(jì)結(jié)果

Fig.9 Optimized design results of long endurance UAV圖9 長航時情況下的優(yōu)化設(shè)計(jì)結(jié)果

4 結(jié)語

本文提出一種現(xiàn)代無人機(jī)的高效氣動彈性優(yōu)化方法，并通過大過載和長航時兩類典型無人機(jī)進(jìn)行了驗(yàn)證，具體如下：

（1）提出基于伴隨理論的氣動彈性優(yōu)化方法，通過氣動渦格法與有限元分析法進(jìn)行氣動彈性松耦合分析，并基于伴隨理論和序列二次規(guī)劃法進(jìn)行非線性敏度優(yōu)化。

（2）針對大過載無人機(jī)進(jìn)行結(jié)構(gòu)優(yōu)化設(shè)計(jì)，以燃油消耗量為目標(biāo)函數(shù)，以機(jī)翼外形參數(shù)和機(jī)翼翼梁結(jié)構(gòu)參數(shù)為設(shè)計(jì)變量，結(jié)合機(jī)動狀態(tài)下氣動、結(jié)構(gòu)兩方面的約束函數(shù)，使無人機(jī)的燃油消耗量分別減少了21.66%、27.93%、32.05%。

（3）針對長航時無人機(jī)進(jìn)行結(jié)構(gòu)優(yōu)化設(shè)計(jì)，以燃油消耗量為目標(biāo)函數(shù)，以機(jī)翼外形參數(shù)和機(jī)翼翼梁結(jié)構(gòu)參數(shù)為設(shè)計(jì)變量，結(jié)合巡航狀態(tài)下氣動、結(jié)構(gòu)兩方面的約束函數(shù)，使無人機(jī)的燃油消耗量減少了22.93%。

綜上所述，本文提出的方法高效、穩(wěn)健、可靠，時間成本和學(xué)習(xí)成本低，可為無人機(jī)氣動彈性優(yōu)化提供理論和方法參考。