具体实施方式

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

本申请提供了一种用于移动风场环境中飞行器气动载荷的确定方法，主要包括：

步骤S1、获取当前时刻体轴系上的飞行器各速度分量；

步骤S2、对飞行器进行网格划分；

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

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

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

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

步骤S7、基于所述压力系数确定飞行器各网格处的气动力及气动力矩，返回步骤S3，进行下一时刻的风场速度计算，直至仿真结束。

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

以下结合算例进行说明。

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

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

1.2、气动载荷计算参数初始化:输入飞行器飞行高度H、飞行器飞行速度V_A(对应飞行马赫数M)、飞行器迎角α、飞行器侧滑角β、飞行器三个欧拉角(滚转角φ、俯仰角θ、偏航角ψ)。例如飞行器飞行高度H为1500m，飞行器飞行速度V_A为200.7m/s(对应飞行马赫数M为0.6)，对飞行器进行1g配平，得到飞行器配平迎角α为1.2°，飞行器侧滑角β为0°、飞行器三个欧拉角(滚转角φ为0°、俯仰角θ为1.2°、偏航角ψ为0°)。

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、获得移动风场时间轴的当前时刻点t。移动风场时间轴的当前时刻点t也就是当前的仿真时刻，根据此时刻，可以知道飞行器在移动风场中的位置。

1.5、飞行器体轴系速度分量计算:根据步骤1.2中的飞行器飞行速度V_A、飞行器迎角α、飞行器侧滑角β，按式(2)计算飞行速度V_A在体轴系X轴上的速度分量u、Y轴上的速度分量v、Z轴上的速度分量w：

根据公式2计算得到t＝0时刻飞行速度V_A在体轴系X轴上的速度分量u为200.66m/s、Y轴上的速度分量v为0m/s、Z轴上的速度分量w为4.1946m/s。

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

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

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

1.7、网格处的移动风场速度插值：任意第i个网格处所受到的水平移动风场速度U_x(t_i)、侧向移动风场速度U_y(t_i)、垂向移动风场速度U_z(t_i)根据移动风场数据插值获得。如果网格未进入风场、或者已经退出风场，则移动风场的速度为零。

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

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

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

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

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

1.12、网格处的压力系数确定：根据步骤1.1中的网格形心无因次化位置步骤1.2中的飞行马赫数M、步骤1.10中的等效气动迎角α_i、等效侧滑角β_i，任意第i个网格处压力系数Cp_i(t)通过对飞行器的原始压力分布数据Cp0(通过风洞试验或CFD仿真获得)插值得到，确定公式如(9)所示：

本实施例中，根据网格的气动状态参数、几何参数，采用多维插值方法得到网格处的压力系数。

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

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

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

通过对F_xS(t)、F_yS(t)、F_zS(t)、M_xS(t)、M_yS(t)、M_zS(t)无因次处理得到变量的无量纲形式；F_zS(t)的无因次化形式为飞行器升力系数CL，M_yS(t)的无因次化形式为飞行器气动俯仰力矩系数Cmy。

1.15、回到步骤1.4，更新计算时间，直至计算结束。

图2是本发明一个实施例的飞行器升力系数的动响应图，横轴为时间，纵轴为飞行器升力系数，带方框符号的曲线表示考虑移动风场包围飞行器过程的升力系数动响应，带三角符号的曲线表示忽略移动风场包围飞行器过程的升力系数动响应，点划线曲线表示阶跃移动风场激励信号；

图3是本发明一个实施例的飞行器气动俯仰力矩系数的动响应图，横轴为时间，纵轴为飞行器气动俯仰力矩系数，带方框符号的曲线表示考虑移动风场包围飞行器过程的气动俯仰力矩系数动响应，带三角符号的曲线表示忽略移动风场包围飞行器过程的气动俯仰力矩系数动响应，点划线曲线表示阶跃移动风场激励信号；

从图2、图3中可以看出对于忽略移动风场包围飞行器过程的计算情况，飞行器的升力、气动俯仰力矩突然变化到对应稳定状态的数值；而考虑移动风场包围飞行器过程的计算过程充分考虑了每个网格进入风场、处于风场中位置等信息细节，因此飞行器的升力和气动俯仰力矩是逐步达到稳定状态的数值；当飞行器全部浸入阶跃移动风场中后，上述两种计算结果的气动力一样。

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