北航飛行力學(xué)理論與應(yīng)用課程大作業(yè)第組

1、1飛行力學(xué)理論及應(yīng)用大作業(yè) 院（系）自動(dòng)化科學(xué)與電氣工程學(xué)院 專 業(yè)控制科學(xué)與工程 學(xué)生姓名 朱付濤 BY1403158李琛 SY1403520左京興 B1月4日飛行力學(xué)理論大作業(yè)大作業(yè)分工情況：李?。猴w機(jī)運(yùn)動(dòng)方程建模左京興：小擾動(dòng)線性化朱付濤：求解飛機(jī)時(shí)域響應(yīng)和附錄第一部分 飛機(jī)運(yùn)動(dòng)建模在平面地球假設(shè)下，飛機(jī)飛行在平面大地上空，不考慮地球曲率和地球自轉(zhuǎn)。首先列出氣動(dòng)力和氣動(dòng)力矩方程：其中，為氣動(dòng)力在飛機(jī)機(jī)體系的分量，坐標(biāo)系選為美式坐標(biāo)系。為飛機(jī)氣動(dòng)力矩在機(jī)體上的分量。為動(dòng)壓，為飛機(jī)的機(jī)翼面積，為翼展，為平均氣動(dòng)弦長(zhǎng)。1.1 飛機(jī)質(zhì)心動(dòng)力學(xué)方程由式(5.1.7)和式(

2、5.4.9)出發(fā)，推導(dǎo)飛機(jī)的質(zhì)心動(dòng)力學(xué)方程。本文推導(dǎo)了兩種情況下的質(zhì)心動(dòng)力學(xué)方程：風(fēng)軸系中的質(zhì)心動(dòng)力學(xué)方程和體軸系中的質(zhì)心動(dòng)力學(xué)方程。(1) 風(fēng)軸系中的質(zhì)心動(dòng)力學(xué)方程式（）在風(fēng)軸系中寫為式（）是一般坐標(biāo)系中慣性加速度的表示，將其寫為地軸系中的表達(dá)，此時(shí)就是質(zhì)心相對(duì)地球的速度，用來表示。假定地軸固定于慣性空間，則。原點(diǎn)的加速度就是與地球轉(zhuǎn)動(dòng)有關(guān)的向心加速度，它與g相比可以略去。其中的向心加速度項(xiàng)通常也可略去。另外令,哥氏加速度為，則有由于相對(duì)于的角速度為，利用將式轉(zhuǎn)換到風(fēng)軸系得到其中最后一項(xiàng)可以化簡(jiǎn)為：于是得到下式：其中利用矩陣表示式可得到和的表示為：將式和式代入到式，并假設(shè)大氣靜止，即W=0，

3、則可得： 將式代入到式中可得到標(biāo)量方程為： 大氣飛行的力矢量由兩部分組成：氣動(dòng)反作用力A和重力mg，即 在風(fēng)軸系中，A的分量可以寫成下面的形式：上式中，D是阻力，C是側(cè)力，L是升力。重力的表達(dá)如下式所示：因此利用風(fēng)軸系的歐拉角，可得 將式、式和式。代入到式中可得如下方程組：當(dāng)忽略地球轉(zhuǎn)動(dòng)時(shí)，和項(xiàng)為零，則有由此得到了風(fēng)軸系下飛機(jī)的質(zhì)心動(dòng)力學(xué)方程。(2) 體軸系中的質(zhì)心動(dòng)力學(xué)方程式（）在體軸系中可寫為設(shè)相對(duì)于的角速度為，則根據(jù)式及式可得：其中最后一項(xiàng)根據(jù)教材中式（）可化簡(jiǎn)為:則式可改為：式矢量的標(biāo)量分量為根據(jù)式可得將式和式代入到式中，并考慮大氣靜止的情況，即W=0，得到如下方程組：類似的可得：考慮

4、到體軸系中的氣動(dòng)力表達(dá)式為：并利用式可得：即將式和式代入到式中得到體軸系下的質(zhì)心動(dòng)力學(xué)方程為：當(dāng)忽略地球轉(zhuǎn)動(dòng)時(shí)，為零，則有：1.2飛機(jī)轉(zhuǎn)動(dòng)動(dòng)力學(xué)方程由式（）出發(fā)，推導(dǎo)飛機(jī)的轉(zhuǎn)動(dòng)動(dòng)力學(xué)方程（考慮發(fā)動(dòng)機(jī)轉(zhuǎn)子產(chǎn)生的角動(dòng)量）。僅在體軸系下考慮：，式中 將式轉(zhuǎn)換成矩陣符號(hào)為：根據(jù)式則可得h在中的分量為根據(jù)教材中式（），有，則上式可以化簡(jiǎn)為：利用可將上式中第二項(xiàng)改寫為則式可重新寫為：在以上分析的基礎(chǔ)上，設(shè)兩個(gè)體軸系和，相對(duì)于有一角速度，并考慮轉(zhuǎn)子產(chǎn)生的角動(dòng)量，則對(duì)這兩個(gè)參考系的角動(dòng)量為：令，則可將其轉(zhuǎn)換成下列形式： 當(dāng)為零，即，則得到了平均軸系下的關(guān)系式：下面再回到(5.4.12)式，它在參考系中的形式為

5、：那么轉(zhuǎn)換到體軸系中為：其中為：當(dāng)采用平均軸系時(shí)，將式代入到式中得到： 對(duì)于剛體，則上式中可以忽略，并進(jìn)行標(biāo)量展開為：將上式寫成方程組的形式為：通常是飛行器的對(duì)稱面，此時(shí),則上式可以簡(jiǎn)化為列出整個(gè)方程的標(biāo)量展開式?jīng)]有什么好處，對(duì)于忽略的且沒有轉(zhuǎn)子項(xiàng)的局限情況，即對(duì)于剛體，式（1.46）標(biāo)量方程為 1.3 飛機(jī)質(zhì)心運(yùn)動(dòng)學(xué)方程由速度矢量三角形出發(fā)，推導(dǎo)飛機(jī)的質(zhì)心運(yùn)動(dòng)學(xué)方程。飛機(jī)的質(zhì)心運(yùn)動(dòng)學(xué)方程在風(fēng)軸系和體軸系展開。(1) 風(fēng)軸系的質(zhì)心運(yùn)動(dòng)學(xué)方程由可寫出風(fēng)軸系中的方程,進(jìn)而可以得到牽連垂直坐標(biāo)系中的分量為：式中 將式展開成標(biāo)量形式為： 當(dāng)忽略地球曲率時(shí)（），即可認(rèn)為平行于，并且假設(shè)大氣相對(duì)于地球靜止

6、，即W=0，而質(zhì)心的位置坐標(biāo)可寫成下面的形式：上式便是風(fēng)軸系中飛機(jī)質(zhì)心運(yùn)動(dòng)學(xué)方程。(2) 體軸系中運(yùn)動(dòng)學(xué)方程由可寫出體軸系中的方程,進(jìn)而得到牽連垂直坐標(biāo)系中的分量：式中當(dāng)忽略地球曲率時(shí)（），即可認(rèn)為平行于，并且假設(shè)大氣相對(duì)于地球靜止，即W=0，而質(zhì)心的位置坐標(biāo)可寫成下面的形式：上式便是體軸系中的飛機(jī)質(zhì)心運(yùn)動(dòng)學(xué)方程。1.4飛機(jī)轉(zhuǎn)動(dòng)運(yùn)動(dòng)學(xué)方程由(5.2.5)式出發(fā)，推導(dǎo)飛機(jī)的轉(zhuǎn)動(dòng)運(yùn)動(dòng)學(xué)方程。(1) 體軸系轉(zhuǎn)動(dòng)運(yùn)動(dòng)學(xué)方程由式(5.2.5)可得所以將上式變換得到歐拉角的變化率的表示為：上式就是體軸系中的飛機(jī)轉(zhuǎn)動(dòng)運(yùn)動(dòng)學(xué)方程。(2) 風(fēng)軸系轉(zhuǎn)動(dòng)運(yùn)動(dòng)學(xué)方程與體軸系中的推導(dǎo)過程完全相同，最終得到的轉(zhuǎn)動(dòng)運(yùn)動(dòng)學(xué)方

7、程為：1.5六自由度全量方程組當(dāng)無風(fēng)時(shí)( W = 0 )，對(duì)于具有對(duì)稱面的剛體飛機(jī)，建立飛機(jī)的6自由度全量運(yùn)動(dòng)方程組飛機(jī)的質(zhì)心動(dòng)力學(xué)方程：對(duì)于具有對(duì)稱面的剛體飛機(jī)的轉(zhuǎn)動(dòng)動(dòng)力學(xué)方程為：飛機(jī)的質(zhì)心運(yùn)動(dòng)學(xué)方程為：飛機(jī)的轉(zhuǎn)動(dòng)運(yùn)動(dòng)學(xué)方程為：第二部分 小擾動(dòng)線性化設(shè)基準(zhǔn)運(yùn)動(dòng)為對(duì)稱定常直線水平飛行，對(duì)飛機(jī)全量運(yùn)動(dòng)方程組進(jìn)行小擾動(dòng)線性化處理。2.1飛機(jī)全量運(yùn)動(dòng)方程組飛機(jī)的全量運(yùn)動(dòng)方程組如下所示，在此基礎(chǔ)上進(jìn)行小擾動(dòng)線性化。由所給的F-16氣動(dòng)參數(shù)表，可以得到氣動(dòng)力和氣動(dòng)力矩的表達(dá)式為：氣動(dòng)力系數(shù)需要按如下方法通過插值求出：以上為飛機(jī)的全量運(yùn)動(dòng)方程組。2.2小擾動(dòng)線性化方程組設(shè)基準(zhǔn)運(yùn)動(dòng)為對(duì)稱定常直線水平飛行，對(duì)

8、剛體飛機(jī)全量運(yùn)動(dòng)方程組進(jìn)行小擾動(dòng)線性化處理。對(duì)如下縱向非線性方程進(jìn)行小擾動(dòng)線性化，有其中上述四個(gè)方程的基準(zhǔn)運(yùn)動(dòng)方程為：用式減去式并化簡(jiǎn)得到：然后將所給的力和力矩的導(dǎo)數(shù)形式代入到式中得到：然后以為狀態(tài)變量，以為控制變量，可將式線化成的形式，即以上是縱向方程組的線化形式，同理可得到橫側(cè)向方程組（其中狀態(tài)變量，控制變量）的線化形式為：2.3求氣動(dòng)導(dǎo)數(shù)在第二部分中已經(jīng)建立了分立的小擾動(dòng)線性化方程，式子中涉及到了無因次量的導(dǎo)數(shù)，需要求出這些導(dǎo)數(shù)的計(jì)算式，計(jì)算過程如下，此過程中。要建立小擾動(dòng)線性化方程，在縱向，需要求出氣動(dòng)力和氣動(dòng)力矩關(guān)于的導(dǎo)數(shù)，在橫側(cè)向，需要求出氣動(dòng)力和氣動(dòng)力矩關(guān)于的導(dǎo)數(shù)。(1) 軸向力

9、的導(dǎo)數(shù)導(dǎo)數(shù)是在基準(zhǔn)狀態(tài)下求得的，所以有軸向力的導(dǎo)數(shù)有、：其中，另外有：(2) 法向力的導(dǎo)數(shù)軸向力的導(dǎo)數(shù)有、：(3) 力矩的導(dǎo)數(shù)力矩的導(dǎo)數(shù)有、： (4) 側(cè)向力的導(dǎo)數(shù)側(cè)向力的導(dǎo)數(shù)有、： (5) 力矩的導(dǎo)數(shù)力矩的導(dǎo)數(shù)有、：(6) 力矩的導(dǎo)數(shù)力矩的導(dǎo)數(shù)有、：第三部分 飛機(jī)時(shí)域響應(yīng)求解已知某飛機(jī)數(shù)據(jù)，建立描述飛機(jī)運(yùn)動(dòng)特性的全量運(yùn)動(dòng)方程。飛機(jī)的慣性數(shù)據(jù)如下表所示：表1 飛機(jī)的慣性數(shù)據(jù)重量(lbs)轉(zhuǎn)動(dòng)慣量(slug-ft2)WIxIyIzIxz2500094965581463100982重量(kg)轉(zhuǎn)動(dòng)慣量(kg-m2)WIxIyIzIxz11339.8093飛機(jī)機(jī)翼幾何尺寸如下表所示：表2 飛機(jī)的幾何

10、尺寸翼展b(ft)機(jī)翼面積S(ft2)平均氣動(dòng)弦(ft)3030011.32翼展b(m)機(jī)翼面積S(m2)平均氣動(dòng)弦(m)由于不同的大氣密度計(jì)算公式算出的大氣密度差別很大，所以采用常用的大氣模型，1976年美國(guó)標(biāo)準(zhǔn)大氣模型USSA76，該模型以地理緯度45°3233地區(qū)海平面為基準(zhǔn)，全年實(shí)際大氣參數(shù)的統(tǒng)計(jì)平均值為標(biāo)準(zhǔn)大氣參數(shù)。大氣基準(zhǔn)值為：溫度，密度，聲速。USSA適用高度范圍為，并分成低層和高層部分。低層高度為，在這里只介紹低層大氣。低層大氣空氣分子量為常數(shù)，大氣是完全氣體，處在靜平衡狀態(tài)，其位勢(shì)高度為，重力加速度為，為海拔，與對(duì)應(yīng)的地球半徑為。于是，。 大氣參數(shù)分段計(jì)算公式：(1

11、). (2). 的大氣參數(shù)在此不列出。聲速。3.1飛機(jī)全量運(yùn)動(dòng)方程組其中，氣動(dòng)力和氣動(dòng)力矩為：氣動(dòng)力系數(shù)需要按如下方法通過插值求出：下面以定直平飛狀態(tài)作為基準(zhǔn)運(yùn)動(dòng)狀態(tài)，完成縱向力和力矩的配平，得到配平迎角、配平升降舵偏角和配平油門的大小。3.2縱向力和力矩的配平地面的重力加速度g=9.8066 m/s2、大氣密度。，定直平飛狀態(tài)下，。縱向方向配平需要用到的力和力矩方程為：縱向氣動(dòng)力和氣動(dòng)力矩的表達(dá)式為又有、和，可將上式進(jìn)一步化簡(jiǎn)為由氣動(dòng)參數(shù)表查表可知，Cz=A*(1-(beta/57.3)2)-0.19*(elevator/25.0)，在定直平飛狀態(tài)中，beta為零，于是配平方程進(jìn)一步寫為：在

12、給定飛行狀態(tài)時(shí)，上式為關(guān)于的非線性方程組。用matlab的fsolve函數(shù)可以先解出，再解出。每一種飛行狀態(tài)對(duì)應(yīng)的速度和高度如下表所示：表3 飛機(jī)的飛行狀態(tài)飛行狀態(tài)速度(m/s)高度(m)飛行狀態(tài)速度(m/s)高度(m)1500131705000260014180100003700151901000048001620010000590017210100006100018220110007110100019230110008120100020240110009130100021250120001014010002226013000111505000232701400012160500024280

13、15000配平攻角、升降舵偏角及油門的值為：圖1 對(duì)應(yīng)所有飛行狀態(tài)的配平結(jié)果我們這一組（第四組）所取飛行狀態(tài)的速度為120m/s，高度為1000m，其配平攻角、舵偏角和油門開度為：6.3296°、-0.5714°和0.169。3.3小擾動(dòng)線性化方程組(1) 縱向小擾動(dòng)方程在第二部分中已經(jīng)建立了分立的小擾動(dòng)線性化方程和氣動(dòng)導(dǎo)數(shù)的求解公式，根據(jù)3.2節(jié)的配平結(jié)果，可以求得每一個(gè)氣動(dòng)導(dǎo)數(shù)，進(jìn)一步求得縱向和橫側(cè)向的小擾動(dòng)線性化方程。其中，縱向的小擾動(dòng)線性化方程為：其中，縱向系統(tǒng)矩陣的特征值和對(duì)應(yīng)的特征向量為：和(2) 橫側(cè)向小擾動(dòng)方程橫側(cè)向的小擾動(dòng)線性化方程為： 橫側(cè)向系統(tǒng)矩陣的特

14、征值和對(duì)應(yīng)的特征向量為：和至此，求出了小擾動(dòng)線性化方程的具體形式。3.4計(jì)算時(shí)域響應(yīng)基于全量運(yùn)動(dòng)方程和小擾動(dòng)線化方程，計(jì)算飛機(jī)各狀態(tài)變量的時(shí)域響應(yīng)。飛機(jī)的狀態(tài)變量為：體軸系速度；角速度；姿態(tài)角；位移。輸入為舵偏角；發(fā)動(dòng)機(jī)油門開度。配平攻角、升降舵和油門為6.3296°、-0.5714和0.169。 (1). 狀態(tài)變量初值為零，外界輸入為零對(duì)于全量方程：初值為零：=120*cos(6.3296/57.3) 0 120*sin(6.3296/57.3) 0 0 0 0 6.3296/57.3 0 0 0 0。外界輸入為零：=0 -0.5714°0 0.169。對(duì)于小擾動(dòng)方程：初

15、值為零：和。 外界輸入為零：。(2). 縱向模態(tài)激勵(lì)下，外界輸入為零縱向特征值為和。對(duì)應(yīng)的特征向量的實(shí)部分別為和。對(duì)于全量方程：對(duì)于長(zhǎng)周期模態(tài)，模態(tài)激勵(lì)為=120*cos(6.3296/57.3) 0 120*sin(6.3296/57.3) 0 0 0 0 6.3296/57.3 0 0 0 0+，外界輸入為=0 -0.5714°0 0.169；對(duì)于短周期模態(tài)，模態(tài)激勵(lì)為=120*cos(6.3296/57.3) 0 120*sin(6.3296/57.3) 0 0 0 0 6.3296/57.3 0 0 0 0+ ，外界輸入為=0 -0.5714° 0 0.169。對(duì)于

16、小擾動(dòng)方程：對(duì)于長(zhǎng)周期模態(tài)，模態(tài)激勵(lì)為，外界輸入為；對(duì)于短周期模態(tài)，模態(tài)激勵(lì)為，外界輸入為。(3). 側(cè)向模態(tài)激勵(lì)下，外界輸入為零側(cè)向特征值為、和，現(xiàn)在取一個(gè)穩(wěn)定模態(tài)和一個(gè)不穩(wěn)定模態(tài)。對(duì)應(yīng)的特征向量分別為和。對(duì)于全量方程：對(duì)于長(zhǎng)周期模態(tài)，模態(tài)激勵(lì)為=120*cos(6.3296/57.3) 0 120*sin(6.3296/57.3) 0 0 0 0 6.3296/57.3 0 0 0 0+ ，外界輸入為=0 -0.5714°0 0.169；對(duì)于短周期模態(tài)，模態(tài)激勵(lì)為=120*cos(6.3296/57.3) 0 120*sin(6.3296/57.3) 0 0 0 0 6.3296

17、/57.3 0 0 0 0+ ，外界輸入為=0 -0.5714° 0 0.169。對(duì)于小擾動(dòng)方程：對(duì)于長(zhǎng)周期模態(tài)，模態(tài)激勵(lì)為，外界輸入為；對(duì)于短周期模態(tài)，模態(tài)激勵(lì)為，外界輸入為。(4). 升降舵單位階躍輸入（1º0.01745 rad）對(duì)于全量方程：初值為零：=120*cos(6.3296/57.3) 0 120*sin(6.3296/57.3) 0 0 0 0 6.3296/57.3 0 0 0 0。外界輸入為=0 -0.5714° 0 0.169+ 0 step 0 0。對(duì)于小擾動(dòng)方程：初值為零：和。 外界輸入為。(5). 副翼單位階躍輸入（1º0.

18、01745 rad）對(duì)于全量方程：初值為零：=120*cos(6.3296/57.3) 0 120*sin(6.3296/57.3) 0 0 0 0 6.3296/57.3 0 0 0 0。外界輸入為=0 -0.5714° 0 0.169+ step 0 0 0。對(duì)于小擾動(dòng)方程：初值為零：和。 外界輸入為。初值為零輸入為零(1) 全量方程模型全量方程模型在MATLAB中的實(shí)現(xiàn)如下圖所示：圖2 零狀態(tài)零輸入的全量方程模型全量方程在初值為零無輸入情況下的響應(yīng)結(jié)果如下圖所示：圖3 初值為零輸入為零全量方程的三軸速度圖4 初值為零輸入為零全量方程的三軸角速度圖5 初值為零輸入為零全量方程的三

19、軸姿態(tài)角(2) 小擾動(dòng)線性化方程模型小擾動(dòng)線性化模型在MATLAB中的實(shí)現(xiàn)如下圖所示：圖6 零狀態(tài)零輸入的小擾動(dòng)線性化模型小擾動(dòng)線性化模型在初值為零無輸入的情況下的響應(yīng)如下圖所示：圖7 初值為零輸入為零小擾動(dòng)方程的三軸姿速度圖8 初值為零輸入為零小擾動(dòng)方程的三軸姿角速度圖9 初值為零輸入為零小擾動(dòng)方程的姿態(tài)角結(jié)果的分析：從圖3-5和圖7-9可以看出，零狀態(tài)零輸入情況下全量方程和小擾動(dòng)方程的響應(yīng)結(jié)果是一致的。小擾動(dòng)方程的響應(yīng)一直為零，這與小擾動(dòng)方程的特性是一致的。全量方程的響應(yīng)與小擾動(dòng)方程的響應(yīng)結(jié)果一致，在開始的地方有振蕩調(diào)整過程，其縱向狀態(tài)的最終值有微小的振蕩?？v向模態(tài)激勵(lì)下全量方程在縱向模態(tài)

20、激勵(lì)下的響應(yīng)，分為2種情況：以對(duì)應(yīng)最大模特征值的特征向量的實(shí)部為狀態(tài)初始值；以對(duì)應(yīng)最小模特征值的特征向量的實(shí)部為狀態(tài)初始值。 (1) 全量方程模型全量方程模型的MATLAB實(shí)現(xiàn)如下圖所示：圖10縱向模態(tài)激勵(lì)全量方程模型以對(duì)應(yīng)最大模特征值的特征向量的實(shí)部為狀態(tài)初始值，縱向模態(tài)激勵(lì)全量模型的響應(yīng)如下圖所示：圖11 初值為零大特征值縱向模態(tài)激勵(lì)全量方程的三軸速度圖12 初值為零大特征值縱向模態(tài)激勵(lì)全量方程的三軸角速度圖13 初值為零大特征值縱向模態(tài)激勵(lì)全量方程的三軸姿態(tài)角以對(duì)應(yīng)最小模特征值的特征向量的實(shí)部為狀態(tài)初始值：圖14 初值為零小特征值縱向模態(tài)激勵(lì)全量方程的三軸速度圖15 初值為零小特征值縱向

21、模態(tài)激勵(lì)全量方程的三軸角速度圖16 初值為零小特征值縱向模態(tài)激勵(lì)全量方程的三軸姿態(tài)角(2) 小擾動(dòng)線性化方程模型小擾動(dòng)線性化模型的MATLAB實(shí)現(xiàn)如下圖所示：圖17縱向模態(tài)激勵(lì)小擾動(dòng)方程模型小擾動(dòng)方程的響應(yīng)，大特征值對(duì)應(yīng)的模態(tài)激勵(lì)下的響應(yīng)： 圖18 初值為零大特征值縱向模態(tài)激勵(lì)小擾動(dòng)方程的三軸速度圖19 初值為零大特征值縱向模態(tài)激勵(lì)小擾動(dòng)方程的三軸角速度圖20 初值為零大特征值縱向模態(tài)激勵(lì)小擾動(dòng)方程的姿態(tài)角小擾動(dòng)方程的響應(yīng)，小特征值對(duì)應(yīng)的模態(tài)激勵(lì)下的響應(yīng)：圖21 初值為零小特征值縱向模態(tài)激勵(lì)小擾動(dòng)方程的三軸速度圖22 初值為零小特征值縱向模態(tài)激勵(lì)小擾動(dòng)方程的三軸角速度圖23 初值為零小特征值縱

22、向模態(tài)激勵(lì)小擾動(dòng)方程的姿態(tài)角分析：對(duì)比圖11-13和圖14-16可知，全量方程模型在兩種縱向激勵(lì)作用下的響應(yīng)差別不大。對(duì)比圖18-20和圖21-23可知，小擾動(dòng)方程模型在長(zhǎng)周期模態(tài)激勵(lì)和短周期激勵(lì)作用下的響應(yīng)差別明顯，短周期響應(yīng)在4秒時(shí)趨于穩(wěn)定，而長(zhǎng)周期響應(yīng)需要近1000秒時(shí)間才能趨于穩(wěn)定。從圖11-23可以看出，小擾動(dòng)線性化方程模型便于分析飛機(jī)的運(yùn)動(dòng)特性，能夠準(zhǔn)確分析飛機(jī)的長(zhǎng)周期模態(tài)和短周期模態(tài)。而全量方程模型雖然能夠體現(xiàn)出一定的短周期特性，但整體上分析不出飛機(jī)的運(yùn)動(dòng)特性。橫側(cè)向模態(tài)激勵(lì)下飛機(jī)橫側(cè)向特征值為、和，現(xiàn)在取一個(gè)穩(wěn)定模態(tài)和一個(gè)不穩(wěn)定模態(tài)來求取響應(yīng)。對(duì)應(yīng)的特征向量分別為和。(1) 全

23、量方程模型全量方程的MATLAB實(shí)現(xiàn)如下圖所示：圖24橫側(cè)向模態(tài)激勵(lì)全量方程模型以大特征值對(duì)應(yīng)的特征向量為初始狀態(tài)求得響應(yīng)如下圖所示：圖25 初值為零大特征值橫側(cè)向模態(tài)激勵(lì)全量方程的速度圖26 初值為零大特征值橫側(cè)向模態(tài)激勵(lì)全量方程的角速度圖27 初值為零大特征值橫側(cè)向模態(tài)激勵(lì)全量方程的姿態(tài)角以小特征值對(duì)應(yīng)的特征向量為初始狀態(tài)求得響應(yīng)如下圖所示：圖28 初值為零小特征值橫側(cè)向模態(tài)激勵(lì)全量方程的速度圖29 初值為零小特征值橫側(cè)向模態(tài)激勵(lì)全量方程的角速度圖30 初值為零小特征值橫側(cè)向模態(tài)激勵(lì)全量方程的姿態(tài)角(2) 小擾動(dòng)線性化方程模型小擾動(dòng)線性化模型在MATLAB中的實(shí)現(xiàn)如下圖所示：圖31橫側(cè)向模

24、態(tài)激勵(lì)小擾動(dòng)方程模型小擾動(dòng)線性化模型在橫側(cè)向模態(tài)激勵(lì)下的響應(yīng)，大實(shí)根的特征向量為狀態(tài)初始值：圖32 初值為零大特征值橫側(cè)向模態(tài)激勵(lì)小擾動(dòng)方程的速度圖33 初值為零大特征值橫側(cè)向模態(tài)激勵(lì)小擾動(dòng)方程的角速度圖34 初值為零大特征值橫側(cè)向模態(tài)激勵(lì)小擾動(dòng)方程的姿態(tài)角小實(shí)根的特征向量為狀態(tài)初始值：圖35 初值為零小特征值橫側(cè)向模態(tài)激勵(lì)小擾動(dòng)方程的速度圖36 初值為零小特征值橫側(cè)向模態(tài)激勵(lì)小擾動(dòng)方程的角速度圖37 初值為零小特征值橫側(cè)向模態(tài)激勵(lì)小擾動(dòng)方程的姿態(tài)角分析：橫側(cè)向長(zhǎng)周期模態(tài)對(duì)應(yīng)的特征值為穩(wěn)定特征點(diǎn)，而短周期模態(tài)對(duì)應(yīng)的特征值是不穩(wěn)定的，這從圖25-30和圖32-37也能看出，無論是全量方程模型還是

25、小擾動(dòng)線性化模型，短周期模態(tài)激勵(lì)下的狀態(tài)響應(yīng)都是不穩(wěn)定的。對(duì)于長(zhǎng)周期模態(tài)激勵(lì)下的響應(yīng)，小擾動(dòng)線性化模型的結(jié)果能夠很快的收斂；而全量方程模型的結(jié)果收斂速度較慢。 升降舵單位階躍輸入全量方程模型在MATLAB中的實(shí)現(xiàn)如下圖所示。其中，配平升降舵偏角為-0.5714，配平油門開度為0.169，升降舵階躍輸入的幅值為0.05。(1) 全量方程模型圖38 零狀態(tài)升降舵階躍輸入全量方程模型全量模型在升降舵階躍輸入下的響應(yīng)結(jié)果如下圖所示：圖39 零狀態(tài)升降舵階躍輸入全量方程的速度圖40 零狀態(tài)升降舵階躍輸入全量方程的角速度圖41 零狀態(tài)升降舵階躍輸入全量方程的姿態(tài)角(2) 小擾動(dòng)線性化方程模型小擾動(dòng)線性化方

26、程模型在MATLAB中的實(shí)現(xiàn)如下圖所示：圖42 零狀態(tài)升降舵階躍輸入小擾動(dòng)方程模型小擾動(dòng)線性化模型在升降舵單位階躍輸入下的響應(yīng)結(jié)果如下圖所示：圖43 零狀態(tài)升降舵階躍輸入小擾動(dòng)方程的速度圖44 零狀態(tài)升降舵階躍輸入小擾動(dòng)方程的角速度圖45 零狀態(tài)升降舵階躍輸入小擾動(dòng)方程的姿態(tài)角分析：從圖39-45可以看出，飛機(jī)狀態(tài)量最終趨于穩(wěn)定，說明飛機(jī)縱向動(dòng)態(tài)穩(wěn)定；對(duì)比圖39-41和圖43-45可以看出全量方程模型和小擾動(dòng)模型的動(dòng)態(tài)調(diào)整時(shí)間都較長(zhǎng)。飛機(jī)的縱向動(dòng)態(tài)穩(wěn)定和2對(duì)復(fù)根的實(shí)部都小于0相對(duì)應(yīng)。副翼單位階躍輸入副翼階躍輸入下全量方程在MATLAB中的實(shí)現(xiàn)如下圖所示，副翼階躍輸入幅值為0.01。(1) 全量

27、方程模型圖46 零狀態(tài)副翼階躍輸入全量方程模型全量模型在副翼階躍輸入下的響應(yīng)如下圖所示：圖47 零狀態(tài)副翼階躍輸入全量方程的速度圖48 零狀態(tài)副翼階躍輸入全量方程的角速度圖49 零狀態(tài)副翼階躍輸入全量方程的姿態(tài)角(2) 小擾動(dòng)線性化方程模型小擾動(dòng)線性化方程在MATLAB中的實(shí)現(xiàn)如下圖所示：圖50 零狀態(tài)副翼階躍輸入小擾動(dòng)方程模型小擾動(dòng)線性化模型在副翼單位階躍輸入下的響應(yīng)如下圖所示：圖51 零狀態(tài)副翼階躍輸入小擾動(dòng)方程的速度圖52 零狀態(tài)副翼階躍輸入全量方程的角速度圖53 零狀態(tài)副翼階躍輸入全量方程的姿態(tài)角分析：從圖47-49和圖51-53可以看出，全量方程模型和小擾動(dòng)線性化模型的副翼階躍響應(yīng)都

28、是不穩(wěn)定的，這與飛機(jī)橫側(cè)向有正實(shí)根是相對(duì)應(yīng)的。附錄5.1 配平程序5.1.1 配平主函數(shù)clear all;close all;clc%-所給數(shù)據(jù)及矩陣將需要插值求出的氣動(dòng)數(shù)據(jù)設(shè)為全局變量-%global given_Alpha given_Elevator given_Cz given_Cx given_Cm given_Cl given_Cn given_Cl_delta_a given_Cl_delta_rglobal given_Cn_delta_a given_Cn_delta_r given_throtte_000 given_throtte_077 given_throtte_1

29、00global given_Cxq given_Cyr given_Cyp given_Czq given_Clr given_Clp given_Cmq given_Cnr given_Cnp given_delta_beta given_beta given_height given_mahglobal m g c_ S b r2d H rou Ix Iy Iz Ixz V0 Mah0global Result_alpha Result_elevator Result_Throtte T resultalfaelevator R g0%-地球半徑，地面的重力加速度-%R = 635676

30、6;g0 = 9.80665;%-飛機(jī)的常數(shù)（國(guó)際單位制）-%m = 9307.8093;c_= 3.4503;S = 27.8709;b = 9.144;r2d=57.29577951;%-飛機(jī)的慣性數(shù)據(jù)（國(guó)際單位制）-%Ix = 12874.8446; Iy = 75673.6077; Iz = 85552.0953; Ixz = 1331.4130; Isum=Ix*Iz-Ixz2;%-F-16的氣動(dòng)參數(shù)數(shù)據(jù)-%given_Alpha = -10 -5 0 5 10 15 20 25 30 35 40 45;given_Elevator = -24 -12 012 24;given_Cz

31、 =0.7700 0.2410 -0.1000 -0.4160 -0.7310 -1.0530 -1.3660 -1.6460 -1.9170 -2.1200 -2.2480 -2.2290;given_Cx = -0.099-0.048-0.022-0.04-0.083 -0.081-0.038-0.02-0.038-0.073 -0.081-0.04-0.021-0.039-0.076 -0.063-0.021-0.004-0.025-0.072 -0.025 0.016 0.032 0.006-0.046 0.044 0.083 0.094 0.062 0.012 0.097 0.127

32、 0.128 0.087 0.024 0.113 0.137 0.13 0.085 0.025 0.145 0.162 0.154 0.1 0.043 0.167 0.177 0.161 0.11 0.053 0.174 0.179 0.155 0.104 0.047 0.166 0.167 0.138 0.091 0.04;given_Cm = 0.205 0.081-0.046-0.174-0.259 0.168 0.077-0.02-0.145-0.202 0.186 0.107-0.009-0.121-0.184 0.196 0.11-0.005-0.127-0.193 0.213 0

33、.11-0.006-0.129-0.199 0.251 0.141 0.01-0.102-0.15 0.245 0.127 0.006-0.097-0.16 0.238 0.119-0.001-0.113-0.167 0.252 0.133 0.014-0.087-0.104 0.231 0.108 0 -0.084-0.076 0.198 0.081-0.013-0.069-0.041 0.192 0.093 0.032-0.006-0.005;given_Cl = 0-0.001-0.003-0.001 0 0.07 0.009 0-0.004-0.009-0.01-0.01-0.01-0

34、.011 0-0.008-0.017-0.02-0.022-0.023-0.023 0-0.012-0.024-0.03-0.034-0.034-0.037 0-0.016-0.03-0.039-0.047-0.049-0.05 0-0.019-0.034-0.044-0.046-0.046-0.047 0-0.02-0.04-0.05-0.059-0.068-0.074 0-0.02-0.037-0.049-0.061-0.071-0.079 0-0.015-0.016-0.023-0.033-0.06-0.091 0-0.008-0.002-0.006-0.036-0.058-0.076

35、0-0.013-0.1-0.014-0.035-0.062-0.077 0-0.015-0.19-0.027-0.035-0.059-0.076;given_Cn = 00.0180.0380.0560.0640.0740.079 00.0190.0420.0570.0770.0860.09 00.0180.0420.0590.0760.0930.106 00.0190.0420.0580.0740.0890.106 00.0190.0430.0580.0730.080.096 00.0180.0390.0530.0570.0620.08 00.0130.030.0320.0290.0490.

36、068 00.0070.0170.0120.0070.0220.03 00.0040.0040.0020.0120.0280.064 0 -0.014 -0.035 -0.046 -0.034 -0.0120.015 0 -0.017 -0.047 -0.071 -0.065 -0.0020.011 0 -0.033 -0.057 -0.073 -0.041 -0.013 -0.001;given_Cl_delta_a = -0.041-0.041-0.042-0.04-0.043-0.044-0.043 -0.052-0.053-0.053-0.052-0.049-0.048-0.049 -

37、0.053-0.053-0.052-0.051-0.048-0.048-0.047 -0.056-0.053-0.051-0.052-0.049-0.047-0.045 -0.05-0.05-0.049-0.048-0.043-0.042-0.042 -0.056-0.051-0.049-0.048-0.042-0.041-0.037 -0.082-0.066-0.043-0.042-0.042-0.02-0.003 -0.059-0.043-0.035-0.037-0.036-0.028-0.013 -0.042-0.038-0.026-0.031-0.025-0.013-0.01 -0.0

38、38-0.027-0.016-0.026-0.021-0.014-0.003 -0.027-0.023-0.018-0.017-0.016-0.011-0.007 -0.017-0.016-0.014-0.012-0.011-0.01-0.008;given_Cl_delta_r = 0.0050.0070.0130.0180.0150.0210.023 0.0170.0160.0130.0150.0140.0110.01 0.0140.0140.0110.0150.0130.010.011 0.01 0.0140.0120.0140.0130.0110.011 -0.0050.0130.01

39、10.0140.0120.010.011 0.0090.0090.0090.0140.0110.0090.01 0.0190.0120.0080.0140.0110.0080.008 0.0050.0050.0050.0150.010.010.01 0 0 -0.0020.0130.0080.0060.006 -0.0050.0040.0050.0110.0080.0050.014 -0.0110.0090.0030.0060.0070 0.02 0.0080.0070.0050.0010.0030.0010 ;given_Cn_delta_a = 0.001 0.002-0.006-0.01

40、1-0.015-0.024-0.022 -0.027-0.014-0.008-0.011-0.015-0.01 0.002 -0.017-0.016-0.006-0.01-0.014-0.004-0.003 -0.013-0.016-0.006-0.009-0.012-0.002-0.005 -0.012-0.014-0.005-0.008-0.011-0.001-0.003 -0.016-0.019-0.008-0.006-0.008 0.003-0.001 0.001-0.021-0.005 0 -0.002 0.014-0.009 0.017 0.002 0.007 0.004 0.00

41、2 0.006-0.009 0.011 0.012 0.004 0.007 0.006-0.001-0.001 0.017 0.015 0.007 0.1 0.012 0.004 0.003 0.008 0.015 0.006 0.004 0.011 0.004-0.002 0.016 0.011 0.006 0.1 0.011 0.006 0.001;given_Cn_delta_r = -0.018-0.028-0.037-0.048-0.043-0.052-0.062 -0.052-0.051-0.041-0.045-0.044-0.034-0.034 -0.052-0.043-0.03

42、8-0.045-0.041-0.036-0.027 -0.052-0.046-0.04-0.045-0.041-0.036-0.028 -0.054-0.045-0.04-0.044-0.04-0.035-0.027 -0.049-0.049-0.038-0.045-0.038-0.028-0.027 -0.059-0.057-0.037-0.047-0.034-0.024-0.023 -0.051-0.052-0.03-0.048-0.035-0.023-0.023 -0.03-0.03-0.027-0.049-0.035-0.02-0.019 -0.037-0.033-0.024-0.04

43、5-0.029-0.016-0.009 -0.026-0.03-0.019-0.033-0.022-0.01-0.025 -0.013-0.008-0.013-0.016-0.009-0.014-0.01;given_throtte_000=1060 670 880 114015001860 635 425 690 101013301700 60 25 345 75511301525 -1020 -710-300 350910 1360 -2700 -1900-1300-247600 1100 -3600 -1400-595-342 -200 700;given_throtte_077=126

44、8091506200395024501400 1268091506313404024701400 1261093126610429026001560 1264098397090466028401660 12390101767750532032501930 1168098488050610038002310;given_throtte_100=200001500010800700040002500 214201570011225732344352600 227001686012250815450002835 242401891013760928557003215 2607021075159751

45、111568603950 2888623319183001348486425057;given_Cxq = -0.267 -0.11 0.308 1.34 2.08 2.91 2.76 2.05 1.5 1.49 1.83 1.21;given_Cyr = 0.882 0.852 0.876 0.958 0.962 0.974 0.819 0.483 0.59 1.21 -0.493 -1.04;given_Cyp = -0.108 -0.108 -0.188 0.11 0.258 0.226 -0.344 0.362 0.611 0.529 0.298 -2.27;given_Czq = -

46、8.8 -25.8 -28.9 -31.4 -31.2 -30.7 -27.7 -28.2 -29 -29.8 -38.3 -35.3;given_Clr = -0.126 -0.126 0.063 0.113 0.208 0.23 0.319 0.437 0.68 0.1 0.447 -0.33;given_Clp = -0.36 -0.359 -0.443 -0.42 -0.383 -0.375 -0.329 -0.294 -0.23 -0.21 -0.12 -0.1;given_Cmq = -7.21 -0.54 -5.23 -5.26 -6.11 -6.64 -5.69 -6 -6.2 -6.4 -6.6 -6;given_Cnr = -0.38 -0.363 -0.378 -0.386 -0.37 -0.453 -0.55 -0.582 -0.595 -0.637 -1.02 -0.804;given_Cnp = 0.061 0.052 0.052 -0.012 -0.013 -0.024 0.05 0.15 0.13 0.158 0