數(shù)值延拓算法應(yīng)用于直升機(jī)配平計(jì)算的研究

萬紹峰,曹 龍,黃俊森,曹義華

(北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院，北京 100191)

?

數(shù)值延拓算法應(yīng)用于直升機(jī)配平計(jì)算的研究

萬紹峰,曹 龍,黃俊森,曹義華

(北京航空航天大學(xué) 航空科學(xué)與工程學(xué)院，北京 100191)

直升機(jī)配平，本質(zhì)上為非線性方程組的求解。以UH-60為例，建立直升機(jī)非線性飛行動(dòng)力學(xué)模型；通過設(shè)置主旋翼軸前傾角為零簡化模型，數(shù)值延拓得到懸停狀態(tài)平衡解；以前飛速度為延拓參數(shù)，數(shù)值延拓完成定直平飛狀態(tài)的配平計(jì)算。使用MATLAB計(jì)算平臺(tái)，配平結(jié)果與參考數(shù)據(jù)吻合較好。結(jié)果表明：數(shù)值延拓方法簡單有效，其結(jié)果具有連續(xù)性和全面性，易于觀察解曲線走向。作為一類計(jì)算方法，數(shù)值延拓適用于直升機(jī)配平計(jì)算。

數(shù)值延拓；直升機(jī)；配平計(jì)算；非線性方程組

0 引言

配平計(jì)算，是一切飛行器建模分析的基礎(chǔ)，也是極其重要的一步。與固定翼飛行器相比，直升機(jī)由于存在旋翼氣動(dòng)力，物理模型通常較為復(fù)雜；表現(xiàn)在數(shù)學(xué)中，為強(qiáng)非線性方程組，各個(gè)狀態(tài)變量之間耦合強(qiáng)烈，一般不易求解。

求解非線性方程組，傳統(tǒng)的數(shù)值方法是牛頓迭代法[1]。因?yàn)榕ｎD迭代法是局部收斂的，所以對于一般的直升機(jī)配平問題[2-3]，往往給不出合適的初值估計(jì)，導(dǎo)致迭代無法收斂。

本文采用動(dòng)力系統(tǒng)理論[4]中的數(shù)值延拓[5-6]方法，以簡化模型為基礎(chǔ)，數(shù)值延拓得到直升機(jī)配平結(jié)果。

1 理論背景

1.1 動(dòng)力系統(tǒng)理論

描述一個(gè)動(dòng)力系統(tǒng)通常需要具備兩個(gè)要素：一是描述系統(tǒng)狀態(tài)的參量X∈Rn(稱為狀態(tài)點(diǎn))，另一個(gè)是給出從一狀態(tài)點(diǎn)到另一個(gè)狀態(tài)點(diǎn)的對應(yīng)法則F。具備以上兩個(gè)要素的系統(tǒng)被稱為動(dòng)力系統(tǒng)。根據(jù)對應(yīng)法則形式的不同，動(dòng)力系統(tǒng)可分為連續(xù)動(dòng)力系統(tǒng)和離散動(dòng)力系統(tǒng)。

連續(xù)動(dòng)力系統(tǒng)通常用一個(gè)常微分方程組的形式來描述：

式中，X∈Rn，為n維狀態(tài)向量；Fa(X)∈Rn，為n維向量函數(shù)；a∈Rk，為k維參數(shù)向量；X′為X的一階時(shí)間導(dǎo)數(shù)。

動(dòng)力系統(tǒng)理論中，分析系統(tǒng)的第一步，是計(jì)算動(dòng)力系統(tǒng)的平衡解。令狀態(tài)變量導(dǎo)數(shù)X′為0，即可得到平衡解，本質(zhì)上為非線性代數(shù)方程組的求解。

1.2 數(shù)值延拓算法

k=1時(shí)計(jì)算系統(tǒng)平衡解。此時(shí)，需要求解方程組：

式中，a∈R，為標(biāo)量參數(shù)。注意到方程組包含n+1個(gè)未知數(shù)，卻只有n個(gè)方程，此時(shí)它的解空間呈現(xiàn)為一條解曲線。

假設(shè)已經(jīng)獲得了解曲線上的某一點(diǎn)(X1,a1)，并能從該點(diǎn)出發(fā)，順序解出曲線上的其他點(diǎn)(X2,a2)、(X3,a3)、…，那么所有點(diǎn)的集合就是我們需要的解曲線。這一類思想誕生出的算法被人們稱為數(shù)值延拓。

大部分?jǐn)?shù)值延拓采用預(yù)報(bào)-校正。為了說明點(diǎn)集合是如何生成的，假設(shè)在解曲線上找到一個(gè)點(diǎn)(Xi,ai)，并且還有點(diǎn)(Xi,ai)處的歸一化切向量Vi。那么下一點(diǎn)(Xi+1,ai+1)的計(jì)算包含下面兩步，見圖1。

1) 預(yù)測一個(gè)新點(diǎn)；

2) 對預(yù)測點(diǎn)的校正。

圖1 預(yù)報(bào)-校正法

2 動(dòng)力學(xué)模型

本文選取UH-60直升機(jī)建模。

文獻(xiàn)[7]詳細(xì)介紹了一類適用于飛行仿真的單旋翼直升機(jī)模型，其中最重要的是旋翼模型。

旋翼作為直升機(jī)升力最重要的來源，因?yàn)閾]舞、擺振等運(yùn)動(dòng)存在，推導(dǎo)其氣動(dòng)方程前需要以下假設(shè)：

1) 槳葉剖面安裝角沿徑向線性變化；

2) 不考慮反流區(qū)和失速的影響；

3) 槳尖損失系數(shù)為1；

4)β(t)=a0(t)-a1(t)cosψ-b1(t)sinψ，揮舞角寫成福氏級數(shù)時(shí)只取到一階。

文獻(xiàn)[7]同時(shí)提及了，對于前進(jìn)比小于0.3的情況，所建模型精度滿足后續(xù)研究的要求。

為了使該模型同時(shí)適用于UH-60直升機(jī)，需要對文獻(xiàn)[7]中的模型做如下修改[8]：

1) 機(jī)身模型

UH-60機(jī)身氣動(dòng)力模型使用文獻(xiàn)[9]中的風(fēng)洞數(shù)據(jù)，通過回歸算法得到力和力矩的方程。機(jī)體力和力矩均為側(cè)滑角和仰角的函數(shù)。

2) 尾槳模型

UH-60尾槳向上斜置20°，不同于常規(guī)直升機(jī)。為了計(jì)算尾槳處的力和力矩，需要引入兩個(gè)坐標(biāo)系：斜置尾槳軸系和斜置尾槳風(fēng)軸系。

尾槳?dú)鈩?dòng)力在斜置尾槳風(fēng)軸系中的表達(dá)式和常規(guī)尾槳相似，只需替換對應(yīng)項(xiàng)就能得到尾槳處產(chǎn)生的力和力矩。

3) 水平安定面模型

水平安定面帶可變?nèi)肓鹘牵@有兩個(gè)目的：一是消除低速時(shí)安定面處下洗沖擊引起的抬頭力矩；二是優(yōu)化爬升、巡航和自旋下降時(shí)的俯仰角。

4) 控制系統(tǒng)模型

UH-60控制系統(tǒng)包括了一個(gè)俯仰角偏斜舵機(jī)，它改變了縱向周期控制和槳盤傾斜角之間的關(guān)系，目的是為了提升飛行器的縱向靜穩(wěn)定性。

將所有力和力矩投影到體軸系，得到飛行器六自由度剛體運(yùn)動(dòng)方程。

式中，CB/E為地軸系到體軸系的轉(zhuǎn)換矩陣；X、Y、Z為合力在體軸系的分量；L、M、N為合力矩在體軸系的分量；uB、vB、wB為空速在體軸系下的分量；pB、qB、rB為角速度在體軸系下的分量；IB為直升機(jī)轉(zhuǎn)動(dòng)慣量矩陣；下標(biāo)B代表體軸系。

3 直升機(jī)配平

3.1 懸停狀態(tài)配平

懸停狀態(tài)平衡解是所有平衡解中最易求也是最重要的。求解懸停平衡解時(shí)，偏航角ψ置為0；空速分量uB、vB和wB，角速度分量pB、qB和rB均為0。將上述條件帶入式(3)、(4)中，得到：

式中，X、Y、Z、L、M和N均為四個(gè)操縱量δe、δa、δc、δp的函數(shù)。

式(5)、(6)聯(lián)立，六個(gè)方程求解六個(gè)未知數(shù)，包括操縱量δe、δa、δc、δp加兩個(gè)姿態(tài)角θ、φ。求解該問題，可選擇牛頓迭代法，但迭代需要給出良好的初值估計(jì)，一般情況下很難實(shí)現(xiàn)。

本文作者在大量嘗試不同初值發(fā)現(xiàn)迭代均無法收斂后，改變思路，選擇簡化物理模型，從簡化模型的解出發(fā)，得到真實(shí)情況下的配平解。

3.1.1 簡化模型求解

一般而言，機(jī)身、水平安定面和垂直安定面貢獻(xiàn)的力和力矩值相對較小，模型簡化時(shí)可考慮置為0。另外，觀察旋翼方程發(fā)現(xiàn)，如果將旋翼軸前傾角is置為O，能大大簡化方程形式。利用上述條件后，式(5)、(6)可簡化為：

上式中，HW、YW、TMR是旋翼氣動(dòng)力分量，均為δe、δa和δc的函數(shù)；TTRCW是尾槳升力，僅為δp的函數(shù)；K是尾槳斜置角；下標(biāo)MR代表旋翼，下標(biāo)TR代表尾槳。

選取δe、δa和δc作為未知數(shù)，通過式(7)可依次得到θ、φ和δp的表達(dá)式，再帶入式(8)中。此時(shí)，問題簡化為三個(gè)方程求解三個(gè)未知數(shù)，求解難度大大降低。另外，當(dāng)未知數(shù)個(gè)數(shù)降至3個(gè)時(shí)，可直接作圖觀察解的大致范圍，給出合理的初值估計(jì)，運(yùn)用牛頓迭代輕松求解，結(jié)果見表1。

表1 is=0°狀態(tài)下的解

3.1.2 以is為參數(shù)延拓

實(shí)際模型中，旋翼軸前傾角值為3°。以表1中is=0°狀態(tài)下的解為起始點(diǎn)，選擇is增大方向?yàn)檠油胤较颍瑪?shù)值延拓，延拓得到的解曲線見圖2。

全部的延拓曲線應(yīng)該有6條，即δe、δa、δc、δp、θ、φ隨is變化的曲線，這里只給出了一條。圖2中每一個(gè)點(diǎn)都對應(yīng)一個(gè)平衡狀態(tài)，只是旋翼前傾角is不同。如果沿反方向延拓，我們就能得到旋翼后傾時(shí)對應(yīng)的平衡狀態(tài)。如果減小步長，理論上我們能精確得到任意is時(shí)對應(yīng)的平衡狀態(tài)。需要注意的是：圖2中每一個(gè)點(diǎn)都是以前一個(gè)點(diǎn)為基礎(chǔ)而得到的，整條曲線的得到只需要知道起始點(diǎn)的值(即表1的值)。

從圖2中可以看出，以is為延拓參數(shù)，總距操縱解曲線先上升后下降，直接得到is=3°狀態(tài)下的解，見表2。

圖2 以is為參數(shù)延拓

狀態(tài)變量值δe-0.0443mδa-0.0404mδc0.1449mδp-0.0742m?-0.0168radθ0.1349rad

回過來考察式(5)、(6)，表2中is=3°狀態(tài)下的解已非常接近實(shí)際懸停配平解，以其為初值估計(jì)，運(yùn)用牛頓迭代輕松求解，得到真實(shí)情況下懸停配平解，見表3。

表3 懸停配平解

3.2 前飛狀態(tài)配平

前飛配平時(shí)以前飛地速uE為延拓參數(shù)，故需補(bǔ)充相應(yīng)方程，體現(xiàn)速度在體軸系和地軸系間的轉(zhuǎn)換：

此時(shí)，以表3懸停配平解為起始點(diǎn)，選擇uE增大方向?yàn)檠油胤较颍瑪?shù)值延拓，延拓得到的解曲線見圖3-8，圖中提及的參考數(shù)據(jù)來自文獻(xiàn)[8]。

圖3 定直平飛縱向周期操縱配平曲線

圖4 定直平飛橫向周期操縱配平曲線

圖5 定直平飛總距操縱配平曲線

圖3-8中，數(shù)值延拓每進(jìn)行一次預(yù)報(bào)-校正，就得到解曲線上的一個(gè)點(diǎn)，以符號“×”顯示。前飛速度0～30m/s階段步長小，曲線顯得密集。由圖3-8可以得到各操縱量和姿態(tài)角隨前飛速度增長的變化趨勢，以及各個(gè)速度下的配平解。

另外，起始點(diǎn)也可以不選擇懸停狀態(tài)。假設(shè)我們知道30m/s時(shí)直升機(jī)的配平解，只需要沿兩個(gè)方向數(shù)值延拓，也能得到整條解曲線。

圖6 定直平飛腳蹬操縱配平曲線

圖7 定直平飛滾轉(zhuǎn)角φ配平曲線

圖8 定直平飛俯仰角θ配平曲線

4 結(jié)論

本文采用數(shù)值延拓方法，完成UH-60直升機(jī)配平。從簡化狀態(tài)解出發(fā)，數(shù)值延拓得到懸停配平解；從懸停配平解出發(fā)，數(shù)值延拓得到任意前飛速度下的配平解。計(jì)算結(jié)果與參考數(shù)據(jù)吻合較好，數(shù)值延拓適用于直升機(jī)配平計(jì)算，得到的配平解連續(xù)、全面，計(jì)算效率高。

另外，因?yàn)閿?shù)值延拓解的全面性，對于飛行器飛行包線、參數(shù)優(yōu)化等研究也有一定的借鑒意義。

[1] Achar N S, Gaonkar G H. Helicopter Trim Analysis by Shooting and Finite Methods with Optimally Damped Newton Iterations[J]. AIAA Journal, 1993, 31(2):225-234.

[2] McVicar J S G, Bradleyt R. Robust and Efficient Trimming Algorithm for Application to Advanced Mathematical Models of Rotorcraft[J]. Journal of Aircraft, 1995, 32(2):439-442.

[3] Chen F, Omri R. A Highly Robust Trim Procedure for Rotorcraft Simulations[C].AIAA Modeling and Simulation Technologies Conference and Exhibit, 2008.

[4] 張錦炎, 馮貝葉. 常微分方程幾何理論與分支問題[M]. 2版. 北京: 北京大學(xué)出版社, 2000:105-109.

Zhang Jinyan, Feng Beiye. Geometric theory of differential equations and bifurcation problems[M]. 2nded. Beijing: Peking University Press, 2000:105-109(in Chinese).

[5] Krauskopf B, Osinga H M, Galan-Vioque J. Numerical Continuation Methods for Dynamical Systems[M]. A A Dordrecht, The Netherlands: Springer, 2007:77-95.

[6] Chow C C, Villac B F. Mapping Autonomous Constellation Design Spaces Using Numerical Continuation[J]. Journal of Guidance, Control, and Dynamics, 2012, 35(5):1426-1434.

[7] Talbot P D, Tinling B E, Decker W A. A Mathematical Model of a Single Main Rotor Helicopter for Polited Simulation[R]. NASA TM-84281,1982.

[8] Kathryn B H. A Mathematical Model of the UH-60 Helicopter[R]. NASA TM-85890,1984.

[9] Holett J J. UH-60A Black Hawk Engineering Simulation Program, Volumes I[R]. NASA CR-166309,1981.

Helicopter Trim Calculation by Numerical Continuation Method

WAN Shaofeng，CAO Long，HUANG Junsen，CAO Yihua

(School of Aeronautic Science and Engineering，Beijing University of Aeronautics and Astronautics，Beijing 100191，China)

Helicopter trim calculation is essentially to solve nonlinear equations. The UH-60 helicopter nonlinear flight dynamics model was built. The model was simplified by setting the main rotor shaft forward tilt angle to zero and equilibrium solutions were found by numerical continuation method in hover state. With the forward velocity as a continuation parameter, trim calculation was accomplished by numerical continuation method at the hover in the forward. The whole process was conducted through using the MATLAB software as the computational platform and trim solutions agree well with the referenced data. The results show that, numerical continuation is simple and effective, the results are continuous and full-scale, and it’s easy to watch the trend of the solution curve. As a computational method, numerical continuation is suitable for helicopter trim calculation.

numerical continuation；helicopter；trim；nonlinear equations

2014-09-15

萬紹峰(1991-)，男，江西南昌縣人，碩士生。

1673-1220(2015)02-011-05

V212.4

A