高温非平衡双原子气体热物性计算方法及数据库生成方法与流程

1.本发明涉及计算空气动力学技术领域，具体涉及一种高温非平衡双原子气体热物性计算及数据库生成方法。


背景技术：

2.高空高速流场的高精度模拟对新型高速飞行器设计、通讯中断预测和未来载人航天器发展有重要意义。真实气体效应和非平衡性是高空高速流动的两大特征，为提升流场模拟精度，有效描述非平衡过程的物理模型和可靠的高温气体数据缺一不可。
3.在非平衡物理模型方面，经过近40余年的发展，目前建立了描述化学非平衡过程的有限速率化学反应模型，如park模型、gupta模型等；描述热非平衡过程的多温度模型，如双温度模型、三温度模型和四温度模型等。尽管现有研究表明，考虑电子能非平衡效应可以提高等离子体流场预测精度、转动非平衡现象在平动温度大于约12000k后凸显，但受限于缺乏可靠高温非平衡气体热物性数据，特别是高温非平衡双原子气体热物性数据，目前工程实践中仅有双温模型得到了广泛的应用。
4.区别于仅激发平动和转动两种能量模态的低温气体，高温双原子气体最多能激发平动、电子、振动和转动四种能量模态，其中电子、振动和转动模态合称为内自由度模态(internal mode)，它们与平动模态相互独立、高温下存在一系列复杂特性，如振动非简谐性、转动非刚性和模态耦合等。正是因为这些特性，高温非平衡气体建模和热物性求解十分困难。


技术实现要素：

5.本发明的目的在于克服上述缺陷，提供一种适用于双温以上热非平衡模型的双原子气体热物性计算方法，并据此建立高温非平衡双原子气体热物性数据库，支撑双温以上热非平衡模型的工程应用。本发明首先根据双原子分子振动-转动耦合能量在转动模态和振动模态间以一定比例分配假设，双原子分子在电子、振动、转动模态上能级分布概率分别符合玻尔兹曼分布假设，建立基于配分函数的高温非平衡气体模型；然后根据双原子分子振动、转动能之和小于分子最小离解能准则，确定双原子分子振动能级和转动能级范围；最后在模态解耦假设t
el
≈t
vib
≈t
rot
的条件下，根据高温非平衡气体模型及双原子分子振动能级和转动能级范围，导出仅包含单一温度变量的非平衡双原子气体电子、振动和转动模态比热计算式。数据库的生成方法中，基于非平衡双原子气体比热计算式，结合分段线性插值和分段二次插值方法，构建高温非平衡双原子气体热物性数据库。
6.其中，确定双原子分子振动能级和转动能级范围的优选方式中，主要步骤包括：根据双原子分子内势能函数确定转动量子数上界j
lim.n
；在j《j
lim,n
范围内根据在极大值点势能函数一阶导数为零、二阶导数小于零的性质，求解势能曲线极大值点对应的r
max,n,j
；r
max,n,j
代入势能函数，得到分子离解所需最小能量求解所有满足约束条件双原子分
子振动、转动能之和小于的振动量子数v和转动量子数j，得到振动能级和转动能级范围。
7.本发明综合考虑高阶电子能激发、振子非简谐性、转子非刚性、转动-振动耦合、离解能约束等高温真实气体效应影响，在合理假设基础上，将双原子气体电子、振动、转动模态完全解耦，获得焓、比热近似值计算方法，能够支撑三、四温模型的工程应用探索，对计算空气动力学技术领域的研究具有重要指导意义。
8.为实现上述发明目的，本发明提供如下技术方案：
9.一种高温非平衡双原子气体热物性计算方法，包括如下步骤：
10.s1根据双原子分子振动-转动耦合能量在转动模态和振动模态间分配假设，以及双原子分子在电子、振动和转动模态上能级分布概率分别符合玻尔兹曼分布假设，建立基于配分函数的高温非平衡气体模型；
11.s2根据双原子分子振动能和转动能之和小于分子最小离解能准则，确定双原子分子振动能级和转动能级范围；
12.s3 t
el
、t
vib
和t
rot
分别为双原子分子电子温度、振动温度和转动温度，在模态解耦假设t
el
≈t
vib
≈t
rot
的条件下，根据高温非平衡气体模型及双原子分子振动能级和转动能级范围，导出仅包含单一温度变量的非平衡双原子气体电子、振动和转动模态比热计算式。
13.进一步的，所述步骤s1中建立的基于配分函数的高温非平衡气体模型如下：
[0014][0015]
其中q
int
是分子内能配分函数，σ是分子对称常数，n、v、j是电子、振动和转动量子数，和分别是电子能、振动能、转动能和振动-转动耦合能，gn是能级简并度，k
bol
是玻尔兹曼常数，f为耦合能量分配控制参数，由振动温度t
vib
、转动温度t
rot
和振动特征温度确定，取值范围为0～1：ε0是零点能；
[0016]nmax
为最大电子能级，v
max,n
和j
max,n,v
分别为根据双原子分子振动能级范围和转动能级范围确定的最大振动能级和最大转动能级。
[0017]
进一步的，所述步骤s1的具体步骤如下：
[0018]
s1.1将双原子分子振动-转动能分为三部分：只受振动模态影响部分，即振动能只受转动模态影响部分，即转动能和同时受转动模态和振动模态影响部分，即振动-转动耦合能
[0019][0020]
[0021][0022]
其中，h
plk
为普朗克常数，c为光速，ω
e,n
、ω
e,n
x
e,n
、ω
e,nye,n
、ω
e,nze,n
、α
e,n
、β
e,n
、b
e,n
、d
e,n
为经由实验测量或理论计算确定的双原子分子光谱常数；
[0023]
s1.2根据双原子分子振动-转动耦合能量在转动模态和振动模态间分配假设引入取值范围0～1的耦合能量分配控制参数f，将分为和两部分，分别分配到振动模态和转动模态中，得到新的振动能与新的转动能
[0024][0025][0026]
其中，耦合能量分配控制参数
[0027]
s1.3根据双原子分子在电子、振动和转动模态上能级分布概率分别符合玻尔兹曼分布假设，以及步骤s1.2所得新的振动能与新的转动能建立基于配分函数的高温非平衡气体模型。
[0028]
进一步的，所述步骤s2的具体步骤为：
[0029]
s2.1根据双原子分子内势能函数确定转动量子数上界j
lim.n
；
[0030]
s2.2根据双原子分子内势能函数，在0≤j《j
lim
范围内确定分子最小解离能
[0031]
s2.3求解所有满足双原子分子振动、转动能之和小于分子最小离解能准则的振动量子数v和转动量子数j组合，即振动能级和转动能级范围。
[0032]
进一步的，所述步骤s2.1中，根据双原子分子内势能函数在极值点处一阶导数为零的性质，得到如下方程，转动量子数上界j
lim.n
为使如下方程无解的最小整数：
[0033][0034]
其中，中，为转动能基态j＝0分子离解所需的最小能量，μ为双原子分子折合质量，r为组成双原子分子的原子间距离，r
e,n
为j＝0对应势能曲线最低点的r，h
plk
为普朗克常数，c为光速，ω
e,n
为经由实验测量或理论计算确定的双原子分子光谱常数，j是转动量子数；
[0035]
所述步骤2.2中，在0≤j《j
lim
范围内，求解满足双原子分子内势能函数一阶导数为零且二阶导数小于零的势能极大值r
max,n,j
，将r
max,n,j
代入双原子分子内势能函数得到
[0036]
进一步的，所述步骤s3所得非平衡双原子气体电子、振动和转动模态比热计算式如下：
[0037][0038][0039][0040]
其中，r为气体常数，q
int
是分子内能配分函数。
[0041]
一种高温非平衡双原子气体热物性数据库生成方法，根据上述非平衡双原子气体电子、振动和转动模态比热计算式，采用分段线性插值和分段二次插值方法，建立包含双原子气体比热和焓的高温非平衡双原子气体热物性数据库。
[0042]
进一步的，一种高温非平衡双原子气体热物性数据库生成方法中，在t
min
～t
max
温度范围内，采用分段线性插值方法得到各温度下双原子气体比热c
p
，采用分段二次插值法得到各温度下双原子气体焓h；
[0043]
具体步骤如下：
[0044]
(1)建立t
min
～t
max
温度范围内焓与比热的插值式为f
h,i
(t)和采用n个点将t
min
～t
max
的区间分为长度相同的n-1段，记相邻插值点间隔为δt，第i个插值点温度ti＝iδt，插值点i到i 1间为第i个插值区域，i＝1，2，3
…
n-1，插值区间ti～t
i 1
内f
h,i
(t)和如下：
[0045]
其中a
i,0
、a
i,1
、a
i,2
为插值系数，t为温度；
[0046]
(2)基于如下条件确定插值式f
h,i
(t)和中的插值系数：
[0047]
条件1：插值式f
h,i
(t)和在t∈[t
min
,t
max
]内连续，在插值点i处等于通过非平衡双原子气体电子、振动或转动模态比热计算式得到的比热值c
p
(ti)；
[0048]
条件2：f
h,i
(t)的一阶导数等于
[0049]
条件3：i＝1时，
[0050]
基于所述条件1～3，根据插值式f
h,i
(t)和建立的用于确定插值系数的方程如下：
[0051][0052]
进一步的，一种高温非平衡双原子气体热物性数据库生成方法中，还包括步骤：
[0053]
(3)在t∈(t
max
, ∞)范围内，定义：
[0054]
h＝f
h,n
(t)，
[0055]
且f
h,n
(tn)＝f
h,n-1
(tn)，
[0056]
其中，f
h,n-1
(tn)表示t
min
～t
max
温度范围内插值式f
h,i
(t)在第n-1个插值区域中第n个插值点的温度tn时的值；tn＝t
max
；
[0057]
(4)在t∈(0,t
min
)范围内，定义t1＝t
min
。
[0058]
进一步的，一种高温非平衡双原子气体热物性数据库生成方法步骤(3)中，在t∈(t
max
, ∞)温度范围内f
h,n
(t)和表达式如下：
[0059]fh,n
(t)＝-a/bexp[-b(t-tn)] c
[0060][0061]
其中，表达式中的参数a、b、c为：
[0062][0063][0064]
c＝f
h,n-1
(tn) a/bexp[-b(t-tn)]。
[0065]
进一步的，一种高温非平衡双原子气体热物性数据库生成方法中，所述t
min
＝50k，t
max
＝50000k，δt＝50k。
[0066]
本发明与现有技术相比具有如下至少一种有益效果：
[0067]
(1)本发明高温非平衡双原子气体热物性计算方法中，综合考虑了高阶电子能激发、振子非简谐性、转子非刚性、转动-振动耦合、离解能约束等高温真实气体效应影响，能够较好地反映高温气体的性质，本发明预测的焓值和比热值准确性明显提高；
[0068]
(2)本发明高温非平衡双原子气体热物性计算方法，在合理假设基础上，将双原子气体电子、振动、转动模态完全解耦，获得焓、比热近似值计算方法，能够支撑三、四温模型的工程应用探索，适用性强；
[0069]
(3)本发明高温非平衡双原子气体热物性数据库生成方法中，根据非平衡双原子气体比热计算式，采用分段线性插值和分段二次插值方法，建立了双原子气体焓与比热数据库，回避了高阶插值带来的runge现象，连续性好，有效提高了计算效率。
附图说明
[0070]
图1为本发明一种高温非平衡气体热物性计算方法和数据库生成方法的步骤图；
[0071]
图2为氮气(n2)在电子能基态(n＝0)不同转动量子数j下分子势能曲线；
[0072]
图3为实施例1中利用本发明方法、低温假设简化方法和理论预测的(a)n2、(b)co、(c)o2和(d)no四种典型双原子分子气体内自由度比热随温度变化曲线。
具体实施方式
[0073]
下面通过对本发明进行详细说明，本发明的特点和优点将随着这些说明而变得更为清楚、明确。
[0074]
在这里专用的词“示例性”意为“用作例子、实施例或说明性”。这里作为“示例性”所说明的任何实施例不必解释为优于或好于其它实施例。尽管在附图中示出了实施例的各种方面，但是除非特别指出，不必按比例绘制附图。
[0075]
目前双温以上热非平衡流场模拟中，低温假设简化方法是获取双原子气体的焓、比热数据的唯一途径，该方法的优点是：计算式具有解析形式，对内自由度模态解耦彻底，适用于三、四温模型；缺点是：考虑能级有限，不考虑模态耦合，对高温气体热物性预测偏差大，低温假设简化方法的假设如下：1)电子能最多激发一阶；2)分子振动简谐；3)分子为刚性转子。基于此得到的焓和比热计算式为：
[0076][0077][0078]hrot
＝rt
rot
ꢀꢀꢀꢀꢀꢀꢀꢀ
[0079]
其中，为电子基态与第一激发态能级简并度，为第一激发态特征温度，为电子基态的振动特征温度，r为气体常数。该方法能较好地反映低温双原子气体性质，但受限于简化假设，其对高温双原子气体热物性预测较大偏差，h代表焓，c
p
代表比热，上标el、rot、vib分别代表双原子分子的电子、转动、振动模态，各电子、振动、转动模态对应温度分别为t
ex
、t
vib
、t
rot
。
[0080]
对此，综合考虑高阶电子能激发、振子非简谐性、转子非刚性、转动-振动耦合、离解能约束等高温真实气体效应影响，完成了高温非平衡气体建模，提出了一种内自由度模态解耦假设，建立了高温非平衡双原子气体热物性计算方法，并据此设计了高温非平衡双原子气体热物性数据库生成方法。
[0081]
如图1所示，本发明一种高温非平衡双原子气体热物性计算方法包括如下步骤：
[0082]
步骤s1：高温非平衡气体建模；
[0083]
双原子分子振动、转动模态间存在强耦合，为建立非平衡双原子气体模型，首先将转动-振动能分为三部分：只受振动模态影响部分只受转动模态影响部分同时
受转动模态和振动模态影响部分
[0084][0085][0086][0087]
其中，h
plk
为普朗克常数，c为光速，n、v、j是电子、振动和转动量子数，是相应的电子、振动和转动能，ω
e,n
、ω
e,n
x
e,n
、ω
e,nye,n
、ω
e,nze,n
、α
e,n
、β
e,n
、b
e,n
、d
e,n
为经由实验测量或理论计算确定的双原子分子光谱常数。引入取值范围0～1的耦合能量分配控制参数f，把分为和两部分，分别分配到振动模态和转动模态中，得到新的振动能与转动能：
[0088][0089][0090]
本发明设计的f形式如下，由振动温度t
vib
、转动温度t
rot
和振动特征温度确定：
[0091][0092]
假设双原子分子在电子、振动、转动模态各能级上概率独立满足玻尔兹曼分布，记各电子、振动、转动模态对应温度分别为t
ex
、t
vib
、t
rot
，得到非平衡形式的配分函数，即基于配分函数的高温非平衡气体模型：
[0093][0094]
其中，q
int
是分子内能配分函数，σ是分子对称常数，k
bol
是玻尔兹曼常数；
[0095]nmax
，v
max,n
和j
max,n,v
分别为最大电子能级，最大振动能级和最大转动能级，分子种类确定时，最大电子能级n
max
为已知量，v
max,n
和j
max,n,v
通过步骤s2确定，v
max,n
为固定双原子分子电子量子数n得到的最大振动能级，j
max,n,v
为固定振动量子数v和电子量子数n得到的最大转动能级j
max,n,v
。
[0096]
步骤s2：振动能级和转动能级范围计算
[0097]
受转动能与振动能之和小于分子最小理解能约束，转动、振动量子数存在有限种组合，下面将以氮气(n2)电子能基态为例，阐述振动能级和转动能级范围确定的具体步骤。考虑转动的影响，n2分子内势能函数及其一阶、二阶导数为：
[0098][0099][0100][0101][0102]
表1 n2气体电子能基态主要光谱数据
[0103][0104]
注*表中除r
e,0
单位为外，其余单位为cm-1
[0105]
上式中，为转动能处于基态时分子离解所需的最小能量，r
e,n
对应转动能基态势能曲线最低点n原子间距离，μ
n2
＝m
nmn
/(mn mn)为折合分子质量，式中mn为单原子分子质量，表1给出了电子能基态下氮气的主要光谱数据。
[0106]
如图2，示出了电子基态下，不同转动量子数j下分子内势能曲线，图中r
min,n,j
与r
max,n,j
分别对应了势能曲线的极小值点与极大值点，r
max,n,j
对应的u(r
max,n,j
)为给定量子数n、j下分子最小离解能，即转动、振动能之和需满足本发明参照下列步骤确定能级范围：
[0107]
1)确定转动量子数上界j
lim.n
。根据函数在极值点处一阶导数为零的性质，得到方程：
[0108][0109]
从j＝0开始逐渐增大j，定义使方程无解的最小整数j为j
min,n
。
[0110]
具体的说，本步骤中，固定双原子分子电子量子数n，根据双原子分子光谱数据绘制分子势能曲线，逐渐增大转动量子数j，直至势能曲线不存在极值点，将此时的转动量子数定义为j
lim,n
。
[0111]
2)在0≤j《j
lim,n
范围内确定离解所需最小能量根据函数在极大值点处一阶导数在为零、二阶导数小于零的性质：
[0112][0113]
求解得到满足上式的r
max,n,j
，代入势能函数得到
[0132][0133]ai,0
、a
i,1
、a
i,2
是插值系数，下面将说明插值系数的确定过程。插值式fh(t)和在50k～50000k范围内满足下列条件1)在50k～50000k范围内连续，在插值点i处等于单一温度变量的非平衡气体比热计算式给出的比热值；2)fh(t)在50k～50000k范围内可导，一阶导数等于基于此建立下列方程：
[0134][0135]
上述方程共包含3n-2项、3(n-1)个未知数，假设在i＝1时得到使方程组定解的补充方程。
[0136]a1,0
a
1,1
t1 0.5a
1,2
t
12
＝(a
1,1
a
1,2
t1)t1[0137]
对数据库覆盖范围以外的情况，本发明采取以下处理方式：
[0138]
0k～50k：将c
p,1
外推作为该区域的比热，焓值在插值点1处连续，由此在该区域内：
[0139]
h＝c
p,1
t,c
p
≡c
p,1
[0140]
》50000k：在该区域内，假设h形式为-a/bexp[-b(t-50000)] c，c
p
形式为aexp[-b(t-50000)]，在插值点n处，h等于f
h,999
(t
1000
)、c
p
等于c
p
的一阶导数等于由此得到：
[0141][0142]
实施例1：
[0143]
本应用实例选取了n2、co、o2和no四种典型双原子分子气体，验证了本发明的高温非平衡气体热物性计算方法(以下简称“新算法”)对高温双原子气体比热预测有效性，并以n2气体的焓与比热计算为例，测试了现有技术中的较低温假设简化方法(以下简称“旧算法”)的效率提升情况。
[0144]
1)双原子分子气体比热计算测试
[0145]
验证测试选取n2、co、o2和no四种典型双原子分子，设定温度范围50k～50000k。
[0146]
如图3，示出了新旧算法对四种双原子气体比热的预测结果，以及分子物理理论计算的总比热数据。在整个测试温度范围内，新算法对总比热预测结果与理论数据相一致，旧算法预测的总比热在温度大于约7000k后显著偏离理论数据。进一步对比新旧算法对各内能模态比热的预测情况：首先对比的预测情况，因为新算法考虑电子能级数多，所以预测
的相对高；再对比的预测情况，在温度不超过约2000k时两种算法预测的相吻合，随温度升高，受振动非简谐性影响，新算法预测的略高于旧算法，随温度进一步上升，受最小离解能约束影响，新算法给出的趋于0，而旧算法预测的趋于气体常数r；最后对比的预测情况，基于刚性转子假设预测的恒等于r，在考虑转子非刚性和最小离解能约束后，新方法预测的在温度大于约15000k后逐渐趋于0。以上结果都说明了新算法显著优于旧算法。
[0147]
2)计算效率测试
[0148]
为说明基于新算法计算效率上的优势，设计了基于n2气体焓、比热计算对比测试。测试共进行五轮，每轮包含新旧算法对三种内能模态焓、比热计算测试，共12个算例，单个算例在温度50k～50000k的范围内，每间隔0.001k布置一个算点，每个算例包含算点49950001个。
[0149]
表3、表4记录了五轮测试结果，以及各参数求解的平均耗时。测试结果表明，新算法解算电子、振动模态焓的平均耗时约是旧算法的四分之一，解算电子、振动模态比热的平均耗时约是旧算法的六分之一，解算转动模态焓、比热的平均耗时略低于旧算法。
[0150]
表3 n2气体焓计算耗时
[0151][0152]
表4 n2气体比热计算耗时
[0153][0154]
以上结合具体实施方式和范例性实例对本发明进行了详细说明，不过这些说明并不能理解为对本发明的限制。本领域技术人员理解，在不偏离本发明精神和范围的情况下，可以对本发明技术方案及其实施方式进行多种等价替换、修饰或改进，这些均落入本发明的范围内。本发明的保护范围以所附权利要求为准。
[0155]
本发明说明书中未作详细描述的内容属本领域技术人员的公知技术。