CN112966401A - 一种热化学非平衡多级气体模型自适应算法 - Google Patents

一种热化学非平衡多级气体模型自适应算法 Download PDF

Info

Publication number
CN112966401A
CN112966401A CN202110513881.7A CN202110513881A CN112966401A CN 112966401 A CN112966401 A CN 112966401A CN 202110513881 A CN202110513881 A CN 202110513881A CN 112966401 A CN112966401 A CN 112966401A
Authority
CN
China
Prior art keywords
temperature
thermochemical
gas
nonequilibrium
flow
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Granted
Application number
CN202110513881.7A
Other languages
English (en)
Other versions
CN112966401B (zh
Inventor
李鹏
丁明松
陈坚强
赵钟
董维中
何先耀
江涛
梅杰
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Original Assignee
Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Computational Aerodynamics Institute of China Aerodynamics Research and Development Center filed Critical Computational Aerodynamics Institute of China Aerodynamics Research and Development Center
Priority to CN202110513881.7A priority Critical patent/CN112966401B/zh
Publication of CN112966401A publication Critical patent/CN112966401A/zh
Application granted granted Critical
Publication of CN112966401B publication Critical patent/CN112966401B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
    • G06F17/12—Simultaneous equations, e.g. systems of linear equations
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2113/00—Details relating to the application field
    • G06F2113/08—Fluids
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2119/00—Details relating to the type or aim of the analysis or the optimisation
    • G06F2119/14—Force analysis or force optimisation, e.g. static or dynamic forces

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • General Engineering & Computer Science (AREA)
  • Data Mining & Analysis (AREA)
  • Operations Research (AREA)
  • Software Systems (AREA)
  • Databases & Information Systems (AREA)
  • Algebra (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)

Abstract

本发明公开了一种热化学非平衡多级气体模型自适应算法，主要用于高超声速热化学非平衡流动的数值模拟过程。该算法基于高温气体热力学模态激发、离解/电离等热化学非平衡特征，从物理机制出发采用由简单到复杂的分段渐次逼近模拟策略，依据压力、离解度、温度和电离度等判据自动完成计算气体模型的演变和逐级流动模拟，实现高超声速热化学非平衡流动的自适应模拟过程。该算法适用于地球大气、火星大气和高温燃气等多种计算介质，应用范围广，计算稳定性好，在保证计算精准度的前提下能加快收敛速度和大幅提升计算效率。

Description

一种热化学非平衡多级气体模型自适应算法
技术领域
本发明涉及空气动力学领域，具体涉及高超声速热化学非平衡流动的数值计算方法。
背景技术
高超声速飞行器在大气中飞行，如果飞行速度较高（马赫数10以上），激波波后将产生数千至上万开尔文的高温，高温流场中将会发生复杂的物理化学现象，混合气体发生离解、复合、交换、电离、多级电离等化学反应，气体分子平转动能、振动能、电子能/电子束缚能等热力学模态激发，形成混合多组分离解-电离-热力学激发气体云，并沿流动方向向后扩展，在飞行器周围形成鞘套，这种现象称为高温气体效应。由于高温化学反应和热力学激发是在流动中进行的，常常处于非平衡态，因此又被称为高温气体热化学非平衡效应。高温气体热化学非平衡效应，会严重影响飞行器的气动力特性、气动热环境、等离子体特性、辐射特性等气动特性，是高超声速飞行器气动操控设计、气动热防护、电磁通讯设计的研究基础。由于其重要性和复杂程度，也一直是高超声速飞行技术研究的前沿课题。
数值模拟是高超声速热化学非平衡流动研究的主要手段之一。高超声速热化学非平衡流动模拟与一般流动模拟存在较大不同，首先是控制方程数目多且形式复杂，计算量巨大。在完全气体流动的数值计算中，控制方程一般最多五个；而热化学非平衡流控制方程可达十几个，精细模拟离解/电离/多次电离、热力学各能量模态多级激发时，方程的数目还会大幅的增加，方程中还包含复杂的热化学源项和组分扩散项以及各种耦合效应项，因此其计算量比完全气体大得多。其次是求解的刚性问题，在高超声速流动过程中，流场各区域流动速度和温度差异较大，流动的特征时间与热化学的特征时间（化学反应的特征时间和振动能量的松弛时间等）常常存在量级上的差别，控制方程的求解出现刚性问题，计算稳定性和收敛性差，计算效率低。
国内外针对高超声速热化学非平衡流动数值模拟方法开展了大量研究，常见的数值模拟方法主要可以分为两类。一类是非耦合（或松耦合）的方法，将流动模拟和热化学物理现象的模拟解耦，分别独立（或近似独立）求解流动控制方程和热化学机制表征方程，可以在一定程度上避开数值模拟的刚性问题，因此计算效率较高。但这种方法与真实的物理耦合机制存在一定差别，其有效性仍有待商榷。
另一类是耦合的计算方法，将流动控制方程和热化学机制表征方程作为整体进行求解，实现流动过程和热化学物理的耦合模拟。为了提高计算稳定性和效率，主要依赖于改造时间推进格式，如全隐式、点隐式和部分隐式等各种隐式计算方法。但这些方法并没有完全解决热化学非平衡流动模拟“计算量大”和“求解刚性”等问题，实际应用过程，仍经常出现“效率低”、“不收敛”“发散”等问题。
因此，仍有必要从热化学非平衡效应的物理机制出发，发展效率更高、稳定性更好的热化学非平衡流动数值模拟方法。
发明内容
本发明的目的在于提供一种热化学非平衡多级气体模型自适应算法，主要用于高超声速热化学非平衡流动的数值模拟过程。该算法基于高温气体热力学模态激发、离解/电离等热化学非平衡特征，从物理机制出发采用由简单到复杂的分段渐次逼近模拟策略，依据压力、离解度、温度和电离度等判据自动完成计算气体模型演变和逐级流动模拟，实现高超声速热化学非平衡流动的自适应模拟过程。
为了实现上述目的，本发明采用如下技术方案：
步骤一：在高超声速流动模拟过程中，首先忽略热化学非平衡效应，数值求解基于完全气体的流动控制方程组，直至满足压力判据获得稳定基础流场；
步骤二：在步骤一的基础上，采用单温度化学非平衡气体模型模拟高温气体离解、复合和置换反应机制，数值求解热力学平衡-化学非平衡的流动控制方程组，直至满足离解度判据使得流场主要化学反应过程模拟稳定；
步骤三：在步骤二的基础上，采用两温度热化学非平衡气体模型模拟高温气体热力学振动激发与松弛机制，数值求解振动非平衡-化学非平衡的流动控制方程组，直至满足温度判据使得流场热力学平转动能-振动能松弛过程模拟稳定；
步骤四：在步骤三的基础上，采用三温度热化学非平衡气体模型模拟高温气体电离反应机制和重粒子热力学电子束缚能激发效应，数值求解三温度热力学非平衡-化学非平衡的流动控制方程组，直至满足电离度判据使得流场主要电离过程模拟稳定；
步骤五：在步骤四的基础上，采用多温度能级激发、多种电离成分的热化学非平衡气体模型模拟热力学多能级/多能态激发效应和重粒子多级电离反应机制，数值求解多能级温度热力学非平衡-多电离成分化学非平衡的流动控制方程组，直至满足计算收敛精度要求，计算流场稳定。
综上所述，由于采用了上述技术方案，本发明的有益效果是：
本发明从热化学非平衡流动特征出发，由主到次地依次考虑流动中的物理化学作用机制，可实现从“流场激波结构生成”、“气体离解/复合”、“分子热力学初步激发”、“气体初步电离”、“气体进一步电离”以及“热力学多能级进一步激发”等过程的逐级逼近模拟，抓住了物理化学现象的主要影响因素，降低了状态参数在计算过程中极度偏离真实物理过程的可能性，因此计算稳定性好。
本发明所采用的分段渐次逼近模拟方法是一种由简单到复杂的模拟过程，兼顾了效率和精度。在数值模拟前期，使用的物理化学模型相对简单，能节约大量的计算时间，得到较好的阶段性流场，后续模拟过程基于前面流动模拟获得的良好流动参数赋初值，收敛所需迭代步数减少，总体计算量和耗时也将大幅降低，因此具有加快收敛速度和提升计算效率的作用。在数值计算的后期，考虑的物理化学机制较为全面，在一定程度上保证了计算的精准度。
本发明适用面广，对气体介质没有限定，可普遍适用于各类气体介质，如地球大气、火星大气、高温燃气等气体介质的高温热化学非平衡流动模拟。
附图说明
本发明将通过例子并参照附图的方式说明，其中：
图1为本方案的算法流程示意图；
图2为实例一OREX返回器外形采用本发明方法和一般方法的计算残差对比；
图3为实例二HEG圆柱实验模型采用本发明方法和一般方法的计算残差对比。
具体实施方式
本说明书中公开的所有特征，或公开的所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以以任何方式组合。
本说明书（包括任何附加权利要求、摘要和附图）中公开的任一特征，除非特别叙述，均可被其他等效或具有类似目的的替代特征加以替换。即，除非特别叙述，每个特征只是一系列等效或类似特征中的一个例子而已。
本发明的算法流程如图1所示，以“地球大气热化学非平衡流动模拟”为实例，说明其具体的实现过程：
步骤S1：在高超声速流动模拟过程中，首先忽略热化学非平衡效应，数值求解基于完全气体的流动控制方程组，直至满足压力判据C1获得稳定基础流场；
步骤S1的主要目的是快速得到基础流场，因此采用忽略热化学非平衡效应的气体模型进行初步流动模拟。忽略热化学非平衡效应的气体模型包括但不限于完全气体模型和平衡气体模型，这里从计算效率考虑主要数值求解基于完全气体的流动控制方程。完全气体流动控制方程组的守恒积分形式为：
Figure 172248DEST_PATH_IMAGE001
其中守恒量
Figure 283423DEST_PATH_IMAGE002
，源项
Figure 559684DEST_PATH_IMAGE003
。完全气体总内能
Figure 730902DEST_PATH_IMAGE004
，比内能
Figure 967980DEST_PATH_IMAGE005
，其定容比热
Figure 882846DEST_PATH_IMAGE006
为常数，t为时间
Figure 279192DEST_PATH_IMAGE007
代表流体微元控制体，V为流体微元控制体的体积，F为对流通量，
Figure 558995DEST_PATH_IMAGE008
为粘性通量，n为流通面的法向量，s（小写）为流通面的面积，S（大写加粗）表示源项，
Figure 673582DEST_PATH_IMAGE009
为来流雷诺数数值，
Figure 858051DEST_PATH_IMAGE010
为密度，u、v、w分别表示流动速度在坐标轴x、y、z方向上的分量速度，E为总内能，e为比内能，T为温度，模拟过程采用LUSGS时间推进，无源项计算。
在S1的数值模拟过程中，压力分布差异是流场结构的主要表征之一，因此选择特征区域的压强相对变化量建立压力判据C1，作为步骤S1流动模拟是否计算稳定的判别条件。压力判据C1的数学表达式为
Figure 515428DEST_PATH_IMAGE011
其中上标
Figure 825187DEST_PATH_IMAGE012
表示时间迭代步序号，下标
Figure 833594DEST_PATH_IMAGE013
表示驻点附近的控制体单元序号，
Figure 824684DEST_PATH_IMAGE014
为压强，阈值
Figure 664464DEST_PATH_IMAGE015
。
步骤S2：在步骤S1的基础上，采用单温度化学非平衡气体模型模拟高温气体离解、复合和置换反应机制，数值求解热力学平衡-化学非平衡的流动控制方程组，直至满足离解度判据C2使得流场主要化学反应过程模拟稳定；
地球大气的主要成分是氮气N2和氧气O2。对于静态的空气，在一个大气压条件下，当环境温度达到2500K左右时，O2首先发生离解反应；在4000K左右时O2几乎完全离解，同时N2开始离解；在9000K左右时N2几乎完全离解。因此，模拟热化学非平衡效应，首先应该考虑的就是随气体温度升高，而发生的离解、复合、置换等反应，即考虑N2，O2，NO，N和O等气体组分相关反应。这里基于热力学平衡-化学非平衡控制方程组，采用单温度5组分（N2，O2，NO，N和O）气体模型进行计算，以离解度判据C2作为离解反应过程模拟是否稳定的依据条件。当满足离解度判据C2时得到离解反应稳定的流场，同时自动切换至步骤S3继续计算过程。
热力学平衡-化学非平衡控制方程组的守恒积分形式见上式，其中守恒量
Figure 82807DEST_PATH_IMAGE016
，源项
Figure 578510DEST_PATH_IMAGE017
，
Figure 170029DEST_PATH_IMAGE018
表示第
Figure 801998DEST_PATH_IMAGE013
个组分的密度，
Figure 656822DEST_PATH_IMAGE019
为第
Figure 436559DEST_PATH_IMAGE013
个组分方程的化学反应生成源项。化学非平衡气体的比内能
Figure 766522DEST_PATH_IMAGE020
，其中
Figure 518577DEST_PATH_IMAGE021
为组分个数，
Figure 341040DEST_PATH_IMAGE022
分别表示第
Figure 545756DEST_PATH_IMAGE013
个组分的质量分数和比焓，
Figure 416760DEST_PATH_IMAGE023
为混合气体的分子量，R是普适气体常量。在化学非平衡气体模型中，气体比热
Figure 616797DEST_PATH_IMAGE024
和
Figure 813423DEST_PATH_IMAGE006
以及比焓
Figure 974277DEST_PATH_IMAGE025
等均是组分
Figure 976868DEST_PATH_IMAGE026
、温度T的函数，数值模拟过程采用LUSGS时间推进，源项采用点隐式计算方法。
在S2数值模拟的初始化阶段，采用步骤S1完全气体模拟获得的参数赋初值，对于增加的组分质量分数
Figure 703516DEST_PATH_IMAGE026
等变量使用（但不限于）以下两种方法：（1）由来流条件赋值；（2）基于完全气体计算结果由平衡气体模型估算的
Figure 539885DEST_PATH_IMAGE026
赋值。为了得到更好的阶段性初场以保证计算的稳定性，这里推荐第（2）种方法。
在S2的数值模拟过程中，驻点区域流动变化的主要趋势是氧气O2和氮气N2等主要分子组分发生离解反应，组分浓度分布随流动发生改变直至稳定，因此这里以氧气O2和氮气N2两种分子组分浓度的相对变化量建立离解度判据C2，作为步骤S2流动模拟是否计算稳定的判别条件。离解度判据C2的数学表达式为
Figure 47090DEST_PATH_IMAGE027
其中上标
Figure 358483DEST_PATH_IMAGE012
表示时间迭代步序号，下标
Figure 408479DEST_PATH_IMAGE013
表示驻点附近的控制体单元序号，
Figure 9225DEST_PATH_IMAGE028
、
Figure 941408DEST_PATH_IMAGE029
分别表示氧气O2和氮气N2的质量分数，阈值
Figure 161168DEST_PATH_IMAGE030
。
步骤S3：在步骤S2的基础上，采用两温度热化学非平衡气体模型模拟高温气体热力学振动激发与松弛机制，数值求解振动非平衡-化学非平衡的流动控制方程组，直至满足温度判据C3使得流场热力学平转动能-振动能松弛过程模拟稳定；
在S2模拟基础上进一步考虑多原子组分的振动能模态激发效应，此时需要用两个热力学温度分别描述组分的平转动能和振动能，这里基于两温度热化学非平衡控制方程组采用两温度5组分（N2，O2，NO，N和O）气体模型进行计算，以驻点温度判据C3作为热力学平转动能-振动能松弛过程是否稳定的依据条件。当满足驻点温度判据C3时得到平转动能-振动能松弛过程模拟稳定的流场，同时自动切换至步骤S4继续计算过程。
两温度热化学非平衡方程组的守恒积分形式见式（1），其中守恒量
Figure 862408DEST_PATH_IMAGE031
，源项
Figure 899634DEST_PATH_IMAGE032
，
Figure 53535DEST_PATH_IMAGE018
表示第
Figure 342565DEST_PATH_IMAGE013
个组分的密度，
Figure 226207DEST_PATH_IMAGE019
为第
Figure 309701DEST_PATH_IMAGE013
个组分方程的化学反应生成源项，
Figure 947968DEST_PATH_IMAGE033
为振动能，
Figure 634165DEST_PATH_IMAGE034
为振动能量方程源项。两温度5组分热化学非平衡气体的比内能
Figure 513259DEST_PATH_IMAGE035
,其中
Figure 829971DEST_PATH_IMAGE021
为组分个数，
Figure 286360DEST_PATH_IMAGE036
分别表示第
Figure 448351DEST_PATH_IMAGE013
个组分的质量分数、平转动能、振动能和电子能。在两温度热化学非平衡气体模型中，平转动能
Figure 916372DEST_PATH_IMAGE037
对应平转动温度
Figure 466303DEST_PATH_IMAGE038
，振动能
Figure 285354DEST_PATH_IMAGE033
和电子能
Figure 516615DEST_PATH_IMAGE039
采用振动温度
Figure 229356DEST_PATH_IMAGE040
计算，数值模拟过程采用LUSGS时间推进，源项采用点隐式计算方法。
在S3数值模拟的初始化阶段，采用S2模拟获得的参数赋初值，增加的振动能变量
Figure 559975DEST_PATH_IMAGE033
使用S2模拟结果中的温度
Figure 990956DEST_PATH_IMAGE041
和组分质量分数
Figure 757399DEST_PATH_IMAGE026
计算，即
Figure 200013DEST_PATH_IMAGE042
。
在S3数值模拟过程中，流动变化的主要趋势是过激波后振动能模态激发，沿流动方向平转动能-振动能松弛直至新的平衡，因此选择特征区域振动温度的相对变化量建立温度判据C3，作为步骤S3流动模拟是否计算稳定的判别条件。驻点温度判据C3的数学表达式为
Figure 91745DEST_PATH_IMAGE043
其中上标
Figure 682127DEST_PATH_IMAGE012
表示时间迭代步序号，下标
Figure 458453DEST_PATH_IMAGE013
表示驻点附近的控制体单元序号，
Figure 614628DEST_PATH_IMAGE040
为振动温度，阈值
Figure 818207DEST_PATH_IMAGE044
。
步骤S4：在步骤S3的基础上，采用三温度热化学非平衡气体模型模拟高温气体电离反应机制和重粒子热力学电子束缚能激发效应，数值求解三温度热力学非平衡-化学非平衡的流动控制方程组，直至满足电离度判据C4使得流场主要电离过程模拟稳定；
在S3模拟基础上进一步考虑重粒子的电子束缚能模态激发效应，需要用三个热力学温度分别描述重粒子组分的平转动能、振动能以及电子能，此时气体的电离过程不能忽略。一般认为，在温度低于8000K时，空气电离以NO+的相关反应为主，因此可考虑N2，O2，NO，N、O、NO+和e-等组分参与的离解、电离等反应。这里基于三温度热化学非平衡控制方程组采用三温度7组分（N2，O2，NO，N、O、NO+和e-）气体模型进行计算，以电离度判据C4作为平转动能-振动能/电子能松弛过程和主要电离反应过程是否稳定的依据条件。当满足电离度判据C4时得到平转动能-振动能/电子能松弛和主要电离反应模拟稳定的流场，同时自动切换至步骤S5继续计算过程。
三温度热化学非平衡控制方程组的守恒积分形式见式（1），其中守恒量
Figure 630305DEST_PATH_IMAGE045
，源项
Figure 69377DEST_PATH_IMAGE046
，
Figure 283321DEST_PATH_IMAGE018
表示第
Figure 657801DEST_PATH_IMAGE013
个组分的密度，
Figure 19512DEST_PATH_IMAGE019
为第
Figure 931449DEST_PATH_IMAGE013
个组分方程的化学反应生成源项，
Figure 203162DEST_PATH_IMAGE047
分别为振动能和电子能，
Figure 607598DEST_PATH_IMAGE048
分别为振动能量方程源项和电子能量方程源项。其比内能
Figure 597551DEST_PATH_IMAGE049
，其中
Figure 50529DEST_PATH_IMAGE050
为组分个数，
Figure 301382DEST_PATH_IMAGE036
分别表示第
Figure 17665DEST_PATH_IMAGE013
个组分的质量分数、平转动能、振动能和电子能。在三温度热化学非平衡气体模型中，平转动能
Figure 26072DEST_PATH_IMAGE037
对应平转动温度
Figure 282741DEST_PATH_IMAGE038
，振动能
Figure 653680DEST_PATH_IMAGE033
和电子能
Figure 275285DEST_PATH_IMAGE039
分别采用振动温度
Figure 98885DEST_PATH_IMAGE040
和电子温度
Figure 156315DEST_PATH_IMAGE051
计算，数值模拟过程采用LUSGS时间推进，源项采用点隐式计算方法。
在S4数值模拟的初始化阶段，采用步骤S3模拟获得的参数赋初值，新增离子组分的质量分数赋初值零，新增振动能
Figure 991547DEST_PATH_IMAGE033
和电子能
Figure 846370DEST_PATH_IMAGE039
变量则使用S3模拟结果中的振动温度
Figure 157266DEST_PATH_IMAGE040
和组分质量分数
Figure 959000DEST_PATH_IMAGE026
重新计算，即
Figure 711055DEST_PATH_IMAGE052
。
在S4数值模拟过程中，流动变化的主要趋势是重粒子束缚电子能模态激发，平转动能-振动能/电子能松弛直至新的平衡，同时进一步考虑主要原子组分发生电离反应，混合气体电离度产生变化直至稳定，因此这里以电离度的相对变化量建立电离度判据C4，作为步骤S4流动模拟是否计算稳定的判别条件。电离度判据C4的数学表达式为
Figure 64676DEST_PATH_IMAGE053
其中上标
Figure 738234DEST_PATH_IMAGE012
表示时间迭代步序号，下标
Figure 140397DEST_PATH_IMAGE013
表示驻点附近的控制体单元序号，
Figure 74855DEST_PATH_IMAGE054
表示混合气体的电离度，阈值
Figure 474743DEST_PATH_IMAGE055
。
步骤S5：在步骤S4的基础上，采用多温度能级激发、多种电离成分的热化学非平衡气体模型模拟热力学多能级/多能态激发效应和重粒子多级电离反应机制，数值求解多能级温度热力学非平衡-多电离成分化学非平衡的流动控制方程组，直至满足计算收敛精度要求，计算流场稳定。
在S4模拟基础上进一步考虑热力学多能态/多能级激发效应或重粒子的多级电离反应过程，此时需要引入多个热力学温度分别描述组分的平转动能、电子能以及其它能态对应的振动能，这里基于多振动温度模型热化学非平衡控制方程组，采用11组分（N2，O2，NO，N、O、O+、O2 +、N+、N2 +、NO+和e-）气体模型进行计算，依据计算收敛精度要求判断计算过程是否需要停止。当满足计算收敛精度要求时，完成计算过程获取流动参数。
多振动温度热化学非平衡控制方程组的守恒积分形式见式（1），其中守恒量
Figure 760231DEST_PATH_IMAGE056
，源项
Figure 697575DEST_PATH_IMAGE057
，
Figure 627485DEST_PATH_IMAGE018
表示第
Figure 995013DEST_PATH_IMAGE013
个组分的密度，
Figure 502217DEST_PATH_IMAGE019
为第
Figure 183866DEST_PATH_IMAGE013
个组分方程的化学反应生成源项，
Figure 358495DEST_PATH_IMAGE058
分别为振动模态
Figure 896924DEST_PATH_IMAGE059
对应的振动能和振动能量方程源项，
Figure 32370DEST_PATH_IMAGE060
分别为电子能和电子能量方程源项。在多振动温度模型中，比内能
Figure 376764DEST_PATH_IMAGE061
，其中
Figure 78003DEST_PATH_IMAGE062
为组分个数，
Figure 990596DEST_PATH_IMAGE063
表示组分
Figure 472393DEST_PATH_IMAGE013
的振动模态个数，
Figure 758493DEST_PATH_IMAGE064
分别表示第
Figure 314239DEST_PATH_IMAGE013
个组分的质量分数、平转动能和电子能，
Figure 522367DEST_PATH_IMAGE065
表示组分
Figure 366826DEST_PATH_IMAGE013
第
Figure 990705DEST_PATH_IMAGE059
个振动模态对应的振动能。因此平转动能
Figure 728854DEST_PATH_IMAGE037
对应平转动温度
Figure 983249DEST_PATH_IMAGE038
，电子能
Figure 111742DEST_PATH_IMAGE039
对应电子温度
Figure 601629DEST_PATH_IMAGE051
，振动能
Figure 335230DEST_PATH_IMAGE066
由多个振动温度
Figure 822843DEST_PATH_IMAGE067
对应描述和计算。数值模拟过程采用LUSGS时间推进，源项采用点隐式计算方法。
计算收敛精度要求包括但不限于残差条件和总迭代步数要求等，根据自身程序设计决定。这里给出一种采用残差条件的计算收敛判据，即
Figure 766528DEST_PATH_IMAGE068
,
其中：
Figure 666964DEST_PATH_IMAGE069
表示时间迭代步序号，
Figure 51809DEST_PATH_IMAGE070
表示全流场控制体单元序号，
Figure 38219DEST_PATH_IMAGE071
为计算残差值，阈值
Figure 78988DEST_PATH_IMAGE072
。
需要特别说明的是，上述分段渐近逼近模拟过程根据计算的实际需要可以整合或删减某些中间步骤，例如前述模拟过程可以删除步骤S4，当S3流动模拟结束后直接切换至步骤S5。再例如对某个工况的模拟只需要用到两温度7组分气体时，可以删除步骤S5，将步骤S4调整为采用两温度7组分气体模拟。
实施例一：OREX返回器高超声速再入流动模拟。计算模型和飞行试验条件参考文献“Yukimitsu Y., Minako Y., CFD and FEM Coupling Analysis of OREXAerothermodynamic Flight Data, AIAA 95-2087, 29th AIAA ThermophysicsConference, June 19-22, 1995/San Diego, CALIFORNIA, 1995:1-14”。计算模拟状态为飞行高度59.6km，来流速度5561.6m/s，来流温度248.12K，壁面温度1519K。数值模拟采用Park化学模型，对流通量选取AUSMDV格式和Minmod限制器，CFL数取10，壁面条件为等温非催化壁。
如图2所示，基于实例一给出了本发明方法与现有方法的计算残差对比，横坐标表示计算时间（twall）单位秒（s），N2为氮气，纵坐标表示平均计算残差（无单位）。现有方法为不采用渐进逼近方式，直接求解完整的热化学非平衡流动控制方程的方法。结果显示，采用一般方法，计算结果发散，而采用本发明，却能较好地收敛。可见，本发明能提升计算稳定性。
实施例二：HEG激波风洞圆柱模型实验流动模拟。计算模型和实验条件参考文献“Hannemann K., High Enthalpy Flows in the HEG Shock Tunnel: Experiment andNumerical Rebuilding, AIAA 2003-0978, 41st Aerospace Sciences Meeting andExhibit 6-9 January 2003, Reno, Nevada, 2003:1-18”。数值模拟采用Dunn-Kang化学模型，对流通量选取AUSMDV格式和Minmod限制器，CFL数取10，壁面条件为等温非催化壁（壁面温度300K）。
如图3所示，基于实例二给出了本发明方法与现有方法的计算残差对比，横坐标表示迭代步数（无单位），纵坐标表示平均计算残差（无单位），图中Tv(K)为振动温度。结果显示，本发明方法计算收敛时的站位较一般方法提前约150s，计算耗时减少约32%。可以看出，本发明具有加快收敛速度和提升计算效率的作用。
本发明并不局限于前述的具体实施方式。本发明扩展到任何在本说明书中披露的新特征或任何新的组合，以及披露的任一新的方法或过程的步骤或任何新的组合。

Claims (2)

1.一种热化学非平衡多级气体模型自适应算法，其特征在于基于高温气体热力学模态激发、离解/电离等热化学非平衡特征，从物理机制出发采用由简单到复杂的分段渐次逼近模拟策略，依据压力、离解度、温度和电离度等判据自动完成计算气体模型的演变和逐级流动模拟，实现高超声速热化学非平衡流动的自适应模拟过程，包括以下步骤：
S1：在高超声速流动模拟过程中，忽略热化学非平衡效应，数值求解基于完全气体的流动控制方程组，直至流场特征区域压强变化相对量的最大值小于阈值
Figure 269334DEST_PATH_IMAGE001
时，基础流场稳定，所述
Figure 276604DEST_PATH_IMAGE002
；
S2：在S1的基础上，采用单温度化学非平衡气体模型模拟高温气体离解、复合和置换反应机制，数值求解热力学平衡-化学非平衡的流动控制方程组，直至流场特征区域主要离解组分质量分数变化相对量的最大值小于阈值
Figure 536684DEST_PATH_IMAGE003
时，流场主要化学反应过程模拟稳定，所述
Figure 287603DEST_PATH_IMAGE004
；
S3：在S2的基础上，采用两温度热化学非平衡气体模型模拟高温气体热力学振动激发与松弛机制，数值求解振动非平衡-化学非平衡的流动控制方程组，直至流场特征区域振动温度变化相对量的最大值小于阈值
Figure 457684DEST_PATH_IMAGE005
时，流场热力学平转动能-振动能松弛过程模拟稳定，所述
Figure 399970DEST_PATH_IMAGE006
；
S4：在S3的基础上，采用三温度热化学非平衡气体模型模拟高温气体电离反应机制和重粒子热力学电子束缚能激发效应，数值求解三温度热力学非平衡-化学非平衡的流动控制方程组，直至特征区域电离度变化相对量的最大值小于阈值
Figure 147346DEST_PATH_IMAGE007
时，流场主要电离过程模拟稳定，所述
Figure 701955DEST_PATH_IMAGE008
；
S5：在S4的基础上，采用多温度能级激发、多种电离成分的热化学非平衡气体模型模拟热力学多能级/多能态激发效应和重粒子多级电离反应机制，数值求解多能级温度热力学非平衡-多电离成分化学非平衡的流动控制方程组，直至流场残差变化相对量的最大值小于阈值
Figure 726543DEST_PATH_IMAGE009
时，满足计算收敛所需精度要求，所述
Figure 341195DEST_PATH_IMAGE011
。
2.根据权利要求1所述的一种热化学非平衡多级气体模型自适应算法，其特征在于在S1-S5中的阀值比较式分别为压力、离解度、温度和电离度判据，它们的数学表达形式为：
Figure 310288DEST_PATH_IMAGE012
，
Figure 901544DEST_PATH_IMAGE013
，
Figure 46218DEST_PATH_IMAGE014
，
Figure 956405DEST_PATH_IMAGE015
,
Figure 553740DEST_PATH_IMAGE016
其中：
Figure 450151DEST_PATH_IMAGE017
表示时间迭代步序号
Figure 682287DEST_PATH_IMAGE018
表示流场特征区域的控制体单元序号，
Figure 763376DEST_PATH_IMAGE019
为压强，
Figure 848007DEST_PATH_IMAGE020
表示混合气体主要离解组分，
Figure 548109DEST_PATH_IMAGE021
为质量分数，
Figure 260850DEST_PATH_IMAGE022
为振动温度,
Figure 122627DEST_PATH_IMAGE023
为混合气体电离度，
Figure 198948DEST_PATH_IMAGE024
为计算残差值。
CN202110513881.7A 2021-05-12 2021-05-12 一种热化学非平衡多级气体模型自适应方法 Active CN112966401B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110513881.7A CN112966401B (zh) 2021-05-12 2021-05-12 一种热化学非平衡多级气体模型自适应方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110513881.7A CN112966401B (zh) 2021-05-12 2021-05-12 一种热化学非平衡多级气体模型自适应方法

Publications (2)

Publication Number Publication Date
CN112966401A true CN112966401A (zh) 2021-06-15
CN112966401B CN112966401B (zh) 2021-07-16

Family

ID=76279723

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110513881.7A Active CN112966401B (zh) 2021-05-12 2021-05-12 一种热化学非平衡多级气体模型自适应方法

Country Status (1)

Country Link
CN (1) CN112966401B (zh)

Cited By (10)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113887063A (zh) * 2021-10-14 2022-01-04 圣泰(北京)工程软件有限公司 特殊相平衡多组分烃类的分离塔模拟计算方法
CN114491962A (zh) * 2021-12-30 2022-05-13 中国航天空气动力技术研究院 高温非平衡双原子气体热物性计算方法及数据库生成方法
CN114880971A (zh) * 2022-07-13 2022-08-09 中国空气动力研究与发展中心计算空气动力研究所 一种计算流体力学软件采用的隐式方法
CN115048775A (zh) * 2022-05-27 2022-09-13 中国空气动力研究与发展中心计算空气动力研究所 一种热化学非平衡流动的组分限制方法
CN115312139A (zh) * 2022-09-23 2022-11-08 中国空气动力研究与发展中心计算空气动力研究所 一种高超声速流动化学反应模型数据的存取与转换方法
CN116090262A (zh) * 2023-04-10 2023-05-09 中国空气动力研究与发展中心计算空气动力研究所 一种有限催化分级隐式数值模拟方法、装置、设备及介质
CN116227388A (zh) * 2023-04-21 2023-06-06 中国空气动力研究与发展中心计算空气动力研究所 高超流动模拟cfl数动态调整方法、系统、设备及介质
CN117059188A (zh) * 2023-10-12 2023-11-14 中国空气动力研究与发展中心计算空气动力研究所 一种化学非平衡气体热力学平衡能量体系改进方法及系统
CN117332511A (zh) * 2023-12-01 2024-01-02 中国空气动力研究与发展中心计算空气动力研究所 面向高超飞行器高温非平衡流的自适应耦合数值模拟方法
CN117393062A (zh) * 2023-12-13 2024-01-12 上海交通大学四川研究院 刚性化学反应流回退自适应半隐半显耦合时间的模拟方法

Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2017014474A (ja) * 2015-07-06 2017-01-19 Ｇｇｉジャパン株式会社 連続式熱化学型バイオマス原料ガス化装置
CN107942309A (zh) * 2017-10-19 2018-04-20 上海无线电设备研究所 一种稀薄大气层内超高速目标电磁散射快速计算方法
CN107992684A (zh) * 2017-12-05 2018-05-04 上海无线电设备研究所 一种时变等离子体等效分层介质模型建模方法
CN109581340A (zh) * 2018-10-11 2019-04-05 上海无线电设备研究所 一种基于时域弹跳射线法的等离子体电磁散射建模方法
CN111222241A (zh) * 2020-01-06 2020-06-02 中国人民解放军国防科技大学 热化学非平衡条件下流场数据的数值计算方法和装置
CN111595490A (zh) * 2020-04-07 2020-08-28 北京空天技术研究所 一种高温气体效应对气动热影响的测量方法
CN112417743A (zh) * 2021-01-25 2021-02-26 中国空气动力研究与发展中心计算空气动力研究所 一种气体能量反演热力学温度的混合迭代方法

Patent Citations (7)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2017014474A (ja) * 2015-07-06 2017-01-19 Ｇｇｉジャパン株式会社 連続式熱化学型バイオマス原料ガス化装置
CN107942309A (zh) * 2017-10-19 2018-04-20 上海无线电设备研究所 一种稀薄大气层内超高速目标电磁散射快速计算方法
CN107992684A (zh) * 2017-12-05 2018-05-04 上海无线电设备研究所 一种时变等离子体等效分层介质模型建模方法
CN109581340A (zh) * 2018-10-11 2019-04-05 上海无线电设备研究所 一种基于时域弹跳射线法的等离子体电磁散射建模方法
CN111222241A (zh) * 2020-01-06 2020-06-02 中国人民解放军国防科技大学 热化学非平衡条件下流场数据的数值计算方法和装置
CN111595490A (zh) * 2020-04-07 2020-08-28 北京空天技术研究所 一种高温气体效应对气动热影响的测量方法
CN112417743A (zh) * 2021-01-25 2021-02-26 中国空气动力研究与发展中心计算空气动力研究所 一种气体能量反演热力学温度的混合迭代方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
梁日辉: "高速飞行器高温流场及其辐射特性数值模拟", 《中国优秀硕士学位论文全文数据库 工程科技Ⅱ辑》 *

Cited By (18)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113887063A (zh) * 2021-10-14 2022-01-04 圣泰(北京)工程软件有限公司 特殊相平衡多组分烃类的分离塔模拟计算方法
CN113887063B (zh) * 2021-10-14 2024-04-26 圣泰(北京)软件科技有限公司 特殊相平衡多组分烃类的分离塔模拟计算方法
CN114491962A (zh) * 2021-12-30 2022-05-13 中国航天空气动力技术研究院 高温非平衡双原子气体热物性计算方法及数据库生成方法
CN115048775B (zh) * 2022-05-27 2024-04-09 中国空气动力研究与发展中心计算空气动力研究所 一种热化学非平衡流动的组分限制方法
CN115048775A (zh) * 2022-05-27 2022-09-13 中国空气动力研究与发展中心计算空气动力研究所 一种热化学非平衡流动的组分限制方法
CN114880971A (zh) * 2022-07-13 2022-08-09 中国空气动力研究与发展中心计算空气动力研究所 一种计算流体力学软件采用的隐式方法
CN114880971B (zh) * 2022-07-13 2022-09-20 中国空气动力研究与发展中心计算空气动力研究所 一种计算流体力学软件采用的隐式方法
CN115312139A (zh) * 2022-09-23 2022-11-08 中国空气动力研究与发展中心计算空气动力研究所 一种高超声速流动化学反应模型数据的存取与转换方法
CN115312139B (zh) * 2022-09-23 2023-01-13 中国空气动力研究与发展中心计算空气动力研究所 一种高超声速流动化学反应模型数据的存取与转换方法
CN116090262A (zh) * 2023-04-10 2023-05-09 中国空气动力研究与发展中心计算空气动力研究所 一种有限催化分级隐式数值模拟方法、装置、设备及介质
CN116227388A (zh) * 2023-04-21 2023-06-06 中国空气动力研究与发展中心计算空气动力研究所 高超流动模拟cfl数动态调整方法、系统、设备及介质
CN116227388B (zh) * 2023-04-21 2023-07-07 中国空气动力研究与发展中心计算空气动力研究所 高超流动模拟cfl数动态调整方法、系统、设备及介质
CN117059188B (zh) * 2023-10-12 2024-01-23 中国空气动力研究与发展中心计算空气动力研究所 一种化学非平衡气体热力学平衡能量体系改进方法及系统
CN117059188A (zh) * 2023-10-12 2023-11-14 中国空气动力研究与发展中心计算空气动力研究所 一种化学非平衡气体热力学平衡能量体系改进方法及系统
CN117332511A (zh) * 2023-12-01 2024-01-02 中国空气动力研究与发展中心计算空气动力研究所 面向高超飞行器高温非平衡流的自适应耦合数值模拟方法
CN117332511B (zh) * 2023-12-01 2024-02-20 中国空气动力研究与发展中心计算空气动力研究所 面向高超飞行器高温非平衡流的自适应耦合数值模拟方法
CN117393062A (zh) * 2023-12-13 2024-01-12 上海交通大学四川研究院 刚性化学反应流回退自适应半隐半显耦合时间的模拟方法
CN117393062B (zh) * 2023-12-13 2024-02-23 上海交通大学四川研究院 刚性化学反应流回退自适应半隐半显耦合时间的模拟方法

Also Published As

Publication number Publication date
CN112966401B (zh) 2021-07-16

Similar Documents

Publication Publication Date Title
CN112966401B (zh) 一种热化学非平衡多级气体模型自适应方法
Wang et al. Unified Navier-Stokes flowfield and performance analysis of liquid rocket engines
Anand et al. An advanced unstructured-grid finite-volume design system for gas turbine combustion analysis
CN111222241A (zh) 热化学非平衡条件下流场数据的数值计算方法和装置
Harris et al. Pulse detonation engine as a ramjet replacement
Yang et al. A mechanism of combustion mode transition for hydrogen fueled scramjet
CN116090110B (zh) 高超声速飞行器高温流场数值仿真的修正方法及相关组件
CN117059188B (zh) 一种化学非平衡气体热力学平衡能量体系改进方法及系统
Li et al. A Unfied computational formulation for multi-component and multi-phase flows
CN113792432A (zh) 基于改进型fvm-lbfs方法的流场计算方法
Cross et al. Conjugate analyses of ablation in rocket nozzles
Miró et al. Ionization and dissociation effects on boundary-layer stability
Longo et al. The challenge of modeling high speed flows
CN115169253A (zh) 一种热力学两温度模型非平衡能量体系改进方法
Cross et al. Conjugate analysis of rocket nozzle ablation
Munafo' et al. Vibrational state to state kinetics in expanding and compressing nitrogen flows
Ashpis et al. DBD plasma actuators for flow control in air vehicles and jet engines-simulation of flight conditions in test chambers by density matching
BOYD Assessment of chemical nonequilibrium in rarefied hypersonic flow
Terzic et al. Prediction of internal ballistic parameters of solid propellant rocket motors
Zeng et al. Computational simulation of reentry flows over hypersonic vehicles using nonlinear coupled constitutive relations
Grover et al. Coupled rotational-vibrational excitation in shock waves using trajectory-based direct simulation Monte Carlo
Ground et al. Computational analysis and characterization of the UTA 1.6 MW arc-heated wind tunnel facility
Zhang et al. Development of the simulation & experiment system for martian atmosphere simulation environment box
Shevyrin et al. Coupled kinetic-continuum modeling of reentry vehicle plasma environment
George et al. Numerical simulation of self-ignition of hydrogen-hydrocarbons mixtures in a hot supersonic air flow

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant