开“窗”结构对环帆伞开伞过程影响

包文龙，贾贺，2，*，薛晓鹏，黄雪姣，高树义，荣伟，王奇，吴壮志

1．北京空间机电研究所，北京 100094

2．南京航空航天大学 航空学院，南京 210016

3．中南大学 航空航天学院，长沙 410083

4．北京航空航天大学 计算机学院，北京 100191

环帆伞因开伞性能优越、抗破坏能力强、工作可靠、阻力系数大和稳定性好，被作为核心减速装置广泛应用于各种航天器的回收着陆系统［1-3］。作为一种依靠其结构形状和性能特征来改变载荷体运动状态的气动减速装置，其气动特性的好坏尤为重要。已有试验表明，适当改变伞衣的透气性可以使阻力效率和开伞可靠性，以及稳定性和开伞动载之间达到一个微妙的平衡［4］。

近年来，随着计算机技术的发展，国内外已有诸多的研究团队利用数值模拟技术对环帆伞透气性进行研究。Tezduyar 等［5］将环帆伞伞衣划分为12 个不同位置的同心圆区域，利用局部变化的织物透气性来模拟整个环帆伞的透气性，从而实现织物透气性与结构透气性统一，该方法有助于减少环帆伞这一复杂伞型数值计算所需的网 格 数 量。 Takizawa 等［6］应 用FSI（Fluid-Structure Interation）方法分析了3 种环帆伞设计构型的性能，结果表明移除靠近伞衣底部的帆能有 效 提 高 静 稳 定 性。Greathouse 和Schwing［7］基于CFD（Computational Fluid Dynamics）方法对不同结构的环帆伞开展性能分析，以研究结构透气性对环帆伞阻力特性和静稳定性的影响，并提出以最小阻力损失最大化静稳定性的设计方案。甘小娇等［8-11］探究了部分结构参数及织物透气性对环帆伞气动性能的影响，但均选择相对简单的稳降阶段。部分学者针对环帆伞充气过程进行了数值模拟，但未能研究透气性对开伞性能的影响［12-13］。

环帆伞的设计改进不可避免地涉及透气性，但当前透气性对其开伞过程的影响研究还不够深入。另外，随着中国航天技术的发展及空间站建设，返回舱的重量及体积将有所增加。而目前采用的单伞系统受工艺、成本、可靠性等因素影响已无法满足工程需求，群伞系统获得关注并已在新一代载人飞船试验船成功应用，但也带来开伞不同步和伞间干扰等问题。

为解决上述问题，国外经多次空投及风洞试验设计修改，采用了开“窗”结构设计，发现该结构有利于提高群伞系统的同步性和稳定性［14-15］，图1 为NASA 猎户座飞船主伞采用的开“窗”结构。中国目前对此结构尚未应用也无相关公开仿真研究成果。

图1 猎户座飞船主伞的开“窗”结构Fig.1 ‘Windows’ structure of NASA Orion main parachute

同时，由于群伞间复杂气动干扰的存在，对其开伞过程数值仿真难度较大。而群伞中各个组成伞结构相同，且很大程度上决定了群伞系统整体性能的优劣。故本文采用流固耦合方法，对不同开“窗”结构的环帆伞单伞在无限质量情况下的开伞过程进行数值模拟，进而分析开“窗”结构对其气动特性的影响，为后续从开“窗”角度设计优化环帆伞，从而提高群伞系统性能提供参考。

1 研究对象

本文以某环帆伞为研究对象，开展开伞过程数值计算。其平面结构及具体尺寸如图2 和表1所示。同时，为方便后续分析，图2（a）对环1～4和帆1～6 进行了序号标注。

图2 环帆伞结构Fig.2 Flat structure of ringsail canopy

表1 环帆伞结构尺寸Table 1 Structural size of ringsail parachute

由于加强带等细化结构会影响流固耦合的仿真结果［16］，故本文的环帆伞模型由环、帆、伞绳、顶孔绳、径向带、纬向带等结构构成，确保仿真模型与真实物理更为贴近。伞衣幅结构和整体绳带结构如图3、图4 所示。

图3 环帆伞伞衣幅结构Fig.3 Structure of ringsail canopy gore

图4 整体绳带结构Fig.4 Integral rope structure

环帆伞的环片和帆片为2 种不同的材料，经测得厚度分别为0.15，0.10 mm。采用Ergun 公式对2 种材料的透气性进行描述［17-18］，可得

式中：a、b分别为描述织物透气性的黏性系数和惯性系数，可经试验获得；e为伞衣厚度；vq为透过伞衣织物的气流速度；Δp为伞衣织物的试验压差。经透气量与压差曲线双点估计反算可得开伞过程2 种织物材料的黏性系数a和惯性系数b。

环片：

帆片：

通常认为在最大投影直径位置以下进行开“窗”对阻力特性影响较小，该伞空投试验中的投影直径Dt与结构直径Dj比为0.82，接近于帆5 位置。数量采用占伞衣幅数的百分比定义，根据图1可知国外已应用的开“窗”数量为20%。考虑到本文基准算例中的环帆伞伞衣幅数为24，且开“窗”数量需为整数，最终在不改变伞衣构型设计的前提下，以25%和帆5 位置为基础选取5 种不同的结构分析开“窗”位置和数量对环帆伞气动特性的影响，各研究对象结构参数如表2 所示。

表2 研究算例Table 2 Research examples

2 控制方程

采用任意拉格朗日-欧拉（ALE）方法对环帆伞开伞过程进行流固耦合仿真［19-24］。

1） 流场域控制方程

式中：ρf为流体密度；t为时间；f为体积力；E为能量；vf和vm分别是流体和网格节点的速度；σf为流 体 的 应 力 张 量，σf=-pδ+μ(grad(vf)+grad(vf)T)；p为压强；δ为Kroneckerδ函数；μ为动力黏度。

2） 结构域控制方程

3） 流固耦合实现

对上述控制方程采用罚函数法进行结构与流场的全耦合计算，主要策略是追踪结构节点与流体节点之间的相对位移d。若以dn表示时间t=tn时的相对位移，并采用中心差分法按时间递增进行求解，则按以下等式对其进行更新：

式中：vr表示结构节点与流体节点之间相对速度，；Δt为时间间隔。

当发生穿透，即在ns·dn＜0 时进行耦合，罚函数耦合力Fp分布到结构节点与流体节点上，ns表示与结构节点相连的结构单元的平均法线方向；若未发生穿透，就不进行任何操作。

罚函数耦合力的大小与穿透深度呈正比，即

式中：i为笛卡尔坐标系中的3 个坐标轴方向；ki为i方向上基于主从节点模型特性的刚度系数；di为i方向上的穿透深度；界面接触耦合力Fi的大小通常与罚函数耦合力的大小相等，但由于本文考虑织物透气性，需基于Ergun 公式对罚函数耦合力进行修正，同时伞衣采用薄壳单元，故i=xn，即仅考虑法向的耦合：

式中：an和bn为法向的黏性系数和惯性系数；vn为法向的流体速度；ε为孔隙率；S为单元表面积。结构节点的界面接触耦合力Fs=-Fxn；而对于流体，界面接触耦合力需基于流体节点j的形函数［25］Nj分配给各个节点，即Fjf=Nj Fxn。

3 仿真模型及验证

3.1 半折叠状态下环帆伞模型

开伞过程流固耦合仿真中，通常需要有一定的进气口面积才能保证伞衣的顺利充气。本文选择具有一定初始展开的伞衣外廓作为初始状态（初始展开直径与名义直径比为0.047），即建立半折叠状态下的伞衣仿真模型，继而进行后续模拟计算。本研究中半折叠状态下的伞衣模型基于数值计算方法获得，主要利用了顶孔绳结点和伞绳结点的同步位移。伞衣采用三角形单元，减小折叠过程中出现网格穿透的概率。数值计算方法相对于直接建模方法计算快，人工干预因素小。最终，经数值计算后生成的伞衣半折叠模型如图5 所示，网格数为26 040。

图5 伞衣半折叠模型Fig.5 Half-folding model of canopy

3.2 流场域模型

通常来讲，环帆伞开伞过程数值仿真中，流场边界至降落伞距离在5D0以上为宜，但考虑到大变形流固耦合仿真所需的计算资源庞大，本文流场域取的相对较小，沿降落伞轴向、法向、侧向的流场计算域尺寸为5D0×3D0×3D0。同时为更好地捕捉流场特征，在综合考虑计算精度与时间成本的基础上，对流场的网格划分采用十字形加密方式，以在伞衣附近区域建立更为细密的网格，确保降落伞伞衣附近的流体网格尺寸接近降落伞伞衣单元尺寸。另外在伞衣加密部分沿流场两边轴向方向的网格采用过渡方式，即网格逐渐稀疏，以避免出现质量较差的网格。最终流场域的六面体网格数为276 万，建立的流场计算域及仿真模型整体网格结构如图6 所示。

图6 仿真模型网格结构Fig.6 Mesh structure of simulation model

仿真中，由于该环帆伞开伞时的来流马赫数为0.216，因此整个流场仿真的速度区间为低速段（可以不考虑空气压缩性影响）。同时，由于伞衣面积较小，开伞过程通常在零点几秒内即已完成，空气组分、密度、压强等参数的时间变化相对平缓，因此可采用恒定大气密度和压强进行仿真。流场顶端和壁面均添加无反射约束，用以模拟远场状态，使仿真结果与环帆伞置于无限大流场的情形更为接近。

3.3 模型验证

计算中开伞速度为70 m/s，开伞高度为4 000 m，初始攻角为0°，仿真时间为3 s，研究对象为基准算例。以伞衣外形、阻力系数及开伞动载系数作为考察对象，比较仿真结果与理论及实际试验之间的异同。

数值模拟得到的伞衣外形如图7 所示。由图7 可以发现：来流经伞衣底边进气口进入后，形成1 个由伞衣包裹的管状区域（0.1 s），之后伞衣顶部的空气开始聚集，张开程度逐渐超过了伞衣底部（0.2 s），但随后受到结构惯性以及张力阻碍减缓了扩张（0.4 s），同时张力使得进气口进一步扩大，引起伞衣的迅速张满（0.6 s）。伞衣外形的变化与降落伞充气理论描述基本吻合。图8 为仿真计算与空投试验获得的环帆伞张满外形对比，由图8 可知，二者基本一致。

图7 开伞过程中的伞衣外形变化Fig.7 Canopy shape changes during inflation

图8 仿真计算结果和空投试验张满外形对比Fig.8 Comparison of canopy shape between numerical simulation and airdrop test

图9为伞衣阻力系数仿真时间历程曲线。以1.5 s 后阻力系数的平均值作为平均阻力系数。由图9 可知：开伞过程中伞衣最大阻力面积为49.371 m2，开伞结束后该环帆伞的平均阻力面积为36.075 m2，对应的平均阻力系数为0.717 2，与空投试验值0.82 较为接近，偏差为12.5%。考虑到影响伞衣阻力系数的因素众多，包括伞衣折叠模型、网格划分质量、空投试验误差、仿真算法适用性以及环境参数真实程度等，认为以上偏差可以接受。另外，根据图9 可得，开伞动载系数kd为1.368，与环帆伞动载系数≥1.1 相符［26］。综上认为，本文建立的仿真模型可以有效模拟环帆伞的开伞过程。

图9 伞衣阻力系数仿真时间历程曲线Fig.9 Time history curve of drag coefficient

4 开“窗”结构影响

4.1 伞衣外形

图10 为各算例开伞过程中的投影面积与名义面积之比随时间变化曲线，同时仿真计算中t= 0.1，0.2，0.4，0.6 s时伞衣外形如图11 所示。就仿真结果整体而言：开“窗”后的环帆伞伞衣投影面积与名义面积比与基准算例相差最大不超过15%；采用不同开“窗”位置及数量的环帆伞开伞状态良好，未出现非对称充气、无法充满等现象。这是由于上述伞衣模型均采用对称开“窗”方式，同时伞衣主要承力部件为加强带，故未对伞衣外形在开伞过程中的变化造成较大影响。

图11 不同开“窗”方式下开伞过程中的伞衣外形变化Fig.11 Canopy shape changes during inflation process of research examples

4.2 气动特性参数

图12为各算例的气动特性参数。根据图12，可以发现开“窗”位置及数量均对该环帆伞的气动特性产生一定影响。其中，图12（d）阻力系数标准差表示各算例的平均阻力系数在1.5～3.0 s 内的波动程度。

图12 各算例的气动特性参数Fig.12 Aerodynamic characteristic parameters of research examples

相比于基准算例，数量不变，不同开“窗”位置变化如下：帆4 和帆5 位置的平均阻力系数分别减小1.56%、2.60%，帆3 位置平均阻力系数增加1.40%；帆 3、帆4、帆5 位置的开伞动载系数依次减小2.5%、6.8%、4.5%；帆4 位置的阻力系数标准差增大27.9%，而帆3 和帆5 位置阻力系数标准差分别减小14.5%、21.1%。

相比于基准算例，位置不变，不同开“窗”数量变化如下：25%、50%、75%数量的平均阻力系数依次减小2.6%、11.6%、12.8%；开伞动载系数依次减小4.5%、4.8%、9.1%；阻力系数标准差依次减小21.1%、41.4%、36.2%。

综上，选择开“窗”位置时，相比于帆3 和帆4，在帆5 位置进行开“窗”可在平均阻力系数略有下降，开伞动载系数降低的同时，最大程度减小阻力系数的波动；选择开“窗”数量时，相比于50%和75%，以25%的数量进行开“窗”可在开伞动载系数降低、阻力系数波动减小的同时，最大程度减小对平均阻力系数的影响。

4.3 摆动角度

图13 为摆动角度仿真时间历程曲线。根据图13 可以看出：总体而言，开“窗”结构会显著减小 最大摆角，减小幅度在5°～10°，其中帆5 开“窗”50%和帆5 开“窗”75%时的减小幅度最大；开“窗”数量不变时，帆3 和帆5 位置的摆动角度最小；开“窗”位置不变时，50%和75%数量的摆动角度最小且基本相同。

受计算资源以及求解能力限制，流固耦合计算对流场中旋涡、压力的捕捉精度有限。为分析上述曲线最大摆动角度的变化原因，基于计算流体力学方法获得如图14、图15 所示摆角20°时的流场速度分布云图及压力云图。

图14 摆角20°时流场速度分布云图Fig.14 Nephogram of flow field velocity distribution at 20°

图15 摆角20°时，压力分布云图Fig.15 Nephogram of pressure distribution at 20°

根据图14、图15 可以看出：各算例的伞衣内外部流场相似，伞衣外部左右两侧均存在一个范围较大的高速区域（低压区）。开“窗”以后，伞衣内部的部分气流从“窗口”流出，一方面，在伞衣左侧与外部气流汇合，阻扰外部高速区域生成并使其远离伞衣，致使从环缝、月牙缝流出的射流偏折角度减小，导致“窗口”上侧伞衣幅外部压力降低；另一方面，在伞衣右侧，破坏了沿顶孔至伞衣底边内部流场的连续性，造成伞衣右侧内部整体压力降低；二者共同作用下，降落伞回正作用增强，最大摆角减小；帆5 较帆4 最大摆动角度较小的原因可能是越靠近伞衣底部，伞衣左侧“窗口”上侧伞衣幅外部压力降低的区域面积越大；而帆3 较帆5 最大摆动角度较小的原因可能是越靠近伞衣顶部，伞衣右侧沿顶孔至伞衣底边内部流场的连续性越差，致使右侧内部整体压力下降越多。

开“窗”数量越多，伞衣内外流场受上述作用影响越大，这可能是造成最大摆角大幅度减小的原因。同时50%和75%开“窗”数量时的伞衣内外部流场基本相同，表明50%以上开“窗”数量已对伞衣随机摆动过程中的流场变化已无明显影响，最大摆角因开“窗”数量增加而减小的程度已达极限。

5 结 论

开“窗”结构有助于提高群伞系统的同步性和稳定性。本文基于流固耦合方法针对单伞开展不同开“窗”结构的环帆伞开伞过程数值计算，探究了开“窗”位置和数量对环帆伞伞衣外形、气动特性参数、摆动角度的影响，所获结论如下：

1） 由于为对称开“窗”，且伞衣主要承力部件为加强带，故不同的开“窗”位置和数量对开伞过程中伞衣外形变化影响较小，未出现非对称充气、无法充满等现象。

2） 开“窗”数量相同时，在帆5 位置开“窗”可在平均阻力系数基本不变的情况下，降低开伞动载系数，最大程度地减小阻力系数的波动。

3） 开“窗”位置相同时，以25%数量开“窗”可在降低开伞动载系数、减小阻力系数波动的同时，最大程度地减小对平均阻力系数的影响。

4） 开“窗”以后，伞衣内部气流从“窗口”流出，一方面在伞衣左侧与外部气流汇合后，阻扰外部高速区域生成并使其远离伞衣，致使从环缝、月牙缝流出的射流偏折角度减小，导致“窗口”上侧伞衣幅外部压力降低；另一方面在伞衣右侧，破坏了沿顶孔至伞衣底边内部流场的连续性，造成伞衣右侧内部整体压力降低。二者共同作用下，环帆伞的最大摆角发生显著减小，且减小幅度随着开“窗”数量的增加而增大，但与开“窗”位置之间并无明显规律。