基于Simulink的飞/发一体化建模技术研究

李 博，田 丰，喻勇涛，郝彬彬

(1. 沈阳航空航天大学自动化学院，辽宁 沈阳 110136；2. 中国航发沈阳发动机研究所，辽宁 沈阳 110015)

1 引言

随着当今科技发展水平的进步，在军事领域，武器装备的科技含量越来越高，尤其体现在航空领域，如机动性强，航程远的固定翼战机；亦如可以垂直起降，可在空中悬停作业的旋翼飞机[1]。随着各国短距起飞/垂直降落类型战斗机的问世，已将上述两种类型飞机的优点融为一体，即能发挥出很强的机动性，又能降低对起飞降落场地的需求，使军事上的战略战术更加灵活多变[2]。

战斗机科技含量提高的同时，也对战斗机设计和制造提出了更高的要求。综合考虑飞机与发动机的性能、架构等方面影响，将飞机与发动机进行一体化设计是一种有效的选择，相比于传统的飞机与发动机通常采用分离建模，再进行组合测试的设计研发体系，飞/发一体化设计能够充分考虑飞机与发动机之间的耦合特性，将飞机与发动机之间的相互影响降至最低，以至于将飞机的整体性能达到最高[3-4]。

本文以F-35B战机配备的F-135发动机为背景，根据飞机与发动机之间的耦合关系，建立基于Simulink的飞/发一体化模型。

2 发动机总体结构

F-35B战机既拥有普通固定翼战机机动性强，航程远的强大战斗力，同时也具备能够短距起飞/垂直降落的特点，这种特点主要源自于F-135发动机独特的结构，如图1。

图1 F-135发动机结构示意图

F-135发动机主要由升力风扇、三轴承喷管以及滚转喷管组成。其中三轴承喷管位于发动机后方，可向下偏转；滚转喷管位于发动机两侧，机翼下方，另外从主发动机低压轴引出驱动轴，用来驱动升力风扇。

F-35B战机拥有三种飞行状态，短距起飞/垂直降落状态、巡航飞行状态和悬停飞行状态。当飞机处于巡航飞行状态时，升力风扇停止工作，三轴承喷管向后偏转，即图1中巡航位置，提供巡航飞行过程中的全部推力；短距起飞/垂直降落状态以及悬停飞行状态时，三轴承喷管的三段筒体通过相互交汇处的机械结构进行旋转，进而改变推力矢量的方向，向下偏折90°，即图1中垂直推力位置，与升力风扇一起提供飞行所需的大部分升力，滚转喷管从主发动机中引气，为飞行提供一部分升力，并且负责调节飞机的横向姿态，保持机体稳定[5]。

3 飞/发一体化建模

飞/发一体化建模能使模型模拟其动态、静态特性，并且能为飞机整体综合控制提供便利[6]。本模型将分为两个部分，分别是发动机动力模块和机体模块。

3.1 发动机动力模块

发动机动力模块包括三轴承喷管、升力风扇和滚转喷管，各部件是否开启及其推力矢量大小和方向均可独立控制。

发动机模块输入量包括：主风扇物理转速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，即为发动机全部输出流量，如式(8)所示，其中一部分为滚转喷管出口流量qm，gz，余下部分为三轴承喷管实际流量qm，z。

(7)

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

(8)

为了方便计算，简化发动机内部散失的热量与压力，令三轴承喷管的总压、总温近似等于主风扇出口。由上文得到的进口气流总温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的升力。当机体产生俯仰运动时，可通过叶栅调节推力方向，平衡机体。

(15)

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

将上述求得的各数值带入式(24)中，即可求得升力风扇推力Flift，其中Alift为升力风扇出口横截面积。

(24)

3.1.3 滚转喷管模型设计

滚转喷管位于机体两侧，机翼正下方，当机体产生滚转运动时，从发动机外涵道引气，为滚转喷管提供动力，通过改变喷管角度，调节机体平衡[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 机体模块

机体模块是飞发一体化建模的重要组成部分，用于模拟在不同大小和方向的矢量推力的作用下，飞机所处的位置及姿态。通过与发动机模块的数据交流，可以计算出机体在各个方向的力和绕各轴旋转的力矩，并采用simulink软件中的6DOF工具计算出机体的运动状态与飞行姿态。

机体模块输入量包括：升力风扇推力Flift、三轴承喷管推力Fzon、滚转喷管推力FgzL、俯仰角θ、滚转角φ；输出量包括：速度Vex、Vey、Vez，位置Xex、Xey、Xez，机体姿态(滚转角φ、俯仰角θ、偏航角ψ)，机体坐标系与地面坐标系变换矩阵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轴正方向；地面坐标系与最初时刻机体坐标系完全重合，当机体运动随机体姿态改变时，地面坐标系不变。

由于整个飞机受力情况复杂，采用综合分析将十分繁琐，因此需要先将各部件推力分解在地面坐标系中分解，再通过simulink中的6DOF工具计算出机体的运动状态与运行姿态。

3.2.1 三轴承喷管推力分解

三轴承喷管在运行中，既可以上下偏转，调节俯仰姿态，也可以左右偏转，调节偏航姿态，因此可以在地面坐标系中将三轴承喷管推力分解为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 滚转喷管推力分解

由于滚转喷管只能调节滚转姿态，因此可以将滚转喷管推力在地面坐标轴中分解为绕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 机体姿态解算

设飞机整体质量为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工具，可以计算出机体的运动状态与飞行姿态，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千米高空时的状态，且姿态角以及各喷管角度全部置零，令初始转速n=7000r/min，在1000秒内，转速匀速提升至10000r/min，观察三轴承喷管、升力风扇和滚转喷管的推力、机体模块中各轴的分力以及绕各轴旋转的力矩的变化，如图3至图7。

图3 三轴承喷管推力仿真曲线

图4 升力风扇推力仿真曲线

图5 滚转喷管推力仿真曲线

图6 分解推力仿真曲线

图7 力矩仿真曲线

随着发动机转速从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并保持稳定的过程进行仿真验证。首先将转速设置为8500 r/min，其它各控制量、飞机位置与运行姿态等输入量的初始值全部设置为0，并开始运行1500s，并观察飞行高度、俯仰角、滚转角以及偏航角来判断飞机是否成功垂直爬升至1000m，并保持机体稳定，如图8至图11。

图8 高度仿真曲线

图9 俯仰角仿真曲线

图10 滚转角仿真曲线

图11 偏航角仿真曲线

飞机高度在250s左右达到1000m，并在750s左右趋于稳定高度；俯仰姿态在起飞时产生微小偏差，仅为-0.3°，并且迅速调整进入平稳状态；滚转姿态由于受左右喷管交替作用，产生持续等幅振荡，但幅度极小，仅为0.75°，并不影响机体整体运行状态；偏航角起飞阶段处于减幅振荡，偏差迅速削减，很快达到平稳状态，完成了使机体平稳上升至1000m高度并保持稳定的目标。验证了本文飞/发一体化建模的有效性。

5 结论

本文基于F-135发动机，进行了针对于短距起飞/垂直降落的飞/发一体化建模研究。建立了由升力风扇、三轴承喷管及滚转喷管组成的发动机模块模型，为飞行提供动力来源；建立了机体模块模型，将发动机模块产生的力进行分解并计算力矩，利用simulink中的6DOF工具，计算机体的运动状态和运行姿态；对模型输出进行仿真，将仿真结果与文献数据进行对比，确保了模型的正确性和有效性。