空间再入充气结构的流固及热固单向耦合研究

张 章, 吴 杰, 侯安平, 王立武, 竺梅芳

(1. 北京空间机电研究所 中国空间技术研究院航天器无损着陆技术核心专业实验室, 北京 100094；2. 北京航空航天大学 能源与动力工程学院 航空发动机气动热力国家级重点实验室, 北京 100191)

0 引 言

随着航天事业的发展，空天之间的物资来往越来越频繁，在当今的航天器返回过程中，普遍使用刚性飞行器回收太空物资，但这一技术正迅速接近其有效载荷极限。刚性飞行器一般被安装在运载火箭内送入太空，因而结构的最大直径尺寸受限，从而影响其承载能力[1-2]。空间再入充气结构是美国宇航局研发的一种新兴的再入与返回式航天器，利用充气薄膜的可压缩性，在发射过程中折叠以占有尽可能小的体积，返回时通过充气使得具有防热性能的柔性结构形成气动外形，并通过这种充气外形得到热防护和气动减速功能[3-4]。具有质量轻、发射体积小、结构简单、灵活机动、有效载荷比大等优点，可以有效降低返回成本，并能够应用于多种复杂的再入返回任务中[5]。空间再入充气结构在返回过程中面临严峻的飞行环境，在进入大气层后需承受超高声速的气流冲击和巨大的气动热载荷[6]。高超声速气动力会导致阻力型面的变形，产生较大的结构应力，影响结构强度和固有振动特性[7-8]；气动热载荷会使得充气结构驻点温度上升、内充压气体受热剧烈膨胀，极易造成柔性充气结构破坏，影响再入返回任务的顺利进行[9-11]。因此，考虑气动力载荷下的非线性结构动力学研究成为空间再入充气结构设计阶段不可忽视的环节[12]。

针对空间再入充气结构的结构动力学特性，国内外学者已开展了相应的研究工作。Lindell采用有限元方法对空间再入充气结构进行了结构分析，研究了在给定内充压情形下的应力分布和模态振型[13]；Kinney在有限元数值建模中将柔性充气薄膜简化为拉伸弹簧和弯曲弹簧两种单元类型，开展了稳态、瞬态结构动力学分析[14]；许家裕在获得高超声速定常气动力结果后进行瞬态结构动力学计算，分析了IRVE系统的瞬态响应情况[15]，但其仅将最大压力加载到迎风面头部，未完整地考虑高超声速流场真实气动力、气动加热对空间再入充气结构的结构动力学特性的影响。

本文以美国第三代充气式返回飞行器IRVE为研究对象，计算了一定初始再入条件下的二维平面弹道方程。以此为出发点，通过CFD数值模拟研究了再入过程中的阻力系数变化曲线，分析了半锥角、初始再入角的改变对飞行过载和结构顶端热流密度的影响关系。同时基于带预紧力的薄膜振动方程和非线性有限元理论，建立了空间再入充气结构的结构动力学模型，研究其在给定充气压力下的结构静力学特性和固有模态特征，并比较了不同充气压力和薄膜厚度对静力学和模态特征的影响。最后本文运用流固与热固单向耦合的方法，建立相应的耦合分析模型，以CFD数值模拟得到的真实气动力作为结构有限元分析的输入，研究了再入过程中各高度处的气动力、气动热载荷作用下的结构热应力及热模态变化，极大程度上描述了真实气动力对空间再入充气结构材料非线性结构动力学行为的影响。

1 数值方法和计算模型

1.1 技术路线

本文针对空间再入充气结构的弹道方程推导、流场数值模拟及结构有限元计算等具体研究内容，采用如图1所示的技术路线。具体方法是：基于空间再入充气结构在返回过程中的运动方程自主编写飞行弹道计算程序，求解飞行轨迹中高度、速度及过载等运动参数随时间的变化关系；以弹道方程所得结果为边界条件，结合国际标准大气表，对空间再入充气结构所处流体域进行CFD数值模拟，得到飞行过程中阻力系数、顶端热流密度的变化，以及结构表面气动压力及温度的分布特征；通过提取CFD计算得到的空间再入充气结构表面气动力及气动热载荷分布，分别加载到有限元模型上进行结构动力学分析，得到考虑真实流场作用下的柔性充气结构热应力分布与热模态振动特征。

图1 计算流程图Fig.1 Flow chart of analysis

1.2 薄膜静力方程与振动方程

根据薄膜静力学方程计算求解预应力和开展静力学分析。由于再入充气结构需要内充压气体产生的张力才能维持一定的型面和刚度，因此一切结构特性计算的前提是具有一定压力的内充压气体。该结构是一种具有预应力的结构，对这一类结构进行模态分析前需进行静力学分析。薄膜静力学方程可以表示如下形式:

(1)

其中，h为薄膜厚度，p为充气压力，E为弹性模量，z为横向位移，z是x、y的函数。通过薄膜自由振动方程及系统的振动微分方程，求解固有频率。微分形式的薄膜自由振动方程可以表示为如下形式:

(2)

式中，ms为单位面积上的质量，p为充气压力，Tx和Ty为由变形引起的内部张力。

依据一定的原则把分布质量换算为集中的质量，接着像处理集中质量一样建立质量矩阵。其中换算的方法有很多，如果质量的分布是均匀的，那么只需要把质量平均分配到所有节点。至于质量分布不均匀的情况，可以预先分配给各个节点所要负责的区域，接着将各个区域的质量分配给这些节点，这样建立起来的质量矩阵具有除对角线元素外其他元素为0的特点。利用此方法，建立再入充气结构的集中质量矩阵，根据多自由度系统自由振动时候的频率行列式

[KA]-ω2[MA]=0

(3)

求解系统各阶固有频率。其中KA为系统刚度矩阵，MA为系统质量矩阵，而ω是所求的固有频率。

1.3 CFD计算方法

根据国际标准大气表以及空间再入充气结构的真实飞行弹道方程，对流场计算设置边界条件，选取合适的湍流模型和求解方程。本文中,针对不同的进口条件相应地选择亚声速或超声速边界条件:对于亚声速飞行，进口给定总温、速度，出口给定静压；对于超声速飞行，进口给定总温、总压和速度，出口不给定边界条件。高超声速气动力和气动热计算采用基于有限体积法的雷诺平均N-S方程进行求解，湍流模型采用SST模型，空间离散采用迎风格式，时间离散采用二阶欧拉后差格式来保证时间方向上的守恒。

1.4 结构热传导及内充压气体热膨胀方程

对具有一定温度边界的结构，加载温度边界条件至有限元节点上，根据材料物理特性确定热膨胀系数和导热系数，根据热传导方程计算稳态温度分布，并采用理想气体状态方程修正内充压气体参数变化，在传热计算与内充压气体参数修正的基础上计算求解热应力。其中，热传导微分方程可以表示为如下形式：

(4)

式中，T为温度，c为比热容，λ为导热系数。

理想气体状态方程可以表示为如下形式：

pM=ρRT

(5)

式中，p为充气压力，M为充压气体摩尔质量，R为气体常数8.314 J/(mol·K)，T是传热计算后得到的内充压气体温度。

热应力计算由热变形方程和胡克定律求解得到:

(6)

式中，α为热膨胀系数，E为弹性模量，σ为应力，ε为应变。

1.5 计算模型

根据空间再入充气结构外形参数，采用建模软件建立几何模型，并涵盖环形气囊、隔层、蒙皮等真实几何特征。本文中所参考的结构实物图为美国IRVE飞行器，如图2所示。

图2 IRVE结构示意图Fig.2 Schematic diagram of IRVE system

图3展示了空间再入充气结构的结构有限元模型，采用Shell 181四节点面单元，网格数约为10万。由于飞行器具有轴对称的外形，所以取飞行器的一半进行计算，并把对称面的边界条件设置为对称边界。选取计算域为直径为30.0 m、高为42.0 m的圆柱形，圆柱里又有一个直径为3.75 m、高为10.5 m的圆柱形区域，该区域的建立是为了进行近壁加密。图4是取一半后的计算流场，采用混合网格进行划分，总网格数约为800万。

图3 有限元模型Fig.3 Finite element model

图4 CFD模型Fig.4 CFD model

2 二维弹道方程计算

2.1 运动学方程推导

通过飞行动力学计算可得到弹道方程、热流密度、阻力系数以及温度分布等参数的工程算法结果，其计算流程如图5所示。具体方法是:建立空间再入充气结构在再入返回过程中的动力学方程，开发飞行弹道计算程序，结合国际标准大气表得到不同再入初始条件下的飞行轨迹及各高度状态参数；并根据气体分子疏密程度，选择不同的驻点热流计算公式，开发热流密度工程算法；采用热流密度工程算法对空间再入充气结构在自由分子流区、过渡流区、连续流区进行外热流数值估算，求出各流态刚性球头驻点的热流值，结合LEES热流分布公式即可得到非驻点区域的热流分布。

图5 轨迹计算流程图Fig.5 Trajectory calculation flow chart

对于弹道式的再入，因飞行时间短，可忽略地球自转的影响，其运动更接近平面运动，因此弹道式再入的质心运动描述，使用二维弹道方程更加简洁有效[16-17]。为了获得简化后的二维弹道运动方程，本文忽略地球自转对弹道的影响并假设飞行器以稳定的姿态飞行，0°攻角，无侧滑。基于此假设，飞行器的运动方向由再入初始条件确定。由再入点的位置矢量和速度矢量就可以唯一确定整个再入弹道所在的平面，在这个平面内可以建立一个再入坐标系，如图6所示。

图6 再入坐标系Fig.6 Re-entry coordinate

图6中oe-xeye就是由再入初始条件确定的再入坐标系，原点oe是再入点在地面的投影。Θ是速度与当地水平线的夹角，这里称为飞行路径角。βe是r和ye之间的夹角，称为射程角。经过推导，可以得到二维运动轨迹方程：

(7)

2.2 半锥角和初始再入角对过载的影响

过载直接反映了飞行器除重力外所受合力的大小，由于结构的强度存在限度，因此飞行过程中对过载峰值也有着明确的要求。对于空间再入充气结构，一般要求过载峰值不能超过10倍重力加速度[18]。本文探讨了半锥角和初始再入角改变对空间再入充气结构飞行过载的影响。图7是再入过程中不同半锥角的过载曲线，可以看出半锥角对再入段过载影响不大。其原因是，半锥角虽然会影响飞行过程中的阻力系数，但阻力的大小还和速度的平方成正比。因为当阻力系数升高的同时，加速度增大导致速度降低更多，从而阻力并没有产生较大的差异，因此在图中表现为不同半锥角的过载曲线差异较小。

图7 不同半锥角下的过载曲线Fig.7 Drag acceleration curve at different half taper angles

图8是再入过程中改变初始再入角时的过载曲线。与改变半锥角不同，初始再入角对过载曲线有显著影响。当初始再入角增大时，最大过载值迅速增加：初始再入角从-1°增加到-6°时，过载峰值从-3.03g迅速增加到-16.24g，此时结构面临严重的不稳定性。

图8 不同初始再入角下的过载曲线Fig.8 Drag acceleration curve at different initial re-entry angles

2.3 半锥角和初始再入角的最优化设计

工程中希望在满足一定的约束条件下，过载峰值能尽可能小，从而为空间再入充气结构的保形设计和结构安全性设计提供有价值的参考，并为其他性能参数的设计选取提供较大的裕度[19]。为探寻能使过载峰值达到最小的再入角和半锥角选择，本文采用遗传算法对其进行优化。遗传算法是进化计算的一个分支，是一种模拟自然界生物进化过程的随机搜索算法，借鉴生物界自然选择原理和自然遗传机制而形成的一种迭代式自适应概率性全局优化搜索算法[20]。

本文问题可表述为：在再入过程的总时间内，求解最优控制变量——半锥角和初始再入角，使给定目标过载峰值在满足状态方程、边界条件、约束条件的前提下达到最小。流程图如图9所示。具体思路是通过一定的编码方式，将半锥角和初始再入角编码成一定的二进制字符串，代表遗传过程中每一代的染色体组成。由此得到个体在问题域中的适应度值，结合从遗传学中借鉴的再造方法进行个体选择，产生一个新的近似解。这个过程使种群中个体不断进化，得到新的个体比原个体具有更小的过载峰值，从而满足研究要求。

图9 遗传算法示意图Fig.9 Genetic algorithm schematic diagram

表1展示了优化后的角度和性能参数，与没有优化相比，返回过程中的过载峰值降低到-1.96g，此时半锥角的取值为58.65°，初始再入角的取值为1.35°，此时再入总时间相比优化前有所延长。

表1 优化前后参数表Table 1 Parameters table before and after optimization

3 高超声速流场计算

3.1 气动压力及表面温度分布

将飞行弹道计算结果作为CFD计算的边界条件，对空间再入充气结构不同高度下的流场进行数值模拟。表2给出了再入过程中的马赫数变化情况，其数值从23.6逐渐下降到0.7。

表2 马赫数分布表Table 2 Mach number distribution table

图10是IRVE系统表面最大气动压力与最高温度随飞行高度的变化曲线。由图可知，随高度的增加，最大气动压力有下降趋势，而表面最大温度有上升趋势。

图10 最大气动压力和表面温度随高度变化曲线Fig.10 Maximum aerodynamic pressure and surface temperature curves vs. height

图11和图12展示了60 km高度处空间再入充气结构迎风面和背风面的气动压力和温度云图。由图可知，最大压力位于迎风面中心，大小为0.591 kPa。表面最大温度同样位于迎风面中心，为640.5 K。

图11 结构迎风面压力分布Fig.11 Pressure distributions in windward side

图12 结构迎风面温度分布Fig.12 Temperature distributions in windward side

3.2 阻力系数求解与反馈修正

为精确计算阻力系数的大小，本文提取CFD模拟后的表面压力分布，分别在迎风面和背风面处进行压力积分，得到空间再入充气结构外表面所受气动阻力的大小，结合结构尺寸及速度大小，计算阻力系数的数值，并反馈至弹道方程中进行修正。图13是不同半锥角构型时，空间再入充气结构的阻力系数随马赫数变化曲线。由图可知，在跨声速区，随着马赫数的增加，各构型的阻力系数先增大后减小；当马赫数继续增加到高超声速区时，50°半锥角构型的阻力系数有明显的下降，而其他构型的阻力系数变化幅度较小；当马赫数大于15时，各构型的阻力系数曲线差异逐渐明显；最后当马赫数超过20时，阻力系数趋于不变。

图13 阻力系数曲线Fig.13 Drag coefficient curve

3.3 热流密度计算

通过对CFD计算结果进行后处理，并改变半锥角和初始再入角，获得飞行过程中空间再入充气结构顶端热流密度的变化曲线。图14是再入过程中不同半锥角对驻点热流密度的影响曲线。结果表明，当初始再入角不变时，随着半锥角的增大，驻点的热流密度呈减小趋势，但驻点最大热流所在高度不断上升。

图15是改变初始再入角时驻点热流密度的变化曲线。结果表明，当固定半锥角时，随着初始再入角的增大，再入过程中驻点热流密度随之增大，且驻点最大热流所在高度呈下降趋势。

图14 热流密度随半锥角变化曲线Fig.14 Heat flux curve at different half taper angles

图15 热流密度随初始再入角变化曲线Fig.15 Heat flux curve at different initial angles

4 考虑内充压作用的静力学和模态分析

4.1 静力学分析

根据之前建立的有限元结构，在其内部气囊的有限元节点上加载一定压力的载荷来等效代替内充压气体。具体方法是打开预应力选项，输入材料参数，进行静力计算，得到一定内充压作用下系统应力分布。图16展示了充气压力20 kPa、薄膜厚度0.5 mm条件下的结构静应力分布，从云图可以看出结构最大应力为15.464 MPa，位于最底层气囊和隔层的连接处；最小应力为61.084 kPa，位于蒙皮上。

图16 结构静应力分布Fig.16 Static stress distribution

4.2 模态分析

根据多自由度系统自由振动时的频率行列式求解系统各阶固有频率，并对求解得到的初始解进行模态扩展，然后在后处理器中查看振型。根据此方法对充气压力20 kPa、薄膜厚度0.5 mm条件下的结构进行模态分析，得到固有频率及模态振型如表3和图17所示。由图可知，前两阶固有频率一致，为106.89 Hz，对应的是扭转振型，分别沿着两个不同的方向；三阶、四阶为频率不同的拉伸振型。

表3 前六阶固有频率Table 3 Natural frequency of first six orders

对内充气压力10～30 kPa、薄膜厚度在0.2～1 mm内的不同结构进行了计算，得到结构最大静应力和一阶固频随内充气压力和薄膜厚度的变化规律分别如图18和图19所示。由图可知，最大静应力和一阶固频都随内充气压力的增加呈线性递增关系；另一方面，随薄膜厚度的增加，最大静应力呈反比下降规律，而一阶固频以近似线性规律下降。

图17 前四阶固有振型Fig.17 Modal shapes of first four orders

图18 最大静应力和一阶固频随内充压变化曲线Fig.18 Maximum stress and first order frequency curves vs. internal pressure

图19 最大静应力和一阶固频随薄膜厚度变化曲线Fig.19 Maximum stress and first order frequency curves vs. film thickness

5 高超声速流场作用下的静力学和模态分析

5.1 考虑气动压力下的静力学和模态分析

在流固单向耦合研究中，提取气动力载荷中表面压力分布信息，通过节点插值的方法加载至有限元结构上进行静力学分析，得到其整体质量矩阵与刚度矩阵，从而求出在一定内充压气体作用下空间再入充气结构的整体静应力分布，再以静力学分析得到的预紧力为基础，对其进行模态分析，得到结构固有频率及各阶振型等模态振动特性，其具体加载示意图如图20所示。

图21展示了考虑气动压力时，结构最大静应力和一阶固频随飞行高度的变化曲线。可以得出，当飞行高度低于40 km时，最大静应力随高度增加而降低，一阶固频随高度增加而增加；当飞行高度超过40 km时，随高度增加，最大静应力变化不明显，第一阶频率略有下降趋势。

图20 表面压力加载示意图Fig.20 Surface pressure loading diagram

图21 考虑气动压力时最大静应力和一阶固频随高度变化曲线Fig.21 Maximum static stress and first order frequency curves vs. height considering aerodynamic pressure

5.2 考虑气动加热下的静力学和模态分析

5.2.1 结构热传导

在热固单向耦合研究中，基于CFD数值模拟得到的气动热载荷，提取可展开式气动减速与热防护结构的表面温度分布信息，通过节点插值的方法将表面温度加载至有限元模型上，进行热传导计算。空间再入充气结构在返回过程中会面临巨大的气动热载荷冲击，因此其结构表面通常配置热防护系统(TPS)以阻挡热流穿透。在热防护结构各功能层瞬态热传导分析中，需推导适用于多功能层的瞬态热传导方程，建立防热层、绝热层、气密层的瞬态热传导模型，设定材料参数及温度、热流边界条件，选用合适的热传导单元进行瞬态热传导的有限元分析，进而得到各功能层温度分布云图。图22展示了返回过程中，TPS最外部防热层和最内部气密层的平均温度随时间的变化，可以看出温度的传递存在明显的延时效应。

图22 TPS最外层及最内层温度变化Fig.22 Temperature variation in outer and inner layers of TPS

5.2.2 内充压气体热膨胀

通过提取最内层气密层节点温度，进行算术平均，近似模拟气动热载荷作用下的内充压气体温度，并结合理想气体状态方程计算，得到温度变化后内充压气体的压力变化情况，从而有效地模拟出内充压气体受热膨胀的物理现象。图23展示了气密层平均温度及受热膨胀后的新充气压力随飞行高度的变化曲线，可见在45 km时，气体受热程度达到最大。

5.2.3 热应力及热模态分析

由于各功能层温度变化及内充压气体的压力变化都能引起结构动力学行为发生较大改变，因此将迎风面和背风面的温度分布加载至结构有限元模型中，并改变内充压气体压力数值，进行结构热应力计算，即可得到在气动热载荷作用下，由于结构温度和内充压气体压力变化产生的结构变形及热应力分布。图24是气动热载荷的加载示意图。

图23 热膨胀下内部气体压力变化Fig.23 Pressure of internal gas after thermal expansion

通过将不同高度气动热载荷依次施加到结构上，得到随高度变化的结构热应力及一阶固频曲线，如图25所示。该曲线仅考虑了高超声速流场中的气动热载荷。与之前考虑气动力载荷相比，可以发现当飞行高度小于40 km时，气动力对结构静应力的影响更大；当飞行高度大于40 km时，气动热对结构静应力的影响更大。原因是，随高度的增加，气动压力呈下降趋势而表面温度呈上升趋势，因此在高度较低时气动压力的作用更为显著，而高度上升时气动热载荷的作用明显提升。

图25 考虑气动加热时最大静应力和一阶固频随高度变化曲线Fig.25 Maximum static stress and first order frequency curves vs. height considering aeroheating

6 结 论

(1) 本文基于空间再入充气结构的二维弹道方程，分析了半锥角和初始再入角对飞行过载和结构顶端热流密度的影响规律，并通过遗传算法对轨迹进行了优化，为该类再入飞行器的轨迹分析提供了一般性的思路。

(2) 建立空间再入充气结构的结构动力学模型，研究了不同充气压力和薄膜厚度对静力学特性和模态特征的影响规律，得出最大应力和固有频率随内充压的增加而增加，随薄膜厚度的增加而减小的结论。

(3) 通过流固及热固单向耦合的方法，研究了气动压力和气动加热对结构静应力及模态的影响。研究成果有望准确描述空间再入充气结构的材料非线性结构动力学行为，为空间充气式再入飞行器的保形设计和结构安全性设计提供有价值的参考。