CN112417743B - 一种气体能量反演热力学温度的混合迭代方法 - Google Patents

一种气体能量反演热力学温度的混合迭代方法 Download PDF

Info

Publication number
CN112417743B
CN112417743B CN202110094711.XA CN202110094711A CN112417743B CN 112417743 B CN112417743 B CN 112417743B CN 202110094711 A CN202110094711 A CN 202110094711A CN 112417743 B CN112417743 B CN 112417743B
Authority
CN
China
Prior art keywords
temperature
energy
iteration
calculation
gas
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.)
Active
Application number
CN202110094711.XA
Other languages
English (en)
Other versions
CN112417743A (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 CN202110094711.XA priority Critical patent/CN112417743B/zh
Publication of CN112417743A publication Critical patent/CN112417743A/zh
Application granted granted Critical
Publication of CN112417743B publication Critical patent/CN112417743B/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
    • G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
    • 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/08—Thermal analysis or thermal optimisation

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Computer Hardware Design (AREA)
  • Evolutionary Computation (AREA)
  • Geometry (AREA)
  • General Engineering & Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Management, Administration, Business Operations System, And Electronic Commerce (AREA)
  • Measuring Temperature Or Quantity Of Heat (AREA)

Abstract

本发明公开了一种气体能量反演热力学温度的混合迭代方法。该方法主要用于高超声速非平衡流动数值模拟过程中，气体模态能量（平动能、转动能、振动能、电子能及其组合的等效能）向对应气体热力学温度的子迭代反演计算过程。该方法在牛顿迭代法的基础上，结合气体能量反演特征，引入局部二分法修正计算判据，将牛顿迭代法和二分法结合，形成混合迭代计算方法。该方法保留了牛顿迭代法“迭代初值趋近于真值时收敛较快”的优点，计算效率较高，又吸收二分法“对于单调函数具有高稳定性”的优点，避免了极端条件下牛顿迭代法的发散问题。

Description

一种气体能量反演热力学温度的混合迭代方法
技术领域
本发明涉及数值模拟计算领域，具体涉及采用各种热力学温度模型时气体各模态能量反演相应热力学温度的迭代计算方法。
背景技术
在高超声速非平衡流动数值模拟过程中，由于涉及多种气体组分混合且每种组分的分子（或原子）的能级、能量特征存在较大差异，因此由混合气体各模态能量（平动能、转动能、振动能、电子能中的单个、多个或多个能组合的等效能）往往无法直接解析计算得到对应模态的等效热力学温度，常依赖迭代计算方法进行反演。
当前，由混合气体各模态能量反演气体温度最常用的迭代方法是牛顿迭代法和二分法，这两种方法，各有优缺点。
牛顿迭代法的优点在于能较好的利用能量函数的微分性质和温度初值。由于较大部分区域能量函数的微分性质变化相对平缓（例如冻结组分的平转动能量），同时非平衡控制方程组每一步时间推进过程中流场绝大部分区域温度初值与真值较为接近，因此牛顿迭代法的计算效率较高。它的缺点在于适用范围受到一定限制，当局部区域能量函数变化较为复杂尖锐（例如重粒子的电子束缚能）或者温度初值远离真值时，可能出现计算结果振荡，不收敛，甚至发散。为了保证牛顿迭代法能较好的收敛，常引入松弛因子等人工参数进行松弛迭代调节，这一方面增加了额外的计算量，另一方面并不能完全保证结果收敛，需要根据实际情况进行人工调节。
二分法的优点在于对于单调函数（能量函数一般都是热力学温度的单调增函数），稳定性好，从而保证能量反演温度的稳定收敛；其缺点在于无法继承“流场绝大部分区域的温度初值接近于真值”的特点，每一次反演计算都需要从计算温度范围的上下限，经过很多次迭代计算，逐步趋近于真值，因而计算效率较低。同时由于采用二分法时，真值必须处于计算温度区间内，对于计算区间的给定必须足够的“宽广”，需要预估全流场中各模态温度的所能达到的上下限，并给出足够的冗余，这将进一步降低计算效率。
发明内容
本发明的目的在于提供一种气体能量反演热力学温度的混合迭代方法。该方法基于牛顿迭代法构建，继承了牛顿迭代法的优点，能较好地利用大部分计算区域“能量函数微分性质变化较为平缓”和“温度初值接近真值”的特点，因而计算效率较高；在此基础上，利用局部二分法修正计算判据，准确捕捉局部区域牛顿迭代出现振荡的过程，引入局部二分法迭代，避免了牛顿迭代法“极端情况下发散”问题，因而稳定性较好；采用局部二分法时，计算温度区间能利用牛顿迭代过程给出，保留“温度初值接近真值”的优点，减少不必要的冗余计算。
为了实现上述目的，本发明采用如下技术方案：
一种气体能量反演热力学温度的混合迭代方法，其特征在于在高超声速非平衡流动数值模拟过程中，基于牛顿迭代法，结合热力学温度模型气体能量特征，引入局部二分法修正计算判据，将牛顿迭代法和二分法结合，形成气体模态能量反演气体热力学温度的混合迭代计算方法，具体过程为：
步骤一：在流动模拟过程中获取需要反演的模态能量分布，及其相关量的初值；
步骤二：判断需要进行迭代计算的流场区域；
步骤三：对于需要进行迭代的流场区域，由牛顿迭代法计算对应模态温度；
步骤四：判断牛顿迭代是否收敛；
步骤五：对于未收敛的区域，基于修正判据，判断牛顿迭代法的不适用的流场区域；
步骤六：对于牛顿迭代法不适用的区域，从牛顿迭代中，获取二分法需要的热力学温度区间；
步骤七 ：由二分法计算对应模态温度，并判断是否收敛；
步骤八：对于未收敛的区域，更新计算温度区间，重复步骤七、八，直至二分法收敛，得到需要反演的模态能量对应的热力学温度分布。
本发明主要用于高超声速非平衡流动数值模拟过程中，气体各模态能量向气体温度的子迭代反演计算过程。这里的气体模态能量可以是热力学温度模型中混合气体平动能、转动能、振动能、电子能中的任意一个、多个或者多个能组合的等效能，气体温度则为对应的平动温度、转动温度、振动温度、电子温度中的一个、多个或者多个组合的等效温度。
综上所述，由于采用了上述技术方案，本发明的有益效果是：
本发明结合热力学温度模型气体各模态能量与对应温度的关系特征，引入局部二分法修正计算判据，将牛顿迭代法和二分法结合，形成混合迭代计算方法，保留了牛顿迭代法与二分法各自的优点，避免了这两者的缺点，计算效率较高、稳定性较好。
本发明适用范围较广：可用于化学非平衡流动、热力学非平衡流动或热化学非平衡流动数值模拟过程；涉及的热力学温度模型可以是热力学一温度模型、热力学两温度模型、热力学三温度模型、热力学多振动温度模型以及更精细的热力学温度模型中的任意一种；气体模态能量可以是气体平动能、转动能、振动能、电子能中的任意一个、多个或者多个能组合的等效能，气体温度则为对应的平动温度、转动温度、振动温度、电子温度中的一个、多个或者多个组合的等效温度。
本发明对气体介质没有限定，可普遍适用于各类气体介质，如地球大气、火星大气、高温燃气等混合或单一组分气体。
附图说明
本发明将通过例子并参照附图的方式说明，其中：
图1为本方案的计算流程示意图；
图2为采用本方案进行不同子迭代方法电子激发温度的迭代结果比较。
具体实施方式
本说明书中公开的所有特征，或公开的所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以以任何方式组合。
本说明书（包括任何附加权利要求、摘要和附图）中公开的任一特征，除非特别叙述，均可被其他等效或具有类似目的的替代特征加以替换。即，除非特别叙述，每个特征只是一系列等效或类似特征中的一个例子而已。
如图1 所示，是本实施例的计算流程，以热力学三温度模型中混合气体振动能量反演计算得到振动温度为实例，具体过程为：
S1: 在高超声速非平衡流动数值模拟过程中，获取网格微元上混合气体振动能量初值
Figure 655917DEST_PATH_IMAGE001
、振动温度初值
Figure 152758DEST_PATH_IMAGE002
以及需要反演的振动能量
Figure 927816DEST_PATH_IMAGE003
等参数的分布。
高超声速非平衡流动数值模拟一般通过流动控制方程组数值迭代计算完成，这里将这一迭代过程称为“外迭代”。在外迭代的每一步迭代开始时，网格微元上混合气体振动能量初值
Figure 179806DEST_PATH_IMAGE001
、振动温度初值
Figure 548470DEST_PATH_IMAGE002
均为已知条件，经流动控制方程组离散计算之后，就能得到“新的”振动能量
Figure 973635DEST_PATH_IMAGE003
，可用于流动控制方程组下一步的求解。在这一过程中，还需要得到与“振动能量
Figure 806462DEST_PATH_IMAGE003
”对应的 “新的”振动温度
Figure 901457DEST_PATH_IMAGE004
，但由于混合气体振动能量与振动温度之间的计算关系复杂，由振动能量无法直接解析计算得到振动温度，因此常常需要进行反演迭代计算，这里将这一反演迭代计算过程称为“子迭代”。本实施例主要针对反演子迭代计算过程，主要用于：由
Figure 882051DEST_PATH_IMAGE001
、
Figure 720694DEST_PATH_IMAGE002
、
Figure 470344DEST_PATH_IMAGE003
等已知条件，通过子迭代反演计算得到
Figure 64137DEST_PATH_IMAGE003
对应的气体振动温度
Figure 141814DEST_PATH_IMAGE004
。
S2:根据网格微元上模态能量值的相对差异与迭代计算精度要求
Figure 174361DEST_PATH_IMAGE005
的大小比值，判定需要进行迭代反演计算的流场区域。
当在网格微元上满足
Figure 653884DEST_PATH_IMAGE006
时，则该网格微元需要进行迭代反演计算，设定
Figure 424437DEST_PATH_IMAGE001
和
Figure 114044DEST_PATH_IMAGE002
为子迭代反演计算的初值，即在迭代次数
Figure 294490DEST_PATH_IMAGE007
时，
Figure 18732DEST_PATH_IMAGE008
；当在网格微元上满足
Figure 626431DEST_PATH_IMAGE009
时，则该网格微元不需要进行迭代反演计算，该网格微元上振动能量
Figure 6597DEST_PATH_IMAGE003
对应的振动温度
Figure 849788DEST_PATH_IMAGE010
。
Figure 303903DEST_PATH_IMAGE011
为振动能量初值
Figure 207137DEST_PATH_IMAGE001
与目标值
Figure 12282DEST_PATH_IMAGE003
之间的相对差异，
Figure 659164DEST_PATH_IMAGE012
和
Figure 967786DEST_PATH_IMAGE013
分别为第
Figure 245183DEST_PATH_IMAGE014
步迭代的振动能量和振动温度。迭代计算精度要求
Figure 131100DEST_PATH_IMAGE005
存在两种形式：对于高超声速非定常流动模拟， 迭代计算精度要求
Figure 457039DEST_PATH_IMAGE015
；对于高超声速定常或准定常流动模拟，
Figure 10380DEST_PATH_IMAGE005
取值动态变化：
Figure 458679DEST_PATH_IMAGE016
Figure 238416DEST_PATH_IMAGE017
n为高超声速非平衡流场控制方程的推进步数，c为计算精度调节因子，
Figure 227101DEST_PATH_IMAGE018
为控制方程总的计算步数，
Figure 572631DEST_PATH_IMAGE019
为高超声速非平衡流场模拟的整体精度要求，由实际的计算任务给定，一般有
Figure 129515DEST_PATH_IMAGE020
。
对于满足
Figure 990023DEST_PATH_IMAGE009
的区域，振动能量初值
Figure 657765DEST_PATH_IMAGE001
与目标值
Figure 920119DEST_PATH_IMAGE003
之间的相对差异很小，因此可认为
Figure 975800DEST_PATH_IMAGE010
，无需子迭代计算；反之，则需要继续进行子迭代计算。之所以在S2进行迭代流场区域判断，是因为在控制方程外迭代过程中，流场部分区域
Figure 933391DEST_PATH_IMAGE003
变化较小（例如在头部主激波之前区域气体的振动能量可能基本保持不变），随着外迭代次数的增多，这样的区域可能进一步增大。这些区域可能直接就满足了
Figure 812615DEST_PATH_IMAGE009
的要求，因此无需进入子迭代过程，从而减少计算量。对于“定常流动状态”、“准定常流动状态”流动控制方程外迭代过程，随外迭代次数n增加，
Figure 804842DEST_PATH_IMAGE021
逐渐趋近于0，因此精度判据
Figure 765845DEST_PATH_IMAGE022
取值考虑了两方面因素：一是
Figure 866525DEST_PATH_IMAGE022
随着外循环迭代次数n增加而减小，从
Figure 469544DEST_PATH_IMAGE023
降低至
Figure 316278DEST_PATH_IMAGE024
，这主要是由于在外循环迭代初期，流场处于初始构建阶段，流场主体结构尚未达到基本稳定，外循环得到的变量波动较大，与真实值的差异较大，其误差占主导，此时其内部的子迭代过程计算，无需过高精度；二是
Figure 713761DEST_PATH_IMAGE022
不小于外循环迭代残差
Figure 770579DEST_PATH_IMAGE011
的10%，这样既能捕捉外循环过程中振动能量的主体变化特征，又在保证了整体计算精度的同时避免了非必要的迭代。对于“非定常状态”基于物理时间的外迭代过程，
Figure 583814DEST_PATH_IMAGE021
反应了真实物理时间变化带来的差异，它可能不随外迭代次数增加而减小，因此不能将
Figure 675267DEST_PATH_IMAGE022
与外迭代步数n相关联，这里取
Figure 650176DEST_PATH_IMAGE015
，即内迭代精度要求高于整体计算精度要求一个量级，从而保证整体精度。
S3:对于需要进行迭代计算的流场网格微元，采用牛顿迭代法，结合混合气体振动能量和振动温度关系，利用切线逼近计算气体振动温度，具体计算式为：
Figure 725448DEST_PATH_IMAGE025
其中：
Figure 76795DEST_PATH_IMAGE026
为第
Figure 22754DEST_PATH_IMAGE027
步迭代的振动温度，
Figure 168565DEST_PATH_IMAGE014
为迭代次数，
Figure 403237DEST_PATH_IMAGE028
为第m步迭代的等效振动能比热，由
Figure 948488DEST_PATH_IMAGE013
通过等效比热表达式计算得到，
Figure 624320DEST_PATH_IMAGE029
。
由混合气体分子振动能理论，可以得到振动能等效比热的解析表达式：
Figure 331245DEST_PATH_IMAGE030
式中，第一个求和符号表示对混合气体中的分子组分求和，第二个求和符合表示对气体第
Figure 725317DEST_PATH_IMAGE031
组分的所有的振动模态求和，
Figure 746363DEST_PATH_IMAGE032
为气体第
Figure 666914DEST_PATH_IMAGE034
组分的振动模态数，
Figure 420107DEST_PATH_IMAGE035
和
Figure 160529DEST_PATH_IMAGE036
分别为第
Figure 922949DEST_PATH_IMAGE031
组分气体的质量分数和分子量，
Figure 698007DEST_PATH_IMAGE037
为普适气体常数，
Figure 622101DEST_PATH_IMAGE038
和
Figure 324521DEST_PATH_IMAGE039
分别为第
Figure 484107DEST_PATH_IMAGE031
组分气体的第
Figure 989037DEST_PATH_IMAGE014
个振动能模态的振动特征温度和简并度。
S4：根据当前网格微元上振动能量值的相对差异与迭代计算精度要求
Figure 474245DEST_PATH_IMAGE022
的大小比值，判断第
Figure 64627DEST_PATH_IMAGE027
步牛顿迭代是否收敛；
如果满足
Figure 27903DEST_PATH_IMAGE040
时，则第
Figure 980816DEST_PATH_IMAGE027
步牛顿迭代收敛，该网格微元上振动能量
Figure 246712DEST_PATH_IMAGE003
对应的振动温度
Figure 449023DEST_PATH_IMAGE041
；如果不满足
Figure 91357DEST_PATH_IMAGE042
，则未收敛，这里
Figure 695514DEST_PATH_IMAGE043
为第
Figure 132311DEST_PATH_IMAGE027
步迭代的振动能量，由
Figure 821919DEST_PATH_IMAGE044
通过振动能量与振动温度关系计算得到。
基于分子振动能理论，得到混合气体振动能解析表达式：
Figure 2364DEST_PATH_IMAGE045
。
S5：对于未收敛的网格微元，结合振动能量迭代变化特征，采用局部二分法修正计算判据捕捉牛顿迭代法的不适用的流场区域，网格微元上满足局部二分法修正计算判据：
当满足
Figure 929869DEST_PATH_IMAGE046
时，该网格微元不适用牛顿迭代法（发散或者收敛速度相对较慢），
当满足
Figure 662202DEST_PATH_IMAGE047
时，该网格微元适用牛顿迭代法，进行下一步牛顿迭代，即
Figure 448892DEST_PATH_IMAGE048
，返回步骤S3。
该过程是牛顿迭代法和二分法结合的关键，对于满足
Figure 823242DEST_PATH_IMAGE049
的流场区域，采用二分法；反之，则继续采用牛顿迭代法。
首先分析牛顿迭代法适用区间，即满足
Figure 605253DEST_PATH_IMAGE050
的区间，亦即
Figure 383853DEST_PATH_IMAGE051
的区间。由于S5执行的前提满足
Figure 782473DEST_PATH_IMAGE052
和
Figure 304722DEST_PATH_IMAGE053
，因此该区间可分为
Figure 737977DEST_PATH_IMAGE054
和
Figure 953058DEST_PATH_IMAGE055
两个区间进行分析。
当
Figure 104553DEST_PATH_IMAGE054
时，由
Figure 492809DEST_PATH_IMAGE056
可知
Figure 921517DEST_PATH_IMAGE057
，即有
Figure 629535DEST_PATH_IMAGE058
，又由
Figure 268327DEST_PATH_IMAGE059
恒大于零可知
Figure 460274DEST_PATH_IMAGE060
与
Figure 743488DEST_PATH_IMAGE061
同正负号，即有
Figure 425005DEST_PATH_IMAGE062
,可见
Figure 160880DEST_PATH_IMAGE044
必然处于
Figure 953255DEST_PATH_IMAGE004
和
Figure 90975DEST_PATH_IMAGE013
之间，也就是说随子迭代进行，计算得到的振动温度将不断渐进逼近真值，不会出现数值振荡，因此牛顿迭代法合理可行。
当
Figure 881077DEST_PATH_IMAGE063
时，可等价为同时满足
Figure 228882DEST_PATH_IMAGE064
和
Figure 700314DEST_PATH_IMAGE065
两个条件。对于
Figure 817175DEST_PATH_IMAGE066
，由
Figure 715861DEST_PATH_IMAGE067
恒大于零可知
Figure 550962DEST_PATH_IMAGE060
与
Figure 560506DEST_PATH_IMAGE061
同正负号、
Figure 735135DEST_PATH_IMAGE068
与
Figure 929356DEST_PATH_IMAGE069
同正负号，即有
Figure 127119DEST_PATH_IMAGE070
，可见
Figure 64988DEST_PATH_IMAGE004
处于
Figure 766228DEST_PATH_IMAGE044
与
Figure 131350DEST_PATH_IMAGE013
之间，即牛顿迭代计算得到振动温度在真值两侧来回振荡。但由于同时满足了
Figure 878727DEST_PATH_IMAGE071
的条件，即每次子迭代，能量差异的幅值衰减50%以上，牛顿迭代法仍能迅速收敛，其收敛速度不小于“2的指数级” (即每次迭代能量差异的幅值减小50%以上)。
反之，则说明牛顿迭代法出现振荡且收敛速度较慢，此时采用二分法具有一定优势。二分法迭代具有较为稳定的收敛速度（“2的指数级”收敛，即每次迭代计算的振动温度区间范围减半），同时由于真值
Figure 230073DEST_PATH_IMAGE004
处于
Figure 441612DEST_PATH_IMAGE044
与
Figure 853002DEST_PATH_IMAGE013
之间，恰好满足了二分法对计算区间的要求，因此可将
Figure 884412DEST_PATH_IMAGE044
与
Figure 39449DEST_PATH_IMAGE013
引入二分法计算温度区间的设定过程，保留“振动温度初值接近真值”的优点，避免二分法“人工设定振动温度范围冗余较大、计算效率低”的缺点。
S6：对于牛顿迭代法不适用的网格微元，则采用二分法进行修正计算，具体为：
从牛顿迭代过程中，获取二分法迭代需要的振动温度区间
Figure 845775DEST_PATH_IMAGE072
，
Figure 755962DEST_PATH_IMAGE073
，函数min为取
Figure 150034DEST_PATH_IMAGE074
中的最小值，
Figure 967817DEST_PATH_IMAGE075
，函数max为取
Figure 763735DEST_PATH_IMAGE076
中的最大值，
Figure 641561DEST_PATH_IMAGE077
为人为设定的流场中振动温度的最小值，
Figure 522929DEST_PATH_IMAGE078
为人为设定的流场中振动温度的最大值。
由
Figure 347666DEST_PATH_IMAGE026
和
Figure 122724DEST_PATH_IMAGE013
计算二分法迭代所需要的振动温度区间，保留了“温度初值接近真值”的优点，避免了二分法“人工设定温度范围冗余较大、计算效率低”的缺点。
S7:以
Figure 46818DEST_PATH_IMAGE079
为计算振动温度区间，计算振动温度中值
Figure 274537DEST_PATH_IMAGE080
，由
Figure 575068DEST_PATH_IMAGE081
通过振动能量关系计算振动能量中值
Figure 470212DEST_PATH_IMAGE082
，判断
Figure 565207DEST_PATH_IMAGE082
是否满足二分法收敛的精度要求；
当满足
Figure 280222DEST_PATH_IMAGE083
时，计算结果收敛，该网格微元上振动能量
Figure 118865DEST_PATH_IMAGE003
对应的振动温度
Figure 71777DEST_PATH_IMAGE084
；反之，则未收敛；这里
Figure 462307DEST_PATH_IMAGE085
为二分法得到的振动能量中值与目标值之间的相对差异。
这里由振动温度
Figure 805564DEST_PATH_IMAGE086
计算
Figure 306952DEST_PATH_IMAGE087
的解析表达式为：
Figure 114371DEST_PATH_IMAGE045
。
S8：对于二分法未收敛的网格微元，以振动能量
Figure 551169DEST_PATH_IMAGE003
为逼近目标，修改二分法迭代计算的振动温度区间，然后在新的温度区间条件下执行S7；
修改振动温度区间的方法为：
当
Figure 240776DEST_PATH_IMAGE088
时，则
Figure 686801DEST_PATH_IMAGE089
，继续采用二分法迭代，返回S7再次计算；
当
Figure 145464DEST_PATH_IMAGE090
时，则
Figure 753163DEST_PATH_IMAGE091
，返回S7再次计算。
通过S7、S8的反复迭代，直至二分法计算收敛，得到网格微元上的振动能量
Figure 930067DEST_PATH_IMAGE092
对应的振动温度
Figure 259677DEST_PATH_IMAGE004
的分布。
实施例一：单网格点氮气的电子束缚能
Figure 448213DEST_PATH_IMAGE093
反演计算电子温度
Figure 617026DEST_PATH_IMAGE094
。
计算场景：针对气体介质为氮气的高超声速化学反应冻结定常流动，其流动控制方程的时间推进过程中，电子束缚能量和电子温度的初值分别为
Figure 156592DEST_PATH_IMAGE095
和
Figure 6737DEST_PATH_IMAGE096
，经流动控制方程组离散计算之后，得到的新的电子束缚能量
Figure 705571DEST_PATH_IMAGE093
，这里主要考虑氮气分子束缚电子第0和第1激发能级的电子束缚能。
计算的目的：由
Figure 655073DEST_PATH_IMAGE097
等，通过迭代过程反演计算得到
Figure 9831DEST_PATH_IMAGE094
，这里假设需要反演的
Figure 194824DEST_PATH_IMAGE098
，即目标真值
Figure 623532DEST_PATH_IMAGE099
。
经测算，采用牛顿迭代法时，当
Figure 134147DEST_PATH_IMAGE100
时，迭代计算收敛；而当
Figure 710622DEST_PATH_IMAGE101
，迭代计算结果发散。可以看出，在较大区间内，牛顿迭代法均能收敛，由于随定常流动时间推进，
Figure 574673DEST_PATH_IMAGE102
逐渐趋近于真值，因此在绝大部分情况下，牛顿迭代法都可以较好地适用。但某些特殊情况下，当流场中由于复杂流动现象（如激波强间断等）或者某种非物理的因素（如流场初始构建过程中数值误差较大等）出现较大波动，造成
Figure 982521DEST_PATH_IMAGE096
远离真值，此时牛顿迭代法就可能不收敛，从而导致整体计算不收敛，甚至发散。
图2为采用不同子迭代方法电子激发温度的迭代结果，这里为了较好的显示电子温度的变化，避免彻底的发散，子迭代过程中将电子激发温度限定在
Figure 539404DEST_PATH_IMAGE103
范围。图中m为子迭代步数，
Case1为牛顿迭代法
Figure 399912DEST_PATH_IMAGE104
时结果；Case2为牛顿迭代法
Figure 67654DEST_PATH_IMAGE105
时结果；
Case3为本发明
Figure 533271DEST_PATH_IMAGE105
的结果；Case4为本发明
Figure 385689DEST_PATH_IMAGE106
的结果；Case5为本发明
Figure 608860DEST_PATH_IMAGE107
的结果。
由图可以看出，采用牛顿迭代法，当子迭代初值
Figure 939347DEST_PATH_IMAGE104
时计算结果收敛，当
Figure 259470DEST_PATH_IMAGE105
时结果不收敛；
而采用本发明
Figure 158156DEST_PATH_IMAGE108
时，计算均能较快地收敛。
这说明本发明稳定性好，能较好地解决极端情况下牛顿迭代法的不收敛问题。
实施例二：全流场空气热离解/电离混合气体，振动-电子能
Figure 258836DEST_PATH_IMAGE109
反演计算振动-电子温度
Figure 2801DEST_PATH_IMAGE110
。
计算场景：采用RAM-C钝锥外形，计算网格约25万，计算飞行高度61km、速度7650m/s；空气化学反应模型采用7组分Park模型，热力学模型采用Park的热力学两温度模型；外迭代求解的流动控制方程组为热化学非平衡N-S方程组；考虑不同的计算时刻，外迭代推进步数n分别为1000、2000和10000步；考虑不同的反演计算方法，
Figure 239747DEST_PATH_IMAGE109
反演
Figure 43755DEST_PATH_IMAGE111
的子迭代分别采用牛顿迭代法、二分法和本发明方法。为了保证牛顿迭代法收敛，引入松弛迭代，松弛因子取为0.5；采用二分法时，振动-电子温度的计算区间设为
Figure 372011DEST_PATH_IMAGE112
。
经重复测试，这三者平均的用时（单位秒）分别为：
Figure 247564DEST_PATH_IMAGE113
可以看出，本发明的计算效率远高于二分法，略优于牛顿迭代法。
本发明并不局限于前述的具体实施方式。本发明扩展到任何在本说明书中披露的新特征或任何新的组合，以及披露的任一新的方法或过程的步骤或任何新的组合。

Claims (7)

1.一种气体能量反演热力学温度的混合迭代方法，其特征在于在高超声速非平衡流动数值模拟过程中，基于牛顿迭代法，结合热力学温度模型气体能量特征，引入局部二分法修正计算判据，将牛顿迭代法和二分法结合，形成气体模态能量反演气体温度的混合迭代计算方法，具体过程为：
S1: 在高超声速非平衡流动数值模拟过程中，获取网格微元上混合气体模态能量初值
Figure DEST_PATH_IMAGE002
、对应模态温度初值
Figure DEST_PATH_IMAGE004
以及需要反演的模态能量
Figure DEST_PATH_IMAGE006
等参数的分布；
S2:根据网格微元上模态能量值的相对差异
Figure DEST_PATH_IMAGE008
与迭代计算精度要求
Figure DEST_PATH_IMAGE010
的大小比值，判定需要进行迭代反演计算的流场区域，
当在网格微元上满足
Figure DEST_PATH_IMAGE012
时，则该网格微元需要进行迭代反演计算，在迭代次数为
Figure DEST_PATH_IMAGE014
时，
Figure DEST_PATH_IMAGE016
；
当在网格微元上满足
Figure DEST_PATH_IMAGE018
时，则该网格微元需要进行迭代反演计算，该网格微元上模态能量
Figure 796230DEST_PATH_IMAGE006
对应的温度
Figure DEST_PATH_IMAGE020
；
S3:对于需要进行迭代计算的流场网格微元，采用牛顿迭代法，结合混合气体各模态能量和温度关系，利用切线逼近计算气体温度，具体计算式为：
Figure DEST_PATH_IMAGE022
其中：
Figure DEST_PATH_IMAGE024
为第
Figure DEST_PATH_IMAGE026
步迭代的模态温度，
Figure DEST_PATH_IMAGE028
为迭代次数，
Figure DEST_PATH_IMAGE030
为第
Figure 274223DEST_PATH_IMAGE028
步迭代的模态能量，
Figure DEST_PATH_IMAGE032
为第
Figure 460485DEST_PATH_IMAGE028
步迭代的等效比热，
Figure 777066DEST_PATH_IMAGE032
由
Figure DEST_PATH_IMAGE034
通过等效比热表达式计算得到，
Figure DEST_PATH_IMAGE036
；
S4：根据当前网格微元上模态能量值的相对差异与迭代计算精度要求
Figure 981651DEST_PATH_IMAGE010
的大小比值，判断第
Figure 272955DEST_PATH_IMAGE026
步牛顿迭代是否收敛；
S5：对于未收敛的网格微元，结合模态能量迭代变化特征，捕捉牛顿迭代法的不适用的流场区域；
S6：对于牛顿迭代法不适用的网格微元，则采用二分法进行修正计算，具体为：
从牛顿迭代过程中，获取二分法迭代需要的温度区间
Figure DEST_PATH_IMAGE038
，
Figure DEST_PATH_IMAGE040
，函数min为取
Figure DEST_PATH_IMAGE042
中的最小值，
Figure DEST_PATH_IMAGE044
，函数max为取
Figure DEST_PATH_IMAGE046
中的最大值，
Figure DEST_PATH_IMAGE048
为人为设定的流场中模态温度的最小值，
Figure DEST_PATH_IMAGE050
为人为设定的流场中模态温度的最大值，
S7:以
Figure 222544DEST_PATH_IMAGE038
为计算温度区间，计算模态温度中值
Figure DEST_PATH_IMAGE052
，由
Figure DEST_PATH_IMAGE054
通过气体能量关系直接计算模态能量中值
Figure DEST_PATH_IMAGE056
，判断能量中值
Figure 749341DEST_PATH_IMAGE056
是否收敛；
S8：对于二分法未收敛的网格微元，以模态能量
Figure 152640DEST_PATH_IMAGE006
为逼近目标，修改二分法的计算的温度区间，然后在新的温度区间条件下执行S7，通过S7、S8的反复迭代，直至满足二分法收敛要求，得到网格微元上的模态能量
Figure 67375DEST_PATH_IMAGE006
对应的气体温度
Figure DEST_PATH_IMAGE058
的分布。
2.根据权利要求1所述的一种气体能量反演热力学温度的混合迭代方法，其特征在于在S4中：
当满足
Figure DEST_PATH_IMAGE060
时，第
Figure 839283DEST_PATH_IMAGE026
步牛顿迭代收敛，该网格微元上模态能量
Figure 45136DEST_PATH_IMAGE006
对应的温度
Figure DEST_PATH_IMAGE062
，
Figure DEST_PATH_IMAGE064
为第
Figure 489893DEST_PATH_IMAGE026
步迭代的模态能量。
3.根据权利要求1所述的一种气体能量反演热力学温度的混合迭代方法，其特征在于在S5中采用局部二分法修正计算判据捕捉牛顿迭代法的不适用的流场区域，具体为：
当满足
Figure DEST_PATH_IMAGE066
时，该网格微元不适用牛顿迭代法，
当满足
Figure DEST_PATH_IMAGE068
时，该网格微元适用牛顿迭代法，进行下一步牛顿迭代，即
Figure DEST_PATH_IMAGE070
，返回步骤S3。
4.根据权利要求1所述的一种气体能量反演热力学温度的混合迭代方法，其特征在于在S7中判断二分法计算是否收敛，具体为
当满足
Figure DEST_PATH_IMAGE072
时，计算结果收敛，该网格微元上模态能量
Figure 716475DEST_PATH_IMAGE006
对应的温度
Figure DEST_PATH_IMAGE074
。
5.根据权利要求4所述的一种气体能量反演热力学温度的混合迭代方法，其特征在于S6中，对于二分法的计算的温度区间：
当
Figure DEST_PATH_IMAGE076
时，则
Figure DEST_PATH_IMAGE078
，继续采用二分法迭代，返回S7再次计算；
当
Figure DEST_PATH_IMAGE080
时，则
Figure DEST_PATH_IMAGE082
，返回S7再次计算。
6.根据权利要求1所述的一种气体能量反演热力学温度的混合迭代方法，其特征在于所述迭代计算精度
Figure 322267DEST_PATH_IMAGE010
具有两种形式，当模拟高超声速非定常流动时
Figure DEST_PATH_IMAGE084
，当模拟高超声速定常或准定常流动时，
Figure 800653DEST_PATH_IMAGE010
的取值动态变化为：
Figure DEST_PATH_IMAGE086
Figure DEST_PATH_IMAGE088
n为高超声速非平衡流场控制方程的推进步数，c为计算精度调节因子，
Figure DEST_PATH_IMAGE090
为控制方程总的计算步数，
Figure DEST_PATH_IMAGE092
为高超声速非平衡流场模拟的整体精度要求。
7.根据权利要求1所述的一种气体能量反演热力学温度的混合迭代方法，其特征在于：
所述气体模态能量由气体平动能、转动能、振动能、电子能中的任意单个、或多个组合的等效能，
所述气体温度为与能量模态对应的平动温度、转动温度、振动温度、电子温度中的任意单个、或多个组合的等效温度，
所述热力学温度模型为热力学一温度模型、热力学两温度模型、热力学三温度模型、热力学多振动温度模型中的任意一种，
所述高超声速非平衡流动模拟为化学非平衡流动、热力学非平衡流动、热化学非平衡流动中的任意一种。
CN202110094711.XA 2021-01-25 2021-01-25 一种气体能量反演热力学温度的混合迭代方法 Active CN112417743B (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110094711.XA CN112417743B (zh) 2021-01-25 2021-01-25 一种气体能量反演热力学温度的混合迭代方法

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110094711.XA CN112417743B (zh) 2021-01-25 2021-01-25 一种气体能量反演热力学温度的混合迭代方法

Publications (2)

Publication Number Publication Date
CN112417743A CN112417743A (zh) 2021-02-26
CN112417743B true CN112417743B (zh) 2021-04-02

Family

ID=74782947

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110094711.XA Active CN112417743B (zh) 2021-01-25 2021-01-25 一种气体能量反演热力学温度的混合迭代方法

Country Status (1)

Country Link
CN (1) CN112417743B (zh)

Families Citing this family (6)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN112966401B (zh) * 2021-05-12 2021-07-16 中国空气动力研究与发展中心计算空气动力研究所 一种热化学非平衡多级气体模型自适应方法
CN115312139B (zh) * 2022-09-23 2023-01-13 中国空气动力研究与发展中心计算空气动力研究所 一种高超声速流动化学反应模型数据的存取与转换方法
CN116090110B (zh) * 2023-04-07 2023-07-07 中国空气动力研究与发展中心计算空气动力研究所 高超声速飞行器高温流场数值仿真的修正方法及相关组件
CN116227388B (zh) * 2023-04-21 2023-07-07 中国空气动力研究与发展中心计算空气动力研究所 高超流动模拟cfl数动态调整方法、系统、设备及介质
CN117059188B (zh) * 2023-10-12 2024-01-23 中国空气动力研究与发展中心计算空气动力研究所 一种化学非平衡气体热力学平衡能量体系改进方法及系统
CN117332511B (zh) * 2023-12-01 2024-02-20 中国空气动力研究与发展中心计算空气动力研究所 面向高超飞行器高温非平衡流的自适应耦合数值模拟方法

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105389439A (zh) * 2015-12-03 2016-03-09 四川大学 天然气关井井筒建模模拟实现方法
CN106125674A (zh) * 2016-08-03 2016-11-16 大连理工大学 一种高精度实时轮廓误差估计方法
CN107122571A (zh) * 2017-06-06 2017-09-01 大连理工大学 一种考虑水合物分解的沉积物多场耦合模型的建模方法
CN111339672A (zh) * 2020-03-02 2020-06-26 上海索辰信息科技有限公司 进气道前缘激波气动热仿真分析方法

Family Cites Families (9)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102422179A (zh) * 2009-04-21 2012-04-18 密歇根宇航公司 大气测量系统
CN102930152B (zh) * 2012-10-26 2016-08-03 中国科学院上海药物研究所 一种模拟配体分子与靶标受体反应并计算预测该反应的热力学与动力学参数的方法和系统
JP2020530159A (ja) * 2017-08-02 2020-10-15 ストロング フォース アイオーティ ポートフォリオ ２０１６，エルエルシー 大量のデータセットを使用する産業用のモノのインターネットのデータ収集環境における検出のための方法及びシステム
CN108389136B (zh) * 2018-01-29 2022-04-22 重庆大学 一种考虑多重不确定性的天然气概率-模糊能流分析方法
JP2021519475A (ja) * 2018-03-19 2021-08-10 ダッソー システムズ ドイチェラント ゲーエムベーハー ＣＯＳＭＯｐｌｅｘ：自己組織化システムの自己一貫性シミュレーション
CN108494021B (zh) * 2018-04-20 2021-06-01 东北大学 电-热-气综合能源系统的稳定评估与静态控制方法
CN109344436B (zh) * 2018-08-28 2023-05-26 中国石油化工股份有限公司 一种大型复杂天然气管网系统在线仿真方法
CN109858129B (zh) * 2019-01-23 2021-04-13 清华大学 一种关于三联供系统的燃气轮机动态仿真方法
CN111476394B (zh) * 2020-02-28 2022-07-22 浙江工业大学 一种适用于电热气等多能源系统的鲁棒运行优化方法

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105389439A (zh) * 2015-12-03 2016-03-09 四川大学 天然气关井井筒建模模拟实现方法
CN106125674A (zh) * 2016-08-03 2016-11-16 大连理工大学 一种高精度实时轮廓误差估计方法
CN107122571A (zh) * 2017-06-06 2017-09-01 大连理工大学 一种考虑水合物分解的沉积物多场耦合模型的建模方法
CN111339672A (zh) * 2020-03-02 2020-06-26 上海索辰信息科技有限公司 进气道前缘激波气动热仿真分析方法

Also Published As

Publication number Publication date
CN112417743A (zh) 2021-02-26

Similar Documents

Publication Publication Date Title
CN112417743B (zh) 一种气体能量反演热力学温度的混合迭代方法
Aftosmis et al. Behavior of linear reconstruction techniques on unstructured meshes
WO2019144337A1 (zh) 一种基于深度学习算法的航空发动机全包线模型自适应修正方法
CN112613250B (zh) 火星进入器表面流动转捩位置预测方法
CN112966401B (zh) 一种热化学非平衡多级气体模型自适应方法
CN109376403B (zh) 一种基于笛卡尔自适应重构技术的二维结冰数值模拟方法
CN111594322B (zh) 一种基于Q-Learning的变循环航空发动机推力控制方法
Gnoffo Updates to multi-dimensional flux reconstruction for hypersonic simulations on tetrahedral grids
Montero Villar et al. Multi-objective optimization of an counter rotating open rotor using evolutionary algorithms
CN114611416A (zh) 导弹非线性非定常气动特性ls-svm建模方法
Kermani et al. Thermodynamically based moisture prediction using Roe’s scheme
Lengyel-Kampmann et al. Generalized optimization of counter-rotating and single-rotating fans
CN117059188A (zh) 一种化学非平衡气体热力学平衡能量体系改进方法及系统
Darbandi et al. Conceptual linearization of Euler governing equations to solve high speed compressible flow using a pressure‐based method
Mist’e et al. Improvements in off design aeroengine performance prediction using analytic compressor map interpolation
Qi et al. Design and experimental calibration of the profile rotating wind tunnel with Mach number varying from 2.0 to 4.0
Skujins et al. Toward an unsteady aerodynamic ROM for multiple mach regimes
Page et al. Inverse design of 3D multi-stage transonic fans at dual operating points
Ren et al. Numerical stability analysis of shock-capturing methods for strong shocks I: second-order MUSCL schemes
Chen et al. On-board Rapid Planning for Booster Trajectory of Multistage Rocket Based on SQP and Surrogate Model
Zhang et al. Optimization of cycle parameters of variable cycle engine based on response surface model
Katzenmeier et al. Correction technique for quality improvement of doublet lattice unsteady loads by introducing CFD small disturbance aerodynamics
CN109816088A (zh) 一种压电精密定位平台参数辨识方法及其应用
Krivcov et al. Account the mutual influence of the simulation components of GTE
VAN DEN BRAEMBUSSCHE Turbomachinery component design by means of CFD

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