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

技术领域

本发明属于稀薄气体动力学技术领域，具体涉及一种用于二维稀薄气体流动的数值计算方法。

背景技术

传统流体力学方程的研究集中于连续区域，即流体分子的平均自由程远小于流动特征尺寸，纳维-斯托克斯-傅里叶方程(NSF方程)是该类方程的典型代表。其通过牛顿粘性应力和傅里叶热传导的本构关系耦合斯托克斯假设来封闭连续、动量、能量三大流动守恒方程。近两个世纪，NSF方程在连续流领域取得了巨大成功，推动了流体力学的发展。然而，当分子的平均自由程逐渐增大，即气体变得稀薄而不再稠密时，连续性假设失效，这时候如果继续采用基于连续介质或平衡态(当系统的宏观热观察量不再随时间改变，此时系统处于平衡态)的NSF方程必然是不合适的。

近年来，随着科学的不断进步和人类对未知领域的探索，使得人们越来越关注稀薄气体领域及其工程应用，如在空气稀薄的大气层(几十公里高空外)中飞行的飞行器的气动力的预测。该条件下气体密度很低，大概为10

根据克努森K

连续流领域(Kn<0.001)，传统的欧拉方程或者NS方程均适用；

滑移流领域(0.001

过渡流领域(0.1

自由分子流(10>Kn)，直接采用分子运动论的方法，如DSMC。

上述各种求解流体问题的方法中，NS方程和分子运动论的方法分别在连续流和自由分子流领域取得了成功，但对于滑移流流域和过渡流流域而言，并没有一个统一的有效的控制方程对该区域的流域进行求解，例如近些年大力发展的临界空间层(海拔高度20-100km)中的气体流动就正好处在这两个流域之中，因此发展一套既能适用于连续流又能适用于稀薄流的气体动力学理论就显得格外重要。

为研究稀薄气体流动，大量数学物理方程被提出，其基本可分为三类：从微观角度进行求解的分子动力模型，包括直接模拟的蒙特卡罗方法(DSMC)、DSMC的信息保存方法(DSMC-IP)和格子玻尔兹曼方法(LBM)；也有对Boltzmann方程中的分布函数直接进行离散建模求解的方法，例如BGK模型方程；以及将宏观的统计表达代入到微观Boltzmann方程中，进而获得宏观量的守恒方程和演化方程，其代表性方法为Grad方法、Chapman-Enskog方法和Eu方法。

对于分子动力模型，由于其采用模拟分子模拟真实分子的运动和气体流的碰撞，因此其在低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)。

因此，如果能够提供一种用于二维稀薄气体流动的数值计算方法，将具有优越的应用前景。

发明内容

为解决上述技术问题，本发明基于Eu方法推导而来的二维NCCR方程，给出了一种适用于稀薄气体流动的数值计算方法，尤其是带有激波现象的超高声速飞行的非平衡热流动分析。

为实现上述发明目的，本发明所采用的技术方案是：一种用于二维稀薄气体流动的数值计算方法，将DG数值方法引入NCCR方程中，进行数值离散计算，且引入限制器进行数值间断侦测和限制。

优选的：所述NCCR方程为通过无量纲参数、相似准则数简化得到的无量纲形式方程。

优选的：通过近似解S

式中，N为基函数的个数；

所述DG-NCCR离散方程为：

/>

优选的：所述限制器为TVB斜率限制器。

优选的：包括如下步骤：

S01、通过网格程序对二维求解流场进行网格划分，并生成网格文件；

S02、读取网格文件，并记录网格节点坐标；

S03、根据节点坐标，获得标准正方形网格单元；

S04、给定压力远场边界条件、Langmuir壁面滑移边界条件；

S05、初始化计算条件，根据CFL条件确定时间步长；根据边界条件给流场参数赋初值x

S06、将x

S07、采用不同的通量计算格式计算所述离散方程中的通量；

S08、根据步骤S06、S07的计算结果，采用三阶TVD-RK对所述离散方程进行求解，得到流场参数x

优选的：所述步骤S01中，生成的二维网格为三角形网格单元；所述步骤S03中，先将步骤S01中的三角形网格单元转换成标准三角形单元，再转换成标准正方形网格单元。

优选的：还包括：

S09、遍历所有求解结果，查找“问题单元”并采用限制器对其单元内的求解参数进行限制，更新x

S10、计算迭代误差x

当x

当x

S11、循环计算结束。

相应的：用于二维稀薄气体流动的数值计算方法在临近空间飞行器飞行中的应用。

相应的：一种电子设备，其特征在于：包括：

一个或多个处理器；

存储装置，用于存储一个或多个程序；

当所述一个或多个程序被所述一个或多个处理器执行，使得所述一个或多个处理器实现用于二维稀薄气体流动的数值计算方法。

相应的：一种计算机可读介质，所述可读介质存储有计算机程序，其特征在于：所述计算机程序被处理器执行时实现用于二维稀薄气体流动的数值计算方法。

与现有技术相比，本发明具有以下有益效果：

本发明将DG方法引入到NCCR方程中进行数值离散计算，同时在数值方法中引入限制器，进行数值间断侦测和限制，使得对于稀薄状态下超高声速的激波问题有了很好的解。该方法解决了传统NS方程无法解决非连续流动、分子动力模型在过渡流和滑移流动中计算量大的技术问题，对于临近空间各种飞行器的超高速流动及热力计算尤其是激波问题求解有着很大的优势。

附图说明

图1为本发明二维DG-NCCR离散方程计算方法流程框图；

图2为本发明实施例中二维圆柱绕流非结构网格图；

图3为本发明任意三角单元Ω与标准三角单元T的坐标转换示意图；

图4为本发明标准三角单元T与标准正方形单元R间的坐标转换示意图；

图5为本发明二维限制单元K及其相邻单元示意图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅是本发明一部分实施例，而不是全部的实施例。若未特别指明，实施例中所用的技术手段为本领域技术人员所熟知的常规手段。

本发明提出了一种适用于二维稀薄气体流动分析的数值计算方法，其核心思想是：将DG数值方法引入NCCR方程中，进行数值离散计算，且引入限制器进行数值间断侦测和限制。该方法所用流体输运方程由Boltzmann方程推导而来，而Boltzmann方程被普遍认为是全流域气体动力学控制方程。故本发明旨在突破传统NS方程无法解决非连续流动，以及分子动力模型在过渡流和滑移流动中计算量大的劣势，提出一套适合于临近空间飞行器飞行过程中的热气动分析方法。

本发明的数值计算方法对二维NCCR方程进行离散过程中采用了DG(间断伽辽金)数值方法，该方法结合了传统有限元(FEM)和有限体积(FVM)的优点，具体包括：(1)对非结构网格有很高的适应能力，适用于复杂构型；(2)加密或者加粗网格时不用考虑连续性，可轻松实现网格加密的自适应技术；(3)每个计算单元是独立的，其具有易编码以及并行的优点；(4)由于单元内采用间断近似的操作，因而DG方法很适合用于求解数值间断问题。DG-NCCR离散方程对于在临近空间的超高声速激波问题有着很好的求解优势，激波是超高声速绕流过程中的一种气流参数发生突跃变化的压缩波，其广泛存在炮弹，火箭，飞机等的超高声速飞行中，由于临近空间中空气稀薄，大气温度压力较低空环境更低，因而临近空间的超高声速问题更难求解分析。

为了更好地解释本发明的数值计算方法，先简要给出二维NCCR方程的来源。对于二维Boltzmann方程，通过引入指数形式的分布函数以及其对应的Rayleigh-Onsager耗散函数，进而得到了满足熵增原理的非线性耦合本构方程暨NCCR本构关系：

其中，

式中，Λ

q(k)＝sinhk/k，k

p为压力，ρ为密度，T为温度，Δ为体积附加正应力，E为单位质量的能量。

此处速度、散度均为二维情形。上述方程中的第一方程为质量、动量和能量的三大守恒方程，第二个方程为本构方程，其给出了粘性应力Π和热通量Q的表达式，用于封闭守恒方程组。

为了简化计算，对上述方程进行无量纲化：

其中，下标r表示参考状态，f

定义下列相似准则数：

Nδ表征粘性应力相对于静压的大小，其大小可以衡量系统远离平衡态的程度，Ma为马赫数，Re为雷诺数，Ec为马赫数Ma的函数，Pr为普朗特数。则通过上述无量纲参数和相似准则数可获得方程组的无量纲矢量形式：

其中，

方程组(2)的解即为本发明方案涉及二维稀薄气体流动的数值计算解。通过求解上述方程的温度、速度、压力以及密度等流动参数，即可完成对二维稀薄流动的温度场、速度场压力场以及密度场的求解，从而实现对稀薄超高声速飞行的热流动分析。

一种用于二维稀薄气体流动的数值计算方法，用于求解DG-NCCR离散方程，具体包括以下步骤：

S01、通过网格程序对二维求解流场进行网格划分，并生成网格文件。可以选用商业网格生成软件生成网格文件，网格生成软件可以为Gambit或者ICEM。该步骤中生成的网格一般为二维三角形网格。

S02、读取网格文件，并记录网格节点坐标参数。二维DG-NCCR程序读取网格文件，二维DG-NCCR程序中写入了本发明公开的用于二维稀薄气体流动的数值计算方法。

S03、根据节点坐标，获得标准正方形网格单元。具体地，先将步骤S01中的三角形网格单元转换成标准三角形单元，再转换成标准正方形网格单元。

S04、给定压力远场边界条件，通过给定的绕流马赫数计算速度，其余参数与周围的大气参数一致。将壁面设置为Langmuir滑移边界条件，此边界条件根据Langmuir吸附等温线，并考虑气体分子与壁面间的相互作用，从而确定壁面上的速度和温度。

S05、进行迭代计算：初始化计算条件，如最大迭代步数M、计算精度要求ε，根据CFL条件确定时间步长；根据边界条件给流场参数赋初值x

S06、将x

S07、采用不同的通量计算格式计算DG-NCCR离散方程中的通量。

S08、根据步骤S06、S07的计算结果，采用三阶TVD-RK对DG-NCCR离散方程进行求解，得到当前求解步骤的流场参数结果x

S09、遍历所有求解结果，查找“问题单元”并采用限制器对其单元内的求解参数进行限制，更新x

S10、计算迭代误差x

S11、循环计算结束。根据每个单元的x

由于NCCR的高度非线性和耦合性，其在一般的离散方法中难以得到很好的解，本发明的数值计算方法将DG方法引入到NCCR方程中进行数值离散计算，同时在数值方法中引入限制器，进行数值间断侦测和限制，使得本发明对于稀薄状态下超高声速的激波问题有着很好的解。

实施例

先确定DG-NCCR离散方程的表达形式。

通过近似解S

式中，

需要说明的是，各个单元中的基函数形式是相同的，，u

N是基函数的个数，其大小受计算阶数I的制约，具体关系为：

为了保证计算精度以及节省计算资源，阶数I的优选值为2，则N的优选值为6。用近似解S

/>

图1中的流程基本都是围绕上式(3)u

S01：采用商业网格生成软件，如Gambit或者ICEM读取待求解流场几何文件，随后在流动区域生成结构或非结构网格，由于二维非结构网格生成简单，且本发明方案中会将非结构网格转换成结构网格，因而该步骤生成的网格为一般的二维三角网格。如图2所示，为针对圆柱形飞行器飞行中的流动域所产生的网格，从图2中可看出，生成的二维网格均为普通的三角形。

S02：根据图2所示的网格图读取其中每个网格的几何信息，主要为每个网格单元的节点坐标。

S03：为了简化计算，基函数采用Dubiner基函数，而其计算需要在标准网格下进行，因而需要通过步骤S03将图2中每个网格单元转换成标准三角单元以及标准四边形单位，从而简化积分计算。具体来说，对于一个普通的非结构三角单元Ω，它的三个顶点1,2,3的坐标分别为(x

三角单元重心的横纵坐标分别为：

有了上述准备，先进行原始三角单元Ω到标准三角单元T的转换，如图3所示，T所在坐标系与笛卡尔坐标系之间的关系为：

Ω→T:

因整个单元内的积分是在标准正方形单元R＝{(a,b),-1≤a,b≤1}内进行的，详细地坐标转换如图4所示，T所在坐标系与R所在坐标系之间的转换关系为：

T→R:

对于方程(3)的边界积分

而对于剩下的单元面积分，则是在标准正方形单元内进行的，其对应雅可比行列式为：

S04：通过给定实际稀薄气体流动时的边界条件来约束方程组(3)将流场外边界设置为压力远场边界，通过给定的绕流马赫数计算速度，其余参数与周围的大气参数一致。将壁面设置为Langmuir滑移边界，此边界条件根据Langmuir吸附等温线，并考虑气体分子与壁面间的相互作用，从而确定壁面上的速度和温度。对于初始条件，根据边界条件进行计算，对于没办法计算出来的参数则将其赋值为0。

S05：通过前面给定的网格坐标转换公式以及雅可比行列式可得出面积分的表达式：

上式中的A

S06：对于边界积分，本发明方案采用数值通量函数来求解边界上的解，现存很多种通量格式，为了使计算中的密度和压力恒正，本发明采用一种恒正的通量格式来计算守恒方程中无粘项的通量：

这里，

上角标-和+分别表示单元交界面内部和外部的值，即计算单元Ω

S07：经过前述步骤中对离散方程(3)积分项的求解可得到其半离散形式：

/>

式中，

其中

同理，对s

这里，CFL为CFL数，其满足条件CFL≤1。步骤S07中的x

S08：由于DG方法在间断解(如激波)附近会产生非物理的数值震荡，使得计算出来的压力和密度出现负值，从而使结果发散，为了解决这一问题，本发明采用限制器技术捕捉“问题单元”(即间断附近的单元)，并对其解进行重构，从而抑制非物理震荡的出现。具体来说，本发明采用的限制器为TVB斜限制器，其核心是引入了rmin mod，从而避免了临界点处的精度弱化到一阶精度，rmin mod函数的具体表达式为：

其中，函数sign(x)的表达式为：

其中，ΔL是判断是否需要限制的参数，其与网格及其周围网格的几何参数和解有关，本发明中其表达式为：

其中,e分别为三角形单元边的总数，

而Δx

Δx

(x

对于二阶情况下(基函数个数为6)任意单元内求解值的平均值，其表达式为：

同理，可求出U

通过rmin mod函数以及单元K(i,j)及其周围单元内求解值的平均值表达式可建立限制的求解系数

/>

上式中，

通过上述式子，即可得到“问题单元”经过限制的近似局部解：

如此，便可对x

以上所述的实施例仅是对本发明的优选方式进行描述，并非对本发明的范围进行限定，在不脱离本发明设计精神的前提下，本领域普通技术人员对本发明的技术方案做出的各种变形、变型、修改、替换，均应落入本发明权利要求书确定的保护范围内。