一种超高速飞行器纵横向耦合运动稳定性分析方法与流程

1.本发明涉及空气动力学及飞行力学领域，尤其涉及一种超高速飞行器纵横向耦合运动稳定性分析方法。


背景技术：

2.超高速飞行器受机体构型及高空高速飞行环境的影响，其横侧向稳定性较差，在机动过程中易发生多通道耦合运动引起的失稳，例如2010年美国htv-2飞行器由于多通道耦合失稳导致试飞失败。在超高速飞行器进行俯仰快速机动时，俯仰通道与滚转通道间存在较强的惯性耦合及气动/运动耦合作用，使得俯仰运动对滚转通道产生干扰，从而影响滚转通道稳定性并可能导致飞行器失稳。但目前针对超高速飞行器纵/横向（即俯仰/滚转）耦合运动稳定性分析的研究较少。
3.传统的稳定性分析方法一般基于小扰动假设，把飞行器的运动分为基准运动和扰动运动。将飞行器的扰动运动用纵向（俯仰）和横航向（滚转及偏航）分开的模型分别描述，使描述飞行器运动的微分方程组线性化，并使用线性稳定性理论作为理论分析方法。但超高速飞行器在俯仰角大范围变化的过程中，俯仰通道与滚转通道存在较强的惯性耦合及气动/运动耦合作用，此时运动方程是非线性的，传统理论方法不再适用，因此需寻求新的稳定性分析方法。
4.除理论分析外，通常可基于cfd计算及风洞试验方法开展飞行器稳定性预测。但是针对超高速飞行器多通道耦合运动工况，直接开展cfd计算需要耗费大量计算资源，且计算周期较长，而进行风洞多自由度动态试验也存在技术复杂、成本高的问题。
5.为了满足超高速飞行器技术的工程应用需求，当前亟需发展一种超高速飞行器纵/横向耦合运动稳定性分析方法，以实现快速高效的稳定性预测。


技术实现要素：

6.本发明意在提供一种超高速飞行器纵横向耦合运动稳定性分析方法，解决了现有的超高速飞行器在多通道耦合失稳分析中效率低的问题。
7.为了实现上述目的，本发明提供的一种技术方案：一种超高速飞行器纵横向耦合运动稳定性分析方法，包括如下步骤：s1、针对需要研究的飞行器运动特性，给定来流工况及飞行器俯仰角运动范围，在俯仰角运动范围选取多个样本点，然后选取俯仰通道正弦振荡运动振幅a和平均俯仰角，从而建立俯仰通道运动方程；s2、通过cfd计算获得样本点处的滚转力矩静导数以及处的滚转力矩动导数；s3、通过多项式拟合获得滚转力矩静导数关于俯仰角的函数表达式，将代入获得函数表达式；
s4、将a、、、代入考虑惯性耦合及气动/运动耦合作用的滚转通道动力学方程，并将方程简化为mathieu方程形式的方程；s5、通过希尔无限行列式法求解方程的稳定界限、；s6、通过、求解临界频率、、、、，当强迫俯仰运动频率满足、、中任意一项时，滚转运动会在惯性耦合和气动/运动耦合的共同作用下发散；s7、通过选取不同的运动振幅a和平均俯仰角， 重复步骤s1-s6，即可获得飞行器在给定运动参数范围内的纵/横向耦合运动稳定性。
8.进一步的，步骤s1中俯仰通道运动方程如下：
ꢀꢀꢀꢀꢀꢀꢀꢀ
（1）其中为飞行器俯仰角，为平均俯仰角，a为振幅，ω为俯仰振荡角频率，为时间。
9.进一步的，步骤s2中的滚转力矩静导数以及滚转力矩动导数的计算方法如下：
ꢀꢀꢀꢀꢀꢀꢀꢀ
（2）
ꢀꢀꢀꢀꢀꢀꢀꢀ
（3）其中为来流动压，为飞行器参考面积，为飞行器参考长度，为俯仰角处的滚转静导数，为平均俯仰角处的滚转动导数，为绕飞行器体轴的转动惯量。
10.进一步的，步骤s3中的多项式拟合方法使用的拟合函数如下：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（4）其中、为拟合系数。
11.进一步的，步骤s4中的滚转通道动力学方程如下：
ꢀꢀꢀꢀꢀꢀꢀ
（5）其中为气动/运动耦合力矩项；为惯性耦合力矩项；为滚转阻尼力矩项；，、分别为绕飞行器体轴、的转动惯量。
12.进一步的，步骤s4中滚转通道动力学方程如下：
ꢀꢀꢀꢀꢀꢀꢀ
（6）其中
；；；
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（7）；
13.进一步的，步骤s5中通过希尔无限行列式法求解方程的稳定界限、的计算方法如下：针对=1附近的稳定界限，将方程以2为周期的解表示如下：
ꢀꢀꢀꢀꢀ
（8）式中为复傅里叶级数的各项系数，将式(8)代入式(6)，并考虑，则有:
ꢀꢀꢀ
（9）式(9)可简化为
ꢀꢀꢀꢀꢀ
（10）仅在的所有系数均为零时，式才能成立，即： （11）要使式构成的关于的无限方程组有非零解，则对应的系数矩阵行列式的值须为零，其中对于任意满足；系数矩阵行列式的形式为：
ꢀꢀꢀꢀꢀꢀꢀ
（12）
其中：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（13）
14.同样地针对附近的稳定界限，将方程以4为周期的解表示如下：
ꢀꢀꢀꢀꢀꢀꢀꢀ
（14）将式(14)代入式(6)则有：
ꢀꢀꢀꢀꢀꢀ
（15）
15.式(15)可简化为 （16）仅在的所有系数均为零时，式(9)才能成立，即：
ꢀꢀꢀꢀꢀ
（17）系数矩阵行列式的形式为：
ꢀꢀꢀ
（18）其中：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（19）
16.将式(12)和式(18)表示为如下形式的方程：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（20）把和看做常数，则式(20)为关于、的隐函数，通过matlab的“fimplicit”命令
求得在=1附近和=0.25附近的ε随λ的变化曲线、，即稳定性界限。
17.进一步的，步骤s6中通过、求解临界频率、、、、的计算方法如下：针对给定转动惯量参数及俯仰运动参数、来建立关于的函数，即：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（21）在稳定性图谱中绘制作为工况曲线，工况曲线依次穿过=0、=0.25和=1附近的稳定性界限形成5个交点，设交点的坐标为，由式(7)可得：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（22）从而求得稳定性临界点对应临界频率,临界频率按照对应的值从小到大的顺序依次定义为、、、、。
18.与现有技术方案相比，本发明的有益效果：本方案针对超高速飞行器纵/横向耦合运动稳定性分析问题，提出了一种构建包含惯性耦合及气动/运动耦合效应的动力学方程并基于mathieu方程稳定性理论进行稳定性分析的方法，能够在保证预测精度的同时降低计算用时及成本，从而实现超高速飞行器纵/横向耦合运动稳定性的定性、定量、快速、高效预测，为超高速飞行器相关理论研究及工程应用提供指导。
附图说明
19.图1为本实施例的流程图；图2为本实施例中超高速飞行器的示意图；图3为本实施例中超高速飞行器拟合函数曲线与cfd结果对比图；图4为本实施例中超高速飞行器工况稳定性图谱；图5为本实施例中超高速飞行器工况稳定性平面图与cfd结果对比图。
具体实施方式
20.下面结合附图和实施方式对本发明作进一步的详细说明。
21.实施例如图1所示，一种超高速飞行器纵横向耦合运动稳定性分析方法，首先需要构建包含惯性耦合及气动/运动耦合效应的动力学方程，然后基于mathieu方程稳定性理论进行稳定性。本实施例具体包括以下步骤：s1、本实施例以图2所示的超高速飞行器为对象研究其纵/横向耦合运动稳定性。给定来流工况为马赫数6，飞行高度28千米，俯仰角运动范围-5~15
°
，在俯仰角运动范围选取样本点。给定俯仰通道正弦振荡运动振幅a=10
°
，平均俯仰角=5
°
，从而建立俯仰通道运动方程。俯仰通道运动方程如下：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（1）
其中为飞行器俯仰角，为平均俯仰角，a为振幅，ω为俯仰振荡角频率，为时间。
22.s2、通过cfd计算获得样本点处的滚转力矩静导数以及处的滚转力矩动导数。其中滚转力矩静导数以及滚转力矩动导数的计算方法如下：
ꢀꢀꢀꢀꢀꢀ
（2）
ꢀꢀꢀꢀꢀ
（3）其中为来流动压，为飞行器参考面积，为飞行器参考长度，为俯仰角处的滚转静导数，为平均俯仰角处的滚转动导数，为绕飞行器体轴的转动惯量，体轴定义如图2所示。
23.s3、通过多项式拟合获得滚转力矩静导数关于俯仰角的函数表达式。通过最小二乘法拟合获得拟合函数曲线与cfd结果对比如图3所示，获得拟合系数，决定系数r2=0.998817，说明采用二阶多项式即可获得较好的拟合效果。将代入获得函数表达式。其中多项式拟合方法使用的拟合函数如下：
ꢀꢀꢀꢀꢀꢀꢀ
（4）其中、为拟合系数。
24.s4、将a、、、代入考虑惯性耦合及气动/运动耦合作用的滚转通道动力学方程，并将方程简化为mathieu方程形式的方程；其中滚转通道动力学方程如下：
ꢀꢀꢀꢀꢀꢀꢀ
（5）其中为气动/运动耦合力矩项；为惯性耦合力矩项；为滚转阻尼力矩项；，、分别为绕飞行器体轴、的转动惯量，体轴、定义如图2所示。
25.滚转通道动力学方程如下：
ꢀꢀꢀꢀꢀꢀ
（6）其中；；；
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（7）
；
26.s5、通过希尔无限行列式法求解方程的稳定界限、，在平面内绘制稳定界限，获得稳定性图谱如图4所示。其中稳定界限、的计算方法如下：如图4所示，由floquet理论可知，mathieu方程的稳定界限能够把平面划分为若干稳定区和不稳定区，从轴上=0处出发有一条稳定性界限，从=0.25和=1附近出发分别有两条稳定性界限。其中=0处出发的稳定性方程在稳定界限上有周期性解，其中从=0.25出发的稳定性界限上有周期为4的解，从=1出发的稳定性界限上有周期为2的解。周期性解可用复傅里叶级数表示。通过希尔无限行列式法能够获得mathieu方程的稳定界限。具体方法如下：针对=1附近的稳定界限，将方程以2为周期的解表示如下：
ꢀꢀꢀꢀꢀꢀꢀ
（8）式中为复傅里叶级数的各项系数，将式(8)代入式(6)，并考虑，则有:
ꢀꢀꢀꢀꢀꢀ
（9）式(9)可简化为
ꢀꢀ
（10）仅在的所有系数均为零时，式才能成立，即：
ꢀꢀ
（11）要使式（11）构成的关于的无限方程组有非零解，则对应的系数矩阵行列式的值须为零，其中对于任意满足；系数矩阵行列式的形式为：
ꢀꢀꢀꢀꢀꢀꢀꢀ
（12）其中：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（13）
27.同样地针对附近的稳定界限，将方程以4为周期的解表示如下：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（14）将式(14)代入式(6)则有：
ꢀꢀꢀꢀꢀꢀꢀ
（15）
28.式(15)可简化为 （16）仅在的所有系数均为零时，式(9)才能成立，即：
ꢀꢀꢀꢀꢀꢀꢀ
（17）系数矩阵行列式的形式为：
ꢀꢀꢀꢀ
（18）其中：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（19）
29.形如式(12)和式(18)的希尔行列式有无限多阶，即式中的n应趋于无穷大，但在实际求解过程中需取有限大的n阶对无限行列式进行截断，可将n取5。
30.将式(12)和式(18)表示为如下形式的方程：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（20）把和看做常数，则式(20)为关于、的隐函数，通过matlab的“fimplicit”命令求得在=1附近和=0.25附近的随的变化曲线、，即稳定性界限。
31.s6、通过、求解临界频率、、、、，当强迫俯仰运动频率满足、、中任意一项时，滚转运动会在惯性耦合和气动/运动耦合的共同作用下发散。
32.其中通过、求解临界频率、、、、的计算方法如下：针对给定转动惯量参数及俯仰运动参数、的特定工况，建立关于的函数，即：
ꢀꢀꢀꢀꢀꢀꢀꢀ
（21）如图4所示，在稳定性图谱中绘制作为工况曲线，工况曲线依次穿过=0、=0.25和=1附近的稳定性界限形成5个交点，设交点的坐标为，由式(7)可得：
ꢀꢀꢀꢀꢀꢀꢀꢀ
（22）从而求得稳定性临界点对应临界频率，临界频率按照对应的值从小到大的顺序依次定义为、、、、。工况曲线与稳定界限交点坐标及对应频率如下表1所示。
33.表1 工况曲线与稳定界限交点坐标及对应频率
34.s7、选取a=5
°
，=5
°
，重复步骤s1~s6，可获得该工况下的临界频率、、、、，基于a=10
°
与a=5
°
工况下的临界频率结果，可绘制稳定性平面图并与cfd结果对比如图5所示。图中阴影区域为理论分析预测的不稳定区，图中每一个点代表一个cfd工况结果，“·”表示该工况下cfd结果为滚转运动收敛，“*”表示该工况下cfd结果为滚转运动发散，“。”表示该工况下cfd结果为极限环运动，即临界稳定状态。图5中可以看到，位于阴影区域外部的工况数值模拟结果均为收敛，位于阴影区域内的工况数值模拟结果基本为发散，表明数值模拟结果与理论预测基本相符，这表明发明公开的稳定性预测方法在本文所研究的工况范围内是有效的，通过本方法可获得超高速飞行器在给定运动参数范围内的纵/横向耦合运动稳定性。
35.以上所述的仅是本发明的实施例，方案中公知的具体技术方案或特性等常识在此未作过多描述。应当指出，对于本领域的技术人员来说，在不脱离本发明技术方案的前提下，还可以作出若干变形和改进，这些也应该视为本发明的保护范围，这些都不会影响本发明实施的效果和专利的实用性。本技术要求的保护范围应当以其权利要求的内容为准，说明书中的具体实施方式等记载可以用于解释权利要求的内容。

技术特征：
1..一种超高速飞行器纵横向耦合运动稳定性分析方法，其特征在于：包括如下步骤：s1、针对需要研究的飞行器运动特性，给定来流工况及飞行器俯仰角运动范围，在俯仰角运动范围选取多个样本点，然后选取俯仰通道正弦振荡运动振幅a和平均俯仰角，从而建立俯仰通道运动方程；s2、通过cfd计算获得样本点处的滚转力矩静导数以及处的滚转力矩动导数；s3、通过多项式拟合获得滚转力矩静导数关于俯仰角的函数表达式，将代入获得函数表达式；s4、将a、、、代入考虑惯性耦合及气动/运动耦合作用的滚转通道动力学方程，并将方程简化为mathieu方程形式的方程；s5、通过希尔无限行列式法求解方程的稳定界限、；s6、通过、求解临界频率、、、、，当强迫俯仰运动频率满足、、中任意一项时，滚转运动会在惯性耦合和气动/运动耦合的共同作用下发散；s7、通过选取不同的运动振幅a和平均俯仰角， 重复步骤s1-s6，即可获得飞行器在给定运动参数范围内的纵/横向耦合运动稳定性。2.根据权利要求1所述的一种超高速飞行器纵横向耦合运动稳定性分析方法，其特征在于：步骤s1中俯仰通道运动方程如下：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（1）其中为飞行器俯仰角，为平均俯仰角，a为振幅，ω为俯仰振荡角频率，为时间。3. 根据权利要求1所述的一种超高速飞行器纵横向耦合运动稳定性分析方法，其特征在于：步骤s2中的滚转力矩静导数以及滚转力矩动导数的计算方法如下：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（2）
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（3）其中为来流动压，为飞行器参考面积，为飞行器参考长度，为俯仰角处的滚转静导数，为平均俯仰角处的滚转动导数，为绕飞行器体轴的转动惯量。4.根据权利要求1所述的一种超高速飞行器纵横向耦合运动稳定性分析方法，其特征在于：步骤s3中的多项式拟合方法使用的拟合函数如下：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（4）其中、为拟合系数。5.根据权利要求1所述的一种超高速飞行器纵横向耦合运动稳定性分析方法，其特征在于：步骤s4中的滚转通道动力学方程如下：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（5）其中为气动/运动耦合力矩项；为惯性耦合力矩项；为滚转阻尼力矩项；，、分别为绕飞行器体轴、的转动惯量。6.根据权利要求1所述的一种超高速飞行器纵横向耦合运动稳定性分析方法，其特征在于：步骤s4中滚转通道动力学方程如下：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（6）其中；；；
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（7）；。7.根据权利要求1所述的一种超高速飞行器纵横向耦合运动稳定性分析方法，其特征在于：步骤s5中通过希尔无限行列式法求解方程的稳定界限、的计算方法如下：针对=1附近的稳定界限，将方程以2为周期的解表示如下：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（8）式中为复傅里叶级数的各项系数，将式(8)代入式(6)，并考虑，则有:
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（9）式(9)可简化为
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（10）
仅在的所有系数均为零时，式才能成立，即：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（11）；要使式构成的关于的无限方程组有非零解，则对应的系数矩阵行列式的值须为零，其中对于任意满足；系数矩阵行列式的形式为：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（12）其中：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（13）；同样地针对附近的稳定界限，将方程以4为周期的解表示如下：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（14）将式(14)代入式(6)则有：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（15）；式(15)可简化为
ꢀꢀꢀꢀꢀ
（16）仅在的所有系数均为零时，式(9)才能成立，即：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（17）系数矩阵行列式的形式为：
ꢀꢀꢀꢀꢀ
（18）其中：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（19）；将式(12)和式(18)表示为如下形式的方程：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（20）把和看做常数，则式(20)为关于、的隐函数，通过matlab的“fimplicit”命令求得在=1附近和=0.25附近的ε随λ的变化曲线、，即稳定性界限。8.根据权利要求1所述的一种超高速飞行器纵横向耦合运动稳定性分析方法，其特征在于：步骤s6中通过、求解临界频率、、、、的计算方法如下：针对给定转动惯量参数及俯仰运动参数、来建立关于的函数，即：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（21）在稳定性图谱中绘制作为工况曲线，工况曲线依次穿过=0、=0.25和=1附近的稳定性界限形成5个交点，设交点的坐标为，由式(7)可得：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
（22）从而求得稳定性临界点对应临界频率，临界频率按照对应的值从小到大的顺序依次定义为、、、、。

技术总结
本发明属于空气动力学及飞行力学领域，公开了一种超高速飞行器纵横向耦合运动稳定性分析方法，包括如下步骤：针对需要研究的飞行器运动特性，建立俯仰通道运动方程；通过CFD计算获得滚转力矩静导数以及滚转力矩动导数；通过多项式拟合获得滚转力矩静导数关于俯仰角的函数表达式；建立考虑惯性耦合及气动/运动耦合作用的滚转通道动力学方程；通过希尔无限行列式法求解方程的稳定界限；通过稳定界限求解导致滚转运动发散的俯仰运动临界频率；通过选取不同的运动振幅和平均俯仰角，重复上述步骤。本发明解决了现有的超高速飞行器在多通道耦合失稳分析中效率低的问题，适用于考虑惯性耦合及气动/运动耦合作用的运动稳定性分析。耦合及气动/运动耦合作用的运动稳定性分析。耦合及气动/运动耦合作用的运动稳定性分析。


技术研发人员：丛戎飞 叶友达 丁智坚 温杰
受保护的技术使用者：中国空气动力研究与发展中心空天技术研究所
技术研发日：2023.06.30
技术公布日：2023/8/1