一种用于二维稀薄气体流动的数值计算方法技术

本发明专利技术属于稀薄气体动力学技术领域，具体涉及一种用于二维稀薄气体流动的数值计算方法。具体技术方案为：将DG数值方法引入NCCR方程中，进行数值离散计算，且引入限制器进行数值间断侦测和限制。该方法解决了传统NS方程无法解决非连续流动、分子动力模型在过渡流和滑移流动中计算量大的技术问题，对于临近空间各种飞行器的超高速流动及热力计算尤其是激波问题求解有着很大的优势。问题求解有着很大的优势。问题求解有着很大的优势。

【技术实现步骤摘要】
一种用于二维稀薄气体流动的数值计算方法


[0001]本专利技术属于稀薄气体动力学
，具体涉及一种用于二维稀薄气体流动的数值计算方法。

技术介绍

[0002]传统流体力学方程的研究集中于连续区域，即流体分子的平均自由程远小于流动特征尺寸，纳维
‑
斯托克斯
‑
傅里叶方程(NSF方程)是该类方程的典型代表。其通过牛顿粘性应力和傅里叶热传导的本构关系耦合斯托克斯假设来封闭连续、动量、能量三大流动守恒方程。近两个世纪，NSF方程在连续流领域取得了巨大成功，推动了流体力学的发展。然而，当分子的平均自由程逐渐增大，即气体变得稀薄而不再稠密时，连续性假设失效，这时候如果继续采用基于连续介质或平衡态(当系统的宏观热观察量不再随时间改变，此时系统处于平衡态)的NSF方程必然是不合适的。
[0003]近年来，随着科学的不断进步和人类对未知领域的探索，使得人们越来越关注稀薄气体领域及其工程应用，如在空气稀薄的大气层(几十公里高空外)中飞行的飞行器的气动力的预测。该条件下气体密度很低，大概为10
‑7到10
‑
10
kg/m3，此时分子间距较大会导致一个分子碰撞两次所走过的距离较大，即分子的平均自由程较大，由于分子的碰撞效应逐渐显现出来，因此传统的宏观方程—NSF方程将不再适用，此时带有碰撞项的Boltzmann方程担任起了解决稀薄流的重任。
[0004]根据克努森K
n
数(衡量气体稀薄程度以及宏观模型正确程度衡量气体稀薄程度以及宏观模型正确程度的重要参数)的大小，现有的气体流动及其相应求解方法可分为四种：
[0005]连续流领域(Kn<0.001)，传统的欧拉方程或者NS方程均适用；
[0006]滑移流领域(0.001<Kn<0.1)，边界外的主流中，NS方程成立，但在边界处要考虑速度滑移和温度跳跃边界条件；
[0007]过渡流领域(0.1<Kn<10)，NS方程不再成立；
[0008]自由分子流(10>Kn)，直接采用分子运动论的方法，如DSMC。
[0009]上述各种求解流体问题的方法中，NS方程和分子运动论的方法分别在连续流和自由分子流领域取得了成功，但对于滑移流流域和过渡流流域而言，并没有一个统一的有效的控制方程对该区域的流域进行求解，例如近些年大力发展的临界空间层(海拔高度20
‑
100km)中的气体流动就正好处在这两个流域之中，因此发展一套既能适用于连续流又能适用于稀薄流的气体动力学理论就显得格外重要。
[0010]为研究稀薄气体流动，大量数学物理方程被提出，其基本可分为三类：从微观角度进行求解的分子动力模型，包括直接模拟的蒙特卡罗方法(DSMC)、DSMC的信息保存方法(DSMC
‑
IP)和格子玻尔兹曼方法(LBM)；也有对Boltzmann方程中的分布函数直接进行离散建模求解的方法，例如BGK模型方程；以及将宏观的统计表达代入到微观Boltzmann方程中，进而获得宏观量的守恒方程和演化方程，其代表性方法为Grad方法、Chapman
‑
Enskog方法和Eu方法。
[0011]对于分子动力模型，由于其采用模拟分子模拟真实分子的运动和气体流的碰撞，因此其在低Kn数(分子浓度较大)下会花费大量的计算资源。而对于BGK模型方程，其简化了Boltzmann方程中的碰撞项，因而其仅在平衡态或者近平衡态时能取得准确的结果，且其对于输运系数的计算也不准确。如通过BGK模型计算出来的单原子气体的普朗特(Pr)数为1，而不是正确值2/3。而在通过Boltzmann方程演化而来的非NS流体控制方程方法中，以Grad提出的矩方法为基础的R
‑
13方程和以Chapman
‑
Enskog展开为基础的Burnett方程，被认为不满足热力学定律的熵条件，即难以满足Boltzmann方程的定解，故此两类方法的发展收到了一定程度的限制。与该两种方法不同，Eu方法严格满足熵条件其从Boltzmann定解条件H定理出发，抓住了非平衡态到平衡态熵增的特点，构造了一种指数形式的分布函数，得到了非平衡态到平衡态的熵增耗散模型，该模型在接近平衡态时会收敛为Rayleigh
‑
Onsager耗散函数，以此为基础，通过构建非平衡态到平衡态统一的非线性熵增模型来处理碰撞项。由此，通过Boltzmann方程导出三大守恒方程和粘性应力与热通量等的本构方程，并将此守恒量和非守恒量的方程耦合在一起从而完成气体动力学方程的封闭形成了Eu统一流体方程。此后，Myong在Eu的基础上重点处理了本构方程的高阶项并以此为基础构建了非线性耦合本构关系(NCCR)。
[0012]因此，如果能够提供一种用于二维稀薄气体流动的数值计算方法，将具有优越的应用前景。

技术实现思路

[0013]为解决上述技术问题，本专利技术基于Eu方法推导而来的二维NCCR方程，给出了一种适用于稀薄气体流动的数值计算方法，尤其是带有激波现象的超高声速飞行的非平衡热流动分析。
[0014]为实现上述专利技术目的，本专利技术所采用的技术方案是：一种用于二维稀薄气体流动的数值计算方法，将DG数值方法引入NCCR方程中，进行数值离散计算，且引入限制器进行数值间断侦测和限制。
[0015]优选的：所述NCCR方程为通过无量纲参数、相似准则数简化得到的无量纲形式方程。
[0016]优选的：通过近似解S
h
、U
h
来表达局部单元Ω内的全局解S和U，
[0017][0018][0019]式中，N为基函数的个数；为基函数，
[0020]所述DG
‑
NCCR离散方程为：
[0021][0022]优选的：所述限制器为TVB斜率限制器。
[0023]优选的：包括如下步骤：
[0024]S01、通过网格程序对二维求解流场进行网格划分，并生成网格文件；
[0025]S02、读取网格文件，并记录网格节点坐标；
[0026]S03、根据节点坐标，获得标准正方形网格单元；
[0027]S04、给定压力远场边界条件、Langmuir壁面滑移边界条件；
[0028]S05、初始化计算条件，根据CFL条件确定时间步长；根据边界条件给流场参数赋初值x0，x
n
‑1＝x0；
[0029]S06、将x
n
‑1代入DG
‑
NCCR离散方程中的数值积分项中进行计算；
[0030]S07、采用不同的通量计算格式计算所述离散方程中的通量；
[0031]S08、根据步骤S06、S07的计算结果，采用三阶TVD
‑
RK对所述离散方程进行求解，得到流场参数x
n
。
[0032]优选的：本文档来自技高网...

【技术保护点】

【技术特征摘要】
1.一种用于二维稀薄气体流动的数值计算方法，其特征在于：将DG数值方法引入NCCR方程中，进行数值离散计算，且引入限制器进行数值间断侦测和限制。2.根据权利要求1所述的数值计算方法，其特征在于：所述NCCR方程为通过无量纲参数、相似准则数简化得到的无量纲形式方程。3.根据权利要求2所述的数值计算方法，其特征在于：通过近似解S
h
、U
h
来表达局部单元Ω内的全局解S和U，Ω内的全局解S和U，式中，N为基函数的个数；为基函数，所述DG
‑
NCCR离散方程为：4.根据权利要求1所述的数值计算方法，其特征在于：所述限制器为TVB斜率限制器。5.根据权利要求3所述的数值计算方法，其特征在于：包括如下步骤：S01、通过网格程序对二维求解流场进行网格划分，并生成网格文件；S02、读取网格文件，并记录网格节点坐标；S03、根据节点坐标，获得标准正方形网格单元；S04、给定压力远场边界条件、Langmuir壁面滑移边界条件；S05、初始化计算条件，根据CFL条件确定时间步长；根据边界条件给流场参数赋初值x0，x
n
‑1＝x0；S06、将x
n
‑1代入DG
‑
NCCR离散方程中的数值积分项中进行计算；S07、采用不同的通量计算格式计算所述离散方程中的通量；S08、根据步骤S06、S07的计算结果，采用三阶TVD
‑
RK对所述离散方程进行求解，得到流场参数x
n
。6.根据权利要求5所述的数值计算方法，其特征在于：所述步骤S01中，生成的二...

【专利技术属性】
技术研发人员：唐轲，肖洪，于艾洋，周磊，陈果，
申请(专利权)人：西北工业大学，
类型：发明
国别省市：