采用贝叶斯估计方法在N-S/DSMC耦合算法中进行信息交换的策略与流程

采用贝叶斯估计方法在n-s/dsmc耦合算法中进行信息交换的策略
技术领域
1.本发明涉及飞行器空气动力技术领域，涉及飞行器高超声速绕流气动力热特性的数值模拟，针对飞行器在过渡流区n-s/dsmc耦合算法中两种计算方法耦合的信息交换的实现，具体涉及一种采用贝叶斯估计方法在n-s/dsmc耦合算法中进行信息交换的策略。


背景技术：

2.在研究飞行器绕流的空气动力学特征时，一般把气体稀薄程度按照kn(knudsen克努森数)(即气体分子平均自由程与流动特征长度比值)数值范围，将航天器经历的绕流问题大致分为：连续流区(kn《10-3
)、过渡流区(10-3
《kn《10)、自由分子流区(kn》10)。其中，以过渡流区的气动问题最为复杂，这一流动区域又可细分为近连续流、滑移流和稀薄过渡流。位于连续流与自由分子流之间的过渡区流动在数值计算方面是难于处理的一种流动。为了研究飞行器跨越各流域的空气动力学特征，传统做法是稀薄流用稀薄流的一套算法，如dsmc(直接模拟monte carlo)方法；而连续流有连续流的一套研究方法，如数值求解euler、n-s(纳维-斯托克斯)方程等；两类方法相差较大，彼此独立，且计算结果很难随高度光滑连接。
3.在过渡流区，n-s/dsmc耦合算法是一种求解途径。在一个流场中，把流场分为连续流区和连续流失效区域，连续流区采用n-s方程求解，连续流方程失效的区域采用dsmc方法，可以拓展n-s方程和dsmc方法的应用范围，提高dsmc方法的计算效率。
4.在两种方法进行耦合计算的过程中，需要进行信息交换。n-s求解区域对dsmc求解区域提供边界信息时，相对简单，根据当地参数抽样出每一个进入dsmc区域的仿真分子的参数。而dsmc为n-s提供边界条件时，由于dsmc瞬时结果存在很大的统计波动，会对n-s求解过程的稳定性产生影响，严重时导致求解发散，不能继续进行。
5.现有技术中一般采用的“亚松弛”技术来实现，其耦合计算结构能满足试验要求，但其缺点在于信息交换的效率太低，影响计算进程。


技术实现要素：

6.本发明的一个目的是解决至少上述问题和/或缺陷，并提供至少后面将说明的优点。
7.为了实现本发明的这些目的和其它优点，提供了一种采用贝叶斯估计方法在n-s/dsmc耦合算法中进行信息交换的策略，在基于贝叶斯估计方法的dsmc计算区域对n-s求解区域进行信息交换的流程被配置为包括：
8.s1、利用dsmc统计分子速度样本，采用正态-逆伽玛共轭先验分布，把n-s上一步计算的速度和移动温度作为先验，估计速度和气体的移动温度t
tr
；
9.s2、利用dsmc统计转动能量样本、振动能量样本，采用正态-逆伽玛共轭先验分布，把n-s上一步计算的转动温度值、振动温度值分别作为先验，估计转动温度t
rot
、振动温度
t
vib
；
10.s3、利用dsmc统计密度样本，采用正态共轭先验分布，把n-s上一步计算的密度值作为先验，估计密度；
11.s4、将s41-s43中得到的估计量作为dsmc的结果与n-s计算进行信息交换。
12.优选的是，在s1中，所述速度和气体的移动温度t
tr
估算流程包括：
13.以x方向为例，将分子的速度整体分布写成如下的常规正态分布：
[0014][0015]
其中，r为气体常数，t
trx
为x方向移动温度，θ为估计均值，σ2为方差，u为x方向的分子速度，为x方向的分子平均速度；
[0016]
在已知某气体组份n个速度样本x(x1，
…
xn)，通过估计均值θ和方差σ2，从而得到t
tr
；
[0017]
θ、σ2的先验分布采用如下的正态—逆伽玛共轭先验分布：
[0018][0019]
其中，k0、υ0为ai中需要通过调试确定的超参数；θ0、为先验均值和方差；
[0020]
(θ，σ2)联合后验密度函数为：
[0021][0022]
其中，
[0023]
vn＝n+υ0[0024]kn
＝n+k0[0025][0026][0027]
其中，n为样本量；通过调试确定k0、后，可得vn，kn、θn、υn；
[0028]
θ、σ2对应的最大后验概率贝叶斯估计分别为：
[0029][0030][0031]
其中，为样本均值：
[0032][0033]
s2为样本方差：
[0034][0035]
通过n-s计算的速度u
n-s
和温度t
trn-s
得到先验均值和方差θ0、
[0036]
90＝u
n-s
[0037][0038]
通过对θ、σ2的估计，可得x方向的和移动温度t
trx
：
[0039][0040][0041]
采用同样的方法对y，z方向进行贝叶斯估计，得到y方向速度和温度t
try
，z方向速度和温度t
trz
，则气体的移动温度t
tr
为：
[0042][0043]
k0、v0作为超参数，通过调试，给出适合高超声速流动的参数设置；
[0044]
优选的是，在s2中，所述转动温度t
rot
的估计流程包括：
[0045]
将转动能量写成如下的常规伽玛分布：
[0046][0047]
其中，γ(r)为伽玛函数，λ、r为伽马函数中的参数，ζ
rot
为转动自由度，k为玻尔兹曼常数，ε为振动能；
[0048]
λ共轭先验分布为γ(α，β)
[0049][0050]
其中，α、β为伽玛函数中的参数；
[0051]
对于气体组份中的n个转动能量样本x(x1，
…
xn)，λ最大后验概率贝叶斯估计为：
[0052][0053]
其中：
[0054][0055][0056]
其中，a为待定超参数，t
rotn-s
为通过n-s计算的转动温度；
[0057]
则转动温度t
rot
通过下式获得：
[0058][0059]
所述振动温度t
vib
的估计方式与转动温度t
rot
估计方式配置为一致。
[0060]
优选的是，在s3中，密度的估计流程被配置为包括：
[0061]
将密度写成如下的常规正态分布：
[0062][0063]
基于气体组份n个密度样本x(x1，
…
xn)，估计均值θ，从而得到密度均值
[0064]
采用s41中相同的共轭先验分布和贝叶斯估计，基于下式得到：
[0065][0066]
其中，θ0为通过n-s计算的密度值ρ
n-s
得到的先验均值，
[0067]
θ0＝ρ
n-s
[0068]
k0作为超参数，通过调试，给出适合高超声速流动的参数设置，这里，k0的值与。
[0069]
优选的是，在采用贝叶斯估计方法的dsmc计算区域对n-s求解区域进行信息交换之前，还包括：
[0070]
步骤一，基于连续流区求解化学非平衡流动n-s方程和基于稀薄流的化学非平衡dsmc计算；
[0071]
步骤二，基于n-s方程计算结果求解当地克努森数，并基于dsmc计算结果把流场分为连续流区域和稀薄流区域，以得到对应的边界面；
[0072]
步骤一，基于n-s方程计算结果对边界面上的非平衡分布进行计算，并据非平衡分布抽样进入dsmc计算区域仿真分子参数，为dsmc计算提供边界条件。
[0073]
本发明至少包括以下有益效果：
[0074]
本发明为了抑制dsmc统计波动对n-s求解过程的影响，采用了一种人工智能的贝叶斯估计算法替代现有技术中一般采用的“亚松弛”技术，能够在保证计算结果、精度具有一致性的情况下，在较少的dsmc采样样本量条件下，给出较好的参数估计，提供给n-s求解，有效提升信息交换效率。
[0075]
(1)本发明采用共轭先验分布，在dsmc向n-s信息交换过程中，避免了贝叶斯估计中的大量计算。
[0076]
(2)本发明把n-s计算的上一步结果作为先验取值，方便了贝叶斯估计过程中对先验的需求。
[0077]
(3)本发明采用统一网格系统，便于计算过程中的自动分区的判断和信息交换过程的计算。
[0078]
(4)本发明采用自动分区技术，在计算过程中判断哪些区域处于连续流方程失效的区域，适应性强。
[0079]
(4)本发明采用非平衡边界，在n-s向dsmc信息交换过程中，考虑非平衡效应的影响。
[0080]
本发明的其它优点、目标和特征将部分通过下面的说明体现，部分还将通过对本发明的研究和实践而为本领域的技术人员所理解。
附图说明
[0081]
图1为本发明中n-s/dsmc耦合计算流程图；
[0082]
图2为本发明中自动分区示意图；
[0083]
图3为本发明中分区判断示意图；
[0084]
图4为本发明分界面判断示意图；
[0085]
图5为本发明与现有技术的曲线式对比示意图。
具体实施方式
[0086]
下面结合附图对本发明做进一步的详细说明，以令本领域技术人员参照说明书文字能够据以实施。
[0087]
为了抑制dsmc统计波动对n-s求解过程的影响，本发明采用了一种人工智能的贝叶斯估计算法替换一般采用“亚松弛”技术，贝叶斯(bayes)算法是基于贝叶斯定理的一种人工智能算法，是基于后验概率分布的一类估计方法，其中后验概率分布中采用了先验信息和样本信息。其基本思想是：把待估计参数θ看作具有先验分布密度l(θ)的随机变量，根据样本集x，给出后验概率分布π(θ|x)，基于后验概率分布，估计最优的θ值。采用贝叶斯估计方法法，在较少的dsmc采样样本量条件下，给出较好的参数估计，提供给n-s求解。
[0088]
实施例：
[0089]
一种适用于n-s/dsmc耦合算法的人工智能贝叶斯方法信息交换技术，如图1，其过程包括：
[0090]
步骤一、基于连续流区求解n-s方程的计算和基于稀薄流的dsmc计算，其包括：
[0091]
s11、基于连续流区的n-s方程数值计算，控制方程为(热)化学非平衡流动n-s方程，采用gordon和mcbride的拟合公式计算气体热力学特性，采用gupta和yos粘性系数拟合公式、eucken关系和wike混合规则计算气体输运特性，采用dunn&kang化学动力学模型、millikan和white的振动松弛公式计算热化学源项。网格系统与dsmc相同，控制方程数值离散采用限体积法、steger-warming分裂、nnd格式、dplur隐格式等数值方法。每个迭代步，n-s方程获得组元质量分数、速度、温度、压力等流场数据后与dsmc进行信息交换。
[0092]
s12、基于稀薄流区的dsmc方法计算，计算中采用变刚球分子模型(vhs)，拉尔森-博格内克(larsen-borgnakke)碰撞模型计算分子碰撞过程的能量传递。采用n-s计算中相同的化学反应模型，计算空气5组元热化学非平衡过程。温度采用三温度模型(移动温度、转动温度、振动温度)，其中移动温度和转动温度按自由度平均后作为一个温度与n-s计算进行信息交换。网格系统采用二级直角网格系统：表面网格a(或称为面元网格)为三角形非结构网格；空间背景网格b(或称为空间网格)为等距笛卡尔网格(如图2)，碰撞子网格根据流场密度自适应，碰撞对的选取限制在碰撞子网格内。
[0093]
步骤二、基于当地克努森数的自动分区，其包括：
[0094]
s21、根据n-s计算结果，计算当地克努森数。当地克努森数表示流场内某网格内稀薄程度
[0095][0096]
其中：λ为当地流场气体分子平均自由程；q为当地流场宏观参数，可以是压力、密度或温度。
[0097]
s22、扫描每个网格，判断是属于哪个分区。分别采用温度、压力、和温度计算当地克努森数。取梯度最大的流动参数来计算kn
l
。当kn
l
≥0.02，认为连续流方程失效，该网格采用dsmc方法进行模拟。在计算过程中，根据n-s计算结果不断调整dsmc方法计算区域，直到流场达到稳定，如图3中c处为n-s区域，d处为界面，e处为dsmc区域。
[0098]
s23、判断边界面。扫描每个网格及其周围的网格，如果一个网格体是dsmc计算网格，而其相邻的网格体为非dsmc网格，则这两个网格中间的界面为分界面(即图4的f处)。
[0099]
步骤三、基于当地弱非平衡态的n-s求解区对dsmc计算区域的信息交换。
[0100]
s31、利用n-s计算的结果，计算边界上非平衡分布。根据连续流失效判断方法，当k
nl
≥0.02时，在耦合边界上，流动已经开始出现非平衡现象，此时速度分布函数已经不再符合maxwell平衡分布，可以采用chapman-enskog分布来描述。
[0101]
对于偏离maxwell分布的非平衡分布，其一阶展开形式为：
[0102]
f(c)＝f0(c)γ(c)
[0103]
其中，c为无量纲化的分子热运动速度
[0104]
c＝c/(2kt/m)
1/2
[0105]
其中，c为分子热运动速度，k为玻尔兹曼常数，t为气体温度，m为气体分子质量。
[0106][0107][0108][0109][0110]
其中，c
x,y,z
为x、y、z三个方向无量纲化的分子热运动速度，τ为剪切应力，q为热流，κ为热传导系数，p为压力，μ为粘性系数。
[0111]
s32、根据非平衡分布确定进入dsmc计算区域仿真分子参数。在边界面上给出符合chapman-enskog分布的分子速度。基于maxwell分布的速度采样方法，根据当地的非平衡度，采用“拒绝-接受”的随机步骤，产生符合chapman-enskog分布的分子速度。步骤如下：
[0112]
(1)根据n-s方程结果，计算出当地的τ
i,j
和qi，取
[0113]
b≡max(|τ
i,j
|,|qi|)
[0114]
(2)设置幅度参数
[0115]
a＝1+30b
[0116]
(3)根据maxwell分布产生一个速度矢量c
try
。
[0117]
(4)产生一个随机数rf，如果
[0118]
arf≤γ(c)
[0119]
接受c
try
，否则，回到(3)产生新的c
try
。
[0120]
(5)给出分子速度
[0121]
c＝(2kt/m)
1/2ctry
+u
[0122]
其中，u为边界面上的宏观速度分量。
[0123]
步骤四、基于贝叶斯估计的dsmc计算区域对n-s求解区域的信息交换。其包括：
[0124]
s41、利用dsmc统计分子速度样本，采用正态-逆伽玛共轭先验分布，把n-s上一步计算的速度和移动温度作为先验，估计速度和移动温度t
tr
。
[0125]
以x方向为例，分子的速度整体分布为正态分布：
[0126][0127]
其中，u为x方向的分子速度，为x方向分子平均速度，r为气体常数，t
trx
为x方向移动温度。
[0128]
写成常规的正态分布：
[0129][0130]
已知某气体组份n个速度样本x(x1，
…
xn)，估计均值θ和方差σ2，从而得到t
tr
。
[0131]
θ、σ2的先验分布采用共轭先验分布：正态—逆伽玛分布：
[0132][0133]
其中，k0、υ0为ai中需要通过调试确定的超参数；θ0、为先验均值和方差；
[0134]
(θ，σ2)联合后验密度函数为
[0135][0136]
其中，
[0137]
vn＝n+υ0[0138]kn
＝n+k0[0139][0140][0141]
(θ，σ2)的最大后验概率贝叶斯估计为：
[0142]
[0143][0144]
其中，样本均值
[0145][0146]
样本方差
[0147][0148]
通过n-s计算的速度u
n-s
和温度t
trn-s
得到先验均值和方差θ0、
[0149]
θ0＝u
n-s
[0150][0151]
通过对θ、σ2的估计，可得x方向的和移动温度t
trx
：
[0152][0153][0154]
其他方向(y，z)采用同样的方法进行贝叶斯估计，得到y方向速度和温度t
try
，z方向速度和温度t
trz
。
[0155]
气体的移动温度t
tr
为
[0156][0157]
k0、υ0作为超参数，通过调试，给出适合高超声速流动的参数设置。
[0158]
s42、利用dsmc统计转动能量样本，采用伽玛共轭先验分布，把n-s上一步计算的转动温度作为先验，估计转动温度t
rot
。
[0159]
转动能量分布满足伽玛分布
[0160][0161]
其中，ε为转动能，ζ
rot
为转动自由度，t
rot
为转动温度。
[0162]
写成常规的伽玛分布
[0163][0164]
λ共轭先验分布为γ(α，β)
[0165]
[0166]
已知某气体组份n个转动能量样本x(x1，
…
xn)
[0167]
可得，λ最大后验概率贝叶斯估计为：
[0168][0169]
其中：
[0170][0171][0172]
其中，a为超参数，通过调试，给出适合高超声速流动的参数设置。t
rotn-s
为通过n-s计算的转动温度。
[0173]
可得转动温度
[0174][0175]
s43、利用dsmc统计振动能量样本，采用伽玛共轭先验分布，把n-s上一步计算的振动动温度作为先验，估计振动温度t
vib
。
[0176]
振动能量分布与转动能量的形式完全一致，满足伽玛分布
[0177][0178]
其中，ε为振动能，ζ
vib
为振动自由度，t
vib
为振动温度。
[0179]
由于振动能量分布与转动能量的形式完全一致，可以采用相同的贝叶斯估计方法，已知某气体组份n个振动能量样本x(x1，
…
xn)
[0180]
可得，λ最大后验概率贝叶斯估计为：
[0181][0182]
其中：
[0183][0184][0185]
其中，a为超参数，通过调试，给出适合高超声速流动的参数设置。t
vibn-s
为通过n-s计算的转动温度。振动自由度通过t
vibn-s
求出，
[0186]
[0187]
其中，θ
vib
为主份的振动特征温度。
[0188]
s44、利用dsmc统计密度样本，采用正态共轭先验分布，把n-s上一步计算的密度作为先验，估计密度。
[0189]
设密度符合正态分布
[0190][0191]
写成常规的正态分布：
[0192][0193]
已知某个气体组份n个密度样本x(x1，
…
xn)，估计均值θ，从而得到
[0194]
采用s41中相同的共轭先验分布和贝叶斯估计，可得
[0195][0196]
其中，通过n-s计算的密度ρ
n-s
得到先验均值θ0[0197]
θo＝p
n-s
[0198]
k0作为超参数，通过调试，给出适合高超声速流动的参数设置，这里，k0的值与s41中的值不同。
[0199]
s45、将以上估计量作为dsmc的结果与n-s计算进行信息交换。
[0200]
如图5所示，通过本方案的人工智能的贝叶斯估计算法替换一般采用“亚松弛”技术，其达到的效果在于：在温度从300k到1300k的变化过程中，贝叶斯估计约30步可以达到目标温度，亚松弛需要80步达到目标温度，提高了信息交换的效率。
[0201]
以上方案只是一种较佳实例的说明，但并不局限于此。在实施本发明时，可以根据使用者需求进行适当的替换和/或修改。
[0202]
这里说明的设备数量和处理规模是用来简化本发明的说明的。对本发明的应用、修改和变化对本领域的技术人员来说是显而易见的。
[0203]
尽管本发明的实施方案已公开如上，但其并不仅仅限于说明书和实施方式中所列运用。它完全可以被适用于各种适合本发明的领域。对于熟悉本领域的人员而言，可容易地实现另外的修改。因此在不背离权利要求及等同范围所限定的一般概念下，本发明并不限于特定的细节和这里示出与描述的图例。

技术特征：
1.一种采用贝叶斯估计方法在n-s/dsmc耦合算法中进行信息交换的策略，其特征在于，在基于贝叶斯估计方法的dsmc计算区域对n-s求解区域进行信息交换的流程被配置为包括：s1、利用dsmc统计分子速度样本，采用正态-逆伽玛共轭先验分布，把n-s上一步计算的速度和移动温度作为先验，估计速度和气体的移动温度t
tr
；s2、利用dsmc统计转动能量样本、振动能量样本，采用正态-逆伽玛共轭先验分布，把n-s上一步计算的转动温度值、振动温度值分别作为先验，估计转动温度t
rot
、振动温度t
vib
；s3、利用dsmc统计密度样本，采用正态共轭先验分布，把n-s上一步计算的密度值作为先验，估计密度；s4、将s41-s43中得到的估计量作为dsmc的结果与n-s计算进行信息交换。2.如权利要求1所述的采用贝叶斯估计方法在n-s/dsmc耦合算法中进行信息交换的策略，其特征在于，在s1中，所述速度和气体的移动温度t
tr
估算流程包括：以x方向为例，将分子的速度整体分布写成如下的常规正态分布：其中，r为气体常数，t
trx
为x方向移动温度，θ为估计均值，σ2为方差，u为x方向的分子速度，为x方向的分子平均速度；在已知某气体组份n个速度样本x(x1，
…
x
n
)，通过估计均值θ和方差σ2，从而得到θ、σ2的先验分布采用如下的正态—逆伽玛共轭先验分布：其中，k0、v0为ai中需要通过调试确定的超参数；θ0、为先验均值和方差；(θ，σ2)联合后验密度函数为：其中，v
n
＝n+υ0k
n
＝n+k
00
其中，n为样本量；通过调试确定、k0、υ0后，可得v
n
，k
n
、θ
n
、υ
n
；θ、σ2对应的最大后验概率贝叶斯估计分别为：
其中，样本均值其中，样本均值样本方差s2：通过n-s计算的速度u
n-s
和温度t
trn-s
得到先验均值和方差θ0、θ0＝u
n-s
通过对θ、σ2的估计，可得x方向的和移动温度t
trx
：：采用同样的方法对y，z方向进行贝叶斯估计，得到y方向速度和温度t
try
，z方向速度和温度t
trz
，则气体的移动温度t
tr
为：k0、υ0作为超参数，通过调试，给出适合高超声速流动的参数设置。3.如权利要求1所述的采用贝叶斯估计方法在n-s/dsmc耦合算法中进行信息交换的策略，其特征在于，在s2中，所述转动温度t
rot
的估计流程包括：将转动能量写成如下的常规伽玛分布：其中，γ(r)为伽玛函数，λ、r为伽马函数中的参数，ζ
rot
为转动自由度，，k为玻尔兹曼常数，ε为振动能；λ共轭先验分布为γ(α，β)：其中，α、β为伽玛函数中的参数；对于气体组份中的n个转动能量样本x(x1，...x
n
)，λ最大后验概率贝叶斯估计为：其中：
其中，a为待定超参数，t
rotn-s
为通过n-s计算的转动温度；则转动温度t
rot
通过下式获得：所述振动温度t
vib
的估计方式与转动温度t
rot
估计方式配置为一致。4.如权利要求1所述的采用贝叶斯估计方法在n-s/dsmc耦合算法中进行信息交换的策略，其特征在于，在s3中，密度的估计流程被配置为包括：将密度写成如下的常规正态分布：x＝ρ基于气体组份n个密度样本x(x1，...x
n
)，估计均值θ，从而得到密度均值采用s41中相同的共轭先验分布和贝叶斯估计，基于下式得到：其中，θ0为通过n-s计算的密度ρ
n-s
得到的先验均值，θ0＝ρ
n-s
k0为待定超参数。5.如权利要求1所述的采用贝叶斯估计方法在n-s/dsmc耦合算法中进行信息交换的策略，其特征在于，在采用贝叶斯估计方法的dsmc计算区域对n-s求解区域进行信息交换之前，还包括：步骤一，基于连续流区求解化学非平衡流动n-s方程和基于稀薄流的化学非平衡dsmc计算；步骤二，基于n-s方程计算结果求解当地克努森数，并基于dsmc计算结果把流场分为连续流区域和稀薄流区域，以得到对应的边界面；步骤一，基于n-s方程计算结果对边界面上的非平衡分布进行计算，并据非平衡分布抽样进入dsmc计算区域仿真分子参数，为dsmc计算提供边界条件。

技术总结
本发明公开了一种采用贝叶斯估计方法在N-S/DSMC耦合算法中进行信息交换的策略，包括：在基于贝叶斯估计方法的DSMC计算区域对N-S求解区域进行信息交换中，利用DSMC统计样本，估计速度气体的移动温度T


技术研发人员：李中华 党雷宁 吴俊林 彭傲平
受保护的技术使用者：中国空气动力研究与发展中心超高速空气动力研究所
技术研发日：2023.05.24
技术公布日：2023/8/21