李 博，田 豐，喻勇濤，郝彬彬

(1. 沈陽航空航天大學(xué)自動化學(xué)院，遼寧 沈陽 110136；2. 中國航發(fā)沈陽發(fā)動機研究所，遼寧 沈陽 110015)

1 引言

隨著當(dāng)今科技發(fā)展水平的進步，在軍事領(lǐng)域，武器裝備的科技含量越來越高，尤其體現(xiàn)在航空領(lǐng)域，如機動性強，航程遠的固定翼戰(zhàn)機；亦如可以垂直起降，可在空中懸停作業(yè)的旋翼飛機[1]。隨著各國短距起飛/垂直降落類型戰(zhàn)斗機的問世，已將上述兩種類型飛機的優(yōu)點融為一體，即能發(fā)揮出很強的機動性，又能降低對起飛降落場地的需求，使軍事上的戰(zhàn)略戰(zhàn)術(shù)更加靈活多變[2]。

戰(zhàn)斗機科技含量提高的同時，也對戰(zhàn)斗機設(shè)計和制造提出了更高的要求。綜合考慮飛機與發(fā)動機的性能、架構(gòu)等方面影響，將飛機與發(fā)動機進行一體化設(shè)計是一種有效的選擇，相比于傳統(tǒng)的飛機與發(fā)動機通常采用分離建模，再進行組合測試的設(shè)計研發(fā)體系，飛/發(fā)一體化設(shè)計能夠充分考慮飛機與發(fā)動機之間的耦合特性，將飛機與發(fā)動機之間的相互影響降至最低，以至于將飛機的整體性能達到最高[3-4]。

本文以F-35B戰(zhàn)機配備的F-135發(fā)動機為背景，根據(jù)飛機與發(fā)動機之間的耦合關(guān)系，建立基于Simulink的飛/發(fā)一體化模型。

2 發(fā)動機總體結(jié)構(gòu)

F-35B戰(zhàn)機既擁有普通固定翼戰(zhàn)機機動性強，航程遠的強大戰(zhàn)斗力，同時也具備能夠短距起飛/垂直降落的特點，這種特點主要源自于F-135發(fā)動機獨特的結(jié)構(gòu)，如圖1。

圖1 F-135發(fā)動機結(jié)構(gòu)示意圖

F-135發(fā)動機主要由升力風(fēng)扇、三軸承噴管以及滾轉(zhuǎn)噴管組成。其中三軸承噴管位于發(fā)動機后方，可向下偏轉(zhuǎn)；滾轉(zhuǎn)噴管位于發(fā)動機兩側(cè)，機翼下方，另外從主發(fā)動機低壓軸引出驅(qū)動軸，用來驅(qū)動升力風(fēng)扇。

F-35B戰(zhàn)機擁有三種飛行狀態(tài)，短距起飛/垂直降落狀態(tài)、巡航飛行狀態(tài)和懸停飛行狀態(tài)。當(dāng)飛機處于巡航飛行狀態(tài)時，升力風(fēng)扇停止工作，三軸承噴管向后偏轉(zhuǎn)，即圖1中巡航位置，提供巡航飛行過程中的全部推力；短距起飛/垂直降落狀態(tài)以及懸停飛行狀態(tài)時，三軸承噴管的三段筒體通過相互交匯處的機械結(jié)構(gòu)進行旋轉(zhuǎn)，進而改變推力矢量的方向，向下偏折90°，即圖1中垂直推力位置，與升力風(fēng)扇一起提供飛行所需的大部分升力，滾轉(zhuǎn)噴管從主發(fā)動機中引氣，為飛行提供一部分升力，并且負責(zé)調(diào)節(jié)飛機的橫向姿態(tài)，保持機體穩(wěn)定[5]。

3 飛/發(fā)一體化建模

飛/發(fā)一體化建模能使模型模擬其動態(tài)、靜態(tài)特性，并且能為飛機整體綜合控制提供便利[6]。本模型將分為兩個部分，分別是發(fā)動機動力模塊和機體模塊。

3.1 發(fā)動機動力模塊

發(fā)動機動力模塊包括三軸承噴管、升力風(fēng)扇和滾轉(zhuǎn)噴管，各部件是否開啟及其推力矢量大小和方向均可獨立控制。

發(fā)動機模塊輸入量包括：主風(fēng)扇物理轉(zhuǎn)速nzon、升力風(fēng)扇物理轉(zhuǎn)速nlift、高度H、速度V0(地面坐標系)；輸出量包括：升力風(fēng)扇推力Flift、三軸承噴管推力Fzon、左滾轉(zhuǎn)噴管推力FgzL、右滾轉(zhuǎn)噴管推力FgzR。

3.1.1 三軸承噴管模型設(shè)計

當(dāng)飛機飛行高度處于11000米以下時，可由飛機當(dāng)前高度H通過式(1)、(2)分別計算出此時的Ts0和Ps0

Ts0=288.15-0.0065H

(1)

Ps0=101325(1-H/44331)5.25588

(2)

結(jié)合當(dāng)前飛機所處高度的大氣壓力Ps0和大氣溫度Ts0，根據(jù)絕熱空氣指數(shù)k1，以及當(dāng)前飛機馬赫數(shù)M0，可通過式(4)、(5)計算出進口氣流總溫T1、進口氣流總壓P1，其中馬赫數(shù)M0由飛機速度V0通過式(3)計算得出

M0=V0/340

(3)

(4)

(5)

由主風(fēng)扇物理轉(zhuǎn)速nzon，進口氣流靜溫T1和懸停設(shè)計點進口總溫T1d，利用式(6)得到主風(fēng)扇相對換算轉(zhuǎn)速nzon，c

(6)

文獻[7]中提及進口導(dǎo)流葉片角度θ1，主風(fēng)扇進出口增壓比πzon和主風(fēng)扇換算轉(zhuǎn)速nzon，c可以通過三維插值法計算求得主風(fēng)扇的效率ηzon以及主風(fēng)扇換算流量qm，a，c。由主風(fēng)扇換算流量qm，a，c，進口氣流總溫T1以及進口氣流總壓P1通過式(7)可以計算得出主風(fēng)扇實際流量qm，a，再加上主燃油流量qm，f，即為發(fā)動機全部輸出流量，如式(8)所示，其中一部分為滾轉(zhuǎn)噴管出口流量qm，gz，余下部分為三軸承噴管實際流量qm，z。

(7)

qm，z=qm，a+qm，f-qm，gz

(8)

為了方便計算，簡化發(fā)動機內(nèi)部散失的熱量與壓力，令三軸承噴管的總壓、總溫近似等于主風(fēng)扇出口。由上文得到的進口氣流總溫T1以及進口氣流總壓P1可以通過式(9)式(10)計算得到主風(fēng)扇出口總溫T2和總壓P2，k1為絕熱空氣指數(shù)，πzon為主風(fēng)扇進出口增壓比，ηzon為主風(fēng)扇的效率。

(9)

P2=P1πzon

(10)

通過式(11)計算得出主風(fēng)扇出口靜壓Ps2，利用所得Ps2通過式(12)得到噴管出口馬赫數(shù)Mazon。馬赫數(shù)Mazon再經(jīng)過式(13)計算得到出口氣流速度Vzon，其中Cp2為比定壓熱容比。

(11)

(12)

(13)

將上述求得的數(shù)值帶入式(14)中，即可求得三軸承噴管推力Fzon，其中Azon為三軸承噴口橫截面積。

Fzon=qm，z(Vzon-V0)+(Ps2-Ps0)Azon

(14)

3.1.2 升力風(fēng)扇模型設(shè)計

升力風(fēng)扇是實現(xiàn)垂直起降的重要升力來源，在懸停或起飛降落過程中可為飛機提供最大約80KN的升力。當(dāng)機體產(chǎn)生俯仰運動時，可通過葉柵調(diào)節(jié)推力方向，平衡機體。

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

將上述求得的各數(shù)值帶入式(24)中，即可求得升力風(fēng)扇推力Flift，其中Alift為升力風(fēng)扇出口橫截面積。

(24)

3.1.3 滾轉(zhuǎn)噴管模型設(shè)計

滾轉(zhuǎn)噴管位于機體兩側(cè)，機翼正下方，當(dāng)機體產(chǎn)生滾轉(zhuǎn)運動時，從發(fā)動機外涵道引氣，為滾轉(zhuǎn)噴管提供動力，通過改變噴管角度，調(diào)節(jié)機體平衡[10]。其中涉及的中間變量包括主風(fēng)扇實際流量qm，a、出口總溫T2、出口總壓P2、靜壓Ps0。Cgz為滾轉(zhuǎn)噴管引氣比，與上文計算出的主風(fēng)扇實際流量即內(nèi)外涵道總流量qm，a通過式(25)可以計算出滾轉(zhuǎn)噴管出口流量[11]。

qm，gz=qm，aCgz

(25)

由于滾轉(zhuǎn)噴管與三軸承噴管高度幾乎相同，由式(2)可知，二者所處位置的大氣靜壓也同為Ps0，取Vgz近似等于Vzon。由于滾轉(zhuǎn)噴管需要從外涵道引氣，因此滾轉(zhuǎn)噴管的進口總壓就等于主風(fēng)扇出口總壓P2。文獻[12]提及σ15可由滾轉(zhuǎn)噴管進口總壓P2和動壓Pd2通過神經(jīng)網(wǎng)絡(luò)映射函數(shù)計算得到，具體過程不做贅述，通過式(26)即可求出滾轉(zhuǎn)噴管出口總壓P3。

P3=σ15P2

(26)

將上述求得的各個數(shù)值帶入式(27)中，即可求得滾轉(zhuǎn)噴管推力FgzL和FgzR，其中Agz滾轉(zhuǎn)噴管出口橫截面積。

Fgz=0.5qm，gz(Vgz-V0)+0.5(P3-Ps0)Agz

(27)

3.2 機體模塊

機體模塊是飛發(fā)一體化建模的重要組成部分，用于模擬在不同大小和方向的矢量推力的作用下，飛機所處的位置及姿態(tài)。通過與發(fā)動機模塊的數(shù)據(jù)交流，可以計算出機體在各個方向的力和繞各軸旋轉(zhuǎn)的力矩，并采用simulink軟件中的6DOF工具計算出機體的運動狀態(tài)與飛行姿態(tài)。

機體模塊輸入量包括：升力風(fēng)扇推力Flift、三軸承噴管推力Fzon、滾轉(zhuǎn)噴管推力FgzL、俯仰角θ、滾轉(zhuǎn)角φ；輸出量包括：速度Vex、Vey、Vez，位置Xex、Xey、Xez，機體姿態(tài)(滾轉(zhuǎn)角φ、俯仰角θ、偏航角ψ)，機體坐標系與地面坐標系變換矩陣DCMbe，機體速度Vbx、Vby、Vbz，機體角速度ωbp、ωbq、ωbr，機體角加速度dωbp/dt、dωbq/dt、dωbr/dt，機體加速度Abx、Aby、Abz；中間量包括沿三個軸的力Fx、Fy、Fz，以及繞三個軸旋轉(zhuǎn)的力矩L、M、N，其中機體坐標軸以機體質(zhì)心為原點，質(zhì)心到機體前方為X軸正方向，質(zhì)心到機體右側(cè)為Y軸正方向，過質(zhì)心且垂直于X軸與Y軸向下為Z軸正方向；地面坐標系與最初時刻機體坐標系完全重合，當(dāng)機體運動隨機體姿態(tài)改變時，地面坐標系不變。

由于整個飛機受力情況復(fù)雜，采用綜合分析將十分繁瑣，因此需要先將各部件推力分解在地面坐標系中分解，再通過simulink中的6DOF工具計算出機體的運動狀態(tài)與運行姿態(tài)。

3.2.1 三軸承噴管推力分解

三軸承噴管在運行中，既可以上下偏轉(zhuǎn)，調(diào)節(jié)俯仰姿態(tài)，也可以左右偏轉(zhuǎn)，調(diào)節(jié)偏航姿態(tài)，因此可以在地面坐標系中將三軸承噴管推力分解為Fx，zon、Fy，zon、Fz，zon三個方向的推力與繞Y、Z軸旋轉(zhuǎn)的力矩Mzon、Nzon。

當(dāng)三軸承噴管左右偏折角度為λzon，ψ時，將推力Fzon在機體坐標軸中分解為Fy，zon、Fxz，zon。

Fy，zon=Fzonsinλzon，ψ

(28)

Fxz，zon=Fzoncosλzon，ψ

(29)

(30)

(31)

(32)

(33)

設(shè)三軸承噴管噴口至質(zhì)心距離為xzon，則繞Z軸力矩Nzon與繞Y軸力矩Mzon如下：

Nzon=Fy，zonxzon

(34)

(35)

3.2.2 升力風(fēng)扇推力分解

升力風(fēng)扇推力Flift方向只受葉柵角度λlift和機體俯仰角θ影響，只在俯仰面內(nèi)，即XOZ面內(nèi)運動，因此可以將推力在地面坐標系中分解為X，Z軸方向的力Fx，lift、Fz，lift，與繞Y軸的力矩Mlift。

(36)

(37)

(38)

(39)

設(shè)升力風(fēng)扇至質(zhì)心距離為xlift，則繞Y軸力矩Mlift如下

(40)

3.2.3 滾轉(zhuǎn)噴管推力分解

由于滾轉(zhuǎn)噴管只能調(diào)節(jié)滾轉(zhuǎn)姿態(tài)，因此可以將滾轉(zhuǎn)噴管推力在地面坐標軸中分解為繞X軸的力矩Lgz和X、Y、Z軸方向的力Fx，gz、Fy，gz、Fz，gz。

假設(shè)機體此時左右滾轉(zhuǎn)噴管的偏轉(zhuǎn)角度分別為λgzL、λgzR，將FgzL與FgzR在機體坐標軸中分解。

Fy，gz=FgzLsinλgzL+FgzRsinλgzR

(41)

(42)

(43)

(44)

設(shè)左右滾轉(zhuǎn)噴管至質(zhì)心距離為xgzL和xgzR，則繞X軸力矩Lgz如下

LgzL=FgzLcosλgzLxgzL

(45)

LgzR=FgzRcosλgzRxgzR

(46)

Lgz=LgzR-LgzL

(47)

3.2.4 機體姿態(tài)解算

設(shè)飛機整體質(zhì)量為m，重力加速度g取9.8，由式(48)可知機體所受重力為Fg，且重力方向不隨機體運動而變化，始終沿地面坐標軸Z軸向下。

Fg=mg

(48)

由各部件推力分解可通過式(49)-(54)求出機體整體受力及繞各軸旋轉(zhuǎn)的力矩。

Fx=Fx，zon-Fx，lift-Fx，gz

(49)

Fy=Fy，gz-Fy，zon

(50)

Fz=Fz，g-Fz，zon-Fz，lift-Fz，gz

(51)

L=Lgz

(52)

M=Mlift+Mzon

(53)

N=Nzon

(54)

將上述Fx、Fy、Fz與L、M、N作為輸入，通過simulink中的6DOF工具，可以計算出機體的運動狀態(tài)與飛行姿態(tài)，6DOF工具如圖2。

圖2 6DOF輸入與輸出

4 仿真研究

4.1 設(shè)計點選取

根據(jù)文獻[12]，以該模型處于H=9.7km處，M0=1.2時的特性作為本模型的設(shè)計點參數(shù)，其對應(yīng)的模型仿真結(jié)果為三軸承噴管推力Fzon=78100N，升力風(fēng)扇推力Flift=81100N，左右滾轉(zhuǎn)噴管推力FgzL=FgzL=8940N，具體部件設(shè)計點參數(shù)見表1。

表1 模型設(shè)計點參數(shù)

4.2 模型仿真

為了確保本模型的正確性與有效性，將對模型進行數(shù)值仿真。將模型輸入調(diào)整為高度H=9.7km，馬赫數(shù)M0=1.2，即機體懸停在9.7千米高空時的狀態(tài)，且姿態(tài)角以及各噴管角度全部置零，令初始轉(zhuǎn)速n=7000r/min，在1000秒內(nèi)，轉(zhuǎn)速勻速提升至10000r/min，觀察三軸承噴管、升力風(fēng)扇和滾轉(zhuǎn)噴管的推力、機體模塊中各軸的分力以及繞各軸旋轉(zhuǎn)的力矩的變化，如圖3至圖7。

圖3 三軸承噴管推力仿真曲線

圖4 升力風(fēng)扇推力仿真曲線

圖5 滾轉(zhuǎn)噴管推力仿真曲線

圖6 分解推力仿真曲線

圖7 力矩仿真曲線

隨著發(fā)動機轉(zhuǎn)速從7000r/min勻速提升至10000r/min，如圖3至圖5，三個部件的推力也明顯提升。各推力經(jīng)過機體模塊分解后，如圖6，由于各部件偏轉(zhuǎn)角度全部為0°，因此升力風(fēng)扇和滾轉(zhuǎn)噴管推力垂直向下，產(chǎn)生沿Z軸負方向的升力Fz，三軸承噴管推力水平向后，產(chǎn)生沿X軸正方向的推力Fx，其中升力風(fēng)扇的升力Flift產(chǎn)生如圖7的繞Y軸旋轉(zhuǎn)的力矩M。

令馬赫數(shù)M0=1.2，轉(zhuǎn)速n=8500r/min固定不變，將高度H在1500秒內(nèi)從0km勻速提升至18km，將三軸承噴管推力Fzon仿真結(jié)果與外文文獻[12]的研究數(shù)據(jù)結(jié)果進行比較，變化趨勢基本相同，確保了動力模型的正確性與有效性，具體數(shù)據(jù)見表2。

表2 三軸承噴管推力數(shù)據(jù)對比

為驗證一體化模型的有效性，將對飛機從0m垂直爬升至1000m并保持穩(wěn)定的過程進行仿真驗證。首先將轉(zhuǎn)速設(shè)置為8500 r/min，其它各控制量、飛機位置與運行姿態(tài)等輸入量的初始值全部設(shè)置為0，并開始運行1500s，并觀察飛行高度、俯仰角、滾轉(zhuǎn)角以及偏航角來判斷飛機是否成功垂直爬升至1000m，并保持機體穩(wěn)定，如圖8至圖11。

圖8 高度仿真曲線

圖9 俯仰角仿真曲線

圖10 滾轉(zhuǎn)角仿真曲線

圖11 偏航角仿真曲線

飛機高度在250s左右達到1000m，并在750s左右趨于穩(wěn)定高度；俯仰姿態(tài)在起飛時產(chǎn)生微小偏差，僅為-0.3°，并且迅速調(diào)整進入平穩(wěn)狀態(tài)；滾轉(zhuǎn)姿態(tài)由于受左右噴管交替作用，產(chǎn)生持續(xù)等幅振蕩，但幅度極小，僅為0.75°，并不影響機體整體運行狀態(tài)；偏航角起飛階段處于減幅振蕩，偏差迅速削減，很快達到平穩(wěn)狀態(tài)，完成了使機體平穩(wěn)上升至1000m高度并保持穩(wěn)定的目標。驗證了本文飛/發(fā)一體化建模的有效性。

5 結(jié)論

本文基于F-135發(fā)動機，進行了針對于短距起飛/垂直降落的飛/發(fā)一體化建模研究。建立了由升力風(fēng)扇、三軸承噴管及滾轉(zhuǎn)噴管組成的發(fā)動機模塊模型，為飛行提供動力來源；建立了機體模塊模型，將發(fā)動機模塊產(chǎn)生的力進行分解并計算力矩，利用simulink中的6DOF工具，計算機體的運動狀態(tài)和運行姿態(tài)；對模型輸出進行仿真，將仿真結(jié)果與文獻數(shù)據(jù)進行對比，確保了模型的正確性和有效性。