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

(1.南京航空航天大學直升機旋翼動力學國家級重點實驗室，江蘇南京210016;2.清華大學 航天航空學院，北京100084)

0 引言

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

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

1 飛行動力學模型

直升機飛行動力學模型是配平算法的基礎，典型的氣動力部件包括旋翼、尾槳、尾翼及機身等部分。

1.1 旋翼動力學模型

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

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

將所有槳葉的葉素氣動載荷沿半徑和方位角積分，并投影到機體坐標軸系下，得到槳榖六力素:

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

1.2 尾槳動力學模型

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

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

1.3 尾翼和機身動力學模型

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

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

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

1.4 全機重力模型

將全機重力投影到機體坐標系下可得fG。

1.5 全機飛行動力學模型

全機動力學方程如下:

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

運動學方程為:

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

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

2 配平算法

2.1 配平定義和對應初始條件

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

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

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

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

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

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

表2 獨立變量和依賴變量表Table 2 Independent variables and dependent variables

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

當J=0時，XI為式(11)的精確解，由此將直升機配平問題轉(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ù)遺傳算法，本文將目標函數(shù)式(12)作為適用度函數(shù)，設定種群數(shù)ΝP=100，交叉概率pc=0.8，變異概率pm=0.2，迭代次數(shù)N=100，數(shù)值精度為雙精度。

計算步驟如下:

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

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

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

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

3 計算結(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;初值Ⅱ均為零。計算結(jié)果如表3和表4所示。

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

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

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

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

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

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

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

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

3.2 算例驗證和分析

協(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)彎進行配平計算，并將4個操縱量轉(zhuǎn)換為操縱桿位置［9］。采用GA/LM混合算法計算的配平曲線如圖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中可以看出，本文計算的總距和腳蹬與GEN HEL［13］計算值比較吻合;由于全動平尾的安裝角是由實驗數(shù)據(jù)給出，因此本文計算的縱向操縱桿和俯仰角比GEN HEL計算值更接近飛測數(shù)據(jù);橫向駕駛桿在較大負滾轉(zhuǎn)角時與試驗數(shù)據(jù)有較大誤差，主要原因是腳蹬位置與試驗數(shù)據(jù)存在一定的靜態(tài)誤差。

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

4 結(jié)束語

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

［1］ 高正，陳仁良.直升機飛行動力學［M］.北京:科學出版社，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］.飛行力學，2010，28(1):24-28.

［6］ 李廣楊，關治，白峰杉.數(shù)值計算原理［M］.北京:清華大學出版社，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)化計算［M］.北京:電子工業(yè)出版社，2012:213-216.

［12］張呈林，郭才根.直升機總體設計［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.