CN102520620A - 单旋翼直升机/涡轴发动机通用综合模型的构建方法 - Google Patents

Abstract

本发明公开了一种单旋翼直升机/涡轴发动机通用综合模型的构建方法。该方法包括：分别建立直升机和发动机的非线性模型；单旋翼直升机/涡轴发动机的稳态配平；单旋翼直升机/涡轴发动机的动态计算；单旋翼直升机/涡轴发动机的实际运行解算。其中稳态配平采用Levenberg-Marquarat算法；动态计算中采用四阶龙格-库塔法求解直升机运动方程组，采用一次通过算法求解涡轴发动机动态平衡方程；实际运行状态解算中采用修正的线性二次型调解器算法进行直升机解算，采用串级PID+扭矩前馈+总距补偿及扭振滤波器方法进行发动机解算及负载滤波。本发明能够满足单旋翼直升机/涡轴发动机通用模型及控制方法的试验需求，对于缩短研制周期、降低试验风险及成本有着积极促进作用。

Description

单旋翼直升机/涡轴发动机通用综合模型的构建方法

技术领域

本发明属于航空宇航推进理论与工程中的系统控制与仿真领域，具体涉及一种单旋翼直升机/涡轴发动机通用综合模型的构建方法。

背景技术

在航空领域，数字仿真已发展为复杂系统设计与验证不可或缺的一环，且为首要的一环，通过数字仿真验证可有效降低物理试验的风险及成本，并大幅缩短研制周期。直升机是一个多自由度、强耦合的高阶系统，包括直升机子系统、发动机子系统以及控制系统等，其稳定性、操纵性和机动性与固定翼飞机相比均较差，随着新一代直升机对机动性、生存性及操控品质等要求的不断提高，传统的将直升机与发动机分开设计的方法已不能满足需求，只有充分考虑其间耦合作用采用综合建模和控制，才能充分发挥发动机的性能潜力，提高直升机的飞行品质。然而若要准确表达真实的飞行状态或参数，体现综合模型性能，首先则需要一个置信度高、性能良好的直升机/发动机综合仿真模型，国内几个主要的直升机设计研究所、航空高等院校均先后提出对该类综合仿真模型建模技术的需求。直升机/涡轴发动机综合仿真模型的建模技术主要包括：直升机非线性建模、涡轴发动机非线性建模、直升机/涡轴发动机稳态配平、直升机/涡轴发动机动态计算以及直升机/涡轴发动机综合模型实际运行解算五个主要环节。

经对现有技术的文献检索发现，王志峰等在《飞行力学》1999年第17卷第4期上发表了“直升机飞行/推进综合模型的建立”，该文中的直升机模型与发动机模型均用线性模型来表征，难以在大包线范围内反映直升机/发动机的非线性动静态特性；南京航空航天大学的旋翼国家重点实验室近年来对直升机/发动机功率匹配作了较深入的研究，如陈仁良1998年博士论文“直升机飞行动力学数学建模及机动性研究”，孙传伟2002年博士论文“直升机飞行动力学模型与飞行品质评估”等，但是这些成果仅实现了直升机建模技术研究，在发动机方面并没有实质性的突破；南京航空航天大学飞行/推进综合控制实验室姚文荣2008年博士论文“涡轴发动机/旋翼综合建模控制与优化研究”实现了带旋翼的涡轴发动机建模，并不包含完整的直升机模型，很难在相关直升机/发动机综合系统设计及研发中发挥出更大效益；国外虽已有成熟的GENHEL、FLYRT和SARC直升机飞行动态仿真模型，但偏重于直升机飞行品质方面的研究，与之匹配发动机模型却采用简化模型无法充分而准确的反映发动机动静态特性。

发明内容

本发明所要解决的技术问题在于克服现有技术的不足，提供一种单旋翼直升机/涡轴发动机通用综合模型的构建方法，用于对不同型号的直升机或发动机综合系统的模型特征进行仿真，可进一步用于综合系统的控制研究(如闭环控制、优化控制、抗扰控制等)。

本发明的单旋翼直升机/涡轴发动机通用综合模型的构建方法，包括以下步骤：

步骤A、分别建立直升机和发动机的非线性模型；

步骤B、单旋翼直升机/涡轴发动机的稳态配平，包括旋翼动力学系统配平、机体动力学模型配平、发动机模型配平；

步骤C、单旋翼直升机/涡轴发动机的动态计算；

步骤D、单旋翼直升机/涡轴发动机的实际运行解算。

所述步骤B具体包括以下步骤：

步骤B1、旋翼动力学系统配平，具体按照以下方法：首先计算每个翼型的迎角，其中旋翼诱导速度采用非均匀一阶谐波入流的计算方法，桨叶挥舞运动则基于桨叶的气流环境分析旋翼桨叶的非定常挥舞运动，以递推形式得到挥舞运动方程，然后根据翼型迎角计算其二元翼型气动力，再通过积分求解每片桨叶的气动力，最终计算出旋翼的气动力，完成直升机旋翼模型的稳态计算；

步骤B2、机体动力学模型配平，具体按照以下方法：进行直升机旋翼、机身、尾桨、水平尾翼和垂直尾翼的动力学解算，分别计算出各部件在直升机重心处产生的力和力矩；根据计算出的各部件在直升机重心处产生的力和力矩，以及直升机重力，解算如下机体平衡方程组：

$\left.\begin{array}{c}{X}_{\mathrm{SUM}}-G{\sin \mathrm{\θ}}_b=0\\ Y_{\mathrm{SUM}}-{G\cos \mathrm{\θ}}_b\sin {\mathrm{\φ}}_b=0\\ Z_{\mathrm{SUM}}+G\cos {\mathrm{\θ}}_b\cos {\mathrm{\φ}}_b=0\\ L_{\mathrm{SUM}}=0\\ M_{\mathrm{SUM}}=0\\ N_{\mathrm{SUM}}=0\end{array}\right\},$

其中，G表示直升机重量，θ_b表示机体的俯仰角，φ_b表示机体的滚转角，下标b表示机体参数，X_SUM＝X_MR+X_WF+X_TR+X_HT+X_VT、Y_SUM＝Y_MR+Y_WF+Y_TR+Y_HT+Y_VT与Z_SUM＝Z_MR+Z_WF+Z_TR+Z_HT+Z_VT分别表示直升机体轴系下X/Y/Z方向的合力，L_SUM＝L_MR+L_WF+L_TR+L_HT+L_VT、M_SUM＝M_MR+M_WF+M_TR+M_HT+M_VT与N_SUM＝N_MR+N_WF+N_TR+N_HT+N_VT分别为其对应的合力矩，下标SUM表示合力或合力矩，MR表示旋翼参量，WF表示机身参量，TR表示尾桨参量，HT表示水平尾翼参量，VT表示垂直尾翼参量；

步骤B3、发动机模型配平，具体按照以下方法：首先进行涡轴发动机气动热力计算；然后根据气动热力计算结果解算如下的发动机稳态平衡方程组：

$\left.\begin{array}{c}E1=\frac{{\stackrel{\‾}{P}}_g}{{\stackrel{\‾}{P}}_c/{\mathrm{\η}}_c+{\stackrel{\‾}{P}}_{\mathrm{ext}}}-1=0\\ E2=\frac{{\stackrel{\‾}{P}}_p}{{\stackrel{\‾}{P}}_l/{\mathrm{\η}}_l}-1=0\\ E3=\frac{W_{41\mathrm{cX}}-W_{41c}}{W_{41c}}=0\\ E4=\frac{W_{44\mathrm{cX}}-W_{44c}}{W_{44c}}=0\\ E5=\frac{P_8^{\′}}{P_8}-1=0\end{array}\right\},$

其中，E1～E5分别表示压气机与燃气涡轮功率平衡、负载与动力涡轮功率平衡、燃气涡轮进口流量连续、动力涡轮进口流量连续、尾喷口压力平衡，

分别为旋翼负载需求功率、压气机功率、燃气涡轮功率、动力涡轮功率、附件功率，η_l，η_c分别为旋翼负载和压气机效率，W_41c，W_44c分别为燃气涡轮和动力涡轮进口换算流量，下标X表示初猜值，P₈，P′₈分别为尾喷口总压和环境压力，压气机相对转速

动力涡轮相对转速

燃气涡轮压比系数Z_G、动力涡轮压比系数Z_P、压气机压比系数Z_C为猜值参数。

优选地，所述机体平衡方程组和发动机稳态平衡方程组的解算均采用Levenberg-Marquarat算法。

所述单旋翼直升机/涡轴发动机的动态计算包括直升机动态计算和涡轴发动机动态计算，其具体方法分别如下：

直升机动态计算：

根据直升机各部件上的力和力矩，解算如下的直升机运动微分方程组：

${\stackrel{\·}{V}}_x=(X_{\mathrm{SUM}}-G{\sin \mathrm{\θ}}_b)g/G+{\mathrm{rV}}_y-q{V}_z$

${\stackrel{\·}{V}}_y=(Y_{\mathrm{SUM}}-G\cos {\mathrm{\θ}}_b\sin {\mathrm{\φ}}_b)g/G+p{V}_z-r{V}_x$

${\stackrel{\·}{V}}_z=(Z_{\mathrm{SUM}}+G\cos {\mathrm{\θ}}_b\cos {\mathrm{\φ}}_b)g/G+q{V}_x-p{V}_y$

$\stackrel{\·}{p}=\frac{I_Z}{I_X{I}_X-I_{\mathrm{XZ}}^2}\{L_{\mathrm{SUM}}-(I_Z-I_Y)\mathrm{qr}+I_{\mathrm{XZ}}\mathrm{pq}\}+\frac{I_{\mathrm{XZ}}}{I_X{I}_Z-I_{\mathrm{XZ}}^2}\{N_{\mathrm{SUM}}-(I_Y-I_X)\mathrm{pq}-I_{\mathrm{XZ}}\mathrm{rq}\}$

$\stackrel{\·}{q}=\frac{1}{I_Y}\{M_{\mathrm{SUM}}-(I_X-I_Z)\mathrm{pr}+I_{\mathrm{XZ}}(r^2-p^2)\},$

$\stackrel{\·}{r}=\frac{I_X}{I_X{I}_Z-I_{\mathrm{XZ}}^2}\{N_{\mathrm{SUM}}-(I_Y-I_X)\mathrm{pq}-I_{\mathrm{XZ}}\mathrm{rq}\}+\frac{I_{\mathrm{XZ}}}{I_X{I}_Z-I_{\mathrm{XZ}}^2}\{L_{\mathrm{SUM}}-(I_Z-I_Y)\mathrm{qr}+I_{\mathrm{XZ}}\mathrm{pq}\}$

${\stackrel{\·}{\mathrm{\θ}}}_b=57.3(q{\cos \mathrm{\φ}}_b-r\sin {\mathrm{\φ}}_b)$

${\stackrel{\·}{\mathrm{\ψ}}}_b=57.3\frac{r\cos {\mathrm{\φ}}_b+q\sin {\mathrm{\φ}}_b}{\cos {\mathrm{\θ}}_b}$

${\stackrel{\·}{\mathrm{\φ}}}_b=57.3p+{\stackrel{\·}{\mathrm{\ψ}}}_b\sin {\mathrm{\θ}}_b$

其中， $V_x=\∫{\stackrel{\·}{V}}_x\mathrm{dt},$ $V_y=\∫{\stackrel{\·}{V}}_y\mathrm{dt},$ $V_z={\∫\stackrel{\·}{V}}_z\mathrm{dt},$ $p=\∫\stackrel{\·}{p}\mathrm{dt},$ $q=\∫\stackrel{\·}{q}\mathrm{dt},$ $r=\∫\stackrel{\·}{r}\mathrm{dt}$ 分别为直升机体轴系下速度和角速度， ${\mathrm{\θ}}_b=\∫{\stackrel{\·}{\mathrm{\θ}}}_b\mathrm{dt},$ ${\mathrm{\φ}}_b={\∫\stackrel{\·}{\mathrm{\φ}}}_b\mathrm{dt},$ ${\mathrm{\ψ}}_b={\∫\stackrel{\·}{\mathrm{\ψ}}}_b\mathrm{dt}$ 分别为机体的俯仰角、滚转角、偏航角，g为重力加速度，G为直升机重量，X_SUM，Y_SUM，Z_SUM，L_SUM，M_SUM，N_SUM分别为直升机体轴系下的合力和合力矩，I_X，I_Y，I_Z，I_XY，I_XZ，I_YZ分别为机体惯性矩和惯性积；涡轴发动机动态计算：

解算如下发动机动态平衡方程组：

$\left.\begin{array}{c}E3=\frac{W_{41\mathrm{cX}}-W_{41c}}{W_{41c}}=0\\ E4=\frac{W_{44\mathrm{cX}}-W_{44c}}{W_{44c}}=0\\ E5=\frac{P_8^{\′}}{P_8}-1=0\\ \frac{{\mathrm{dn}}_g}{\mathrm{dt}}=(W_g*{\mathrm{\η}}_g-W_c)/[n_g*J_g*{\left(\frac{\mathrm{\π}}{30}\right)}^2]\\ \frac{{\mathrm{dn}}_p}{\mathrm{dt}}=(W_p*{\mathrm{\η}}_p-W_l)/[n_p*J_p*{\left(\frac{\mathrm{\π}}{30}\right)}^2]\end{array}\right\},$

其中，E3～E5分别表示燃气涡轮进口流量连续、动力涡轮进口流量连续、尾喷口压力平衡，W_41c，W_44c分别为燃气涡轮和动力涡轮进口换算流量，下标X表示初猜值，P₈，P′₈分别为尾喷口总压和环境压力；燃气涡轮压比系数Z_G、动力涡轮压比系数Z_P、压气机压比系数Z_C为猜值参数，J_g，J_p分别表示燃气涡轮、动力涡轮的转子转动惯量，η_g，η_p表示燃气涡轮、功率涡轮的转子机械效率，n_g，n_p燃气涡轮和动力涡轮转子转速，m_g，m_c，m_p，m_l表示燃气涡轮、压气机、动力涡轮、旋翼负载扭矩，W_g，W_p分别表示燃气涡轮和动力涡轮功率。

优选地，所述直升机运动微分方程组的解算采用四阶龙格-库塔法。

优选地，所述发动机动态平衡方程组的解算采用一次通过算法。

所述单旋翼直升机/涡轴发动机的实际运行解算包括直升机实际运行解算和发动机实际运行解算，其具体方法分别如下：

直升机实际运行解算：

采用修正的线性二次型调节器方法，具体如下：

直升机的参考指令为常值信号r，则误差为e＝r-y，对误差信号微分可得增广系统：

$\left[\begin{array}{c}\stackrel{\·\·}{x}\\ \stackrel{\·}{e}\end{array}\right]=\left[\begin{array}{c}\stackrel{\·\·}{x}\\ \stackrel{\·}{r}-\stackrel{\·}{y}\end{array}\right]=\left[\begin{array}{cc}A& 0\\ -C& 0\end{array}\right]\left[\begin{array}{c}\stackrel{\·}{x}\\ e\end{array}\right]+\left[\begin{array}{c}B\\ -D\end{array}\right]\stackrel{\·}{u}$

其中，x为直升机体轴系下三个方向的速度、角速度以及机体的俯仰角、滚转角、偏航角，共9个状态量；u为总矩、横向周期变距、纵向周期变距和尾桨总矩4个控制输入；A，B，C，D为适维系统状态矩阵；取增广状态量为： $\stackrel{\‾}{x}={\left[\begin{array}{cc}{\stackrel{\·}{x}}^T& e^T\end{array}\right]}^T,$ 控制量为：

则 $\stackrel{\·}{\stackrel{\‾}{x}}=\stackrel{\‾}{A}\stackrel{\‾}{x}+\stackrel{\‾}{B}\stackrel{\‾}{u},$ 其中 $\stackrel{\‾}{A}=\left[\begin{array}{cc}A& 0\\ -C& 0\end{array}\right],$ $\stackrel{\‾}{B}=\left[\begin{array}{c}B\\ -D\end{array}\right];$ 对该增广系统进行线性二次型调节器设计，使其所有状态变量都保持为零： $\stackrel{\‾}{x}={\left[\begin{array}{cc}{\stackrel{\·}{x}}^T& e^T\end{array}\right]}^T=0,$ 且参考指令与系统输出之差e为零；而后用线性二次型调节器法求出的解算公式为

将

按

和e表示为分块矩阵并取拉氏变换，可得实际运行系统控制输入u解算公式为：

$u=-K_{\stackrel{\·}{x}}x-\frac{K_{e}e}{s};$

发动机实际运行解算：

解算如下公式并接入扭振滤波器进行负载滤波：

$\mathrm{\Δ}{W}_f={\mathrm{\Σ}}_{m=1}^2(K_{\mathrm{pm}}{e}_m+K_{\mathrm{im}}\mathrm{\Σ}{e}_m)+K_1\mathrm{\ΔT}+K_2\mathrm{\Δ}{\mathrm{\θ}}_0,$

其中，K_p1＝0.3，K_i1＝0.29为串级PID外环参数，K_p2＝0.002，K_i2＝0.000为内环参数，e₁为动力涡轮转速偏差量，e₂为燃气涡轮转速偏差量，ΔT为扭矩偏差量，Δθ₀为总距偏差量，K₁、K₂为可调系数，K₂根据总距杆角度的大小插值得到，K₁为一给定的参数。

优选地，所述扭振滤波器采用带通滤波与均值滤波串联的形式。

相比现有技术，本发明具有以下有益效果：

(1)模型具有通用性及可移植性。单旋翼直升机与涡轴发动机模型均是基于特性数据建立的非线性实时模型，具有通用性，仅需将特性数据作简单替换，即可实现不同型号的直升机或发动机的建模，而且该模型可与半物理仿真试验完成无缝结合，达到很好的仿真效果。

(2)模型具有初始配平任意性。采用L-M算法进行综合系统稳态配平的求解，对初猜解不敏感，无需求逆，非线性适应性能力强，从设计点出发就可求解综合系统任意状态的配平解，可在包线范围内任意点进行功率解算和初始化。

(3)模型具有较好动态计算精度及实时性。采用四阶龙格-库塔法求解直升机运动方程组，提高了模型计算精度，并且采用一次通过算法求解涡轴发动机动态工作方程组，改善了模型实时性。

(4)直升机/发动机实际运行解算方法的先进性。采用引入积分项可以实现无静差跟踪的修正线性二次调节器方法进行直升机实际运行状态解算，实现较优的控制效果及鲁棒性；采用串级PID+扭矩前馈+总距补偿方法进行涡轴发动机实际运行状态解算，结构简单且具有很强工程实用性，并接入扭振滤波器对旋翼负载进行滤波以降低对发动机控制器带宽的要求，综合作用下提升直升机的飞行品质。

附图说明

图1是直升机/涡轴发动机通用仿真模型结构图；

图2是直升机/发动机稳态配平流程图；

图3是旋翼非定常气动力特性原理图；

图4是直升机/发动机动态计算流程图；

图5是扭振滤波器结构图；

图6是直升机/发动机稳态配平效果图，其中图(a)为总距θ₀配平图，图(b)为横向周期变距A_1C配平图，图(c)为纵向周期变距B_1S配平图，图(d)为尾桨总距θ_T配平图，图(e)为俯仰角θ配平图，图(f)为滚转角□配平图；

图7是直升机/发动机动态计算效果图，其中，图(a)为总距θ₀+1°阶跃响应图，图(b)为横向周期变距A_1C+0.5°阶跃响应图，图(c)为纵向周期变距B_1S+0.5°阶跃响应图，图(d)为尾桨总矩θ_T+1°阶跃响应图；

图8是直升机/发动机综合系统半物理仿真平台结构图；

图9是直升机/发动机综合系统半物理仿真效果图。

具体实施方式

下面结合附图对本发明的技术方案进行详细说明：

本具体实施方式以构建黑鹰UH-60A直升机/T700涡轴发动机综合系统为例，如图1为直升机/涡轴发动机通用仿真模型结构图所示，主要包括相互耦合的直升机子系统和发动机子系统，直升机模型是依照美国NASA提出的直升机通用实时非线性模型建模方法建立的包含主旋翼、机身、尾桨、平尾和垂尾的非定常、非线性气动力模型。其中旋翼模型建立根据叶素理论，并采用了非均匀一阶谐波入流模型，能准确反映旋翼的挥舞、摆振等非线性非定常动力学行为。涡轴发动机部件级模型是用部件法建立的非线性实时数学模型，且能够模拟压气机导叶可调功能。本发明主要研究直升机/涡轴发动机的稳态配平、动态计算及直升机/涡轴发动机的实际运行解算方法设计。

本发明的直升机/涡轴发动机稳态配平包括旋翼动力学系统配平、机体动力学模型配平及发动机模型配平三个部分，参照附图2所示的直升机/发动机稳态配平流程图，其实现过程如下：

步骤B1、旋翼动力学系统配平；

直升机旋翼模型是根据叶素理论以等环形面积作为桨叶微段划分准则而建立，对照附图3所示的旋翼非定常气动力特性原理图，首先计算每个翼型的迎角，其中旋翼诱导速度采用非均匀一阶谐波入流的计算方法，桨叶挥舞运动则基于桨叶的气流环境分析旋翼桨叶的非定常挥舞运动，以递推形式得到挥舞运动方程，然后根据翼型迎角计算其二元翼型气动力，再通过积分求解每片桨叶的气动力，最终计算出旋翼的气动力，完成直升机旋翼模型的稳态计算。

步骤B2、机体动力学模型配平；

机体动力学模型配平包括旋翼总矩、纵横向周期变距、尾桨总矩和平尾安装角的配平，本具体实施方式中采用L-M算法，其实现过程如下：

(a)机身模型动力学计算：主要考虑旋翼下洗的影响，根据当地气流速度与旋翼尾迹在机身处干扰速度的叠加求取机身处的气流速度，进而可计算机身的攻角、侧滑角及动压，机身的力和力矩是根据风洞试验数据通过机身攻角和侧滑角插值获得。

(b)尾桨模型动力学计算：只计算尾桨拉力，忽略对直升机作用很小的尾桨后向力及反扭矩。首先计算旋翼与机身在尾桨处的干扰速度求取尾桨处的合速度，而后采用简化的封闭形式的Bailey解(即忽略了所有前进比二次方及以上的项)来解算尾桨拉力。

(c)尾翼模型动力学计算：包括水平尾翼及垂直尾翼，两者计算方法相同。首先通过主旋翼和机身在水平尾翼(或垂直尾翼)处的干扰速度，求取水平尾翼(或垂直尾翼)处的合速度，进而计算出水平尾翼(或垂直尾翼)的攻角和侧滑角，并由其插值出水平尾翼(或垂直尾翼)的气动力和力矩系数计算各自的力和力矩。

(d)机体稳态配平：通过上述直升机旋翼、机身、尾桨、水平尾翼和垂直尾翼的动力学解算，分别计算出各部件在直升机重心处产生的力和力矩，再加上直升机本身的重力，要保证直升机稳定平衡，就是令机体加速度、角速度、角加速度为零，通过优化算法求解机体平衡方程，优化变量为旋翼总距操纵量、横向周期变距、纵向周期变距、尾桨总距、机体俯仰角、机体滚转角，最终发动机模型根据直升机的需用功率解算出发动机的工作状态。即保证以下的机体平衡方程组成立：

$\left.\begin{array}{c}{X}_{\mathrm{SUM}}-G{\sin \mathrm{\θ}}_b=0\\ Y_{\mathrm{SUM}}-{G\cos \mathrm{\θ}}_b\sin {\mathrm{\φ}}_b=0\\ Z_{\mathrm{SUM}}+G\cos {\mathrm{\θ}}_b\cos {\mathrm{\φ}}_b=0\\ L_{\mathrm{SUM}}=0\\ M_{\mathrm{SUM}}=0\\ N_{\mathrm{SUM}}=0\end{array}\right\}---\left(1\right)$

其中，G表示直升机重量，θ_b表示机体的俯仰角，φ_b表示机体的滚转角，下标b表示机体参数，X_SUM＝X_MR+X_WF+X_TR+X_HT+X_VT、Y_SUM＝Y_MR+Y_WF+Y_TR+Y_HT+Y_VT与Z_SUM＝Z_MR+Z_WF+Z_TR+Z_HT+Z_VT分别表示直升机体轴系下X/Y/Z方向的合力，L_SUM＝L_MR+L_WF+L_TR+L_HT+L_VT、M_SUM＝M_MR+M_WF+M_TR+M_HT+M_VT与N_SUM＝N_MR+N_WF+N_TR+N_HT+N_VT分别为其对应的合力矩，下标SUM表示合力或合力矩，MR表示旋翼参量，WF表示机身参量，TR表示尾桨参量，HT表示水平尾翼参量，VT表示垂直尾翼参量。

上述机体平衡方程组(1)可采用现有的多种方法解算，本具体实施方式中采用L-M算法进行解算，具体如下：

Step1：选取β∈(0，1)，初始正则化实参数μ₀，增长因子v＞1；

Step2：给定初猜值x₀，设置迭代次数k＝0，其中x为直升机机体模型初猜变量；

Step3：计算f_k＝f(x_k)，

其中f为机体平衡方程组(1)；

Step4：计算A_k＝Df(x_k)，

其中D为雅克比矩阵；

Step5：计算中间参数

Step6：解方程(Q_k+μ_kI)p_k＝-g_k，其中p_k为当前下降方向，μ_k为修正实因子

Step7：计算x_k+1＝x_k+p_k(μ_k)，f_k+1＝f(x_k)，

Step8：若满足H-收敛准则(或Himmelblau终止准则)，则输出x_k+1，S_k+1；

Step9：检验：若成立，则令μ_k+1＜μ_k/v，转Step10，否则令μ_k+1＝vμ_k，转Step6；

Step10：k＝k+1，转Step3。

步骤B3、发动机模型配平；

发动机模型配平就是根据旋翼的需用功率进行涡轴发动机功率解算模型的配平计算，求解出涡轴发动机当前的工作状态。涡轴发动机的功率解算模型的配平包括两个部分：一是沿气流方向的发动机各部件气动热力计算，二是建立反映稳态时各部件共同工作关系的平衡方程组并解算，其实现过程如下：

(a)涡轴发动机气动热力计算：根据气动热力学特性及压气机、燃气涡轮、功率涡轮等的部件特性，可建立包括进气道、压气机、燃烧室、燃气涡轮、动力涡轮及尾喷口六大部件的涡轴发动机数学模型，而后按进气流程进行各部件的气动热力计算，该部件法建模业内比较成熟，在此不加详述。

(b)解算发动机稳态平衡方程组：涡轴发动机功率解算模型也应该满足各转动部件的功率平衡，各截面的质量流量连续条件，这样就得到如下五个平衡方程：

$\left.\begin{array}{c}E1=\frac{{\stackrel{\‾}{P}}_g}{{\stackrel{\‾}{P}}_c/{\mathrm{\η}}_c+{\stackrel{\‾}{P}}_{\mathrm{ext}}}-1=0\\ E2=\frac{{\stackrel{\‾}{P}}_p}{{\stackrel{\‾}{P}}_l/{\mathrm{\η}}_l}-1=0\\ E3=\frac{W_{41\mathrm{cX}}-W_{41c}}{W_{41c}}=0\\ E4=\frac{W_{44\mathrm{cX}}-W_{44c}}{W_{44c}}=0\\ E5=\frac{P_8^{\′}}{P_8}-1=0\end{array}\right\}---\left(2\right)$

其中，E1～E5分别表示压气机与燃气涡轮功率平衡、负载与动力涡轮功率平衡、燃气涡轮进口流量连续、动力涡轮进口流量连续、尾喷口压力平衡，

分别为旋翼负载需求功率、压气机功率、燃气涡轮功率、动力涡轮功率、附件功率，η_l，η_c分别为旋翼负载和压气机效率，W_41c，W_44c分别为燃气涡轮和动力涡轮进口换算流量，下标X表示初猜值，P₈，P′₈分别为尾喷口总压和环境压力。选压气机相对转速

动力涡轮相对转速

燃气涡轮压比系数Z_G、动力涡轮压比系数Z_P、压气机压比系数Z_C为猜值参数，在动态时选Z_G，Z_P，Z_C为猜值参数。

上述发动机稳态平衡方程组(2)可采用现有的多种方法解算，本具体实施方式中采用L-M算法进行解算，具体如下：

Step1：选取β∈(0，1)，初始参数μ₀，增长因子v＞1；

Step2：给定初猜值x₀，设置迭代次数k＝0，其中x为发动机模型猜值参数；

Step3：计算f_k＝f(x_k)，

其中f为发动机稳态平衡方程组(2)；

Step4：计算A_k＝Df(x_k)，

其中D为雅克比矩阵；

Step5：计算 $g_k=A_k^T{f}_k;$

Step6：解方程(Q_k+μ_kI)p_k＝-g_k，其中p_k为当前下降方向，μ_k为修正实因子

Step7：计算x_k+1＝x_k+p_k(μ_k)，f_k+1＝f(x_k)，

Step8：若满足H-收敛准则(或Himmelblau终止准则)，则输出x_k+1，S_k+1；

Step9：检验：

若成立，则令μ_k+1＜μ_k/y，转Step10，否则令μ_k+1＝vμ_k，转Step6；

Step10：k＝k+1，转Step3。

本发明的直升机/涡轴发动机动态计算包括直升机的动态计算及发动机的动态计算两个部分，参照附图4所示的直升机/发动机动态计算流程图，其实现过程如下：

步骤C1、直升机的动态计算；

直升机的动态计算就是稳态配平后，进行动态仿真，根据主旋翼、机身、尾翼、尾桨等各部件力和力矩的改变求解直升机运动微分方程，解算出直升机的飞行姿态(姿态角)、飞行状态(高度、速度，加速度等)变化，其实现过程如下：

(a)直升机各部件上的力和力矩在体轴系的重心处叠加得到总外力和总外力矩，将上述各力和力矩引入直升机通用运动方程(欧拉方程)，得到描述直升机运动的微分方程组：

${\stackrel{\·}{V}}_x=(X_{\mathrm{SUM}}-G{\sin \mathrm{\θ}}_b)g/G+{\mathrm{rV}}_y-q{V}_z$

${\stackrel{\·}{V}}_y=(Y_{\mathrm{SUM}}-G\cos {\mathrm{\θ}}_b\sin {\mathrm{\φ}}_b)g/G+p{V}_z-r{V}_x$

${\stackrel{\·}{V}}_z=(Z_{\mathrm{SUM}}+G\cos {\mathrm{\θ}}_b\cos {\mathrm{\φ}}_b)g/G+q{V}_x-p{V}_y$

$\stackrel{\·}{p}=\frac{I_Z}{I_X{I}_X-I_{\mathrm{XZ}}^2}\{L_{\mathrm{SUM}}-(I_Z-I_Y)\mathrm{qr}+I_{\mathrm{XZ}}\mathrm{pq}\}+\frac{I_{\mathrm{XZ}}}{I_X{I}_Z-I_{\mathrm{XZ}}^2}\{N_{\mathrm{SUM}}-(I_Y-I_X)\mathrm{pq}-I_{\mathrm{XZ}}\mathrm{rq}\}$ (3)

$\stackrel{\·}{q}=\frac{1}{I_Y}\{M_{\mathrm{SUM}}-(I_X-I_Z)\mathrm{pr}+I_{\mathrm{XZ}}(r^2-p^2)\}$

$\stackrel{\·}{r}=\frac{I_X}{I_X{I}_Z-I_{\mathrm{XZ}}^2}\{N_{\mathrm{SUM}}-(I_Y-I_X)\mathrm{pq}-I_{\mathrm{XZ}}\mathrm{rq}\}+\frac{I_{\mathrm{XZ}}}{I_X{I}_Z-I_{\mathrm{XZ}}^2}\{L_{\mathrm{SUM}}-(I_Z-I_Y)\mathrm{qr}+I_{\mathrm{XZ}}\mathrm{pq}\}$

${\stackrel{\·}{\mathrm{\θ}}}_b=57.3(q{\cos \mathrm{\φ}}_b-r\sin {\mathrm{\φ}}_b)$

${\stackrel{\·}{\mathrm{\ψ}}}_b=57.3\frac{r\cos {\mathrm{\φ}}_b+q\sin {\mathrm{\φ}}_b}{\cos {\mathrm{\θ}}_b}$

${\stackrel{\·}{\mathrm{\φ}}}_b=57.3p+{\stackrel{\·}{\mathrm{\ψ}}}_b\sin {\mathrm{\θ}}_b$

其中， $V_x=\∫{\stackrel{\·}{V}}_x\mathrm{dt},$ $V_y=\∫{\stackrel{\·}{V}}_y\mathrm{dt},$ $V_z={\∫\stackrel{\·}{V}}_z\mathrm{dt},$ $p=\∫\stackrel{\·}{p}\mathrm{dt},$ $q=\∫\stackrel{\·}{q}\mathrm{dt},$ $r=\∫\stackrel{\·}{r}\mathrm{dt}$ 分别为直升机体轴系下速度和角速度， ${\mathrm{\θ}}_b=\∫{\stackrel{\·}{\mathrm{\θ}}}_b\mathrm{dt},$ ${\mathrm{\φ}}_b={\∫\stackrel{\·}{\mathrm{\φ}}}_b\mathrm{dt},$ ${\mathrm{\ψ}}_b={\∫\stackrel{\·}{\mathrm{\ψ}}}_b\mathrm{dt}$ 分别为机体的俯仰角、滚转角、偏航角，g为重力加速度，G为直升机重量，X_SUM，Y_SUM，Z_SUM，L_SUM，M_SUM，N_SUM分别为直升机体轴系下的合力和合力矩，I_X，I_Y，I_Z，I_XY，I_XZ，I_YZ分别为机体惯性矩和惯性积。

(b)对上述直升机运动微分方程组(3)的求解，可采用多种方法，本具体实施方式中采用四阶的龙格-库塔法进行解算，具体如下：

具体的龙格-库塔法已比较成熟，在此不加详述。本具体实施方式中采用四阶的龙格-库塔法是在一阶欧拉法基础上进行改进，由于超过四阶的龙格-库塔公式效率并不高，故最终选择四阶，引入四个点的斜率值，以减少累积误差，提高仿真精度。本具体实施方式中在选取直升机仿真步长h时，结合了直升机系统本身的特性，采用调试法，最后使得累积误差和截断舍入误差达到一个均衡，最终确定h＝5ms。

步骤C2、发动机的动态计算；

发动机的动态计算就是在直升机动态计算后，根据直升机功率变化，计算出所引起的发动机动态响应，即根据直升机飞行状态及需用功率的改变建立涡轴发动机动态平衡方程组并进行解算，其实现过程如下：

(a)发动机动态平衡方程组：动态过程中，由于发动机工作在非平衡状态，各截面的流量仍满足连续条件，但燃气涡轮和动力涡轮的转子功率不再平衡，这样根据式(2)就可得到如下的3个平衡方程：

$\left.\begin{array}{c}E3=\frac{W_{41\mathrm{cX}}-W_{41c}}{W_{41c}}=0\\ E4=\frac{W_{44\mathrm{cX}}-W_{44c}}{W_{44c}}=0\\ E5=\frac{P_8^{\′}}{P_8}-1=0\end{array}\right\}---\left(3\right)$

其中参数解释参见公式(2)，而公式(2)中相应的功率平衡方程由如下转子动力学方程取代：

$\left.\begin{array}{c}{m}_g*{\mathrm{\η}}_g-m_c=J_g*\frac{\mathrm{\π}}{30}*\frac{{\mathrm{dn}}_g}{\mathrm{dt}}\\ m_p*{\mathrm{\η}}_p-m_l=J_p*\frac{\mathrm{\π}}{30}*\frac{{\mathrm{dn}}_p}{\mathrm{dt}}\end{array}\right\}---\left(4\right)$

其中，J_g，J_p分别表示燃气涡轮、动力涡轮的转子转动惯量，η_g，η_p表示燃气涡轮、功率涡轮的转子机械效率，n_g，n_p燃气涡轮和动力涡轮转子转速，m_g，m_c，m_p，m_l表示燃气涡轮、压气机、动力涡轮、旋翼负载扭矩。将功率与扭矩的关系代入，可得转加速度与功率的关系：

$\left.\begin{array}{c}\frac{{\mathrm{dn}}_g}{\mathrm{dt}}=(W_g*{\mathrm{\η}}_g-W_c)/[n_g*J_g*{\left(\frac{\mathrm{\π}}{30}\right)}^2]\\ \frac{{\mathrm{dn}}_p}{\mathrm{dt}}=(W_p*{\mathrm{\η}}_p-W_l)/[n_p*J_p*{\left(\frac{\mathrm{\π}}{30}\right)}^2]\end{array}\right\}---\left(5\right)$

其中，W_g，W_p分别表示燃气涡轮和动力涡轮功率。

(b)发动机动态平衡方程组公式(3)和公式(5)可采用多种方法解算，本具体实施方式中采用一次通过算法进行解算，具体如下：

首先根据飞行马赫数、飞行高度、主燃油量进行发动机各部件模型的计算，而后根据公式(5)转子动力学方程计算燃气涡轮及动力涡轮转子的加速度，结合上一时刻的转速求得当前时刻转速，并求得公式(3)动态平衡方程E3、E4和E5的值，若其值不在允许范围时，则通过一次通过算法迭代修正初猜值，得到该时刻发动机的各截面的参数值，如此反复计算下去，直到转子的加速度为零，发动机再次进入稳定状态。其中，一次通过算法是在Newton-Raphson算法上进行改进，不同之处在于两次相邻计算时刻之间不再进行反复迭代，而只用雅可比矩阵对初猜值进行一次修正，便进入下一个时刻的计算，几乎不需要进行迭代，计算速度提高了很多，实时性有明显改善。

本发明的直升机/涡轴发动机实际运行状态解算方法，包括直升机实际运行解算方法以及涡轴发动机实际运行状态解算方法两部分，其实现过程如下：

直升机实际运行解算方法设计：

直升机实际运行解算方法采用修正的线性二次型调解器算法，该算法是基于状态变量小偏差模型，具体实现步骤如下：

(a)求取直升机状态变量模型：采用最小二乘拟合法求取直升机9阶状态变量模型，由于最小二乘拟合法是控制领域熟知的方法，在此不作介绍。在设计过程中应用的名义对象是基于考虑了动力入流的六自由度直升机非线性模型的海平面悬停点小偏差线化模型，其状态空间表达式为：

$\left\{\begin{array}{c}\stackrel{\·}{x}=\mathrm{Ax}+\mathrm{Bu}\\ y=\mathrm{Cx}+\mathrm{Du}\end{array}\right.---\left(6\right)$

其中，x为直升机体轴系下三个方向的速度、角速度以及机体的俯仰角、滚转角、偏航角，共9个状态量；u为总矩、横向周期变距、纵向周期变距和尾桨总矩4个控制输入；y为直升机三个方向的速度及偏航角共4个轨迹输出量，其中各量均为偏差量，省去了偏差符号Δ；A为9X9矩阵，B为9X4矩阵，将二者总共117个元素设定成为决策向量；C为4X9矩阵，D为4X4矩阵。

(b)修正的线性二次型调解器算法设计：对于直升机系统即式(6)，设直升机的参考指令为常值信号r，则误差为e＝r-y，对误差信号微分可得：

$\left[\begin{array}{c}\stackrel{\·\·}{x}\\ \stackrel{\·}{e}\end{array}\right]=\left[\begin{array}{c}\stackrel{\·\·}{x}\\ \stackrel{\·}{r}-\stackrel{\·}{y}\end{array}\right]=\left[\begin{array}{cc}A& 0\\ -C& 0\end{array}\right]\left[\begin{array}{c}\stackrel{\·}{x}\\ e\end{array}\right]+\left[\begin{array}{c}B\\ -D\end{array}\right]\stackrel{\·}{u}---\left(7\right)$

这里取增广状态量为： $\stackrel{\‾}{x}={\left[\begin{array}{cc}{\stackrel{\·}{x}}^T& e^T\end{array}\right]}^T,$ 控制量为：

则其中 $\stackrel{\‾}{A}=\left[\begin{array}{cc}A& 0\\ -C& 0\end{array}\right],$ $\stackrel{\‾}{B}=\left[\begin{array}{c}B\\ -D\end{array}\right],$ 对该增广系统可以设计线性二次型调节器，使其所有状态变量都保持为零： $\stackrel{\‾}{x}={\left[\begin{array}{cc}{\stackrel{\·}{x}}^T& e^T\end{array}\right]}^T=0,$ 亦即系统式(6)的状态保持不变，且参考指令与系统输出之差e为零，这就满足了消除稳态误差的要求。而后用LQR方法求出的解算公式为将按

和e表示为分块矩阵并取拉氏变换，可得系统式(6)的解算公式为：

$u=-K_{\stackrel{\·}{x}}x-\frac{K_{e}e}{s}---\left(8\right)$

其中，K_x阵负责原始状态量的解耦和状态调节增稳，K_e负责实现前飞、侧飞、高度、航向通道的指令跟踪控制。到此，便得到了系统式(6)的增广LQR解算方法，最终，在直升机设计点得到的解算相关参数如下：

$K_x=\left[\begin{array}{ccccccccc}-0.4001& -0.0467& -1.1564& -0.0797& 1.0657& 1.8002& -0.3976& 1.4688& 4.2464\\ -0.4274& 1.6979& 0.0909& 2.4088& 2.6660& 1.9653& 13.4524& 0.8429& 5.1153\\ 2.3167& 0.2178& 0.0080& 0.3494& -13.4887& 1.9174& 2.1999& 1.0919& -26.3343\\ 0.4334& 0.1857& -0.8609& 0.2838& -1.3497& -7.4928& 0.7367& -4.4255& -4.7084\end{array}\right]$ $K_e=\left[\begin{array}{cccc}0.1333& 0.0479& 0.8483& -0.4623\\ 0.1608& -0.9210& -0.0461& -0.1782\\ -0.9102& -0.1216& 0.0390& -0.2712\\ -0.1668& -0.1864& 0.4203& 0.9386\end{array}\right]$

在上述修正线性二次型调节器算法解算过程中，选取了R＝100*diag([1，1，1，1.5])，Q＝100*diag([1，1，1，1，1，1，1，1，1，0.9，0.9，0.9，1.5])，其中对角阵Q前面9个元素分别用于调节与其对应的状态变量的作用程度，新增广的后四个元素分别表征了控制器作为指令跟踪器使用时对于指令偏差的敏感程度。

发动机实际运行解算方法设计：

涡轴发动机实际运行解算方法采用串级PID+扭矩前馈+总距补偿方法，并接入扭振滤波器进行负载滤波，其实现步骤如下：

(a)串级PID+扭矩前馈+总距补偿解算方法：国内通常采用总距前馈+串级PID的方法，虽然总距前馈一定程度上反映了直升机功率需求的变化，但是由于直升机功率需求是飞行状态(飞行高度、飞行速度)与操纵量(总距、纵横周期变距以及尾桨总距)的复杂组合，因此总距并不能很好反映直升机当前的飞行状况或功率需求，而扭矩信息直接反映了系统不同条件下的负载变化。直升机扭矩负载的变化会引起涡轴发动机动力涡轮转速的强扰动，采用扭矩前馈+总距补偿的方法可以有效地抑制发动机子系统的强扰动，提高直升机飞行操纵品质。串级PID+扭矩前馈+总距补偿解算方法可如下表达：

$\mathrm{\Δ}{W}_f={\mathrm{\Σ}}_{m=1}^2(K_{\mathrm{pm}}{e}_m+K_{\mathrm{im}}\mathrm{\Σ}{e}_m)+K_1\mathrm{\ΔT}+K_2\mathrm{\Δ}{\mathrm{\θ}}_0---\left(9\right)$

其中，K_p1＝0.3，K_i1＝0.29为串级PID外环参数，K_p2＝0.002，K_i2＝0.000为内环参数，e₁为动力涡轮转速偏差量，e₂为燃气涡轮转速偏差量，K₁、K₂为可调系数，ΔT为扭矩偏差量，Δθ₀为总距偏差量。K₂根据总距杆角度的大小插值得到，K₁为一给定的参数。

(b)扭振滤波器设计：由于旋翼本身特性，旋翼需用功率是一个周期振荡的载荷，并且随着飞行速度的增加，其振荡幅值也随之增加，这就对发动机控制器的带宽提出了非常高的要求。本具体实施方式中，提出了一种基于带通滤波和均值滤波的扭振滤波器设计方法，参照附图5所示的滤波器结构图，其结构简单并具有很好的精度，能够有效对旋翼负载进行滤波以降低对发动机控制器带宽的要求，其具体形式如下：

带通滤波器的表达形式为：

$G\left(s\right)=(\frac{1}{{\mathrm{\ω}}^2}{s}^2+\frac{T_1}{\mathrm{\ω}}s+1)/(\frac{1}{{\left(\mathrm{k\ω}\right)}^2}{s}^2+\frac{T_2}{\mathrm{\ω}}s+1)---\left(10\right)$

Y(s)＝G(s)·X(s)

均值滤波器的表达形式为：

$Y_n=\frac{1}{N}\underset{i=0}{\overset{N}{\mathrm{\Σ}}}{X}_i---\left(11\right)$

其中，N是一个振荡周期内数值仿真计算的次数。该滤波器的准稳态判断依据是当前振荡周期中带通滤波器的输出均值与上一个振荡周期的输出均值之差小于ε，ε为事先设定的较小正数。

最后，通过数字仿真检验本发明的一种单旋翼直升机/涡轴发动机通用仿真模型的稳态配平与动态计算的精度，并通过半物理试验检验本发明方法的可行性，如下：

(a)直升机/发动机稳态配平效果：如附图6所示，以飞行高度ΔH＝600m，前飞速度ΔVx＝5m/s划分飞行包线，对每个包线点进行直升机稳态配平解算，经过与UH-60A直升机飞行试验数据对比分析，关键参数曲线基本一致，反映出了本文所采用的建模机理是合理的。其中H为飞行高度，图(a)为总距θ₀配平图，图(b)为横向周期变距A_1C配平图，图(c)为纵向周期变距B_1S配平图，图(d)为尾桨总距θT配平图，图(e)为俯仰角θ配平图，图(f)为滚转角□配平图。

(b)直升机/发动机动态仿真效果：在飞行高度H＝3km，前飞速度Vx＝30m/s状态下，将主旋翼轴转速限定在27rad/s对直升机动态模型进行开环校验，如附图7所示为直升机四个主控通道的响应情况，响应曲线仅画出前飞速度Vx、侧飞速度Vy、爬升速度Vz和偏航角ψ，从响应曲线可以看出直升机开环操纵响应结果基本与UH-60A试飞结果吻合。其中，图(a)为总距θ₀+1°阶跃响应图，图(b)为横向周期变距A_1C+0.5°阶跃响应图，图(c)为纵向周期变距B_1S+0.5°阶跃响应图，图(d)为尾桨总矩θ_T+1°阶跃响应图。

(c)半物理仿真验证效果：对照附图8的直升机/发动机综合系统半物理仿真平台结构图，将黑鹰UH-60A直升机/T700涡轴发动机综合模型与半物理仿真平台无缝对接，仿真步长为20ms，发动机控制器采用与实物电子控制器性能相当的工业控制计算机，控制算法中加入了信号的噪声及滞后以模拟真实信号，信号的输入输出均通过仿真接口模拟器经AD/DA数模转换得到，采用多线程程序设计方法设计系统的串口通讯，保证信号传输的实时性与同步性。如附图9所示：(1-2)阶段模拟发动机地面起动；(2-3)阶段模拟发动机从起动到慢车；(3-4)阶段模拟发动机从地面慢车到空载飞行；(4-5)阶段通过推拉控制台上的负载杆，对发动机模型施加负载扭矩，当负载扭矩超过某一值后接通直升机模型；(5-6)阶段接通直升机模型后，通过直升机模型界面给定飞行操作指令，使直升机模型以爬升速度Vz＝6m/s爬升到高度H＝2000m；(6-7)阶段直升机模型爬升到H＝2000m后爬升速度Vz降为0m/s；(7-8)阶段直升机模型以前飞速度Vx＝75m/s前飞；(8-9)阶段直升机模型前飞速度Vx到达75m/s后，等待10s，以爬升速度Vz＝3m/s爬升到高度H＝2500m；(9-10)阶段直升机模型爬升到H＝2500m后，爬升速度Vz降到0m/s，前飞速度Vx降到0m/s；(10-11)阶段直升机模型以爬升速度Vz＝-4m/s速度下降，降到高度H＝0m处悬停；(11-12)阶段模拟发动机慢车；(12之后)阶段模拟发动机停车。附图中动力涡轮转速PNP、燃气涡轮转速PNC、燃油流量WFB、动力涡轮功率即发动机输出功率HPP、前飞速度Vx、飞行高度H及旋翼扭矩Trotor的动态响应曲线可综合看出，根据本发明方法构建的单旋翼直升机/涡轴发动机通用仿真模型具有可靠的置信度。

Claims (9)

1.一种单旋翼直升机/涡轴发动机通用综合模型的构建方法，其特征在于，包括以下步骤：

步骤A、分别建立直升机和发动机的非线性模型；

步骤B、单旋翼直升机/涡轴发动机的稳态配平，包括旋翼动力学系统配平、机体动力学模型配平、发动机模型配平；

步骤C、单旋翼直升机/涡轴发动机的动态计算；

步骤D、单旋翼直升机/涡轴发动机的实际运行解算。

2.如权利要求1所述单旋翼直升机/涡轴发动机通用综合模型的构建方法，其特征在于，所述步骤B具体包括以下步骤：

步骤B1、旋翼动力学系统配平，具体按照以下方法：首先计算每个翼型的迎角，其中旋翼诱导速度采用非均匀一阶谐波入流的计算方法，桨叶挥舞运动则基于桨叶的气流环境分析旋翼桨叶的非定常挥舞运动，以递推形式得到挥舞运动方程，然后根据翼型迎角计算其二元翼型气动力，再通过积分求解每片桨叶的气动力，最终计算出旋翼的气动力，完成直升机旋翼模型的稳态计算；

步骤B2、机体动力学模型配平，具体按照以下方法：进行直升机旋翼、机身、尾桨、水平尾翼和垂直尾翼的动力学解算，分别计算出各部件在直升机重心处产生的力和力矩；根据计算出的各部件在直升机重心处产生的力和力矩，以及直升机重力，解算如下机体平衡方程组：

$\left.\begin{array}{c}{X}_{\mathrm{SUM}}-G{\sin \mathrm{\θ}}_b=0\\ Y_{\mathrm{SUM}}-{G\cos \mathrm{\θ}}_b\sin {\mathrm{\φ}}_b=0\\ Z_{\mathrm{SUM}}+G\cos {\mathrm{\θ}}_b\cos {\mathrm{\φ}}_b=0\\ L_{\mathrm{SUM}}=0\\ M_{\mathrm{SUM}}=0\\ N_{\mathrm{SUM}}=0\end{array}\right\},$

其中，G表示直升机重量，θ_b表示机体的俯仰角，φ_b表示机体的滚转角，下标b表示机体参数，X_SUM＝X_MR+X_WF+X_TR+X_HT+X_VT、Y_SUM＝Y_MR+Y_WF+Y_TR+Y_HT+Y_VT与Z_SUM＝Z_MR+Z_WF+Z_TR+Z_HT+Z_VT分别表示直升机体轴系下X/Y/Z方向的合力，L_SUM＝L_MR+L_WF+L_TR+L_HT+L_VT、M_SUM＝M_MR+M_WF+M_TR+M_HT+M_VT与N_SUM＝N_MR+N_WF+N_TR+N_HT+N_VT分别为其对应的合力矩，下标SUM表示合力或合力矩，MR表示旋翼参量，WF表示机身参量，TR表示尾桨参量，HT表示水平尾翼参量，VT表示垂直尾翼参量；

步骤B3、发动机模型配平，具体按照以下方法：首先进行涡轴发动机气动热力计算；然后根据气动热力计算结果解算如下的发动机稳态平衡方程组：

$\left.\begin{array}{c}E1=\frac{{\stackrel{\‾}{P}}_g}{{\stackrel{\‾}{P}}_c/{\mathrm{\η}}_c+{\stackrel{\‾}{P}}_{\mathrm{ext}}}-1=0\\ E2=\frac{{\stackrel{\‾}{P}}_p}{{\stackrel{\‾}{P}}_l/{\mathrm{\η}}_l}-1=0\\ E3=\frac{W_{41\mathrm{cX}}-W_{41c}}{W_{41c}}=0\\ E4=\frac{W_{44\mathrm{cX}}-W_{44c}}{W_{44c}}=0\\ E5=\frac{P_8^{\′}}{P_8}-1=0\end{array}\right\},$

其中，E1～E5分别表示压气机与燃气涡轮功率平衡、负载与动力涡轮功率平衡、燃气涡轮进口流量连续、动力涡轮进口流量连续、尾喷口压力平衡，

分别为旋翼负载需求功率、压气机功率、燃气涡轮功率、动力涡轮功率、附件功率，η_l，η_c分别为旋翼负载和压气机效率，W_41c，W_44c分别为燃气涡轮和动力涡轮进口换算流量，下标X表示初猜值，P₈，P′₈分别为尾喷口总压和环境压力，压气机相对转速

动力涡轮相对转速

燃气涡轮压比系数Z_G、动力涡轮压比系数Z_P、压气机压比系数Z_C为猜值参数。

3.如权利要求2所述单旋翼直升机/涡轴发动机通用综合模型的构建方法，其特征在于，所述机体平衡方程组和发动机稳态平衡方程组的解算均采用Levenberg-Marquarat算法。

4.如权利要求2所述单旋翼直升机/涡轴发动机通用综合模型的构建方法，其特征在于，所述单旋翼直升机/涡轴发动机的动态计算包括直升机动态计算和涡轴发动机动态计算，其具体方法分别如下：

直升机动态计算：

根据直升机各部件上的力和力矩，解算如下的直升机运动微分方程组：

${\stackrel{\·}{V}}_x=(X_{\mathrm{SUM}}-G{\sin \mathrm{\θ}}_b)g/G+{\mathrm{rV}}_y-q{V}_z$

${\stackrel{\·}{V}}_y=(Y_{\mathrm{SUM}}-G\cos {\mathrm{\θ}}_b\sin {\mathrm{\φ}}_b)g/G+p{V}_z-r{V}_x$

${\stackrel{\·}{V}}_z=(Z_{\mathrm{SUM}}+G\cos {\mathrm{\θ}}_b\cos {\mathrm{\φ}}_b)g/G+q{V}_x-p{V}_y$

$\stackrel{\·}{p}=\frac{I_Z}{I_X{I}_X-I_{\mathrm{XZ}}^2}\{L_{\mathrm{SUM}}-(I_Z-I_Y)\mathrm{qr}+I_{\mathrm{XZ}}\mathrm{pq}\}+\frac{I_{\mathrm{XZ}}}{I_X{I}_Z-I_{\mathrm{XZ}}^2}\{N_{\mathrm{SUM}}-(I_Y-I_X)\mathrm{pq}-I_{\mathrm{XZ}}\mathrm{rq}\}$

$\stackrel{\·}{q}=\frac{1}{I_Y}\{M_{\mathrm{SUM}}-(I_X-I_Z)\mathrm{pr}+I_{\mathrm{XZ}}(r^2-p^2)\},$

$\stackrel{\·}{r}=\frac{I_X}{I_X{I}_Z-I_{\mathrm{XZ}}^2}\{N_{\mathrm{SUM}}-(I_Y-I_X)\mathrm{pq}-I_{\mathrm{XZ}}\mathrm{rq}\}+\frac{I_{\mathrm{XZ}}}{I_X{I}_Z-I_{\mathrm{XZ}}^2}\{L_{\mathrm{SUM}}-(I_Z-I_Y)\mathrm{qr}+I_{\mathrm{XZ}}\mathrm{pq}\}$

${\stackrel{\·}{\mathrm{\θ}}}_b=57.3(q{\cos \mathrm{\φ}}_b-r\sin {\mathrm{\φ}}_b)$

${\stackrel{\·}{\mathrm{\ψ}}}_b=57.3\frac{r\cos {\mathrm{\φ}}_b+q\sin {\mathrm{\φ}}_b}{\cos {\mathrm{\θ}}_b}$

${\stackrel{\·}{\mathrm{\φ}}}_b=57.3p+{\stackrel{\·}{\mathrm{\ψ}}}_b\sin {\mathrm{\θ}}_b$

其中， $V_x=\∫{\stackrel{\·}{V}}_x\mathrm{dt},$ $V_y=\∫{\stackrel{\·}{V}}_y\mathrm{dt},$ $V_z={\∫\stackrel{\·}{V}}_z\mathrm{dt},$ $p=\∫\stackrel{\·}{p}\mathrm{dt},$ $q=\∫\stackrel{\·}{q}\mathrm{dt},$ $r=\∫\stackrel{\·}{r}\mathrm{dt}$ 分别为直升机体轴系下速度和角速度， ${\mathrm{\θ}}_b=\∫{\stackrel{\·}{\mathrm{\θ}}}_b\mathrm{dt},$ ${\mathrm{\φ}}_b={\∫\stackrel{\·}{\mathrm{\φ}}}_b\mathrm{dt},$ ${\mathrm{\ψ}}_b={\∫\stackrel{\·}{\mathrm{\ψ}}}_b\mathrm{dt}$ 分别为机体的俯仰角、滚转角、偏航角，g为重力加速度，G为直升机重量，X_SUM，Y_SUM，Z_SUM，L_SUM，M_SUM，N_SUM分别为直升机体轴系下的合力和合力矩，I_X，I_Y，I_Z，I_XY，I_XZ，I_YZ分别为机体惯性矩和惯性积；涡轴发动机动态计算：

解算如下发动机动态平衡方程组：

$\left.\begin{array}{c}E3=\frac{W_{41\mathrm{cX}}-W_{41c}}{W_{41c}}=0\\ E4=\frac{W_{44\mathrm{cX}}-W_{44c}}{W_{44c}}=0\\ E5=\frac{P_8^{\′}}{P_8}-1=0\\ \frac{{\mathrm{dn}}_g}{\mathrm{dt}}=(W_g*{\mathrm{\η}}_g-W_c)/[n_g*J_g*{\left(\frac{\mathrm{\π}}{30}\right)}^2]\\ \frac{{\mathrm{dn}}_p}{\mathrm{dt}}=(W_p*{\mathrm{\η}}_p-W_l)/[n_p*J_p*{\left(\frac{\mathrm{\π}}{30}\right)}^2]\end{array}\right\},$

其中，E3～E5分别表示燃气涡轮进口流量连续、动力涡轮进口流量连续、尾喷口压力平衡，W_41c，W_44c分别为燃气涡轮和动力涡轮进口换算流量，下标X表示初猜值，P₈，P′₈分别为尾喷口总压和环境压力；燃气涡轮压比系数Z_G、动力涡轮压比系数Z_P、压气机压比系数Z_C为猜值参数，J_g，J_p分别表示燃气涡轮、动力涡轮的转子转动惯量，η_g，η_p表示燃气涡轮、功率涡轮的转子机械效率，n_g，n_p燃气涡轮和动力涡轮转子转速，m_g，m_c，m_p，m_l表示燃气涡轮、压气机、动力涡轮、旋翼负载扭矩，W_g，W_p分别表示燃气涡轮和动力涡轮功率。

5.如权利要求4所述单旋翼直升机/涡轴发动机通用综合模型的构建方法，其特征在于，所述直升机运动微分方程组的解算采用龙格-库塔法。

6.如权利要求5所述单旋翼直升机/涡轴发动机通用综合模型的构建方法，其特征在于，所述直升机运动微分方程组的解算采用四阶龙格-库塔法。

7.如权利要求3所述单旋翼直升机/涡轴发动机通用综合模型的构建方法，其特征在于，所述发动机动态平衡方程组的解算采用一次通过算法。

8.如权利要求3所述单旋翼直升机/涡轴发动机通用综合模型的构建方法，其特征在于，所述单旋翼直升机/涡轴发动机的实际运行解算包括直升机实际运行解算和发动机实际运行解算，其具体方法分别如下：

直升机实际运行解算：

采用修正的线性二次型调节器方法，具体如下：

直升机的参考指令为常值信号r，则误差为e＝r-y，对误差信号微分可得增广系统：

$\left[\begin{array}{c}\stackrel{\·\·}{x}\\ \stackrel{\·}{e}\end{array}\right]=\left[\begin{array}{c}\stackrel{\·\·}{x}\\ \stackrel{\·}{r}-\stackrel{\·}{y}\end{array}\right]=\left[\begin{array}{cc}A& 0\\ -C& 0\end{array}\right]\left[\begin{array}{c}\stackrel{\·}{x}\\ e\end{array}\right]+\left[\begin{array}{c}B\\ -D\end{array}\right]\stackrel{\·}{u}$

其中，x为直升机体轴系下三个方向的速度、角速度以及机体的俯仰角、滚转角、偏航角，共9个状态量；u为总矩、横向周期变距、纵向周期变距和尾桨总矩4个控制输入；A，B，C，D为适维系统状态矩阵；取增广状态量为： $\stackrel{\‾}{x}={\left[\begin{array}{cc}{\stackrel{\·}{x}}^T& e^T\end{array}\right]}^T,$ 控制量为：则 $\stackrel{\·}{\stackrel{\‾}{x}}=\stackrel{\‾}{A}\stackrel{\‾}{x}+\stackrel{\‾}{B}\stackrel{\‾}{u},$ 其中 $\stackrel{\‾}{A}=\left[\begin{array}{cc}A& 0\\ -C& 0\end{array}\right],$ $\stackrel{\‾}{B}=\left[\begin{array}{c}B\\ -D\end{array}\right];$ 对该增广系统进行线性二次型调节器设计，使其所有状态变量都保持为零： $\stackrel{\‾}{x}={\left[\begin{array}{cc}{\stackrel{\·}{x}}^T& e^T\end{array}\right]}^T=0,$ 且参考指令与系统输出之差e为零；而后用线性二次型调节器法求出的解算公式为

将

按

和e表示为分块矩阵并取拉氏变换，可得实际运行系统控制输入u解算公式为：

$u=-K_{\stackrel{\·}{x}}x-\frac{K_{e}e}{s};$

发动机实际运行解算：

解算如下公式并接入扭振滤波器进行负载滤波：

$\mathrm{\Δ}{W}_f={\mathrm{\Σ}}_{m=1}^2(K_{\mathrm{pm}}{e}_m+K_{\mathrm{im}}\mathrm{\Σ}{e}_m)+K_1\mathrm{\ΔT}+K_2\mathrm{\Δ}{\mathrm{\θ}}_0,$

其中，K_p1＝0.3，K_i1＝0.29为串级PID外环参数，K_p2＝0.002，K_i2＝0.000为内环参数，e₁为动力涡轮转速偏差量，e₂为燃气涡轮转速偏差量，ΔT为扭矩偏差量，Δθ₀为总距偏差量，K₁、K₂为可调系数，K₂根据总距杆角度的大小插值得到，K₁为一给定的参数。

9.如权利要求8所述单旋翼直升机/涡轴发动机通用综合模型的构建方法，其特征在于，所述扭振滤波器采用带通滤波与均值滤波串联的形式。

Images