随着“大鹏”、“蓝鲸”等无人机的问世, 多旋翼固定翼无人机已成为目前研究的热点, 由于其兼顾垂直起降和前飞巡航性能, 具有较大市场前景。然而对于多旋翼固定翼无人机的动力学建模问题到目前为止仍然没有解决。国外^[1-2]和国内^[3]分别采用牛顿-欧拉方程和拉格朗日方程建立固定翼和旋翼无人机的动力学模型, 以上文献均是将无人机看成完整刚体进行分析。对于多旋翼固定翼这类特殊无人机, 其在运动过程当中与传统无人机存在以下区别, 首先其旋翼数量较多且旋翼叶片较大, 因此旋翼对于无人机的影响不可忽略; 其次, 对于该无人机的动力学建模方法而言, 拉格朗日方程与牛顿-欧拉方程相比, 其主要采用能量交换的方式建立动力学模型, 建模更为系统, 且拉格朗日方程直接通过旋翼的转动动能引入旋翼对无人机整体的影响, 因此选用拉格朗日方程作为此类无人机的建模基本方程具有优势。为了解决上述问题, 建立准确的无人机动力学模型, 采用多体动力学思路, 将该无人机划分为左右机翼、机体、左右平尾、旋翼、垂尾、舵面共计16个刚体, 针对每个刚体单独进行动力学分析, 通过虚功和第二类拉格朗日方程^[4], 建立该无人机动力学模型, 该建模方法的主要创新点在于, 引入圆柱扰流和平板扰流建立垂直起降过程的无人机非定常动力学模型, 采用四元数描述无人机质心姿态引入拉格朗日乘子表征各刚体之间的相对约束关系, 基于多体动力学思路和拉格朗日方程建立无人机动力学模型。

1 多旋翼固定翼无人机运动学分析 1.1 坐标系建立

建模整体思路如图 1所示:

该无人机具体构型如图 2所示, 采用多旋翼与固定翼相结合的方式, 其中有4个垂起旋翼和2个前拉旋翼, 无人机采用“T”形尾。基于多体动力学^[5-7], 将该无人机划分为左机翼、右机翼、机体、左平尾、右平尾、垂尾、升降舵、左副翼、右副翼、方向舵、1号旋翼、2号旋翼、3号旋翼、4号旋翼、5号旋翼、6号旋翼共计16个刚体, 分别以每个刚体的质心为坐标原点建立右手坐标系。选取F_e作为基准坐标系, F_l为气流坐标系。

1.2 广义坐标选取

由于四元数没有奇点且其自带约束^[8]有利于消去拉格朗日乘子, 因此选取四元数作为描述无人机质心姿态的广义坐标, 其形式如下

(1)

式中，q为无人机质心广义坐标, r^be用于描述无人机质心的位置, q^be为用于描述无人机质心姿态的四元数, 四元数q^be的具体含义为:基准坐标系F_e通过绕向量a旋转δ角度转换为机体坐标系F_b, a₀, a₁, a₂, a₃为四元数q^be的分量。

2 多旋翼固定翼无人机动力学分析 2.1 机翼、舵面、平尾动力学模型

机翼、舵面、平尾的建模方式相同, 在此仅以左机翼为例给出建模思路, 左机翼气动力模型由两部分构成, 即正常飞行的准定常气动力模型和垂直起降的非定常气动力模型, 先给出左机翼准定常气动力模型

(2)

式中，C_el为气流坐标系向基准坐标系转换矩阵, V_z为左机翼质心处的合速度, S_z为左机翼面积, C_D, C_y, C_L均为力系数, 采用牛顿插值方法求得。

左机翼非定常气动力模型

(3)

式中，w为垂起速度, C_{L_f}和C_{D_f}为无人机非定常气动力系数, 采用平板扰流的方式给出, 一般只有当无人机垂起或者垂降时, 其值才不为零。因此左机翼总力模型为

(4)

对于无人机左机翼所产生的力矩同样由两部分产生, 即无人机左机翼质心处气动力相对于无人机机体质心所产生的力矩和无人机左机翼在运动过程当中所产生的非定常气动力矩。由无人机左机翼质心气动力所产生的力矩主要在下文通过虚功的形式引入。无人机左机翼非定常气动力矩主要为横向阻尼力矩、航向交叉力矩、航向阻尼力矩, 横向交叉力矩, 其具体形式如下

(5)

式中，c为无人机左机翼平均气动弦长, p, r为无人机角速度分量, C_lp, C_np, C_nr, C_lr为无人机动导数。

2.2 机体动力学模型

无人机机体在飞行时主要受到非定常气动力, 采用圆柱扰流近似求得无人机机体在运动过程中的升力系数与阻力系数, 并通过力系数与迎角和机体质心合速度的耦合得到无人机机体的非定常气动力模型, 其具体形式如下所示

(6)

式中，C_{D_y}为圆柱扰流阻力系数, 当无人机垂起时S_{j_1}为零, 当无人机前飞时其值为机体纵向截面积, 当无人机垂起时S_{j_2}为机体横向截面积, 当无人机前飞时其值为零, 由于无人机机体所产生的力矩较小, 因此在建模过程当中予以忽略。

2.3 旋翼动力学模型

采用叶素理论所推公式作为旋翼动力学模型的基本公式, 将滑流理论作为动态入流模型, 建立旋翼动力学模型, 其具体形式如下

(7)

(8)

其中各力和力矩系数形式如下

(9)

式中，C_ea为旋翼坐标系向基准坐标系的转换矩阵, C_ba为旋翼坐标系向机体坐标系的转换矩阵, k为桨叶片数, a_∞为桨叶翼型升力系数, μ为前进比, θ₀为桨叶根部安装角, θ₁为桨叶扭度, v₀为旋翼平均诱导速度, λ₀为旋翼流入比, C_x为桨叶翼型阻力系数, m_c为旋翼质量, ω为旋翼转速, R为桨叶半径, M_x, M_y分别为螺旋桨桨毂力矩在其X轴和Y轴中的投影。

为了方便后续进行无人机虚功的求解, 首先求得无人机整机非定常力矩之和, 应当注意的是各刚体的力模型都是基于基准坐标系建立的, 力矩模型都是基于无人机机体坐标系而建立的。

(10)

式中, M_sum为无人机非定常合力矩, M_z, M_y分别为无人机左右机翼质心处的非定常气动力矩模型, M_pz, M_py分别为无人机左右平尾质心处的非定常气动力矩模型, 为无人机旋翼1到6号质心处的非定常气动力矩模型。

3 拉格朗日方程 3.1 虚功求解

通过虚位移的方法求得虚功, 将各刚体虚功叠加得到无人机广义力和广义力矩模型, 在此以无人机机体和1号旋翼为例进行虚功求解, 其余刚体类似求得即可。

对于无人机机体而言, 其虚功分为两部分, 即无人机机体气动力模型所做虚功和无人机非定常合力矩模型所做虚功, 其中无人机机体质心到基准坐标系原点的实位移为r^be, 求得其虚位移为δr^be, 通过虚位移求得其气动力所做虚功为

(11)

同理求得其合力矩所做虚功为

(12)

式中，S_b^be为四元数求偏导后的分量, 对于1号旋翼而言, 其非定常气动力矩所做虚功已在M_sum中求解, 其气动力模型所做虚功按上文方式进行求解, 旋翼质心到基准坐标系原点的实位移为r^ae, 其虚位移为δr^ae, 所做虚功具体形式如下

(13)

将所有虚功进行求和运算, 便可以得到无人机质心处的广义力和广义力矩模型为

(14)

式中，F为无人机整机广义力和力矩模型, δW_z, δW_y分别为无人机左右机翼质心气动力所做虚功, δW_p, δW_q分别为左右平尾质心气动力所做虚功, 为无人机1到6号旋翼所做虚功, δW_k, δW_j分别为方向舵和升降舵质心气动力所做虚功, δW_h, δW_i分别为无人机左右副翼质心气动力所做虚功, δW_m为垂尾质心气动力所做虚功。

3.2 拉格朗日函数

为了建立该无人机的动力学模型, 首先应该进行拉格朗日函数的求解, 其主要由无人机的动能与势能所构成, 其中无人机的动能主要为旋翼动能和无人机整机动能之和, 旋翼动能的求解可以分为两部分, 即跟随无人机一起运动时的动能和单独旋转时的转动动能^[9], 无人机动能的具体形式如下

(15)

(16)

(17)

(18)

式中，v^be为无人机质心速度, ω^be为无人机质心角速度, m_sum为无人机各刚体质量和, J_sum为无人机各刚体转动惯量之和, J_a, J_c, J_d, J_n, J_f, J_g分别为各旋翼转动惯量, ω_a, ω_c, ω_d, ω_n, ω_f, ω_g分别为各旋翼自身转速, A为6个旋翼转动动能之和。

对于该无人机, 其势能直接进行求解即可

(19)

因此拉格朗日函数L为

(20)

3.3 动力学模型建立

将上述求得的各项带入下面的拉格朗日方程当中

(21)

式中，λ为拉格朗日乘子, 带入四元数约束后可以直接消掉, F为无人机质心的广义力和广义力矩矩阵, Ξ为四元数自带约束, 其中拉格朗日乘子主要表征各刚体之间的约束。

4 仿真与实验

采用MATLAB基于四阶龙格库塔方法对该无人机的动力学模型进行数值仿真, 将仿真结果与牛顿-欧拉方程六自由度模型进行对比, 仿真初始条件为:旋翼1-6号的转速分别为1 500, 1 000, 1 500, 1 500, 1 000, 1 500 r/min, 无人机飞行高度为100 m, 前飞速度为6 m/s, 初始迎角为3°, 其余量的初始值均为零, 其仿真结果如下。在该飞行过程当中无人机没有侧向运动, 因此其侧滑角、侧向速度等量均为零。

图 3主要描述了无人机迎角的变化过程, 无人机迎角初始时刻为3°, 并未达到稳定, 因此无人机迎角一开始便快速由3°变化到5°, 而后再逐步趋于稳定状态; 图 7主要表示无人机总能量变化曲线, 拉格朗日方程与牛顿-欧拉方程仿真结果差异不大, 符合逻辑。通过图 3~6数值仿真结果的对比, 可以发现拉格朗日方程与牛顿-欧拉方程仿真结果在120 s之前两者变化趋势一致, 但动态变化过程不一致, 即拉格朗日方程以一种波动的形式趋于稳定, 而牛顿-欧拉方程则不存在波动, 在120 s之后两者仿真曲线几乎相同。接下来对这种现象进行分析, 首先牛顿-欧拉方程与拉格朗日方程相比, 两者在建模过程当中所忽略和着重考虑的因素略有不同, 但两者最终都是回归到牛顿第二定律和转动惯量定律上, 因此最终两者仿真结果趋于一致。对于两者仿真结果为何在前120 s趋势一致而动态变化过程不一致, 其主要原因在于以下三方面, 从建模基本方程而言, 拉格朗日方程是采用四元数描述无人机姿态, 由于四元数自带约束将产生拉格朗日乘子, 而牛顿-欧拉方程则是采用欧拉角描述无人机姿态, 欧拉角不存在约束, 因此由于处理四元数约束将导致两者动态变化过程有所不同; 其次拉格朗日方程在处理无人机力矩时, 是将其化为两部分, 即非定常气动力矩直接给出, 气动力所产生力矩通过虚功形式引入, 而牛顿-欧拉方程更倾向于直接将力矩引入无人机质心, 因此拉格朗日方程在考虑无人机力矩时更为全面; 最后, 从旋翼对无人机整机的影响考虑, 拉格朗日方程以动能的形式引入旋翼对无人机整机的影响, 而牛顿-欧拉方程则是通过力和力矩的形式引入旋翼的影响。因此拉格朗日方程对于此类特殊无人机的动力学建模问题具有较大优势。

对多旋翼固定翼无人机而言, 更为关注的过程应当是垂起改平飞过程, 接下来采用如下方式进行该无人机的垂起改平飞实验, 并将仿真结果与实验进行对比, 过渡策略为多旋翼无人机0~28 s, 4个垂起旋翼转速由1 500 r/min线性增加到2 000 r/min, 2个前拉旋翼无转速, 28 s开始进行过渡, 从28~40 s, 4个垂起旋翼转速线性减少为0。2个前拉旋翼转速线性增加到2 000 r/min, 40~50 s, 2个前拉旋翼转速保持不变, 4个垂起旋翼转速均为零, 多旋翼固定翼无人机前飞速度达到15 m/s, 垂起改平飞顺利结束, 其具体过程如下:

图 8~11描述了无人机由垂起改平飞的整个过程, 其中无人机在前28 s内进行垂直起飞, 当飞行高度达到171 m后, 无人机由垂直起飞开始改向前飞状态; 由图 9可以观察得到无人机垂起改平飞过程高度损失约2 m, 高度损失不大, 能够保证飞行的稳定性; 由图 10可知, 无人机在前30 s内, 其迎角均为-90°, 在30 s后由于2个前拉旋翼转速的增加和4个垂起旋翼转速的减小, 其迎角变化曲线有凸起, 最后无人机迎角再次慢慢趋于稳定, 无人机改平飞过程顺利完成; 图 11表示了无人机垂起改平飞过程当中总能量的变化, 其中总能量曲线变化较为平缓, 说明该无人机垂起改平飞过程较为顺利, 所制定的过渡策略能完成目标。

图 8~11同时也描述了拉格朗日方程与牛顿-欧拉方程在无人机垂起改平飞过程当中仿真结果的差异, 通过图 8, 10, 11的对比可以发现, 拉格朗日方程与牛顿-欧拉方程仿真结果具有较好的一致性, 能够较为准确地描述无人机垂起改平飞的过渡过程; 通过图 9可以发现, 拉格朗日方程、牛顿-欧拉方程以及实验数据, 三者具有一定差异, 产生差异的主要原因在于真实情况下无人机由垂起改平飞的过程中将受到不同旋翼之间的干扰, 由于无人机旋翼数量较多, 因此忽略旋翼间的干扰阻力将导致数值仿真结果与实验结果的差异。

5 结论

通过该无人机动力学建模与仿真结果分析可以得到如下结论：

(1) 本文主要对拉格朗日方程动力学建模的整体思路和建模方法进行归纳, 其中拉格朗日方程与牛顿欧拉方程相比, 其主要是通过能量交换的方式进行动力学建模, 通过虚功引入广义力和力矩, 通过拉格朗日乘子处理刚体之间的约束, 这种建模方法能够更为全面地考虑多旋翼对于无人机整体的影响, 更适合于多旋翼固定翼这类结构复杂的无人机；

(2) 通过图 3~7拉格朗日方程多体动力学建模与牛顿-欧拉方程动力学建模仿真结果的对比可以发现, 两者的最终结果是一致的, 但两者由仿真初始条件过渡到稳定状态的过程不一致, 拉格朗日方程主要是通过波动的形式变化到稳定, 而牛顿-欧拉方程不存在波动, 其主要原因在于两类方程的基准方程都是牛顿第二定律和转动惯量定律, 但在方程推导和建模的过程当中, 拉格朗日方程基于能量交换的方式更为全面地考虑了旋翼和非定常气动力矩对于无人机的影响, 这将导致其存在上述差异。

(3) 通过图 8~11无人机垂起改平飞过渡阶段实验与仿真结果对比, 上文所提过渡策略能完成无人机垂起改平飞过程, 且该动力学模型能够准确描述无人机动态变化过程, 通过仿真曲线发现过渡阶段无人机高度下降约2 m, 能量和迎角曲线变化相对光滑, 如果对无人机高度有严格要求, 可以进一步优化过渡策略, 减小过渡过程的高度下降。