基于凯恩方程的无人机伞降回收动力学建模与仿真

吴翰，王正平，周洲，王睿

（西北工业大学 无人机特种技术国防科技重点实验室，西安710065）

随着美国 X47b、欧洲“神经元”等无人机的问世，飞翼类[1-2]无人机已逐步成为国内外研究的热点，飞翼类无人机具有升阻比高等特性，这些特性在无人机起飞时具有较大优势，但飞翼无人机在着陆时尤其对于低翼载飞翼无人机而言，其着陆轨迹容易受到风力的影响，而飞翼无人机本身又不是很稳定，因此将会对此类无人机的安全回收造成阻碍。目前无人机常采用3种回收方式，即伞降回收、滑跑着陆回收和拦网回收，为了将飞翼无人机机翼面积较大的特性转化为其回收时的优势，即将飞翼无人机的机翼当成阻力板以减小飞翼无人机的着陆速度，选取伞降方式对此类无人机进行回收。

降落伞在工程领域应用较多，目前已有较为成熟的气动研究[3]和附加质量[4]研究，但是关于伞降回收系统动力学建模的研究更多集中在航天领域，国外文献[5-7]和国内文献[8]提出了两种建立航天器伞降回收动力学模型的方法。由于无人机的伞降回收与航天器的伞降回收具有较大差异，因此上述模型对于无人机的伞降回收有不适用之处。针对这样的情况提出以下需要解决的问题：①无人机在伞降回收过程中和降落伞一直处于动态平衡状态，由于两者之间的耦合关系较为复杂，因此很难建立精确的无人机伞降回收动力学模型，如何采用较为简单的方式得到能表征无人机伞降回收运动趋势的动力学模型，值得研究。②目前大部分关于降落伞回收的动力学模型，都是将回收物看成单刚体进行处理，对于低翼载飞翼无人机的伞降回收，飞翼无人机对整个伞降回收系统有一定的影响，因此能否采用新的建模方法，更为精确地引入无人机的影响，值得研究。③大型降落伞在充气阶段，其质心会发生较为明显的变化，因此在建立该过程降落伞动力学模型时需要大量的实验数据进行拟合，较为复杂；对于无人机伞降回收，特别是小型无人机，其采用的降落伞一般较小，因此降落伞在拉直以后的充气过程中质心位置变化不大，能否采用更为简单的方式建立此类降落伞充气阶段的动力学模型，值得研究。④对于无人机伞降回收，无人机由稳定平飞过渡到稳定下降，其实是一个迎角变化很大的运动过程，如果直接采用CFD进行数值模拟得到该阶段无人机的非定常气动力模型将是较为复杂的，能否基于合理的假设快速建立该阶段的非定常气动力模型，以满足工程的需要，值得研究。关于以上问题的解决方法与思路将是本文的创新点所在。

为了解决上述问题，建立准确的无人机伞降回收动力学模型，首先，采用平板绕流与迎角之间的耦合关系建立无人机从平飞过渡到稳定下降阶段的非定常气动力和力矩模型；其次，通过降落伞的阻力面积随其充气时间的变化，建立降落伞充气过程中的动力学模型；最终，基于多体动力学和凯恩方程[9-10]建立低翼载飞翼无人机伞降回收六自由度模型，并分析了伞降回收位置海拔高度和风力对伞降回收的影响。需要强调的是，本模型是将各系统的力和力矩引入整个伞降回收系统的质心建立动力学模型，该模型能够描述无人机伞降回收的趋势，但无法描述降落伞与无人机之间的相对运动，该模型主要描述的是降落伞从完全拉直到和无人机一起稳定下降的运动过程。

1 伞降回收系统运动学描述

1.1 坐标系建立

低翼载飞翼无人机伞降回收系统主要由低翼载飞翼无人机和回收降落伞两部分组成，如图1所示。

图1 伞降回收坐标系示意图Fig.1 Schematic diagram of parachute recovery coordinate systems

将伞降回收系统分为降落伞和无人机，基于多体动力学的思路[11-12]将无人机离散为无人机左机翼d、右机翼c和机身 a的多体系统，其中由于垂尾、舵面以及动力装置在回收过程中影响不大，将其并入无人机机身进行分析。由于无人机伞降回收的降落伞较小，因此可以假设降落伞拉直后，这一时变系统在充气过程中的质心与其气动合力作用点位于同一点。以各体质心为原点建立右手坐标系。首先基于机身坐标系Fa求得整个伞降回收系统质心b的位置，然后基于整个伞降系统质心b进行建模。伞降回收系统质心b的位置为

式中：lx、ly和lz分别为降落伞质心 g到机身质心a的位移在机身坐标系X、Y和Z轴方向的投影；xda、xca分别为左右机翼质心 d、c到机身质心 a的位移在机身坐标系X轴方向的投影；yda和yca分别为左右机翼质心d和c到机身质心a的位移在机身坐标系Y轴方向的分量；zda和zca分别为左右机翼质心d和c到机身质心a的位移在机身坐标系Z轴方向的分量；mg、ma、md和 mc分别为降落伞、机身和左右机翼的质量。需要强调的是，伞降回收系统质心位置b会随着降落伞的充气而发生变化，在实际情况中由于降落伞质量相对于无人机而言较小，因此降落伞质量的变化对全机质心的影响较小。本文最终采用迭代的方式进行动力学模型的仿真，因此可以在每一时刻采用式（1）得到伞降回收系统质心，然后将力和力矩引入该质心处进行动力学建模。

1.2 广义坐标与速度矩阵建立

选取四元数[13]作为描述伞降回收系统质心姿态的广义坐标，具体形式为

式中：φ、θ和ψ为伞降回收系统质心处的姿态角；a0、a1、a2和 a3为四元数分量；xb、yb和 zb为伞降回收系统质心位置。

接下来1.1节文求得的伞降系统质心b为基准建立降落伞以及各体质心处的线速度、角速度、线加速度和角加速度矩阵，各体质心处的线速度和角速度矩阵具体形式为

式中：r为体与体质心间的位移向量；V为各体质心处的线速度；ω为各体质心处的角速度；下标c、d、a和g分别代表右机翼、左机翼、机身以及降落伞，无特别说明，下文一致。除此之外需要强调的是，各体质心的角速度近视是与伞降回收系统质心处角速度一致的，但各体质心处的角加速度不一定相同，具体公式由下文给出。

同理各体质心处的线加速度和角加速度矩阵分别为

1.3 质量和附加质量描述

接下来对伞降回收系统各体的质量、转动惯量、附加质量和附加质量惯性矩进行介绍。

图2 降落伞充气过程容积变化示意图Fig.2 Schematic diagram of parachute volume variation during inflation process

图2展示了降落伞整个充气过程中的容积变化，其实降落伞在充气过程当中是一个时变系统，其时变性主要由两部分体现，即其质心处的力和力矩随充气时间而发生变化；其质量、转动惯量、附加质量和附加质量惯性矩随充气时间而变化。其质心处的力和力矩随充气时间的变化公式将在后文给出，在此先给出其质量、转动惯量、附加质量和附加质量惯性矩随充气时间的变化公式。IR为降落伞与无人机的连接点到降落伞伞衣边缘的距离；R0为降落伞伞衣顶端到伞衣底部的弧线距离。

降落伞容积、质量和转动惯量随充气时间的变化公式为[14]

式中：Δ为降落伞的容积；ρ为大气密度；Dp为降落伞投影直径；t为降落伞充气时间；tm1为降落伞初始充气阶段结束时的充气时间；tm为降落伞主充气阶段结束时的充气时间；ht为降落伞伞衣未充满部分的高度；din为降落伞伞衣未充满部分底部的直径；m1和m2分别为降落伞伞衣和伞绳的质量；mg为降落伞质量矩阵；Jg为降落伞转动惯量矩阵；Jx1、Jy1和Jz1为降落伞伞衣转动惯量矩阵的分量；Jx2、Jy2和Jz2为降落伞伞绳转动惯量的矩阵分量。

降落伞附加质量和附加质量惯性矩随充气时间的变化关系式为

式中：ρ为大气密度；kii和 kjj均为降落伞附加质量和附加质量惯性矩系数；If为降落伞所带动周围运动空气的体积；Ng和 Ig分别为降落伞的附加质量和附加质量惯性矩矩阵；Aii和 Ajj分别为降落伞附加质量和附加质量惯性矩矩阵分量。无人机各体的质量和转动惯量矩阵与常规形式一样，因此不再进行介绍。

1.4 广义惯性力和广义惯性力矩模型建立

将上述得到各体的线加速度矩阵、角加速度矩阵以及各体的质量和转动惯量矩阵代入广义惯性力和广义惯性力矩方程当中，可以得到无人机伞降回收系统的广义惯性力和广义惯性力矩模型[15]为

式中：mj为各体的质量矩阵；Nj和 Ij分别为降落伞附加质量矩阵和附加质量惯性矩矩阵，只有当j取g时，其值才不为零；Jj为各体转动惯量矩阵；分别为各体质心处的偏线速度和偏角速度矩阵。

2 伞降回收系统动力学描述

2.1 降落伞动力学模型建立

对于降落伞充气过程的气动力和力矩模型，目前大多采用半实验半理论的方法进行建模。为了建立降落伞在充气过程中的动力学模型，首先将降落伞从拉直到稳定下降划分为3个阶段，如图3所示，即初始充气阶段、主充气阶段和稳定下降阶段，求得降落伞每个阶段所需充气时间，根据降落伞充气时间表征降落伞阻力面积，最终基于降落伞阻力面积随充气时间的变化公式来建立降落伞的动力学模型。

图3 无人机伞降回收示意图Fig.3 Schematic diagram of UAV parachute recovery

采用充气时间法[16]得到降落伞前2个充气阶段所用时间分别为

式中：λ和k分别为初始充气阶段与主充气阶段的修正因子；D0为降落伞名义直径；Vl为伞衣拉直速度；Vp为主充气过程中的平均速度。

通过式（14）得到的充气时间可建立降落伞阻力面积随充气时间的变化公式为

式中：（CA）1为降落伞初始充气阶段完成时的阻力面积；（CA）s为降落伞主充气阶段结束时的阻力面积；CT为降落伞轴向力系数；A0为降落伞名义面积。

式中：Cbg表示坐标系 g向坐标系 b的转换矩阵；vg、ug和wg为降落伞沿 X轴、Y轴和 Z轴的速度分量；Vg为降落伞质心合速度；CN为降落伞法向力系数；mg为降落伞质量。

2.2 无人机动力学模型建立

基于多体动力学，将无人机进行离散化处理，由于轻型飞翼无人机的特殊性，左右机翼以及机身的建模思路一致只是力和力矩系数稍有不同，下面以右机翼为例建立其动力学模型。

首先按照无人机的失速迎角将其力和力矩划分为两部分，当无人机迎角小于失速迎角时，其气动力和力矩模型与正常飞行时一样，本文不再给出；当其迎角大于失速迎角时，由于降落伞的作用，无人机的迎角会很快达到大迎角状态，此时流过飞翼无人机机翼的气流将会发生流动分离，故可以采用平板扰流系数[17-18]与此阶段迎角的关系来建立这一阶段的无人机气动力模型，具体右机翼力模型为

式中：Cp为平板扰流系数；Cbc和Cbe分别为右机翼坐标系和地面坐标系向系统质心坐标系的转换矩阵；αc为无人机右机翼质心处的迎角；Sc为右机翼面积；mc为右机翼质量；gc为重力加速度。

无人机在运动过程当中的力矩模型主要由3部分叠加产生，即无人机运动过程当中各体所产生的非定常力矩、各体质心处气动力相对伞降回收系统质心所产生的力矩和各体气动合力作用点处的气动力相对于各体质心所产生的力矩。无人机各刚体质心处的气动力相对伞降回收系统质心的力矩，由其力模型乘以其质心偏速度矩阵得到，无人机右机翼非定常力矩模型和右机翼气动合力作用点处气动力相对于右机翼质心所产生的力矩之和为

式中：cc和l分别为右机翼展长与弦长；pc和rc为右机翼角速度分量；Clp、Clc、Cnp、Cnr、Cmθ和Cmθ·均为无人机动导数，其表征的是无人机横向、航向和纵向阻尼特性。同理可得无人机机身质心处的力和非定常力矩模型分别为Fa和Ma；左机翼质心处的力和非定常力矩模型为 Fd和 Md；Fc_q为无人机右机翼的气动合力作用点处的气动力，（xc_q，yc_q，zc_q）为无人机右机翼气动合力作用点到无人机右机翼质心处位移的分量。

2.3 广义主动力和广义主动力矩模型建立

将上述得到各体的力和非定常力矩模型导入广义主动力和广义主动力矩方程，可得低翼载飞翼无人机伞降回收系统的广义主动力和广义主动力矩模型为

式中：Fj为伞降系统各体质心处的力模型；Mj为伞降系统各体质心处的力矩模型；需要注意的是伞降回收系统的广义主动力矩模型是由两部分构成的，一部分是各体自身运动绕其自身质心所产生的非定常力矩和各体气动合力作用点处气动力相对于各体质心所产生的力矩之和，已由式（20）和式（21）直接给出。另一部分为各体质心广义主动力相对于伞降回收系统质心所产生的力矩，该部分由这一项引入。

3 无人机伞降回收动力学模型建立

凯恩方程由动力学普遍方程演化而来，其形式为

式中：Fi为系统广义主动力；mi为各体质量；δri为力作用点到所选参考点的虚位移；ai为各体质心处的加速度，将式（23）化简可得

因此可以得到凯恩方程的具体形式为

式中：K*称为广义惯性力和广义惯性力矩模型；K称为广义主动力和广义主动力矩模型。凯恩方程建模的基准为：系统的广义主动力、广义主动力矩与系统的广义惯性力、广义惯性力矩之和为零。

将2.1节和2.2节所求伞降回收系统广义主动力、广义主动力矩模型和广义惯性力、广义惯性力矩模型代入凯恩方程[19]可得

式（26）即为飞翼无人机伞降回收系统六自由动力学模型，式中共有6个方程，13个变量，再加上伞降回收系统质心角速度与四元数的关系式以及质心位移与速度的关系式，共计13个方程，13个变量，可直接采用MATLAB进行数值仿真。

4 伞降回收仿真实验

基于四阶龙格库塔方法进行数值仿真并与实验结果进行对比。本次实验过程为：海拔高度1 300 m，无人机完成由稳定前飞到开伞，再到无人机和降落伞共同稳定下降的实验。为了验证本文模型的准确性，所选取实验数据的初始时刻为：降落伞已经被拉直并且刚好开始充气，无人机飞行高度56 m，其前飞速度11 m／s，同时具有1 m／s的上升速度。在此对仿真结果与实验结果的对比做以下几点说明：①由于建模所用坐标系为欧美坐标系，为了与实验结果进行对比，沿z轴方向的物理量在仿真结果展示时均乘以-1。②由于实验仅得到无人机数据，因此仿真结果已经通过式（4）和式（5），由伞降回收系统质心转换到无人机质心处，方便进行对比。③由于实验得到的数据是基于地面坐标系的，因此仿真数据以由无人机机体坐标系转换到地面坐标系。具体仿真结果与实验结果对比如图4所示。

通过图4的对比可以发现：数值仿真与实验结果吻合较好，验证了本文模型的准确性。但是由开伞到稳定下降阶段的动态变化过程中，两者略有差异，归结其原因主要在于：①采用平板扰流建立无人机该阶段的非定常气动力模型，虽然较为便利，但是精度却相对不足。②所建模型无法表征该阶段降落伞与无人机之间的相对运动关系，而降落伞与无人机之间的相对运动关系，会进一步影响降落伞和无人机所受气动力的大小，进而影响该阶段的仿真精度。正是上述原因导致了该阶段的数值仿真与实验结果的差异，进一步优化模型时可从这两方面入手。

图4 无人机俯仰角、前飞速度、下降速度和飞行高度仿真结果与实验结果对比Fig.4 Comparison of UAV angle of pitch，forward velocity，descent velocity and flight height between simulation and experimental results

5 海拔高度与风力对伞降回收影响

无人机伞降回收位置的海拔高度以及伞降过程中的风力都将对无人机的伞降回收产生一定的影响，接下来引入海拔高度和风力对无人机伞降回收系统的影响，并对其进行分析。

为了得到无人机稳定下降速度随回收位置海拔高度的变化曲线，从海平面开始，海拔高度每增加500 m进行一次伞降仿真，通过数值仿真得到无人机最终稳定下降速度随海拔高度的变化曲线如图5所示。

通过图5可以发现：无人机伞降回收系统的稳定下降速度随回收位置海拔高度的增加而增大，其主要原因是海拔的变化将影响密度的变化，进而影响降落伞与无人机的气动力。海拔5 500 m与海拔0 m相比，其稳定下降速度由4.7 m／s变为了 5.51 m／s，增加了约 17%，这意味着随着海拔高度的增加，将对无人机伞降回收的安全性提出更高要求。

无人机在伞降回收过程中常会受到风力的影响，因此通过数值仿真得到无人机伞降回收系统在降落伞充气完成后（4 s后），分别受到1 m／s的顺风、逆风、正侧风和逆侧风的影响，无人机伞降回收速度变化的曲线，该仿真的初始条件为：海拔2 500 m，开伞高度 100 m，无人机前飞速度9 m／s，同时具有 1 m／s的上升速度，0 s时刻降落伞已经被拉直并且开始充气，其中规定沿Y轴正方向的风为正侧风，沿 Y轴反方向的风为逆侧风，沿X轴正向的为顺风，沿X轴反向为逆风，具体仿真结果如图6所示。

图5 无人机稳定下降速度随海拔高度的变化Fig.5 Variation of UAV steady descent velocity with altitude

图6 受1 m／s顺风、逆风、正侧风和逆侧风时，无人机前飞速度、侧向速度、下降速度和俯仰角的变化Fig.6 Variation of UAV forward velocity，side velocity，descent velocity and angle of pitch under 1 m／s down wind，against wind，positive side wind and reverse side wind

通过图6（c）以及无人机的开伞高度可以得到无人机由开伞到着陆的总时间，通过图6（a）、（b）以及得到的无人机由开伞到着陆的总时间可以得到无人机在伞降过程中沿其X轴方向的最大位移和沿Y轴方向的最大位移，具体结果如表1所示，选取无人机的开伞点作为参考点。

通过图6可以发现，无人机伞降回收系统速度会受到风力的影响而发生变化。以图6（c）为例，当无人机伞降回收系统受到1 m／s的顺风或逆风时，其稳定下降速度由 5.05 m／s变为了5.16 m／s和 5.22 m／s，分别增加了 2.2%和3.4%，这是由于顺风或者逆风导致无人机与降落伞间夹角减小，虽然无人机迎角增加将导致无人机升力增大，但是该过程中降落伞阻力减小占主导地位，因此其稳定下降速度会增大；受到1 m／s的正侧风或者逆侧风时其稳定下降速度均由5.05 m／s变为了4.23 m／s，均减小约16%，这是由于在侧风作用下无人机与降落伞共同发生侧向运动，虽然降落伞竖直方向阻力分量减少，但在该过程中无人机竖直方向升力增大占主导地位，因此其下降速度减小。图6（d）展示了无人机的俯仰角随降落伞充气完成后风力变化的曲线。通过曲线可以看出，无人机的俯仰角在降落伞充气完成后风力的作用下仍然处于稳定的状态，其中顺风与逆风对无人机的俯仰角影响较大，正侧风与逆侧风对无人机俯仰角影响较小，可以通过本文模型预估无人机在着陆时的姿态进而确保无人机在着陆时不会损伤机头、机翼等部件，保证无人机的安全回收。除此之外，无人机伞降受风力影响后的相关计算结果如表1所示。

通过表1可以发现，本文模型及仿真曲线能够得到无人机的稳定下降速度，进而确保无人机在着陆时不会由于冲击载荷太大而破坏无人机的结构；本文模型能够得到无人机在伞降过程中沿其X轴、Y轴的最大位移（相对于开伞点），这些量可以作为参考值进而为无人机回收着陆地点的选取提供指导。

表1 受1 m／s顺风、逆风、正侧风和逆侧风时，无人机伞降回收相关结果Table 1 Related results of UAV parachute recovery under 1 m／s down wind，against wind，positive side wind and reverse side wind

综上所述，风力和回收位置的海拔高度对于无人机的伞降回收有一定影响，应当实时考虑风力和海拔高度影响，进而保证无人机伞降回收的安全。

6 结 论

通过无人机伞降回收动力学建模与仿真实验可以得到如下结论：

1）本文采用多体动力学的方法引入了无人机对伞降回收系统的影响，基于平板扰流与迎角的耦合关系建立了无人机大迎角非定常气动力模型，采用凯恩方程参照无人机伞降回收系统质心建立了无人机伞降回收六自由度模型。通过仿真结果与实验数据的对比可以发现，该模型能够较好地反映无人机伞降回收过程中的运动趋势，能够为无人机的伞降回收提供指导。

2）分析了回收位置海拔高度和伞降回收过程当中风力对无人机伞降回收的影响，通过相应仿真曲线可以得出结论：风力与回收位置海拔高度对无人机伞降回收有一定影响，应当在伞降回收初期予以考虑，以保证伞降回收过程的稳定性和安全性。

3）在无人机伞降回收过程中，无人机的稳定下降速度、无人机所承受的最大冲击载荷、无人机的姿态变化以及无人机开伞后的最大水平位移均对无人机伞降回收的设计有重要意义。无人机所承受的冲击载荷在伞降过程中可分为两部分：即射伞过程中的冲击载荷以及无人机着陆时的冲击载荷，射伞过程中的冲击载荷主要由相应的射伞装置以及降落伞的大小决定，无人机着陆时的冲击载荷主要由无人机的稳定下降速度决定；通过无人机伞降回收过程中的姿态变化曲线能够推断无人机在着陆时的姿态，进而保证无人机在着陆时不会损伤无人机的机头、机翼等部位；开伞后无人机的最大水平位移将受到无人机所选用降落伞的大小、伞降过程中的风力、回收位置的海拔高度等因素的影响，开伞后无人机的最大水平位移对无人机伞降回收着陆地点的选取具有指导意义；本文所建模型能够仿真得到无人机伞降回收过程中的稳定下降速度、姿态变化、无人机回收位置海拔高度以及无人机伞降回收过程中风力变化对无人机伞降回收的影响，进而通过速度仿真曲线得到无人机伞降回收过程中的最大水平位移，通过稳定下降速度确定无人机伞降过程中的冲击载荷，在无人机的伞降回收设计初期可采用该模型进行大量的数值仿真以确定在不同的情况下、对于不同的无人机应如何选取降落伞才能达到设计要求，进而保证无人机的安全回收。

针对接下来的研究，其主要可以分为两部分：①继续完善无人机和降落伞动态过渡过程的动力学模型，以便能够更为准确地描述无人机伞降系统由降落伞完全拉直到降落伞和无人机一起稳定下降阶段的运动过程；②针对无人机伞降系统，可以将降落伞和无人机单独拿出来进行研究，建立无人机伞降回收系统更高自由度模型，以描述无人机与降落伞的相对运动，为无人机的伞降回收提供更为精确的指导。