一种考虑全机气动干扰的直升机飞行特性计算方法

摘要

本发明公开的一种考虑全机气动干扰的直升机飞行特性计算方法，包括耦合气动干扰的直升机配平计算方法、耦合气动干扰的直升机操纵响应计算方法和耦合气动干扰的直升机稳定性计算方法。采用涡面元法和粘性涡粒子法建立的直升机全机气动干扰计算模型，并利用多极展开和涡元重构两种加速技术提高计算效率，完成直升机全机干扰流场计算，同时解决了直升机全机干扰流场计算效率和精度问题。将该全机气动干扰计算模型与直升机飞行动力学模型耦合，建立耦合计算方法，形成一套彻底摆脱气动干扰风洞试验的高效高精度飞行特性计算方法，从而达到提高直升机飞行特性计算精度，缩短气动布局设计周期的目的。本计算方法适用于直升机可用包线内的所有飞行状态。

说明书

技术领域

本发明涉及直升机总体设计技术领域，具体涉及一种考虑全机气动干扰的直升机 飞行特性计算方法，可以应用于直升机概念和初步设计阶段的直升机气动布局设计过程中 飞行性能和飞行品质计算分析。

背景技术

直升机全机气动布局设计是总体设计中的主要内容之一，它包括直升机各部件的 气动外形、参数和相对位置的确定。直升机是否具有良好的飞行特性(包括配平特性、稳定 性、操纵性及飞行性能等)在很大程度上决定于气动布局，因此气动布局设计在直升机总体 设计中处于核心地位，而气动布局设计过程中飞行特性的准确计算是关键。由于直升机气 动布局紧凑，机身、尾翼等处于旋翼尾流当中，气动干扰问题严重，对直升机的飞行特性具 有重要影响，因此为获得最佳气动布局，在飞行特性计算中必须考虑各部件之间的气动干 扰。

现有考虑气动干扰的直升机飞行特性计算方法主要分为两类，第一类为目前主要 采用的方法。第一类方法先通过组合部件风洞试验确定气动干扰因子，然后在计算机身、尾 翼等部件气动力时利用试验获得的干扰因子计入干扰速度的影响。该类方法的主要缺点在 于：

(1)依赖于风洞试验，效率低、经济性差。直升机气动干扰与许多飞行状态参数密 切相关，必须进行大量不同状态的试验；其次，在迭代设计过程中直升机气动布局可能不断 变化，需进行各种变布局的风洞试验，而风洞试验本身周期长，费用高，因此这类计算方法 效率低、经济性差。

(2)依赖于经验公式，难以包括所有飞行状态，人为因素较大。直升机气动干扰与 许多飞行状态参数密切相关，如飞行重量、高度、空速、姿态、旋翼桨盘角度等等，而风洞试 验状态并不能覆盖所有飞行状态，因此通常的方法是利用试验数据根据人为经验建立气动 干扰因子与飞行状态参数的函数关系，进而再应用于飞行特性计算中。

第二类方法是采用旋翼固定或自由尾迹模型计算旋翼尾迹对各部件的气动干扰 因子，进而在机身、尾翼等部件气动力计算时考虑旋翼的干扰作用。这类方法的主要缺点是 气动干扰计算精度差，主要原因在于：

(1)旋翼尾迹模型基于势流理论，忽略粘性耗散，人为设定参数对结果影响大，如 涡核半径等。

(2)只考虑旋翼尾流对其它部件的单向干扰作用，未考虑全机气动干扰，即所有气 动部件及尾流的相互干扰作用。

另外，现有的直升机计算流体力学(CFD)数值计算方法尽管可较准确的计算各部 件之间的气动干扰，但是由于直升机全机气动外形复杂、运动部件较多、涡流特征显著，CFD 计算方法依赖于大规模、精细的、分块的复杂网格系统和各种尾涡捕捉方法，计算量巨大， 即使采用大规模并行计算方法，其计算周期仍以周甚至月为度量单位，难以满足日常工程 设计中飞行特性计算的需求。

发明内容

本发明的目的是为了解决直升机气动布局工程设计中现有的飞行特性计算方法 的不足，提出了一种不依赖于气动干扰风洞试验数据，且计算效率和精度较高的直升机飞 行特性计算方法。

本发明的方法采用基于涡面元法和粘性涡粒子法的直升机全机气动干扰计算模 型(为简化描述，本发明后续部分将该模型简称为直升机全机气动干扰计算模型)，进行直 升机全机干扰流场计算。直升机全机气动干扰计算模型中包括了多极展开和涡元重构两种 加速技术提高计算效率，同时解决了直升机全机干扰流场计算效率和精度问题。将直升机 全机气动干扰计算模型与直升机飞行动力学模型耦合，建立耦合计算方法，形成了一套彻 底摆脱气动干扰风洞试验的高效高精度飞行特性计算方法，从而达到提高直升机飞行特性 计算精度，缩短气动布局设计周期的目的。

本发明的一种考虑全机气动干扰的直升机飞行特性计算方法，通过如下步骤实 现：

11)数据准备：对直升机各孤立部件三维几何模型物面进行网格划分(物面网格为 三角形或四边形平面网格)，导出各孤立部件网格单元节点坐标信息数据表，制作直升机设 计参数以及孤立部件气动特性数据表，设定直升机稳定飞行状态参数；

12)数据输入：将直升机设计参数数据表和孤立部件气动特性数据表输入直升机 飞行动力学模型，将各孤立部件网格单元节点坐标信息数据表输入直升机全机气动干扰计 算模型；

13)耦合气动干扰的直升机配平计算：首先进行直升机飞行动力学模型配平计算 输出配平值，将配平值传递给直升机全机气动干扰计算模型，然后进行直升机全机气动干 扰计算模型配平计算，并将直升机全机气动干扰计算模型配平计算输出的旋翼和尾桨桨盘 诱导速度分布、各部件受到其它部件干扰速度量传递给直升机飞行动力学模型，重复上述 计算过程，直到直升机飞行动力学模型配平计算输出的配平值收敛，得到直升机稳定飞行 状态配平值；

14)耦合气动干扰的直升机操纵响应计算：根据步骤13)输出的配平值对直升机飞 行动力学模型和直升机全机气动干扰计算模型进行初始化，设定操纵响应计算时间、时间 步长以及直升机操纵输入变化量时间历程，首先进行1个时间步的直升机飞行动力学模型 时间积分，将时间积分得到的直升机运动状态量信息传递给直升机全机气动干扰计算模 型，然后进行直升机全机气动干扰计算模型时间积分，将全机气动干扰计算模型时间积分 后输出的旋翼和尾桨桨盘诱导速度分布、各部件受到其它部件干扰速度量传递给直升机飞 行动力学模型，重复上述操纵响应计算过程，直到达到设定的计算时间，得到直升机的操纵 响应；

15)耦合气动干扰的直升机稳定性计算：根据步骤13)输出的配平值对直升机飞行 动力学模型和直升机全机气动干扰计算模型进行初始化，设定直升机运动状态量和操纵输 入量的小扰动量，采用数值差分求导的方法，对直升机飞行动力学模型进行线化，获得直升 机气动导数矩阵和操纵导数矩阵，对直升机气动导数矩阵进行特征值和特征向量计算，得 到直升机各运动模态的稳定根。

作为上述技术方案的进一步改进，所述步骤13)耦合气动干扰的直升机配平计算 具体过程为：

43)配平初值计算：将直升机稳定飞行状态参数输入直升机飞行动力学模型，进 行稳定飞行配平计算，输出配平初值；

44)全机气动干扰计算模型配平计算：将步骤43)输出的配平初值或步骤45)输出 的配平值输入直升机全机气动干扰计算模型，进行全机气动干扰计算模型配平计算，输出 旋翼和尾桨桨盘诱导入流分布、各部件受到其它部件干扰速度量；

45)飞行动力学模型配平计算：将步骤44)输出旋翼和尾桨桨盘诱导入流分布和各 部件受到其它部件干扰速度量输入直升机飞行动力学模型，进行稳定飞行配平计算，输出 配平值；

46)耦合配平计算收敛判断：若步骤45)仅执行过一次，计算步骤43)输出的配平初 值与步骤45)输出的配平值之间的相对残差值，若步骤45)执行次数大于1，计算前后两次步 骤45)输出配平值之间的相对残差值，当相对残差值大于某小量(通常设定为10^-3量级)，配 平计算未收敛，重复所述步骤44)和45)，当相对残差值小于等于某小量，配平计算已收敛， 执行步骤47)；

47)耦合配平计算结果输出：输出最后一次执行步骤45)配平值和最后一次执行步 骤44)过程中各部件涡面元强度、各部件尾迹结构中涡粒子位置和强度数据信息以及旋翼 和尾桨桨盘诱导入流分布、各部件受到其它部件干扰速度量，耦合气动干扰的直升机配平 计算结束。

作为上述技术方案的进一步改进，所述步骤14)耦合气动干扰的直升机操纵响应 计算具体过程为：

52)耦合配平计算：执行耦合气动干扰的直升机配平计算方法所有步骤，即步骤 41)～步骤46)，直到耦合气动干扰的直升机配平计算结束；

53)初始化一：将步骤52)输出的配平值转化为直升机飞行动力学模型中各状态量 和操纵输入量的初始值，将步骤52)输出的旋翼和尾桨桨盘诱导入流分布、各部件受到其它 部件干扰速度量输入直升机飞行动力学模型，设定当前时刻t＝0，完成直升机飞行动力学 模型初始化；

54)初始化二：将步骤52)输出的配平值、各部件涡面元强度、各部件尾迹结构中涡 粒子位置和强度数据信息输入直升机全机气动干扰计算模型，完成直升机全机气动干扰计 算模型初始化；

55)直升机飞行动力学模型时间积分：根据直升机操纵输入变化量时间历程数据 表，读入t时刻直升机操纵输入变化量，采用DE/STEP常微分方程积分算法，对直升机飞行动 力学模型进行1个时间步长积分，得到t+△t时刻直升机运动状态量信息，并将此运动状态 量信息传递给直升机全机气动干扰计算模型；

56)直升机全机气动干扰计算模型时间积分：对直升机全机气动干扰计算模型进 行1个时间步长积分，得到t+△t时刻旋翼和尾桨桨盘诱导入流分布、各部件受到其它部件 干扰速度量、各部件涡面元强度、各部件尾迹结构中涡粒子位置和强度，并将t+△t时刻旋 翼和尾桨桨盘诱导入流分布、各部件受到其它部件干扰速度量传递给直升机飞行动力学模 型；

57)时间积分判断和输出：令t＝t+△t，输出步骤55)得到的直升机运动状态量信 息，若t小于T，重复步骤55)和56)，否者，直升机操纵响应计算结束。

作为上述技术方案的进一步改进，所述步骤15)耦合气动干扰的直升机稳定性计 算具体过程为：

62)耦合配平计算：执行耦合气动干扰的直升机配平计算方法所有步骤，即步骤 41)～步骤46)，直到耦合气动干扰的直升机配平计算结束；

63)初始化一：将步骤62)输出的配平值转化为直升机飞行动力学模型中各状态量 和操纵输入量的初始值，将步骤62)输出的旋翼和尾桨桨盘诱导入流分布、各部件受到其它 部件干扰速度量输入直升机飞行动力学模型，设定当前时刻t＝0，完成直升机飞行动力学 模型初始化；

64)初始化二：将步骤62)输出的配平值、各部件涡面元强度、各部件尾迹结构中涡 粒子位置和强度数据信息输入直升机全机气动干扰计算模型，完成直升机全机气动干扰计 算模型初始化；

65)t时刻气动导数和操纵导数矩阵计算：采用数值求导方法利用直升机飞行动力 学模型，计算t时刻气动导数和操纵导数矩阵，具体步骤为，先计算t时刻直升机运动状态量 时间变化率，在t时刻直升机运动状态量和操纵输入量的基础上，根据设定的直升机运动状 态量和操纵输入量的小扰动量，分别对每一个直升机运动状态量和操纵输入量进行扰动， 计算t时刻直升机运动状态量扰动后的时间变化率，将扰动后和扰动前的每个一个运动状 态量时间变化率与每一个直升机运动状态量和操纵输入量的小扰动量进行数值求导计算， 获得t时刻直升机的气动导数矩阵和操纵导数矩阵，并记录保存；

66)直升机飞行动力学模型时间积分：保持步骤62)输出的配平值中的操纵输入 量，采用DE/STEP常微分方程积分算法，对直升机飞行动力学模型进行1个时间步长积分，得 到t+△t时刻直升机运动状态量信息，并将此运动信息传递给直升机全机气动干扰计算模 型；

67)直升机全机气动干扰计算模型时间积分：对直升机全机气动干扰计算模型进 行1个时间步长积分，得到t+△t时刻旋翼和尾桨桨盘诱导入流分布、各部件受到其它部件 干扰速度量、各部件涡面元强度、各部件尾迹结构中涡粒子位置和强度，并将t+△t时刻旋 翼和尾桨桨盘诱导入流分布、各部件受到其它部件干扰速度量传递给直升机飞行动力学模 型；

68)气动导数和操纵导数矩阵计算判断：令t＝t+△t，若t小于旋翼旋转周期，重复 步骤65)和67)，若t等于旋翼旋转周期，转入步骤69)；

69)稳定性计算：将执行N次步骤65)保存的N个气动导数矩阵和操纵导数矩阵进行 代数平均，获得设定稳定飞行状态下直升机的气动和操纵导数矩阵，对气动导数矩阵进行 特征值和特征向量计算，获得直升机各运动模态的稳定根，耦合气动干扰的直升机稳定性 计算结束。

作为上述技术方案的进一步改进，所述孤立部件包括直升机旋翼、机身、平尾、垂 尾以及尾桨。

作为上述技术方案的另一种改进，所述稳定飞行状态参数包括直升机的飞行高 度、飞行速度、爬升角和偏航角速率。

作为上述技术方案的另一种改进，所述配平初值和配平值包括直升机体轴系下机 体的3个线速度、3个姿态角、3个旋翼操纵输入量、1个尾桨操纵输入量、旋翼桨盘锥度角、旋 翼桨盘后倒角和旋翼桨盘侧倒角，若考虑了旋翼桨叶摆振运动，则还包括旋翼集合型、前进 型和后退型摆振角。

作为上述技术方案的另一种改进，所述直升机运动状态量包括直升机体轴系下机 体的3个线速度、3个角速度、3个姿态角、旋翼参考桨叶方位角、旋翼各片桨叶挥舞角、旋翼 各片桨叶挥舞角速度、旋翼各片桨叶变距角、旋翼各片桨叶变距角速率，若考虑了旋翼桨叶 摆振运动，则还包括旋翼各片桨叶摆振角、旋翼各片桨叶摆振角速度。

作为上述技术方案的另一种改进，所述各部件受到其它部件干扰速度量包括机身 气动中心受到旋翼、平尾、垂尾和尾桨涡面元及其尾迹诱导的气流速度量，平尾气动中心受 到旋翼、机身、垂尾和尾桨涡面元及其尾迹诱导的气流速度量，垂尾气动中心受到旋翼、机 身、平尾和尾桨涡面元及其尾迹诱导气流速度量。

作为上述技术方案的另一种改进，所述旋翼和尾桨桨盘诱导入流分布包括旋翼和 尾桨桨盘平面受到旋翼、机身、平尾、垂尾和尾桨涡面元及其尾迹诱导的气流速度分布。

与现有方法相比，本发明公开的方法的效果和优点如下：

(1)彻底摆脱直升机飞行特性计算依赖气动干扰试验数据和气动干扰经验修正的 现状。

(2)计算精度高：在低速和机动飞行状态，直升机配平值、异轴响应计算精度比已 有理论计算方法提高10％以上。

(3)计算效率高：采用普通个人电脑，稳定平飞状态配平计算时间约1－3小时。

(4)本计算方法适用于直升机可用包线内的所有飞行状态。

附图说明

图1是直升机全机气动干扰计算模型配平计算步骤流程图；

图2是直升机全机气动干扰计算模型时间积分步骤流程图；

图3是耦合气动干扰的直升机配平计算方法步骤流程图；

图4是耦合气动干扰的直升机操纵响应计算方法步骤流程图；

图5是耦合气动干扰的直升机稳定性计算方法步骤流程图；

图6是某直升机各部件物面网格图；

图7.1-7.4是某直升机配平计算结果图，其中，图7.1是旋翼总距计算结果图，图 7.2是旋翼需用功率计算结果图，图7.3是纵向周期变距杆量计算结果图，图7.4是俯仰姿态 角计算结果图；

图8.1-8.4是某直升机操纵输入变化量时间历程及操纵响应计算结果图，其中，图 8.1是操纵输入计算结果图，图8.2是滚转角速度响应计算结果图，图8.3是俯仰角速度响应 计算结果图，图8.4是偏航角速度响应计算结果图。

具体实施方式

本发明的直升机飞行动力学模型，是在已有通用的常规直升机飞行动力学模型基 础上改造而来的。通用常规直升机飞行动力学模型一般可表示为：$\stackrel{·}{y}=f(y,u)---\left(1\right)$

其中，y＝(y_f,y_m,y_tr,y_mi,y_tri,y_e)为直升机状态量向量， 为由直升机旋翼横向周期变距角B₁、旋翼纵向周期变距角A₁、旋 翼总距角θ₀、尾桨总距角θ_tr以及它们的时间变化率所组成的操纵输入向量。直升机状态量 向量由机体运动状态量向量y_f、旋翼运动状态量向量y_m、尾桨运动状态量向量y_tr、旋翼桨盘 入流状态量向量y_mi、尾桨桨盘入流状态量向量y_tri以及发动机状态量向量y_e组成。其中y_f中 包括机体3个线速度分量、3个角速度分量、3个姿态角；y_m中一般需包括旋翼参考桨叶方位 角、旋翼转速、旋翼各片桨叶的挥舞角和角速度、旋翼各片桨叶的摆振角和角速度(若考虑 摆振运动)；y_tr中一般需包括尾桨参考桨叶方位角、尾桨各片桨叶的挥舞角和角速度(若考 虑挥舞运动)；若采用动态入流模型，y_mi和y_tri应分别包括旋翼和尾桨桨盘各阶入流系数；y_e与所选择的发动机模型相关。

本发明对通用常规直升机飞行动力学模型进行改造，以满足直升机全机气动干扰 计算模型与直升机飞行动力学模型耦合求解的需要。改造后的直升机飞行动力学模型表示 为：

$\stackrel{·}{y}=f(y,u,[{\upsilon }_{mi}],[{\upsilon }_{tri}],[{\upsilon }_{ia}],t)---\left(2\right)$

其中t为时间，[υ_mi]和[υ_tri]分别表示t时刻旋翼和尾桨桨盘平面诱导入流速度分 布，[υ_ia]表示t时刻机身、平尾、垂尾等部件气动中心分别受到其它部件引起的干扰气流速 度，而改造后的直升机飞行动力学模型中的直升机状态量向量y＝(y_f,y_m,y_tr,y_e)中不再包 括旋翼和尾桨桨盘入流状态量。

本发明中所述直升机飞行动力学模型均指改造后的直升机飞行动力学模型，所述 直升机飞行动力学模型配平计算和时间积分计算，均指通用常规直升机飞行动力学模型 已有的配平和时间积分算法。

本发明的直升机全机气动干扰计算模型基于涡面元法和粘性涡粒子法而建立，它 由直升机各部件物面网格(各部件主要包括旋翼、机身、平尾、垂尾、尾桨等气动面)、各部件 物面涡面元、各部件尾迹涡粒子构成。

如图1所示，直升机全机气动干扰计算模型配平计算包括以下步骤：

21)数据准备：对直升机各孤立部件三维几何模型物面进行网格划分(物面网格为 三角形或四边形平面网格)，导出各孤立部件网格单元节点坐标信息数据表,设定直升机稳 定飞行配平值；

22)数据输入：将各孤立部件网格单元节点坐标信息数据表输入直升机全机气动 干扰计算模型；

23)初始化：采用涡面元法中面元影响系数矩阵计算方法，计算直升机各孤立部件 涡面元影响系数矩阵，设定当前时刻t＝0，时间步长△t，要求△t为旋翼旋转周期的1/N，将 直升机稳定飞行配平值转化旋翼旋转周期内以△t为步长N组离散的直升机运动状态量信 息；

24)执行N步全机气动干扰气动模型时间积分，记录每步执行后各部件尾迹涡粒子 位置和强度；

25)直升机全机气动干扰计算模型配平收敛判断：若步骤24)仅执行过1次，则执行 步骤24)，若步骤24)执行次数大于1，判断步骤24)前后两次执行所记录的各部件尾迹涡粒 子位置和强度的相对残差值，当相对残差值大于某小量(通常设定为10^-3量级)，配平计算未 收敛，重复所述步骤14)，当相对残差值小于等于某小量，配平计算已收敛，执行步骤26)；

26)输出结果：计算旋翼和尾桨桨盘诱导入流分布，计算各部件受到其它部件干扰 速度量，并输出，结束。

如图2所示，所述直升机全机气动干扰计算模型时间积分由t时刻到t+△t时刻一 步积分包括以下步骤：

31)输入t时刻直升机运动状态量信息；

32)计算各部件涡面元强度：根据t时刻直升机运动状态量信息和各部件尾迹对各 部件涡面元控制点诱导速度，利用各部件影响系数矩阵，计算t时刻各部件涡面元强度；

33)各部件生成新的尾迹涡粒子：根据t时刻直升机运动状态量信息和t-△t时刻 及t时刻各部件涡面元强度，确定各部件新生成尾迹涡粒子的位置和强度；

34)计算各部件尾迹涡粒子运动速度和梯度：基于多极展开加速技术计算t时刻各 部件尾迹涡粒子运动速度和梯度；

35)各部件尾迹涡粒子位置和涡强时间积分；根据t时刻各部件尾迹涡粒子运动速 度和梯度，采用数值积分算法，获得t+△t时刻各部件尾迹涡粒子的位置和强度；

36)涡粒子重构：采用涡粒子重构加速技术，对t+△t时刻各部件尾迹涡粒子进行 合并和分裂；

37)计算各部件涡面元控制点诱导速度：根据t+△t时刻各部件尾迹涡粒子的位置 和强度，计算t+△t时刻各部件尾迹对各部件涡面元控制点诱导速度，令t＝t+△t，则t时刻 到t+△t时刻直升机全机气动干扰计算模型时间积分结束。

如图3所示，耦合气动干扰的直升机配平计算方法包括以下步骤：

41)数据准备：对直升机各孤立部件三维几何模型物面进行网格划分(物面网格为 三角形或四边形平面网格)，导出各孤立部件网格单元节点坐标信息数据表,制作直升机设 计参数以及孤立部件气动特性数据表,设定直升机稳定飞行状态参数；

42)数据输入：将直升机设计参数数据表和孤立部件气动特性数据表输入直升机 飞行动力学模型，将各孤立部件网格单元节点坐标信息数据表输入直升机全机气动干扰计 算模型；

43)配平初值计算：将直升机稳定飞行状态参数输入直升机飞行动力学模型，进行 稳定飞行配平计算，输出配平初值；

44)全机气动干扰计算模型配平计算：将步骤43)输出的配平初值或步骤45)输出 的配平值输入直升机全机气动干扰计算模型，进行全机气动干扰计算模型配平计算，输出 旋翼和尾桨桨盘诱导入流分布、各部件受到其它部件干扰速度量；

45)飞行动力学模型配平计算：将步骤44)输出旋翼和尾桨桨盘诱导入流分布和 各部件受到其它部件干扰速度量输入直升机飞行动力学模型，进行稳定飞行配平计算，输 出配平值；

46)耦合配平计算收敛判断：若步骤45)仅执行过一次，计算步骤43)输出的配平初 值与步骤45)输出的配平值之间的相对残差值，若步骤45)执行次数大于1，计算前后两次步 骤45)输出配平值之间的相对残差值，当相对残差值大于某小量(通常设定为10^-3量级)，配 平计算未收敛，重复所述步骤44)和45)，当相对残差值小于等于某小量，配平计算已收敛， 执行步骤47)；

47)耦合配平计算结果输出：输出最后一次执行步骤45)配平值和最后一次执行步 骤44)过程中各部件涡面元强度、各部件尾迹结构中涡粒子位置和强度数据信息以及旋翼 和尾桨桨盘诱导入流分布、各部件受到其它部件干扰速度量，耦合气动干扰的直升机配平 计算结束。

如图4所示，耦合气动干扰的直升机操纵响应计算方法包括以下步骤：

51)数据准备：制作直升机操纵输入变化量时间历程数据表，设定操纵响应计算时 间T和时间步长△t；

52)耦合配平计算：执行耦合气动干扰的直升机配平计算方法所有步骤，即步骤 41)～步骤46)，直到耦合气动干扰的直升机配平计算结束；

53)初始化一：将步骤52)输出的配平值转化为直升机飞行动力学模型中各状态量 和操纵输入量的初始值，将步骤52)输出的旋翼和尾桨桨盘诱导入流分布、各部件受到其它 部件干扰速度量输入直升机飞行动力学模型，设定当前时刻t＝0，完成直升机飞行动力学 模型初始化；

54)初始化二：将步骤52)输出的配平值、各部件涡面元强度、各部件尾迹结构中涡 粒子位置和强度数据信息输入直升机全机气动干扰计算模型，完成直升机全机气动干扰计 算模型初始化；

55)直升机飞行动力学模型时间积分：根据直升机操纵输入变化量时间历程数据 表，读入t时刻直升机操纵输入变化量，采用DE/STEP常微分方程积分算法，对直升机飞行动 力学模型进行1个时间步长积分，得到t+△t时刻直升机运动状态量信息，并将此运动状态 量信息传递给直升机全机气动干扰计算模型；

56)直升机全机气动干扰计算模型时间积分：对直升机全机气动干扰计算模型进 行1个时间步长积分，得到t+△t时刻旋翼和尾桨桨盘诱导入流分布、各部件受到其它部件 干扰速度量、各部件涡面元强度、各部件尾迹结构中涡粒子位置和强度，并将t+△t时刻旋 翼和尾桨桨盘诱导入流分布、各部件受到其它部件干扰速度量传递给直升机飞行动力学模 型；

57)时间积分判断和输出：令t＝t+△t，输出步骤55)得到的直升机运动状态量信 息，若t小于T，重复步骤55)和56)，否者，直升机操纵响应计算结束。

如图5所示，耦合气动干扰的直升机稳定性计算方法包括以下步骤：

61)数据准备：设定时间步长△t，要求△t为旋翼旋转周期的1/N,设定直升机运动 状态量和操纵输入量的小扰动量；

62)耦合配平计算：执行耦合气动干扰的直升机配平计算方法所有步骤，即步骤 41)～步骤46)，直到耦合气动干扰的直升机配平计算结束；

63)初始化一：将步骤62)输出的配平值转化为直升机飞行动力学模型中各状态量 和操纵输入量的初始值，将步骤62)输出的旋翼和尾桨桨盘诱导入流分布、各部件受到其它 部件干扰速度量输入直升机飞行动力学模型，设定当前时刻t＝0，完成直升机飞行动力学 模型初始化；

64)初始化二：将步骤62)输出的配平值、各部件涡面元强度、各部件尾迹结构中涡 粒子位置和强度数据信息输入直升机全机气动干扰计算模型，完成直升机全机气动干扰计 算模型初始化；

65)t时刻气动导数和操纵导数矩阵计算：采用数值求导方法利用直升机飞行动力 学模型，计算t时刻气动导数和操纵导数矩阵，具体步骤为，先计算t时刻直升机运动状态量 时间变化率，在t时刻直升机运动状态量和操纵输入量的基础上，根据设定的直升机运动状 态量和操纵输入量的小扰动量，分别对每一个直升机运动状态量和操纵输入量进行扰动， 计算t时刻直升机运动状态量扰动后的时间变化率，将扰动后和扰动前的每个一个运动状 态量时间变化率与每一个直升机运动状态量和操纵输入量的小扰动量进行数值求导计算， 获得t时刻直升机的气动导数矩阵和操纵导数矩阵，并记录保存；

66)直升机飞行动力学模型时间积分：保持步骤62)输出的配平值中的操纵输入 量，采用DE/STEP常微分方程积分算法，对直升机飞行动力学模型进行1个时间步长积分，得 到t+△t时刻直升机运动状态量信息，并将此运动信息传递给直升机全机气动干扰计算模 型；

67)直升机全机气动干扰计算模型时间积分：对直升机全机气动干扰计算模型进 行1个时间步长积分，得到t+△t时刻旋翼和尾桨桨盘诱导入流分布、各部件受到其它部件 干扰速度量、各部件涡面元强度、各部件尾迹结构中涡粒子位置和强度，并将t+△t时刻旋 翼和尾桨桨盘诱导入流分布、各部件受到其它部件干扰速度量传递给直升机飞行动力学模 型；

68)气动导数和操纵导数矩阵计算判断：令t＝t+△t，若t小于旋翼旋转周期，重复 步骤65)和67)，若t等于旋翼旋转周期，转入步骤69)；

69)稳定性计算：将执行N次步骤65)保存的N个气动导数矩阵和操纵导数矩阵进行 代数平均，获得设定稳定飞行状态下直升机的气动和操纵导数矩阵，对气动导数矩阵进行 特征值和特征向量计算，获得直升机各运动模态的稳定根，耦合气动干扰的直升机稳定性 计算结束。

应用例：

本实例中，将本发明提出的一种考虑全机气动干扰的直升机飞行特性计算方法应 用于实际的某直升机配平、操纵响应和稳定性计算中。

71)配平计算

如图6所示，应用的模拟直升机包括机身1、旋翼2、平尾3、垂尾4和尾桨5。照具体实 施方式中耦合气动干扰的直升机配平计算方法实施步骤，首先采用CATIA或Gabit等商用软 件中的网格划分功能对各孤立部件三维几何模型物面进行网格划分(如图6所示)，导出各 孤立部件网格单元节点坐标信息数据表，制作直升机设计参数表(如表1所示)和各孤立部 件气动特性数据表(数据量太大，在本应用例中未给出)，设定8组直升机稳定平直飞行状态 参数(如表2所示)。然后，按照耦合气动干扰的直升机配平计算方法余下各步骤对8组稳定 飞行状态参数分别进行计算，得到了对应稳定飞行状态的配平值，主要配平值如图7.1至 7.4所示，图中同时给出了飞行试验结果和不考虑全机气动干扰的飞行动力学模型计算结 果，由图中对比可知本本发明提出的方法，大幅提高了直升机配平计算精度，尤其是在低 速飞行状态。

表1

表2

72)操纵响应计算

按照具体实施方式中耦合气动干扰的直升机操纵响应计算方法实施步骤，首先设 定直升机前飞右压杆操纵输入变化量时间历程(如图8.1)，设定操纵响应计算时间T＝6s和 时间步长△t＝0.01s，选择表2中的第8组稳定飞行状态参数，按照耦合气动干扰的直升机 操纵响应计算方法余下各步骤进行操纵响应计算，得到直升机的操纵响应(如图8.2-图 8.4)，图中同时给出了飞行试验结果和不考虑全机气动干扰的飞行动力学模型计算结果， 由图中对比可知本发明提出的方法有效提高了机动飞行操纵响应计算精度，尤其是异轴响 应的计算精度。

73)稳定性计算

按照具体实施方式中耦合气动干扰的直升机稳定性计算方法实施步骤，首先设定 直升机运动状态量和操纵输入量的小扰动量，对于速度量纲设置扰动量为6.5×10^-3，对于 角度和角速度量纲扰动量设置为6.5×10^-2，设定时间步长△t＝0.00727，则N＝32，选择表2 中的第8组稳定飞行状态参数按照耦合气动干扰的直升机操纵响应计算方法余下各步骤进 行稳定性计算，得到直升机各运动模态的稳定根，如表3所示。

表3

运动模态 特征根实部 特征根虚部 俯仰角速度模态 -1.93 0 垂向运动模态 -0.243 0 沉浮运动模态 -0.0985 ±0.274

滚转角速度模态 -3.655 0 螺旋运动模态 -0.191 0 荷兰滚运动模态 -0.53 ±1.53 偏航运动模态 0 0

本发明具体应用途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于 本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进，这 些改进也应视为本发明的保护范围。