符号说明
翼伞具有良好的滑翔性能、操纵性和稳定性，能够实现“雀降”和无损着陆，克服了传统降落伞减速系统飞行轨迹“随风飘”、落点散布大的缺点，在航空航天飞行器的精确定点回收、军民用物资精确空投和人员装备空降等方面有着很高的应用价值和广阔的发展前景。载荷平台摄像机拍摄的相对偏航运动如图 1所示。
在翼伞系统的控制系统和控制算法设计时，首先要对翼伞系统的运动特性进行深入研究。但翼伞系统与传统航空器不同，翼伞与有效载荷之间一般存在相对运动，因此要从多体系统动力学的角度建立翼伞-载荷系统动力学模型。此外，翼伞是一种超轻结构，在分析翼伞系统的运动特性时，需要考虑翼伞的附加质量影响[1-2]。
对于翼伞系统的多体动力学研究，最早出现在1968年Wolf翼伞系统稳定性研究论文[3]。目前翼伞系统的多体动力学分析模型归纳起来主要分为3类：
(1) “两体+弹簧”，即将翼伞和有效载荷分别表示为两个独立的刚体，两体之间通过弹簧模型相连，弹簧参数设置依赖于设计和经验，而两体上的连接点的空间位置显式求解，这增加了求解自由度，以Vishnyak[4]、国防科技大学熊菁[5]等为代表。
(2) “两体+约束方程”，即将伞体和有效载荷分别表示为两个独立的刚体，其连接点通过约束方程建立关系，这保证了连接点空间位置的一致性，且建模可以程式化，但增加了求解的自由度，以Wolf[3]、Pillasch[6]、 Wise[7]、 Stricker[8]和Christiaan[9]等为代表。
(3) 采用“共铰点两体+扭簧”模型，是最小解集多体系统，不增加求解自由度，但扭簧参数设置依赖于实际系统设计和仿真经验，以Barrows[10]、 Prakash[11]、 Muller[12]、 Mooij[13]、 Gorman[14]、 Slegers[14]和Costello[15]等为代表。
本文采用拉格朗日乘子法对一般可控翼伞-载荷系统建立了8自由度多体动力学仿真模型，采用Barrows表观质量估算方法[1-2]进行翼伞质量计算，对3种工况进行了仿真分析，并与相应的试验结果进行了对比分析，初步验证了动力学仿真模型的有效性。
1 翼伞系统动力学模型
1.1 研究对象描述
本文主要开展翼伞系统滑翔飞行阶段直至雀降过程的飞行动力学行为研究。研究对象如图 2所示，由翼伞(主要包括伞衣、伞绳和吊带)和有效载荷两大部分组成，控制系统设备通常布置在有效载荷上。根据该阶段的运动特性，动力学建模做如下基本假设：
(1) 大地是水平的，忽略地球自转，且重力加速度恒定。
(2) 翼伞完成充气张满后几何形状不变，左右对称，后缘操纵只影响气动力。
(3) 翼伞和有效载荷均看作6自由度刚体，两体间通过刚性吊带连接约束。
(4) 有效载荷的质量特性始终不变，但翼伞的附加质量特性随大气密度变化。
(5) 不考虑有效载荷的气动力。
(6) 忽略翼伞操纵过程的响应延迟。
1.2 动力学方程 本文动力学模型涉及到惯性坐标系、各物体本体坐标系及气流坐标系，坐标系之间的关系可用z-y-x转序的转换矩阵描述，可参见飞行力学或多体系统动力学教科书。一般翼伞系统动力学模型见图 3。如图 3所示，翼伞和载荷可分别看作6自由度刚体，两体间通过约束建立耦合运动关系，翼伞系统的拉格朗日动力学方程[16]为
$\frac{d}{dt}\frac{\partial T}{\partial \dot{q}}-\frac{\partial T}{\partial q}+{{\left( \frac{\partial f}{\partial q} \right)}^{T}}\lambda =Q$
(1)
$q={{\left( {{x}_{p}},{{\theta }_{p}},{{x}_{C}},{{\theta }_{C}} \right)}^{T}}$
(2)
式中q为广义坐标，为12×1列阵，不能直接积分。从本体系到惯性系的坐标转换矩阵为S，则有
$\begin{align}
& {{{\dot{x}}}_{p}}={{S}_{p}}{{v}_{p}} \\
& {{{\dot{x}}}_{C}}={{S}_{C}}{{v}_{C}} \\
\end{align}$
(3)
物体(翼伞或载荷)的动能T可以表示为
$2T={{{\dot{x}}}^{T}}SM{{S}^{T}}\dot{x}+{{\omega }^{T}}I\omega +{{{\dot{x}}}^{T}}SH\dot{\omega }$
(4)
式中M为物体的真实物理质量矩阵，对于翼伞还应再加上附加质量项[5]，后文说明其计算方法。
根据式(1)，可以推导得到
$\frac{d}{dt}\left( \frac{\partial T}{\partial \dot{q}} \right)-\frac{\partial T}{\partial q}=\left. \left| \begin{align}
& {{S}_{p}}\left( {{M}_{p}}{{{\dot{V}}}_{p}}+{{\omega }_{p}}\times \left( {{M}_{p}}{{V}_{p}}+{{H}_{p}}{{\omega }_{p}} \right)+{{H}_{p}}{{\omega }_{p}} \right) \\
& H_{p}^{T}{{{\dot{V}}}_{p}}+{{V}_{p}}\times ({{M}_{p}}{{V}_{p}}+{{H}_{p}}{{\omega }_{p}})+{{I}_{p}}{{{\dot{\omega }}}_{p}}+{{\omega }_{p}}\times ({{I}_{p}}{{\omega }_{p}}+H_{p}^{T}{{V}_{P}}) \\
& {{S}_{C}}({{M}_{C}}{{{\dot{V}}}_{C}}+{{\omega }_{C}}\times {{M}_{C}}{{V}_{C}}+{{H}_{C}}{{{\dot{\omega }}}_{C}}) \\
& H_{C}^{T}{{{\dot{V}}}_{C}}+{{V}_{C}}\times {{H}_{C}}{{\omega }_{C}}+{{I}_{C}}{{{\dot{\omega }}}_{C}}+{{\omega }_{C}}\times H_{C}^{T}{{V}_{C}} \\
\end{align} \right. \right\}$
(5)
通过如下3个约束模型假设，建立翼伞与载荷间的4个约束方程：
(1) 翼伞-平台连接特征线A0B0与载荷平台B1B2正交，即约束f1，表达式为
${{S}_{C}}{{e}_{2}}\cdot {{S}_{p}}\varepsilon =0$
(6)
式中ε=(-sinθr,0,-cosθr)T。
(2) 翼伞和载荷连接特征线A0B0的A0点在翼伞纵向对称面内，且相对翼伞坐标系z方向允许相对位移，而x和y方向限制其位移，从而f2和f3分别为
${{\left( {{x}_{C}}+{{S}_{C}}{{r}_{B0}}-{{x}_{p}} \right)}^{T}}{{S}_{p}}{{e}_{1}}=0$
(7)
${{({{x}_{C}}+{{S}_{C}}{{r}_{B0}}-{{x}_{p}})}^{T}}{{S}_{p}}{{e}_{2}}=0$
(8)
式中rB0为B0点在载荷坐标系下的位置向量。
(3) 翼伞单个吊带长度假设不发生变化，即f4表示为
$\begin{align}
& {{({{x}_{B1}}+{{S}_{C}}{{r}_{B1}}-{{x}_{A1}}-{{S}_{p}}{{r}_{A1}})}^{T}}({{x}_{B1}}+ \\
& {{S}_{C}}{{r}_{B1}}-{{x}_{A1}}-{{S}_{p}}{{r}_{A1}})-{{l}^{2}}=0 \\
\end{align}$
(9)
式中rA1和rB1分别为单根吊带两端点在载荷和翼伞本体坐标系下的位置向量。
通过以上约束，翼伞与载荷间允许相对俯仰和偏航，但偏航运动耦合着相对重心位置运动和相对滚动运动。
4个约束方程的一阶和二阶导数可以统一表示为
$\left\{ \begin{align}
& \dot{f}={{f}_{q}}\dot{q}=B\dot{q}=0 \\
& \ddot{f}={{f}_{q}}\ddot{q}+{{{\dot{f}}}_{q}}\dot{q}=B\ddot{q}+R=0 \\
\end{align} \right.$
(10)
式中：为本体系下表达的拉格朗日方程广义坐标；${\dot{q}}$=(Vp,ωp,VC,ωC)T为本体系下表达的广义速度；f=(f1,f2,f3,f4)T为约束方程；B为雅可比矩阵；R为余项列向量，分别表示为
$B=\left[ \begin{matrix}
0 & {{(S_{\Delta }^{T}{{e}_{2~}}\times \varepsilon )}^{T}} & 0 & -{{({{S}_{\Delta }}\varepsilon \times {{e}_{2~}})}^{T}} \\
-e_{1}^{T} & -{{(S_{p}^{T}({{x}_{C~}}+\text{ }{{S}_{C~}}{{r}_{B0~}}-{{x}_{p~}})\times {{e}_{1~}})}^{T}} & {{({{S}_{\Delta }}{{e}_{1~}})}^{T}} & {{({{S}_{\Delta }}{{e}_{1~}}\times {{r}_{B0~}})}^{T}} \\
-e_{2}^{T} & -(S_{p}^{T}{{x}_{C~}}+\text{ }{{S}_{C~}}{{r}_{B0~}}-{{x}_{p~}})\times {{e}_{2~}}{{)}^{T}} & {{({{S}_{\Delta }}{{e}_{2~}})}^{T}} & {{({{S}_{\Delta }}{{e}_{2~}}\times {{r}_{B0~}})}^{T}} \\
-2{{\rho }^{T}}{{S}_{p~}} & -2{{\rho }^{T}}{{S}_{p~}}({{\omega }_{p~}}\times {{r}_{p~}}) & 2{{\rho }^{T}}{{S}_{C~}} & 2{{\rho }^{T}}{{S}_{C~}}({{\omega }_{C~}}\times {{r}_{C~}}) \\
\end{matrix} \right]$
(11)
$R=\left[ \begin{matrix}
-{{(\dot{S}_{\Delta }^{T}{{e}_{2}}\times \varepsilon )}^{T}}{{\omega }_{p}}-{{(\dot{S}_{\Delta }^{T}\varepsilon \times {{e}_{2}})}^{{}}}T{{\omega }_{C}} \\
{{R}_{2}} \\
{{R}_{3}} \\
2{{{\dot{\rho }}}^{T}}\dot{\rho }+2\rho _{2}^{T}{{S}_{p}}\left( {{\omega }_{p}}\times ({{v}_{p}}+{{\omega }_{p}}\times {{r}_{A0}}) \right)-{{S}_{C}}{{\omega }_{C}}\times ({{v}_{C}}+{{\omega }_{C}}\times {{r}_{B0}}) \\
\end{matrix} \right]$
(12)
式中ρ=xB1+SCrB1-xA1-SprA1，${\dot{\rho }}$=SC(vC+ω×rB1)-Sp(vp+ω×rA1)，SΔ=SCTSp，${\dot{S}}$Δ=SΔωp×-ωC×SΔ。
$\begin{align}
& {{R}_{2}}={{({{S}_{\Delta }}{{e}_{2}})}^{T}}({{\omega }_{C}}\times {{v}_{C}})+2{{(S_{\Delta }^{T}T{{v}_{C}}-{{v}_{p}})}^{T}}({{\omega }_{p}}\times {{e}_{2}})-{{({{x}_{C}}-{{x}_{p}})}^{T}}{{S}_{p}}({{\omega }_{p}}\times ({{e}_{2}}\times {{\omega }_{p}}))- \\
& {{e}_{2}}({{\omega }_{p}}\times {{v}_{p}})+{{({{({{\omega }_{C}}\times {{r}_{B0}})}^{T}}\times {{S}_{\Delta }}{{e}_{2}})}^{T}}{{\omega }_{C}}+{{({{r}_{B0}}\times {{S}_{\Delta }}({{\omega }_{p}}\times {{e}_{2}}))}^{T}}{{\omega }_{C}}- \\
& {{({{e}_{2}}\times S_{\Delta }^{T}({{\omega }_{C}}\times {{r}_{B0}}))}^{T}}{{\omega }_{p}}+{{(({{\omega }_{p}}\times {{e}_{2}})\times S_{\Delta }^{T}{{r}_{B0}})}^{T}}{{\omega }_{p}} \\
& {{R}_{3}}={{({{S}_{\Delta }}\varepsilon )}^{T}}({{\omega }_{C}}\times {{v}_{C}})+2{{(S_{\Delta }^{T}{{v}_{C}}-{{v}_{p}})}^{T}}({{\omega }_{p}}\times \varepsilon )-{{({{x}_{C}}-{{x}_{p}})}^{T}}{{S}_{p}}({{\omega }_{p}}\times (\varepsilon \times {{\omega }_{p}}))- \\
& \varepsilon ({{\omega }_{p}}\times {{v}_{p}})+{{({{({{\omega }_{C}}\times {{r}_{B0}})}^{T}}\times {{S}_{\Delta }}\varepsilon )}^{T}}{{\omega }_{C}}+{{({{r}_{B0}}\times {{S}_{\Delta }}({{\omega }_{p}}\times \varepsilon ))}^{T}}{{\omega }_{C}}- \\
& {{(\varepsilon \times S_{\Delta }^{T}({{\omega }_{C}}\times {{r}_{B0}}))}^{T}}{{\omega }_{p}}+{{(({{\omega }_{p}}\times \varepsilon )\times S_{\Delta }^{T}{{r}_{B0}})}^{T}}{{\omega }_{p}} \\
\end{align}$
则二体系统的动力学方程可以表示为
$\underbrace{\left[ \begin{matrix}
{{M}_{p}} & {{H}_{p}} & 0 & 0 \\
H_{p}^{T} & {{J}_{p}} & 0 & 0 \\
0 & 0 & {{M}_{C}} & {{H}_{C}} \\
0 & 0 & H_{C}^{T} & {{J}_{C}} \\
\end{matrix} \right]}_{N}\underbrace{\left[ \begin{align}
& {{{\dot{V}}}_{p}} \\
& {{\omega }_{p}} \\
& {{{\dot{V}}}_{C}} \\
& {{\omega }_{C}} \\
\end{align} \right]}_{{\ddot{q}}}=\underbrace{\left[ \begin{align}
& {{F}_{p}}-{{\omega }_{p}}\times ({{M}_{p}}{{V}_{p}}+{{H}_{p}}{{\omega }_{p}}) \\
& {{L}_{p}}-{{\omega }_{p}}\times ({{I}_{p}}{{\omega }_{p}}+{{H}_{p}}{{V}_{p}})-{{V}_{p}}\times ({{M}_{p}}{{V}_{p}}+{{H}_{p}}{{\omega }_{p}}) \\
& {{F}_{C}}-{{\omega }_{C}}\times ({{M}_{C}}{{V}_{C}}+{{H}_{C}}{{\omega }_{C}}) \\
& {{L}_{C}}-{{\omega }_{C}}\times ({{H}_{C}}{{V}_{C}})-{{V}_{C}}\times ({{H}_{C}}{{\omega }_{C}})F \\
\end{align} \right]}_{F}-{{S}^{T}}f_{q}^{T}\lambda $
(13)
式(13)可以简写为
$M\ddot{q}=F-{{B}^{T}}\lambda $
(14)
式中
$S=\left[ \begin{matrix}
{{S}_{p}} & 0 & 0 & 0 \\
0 & E & 0 & 0 \\
0 & 0 & {{S}_{C}} & 0 \\
0 & 0 & 0 & E \\
\end{matrix} \right]$
(15)
${{S}^{T}}f_{q}^{T}={{B}^{T}}$
(16)
将式(14)解出${\ddot{q}}$，并代入式(10)，可以显式求解出拉格朗日乘子
$\lambda ={{(B{{M}^{-1}}{{B}^{T}})}^{-1}}(B{{M}^{-1}}F+R)$
(17)
拉格朗日乘子具有物理意义，其模值表示约束力的大小。从而得到系统动力学方程为
$M\ddot{q}=F-{{B}^{T}}{{(B{{M}^{-1}}{{B}^{T}})}^{-1}}(B{{M}^{-1}}F+R)$
(18)
本文采用亚当斯积分法[16]进行数值求解，采用Baumgarte违约修正方法进行修正。
1.3 附加质量及气动模型 本文翼伞的附加质量采用Barrows估算方法[1]得到，如式(19)所示，翼伞附加质量矩阵中各项计算见表 1，下标c,p,r分别表示翼伞弦长/伞绳汇交点，翼伞俯仰运动中心，翼伞滚动中心。
$\left[ \begin{matrix}
M & -Mu_{rc}^{\times }+u_{pr}^{\times }{{S}_{2}} \\
u_{rc}^{\times }+{{S}_{2}}u_{pr}^{\times }M & {{J}_{c}} \\
\end{matrix} \right]$
(19)
表 1 附加质量矩阵各项
Table 1 Items of apparent mass matrix
翼伞气动力计算式及相关气动系数参考了X-38翼伞系统数据[17-18]，如表 2所示。
表 2 翼伞气动力模型
Table 2 Aerial model of parafoil
本文采用简单线性下偏操纵归航控制律，即翼伞后缘下偏输出量δ与航向偏差量Δα成正比关系，如表 3所示，航向偏差角是指测量的航向角与目标航向角(地面投影位置-目标点连线的夹角)。其中H(t)和R(t)分别为任意时刻t，系统离目标点高度和水平距离。
表 3 各飞行模式归航控制算法
Table 3 Homming control algorithm in flight modes
2 动力学仿真研究
2.1 仿真工况
翼伞的主要参数：翼型Clark-Y，c=5.55 m，b=14.4 m，e=0.15，R=11.5 m，θr=6°(不含剖面4°)；翼伞净重35 kg，落点目标均为地面惯性坐标系下(0,0,1 244)，各分析工况的其他初值输入条件如表 4所示。
表 4 各仿真工况初值输入条件
Table 4 Initial inputs in all simulation cases
2.2 仿真结果及分析 图 4～6分别给出了工况1～3的试验数据与仿真数据的对比分析结果，根据仿真结果，可以得出以下结论：
(1) 3个工况的系统高度随时间基本都是直线下降的，即下降速度基本维持在一个定值，伴随着不同程度的波动。这一特性可以应用于雀降机动的时机判断。
(2) 尽管仿真所用气动参数可能与实际有偏差，但3个工况的仿真结果与试验结果均基本吻合，说明了仿真模型、控制律逻辑和仿真程序具有一定的可信度。
(3) 工况1和3虽然初始条件接近，但试验和仿真数据显示飞行下降速度、水平速度都有较大差别，即图 4(c)中下降速度约为5 m/s，而图 6(c)中则约4.4 m/s。分析原因可能是操纵状态差异导致，工况1中因有航向偏差，不断执行较大的单边操纵，而工况3中则大部分时间操纵量较小，几乎处于自由直飞状态，如图 7所示。
(4) 工况2仿真与试验对比图中，航迹和航向角随时间的变化有偏差，试验结果不符合控制逻辑。分析原因是，图 5(d)显示试验第7~20 s出现导航数据中断，翼伞向-X向飞行一段(图 5(b)中(1 150,2 000)~(1 050,2 000))后才转弯向目标飞行，因此翼伞一直在初始点-目标点连线的右侧。
(5) 图 8为仿真得到的工况1和工况2中翼伞姿态的动态变化情况，结果表明翼伞在工况1条件下比较稳定，而在工况2中出现较长时期的俯仰振荡运动，说明了有效载荷重量影响翼伞的配平攻角，并且翼伞存在多个平衡点，而有些平衡点不稳定。
3 结束语 本文采用拉格朗日乘子法对一般翼伞-载荷系统建立了两体8DOF动力学仿真模型，考虑了翼伞附加质量特性。重建了某小型翼伞系统的3个工况飞行过程，并与相应的空投试验数据进行了对比，吻合度较好，发现并解释了试验过程中的某些运动特征，初步验证了仿真模型和仿真程序的有效性，说明仿真模型具有一定的可信度。获得更准确的气动参数将有助于提高仿真模拟的准确性。
