民机弹射降落伞离机充气过程和安全分析

2022-05-10张延泰孙建红许常悦冯传奇

空气动力学学报 2022年2期

张延泰，孙建红，2，*，侯斌，许常悦，冯传奇

(1. 南京航空航天大学飞行器环境控制与生命保障工业和信息化部重点实验室，南京 210016；2. 南京航空航天大学民航学院，南京 211106；3. 中国飞行试验研究院，西安 710089)

0 引言

传统飞行数据记录器（flight data recorder, FDR）都是固定在机身中。如果客机在跨洋飞行中发生空难，FDR会随飞机残骸沉入海底。马航370航班失联后，南京航空航天大学民机救生团队与中国商飞针对民航客机联合研发了分离式应急记录跟踪系统，该系统又被称为报信者（Harbinger, HBG）系统[1-4]。系统中装备了可快速弹射的小型降落伞，可以依靠气动力在飞机处于紧急状态时将应急飞行数据记录器（emergency flight data recorder, EFDR）拉出飞机。随后利用降落伞和气囊气动阻力，降低EFDR的坠落速度，减小坠落冲击。如果降落伞在充气过程中剐蹭机身，则可能导致伞衣破损，气囊漏气，引起气囊降落伞失效，设备损坏。研究不同工况下的降落伞充气过程，分析伞的运动轨迹，判断其工作过程中的安全性对设计理论和工程应用都有重要意义。

在空降空投领域，研究手段主要是空投试验、风洞试验和数值模拟。但在涉及降落伞的离机过程时，空投试验不仅费用昂贵，且载机还有一定安全风险，所以在设计研究阶段往往以风洞试验和数值模拟为主。

在人员离机过程的研究中，刘洋等[5]利用风洞试验研究了民航飞机试飞员应急离机的运动轨迹，分析了离机方案的安全性。刘晓宇等[6]通过动态嵌套网格技术模拟了人员跳伞离机的过程。在运输机空投货物的研究中，Bergmann等[7]针对降落伞的离机过程进行了风洞试验，试验中将降落伞简化为刚性半球壳，测量了降落伞的运动轨迹。Schade等[8]采用嵌套网格技术模拟了试验[7]中的空投过程，但发现过程中降落伞离机的初始阶段误差较小，随后误差逐渐增大，运动轨迹和俯仰角的最大误差约为18%。Roosenboom等[9]进一步完善了空投过程的风洞实验，试验中同样采用了刚性降落伞模型，应用PIV技术测量了飞机尾流场。Geisbauer等[10]对比了数值模拟与试验[9]中的流场结构，验证了数值模拟方法的可靠性。为了研究飞机尾流对降落伞气动力的影响，Tezduyar等[11]采用多计算域（multi-domain method, MDM）方法，分别将飞机和降落伞放置于不同的计算域内，模拟了柔性降落伞在飞机远尾流区中的运动。Serrano等[12]将刚性牵引伞固定在运输机近尾流区中，采用数值模拟方法分析了运输机尾流对牵引伞气动力的影响。除此以外，近年来美国陆军纳蒂克士兵研发与工程中心（NSRDEC）和美国空军学院高性能计算中心合作[13-14]，采用大规模数值模拟的方法（最大网格数约为2亿），研究了空投空降过程中飞机附近的流场结构，并分析了飞机尾流对伞兵、货物和牵引伞气动特性的影响。其中Bergeron等[14]将刚性牵引伞模型固定在C17运输机尾部舱门后方，分析了牵引伞前方来流的特点，研究了运输机尾流对牵引伞的气动力影响。

可以看到，研究中往往将人员、货物和降落伞等简化为刚性物体。虽然刚性物体从机身离机研究较多，但柔性、随机、小质量载荷、大变形的降落伞在民机弹射离机研究还不多见。本文针对从机身下方弹射离机的降落伞充气过程进行了数值研究，分别模拟了有/无机身影响下的伞衣充气过程，分析了伞衣剪切层和近尾迹区的流场特征，以及机身和来流角度对伞衣充气过程的影响，给出了降落伞安全离机的临界迎角。

1 数值计算方法与验证算例

针对伞衣充气过程中的柔性体大变形流固耦合问题，采用ALE方法进行数值模拟。HBG系统降落伞最大开伞速度为160 m/s，马赫数约为0.46，可压缩流体控制方程为：

其中，t为时间，ρ 为流体密度，ui为流体速度张量，p为流体压力，fi为体积力张量，ε为流体比内能， Φ为耗散函数，T为流体温度。

由于ALE方法在欧拉和拉格朗日坐标之外引入了参考坐标，所以在ALE方法中，材料域的节点信息需要被映射到空间域内，其流场变量的输运方程为：

其中，φ为流场变量， φ0为当前时间步流场变量的初始值，为材料速度张量，当=0 时对应欧拉方法，当=ui时对应拉格朗日方法。

本文采用罚函数方法进行流固耦合计算。罚函数算法会时刻监控流体与结构体之间是否发生穿透，如果流体与结构体没有发生穿透，则流固耦合力为零；如果流体与结构体发生穿透，则同时对流体和结构体施加大小相同方向相反的耦合力。

为了验证降落伞展开过程，采用Heinrich等[15]的风洞试验对数值模拟方法进行验证。试验中采用了C9伞的缩比模型，伞衣名义直径D0为0.914 m，前置体重0.251 kg，来流速度15.3 m/s。降落伞折叠模型如图1所示，计算域尺寸为 6D0×6D0×9D0。在风洞试验中，前置体悬挂于风洞外部，通过滑轮绳索系统与降落伞相连，无前置体尾流影响。因此，数值模拟中前置体简化为质点。

图1 折叠后的降落伞网格Fig. 1 The folded parachute mesh

图2给出了当前计算的开伞动载gn与三次试验数据的对比。可以看出，在充气初始阶段，计算和试验数据趋势一致但有差异，主要缘于随机性影响[16]。在充气过程后期，伞衣接近充满状态，这时随机因素的影响会明显减弱，计算结果和试验数据吻合更好。

图2 试验和数值模拟中的开伞动载Fig. 2 The parachute opening loads in experiments and numerical simulation

2 计算模型

机身和降落伞模型如图3所示，伞型为平面圆形伞，伞衣名义直径为1.9 m。坐标系x轴为流向，y轴为展向，z为横向。离机后降落伞位于机身下方贴近蒙皮的位置，伞衣初始状态如图1所示。因为在开伞过程中前置体体积较小，所以将前置体简化为一个质点，质量为4 kg。计算域采用O形网格（图4），计算域直径为 26L，其中L为机身长度，在伞衣充气区域进行局部加密，总网格数约为5.4×106，为保证质量守恒，计算域边界采用速度入口。

图3 机身与降落伞模型示意图Fig. 3 Schematic model of the fuselage and parachute

图4 机身对称面局部网格Fig. 4 The computational mesh at the symmetrical plane

根据美国国家运输安全委员会（National Transportation Safety Board, NTSB）在1990年～2020年的事故报告，取有人员死亡的大型客机空难案例进行分析。在83%的案例中，飞机坠毁前的空速处于40～160 m/s的区间。模拟工况来流速度分别为40、70、100、130、160 m/s。初始时刻，降落伞弹射离机速度10 m/s，方向为z轴负方向，本文所述迎角均为飞机来流迎角。

3 结果和分析

3.1 流场分析

为方便分析，定义在开伞动载峰值时刻伞衣底边中心点A沿降落伞径向与飞机机身的距离为δ，以及伞衣底边投影半径，其中Ap为伞衣底边投影面积，如图5所示。为获得迎角的影响，分析了来流速度100 m/s，来流迎角分别为0°～30°的5种工况。同样，为研究来流速度影响，分析了迎角30°，来流速度U∞的5种工况，如图6所示。结果表明，伞衣与机身的距离随着来流迎角和速度的上升而下降。根据伞衣与机身的相对位置，伞衣充气过程可以被分为3中典型状态：当 δ>R时，伞衣远离；当 δ≈R时，伞衣贴近机身；当 δ

图5 降落伞与机身尾部相对位置示意图Fig. 5 Schematic diagram of the relative position between a parachute and a fuselage tail

图6 动载峰值时刻伞衣与机身的距离Fig. 6 The distance between a parachute and a fuselage when the opening load reaches its maximum

以往研究表明，伞衣外形变化过程与剪切层的演化密切相关[17]。图7给出了来流速度160 m/s、迎角30°、来流速度100 m/s、迎角30°和来流速度100 m/s、迎角0°这3种工况计算结果。图7（a～c）为不考虑机身影响的工况，如图所示，伞衣底边分离点后形成了相对稳定的剪切层。而机身影响会改变流场特性，如图7（d～f）：当伞衣贴近机身或与机身发生挤压时（图7（d、e）），伞衣靠近机身位置处剪切层消失，伞衣近尾流区的压力没有出现急剧下降；当伞衣与机身距离较远时（图7（f）），伞衣靠近机身一侧又出现了剪切层结构，机身影响变得不明显。

图7 动载峰值时刻伞衣对称面的压力云图，等值线为流向速度Fig. 7 Flow fields at the parachute symmetrical plane (Shade and line contours denote pressure and streamwise velocity, respectively)

Lamb矢量散度与流场中的动量输运过程相关，可被用于研究流场演化的动力学过程[18]。Lamb矢量散度定义为 ▽·L。这里，L=ω×u为Lamb矢量；ω为伪涡矢量；u为流体速度矢量。因为伞衣充气过程主要受剪切层和涡结构运动影响，故 ▽·L可被用于描述伞衣附近流场特征[19]。

图8给出了利用Lamb矢量散度描述的流场拓扑结构，其中图8（a～c）为不考虑机身影响的工况。可以看出，在自由来流中，剪切层在伞衣底边分离后，▽·L产生了负值包裹正值的三层结构，类似结构可以减缓剪切层失稳[20]。在伞衣的近尾迹区中， ▽·L演化为正负交替结构（图8（a～c））。此类结构一般出现在钝体绕流的剪切层和近尾迹区[21]，意味着此处存在剧烈的高低速流体动量交换。从图8（d）中可以看出，伞衣挤压机身时不会明显减弱伞衣后方的动量交换。另外，从图8（e）中还可以看到，当伞衣与机身存在较小间距时，伞衣靠近机身侧与机身之间形成通道流动，使得该处的动量交换加剧。因此，我们有理由认为伞衣与机身之间间距较小时，靠近机身的伞衣后方不会出现严重低压区（见图7），即伞衣被吸在机身上的风险较低。

图8 动载峰值时刻伞衣对称面的Lamb矢量散度云图Fig. 8 Divergence of the Lamb vector at the parachute symmetrical plane

涡结构的动力学特征与流场演化密切相关，涡量场的演化被涡量动力学方程所制约[22]：

其中 (ω·▽)u表示涡线发生拉伸扭曲导致的涡量场的变化； −ω(▽·u)被称为散度项，表示流体体积膨胀（或压缩）导致的涡量变化； ρ−2(▽ρ×▽p)被称为斜压项，表示流体等密度面和等压面不重合导致的涡量变化；方程右侧第四项（黏性扩散项）和第五项（体积力项）幅值较小，涡量的影响主要集中在前三项。

图9～图11给出了伞衣流场中的涡量输运项分布。不难得知，涡量的生成和输运主要发生在伞衣的近尾迹区。在三种工况下，涡拉伸/扭曲项对涡量的生成明显占主导地位，其次是散度项的影响，斜压项对涡量生成的贡献最小，这与钝体可压缩绕流的现象类似[23]。当伞衣靠近机身时，由于通道流动的存在和机身“地面效应”的影响，使得伞衣后方靠近机身侧的涡运动更加活跃。此时，涡量生成的加剧主要来源于涡拉伸/扭曲以及流体的胀压效应。这种涡运动没有导致伞衣近尾迹区产生严重低压区。Tsutsui[24]在圆球地面效应的研究中也观测到了类似的现象。由此可以得知，当伞衣靠近机身时，伞衣后方存在明显的涡运动，不会导致伞衣被吸在机身上。

图9 伞衣对称面的涡量和涡量动力学方程中前三项的幅值(U ∞=160m/s,α=30°)Fig. 9 The magnitude of vorticity and the first three terms in the vorticity dynamics equation at the parachute symmetrical plane (U∞=160m/s,α=30°)

图11 伞衣对称面的涡量和涡量动力学方程中前三项的幅值(U ∞=100m/s,α=0°)Fig. 11 The magnitude of vorticity and the first three terms in the vorticity dynamics equation at the parachute symmetrical plane (U ∞=100m/s,α=0°)

图10 伞衣对称面的涡量和涡量动力学方程中前三项的幅值(U ∞=100m/s,α=30°)Fig. 10 The magnitude of vorticity and the first three terms in the vorticity dynamics equation at the parachute symmetrical plane (U ∞=100m/s,α=30°)

3.2 伞衣充气过程

为了分析来流迎角对伞衣充气过程的影响，图12给出了来流速度100 m/s，迎角0°和30°时，伞衣在充气过程中的动载与投影面积变化。当迎角为0°时，伞衣动载和投影面积的变化曲线非常相近。说明该工况下机身对伞衣充气过程的影响较弱，伞衣充气过程与自由来流工况较为相近。当迎角为30°时，机身下方的伞衣充气速度较慢。其原因是，在充气过程的第一阶段结束时（t= 0.06 s），机身下方降落伞的底边进气口面积较小（见图13）。由于机身影响了伞衣剪切层和尾流的演化，从图13中还可以看到，机身下方的降落伞在充气过程的第三阶段中没有发生伞顶局部塌陷。这导致了伞衣投影面积和动载峰值较小。

图12 伞衣充气过程中的动载与投影面积(U ∞=100m/s)Fig. 12 The opening load and projected area during the parachute inflation process (U ∞=100m/s)

图13 伞衣在充气过程中的外形变化(U ∞=100m/s,α=30°)Fig. 13 The canopy deformation during the parachute inflation process (U ∞=100m/s,α=30°)

3.3 降落伞离机安全分析

降落伞安全离机主要指，1）伞衣不会受到机身影响导致开伞失败；2）降落伞离机过程中，伞衣不会剐蹭机身导致破损失效；3）降落伞离机过程中，飞机机身结构不会因为伞衣伞绳钩挂或前置体其他部件撞击导致损伤。

在伞衣充气的初始时刻，降落伞有沿z轴负方向的弹射速度Ueject。因此在伞绳拉直的过程中，前置体质点会向z轴负方向运动。如图14（a）所示，当来流速度与弹射速度比值V=U∞/Ueject较小时，在充气过程中，伞衣会在机身对称面内向机身斜下方摆动并逐渐远离机身。如图14（b）所示，当来流速度增大为160 m/s时，在伞绳拉直时伞衣与机身距离较小。在充气过程中，伞衣主要沿来流方向运动更容易挤压剐蹭机身。如图14（c）所示，当弹射速度增大为20 m/s时，降落伞在充气过程中会向机身斜下方产生更大幅度的摆动，伞衣在充气过程中不会接触机身。因此可以采用增大弹射速度的方法提高降落伞离机过程的安全性。

图14 不同来流速度和弹射速度下，伞衣在充气过程中的外形变化Fig. 14 The canopy deformation during the parachute inflation process at different free-stream velocities and eject velocities

为了进一步分析来流速度和弹射速度对伞衣充气过程的影响，图15（a）中给出了来流迎角30°，来流速度为40 ～160 m/s时，充气过程中的动载变化。可以看到，在40 m/s和100 m/s的来流速度下，机身不会对动载峰值产生明显影响；在160 m/s的来流速度下，由于伞衣与机身发生挤压剐蹭，动载峰值明显降低。如图15（b）所示，当弹射速度增加到20 m/s时，由于伞衣不再接触机身，动载峰值与自由来流的工况几乎一致。

图15 伞衣充气过程中的动载幅值Fig. 15 The opening load during the parachute inflation process

现代大型客机具有飞行包线保护功能，其中迎角限制是该功能中最重要的一个部分[25]。HBG系统可将超出飞行包线的状态识别为紧急状态。为避免伞衣剐蹭机身，HBG系统需要在飞机达到临界迎角前弹射降落伞。以B737的最大失速迎角作为参考[26]，以起飞安全速度（V2= 59 m/s）为最小客机飞行速度。如图16所示，随着来流速度增加，HBG系统安全离机的临界迎角下降。飞机失速迎角与HBG系统临界迎角共同组成了HBG系统的工作包线（图16中黑色虚线区域）。随着来流速度增加，可供HBG系统工作的来流迎角范围不断减小。考虑到在相同来流工况下，降落伞弹射离机的速度越大，在伞绳拉直时降落伞与机身的距离就越远。所以为了增大HBG系统的适用范围，可以采用提高降落伞弹射速度的方法。当弹射速度从10 m/s提高到20 m/s时，临界迎角大幅提高，使HBG系统获得了更大的适用范围。