基于Simulink的飛/發(fā)一體化建模技術研究

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

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

1 引言

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

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

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

2 發(fā)動機總體結構

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

圖1 F-135發(fā)動機結構示意圖

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

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

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

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

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

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

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

3.1.1 三軸承噴管模型設計

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

Ts0=288.15-0.0065H

(1)

Ps0=101325(1-H/44331)5.25588

(2)

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

M0=V0/340

(3)

(4)

(5)

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

(6)

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

(7)

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

(8)

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

(9)

P2=P1πzon

(10)

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

(11)

(12)

(13)

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

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

(14)

3.1.2 升力風扇模型設計

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

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

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

(24)

3.1.3 滾轉噴管模型設計

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

qm，gz=qm，aCgz

(25)

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

P3=σ15P2

(26)

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

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

(27)

3.2 機體模塊

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

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

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

3.2.1 三軸承噴管推力分解

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

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

Fy，zon=Fzonsinλzon，ψ

(28)

Fxz，zon=Fzoncosλzon，ψ

(29)

(30)

(31)

(32)

(33)

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

Nzon=Fy，zonxzon

(34)

(35)

3.2.2 升力風扇推力分解

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

(36)

(37)

(38)

(39)

設升力風扇至質心距離為xlift，則繞Y軸力矩Mlift如下

(40)

3.2.3 滾轉噴管推力分解

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

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

Fy，gz=FgzLsinλgzL+FgzRsinλgzR

(41)

(42)

(43)

(44)

設左右滾轉噴管至質心距離為xgzL和xgzR，則繞X軸力矩Lgz如下

LgzL=FgzLcosλgzLxgzL

(45)

LgzR=FgzRcosλgzRxgzR

(46)

Lgz=LgzR-LgzL

(47)

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

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

Fg=mg

(48)

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

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 設計點選取

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

表1 模型設計點參數

4.2 模型仿真

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

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

圖4 升力風扇推力仿真曲線

圖5 滾轉噴管推力仿真曲線

圖6 分解推力仿真曲線

圖7 力矩仿真曲線

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

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

表2 三軸承噴管推力數據對比

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

圖8 高度仿真曲線

圖9 俯仰角仿真曲線

圖10 滾轉角仿真曲線

圖11 偏航角仿真曲線

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

5 結論

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