火星进入器作强迫震荡运动壁面脉动压力数值模拟

2019-03-14石小潘荣吉利

宇航学报 2019年2期

石小潘，赵瑞，荣吉利，袁武，李齐

(1.北京理工大学飞行器动力学与控制教育部重点实验室，北京 100081；2. 中国航空救生研究所，襄阳 441003；3. 中国科学院计算机网络信息中心，北京 100190；4. 北京空间飞行器总体设计部，北京 100094)

0 引言

长期以来，世界各国先后展开了多次的深空探测任务。其中包括火星、金星、木星、土卫六等有大气天体的无人进入或软着陆[1]。火星作为距离地球最近的地外行星，和其他行星相比其自然环境更接近地球，因此，火星成为人类探索宇宙重要的一部分。由于火星大气比较稀薄，密度只有地球大气的1%，几乎所有的火星进入器均采用钝体设计减速，沿双曲轨道直接进入火星大气[2]。进入器进入火星大气层之后，由于进入器外形几何剖面不连续过渡，流动绕过肩部时出现流动分离，在进入器尾段形成分离区，分离区内旋涡流动的不稳定性及强脱体涡的动力学特性使得作用于进入器上气动力的脉动剧烈[3-4]，致使进入器的动态稳定性遭到破坏，进入器可能会出现极限环震荡的问题。进入器的非定常运动将会加剧进入器壁面的脉动压力环境(或称气动噪声环境)，使得进入器表面出现较大的局部载荷，引起进入器结构的抖振响应，极大地缩短材料的使用寿命。另外，脉动压力信号以噪声的形式通过透射及结构共振转变为舱内噪声，极易引起进入器内部仪器失效，导致飞行任务的失败[5-6]。因此，研究火星进入器作非定常运动时壁面所受到的脉动压力环境，对火星进入器的结构设计及舱内声振环境的预测都具有重要的意义。

脉动压力来源于飞行器边界层的转捩/湍流、分离/再附及激波震荡，具有显著的随机特性[7]。对于脉动压力的研究，目前常用的预示方法包括风洞实验、经验公式和数值计算。自20世纪50年代以来，国内外研究机构对飞行器壁面脉动压力作了诸多的研究。Plotkin等[8]通过对前人的研究成果进行分析归纳，整理了一套适用于预测航天飞机壁面脉动压力的经验算法，给出了均方根脉动压力系数、功率谱及空间相关性的计算方法。Laganelli等[9-10]分析研究了光滑/粗糙表面上附体和分离的边界层流动，归纳了一套预测壁面脉动压力分布的工程算法。徐立功等[11]对机动再入飞行器壁面所受的空气动力脉动环境进行了分析，用统一的参变量给出了高超声速范围内飞行器壁面脉动压力统计特性的工程计算公式。张志成等[12]使用工程估算公式分析了来流马赫数和来流攻角对球头双锥体表面脉动压力分布的影响。王娜等[13]通过实验对弹体模型表面的脉动压力进行研究，分析了弹体表面上脉动压力与马赫数、迎角的关系以及脉动压力的频谱特征。赵瑞等[5,14]使用隐式大涡模拟方法(ILES)对跨声速旋成体火箭外壁面的脉动压力环境进行数值模拟，验证了该方法的可行性，并对外壁面的流场结构进行分区归纳，整理改进了各分区脉动压力计算的经验公式。

对于返回舱这类飞行器，为了再入大气层时减速和防热的需要，一般采用大钝头倒锥外形[15]，其再入过程的稳定性一直以来都是国内外研究者首要关注的问题。根据已有的研究成果[16]，返回舱在亚跨超声速阶段可能会出现极限环震荡。因此，对于这类飞行器壁面脉动压力环境的研究，需要考虑飞行器非定常运动带来的影响。Ross等[17]通过风洞实验分析了雷诺数、来流马赫数、攻角等对标准返回舱外形脉动压力的影响。石小潘等[18]针对火星进入器在开伞之前的跨超声速阶段，采用DES方法分析了马赫数、攻角及配平翼展开角对壁面脉动压力环境的影响。以上对于标准返回舱及火星再入器壁面脉动压力环境的研究均是基于飞行器保持相对静止，由非定常流动诱导产生的脉动压力，没有考虑飞行器非定常运动带来的影响。

综上所述，国内外对飞行器壁面脉动压力环境作了较多的研究，并取得了显著的成果。但对于钝体外形的飞行器由于自身非定常运动对壁面脉动压力环境的影响研究相对较少。本文基于DES方法，耦合进入器运动方程和流体力学方程，利用刚性动网格技术，对火星大气下进入器作强迫俯仰震荡运动进行非定常数值模拟，综合分析进入器在开伞之前的超声速气动减速阶段壁面的脉动压力环境。

1 计算方法和模型

1.1 控制方程

流动控制方程为三维非定常可压缩N-S方程组，其在直角坐标系下的守恒积分形式可表示为：

(1)

式中：U为控制体，∂U为控制面，Q为守恒变矢量，F为矢通量，n为边界的外法向量。

采用有限体积法对控制方程进行数值离散，在数值离散过程中，空间对流项离散使用5阶的WENO格式，黏性项使用4阶的中心差分方法[19]。时间推进采用二阶隐式双时间法进行迭代[20]，内迭代选择常用的LU-SGS(lower-upper symmetric-Gauss-Seidel)方法[21]。定常计算采用一方程S-A湍流模型[22]并将其计算结果作为非定常DES计算的初始流场，DES方法的详细介绍参考文献[18]。

目前对于火星大气介质的模拟，有两种方法：一种方法是采用真实气体模型[23]，考虑气体介质中的各种化学组分；另一种方法是等效比热比模型[24]。本文使用的是基于等效比热比方法进行数值计算，其等效比热比取为1.17[24]。另外，对于火星大气环境，以CO2为主，其黏性系数随温度而改变，使用Sutherland公式来表征火星大气以CO2为主，其黏性系数随温度而改变，使用Sutherland公式来表征

(2)

式中：参考黏性系数μ0=1.48×10-5kg/m·s，参考温度T0=293.15 K。和黏性力一样，气体微团之间的热传导也是由分子运动造成，导热系数和黏性系数之间存在相似关系，在处理高速黏性流动问题，导热系数由Pr=μcp/k确定，Pr=0.71。

进入器强迫震荡运动可认为是6自由度运动的简化，忽略质心的平动，仅考虑某一方向绕质心的转动。对于单自由度强迫俯仰震荡，给定简谐振动形式为：

α=αm+α0sin(kt)=αm+θ

(3)

式中：α为攻角，αm为初始攻角，α0为攻角振幅，θ为俯仰角，k=ωLref/2V∞为缩减频率。本文计算和风洞实验条件保持一致，强迫振动的缩减频率为0.5 Hz，最大振幅2°。

1.2 动网格技术

由于进入器作强迫俯仰震荡运动，网格模型将随着改变，不同时刻必须重新生成新的网格，这将极大地增加计算的时间。采用刚性动网格技术，即令计算网格与物面刚性固连，整个非定常计算过程中，网格无须重新生成，这种方法动网格的质量与静网格质量完全一致，精确地满足几何守恒定律。实现该方法较为简单，在计算过程中对位置矢量乘以一个坐标转化矩阵，n+1时刻的网格坐标：

Gn+1=Ln+1G0

(4)

式中：Ln+1为由n+1时刻姿态角确定的转换矩阵；G0为姿态角全为0时的网格坐标。

1.3 计算模型和网格

本文研究进入器进入火星大气层之后由于自身非定常运动诱发的壁面脉动压力，同时作为对比，本文计算了进入器保持相对静止仅由非定常分离流动引起的脉动压力。计算工况如表1所示。

表1 进入器计算条件Table 1 Computation condition for the capsule

计算模型为类似文献[25]中所示的具有配平翼的进入器外形，基本尺寸如图1所示。配平翼展开后计算网格模型如图2所示，网格量总计600万，并在分离区附近(配平翼及舱体后方)进行加密，用于捕捉脱体涡结构。壁面第一层网格保证y+<1，用于捕捉黏性边界层流动。

图1 进入器模型Fig.1 Physical model of the capsule

图2 进入器计算网格示意图Fig.2 Schematic of the capsule computational grid

对进入器z=0截面进行测点采样，采样点共计29个，舱体迎风面测点为1～14，舱体背风面测点为15～29。壁面脉动压力的采样频率为200 Hz，从瞬时流场中提取壁面特征点压力信号随时间变化历程曲线，获得均方根脉动压力系数，并采用快速傅里叶变换(FFT)，研究频域空间中脉动压力的功率谱密度分布。

1.4 数值方法验证

为验证本文使用的计算方法，以具有风洞实验数据和计算数据的再入返回舱模型[26]作为验证算例。风洞实验模型及壁面测点如图3所示，风洞实验测试和计算条件见表2。针对该再入返回舱模型，基于SA-DES方法进行数值计算。

图3 风洞实验模型测点布置Fig.3 Observation points of wind tunnel model

MaU∞/(m·s-1)ρ∞/(kg·m3)T∞Re1.05341.3080.659263.031.35×1071.1326.080.0636218.671.21×1061.2380.3850.5752501.37×1071.4392.190.0428221.131.03×1062600.930.0246224.628.45×105

图4和图5分别给出了计算所得的壁面测点L1和W1的均方根脉动压力系数、风洞实验结果及Fujimoto等[26]的计算结果，通过对比可以看出数值计算所得的壁面测点的均方根脉动压力系数与风洞实验测量值及文献中Fujimoto等[26]的计算值误差相对较小，可以验证本文所采用的数值计算方法具有一定的准确性。

图4 L1测点计算与实验均方根脉动压力系数对比Fig.4 Surface pressure coefficient comparison with experiments at L1 position

图5 W1测点计算与实验均方根脉动压力系数对比Fig.5 Surface fluctuating pressure coefficient comparison with experiments at W1 position

2 结果与分析

2.1 进入器非定常运动脉动压力特征

图6 来流马赫数1.2不同状态的各测点脉动压力系数对比Fig.6 The coefficient of fluctuating pressure of different states under inflow Mach number is 1.2

在配平翼迎风面，进入器作强迫震荡运动时，攻角震荡加剧了分离区再附点的脉动，相比进入器保持相对静止时，配平翼迎风面分离区的脉动压力环境显著增强。图7和图8分别为来流马赫数为1.2时进入器保持相对静止和强迫震荡运动时不同时刻流场结构示意图。对于进入器背风面的脉动压力环境，分离区内涡流脉动是壁面压力脉动的主要原因[18]。从图7可以看出，进入器保持相对静止时，分离区剪切层位置基本一致，流场内大尺度结构变化较小，脉动压力主要由分离区内的小尺度涡流脉动所致。而进入器作强迫震荡运动时，流动分离区内大尺度结构非定常效应显著、变化剧烈(如图8所示)，致使涡流脉动更加剧烈。因此，进入器作强迫震荡运动时加剧了进入器背风壁面的脉动压力环境，使得均方根脉动压力系数为0.04～0.07，增加了近一倍。

图7 Ma=1.2进入器相对静止时不同时刻流场结构示意图Fig.7 The sketch of flowfield structure of different times at inflow Mach number 1.2 under relatively static condition

图8 Ma=1.2进入器强迫震荡时不同时刻流场结构示意图Fig.8 The sketch of flowfield structure of different times at inflow Mach number 1.2 under forced oscillation

图9为来流马赫数为3时进入器作强迫震荡运动与进入器保持相对静止时壁面测点均方根脉动压力系数对比曲线。从图9可以看出，进入器作强迫震荡运动使得舱体迎风面的均方根脉动压力系数明显增加。

图9 来流马赫数3不同状态的各测点脉动压力系数对比Fig.9 The coefficient of fluctuating pressure of different states under inflow Mach number is 3

图10和图11分别为来流马赫数为3时进入器保持相对静止和强迫震荡运动时不同时刻流场结构示意图。通过对比可以发现，进入器保持相对静止脱体激波形态和位置保持稳定时，对迎风面的压力脉动起到抑制作用，均方根脉动压力系数约为0。进入器作强迫震荡运动时，迎风面脱体激波形态和位置发生改变，脱体激波震荡将会在舱体迎风面诱导产生均方根脉动压力系数在0.03～0.09的脉动压力环境。对于进入器背风面，相比进入器保持相对静止，进入器强迫震荡运动时迎风面脱体激波震荡使得背风面高速区范围增大，低速分离区范围减小(如图11所示)。分离区范围减小使得涡流脉动减弱，但背风面分离区脉动压力环境同样受到迎风面脱体激波震荡的影响，两者共同作用下，使得进入器背风面脉动压力环境和进入器保持相对静止时差异不大，均方根脉动压力系数在0.01以下。

图10 Ma=3进入器相对静止时不同时刻流场结构示意图Fig.10 The sketch of flowfield structure of different times at inflow Mach number 3 under relatively static condition

图11 Ma=3进入器强迫震荡时不同时刻流场结构示意图Fig.11 The sketch of flowfield structure of different times at inflow Mach number 3 under forced oscillation

图12为来流马赫数为1.2和3时进入器舱翼连接段的局部流场图，可以看出在配平翼迎风面靠近翼根区存在小分离区，进入器保持相对静止时小分离区非定常效应显著，再附点剧烈脉动会在配平翼迎风面诱导产生均方根脉动压力系数为0.022(Ma=1.2)和0.002(Ma=3)的脉动压力环境[18]，而进入器作强迫震荡运动时，在配平翼迎风面的诱导产生均方根脉动压力系数高达0.06(Ma=1.2)和0.12(Ma=3)的脉动压力环境。通过对比可以看出，低马赫数时，进入器强迫震荡运动使得均方根脉动压力系数增大约3倍；高马赫数时，进入器作强迫震荡运动使得均方根脉动压力系数增大约60倍。进入器强迫震荡运动时的攻角震荡和激波震荡加剧了翼根分离区的非定常脉动，造成更为恶劣的脉动压力环境，在结构设计中需要尤为关注。

图12 进入器作强迫震荡时翼根区局部流场图Fig.12 The local flowfield of the flap-root under forced oscillation

图13 进入器强迫震荡时迎风面三个测点的脉动压力变化图Fig.13 The fluctuating pressure at three points around the windward of the capsule under forced oscillation

图14 进入器强迫震荡时配平翼迎风面测点脉动压力变化图(测点11,Ma=3)Fig.14 The fluctuating pressure for the windward side of the capsule under forced oscillation (Point 11,Ma=3)

2.2 进入器非定常运动脉动压力的功率谱密度特征

功率谱密度表征脉动能量随频率的分布特性。通过对火星进入器脉动压力的频谱特性分析，可以确定脉动能量集中的频率区间，进而对进入器进行针对性的结构设计。将进入器分为三个区域：舱体迎风面，配平翼迎风面，进入器背风面，选取各区域特征点进行功率谱密度分析。图15为上述3个区域特征测点的功率谱密度频谱在不同来流马赫数下的变化曲线。测点5位置(图15(a))位于舱体迎风面大底中心，来流马赫数较小时，脱体激波较弱且距离进入器较远，脉动能量集中在0.5 Hz，和攻角震荡频率一致；来流马赫数较大时，进入器迎风面脉动压力环境受攻角震荡和激波震荡共同影响，且激波震荡的影响要大于攻角震荡的影响，脉动能量密度在30 Hz左右达到峰值，激波震荡诱导的脉动压力在频域表现出低频高幅的特性，在结构设计中要着重注意。测点11位置(图15(b))位于配平翼迎风面分离区，来流马赫数较小时，脉动能量主要集中在0.5 Hz，但相比大底迎风面，脉动能量幅值要高出几个量级。来流马赫数较大时，功率谱密度和大底迎风面类似，但由于该区域存在激波与分离区再附点相互作用，导致脉动能量幅值比大底迎风面要大。测点17位置(图15(c))位于进入器背风面分离区，脉动压力主要受舱体分离区内涡流脉动和迎风面激波震荡双重影响，低马赫数时背风面分离区较大，涡流脉动较强，致使10 Hz左右的脉动能量密度幅值较高；高马赫数时，进入器背风面分离区减小，涡流脉动较弱，致使10 Hz左右的脉动能量密度有所减弱，但是激波震荡对于背风面分离区影响增强，使得脉动能量密度在30～40 Hz范围有所增强。