深空行星探测实力是一个国家航天高科技发展水平的象征,也是国家提高在国际社会地位的一种重要途径,因此要发展成为如美国、俄罗斯、日本以及欧洲那些拥有空间技术实力的强国,是我国在21世纪中的重要战略目标之一.由于火星是我们的近邻,其特征在很多方面与我们居住的地球很相似,特别是美国的多项火星探测计划,使我们对火星有了更进一步的认识.探测结果表明:火星可能具有生命存在的某些必要条件,极有可能成为人类未来开发的理想星球[1].因此,探测火星成为当今和未来很长一段时间内人类进行空间探索和开发的首要目标,而中国在2020年7月23日也发射了首次自主研发的火星探测器——“天问一号”.在火星探测任务中,空气动力学减速器(一般简称为气动减速器)在探测器进行再入、降落以及着陆过程中起到了至关紧要的作用.气动减速器一般定义为是一种工作环境在5Ma速度以下的特纺材料制成的减速设备装置,其中最常见的当属于降落伞系统和充气式气动减速器(又称为减速气球).采用柔性织物材料制成的降落伞系统的伞衣具有良好的透气性以及超高的柔韧性等特点,而且它可以紧密的包装在一起,在气流的作用下可以迅速展开到比原来折叠状态大几十倍到几百倍的阻力面积的工作状态,使得单位质量产生的阻力效率远高于其他气动减速器,从而能够帮助目标载荷减速到安全的着陆速度[2].鉴于降落伞系统优良的气动减速性能、可靠性以及质量轻等特点,使其成为各项深空行星探测计划中气动减速器的首选.
在火星探测任务中,虽然刚性热防护罩能够在探测器进行再入阶段过程中保护其承受住热辐射的侵蚀,但是刚性热防护罩没有足够的阻力面积来使得探测器充分减速安全着陆到火星表面上.尽管火星表面的大气压力不及地球表面的1%,但由于降落伞系统优良的减速性能,使得降落伞系统仍是火星探测着陆过程中最经济、最有效的一种气动减速工具[3-4].降落伞系统一般的工作过程主要包括:自由降落阶段、伞绳拉直阶段、伞衣充气展开阶段和稳定下降阶段.其中,伞绳拉直阶段过程时间非常短暂,而伞衣充气展开阶段是一个降落伞系统周围流体与降落伞系统结构相互耦合、相互强烈作用的过程,有一个几何非线性与材料非线性并存的瞬态大变形结构动力学现象[5].早期的科研学者们采用理论方法、风洞实验以及飞行空投试验等对降落伞系统的动力学特性以及气动性能等进行了研究[6-8].在超音速流场中,降落伞系统伞衣与其周围流场中的非定常激波、湍流以及分离流等多物理相互作用更进一步加剧了降落伞系统伞衣周围流场物理量的剧烈变化,使得难以准确得到流场的流体结构特征[9-10].非周期性的流场流体与降落伞系统的相互作用会导致降落伞系统产生剧烈的“呼吸运动”,这将使得降落伞系统的气动减阻特性发生改变,同时也会导致伞衣材料的疲劳破坏,以至于整个降落伞系统工作的失败,所以对超音速下探测器-降落伞系统进行深入研究是航空航天领域的迫切需求,具有普适性、前瞻性和实用性,特别是对于火星探测等深空行星探测领域,更具有基础支撑作用[11].基于此,本文将研究不同流场块结构自适应网格加密精度对探测器-刚性盘-缝-带型降落伞系统的气动减速性能以及流场流体结构特性的影响,验证流场自适应网格的收敛性,为后续精确模拟计算超音速探测器-降落伞系统复杂流场流体与柔性降落伞系统的大变形相互耦合作用作准备,以期找到影响降落伞系统气动减速性能的机理原因,为未来火星探测降落伞系统的设计和制造提供基础数据支撑.
本文研究的是超音速流场中的探测器-刚性盘-缝-带型降落伞系统模型,它们的位置都是固定在物理空间中的,且没有发生位移以及变形等情形.流场流体属于三维非定常可压缩流动的流体,在超音速流场中涉及到激波或接触不连续间断等情形,所以流场流体的控制方程选择的是三维空间位置固定的极小有限控制体模型的守恒积分形式的Navier-Stokes方程,这是捕捉求解流场中存在间断不连续处的速度、压强等参量的关键.流体的通用控制方程可以表示为
(1)
式中,Q,Fc,Fv分别表示流场流体单元的守恒变量向量、对流矢通量以及黏性矢通量,源项S涉及到流体单元所受体积力以及体积加热的作用,具体形式参考文献[12-14],这里不再赘述.
采用密度加权形式(Favre过滤器)过滤方法得到的大涡数值模拟方程中的亚格子应力和亚格子标量通量统称为亚格子通量.亚格子通量是大涡数值模拟方程中的不封闭量,需要建立模型予以封闭.通常的亚格子应力的表达式为
(2)
亚格子模型的构造有很多类型,这里我们采用的是文献[15-16]提出的拉伸涡亚格子模型(stretched-vortex subgrid model).在该模型中,它假设计算网格单元内的亚格子尺度流动是由当地可分辨尺度涡和近轴对称涡的集合产生的[17].
亚格子应力可以表示为
(3)
式中,为流体运动黏性系数,eν为亚格子涡轴单元向量,为亚格子湍流能量,其表达式为
(4)
其中,kc为湍流流场最大可解波数,kc=π/Δ.假设亚格子运动为拉伸螺旋涡[18],其能量谱可以表示为
(5)
其中,K0为Kolmogorov前因子,为当地网格单元平均耗散量,为亚格子涡轴方向上的应变,
(6)
为了封闭模型,前因子组K02/3需要在流场中的每一个单元里求解计算,然后通过结构函数匹配关系[19-20],在半径为Δ的球面上的平均二阶速度结构函数F2(r)表达式为
(7)
假设,在过滤尺度Δ下,式(5)中的指数可以忽略.可以得到前因子组的表达式为
(8)
其中
(9)
采用有限体积法对流体控制方程进行离散,其中对流项采用WENO和TCD混合格式来分别处理激波/非连续流场区域和光滑连续的流场区域;黏性项采用的是中心二阶三点格式;而时间项离散采用TVD(total variation diminishing)型三步三阶强稳定性的Runge-Kutta法.在激波附近采用WENO格式计算通量其他流场区域则采用TCD格式计算通量根据压强和密度的相关曲率阈值来判断是否有激波的存在,选择对应的计算格式.其定义域C为
C={x∈R3: |αp|>cΔx2, |αρ|>cΔx2, αpαρ>0},
(10)
其中
(11)
(12)
在全部的C区域及其周围nΔx半径范围内属于不连续的区域采用WENO格式计算通量,其他区域采用TCD格式计算通量,通常c=2.5×103,n=3,则有
(13)
其中,表示C之外的区域.
本文研究的三维探测器-刚性盘-缝-带型降落伞系统,探测器的选择为类似海盗号(viking-type)探测器形状:其中前端倾角为70°,后端倾角为40°.图1为探测器-刚性盘-缝-带型降落伞系统模型二维结构示意图.可以看到,在沿着X轴正向的主流动方向上,探测器位于流场流动的上游,固定在图1的左侧;而刚性盘-缝-带型降落伞系统则固定在图1的右侧.我们将探测器最大截面形心处作为X轴方向上的原点,而该点到降落伞系统的带结构边缘的距离长度记为H;另外,将刚性盘-缝-带型降落伞系统带结构的宽度记为LB、缝结构宽度记为LG.另外,降落伞系统的伞衣厚度为5×10-4 m,降落伞系统是通过多根柔性伞绳连接固定在A点的吊绳上与探测器组成统一系统,由于伞绳的长度特征尺度与降落伞系统的长度特征尺度相比可以忽略不计,在数值模拟计算中将不考虑伞绳对探测器-刚性盘-缝-带型降落伞系统流场结构的影响.
图1 探测器-刚性盘-缝-带型降落伞系统的二维结构示意图
Fig. 1 Configuration of the 2D capsule-rigid disk-gap-band parachute system
图2是探测器以及刚性盘-缝-带型降落伞系统盘结构的几何结构示意图.在图2(a)中探测器的最大截面直径记为d,探测器的宽度记为h;而图2(b)中所示的刚性盘-缝-带型降落伞系统盘结构的公称直径记为Dp,降落伞系统盘结构顶端通气孔直径记为Dv.
(a) 探测器 (b) 降落伞系统盘结构
(a) The capsule (b) The disk geometry of the parachute system
图2 探测器以及降落伞系统盘结构的几何尺寸
Fig. 2 Sizes of the capsule geometry and the disk geometry of the parachute system
探测器以及刚性盘-缝-带型降落伞系统的模型尺寸大小参数如表1所示.
表1 探测器-刚性盘-缝-带型降落伞系统模型尺寸参数表
Table 1 Size parameters of the capsule-rigid disk-gap-band parachute system model
H/mLB/mLG/md/mh/mDp/mDv/m0.90.0630.0240.0900.0510.3080.033
本文将研究流场不同块结构自适应网格加密精度对探测器-刚性盘-缝-带型降落伞系统的气动减速性能以及流场结构特性的影响.其中,所有算例中公共不变的流场参数如下:三维的流场计算域范围为[-0.5 m,2.5 m]×[-0.75 m,0.75 m]×[-0.75 m,0.75 m],流场区域网格的划分采用的是多层块结构自适应网格加密技术,且各层网格之间的加密因子f=2,流体动力黏性系数μ=5.325 3×10-5 kg/(m·s).数值模拟流场参数如表2所示.
表2 数值模拟流场参数表
Table 2 Flow parameters for simulation
initial Mach numberMa∞Reynolds number Reinitial density ρ/(kg/m3)initial pressure q∞/Painitial velocity u∞/(m/s)initial temperature T/K2195 9290.219 910 914.1527.2172.8
数值模拟计算的初始边界条件:在流体沿X轴主流动方向的上游边界为流体流动进口边界,设置了流体的初始密度ρ=0.219 9 kg/m3,初始压强q∞=10 914.1 Pa以及初始Mach数大小为2.计算区域的其他边界为流动出口边界,即流体流动是充分发展状态的,流场流体变量沿流动方向上的梯度为零.在探测器以及刚性盘-缝-带型降落伞系统的固体壁面上,则设置为无滑移边界条件,即壁面上的速度分布为零.
探测器-刚性盘-缝-带型降落伞系统流场流体Reynolds数的表达式为
(14)
其中,d为探测器最大截面处的直径.
对于探测器-刚性盘-缝-带型降落伞系统来说,网格划分采用的是非结构类型的三角形单元表面网格,使用ICEMCFD以及GMSH软件联合进行生成.探测器-刚性盘-缝-带型降落伞系统的网格节点总数为73 626,三角形表面网格单元总数为146 064,如图3所示.
图3 探测器-刚性盘-缝-带型降落伞系统表面网格
Fig. 3 Surface meshing of the capsule-rigid disk-gap-band parachute system
表3 数值模拟算例网格参数表
Table 3 Parameters of the simulation grid
exampleMach numberadaptive basic grid X×Y×Znetwork layer nmagnitude of adaptive basic grid total number NB2120×60×603 3 000 000B-12120×60×602700 000B-22140×70×7034 200 000B-32160×80×8034 800 000
计算域中流场区域的网格划分则采用的是Berger-Oliger类型的多层块结构自适应网格加密技术.为了研究流场不同块结构自适应网格加密精度对探测器-刚性盘-缝-带型降落伞系统的气动减速性能以及流场流体结构特性的影响,本文将进行的四组算例分别记为B、B-1、B-2以及B-3,由于流场块结构多层自适应网格总数是动态变化的,因此对流场中的自适应网格总数采用大致量级来进行统计比较分析,四组数值模拟计算算例参数如表3所示.
当探测器-刚性盘-缝-带型降落伞系统的流场区域采用两层块结构自适应网格加密技术时,其示意图如图4所示.其中,图4(a)为流场两层自适应网格沿着X轴方向上的二维切面示意图,图4(b)对应的是自适应网格沿主流动方向上的三维视图.此时我们可以发现,流场中的网格分布只是在探测器以及降落伞系统的周围附近生成了精度更高的细网格,而对于超音速下的探测器前方弓形激波区域以及降落伞系统后端的尾迹复杂湍流流场区域等都没有进行相应的网格加密.
(a) 2D (b) 3D
图4 流场两层块结构自适应网格加密示意图
Fig. 4 An overall view of 2-layer block-structured adaptive mesh refinement in the flow field
注 为了解释图中的颜色,读者可以参考本文的电子网页版本,后同.
以下所有算例都将采用32核数(8*AMD、CPU Opteron 8350四核、主频2.0 GHz)并行计算集群进行数值模拟计算.其中探测器以及刚性盘-缝-带型降落伞固体部分分配的核数为2,而流场流体运动部分的计算使用的核数为30. 数值模拟计算初始时间步长为1.0×10-5 s,时间推进过程中的时间步长量级均为10-6 s.为了保证数值计算过程中的稳定性以及收敛性条件,数值模拟计算中设置的最大的CFL数为0.95,数值模拟计算设置的总计算时间为200 ms,输出的结果数为400,因此所得到的流场流体结果时间间隔为0.5 ms.
图5 两层自适应网格加密精度下的探测器-刚性盘-缝-带型降落伞系统的气动阻力变化
Fig. 5 Variation of the aerodynamic drag of the capsule-rigid disk-gap-band parachute system under 2-layer adaptive mesh refinement
首先对于算例B-1来说,此时流场动态的自适应网格总数量级约为70万,数值模拟计算总时长为1.14 d.图5为流场流体初始Mach数为2时,探测器-刚性盘-缝-带型降落伞系统的气动阻力随计算时间的变化示意图.考虑到数值模拟计算开始时的流场是不充分发展的状态,导致了降落伞系统的气动阻力出现了剧烈的振荡变化.降落伞系统所受的气动阻力最大载荷出现在最初的流场不充分发展过程中,经过短暂的振荡后从T=4 ms时开始,降落伞系统的气动阻力变化振荡波动趋于平稳,此时的流场流体逐渐过渡到了充分发展状态.通常降落伞系统的气动阻力系数CD的通用表达式为
(15)
式中,D为降落伞系统的气动阻力,C为降落伞系统在X轴流体主流动方向上的投影面积.
当流场只采用两层块结构自适应网格加密技术时,数值模拟计算得到的降落伞系统所受的气动阻力最大载荷为2 592.43 N;当流场流体过渡到充分发展状态以后,此时的探测器-刚性盘-缝-带型降落伞系统时间平均的气动阻力为766.25 N.气动阻力变化的标准差为67.76 N,即降落伞系统气动阻力振荡脉动幅值占时间平均气动阻力的8.84%.此时将表1以及表2中相应的参数代入式(15)中,计算获得的探测器-刚性盘-缝-带型降落伞系统在流场流体充分发展状态下的时间平均气动阻力系数CD为0.337.
(a) 2D (b) 3D
图6 流场三层块结构自适应网格加密示意图
Fig. 6 An overall view of the 3-layer block-structured adaptive mesh refinement in the flow field
在算例B、B-2以及B-3中,流场网格都是采用的三层块结构自适应网格加密技术,它们之间只是对应的流场底层基础网格数有所不同,如表3所示.图6所示的是算例B中流场对应的三层块结构自适应网格加密示意图.其中图6(a)为三层自适应网格沿X轴流体主流动方向上的二维切面视图,而图6(b)则对应的是该时刻三层自适应网格沿主流动方向上的三维示意图.与图4中的流场两层自适应网格精度相对比,可以明显地观察到,此时在刚性盘-缝-带型降落伞系统伞衣后端复杂尾迹流场区域的网格已经进行了较为充分地网格加密.同时,在探测器的周围也进行了较大区域更为细致地网格加密,这将大大提高该时刻下的数值模拟计算结果中网格单元内的流场物理量的精度,使得数值模拟计算探测器-刚性盘-缝-带型降落伞系统的气动减速性能以及流场流动特征等可靠性增强,数值模拟计算结果更为精确.
在算例B中,探测器-刚性盘-缝-带型降落伞系统的气动阻力的变化如图7所示.此时三维流场自适应底层基础网格数为120×60×60,流场中动态的三层自适应网格总数量级达到了300万,相比于流场两层自适应网格总数量级的70万,约增加到了原来的4.28倍,而相应的并行计算总时长也增加到了4.50 d.从图7中可以发现,在数值模拟计算的开始,降落伞系统的气动阻力振荡波动也是剧烈的,此时的降落伞系统所受的气动阻力最大载荷为6 133.26 N,当流场流体达到充分发展状态以后,降落伞系统的气动阻力较算例B中有明显的提高,且气动阻力的变化呈现出了一定的周期性规律,振荡脉动也相应地增加了.此时降落伞系统时间平均的气动阻力为1 285.23 N,其变化标准差为166.19 N,气动阻力振荡脉动幅值占时间平均气动阻力的12.93%.与算例B-1类似,计算得到的探测器-刚性盘-缝-带型降落伞系统在流场流体充分发展状态后的时间平均气动阻力系数CD为0.564.
在算例B-2中,此时的三维流场自适应底层基础网格数为140×70×70.在算例B-2中,流场动态的自适应网格总数量级达到了约为420万,相应的数值模拟计算并行计算时长也增加到了6.91 d. 图8为刚性盘-缝-带型降落伞系统的气动阻力的变化.可以观察到,降落伞系统的气动阻力变化趋势及振荡脉动幅值与算例B中降落伞系统的气动阻力变化相类似.此时刚性盘-缝-带型降落伞系统气动阻力最大的载荷为6 693.93 N,经过计算时间T=4 ms之后,流场流体过渡到了充分发展状态,此时降落伞系统的时间平均气动阻力有小幅度的增加, 大小为1 447.73 N, 其变化标准差为185.63 N, 振荡脉动幅值占平均气动阻力的12.82%.同理,对应得到的探测器-刚性盘-缝-带型降落伞系统在流场流体充分发展状态后的时间平均气动阻力系数CD为0.636.
图7 三层自适应网格加密的探测器-刚性盘-缝-带型降落伞系统的气动阻力变化(基础网格为120×60×60)
Fig. 7 Variation of the aerodynamic drag of the capsule-rigid disk-gap-band parachute sytem under 3-layer adaptive mesh refinement(the basic grid is 120×60×60)
图8 三层自适应网格加密的探测器-刚性盘-缝-带型降落伞系统的气动阻力变化(基础网格为140×70×70)
Fig. 8 Variation of the aerodynamic drag of the capsule-rigid disk-gap-band parachute system under 3-layer adaptive mesh refinement(the basic grid is 140×70×70)
为了更进一步研究流场不同块结构自适应网格加密精度对刚性盘-缝-带型降落伞系统的气动阻力的影响.在表3所示的算例B-3中,我们进一步提高了流场三层块结构自适应网格加密的底层基础网格的分辨率,考虑到计算机并行集群软硬件的条件,此时流场的自适应底层基础的网格数为160×80×80,流场动态的自适应网格总数量级达到了约为480万,此时并行计算所需的总时长为9.28 d.探测器-刚性盘-缝-带型降落伞系统的气动阻力的变化如图9所示.降落伞系统所受气动阻力最大载荷为6 919.74 N,与算例B以及算例B-2相比,降落伞系统的气动阻力有所降低,伴随着气动阻力振荡脉动幅值的加剧,此时刚性盘-缝-带型降落伞系统的时间平均气动阻力为1 187.01 N,对应的气动阻力变化标准差为255.64 N,其振荡脉动幅值占时间平均气动阻力的21.54%,此时得到相应的探测器-刚性盘-缝-带型降落伞系统的时间平均气动阻力系数CD为0.521.
前面,我们研究讨论了表3中各算例中的探测器-刚性盘-缝-带型降落伞系统的气动阻力变化以及其特点.图10中展示的是流场流体初始Mach数为2时,四组流场不同块结构自适应网格加密精度下的探测器-刚性盘-缝-带型降落伞系统的气动阻力变化对比情况.为了能突出流场块结构自适应网格精度对刚性盘-缝-带型降落伞系统的气动阻力变化等的影响,我们暂不考虑数值模拟计算开始时产生的剧烈振荡,只考虑流场流体过渡到充分发展状态之后的降落伞系统的气动阻力变化对比.因此,在图10中的计算时间是从T=4 ms开始的.从图中可以明显发现,所有算例中的降落伞系统的气动阻力,在经过了初始阶段振荡波动之后都能够保持小幅度脉动且稳定的状态,此时的流场流体处于充分发展的状态,并且刚性盘-缝-带型降落伞系统的气动阻力变化都呈现出了一定的周期性变化的特点.
图9 三层自适应网格加密的探测器-刚性盘-缝-带型降落伞系统的气动阻力变化(基础网格为160×80×80)
Fig. 9 Variation of the aerodynamic drag of the capsule-rigid disk-gap-band parachute system under 3-layer adaptive mesh refinement(the basic grid is 160×80×80)
图10 流场不同自适应网格加密精度下的探测器-刚性盘-缝-带型降落伞系统的气动阻力变化
Fig. 10 Variation of the aerodynamic drag of the capsule-rigid disk-gap-band parachute system for different adaptive mesh refinement accuracies in the flow field
数值模拟计算结果表明:探测器-刚性盘-缝-带型降落伞系统在流场两层块结构自适应网格加密精度下的气动阻力要明显小于在流场三层块结构自适应网格加密精度下的气动阻力,且气动阻力变化的脉动幅值也较小,脉动幅值大小只占时间平均气动阻力的8.84%;而在流场三层块结构自适应网格加密精度下的降落伞系统的气动阻力变化趋势大致相同,随着流场自适应网格底层基础网格分辨率的增加,可以发现,所需要的数值模拟计算并行计算时间是逐渐增加的;探测器-刚性盘-缝-带型降落伞系统所受的气动阻力最大载荷也是逐渐增加的,在流场流体过渡到充分发展状态后的降落伞系统的时间平均气动阻力先增大后减小,但降落伞系统气动阻力变化的振荡脉动幅值是先减小后增大的.在流场自适应底层基础网格分辨率最高的算例B-3中,可以发现此时的气动阻力呈现出一定的周期性变化,流场流体表现出类似于降落伞系统“呼吸运动”的类周期性变化.其归根到底的原因是,随着流场块结构自适应网格总数的增加,数值模拟计算对探测器-刚性盘-缝-带型降落伞系统中的流场流体流速、压强等物理量的计算精确度有一定的提高,流场网格分辨率越是精细,得到的流场流体细节就越精确,但是越是精细的自适应网格分辨率也会使得物理并行计算时长相应增加,同时需要消耗更多的计算内存和存储资源等.
在流场流体初始Mach数为2时,以上四组流场不同块结构自适应网格加密精度下的探测器-刚性盘-缝-带型降落伞系统的气动阻力最大载荷、计算总时长以及降落伞系统的时间平均气动阻力和对应的气动阻力系数的结果对比如表4所示.
表4 流场不同自适应网格加密精度下的降落伞系统气动性能对比表
Table 4 The aerodynamic performance comparison for the parachute system at different adaptive mesh refinement accuracies
examplemaximum load Fmax/Ntotal time Ttotal/daerodynamic drag FX/Noscillation amplitude ratio η/%aerodynamic drag coefficient CDB6 133.264.501 285.2312.930.564B-12 592.431.14766.258.84 0.337B-26 693.936.911 447.7312.820.636B-36 919.749.281 187.0121.540.521
前面我们对在流场流体初始Mach数为2时,流场不同块结构自适应网格加密精度下的探测器-刚性盘-缝-带型降落伞系统的气动阻力的变化进行了对比分析.现在我们将讨论研究流场不同块结构自适应网格加密精度对探测器-刚性盘-缝-带型降落伞系统流场流体结构特性的影响.在流场流体初始Mach数为2、计算时间T=113.0 ms时,流场在两层和三层块结构自适应网格加密精度下,探测器-刚性盘-缝-带型降落伞系统沿X轴主流动方向上的流体Mach数分布的切面示意图如图11所示.
(a) 两层自适应网格 (b) 三层自适应网格
(a) The 2-layer adaptive grid (b) The 3-layer adaptive grid
图11 计算时间T=113.0 ms时, 流场不同自适应网格加密精度下的探测器-刚性盘-缝-带型降落伞系统的Mach数大小分布切面示意图
Fig. 11 Schematic of the Mach number of the capsule-rigid disk-gap-band parachute system for T=113.0 ms and different adaptive mesh refinement accuracies in the flow field
如图11(a)中所示,流场采用的是两层块结构自适应网格加密技术,此时的探测器和刚性盘-缝-带型降落伞系统的前端都形成了明显的弓形激波,在探测器后端以及刚性盘-缝-带型降落伞伞衣内部,流体Mach数是急剧减小的.探测器后端的尾流向下游运动发展与降落伞系统伞前的弓形激波相互作用,同时降落伞系统伞衣内逆向溢出的流体使得降落伞系统伞衣前端弓形激波形状和位置都发生了改变,继而产生了一系列Mach数较小的膨胀波,沿着降落伞系统伞衣两端方向发展;刚性盘-缝-带型降落伞系统伞衣后端尾流流场是一个极度复杂的流场,在这里从降落伞系统带结构边缘以及缝结构溢出的流体与伞衣前端激波以及膨胀波等相互作用,使得流场流体发生分离,再加上从降落伞系统伞衣边缘脱落的漩涡结构以及盘结构顶端通气孔直接流过的超音速状态流体,使得降落伞系统后端尾流流场是一个湍流、漩涡以及分离流等多物理相互耦合作用的流场.而从图11(a)中可以发现,由于流场采用的两层块结构自适应网格加密技术的分辨率不高,使得降落伞系统伞衣后端的复杂流场流体细节难以充分表现,这也造成了探测器-刚性盘-缝-带型降落伞系统的气动阻力变化的模拟计算不够准确,算例B-1中的降落伞系统的气动阻力系数如表4所示,相应也是最低的.
而在图11(b)中,此时的流场采用的是三层块结构自适应网格加密技术,与图11(a)中的流场流体结构相比,此时的流场流体细节更加丰富,探测器和刚性盘-缝-带型降落伞系统前端的弓形激波更加光滑以及稳定.可以观察到,刚性盘-缝-带型降落伞系统伞衣内的部分流体从带和缝结构的边缘溢出后与伞前弓形激波相互作用,形成了Mach数较低的膨胀波,并逐渐与伞前的弓形激波相混合;另外,部分溢出的流体经过降落伞系统的盘结构边缘以及盘结构顶端通气孔与降落伞系统后端的尾流区流体相互作用,形成了大量的漩涡结构,漩涡运动发展使得尾流区里的流体流速降低了,同时也使得部分流体流向发生了改变,相应的流场流体动压也降低了,这使得刚性盘-缝-带型降落伞系统伞衣内外压强差增大了,从而使得降落伞系统的气动阻力系数增加了.在图10中,在算例B、B-2以及B-3中,流场都是采用的三层块结构自适应网格加密技术,可以发现,刚性盘-缝-带型降落伞系统的气动阻力大小变化趋势大致是相同.
(a) 两层自适应网格 (b) 三层自适应网格
(a) The 2-layer adaptive grid (b) The 3-layer adaptive grid
图12 计算时间T=113.0 ms时, 流场不同自适应网格加密精度下的探测器-刚性盘-缝-带型降落伞系统的Mach数大小分布切面网格线框示意图
Fig. 12 The wire frame of the Mach number of the capsule-rigid disk-gap-band parachute system for T=113.0 ms and different adaptive mesh refinement accuracies in the flow field
为了能够更为直观地观察研究流场块结构自适应网格加密分辨率对超音速下的探测器-刚性盘-缝-带型降落伞系统流场流体结构特性的影响,我们将研究讨论流场在两层和三层块结构自适应网格加密精度下,沿着X轴主流动方向上的流场流体Mach数分布的切面网格线框示意图.此时的流场流体初始Mach数为2且计算时间T=113.0 ms,如图12所示.
在流场自适应网格加密的线框图中,对比图12中的(a)和(b)可以得到,随着流场块结构自适应网格加密层数的增加,在三层块结构自适应网格加密的图12(b)中,流场中流体Mach数分布状况更加细致准确,其中探测器前端的弓形激波分布更加精细,激波区域处的流体Mach数分布更加细致均匀,同时在其后端的尾流低Mach数区域也可以表现得很好.从图12(b)中可以观察到,流场流体的Mach数分布范围更广泛,流体最大的Mach数为2.48,同时最小的Mach数为0.018 2,而在图12(a)中的流场流体的Mach数分布情形,其中流体最大的Mach数为2.37,对应的最小Mach数为0.034 5.采用三层块结构自适应网格加密分辨率时,更为突出的特点是在刚性盘-缝-带型降落伞系统的伞衣内部以及伞衣后端的复杂尾流区域,得到了流场流体细节特点.此时的降落伞系统伞衣内部流场流体是亚音速流动的,随着沿着X轴主流动方向上的流体不断涌进降落伞系统伞衣内部,使得伞衣内部的压强不断增加,而流场流体的Mach数是减小的.当降落伞系统伞衣内部流体达到饱和状态以后,伞衣内部流体将不断沿着降落伞系统的带边缘以及缝结构往外溢出,此时伞衣内部的流体Mach数较低,压强较高,同时流场中还存在着涡旋回流运动.在降落伞系统伞衣后端的复杂尾流区域,是沿降落伞系统伞衣边缘脱落的漩涡结构、湍流以及分离流等多物理相互作用的流场区域.其中流场中既有亚音速区也存在着超音速区,从降落伞系统盘结构顶端通气孔通过的流体还是保持着超音速状态,随后该流体与湍流、漩涡等相互作用使得流场区域流体流速降低,形成了大范围的低压区,这种相互耦合的复杂尾流逐渐向下游运动发展,使得流场流体Mach数有所增加.在图12(a)中,由于流场采用的两层块结构自适应网格加密分辨率较低,导致了在刚性盘-缝-带型降落伞系统伞衣内部以及伞衣后端复杂的尾流区域流场流体Mach数分布基本是空白的,这就难以表现突出尾流区的复杂湍流流动以及流场低压区的特性,从而造成了探测器-刚性盘-缝-带型降落伞系统的气动阻力变化振荡脉动幅值较小,且对应的流场流体时间平均气动阻力系数CD也偏小.
(a) 两层自适应网格 (b) 三层自适应网格
(a) The 2-layer adaptive grid (b) The 3-layer adaptive grid
图13 计算时间T=113.0 ms时, 流场不同自适应网格加密精度下的探测器-刚性盘-缝-带型降落伞系统流场流体速度分量v的大小分布切面示意图
Fig. 13 Schematic of velocity v of the capsule-rigid disk-gap-band parachute system for T=113.0 ms and different adaptive mesh refinement accuracies in the flow field
图13展示的是流场流体初始Mach数为2且计算时间T=113.0 ms时,在流场不同块结构自适应网格加密精度下的流场流体速度分量v沿X轴主流动方向上的分布切面示意图.其中,图13(a)为流场采用两层块结构自适应网格加密精度下的流场流体速度分量v的分布切面示意图.从图上可以观察到,流场流体在探测器的后端尾流区域,大部分流体还是保持着沿X轴主流动方向上流动的,而沿着横向Y轴方向上的流体流动较少.在刚性盘-缝-带型降落伞系统伞衣内部以及伞衣后端的尾流区域,可以发现,流场中存在着多块较大区域的流动分布,而同时存在着降落伞系统的伞衣内外两面的流场流体速度分量v大小和方向分布不合理情形,流体速度分量v出现跨伞衣连续分布的状况.同时在降落伞系统伞衣后端的复杂尾流区域,由于流场块结构自适应网格加密精度的不足,也只能粗略地计算显示流场流体速度分量v分布的大体不同区域以及流体运动方向,这对于降落伞系统的复杂湍流流场的精确模拟计算是不够理想的.
在图13(b)中可以看到,在探测器的尾流区域,流场流体也是保持着沿X轴主流动方向上的流动,而沿横向Y轴方向的流体运动几乎也是可以忽略的.当流体流动到达降落伞系统伞衣内部时,此时伞衣内部的压强急剧增大,当伞衣内部流体达到饱和状态后,流体开始从降落伞系统伞衣边缘溢出,此时流场流体速度分量v是增加的.可以看到,此时的降落伞系统伞衣内部出现流场流体速度分量v增大的区域.此时的流体流动将不断从降落伞系统的带结构边缘以及缝结构溢出,在降落伞系统伞衣带结构边缘,可以观察到较大的流体速度分量v的分布,此时伞衣内部存在着沿流动方向逆向的流体流动情况,而形成漩涡结构,同时漩涡伴随着沿其他方向的流动.而在降落伞系统伞衣后端的尾流区域,因为降落伞系统结构为非流线型结构,因此在其尾流区存在着大量涡旋结构,此时的流场是分离流、漩涡等多物理结构相互作用的区域,流场细节丰富.在沿X轴主流动方向上流场流体速度分量v的大小以及方向是交替变换的,这表明该区域是存在类似卡门涡街型的漩涡运动的,伴着这些漩涡的流动发展.在沿Y轴方向上的流场流体速度分量v的方向分布也是交替变化的,此时的流体速度分量v的最大值达到了219.8 m/s.而在图13(a)中的流体速度分量v的最大值为121.7 m/s,只有流场三层块结构自适应网格加密精度下的流场流体速度分量v的最大值的55.4%.
图14为流场不同块结构自适应网格加密精度下的探测器-刚性盘-缝-带型降落伞系统表面单元所受压力分布示意图.此时的流场流体初始Mach数大小为2,计算时间T=113.0 ms.图14(a)为流场采用两层块结构自适应网格精度,可以观察到,整个探测器-刚性盘-缝-带型降落伞系统表面所受压力分布较为均匀,与流场进口边界处的流体初始压强相差不大;而在图14(b)中可以发现,在探测器的前端表面所受的压力明显增加了,同时在降落伞系统伞衣内部,此时所受到的压强增加得更为明显,特别是在降落伞系统盘结构的边缘所受的压力值更高,这与以前的科研学者们在超音速流场中所得到的降落伞系统的伞衣边缘容易振荡是一致的.
(a) 两层自适应网格 (b) 三层自适应网格
(a) The 2-layer adaptive grid (b) The 3-layer adaptive grid
图14 计算时间T=113.0 ms时, 流场不同自适应网格加密精度下的探测器-刚性盘-缝-带型降落伞系统表面单元所受压力大小分布示意图
Fig. 14 Pressure distributions on the surface of the capsule-rigid disk-gap-band parachute system for T=113.0 ms and different adaptive mesh refinement accuracies in the flow field
(a) 算例B (b) 算例B-3
(a) Case B (b) Case B-3
图15 计算时间T=113.0 ms时, 流场三层自适应网格加密精度下的探测器-刚性盘-缝-带型降落伞系统流场流体流体Mach数大小分布切面示意图
Fig. 15 Schematic of the Mach number of the capsule-rigid disk-gap-band parachute system for T=113.0 ms and 3-layer adaptive mesh refinement accuracy in the flow field
为了对比分析流场在不同的三层块结构自适应网格精度下的流场流体结构特性之间的区别差异.图15所示的是在不同三层块结构自适应网格加密分辨率下的两算例的探测器-刚性盘-缝-带型降落伞系统的流场流体Mach数分布切面示意图.此时的流场流体初始Mach数为2,对应的计算时间T=113.0 ms.其中图15(a)中为算例B中自适应底层基础网格为120×60×60;而图15(b)则对应的是算例B-3中自适应底层基础网格为160×80×80的情形.从图15中可以观察到,两算例中的探测器-刚性盘-缝-带型降落伞系统的流场结构很相似,无论是探测器前端的弓形激波形状以及位置,还是探测器后端的尾流区域流场,算例B以及算例B-3数值模拟计算的结果是一致的;而对于刚性盘-缝-带型降落伞系统的伞衣内部以及伞衣后端复杂的尾流流场区域,都存在着探测器后端尾流与降落伞系统伞衣前端弓形激波相互作用的情形,伞衣内溢出的流体与降落伞系统伞前弓形激波相互作用形成了膨胀波向下游发展,同时,从伞衣边缘脱落的漩涡结构与降落伞系统后端的复杂尾流流场相互作用,流场流体结构分布类似.只是在流场块结构自适应网格加密分辨率更高的算例B-3中,降落伞系统伞衣内部的流体结构以及后端的尾流区域流场流体结构细节更为精细,流场流体单元的物理量精度更高,同时整个流场流体的Mach数分布更加广泛.
本文研究讨论了四组流场不同块结构自适应网格加密分辨率对流场流体初始Mach数为2时的探测器-刚性盘-缝-带型降落伞系统的气动减速性能以及流场流体结构的影响.结果表明:在较低的流场块结构自适应网格分辨率下,例如流场采用两层块结构自适应网格加密精度时,是难以准确模拟计算得到降落伞系统重要的气动阻力性能以及流场的流动特征细节的.本文随后验证了流场自适应网格的收敛性,对超音速下探测器-刚性盘-缝-带型降落伞系统的精确数值模拟计算,需要确保流场网格分辨率.这对于研究超音速下探测器-刚性盘-缝-带型降落伞系统的气动减速性能以及流场中的激波、湍流以及分离流等多物理相互作用的机理是很有必要的,也为后续对超音速下探测器-降落伞系统复杂的流场流体运动与柔性降落伞系统的大变形相互耦合作用提供了数据支撑.
[1] O’FARRELL C, MUPPIDI S, BROCK J M. Development of models for disk-gap-band parachutes deployed supersonically in the wake of a slender body[C]//2017 IEEE Aerospace Conference. 2017.
[2] 王利荣. 降落伞理论与应用[M]. 北京: 宇航工业出版社, 1997.(WANG Lirong. Parachute Theory and Application[M]. Beijing: Aerospace Industry Press, 1997.(in Chinese))
[3] BAYLE O, LORENZONI L, BLANCQUAERT T, et al. Exomars entry descent and landing demonstrator mission and design overview[C]//European Space Agency. 2015.
[4] SENGUPTA A, WITKOWSKI A, ROWAN J, et al. An overview of the Mars science laboratory parachute decelerator system[C]//19th Aerodynamic Decelerator Systems Technology Conference and Seminar. Williamsburg, VA: AIAA, 2007.
[5] PETERSON C W, STRICKLAND J H. The fluid dynamic of parachute inflation[J]. Annual Reviews Fluid Mechanics, 1961, 28(1): 361-387.
[6] 余莉, 李水生, 明晓. 降落伞弹性现象对伞衣载荷的影响[J]. 宇航学报, 2008, 29(1): 381-385.(YU Li, LI Shuisheng, MING Xiao. Influence of the parachute elastic behavior on the canopy payload[J]. Journal of Astronautics, 2008, 29(1): 381-385.(in Chinese))
[7] JIN Z Y, PASQUALINI S, QIN B. Experimental investigation of the effect of Reynolds number on flow structures in the wake of a circular parachute canopy[J]. Acta Mechanica Sinica, 2014, 30: 361-369.
[8] JOHARI H, DESABRAIS K J. Vortex shedding in the near wake of a parachute canopy[J]. Journal of Fluid Mechanics, 2005, 536: 185-207.
[9] XUE X P, KOYAMA H, NAKAMURA Y. Numerical simulation on supersonic aerodynamic interference for rigid and flexible parachute[C]//AIAA Fluid Dynamics Conference and Exhibits. 2013.
[10] KARAGIOZIS K, KAMAKOTI R, CIRAK F, et al. A computational study of supersonic disk-gap-band parachutes using large-eddy simulation coupled to a structural membrane[J]. Journal of Fluids and Structures, 2011, 27(2): 175-192.
[11] BARNHARDT M, DRAYNA T, NOMPELIS I, et al. Detached eddy simulations of the MSL parachute at supersonic conditions[C]//19th Aerodynamic Decelerator Systems Technology Conference and Seminar. Williamsburg, VA: AIAA, 2007.
[12] BLAZEK J. Computational Fluid Dynamic Principles and Application[M]. USA: Elsevier, 2015.
[13] 阎超. 计算流体动力学方法与应用[M]. 北京: 北京航空航天大学出版社, 2006.(YAN Chao. Computational Fluid Dynamics Methods and Application[M]. Beijing: Beihang University Press, 2006.(in Chinese))
[14] 李涛, 随晶侠, 吴锤结. 超声速流场中6自由度物体运动的模拟研究[J]. 应用数学和力学, 2016, 37(1): 33-53.(LI Tao, SUI Jingxia, WU Chuijie. Simulation of 6-DOF rigid bodies moving in supersonic flow[J]. Applied Mathematics and Mechanics, 2016, 37(1): 33-53.(in Chinese))
[15] HILL D J, PANTANO C, PULLIN D I. Large-eddy simulation and multiscale modelling of a Richtmyer-Meshkov instability with reshock[J]. Journal of Fluid Mechanics, 2006, 557: 29-61.
[16] B, PULLIN D I, SAMTANEY R. Subgrid-scale modelling for large-eddy simulations of compressible turbulence[J]. Physics of Fluids, 2002, 14(4): 1511-1522.
[17] MISRA A, PULLIN D I. A vortex-based subgrid stress model for large-eddy simulation[J]. Physics of Fluids, 1997, 9(8): 2443-2454.
[18] LUNDGREN T S. Strained spiral vortex model for turbulence fine structure[J]. Physics of Fluids, 1982, 25(12): 2193-2203.
[19] VOELKL T, PULLIN D I, CHAN D C. A physical-space version of the stretched-vortes subgrid-stress model for large-eddy simulation[J]. Physics of Fluids, 2001, 12(7): 1810-1825.
[20] PULLIN D I. A vortex-based model for the subgrid flux of a passive scaler[J]. Physics of Fluids, 2000, 12(9): 2311-2319.