[发明专利]一种针对弹性飞机阵风响应的仿真分析方法和系统

权利要求书

1.一种弹性飞机阵风响应的仿真分析方法，其特征在于，具体步骤包括：(一)对飞机本体飞行状态建模，(二)对飞机飞行环境建模，(三)对飞机所受气动力以及飞机气动力的弹性影响建模；(四)将这三个模型连成通路，建立模拟系统，最终在Simulink平台上得到弹性飞机的阵风相应；

(一)对飞机本体飞行状态建模，具体流程为：

飞机本体飞行状态为六自由度运动状态，为了对飞机本体飞行状态建模，首先定义飞机飞行的运动坐标系；飞机运动坐标系用到以下几种坐标系：地面坐标轴系(S_g)、机体坐标轴系(S_b)、气流坐标轴系(S_a)、稳定性坐标轴系(S_S)；坐标系按照欧美惯用形式定义，其中，地面坐标系的原点为地面上任选的一点(O_g)，其它坐标系的原点都取在飞机质心(O)处且坐标系与飞机固连；各坐标轴的指向定义为x_g轴在水平面内指向某一方向，z_g轴垂直于地面并指向地心，y_g轴也在水平面内并垂直于x_g轴，指向按照右手定则确定；x_b轴在飞机对称平面内并平行于飞机的设计轴线且指向机头，y_b轴垂直于飞机对称平面指向机身右方，z_b轴在飞机对称平面且与x_b轴垂直并指向机身下方；x_a轴与飞行速度V重合一致，z_a对称平面内与x_a轴垂直并指向机身下方，y_a轴垂直于Ox_az_a平面并指向机身右方；x_s轴与飞行速度V在飞机对称平面内的投影重合一致，z_s轴在飞机对称平面内与x_s轴垂直并指向机身下方，y_s轴与机体轴y_b重合一致；

几个坐标系的转换矩阵如下：

体轴到风轴：

地轴到体轴：

地轴到风轴：

其中，α迎角；β侧滑角；θ俯仰角；ψ偏航角；φ滚转角；γ航迹倾斜角；χ航迹方位角；μ航迹滚转角；

在坐标系中对飞机质量特性建模：

飞机操稳估算中的质量特性主要指转动惯量参数，也就是相对于飞机体轴系三个坐标轴的惯性矩I_x、I_y、I_z以及相对于任意两个坐标轴的惯性积I_xy、I_yz、I_xz；在研究飞机的运动特性时作出基本假设，飞机相对于对称面几何外形对称而且内部质量分布也对称，因此I_xy＝I_yz＝0；

转动惯量采用两种坐标系给出，即体轴坐标系和中心惯性主轴坐标系；在主轴坐标系中，相对于任意两个坐标轴的惯性积均为0，主轴坐标系与体轴坐标系的y轴重合，两坐标系的x轴与z轴相差一个角度，这个角度称为主轴方位角，用a表示，这个角度由公式(1)求出：

对飞机转动惯量的估算公式都是相对飞机中心惯性主轴系的，并且不考虑机落架的影响；由一组估算公式(5)和(6)组成：

公式(5)和(6)中，I_xp、I_yp、I_zp为相对飞机中心惯性主轴系的惯性矩，单位kg-m2；W为飞机重量，单位kg，b为飞机机翼展长，单位m，L_o为不含空速管的飞机总长，单位m，L_F为飞机机身总长，单位m；k_x、k_y、k_z为统计系数，对于W≤20t的飞机有k_x＝0.10～0.12，k_y＝0.18～0.19，k_z＝0.26；对于W20t的飞机有k_x＝0.13～0.15，k_y＝0.18～0.19，k_z＝0.27；公式(2)适用于小展弦比飞机的惯性矩估算，公式(3)适用于大展弦比飞机的惯性矩估算；

在确定了飞机主轴方位角的条件下，根据飞机中心惯性主轴系的惯性矩，利用公式(7)求出飞机体轴系中的惯性矩和惯性积；

建立六自由度运动的动力学方程：

把大地视为平面，把飞机看成刚体且有纵向对称平面，飞机在运动中，既有质点的相对运动，又有刚体的牵连运动，根据牛顿第二定律，得到力和力矩的两个运动方程：

其中，m为飞机质量；为速度矢量；为力矢量；为力矩矢量；为转动矢量；

经推导得到飞机的基本运动方程：

其中，v_x、v_y、v_z为x、y、z方向的速度；

X、Y、Z分别为x、y、z方向的合力；

ω_x、ω_y、ω_z为x、y、z方向的角速度；

M_x、M_y、M_z为x、y、z方向的力矩；飞机的转动惯量矩阵为：

在飞行载荷计算中，假设飞机的飞行速度V为常数，V_x≈V，阻力等于推力，得到如下五自由度方程：

其中：

根据几何关系得到如下三个辅助方程：

将气动力系数代入上述方程得到如下方程：

其中：

飞机姿态角与速度之间的关系如下：

(二)对飞机飞行环境建模

环境模型包括切变风模型、紊流模型、离散突风模型、大气参数模型和与飞行高度相关的重力加速度模型；

(1)建立风切变模型

风切变的概念：风切变指的是风矢量在竖直方向上的变化，以及飞机相对与风矢量及其变化的各种情况，按航迹，风切变分为下列4种表现形式：

顺风切变，是指水平风的变量对飞机来说是顺风；

逆风切变，是指水平风的变量对飞机来说的逆风；

侧风切变，是指飞机从一种侧风或无风状态进入另一种明显不同侧风状态，它有左右之分，使飞机发生侧滑、滚转或偏转；

垂直切变，是指飞机从无明显的升降气流区进入强烈的升降气流区的情况；

风切变的天气影响因素，包括：雷暴，锋面，超低空急流，低空逆温层；

风切变的数学描述:

以地面坐标系下定义风切变，即认为：

每个风速分量定义为：水平风速u_Wd，顺风为正；侧风分量v_Wd，向左为正；铅锤风分量w_Wd，下降风为正；

风速梯度表示为：

其中，u_Wx,v_Wx,w_Wx为纵向风、侧向风和铅锤风的水平梯度，u_Wy,v_Wy,w_Wy为纵向风、侧向风和铅锤风的侧向梯度，u_Wz,v_Wz,w_Wz为纵向风、侧向风和铅锤风的垂直梯度；

建立风切变模型：

风切变模型是将风矢量的变化加入到系统仿真的模型中；基于MilitarySpecification MIL-F-8785C中提出的算法，认为风速分量的大小仅仅是高度h的函数，纵向水平梯度u_Wx,v_Wx,w_Wx和侧向水平梯度u_Wy,v_Wy,w_Wy均为零，即：

使用如下的方程来表示切变风的大小：

其中，u_w是切变风在各个方向下的平均风速，W₂₀是离地6m(20ft)的测量风速，h是海拔高度，z0是一个常量，C种飞行阶段时，值为0.0612m；其他飞行阶段时，值为0.816m；

这样得到的切变风是在地面坐标系下的值，与六自由度非线性运动模块得到的DCM矩阵相乘，得到体轴下的风速代入仿真模型；

(2)建立紊流模型，为Dryden连续紊流模型

大气紊流是一个随机函数，它与时间和位置有着紧密的关系，这种函数关系是基于大量的测量和统计数据来进行构建的；假设:紊流是平稳的、均匀的；在这种假设条件下可以满足对飞机品质特性分析的需要；

Dryden紊流模型通过有限差分方程的滤波器，使用有界的白噪声，把紊流对系统的影响施加进仿真过程，

大气紊流模型，其对应的指数型纵向相关函数为:f(ξ)＝e^－ξ/L；基于测量和统计数据求得f(ξ)后，对得到的结果通过Fourier变换，求出Dryden模型的纵向和横向频谱函数如下：

其中，Ω为空间频率，单位为弧度/米，L_u，L_v，L_w为湍流尺度；σ_u，σ_v，σ_w为湍流强度；

使用成型滤波器进行计算，生成大气紊流数据；以伪随机信号作为输入，成型滤波器的设计必须依据于理论频谱函数，这样才能够保证所生成的大气紊流数据频谱特性的正确性；成型滤波器就是通过滤波环节将白噪声转化为期望的有色噪声；

根据公式：

其中，ω为时间频率，得出：

Φ_xx(ω)＝G(jω)²Φ_rr(ω)         (40)

把给定的输出频谱按照上式分解，得到成型滤波器的传递函数G(s)；低空扰动模型使用的是Dryden紊流扰动模型，其空间域频谱表达式，再根据方向上飞机分速度V和频率的关系：

ω＝ΩV                     (41)

得到模型的时间域频谱表达式为：

对上述各式进行分解，得到产生给定频谱Φ_x(ω)所需要的成型滤波器的传递函数G_x(s)；对于3个分量，求出所需的传递函数如下：

以上得到的就是可以进行大气紊流数值仿真计算所需要的成型滤波器所对应的传递函数的数学表达式；为能直接在程序中运算，将滤波器所对应的传递函数的形式转化为可以在程序中以数值方法进行运算的递推公式，对于低空扰动的Dryden紊流扰动模型，有如下递推公式:

(3)离散突风模型

对离散突风来讲，突风的垂直速度变化是确定的，属于确定性突风；通常离散突风模型有两种：简单点的锐边突风模型以及复杂点的1-cos离散突风模型；

使用在Military Specification MIL-F-8785C中提出的1-cos形状的阵风模型，阵风单独作用于每一个轴或者同时作用与三个轴；指定阵风强度，阵风长度，以及阵风开始时间；离散阵风模型以英尺/秒和米/秒表示风速；离散的阵风可以一次或者多次应用，以评估飞机多大风干扰的响应；离散阵风的数学表达式为：

其中，V_m为阵风强度，d_m为阵风长度，V_wind为作用在体轴轴线方向上的最终风速；

(4)大气参数模型

采用1976年美国COESA机构发布的美国标准大气参数，这一参数仅与飞行器所在高度有关，包含了大气绝对温度、大气压力、大气密度和声速，在32000米以下；

飞行高度大于84852米时，温度值是由线性外推得到，压力值是由对数外推得到，密度和声速的值则由理想气体方程得到；

(5)重力加速度模型

采用WGS-84坐标系，其原点是地球的质心，空间直角坐标系的Z轴指向BIH(1984.0)定义的地极(CTP)方向，即国际协议原点CIO，它由IAU和IUGG共同推荐；X轴指向BIH定义的零度子午面和CTP赤道的交点，Y轴和Z,X轴构成右手坐标系；WGS-84椭球采用国际大地测量与地球物理联合会第17届大会测量常数推荐值，采用的两个常用基本几何参数；

WGS-84坐标系是一个地固坐标系；有公式：

长半径：a＝6378137±2(m)；

地球引力常数和地球质量的乘积：GM＝3986005×10⁸m³s-2±0.6×108m³s^-2；

正常化二阶带谐系数：C₂₀＝-484.16685×10-6±1.3×10-9；

地球重力场二阶带球谐系数：J₂＝108263×10-8；

地球自转角速度：ω＝7292115×10-11rads^-1±0.150×10-11rads^-1；

扁率f＝0.003352810664；

在WGS84中，根据经度、纬度和高度，在地球的等位椭球上确定重力加速度；

(三)对飞机所受气动力以及飞机气动力的弹性影响建模

气动系数计算，是根据飞行器的基本构型，分别计算气流坐标系下纵向基本气动系数、横航向基本气动系数、升降舵产生的气动系数、副翼产生的气动系数、方向舵产生的气动系数、襟翼产生的气动系数、扰流板产生的气动系数、动导数产生的气动系数，然后计算各部分的气动系数之和作为飞机整体的气动系数；

计算出气流坐标系下的气动系数后，根据气流角将气动系数转化到机体坐标系下；

根据飞行器当前的动压、重心位置计算气动力和力矩沿机体轴的分量；

构建气动力系数数据库：

先建立相关参数关于气动系数及气动导数的多维数据库，利用参数辨识校核风洞实验数据，相关结果与实验值的对比如下：其中，升力计算结果在实验点满足5％的误差要求；阻力系数满足误差范围要求；除了在小攻角0°附近由于绝对值很小导致相对误差略有超出；

在仿真运行时，对所有气动参数进行插值，进而直接得到当前迎角、操纵面下机翼的气动参数，从而输入下一模块进行力和力矩的设计，最后与六自由度非线性运动模块相耦合，构建信息反馈；

分析地面效应对气动力的影响：

由镜象涡理论配合经验修正的综合方法出发,考虑到飞行力学中的力矩配平关系,形成考虑地面效应时全机的附加升力、阻力、迎角和下洗角的工程估算方法，在给定升力系数C_L和离地高度h后，由尾涡引起的迎角减少量弧度和附着涡引起的升力系数变化量，具体为：

其中，各参量的计算公式为：

其中，h为平均气动弦1/4处距地高度，b为翼展，为平均气动弦长

C_m1/4为翼身组合体1/4弦点的力矩系数，同样迎角而襟翼收起时的力矩系数C'_m1/4，N'为同样迎角下襟翼副翼收起的值，近似为：

构建气动弹性模型：

在气动弹性的计算中，舍去对计算结果影响不大的部分，首先根据当前的气动参数，计算机翼上的气动载荷，再将集中载荷根据经验公式分布在机翼表面，从而计算机翼的二维气动-结构耦合的结果，再将机翼变形量输回到气动力模块，得到当前状态下的弹性影响的气动力；

大展弦比的升力面，包括直机翼、后掠翼；当其结构的弹性特性可用工程梁理论表达时，则是一维静气动弹性问题；其定常气动力用修正片条理论，如考虑沿翼展下洗流变化，则用L.Prandtl和J.Weissinger的升力线理论；

解一维静气动弹性问题，用积分方程矩阵数值解法；根据机翼几何模型，积分方程为：

对直机翼：

对后掠翼，公式调整为：

而：

C^θθ(y,η)和C^θz(y,η)均分别与盒式梁的弯曲、扭转刚度有关；上面各式中：

q：远方来流速压δ：操纵面偏转角,rad；θ(y)：y处剖面顺航向弹性变形角，rad；C^θθ(y,η)：结构影响系数，η剖面作用单位力矩θ方向引起y剖面的弹性变形角；C^θz(y,η)：结构影响系数，η剖面作用单位力z方向引起y剖面的弹性变形角θ方向；a,c,e,d,m,n_z,g：分别指η站位的几何迎角、弦长、气动中心至刚轴的距离、翼段重心至刚轴的距离，单位展长y方向质量，z向过载、重力加速度；C_M,AC：气动力系数，下标L表示升力，M表示力矩，AC表示对气动中心，上标r表示由刚性，e表示由弹性变形引起的；

定义新的柔度影响系数：

则直机翼、后掠翼的弹性变形角用统一的矩阵方程表示：

式中：

其中为数值积分加权矩阵，如为等间距间隔，用Simpson或梯形法则；用升力线理论时，展向均非等间隔，用Multhopp求积方程，把区间(a,b)划分出Multhopp站位，则按Multhopp求积方程:

如为对称升力分布，按Multhopp法求半翼展的积分:

即可计算得到矩阵。