吴翰 王正平 周 洲 王睿

摘 要：

由于垂起固定翼无人机结构复杂， 无论是采用数值模拟得到其非定常气动力还是采用单刚体建模方法建立其动力学模型都较为困难。 为了解决该问题， 通过CFD数值模拟機翼翼型的迟滞环线， 引入无人机俯仰运动所产生的迟滞效应， 通过ONERA方程建立无人机垂起改平飞过程的非定常气动力模型， 最终基于多体动力学， 将无人机划分成机翼、 机身、 旋翼、 舵面的多刚体系统， 通过凯恩方程推导并建立无人机动力学模型， 通过数值仿真可以发现： 多体动力学的方法适用于这类结构复杂无人机的建模， 引入机翼的迟滞环对无人机的俯仰运动起阻尼作用; 引入非定常气动力与未引入非定常气动力相比垂起速度存在8%的差异， 垂起改平飞过程结束时两者上升速度存在40%的差异。

关键词：       垂起固定翼无人机; 迟滞环; ONERA方程; 多体动力学; 凯恩方程

中图分类号：      TJ011;  V212   文献标识码：     A   文章编号：      1673-5048（2019）02-0021-08

0 引  言

垂起固定翼无人机一般由旋翼和固定翼所组成， 兼具垂直起降和高巡航效率的性能， 具有很大的研究价值和市场前景， 然而对于这类无人机的动力学建模仍然是个难点。 文献[1-2]分别采用牛顿-欧拉方程建立旋翼无人机动力学模型， 文献[3-4]分别采用拉格朗日方程[5]建立固定翼无人机动力学模型， 文献[6]采用牛顿-欧拉方程建立倾转旋翼无人机的动力学模型， 以上文献是将无人机的气动力看成准定常模型， 未系统地考虑无人机的非定常气动力模型。 本文主要对垂起固定翼无人机的动力学建模与非定常气动力进行了研究， 针对垂起固定翼这类特殊布局的无人机， 提出三个值得研究的问题： （1）现有单刚体建模方法， 在建立复杂结构无人机的动力学模型时较为繁琐， 能否采用新的建模方法， 以便更合理、 更快速地建立这类结构复杂的无人机的动力学模型; （2）以前的动力学模型都是将无人机看成准定常模型进行处理， 在实际情况中， 由于垂起固定翼无人机的机翼面积较大， 机翼的迟滞效应将对无人机的飞行动力学产生较大影响， 因此， 如何将迟滞效应与无人机的动力学模型耦合起来， 值得研究; （3）垂起固定翼无人机最具有研究价值的环节为垂起改平飞的过渡过程。 在该过程中无人机的迎角变化幅度较大， 原有气动模型不再适用， 而如果直接采用CFD进行该阶段的数值模拟， 则运算量太大不适用于工程领域， 因此， 能否采用更为便利的方法建立这一过程无人机的非定常气动力模型。

为了解决上述的三个问题， 对比现有的动力学建模方法， 选取凯恩方程[7]基于多体动力学[8] 建立该垂起固定翼无人机的动力学模型;  采用CFD数值模拟[9]该无人机机翼力矩系数的迟滞环线， 采用积分法得到无人机纵向动导数， 通过纵向动导数引入无人机机翼的迟滞效应;  采用ONERA

非线性气动力方程建立无人机垂起改平飞过渡过程中的非定常气动力模型， 完善垂起改平飞过渡阶段的动力学模型。

1 垂起固定翼无人机多刚体划分

该垂起固定翼无人机采用常规布局， “T”形尾， 其主要由四个垂起旋翼和两个前拉螺旋桨组成。 按照多体动力学的思路[10]将其划分为左右机翼、 机身、 旋翼、 螺旋桨、 垂尾、 左右平尾、 电机等18个刚体， 其中机翼、 垂尾、 电机与机身固结， 平尾与垂尾固结， 各舵面与翼面铰接。 忽略无人机在空中飞行时的柔性变形， 其广义坐标为 q =[ x y z φ θ ψ ]T， 表示无人机质心处的位置和姿态， 以各刚体质心为原点建立右手坐标系， 具体坐标系建立如图1所示。

2 机翼迟滞效应分析

2.1 无人机翼型迟滞环模拟

当无人机在做俯仰运动时， 机翼所产生的气动力为非定常气动力， 其控制方程为描述气体非定常运动的NS方程[11]：

由于该无人机展弦比大且是平直机翼， 因此， 可以用机翼二维翼型的数值模拟结果来代替整个三维机翼， 该机翼翼型为FX-60-126-1， 划分网格如图2所示。 强迫其绕刚心做0°到15°的俯仰运动， 得到不同飞行速度下力矩系数的迟滞环线， 如图3所示。

2.2 无人机机翼动导数求解

以迟滞环为基准， 推导无人机机翼俯仰运动动导数， 动导数是气动载荷对状态变化率（速度大小）的导数， 是描述无人机动态变化过程的重要参数， 应用积分法求取无人机动导数。 假设机翼在俯仰自由度下进行小振幅强迫振荡运动， 其迎角为

式中： k=（ωb/v）， 为振荡运动的减缩频率;  t为无量纲时间; θ=Δα=α-α0为机翼俯仰角。

式中： Cmθ0为俯仰刚度;  Cmθ  · 0为俯仰阻尼导数， 利用上文得到的周期性简谐强迫运动的气动力矩， 通过数值辨识即可直接获得俯仰阻尼导数。 当强迫运动为周期性正弦运动时， 对式（6）两边同时乘以sin（kt）和cos（kt）， 再进行积分可得

式中： ts为简谐强迫运动达到周期性解后的任意时刻; T为简谐强迫运动周期; Cm（t）的数据由CFD算得;  Cmθ， Cmθ  · 分别为无人机俯仰运动所产生的动导数。

3 ONERA非线性气动力方程分析

ONERA方程[13]目前主要被用于解决机翼的气动弹性问题。 该方程为非线性气动力方程， 适用于大迎角的情形。 ONERA方程的线性部分是对Theodorsen气动力的模拟， 其非线性部分则考虑了静态失速的影响。 通过该方程的推导以及推导的前提[14]可以发现该方程适用于无人机垂起改平飞的过渡过程， 选取ONERA方程建立该过程的非定常气动力模型。

首先给出升力系数为

CL=CLa+CLbCLa=Sz1×α  · +Sz2×α  ¨ +Sz3×α  · +CLγC  · Lb+rz1×C  · Lb+rz2×CLb=-rz2×ΔCL-

rz3× ΔCL α ×α  ·    （10）

式中：α为瞬时迎角; CLa， CLb为升力系数的组成部分; ΔCL为ONERA方程线性部分与非线性部分的差值; CLγ为非线性部分的升力系数。 式中rz1， rz2， rz3均与雷诺数有关， 需要通过雷诺数进行确定。 由于ONERA方程的线性部分是对经典Theodorsen气动力模型的拟合， 因此可以得到如下参数值：

sz1=π， sz2=π/2， sz3=0 （11）

同样可以得到ONERA方程线性部分与非线性部分的差值ΔCL为

ΔCL=

0 0≤α≤0.139 6 rad6.322 84×（α-0.139 6）  0.139 6 rad<

α≤0.314 2 rad6.322 84×（α-0.139 6）-0.422 84×

（α-0.314 2）  0.314 2 rad<α   ΔCL= 0 -0.139 6 rad≤α≤06.322 84×（α-0.139 6）  -0.314 2 rad≤

α<-0.139 6 rad6.322 84×（α-0.139 6）-0.422 84×

（α-0.314 2）   α<-0.314 2 rad

（12）

通过上式可以得到 ΔCL α ， 将等式（12）和 ΔCL α 带入式（10）中便可以得到相应情况下的升力系数。

ONERA方程的阻力系数为

CD=CDa+CDbC  ¨ Db+rD1×C  · Db=-r2D×ΔCD-rD3×α  · ΔCD=-0.042×α-0.147 3×α2-4.923×α3   （13）

式中： CDa， CDb为阻力系数的组成部分， ΔCD为线性部分阻力系数与非线性部分阻力系数的差值。 rD1， rD2， rD3同样由雷诺数确定。 由于无人机在过渡阶段对其侧向力系数Cyβ影响较小， 因此侧向力系数采用准定常假设进行求解即可。

4 机翼非定常动力学模型建立

4.1 机翼非定常气动力模型

以左机翼为例基于上文得到的动导数和ONERA方程建立其非定常气动力模型为

f z= C zl - 1 2 ρV2SCD（Re， α， θ）- 1 2 ρV2SCyβ（Re， α， θ）- 1 2 ρV2SCL（Re， α， θ）（14）

式中：  C zl为气流向左机翼的坐标系转换矩阵; V为左机翼质心处的合速度， CD， Cyβ， CL分别为左机翼的气动力系数。

左机翼质心处的力模型为

F z= f z+  00mzg （15）

式中：  F z为左机翼的力模型; mz为左机翼的质量;  f z为左机翼的非定常气动力模型。

4.2 机翼非定常力矩模型

垂起固定翼无人机的力矩主要由两部分产生， 一部分是各刚体质心处的广义主动力相对于无人机质心所产生的力矩， 该部分力矩將通过凯恩方程的偏速度矩阵进行引入。 另一部分是由于无人机在飞行过程中， 滚转、 偏航以及俯仰运动所引起的非定常气动力矩。 直接给出左机翼非定常力矩模型如下：

M z=   Clpp +Clrr   1 2 ρV 2Sc

（Cmθθ+Cmθ  · q  - ） 1 2 ρV 2Sl Cnpp +Cnrr   1 2 ρV 2Sc （16）

式中： p  - = pc 2V ， q  - = qc 2V ， r = rc 2V ， Clp， Clr， Cnp， Cnr为横航向动导数， 表示该垂起固定翼无人机滚转和偏航时左机翼所产生的非定常气动力矩;  Cmθ， Cmθ  · 为无人机做俯仰运动时所产生非定常气动力矩系数;  S为机翼面积;  l， c分别为左机翼展长和弦长;  p， q， r为机翼的角速度分量。 需要特别注意的是主要通过无人机纵向动导数Cmθ， Cmθ  · 引入无人机俯仰运动中机翼所产生的迟滞效应， 而且迟滞效应仅在无人机迎角小于15°时才具有。

5 垂起固定翼无人机动力学模型建立

多体动力学建模方法主要应用于机器人[15]等结构较为复杂的对象， 这种建模方法是将建模对象划分为多刚体系统， 针对每个刚体建立其动力学模型， 然后通过方程自身的约束条件， 将各刚体的动力学模型引入整个对象的质心， 进而完成动力学建模。 对于该垂起固定翼无人机， 这种建模方法较为适用。

基于多体动力学中的凯恩方程建立该垂起固定翼无人机动力学模型， 凯恩方程所建模型为一阶偏微分方程组， 可减少仿真计算量， 其建模的基准为系统广义主动力与系统广义惯性力之和为零， 系统广义主动力矩与系统广义惯性力矩之和为零。

首先建立无人机各刚体的广义惯性力和广义惯性力矩模型， 依据无人机各刚体质心的线加速度和角加速度， 建立各刚体广义惯性力和广义惯性力矩模型如下：

R i=- m i a i，   R *i=- J i ω   · i （17）

式中：  m i，  J i分别为刚体i的质量矩阵和转动惯量矩阵;   a i，  ω   · i分别为刚体i质心处的线加速度和角加速度矩阵;   R i，  R *i分别为刚体i质量所产生的广义惯性力和转动惯量所产生的广义惯性力矩。

无人机整机的广义惯性力和广义惯性力矩模型即为各刚体广义惯性力和广义惯性力矩模型的叠加， 具体形式如下[6]：

K =∑ i  V iq   ·   R i+ω iq   ·   R *i（18）

式中： i代表刚体的编号， 由于无人机被划分为18个刚体， 因此该方程为18个子方程的叠加;  V iq   ·  和ω iq   ·  分别为刚体i的偏线速度矩阵和偏角速度矩阵。

无人机的广义主动力为各刚体质心力模型的叠加， 无人机系统的广义主动力矩由两部分组成， 一部分是各刚体质心力模型相对于全机质心所产生的力矩; 另一部分是各刚体自身运动时绕其自身质心所产生的力矩。 采用左机翼动力学模型， 建立无人机右机翼、 机身、 左右平尾、 垂尾的力模型和非定常力矩模型， 分别为（ F y，  F b，  F p，  F q，  F r）和 （M y，  M b，  M p，  M q，  M r）; 对于无人机的前拉螺旋桨和垂起旋翼， 采用文献[1]中的方法通过流入比和前进比建立其力模型和力矩模型， 分别为（ F c，  F f，  F a，  F n，  F d，  F g）和（ M c，  M f，  M a  M n，  M d，  M g）; 对于无人机的电机采用文献[16]中的方法建立其力和力矩模型， 分别为（ F 1，  F 2，  F 3，  F 4，  F 5，  F 6）和（ M 1，  M 2，  M 3，  M 4，  M 5，  M 6）。 將上述的力和力矩模型以及左机翼模型进行叠加， 可得到该垂起固定翼无人机的广义主动力和广义主动力矩模型， 其具体形式如下[17]：

K *=∑  i   V iq   ·   （ F i）+ ω iq   ·   （ M i）（19）

式中： i表示各个刚体的编号;  F i中包含刚体i自身的重力;   M i为刚体i在机体坐标系下的合力矩， 其包含刚体i质心气动力所产生的力矩以及刚体i的自身运动所产生的非定常力矩。 各刚体质心气动力和各刚体重力相对于无人机质心所产生的力矩， 由偏线速度矩阵与力 F i的乘积引入。

将式（18）和式（19）叠加便可得到该无人机的动力学模型[18]为

K + K *=∑  i   V iq   ·   （ F i+ R i）+ ω iq   ·   （ M i+ R *i） =0  （20）

式（20）即为该垂起无人机的六自由度动力学模型， 其中 V iq   ·   ，    ω iq   ·   均为6×3的矩阵，F i，  M i，  R i，  R *i均为3×1矩阵， 因此通过式（20）可以得到6个子方程， 再加上位移与速度和姿态角与角速度关系式， 共计12个未知量12个方程， 可直接采用Matlab进行数值仿真。

6 数值仿真

6.1 俯仰运动仿真

为了研究无人机机翼的迟滞环对无人机运动的影响， 对无人机小迎角俯仰运动进行仿真。 该仿真初始条件为： 垂起固定翼无人机飞行高度120 m，  前飞速度为10 m/s， 四个垂起旋翼转速均为2 000 r/min， 前拉螺旋桨转速为1 500 r/min， 初始迎角为4°， 其余状态量初始时刻均为0。 有迟滞环修正的模型仿真曲线对应有迟滞环， 无迟滞环修正的对应无迟滞环， 在该飞行过程当中无横航向运动。

6.2 垂起改平飞过渡阶段仿真

对于垂起固定翼无人机而言， 其垂起改平飞的过渡阶段是数值仿真和研究的难点。 采用非定常气动力模型和未引入非定常气动力（即迎角小于15°时为定常气动模型， 迎角大于15°时， 机翼不产生气动力）的动力学模型进行对比。 仿真过程為， 在前40 s无人机进行垂直起飞， 四个垂起旋翼转速均为2 000 r/min， 并且保持不变， 40～45 s， 无人机四个垂起旋翼转速逐渐减为零， 两个前拉螺旋桨转速均由零增加为2 000 r/min， 45 s后， 四个垂起旋翼转速为零， 两个前拉螺旋桨转速保持不变， 直到60 s无人机前飞速度达到15 m/s， 无人机完成垂起改平飞的过渡过程。

6.3 盘旋运动仿真

无人机的盘旋运动是一个横航向与纵向耦合的运动， 无人机在空中执行任务时， 常会发生这样的耦合运动。 对垂起固定翼无人机的盘旋运动进行数值仿真以研究非定常气动力对其动力学的影响。 仿真初始条件为： 无人机前飞速度8 m/s， 迎角为2°， 无人机完成垂起改平飞过程后进行盘旋， 四个垂起旋翼转速为零， 两个前拉螺旋桨转速分别为1 000 r/min和2 000 r/min， 通过两个前拉螺旋桨的拉力差形成力矩进而发生盘旋运动， 其余状态量初始时刻均为零。

7 结  论

通过对垂起固定翼无人机动力学建模和非定常气动力的研究， 可以得到以下结论：

（1） 采用多体动力学基于凯恩方程能够系统和方便地建立垂起固定翼无人机的动力学模型， 通过迟滞环和ONERA方程能够建立垂起固定翼无人机的非定常气动力模型。

（2） 通过仿真曲线图4～6可以发现， 引入迟滞环将导致无人机俯仰振荡速率减小， 其主要原因在于： 由于无人机机翼的气动中心位于其重心之后， 当无人机做俯仰运动， 俯仰角速度为正值时，  机翼气动力起阻尼作用，  俯仰角速度为负值

时， 机翼气动力起促进作用， 其产生的阻尼作用比促进作用大， 因此， 整体呈现阻尼效果， 对于该类无人机进行俯仰运动时引入迟滞环进行修正是合理的。 引入迟滞环的本质其实是引入气流在流动过

程中所产生的非定常效应， 这种非定常效应不

仅在无人机做俯仰等纵向运动时产生阻尼效果， 而且当无人机进行盘旋等横航向与纵向耦合的运动时同样也将产生阻尼效果， 其主要原因在于引入迟滞环将导致无人机的力矩系数随运动而发生变化， 进而使得无人机在运动过程中产生较大的阻尼力矩。

（3） 通过仿真曲线图7～10可以发现， 两者仿真结果整体变化趋势一致， 但是变化过程有差异。 从图9可以明显发现， 无人机垂起速度引入非定常气动力与未引入有8%的差异， 引入非定常气动力模型的垂起速度更小， 在图9的动态过渡过程中， 即40～60 s， 两者也具有较大差异， 引入非定常气动力模型与未引入非定常气动力模型两者的最终上升速度存在40%的差异， 而上升速度是垂起改平飞过程中需要重点关注的物理量， 因此在控制

和设计过渡阶段时， 应以引入非定常气动力的动力学模型为基准， 这样才能保证无人机的安全性。 图10展示了引入非定常气动力和未引入非定常气动力模型的总能量曲线对比， 可以发现引入非定 常气动力模型的总能量更小，  其主要原因在于无

式中各变量的定义和具体的非定常动态流场数值方法可参考文献[11-12]。

人机在过渡阶段需要克服非定常气动力做功， 进而导致其总能量减小。 垂起改平飞过渡过程的总能量曲线整体变化光滑， 该无人机能够顺利完成垂起改平飞过程。

（4） 通过仿真曲线图11～12的对比可以得到两个结论： a. 垂起固定翼无人机在盘旋运动中引入非定常气动力与未引入非定常气动力仿真结果最终趋于一致， 证明了本文采用ONERA方程的准确性， 主要在于ONERA方程在小迎角情况下为线性函数， 其与传统气动力模型相差不大， 两者仿真结果在小迎角情况下是一致的; b. 通过图12可以发现， 垂起固定翼无人机在大侧滑角的情况下引入非定常气动力后其波动的振幅更小， 收敛速度更快。 盘旋运动是一个耦合运动， 在大侧滑角情

况下采用ONERA非线性方程建立了升力系数与阻力系数模型，  由于这种情况下的升力系数与阻力系数与未引入非定常气动力的情况相比， 其值更小， 因此其波动的振幅会更小。 通过图11与图12的对比可以发现， 盘旋运动过程中大侧滑角情况其实就对应着大迎角情况。

参考文献：

[1]  Kim S K，  Tilbury D M. Mathematical Modeling and Experimental Identification of a Model Helicopter[C]∥AIAA Modeling and Simulation Technologies Conference  and Enhibit，  AIAA-98-4357，  Boston， 1998，  203-213.

[2]  DeTore J，  Conway S. Technology Needs for High Speed Rotorcraft（3）， NASA-CR-177592[R].   1991.

[3]  李家乐，  王正平. 基于Lagrange方法的单旋翼飞行器动力学建模[J]. 飞行力学，  2016，  34（4）： 15-18.

Li Jiale，  Wang Zhengping. Dynamics Modeling for Monowing Rotorcraft Using Lagrange Method[J]. Flight Dynamics，  2016，  34（4）： 15-18.（in Chinese）

[4]  李家乐. 弧翼飞行器飞行控制与仿真[D]. 西安： 西北工业大学，  2016： 32-43.

Li Jiale. Research on Control and Simulation for a Cambered Wing Vehicle[D]. Xian： Northwestern Polytechnical University，  2016： 32-43.（in Chinese）

[5]  Hogan F R，  Forbes J R. Modeling of Spherical Robots Rolling on Generic Surfaces[J]. Multibody System Dynamics，  2015，  35（1）： 91-109.

[6]  曹芸芸. 倾转旋翼飞行器飞行动力学数学建模方法研究[D]. 南京： 南京航空航天大学，  2012： 29-63.

Cao Yunyun. Research on Mathematical Modeling Method for Tilt Rotor Aircraft Flight Dynamics[D]. Nanjing： Nanjing University of Aeronautics and Astronautics，  2012： 29-63.（in Chinese）

[7]  Kane T R，  Likins P W，  Levinson D A. Spacecraft Dynamics[M]. New York： McGrawHill Book Company，  1983.

[8]  Zhao Zhenjun，  Ren Gexue. Multibody Dynamic Approach of Flight Dynamics and Nonlinear Aeroelasticity of Flexible Aircraft [J]. AIAA Journal，  2011，  49 （1）： 41-54.

[9]  喬宇航，  马东立，  李陟. 螺旋桨/机翼相互干扰的非定常数值模拟[J]. 航空动力学报，  2015， 30（6）： 1366-1373.

Qiao Yuhang，  Ma Dongli，  Li Zhi. Unsteady Numerical Simulation of Propeller/Wing Interaction[J]. Journal of Aerospace Power，  2015， 30（6）： 1366-1373.（in Chinese）

[10]  袁先旭，  张涵信，  谢昱飞. 基于CFD方法的俯仰静、 动导数数值计算[J]. 空气动力学学报，  2005，  23（4）： 458-463.

Yuan Xianxu，  Zhang Hanxin，  Xie Yufei. The Pitching Static/Dynamic Derivatives Computation Based on CFD Methods[J]. Acta Aerodynamica Sinica，  2005，  23（4）： 458-463. （in Chinese）

[11]  East R A，  Hutt G R. Comparison of Predictions and Experimental Data for Hypersonic Pitching Motion Stability [J]. Journal of Spacecraft and Rockets，  1998，  25（3）：225-233.

[12]  刘罗成，  徐敏， 姚伟刚， 等. 翼型纵向俯仰非定常气动力研究[J]. 飞行力学，  2011，  29（4）： 28-31.

Liu Luocheng，  Xu Min，  Yao Weigang，  et al. Study of Hypersonic Airfoil Pitch Oscillation Unsteady Aerodynamic Force[J]. Flight Dynamics，  2011，  29（4）： 28-31.（in Chinese）

[13]  Kim T， Dugundji J. Nonlinear Large Amplitude Aeroelastic Behavior of Composite Rotor Blades[J]. AIAA Journal，  1993，  31（8）： 1489-1497.

[14]  林學海，  任勇生. 基于ONERA气动力模型的风力机叶片颤振时域分析[J]. 山东科技大学学报： 自然科学版，  2009，  28（3）： 56-60.

Lin Xuehai，  Ren Yongsheng. Analysis of Flutter TimeDomain for Blades of Wind Turbine Based on the ONERA Aerodynamic Model[J]. Journal of Shandong University of Science and Technology： Natural Science，  2009，  28（3）： 56-60.（in Chinese）

[15]  Forbes J R，  Barfoot T D，  Damaren C J. Dynamic Modeling and Stability Analysis of a PowerGenerating Tumbleweed Rover[J]. Multibody System Dynamics，  2010，  24（4）： 413-439.

[16]  蒋达. 新型无刷直流电动机建模与仿真研究[D]. 南宁： 广西大学，  2012：  17-19.

Jiang Da. Research on Modeling and Simulation of New Type Brushless DC Motor[D]. Nanning：  Guangxi University，  2012：  17-19.（in Chinese）

[17]  夏丹， 陈维山， 刘军考， 等. 基于Kane方法的仿鱼机器人波状游动的动力学建模[J]. 机械工程学报，  2009， 45（6）： 41-49.

Xia Dan，  Chen Weishan，  Liu Junkao，  et al. Dynamic Modeling of a Fishlike Robot with Undulatory Motion Based on Kanes Method[J]. Journal of Mechanical Engineering，  2009，  45（6）： 41-49.（in Chinese）

[18]  Tarn T J，  Shoults G A，  Yang S P. A Dynamic Model of Underwater Vehicle with a Robotics Manipulator Using Kanes Method[J]. Autonomous Robots，  1996，  3（2-3）： 269-283.