随着航天事业的不断发展，天地之间的物资和人员往来也越来越密切。由于传统刚性返回器体积大、有效载荷低，迫切需要发展新型充气式返回技术用于天地往返^[1-2]。空间充气式返回器集成了回收着陆过程中的气动减速与热防护等功能，能够有效降低返回成本并应用于多种返回任务中^[3-4]。然而，空间充气式返回器在返回过程中将承受巨大的气流冲击，产生严重的动力响应并影响结构安全^[5-7]。这一现象在超声速段尤为显著，这是由于气流的黏性分离及涡运动效应显著增大，形成了不同于高超声速流动中的后体流动特性及物面压力分布，极大地影响了充气式返回器的气动载荷与振动特性^[8]。如日本宇宙航空研究开发机构(Japan Aerospace Exploration Agency，JAXA)于2009年开展了单充气环薄膜型减速器的飞行试验，起初自由飞行稳定，但是当飞行至低空时充气环内压力仅为5.2 kPa，在非定常气动力的作用下出现塌陷失稳^[9-10]。由于空间再入柔性充气结构在气动力作用下会产生大变形，并反过来影响气动载荷的分布，因此这是一个复杂的气动弹性系统^[11-12]。

传统的流固耦合问题大多针对常规弹性体，如机翼、舵面及叶轮机械等，对于柔性充气结构则较少涉及^[13-14]。与常规弹性体不同，空间充气式返回器具有独特的柔性材料属性，导致其结构动力学呈现强烈的几何及材料非线性特点，这给结构动力学与超声速流场之间的耦合增加了难度。在过去的研究中，空间充气式返回器的动力响应大多通过单向耦合进行研究，如许家裕^[15]在获得稳态气动力后研究了充气式再入飞行器(inflatable reentry vehicle experiment，IRVE)的动力响应并校核了结构应力；Lin等^[16]在考虑材料非线性特征后对空间再入充气系统的瞬态结构动力学进行了研究，分析了给定气动力激励下的瞬态响应情况及应力分布。然而，现有研究方法仍存在较多缺陷，无法合理描述空间充气式返回器在超声速工况下的气动弹性特征。首先，众多试验与计算数据表明：结构变形对流场的影响改变较大，因此单向耦合的研究方法会产生较大的偏差，已不能完全适用^[17-18]。其次，现有研究对内充压气体如何参与流固耦合进程中尚不明晰，并且柔性充气结构动力学模型也缺乏试验的有效验证^[19]。最后，飞行器在超声速阶段所受气动力较大，姿态运动特征显著，仅考虑三自由度质点的运动已不适用于飞行动力学建模与求解。

因此，本文建立了考虑内充压气体作用的流固双向耦合模型，能够揭示更符合实际的物理过程；同时，采用地面模态振动试验来验证并修正结构动力学方法，并采用大涡模拟来刻画流场细节。基于此模型，研究了空间充气式返回器在超声速流场中的气动弹性动力响应问题，合理揭示了外部流场与内充压气体共同作用下的结构动力学机理。

1 计算方法及数值模型 1.1 计算流程

计算流程分为2部分，如图 1所示。一部分为超声速段的轨迹修正。由于超声速段结构姿态运动显著，因此需要对超声速阶段进行六自由度飞行动力学修正。在六自由度飞行动力学中，考虑了绕3根轴的转动，因此将其作为计算流体动力学(computational fluid dynamics，CFD)的动边界条件时，也需引入攻角和侧滑角等参数进行数值模拟。计算选取Ma为0.8~5.0，这一过程将不断重复瞬态CFD计算以及气动力系数的反馈过程，直到2次气动力系数的残差收敛到指定范围，则认为超声速区域的六自由度飞行动力学计算完成。另一部分为典型工况下的气动弹性动力响应特征研究。针对飞行动力学中的典型工况进行更细致的流固双向耦合研究，以观察更细微的气动弹性动力响应现象及发生机理。

1.2 流固双向耦合的实现

研究方法采用时域推进的流固耦合，计算流程如图 2所示。实际计算过程是一个3层嵌套循环结构，在同一时间步内，流体域和固体域先后分别进行流体控制方程和结构控制方程的迭代求解，当达到各自的收敛目标后，跳出分域求解循环，并通过流固耦合交界面将位移和压力数据进行流体域和固体域间的传递，并进入下一轮的分域求解循环。进行若干次数据传递，当传递的节点位移和节点压力的残差达到收敛目标后，则跳出耦合求解循环，完成该时间步的流固耦合计算，并推进至下一时间步，进行下一时刻的流固耦合迭代计算。

在流固耦合计算中，气流对于叶片结构的影响通过叶片表面非定常气动力F(t)来体现，叶片固体域的结构方程表示如下：

$ \boldsymbol{M} \ddot{\boldsymbol{x}}+\boldsymbol{C} \dot{\boldsymbol{x}}+\boldsymbol{K} \boldsymbol{x}=\boldsymbol{F}(t). $

(1)

其中：M为质量矩阵，由结构质量矩阵和气体附加质量矩阵构成；C为Rayleigh矩阵，由机械Rayleigh矩阵和气动Rayleigh矩阵构成；K为刚度矩阵，由结构刚度矩阵和气体附加刚度矩阵构成；x为固体网格各单元点的瞬时变形向量；t为时间。

每一次数据传递后，流体域和固体域就要根据节点位移和速度信息进行局部变形，流固耦合交界面的形状也随之进行改变，这一过程通过动网格技术实现。动网格技术在每完成一次耦合求解循环后都会重新求解并生成与新的流体/固体域几何相适应的网格，避免因交界面的变形而挤压网格甚至产生负网格，导致计算崩溃。

1.3 结构有限元模型及验证

本文研究对象参考美国IRVE-3飞行器^[20]，该飞行器包含蒙皮、气囊及隔层3部分柔性结构，如图 3所示。其中：D_A为最大直径；D_T为气囊直径；α为半锥角；P为内充压；q为气动压力；D_C为中心体直径；H为隔层高度；W为隔层间距，气囊内部独立充气以提供刚度，隔层可有效提高结构强度，蒙皮用以承受超声速气流的冲击。柔性薄膜选择KEVLAR纤维膜，密度为1 440 kg/m³，弹性模量为131 GPa，Poisson比为0.35。图 4展示了空间再入充气结构的结构有限元模型，采用Shell 181四节点面单元，网格数10万。Shell 181单元是一种只能承受拉应力的非线性单元，其刚度仅由内应力产生，符合柔性纤维膜的特性。在结构动力学模型中，充分考虑内充压气体的刚度效应，通过加载压力载荷实施。

为验证结构动力学方法，本文开展了基于锤击法的模态振动试验，试验对象为简化的3层气囊结构，如图 5所示。模态方法验证结果如图 6所示，结果表明：尽管充气压力增加，但是第1阶(摇摆)和2阶(横向伸缩)频率变化不大，数值与试验值的偏差约为9%。第3阶(纵向伸缩)频率随充气压力呈现一定的非线性变化，数值与试验值的偏差约为7%，从而验证了充气式再入飞行器具有的材料非线性特征。考虑到销、充气阀等连接装置会产生附加质量，从而降低固有频率，因此该偏差在合理范围内。

1.4 CFD模型及验证

采用基于有限体积法的大涡模拟(large eddy simulation，LES)进行流场CFD模拟，其中空间离散采用迎风格式，时间离散采用二阶Euler格式，以保证时间方向上的守恒。此外，在求解气动热时需考虑边界层内高焓气流与结构之间的对流换热，以及壁面产生的辐射。设σ为黑体辐射系数，结构发射率ε为0.89，热辐射方程可表示如下：

$ q=\varepsilon \sigma T_{\mathrm{w}}^4. $

(2)

其中：q为热流密度，T_w为壁面温度。流体计算域如图 7所示，采用结构与非结构混合网格进行划分，总网格数为1 000万。本文选取计算域为直径30 000 mm、高42 000 mm的圆柱形，圆柱形中有一个直径3 750 mm、高10 500 mm的圆柱形区域，该区域的建立是为了进行近壁加密。结构表面第1层网格高度为0.8 mm，以保证无量纲壁面距离y⁺的值在1左右。利用文[21]的结果对CFD方法进行验证，如图 8所示。驻点动压的最大偏差发生在峰值处，此时文献和本文计算值分别为1 144和1 077 Pa，偏差约为5.8%，在合理误差范围内。

1.5 六自由度飞行动力学建模

1) 飞行动力学方程。

六自由度飞行动力学方程包含质心平动和绕质心转动2部分，分别推导如下：

(1) 飞行器质心移动的动力学方程。根据动量定理，质心平动可表示为

$ m \frac{\mathrm{d} \boldsymbol{V}}{\mathrm{d} t}=\boldsymbol{F}. $

(3)

其中：m为质量，V为速度矢量，t为时间，F为合力矢量。在动坐标系中的质心动力学方程可表示为

$ m\left(\frac{\delta \boldsymbol{V}}{\delta t}+\boldsymbol{\omega} \times \boldsymbol{V}\right)=\boldsymbol{F}. $

(4)

其中ω为角速度矢量。

(2) 飞行器绕质心转动的动力学方程。根据动量矩定理，飞行器绕质心的转动可表示为

$ \frac{\mathrm{d} \boldsymbol{h}}{\mathrm{d} t}=\boldsymbol{M} . $

(5)

其中：h为飞行器对原点的动量矩，M为外力对原点的合力矩。

(3) 飞行器质心的运动学方程。为得到质心平动的运动学方程，将飞行速度在地面坐标系上投影，由于地轴系速度(V_x, V_y, V_z)是质心空间坐标(x_g, y_g, z_g)的微分，如$ \frac{\mathrm{d} x_{\mathrm{g}}}{\mathrm{d} t}=V_x $等，可得

$ \left\{\begin{array}{l} \frac{\mathrm{d} x_{\mathrm{g}}}{\mathrm{d} t}=V \cos \gamma \cos \chi, \\ \frac{\mathrm{d} y_{\mathrm{g}}}{\mathrm{d} t}=V \cos \gamma \sin \chi , \\ \frac{\mathrm{d} z_{\mathrm{g}}}{\mathrm{d} t}=-V \sin \gamma . \end{array}\right. $

(6)

其中：γ是速度倾角，χ是速度方向角。

(4) 飞行器绕质心转动的运动学方程。通过Euler角(ψ, θ, φ)来描述飞行器在飞行过程中的姿态变化，其中：ψ角沿z_g轴旋转形成，θ角沿y_g轴旋转形成，φ角沿x_g轴旋转形成。因此，可将旋转角速度在坐标系上的投影表示为

$ \left\{\begin{array}{l} \frac{\mathrm{d} \varphi}{\mathrm{d} t}=p+\tan \theta(q \sin \varphi+r \cos \varphi), \\ \frac{\mathrm{d} \theta}{\mathrm{d} t}=q \cos \varphi-r \sin \varphi, \\ \frac{\mathrm{d} \psi}{\mathrm{d} t}=\frac{1}{\cos \theta}(q \sin \varphi+r \cos \varphi) . \end{array}\right. $

(7)

其中p、q、r分别为角速度在3个坐标方向的投影。

2) 航天器模型及初始飞行参数。

空间再入充气系统结构相关参数如表 1所示。根据美国宇航局发射IRVE-3号飞行器的任务要求，其需要携带一定质量的载荷返回地面。由于结构具有轴对称的特点，因此3个坐标轴的惯性积I_xy=I_yz=I_zx=0。

1.6 边界条件

为获得超声速流场的一般特性，选取3个工况点进行研究，如表 2所示，其中，Ma从5.0变化至0.8的过程中，高度从48.8 km下降至28.8 km，大气压力从92 Pa增大至1 426 Pa，而Re从3.4×10⁵增加至1.08×10⁶。在时间步的选择上，综合考虑结构尾缘非定常湍流尾迹的脉动频率以及空间再入充气系统的固有振动频率，确定每个流固耦合算例的总计算时长为10 s，时间步长为0.001 s。同时，针对各个工况点设定了0°和10°攻角，以分析非对称来流产生的效应。

2 飞行动力学

通过CFD计算得到的气动力系数如图 9所示，此处选取Ma为0.8、1.2、2.0、3.0、4.0、5.0共5种工况下不同攻角的情形进行展示。分析图中气动力系数的变化规律可以发现，对于所有的Ma工况，轴向力系数随攻角的增大均呈单调下降趋势；超声速的轴向力系数大于跨声速的轴向力系数，而在超声速区域内随着Ma的增加，轴向力系数又呈下降趋势。法向力系数的变化则有所区别，最初，对于不同的Ma，法向力系数都为0；之后，随着攻角的增加，法向力系数呈先增加后降低的趋势，极值点不尽相同，但都出现在50°攻角之后。俯仰力矩系数的变化趋势与法向力系数类似，两者本质上都反映了俯仰力或力矩大小。由于方向定义的原因，俯仰力矩的变化规律为先减小后增加，但是达到极值时的攻角大小和法向力系数曲线二者保持一致。

通过俯仰力矩分布可以初步判别其纵向稳定性。曲线在极值点处的攻角为临界攻角，那么在临界攻角之前，曲线单调下降，为负值，飞行器具有静稳定性；而当攻角大于临界攻角后，曲线单调上升，为正值，飞行器具有静不稳定性。例如，对于Ma=5.0工况，临界攻角为52°，因此在52°之前，飞行器是静稳定的，即在受到扰动使攻角发生变化时，能够有自行恢复至原平衡位置的能力。随着Ma下降，临界攻角也不断增加，因此当Ma降低时，空间再入充气系统静稳定的攻角裕度不断变大，能够在更大范围内达到静稳定状态。

在得到各气动力系数的基础上，开展飞行动力学和瞬态CFD的耦合计算，并迭代数轮，现分析迭代完成后的飞行动力学结果。图 10分别展示了超声速区域飞行过程中气动力和Re随Ma变化的曲线，以及攻角和俯仰角随时间的振荡历程。气动力衡量的是合力大小，在超声速区域内，随着Ma的减小，气动力整体呈增大趋势，最大值位于Ma= 1.2工况处。从Re曲线可以看出，随着Ma的增加，Re呈单调下降趋势，逐渐从1.08×10⁶减小至3.4×10⁵，Re的变化反映了流动现象的区别。值得注意的是，攻角和俯仰角的变化曲线都呈现明显的振荡趋势。对于攻角曲线，从17°攻角释放后围绕0°平衡位置作摇摆振动，其间在气动力的作用下不断收敛，摆动周期约为4 s。俯仰角曲线和攻角类似，只是值域上存在差异，在此不再赘述。攻角和俯仰角收敛的结果也证实了之前稳态CFD结果得到的气动力系数规律，即在50°攻角之前，俯仰力矩系数的导数为负，结构具有静稳定的特性，因此攻角和俯仰角在释放后不会发散，而是以一定的速率逐渐衰减。

3 非定常流场

相比于高超声速飞行，超声速下气流分离现象严重，极大地影响飞行器的结构稳定性，因此需要重点进行分析。以Ma为0.8和2.0为例，图 11展示了子午面流场Ma分布。在跨声速下，流动典型特征包括结构前方高压区、迎风面肩部膨胀波、背风面低速流以及湍流尾迹。来流在刚性头锥处滞止，并在迎风面重新加速，而后在肩部发生流动分离，在背风面形成较大的分离涡。而在结构后方，存在明显的非定常湍流尾迹，影响宽度甚至超过飞行器直径。

在超声速下，气流在钝体前方产生强烈的弓形激波，离体激波中间一部分很接近正激波，波后压强急剧升高，激波阻力很大。波后低速区及湍流尾迹同样存在，但尾迹宽度较跨声速有所减小。

为描述湍流尾迹的形成及发展过程，图 12选取了不同时刻的流动图像进行分析。可以看出，为了定量分析非定常气动力的变化特征，本文对飞行器轴向力进行观测，如图 13a所示。在湍流尾迹的作用下，气动力呈现明显的波动，其中轴向力幅值达到323 N。进一步探究气动力波动的频谱特征，对其进行Fourier变换，如图 13b所示。轴向力和法向力频谱带范围基本一致，都集中在100 Hz以下，其中存在15、43和61 Hz等占优频率。这表示Ma=0.8工况下湍流尾迹属于低频高幅脉动。而对于空间再入充气系统，其结构刚度较低，固有频率较低，因此有诱发共振的可能性，容易对结构造成较大破坏。

4 结构动力学

空间再入充气系统在超声速下会产生明显的振动，影响结构稳定性及各部件的应力安全。本文以Ma=0.8工况为例对结构振动特性进行具体分析。

柔性结构具有很强的材料非线性特征，所建立的流固耦合模型必须能够模拟气动力和内充压气体的共同影响。首先分析空间再入充气系统的结构变形，为了有效反映变形特征，图 14对变形前后的子午线进行了比较。可以看出，在气动力作用下，结构发生2种形式的形变，分别是气囊弯曲和蒙皮凹陷。气囊弯曲是以刚性头锥为圆心的整体向后变形，会导致结构半锥角减小；而蒙皮凹陷则是由于气囊与蒙皮间的间隙产生的局部变形。

气囊变形表征了结构整体振动的趋势，对图 14中最外层气囊上监测点A进行分析，获取轴向振动位移，如图 15a所示。轴向振动均值为54 mm，平均幅值为40 mm。对振动位移进行频谱分析，如图 15b所示。可以看出，振动频谱分布和气动力分布基本一致，占优频率带都集中在100 Hz以下，属于低频高幅振动。这表明结构振动本质上是流动分离产生的大尺度湍流尾迹诱发的抖振效应。

值得注意的是，Ma=0.8工况下气动力诱导结构振动频率较低，振幅偏高，为判断其与结构模态的耦合情况，本文对内充压作用下的模态特征进行分析。图 16中，第1阶和第2阶模态相同，固有频率均为16.434 Hz，这是由于前两阶模态对应扭转振型，而由于结构具有周向对称性，故在2个垂直方向上均有所体现。第3阶和第4阶均为伸缩模态，但伸缩程度有所区别，固有频率也不相同，分别为54.532和59.211 Hz。

可以看出，这一频率带覆盖了空间再入充气系统的前4阶固有频率，表明非定常气流脉动容易与低阶模态相耦合，容易诱发共振。这在设计过程中需格外关注。

5 飞行Mach数及攻角的影响

探究不同飞行Ma下气动弹性特征的变化，图 17给出了不同超声速条件下的流场分布，其中湍流尾迹已发展完全。由图可知，随着Ma由0.8增加至5.0，湍流尾迹尺度不断减弱，宽度及影响范围不断减小。

如图 18所示，轴向力的变化证实了这一现象。Ma=2.0工况下轴向力平均幅值仅为157 N，而Ma=5.0工况下轴向力平均仅为37 N。法向力分布类似，不再赘述。

在实际飞行中，由于气流不规则或姿态的主动调整，攻角可能会发生变化。对于空间再入充气系统，攻角变化较大，有必要进行研究。本文选取10°攻角进行初步探索。图 19展示了10°攻角下不同Ma的流动图像。当攻角存在时，尾迹发生明显变化，非对称来流会显著增加湍流尾迹的非定常效应。相比对称来流，10AOA的湍流尾迹宽度更大，影响长度更远，即使在Ma=5.0工况下，尾迹也存在较大的非定常特性。

综合各工况下结构振动特性，如表 3所示。可以看出，在超声速区域内，随着Ma的增加，湍流尾迹的尺度减弱，频率增大，导致结构振动幅值减弱。在Ma=0.8工况处，尾流频率带范围和结构前3阶固有频率较为一致，容易使结构产生共振。而攻角存在时，由于非对称来流作用下的湍流尾迹宽度更大，波及范围更远，对结构振幅影响更大。

6 结论

1) 本文建立的流固双向耦合模型能够有效解释超声速工况下空间充气式返回器的气动弹性动力响应现象。较已有方法更全面地考虑了结构变形对流场的影响，并且LES在描述具有强黏性分离特性的流场时具有良好的能力。

2) 超声速工况下，空间充气式返回器在小于50°攻角时的气动力矩导数为负，结构有维持静稳定状态的能力。同时从气动弹性动力响应计算结果来看，结构在17°攻角释放后振幅也不断收敛，进一步证实了不会发生严重的气动失稳现象。

3) 空间充气式返回器在超声速流场中会产生剧烈的振动，其本质为大尺度湍流尾迹作用下的抖振效应。在Ma为0.8工况下，结构产生的轴向振动幅值达到40 mm，而频率和前3阶固频接近，有诱发共振的潜在风险。

4) 在超声速区域内，随着Ma的增加，湍流尾迹的尺度减弱，频率增大，导致结构振动幅值减弱；而当攻角存在时，非对称来流作用下的湍流尾迹宽度更大，波及范围更远，对结构振幅影响更大。