CN110941920B - 一种用于无人机飞行载荷数据计算与后处理的方法 - Google Patents

Abstract

本发明提供了一种用于无人机飞行载荷数据计算与后处理的方法，本方法以流体网格作为计算单元，通过对每个流体网格气动载荷与惯性载荷数据的计算与积分，得到整机飞行载荷分布。并基于无人机强度专业载荷筛选与计算需求，计算无人机全机轴向任意位置剪力图与弯矩图，并对飞行载荷进行配平处理，最终处理成为能够输入有限元软件进行强度计算的载荷数据。本发明适用于飞行任务马赫数较高，飞行载荷较大，对结构强度有设计需求的无人机型号。

Description

一种用于无人机飞行载荷数据计算与后处理的方法

技术领域

本发明属于飞行载荷数据计算的建立方法，具体的涉及一种用于无人机飞行载荷数据计算与后处理的建立方法。

背景技术

无人机飞行载荷的计算作为无人机飞机设计中的关键流程之一，对无人机全机结构设计，结构尺寸重量优化，强度校核等具有非常重要的影响。可以断定的是，无人机飞行载荷能否准确计算，直接关系到无人机结构设计是否安全合理，重量优化是否准确可靠等。无人机的飞行载荷计算的流程主要是基于无人机气动专业机体气动数据与无人机总体专业提供的飞行姿态数据，进行飞机气动载荷与惯性载荷的解算并进行耦合叠加，得到不同飞行姿态下飞机飞行载荷数据，之后通过载荷数据筛选最终得到无人机危险状态点，最终基于危险状态点飞行载荷输入进行全机结构详细设计与优化及零部件强度校核。

基于飞行载荷计算相关已发表论文及《无人机设计手册》等专著内容，现有无人机飞行载荷计算主要通过载荷专业进行气动数据与惯性载荷工程近似计算与叠加，危险载荷通过设计状态点进行手动筛选。这一传统算法存在的问题主要有：

(1)气动载荷与惯性载荷计算单元不匹配

在传统算法中，气动载荷输出主要通过人工近似计算或者计算流体力学软件分析计算，前者计算单元主要是近似剖面，后者计算单元主要是机体表面流体网格，而惯性载荷依托于全机飞行参数，计算单元为全机整体，因此传统算法中将气动载荷与惯性载荷耦合需要将气动载荷与惯性载荷共同简化为同一计算单元，再通过人工分析计算叠加。这一流程放大了气动载荷与惯性载荷耦合叠加的误差，降低了飞行载荷计算的精确性和可靠性。

(2)飞行载荷输出与有限元软件输入格式不匹配

在传统算法中，飞行载荷输出主要通过耦合近似模型决定，具有较大不确定性。而有限元结构分析软件基于有限单元法，分析输入类型主要有节点加载，线加载，面加载等加载方式。耦合模型与有限元输入往往不能直接匹配，需要通过进一步后处理后才能作为有限元软件结构分析输入。而耦合模型的不确定性导致后处理步骤往往难以标准化处理，造成了飞行载荷计算难以流程化和规范化的问题。

(3)危险载荷状态点筛选时间成本大

在传统算法中，危险载荷筛选主要通过选择控制点，绘制不同控制点在各种飞行模态下剪力-弯矩包线图。这一做法的准确性高低直接取决于控制点选取多少。传统做法为了保证危险载荷筛选全面，通常选取较多典型状态点在不同飞行状态的剪力数值与弯矩数值，这一流程耗时较久，往往与无人机设计的短周期，低成本要求相悖。

(4)飞行载荷配平方法随机性大

在传统算法中，飞行载荷配平方法主要通过改变机体质量分布调节惯性载荷数值，最终使加载满足空间静力学三方向力与弯矩平衡。这一方法随机性较大，原因是满足空间静力学六个方程的质量分布并不是唯一的(理论上有无穷大种可能性的质量分布均可以满足配平要求)，质量分布的不唯一性直接导致惯性载荷输出不唯一，导致最终配平结果缺乏稳定性与可靠性。

综上所述，目前针对无人机飞行载荷传统计算流程具有精确度不高，随机性大，无法规范化等特点。因此开发一种具有规范化和流程化的无人机飞行载荷计算和后处理方法能够提高无人机飞行载荷计算精度，缩短无人机研制周期，具有重要实际意义。

发明内容

鉴于上述技术问题，本发明提供了一种用于无人机飞行载荷数据计算与后处理的建立方法。

本发明的技术方案如下：

一种用于无人机飞行载荷数据计算与后处理的方法，

第一步，建立飞行载荷数据计算标准坐标系，在本方法中，该标准坐标系与飞机坐标系一致；

第二步，针对气动载荷网格数据，建立对应的网格气动载荷计算公式：

F_i＝(F_x，F_y，F_z)^T＝F·(n_x，n_y，n_z)^T

其中：F为气动载荷数值，n_i为气动载荷在i方向标准方向向量；

第三步，基于流体网格做为计算单元，计算每个流体网格惯性载荷，其表达式为：

F₁＝F_a+F_ω+F_α

其中：F₁为网格惯性载荷，F_a为无人机飞行线加速度产生惯性载荷，F_ω为无人机飞行角速度产生惯性载荷，F_α为无人机飞行角加速度产生载荷；

第四步，将飞行载荷线性耦合，得到飞行载荷计算公式：

K＝F+F₁+m_ig

其中：K为网格飞行载荷，F为网格气动载荷，F₁为网格惯性载荷，m_ig为网格质量产生重力。

第五步，定义无人机飞行载荷筛选区域方法如下，给定的封闭空间Φ，若计算单元载荷作用点落在该封闭空间内，则判定该网格符合封闭空间筛选要求。

第六步，根据飞行载荷的耦合，进一步给出飞机沿轴向方向剪力和弯矩的计算公式，假定该数值位于轴向x₀位置处，并基于飞行载荷剪力与弯矩数值进行危险载荷筛选，

F_q＝∑m_i+F_1x+F_x

M_q＝∑(m_ix_q+F_1xx_q+F_xx_q)+J_yα_y

其中加和符号表示筛选累加位于0＜x＜x₀筛选条件下的所有网格，m_i表示单一网格质量，F_1x表示单一网格惯性载荷在x方向分量，F_x表示单一网格气动载荷在x方向分量，x_q表示单一网格在x方向距机体质心距离，J_y表示y方向无人机转动惯量数值。α_y表示无人机y方向角加速度。

第七步，将选取飞行载荷进行配平处理。具体的，将α_x，α_y，α_z，ω_x，ω_y，ω_z六个参数设为未知参数，连立解算下列线性方程组，

∑F_x+F_1x＝0

∑F_y+F_1y+m_i＝0

∑F_z+F_1z＝0

∑F_xk_ix+F_1xk_ix+m_ik_iy+J_yα_y＝0

∑F_zk_iz+F_1zk_iz+J_xα_x＝0

其中k_ij表示该网格在j方向投影到机头参考点的距离，F_i表示网格气动载荷在i方向数值，F_1i表示网格惯性载荷在i方向数值，J_i表示机体关于i轴对应转动惯量，α_i表示机体在i轴方向上角加速度数值。

当方程组对应矩阵为非奇异矩阵时，方程组具有唯一解，解算变量α_x，α_y，α_z，ω_x，ω_y，ω_z后迭代入第二步进行惯性载荷重新解算，解算出配平后的飞行载荷数值用于有限元加载；注意到，工程上无人机主要受到对称载荷影响，一般在只考虑对称载荷的影响下，原方程组退化为仅含α_x，ω_x两个未知量的方程组，给出该方程组如下：

∑F_x+F_1x+m_i＝0

进一步地，所述第二步中，

F＝PS

式中，P为网格对应压强数值，S为三角形网格面积，x_ij，y_ij，z_ij分别为网格参数。

进一步地，所述第三步中，

F_a＝(F_ax，F_ay，F_az)^T＝-m_i·(a_x，a_y，a_z)^T

其中，m_i为网格质量，a_i为线加速度方向，F_ωi ^T，F_αi ^T为角速度，角加速度对应惯性力数值矩阵，n，l为其对应方向矩阵。

进一步地，所述第五步中，若该计算单元载荷作用点位于封闭空间对应边界曲面Γ上，此时也判定该计算单元属于封闭空间内。

进一步地，在工程上，为了简化积分区域，定义该封闭空间Φ为长方体，且长方体互相垂直的三边与飞机坐标系相平行，即给定空间范围为： x₀＜x＜x₁，y₀＜y＜y₁，z₀＜z＜z₁。

本发明产生的有益效果在于：本发明采取流体网格作为计算单位，基于流体网格计算的无人机飞行载荷相对传统方法选取机身剖面来看，细分了计算尺度，从而提高了计算精度与准确性；本方法定义了一种载荷积分方法和危险载荷筛选方法，该方法流程性较好，能够通过程序语言实现筛选，相对传统手工筛选危险点方法能够提高筛选范围，并且降低设计时间，缩短了无人机设计周期；本方法定义了一种基于解算线性方程组的无人机载荷配平方法，这种方法能够统一无人机配平后修正惯性载荷数值，避免了传统修改配重方法配平时全机重量无法闭合的问题。总体来说，本发明基于无人机气动专业的飞行数据库基础，输出数据能够直接作为无人机设计专业和强度专业作为参考，具有较强实际应用价值。本发明中提出的一系列创新方法能够有效提高无人机载荷计算精度，缩短设计周期，从而节约无人机研发人力成本与物力成本。

附图说明

图1本发明计算标准坐标系示意图；

图2本发明计算流程图；

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。

在本发明的描述中，需要理解的是，术语“中心”、“长度”、“宽度”、“厚度”、“上”、“下”、“前”、“后”、“左”、“右”、“竖直”、“水平”、“顶”、“底”“内”、“外”、“轴向”、“周向”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此不能理解为对本发明的限制。

在本发明中，除非另有明确的规定和限定，术语“设置”、“安装”、“相连”、“连接”、“固定”等术语应做广义理解，例如，可以是固定连接，也可以是可拆卸连接；可以是机械连接；可以是直接相连，也可以通过中间媒介间接相连。对于本领域的普通技术人员而言，可以根据具体情况理解上述术语在本发明中的具体含义。

此外，术语“第一”、“第二”等仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括一个或者更多个该特征。

本发明为了解决上述技术问题，提出了如下解决方案

具体方案如下：

一种用于无人机飞行载荷数据计算与后处理的方法，

第一步，建立飞行载荷数据计算标准坐标系，在本方法中，该标准坐标系与机体坐标系一致。机体坐标系规定飞机右滚、抬头、右偏航为正，另规定坐标系原点为机头。如图1所示。

第二步，建立基于无人机表面流体网格输出格式的网格物理参数储存矩阵。

流体网格解算后输出数据格式为(p₁，p₂，p₃，P)及(p_i，x_i，y_i，z_i)。其中p₁， p₂，p₃表示网格逆时针排列的节点序号，P表示该三角形网格解算压强，x_i，y_i， z_i表示对应标号节点坐标，建立矩阵如下：

其中P表示该网格压强参数，S表示该网格面积参数，F表示该网格气动力参数。考虑到流体网格尺寸和曲率较小，近似的用平面网格面积代替曲面网格面积，此时曲面网格几何形状退化为三角形形状。

基于海伦公式，S可由以下方程得出。

其中a，b，c为网格三角形对应边长。

b，c取值循环下标。

网格气动力参数数值F由以下公式计算。

F＝PS (4)

其中S由式(2)计算得出。

推导基于单个网格下的气动载荷计算方法，首先推导单个网格下气动载荷作用点。考虑到流体网格实际尺寸非常小，因此可以近似认为以流体网格为计算单位的前提下，单个流体网格是均质的。此时假定气动载荷作用点为网格质心，网格质心计算公式如下。

推导单个网格下气动载荷法向公式，同样的，考虑到网格尺寸关系，近似地认为单个曲面网格为平面网格，此时气动载荷方向与网格平面法向量平行。气动载荷方向向量计算公式如下。

a_i1＝(x₁-x₂，y₁-y₂，z₁-z₂)＝(x_i1，y_i1，z_i1)

a_i2＝(x₂-x₃，y₂-y₃，z₂-z₃)＝(x_i2，y_i2，z_i2)

将方向向量标准化后，得到气动载荷方向向量为

从而气动载荷在三个基准方向分量数值大小计算公式如下：

F_i＝(F_x，F_y，F_z)^T＝F·(n_x，n_y，n_z)^T (10)

其中n_i由式(9)计算得出。

第三步，计算基于该流体网格计算单元下惯性载荷F₁数值，首先基于刚体动力学基本原理，以下方程成立。

F₁＝F_a+F_ω+F_α (11)

其中：F_a为无人机飞行线加速度产生惯性载荷，F_ω为无人机飞行角速度产生惯性载荷，F_α为无人机飞行角加速度产生载荷。

基于达朗贝尔原理，F_a可由以下方程推知

F_a＝(F_ax，F_ay，F_az)^T＝-m_i·(a_x，a_y，a_z)^T (12)

其中，m_i为单个流体网格质量，a_i为无人机给定飞行姿态线加速度方向无人机飞行角速度产生惯性载荷由以下方程推知

F_ω＝F_ω′·n^T

该方程中包含角速度产生惯性载荷数值矩阵与方向矩阵两项，由以下方程推知

F_ω′＝(F_ωx′，F_ωy′，F_ωz′)^T＝-m_i·(ω_x ²r₁，ω_y ²r₂，ω_z ²r₃)^T

其中Q_i表示机体在i方向重心数值，ω_x，ω_y，ω_z表示机体在x，y，z方向上各自飞行姿态角速度数值。方向矩阵n_ij中i依次表示气动载荷角速度关于 x，y，z三轴旋转时的方向特征，j依次表示气动载荷惯性力方向依次分解到 x，y，z轴后标准化的方向特征向量。

无人机飞行角加速度惯性载荷由以下方程推知

该方程中包含角加速度产生惯性载荷数值与方向两项，由以下方程推知

其中l_ij中i依次表示气动载荷角速度关于x，y，z三轴旋转时的方向特征， j依次表示气动载荷惯性力方向一次分解到x，y，z轴后标准化的方向特征向量。

第四步，将气动载荷与惯性载荷线性耦合为飞行载荷。基于广义力学原理，对于全机系统，所有惯性力及施加的外力，经过符合约束条件的虚位移，所作的虚功总和为零。也就是说，对于无人机全机而言，以下力学关系式是恒成立的。

F+G+F₁＝0 (20)

其中，F表示机体受到气动力，G表示机体自重，F1表示机体所受惯性力。

基于以上方程，可推导单个网格飞行载荷K计算方程如下

K＝F+F₁+m_i

K_x＝F_nx+F_1x+m_i

F_y＝F_ny+F_1y

K_z＝F_nz+F_1z (21-24)

第五步，定义网格飞行载荷累加积分筛选方法。，给定的封闭空间Φ，若该计算单元载荷作用点落在该封闭空间内，则判定该网格符合封闭空间筛选要求。特别的，若该计算单元载荷作用点位于封闭空间对应边界曲面Γ上，此时也判定该计算单元属于封闭空间内。

在工程上，为了简化积分区域，往往定义该封闭空间Φ为长方体，且长方体互相垂直的三边与飞机坐标系相平行。即给定空间范围为：

x₀＜x＜x₁，y₀＜y＜y₁，z₀＜z＜z₁。

第六步，推导无人机轴向在给定飞行姿态下弯矩数值与剪力数值计算公式。

假定该数值位于轴向x₀位置处，此位置剪力和弯矩计算公式如下给出

F_q＝∑m_i+F_1x+F_x (25)

M_q＝∑(m_ix_q+F_1xx_q+F_xx_q)+J_yα_y (26)

其中加和符号表示筛选累加位于0＜x＜x₀筛选条件下的所有网格，J_y表示y 方向无人机转动惯量数值。α_y表示无人机y方向角加速度。

基于上述公式，显然地，无人机弯矩数值和剪力数值直接与飞行载荷相关，也就是与飞行姿态相关。因此有必要定义危险飞行姿态的筛选方法，通过弯矩和剪力计算数值判定在各类飞行姿态中对飞机结构可能造成较大影响的姿态，定义如下：

对每种飞行姿态，将x₀遍历全机轴向长度，计算每种飞行姿态下任意位置的弯矩F_q(x)和M_q(x)，判断F在x取任意位置时，所有飞行姿态F_q(x)和M_q (x)最大值。定义取最大值时对应的飞行姿态为危险姿态，将所有危险姿态进行后处理后进行有限元解算。

第七步，将危险飞行载荷进行配平处理。由于飞行载荷在计算过程中，总体重心，转动惯量等参数与基于机体实际重量分布的网格质量间存在偏差，因此危险载荷无法满足三方向静力方程，需要对飞行载荷作相应配平处理。具体的，将原本已知的α_x，α_y，α_z，ω_x，ω_y，ω_z六个参数令为未知参数，连立解算下列线性方程组。

∑F_x+F_1x＝0

∑F_y+F_1y+m_i＝0

∑F_z+F_1z＝0

∑F_xk_ix+F_1xk_ix+m_ik_iy+J_yα_y＝0

∑F_zk_iz+F_1zk_iz+J_xα_x＝0 (27-33)

其中k_ij表示对应网格在j方向投影到机头参考点的距离，当方程组对应矩阵为非奇异矩阵时，方程组具有唯一解，解算变量α_x，α_y，α_z，ω_x，ω_y，ω_z后迭代入第二步进行惯性载荷重新解算。由于方程(27-33)实质上就是飞机三方向静力平衡方程，因此解算所得载荷必然能够满足该方程。因此可以将解算出配平后的飞行载荷数值用于有限元加载。注意到，工程上无人机主要受到对称载荷影响，一般在只考虑对称载荷的影响下，原方程组退化为仅含α_x，ω_x两个未知量的方程组，给出该方程组如下：

∑F_x+F_1x+m_i＝0

本发明产生的有益效果在于：本发明采取流体网格作为计算单位，基于流体网格计算的无人机飞行载荷相对传统方法选取机身剖面来看，细分了计算尺度，从而提高了计算精度与准确性；本方法定义了一种载荷积分方法和危险载荷筛选方法，该方法流程性较好，能够通过程序语言实现筛选，相对传统手工筛选危险点方法能够提高筛选范围，并且降低设计时间，缩短了无人机设计周期；本方法定义了一种基于解算线性方程组的无人机载荷配平方法，这种方法能够统一无人机配平后修正惯性载荷数值，避免了传统修改配重方法配平时全机重量无法闭合的问题。总体来说，本发明基于无人机气动专业的飞行数据库基础，输出数据能够直接作为无人机设计专业和强度专业作为参考，具有较强实际应用价值。本发明中提出的一系列创新方法能够有效提高无人机载荷计算精度，缩短设计周期，从而节约无人机研发人力成本与物力成本。

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

Claims (4)

1.一种用于无人机飞行载荷数据计算与后处理的方法，其特征在于：

第一步，建立飞行载荷数据计算标准坐标系，在本方法中，该标准坐标系与飞机坐标系一致；

第二步，针对气动载荷网格数据，建立对应的网格气动载荷计算公式：

F₁＝(F_x，F_y，F_z)^T＝F·(n_x，n_y，n_z)^T

其中：F为气动载荷数值；

第三步，基于流体网格做为计算单元，计算每个流体网格惯性载荷，其表达式为：

F₂＝F_w+F_a+F_α

其中：F₂为网格惯性载荷，F_a为无人机飞行线加速度产生惯性载荷，F_w为无人机飞行角速度产生惯性载荷，F_α为无人机飞行角加速度产生载荷；

第四步，将飞行载荷线性耦合，得到飞行载荷计算公式：

K＝F₁+F₂+m_ig

其中：K为网格飞行载荷，F₁为网格气动载荷，F₂为网格惯性载荷，m _ig为网格质量产生重力；

第五步，定义无人机飞行载荷筛选区域方法如下，给定的封闭空间

若计算单元载荷作用点落在该封闭空间内，则判定该网格符合封闭空间筛选要求；

第六步，根据飞行载荷的耦合，进一步给出飞机沿轴向方向剪力和弯矩的计算公式，假定数值位于轴向x₀位置处，并基于飞行载荷剪力与弯矩数值进行危险载荷筛选，

F_q＝∑m_i+F_1x+F_x

M_q＝∑(m_ix_q+F_1xx_q+F_xx_q)+J_yα_y

其中加和符号表示筛选累加位于0＜x＜x₀筛选条件下的所有网格，m_i表示单一网格质量，F_1x表示单一网格惯性载荷在x方向分量，F_x表示单一网格气动载荷在x方向分量，x_q表示单一网格在x方向距机体质心距离，J_y表示y方向无人机转动惯量数值；α_y表示无人机y方向角加速度；

第七步，将选取飞行载荷进行配平处理；具体的，将α_x，α_y，α_z，J_x，J_y，J_z六个参数设为未知参数，联立 解算 下列线性方程组，

∑F_x+F_1x＝0

∑F_y+F_1y+m_i＝0

∑F_z+F_1z＝0

∑F_xk_ix+F_1xk_ix+m_ik_iy+J_yα_y＝0

∑F_yk_iy+F_1yk_iy+J_zα_z＝0

∑F_zk_iz+F_1zk_iz+J_xα_x＝0；

当方程组对应矩阵为非奇异矩阵时，方程组具有唯一解，解算 变量α_x，α_y，α_z，J_x，J_y，J_z后迭代入第二步进行惯性载荷重新解算 ，解算 出配平后的飞行载荷数值用于有限元加载；原方程组退化为仅含α_y，J_y两个未知量的方程组，给出该方程组如下：

∑F_x+F_1x+m_i＝0

∑F_xk_ix+F_1xk_ix+m_ik_ix+J_yα_y＝0。

2.如权利要求1所述的一种用于无人机飞行载荷数据计算与后处理的方法，其特征在于：所述第二步中，

F＝PS

式中，P为网格对应压强数值，S为三角形网格面积。

3.如权利要求1所述的一种用于无人机飞行载荷数据计算与后处理的方法，其特征在于：所述第三步中，

F_a＝(F_ax，F_ay，F_az)^T＝-m_i·(a_x，a_y，a_z)^T

其中，m_i为单一网格质量。

4.如权利要求1所述的一种用于无人机飞行载荷数据计算与后处理的方法，其特征在于：所述第五步中，若该计算单元载荷作用点位于封闭空间对应边界曲面Γ上，此时也判定该计算单元属于封闭空间内。

Images