伞翼无人机精确建模与控制

朱虹，孙青林，邬婉楠，孙明玮，陈增强

南开大学 人工智能学院，天津 300350

近年来，伞翼无人机凭借其载荷比高，成本低廉，可低空低速飞行以及可实现雀降无损着陆等优势在航空航天等领域得到了广泛应用［1－3］。翼伞的伞衣由柔性织物构成，展开过程中，空气从伞衣的前缘切口进入，在伞衣内滞留从而维持伞衣的气动外型［4－5］。充气完成后，伞衣类似一个带有前缘切口的小展弦比翼型［6］。伞翼无人机的操纵方式主要为通过伞衣后缘的单侧下偏实现转弯，通过操纵绳双侧下偏实现着陆。伞翼无人机的精确建模对归航与控制起着至关重要的作用。在建模的过程中，关键是伞衣气动力的计算。气动力的计算有两个难点，首先是伞衣的柔性导致伞衣与流场之间的强非线性的耦合很难用机理建模的方法进行分析，其次，前缘切口以及后缘下偏对气动力的影响无法直接从理论的角度计算。以上两个因素对翼伞的精确建模以及控制提出了挑战。

现有基于CFD对翼伞的研究以定性分析为主，针对前缘切口和后缘下偏对气动力的影响进行定量分析的研究较少。朱旭和曹义华［4］研究了翼伞弧面下反角、前缘切口和翼型对翼伞气动性能的影响。黄炎等［7］研究了不同前缘切口参数下的大型冲压翼伞的气动性能随迎角的变化趋势。谢亚荣［8］、Muller［9］和胡文治［10］等分别建立了六自由度、八自由度和九自由度翼伞模型，其中翼伞气动力的计算都采用传统的基于升力线理论的方法。梁炜等［11－12］综合风雨环境与襟翼偏转的影响，采用最小二乘法进行模型参数辨识，实现了翼伞气动方程的修正；陶金等［13］结合CFD技术和Kirchhoff运动方程，提出了一种更加精确的翼伞动力学模型。两者对翼伞的气动模型做了简化，均未考虑前缘切口和后缘下偏对翼伞动力学模型精度的影响。

本文基于上述研究成果，针对后缘下偏和前缘切口对翼伞气动性能的影响展开研究，并将翼伞的肋片考虑在内［14］。通过求解 Navier-Stokes方程对气动模型的参数进行辨识，建立了更加精确的伞翼无人机的六自由度动力学模型，并通过空投试验数据对本文改进的动力学模型的有效性进行了验证，有效提高了翼伞系统的控制精度。

1 伞翼无人机动力学模型

为简化分析，建模之前首先引入3个必要的坐标系描述伞体运动。

1）大地坐标系OdXdYdZd：主要用于描述伞体的姿态、位置等信息，Od为空间中某一固定点，OdZd轴铅垂向下，ΟdXdYd与水平面平行，其中OdXd指向翼伞初始运动方向。

2）伞体坐标系OsXsYsZs：主要用于伞衣、伞绳和载荷（以下简称伞体）建立动力学方程，坐标原点位于翼伞系统的质心，OsZs轴位于伞体纵向平面内指向负载，ΟsXs同样位于伞体纵向平面内指向伞衣前缘，ΟsYs轴与其他两轴构成右手坐标系。

3）气流坐标系OqXqYqZq：用于气动力计算，原点位于翼伞的压心，ΟqZq轴位于伞体纵向平面内指向下方，ΟqXq为相对气流的方向，ΟqYq轴与其他两轴构成右手坐标系。

伞体坐标系到大地坐标系的转换矩阵及以上坐标系的具体定义参考文献［15］。由于本文的重点是气动力的计算及动力学建模，只考虑伞体沿自身坐标系3个轴的平动和绕轴的转动6个自由度，忽略伞衣和载荷之间的相对运动，建立如下伞翼无人机的六自由度动力学方程：

式中：Ar、Aa分别为真实质量和附加质量的惯性矩阵；FG、MG分别为重力及其力矩；Vo、W 分别为伞体质心的速度和角速度的一阶微分形式；Fr、Fa、Mr、Ma由 动 量 和 动 量 矩 展 开 式 得 到；Faero、Maero分别为伞翼无人机的气动力和气动力矩。

传统的翼伞模型在进行气动力计算时通常将伞衣沿展向均匀分成8片，不考虑前缘切口和后缘下偏的影响，限制了模型的精度。计算气动力时每片伞衣的升力系数从外到内分别乘以修正因子0.6、1.0、1.16、1.24，根据各自对应的速度和迎角，计算每片的气动力和力矩，对其求和即可得到伞衣整体的气动力和力矩：

式中：Ti－Os为第i片伞衣自身的体坐标系相对于伞体坐标系的转换矩阵；LOs－i为每片伞衣的质心在伞体坐标系下的位置矢量；FLi和FDi分别为分片i的升力与阻力，一般形式为

其中：i（i＝1，2，…，8）表示每片伞衣的序号；k为乘积因子；S为翼伞的特征面积；ρ为空气密度；u、v、w分别为速度在以压力中心为原点的坐标系中的分量；CL和CD分别为升力、阻力系数，传统的基于升力线理论的计算公式为

式中：CLα为升力线斜率；α为迎角；α0为零升力迎角；φ为弧面下反角；CD0为翼伞的型阻；式（7）右端第2项为诱导阻力；δ为翼型影响因子；AR为伞衣的展弦比；j为关于展弦比的函数，经验值为

升力、阻力系数CL和CD的计算是伞翼无人机动力学模型中气动力计算的重点和难点，以上传统方法中存在较多的人为因素，如修正因子的设置忽略了伞衣的形变，这对于模型的精度有很大影响。要想提高伞翼无人机动力学模型的精度，主要从影响气动力的因素着手，修正和完善气动方程。

从以上的分析中可以看出，虽然升力线理论可以用来预测伞衣的气动特性，但却忽略了前缘切口对翼伞气动性能的影响，并且未给出关于后缘下偏影响的具体表现形式。由文献［16］可知，前缘切口对翼伞的气动性能有重要影响，仅对前缘切口和后缘下偏的影响做定性分析是远远不够的，定量分析才是伞翼无人机精确建模的关键。

2 气动力计算

2.1 气动模型和流场设置

基于CFD对前缘切口和后缘下偏的影响进行定量分析，从而更加准确地计算气动力参数，实现伞翼无人机的精确建模。在保证计算精度的前提下，做出如下假设：翼伞近似为刚体，主要研究伞衣稳态时的特性；忽略伞绳及负载对气动力的影响；伞衣材料不透气。

根据空投试验中所用翼伞的物理模型给出CFD模型，如图1所示。该翼伞模型由15个气室组成，各气室相连的肋片上开有气孔以保证各气室压力平衡。翼伞模型的详细参数见表1，其中c表示伞衣的弦长。流场边界条件应尽量远离扰动源，但实际计算中流场尺寸不宜过大，否则网格数量过多会造成误差积累，根据计算经验将流场边界与翼伞的距离设置为10c、5c、10c，如图2（a）所示。

翼伞工作的雷诺数为1.4×106，因此伞衣周边呈湍流形态，需要启用湍流模型，这里使用Spalart-Allmaras（SA）湍流模型来计算气动参数。

流场和伞衣的网格如图2所示。所有流场的计算都是在软件Fluent17.2中完成的。将流场左侧边界和底面作为速度入口，前侧边界和后侧边界设置为非滑移壁面，其余两个面为压力出口，流体为不可压空气，来流速度为10m／s。控制方程采用二阶精度，利用SIMPLIC算法求解速度和压力之间的耦合。

图1 空投试验中所用翼伞的物理模型Fig.1 Physical model for parafoil used in airdrop experiments

表1 翼伞模型参数Table 1 Parameters of parafoil models

图2 流场和伞衣模型网格划分Fig.2 Mesh generation for flow field and parafoil model

2.2 网格无关性验证

基于CFD的数值模拟结果与网格密度紧密相关，只有当计算域的网格数量足够多时结果才可靠，为了保证CFD模拟的准确性，首先进行网格无关性验证。表2给出了模型1不同网格密度的详细数据，通过调整网格尺寸来改变网格密度。

表2 翼伞模型1的详细网格参数Table 2 Detailed grid parameters of parafoil model 1

图3 网格无关性验证Fig.3 Verification of mesh independence

图3 为4个模型使用SA湍流模型数值模拟得到的升力系数CL和阻力系数CD曲线，由图中可以看出，中等网格和小网格的计算结果基本吻合，而大网格的计算结果只在部分范围内匹配，且前缘切口越大，仿真结果对网格密度越敏感。由于采用中等网格和小网格的计算结果基本一致，因此为了节省时间，以下仿真中全部采用中等网格。

2.3 CFD仿真和风洞试验对比

图4 风洞试验［18］和CFD仿真在无后缘下偏时的结果对比Fig.4 Comparison of wind tunnel test［18］and CFD simulation without flap deflection

图5 风洞试验［18］和CFD仿真在后缘下偏时的结果对比Fig.5 Comparison of wind tunnel test［18］and CFD simulation with flap deflection

将数值模拟结果与风洞试验数据做对比来验证数值模拟方法的有效性［17］。由于空投试验所采用的模型（模型1）没有风洞试验数据，故使用文献［18］中风洞试验的模型1（并非本文模型）作为验证模型，结果如图4和图5所示。在小迎角的情况下CFD的计算结果与风洞试验数据基本一致，随着迎角的增加，仿真结果与风洞试验数据相比逐渐偏大，这是由于忽略了伞衣气室的“鼓包”现象以及伞衣的柔性。在数值模拟的过程中这些偏差的存在是合理的，从而验证了CFD数值模拟的有效性。

2.4 气动力参数辨识

考虑前缘切口和后缘下偏，将伞衣的气动方程修正为

式中：γ、h分别为前缘切口的角度和长度；δe为下偏量。

式（8）中前缘切口系数CL，γ、CL，h与迎角之间是非线性关系，因此所有与迎角有关的影响因子使用式（9）所示的二阶模型即可满足，其中X∈｛L，D｝，采用最小二乘法对以上参数进行辨识［19］。

辨识所需的气动力数据为离线的数据［20］，根据最小二乘法得到待辨识参数的形式为

式中：θ为待辨识的参数矩阵；YN、ΦN分别为气动力系数矩阵和与迎角、切口因子、偏转因子相关的输入矩阵。

由于偏转系数和前缘切口系数属于静态参数，因此可以利用最小二乘法对参数进行离线辨识，得到的详细参数值如表3所示。

表3 与气动力有关的参数辨识结果Table 3 Identification parameters related to aerodynamics

3 伞翼无人机自主归航轨迹跟踪控制

将以上气动力参数的辨识部分加入伞翼无人机的动力学模型中，考虑到翼伞在飞行过程中不可预知的干扰，采用线性自抗扰控制器（Linear Active Disturbance Rejection Controller，LADRC）对伞翼无人机的位置进行跟踪控制［21］，实现沿目标轨迹的自主归航运动。控制器结构如图6所示。图中：y为系统输出；vo为给定值；sat（ ）为限幅函数；kp、kd分别为控制器比例和微分系数；w为系统的外部扰动；b为系统参数。反馈控制率为

观测器ESO的形式为

式中：z＝［z1z2z3］T分别为系统输出及其微分y、y和系统总扰动f的观测值；L为ESO的增益向量；系数矩阵分别为

式（11）～式（14）共同构成了 LADRC，对伞翼无人机进行自主归航轨迹跟踪控制。

图6 LADRC结构示意图Fig.6 Structural diagram of LADRC

4 仿真分析

4.1 前缘切口的影响

图7 CFD仿真得到的升力、阻力系数Fig.7 Lift and drag coefficients based on CFD simulation

Fig.8 Lift-to-drag ratio curves based on CFD simulation

图7 和图8分别为模型1～模型4在无后缘下偏条件下采用CFD数值模拟方法所得的升力、阻力系数和升阻比K的变化曲线，分析发现：前缘切口导致升力系数降低，阻力系数增大，且迎角越大，对翼伞气动性能的影响也越大；前缘切口对失速迎角的影响较小，且切口长度相同的条件下升力曲线斜率基本一致；随着切口长度的增加和角度的减小，在相同迎角下伞衣的升力系数越来越小，阻力系数越来越大，切口长度比角度的影响更大；前缘切口的尺寸越小，升阻比变化越缓慢，越有利于伞翼无人机的稳定性。

为进一步分析前缘切口对翼伞气动性能的影响，图9和图10给出了模型1～模型4在迎角α＝10°时翼伞上、下表面的压力分布，由图可知，前缘切口会影响伞衣表面的压力分布，切口长度越大，伞衣上表面的低压区越靠近前缘，且面积逐渐减小，同时下表面的高压区域减小。此外，图11给出了模型1～模型4翼伞中部弦向剖面的速度分布云图：切口越长，流经前缘的气流速度损失越严重，同时切口角度的减小会导致后缘处的气流分离，降低翼伞的气动性能。

Fig.9 Pressure distribution on parafoil upper surface

Fig.10 Pressure distribution on parafoil bottom surface

图11 伞衣中部弦向剖面的速度分布云图Fig.11 Velocitydistribution contour on central crosssection of parafoil

4.2 后缘下偏的影响

文献［18］中NASA风洞试验和空投试验表明，在制动偏转期间，只有伞衣后缘1／4的翼尖向下弯折，整个偏转角为75°，以模型1为例，下偏量为1时的气动模型如图12所示。图13为模型1～模型4在下偏量δ分别为1／3（25°）、2／3（50°）和1（75°）时的气动力参数变化情况。无论前缘切口的大小如何，后缘下偏都可以使升力、阻力系数增加。小迎角条件下升力系数随着δ的增加而增加，大迎角时，升力系数反而减小，阻力系数始终呈上升趋势，而且下偏量相同时，不同模型的气动特性不同。

模型1～模型3在下偏量δ＝2／3时的压力分布云图如图14所示，在单侧下偏的情况下，伞衣表面的压力呈非对称分布［22］。由图14可以看出，随着切口长度的增加，伞衣下表面的高压区域逐渐减小，上表面的低压区域逐渐增大，这是因为切口越长，气流经历的减速过程就越长，因此上表面的低压区域随之增大，而且随着偏转量的增加，该趋势越来越明显。

图12 模型1后缘下偏量δ＝1的CFD模型Fig.12 CFD model with trailing edge deflectionδ＝1of model 1

图13 4组模型不同后缘下偏量对应的升力、阻力系数Fig.13 Lift and drag coefficients of four models under different trailing edge deflections

图14 δ＝2／3时的伞衣上下表面压力分布云图Fig.14 Pressure distribution contour on upper and bottom surfaces of parafoil whenδ＝2／3

4.3 基于LADRC轨迹跟踪控制的仿真结果

基于以上气动力参数的辨识对伞翼无人机动力学模型的改进以及LADRC控制器，对伞翼无人机在水平面内进行沿圆心在［1 500，1 500］m、半径为800m的圆形轨迹的运动仿真，实现沿目标轨迹的自主归航运动［23］。仿真环境设置为：伞翼无人机的初始位置为［500，3 000，2 000］m，控制器的参数设置为kp＝0.42，kd＝2.9，b＝0.02。仿真时间为600s，仿真结果如图15～图17所示。

从仿真结果可以看出，LADRC能够很好地控制伞翼无人机沿目标圆的运动，从初始位置到达目标圆轨迹需要约100s，随后的500s沿目标 圆运动的误差几乎恒为零，下偏控制量维持在合理范围内，验证了伞翼无人机改进后的动力学模型的有效性。

图15 水平面内的自主归航轨迹跟踪曲线Fig.15 Autonomous homing trajectory curves in horizontal plane

图16 归航轨迹跟踪误差Fig.16 Tracking error of homing trajectory

图17 控制量变化曲线Fig.17 Control quantity variation curve

5 空投试验验证

本课题组利用模型1于2018年在河南省许昌市进行了一次空投试验，用于验证以上所建立的伞翼无人机动力学模型的准确性，翼伞模型的基本参数见表1。伞翼无人机被带到距海平面400m的高空，从热气球上释放，如图18所示。进行空投试验的地面海拔是90m，释放位置的海拔是403m，降落伞打开瞬间的下落高度是40m，左侧操纵线的偏转量δ＝1／3。半圆轨迹的仿真和空投试验的结果对比如图19所示。

图18 空投试验Fig.18 Airdrop experiment

图19 空投试验和仿真结果对比Fig.19 Comparison of airdrop experiment and simulation results

图19 （b）和图19（c）中，模型的垂直速度和水平速度与空投试验的结果基本吻合，误差在允许范围内。由于大气的水平运动比垂直运动对翼伞的影响更大［6］，所以水平方向的误差比垂直方向大。试验当天，风速计测得的风速比较小，后缘下偏量设置为常值，因此垂直方向速度变化较小。误差一方面由于未考虑翼伞运动过程中气室“鼓包现象”对气动力的影响，另一方面来自于飞行环境中的复杂噪声。通过以上分析，基于CFD的气动力参数计算能够精确描述实际伞翼无人机的特征。

6 结 论

1）研究了前缘切口和后缘下偏对翼伞气动力的影响：切口长度和角度的增加会使翼伞的气动性能变差，对升阻比有很大影响；在单侧后缘下偏的情况下，伞衣表面的压力分布不对称。综合考虑前缘切口和后缘下偏可以提高伞翼无人机动力学模型的精度。

2）结合前缘切口和后缘下偏对伞衣气动力的影响，修正了气动力的计算。将升力、阻力系数的辨识加入伞翼无人机的动力学模型中，提高了动力学模型的精度。

3）基于改进的伞翼无人机的动力学模型采用LADRC控制，进行了自主归航轨迹跟踪仿真，仿真结果和与空投试验的对比验证了气动参数辨识方法和改进的动力学模型的有效性，对于伞翼无人机控制器的设计和仿真具有重要意义。