基于伞载多体动力学的降落伞开伞过程研究

贾华明 王文强 刘乃彬 林汝领 李少腾

(1.北京空间机电研究所, 北京 100094;2.中国航天科技集团有限公司航天进入、减速与着陆技术实验室, 北京 100094)

1 引言

降落伞作为一种典型的气动减速装置,在航空航天领域应用广泛。 随着中国载人航天、火星探测等重大专项工程的持续推进,对降落伞的减速性能提出了更高的要求[1-4],同时也给降落伞的理论研究和工程应用带来了新的挑战。

目前降落伞理论研究上的不足造成设计分析中某些缺陷无法避免,难点主要集中在降落伞开伞过程研究[5-7]。 降落伞从拉直到充气,伞衣会经历急剧的结构大变形,而且伞衣内外的气流流动状况非常复杂,材料属性和外界环境对降落伞的开伞过程都有重要影响,存在多种因素剧烈耦合、学科交叉的鲜明特点,这是降落伞开伞过程研究的重点和难点[8]。 目前降落伞开伞过程的研究方法有流固耦合力学方法(Fluid Structure Interaction,FSI)、计算流体力学方法(Computational Fluid Dynamics,CFD)和伞载多体动力学方法(Parachute Payload Dynamics,PPD)。 流固耦合力学方法是目前研究的热点,它同时考虑伞衣变形和周围流场的变化,可得到流场、结构、运动等多场参数,但模型建立复杂,计算消耗量巨大,依然存在诸多技术难题,其中解决负体积问题,提高计算效率是当前需要重点关注的问题,一般用于短时的瞬态分析[9-11]。 计算流体力学方法能够得到内外流场结构细节,计算消耗量中等,但不考虑伞衣变形,一般用于稳降状态下的气动性能和静稳定性能研究,其中外形、透气、是否考虑尾流是该方法分析问题的关键点[12-15]。 伞载多体动力学方法模型建立较简单,计算消耗量最小,适用于长时间的运动计算[16-19]。

由于流固耦合方法和流体力学方法建模复杂、计算量极大,尤其是大型伞系统,不能够满足快速设计迭代、性能分析和预测评估的工程需要。本文基于伞载多体动力学,以充气时间为自变量,将轴向动量方程作为参考模型,进行了降落伞开伞过程研究,并将计算结果与飞行试验数据进行了对比,用以验证降落伞开伞特性。

2 典型的降落伞开伞过程

如图1 所示,降落伞的开伞过程包括拉直、充气和稳降3 个阶段。 一般情况下,包装于伞包内的主伞是由引导伞将其拉出并拉直。 目前,航天器回收领域降落伞基本采用倒拉法,即先拉伞绳的拉直方法,这是因为在相同条件下,倒拉法的拉直力会小于顺拉法。 在倒拉法中首先由引导伞将伞包提起,使其与载荷分离,然后从伞包中先拉出伞绳,再逐步拉出伞衣。 从伞系统全长被拉直,进入充气阶段,空气团从伞衣底边进入伞衣内部,引起伞衣迅速充满,在工程上,为了降低开伞载荷,通常采用底边收口方式,这样在充气过程中就增加了收口工作阶段,由火工切割器切断收口绳后,伞衣底边快速张开,伞衣底边进气口迅速增大,直至伞衣完全充满,之后载荷在降落伞的作用下稳定下降[20]。

图1 降落伞开伞过程Fig.1 The typrocess of parachute opening

3 动力学模型

3.1 拉直阶段

3.1.1 物伞系统动力学方程

整个伞系统的拉直过程是引导伞和载荷二者相对运动的结果,从理论上看,是2 个变质量体的相对运动。 图2 所示为倒拉法开伞拉直过程的动力学模型,图中dm为拉直中的主伞系统微元质量,mw+me为载荷连同已拉直主伞系统的质量;mys+mv为引导伞、伞包连同未拉出的主伞系统的质量;Fsh为伞绳拉出阻力,FL为拉直力,θ为航迹倾斜角,L为伞系拉直距离。

图2 倒拉法开伞拉直过程Fig.2 lines-first deployment

为简化模型,做如下假设:

1)拉直过程中,载荷和引导伞每一瞬时的运动遵循同一条轨迹,航迹倾斜角θ按载荷计算;

2)物伞系统没有升力,在二维平面内运动;

3)在拉直过程中,不考虑伞绳伸长;

4)把引导伞、载荷和拉直中的主伞系统微元质量dm作为3 个质点处理。

拉直过程中运动方程如公式(1)所示:

式中,Dw为载荷的气动阻力,De为主伞系统已拉出部分的气动阻力,Dys为引导伞的气动阻力,Dv为主伞包的气动阻力,m＇为拉动的伞系统单位长度的质量,vw为载荷速度,vys为引导伞速度,t为时间,g为重力加速度。

拉直阶段的轨迹如公式(2)所示:

式中,xd为载荷的水平位移,yd为载荷的竖直位移。

拉直力计算公式如公式(3)所示:

为使计算结果接近实际情形,动力学模型中加入了参考风场数据模型,其中风向定义为逆风速矢量(即风的来向)与正北方向的夹角,以顺时针方向为正,并进行三次样条插值,获得任意高度上的风速和风向。

3.1.2 伞系质量变化规律

在拉直过程中,已拉出主伞系统质量me和未拉出主伞系统质量mv是随时间而变化的,倒拉法降落伞质量变化规律如公式(4)、(5)所示。

3.1.3 拉直阶段计算

根据式(1)、(2)的5 个微分方程组以及式(4)、(5)的伞系质量变化规律,以时间t为计算步长,可以计算出各个时间t内的各种运动学参数vw、vys、θ、xd、yd、L。 根据计算得出的各个时间t内的vw、vys值,代入公式(5)中可求出各个时间t内的伞绳张力,而最大的张力即为拉直力。

3.2 充气阶段

3.2.1 伞载组合体动力学方程

伞载组合体充气阶段如图3 所示,其中,Ds是降落伞阻力,GS是降落伞重力,FS是开伞载荷,DW是载荷阻力,GW为载荷重力,θ为伞载组合体航迹倾斜角。

图3 降落伞充气过程Fig.3 Inflation phase of parachute

为简化模型,做如下假设:

1)不考虑伞载组合体伸长;

2)不考虑伞载组合体的升力,在二维平面内运动;

3)伞载组合体是两质点运动,载荷的质量mw集中在载荷重心位置,伞的质量ms集中在伞衣底边中心位置。 充气阶段,伞的质心相对底边保持位置不变;

4)伞载组合体在充气阶段轴线呈直线。

充气过程中运动学方程如公式(6)所示:

其中,mw是载荷质量,ms是降落伞质量,mf是伞附加质量,dmf/dt是伞附加质量变化率,vwx是载荷水平速度,vwy是载荷竖直速度,t是时间,g是重力加速度。

开伞载荷算计算如公式(7)所示:

3.2.2 降落伞单级收口阻力面积

图4 是降落伞单级收口阻力面积示意图,(CA)sk是收口状态下阻力面积,(CA)S是全张满状态下阻力面积,tm1、tm2、tm3为充气阶段的3 个过程,ty是收口过程,其中ty=tm1+tm2。

图4 充气阶段阻力面积变化规律Fig.4 Variation of drag area in inflation

降落伞单级收口充气阶段分为3 个过程:tm1过程是伞开始充气到伞衣呈乌贼状外形,这个过程阻力面积与时间为线性;tm2过程阻力面积维持恒定值(CA)sk;tm3过程是结束收口至伞全张满,这时候阻力面积维持在另一个恒定值(CA)s,这个过程中阻力面积是时间的n次方,用公式(8)表示充气阶段阻力面积变化。

其中,(CA)是阻力面积;n是充气参数,一般在2～3 之间。

公式(9)为tm1过程中阻力面积与时间呈现的线性关系。

其中,λm1是设定参数,vL是伞开始充气速度。

公式(10)为tm3过程中阻力面积与时间的关系。

其中,λm3是设定参数,vDR是收口完成时伞载组合体速度。

3.2.3 充气阶段附加质量变化

充气阶段会受到降落伞附加质量的影响,附加质量由降落伞特征和大气参数等决定。 公式(11)为降落伞附加质量的计算公式。

其中,kf是设定参数,ρ是空气密度。

降落伞充气阶段附加质量变化率dmf/dt如公式(12)所示:

3.2.4 充气阶段计算

使用公式(8)～(10)的伞衣阻力面积变化和公式(11)、(12)的附加质量变化,按照公式(6)、(7)能够计算降落伞单级收口充满阶段所有运动参数和开伞动载。

3.3 稳定阶段

在伞衣充满后,物伞系统在气动阻力的作用下继续减速,直至系统达到稳定下降。 稳定下降阶段,伞衣已经充满后,降落伞的阻力面积和附加质量不再发生变化,即d(CA)/dt=0,dmf/dt=0。继续使用公式(6)、(7)、(11)就可以完成稳定下降阶段计算。

3.4 开伞过程计算

以时间t为自变量,采用Rouge-Kutta 法,可以进行开伞过程拉直阶段、充气阶段和稳定阶段的物伞系统全部运动学和动力学计算,充气阶段用到的3 个设定参数λm1、λm3、kf,可通过已获取的试验数据进行归纳整理得出。 某十字型降落伞充气过程动力学模型的3 个设定参数分别取为λm1=23,λm3=8,kf=0.5[21],本文建立的模型增加了多种伞型,如环帆伞、带条伞、平面圆伞和盘缝带伞,这大大提高了动力学模型的适用范围。

设定参数λm1和λm3与伞型密切相关,美国的标准环帆伞λm1=37.3,λm3=4.42;美国的改型环帆伞λm1=26.0,λm3=6.80;中国的某环帆伞λm1=23.5,λm3=4.3[22]。

设定参数kf与伞型和透气量密切相关,美国阿波罗环帆伞kf=0.66;中国某平面圆伞kf=0.41[23]。

4 验证与分析

某航天器返回舱采用了带条伞+环帆伞降落伞方案,飞行试验后通过对数据的分析整理,分别拟合出带条伞和环帆伞3 个参数λm1、λm3、kf的取值,将其固化在建立的开伞过程的动力学模型中,进行编程计算,根据多次飞行试验数据对模型中的一些参数进行反复修正优化,计算结果精度,尤其是开伞动载得到了很大提高,更加接近于真实情况。

利用参数固化后的动力学模型对某型号降落伞某次飞行试验开伞过程进行计算,降落伞特性参数见表1 所示,开伞工况见表2 所示。

表1 降落伞参数Table 1 Parachute parameters

表2 计算工况Table 2 Simulation status

计算结果与试验数据的对比曲线见图5～8所示,其中图5 为弹道轨迹的对比曲线,图6 为高度变化对比曲线,图7 为速度变化对比曲线,图8为航迹倾斜角变化对比曲线。

图6 高度对比Fig.6 Comparison of altitude

图8 航迹倾斜角对比Fig.8 Comparison of track angle

图5 中,将开伞时刻的载荷经纬度与高度置为零,减速伞工作阶段轨迹较陡,试验数据曲线与计算数据曲线一致性较好。 主伞工作阶段,轨迹变得较平滑,试验数据与计算数据出现明显偏差,这是因为仿真计算使用的是参考大气数据,与试验场当天真实风场数据偏差较大,主伞工作阶段的载荷速度已经降到很低,风的影响会产生较大的水平位移,从而造成了试验数据与计算数据的偏差。

在图6、图7 的高度和速度的对比曲线中,试验与计算数据的变化趋势一致性很好。 计算曲线较为平滑,而试验数据曲线有一定的峰值变化,这主要是由于大气参数偏差和气象条件的影响。

图7 速度对比Fig.7 Comparison of velocity

在图8 航迹倾斜角对比曲线中,试验数据曲线和计算数据曲线存在一些差异,减速伞及主伞收口阶段,二者的差异不大,在-55°～-70°范围内波动,主伞解除收口后,载荷处于稳定下降阶段,计算数据在参考大气数据的情况下,航迹倾斜角保持在-35°,而试验数据在实际风场的情况下,载荷随风飘动,航迹倾斜角在-20°左右波动。

开伞动载是降落伞开伞过程中的一个重要性能指标参数,是降落伞结构设计的依据,也是回收系统接口设计的输入。 表3 分别对减速伞工作阶段和主伞工作阶段的开伞动载进行了试验数据和计算数据的比较,计算偏差均小于10%,计算结果很好地体现了开伞特性,与试验结果的一致性较好。 其中减速伞工作段的开伞载荷偏差要高于主伞工作段,最大为减速伞收口段,偏差率为8.6%。 最小为主伞全展开段,偏差率为3.9%。减速伞和主伞分别采用了带条伞和环帆伞2 种不同的伞型,所以计算过程中的λm1、λm3、kf也分别采用不同的取值。

表3 开伞动载对比Table 3 Comparison of opening force

5 结论

本文采用伞载多体动力学方法对某型号降落伞的开伞过程进行了计算,并将计算数据与试验数据进行了对比分析,得出如下结论:

1)减速伞工作段,试验数据和计算数据的弹道和航迹倾斜角一致性较好,主伞工作段,尤其是主伞全展开后,载荷速度明显降低,由于参考大气数据与真实风场数据偏差较大,造成试验数据和计算数据产生明显偏差,尤其是航迹倾斜角。

2)在速度和高度的对比曲线中,试验数据和计算数据的变化趋势一致性较好。

3)减速伞和主伞工作阶段计算得到的开伞动载与试验数据的偏差均小于10%,计算结果很好地体现了开伞特性,主伞采用的环帆伞伞型的计算精度要高于减速伞采用的带条伞伞型。

4)通过计算数据与试验数据的对比,证明基于伞载多体动力学的开伞过程的动力学模型是正确的,λm1、λm3、kf3 个经验参数的取值是可信的,由于参考大气数据与真实风场数据存在偏差,使得降落伞开伞过程计算方法有一定的局限性。

5)鉴于降落伞开伞过程的复杂性,尤其是大型降落伞,采用流固耦合和流体力学的数值仿真计算方法,计算量极大、耗时长、可行性低,工程上不易实现。 而本文基于伞载多体动力学方法,简单方便、用时短,能够快速指导降落伞的设计迭代。