CN110309579B - 一种针对弹性飞机阵风响应的仿真分析方法和系统 - Google Patents

Abstract

本发明属于飞机仿真技术领域，具体为针对弹性飞机阵风响应的仿真分析方法和系统。本发明根据弹性飞机在阵风激励下的动力学及结构相应，首先根据飞行状态，对飞机六自由度系统、大气湍流、阵风激励以及气动力数据进行建模；其次，针对在静弹性条件下，对飞机机翼在气动力载荷下的变形，以及造成的迎角变化及其对气动力系数造成的影响进行分析，并对气动力系数进行修正，使之更加符合弹性飞机飞行的真实状态，最后，将以上各个模块连接成系统，实现对弹性飞机阵风激励下的相应的分析。该系统在某通用飞机和某民用大型客机的案例中取得了很好的效果。

Description

一种针对弹性飞机阵风响应的仿真分析方法和系统

技术领域

本发明属于飞机仿真技术领域，具体涉及弹性飞机阵风响应的仿真分析方法和系统。

背景技术

准确的飞机仿真模型对飞机控制系统的设计分析有着重要的作用。飞机实际在空中的运动 很复杂，不会像教科书上为了说明问题方便而简化的运动模式那么纯粹直接，会受到气流、速 度、温度、空气压缩性等多种因素的影响。仿真分析是无破坏性的、可以个性化控制的、能够 多次重复的、不受实际工作条件限制，可以实现的功能非常齐全的一种实验手段。利用仿真实 验，可以加深研究人员对实际系统的理解，快速找出研制过程中出现的问题的原因并提供解决 方法。

仿真工具

本课题在对飞机燃油系统分析计算的基础上，将采用MATLAB/SIMULINK软件进行建模仿 真。MATLAB是矩阵实验室(Matrix Laboratory)的简称，是一种用于算法开发、数据可视化、 数据分析及数值计算的高级技术计算语言和交互式环境。MATLAB作为三大数学软件之一，强 项就是其强大的矩阵计算及仿真能力。MATLAB的应用范围非常广，包括信号和图像处理、通 信、控制系统设计、测试和测量、财务建模和分析，以及计算生物学等众多应用领域。

MATLAB软件中集成的Simulink工具包极大地方便了系统的仿真建模过程，Simulink提供 了一种图形化的仿真界面，可以进行动态系统的建模、仿真以及综合性能分析，能够对线性和 非线性系统包括离散、连续和混合系统进行分析。SIMULINK提供了一个交互式图形环境， SIMULINK是一个进行动态系统建模、仿真和综合分析的集成软件包。它可以处理的系统包括： 线性、非线性系统；离散、连续及混合系统；单任务、多任务离散事件系统。Simulink的开 放性非常高，用户可以将系统的模型通过框图和连线的形式表示出来，和使用者有很高的交互 性，允许用户根据需要修改仿真模块的参数，而且在Simulink环境下可以直接使用MATLAB软 件中的各种分析工具，为仿真模型的性能分析带来方便。软件可以对模型仿真的结果进行分析， 并可以实时地显示在可视化窗口中。Simulink内含大量的模块库，使用户能够快速便捷地建 立动态的系统模型，只需要使用鼠标简单的拖放和添加模块间的连线，就能够建立基于实际系 统的非常复杂的仿真模型。一般的仿真模型要受到线性系统的限制，而使用Simulink工具包， 能够建立更加贴近实际系统的非线性模型。一个自定义模块库，并可针对特定应用加以扩展， 可应用于控制系统设计、信号处理和通信及图像处理等众多领域。在欧美国家，很多大公司将 产品投入实际使用之前都会进行仿真试验，他们主要使用的仿真软件就是SIMULINK。MATLAB 提供了自己的编译器，全面兼容C++及Fortran两大语言。因此MATLAB成为工程师、科研工 作者最好的语言、最好的工具。

为了丰富Simulink创建仿真模型的功能，MathWorks公司陆续开发了许多有特殊功能的 模块库，使得用户可以非常方便地建立仿真模型或者完成系统分析。其中主要用于航空航天仿 真分析的模块包Aerospace Blockset提供了大量的在Simulink环境中使用的航空航天模块， 涵盖了航空航天飞行器设计中的主要应用模型。Aerospace Blockset集成了航空航天领域的 通用的标准模块，使用户可以方便地建立航空航天飞行器的仿真模型。研究者可以在仿真模型 的基础上进行全面的系统开发和测试工作，实现航空航天飞行器系统的概念设计。Aerospace Blockset通过RTW能够自动生成实时代码，开展半物理仿真研究。Aerospace Blockset模块 库提供了航空航天动力学和舵机模块、飞行器的推进系统、控制仿真系统，提供了三自由度和 六自由度运动模型库，可以对固定质量或者变质量的飞行器系统进行建模与仿真，包含了地球 的重力场、磁场、大气环境和风场的标准模型，能够帮助用户在设计中加入航空飞行模拟的关 键特性。重力场模型为1984年全球测地系统重力模型WGS1984和WMM2000地球磁场模型；大 气模型包括1976COESA和标准大气模型，以及MIL-HDBK-310和MIL-STD-210C规范定义的大气 模型；风模型包括满足MIL-F-8785C和MIL-HDBK-1797的离散阵风模型、Dryden紊流模型和 风切变模型。用户可以使用各种单位转换模块和坐标轴转换模块以满足研究的需要，同时还包 含动画模块，在可视化窗口中可以直观显示仿真模型的动画。通过在Simulink环境中建模， 使得系统结构非常清晰。与控制系统工具箱Control System Toolbox一起完成飞行控制系统 的设计，利用软件中的虚拟现实工具箱Virtual Reality Toolbox可以在三维空间中实时地显 示飞行器的飞行仿真姿态。

SIMULINK仿真步骤

在SIMULINK提供的图形用户界面上，只要进行鼠标的简单拖拉操作就可以构造出复杂的 仿真模型。它外表以方块图形呈现，且采用分层结构。从建模角度讲，这既适用于自上而下的 设计流程，又适用于自下而上的逆程设计。从分析研究角度讲，这种SIMULINK模型不仅能让 用户知道具体环节的动态细节，而且能让用户清晰地了解各器件、各子系统、各系统间的信息 交换，掌握各部分之间的交互影响。

SIMULINK的每个模块对于用户来说都是相当于一个黑匣子，用户只需要知道模块的输入 和输出及模块功能即可，而不必管模块是怎么实现的。因此用户使用SIMULINK进行建模系统 建模时的任务就是：

(a)分析计算所需建模的系统的各个物理模型的功能及原理，并找出不同元件之间的数 据关系；

(b)如果必要，根据任务要求对所需仿真系统进行适当简化，找出物理系统必须建模的 模块和它们彼此的内在逻辑联系；

(c)在SIMULINK中选择合适的模块并把他们按照自己的模型连接起来。当SIMULINK模 型库中的模型不符合自己的特定要求时，用户还可以根据特定需求建立自己定义的模块，并把 它封装成一个“黑匣子”；

(d)进行调试和仿真，若仿真结果不满足要求，用户还可以改变模块的相关参数并再次 运行仿真，直至结果符合要求。

气动弹性分析

实际情况中，飞行器的结构并不是绝对刚性的，飞行器外形结构受气动力作用会发生一定 程度的弹性变形，这种结构的变形会反作用于气动力，引起气动力的改变，从而导致结构外形 进一步的弹性变形，这样的结构变形与气动力交互作用的现象称为气动弹性现象。气动弹性效 应是飞行器研究过程中不可忽视的重要部分。对弹性体飞行器的研究涉及到空气动力作用与弹 性结构变形之间相互影响的问题，对飞行器的安全与性能起着关键的作用。随着航空航天领域 的发展，气动弹性问题对飞行器的安全和发展的威胁逐渐突显，历史上就曾发生多起事故：1903 年，Lanley的单翼机首次作有动力的飞行试验时，就由于气动弹性效应而发生了机翼的断裂， 并最终坠毁；1912年，英国的Handly Page轰炸机在作战时由于尾翼颤振而坠毁；20世纪 30年代，英国的“蛾号”与“鸽号”飞机都因发生颤振而失事；40年代，德国的V-2火箭 曾遭遇颤振而导致结构破坏。弹性效应造成的严重损失使得人们逐步认识到结构设计环节考虑 气动弹性效应的重要性。二十世纪中叶以来，研究者们将气动弹性力学定义成为一门独立的科 学分支并不断发展，其研究的结果虽然使弹性问题造成的事故次数有所减少，但依然存在一些 损失严重的事故：20世纪60年代，北京航空学院研制的无人靶机，因气动弹性造成机体失 稳而坠毁；2001年，美国的高超声速验证机X-43因控制舵面颤振而首飞失败，导致整个 Hyper-X计划遭遇重创。由此可见，对飞行器弹性问题的研究已经成为飞行器研究的重要环节。 现代飞行器的性能设计日益追求高速度、高机动性，结构设计致力于减重增效，而飞行速度不 断提高的同时也带来了复杂的气动载荷，飞行器越来越呈现出轻结构、高弹性、高频响的控制 系统特点，这些都会导致结构、空气动力和控制系统等方面的强烈非线性及其相互的强耦合作 用，因而气动弹性与控制系统相耦合的问题在设计和飞行环节中广泛存在。由于通常情况下气 动弹性现象会造成设计失效、结构破坏的有害影响，工程界常将其视为不利因素，这一矛盾突 显了气动弹性力学与控制系统耦合研究的重要性和必要性。研究这类问题，涉及到结构动力学、 非定常空气动力学和自动控制系统动力学等多学科领域的交叉，具有较强的综合性和复杂性。三者的关系可以用气动伺服弹性力学三角形来表述，如图1所示。

目前，气动弹性研究的主要方法有两种：风洞实验和数值模拟。风洞实验在揭示现象、探 索新课题方面准确可靠，因此，风洞实验在空气动力学的研究、各种飞行器的设计研制方面， 以及其他流体有关的领域中，都有广泛应用，始终是气动弹性研究的主要分支。气动弹性实验 大体上可以分为两大类：一类实验是用来为气动弹性分析取得原始数据，如飞行器结构的刚度 试验、结构地面振动试验以及地面伺服弹性试验等；另一类实验用来获取数据分析气动弹性特 性，如颤振模型风动试验、颤振飞行试验以及抖振试验等。然而需要指出的是，气动弹性问题 交叉涉及众多学科，试验研究的开展面临较大的困难。例如静气动弹性实验中，测量模型的刚 度会对气动特性产生一定的影响，通常的风洞实验采用的模型由于选材问题都具有较大的强度 和刚度，而真实的飞行器其刚度比模型低得多。因此，风动实验中必需制造出一种适用的弹性 模型，该模型能够准确模拟飞行器各部件弯曲和扭转刚度，用它进行风洞高动压模拟实验，用 以测量模型刚度对气动特性的影响。首先，弹性模型的建立本身与实际飞行器相比很难保证弹 性变形的模态完全一致；其次，当发生颤振导致模型结构破坏时，容易对风洞设备甚至实验人 员造成伤害，具有相当大的危险性。风洞实验的基础原理是相似原理，这就要求风洞流场与真 实飞行状态之间或能满足所有的相似准则，或能保证两个流场所对应的所有相似准则系数相 等。通常的风洞实验很难完全满足上述两个条件，最常见的不满足相似准则的情况是亚跨声速 风洞的雷诺数不符。以波音737为例，该飞机在巡航高度九千米以上，以每小时927公里的巡 航速度飞行，雷诺数达Re＝2.4×107，而在3米亚声速风洞实验中，风速每秒100米的情况下 雷诺数仅约为Re＝1.4×106，相差较大，实验结果则会受到影响。

与此同时，随着计算机技术和数值计算方法的发展，数值模拟方法的研究范围以及其结果 的精确度均得到了突破性的提高，因此计算气动弹性也成为了气动弹性研究中重要的发展方 向。Schuster和Liu等给出了较为通用的计算气动弹性的定义，将采用数值分析手段、包含 线性和非线性在内各种精度层次的气动弹性研究均列入了计算气动弹性的范畴^[1]。依据研究的 问题同结构、气动方面的非线性因素是否相关，计算气动弹性又可以分为经典线性气动弹性分 析与比较复杂的非线性气动弹性分析。得益于空气动力学的快速发展，计算气动弹性得到了进 一步的发展。CFD方法采用了精确形式的流动控制方程，并且数值格式、离散方法、网格策略、 并行计算技术等也逐渐成熟，目前，气动弹性计算方法正不断完善。现代飞行器设计过程中可 以用商业软件MSC.Nastran或自行编写的求解器进行求解，将气动弹性变形的影响引入到飞行 器设计的方法和过程之中。提高计算精度、扩展分析能力是气动弹性研究不懈追求的目标。因 此，面对愈加复杂的飞行器外型、更加优越的飞行性能，研究包括非定常气动力、结构动力学 和控制方程等相互耦合的多场耦合分析方法，是当前研究气动弹性问题的一项重要任务。

近年来，计算机技术和计算力学方法领域的快速发展使得采用计算流体力学进行气动弹性 动力学仿真成为可能。研究气动弹性是研究固体形变在流场作用下的各种行为以及固体形变对 流场影响这二者相互作用的一门学科。1946年，英国人Collar用一个力的三角形对气动弹 性问题作了形象的分类，如图2所示，三角形三个顶点分别代表空气动力、弹性力和惯性力， 从而区分了各学科的研究范畴；三角形中心的动气动弹性力学需要对这三种力进行耦合分析， 是飞行力学、振动和静气弹等多学科领域的交叉，具有较强的综合性和复杂性^[2]。

早期的气动弹性问题主要考虑的是机翼的颤振边界问题。第一次世界大战初期，英国的 Handley Page轰炸机发生剧烈的尾翼颤振而坠毁，这一事故使得Lanchester、Bairstow和 Fage等人开始研究飞行器的颤振问题，直至1922年，他们研究发现，质量平衡可以有效地 消除舵面颤振。1934年，Theodorsen成功地得到了低速翼型在简谐振动情况下的非定常气动 力精确解，这一理论被称为Theodorsen理论，对气动弹性力学而言具有划时代的意义^[6]。 Dugundji综合叙述了飞行器的气动弹性的研究进展，并且提出了研究气动弹性问题的四个主 要手段:l)有效的数值和理论分析；2)风洞试验；3)地面振动结构测试；4)真实条件的飞行测 试^[7]。

在经典的气动弹性问题中，早期计算气动力采用活塞、片条和准定常理论等比较简单的方 法,文献^[5]叙述了这一阶段国际上的重大研究成果。从二十世纪七十年代到九十年代，核函数 法、小扰动位势理论和偶极子格网法开始应用于亚音速、跨音速和低超音速非定常气动力的分 析中^[6]，这些计算气动力的方法随后均有效应用在了气动弹性的分析上。到了二十世纪九十年 代的中后期，随着计算机技术的快速发展，研究者们把Euler方程和N-S方程求解非定常气 动力的技术运用到气动弹性分析方面，特别是跨音速和大攻角气动弹性分析方面，取得了显著 效果，促使了气动弹性分析的飞速发展^[7][8]。

研究气动弹性问题，需要对非定常空气动力学和结构动力学(CSD)等学科进行交叉研究， 其中包括了CFD与CSD之间的耦合；CFD计算与刚体运动涉及动边界的耦合；以及CFD计算与 结构、运动三者之间的耦合等。其中，求解结构动力学系统主要有以下六种数值方法:l)振型 叠加法；2)直接积分法；3)状态空间法；4)时间有限元法；5)一阶常微分方程组初值问题的数 值解法；6)复模态法。其中直接积分法常用的有线性加速度法、Newmark方法、Wilson方法等。 在研究直接模态叠加法过程中，邢誉峰等人利用DMSM策略，分析了等截面杆、梁的碰撞问题。 文献^[9]指出:这种方法可以用来获得结构弹性碰撞问题的解析解，而且这种方法不但可以用来 分析平动结构的碰撞问题，还可以用来分析机构的各种弹性锁定问题；不但可以用来分析结构 的点碰撞问题，还可以有效分析结构的线、面接触和碰撞等问题。关于CFD/CSD的耦合方式， 流体和结构系统进行独立求解，采用合适的交替算法完成两个计算场之间的耦合，降低了各时 间步下积分的复杂度，同时也便于程序实行模块化。当前流固耦合模型主要分为紧耦合与松耦 合，采用紧耦合能够减少耦合响应过程的时间延迟，时间精度得到了提高同时增大了时间步长， 但计算量也大大增加；采用松耦合在每个物理时间步进行一次数据交换，计算量与正常的非定 常流场分析相当。此外，Farhat和Zee等人[10],叶正寅和蒋跃文等人^[11]还对二阶时间精度 松耦合格式进行了研究。气动弹性数值模拟过程中，力与位移将通过插值技术在CFD和CSD之 间进行传递。因此存在两套网格：结构网格和气动网格。非定常气动力与结构位移等参数需要 在结构/气动网格上进行交换，但是结构/气动网格两者的拓扑方式并不一致。因此必须要采用 合适的插值方法完成数据的交换，插值方法需要满足以下四点要求：1)确保插值的光顺；2) 能够精确地转换两套网格之间的信息；3)在网格较稀疏的情况下也能够完成插值；4)能够适用 于复杂外形。目前比较常用的气动/结构网格间插值方法有常体积转换法(constant volumetetrahedron,CVT)、无限平板样条法(infinite plate spline,IPS)和薄平板样条法(thinplate spline,TPS)等。

发明内容

本发明的目的在于提供一种计算量小、仿真分析性能优异的弹性飞机阵风响应的仿真分析 方法和系统。

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

(一)对飞机本体飞行状态建模

飞机的气动数据和非线性运动方程是飞机非线性动力学模型的基础。飞机的气动力和力矩 系数、静稳定性和动稳定性导数是由飞机的基本构型和飞行状态确定的，操纵导数或操纵引起 的气动系数变化与操纵面的设计和操纵输入量密切相关。在本节中，由于飞机本体飞行状态为 六自由度运动状态，为了对飞机本体飞行状态建模，因此首先要定义飞机飞行的运动坐标系， 并在机体坐标系中，对飞机质量特性分布和飞机转动惯量进行建模，最后推导飞机六自由度动 力学方程，最后在Simulink上得以实现。

飞机运动坐标系的定义

作用在飞机上的各种力和力矩的产生原因各不相同，因此，选择合适的坐标系来方便确切 的描述飞机的空间运动状态是非常必要的。本文中主要用到以下几种坐标系：地面坐标轴系 (S_g)、机体坐标轴系(S_b)、气流坐标轴系(S_a)、稳定性坐标轴系(S_S)。坐标系按照欧美 惯用形式定义，如图3所示。其中，地面坐标系的原点为地面上任选的一点(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；对于W>20t的飞机有k_x＝0.13～0.15，k_y＝0.18～0.19，k_z＝0.27。公式 (2)适用于战斗机等小展弦比飞机的惯性矩估算，公式(3)适用于中型运输机、旅客机等大 展弦比飞机的惯性矩估算。

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

在早期设计阶段，飞机的主轴方位角往往难以确定，对于轻型飞机或展弦比较小的战斗机 通常可以忽略，即将主轴方位角假设为0，此时飞机相对于中心惯性主轴系与机体轴系的转动 惯量相同，且惯性积的值为0；对于大型飞机，I_xz的作用较明显，必须考虑主轴方位角，初步 设计中通常取某一近似值。

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

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

其中，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，阻力等于推力，得到如下五 自由度方程：

其中：

/>

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

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

其中：

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

基于以上方程，在Simulink中建立飞机的非线性动力学模型，然后加入输入输出接口就 可以作为飞机的本体非线性模型在飞行仿真中使用。这样建立的模型还没有加入反馈控制环 节，只有从输入到输出的单向传递关系。而这里非线性动力学模型全部采用Simulink中已经 建立运算关系的封装模块，联合其他动力学及模块进行有反馈的非线性动力学仿真，质量模块 的计算结果也作为输入参数参与其中。

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

在建立飞机本体非线性动力学模型时，为了提高仿真分析的逼真度，在全机仿真模型中需 要建立较为真实的环境模型，环境模型包括切变风模型、紊流模型、离散突风模型、大气参数 模型和与飞行高度相关的重力加速度模型。

风切变模块和紊流模块参见图4所示。

(1)风切变模型

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

顺风切变(TAIL WIND SHEAR)指的是水平风的变量对飞机来说是顺风。顺风切变使飞 机空速减小，升力降低。飞机下沉，这是一中比较危险的风切变形式。

逆风切变(HEAD WIND SHAEAR)指的是水平风的变量对飞机来说的逆风。逆风切变使飞 机空速突然增大，升力突然增大。飞机抬升，危害相对轻些。

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

垂直切变(VERTICAL WIND SHEAR)指的是飞机从无明显的升降气流区进入强烈的升降 气流区的情况，特别的强烈的下降气流，往往有很强的猝发性，强度大，使飞机突然下沉，危 险很大。参见图5。

天气因素的影响：

雷暴

雷暴是产生强烈低空风切变的重要天气。雷暴单体下方的下曳气流在相当宽的范围内可以 造成由下击暴流和雷暴外流组成的两种不同的风切变。一种是雷暴单体中心附近下方的下击暴 流切变，他表现为范围小，生命期短，强度大的特征。另一种是下击暴流接近地面时转化为强 烈的冷空气外流，它向四外传播，沿雷暴单体前进方向伸出15—25km，并使暖湿空气入流抬 升形成阵风锋，在雷暴下方大范围内引起180度的风向变化，表现为强顺风切变和强逆风切变。 由于其中一部分强风切变区向前伸展远离雷暴主体，不易察觉，对飞行安全威胁很大。当飞机 穿越这类风切变时，先是逆风增大，而后逆风对飞机速度减小至零，并有强的下降气流，紧跟 着又是大顺风，这种风切变对飞行危害是最大的，飞机在这一过程中先是空速增大，使飞机有 利于获得高，然后空速迅速减小，使飞机迎角减小也造成升力降低，紧跟的大顺风会使风速再 次减小，如果空速减到一定程度会使飞机掉高度甚至坠地。

锋面

锋面是产生风切变最多的天气系统。锋面两侧气象要素有很大差异，锋面过渡区的垂直结 构，是产生风切变的重要条件。一般当锋面两侧的温度差大于等于5度，锋面移动速度大于等 于15m/s时，都会在锋面附近产生对飞行有影响的低空风切变，其中尤其以冷锋型低空风切变 危害较明显，但这种低空风切变一般持续时间较短。在冷风后的偏北大风区内往往也存在较严 重的低空风切变。如1973年在波士顿罗津国际机场一架坠毁的DC-10就与锋面有关，另外冷 锋附近会经常形成雷暴，也会产生强的低空风切变，因此我们不能对它忽视。

超低空急流

在1500米以下，中心位置常见于120--160米之间的超低空急流常伴随产生低空风切变， 最大风切变位于急流轴之下，我国曾观测过到的最低四十七米的超低空风切变。超低空急流出 现时，一定伴有逆温。正是逆温层阻挡了在其上的大尺度运动与地面之间的动量交换，为逆温 层以上动量的累积和储存提供了条件，同时减缓了动量的耗散，利于逆温层上急流的形成和维 持。此时天气晴朗，地面静风。

低空逆温层

晴夜，在低空易存在辐射逆温层，并常伴有低空急流，由于稳定的逆温层阻碍了上层大风 向下层的传输，使地面风很弱，在逆温层上下形成风的垂直切变。单纯辐射逆温层引起的低空 风切变的强度较雷暴和锋面引起的风切变弱的多，也比超低空急流引起的风切变稍弱些。这种 风切变存在时往往逆温层，下层风平浪静，而上层狂风大做，飞机穿越逆温层时会造成上升或 损失高度，如果逆温层很低，在飞机着陆时，飞行员处置不及时，很容易造成飞机坠地，或冲 出跑道。这种风切变存在时往往逆温层，下层风平浪静，而上层狂风大做，飞机穿越逆温度时 会造成上升或损失高度，如果逆温度很低，在飞机着陆时，飞行员处置不及时，很容易造成飞 机坠地，或冲出跑道。此外当机场周围的环境和地形比较复杂时，也会产生对飞机的起飞，着 陆有影响的低空风切变。如地形波，较大水陆界面等，一般山地高差大，水域面积大，机场附 近高大建筑物群等，均容易产生风切变，尤其在较强阵风条件下。它出现的高度一般较低，并 且通常伴随高度增加而风速减小，飞机在降过程中容易遭遇这种风切变，因在起降过程中飞机 高度低，速度小，而我们通常是逆风起落，逆风的减小使飞机空速减小，如果无法使飞机速度 加速来应付风对飞机空速的影响，轻则影响目测，重则导致飞机坠地，造成严重后果。

风切变的数学描述：

为了研究方便，这里在地面坐标系下定义风切变，即认为：

每个风速分量定义为：水平风速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。

切变风在竖直方向上的分布，参见图6所示。

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

(2)Dryden连续紊流模型

实际的大气紊流是十分复杂的物理现象。为了使飞机响应问题的研究不至于过分复杂，在 保证模型准确度的基础上可进行适当的假设。一般来说，大气紊流是一个随机函数，它与时间 和位置有着紧密的关系，这种函数关系是基于大量的测量和统计数据来进行构建的。在航空工 程应用中，可以假设:大气紊流的统计特征(即平均值和均方差，以及相关函数和频谱函数)既 不随时间而变(认为紊流是平稳的)，也不随位置而变(认为紊流是均匀的)，在这种假设条件下 可以满足对飞机品质特性分析的需要。

Dryden紊流模型通过有限差分方程的滤波器，使用有界的白噪声，把紊流对系统的影响 施加进仿真过程，这种算法首先在Military Specification MIL-F-8785C和Military Handbook MIL-HDBK-1797中提出。

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

/>

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

得到的频谱函数是基于“平稳”和“均匀”两个假设条件得到的，这就会使得所得到的这 个模型在无穷远处的渐进性质是不符合实际紊流理论的，但这并不影响实际的工程应用。该 Dryden模型的优点在于:相对来说频谱的形式较为简单，可以通过常规的数学方法进行处理， 而这一点对于紊流数值仿真是非常必要的

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

根据公式：

可以得出：

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

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

ω＝ΩV (41)

其中，ω为时间频率，得到模型的时间域频谱表达式为：

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

以上所得到的就是可以进行大气紊流数值仿真计算所需要的成型滤波器所对应的传递函 数的数学表达式。对于成型滤波器的形式还不能直接在程序中运算，必须要由滤波器所对应的 传递函数的形式得到可以在程序中以数值方法进行运算的递推公式，这也就是滤波器在程序中 的实现方式。对于低空扰动的Dryden紊流扰动模型，则有如下递推公式:

(3)离散突风模型

对离散突风来讲，突风的垂直速度变化是确定的，属于确定性突风，也就是说其突风速度 不随时间变化，是定常的。通常我们用到的离散突风模型有两种：简单点的锐边突风模型以及 复杂点的1-cos离散突风模型。在实际应用中，由于1-cos突风模型考虑了突风场的空间尺度， 在评价飞行品质和飞控系统设计时，广泛采用此模型

在Military Specification MIL-F-8785C中提出了一直种1-cos形状的阵风模型(参见 图8所示)。阵风可单独作用于每一个轴或者同时作用与三个轴。在模型中，可以指定阵风强 度(由阵风产生的风速增加)，阵风长度(以米为单位的阵风长度)，以及阵风开始时间。离 散阵风模型可以以英尺/秒和米/秒表示风速。离散的阵风可以一次或者多次应用，以评估飞机 多大风干扰的响应。离散阵风的数学表达式为：

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

(4)大气参数模型

1976年美国Committee on Extension to the Standard Atmosphere(COESA)机构发布 了美国标准大气参数，这一参数仅与飞行器所在高度有关，包含了大气绝对温度、大气压力、 大气密度和声速，在32000米以下，美国标准大气与国际民航组织定义的标准大气相同。参见 图9所示。

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

(5)重力加速度模型

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

WGS-84是修正NSWC9Z-2参考系的原点和尺度变化，并旋转其参考子午面与BIH定义的零 度子午面一致而得到的一个新参考系，WGS--84坐标系的原点在地球质心，Z轴指向BIH1984.0 定义的协定地球极(CTP)方向，X轴指向BIH1984.0的零度子午面和CTP赤道的交点，Y轴和 Z、X轴构成右手坐标系。它是一个地固(地心固连)坐标系。有公式：

长半径：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

在World Geodetic System(WGS84)中，根据经度、纬度和高度，在地球的等位椭球上 确定重力加速度。

重力加速度和飞机所在高度的关系，参见图11所示。

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

准确的气动力分布是进行准确静气动弹性特性分析和飞机型架设计的基础。飞机不同设计 阶段有不同的精度要求，需要采用不同的气动分析方法。在飞机设计初期，快速高效地形成设 计方案，确定气动主要参数；飞机详细设计阶段需要计算准确的弹性气动力，为相关专业提供 气动数据输入。因此，采用满足飞机所在设计阶段精度需求的弹性气动流场求解方法十分必要。

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

计算出气流坐标系下的气动系数后，然后根据气流角(迎角、侧滑角)将气动系数转化到 机体坐标系下。

根据飞行器当前的动压(跟飞行速度、高度相关)、重心位置计算气动力和力矩沿机体轴 的分量。

气动力系数数据库

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

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

地面效应对气动力的影响

从简化的空气动力学理论看,地面对于单独机翼或尾翼相当于一个镜象涡系〕对此简化模 型的理论分析和相应的实验研究早在40-50年代即已完成，然而作为机翼、机身和平尾的组合 体,并考虑到地面对大偏度襟翼和干扰流场及配平力矩的综合影响,情况就复杂得多,目前尚无 完整和可靠的理论计算方法，由于雷诺数及其它实验条件(如活动地板和配平舵偏角)的限制, 风洞实验的精度也受到限制。

这里由镜象涡理论配合经验修正的综合方法出发,考虑到飞行力学中的力矩配平关系,形 成了计及地效时全机的附加升力、阻力、迎角和下洗角的工程估算方法，在给定升力系数CL 和离地高度h(平均气动弦1/4处距地高度)后，由尾涡引起的迎角减少量(弧度)和附着涡引 起的升力系数变化量各为：

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

其中，h为平均气动弦1/4处距地高度，b为翼展，c为平均气动弦长。N与翼身组合体1/4弦点的力矩系数Cm1/4和在同样迎角而襟翼收起时的力矩系数C'm1/4；

N'为同样迎角下襟翼副翼收起的值，近似为：

/>

地面效应影响下的升力曲线如图14所示。

考察地面效应导致的飞机升力系数增量，与高度的关系如下图15所示。

气动弹性模块设计

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

大展弦比的升力面(直机翼，后掠翼)，当其结构的弹性特性可用工程梁理论表达时，则 是一维静气动弹性问题。其定常气动力可以用修正片条理论，如考虑沿翼展下洗流变化，则可 以用L.Prandtl和J.Weissinger的升力线理论。当然把工程问题简化为一维问题处理，无 论气动力和弹性力(变形)，均做一定简化假设，为求解更精确，大展弦比升力面也可以作为 二维问题处理。

解一维静气动弹性问题可用微分方程和积分方程，由于微分方程时其气动力只能用片条理 论。而积分方程可以用片条或升力线理论，故工程上认为用积分方程矩阵数值解法更加合理。

根据图16所示的机翼几何模型，经推导，积分方程为：

对直机翼：

对后掠翼，公式调整为：

而：

后掠翼与直机翼不同，弯曲变形可以引起顺流向气流迎角的变化，所以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求积方程，按图17把区 间(a,b)划分出Multhopp站位，则按Multhopp求积方程:

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

即可计算得到

矩阵。

根据以上公式，迎角2°，马赫数为0.8时飞机升力系数沿展向分布如图18所示。可以 看出，当飞机机翼为弹性体时，机翼沿展向发生扭转，飞机升力系数较刚性状态下略有偏大。

(四)将上述三个模型连成通路，建立模拟系统，最终得到弹性飞机的阵风相应。

基于所建立的飞机本体、环境和气动力模型，将以上模型按模拟过程中的物理规则连接起 来，首先，在每次模拟中，根据六自由度模块得到的飞机飞行状态量，经环境模块计算得到飞 机当地所处的风速、温度、重力加速度等物理变量；将这些物理变量与飞机飞行状态量结合， 得到气流坐标系中的飞机速压、迎角、侧滑角等信息；把这些信息放入气动力数据库进行插值， 得到飞机的气动系数；在地面效应和气动弹性模块的作用下，气动系数的值得到修正；最后和 重力一起，计算飞机飞行的合力以及合力矩，传回六自由度模块，计算下一时刻的飞机飞行状 态，形成通路。如图19所示为Simulink模块的连接关系，

本发明的特点在于把飞行动力学、飞机弹性变形和空气动力三者结合，分别建立模型，满 足实时仿真要求，建立统一的弹性飞机飞行动力学仿真平台，针对任意弹性飞机和任意飞行状 态，得到对应的全波长和半波长阵风相应，为飞机的柔性设计提供了统一的仿真平台，减少设 计投入，缩短设计周期。

附图说明

图1为气动伺服弹性力学三角形图示。

图2为Collar气动弹性三角形图示。

图3为飞机运动坐标系图示。

图4为风模块和紊流模块。

图5为风矢量在垂直方向上的变化。

图6为切变风在竖直方向上的分布。

图7为紊流强度的均方差和高度的关系。

图8为1-cos突风模型。

图9为美国标准大气1976。

图10为WGS-84坐标系图示。

图11为重力加速度和飞机所在高度的关系。

图12为升力与迎角的关系曲线。

图13为阻力与迎角的关系曲线。

图14为地面效应影响下的升力曲线。

图15为地面效应导致升力系数增量与高度关系。

图16为机翼几何模型。

图17为把区间(a,b)划分出Multhopp站位。

图18为弹性飞机与刚性飞机升力系数沿展向分布对比。

图19为搭建的整体仿真系统

图20为某小型通用飞机几何模型。

图21为某小型通用飞机中作用在副翼重心上的载荷随时间变化曲线。

图22为某小型通用飞机中作用在飞机重心上的载荷随时间变化曲线。

图23为某大型客机几何模型。

图24为某大型客机中翼梢位移随时间变化曲线。

图25为某大型客机迎角随时间变化曲线。

具体实施方式

下面对该阵风相应分析方法进行案例验证。

(一)某小型通用飞机弹性阵风相应的模拟验证

在某小型通用飞机的弹性阵风相应的分析中，飞机几何外形如图20所示，机翼为平直下 单翼，由于飞机较小，机翼、尾翼的弹性变形相差不大，因此需要同时对机翼、尾翼进行弹性 建模，考虑弹性影响，在全波长离散阵风的激励下，飞机状态经过1秒左右的震荡逐渐回到稳 定位置，其中，作用在副翼重心上的载荷随时间变化曲线以及作用在飞机重心上的载荷随时间 变化曲线如图21、图22所示。

(二)某大型客机弹性阵风相应模拟验证

在某大型客机的弹性阵风相应的分析中，飞机几何外形如图23所示，机翼为后掠下单翼， 由于飞机较大，机翼的弹性变形显著大于全机其他位置，因此仅需要对机翼进行弹性建模，考 虑其弹性影响，在半波长离散阵风的激励下，飞机状态经过8秒左右的震荡逐渐回到稳定位置， 其中，翼梢位移随时间变化曲线以及飞机迎角随时间变化曲线如图24、图25所示。

参考文献

[1]Schuster D M,Liu D,Huttsell L J.Computational Aeroelasticity:Success,Progress,Challenge[J].Journal of Aircraft,2003,40(5):843～856.

[2]Friedmann P,McNamara J.Aeroelastic and Aerothermoelastic AnalysisofHypersonic Vehicles:Current Status and Future Trends[R].AIAA-2007-2103,2007.

[3]Theodorsen T.General Theory of Aerodynamic Instability andtheMechanism of Flutter[J].Technical Report Archive&Image Library,1935:413～433.

[4]Dugundji J.Personal Perspective of Aeroelasticity During theYears1953-1993[J].Journal of Aircraft,2015,40(40):809～812.

[5]专题译丛：颤振问题(内部资料)[M].国防科技资料编辑部第四室,1961.

[6]Landahl M.Kernel Function for Nonplanar Oscillating Surfaces in aSubsonic Flow[J].AIAA Journal,1967,5(5):1045～1046.

[7]邹丛青,陈桂彬.机翼/外挂颤振主动抑制的控制率研究[J].力学学报,1991,23(3):274～281.

[8]伏欣沈克杨.气动弹性力学原理[M].上海科学技术文献出版社,1982.

[9]邢誉峰,诸德超.杆和梁在锁定过程的响应[J].计算力学学报,1998,15(2):192～196.

[10]Tang D,Henry J K,Dowell E H.Limit Cycle Oscillations of DeltaWing Models in Low Subsonic Flow[J].AIAA Journal,2012,37(37):1355～1362.

[11]蒋跃文,张伟伟,叶正寅,等.基于CFD技术的流场/结构时域耦合求解方法研究 [J].振动工程学报,2007,20(4):396～400.。

Claims (3)

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；对于W>20t的飞机有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法求半翼展的积分:

即可计算得到

矩阵。

2.根据权利要求1所述的弹性飞机阵风响应的仿真分析方法，其特征在于：

(四)将三个模型连成通路，建立模拟系统，最终得到弹性飞机的阵风相应，具体流程为基于所建立的飞机本体、环境和气动力模型，将以上模型按模拟过程中的物理规则连接起来，首先，在每次模拟中，根据六自由度模块得到的飞机飞行状态量，经环境模块计算得到飞机当地所处的风速、温度、重力加速度；将这些物理变量与飞机飞行状态量结合，得到气流坐标系中的飞机速压、迎角、侧滑角等信息；把这些信息放入气动力数据库进行插值，得到飞机的气动系数；在地面效应和气动弹性模块的作用下，气动系数的值得到修正；最后和重力一起，计算飞机飞行的合力以及合力矩，传回六自由度模块，计算下一时刻的飞机飞行状态，形成通路。

3.一种基于权利要求1或2所述方法的弹性飞机阵风响应的仿真分析系统，由构建的三个模型：飞机本体飞行状态模型，飞机飞行环境模型，飞机所受气动力以及飞机气动力的弹性影响模型，连成通路而构成；在Simulink平台上得到弹性飞机的阵风相应。

Images