降落伞具有体积小、质量轻、可靠性高等优点，是目前航天器进入其他具有大气层的天体或返回地球过程的主要气动减速装置之一^[1]。降落伞工作过程一般包括拉直、充气和稳降等3个阶段^[2]，整个过程一般处于飞行器的尾流环境。其中，降落伞的弹射拉直过程一般通过火工作动装置实现，该过程为后续降落伞充气张满创造条件。因此降落伞的弹射拉直对于伞系统工作的可靠性具有重要意义。

降落伞的弹射拉直过程研究一般可以通过试验和仿真技术实现。其中，试验可以分为空投试验和地面试验，空投试验可以模拟真实的飞行环境，但是周期长、成本高，而且可验证工况有限；地面试验很难模拟真实的飞行环境，尤其是飞行器的尾流环境，而且周期也相对较长。因此，目前国内外对降落伞的弹射拉直研究较多集中于仿真模拟。

国外针对降落伞的弹射拉直过程的研究起步较早，20世纪70年代，Mcvey等^[3]基于拉直过程中伞绳始终保持直线状态的假设提出了简单连续拉直模型，采用动力学方法用于研究降落伞的弹射拉直过程。随后动力学模型方法仿真被广泛应用于降落伞弹射拉直仿真模拟^[4-8]。动力学模型方法在一定程度取得了成功，但是也有一定自身局限：已有研究基本没有考虑飞行器尾流对弹射拉直过程的影响^[9]。王海涛等^[9]采用计算流体力学(computational fluid dynamics，CFD)将弹射拉直时刻计算的航天器尾流区速度场叠加于弹射分离伞包空速上，通过动力学仿真研究尾流对弹射拉直的影响。该方法在亚声速0°攻角工况下较为适用，但在超声速非0°攻角工况下尚显不足：一方面，超声速工况下尾流与伞包等会产生尾流及尾流-激波干扰等复杂的流动现象；另一方面，非0°攻角会产生气动力矩，从而使分离体姿态发生变化，上述现象均会降低降落伞的弹射拉直过程仿真精度。美国国家航天局(National Aeronautics and Space Administration，NASA) 的设计人员提出，攻角是一个关键参数，其影响降落伞正常拉直和充气、着陆器的轴向过载，以及开伞后舱伞系统的动力学特性。火星探测器“海盗号”为保证降落伞顺利弹射拉直及充气，要求进入器进入攻角不超过13°^[10]；新一代火星探测器“火星科学实验室”更是以不惜抛掉6个质量块调节配平攻角的方式保证降落伞安全工作^[11]。因此，攻角效应对降落伞弹射拉直过程具有重要研究意义。

降落伞的弹射拉直过程实质为变质量多体分离问题，物理过程十分复杂，计算难度较大。当前处理多体分离问题的技术中嵌套了网格，因使用简单、技术成熟，被广泛应用于航空航天领域^[12]。

本文基于嵌套网格技术，采用CFD技术耦合六自由度(six degree of freedom，6DoF)运动方程的方法研究了攻角效应对超声速降落伞弹射拉直过程的影响。

1 数值方法

本文采用非定常Reynolds平均N-S(unsteady Reynolds averaged Navier-Stokes, URANS)方程耦合6DoF运动方程的方法对超声速降落伞的弹射拉直动态过程进行数值仿真。

1.1 流体控制方程

采用URANS方程，其守恒形式在直角坐标系中可表示为^[13]

$ \frac{\partial \boldsymbol{Q}}{\partial t}+\frac{\partial \boldsymbol{F}}{\partial x}+\frac{\partial \boldsymbol{G}}{\partial y}+\frac{\partial \boldsymbol{H}}{\partial z}=\frac{\partial \boldsymbol{F}_{\mathrm{v}}}{\partial x}+\frac{\partial \boldsymbol{G}_{\mathrm{v}}}{\partial y}+\frac{\partial \boldsymbol{H}_{\mathrm{v}}}{\partial z} $

(1)

其中：Q为守恒变量向量；F、G、H为3个方向的无黏通量；x、y、z为直角坐标系中的3个坐标；F_v、G_v、H_v为3个方向的黏性通量。

1.2 湍流模式

SST(shear stress transport)湍流模式^[13]是一个依赖于切换函数F₁的分区模型：在靠近壁面的边界层内为k-ω模型，而在远离壁面的主流区则切换为k-ε模型。它克服了k-ω模型对来流湍流度过于敏感以及k-ε模型在前缘滞止点湍动能过大的问题，计算鲁棒性好，在工程中应用广泛。

SST湍流模式的控制方程^[14]可表示为

$ \begin{aligned} \frac{\partial(\rho k)}{\partial t}+\frac{\partial\left(\rho u_i k\right)}{\partial x_i}= & \frac{1}{R e_{\infty}} \frac{\partial}{\partial x_i}\left[\left(\mu+\sigma_k \mu_t\right) \frac{\partial k}{\partial x_i}\right]+P_k- \\ D_k \frac{\partial(\rho \omega)}{\partial t}+\frac{\partial\left(\rho u_i \omega\right)}{\partial x_i} & =\frac{1}{R e_{\infty}} \frac{\partial}{\partial x_i}\left[\left(\mu+\sigma_\omega \mu_t\right) \frac{\partial \omega}{\partial x_i}\right]+P_\omega- \\ D_\omega & +\left(1-F_1\right) \mathrm{CD}_{k \omega} . \end{aligned} $

(2)

其中：ρ为气体密度；t为时间；u_i为速度分量；x_i为坐标分量；Re_∞为远场单位Reynolds数；μ为分子黏性系数；σ_k与σ_ω为常系数；μ_turb为涡黏系数；P_k为湍动能k的生成项；D_k为湍动能k的破坏项；P_ω为耗散率生成项，D_ω为耗散率破坏项；CD_kω为交叉系数。

1.3 6DoF运动计算

根据刚体动力学，刚体运动可分解为平动和转动，对应的6DoF控制方程^[15]可表示为

$ \begin{gathered} \boldsymbol{F}=m \frac{\mathrm{d} \boldsymbol{v}}{\mathrm{d} t}, \end{gathered} $

(3)

$ \boldsymbol{M}=\frac{\partial \boldsymbol{h}}{\partial t}+\boldsymbol{\varOmega} \cdot \boldsymbol{h} . $

(4)

其中：v为速度，h为角动量，M为刚体受到的力矩，Ω为角速度。

该仿真模拟基本计算过程如图 1所示。

1.4 其他数值方法

本模拟以CFD++商业软件为仿真平台，使用有限体积求解器；空间离散格式为总变差不增(total variation diminishing，TVD)格式；湍流模式采用SST湍流模型；时间推进格式采用隐式LU-SGS(lower-upper symmetric Gauss-Seidel)方法。

2 数值方法标模验证

三维机翼下外挂物投放是验证6-DoF运动物体动态分离的标准验证算例，具有详细的实验设计参数和风洞实验测量数据，一般采用该模型验证动态分离问题的正确性，因此，本研究也选用该模型作为验证算例。

该模型由机翼、挂架及带翼弹体共3部分组成，弹体与挂架之间存在3.66 cm的缝隙，具体几何参数可参见文[16]。

2.1 仿真工况

分离参数：分离时刻来流Ma为1.2，攻角与侧滑角均为0，分离高度H为11.6 km，来流压强p为20 659 Pa，来流密度ρ为0.332 16 kg/m³。在分离指令下达0.054 s内，在弹体上施加2个数值向下的弹射力以加速弹体的分离。弹体的质量特性参数及弹射力相关作用参数参见文[17]，弹射力的位置如图 2所示，其中F_FE为前弹射力，F_AE为后弹射力，M_x为前弹射力产生的相对于质心的力矩，M_n为后弹射力产生的相对于质心的力矩。弹体在气动力、弹射力及重力的共同作用下向下运动。同时，由于弹体与机翼、挂架等存在相互干扰，弹体姿态会在下落过程中发生较为明显的变化。

2.2 仿真网格

该仿真网格采用嵌套网格，网格总数约600万，如图 3所示。图 3a为弹体对称面位置网格，在弹体的运动轨迹区域进行了网格加密处理，弹体网格通过嵌套边界与机翼背景网格传递流场信息；图 3b为图 3a中弹体周围网格(红圈位置处)局部放大。

2.3 仿真结果对比分析

本算例仿真模拟时间设定弹射初始时刻t=0.0 s，总计算时间0.6 s，图 4给出了不同分离时刻机翼与弹体的相对位置，从图中可以看出，随着时间推进，弹体在远离机翼且姿态发生明显变化。

为进一步分析弹体与机翼分离过程，图 5给出了弹体在分离过程中速度、位移及Euler角随时间变化的曲线，并与文[17]实验和文[18]仿真结果进行了对比，其中CFD++为本研究仿真结果，Fluent为文[18]计算结果。可以看出，本文仿真结果中的速度和位移随时间变化的曲线与实验测量数据吻合；俯仰角和偏航角变化也吻合较好，但滚转角在t=0.2 s后，本文和Fluent结果均与实验出现了较大偏差，相比之下，本文计算结果更趋近实验值。出现这种偏差主要有2个原因^[17-18]：1) 仿真模型中无黏假设会导致对阻力的低估，实验过程中准静态特性与实际动态过程也存在差别；2) 相比其余2个轴，弹体绕中心轴的转动惯量I_xx相对较小，滚转角对气动力误差较为敏感。

总体来看，本文数值方法可用于物体动态分离仿真模拟。

3 仿真模拟的建立 3.1 仿真模型与工况

本研究对象为返回舱在超声速状态下与伞包-舱盖组合体(简称“分离体”)弹射分离过程。返回舱与分离体的模型如图 6和7所示，其中D为返回舱最大直径。

仿真工况如表 1所示，3个攻角工况均为超音速来流工况，其余参数完全相同。

3.2 计算网格

降落伞弹射拉直过程中，伞绳处于松弛状态且质量较轻，因此该仿真模拟中暂未考虑伞绳的影响。本文仿真分别单独生成返回舱网格、分离体网格以及背景网格，3套网格通过嵌套边界生成计算网格，网格规模约1 400万，如图 8所示，在物体的运动轨迹区域进行了网格加密，保证了计算精度。

3.3 初始边界条件

1) 边界条件。

返回舱网格外边界为嵌套边界，内边界为物面边界；分离体网格外边界为嵌套边界，内边界为物面边界；背景网格为远场边界。

2) 初始条件。

本仿真涉及来流速度与弹射分离速度。来流速度反映分离体分离前的飞行速度；分离速度体现分离时分离体与返回舱的相对速度，其方向沿返回舱底部外法向方向。本仿真中初始速度通过经验公式^[8]计算如下：

$ v_{\mathrm{r}}=\sqrt{2 n g L} \text {. } $

(5)

其中：v_r为分离体相对于返回舱的分离速度，n为过载系数，g为重力加速度，L为有效尾流长度。由式(4)可得v_r=21 m/s。

4 仿真结果分析与讨论

本文仿真分为2个阶段，第1阶段为静态流场求解；第2阶段为弹射拉直的动态流场求解，该阶段通过在静态流场基础上续算得到。

4.1 时间步长

为提高仿真结果的可靠性，本研究首先开展时间步长对仿真结果的影响。时间步长根据典型工况圆柱的St=0.2进行估算，Δt=0.000 5 s，因此，分别选取Δt为0.001、0.000 5和0.000 25 s的情形进行研究。

$ S t=f l / U. $

(6)

其中：f为涡脱频率，l为特征长度，U为流体速度。

采用不同时间步长获得的返回舱与分离体受到的气动力随时间变化的曲线如图 9和10所示。由曲线可以看出，Δt=0.001 s时得到的返回舱与分离体x向气动力比Δt=0.000 5 s与Δt= 0.000 25 s的仿真结果略大，y和z向气动力基本一致；而Δt=0.000 5 s与Δt=0.000 25 s仿真结果基本一致，因此本文选用Δt=0.000 5 s作为时间步长进行仿真模拟。

4.2 静态流场

图 11给出了3种攻角工况下对称面Ma云图分布，呈现如下流动特征：由于来流为超声速流动，在舱体头部会形成一道明显的弓形激波，波后Ma明显降低；流动在经过返回舱肩部时，形成膨胀波，Ma升高；流动在返回舱后体开始出现流动分离；在返回舱底部形成大范围的回流区域，回流区内流动与分离体相互作用，流动更为复杂；回流区下游形成较长距离的低能尾迹区。

除此之外，3种攻角工况也表现出一些不同的流动特征：由于来流攻角不同，返回舱头部激波形态、分离点位置、回流区形态以及尾迹等流动结构明显不同，即所谓的“攻角效应”。

4.3 弹射拉直动态过程流场

1) 不同时刻流场对比分析。

降落伞弹射拉直过程是一个动态过程，图 12给出了上述3种工况在不同时刻对称面Ma云图。t=0.0 s时刻与上述静态流场类似，这里不再赘述。t=0.15 s时对称面Ma云图如图 13所示，该时刻分离体已经完全离开返回舱后方回流区。

(1) 在0°攻角工况下，分离体的弹射方向与返回舱尾流方向一致，该时刻分离体仍完全处于返回舱尾迹区域，形成了明显的尾流-激波干扰现象。

(2) 相比于0°攻角工况，10°攻角工况下返回舱尾流发生了明显的偏转，该时刻分离体部分已经离开尾流区，分离体在力矩作用下产生明显的姿态变化。

(3) 20°攻角工况下，返回舱后方尾流进一步偏转，此时分离体完全处于返回舱尾流区外，但是分离体产生的激波仍与尾流发生干扰现象，分离体的姿态也发生了较大变化。

图 14为3种攻角工况下降落伞完成弹射拉直过程(返回舱与分离体之间距离达到10 D)时对称面Ma云图，从图中可以看出：0°攻角工况下，分离体仍完全处于返回舱尾流中，且姿态变化较小；10°工况下，分离体在0.28 s时刻分离体再次进入返回舱的尾流区域；20°攻角工况下分离体基本离开了尾流区域。

2) 物体运动轨迹与姿态分析。

为进一步定量研究弹射拉直过程中物体运动轨迹与姿态，本文给出了返回舱与分离体受到的气动力、质心位移、力矩以及Euler角等随时间的变化曲线。

返回舱所受气动力随时间的变化曲线如图 15所示，气动力x向分量在0°攻角工况下始终最大，而气动力y向分量则在20°攻角工况时最大。返回舱质心位移如图 16所示，上述3种工况质心位移主要为正x向位移，表明返回舱在气动力与重力合力作用下作减速运动，其余2个方向质心位移基本可忽略。

返回舱在弹射拉直过程中还将受到气动力矩的作用，从而引起返回舱姿态变化。返回舱所受力矩随时间的变化曲线如图 17所示，返回舱在3种工况下受到的x和y向力矩均较小，由于攻角效应，除0°工况外，z向均受到明显的力矩作用，在z向力矩作用下，返回舱的俯仰姿态也随之发生变化，如图 18所示，返回舱姿态的改变又会影响尾流形态。

分离体在火工装置的作用下与返回舱产生相对速度，其弹射拉直过程中受到的气动力及位移随时间的变化曲线如图 19和20所示，从图中可以看出：

(1) 0°攻角工况下，由于返回舱底部回流区的作用，t≤0.02 s时刻分离体在x向受到的气动力方向与来流方向相反，随着分离体远离返回舱，x向气动力反向，且随着时间推移气动力不断升高，相比之下，y与z向气动力相对较小。

(2) 对于10°和20°攻角工况，由于回流区作用而产生的与来流速度方向相反的气动力已经不明显，由于攻角效应，分离体受到的气动力在x向和y向分量不断变化。

由于攻角效应，除0°攻角工况外，返回舱下游尾流均表现出非对称流动结构，非对称的流动作用于分离体，从而产生气动力矩，最终导致分离体的姿态变化。图 21和22分别给出了分离体气动力矩与Euler角随时间的变化曲线，0°攻角工况下，分离体受到的气动力矩较小；10°和20°攻角工况下，由于在该过程受到非对称尾流、尾流-激波干扰等复杂的流动与分离体之间相互作用，分离体在z向受到明显的气动力矩作用，该力矩呈现一定周期特性，在该力矩作用下分离体发生明显的俯仰姿态变化。

图 23给出了返回舱与分离体之间距离随着时间的变化曲线，其中，纵轴为基于返回舱最大直径的无量纲距离。可以看出，随着攻角的增加，降落伞弹射拉直时间减少。攻角效应造成尾流方向与初始分离速度方向不一致。在0°攻角工况下，分离体在整个分离过程中处于低能的尾流环境中，因此弹射分离时间相对较长；而对于10°和20°工况而言，分离体在分离过程中部分时间与含能更高的来流相互作用，从而缩短了分离时间。

5 结论

本文基于嵌套网格技术，采用CFD耦合六自由度运动方程的方法研究了攻角效应对超声速降落伞弹射拉直过程的影响，通过仿真结果的对比分析得出如下结论：

1) 对于0°、10°及20°攻角3种工况，由于攻角效应的影响，分离体在与返回舱相对运动分离的过程中，分离体与尾流的相对位置、尾流-激波干扰形态出现了明显差异。

2) 0°攻角工况下，分离体在回流区会受到与来流反向的气动力的作用，随着攻角的增加，分离体在回流区内受到的反向气动力作用逐渐减弱。

3) 对于10°与20°攻角工况，由于分离体初始分离速度方向与尾流方向不一致，从而导致分离体受到明显的气动力矩作用，该力矩使分离体姿态产生周期性变化。

4) 攻角效应导致分离体尾流与分离体的相对位置以及分离体姿态变化，进而影响了整个弹射拉直时间，随着开伞时飞行器攻角的增加，降落伞弹射拉直过程的时间减少。