基于DSMC方法的再入飞行器微烧蚀研究.pdf

Mechanics,2023,55(9):1847-1857Gan Chi,Chen Song,Zhang Jun,DSMC study of micro-ablation for reentry vehicles.Chinese Journal of Theoretical and Applied引用格式：甘驰，陈松，张俊.基于DSMC方法的再入飞行器微烧蚀研究.力学学报，2 0 2 3，55(9)：18 4 7-18 57流体力学Sep.,2023Chinese Journal Theoretical and Applied Mechanics2023年9 月Vol.55,No.9力第55卷第9 期报学学基于DSMC方法的再入飞行器微烧蚀研究甘驰*陈松*,2)张俊十*(北京航空航天大学中法工程师学院/国际通用工程学院，北京10 0 19 1)t(北京航空航天大学航空科学与工程学院，北京10 0 19 1)摘要高超声速再入飞行器面临着严峻的气动热环境，准确预测微烧蚀过程对热防护系统的设计至关重要.由于烧蚀会改变气动外形，进而影响周围的气动热环境以及烧蚀过程本身，因此需要将烧蚀过程与流场变化进行耦合计算.文章采用基于直接模拟蒙特卡洛(DSMC)方法的开源程序SPARTA,对高超声速条件下的再入飞行器表面微烧蚀问题展开研究.为构建并测试通用的耦合算法，通过典型的一维烧蚀模型改进了SPARTA的动能烧蚀模型并采用烧蚀表面热平衡模型计算烧蚀速率，结合Marching Square算法的特点修改了网格节点的计算方法.针对柱体、球锥以及带有微小粗糙元的斜楔体等典型外形，文中计算了二维条件下不同气动外形的烧蚀过程并进行了详细分析.其中球锥截面烧蚀预测结果中沿驻点线的烧蚀面呈现出较快的衰退速率，并且与文献中驻点附近的结果吻合情况较好.斜楔体的烧蚀结果表明，微小粗糙元的附近存在着非常稀薄的流场区域，并且其与头部驻点区域会率先发生烧蚀，反映了再入飞行器表面的微烧蚀特征.烧蚀结果对高超声速下微烧蚀机理的研究以及热防护系统的设计具有参考意义.关键词烧蚀，DSMC方法，SPARTA，高超声速再入，化学非平衡流动中图分类号：V411.8文献标识码：Adoi:10.6052/0459-1879-23-112DSMC STUDY OF MICRO-ABLATION FOR REENTRY VEHICLESI)Gan Chi*Chen Song*2)Zhang Jun t(Sino-French Engineer School/School of General Engineering,Beihang University,Beijing 100191,China)(School of Aeronautic Science and Engineering,Beihang University,Bejing 100191,China)Abstractt Hypersonic reentry vehicles are subjected to a severe aerodynamic thermal environment during a high-speedflight.Accurate prediction of micro-ablation process is quite crucial for the design of thermal protection systems(TPS).Since the aerodynamic shape changes due to the ablation process,which in turn will affect the aerodynamic thermalenvironment around the aircraft and then the ablation process itself,coupled calculation between the ablation and theflow field is required.This study uses the open source DSMC kernel SPARTA to study the micro-ablation phenomenaunder extreme conditions for the surfaces of hypersonic reentry vehicles.In order to construct and test a general couplingalgorithm,significant improvements have been made to the program.The kinetic ablation model in SPARTA is improvedbased on the one-dimensional ablation model.Besides,the computing method of the ablation rate has been adapted to bemore accurate after employing the ablation surface thermal equilibrium model,and the grid corner values are modifiedcombined with the basic characteristics of Marching Square algorithm and the ablation rate.In this paper,the ablation2023-03-29收稿，2 0 2 3-0 8-12 录用，2 0 2 3-0 8-13网络版发表.1)国家自然科学基金资助项目（12 2 7 2 0 2 8).2)通讯作者：陈松，副教授，主要研究方向为稀薄气体动力学与高超声速烧蚀.E-mail:力18482023年第55卷报学学processes oftypical aerodynamic shapes such as cylinders,a blunt cone,and a wedge with a small obstacle are simulatedand analyzed under two-dimensional conditions.The predicted surface of cone section along the stagnation streamlineshows a faster recession rate,which achieves a good agreement with the result in the literature.In addition,the ablationresults of the wedge indicate the existence of an extremely rarefied flow region near the obstacle,which undergoes adistinct shape change rapidly along with the head stagnation point after the ablation process happens.These findings haveilluminated the characteristics of the micro-ablation process for the hypersonic reentry vehicles,offering valuableguidance for the design of thermal protection systems and advancing our understanding of the ablation mechanismspertinent to the hypersonic flight vehicles.Key wordsablation,DSMC method,SPARTA,hypersonic reentry,chemical nonequilibrium flow引言飞行器在以高超声速再入大气层时将会产生一道强激波，使得飞行器壁面周围温度急剧升高，导致空气发生振动激发、离解甚至电离等一系列物理化学变化!.针对这种极端的气动热环境，常常需要进行热防护系统（TPS）设计.相比于可反复使用的热防护系统而言，高超声速条件下更需要烧蚀性的热防护系统.这种烧蚀性TPS一般是多孔设计，通过材料质量的损失例如热解、氧化和升华等一系列复杂的物理化学反应来为飞行器缓解极端的热载荷.隔热材料热解后产生的气体会沿着材料向外移动到边界层区域，由温度相对较低的热解气体带走部分热量.同时由于TPS材料的多孔设计，高温的边界层气体也可能会向内流入材料微观结构中，引起更深层次的化学反应，影响结构内部的热流.因此研究局部区域的微烧蚀问题对于整个热防护系统的设计显得十分重要.当飞行器因强烈气动加热而发生烧蚀时，烧蚀外形的变化会进一步影响流场情况，因此在模拟烧蚀过程时需要进行耦合计算 2 .为此，除了进行气动热环境分析外，还需要建立准确的烧蚀模型，以及如何在烧蚀过程中高效地处理烧蚀移动边界的后退问题.Liang等 3 针对“天宫一号 的再入过程建立了一维烧蚀模型和传热模型；周印佳等 4 采用Sutton-Graves和Tauber-Sutton理论计算驻点热流以及Landau坐标系实现烧蚀动网格，通过求解瞬态有限差分热传导方程实现了烧蚀热防护一体化计算方法；李伟等 5 采用有限体积法建立了碳纤维增强复合材料的微观烧蚀数值模型，并引入了固体相体积分数与界面重构技术；胥建宇 6 提出了基于比色测温的烧蚀模型表面温度测量方法，并建立了基于LSTM神经网络的烧蚀量预测模型；王晓婕 7 组合使用生死单元法、网格自适应技术和网格重划分技术建立了气动热环境-固体导热一烧蚀一体化耦合计算方法；杨凯威等 8 针对有限元软件采用了自编热流加载和烧蚀移动边界用户子程序方法，建立了高速空气舵前缘三维烧蚀和温度场耦合计算方法.然而，在相当大的高度时，飞行器所处的气体环境十分稀薄，克努森数Kn0.01,这将使得传统的计算流体方法难以给出准确的模拟结果.在高度不断下降的过程中，虽然飞行器外部周围的流体可能整体上是连续的，但由于材料尺度较小，并且飞行器头部存在各种端口、紧固件等微小粗糙元，局部的克努森数可能仍处于过渡区甚至稀薄区.因此，在这种情况下需要采用适用于稀薄和近连续领域的具有高保真度的气体动力学方法.直接模拟蒙特卡洛方法(DSMC)于2 0 世纪6 0年代由学者Bird提出 9 ,从统计学的角度求解玻尔兹曼方程，通过模拟大量微观粒子的运动和碰撞来反映真实的流场情况.DSMC方法适用于模拟稀薄气体流动，并且发展至今已能够应用于许多领域 10 .现有的众多程序都能够很好地实现DSMC方法，如Bird开发的代码DS2V11和DS3V12,美国国家航空航天局NASA开发的代码DAC13-14,Dietrich 等 15开发的代码MONACO,Scamlon等 16-17 开发的并行代码dsmcFoam等.近年来,Plimpton等 18 基于DSMC方法开发出了开源并行程序SPARTA.SPARTA采用多层笛卡尔网格对粒子进行分组，能够高效地模拟粒子间的碰撞以及化学反应，同时允许在流体网格中加入固体壁面对象，并模拟粒子与壁面之间的碰撞与反射过程.程序采用了MarchingSquare算法来生成以及实时更新烧蚀外形，并记录粒子与壁面碰撞时传递1849甘驰等：基于方法的再入飞行器微烧蚀研究第9 期给壁面的能量，模拟烧蚀面吸收的热量来扣除网格的节点值，进而模拟烧蚀过程.该程序能够同时考虑多组分化学反应，具备模拟高超声速化学非平衡流动的能力.流场初始条件以及网格信息通过读取外部输入文件来设置，随着统一随机粒子(USP)方法的提出，SPARTA还在不断发展之中 19 .针对于飞行器头部局部粗糙区域的烧蚀问题，本文采用最新的标准SPARTA源码，并在此基础上进行了改进.首先对粒子与壁面碰撞时动能的记录规则进行了修改，避免了负能量出现的情况；其次添加了缩放系数使网格节点值与粒子动能值处于同一量级，减少一个烧蚀循环内的溢出损失；最后对程序中烧蚀速率的计算方法进行了改进，以更准确地描述烧蚀物理过程.本文模拟了几种典型几何外形的烧蚀过程，并给出了最终的烧蚀外形以及流场图像.在不考虑转换现象 2 0 以及气体离解反应和电离反应的情况下，球锥烧蚀的模拟结果与参考文献 2 1 的烧蚀结果在驻点线附近取得了较好的一致性，验证了DSMC方法以及SPARTA程序对于烧蚀问题研究的可行性.同时，通过斜楔体的烧蚀计算本文对带有微小粗糙元的飞行器表面的烧蚀问题进行了探讨.如何将更加复杂的因素加入到程序中从而给出更加准确的烧蚀计算结果也是今后主要的研究方向1研究方法1.1分子动理论与DSMC方法分子动理论是物理学中研究物质微观粒子运动行为的理论.该理论认为分子的运动是随机的，并遵循牛顿力学和统计物理学的原理.所有微观分子的运动和相互作用共同影响了气体的宏观性质，例如密度、温度和压强等性质.当气体密度很低时，分子平均自由程远大于分子自身的尺寸，连续介质假设失效，此时的气体被称之为稀薄气体.在稀薄气体领域，分子间的相互作用可以忽略不计，分子个体的行为所产生的影响更加显著.在这种情况下，分子被认为是非常小的硬球，其主要行为包括平动和碰撞.分子在平均自由程内可视为作直线运动.稀薄气体动力学的核心是玻尔兹曼方程.该方程描述了分子在流场中的运动和碰撞.然而由于玻尔兹曼方程形式复杂，直接求解相当困难.因此，许多数值模拟方法被发展出来用以求解稀薄气体问题.其中最常用的方法是DSMC方法.DSMC方法是基于统计物理学的数值模拟方法，最初被用来模拟稀薄气体的动力学行为，经过数十年的发展已经在烧蚀问题以及与流场的耦合计算等研究中具有重要作用.在该方法中，每个模拟粒子代表了大量相同的真实气体分子.粒子的运动和碰撞在一个时间步内可认为是解耦的.粒子在该步长下沿直线运动，并被随机选取来模拟分子间的碰撞.常用的分子碰撞模型有VHS模型、VSS模型.此外，碰撞过程中的化学反应由热化学模型来模拟，例如TCE模型和QK模型等.在烧蚀问题中，粒子会与固体表面发生碰撞并反射或被吸收从而反映气一固间的相互作用过程1.2SPARTA程序SPARTA是基于DSMC方法的开源并行程序，其运行效率高，用户可以方便地添加新功能.该程序中的对象可分为3类：模拟粒子(simulatedparticles)、计算网格（gridcells）以及流体中的壁面（surfaceelements).计算网格始终是结构化网格，可用于二维和三维流动的模拟.在二维流动中，网格为与两条轴相互平行的矩形，壁面外形由数个点两两连接而成的线段组成；而在三维流动中，网格则为长方体形状，壁面由众多点三三连接形成的三角形组成.此外，网格可以根据需要进行多层加密，以及在不同区域进行不同程度的加密.SPARTA的输入文件包括了许多计算参数的设置，例如计算域大小、网格尺寸、来流条件、采样设置以及需要计算的流体性质等.在流场中，粒子的性质也需要进行设置，如粒子的质量、转动自由度、振动自由度等.SPARTA使用VHS或VSS等碰撞模型来模拟粒子之间的相互作用.用于烧蚀的隐式壁面(implicit surface)信息由二进制文件存储.SPARTA工作流程示意图如图1所示.在开始计算过程之后，程序首先会将输入文件中的命令全部读入，之后再开始执行相应的命令.程序会根据输入的网格信息自动生成结构化网格并初始化.接着，程序会在指定的边界处生成模拟粒子并为其赋予初始速度等信息，然后按照该速度来移动粒子.在每个时间步内，计算域外的粒子会被移除，而与壁面碰撞的粒子会发生反射或被吸收.此外，程序会随机选取粒子对并模拟碰撞过程，碰撞后粒子的速度会发生改变.循环内程序会每隔一定时间步进行采样并在2023年第55卷力1850报学学startcalculationread all the inputsgenerate grid cells andinitial geometrygenerate particlesparticles move&collisioncomputeflowpropertiesnoNNoal?tyesend图1SPARTA基本工作流程图Fig.1The basic flowchart of SPARTA循环的最后输出计算结果，如包括粒子和壁面在内的流场图像、每个网格的坐标以及网格内流场的温度、压强、密度等信息.在输入文件中会规定总共模拟的步数Ntotal，如果当前模拟步数N达到了总共的模拟步数Ntotal，程序就会停止运算.1.3烧蚀模型目前，计算烧蚀速率的方法主要包括解析烧蚀外形方程 2 、求解烧蚀边界热平衡方程 3,2 2-2 3、采用气动热化学烧蚀方法 2 4-2 6 等.SPARTA采用类似于热平衡边界假设的烧蚀原理，认为气动加热量与烧蚀吸热量达到平衡.SPARTA采用MarchingSquare算法对烧蚀后的壁面外形进行实时更新.该算法与Marchingcubes方法类似，能够根据输入的二维标量场生成云图或等值线.该算法一般用于作地形图上的等高线或者气象图中的等压线，而在SPARTA中用于生成烧蚀外形.算法的输入可以看作是一个二维矩阵，矩阵中的每一个元素是一个节点，4 个节点可以连接形成一个矩形，相当于一个计算网格.因此每个网格都具有4 个节点值.通常情况下，该算法会被给定一个阈值，大于该阈值的节点值被视为1，小于则视为0,因此每个网格的4 个节点分别可能为0 或者1,总共有16 种不同的情况.对应这16 种情况，算法给定一张查找表，不同的情况下网格会有不同的线段几何外形，并且每个网格一定能在该索引表中寻找到相应的线段几何外形，如图2 所示.Marching Square算法的基本步骤如下：(1)遍历每一个网格，获取其4 个节点的数值；(2)根据节点数值的分布情况确定该网格对应的索引值，采用二进制编码的方式,将4 个节点的数值映射到一个4 位二进制数上，例如0 0 0 0，0 0 0 1,0010等;(3)利用预先建立的查找表,在查找表中找到对应的线段几何外形，例如水平线段、垂直线段和斜线段等；(4)最终将所有网格的线段合并起来构成整体壁面外形.在SPARTA中，被线段所围的区域被视为烧蚀体，能够与模拟粒子发生碰撞并烧蚀.主程序根据输入条件如计算域尺寸、网格尺寸和加密信息等自动生成结构化隐式网格.其中每个网格为正方形(二维计算)或立方体(三维计算),每个二维正方形网格具有4 个节点，三维立方体网格具有8 个节点.在给定阈值之后，根据每个网格的4 个节点值，使用MarchingSquare算法高效地划分出壁面外形.以二维计算为例，初始烧蚀外形的生成如图3所示.该算法能够根据实时的节点数值信息来高效地更新烧蚀外形.在整个烧蚀模拟的过程中，模拟粒子不断地在来流边界生成并以指定的速度进行运动.在每一个时间步内，粒子进行直线运动或与其他粒子发生碰撞.同时部分粒子可能会与烧蚀壁面发生碰撞并依0000000000000000000001000000000100图2 MarchingSquare算法的填充规则Fig.2 Look up table of Marching Square algorithm1851第9 期甘驰等：基于DSMC方法的再入飞行器微烧蚀研究0000000000000000111surface0generation000000111000000011100000000000图3SPARTA烧蚀外形划分示意图Fig.3The schematic of how SPARTA generates ablation geometry据壁面的温度以一定的速度发生漫反射或镜面反射.因此在每一时间步内，烧蚀外形上的不同部分会经受不同程度的粒子碰撞，对应于不同程度的气动热环境.SPARTA会记录下每个烧蚀循环内每一时间步内壁面上发生碰撞的信息，例如与粒子发生碰撞的次数，粒子碰撞前后的速度、动量等，并根据记录下来的物理量对壁面进行烧蚀，例如计算粒子传递给壁面的总动能Efluxx=jBm(crere-Vo(1)prepost2式中，vore和vost分别为粒子与壁面碰撞前后各方向post上的速度的平方之和，为用于转化动能数值的标准化系数.此外，SPARTA的烧蚀外形变化是在每次烧蚀循环后进行更新的，而每一个烧蚀循环由若干时间步组成.SPARTA烧蚀循环流程图如图4 所示.图4中N为当前模拟步数，Nac为一个烧蚀循环所包含的模拟步数，k为正整数.以使用粒子动能变化作为烧蚀判据为例，每一个时间步内粒子传递给壁面的动能都会被记录下来，直到该烧蚀循环结束.烧蚀循环结束时，每一个壁面单元所吸收的能量将会被累加起来形成一个“烧蚀量”程序会读取该壁面单元所在网格的4 个节点值，从节点值中扣除该“烧蚀量”在该烧蚀循环的最后，程序会根据更新后的节点值来生成新的烧蚀外形.例如当烧蚀循环步数Nac为100时，每10 0 个时间步内壁面外形的“耐久值”都在不断减小，最后在第10 0 步时壁面外形会被更新.SPARTA会将动能值转化成与节点值相当量级的数值并从对应的网格节点值中扣除，直至该节点值为0，然后从该网格的其他节点值中继续进行扣除.当一个网格的4 个节点值均为0 后，该网格内的壁面单元即为“完全烧蚀”，此时该网格中不再存在壁面单元，粒子则可以进入该网格与“更里面 的壁面再次发生碰撞，并将该过程循环下去直至计算结束.因此在每一个烧蚀循环内，计算网格中节点值的分布会不断发生变动，在循环结束时对壁面的外形进行更新.图5为SPARTA模拟烧蚀过程的示意图.图5中，左侧为初始的壁面外形，其内部的节点值均大于阈值，程序根据MarchingSquare算法生成了该近似于正方形的外形.经过一个烧蚀循环后，每个壁面单元经受了不同程度的粒子碰撞，受到了不同程度的气动加热，因此节点值分布发生了变化，如图5右侧所示.程序会依照新的网格节点值更新烧蚀后的壁面外形，并将该循环继续进行下去.在SPARTA的烧蚀循环中，最外层的壁面单元总是会最先被烧蚀.并且粒子碰撞最频繁、碰撞粒子动能最高的区域烧蚀量最大，以此来反映真实、物理的烧蚀过程.在SPARTA中，粒子与壁面碰撞后反射的速度由基于壁面温度的分布函数给出.在这种条件下，温度较高的区域粒子反射的速度较大，传递给壁面的能量较小，因此反而烧蚀的速率较慢，与真实的情况相背；甚至当壁面温度较高时，可能会出现碰撞后速external cyclegenerate particlesparticles move&collisionupdatesurfacecompute heat flux andgeometryreduce corner pointvaluescomputeflowpropertiesN=kNa?yesnoexternal cycle图4 SPARTA烧蚀循环示意图Fig.4Flowchart of ablation cycle in SPARTA00000000000ablation11100cycle011100.111001110011100000000000图5SPARTA模拟烧蚀过程示意图Fig.5Simulated ablation process of SPARTA1852力20233年第55卷报学学度大于碰撞前的情况，该能量值成为负数.因此本文对此项进行了优化，将式（1)中碰撞后的动能项舍弃，只考虑入射的速度，从而保证温度越高的区域烧蚀速率越快.另外，由于SPARTA中烧蚀外形的更新是在每次烧蚀循环结束之后进行，因此每个循环内只有最外层的网格承受热量，而内层的网格并不会被烧蚀.在某些情况下，粒子传递给外层网格的能量可能足以在循环结束前烧蚀掉外层网格.然而，由于外形没有及时更新，最外层的网格仍会存在并继续与粒子碰撞.在这种情况下，本应烧蚀内层网格的能量在烧蚀实际并不存在的外层网格.为了解决这个问题，本文在记录能量时添加了一个缩放系数，以确保在一个烧蚀循环内无法完全烧蚀掉一个完整的网格，从而减少能量溢出的损失 2 7=frumEmax(2)Aceldt式中，num为真实分子与模拟粒子的比值，即一个模拟粒子代表多少真实的分子；Acel为网格的尺寸，如二维方形网格的边长；dt为时间步长；Emax为模拟粒子最高的动能.此外，当前模拟中计算网格的尺寸会影响烧蚀速率.fnum中真实分子的数量与气体的数密度有关，如下式所示lxlylzPnumJrum=Pum:Vell,(3)npercellnpercell为了保证计算效率，每个网格单元中模拟粒子的数目nperel 通常是个定值.如果其值过大，模拟过程将耗费大量时间；而如果太小,DSMC方法的统计特性将会引起较大的误差.另外，在同一算例中，体积数密度pnum始终保持不变.在二维条件下，沿z方向的网格单元长度l，也是常数，可以忽略.因此，frum与x和y方向的网格单元大小l和l,的乘积有关.通过代入和简化，并假设用k表示其中的常数项，则式(1)中的系数变为IxlylzPnum=hum(4)dtadt anpercella因此，正比于和l,的乘积与壁面元素长度a的比值.而a的大小与网格尺寸非常接近，因此l和a可以同时消去，还剩下l.这表明烧蚀速率实际上受到网格尺寸的影响.较大的网格意味着更快的烧蚀速率，而较小的网格则会导致较慢的烧蚀速率.烧蚀过程本身的速率受到了模拟参数的影响，将会产生误差.为了解决这一问题，对于二维模拟，系数的分母被修改为了?.修改后保证了烧蚀速率不再受网格单元大小所影响，而是取决于烧蚀过程本身.由于SPARTA程序受到文件数据格式的限制，此前节点最大值被设定为2 55.为了考虑材料的传热特性，本文采用一维热传导方程以及烧蚀表面的热平衡方程计算烧蚀速率 3.在烧蚀表面的热交换过程包括高超声速来流的气动加热热流qaero、材料结构吸热热流qcond以及烧蚀表面的辐射热流qradqaero=qrad+qcond(5)在DSMC方法中，气动加热是由粒子与壁面碰撞前后传递的动能所实现的，因此气动加热热流为粒子传递的动能流量之和Aaro=Z;Eflux(6)烧蚀表面的辐射热流与壁面的温度相关qrad=G:T.w4(7)其中为Stephan-Boltzmann常数，其值为5.6 7 10-8W/(mK4);为发射率,介于0 1之间.当壁面的温度到达最高点后，结构通过热传导方式所吸收的热量可以认为完全用于结构的烧蚀，即qcond=pVQ(8)式中，V为烧蚀表面的后退速率，Q为壁面单位质量的材料烧蚀时所吸收的热量.因此通过烧蚀表面的热平衡方程可以求出烧蚀进行的速率V为V=qaero-qcondlaero-OsTw4(9)PQPQ在以Marching Square算法实现的烧蚀问题中，网格的节点值G将根据该烧蚀速率来设置1xGZ,Emux(10)Vdt以此来实现更加物理的烧蚀过程2典型外形烧蚀模拟2.1二维方柱烧蚀模拟基于SPARTA的烧蚀模型，本文就一些简单的外形进行了烧蚀模拟.表1展示了来流条件以及1853甘驰等：基于DSMC方法的再入飞行器微烧蚀研究第9 期SPARTA中的网格尺寸、时间步长等条件.烧蚀材料采用典型的碳/碳复合材料的性能参数，其密度为18 0 0 kg/m3,表面发射率为0.8.材料在高空稀薄大气环境下发生化学烧蚀，大气的组成为77%的N2和2 3%的0 2，反应体系中包括碳的燃烧、碳氮反应以及碳的升华等反应，在该计算条件下Q=36.95MJ/kg.来流的气动加热热流qaero由SPARTA计算得到，在该条件下驻点处的气动加热热流为6.58 6 MW/m，除去热辐射所耗散热流后，剩余的热量完全用于材料结构的烧蚀，由此可以计算出烧蚀速率约为0.0 9 8 3mm/s，以及网格的节点数值，最后对SPARTA的烧蚀过程进行“标定”SPARTA输出了整个烧蚀过程中烧蚀外形的变化图6 和某一时刻流场的速度分布图7.高超声速来流会在方柱前方形成脱体弓形激波，激波正后方温度急剧上升.在速度分布图中，方柱前缘有较大一块区域内粒子速度明显较低，而模拟粒子在方柱的棱角处的平均速度要高于前缘中心附近.因此棱角处表现出的平均温度更高，烧蚀后退位移多于前缘中心，前缘将会逐渐烧蚀成弧形.与方柱侧边和后缘碰撞的粒表1SPARTA计算参数Table 1Simulation parameters in SPARTAParameterValueDescriptionPo1.433 1020number density/m-3To187particletemperature/KVoo6813stream velocity/(ms-l)dx,dy0.0025grid size/mdt4 10-7time steplength/sTw1000wall temperature/K图6 二维方柱烧蚀过程Fig.6The ablation process of the rectangular cylinder0.5W/0-0.5-0.6-0.4-0.200.20.4X/mvelocity/(msrl)100020003000400050006000图7 二维方柱绕流速度场云图Fig.7The velocity contour of the rectangular cylinder子数量远小于前缘与棱角，并且碰撞时的入射速度较小，所以靠近后缘处的烧蚀量较少，基本维持着原来的外形.根据烧蚀外形的变化图6 可以发现，烧蚀过程与速度场云图所展示的情况基本符合，方柱截面逐渐变成了“子弹”形状，前缘形成了圆弧形状，反映了SPARTA烧蚀模型的有效性.2.2二维圆柱烧蚀模拟本文同时模拟了二维圆柱的烧蚀，其计算参数与方柱烧蚀相同.在圆柱烧蚀模拟中激波的形状更加贴近前缘圆弧形，因此圆柱表面的温度分布更加均匀.圆柱前缘驻点处承受的温度最高，并且沿着两侧弧形温度逐渐下降.肩部的温度要低于头部驻点，所以在圆柱表面上头部驻点处的烧蚀量最大，离驻点越远的地方烧蚀量越小，整体的烧蚀外形表现出弧形逐渐“变平”的特征.图8 所展示的白色区域为圆柱截面的外形变化.图9 展示了烧蚀过程中某一时刻的温度场分布情况.在不考虑化学反应以及电离、离解反应的情况图8 二维圆柱烧蚀过程Fig.8The ablation of the two-dimensional cylinder1854力2023年第55卷报学学0.50.40.30.20.10-0.1-0.2-0.3-0.4-0.5-0.6-0.4-0.200.20.4X/mtemperature/K0.51.01.52.02.5104图9 二维圆柱绕流温度场云图Fig.9The temperature contour of the two-dimensional cylinder下，激波后驻点线上的最高温度达到了2 50 0 0 K.前缘表面上的温度远高于侧后方，因此在整个过程中后退位移要远多于侧后方，使得圆柱的前缘表面呈现明显的扁平形状，而后缘表面仍基本保持圆弧形状.图8 的烧蚀结果和温度场的结果基本吻合.2.3二维球锥烧蚀模拟本文就典型的球锥气动外形 2 1,2 8 的高超声速再入过程进行了模拟.该球锥头的截面外形的头部半径为0.0 32 m，锥角9.8，其示意图如图10 所示.因为球锥的外形为三维轴对称外形,因此在模拟过程中选取其中一个对称截面进行计算.球锥头从70km高度的轨道以58 0 0 m/s的速度再入大气层在该高度下，大气的密度为8.310-5kg/m，气体十分稀薄，DSMC方法能够很好地适用.表2 给出了模拟再入时的稀薄气体环境条件以及SPARTA模拟时的计算网格设置由于选取的球锥外形的尺寸较小，相应地，为保证计算的准确度，网格尺寸设置得足够小，并且小于分子的平均自由程；时间步长的选取同样小于分子平均自由时间.壁面采用恒温条件，在烧蚀过程中壁6R=0.032m图10 球锥截面气动外形Fig.10Aerodynamics geometry of the blunt cone section表2 来流条件以及SPARTA设置Table2Free stream conditions and SPARTA parameters of thesimulated blunt coneParameterValueDescriptionP1.7 1021number density/m-3h70height/kmTo219.58particle temperature/KVoo5800stream velocity/(ms-l)Twll4000temperatureofthewall/Kdx,dy0.0005grid size/mdt4 10-7timestep length/s面的平均温度可达到4 0 0 0 K左右，因此认为在烧蚀过程中壁面温度保持在该水平上.来流的气动加热热流最高为12 6.0 0 9 MW/m，此时沿驻点线的烧蚀速率约为1.7 2 mm/s.烧蚀过程中某时刻的温度云图如图11所示.SPARTA计算的最终烧蚀外形以及参考文献 2 1 中的原始外形和烧蚀外形由图12 给出.烧蚀过程中初步考虑了O2,N2,O,N,NO五组分的化学反应.不过图11的结果中显示的温度最高达到了14 0 0 0 K以上，事实上在该温度下气体已经发生了电离和离解反应，反应过程会吸收很大一部分热量，因此真实情况下的最高温度将低于14 0 0 0 K.由图12 可以看出，在头部驻点线附近，SPARTA模拟得到的烧蚀外形与参考文献结果相符，驻点后退量为9.50 mm，而文献中的烧蚀后退量为9.7 9 mm，误差在3%以内.在再入过程中，球锥前方形成了弓0.080.060.040.020-0.02-0.04-0.06-0.08-0.10-0.0500.050.10X/mtemperature/K200060001000014000图11球锥截面烧蚀温度场云图Fig.11The temperature contour of the ablated blunt cone section1855第9 期甘驰等：基于DSMC方法的再入飞行器微烧蚀研究1.21.00.80.60.4.originalaerodynamicshape.surface simulatedbySPARTA0.2resultsgivenbyRef.21000.51.01.52.02.5X/R图12 SPARTA预测烧蚀外形与文献 2 1 结果对比Fig.12 The surface recession predicted by SPARTA and thecomparison with the result in Ref.21形脱体激波，激波后头部承受的热流最大，因此烧蚀量最大.肩部区域所承受的热流低于头部区域，因此烧蚀后退位移少于头部.对比参考文献中的烧蚀外形，SPARTA计算的烧蚀外形整体呈现较为平滑的弧线型，而文献中的结果在肩部附近呈现“削平”的状态.本文结果中肩部区域的烧蚀量明显较少，这是由于参考文献中考虑了转换现象的影响.肩部的转换点后的热流会迅速升高，局部的烧蚀量较大，因此肩部区域会产生凹陷.并且由于外形的变化，转换点不断沿肩部移动，使得整个肩部区域的烧蚀量要多于不考虑转换现象的情况，导致外形越来越尖.所以SPARTA计算的外形仍呈圆弧形，而参考文献为典型的“双锥形”.SPARTA计算中暂未考虑转换现象的影响，因此导致了部分误差.如何将复杂的转换机制加入SPARTA程序中也是今后的改进之处.3带有微小粗糙元的二维斜楔体烧蚀模拟本文模拟了二维条件下带有微小的方块状粗糙元的斜楔体烧蚀 2 9