CN111859540A

CN111859540A - 一种大气扰动中飞机颠簸响应的计算方法

Info

Publication number: CN111859540A
Application number: CN202010644869.5A
Authority: CN
Inventors: 高振兴; 王得宝; 向志伟; 朱园
Original assignee: Nanjing University of Aeronautics and Astronautics
Current assignee: Nanjing University of Aeronautics and Astronautics
Priority date: 2020-07-07
Filing date: 2020-07-07
Publication date: 2020-10-30
Anticipated expiration: 2040-07-07
Also published as: CN111859540B

Abstract

本发明公开了一种大气扰动中飞机颠簸响应的计算方法，属于计算、推算或计数的技术领域。本发明的目的是准确计算大气扰动飞行中飞机的法向加速度，为进一步评估飞机颠簸严重程度提供定量方法。该方法包括以下步骤：对飞机机翼和平尾的中弧面进行网格划分；采用非平面非定常涡环法计算飞机在大气扰动中的气动特性；计算扰动风增量气动导数并加入到飞行动力学模型中；实时计算飞机机身各处的法向加速度，即飞机颠簸响应，同时通过飞行动力学模型的数值仿真获得非平面非定常涡环法计算所需的各项飞行参数。

Description

一种大气扰动中飞机颠簸响应的计算方法

技术领域

本发明公开了一种大气扰动中飞机颠簸响应的计算方法，涉及民航安全技术应用领域，具体涉及一种使用非平面非定常涡环法和动力学建模方法计算飞机法向加速度的方法，属于计算、推算或计数的技术领域。

背景技术

风切变、大气湍流等大气扰动现象严重影响民航飞机的飞行品质、乘坐品质和飞行安全。在大气扰动飞行中，受到扰动风的影响，飞机加速度、角速度等状态参数发生变化，产生颠簸。轻度的颠簸会降低乘坐品质，严重的颠簸会造成人员伤亡乃至严重飞行事故。

飞机颠簸主要由飞机法向加速度和俯仰角速度变化引起。目前，计算飞行中飞机的颠簸响应主要采用小扰动线化模型和函数拟合模型。采用小扰动线化模型计算飞机颠簸响应是将飞机纵向小扰动线化方程推导成传递函数形式的滤波器，扰动风作为传递函数的输入，经计算获得飞机法向加速度，这种方法可以实时计算飞机响应，但小扰动线化模型只在偏离平衡点的小范围内有效，无法获得高频瞬时动力学响应，因此对高频、剧烈的加速度颠簸响应计算并不精确。采用函数拟合模型计算飞机颠簸响应则是一种简化方法，没有考虑具体的机型，精度较低。由于不同重量和飞行状态的飞机对大气扰动的响应不同，上述两种方法都不能准确计算飞机颠簸响应。此外，这两种方法都不能有效区分飞机机动飞行与大气扰动引起的法向加速度变化。

发明内容

本发明的目的是针对上述背景技术的不足，提供了一种大气扰动中飞机颠簸响应的计算方法，实现了大气扰动飞行中机身各处法向加速度的精确计算，解决了现有颠簸响应计算精度不足、无法分辨大气扰动和机动飞行诱发的法向加速度变化的技术问题。

本发明为实现上述发明目的采用如下技术步骤：

步骤一：根据目标飞机机翼和平尾的翼型、翼展、前缘后掠角、翼根和翼梢弦长等参数，对机翼和平尾的中弧面进行网格划分并布置涡环；

步骤二：使用非平面非定常涡环法计算目标飞机在瞬时三维扰动风下的气动特性，分别获得扰动风影响下任意涡环控制点处的气动力和无扰动风影响的气动力；

步骤三：计算飞机机翼和平尾在大气扰动下的气动力和俯仰力矩形成增量气动导数项，将增量气动导数项加入到飞机气动模型中，并通过飞行动力学模型解算获得飞机机身任意点处的法向加速度；

步骤四：更新各项飞行状态参数，从而可连续计算飞机在大气扰动中的法向加速度颠簸响应。

本发明采用上述技术方案，具有以下有益效果：针对现有飞机颠簸响应一般采用线化传递函数或函数拟合模型进行计算所存在的缺陷，在布局涡环时充分考虑翼型中弧面的构造特点在网格分界点处将涡丝分割为两段，通过多段涡丝诱导速度矢量和的方式准确描述完整涡环对空间任意点的诱导速度，再结合非平面非定常涡环法分别计算有无扰动风影响的气动特性，将扰动风导致的气动特性变化引起的升力增量导数和俯仰力矩增量导数加入到飞行动力学模型中，进而通过气动特性变化引起的气动力变化和俯仰力矩变化准确计算飞机机身各处的法向加速度，并分辨出大气扰动引起的飞机颠簸，有效地避免了构建线化传递函数或拟合函数过程中产生的误差。本发明对飞机在大气扰动中的颠簸响应计算、飞行品质监控、大气扰动诱发的飞行事故与安全性分析具有积极意义。

附图说明

图1为飞机机翼和平尾网格划分及涡环布置的示意图。

图2为非平面非定常涡环法的示意图。

图3为相邻涡环的环量分布示意图。

图4为本申请计算大气扰动中飞机颠簸响应的流程图。

具体实施方式

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

本发明公开的一种大气扰动中飞机颠簸响应的计算方法如图4所示，包括以下四个步骤。

步骤一：根据目标飞机的机翼和平尾的翼型、翼展、前缘后掠角、翼根和翼梢弦长等参数，对机翼和平尾的中弧面进行网格划分。

以毛机翼的顶部定义原点，建立空间坐标系，如图1所示。在机翼的中弧面上布置N_RW行、N_CW列的网格；在平尾的中弧面上布置N_RT行、N_CT列的网格，总的网格数为N个。由于翼梢处的环量分布沿展向变化较大、前缘处的环量分布沿弦向变化较大，因此按半圆法来进行网格划分。考虑到高空飞行过程中，平尾上的升降舵、机翼上的扰流板会根据飞行控制系统指令发生偏转。因此将这类控制面的结构边界作为网格划分的自然边界。确定各个网格点的空间几何坐标。针对控制面的偏转，对相应网格点的空间几何坐标进行修正，从而准确反映由于舵面偏转带来的坐标变化，方便后续进行准确的气动计算。

在网格上布置涡环。如图1所示，涡环的前端放置于所划网格1/4弦线位置，后端放置于后一行网格的1/4弦线处，控制点位于网格3/4弦线中点。对于位于翼面后缘的涡环，其后端拖出后缘1/4弦长。为更加逼近翼型的中弧面，将AE 段涡丝在分界点F处分为AF、FE两段，BD段涡丝在分界点C处分为BC、CD 段。物面法向量n位于控制点处。

根据毕奥-萨瓦尔定律，涡环上某段涡丝在空间中任意一点处产生的诱导速度矢量为：

式(1)中，V是涡丝所诱导产生的速度矢量，Γ为涡环环量且满足右手定则，r₀是涡丝起点到终点的位置矢量，r₁是涡丝起点到诱导点的位置矢量，r₂是涡丝终点到诱导点的位置矢量。r₁，r₂分别是矢量r₁，r₂的模，引入涡丝诱导速度系数矢量k，令，

从而V＝Γk。因此，对于图1中任意一个完整矩形涡环，其对空间任意点的诱导速度可以表示为6段涡丝诱导速度的矢量和：

V＝Γ(k_AB+k_BC+k_CD+k_DE+k_EF+k_FA)＝ΓK (3)，

式(3)中，K为完整的矩形涡环对空间任意点的诱导速度系数矢量。

步骤二：在每个时间步长内，根据瞬时三维扰动风和空速矢量、迎角、侧滑角、姿态角等飞行状态参数，基于非平面非定常涡环法计算飞机气动特性。

基于非平面非定常涡环法计算飞机气动特性的过程包括：步骤A、步骤B、步骤C。

步骤A：计算飞机姿态运动和瞬时三维扰动风诱发的当地速度变化。飞机姿态运动和瞬时三维扰动风诱发的当地速度由自由来流速度分量 V_∞(t_k)＝[V_∞x(t_k),V_∞y(t_k),V_∞z(t_k)]^T、飞机姿态角速度分量[p(t_k),q(t_k),r(t_k)]^T和瞬时扰动风分量W(t_k)＝[W_x(t_k),W_y(t_k),W_z(t_k)]^T三部分组成。在t_k时刻，翼面任意点的当地速度可以表示为：

式(4)中，[x,y,z]^T为翼面任意点相对于飞机质心的当地坐标，α(t_k),β(t_k) 分别表示t_k时刻的迎角、侧滑角。

步骤B：计算物面边界条件，获得翼面上各个涡环的强度。

在以机翼翼根前缘顶点为原点的体轴系中，翼面任意控制点处满足流场法向合速度为零的无穿透物面边界条件：

式(5)中，Φ为翼面任意控制点处的流场位函数，

为控制点处的流场三轴速度分量，i、j、k为以机翼翼根前缘顶点为原点的坐标系的 x轴、y轴、z轴的三轴单位向量。

在t_k时刻，如图2所示，飞机机翼和平尾后缘的涡环以当地速度脱落形成一系列尾涡。按非定常涡环方法，首先按下式确定后缘网格分界点m在t_k-1时刻的当地速度：

式(6)中，

为t_k-1时刻翼面附着涡系在m点处的诱导速度，Kj

m 为翼面上第j个涡环对分界点m的诱导速度矢量，N为翼面上布局的涡环总数，Γ_j(t_k-1)为t_k-1时刻第j个涡环的涡环环量；

为t_k-1时刻m点处飞机姿态运动和瞬时三维扰动风诱发的当地速度，可由(4)式求得；

为t_k-1时刻尾涡系在分界点m处的诱导速度，可由(1)式求得。t_k时刻尾涡涡环拐点系相较于t_k-1时刻顺坐标系方向前移了 V_m(t_k-1)Δt，从而确定一对新的涡环拐点，并构建新的尾涡涡环。最新尾涡涡环的强度与t_k时刻的后缘涡环强度相同。

进一步修正尾涡系任意拐点n的位置：

式(7)中，

为t_k-1时刻翼面附着涡在拐点n的诱导速度，

为翼面上第j个涡环对拐点n的诱导速度矢量；[V_Wx,n(t_k-1),V_Wy,n(t_k-1),V_Wz,n(t_k-1)]^T为t_k-1时刻尾涡系在拐点n处的诱导速度，可由(1)式求得；[W_x(t_k-1),W_y(t_k-1),W_z(t_k-1)]^T为t_k-1时刻的扰动风速。

从而获得t_k时刻尾涡系在控制点i处的诱导速度。控制点处的无穿透物面边界条件为：

求解该代数方程，即可获得各涡环强度Γ₁(t_k),Γ₂(t_k),...,Γ_N(t_k)。

在t＝0时刻，飞机开始以速度V_∞(0)运动，此时无涡环存在。在第一个时间步长t₁＝Δt时，仍可按(8)式计算获得各面元的涡环强度，但此时只存在翼面附着涡环，无尾涡影响。

步骤C：根据库塔-茹科夫斯基定理进行空气动力计算。如图3所示，任意涡环环量与其相邻涡环环量需进行叠加以得到实际环量分布。以右半对称翼面为例，任意涡环控制点i上的气动力为：

式(9)中，F_i(t_k)＝[F_x,i(t_k),F_y,i(t_k),F_z,i(t_k)]^T表示t_k时刻扰动风影响下第n_RW行n_CW列涡环控制点i上的气动力，ρ为空气密度，

为t_k时刻第n_RW行 n_CW列的控制点i对应涡环的环量，

分别为t_k时刻第n_RW-1 行n_CW列、第n_RW行n_CW+1列对应涡环的环量，V_b,i(t_k)为t_k时刻翼面附着涡系在控制点i处的诱导速度，

为t_k时刻控制点i处飞机姿态运动和瞬时三维扰动风诱发的当地速度，

为t_k时刻尾涡系在控制点i处的诱导速度，r_AB、r_BC、r_CD分别为控制点i对应涡环的AB、BC、 CD段涡丝向量；对于翼面前缘，上方不存在涡环，(9)式中

对于后缘位置，后缘涡环的后半部分并不处于翼面内，(9)式中无r_CD项。

上述计算过程中，采用同样的方法，令瞬时扰动风分量W(t_k)＝0，获得无风扰动下的各面元涡环强度和任意涡环控制点i上的气动力

该气动力包含了飞机舵面偏转以及机动飞行对气动力的影响。

步骤三：计算飞机机翼和平尾在大气扰动下的气动力和力矩，形成增量气动导数项，加入到飞机气动模型中，并通过飞行动力学模型解算获得质心处法向加速度a_n和机身任意点处法向加速度a′_n

将步骤二获得的体轴系下的含扰动风的气动力F_i(t_k)和不含扰动风的气动力F_i ⁰(t_k)相减，获得完全由扰动风引起的气动力变化，并将气动力变化转换到风轴系中：

式(10)中，ΔF_x(t_k),ΔF_y(t_k),ΔF_z(t_k)分别为扰动风导致的阻力增量、侧力增量和升力增量。相应地，t_k时刻绕飞机重心[x_cg,y_cg,z_cg]^T的俯仰力矩增量ΔM_y(t_k)为：

式(11)中，x_i、z_i为控制点i的x轴坐标和z轴坐标。

从而，计算由大气扰动引起的升力和俯仰力矩增量导数为：

式(12)中，

分别为t_k时刻升力和俯仰力矩的增量导数，S 为翼面有效面积，c为平均弦长，ρ是大气密度，V(t_k)是t_k时刻的空速。该增量气动导数只包含扰动风影响，与飞机舵面偏转和机动飞行无关。将计算获得的增量气动导数加入到飞机气动模型中，获得实时气动力F_Az和俯仰力矩M_Ay并代入动力学方程：

式(13)中，

为t_k时刻飞机地速的x轴分量和y轴分量，θ(t_k) 为t_k时刻的俯仰角，φ(t_k)为t_k时刻的滚转角，m为飞机质量，

为飞机惯性矩阵，I_x、I_y、I_z为三轴惯性矩，I_xz为x轴z轴惯量积。t_k时刻飞机质心处的法向加速度a_n(t_k)通过(13)式的第1个等式计算获得，t_k时刻俯仰角速度q(t_k)由(13)式的第2个等式计算获得。飞机机身上位于任意位置[x,y,z]^T处的法向加速度a′_n可由下式计算：

步骤四：对完整的动力学方程进行数值积分获得空速、迎角、侧滑角、姿态角等若干飞行状态或通过外部传感量测值更新飞行状态参数，将更新后的飞行状态参数作为步骤二的输入条件，循环计算，从而可连续计算飞机在大气扰动中的法向加速度颠簸响应。

Claims

1.一种大气扰动中飞机颠簸响应的计算方法，其特征在于，

根据目标飞机机翼和平尾的参数对机翼和平尾的中弧面进行网格划分，按照在当前网格下移1/4弦线位置布局涡环的方式布置涡环；

采用非平面非定常涡环法计算目标飞机在瞬时三维扰动风下的气动特性获得扰动风影响下任意涡环控制点处的气动力和无扰动风影响下任意涡环控制点处的气动力；

根据扰动风影响下任意涡环控制点处气动力与无扰动风影响下任意涡环控制点处气动力的差值得到仅由扰动风引起的气动力变化，根据仅由扰动风引起的气动力变化计算目标飞机机翼和平尾在大气扰动下的气动力和俯仰力矩形成增量气动导数项，将增量气动导数项加入飞机气动模型后求解实时气动力及俯仰力矩，再结合飞行动力学方程获得目标飞机机身任意点处的法向加速度；

更新目标飞机的飞行状态参数，通过连续计算目标飞机法向加速度得到大气扰动中的飞机颠簸响应。

2.根据权利要求1所述一种大气扰动中飞机颠簸响应的计算方法，其特征在于，所述涡环为由沿相邻上网格分界线方向的AB段涡丝、沿弦长方向的BD段涡丝、沿相邻下网格分界线方向的DE段涡丝、沿弦长方向的EA段涡丝组成的矩形涡环，沿弦长方向的BD段涡丝在网格分界线处分割为BC段涡丝和CD段涡丝，沿弦长方向的EA段涡丝在网格分界线处分割为EF段涡丝和FA段涡丝。

3.根据权利要求2所述一种大气扰动中飞机颠簸响应的计算方法，其特征在于，所述涡环对空间任意点的诱导速度为各段涡丝诱导速度的矢量和，V＝Γ(k_AB+k_BC+k_CD+k_DE+k_EF+k_FA)＝ΓK，V为涡环对空间任意点的诱导速度，Γ为涡环环量，k_AB、k_BC、k_CD、k_DE、k_EF、k_FA分别为AB段、BC段、CD段、DE段、EF段、FA段涡丝诱导速度系数矢量，K为整个涡环对空间任意点的诱导速度系数矢量。

4.根据权利要求3所述一种大气扰动中飞机颠簸响应的计算方法，其特征在于，扰动风影响下任意涡环控制点处的气动力为

F_i(t_k)＝[F_x,i(t_k),F_y,i(t_k),F_z,i(t_k)]^T表示t_k时刻扰动风影响下第n_RW行n_CW列的控制点i上的气动力，ρ为空气密度，

为t_k时刻第n_RW行n_CW列的控制点i对应涡环的环量，

分别为t_k时刻第n_RW-1行n_CW列、第n_RW行n_CW+1列涡环的环量，V_b,i(t_k)为t_k时刻翼面附着涡系在控制点i处的诱导速度，[V_x,i(t_k),V_y,i(t_k),V_z,i(t_k)]^T为t_k时刻控制点i处飞机姿态运动和瞬时三维扰动风诱发的当地速度，

为t_k时刻尾涡系在控制点i处的诱导速度，r_AB、r_BC、r_CD分别为控制点i对应涡环的AB、BC、CD段涡丝向量；计算翼面前缘涡环控制点i在扰动风影响下的气动力时，

计算翼面后缘涡环控制点i在扰动风影响下的气动力时，省略r_CD项。

5.根据权利要求4所述一种大气扰动中飞机颠簸响应的计算方法，其特征在于，计算无扰动风影响下任意涡环控制点处的气动力的方法为：将[V_x,i(t_k),V_y,i(t_k),V_z,i(t_k)]^T中的扰动风分量置零后按照扰动风影响下任意涡环控制点处气动力的计算方法进行计算，记t_k时刻无扰动风影响下任意涡环控制点i上的气动力为

6.根据权利要求5所述一种大气扰动中飞机颠簸响应的计算方法，其特征在于，根据仅由扰动风引起的气动力变化计算目标飞机机翼和平尾在大气扰动下的气动力和俯仰力矩的方法为：将气动力变化转换到风轴系中得到扰动风导致的阻力增量、侧力增量和升力增量以及绕飞机重心的俯仰力矩增量，

ΔF_x(t_k),ΔF_y(t_k),ΔF_z(t_k)分别为t_k时刻扰动风导致的阻力增量、侧力增量和升力增量，α(t_k)、β(t_k)分别为t_k时刻的迎角、侧滑角，[x_cg,y_cg,z_cg]^T为飞机重心坐标，ΔM_y(t_k)为t_k时刻绕飞机重心的俯仰力矩增量。

7.根据权利要求6所述一种大气扰动中飞机颠簸响应的计算方法，其特征在于，形成的增量气动导数项为

分别为t_k时刻升力和俯仰力矩的增量导数，S为翼面有效面积，c为平均弦长，ρ是大气密度，V(t_k)是t_k时刻的空速。

8.根据权利要求7所述一种大气扰动中飞机颠簸响应的计算方法，其特征在于，再结合飞行动力学方程获得目标飞机机身任意点处的法向加速度的方法为：将实时气动力和俯仰力矩带入飞行动力学方程：

求解飞行动力学模型得到飞机质心处的法向加速度和俯仰角速度，再由表达式：

计算目标飞机机身任意点处的法向加速度，其中，

为t_k时刻目标飞机地速的x轴分量和y轴分量，θ(t_k)为t_k时刻的俯仰角，φ(t_k)为t_k时刻的滚转角，m为飞机质量，

为飞机惯性矩阵，I_x、I_y、I_z为三轴惯性矩，I_xz为x轴z轴惯量积，a_n(t_k)为t_k时刻飞机质心处的法向加速度，q(t_k)为t_k时刻俯仰角速度，[x,y,z]^T为目标飞机机身任意点处的位置坐标，a′_n(t_k)为t_k时刻飞机机身任意点处的法向加速度。

9.根据权利要求1至8中任意一项所述一种大气扰动中飞机颠簸响应的计算方法，其特征在于，通过飞行动力学模型的数值仿真或外部传感量测数据更新目标飞机的飞行状态参数。