动力翼伞系统拟坐标形式的多体动力学建模

张青斌， 高峰， 郭锐， 丰志伟， 葛健全

(1.国防科技大学 空天科学学院， 湖南 长沙 410073； 2.96796部队， 吉林 白山 135200)

0 引言

翼伞系统不仅具有折叠体积小、质量轻的优点，而且拥有优良的操稳和滑翔性能，可以在一定初始偏差和风场等干扰条件下实现定点无损着陆，广泛应用于高空投放、远距离渗透和精确着陆等任务中，具有极大的军事和民用价值。目前国外研发了一系列颇具代表性的空投系统，例如美国Atair Aerospace 公司研制的自主导航空投系统、美国Strong Enterprises公司研制的Screamer空投系统、美国Alabama大学和海军研究院联合研发的Snowflake空投系统，以及加拿大米斯特机动综合系统技术公司研制的Sherpa系统和德国Cassidian公司研制的Paralander空投系统等[1]。上述空投系统能完成从数千克到数十吨的货物精确空投，落点偏差甚至在几十米以内。翼伞空投系统已成为世界各军事强国空投技术的发展方向，但是由于理论研究难度大且技术复杂，目前世界上掌握精确空投核心技术并有产品正式列装的国家较少。

翼伞系统飞行动力学是开展翼伞系统设计研制的理论基础。文献[2-3]建立了多种自由度的翼伞多体动力学模型，分析了展弦比、后缘下拉方式、安装角及系统长度对雀降性能的影响。Goodrick[4]提出并建立了翼伞和回收物刚性连接的6自由度模型，研究了翼伞系统纵向运动的稳定性和飞行性能。Zhang等[5]在建立翼伞6自由度模型基础上，通过实验验证了6自由度模型的可行性，用该模型仿真得到了翼伞各项参数变化对翼伞飞行状态的影响。Slegers[6-7]建立了9自由度多体动力学模型，分析了安装角、伞衣弯曲度及后缘下偏量对翼伞侧向运动的影响。蒋万松等[8]采用Lagrange乘子法建立了两体8自由度动力学仿真模型，并通过空投试验验证了仿真模型的有效性。Meirovitch[9]利用分布系统的扩展哈密顿原理，导出混合系统的标准Lagrange方程，随后以拟坐标表示刚体运动，将混合系统的2阶Lagrange方程转化为一组适用于控制的状态方程。陈建平等[10]通过引入Lagrange乘子，获得了大型翼伞系统封闭的飞行动力学方程组，对某大型翼伞系统飞行动力学过程进行了数值仿真。杨海鹏等[11]假定伞绳和吊带分别通过铰链与翼伞和载荷相连，考虑了伞绳和吊带的弹性以及伞绳相交点与伞衣的相对运动，通过建模仿真得到了吊带、伞绳交点运动的情况。

但是在目前文献中，还很少有关于载荷连接模式和稳定性的研究成果。为此，本文借鉴航天器多柔体动力学建模，以某动力的弹载翼伞系统为工程背景(见图1)，考虑伞绳、载荷吊带的连接方式，针对翼伞稳定飞行的运动特性，将伞衣和载荷简化为刚体，伞绳和吊带交汇点处理为质点；分析翼伞系统各体的受力特性，给出了翼伞多项式形式的气动表达式；将各体平动速度、角速度作为18自由度拟坐标，采用Lagrange拟速度形式的建模原理[12]，建立了18自由度的动力学模型。利用空投试验证明了模型的有效性，仿真分析下拉后缘、突风扰动和不同吊挂模式下的相对运动特性，对比了6自由度单刚体和18自由度多体动力学模型的仿真结果。

图1 某小型翼弹载系统Fig.1 A small wing projectile-borne parafoil system

1 模型描述

如图2所示，翼伞系统4体模型主要包括伞衣、伞绳、吊带和载荷部分。本文将稳定飞行状态下的伞衣、载荷以及2个连接点作为4个刚体，伞衣和载荷分别记为刚体A和刚体B；伞绳和吊带的交汇点简化为质点，记为L和R；伞绳简化为4根绳索LA1、RA2、RA3和LA4，吊带简化为4根绳索LB1、RB2、RB3和LB4. 上述多刚体模型中刚体A和B分别有6个自由度，质点L和质点R分别有3个自由度，总共18个自由度。该模型可以较完整地描述伞衣与载荷的相对运动，反映了实际的空投工况和载荷吊挂模式。

图2 翼伞系统4体模型Fig.2 4-body model of parafoil system

图3 翼伞系统4体模型坐标系示意图Fig.3 Schematic diagram of 4-body model coordinate system of parafoil system

图4 侧面示意图Fig.4 Side view of parafoil system

2 翼伞系统受力分析

当翼伞系统在其稳定状态飞行时，主要受力包括重力、气动力和螺旋桨推力等。伞衣具有柔性大的变形特性，而且伞绳和载荷的连接吊带均为柔性绳索，很难获得其翼伞系统精确的气动特性。为此，本文不考虑翼伞下拉后缘过程中的迟滞效应，并忽略各部件的气动耦合作用，采用工程简化模型表示各体的气动力。

2.1 伞衣受力分析

关于伞衣气动模型的建立，刘媛等[13]通过有限体积元法求解Navier-Stokes湍流模型的控制方程，并对美国国家航空航天局兰利研究中心的三维翼伞进行了气动性能分析，验证了该方法的可行性；进而利用该方法分析得到了不同迎角来流条件以及伞衣后缘不同程度下偏的气动力数据。汪龙芳等[14]对伞衣透气性进行了二维和三维数值模拟，得到了透气性和迎角对翼伞升阻力系数变化的影响。为了便于后续归航控制设计，本文在计算流体力学基础上，采用多项式表示伞衣质心处的气动力和气动力矩[1]，即

(1)

(2)

式中：ρ为大气密度；v为伞衣质心相对于气流的速度大小；SP为伞衣面积；δS和δa为对称舵偏和反对称舵偏；p、q和r分别为伞衣滚转角速度、俯仰角速度和偏航角速度；b为伞衣展长；c为伞衣弦长；其他参数k、CD0、CL、CDδS、CDδa、CL0、CLα、CLδS、CLδa、CYβ、CYp、CYr、CYδa、Clβ、Clp、Clr、Cl、Clδa、Cm0、Cmα、Cmq、CmδS、Cmδa、Cnβ、Cnp、Cnr、Cnδa均为气动力常系数，由试验得到。

由于伞衣充气张满后，其平均密度与空气密度相当，必须考虑附加质量的影响。本文基于文献[15-16]，将伞衣附加质量等价地处理为伞衣质心的作用力和力矩，并在伞衣坐标系下表示为

FPM=-ωP×(MavP)，

(3)

MPM=-ωP×(JaωP)，

(4)

式中：ωP和vP为伞衣角速度和质心相对气流速度；Ma和Ja分别为附加质量相对伞衣质心的质量矩阵和惯量矩阵，即

(5)

(6)

ma,ii和Ia,ii(i=1，2，3)分别为附加质量相对伞衣质心的平动和转动分量。

作用于刚体A质心的合力，可在刚体A体坐标系下表示为

(7)

作用于刚体A的合力矩可在刚体A体坐标系下表示为

(8)

式中：MrtA和MrqA分别为各伞绳张力和气动力对刚体A作用的力矩矢量和。

2.2 载荷受力分析

考虑到坐标变换，作用于刚体B质心处的合外力在其体坐标系下表示为

(9)

式中：Ft为作用于载荷的螺旋桨推力；GB为刚体B所受重力；FrtB和FrqB分别为与载荷相连的各吊带张力总和与气动力矢量和；FqB为载荷所受气动力，忽略螺旋桨对载荷气动力的影响，则载荷所受气动力为

(10)

式中：vB和Sd分别为载荷质心速度和阻力面积。

作用于刚体B的合力矩在体坐标系下表示为

(11)

式中：MrtB和MrqB分别为与载荷相连的各吊带张力和气动力对载荷的作用力矩矢量和，在刚体B体坐标系下分别表示为

(12)

(13)

式中：rOBi为载荷体坐标系下吊带端点Bi的位置矢量，i=1，2，3，4；FrtLB1、FrtLB4、FrtRB2、FrtRB3分别为吊带LB1、LB4、RB2、RB3的张力，方向从载荷接触点指向吊带交汇点；FrqLB1、FrqLB4、FrqRB2、FrqRB3分别为吊带LB1、LB4、RB2、RB3受到的气动力在B1、B4、B2、B3上的集中力。

2.3 伞绳和吊带交汇点的受力分析

作用于连接点L的合力可表示为

FL=GL+FrtL+FrqL,

(14)

式中：GL为质点L所受重力；FrtL和FrqL分别为与质点L相连的各伞绳(吊带)张力和气动力。类似地，作用于连接点R的合力在惯性坐标系下可表示为

FR=GR+FrtR+FrqR,

(15)

式中：GR为质点R所受重力；FrtR和FrqR分别为与质点R相连的各伞绳(吊带)张力和气动力。另外，在具体仿真计算中将作用在绳索上的气动力等效为在两端点的集中力。

3 基于拟坐标形式的翼伞动力学方程

(16)

本文采用拟坐标形式Lagrange方程进行动力学建模，经过复杂推导可得到如(17)式的动力学方程：

(17)

式中：相关状态变换矩阵

(18)

(19)

广义质量矩阵

(20)

I3×3为3阶单位矩阵，03×3为3阶0矩阵；广义力的相关表达式为

(21)

(22)

(23)

(24)

(25)

(26)

(27)

4 仿真分析

某翼伞系统的伞衣展长b=3.404 m，弦长c=0.920 m，伞衣和载荷安装角均为0°，伞衣质量mP为0.5 kg、载荷质量mS为7.0 kg，初始高度和系统速度分别为1 000.0 m、v0=[7.5 0 3.8]Tm/s，其他状态量的初始值均为0. 下面分别基于单刚体和多刚体模型进行飞行特性分析。需要说明的是，本文基于文献[3]中的研究成果，不考虑载荷的姿态变化和吊挂模式，将R点和L点的载荷等效为质点，建立了翼伞系统6自由度单刚体动力学模型。

1)单刚体和多刚体动力学模型的对比分析。图5和图6所示分别为无控状态下(无推力、无后缘下拉偏量)翼伞姿态和伞绳连接点位移的变化曲线。由图5和图6可以看出，因为初始状态值为非稳定状态，所以仿真初始时刻出现一定周期的震荡才能达到稳定飞行状态。仿真结果表明，尽管两种模型仿真的轨迹曲线几乎一致、衰减趋势基本一致，但是伞衣姿态角度稍有差别。在单刚体动力学模型中，由于将载荷简化为质点、忽略了载荷在物- 伞运动过程中的相对运动，其稳定过程中主要依靠伞衣来消耗能量，因此伞衣摆动幅度相对结果略大。

图5 无控时伞衣质心的位移变化曲线Fig.5 Changing curves of center of mass of canopy without control

图6 无控时伞衣姿态变化曲线Fig.6 Changing curves of attitude of canopy without control

2)转弯特性分析。图7和图8所示为基于18自由度多体模型、单侧下偏量100%的载荷质心轨迹和角速度变化曲线。从图7和图8中可以看出：在恒定下偏量作用下翼伞系统具有稳定的转弯半径(本例最小转弯半径约52.0 m)；同时在翼伞系统稳定转弯下滑过程中，载荷也具有稳定的角速度。

图7 单侧下拉时载荷质心轨迹Fig.7 Centroid locus of payload when pulling down on the side

图8 转弯时载荷角速度变化曲线Fig.8 Changing curves of angular velocity of payload during turning

3) 带动力飞行时的运动特性分析。图9和图10所示为30 s时在载荷上施加50 N螺旋桨推力情况下载荷的轨迹和姿态角变化曲线。从图9和图10中可以看出，在一定推力作用下，翼伞系统可以实现平飞，但是伞衣和载荷都具有较大的俯仰角度。

图9 平飞时的载荷质心位移Fig.9 Centroid displacement of payload during level flight

图10 平飞时的伞衣和载荷俯仰角Fig.10 Attitude angles of canopy and payload during level flight

4)突风作用的影响分析。图11所示为在突风2 m/s情况下，设计状态(工况1)、伞绳长度增加20%(工况2)、吊带增加30%(工况3)3种工况下载荷的角速度变化曲线。由图12可以看出，在突风影响下，适当增加吊带长度有利于载荷姿态的稳定。

图11 突风作用下的载荷姿态变化曲线Fig.11 Changing curves of payload attitude under wind gust

5)空投试验分析。采用直升机进行了空投试验(天气良好)，测得了某次轨迹曲线(单侧下拉后缘50%)。图12～图14为载荷位置随时间的变化曲线。由图12～图14可见，理论仿真数据与试验数据基本一致，证明了模型的有效性。二者的差异性主要是受高空风场的作用。

图12 载荷xB轴方向位移随时间变化曲线Fig.12 Payload displacement versus time in xB direction

图13 载荷yB轴方向位移随时间变化曲线Fig.13 Payload displacement versus time in yB direction

图14 载荷zB轴方向位移随时间变化曲线Fig.14 Payload displacement versus time in zB direction

5 结论

本文针对翼伞稳定飞行时的运动状态，综合分析各体受力特点和运动特性，利用拟坐标形式的Lagrange原理建立了反映伞绳构型和载荷吊带连接方式的4体18自由度动力学模型。以某小型动力翼伞系统为例，对比研究单刚体和多体动力学模型的仿真特点，分析了多体模型下单侧下拉、突风作用和各种吊挂模式下的伞衣和载荷运动特性，以及外加推力后的伞衣姿态角和载荷位移变化情况，并与试验飞行数据做比较。主要得出如下结论：

1)相比于6自由度模型，18自由度模型能较真实地反映翼伞系统内部的相对运动，更接近真实情况。

2)在转弯过程中，伞衣和载荷均有稳定的角速度。

3)载荷施加推力作用，伞衣将具有较大的俯仰角。

4)适当增加吊带的长度，可增大系统受突风影响的鲁棒性。

5)建立的模型仿真结果与试验数据比较，拟合度较好，证明了该模型的准确性。