一种热化学非平衡多级气体模型自适应算法的制作方法

本发明涉及空气动力学领域，具体涉及高超声速热化学非平衡流动的数值计算方法。

背景技术：

高超声速飞行器在大气中飞行，如果飞行速度较高（马赫数10以上），激波波后将产生数千至上万开尔文的高温，高温流场中将会发生复杂的物理化学现象，混合气体发生离解、复合、交换、电离、多级电离等化学反应，气体分子平转动能、振动能、电子能/电子束缚能等热力学模态激发，形成混合多组分离解-电离-热力学激发气体云，并沿流动方向向后扩展，在飞行器周围形成鞘套，这种现象称为高温气体效应。由于高温化学反应和热力学激发是在流动中进行的，常常处于非平衡态，因此又被称为高温气体热化学非平衡效应。高温气体热化学非平衡效应，会严重影响飞行器的气动力特性、气动热环境、等离子体特性、辐射特性等气动特性，是高超声速飞行器气动操控设计、气动热防护、电磁通讯设计的研究基础。由于其重要性和复杂程度，也一直是高超声速飞行技术研究的前沿课题。

数值模拟是高超声速热化学非平衡流动研究的主要手段之一。高超声速热化学非平衡流动模拟与一般流动模拟存在较大不同，首先是控制方程数目多且形式复杂，计算量巨大。在完全气体流动的数值计算中，控制方程一般最多五个；而热化学非平衡流控制方程可达十几个，精细模拟离解/电离/多次电离、热力学各能量模态多级激发时，方程的数目还会大幅的增加，方程中还包含复杂的热化学源项和组分扩散项以及各种耦合效应项，因此其计算量比完全气体大得多。其次是求解的刚性问题，在高超声速流动过程中，流场各区域流动速度和温度差异较大，流动的特征时间与热化学的特征时间（化学反应的特征时间和振动能量的松弛时间等）常常存在量级上的差别，控制方程的求解出现刚性问题，计算稳定性和收敛性差，计算效率低。

国内外针对高超声速热化学非平衡流动数值模拟方法开展了大量研究，常见的数值模拟方法主要可以分为两类。一类是非耦合（或松耦合）的方法，将流动模拟和热化学物理现象的模拟解耦，分别独立（或近似独立）求解流动控制方程和热化学机制表征方程，可以在一定程度上避开数值模拟的刚性问题，因此计算效率较高。但这种方法与真实的物理耦合机制存在一定差别，其有效性仍有待商榷。

另一类是耦合的计算方法，将流动控制方程和热化学机制表征方程作为整体进行求解，实现流动过程和热化学物理的耦合模拟。为了提高计算稳定性和效率，主要依赖于改造时间推进格式，如全隐式、点隐式和部分隐式等各种隐式计算方法。但这些方法并没有完全解决热化学非平衡流动模拟“计算量大”和“求解刚性”等问题，实际应用过程，仍经常出现“效率低”、“不收敛”“发散”等问题。

因此，仍有必要从热化学非平衡效应的物理机制出发，发展效率更高、稳定性更好的热化学非平衡流动数值模拟方法。

技术实现要素：

本发明的目的在于提供一种热化学非平衡多级气体模型自适应算法，主要用于高超声速热化学非平衡流动的数值模拟过程。该算法基于高温气体热力学模态激发、离解/电离等热化学非平衡特征，从物理机制出发采用由简单到复杂的分段渐次逼近模拟策略，依据压力、离解度、温度和电离度等判据自动完成计算气体模型演变和逐级流动模拟，实现高超声速热化学非平衡流动的自适应模拟过程。

为了实现上述目的，本发明采用如下技术方案：

步骤一：在高超声速流动模拟过程中，首先忽略热化学非平衡效应，数值求解基于完全气体的流动控制方程组，直至满足压力判据获得稳定基础流场；

步骤二：在步骤一的基础上，采用单温度化学非平衡气体模型模拟高温气体离解、复合和置换反应机制，数值求解热力学平衡-化学非平衡的流动控制方程组，直至满足离解度判据使得流场主要化学反应过程模拟稳定；

步骤三：在步骤二的基础上，采用两温度热化学非平衡气体模型模拟高温气体热力学振动激发与松弛机制，数值求解振动非平衡-化学非平衡的流动控制方程组，直至满足温度判据使得流场热力学平转动能-振动能松弛过程模拟稳定；

步骤四：在步骤三的基础上，采用三温度热化学非平衡气体模型模拟高温气体电离反应机制和重粒子热力学电子束缚能激发效应，数值求解三温度热力学非平衡-化学非平衡的流动控制方程组，直至满足电离度判据使得流场主要电离过程模拟稳定；

步骤五：在步骤四的基础上，采用多温度能级激发、多种电离成分的热化学非平衡气体模型模拟热力学多能级/多能态激发效应和重粒子多级电离反应机制，数值求解多能级温度热力学非平衡-多电离成分化学非平衡的流动控制方程组，直至满足计算收敛精度要求，计算流场稳定。

综上所述，由于采用了上述技术方案，本发明的有益效果是：

本发明从热化学非平衡流动特征出发，由主到次地依次考虑流动中的物理化学作用机制，可实现从“流场激波结构生成”、“气体离解/复合”、“分子热力学初步激发”、“气体初步电离”、“气体进一步电离”以及“热力学多能级进一步激发”等过程的逐级逼近模拟，抓住了物理化学现象的主要影响因素，降低了状态参数在计算过程中极度偏离真实物理过程的可能性，因此计算稳定性好。

本发明所采用的分段渐次逼近模拟方法是一种由简单到复杂的模拟过程，兼顾了效率和精度。在数值模拟前期，使用的物理化学模型相对简单，能节约大量的计算时间，得到较好的阶段性流场，后续模拟过程基于前面流动模拟获得的良好流动参数赋初值，收敛所需迭代步数减少，总体计算量和耗时也将大幅降低，因此具有加快收敛速度和提升计算效率的作用。在数值计算的后期，考虑的物理化学机制较为全面，在一定程度上保证了计算的精准度。

本发明适用面广，对气体介质没有限定，可普遍适用于各类气体介质，如地球大气、火星大气、高温燃气等气体介质的高温热化学非平衡流动模拟。

附图说明

本发明将通过例子并参照附图的方式说明，其中：

图1为本方案的算法流程示意图；

图2为实例一orex返回器外形采用本发明方法和一般方法的计算残差对比；

图3为实例二heg圆柱实验模型采用本发明方法和一般方法的计算残差对比。

具体实施方式

本说明书中公开的所有特征，或公开的所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以以任何方式组合。

本说明书（包括任何附加权利要求、摘要和附图）中公开的任一特征，除非特别叙述，均可被其他等效或具有类似目的的替代特征加以替换。即，除非特别叙述，每个特征只是一系列等效或类似特征中的一个例子而已。

本发明的算法流程如图1所示，以“地球大气热化学非平衡流动模拟”为实例，说明其具体的实现过程：

步骤s1：在高超声速流动模拟过程中，首先忽略热化学非平衡效应，数值求解基于完全气体的流动控制方程组，直至满足压力判据c1获得稳定基础流场；

步骤s1的主要目的是快速得到基础流场，因此采用忽略热化学非平衡效应的气体模型进行初步流动模拟。忽略热化学非平衡效应的气体模型包括但不限于完全气体模型和平衡气体模型，这里从计算效率考虑主要数值求解基于完全气体的流动控制方程。完全气体流动控制方程组的守恒积分形式为：

其中守恒量，源项。完全气体总内能，比内能，其定容比热为常数，t为时间代表流体微元控制体，v为流体微元控制体的体积，f为对流通量，为粘性通量，n为流通面的法向量，s（小写）为流通面的面积，s（大写加粗）表示源项，为来流雷诺数数值，为密度，u、v、w分别表示流动速度在坐标轴x、y、z方向上的分量速度，e为总内能，e为比内能，t为温度，模拟过程采用lusgs时间推进，无源项计算。

在s1的数值模拟过程中，压力分布差异是流场结构的主要表征之一，因此选择特征区域的压强相对变化量建立压力判据c1，作为步骤s1流动模拟是否计算稳定的判别条件。压力判据c1的数学表达式为

其中上标表示时间迭代步序号，下标表示驻点附近的控制体单元序号，为压强，阈值。

步骤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继续计算过程。

热力学平衡-化学非平衡控制方程组的守恒积分形式见上式，其中守恒量，源项，表示第个组分的密度，为第个组分方程的化学反应生成源项。化学非平衡气体的比内能，其中为组分个数，分别表示第个组分的质量分数和比焓，为混合气体的分子量，r是普适气体常量。在化学非平衡气体模型中，气体比热和以及比焓等均是组分、温度t的函数，数值模拟过程采用lusgs时间推进，源项采用点隐式计算方法。

在s2数值模拟的初始化阶段，采用步骤s1完全气体模拟获得的参数赋初值，对于增加的组分质量分数等变量使用（但不限于）以下两种方法：（1）由来流条件赋值；（2）基于完全气体计算结果由平衡气体模型估算的赋值。为了得到更好的阶段性初场以保证计算的稳定性，这里推荐第（2）种方法。

在s2的数值模拟过程中，驻点区域流动变化的主要趋势是氧气o2和氮气n2等主要分子组分发生离解反应，组分浓度分布随流动发生改变直至稳定，因此这里以氧气o2和氮气n2两种分子组分浓度的相对变化量建立离解度判据c2，作为步骤s2流动模拟是否计算稳定的判别条件。离解度判据c2的数学表达式为

其中上标表示时间迭代步序号，下标表示驻点附近的控制体单元序号，、分别表示氧气o2和氮气n2的质量分数，阈值。

步骤s3：在步骤s2的基础上，采用两温度热化学非平衡气体模型模拟高温气体热力学振动激发与松弛机制，数值求解振动非平衡-化学非平衡的流动控制方程组，直至满足温度判据c3使得流场热力学平转动能-振动能松弛过程模拟稳定；

在s2模拟基础上进一步考虑多原子组分的振动能模态激发效应，此时需要用两个热力学温度分别描述组分的平转动能和振动能，这里基于两温度热化学非平衡控制方程组采用两温度5组分（n2，o2，no，n和o）气体模型进行计算，以驻点温度判据c3作为热力学平转动能-振动能松弛过程是否稳定的依据条件。当满足驻点温度判据c3时得到平转动能-振动能松弛过程模拟稳定的流场，同时自动切换至步骤s4继续计算过程。

两温度热化学非平衡方程组的守恒积分形式见式（1），其中守恒量，源项，表示第个组分的密度，为第个组分方程的化学反应生成源项，为振动能，为振动能量方程源项。两温度5组分热化学非平衡气体的比内能,其中为组分个数，分别表示第个组分的质量分数、平转动能、振动能和电子能。在两温度热化学非平衡气体模型中，平转动能对应平转动温度，振动能和电子能采用振动温度计算，数值模拟过程采用lusgs时间推进，源项采用点隐式计算方法。

在s3数值模拟的初始化阶段，采用s2模拟获得的参数赋初值，增加的振动能变量使用s2模拟结果中的温度和组分质量分数计算，即。

在s3数值模拟过程中，流动变化的主要趋势是过激波后振动能模态激发，沿流动方向平转动能-振动能松弛直至新的平衡，因此选择特征区域振动温度的相对变化量建立温度判据c3，作为步骤s3流动模拟是否计算稳定的判别条件。驻点温度判据c3的数学表达式为

其中上标表示时间迭代步序号，下标表示驻点附近的控制体单元序号，为振动温度，阈值。

步骤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），其中守恒量，源项，表示第个组分的密度，为第个组分方程的化学反应生成源项，分别为振动能和电子能，分别为振动能量方程源项和电子能量方程源项。其比内能，其中为组分个数，分别表示第个组分的质量分数、平转动能、振动能和电子能。在三温度热化学非平衡气体模型中，平转动能对应平转动温度，振动能和电子能分别采用振动温度和电子温度计算，数值模拟过程采用lusgs时间推进，源项采用点隐式计算方法。

在s4数值模拟的初始化阶段，采用步骤s3模拟获得的参数赋初值，新增离子组分的质量分数赋初值零，新增振动能和电子能变量则使用s3模拟结果中的振动温度和组分质量分数重新计算，即。

在s4数值模拟过程中，流动变化的主要趋势是重粒子束缚电子能模态激发，平转动能-振动能/电子能松弛直至新的平衡，同时进一步考虑主要原子组分发生电离反应，混合气体电离度产生变化直至稳定，因此这里以电离度的相对变化量建立电离度判据c4，作为步骤s4流动模拟是否计算稳定的判别条件。电离度判据c4的数学表达式为

其中上标表示时间迭代步序号，下标表示驻点附近的控制体单元序号，表示混合气体的电离度，阈值。

步骤s5：在步骤s4的基础上，采用多温度能级激发、多种电离成分的热化学非平衡气体模型模拟热力学多能级/多能态激发效应和重粒子多级电离反应机制，数值求解多能级温度热力学非平衡-多电离成分化学非平衡的流动控制方程组，直至满足计算收敛精度要求，计算流场稳定。

在s4模拟基础上进一步考虑热力学多能态/多能级激发效应或重粒子的多级电离反应过程，此时需要引入多个热力学温度分别描述组分的平转动能、电子能以及其它能态对应的振动能，这里基于多振动温度模型热化学非平衡控制方程组，采用11组分（n2，o2，no，n、o、o⁺、o2⁺、n⁺、n2⁺、no⁺和e^-）气体模型进行计算，依据计算收敛精度要求判断计算过程是否需要停止。当满足计算收敛精度要求时，完成计算过程获取流动参数。

多振动温度热化学非平衡控制方程组的守恒积分形式见式（1），其中守恒量，源项，表示第个组分的密度，为第个组分方程的化学反应生成源项，分别为振动模态对应的振动能和振动能量方程源项，分别为电子能和电子能量方程源项。在多振动温度模型中，比内能，其中为组分个数，表示组分的振动模态个数，分别表示第个组分的质量分数、平转动能和电子能，表示组分第个振动模态对应的振动能。因此平转动能对应平转动温度，电子能对应电子温度，振动能由多个振动温度对应描述和计算。数值模拟过程采用lusgs时间推进，源项采用点隐式计算方法。

计算收敛精度要求包括但不限于残差条件和总迭代步数要求等，根据自身程序设计决定。这里给出一种采用残差条件的计算收敛判据，即

,

其中：表示时间迭代步序号，表示全流场控制体单元序号，为计算残差值，阈值。

需要特别说明的是，上述分段渐近逼近模拟过程根据计算的实际需要可以整合或删减某些中间步骤，例如前述模拟过程可以删除步骤s4，当s3流动模拟结束后直接切换至步骤s5。再例如对某个工况的模拟只需要用到两温度7组分气体时，可以删除步骤s5，将步骤s4调整为采用两温度7组分气体模拟。

实施例一：orex返回器高超声速再入流动模拟。计算模型和飞行试验条件参考文献“yukimitsuy.,minakoy.,cfdandfemcouplinganalysisoforexaerothermodynamicflightdata,aiaa95-2087,29thaiaathermophysicsconference,june19-22,1995/sandiego,california,1995:1-14”。计算模拟状态为飞行高度59.6km，来流速度5561.6m/s，来流温度248.12k，壁面温度1519k。数值模拟采用park化学模型，对流通量选取ausmdv格式和minmod限制器，cfl数取10，壁面条件为等温非催化壁。

如图2所示，基于实例一给出了本发明方法与现有方法的计算残差对比，横坐标表示计算时间（twall）单位秒（s），n2为氮气，纵坐标表示平均计算残差（无单位）。现有方法为不采用渐进逼近方式，直接求解完整的热化学非平衡流动控制方程的方法。结果显示，采用一般方法，计算结果发散，而采用本发明，却能较好地收敛。可见，本发明能提升计算稳定性。

实施例二：heg激波风洞圆柱模型实验流动模拟。计算模型和实验条件参考文献“hannemannk.,highenthalpyflowsinthehegshocktunnel:experimentandnumericalrebuilding,aiaa2003-0978,41staerospacesciencesmeetingandexhibit6-9january2003,reno,nevada,2003:1-18”。数值模拟采用dunn-kang化学模型，对流通量选取ausmdv格式和minmod限制器，cfl数取10，壁面条件为等温非催化壁（壁面温度300k）。

如图3所示，基于实例二给出了本发明方法与现有方法的计算残差对比，横坐标表示迭代步数（无单位），纵坐标表示平均计算残差（无单位），图中tv(k)为振动温度。结果显示，本发明方法计算收敛时的站位较一般方法提前约150s，计算耗时减少约32%。可以看出，本发明具有加快收敛速度和提升计算效率的作用。

本发明并不局限于前述的具体实施方式。本发明扩展到任何在本说明书中披露的新特征或任何新的组合，以及披露的任一新的方法或过程的步骤或任何新的组合。