高超高焓非平衡流动数值模拟方法研究综述.pdf

研究综述力展进学2023年9 月第3 期第53 卷高超高烩非平衡流动数值模拟方法研究综述高振勋*蒋崇文李椿萱北京航空航天大学航空科学与工程学院，北京10 0 19 1摘要高性能计算流体力学（computational fluid dynamic，CFD）模拟可以与高超飞行试验、高恰地面设备实验研究相互印证，在热化学非平衡效应研究以及未来高超声速飞行器研制中将发挥更重要的作用.本文回顾了国内外在热化学非平衡流动CFD研究方面的进展，概述了相关热化学模型、数值格式研究以及CFD软件研制方面的现状和发展趋势，最后指出了今后在基础研究、软件开发、模拟应用等方面需要关注的问题.（1）在热化学模型方面，常用温度模型并不完全精确，多振动温度模型具有发展潜力但工程应用受限，态-态模型更精确但模拟技术尚不成熟，更为精确的热力学输运模型、有限速率化学反应模型、振动-离解耦合模型以及表面效应模型等是提升热化学非平衡模拟精度的重要物理模型，值得深入研究；（2）在数值方法方面，多物理场耦合模拟是高超热化学非平衡流动CFD研究的热点和趋势，对CFD方法的鲁棒性和收敛性提出了更高的要求，值得重点关注和研究，此外常用数值格式需要针对热化学非平衡流动特征进行适应性改造，RANS方法在热化学非平衡流模拟中的计算可靠性仍有待验证；（3）在数值软件方面，基于结构/非结构混合网格的数值求解器更加符合工业应用需求，未来高超数值软件需要具备稳定、鲁棒的多学科、多物理场耦合求解功能，且能够适应更大网格规模大尺度复杂外形的模拟需求；（4）可综合应用多种加速技术手段提升热化学非平衡流动数值模拟的计算效率，计算刚性是热化学非平衡流动数值模拟方法研究的共性基础问题，刚性消除方法仍需进一步研究和发展.关键词高超声速，热化学非平衡，高恰，计算流体力学（CFD)，计算刚性中图分类号：0 3文献标识码：AD0I:10.6052/1000-0992-22-051收稿日期：2 0 2 2-12-2 0；录用日期：2 0 2 3-0 4-2 7；在线出版日期：2 0 2 3-0 4-2 8E-mail:引用方式：高振勋，蒋崇文，李椿萱.高超高恰非平衡流动数值模拟方法研究综述.力学进展，2 0 2 3，53(3)：56 1-59 1Gao Z X,Jiang C W,Li C X.Review of numerical simulation methods for hypersonic and high-en-thalpy non-equilibrium flow.Advances in Mechanics,2023,53(3):561-5912 0 2 3力学进展版权所有力56220238年展学第53 卷进1引言临近空间飞行器在高超声速（Ma5）飞行条件下，在绕流中产生激波后空气温度急剧升高，驻点区域温度一般可以达到几千开甚至更高.在这一流动过程中，空气在高温条件下将出现一系列复杂的热化学现象，例如分子的内能模态被激发以及发生离解、交换、复合、电离等化学反应，使热力学性质发生显著改变，呈现非完全气体性质，Anderson（2 0 0 6）称之为高温气体效应或真实气体效应.由于上述热力学内能激发与化学反应过程伴随高速流动进行，存在耦合相互作用，且各能量模态松弛、化学反应与流动的特征时间尺度通常相近可比(Kane&Peetala2022)，根据Damkohler数定义往往处于非平衡条件，因此文献中也称这种高温气体效应下的流动耦合关系为热化学非平衡效应（董维中19 9 6，王京盈2 0 17，郝佳傲2 0 18).热化学平衡、热化学非平衡和冻结三种状态在高超绕流流场中往往同时存在，但很难严格区分，研究中普遍将平衡和冻结流视为非平衡流动的两种极限状态，并采用热化学非平衡流动统一描述全流域流动状态(李海燕2 0 0 7，丁明松2 0 19).高温气体热化学非平衡效应是高超飞行器绕流和高恰地面风洞实验流动（李海燕2 0 0 7，曾明等2 0 0 6，姚轩宇等2 0 19,Shangetal.2020)的典型特征，对高超声速飞行器气动性能会产生重要影响，主要体现在4个方面（丁明松2 0 19，莫凡等2 0 2 1)：（1）高温下非完全气体性质使流场中的激波形状、分离区大小等流动性状发生改变，进而改变了力、力矩等气动力特性，对飞行器气动布局操控设计产生重要影响；（2）分子热力学内能模态的激发与松弛以及化学反应使流场中的温度、能量分布发生改变，从而影响飞行器的气动热环境，催化复合等放热反应加剧了飞行器表面的受热情况，对飞行器的热防护和结构安全设计提出考验；（3）高温下大量粒子电离反应产生的离子和电子，在飞行器绕流中形成等离子体鞘套，使飞行器绕流流场的电磁传输特性受到影响，对飞行器的无线电通信产生屏蔽作用，导致通信黑障；（4）高温下粒子的能级跃迁产生光子辐射，使飞行器绕流环境的光学特性发生改变，严重的光辐射可以进一步加剧飞行器表面的受热情况.因此，准确预测高超声速飞行器的热化学非平衡流动特征及其相关的气动力特性、气动热环境、气动物理特性等，是开展飞行器气动布局操控、气动热防护、电磁通信设计等的关键基础研究（丁明松2 0 19)目前，高超飞行试验、高地面设备实验以及CFD数值模拟是研究热化学非平衡效应的主要手段.为检验不同高烩地面设备的测试性能以及开展天地关联性研究，欧洲宇航局于19 7 1年和19 7 3年两次开展了Electre球钝锥标准模型的飞行试验（Muylaertetal.1992).为研究高超声速飞行器再入过程的通信“黑障”问题，美国国家航天局（NationalAeronautics SpaceAdministra-tion，NA SA）兰利研究中心进行了RAM-C系列飞行试验（Grantbam1970,Jones&Cross1972).为支持月球、火星等深空探测活动的星座计划，美国航空航天学会（AmericaInstitute of Aero-nautics and Astronautics，A IA A)开展了Orion和FIRE-II等CEV系列飞行试验（MacLean et al.2008,Hashetal.2007，M a c Le a n e t a l.2 0 15).此外，美国进行了多次航天飞机飞行试验活动，完成了STS-1,STS-2,STS-3,STS-5等高超声速气动力和气动热相关的技术验证试验（Zoby1983a,1983b；Sc o t t 19 8 5).其他著名的飞行试验还有欧洲宇航局的EXPERT飞行试验（Morsaetal.563高振勋，蒋崇高超高非平衡流动数值模拟方法研究综述第3 期2014)、日本HOPE计划的OREX飞行试验（(Guptaetal.1996，Vl a s o v&G o r s h k o v 2 0 0 1）等.这些飞行试验为开展热化学非平衡效应相关的气动热力学问题研究提供了宝贵的试验数据支撑，推进了航空航天科学技术的极大发展.能够开展高超声速热化学非平衡效应研究的主要高恰地面设备（王京盈2 0 17)有美国NASA兰利研究中心的T5自由活塞激波风洞（Maus etal.1992)、美国卡尔斯潘大学布法罗研究中心的LENS系列膨胀反射激波风洞（MacLean etal.2008，M a c Le a n e t a l.2 0 15)、法国宇航局的F4电弧加热风洞（Muylaertet al.1992)、德国宇航局的HEG自由活塞激波风洞（Han-nemann2003)、俄罗斯的Y-12激波风洞(Reynier2016)、日本的HIEST自由活塞激波风洞(Ishiharaetal.2016)，以及中国科学院力学研究所的JF10和F12爆轰驱动激波风洞（姚轩宇等2019，唐志共，许晓斌等2 0 15)等.高恰地面测试设备的建立为开展高超声速飞行器研制及发展相关技术提供了可靠验证，其提供的测量数据被广泛应用于飞行试验、数值方法的验证考核.高超声速飞行试验风险较大，成本很高，数据较少，例如美国的X-51项目，2 0 0 5年至2 0 13年共进行了4次飞行试验，耗费近3.4亿美元，却仅成功试飞2 次（Cummings2022).相对而言，风洞实验可重复，成本较低，能够获取大量测量数据，因此被广泛采用.但是受运行成本和技术条件的制约，高恰测试设备的尺寸很难做大，有效运行时间也较短，无法完全复现真实飞行条件下的所有高超声速飞行状态，流场分析也存在一定的局限性（王京盈2 0 17).因此，CFD数值模拟方法成了继飞行试验测试、风洞实验测量的重要补充手段，可以与之相互印证（王京盈2 0 17).美国空军学院于2 0 18 年组建了“高超声速飞行器仿真研究所”，针对高超声速飞行器模拟研究中的诸多不足，联合相关科研机构研究提升数值模拟软硬件的能力，由此可见美国对高超声速飞行器相关CFD数值模拟技术的重视（Cummings2022).随着高性能计算机的快速发展，以及计算流体力学方法的应用成熟，CFD数值模拟在高超声速飞行器及其相关技术预研中的作用日益突出，甚至可以取代风洞实验开展某些有关高超声速物理现象的研究，例如有研究比较了地面风洞实验与CFD数值模拟对各类高超声速飞行物理现象的预示能力，结果显示CFD数值模拟在热化学非平衡效应研究方面具有更好的优势（Hirschel2005，王京盈2 0 17).随着任务要求提高，速度域、高度范围等飞行包线极度增大，高超飞行器绕流条件更加极端、物理现象愈加复杂，现有的高超飞行试验、高地面设备实验技术已无法满足下一代新型高超飞行器设计研制的需求，数值模拟技术在高超热化学非平衡相关研究中将发挥更重要的作用和更大的价值（唐志共，张益荣等2 0 15).本文从热化学模型、数值格式、数值模拟软件等方面分别详细介绍国内外在热化学非平衡流动CFD方法研究方面的进展，最后从现状出发总结和讨论今后在基础研究、软件开发、模拟应用等方面需要关注的问题，2热化学模型研究进展热化学模型是热化学非平衡流动控制方程建模的物理基础，对正确描述流动过程中的物理机理和提高非平衡流动的模拟精度具有重要意义.下文从控制方程建模过程中涉及的热力学温度模型、热力学参数模型、化学动力学模型及表面效应模型4个方面进行介绍.564力2023年展第53 卷进学2.1热力学温度模型选取合适的热力学温度模型对精确预测飞行器的气动热环境至关重要.常用的热力学温度模型包括一温度模型、两温度模型以及三温度模型.此外，还有更加复杂的多振动温度模型（王源杰2 0 16)和态-态模型（徐丹等2 0 14).图1展示了高超飞行器在不同高度和速度条件下驻点区流动的空气热化学状态以及高温空气发生典型化学反应的规律特点（Anderson2006，G u p t a e t a l.19 9 0).一般而言，一温度模型常用于模拟化学非平衡流动，如燃烧性质的化学反应流动、高度较低（空气较稠密）的高超飞行绕流等.两温度模型常用于模拟热化学非平衡流动，如高度较高（空气较稀薄）的高超飞行绕流、火星大气再入绕流（CO2第一振动能级特征温度低、火星大气相对稀薄）（刘庆宗2 0 16）等.三温度模型常用于模拟电离反应显著、辐射效应较强的热化学非平衡流动，例如飞行速度达到第二宇宙速度及以上的高超飞行绕流（王京盈2 0 17，郝佳傲2 0 18）等.然而，影响流动热化学状态的因素很多，除飞行速度、高度以外，还包括飞行器几何尺度、流场压强、温度等，因此实际模拟中还须依据流动特点进行合理近似，采取对应的热力学温度模型。目前应用最广的是Park两温度模型，该方法在数值模拟中取得了很大的成功，但在模拟某些高超流动现象时无法得出合理的计算结果，如高超声速激波脱体距离等.图2 显示了不同温度模型对EXPERT热化学非平衡绕流的影响（王源杰2 0 16).不少研究学者认为一个振动温度的假设并不十分合理，因为不同类型分子振动能激发时的特征温度不同，对化学反应速率的影响也应区别对待，为此提出了更为复杂的多振动温度模型和态-态模型，多振动温度模型（徐丹等2014)常用于模拟高恰风洞设备试验段的高超声速流动，以便考察不同分子振动能冻结对模型气动特性的影响.该方法在描述化学反应机理时仍然存在局限性，其暗含振动能级分布满足振动温度下Boltzmann分布假设，这个假设是否合理需要进一步研究.同时，随着混合气体组元个数增加，计算量呈倍数增长，对计算资源要求较高，现阶段在工程型号任务中应用有限.态-态模型在多振动温度模型的基础上进一步细化，将位于不同振动能级上的分子视为不同组元进行研究该方法在理论上展现出诸多优点，但各种内能模式能级间的跃迁过程十分复杂，因此还无法实现态一态模型与多维流场程序的耦合计算.目前，态-态模型研究主要集中在不涉及流动的零维问题，或是正激波后、边界层和准一维喷管非平衡流动问题，相关模拟技术尚不成熟（徐丹等2 0 14).2.2热力学参数模型高温气体的热力学参数、输运系数等是热化学非平衡气动热问题研究的重要参数（刘庆宗2016，Ca n d l e r 2 0 19).在热化学非平衡条件下，通常采用Dalton分压定理描述高温气体的状态方程，即混合气体中每个粒子的压强、温度和密度仍满足完全气体状态方程.对于一温度模型或者化学非平衡模型，混合气体中粒子的内能、烩、比热等关键参数一般采用分段拟合多项式方法计算，最常用的便是Chemkin模型（Guptaetal.1990,Kee1996)，以及针对火星大气的ESA数据模型（Capitelli etal.2005）等.对于两温度模型及其他多温度模型，通常基于Born-Oppen-heimer近似分解和统计热力学方法计算相关温度模态的能量项，基于所有能量模态平衡的假设，一温度模型也可以采用这种方法计算内能等热力学参数.565高振勋，蒋崇文高超高烩非平衡流动数值模拟方法研究综述第3期aIV11组分b115组分120I17组分组分N2,O2,N,O,NO,N2+2N2,O2,NN2,O2,N,O,NO,NO+,1OLN+,O+,NO+,eN2O2ONO1激波/壁面电离范围典型上o-o+e-1化学、滑移效应390N-N+e-升轨迹1热非平衡N2完全离解个9000KuX/H1开始电离振动能60激发发典型再化学非平衡离解范围N22N入轨迹壁面热辐射主导完全流峰值&热平衡4000KN2开始离解气体O2几乎完全离解1300220花学O2开始离解氧气、二次电离2500K氮气离解电离热平衡无化学反应800K振动激发OK03691215U/(kms-1)高超飞行器不同高度和速度条件下驻点区空气热化学状态标准大气压下高温空气化学反应温度范围图1不同高度和速度条件下驻点区空气热化学状态以及典型化学反应对应的温度范围（Anderson2006,Gupta et al.1990)ab多温度(No-Pref.Diss)多温度振动温度(Park)Park两温度模型催化跳跃头部区域平动温度分布云图头部区域原子氧分布云图图2不同温度模型对EXPERT热化学非平衡绕流的影响（王源杰2 0 16)目前，高温真实气体的输运模型可分为碰撞积分方法和拟合近似方法两类（郝佳傲2 0 18)碰撞积分方法（或SCEBD模型方法）物理机制清晰，但计算相对复杂，因此基于统计学方法和数学拟合方法的近似计算式应用广泛.单个粒子黏性系数的常用近似方法有Sutherland公式（李海燕2 0 0 7)、Blottner拟合式(3系数（Gupta etal.1990)和5系数（董维中19 9 6)两种）Lennard-Jones模型公式及其简化式（Gupta etal.1990）等；单个粒子热导系数的近似方法有Eu-cken经验式（王京盈2 0 17，郝佳傲2 0 18）、Gupta拟合式（Guptaetal.1990）等；单个粒子扩散系数的近似方法包含基于Lewis数的双极性扩散近似模型（董海波2 0 18)、基于Schmidt数的双极性扩散近似模型（王京盈2 0 17)等.在求得单个粒子的输运系数之后，混合气体的输运系数则采用Wilke混合律或Gupta-Yos混合律（Adhikarietal.2022）加权计算，Gupta-Yos方法被认为在高温条件下比Wilke方法更为精确最后，在碰撞积分模型中，扩散碰撞积分元2,-1)、粘性碰撞积分元,2)等关键参数也采用拟合公式(Hashet al.2007,Gupta et al.1990）计算，碰撞积分参力566展学2023年第53 卷进数是实验和数值研究的重要参数，在更高温度条件下，需要构建更为精确的计算模型（Candler2015).此外，振动能量（或电子能量）松弛过程是热力学模型的重要研究内容.在热化学非平衡流动的许多CFD研究中，电子能模态激发通常被忽略，振动能激发以及振动-离解耦合效应通常被认为是最重要的物理现象（Candler 2019).目前，平动-振动能量松弛过程通常采用Landau-Tell-er模型计算，其中振动松弛时间主要采用Millikan-White关系式及Park高温修正式计算（Hashetal.2007）针对不同粒子碰撞对给出了多种修正系数.当计及电子能量模态激发时，受制于实验数据，目前仅考虑氮气分子的振动-电子能量松弛过程（郝佳傲2 0 18).各能量模态之间的松弛速率是热力学建模的关键基础问题，需要得到更多的实验支持和验证，是今后的重要研究方向(Candler 2019).2.3化学动力学模型化学动力学模型直接影响高超声速流动中混合气体热力学性质、组元构成、温度/能量分布等，是高超飞行器气动热环境精准预测的关键模型，其中化学反应方程和化学反应速率系数是化学动力学模型研究的主要方向.常用的地球空气化学动力学模型有Gupta等的11组元2 0 反应方程模型（Gupta etal.1990)、Dunn和Kang的11组元2 6 反应方程模型(Dunn&Kang1973)、Pa r k 的系列化学反应模型（包括Park 1985/1990/2001等)(Dunn&Kang 1973,Park 1985,Park 1993)、M a r t i n-Bo r t n e r 化学反应模型（SurzhikovSh a n g 2 0 0 8)，以及基于它们进行修正的简化模型.对于火星大气，常用的化学动力学模型主要有McKenzie的CO2-N2化学反应体系模型（McKenzie1966)、Pa r k 等的CO2-N2-Ar化学反应体系模型(Parket al.1991)，以及基于两者进行截取的简化模型（刘庆宗2016).在化学动力学模型中，正反应速率系数采用Arrhenius拟合式计算，但逆反应速率系数计算包含两种方式：一种是同样采用Arrhenius拟合式，另一种是采用基于平衡常数的计算式.平衡常数的计算方法又包含Gibbs自由能方法（Anderson2006）和多项式拟合方法（Guptaetal.1990)等.完整的化学动力学模型应用范围更广，但也应根据计算条件变化进行简化应用（刘庆宗2016).依据图1和图2 所示飞行速度-高度包线和高温空气化学反应特点，一般来说，空气5组元模型主要用于模拟飞行速度不太高（飞行速度U1.55k m/s）、以离解反应为主的化学非平衡流动，空气7 组元模型主要用于模拟温度在9 0 0 0 K以下（飞行速度U511k m/s）、包含弱电离反应的化学非平衡或热化学非平衡流动，空气11组元模型用于模拟温度在9 0 0 0 K以上（飞行速度U11km/s）、包含强电离反应或强辐射效应的热化学非平衡流动.对于火星大气模拟，8 组元模型足以表征火星流场信息，因此被广泛采用（刘庆宗2 0 16).在纯CO2的风洞实验条件下，5组元模型已经足够.在考虑通信中断问题或辐射效应较强时，才需要使用包含更多组元的电离反应模型或完整的火星气体化学动力学模型（刘庆宗2 0 16).数值实验显示不同化学动力学模型对热化学非平衡流动模拟的结果存在较大差异（图3），并且指出空气Park化学动力学模型更为准确，其计算的电子数密度等与飞行试验数据更相符（丁明松2 0 19，赵法明2 0 19).对纯土源仙2U10)木用不问化子风应模坐订异的电于数岔皮567高振勋，蒋崇文高超高恰非平衡流动数值模拟方法研究综述第3期a1015b10167组分D-K非催化11组分D-K非催化D-K完全催化D-K完全催化1014Gupta非催化1015Gupta.完全催化Gupta非催Gupta.完全催化Park非催化Park非催化1014Park完全催化1013Park完全催化试验测量值试验测量值101310121012101110111010101000.10.20.30.40.500.10.20.30.40.5法向无量纲位置点（/Rn）法向无量纲位置点（y/Rn）空气7 组分模型空气11组分模型图3CO2气体的化学反应机理研究和对比试验显示火星Park化学动力学模型优于McKenzie化学动力学模型(刘庆宗2 0 16).由于化学反应对激波脱体距离、组分浓度分布、壁面热流分布等热化学非平衡流动参数的影响较大，因此工程应用中需要一种适用范围广的通用化学反应模型.反应速率系数在不同化学模型中存在较大差异，后来学者也不断对反应速率系数进行修正，但均未完全解决该问题(Candler2019，Si n g h&Sc h w a r t z e n t r u b e r 2 0 2 2).此外，伴随新的工程应用需求，例如极端高温条件电离反应流动、高温空气/燃烧反应内外混合流动、高温空气/燃气喷流混合流动以及高温空气/烧蚀产物引射混合流动等模拟问题，现有化学动力学模型均存在应用局限，化学反应模型及化学反应速率系数必须修正（胡雨濛2 0 18).由于混合气体的组分构成、热力学状态、稀薄状态、辐射效应以及流动的形态等因素均会对化学反应速率系数产生影响，因此适应复杂流动条件的化学动力学模型不是基于现有化学模型的简单迭加和筛选，且化学反应速率系数必须结合理论和实验加以校正.如图4所示，Chaudhry等（2 0 2 0 采用高精度量子化学数据及准经典轨道理论计算分子微观能级跃迁率，提出了化学反应速率更精准的MMT化学反应模型，该方法对开展相关基础研究具有一定参考价值.最后，化学反应与热力学松弛存在密切的耦合作用关系，这种耦合效应是影响化学反应速率的重要因素，也是热化学非平衡效应物理建模的重要研究内容，包括振动-离解耦合、电子-电离耦合等（郝佳傲2 0 18).目前，多数数值模拟研究主要关注振动-离解耦合效应，通常采用模化方法，如Park振动-离解耦合模型、Marrone-Treanor模型（Saboetal.2022）以及Machere-Fridman模型（董维中19 9 6).其中，Park振动-离解耦合模型采用平动温度与振动温度的算术平均值以描述化学反应的控制温度，这种方法简单易实现，且能够在大多数情况下得到合理的预测结果，因此被数值模拟软件广泛采用.对于平动温度的指数有不同取值方法，如DPLR取0.5力568展2023年第53 卷进学ab1500070Park化学模型Park化学模型MMT化学模型MMT化学模型口实验值60口实验值50(z-.MX)/洋10000/4030中50002010000.050.100.1500.050.100.15轴向位置点/m轴向位置点/m壁面压力分布壁面热流分布图4不同化学反应模型对高烩双锥流动壁面参数的影响比较（Chaudhryetal.2020)(Hashetal.2007)，Sh a r ma 等取0.6,LAURA取0.7，Jaffe等人认为振动温度对离解化学反应正反应速率的影响较小，因此建议采用0.5的取值方法（Guptaetal.1990，李鹏等2 0 2 1a).尽管Park振动-离解耦合模型在CFD中得到了广泛应用，但更为精确的模型方法仍需进一步研究.关于振动-离解耦合、电子-电离耦合效应的更多专业研究可参考郝佳傲（2 0 18）的文章2.4表面效应模型表面效应主要包括壁面催化、烧蚀反应和稀薄滑移等，也是影响热化学非平衡气动热环境预示精准性的重要因素，在当前工程应用研究中非常受关注（Candler2019).壁面催化条件用于预测近壁处催化复合反应对气动热环境的影响大小，包含完全催化、有限催化和完全非催化三种模型（李鹏等2 0 2 1a).完全催化和完全非催化条件是两种理论上的极限情况，有限催化条件最符合真实情况，研究最多.催化复合系数是有限催化模型的关键参数，目前主要采用两种方法确定：一种方法是在0，1 区间取常数值，另一种方法是按照基于壁温的模型函数取值.催化复合系数的半经验预测模型主要包括Stewart-RCG，Zo b y-R CG，Sc o t t-R CG,RCC-LVP，CCA T-A CC以及CSiC-SENCMA等（原志超2 0 17)，但并不完全准确，普遍认为催化复合系数与壁面温度、反应机理、结构材料以及表面光滑度等因素有关，因此更为精确的预测模型条件值得研究.近年来，国内外学者就壁面催化反应机理进行了大量的深入研究，提出了多种催化模型条件（Lietal.2015，杨肖峰等2 0 2 1,李芹等2 0 2 1).结构烧蚀效应通过烧蚀反应改变近壁处组元构成、温度分布，从而影响壁面热流分布。壁面结构烧蚀包括微烧蚀和严重烧蚀等情形，其中严重烧蚀会引起结构变形，进而影响流动形态，因此通常需要采用流-固耦合模拟技术进行研究（Lee&Kim2021)，这里仅简要提及.烧蚀边界条件基于烧蚀反应机理，在壁面组分质量守恒方程中增加质量引射条件进行计算.现有烧蚀混合物569高振勋，蒋崇超高烩非平衡流动数值模拟方法研究综述第3期模型包括JANAF升华碳组元模型（Chaseetal.1975)、高温空气/碳-碳复合材料烧蚀反应模型（董维中等2 0 15)、高温空气/碳-碳化硅复合材料烧蚀反应模型（董维中等2 0 10)、以及高温空气/碳-酚醛复合材料烧蚀反应模型（Park2001,Martinetal.2010）等.烧蚀产物/流场气体混合是烧蚀化学模型的研究重点，针对不同热防护材料，需要发展更多新的烧蚀混合物化学模型。稀薄滑移效应通常采用滑移边界条件以描述速度滑移、温度跳跃和热滑移等现象对高超飞行器气动热的影响（Lofthouse2008，丁明松2 0 19).常用的速度和温度滑移条件有经典的Max-well模型、Gokcen模型、Lockerby函数修正滑移边界模型(Lockerby2005)，以及Kogan简化模型（Vottaet al.2009）等，按照偏导数近似方法又分为一阶、二阶模型等.平均自由程是滑移边界条件的关键参数，它的计算方法主要分为基于硬球模型（HS）假设的方法、基于变径硬球模型(VHS）假设的方法以及基于变径软球模型（VSS）假设的方法，常用的基于热力学参数的计算方法本质是基于硬球模型假设的方法（沈青2 0 0 3.在热化学非平衡条件下，滑移边界还需考虑组分的滑移条件.组分滑移条件包括一阶滑移模型、二阶滑移模型等（丁明松2 0 19).壁面组分滑移与催化效应相互影响，两者的耦合计算较为困难，目前在公开文献中少见记载.3CFD数值方法研究进展高温气体热化学非平衡效应是一种高马赫数高温极端条件下的复杂多尺度、多物理效应耦合流动问题，其CFD数值方法研究面临精准、高效和稳定等方面的考验，需要解决以下主要问题：（1）热-化学耦合机理复杂，激波强间断问题易引起数值振荡，计算精度尤其是气动热预测精度对网格质量、数值格式敏感；（2）方程数量众多，热力学能量体系形式复杂，计算量巨大，计算耗时；(3）时间尺度分布范围广，最大和最小特征时间存在数量级差异（10 以上）（Bussing&Murman1988)，迭代矩阵刚性强，数值稳定性和收敛性较差.因此，下文首先介绍热化学非平衡流动主控方程及相关耦合方程，然后从CFD数值计算的共性问题出发分别对数值格式、高效求解技术、刚性处理方法等进行概述.3.1流动主控方程忽略端流效应，综合考虑完全气体模型、一温度模型、两温度模型、三温度模型以及多振动温度模型，完整热化学非平衡气体Navier-Stokes（N-S）方程的微分守恒形式（董维中19 9 6，王京盈2 0 17，郝佳傲2 0 18，李海燕2 0 0 7）为(p)(pui)=0(1)t8(pui)Ti(puiuj+dijp)=(2)taTv.km8(pE)8T8TEns(pHui）入k+入EpDshhsoj(3)Tiiot8jk=1S-18(pcs)pcsujPDs=1,2,.,n(4)一t力57020233年第53 卷展进学(pev,k)OTv.knsOcspev.kZpDev.kaj+Ov,k,k=1,2,.,m(5)tS=1(peE)8TEnsOcs入EpDseE(6)tS=1其中式(1)为总质量连续性方程，式(2)为动量方程，式(3)为总能量方程，式（4)为组分输运方程或组元的质量连续性方程，式（5）为振动能量方程，式（6)为电子能量方程.此外，完整的完全气体N-S方程在简化能量项的基础上，由式(1)、式(2)和式(3)紧耦合得到.利用高斯散度定理，微分形式的N-S方程可转换为积分守恒形式（李鹏等2 0 2 1a，2 0 2 1b)QdV+Fds=Fuds+SdV(7)t82式中，Q为守恒量向量，F和F，分别为对流项通量和黏性项通量，S表示源项.在不同气体模型和温度模型中，Q和S的具体形式如图5所示.质量连续性方程目前有三种构建方式（李鹏等2 0 2 1a).（1）保留总质量连续性方程，减少组分输运方程.当混合气体不包含离子组元时，n=ns-1（n s 为混合气体组元个数），结合质量守恒条件封闭方程；当混合气体包含离子组元时，n=ns-2,结合质量守恒和电荷数守恒条件封闭方程.这样处理的优点在于可较好地继承完全气体模型方程形式，方便功能拓展和编程，缺点在于相关非平衡偏导数形式较为复杂.(2）去除总质量连续性方程，保留全部组分输运方程，混合气体密度由质量守恒条件计算.该方法不易退化为完全气体模型形式，对底层模型兼容性较差.(3）同时采用总质量连续性方程和全部组分输运方程.该方法相关偏导数形式简单，编程容易，但相对计算量较大，同时方程超定，易受数值格式影响出现非物理值问题，需要额外步数检验物理量的合理性，对组元质量分数等作归一化处理，目前，热化学非平衡流动的CFD数值方法可以大体分为紧耦合形式和非耦合形式.紧耦合方式是联合上述式（1）(6)方程进行整体计算和时间推进的方法，而非耦合方式通常将式（1）(3)作为一个整体与其余方程分开交替计算，或者将源项从方程中分离出来单独作为一个整体与不含源项的方程分开交替计算（刘君等2 0 18).非耦合方式在消除计算刚性方面具有较大优势，但与真实的流动/热化学反应耦合机制存在差距，尤其在强热化学非平衡效应情况下，其计算误差和有效性需要严格考证.紧耦合方式遵循真实的多物理效应耦合机理，因此被数值软件广泛采用，但需要解决刚性问题，工程上对数值稳定性要求较高3.2多物理场耦合方程在气动热环境数值研究方面，多物理效应耦合模拟技术是国内外的发展趋势（董维中等2015，Zi b i t s k e r e t a l.2 0 2 2).主要应用研究包括热化学非平衡端流、高超声速磁流体力学控制、结构热防护、辐射加热等，涉及热化学非平衡等离子体流场/流模式方程耦合模拟、热化学非平衡等离子体流场/电磁场耦合模拟、热化学非平衡等离子体流场/结构温度场耦合模拟，以及热化学非平衡等离子体流场/辐射传输场耦合模拟等.这些耦合问题研究是热化学非平衡流动模571高振勋，蒋崇高超高烩非平衡流动数值模拟方法研究综述第3期多振动温度热化学非平衡气体三温度热化学非平衡气体两温度热化学P平衡气体或非平衡气体单温度化学puQ非平衡气体puPpw完全气体SPpupEPpupupcipupupwPC2SpupwpEpwpESTpCnsnsPpESpcipC2Tpeu.1pupcpC2：T2peu.22pC22：TpCnsnpwTpCnsnsTpeTpeu.mTpETpCnsnsT.peTpee。T.pee图5热力学温度模型的数组拓展与简化示意图（李鹏等2 0 2 1a)拟方法的延伸拓展，与热化学非平衡N-S方程数值模拟方法具有共性问题，因此以下仅从耦合方程形式等方面进行简单介绍.在热化学非平衡流模拟方面，主要数值方法包括DNS，R A NS和LES(Cummings2022，唐志共，张益荣等2 0 15,Caryetal.2021,刘朋欣等2 0 2 2).DNS方法主要用于基础理论研究，例如端流效应与热化学耦合机理研究(Cummings2022)、热化学非平衡端流边界层特性分析(Candler2019，刘朋欣等2 0 2 2）等.RANS和LES主要用于工程流模拟，当前RANS方法应用最多.RANS方法目前直接沿用由低速无化学反应流动发展并经可压缩性修正得到的方程模型进行求解，其在形式上与式（1）（7）相同，其中湍流封闭方程直接采用现有的一方程模式、两方程模式以及雷诺应力模式，如S-A，k-s，k-，k-SST 等，以及它们的各种修正模型（李海燕2 0 0 7,Cummings2022).尽管RANS方法在热化学非平衡流的工程模拟中成功取得了广泛应用，但其计算有效性值得商榷（Candler 2015).由于当前使用的RANS方程并不是严格由热化学非平衡N-S方程经过Favr平均推导，并经过端流/热化学耦合效应机理分析或实验数据封闭端流模式相关项得到的，其应用于高超热化学非平衡流动模拟并不十分合理（Cummings2022).因此湍流/热化学非平衡效应的耦合机理及物理建模是国内外CFD研究者最为关注的课题（Cum-mings 2022,Cary et al.2021).在高超声速磁流体力学（magnetohydrodynamic，M H D）控制方面，数值模拟方法包括全MHD方法和低磁雷诺数MHD方法（丁明松2 0 19).全MHD方法在方程形式上以非耦合方式分开计算热化学非平衡N-S方程（流动主控方程）与电磁场控制方程。其中流动主控方程中动量与能量方程需要增加电磁场带来的动量和能量项，即电磁作用源项；电磁场控制方程由Maxwell电磁方程组和广义欧姆定律构成.两组方程尽管在数学形式上分开计算，但在流场推进中表现为强耦合关系，因此计算刚性问题不可忽略.低磁雷诺数MHD方法在全MHD模型基础上，保持流力20233年第53卷572展进学动主控方程不变，但忽略感应磁场，电磁相互作用以源项形式出现在主控方程中，同时耦合电场感应方程得到低磁雷诺数MHD模型.MHD方法物理意义明确，与一般N-S方程形式类似，在N-S方程中使用的方法在低磁雷诺数MHD模型都能较好适用.MHD相关方法研究可参考丁明松(2019)的文章.在结构热防护研究方面，热化学非平衡流场/结构温度场耦合模拟目前多采用非耦合的方式构建计算方程组，并将表面效应作为一种边界条件分别作用于非平衡流场和结构温度场，结构温度场主要由热传导方程计算，必要时需考虑烧蚀问题.热化学非平衡流场/表面效应/结构温度场的耦合模拟是国内外学者最为关注的研究课题之一（Cummings 2022)，是数值模拟软件必须具备的能力，也是非平衡气动热、表面效应、耦合计算策略等数值模拟方法研究的重点（董维中等2 0 15).在辐射加热特性研究方面，热化学非平衡等离子流场/辐射传输场耦合模拟基于热化学非平衡N-S方程采用非耦合的方式联合辐射传输方程计算，在N-S总能量方程中增加辐射热流项，同时修改相应能量方程的源项以描述热力学/强电离化学反应的耦合效应（王京盈2 0 17).辐射传输方程的计算方法有光学薄/厚近似、球谐函数近似（PN近似）、两热流法、离散坐标法、有限体积法、蒙特卡洛模拟法、离散传递法等，具体可参考王京盈（2 0 17）的文章，这里不作赘述，3.3数值格式有限差分方法和有限体积方法是数值求解热化学非平衡N-S方程的两种主要CFD离散方法（董维中