翼伞具有极优的操纵能力，在精确空投、伞兵空降及运动表演等领域都得到了广泛应用^[1-2]。翼伞通过下拉双侧操纵绳使其气动外形发生变化，升力、阻力系数迅速增大，前进速度和垂直速度在着陆前迅速减小^[3]，使之实现安全着陆的过程称为雀降过程。雀降操纵本质上是一种小心操纵的动力失速 ^[4]。

目前，国内外学者多采用单体六自由度(6 degree of freedom，6DOF)模型开展翼伞雀降着陆阶段的研究^[5-6]。早期开展此项工作代表人物Goodrick等^[7-8]提出并建立了系统平动及转动的6DOF模型，对翼伞系统雀降过程进行了仿真研究，通过与试验翼伞下降率、速度等对比，验证了模型的准确性。国内一些学者利用6DOF模型开展了翼伞雀降着陆阶段的研究。卞琨^[9]基于雀降着陆操纵时长及平均垂直速度，得到了翼伞进入雀降着陆段的平均高度；胡容^[10]针对2种雀降实施高度下着陆速度的对比，得到了相对优化的雀降高度；张兴会等^[11]探究了逆风对雀降高度的影响；熊菁^[12]则针对翼伞部分结构参数对雀降性能的影响进行了较为深入的研究。

6DOF模型将翼伞与载荷视为固连的刚体，没有考虑到伞-载两体之间的相对运动状态变化。为更好地探究翼伞-载荷体之间的相对运动变化对翼伞雀降着陆阶段的影响规律，一些学者开始进行翼伞系统多体动力学模型下的雀降研究^[13-16]。如杨华等^[17]基于翼伞系统纵向对称面内的4DOF模型，根据载荷体接地机械能及速度变化得到翼伞安全雀降的高度及时间范围；李东^[18]基于9DOF模型，根据速度的合理范围得到相应的雀降高度范围。上述研究可以得到雀降操纵实施的时间或高度范围，但鲜有对最优雀降实施高度方面的研究报道。

本文根据文献[19]建立了翼伞-载荷两体系统9DOF动力学模型，以垂直着陆速度最小为优化目标，采用时间二分法开展了雀降实施高度的优化研究，探究了初始工况、翼载荷、着陆高度及安装角对于最优雀降实施高度、着陆速度的影响，并完成公式拟合，为工程设计提供了理论依据。

1 翼伞系统动力学建模 1.1 坐标系定义

本文动力学方程、运动学方程之间的关系涉及到3个坐标系，包括惯性坐标系O_eX_eY_eZ_e和两类体轴坐标系(伞体坐标系O_pX_pY_pZ_p、载荷坐标系O_bX_bY_bZ_b)，如图 1所示。

惯性坐标系到体轴系的转换矩阵为

$ \boldsymbol{T}_{i}=\left(\begin{array}{ccc} C_{\theta_{i}} C_{\varphi_{i}} & C_{\theta_{i}} S_{\varphi_{i}} &-S_{\theta_{i}} \\ S_{\phi_{i}} S_{\theta_{i}} C_{\varphi_{i}}-C_{\phi_{i}} C_{\varphi_{i}} & S_{\phi_{i}} S_{\theta_{i}} S_{\varphi_{i}}+C_{\phi_{i}} C_{\varphi_{i}} & S_{\phi_{i}} C_{\theta_{i}} \\ C_{\phi_{i}} S_{\theta_{i}} C_{\varphi_{i}}+S_{\phi_{i}} S_{Q_{i}} & C_{\phi_{i}} S_{\theta_{i}} S_{\varphi_{i}}-S_{\phi_{i}} C_{\varphi_{i}} & C_{\phi_{i}} C_{\theta_{i}} \end{array}\right) $

(1)

式中：C_*=cos(*)，S_*=sin(*)；ϕ_i为滚转角，θ_i为俯仰角，φ_i为偏航角，其中下标i为不同的体轴系。

1.2 9DOF动力学模型

本文将翼伞系统中的翼伞、载荷看作是两个中间由球铰连接的刚体，球铰只传递力不传递力矩，翼伞及载荷之间不发生相对移动。因此，翼伞系统的运动模型可简化为系统在连接点(图 1中C点)处的3自由度平动以及翼伞、载荷绕C点的3自由度转动。由此建立刚柔组合系统9DOF动力学模型，动力学方程为

$ \left[\begin{array}{cccc}-\boldsymbol M_{\mathrm{b}} \tilde{\boldsymbol R}_{\mathrm{cb}} & 0 & \boldsymbol M_{\mathrm{b}} T_{\mathrm{b}} & \boldsymbol{T}_{\mathrm{b}} \\ 0 &-\left(\boldsymbol M_{\mathrm{p}}+\boldsymbol M_{\mathrm{F}}\right) \tilde{\boldsymbol{R}}_{\mathrm{cp}}&\left(\boldsymbol M_{\mathrm{p}}+\boldsymbol M_{\mathrm{F}}\right) \boldsymbol T_{\mathrm{p}} &-\boldsymbol T_{\mathrm{p}} \\ \boldsymbol I_{\mathrm{b}} & 0 & 0 &-\tilde{\boldsymbol R}_{\mathrm{ch}} \boldsymbol T_{\mathrm{b}} \\ 0 & \boldsymbol I_{\mathrm{p}}+\boldsymbol I_{\mathrm{F}} & 0 & \tilde{\boldsymbol R}_{\mathrm{cp}} \boldsymbol T_{\mathrm{p}}\end{array}\right]\left[\begin{array}{l}\dot{\boldsymbol{\varOmega}}_{\mathrm{b}} \\ \dot{\boldsymbol{\varOmega}}_{\mathrm{p}} \\ {\dot{\boldsymbol V}_{\mathrm{c}}} \\ \boldsymbol F_{\mathrm{c}}\end{array}\right]=\left[\begin{array}{l}\boldsymbol{B}_{1} \\ \boldsymbol{B}_{2} \\ \boldsymbol{B}_{3} \\ \boldsymbol{B}_{4}\end{array}\right] $

(2)

式中$\boldsymbol{B}$为翼伞系统外载荷矩阵，具体为

$\left[\begin{array}{l}\boldsymbol B_{1} \\ \boldsymbol B_{2} \\ \boldsymbol B_{3} \\ \boldsymbol B_{4}\end{array}\right]=\left[\begin{array}{c}\left[\begin{array}{c} \boldsymbol{F}_{\mathrm{b}}^{\mathrm{A}}+\boldsymbol{F}_{\mathrm{b}}^{\mathrm{G}}-\tilde{\boldsymbol{\varOmega}}_{\mathrm{b}} \boldsymbol{M}_{\mathrm{b}} \tilde{\boldsymbol{\varOmega}}_{\mathrm{b}} \boldsymbol{R}_{\mathrm{cb}} \\ \boldsymbol{F}_{\mathrm{p}}^{\mathrm{A}}+\boldsymbol{F}_{\mathrm{p}}^{\mathrm{G}}-\tilde{\boldsymbol{\varOmega}}_{\mathrm{p}}\left(\boldsymbol{M}_{\mathrm{p}}+\boldsymbol{M}_{\mathrm{F}}\right) \tilde{\boldsymbol{\varOmega}}_{\mathrm{p}} \boldsymbol{R}_{\mathrm{cp}}+\boldsymbol{M}_{\mathrm{F}} \tilde{\boldsymbol{\varOmega}}_{\mathrm{p}} \boldsymbol T_{\mathrm{p}} \boldsymbol{V}_{\mathrm{c}}-\tilde{\boldsymbol{\varOmega}}_{\mathrm{p}} \boldsymbol{M}_{\mathrm{F}} \boldsymbol{T}_{\mathrm{p}} \boldsymbol{V}_{\mathrm{c}} \\ \boldsymbol{M}_{\mathrm{b}}^{\mathrm{A}}-\tilde{\boldsymbol{\varOmega}}_{\mathrm{b}} \boldsymbol I_{\mathrm{b}} \boldsymbol{\varOmega}_{\mathrm{b}} \\ \boldsymbol{M}_{\mathrm{p}}^{\mathrm{A}}-\tilde{\boldsymbol{\varOmega}}_{\mathrm{p}}\left(\boldsymbol I_{\mathrm{p}}+\boldsymbol I_{\mathrm{F}}\right) \boldsymbol{\varOmega}_{\mathrm{p}} \end{array}\right]\end{array}\right] $

(3)

运动学关系式为：

$ \left[\begin{array}{c} \dot{x}_{\mathrm{e}} \\ \dot{y}_{\mathrm{e}} \\ \dot{z}_{\mathrm{e}} \end{array}\right]=\left[\begin{array}{l} \dot{x}_{\mathrm{c}} \\ \dot{y}_{\mathrm{c}} \\ \dot{z}_{\mathrm{c}} \end{array}\right]=\left[\begin{array}{c} u_{\mathrm{c}} \\ v_{\mathrm{c}} \\ w_{\mathrm{c}} \end{array}\right] $

(4)

$ \left[\begin{array}{l} u_{\mathrm{p}} \\ v_{\mathrm{p}} \\ w_{\mathrm{p}} \end{array}\right]=\boldsymbol{T}_{\mathrm{p}}\left[\begin{array}{c} u_{\mathrm{c}} \\ v_{\mathrm{c}} \\ w_{\mathrm{c}} \end{array}\right]+\dot{\boldsymbol{\varOmega}}_{\mathrm{p}} \boldsymbol{R}_{\mathrm{cp}} $

(5)

$ \left[\begin{array}{l} \dot{\phi}_{\mathrm{b}} \\ \dot{\theta}_{\mathrm{b}} \\ \dot{\varphi}_{\mathrm{b}} \end{array}\right]=\left[\begin{array}{lll} 1 & S_{\phi_{\mathrm{b}}} t_{\theta_{\mathrm{b}}} & C_{\phi_{\mathrm{b}}} t_{\theta_{\mathrm{b}}} \\ 0 & C_{\phi_{\mathrm{b}}} &-S_{\phi_{\mathrm{b}}} \\ 0 & S_{\phi_{\mathrm{b}}} / C_{\theta_{\mathrm{b}}} & C_{\phi_{\mathrm{b}}} / C_{\theta_{\mathrm{b}}} \end{array}\right]\left[\begin{array}{l} p_{\mathrm{b}} \\ q_{\mathrm{b}} \\ r_{\mathrm{b}} \end{array}\right] $

(6)

$ \left[\begin{array}{l} \dot{\phi}_{\mathrm{p}} \\ \dot{\theta}_{\mathrm{p}} \\ \dot{\varphi}_{\mathrm{p}} \end{array}\right]=\left[\begin{array}{lll} 1 & S_{\phi_{\mathrm{p}}} t_{\theta_{\mathrm{p}}} & C_{\phi_{\mathrm{p}}} t_{\theta_{\mathrm{p}}} \\ 0 & C_{\phi_{\mathrm{p}}} &-S_{\phi_{\mathrm{p}}} \\ 0 & S_{\phi_{\mathrm{p}}} / C_{\theta_{\mathrm{p}}} & C_{\phi_{\mathrm{p}_{\mathrm{p}}}} / C_{\theta_{\mathrm{p}}} \end{array}\right]\left[\begin{array}{l} p_{\mathrm{p}} \\ q_{\mathrm{p}} \\ r_{\mathrm{p}} \end{array}\right] $

(7)

式中：$\boldsymbol{F}_{*}^{\mathrm{A}}$为气动力矢量；$\boldsymbol{F}_{*}^{\mathrm{G}}$为重力矢量；$\boldsymbol{M}_{*}^{\mathrm{A}}$为气动力矩矢量；$\boldsymbol{M}_{*}$为质量矢量；$\boldsymbol{I} *$为转动惯量矢量；$\boldsymbol{V}_{*} $为速度矢量；$\boldsymbol{\varOmega}_{*}$为角速度矢量；$\boldsymbol{R} _* $为距离矢量；$\tilde{\boldsymbol{\varOmega}}_{*}、\tilde{\boldsymbol{R}}_{*}$分别为各自的叉乘矩阵；t_*=tan(*)；下标p为翼伞；下标b为载荷；下标c为伞载系统连接点；下标e为地面。翼伞结构示意如图 2所示。

1.3 气动力模型

雀降操纵会引起翼伞气动外形改变，从而使翼伞气动特性发生改变，因此气动力系数和气动力矩系数计算过程中需加入双侧操纵量的影响，计算公式如下：

$ \left\{\begin{array}{l} \mathrm{C}_{L}=\mathrm{C}_{L}\left(\alpha_{p}, \delta_{e}\right) \\ \mathrm{C}_{D}^{p}=\mathrm{C}_{D}^{p}\left(\alpha_{p}, \delta_{e}\right) \end{array}\right. $

$ \begin{align*} & \left\{\begin{array}{l} C_{X}=\left(-C_{\mathrm{D}}^{\mathrm{p}} u_{\mathrm{p}}+C_{\mathrm{L}} w_{\mathrm{p}}\right) / V_{\mathrm{p}} \\ C_{Y}=C_{Y_{\beta}} \beta+C_{Y_{r}}\left(r_{\mathrm{p}} b / 2 V_{\mathrm{p}}\right) \\ C_{Z}=\left(-C_{\mathrm{D}}^{\mathrm{p}} w_{\mathrm{p}}-C_{\mathrm{L}} u_{\mathrm{p}}\right) / V_{\mathrm{p}} \end{array}\right. \\ & \left\{\begin{array}{l} C_{l}=C_{l_{\beta}} \beta+C_{l_{p}} \frac{p_{\mathrm{p}} b}{2 V_{\mathrm{p}}}+C_{l_{r}} \frac{r_{\mathrm{p}} b}{2 V_{\mathrm{p}}} \\ C_{m}=C_{m c / 4}\left(\alpha_{\mathrm{p}}, \delta_{\mathrm{e}}\right)+C_{m_{q}} \frac{p_{\mathrm{p}} c}{2 V_{\mathrm{p}}} \\ C_{n}=C_{n_{\beta}} \beta+C_{n_{\mathrm{p}}} \frac{p_{\mathrm{p}} b}{2 V_{\mathrm{p}}}+C_{n_{r}} \frac{r_{\mathrm{p}} b}{2 V_{\mathrm{p}} b} \end{array}\right. \end{align*} $

(8)

式中：δ_e为双侧下偏量，C_X、C_Y、C_Z为翼伞体轴坐标系下的气动力系数，C_l、C_m、C_n为翼伞的滚转、俯仰、偏航力矩系数，C_L、Cp_D、C_mc/4为采用插值拟合法拟合得到和攻角、双侧下偏量同时相关的升力、阻力、1/4弦线处俯仰力矩系数，其他气动数据取值为常数。

2 翼伞系统动力学仿真 2.1 研究对象及计算工况

为验证本文模型及方法的准确性，以文献[19]中翼伞系统为研究对象采用龙格库塔方法进行了编程计算。其中，翼伞系统结构参数见表 1，翼伞、翼载荷质量分别为5、135 kg，式(7)中的气动参数和文献[19]一致，初始运动参数见表 2。为更好模拟雀降操纵，本文设计如下仿真工况：在5 s时对翼伞施加双侧下拉操纵，操纵下偏量为100%，30 s时操纵下偏量回复至0，操纵实施时间均为1 s，计算总时长60 s。

2.2 结果分析

翼伞速度和伞载姿态的动态变化如图 3所示。由图 3可知，本文的仿真结果与文献结果规律一致，翼伞的运动特性变化与操纵绳下拉相互对应，翼伞系统工作过程分为无操纵稳定滑翔阶段、双侧全下拉减速雀降阶段、操纵绳回复至稳定滑翔3个阶段。

操纵绳在第5 s下拉后，翼伞的升力、阻力、俯仰力矩系数均增加，从而导致翼伞的前向速度和垂相速度减小，俯仰角增加，翼伞垂向速度在第6 s达到最小值，该值被称作最优着陆速度，此时被认为是雀降操纵的最佳着陆点，也是确定最优雀降高度的依据；之后由于翼伞的攻角急速增加，升力系数减小，翼伞的垂向速度逐渐增加，俯仰力矩系数逐渐减小从而俯仰角减小。操纵绳在第30 s恢复后，翼伞的升力系数增加，阻力、俯仰力矩系数减小，从而导致翼伞的前向速度增加、垂相速度减小、俯仰角减小，之后由于俯仰力矩系数增加，俯仰角增加，最终恢复至初始值。

表 3为稳定时各阶段运动特征量和文献的对比，最大误差为8.20%。上述结果验证了本文雀降操纵两体动力学模型的准确性。

3 雀降操纵优化高度确定方法 3.1 优化目标及求解方法

根据翼伞系统动力学仿真的分析结果，雀降操纵后系统的垂直速度会迅速降低到最小值随后增加，为保证着陆安全，应该确保翼伞系统以载荷体最小垂直速度着陆。本文以载荷体垂直着陆速度最小为优化目标，采用时间二分法获得最优雀降实施高度，技术方案如图 4所示，若载荷体垂直速度达到最小时的高度和预设着陆高度之间差值的误差小于0.01 m则可认为符合优化目标。

翼伞雀降采用一次操纵、“匀速下拉”模式^[20]，操纵总时长为3 s，操纵量变化满足:

$ \begin{equation*} \delta_{s}=\frac{t-t_{k}}{3}, \quad t_{k} \leqslant t \leqslant t_{k}+3 \end{equation*} $

(9)

3.2 雀降操纵实施高度影响分析

为探究翼伞系统滑翔运动状态对雀降实施高度的影响，本文以翼伞系统为对象，设计了10种计算工况，采用稳定滑翔阶段初始工况作为稳定状态计算工况，同时设计初始非稳定状态，具体见表 4。

通过计算得到了最优雀降实施高度随不同初始高度H₀的变化情况如图 5所示。随着高度的变化，初始稳定工况和非稳定工况的雀降高度集中在7.02 m左右；同时仿真结果中不同初始工况下的雀降时长、着陆速度均相同，分别为2.63 s、1.95 m/s。对同一对象无外在干扰的情况下，不同的初始运动参数最终会稳定至同一滑翔状态，本文由于算法的影响，不同初始工况下的雀降高度存在较小的波动，平均保持在7.02 m左右。

假设翼伞系统雀降操纵前处于稳定滑翔状态，初始高度为1 km开始计时，若雀降提前实施或延迟实施3 m，操纵时刻虽然相差不大但着陆速度却相差较大(如表 5及图 6所示)。由此可知，雀降操纵时机对着陆速度的影响较大，需合理判断雀降实施高度以更好地保证载荷体的无损着陆。

4 雀降操纵影响因素研究 4.1 翼载荷影响

为探究翼载荷m_b对雀降操纵的影响，设置翼载荷分别为55~215 kg的5组对比工况，初始工况见表 2，得到结果如图 7所示。

随着翼载荷的增加，滑翔速度增加，雀降操纵时的高度损失增加，因此需要在更高的高度实施雀降操纵，相对雀降高度随之增加(如图 7所示)，且相对雀降高度与翼载荷的1/2次方成正比，采用最小二乘法进行拟合，得到如式(10)所示的计算模型；当翼载荷越大时，系统重力越大，即翼伞竖直向上的气动力越大，着陆速度增加，且着陆速度与翼载荷的1/3次方成正比，采用最小二乘法进行拟合，得到如式(11)所示的计算模型。

$ h_{\text {land }}=0.6181 m_{\mathrm{b}}^{1 / 2}, R^{2}=0.998 $

(10)

$w_{\text {land }}=0.3784 m_{\mathrm{b}}^{1 / 3}, R^{2}=0.998 $

(11)

为验证相对雀降高度和着陆速度模型的准确性和适用性，本文选取大(250 kg)、中(100 kg)、小(20 kg) 翼载荷进行雀降数值计算，所得到的结果与式(10)、(11)计算结果进行对比，结果如表 6及图 8所示。从图 8中可以看出，本文提出的相对雀降高度和着陆速度计算模型与仿真计算结果较为符合，最大误差小于4.00%，表明该模型在全翼载荷范围内均有良好的适用性。

4.2 着陆高度影响

考虑到空气密度的变化，设置5组不同着陆高度(0~4 km)完成着陆高度对雀降性能的影响，仿真初始高度为5 km，初始运动参数见表 2，得到翼伞系统相对雀降高度和着陆速度变化如图 9所示。

随着着陆高度减小，空气密度增大，近地面的滑翔速度减小，因此雀降操纵时的高度损失减小，雀降实施高度减小；同时由于越小的着陆海拔高度雀降操纵前滑翔速度越小，操纵后着陆速度越小。分析得到相对雀降高度、着陆速度分别和密度的-1/2、-1/3次方成正比，计算模型如式(12)、(13)所示。

$ h_{\text {land }}=7.8590 \rho^{-1 / 2}, R^{2}=0.987 $

(12)

$ w_{\text {land }}=2.0800 \rho^{-1 / 3}, R^{2}=0.999 $

(13)

选取大(4.5 km)、中(1.5 km)、小(-0.2 km)着陆高度所对应的密度进行雀降数值计算，所得到的结果对比如表 7及图 10所示。从图 10中可以看出，本文提出的相对雀降高度和着陆速度计算模型与仿真计算结果最大误差小于2.00%，表明该模型在全着陆高度范围内均有良好的适用性。

4.3 安装角影响

翼伞安装角是定位翼伞相对于载荷体位置的重要参数，工程中其大小可通过改变伞绳长度来调节。采用初始工况见表 2，通过仿真计算得到不同安装角下翼伞系统着陆速度、相对雀降高度的变化如图 11所示，高度随攻角和气动力的变化曲线如图 12、13所示。

随着安装角的增加，稳定滑翔阶段平衡攻角减小(如图 12所示)，从而引起翼伞气动特性改变，升、阻力系数减小，由升、阻力计算公式及系统平衡可知稳定阶段滑翔速度增大，同理可知雀降实施操纵提前，相对雀降高度增加(如图 11所示)，且相对雀降高度与平衡攻角的-1/2次方成正比，得到如式(14)所示的计算模型。雀降操纵后，不同安装角下气动力变化幅度接近(如图 13所示)，着陆速度随安装角的增加先减小后增加，在此范围内差值小于0.1 m/s，几乎可以忽略不计。

$\begin{equation*} h_{\text {land }}=17.8000 \alpha^{-1 / 2}, R^{2}=0.983 \end{equation*} $

(14)

选取大(15°)、中(10°)、小(5°)安装角进行雀降数值计算，所得到的结果对比如表 8及图 14所示。从图 14中可以看出，该计算模型与数值计算结果较为符合，最大误差小于4.00%，表明该模型在全安装角范围内均有良好的适用性。

5 结论

1) 翼伞系统雀降高度需准确实施，提前或延迟实施均会引起着陆速度的增加。

2) 初始运动参数对最优雀降高度和着陆速度几乎没有影响；随着翼载荷和着陆高度的增加，系统雀降实施高度和着陆速度增加，安全性降低；随着安装角的增加，雀降实施高度增加。

3) 相对雀降高度及着陆速度的拟合计算模型结果和仿真结果最大误差小于4.00%，计算模型具有良好的适用性。