CN113761646A - 一种用于移动风场环境中飞行器动响应的确定方法 - Google Patents

Abstract

本申请属于飞行器仿真技术领域，特别涉及一种用于移动风场环境中飞行器动响应的确定方法。该方法包括气动载荷计算参数初始化；更新飞行器表面各个网格当地的“等效时刻”；插值得到各个网格所受到的波阵风速度；计算各个网格处的迎角、侧滑角、动压；插值得到各个网格形心处的压力系数；对所有网格进行积分，得到飞行器的外部气动载荷。本发明提供了一种改进的移动风场环境中飞行器外部气动载荷的确定方法，考虑了移动风场包围飞行器过程对气动载荷的影响，简化了气动载荷的确定流程，提高了气动载荷的确定精度。

Description

一种用于移动风场环境中飞行器动响应的确定方法

技术领域

本申请属于飞行器仿真技术领域，特别涉及一种用于移动风场环境 中飞行器动响应的确定方法。

背景技术

动态响应一般是指控制***在典型输入信号的作用下，其输出量从 初始状态到最终状态的响应。

移动风场是一种特殊的突风，与普通风场模型相比有以下3点差异： (1)移动风场具有高速传播的特点，一般以超音速或音速扩张，而普通 风场一般不考虑自身的传播速度；(2)移动风场携带的扰动能量更大， 这表现在它的风场速度较大；(3)移动风场的速度是空间的，可以同时 沿着飞机的航向和法向扰动。当移动风场以音速(或略高于音速)从后方包围飞机时，相对速度较小，尾翼会先受到移动风场的作用，然后才 是机翼受到扰动。

目前，移动风场环境中飞行器动响应没有较好的方法。移动风场环 境中飞行器外部气动载荷的确定方法普遍采用面元法，气动载荷计算技 术及流程复杂，计算精度较低，尤其是在气动非线性区偏差将更大。

发明内容

为了解决上述技术问题，本申请提供了一种用于移动风场环境中飞 行器动响应的确定方法，包括：

步骤S1、对划分为网格的飞行器进行所有网格的初始气动力计算， 确定飞行器外部的气动力及气动力矩；

步骤S2、基于所述气动力及气动力矩进行动响应求解，获得在体轴 系上投影的飞行器各速度分量；

步骤S3、对给定的移动风场数据进行插值，获得飞行器各网格处受 到的与所述体轴系各轴对应的各方向的风场速度；

步骤S4、根据飞行器各速度分量及各方向的风场速度，确定飞行器 各网格处的等效速度分量；

步骤S5、根据所述等效速度分量确定飞行器各网格处的等效气动迎 角及等效侧滑角；

步骤S6、基于所述等效气动迎角及等效侧滑角，对压力分布数据进 行插值，获得飞行器各网格处的压力系数；

步骤S7、基于所述压力系数确定飞行器各网格处的气动力及气动力 矩，返回步骤S1，采用新的气动力及气动力矩替换初始的气动力及气动 力矩，进行下一时刻的飞行器动响应计算，直至仿真结束。

优选的是，步骤S2中，采用四阶龙格-库塔法对飞行器进行动响应求 解。

优选的是，步骤S3进一步包括：

步骤S11、确定仿真时间t；

步骤S12、基于移动风场包围飞行器的速度V_Δ确定飞行器各网格在移 动风场中的等效时刻t_i：t_i＝t-Δx_i/V_Δ，其中，Δx_i为第i个网格形心与飞行 器尾翼后边界的水平距离；

步骤S13、根据移动风场数据插值获得各网格处所受到的水平移动风 场速度U_x(t_i)、侧向移动风场速度U_y(t_i)、垂向移动风场速度U_z(t_i)。

优选的是，步骤S4进一步包括：

步骤S41、获取大地坐标系到飞机体轴系的变换矩阵L_bg；

步骤S42、基于所述变换矩阵确定飞行器各网格处的等效速度分量 u_i(t)、v_i(t)、w_i(t)：

其中，u(t)、v(t)、w(t)为飞行器各速度分量，U_x(t_i)、U_y(t_i)、U_z(t_i)为各 方向的风场速度。

优选的是，所述变换矩阵根据步骤S2动响应求解过程中的气动特性 数据计算。

优选的是，步骤S6中，所述压力分布数据包括飞行器的原始基础压 力分布数据及升降舵偏转增量压力分布数据，所述飞行器的原始基础压 力分布数据及升降舵偏转增量压力分布数据均通过风洞试验或CFD仿真 获得。

优选的是，在步骤S6中进行插值过程前，进一步包括：

获取飞行器的飞行马赫数；

进行飞行器1g配平，获得飞行器配平迎角及配平升降舵偏度；

获取进行飞行器表面网格划分时的各网格的形心无因次化位置。

优选的是，步骤S7进一步包括：

步骤S71、根据飞行器各网格处的压力系数、等效速压、物面法向矢 量及网格面积确定飞行器各网格处的气动力；

步骤S72、根据飞行器各网格处的气动力及形心位置矢量确定飞行器 各网格处的气动力矩；

步骤S73、根据飞行器各网格处的气动力及气动力矩确定飞行器外部 的气动力及气动力矩；

其中，物面法向矢量、网格面积及形心位置矢量均在网格划分时确 定。

优选的是，所述飞行器表面划分的网格为三角形网格。

本申请考虑移动风场包围飞行器过程对气动载荷的影响，采用风洞 试验及CFD仿真得到的压力分布数据，简化了气动载荷的确定流程，提高 了飞行器气动载荷确定精度及动响应仿真精度。

附图说明

图1是本申请用于移动风场环境中飞行器动响应的确定方法的一优 选实施例的流程图。

图2是本申请图1所示实施例的俯仰角加速度动响应图。

图3是本申请图1所示实施例的飞行器法向载荷系数动响应图。

具体实施方式

为使本申请实施的目的、技术方案和优点更加清楚，下面将结合本 申请实施方式中的附图，对本申请实施方式中的技术方案进行更加详细 的描述。在附图中，自始至终相同或类似的标号表示相同或类似的元件 或具有相同或类似功能的元件。所描述的实施方式是本申请一部分实施 方式，而不是全部的实施方式。下面通过参考附图描述的实施方式是示 例性的，旨在用于解释本申请，而不能理解为对本申请的限制。基于本 申请中的实施方式，本领域普通技术人员在没有作出创造性劳动前提下 所获得的所有其他实施方式，都属于本申请保护的范围。下面结合附图 对本申请的实施方式进行详细说明。

本申请提供了一种用于移动风场环境中飞行器动响应的确定方法， 主要包括：步骤S1、对划分为网格的飞行器进行所有网格的初始气动力 计算，确定飞行器外部的气动力及气动力矩；步骤S2、基于所述气动力 及气动力矩进行动响应求解，获得在体轴系上投影的飞行器各速度分量； 步骤S3、对给定的移动风场数据进行插值，获得飞行器各网格处受到的 与所述体轴系各轴对应的各方向的风场速度；步骤S4、根据飞行器各速 度分量及各方向的风场速度，确定飞行器各网格处的等效速度分量；步 骤S5、根据所述等效速度分量确定飞行器各网格处的等效气动迎角及等 效侧滑角；步骤S6、基于所述等效气动迎角及等效侧滑角，对压力分布 数据进行插值，获得飞行器各网格处的压力系数；步骤S7、基于所述压 力系数确定飞行器各网格处的气动力及气动力矩，返回步骤S1，采用新 的气动力及气动力矩替换初始的气动力及气动力矩，进行下一时刻的飞 行器动响应计算，直至仿真结束。

图1给出了一个具体的实施流程，参考图1，本申请的工作原理是： 将整个飞行器按三角网格进行划分，计算得到所有网格的形心位置矢量、 物面法向矢量、面积、形心无因次化位置；按照移动风场时间轴线推进， 更新飞行器表面各个网格当地的“等效时刻”，按照“等效时刻”插值 获得各网格处的移动风场速度，可以详细地考虑移动风场对不同位置网 格的影响，进而考虑网格进出移动风场的动态效应；针对不同网格，分 别计算网格当地的迎角、侧滑角、速压等参数，进而插值获得网格处的 压力系数，这对提高飞行器外部气动载荷的计算精度有积极作用；飞行 器的原始压力分布数据通过风洞试验或CFD仿真获得，数据精度高，计算 结果精度高，这对提高移动风场环境中的飞行器外部气动载荷的计算精度有积极作用；将所有网格的气动载荷累加起来得到飞行器的气动载荷， 代入飞行器动力学方程进行动响应求解。

以下结合算例进行说明。

1.1、飞行器表面三角网格划分：将整个飞行器表面划分为三角形网 格，任意第i个网格的主要信息有：形心位置矢量r_i、物面法向矢量n_i、面 积s_i、形心无因次化位置

飞行器表面三角网格的划分可以采用CATIA等商业软件完成，网格 的数量要尽量多、以提高飞行器外部气动载荷的计算精度；通过CATIA 等商业软件将飞行器表面划分为三角形网格，需要导出网格节点坐标、 网格组成节点编号等信息，根据右手法则计算网格的上述5个参量。在一 个具体实施例中，将某飞行器表面划分三角形网格，全机共得到26892个 网格。

1.2、计算参数初始化：输入飞行器飞行高度H、飞行器飞行速度V_A(对 应飞行马赫数M)、飞行器质量m、飞行器气动特性数据。例如飞行器 飞行高度H为1500m，飞行器飞行速度V_A为200.7m/s(对应飞行马赫数M 为0.6)，飞行器质量为80t。

1.3、移动风场包围飞行器速度计算：根据步骤1.2中的飞行器飞行高 度H计算出当前高度上的音速V_S、大气密度ρ，移动风场包围飞行器的速 度V_Δ计算公式如(1)所示：

V_Δ＝V_S-V_A (1)

依据飞行器飞行高度，采用标准大气计算方法得到当前高度上的音 速V_S、大气密度ρ；当移动风场从后方包围飞行器时，包围速度需按照上 式进行计算。例如，飞行高度H为1500m处的大气密度ρ为1.0581kg/m³； 音速V_S为334.5m/s，根据公式计算移动风场包围飞行器的速度V_Δ为 133.8m/s。

1.4、飞行器1g配平计算：采用步骤1.2中的飞行器飞行高度H、飞行 器质量m、飞行器气动特性数据，对飞行器进行1g配平，得到飞行器配平 迎角α^trim、配平升降舵偏度δe^trim。本实施例得到飞行器配平迎角α^trim为 1.2°、配平升降舵偏度δe^trim为-2.3°。

1.5、计算飞行器全机气动力：对飞行器所有网格的气动力进行计算， 得到飞行器外部的气动力F_S(t)、气动力矩M_S(t)；F_S(t)在X轴、Y轴、Z轴 上的投影分量分别为F_xS(t)、F_yS(t)、F_zS(t)，M_S(t)在X轴、Y轴、Z轴上的 投影分量分别为M_xS(t)、M_yS(t)、M_zS(t)。

1.6、飞行器动响应求解：根据步骤1.5中的F_xS(t)、F_yS(t)、F_zS(t)、M_xS(t)、 M_yS(t)、M_zS(t)，采用四阶龙格-库塔法对飞行器进行动响应求解，计算公 式如(2)和(3)所示：

其中，g为重力加速度；u、v和w分别为飞行器速度在体轴系X轴、 Y轴、Z轴上的投影；p、q和r分别为飞行器角速度在体轴系X轴、Y轴、 Z轴上的投影；I_xx、I_yy和I_zz分别为飞行器相对体轴系X轴、Y轴、Z轴的 惯性矩；I_zx为飞行器的惯性积。

1.7、仿真时间更新：按照移动风场的时间轴线，更新动响应仿真时 间t。

1.8、网格“等效时刻”的更新：任意第i个网格在移动风场中的“等 效时刻”t_i计算如公式(4)所示；

t_i＝t-Δx_i/V_Δ (4)

其中，Δx_i为第i个网格形心与飞行器尾翼后边界的水平距离。

需要说明的是，网格的位置差异引起的时间滞后影响，越是靠近飞 行器头部方向的网格、越晚进入移动风场，相对尾翼上的网格具有时间 滞后。

1.9、网格处的移动风场速度插值：任意第i个网格处所受到的水平移 动风场速度U_x(t_i)、侧向移动风场速度U_y(t_i)、垂向移动风场速度U_z(t_i)根据 移动风场数据插值获得。

需要说明的是，如果网格未进入风场、或者已经退出风场，则移动 风场的速度为零。

1.10、网格当地等效速度更新：任意第i个网格处的等效速度分量u_i(t)、 v_i(t)、w_i(t)计算如公式(5)所示：

其中，U_x(t_i)、U_y(t_i)、U_z(t_i)分别表示移动风场的水平分量、侧向分量、 垂向分量在时刻处t_i的风速；L_bg为大地坐标系到飞机体轴系的变换矩阵， L_bg计算如公式(6)所示：

1.11、网格当地等效空速确定：根据步骤1.10中的u_i(t)、v_i(t)、w_i(t)， 任意第i个网格处的等效空速V_i计算公式如(7)所示：

1.12、网格当地等效气动迎角α_i确定：根据步骤1.10中的u_i(t)、v_i(t)、 w_i(t)，根据步骤1.11中的V_i，任意第i个网格处的等效气动迎角α_i、等效侧 滑角β_i计算公式如(8)所示：

1.13、网格当地等效速压确定：根据步骤1.3中的大气密度ρ、步骤1.11 中的V_i，任意第i个网格处的等效速压Q_i计算公式如(9)所示：

1.14、网格处的压力系数确定：根据步骤1.1中的网格形心无因次化 位置

步骤1.2中的飞行马赫数M、步骤1.4中的δe^trim、步骤1.12 中的等效气动迎角α_i、等效侧滑角β_i，任意第i个网格处压力系数Cp_i(t)通 过对飞行器的原始基础压力分布数据Cp0(通过风洞试验或CFD仿真获得) 和升降舵偏转增量压力分布数据ΔCp0(通过风洞试验或CFD仿真获得) 插值得到，确定公式如(10)所示：

1.15、网格气动载荷的确定：根据步骤1.1中的形心位置矢量r_i、物面 法向矢量n_i、面积s_i，根据步骤1.13中的等效速压Q_i、步骤1.14中的压力 系数Cp_i(t)，任意第i个网格的气动力f_i和气动力矩m_i的计算公式如(11)所 示：

1.16、飞行器外部气动力的确定：根据步骤1.13中的气动力f_i和气动 力矩m_i，通过对全部网格积分得到飞行器外部的气动力F_S(t)、气动力矩 M_S(t)，具体的计算公式分别如(12)、(13)所示：

其中，N为飞行器外表面的所有网格数目；F_S(t)在X轴、Y轴、Z轴上 的投影分量分别为F_xS(t)、F_yS(t)、F_zS(t)，M_S(t)在X轴、Y轴、Z轴上的投 影分量分别为M_xS(t)、M_yS(t)、M_zS(t)。

1.17，回到步骤1.5，进行迭代计算，直至仿真结束，输出仿真结果。

图2是本发明一个实施例的俯仰角加速度动响应图，横轴为时间，纵 轴为俯仰角加速度，带方框符号的曲线表示考虑移动风场包围飞行器过 程的俯仰角加速度动响应，带三角符号的曲线表示忽略移动风场包围飞 行器过程的俯仰角加速度动响应，虚线曲线表示“1-cos”移动风场激励 信号。

图3是本发明一个实施例的飞行器法向载荷系数动响应图，横轴为时 间，纵轴为飞行器法向载荷系数，带方框符号的曲线表示考虑移动风场 包围飞行器过程的法向载荷系数动响应，带三角符号的曲线表示忽略移 动风场包围飞行器过程的法向载荷系数动响应，虚线曲线表示“1-cos” 移动风场激励信号。

从图2、图3中可以看出：对于考虑移动风场包围飞行器过程的计算 情况，飞行器的俯仰角加速度动响应要更为剧烈一些，且法向载荷系数 出现了明显的双峰值现象。

以上所述，仅为本申请的具体实施方式，但本申请的保护范围并不 局限于此，任何熟悉本技术领域的技术人员在本申请揭露的技术范围内， 可轻易想到的变化或替换，都应涵盖在本申请的保护范围之内。因此， 本申请的保护范围应以所述权利要求的保护范围为准。

Claims (9)

1.一种用于移动风场环境中飞行器动响应的确定方法，其特征在于，包括：

步骤S1、对划分为网格的飞行器进行所有网格的初始气动力计算，确定飞行器外部的气动力及气动力矩；

步骤S2、基于所述气动力及气动力矩进行动响应求解，获得在体轴系上投影的飞行器各速度分量；

步骤S3、对给定的移动风场数据进行插值，获得飞行器各网格处受到的与所述体轴系各轴对应的各方向的风场速度；

步骤S4、根据飞行器各速度分量及各方向的风场速度，确定飞行器各网格处的等效速度分量；

步骤S5、根据所述等效速度分量确定飞行器各网格处的等效气动迎角及等效侧滑角；

步骤S6、基于所述等效气动迎角及等效侧滑角，对压力分布数据进行插值，获得飞行器各网格处的压力系数；

步骤S7、基于所述压力系数确定飞行器各网格处的气动力及气动力矩，返回步骤S1，采用新的气动力及气动力矩替换初始的气动力及气动力矩，进行下一时刻的飞行器动响应计算，直至仿真结束。

2.如权利要求1所述的用于移动风场环境中飞行器动响应的确定方法，其特征在于，步骤S2中，采用四阶龙格-库塔法对飞行器进行动响应求解。

3.如权利要求1所述的用于移动风场环境中飞行器动响应的确定方法，其特征在于，步骤S3进一步包括：

步骤S11、确定仿真时间t；

步骤S12、基于移动风场包围飞行器的速度V_Δ确定飞行器各网格在移动风场中的等效时刻t_i：t_i＝t-Δx_i/V_Δ，其中，Δx_i为第i个网格形心与飞行器尾翼后边界的水平距离；

步骤S13、根据移动风场数据插值获得各网格处所受到的水平移动风场速度U_x(t_i)、侧向移动风场速度U_y(t_i)、垂向移动风场速度U_z(t_i)。

4.如权利要求1所述的用于移动风场环境中飞行器动响应的确定方法，其特征在于，步骤S4进一步包括：

步骤S41、获取大地坐标系到飞机体轴系的变换矩阵L_bg；

步骤S42、基于所述变换矩阵确定飞行器各网格处的等效速度分量u_i(t)、v_i(t)、w_i(t)：

其中，u(t)、v(t)、w(t)为飞行器各速度分量，U_x(t_i)、U_y(t_i)、U_z(t_i)为各方向的风场速度。

5.如权利要求4所述的用于移动风场环境中飞行器动响应的确定方法，其特征在于，所述变换矩阵根据步骤S2动响应求解过程中的气动特性数据计算。

6.如权利要求1所述的用于移动风场环境中飞行器动响应的确定方法，其特征在于，步骤S6中，所述压力分布数据包括飞行器的原始基础压力分布数据及升降舵偏转增量压力分布数据，所述飞行器的原始基础压力分布数据及升降舵偏转增量压力分布数据均通过风洞试验或CFD仿真获得。

7.如权利要求6所述的用于移动风场环境中飞行器动响应的确定方法，其特征在于，在步骤S6中进行插值过程前，进一步包括：

获取飞行器的飞行马赫数；

进行飞行器1g配平，获得飞行器配平迎角及配平升降舵偏度；

获取进行飞行器表面网格划分时的各网格的形心无因次化位置。

8.如权利要求1所述的用于移动风场环境中飞行器动响应的确定方法，其特征在于，步骤S7进一步包括：

步骤S71、根据飞行器各网格处的压力系数、等效速压、物面法向矢量及网格面积确定飞行器各网格处的气动力；

步骤S72、根据飞行器各网格处的气动力及形心位置矢量确定飞行器各网格处的气动力矩；

步骤S73、根据飞行器各网格处的气动力及气动力矩确定飞行器外部的气动力及气动力矩；

其中，物面法向矢量、网格面积及形心位置矢量均在网格划分时确定。

9.如权利要求1所述的用于移动风场环境中飞行器动响应的确定方法，其特征在于，所述飞行器表面划分的网格为三角形网格。

Images