吳超，譚劍鋒，王浩文，倪先平

(1.南京航空航天大學(xué)直升機(jī)旋翼動(dòng)力學(xué)國家級重點(diǎn)實(shí)驗(yàn)室，江蘇南京210016;2.清華大學(xué) 航天航空學(xué)院，北京100084)

0 引言

直升機(jī)配平是空氣動(dòng)力學(xué)、飛行動(dòng)力學(xué)和旋翼動(dòng)力學(xué)分析的基礎(chǔ)，對直升機(jī)姿態(tài)響應(yīng)計(jì)算、穩(wěn)定性分析及飛控系統(tǒng)的設(shè)計(jì)與開發(fā)都具有重要作用。直升機(jī)配平問題為多元非線性方程組求解問題，通常需通過數(shù)值求解得到。直升機(jī)配平傳統(tǒng)采用牛頓法［1］，該方法簡單、高效，但需根據(jù)飛行狀態(tài)采用不同的簡化模型計(jì)算初值。文獻(xiàn)［2］根據(jù)初始預(yù)測值將配平分為縱向、橫向和旋翼轉(zhuǎn)速三部分，該算法比較復(fù)雜，且需引入阻尼項(xiàng)構(gòu)建一個(gè)依賴工程經(jīng)驗(yàn)的算法。文獻(xiàn)［3］采用廣義簡約梯度優(yōu)化法，但該算法配平結(jié)果依賴于初值，且易收斂到局部最優(yōu)解。文獻(xiàn)［4］綜合了梯度下降法和擬牛頓法，一定程度上降低了配平結(jié)果對初值的依賴，但該算法涉及大量偏導(dǎo)數(shù)，因而計(jì)算量偏大。文獻(xiàn)［5］將浮點(diǎn)遺傳算法計(jì)算的全局較優(yōu)解作為擬牛頓法的初值，而后計(jì)算最優(yōu)配平解，雖不依賴于初值，但解空間偏大，因而以一定概率收斂到非物理配平解。

非線性方程組解的存在性及唯一性仍未完全解決［6］，雖然直升機(jī)配平全局?jǐn)?shù)值解并不唯一，但直升機(jī)穩(wěn)定飛行狀態(tài)是存在的，因而直升機(jī)優(yōu)化配平解必須為物理解。本文在取值范圍約束和非線性約束的條件下，將遺傳算法(GA)計(jì)算的全局較優(yōu)解作為LM(Levenberg-Marquardt)算法的初值，利用LM算法快速計(jì)算物理最優(yōu)解，由此建立基于GA/LM混合優(yōu)化算法的直升機(jī)全機(jī)配平方法。

1 飛行動(dòng)力學(xué)模型

直升機(jī)飛行動(dòng)力學(xué)模型是配平算法的基礎(chǔ)，典型的氣動(dòng)力部件包括旋翼、尾槳、尾翼及機(jī)身等部分。

1.1 旋翼動(dòng)力學(xué)模型

為突出配平算法研究，本文采用均勻入流模型和剛性槳葉假設(shè)。由揮舞鉸力矩平衡得到旋翼揮舞動(dòng)力學(xué)方程［7］:

式中，a=［a0，a1，b1］T為一階諧波形式的揮舞運(yùn)動(dòng)系數(shù);D，K分別為阻尼和剛度矩陣;f為外力向量，是前進(jìn)比μ、入流比λ、機(jī)體角速度、角加速度以及操縱輸入的函數(shù)。

將所有槳葉的葉素氣動(dòng)載荷沿半徑和方位角積分，并投影到機(jī)體坐標(biāo)軸系下，得到槳榖六力素:

式中，θ0，θ1c和 θ1s分別為總距和兩個(gè)周期變距;λ 隱含在旋翼入流方程中。

1.2 尾槳?jiǎng)恿W(xué)模型

尾槳無周期變矩且錐度角為常值［8］。與主旋翼槳轂六力素相似，尾槳六力素:

式中，aTR，θ0TR，μTR和 λTR分別為尾槳揮舞運(yùn)動(dòng)向量、尾槳槳距、尾槳前進(jìn)比和入流比;λTR隱含在尾槳入流方程中，尾槳揮舞運(yùn)動(dòng)向量可通過尾槳揮舞運(yùn)動(dòng)方程來計(jì)算。

1.3 尾翼和機(jī)身動(dòng)力學(xué)模型

平尾、垂尾的升力和阻力分別為:

式中，qHS，qVF分別為平尾和垂尾處的動(dòng)壓;SHS，SVF分別為平尾和垂尾的面積;CLHS，CLVF分別為平尾和垂尾的升力系數(shù);CDHS，CDVF分別為平尾和垂尾的阻力系數(shù)。對于采用全動(dòng)平尾的UH-60A黑鷹直升機(jī)，假設(shè)平尾安裝角僅是前進(jìn)比的函數(shù)，并由飛測數(shù)據(jù)［9］擬合得到，垂尾安裝角為零。

由于機(jī)身沒有精確的理論模型，因此本文根據(jù)風(fēng)洞吹風(fēng)試驗(yàn)數(shù)據(jù)［9］計(jì)算機(jī)身六力素fF。

1.4 全機(jī)重力模型

將全機(jī)重力投影到機(jī)體坐標(biāo)系下可得fG。

1.5 全機(jī)飛行動(dòng)力學(xué)模型

全機(jī)動(dòng)力學(xué)方程如下:

式中，v=［u，v，w］T和 ω =［p，q，r］T分別為線速度和角速度;f和m分別為合外力和合外力矩;I為全機(jī)慣性矩陣;Ω為叉乘向量算子。定義如下:

運(yùn)動(dòng)學(xué)方程為:

式中，α=［φ，θ，ψ］T為歐拉角;E 為角速度至歐拉角速度的轉(zhuǎn)換矩陣，定義如下:

式(1)和式(5)～式(7)構(gòu)成了直升機(jī)的非線性狀態(tài)方程組:

2 配平算法

2.1 配平定義和對應(yīng)初始條件

配平狀態(tài)對應(yīng)于狀態(tài)量導(dǎo)數(shù)為常數(shù)或零且操縱輸入為固定值的狀態(tài)。直升機(jī)穩(wěn)定飛行時(shí)為零為定值。若為非零常數(shù)，則α將隨時(shí)間線性變化，由此直升機(jī)穩(wěn)定狀態(tài)包含了懸停、勻速平飛、勻速直線爬升、螺旋爬升及穩(wěn)定自轉(zhuǎn)下滑等狀態(tài)，具有廣泛的意義。

由于穩(wěn)定飛行時(shí)的機(jī)體角加速度為零，式(1)退化為:

式中，下標(biāo) e表示平衡狀態(tài);Xe=［ve，ωe，αe，ae］T為配平狀態(tài)量;Ue為配平操縱量。

狀態(tài)方程式(10)由12個(gè)非線性方程構(gòu)成，而待求變量包含12個(gè)狀態(tài)量和4個(gè)控制量，為得到唯一確定解，需給定4個(gè)運(yùn)動(dòng)約束量［2］:飛行速度VT、爬升角γ、轉(zhuǎn)彎速率和側(cè)滑角βe。這4個(gè)量的不同組合對應(yīng)了不同的穩(wěn)定飛行狀態(tài)，如表1所示。

本文將16個(gè)未知量分為9個(gè)獨(dú)立變量XI和7個(gè)依賴變量 XD(ve，ωe，θe)。依賴變量是獨(dú)立變量和運(yùn)動(dòng)約束量的函數(shù)，將錐度角和周期揮舞系數(shù)作為獨(dú)立變量，如表2所示。

表1 穩(wěn)定飛行狀態(tài)和運(yùn)動(dòng)約束量對照表Table 1 List of motion constraints corresponding to steady state

表2 獨(dú)立變量和依賴變量表Table 2 Independent variables and dependent variables

構(gòu)造目標(biāo)函數(shù):

當(dāng)J=0時(shí)，XI為式(11)的精確解，由此將直升機(jī)配平問題轉(zhuǎn)化為優(yōu)化問題，即尋找XI使得J最小。

2.2 GA/LM混合優(yōu)化算法

2.2.1 解空間的確定

為避免優(yōu)化算法收斂到非物理配平解，本文將構(gòu)建配平解空間約束條件。對UH-60A而言，旋翼總距范圍為［9.9°，25.9°］，縱向周期變距范圍為［－8.0°，8.0°］，橫向周期變距范圍為［－ 12.5°，16.3°］，尾槳范圍為［4.5°，36.5°］［9］，旋翼揮舞角范圍為［－ 6°，25°］［10］，俯仰角和滾轉(zhuǎn)角范圍為［－90°，90°］，則配平解空間約束條件可表示為:

2.2.2 GA/LM 混合算法

依據(jù)遺傳算法，本文將目標(biāo)函數(shù)式(12)作為適用度函數(shù)，設(shè)定種群數(shù)ΝP=100，交叉概率pc=0.8，變異概率pm=0.2，迭代次數(shù)N=100，數(shù)值精度為雙精度。

計(jì)算步驟如下:

(1)在式(13)解空間范圍內(nèi)隨機(jī)產(chǎn)生ΝP個(gè)個(gè)體;(2)將十進(jìn)制個(gè)體轉(zhuǎn)變?yōu)槎M(jìn)制;(3)選擇操作:利用適應(yīng)度優(yōu)的生存概率大，采取輪盤賭方式選擇適應(yīng)度優(yōu)的個(gè)體(適用度函數(shù)小);(4)交叉操作:在種群中隨機(jī)選出兩個(gè)個(gè)體，隨機(jī)確定交換位置交換基因，得到新一代個(gè)體;(5)變異操作:個(gè)體中的某一個(gè)基因以一定概率發(fā)生變異;(6)根據(jù)式(12)計(jì)算十進(jìn)制個(gè)體的適應(yīng)度，若達(dá)到收斂精度，則將最新一代適應(yīng)度最優(yōu)的個(gè)體作為LM算法的初值，否則跳回到步驟(2)繼續(xù)計(jì)算。

在GA給出的較優(yōu)解基礎(chǔ)上求解式(12)，目標(biāo)函數(shù)為平方和函數(shù)，可得到雅可比矩陣，將其代入牛頓法得到高斯牛頓法，但高斯牛頓法需計(jì)算在當(dāng)前步k的逆，可能出現(xiàn)不可逆現(xiàn)象。為此LM算法增加了系數(shù)阻尼矩陣［11］，求解得出自變量的增量，其中λk為阻尼系數(shù)，I為單位矩陣。由此計(jì)算出，依次迭代直至達(dá)到收斂精度。LM算法的搜索方向是牛頓高斯和最速下降法方向的混合，當(dāng)λk=0時(shí)搜索方向ΔXIk與高斯牛頓法相同;當(dāng)λk→∞時(shí)搜索方向ΔXIk與最速下降法相同。

GA算法全局搜索性能好，不易陷入局部最優(yōu)解，且易并行化［12］，但當(dāng)平均適應(yīng)度函數(shù)和最佳適應(yīng)度函數(shù)相差較小時(shí)，GA算法尋優(yōu)速度非常緩慢。為提高收斂速度，將GA給出的較優(yōu)解作為LM的初值，利用LM算法快速收斂性提高直升機(jī)優(yōu)化配平速度。由此，本文通過結(jié)合GA和LM算法各自優(yōu)點(diǎn)建立了適用于直升機(jī)配平的高效全局優(yōu)化算法，計(jì)算流程如圖1所示。

圖1 主程序流程圖Fig.1 Main program flow chart

3 計(jì)算結(jié)果與分析

3.1 物理最優(yōu)解和常用優(yōu)化算法對比

對式(12)無約束優(yōu)化的求解，采用了8種單一優(yōu)化算法。UH-60A以20 m/s作無側(cè)滑的平飛，初值Ⅰ:θ0=22.92°，θ1c= － 0.57°，θ1s=4.58°，θ0TR=21.03°，其余初值為0;初值Ⅱ均為零。計(jì)算結(jié)果如表3和表4所示。

從表3和表4中可以看出，單形調(diào)優(yōu)法總迭代步數(shù)較多，且該算法收斂速度依賴于初值選取;擬牛頓法和最速下降法收斂速度極快，但均對初值要求較高;BFGS擬牛頓法對初值要求較低;LM算法收斂速度最快，信賴域反射法收斂速度次之;模擬退火算法和無解空間約束的遺傳算法均不收斂，不適用于直升機(jī)配平問題。

表5中給出了達(dá)到收斂精度算法的獨(dú)立變量值XI。對于初值Ⅰ，單形調(diào)優(yōu)法、BFGS擬牛頓法、LM法和信賴域反射法均收斂到表5中的解1;對于初值Ⅱ，BFGS擬牛頓法、LM法和信賴域反射法數(shù)值均收斂到表5中的解2。從表5中可以看出，存在兩個(gè)數(shù)值解，因此目標(biāo)函數(shù)收斂是直升機(jī)配平的必要條件而非充要條件。根據(jù)式(13)的解空間易知解1是配平解。

表3 優(yōu)化計(jì)算結(jié)果對比1Table 3 Results of optimization algorithm part 1

表4 優(yōu)化計(jì)算結(jié)果對比2Table 4 Results of optimization algorithm part 2

表5 配平值(°)Table 5 Results of trim

圖2和圖3為隨機(jī)一次GA/LM混合算法的初值和快速收斂計(jì)算的迭代曲線圖，混合算法在解空間范圍內(nèi)經(jīng)多次計(jì)算均收斂到解1。圖2反映GA算法在解空間內(nèi)能很快計(jì)算出給定精度的初值，解空間約束的引入改善了無約束GA算法的收斂性，加快了收斂速度。圖3說明了LM算法的快速收斂性，最終目標(biāo)函數(shù)J=6.6×10－30。顯然可以看出，采用GA/LM混合算法進(jìn)行約束優(yōu)化求解直升機(jī)的配平問題是有效的。

圖2 初值計(jì)算迭代曲線Fig.2 Iterative curve of initial value

圖3 快速收斂迭代曲線Fig.3 Iterative curve of fast convergence

3.2 算例驗(yàn)證和分析

協(xié)調(diào)轉(zhuǎn)彎是最典型的穩(wěn)定飛行狀態(tài)。本文對UH-60A以51.44 m/s在1600 m高度的水平面內(nèi)作不同轉(zhuǎn)彎速率的協(xié)調(diào)轉(zhuǎn)彎進(jìn)行配平計(jì)算，并將4個(gè)操縱量轉(zhuǎn)換為操縱桿位置［9］。采用GA/LM混合算法計(jì)算的配平曲線如圖4～圖8所示。

圖4 總距桿協(xié)調(diào)轉(zhuǎn)彎配平曲線Fig.4 Trim results of collective lever for coordinated turn

圖5 橫向駕駛桿協(xié)調(diào)轉(zhuǎn)彎配平曲線Fig.5 Trim results of lateral stick for coordinated turn

圖6 縱向操縱桿協(xié)調(diào)轉(zhuǎn)彎配平曲線Fig.6 Trim results of longitudinal stick for coordinated turn

圖7 腳蹬協(xié)調(diào)轉(zhuǎn)彎配平曲線Fig.7 Trim results of pedal for coordinated turn

圖8 俯仰角協(xié)調(diào)轉(zhuǎn)彎配平曲線Fig.8 Trim results of pitch angle for coordinated turn

從圖4～圖8中可以看出，本文計(jì)算的總距和腳蹬與GEN HEL［13］計(jì)算值比較吻合;由于全動(dòng)平尾的安裝角是由實(shí)驗(yàn)數(shù)據(jù)給出，因此本文計(jì)算的縱向操縱桿和俯仰角比GEN HEL計(jì)算值更接近飛測數(shù)據(jù);橫向駕駛桿在較大負(fù)滾轉(zhuǎn)角時(shí)與試驗(yàn)數(shù)據(jù)有較大誤差，主要原因是腳蹬位置與試驗(yàn)數(shù)據(jù)存在一定的靜態(tài)誤差。

本算例計(jì)算中假設(shè)側(cè)滑角為零，而飛行測試時(shí)飛行員很難保證，因此在受側(cè)滑角影響明顯的橫向方向上，本文計(jì)算結(jié)果與飛行測試結(jié)果存在一定偏差。

4 結(jié)束語

在直升機(jī)動(dòng)力學(xué)模型基礎(chǔ)上，通過構(gòu)造優(yōu)化目標(biāo)函數(shù)，將直升機(jī)全機(jī)配平問題轉(zhuǎn)化為優(yōu)化問題。針對單一優(yōu)化算法對初值依賴高，易收斂到非物理解問題，構(gòu)建了解空間約束條件，提出了一種基于GA/LM混合優(yōu)化的直升機(jī)配平算法。通過UH-60A直升機(jī)對協(xié)調(diào)轉(zhuǎn)彎的計(jì)算，表明本文所建立的混合優(yōu)化算法準(zhǔn)確高效，且不依賴于初值，適用于直升機(jī)配平分析。

［1］ 高正，陳仁良.直升機(jī)飛行動(dòng)力學(xué)［M］.北京:科學(xué)出版社，2003:65-71.

［2］ Padfield G D.Helicopter flight dynamics:the theory and application of flying qualities and simulation modeling［M］.Oxford:Blackwell Publishing Ltd，2007:192-204.

［3］ Schank T C.Optimal aeroelastic trim for rotorcraft with constrained non-unique trim solutions［D］.Atlanta:Georgia Institute of Technology，2008.

［4］ Subramanian S，Gaonkar G H，Maier T H.A theoretical and experimental investigation of hingeless-rotor stability and trim［C］//Twenty-Third European Rotorcraft Forum.Dresden，Germany，1997:10.

［5］ 代冀陽，吳國輝，朱國民.適于直升機(jī)配平計(jì)算的混合遺傳算法［J］.飛行力學(xué)，2010，28(1):24-28.

［6］ 李廣楊，關(guān)治，白峰杉.數(shù)值計(jì)算原理［M］.北京:清華大學(xué)出版社，2000:243.

［7］ Talbot P D，Tinling B E，Decker WA，et al.A mathematical model of a single main rotor helicopter for piloted simulation［R］.NASA-TM-84281，1982.

［8］ Hilbert K B.A mathematical model of the UH-60 helicopter［R］.NASA-TM-85890，1984.

［9］ Howlett J J.UH-60A black hawk engineering simulation program［R］.NASA CR-166309，1981.

［10］ Cross J，Brilla J，Kufeld R，et al.The modern rotor aerodynamic limits survey:a report and data survey［R］.NASATM-4446，1993.

［11］龔純，王正林.精通MATLAB最優(yōu)化計(jì)算［M］.北京:電子工業(yè)出版社，2012:213-216.

［12］張呈林，郭才根.直升機(jī)總體設(shè)計(jì)［M］.北京:北京國防工業(yè)出版社，2006:141-143.

［13］ Ballin M G.Validation of a real-time engineering simulation of the UH-60A helicopter［R］.NASA-TM-88360，1987.