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

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

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

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

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

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

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

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

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

9、數(shù)，在橫側(cè)向，需要求出氣動力和氣動力矩關(guān)于的導(dǎo)數(shù)。(1) 軸向力的導(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ī)時域響應(yīng)求解已知某飛機(jī)數(shù)據(jù)，建立描述飛機(jī)運(yùn)動特性的全量運(yùn)動方程。飛機(jī)的慣性數(shù)據(jù)如下表所示：表1 飛機(jī)的慣性數(shù)據(jù)重量(lbs)轉(zhuǎn)動慣量(slug-ft2)WIxIyIzIxz2500094965581463100982重量(kg)轉(zhuǎn)動慣量(kg-m2)WIxIyIzIx

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

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

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

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

14、2) 橫側(cè)向小擾動方程橫側(cè)向的小擾動線性化方程為： 橫側(cè)向系統(tǒng)矩陣的特征值和對應(yīng)的特征向量為： 和至此，求出了小擾動線性化方程的具體形式。3.4計算時域響應(yīng)基于全量運(yùn)動方程和小擾動線化方程，計算飛機(jī)各狀態(tài)變量的時域響應(yīng)。飛機(jī)的狀態(tài)變量為：體軸系速度；角速度；姿態(tài)角；位移。輸入為舵偏角；發(fā)動機(jī)油門開度。配平攻角、升降舵和油門為6.3296°、-0.5714和0.169。 (1). 狀態(tài)變量初值為零，外界輸入為零對于全量方程：初值為零：=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。外界輸入

15、為零：=0 -0.5714°0 0.169。對于小擾動方程：初值為零：和。 外界輸入為零：。(2). 縱向模態(tài)激勵下，外界輸入為零縱向特征值為和。對應(yīng)的特征向量的實(shí)部分別為和。對于全量方程：對于長周期模態(tài)，模態(tài)激勵為=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；對于短周期模態(tài)，模態(tài)激勵為=120*cos(6.3296/57.3) 0 120*sin(6.3296/57.3) 0 0 0 0 6.3296/57.3 0 0 0

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

17、 120*sin(6.3296/57.3) 0 0 0 0 6.3296/57.3 0 0 0 0+ ，外界輸入為=0 -0.5714° 0 0.169。對于小擾動方程：對于長周期模態(tài)，模態(tài)激勵為，外界輸入為；對于短周期模態(tài)，模態(tài)激勵為，外界輸入為。(4). 升降舵單位階躍輸入（1º0.01745 rad）對于全量方程：初值為零：=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。對于小擾動方程：初

18、值為零：和。 外界輸入為。(5). 副翼單位階躍輸入（1º0.01745 rad）對于全量方程：初值為零：=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。對于小擾動方程：初值為零：和。 外界輸入為。3.4.1初值為零輸入為零(1) 全量方程模型全量方程模型在MATLAB中的實(shí)現(xiàn)如下圖所示：圖2 零狀態(tài)零輸入的全量方程模型全量方程在初值為零無輸入情況下的響應(yīng)結(jié)果如下圖所示：圖3 初值為零輸入為零全量方程的三

19、軸速度圖4 初值為零輸入為零全量方程的三軸角速度圖5 初值為零輸入為零全量方程的三軸姿態(tài)角(2) 小擾動線性化方程模型小擾動線性化模型在MATLAB中的實(shí)現(xiàn)如下圖所示：圖6 零狀態(tài)零輸入的小擾動線性化模型小擾動線性化模型在初值為零無輸入的情況下的響應(yīng)如下圖所示：圖7 初值為零輸入為零小擾動方程的三軸姿速度圖8 初值為零輸入為零小擾動方程的三軸姿角速度圖9 初值為零輸入為零小擾動方程的姿態(tài)角結(jié)果的分析：從圖3-5和圖7-9可以看出，零狀態(tài)零輸入情況下全量方程和小擾動方程的響應(yīng)結(jié)果是一致的。小擾動方程的響應(yīng)一直為零，這與小擾動方程的特性是一致的。全量方程的響應(yīng)與小擾動方程的響應(yīng)結(jié)果一致，在開始的地

20、方有振蕩調(diào)整過程，其縱向狀態(tài)的最終值有微小的振蕩。3.4.2縱向模態(tài)激勵下全量方程在縱向模態(tài)激勵下的響應(yīng)，分為2種情況：以對應(yīng)最大模特征值的特征向量的實(shí)部為狀態(tài)初始值；以對應(yīng)最小模特征值的特征向量的實(shí)部為狀態(tài)初始值。 (1) 全量方程模型全量方程模型的MATLAB實(shí)現(xiàn)如下圖所示：圖10縱向模態(tài)激勵全量方程模型以對應(yīng)最大模特征值的特征向量的實(shí)部為狀態(tài)初始值，縱向模態(tài)激勵全量模型的響應(yīng)如下圖所示：圖11 初值為零大特征值縱向模態(tài)激勵全量方程的三軸速度圖12 初值為零大特征值縱向模態(tài)激勵全量方程的三軸角速度圖13 初值為零大特征值縱向模態(tài)激勵全量方程的三軸姿態(tài)角以對應(yīng)最小模特征值的特征向量的實(shí)部為狀

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

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

23、機(jī)橫側(cè)向特征值為、和，現(xiàn)在取一個穩(wěn)定模態(tài)和一個不穩(wěn)定模態(tài)來求取響應(yīng)。對應(yīng)的特征向量分別為和。(1) 全量方程模型全量方程的MATLAB實(shí)現(xiàn)如下圖所示：圖24橫側(cè)向模態(tài)激勵全量方程模型以大特征值對應(yīng)的特征向量為初始狀態(tài)求得響應(yīng)如下圖所示：圖25 初值為零大特征值橫側(cè)向模態(tài)激勵全量方程的速度圖26 初值為零大特征值橫側(cè)向模態(tài)激勵全量方程的角速度圖27 初值為零大特征值橫側(cè)向模態(tài)激勵全量方程的姿態(tài)角以小特征值對應(yīng)的特征向量為初始狀態(tài)求得響應(yīng)如下圖所示：圖28 初值為零小特征值橫側(cè)向模態(tài)激勵全量方程的速度圖29 初值為零小特征值橫側(cè)向模態(tài)激勵全量方程的角速度圖30 初值為零小特征值橫側(cè)向模態(tài)激勵全量方

24、程的姿態(tài)角(2) 小擾動線性化方程模型小擾動線性化模型在MATLAB中的實(shí)現(xiàn)如下圖所示：圖31橫側(cè)向模態(tài)激勵小擾動方程模型小擾動線性化模型在橫側(cè)向模態(tài)激勵下的響應(yīng)，大實(shí)根的特征向量為狀態(tài)初始值：圖32 初值為零大特征值橫側(cè)向模態(tài)激勵小擾動方程的速度圖33 初值為零大特征值橫側(cè)向模態(tài)激勵小擾動方程的角速度圖34 初值為零大特征值橫側(cè)向模態(tài)激勵小擾動方程的姿態(tài)角小實(shí)根的特征向量為狀態(tài)初始值：圖35 初值為零小特征值橫側(cè)向模態(tài)激勵小擾動方程的速度圖36 初值為零小特征值橫側(cè)向模態(tài)激勵小擾動方程的角速度圖37 初值為零小特征值橫側(cè)向模態(tài)激勵小擾動方程的姿態(tài)角分析：橫側(cè)向長周期模態(tài)對應(yīng)的特征值為穩(wěn)定特征

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

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

27、3.4.5副翼單位階躍輸入副翼階躍輸入下全量方程在MATLAB中的實(shí)現(xiàn)如下圖所示，副翼階躍輸入幅值為0.01。(1) 全量方程模型圖46 零狀態(tài)副翼階躍輸入全量方程模型全量模型在副翼階躍輸入下的響應(yīng)如下圖所示：圖47 零狀態(tài)副翼階躍輸入全量方程的速度圖48 零狀態(tài)副翼階躍輸入全量方程的角速度圖49 零狀態(tài)副翼階躍輸入全量方程的姿態(tài)角(2) 小擾動線性化方程模型小擾動線性化方程在MATLAB中的實(shí)現(xiàn)如下圖所示：圖50 零狀態(tài)副翼階躍輸入小擾動方程模型小擾動線性化模型在副翼單位階躍輸入下的響應(yīng)如下圖所示：圖51 零狀態(tài)副翼階躍輸入小擾動方程的速度圖52 零狀態(tài)副翼階躍輸入全量方程的角速度圖53 零

28、狀態(tài)副翼階躍輸入全量方程的姿態(tài)角分析：從圖47-49和圖51-53可以看出，全量方程模型和小擾動線性化模型的副翼階躍響應(yīng)都是不穩(wěn)定的，這與飛機(jī)橫側(cè)向有正實(shí)根是相對應(yīng)的。附錄5.1 配平程序5.1.1 配平主函數(shù)clear all;close all;clc%-所給數(shù)據(jù)及矩陣將需要插值求出的氣動數(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_C

29、n_delta_r given_throtte_000 given_throtte_077 given_throtte_100global 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

30、_Throtte T resultalfaelevator R g0%-地球半徑，地面的重力加速度-%R = 6356766;g0 = 9.80665;%-飛機(jī)的常數(shù)（國際單位制）-%m = 11339.8093;c_= 3.4503;S = 27.8709;b = 9.144;r2d=57.29577951;%-飛機(jī)的慣性數(shù)據(jù)（國際單位制）-%Ix = 12874.8446; Iy = 75673.6077; Iz = 85552.0953; Ixz = 1331.4130; Isum=Ix*Iz-Ixz2;%-F-16的氣動參數(shù)數(shù)據(jù)-%given_Alpha = -10 -5 0 5 10

31、 15 20 25 30 35 40 45;given_Elevator = -24 -12 012 24;given_Cz =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.0

32、16 0.032 0.006-0.046 0.044 0.083 0.094 0.062 0.012 0.097 0.127 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.18

33、6 0.107-0.009-0.121-0.184 0.196 0.11-0.005-0.127-0.193 0.213 0.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

34、-0.001-0.003-0.001 0 0.07 0.009 0-0.004-0.009-0.01-0.01-0.01-0.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

35、-0.023-0.033-0.06-0.091 0-0.008-0.002-0.006-0.036-0.058-0.076 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.

36、096 00.0180.0390.0530.0570.0620.08 00.0130.030.0320.0290.0490.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-

37、0.043-0.044-0.043 -0.052-0.053-0.053-0.052-0.049-0.048-0.049 -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.

38、036-0.028-0.013 -0.042-0.038-0.026-0.031-0.025-0.013-0.01 -0.038-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.01

39、30.010.011 0.01 0.0140.0120.0140.0130.0110.011 -0.0050.0130.0110.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.

40、0050.0010.0030.0010 ;given_Cn_delta_a = 0.001 0.002-0.006-0.011-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

41、-0.021-0.005 0 -0.002 0.014-0.009 0.017 0.002 0.007 0.004 0.002 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.06

42、2 -0.052-0.051-0.041-0.045-0.044-0.034-0.034 -0.052-0.043-0.038-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 -

43、0.03-0.03-0.027-0.049-0.035-0.02-0.019 -0.037-0.033-0.024-0.045-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-

44、247600 1100 -3600 -1400-595-342 -200 700;given_throtte_077=1268091506200395024501400 1268091506313404024701400 1261093126610429026001560 1264098397090466028401660 12390101767750532032501930 1168098488050610038002310;given_throtte_100=200001500010800700040002500 214201570011225732344352600 2270016860

45、12250815450002835 242401891013760928557003215 2607021075159751111568603950 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 = -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