英语原文共 10 页

用于光滑流动的半隐式气体动理学方法

王朋，郭照立. 计算机物理通信，卷205,2016年8月,页面22-31

摘要

在本文中，基于Bhatnagar-Gross-Krook（BGK）方程推导出了用于光滑流动的半隐式气体动理学方案（SIGKS）。作为一种有限体积格式，控制体积中平均流量变量的演化在欧拉框架下进行，而单元界面上的数值通量的构造则是从拉格朗日角度出发的。拉格朗日相位的采用使得碰撞和输运机制在通量评估中内在耦合在一起。因此，时间步长与颗粒碰撞时间无关，仅由Courant-Friedrichs-Lewy（CFL）条件决定。对单元界面上的重构分布函数的分析表明，SIGKS可以被视为具有附加项的改进型Lax-Wendroff格式。此外，来自重构中的隐含性产生的附加项能够增强该方案的数值稳定性。为了SIGK进行基准测试，我们对低马赫数和中等马赫数的光滑流动进行了大量的数值试验。结果表明，该方法具有二阶空间精度，与基准结果相比，可以给出精确的数值解。结果还表明，对于光滑流动，该格式的数值稳定性优于原GKS格式。

关键词：半隐式方案；数值耗散；气体动理学方案

1.简介

近年来，动理学方法作为新开发的计算流体动理学（CFD）技术引起了特别的关注。与基于直接离散化Navier-Stokes（NS）方程的传统CFD方法不同，动理学方法基于动理学理论或微粒动理学，它提供了流体动理学与潜在微观物理学之间的理论联系，从而为多尺度流动提供了有效的工具。迄今为止，已经提出了各种介观方法，例如格子气自动机（LGCA）、格子Boltzmann方程（LBE），气体动理学方案（GKS）和光滑粒子流体动理学（SPH），其中LBE和GKS是专门为CFD设计的。 LBE和GKS的动理学特性带来了许多独特的优势，使它们与经典的CFD方法区别开来。具体来说，对于LBE方法，相空间中的对流算子（或流动过程）是线性的，与Boltzmann方程相同，并与经典CFD方法中的非线性对流项形成对比;此外，LBE方法中的压力是使用状态方程计算的，而在经典CFD方法中，为得到压力，应求解具有速度应变的泊松方程，这经常产生需要特殊处理的数值困难，例如迭代或松弛。此外，LBE的动理学特性保证了LBE的内在并行性特征，这使得它比传统的CFD方法更直接地实现消息传递接口（MPI）和图形处理单元（GPU）以提高计算效率。与旨在解决几乎不可压缩流动的LBE方法不同，GKS主要设计用于可压缩流动，特别是涉及不连续和连续状态的冲击问题，这是经典CFD方法的数值挑战。在GKS中，通过基于某些动理学方程构造界面通量，可以实现逆流向和中心差异力学之间的光滑过渡，这确保了该方案同时捕获冲击不连续性和光滑流动的能力.特别是,对解决Navier-Stokes方案中,GKS已得到很好的发展，并成功应用于各种流体问题。

作为一种有限体积方案，GKS中NS解决方案的关键在于单元界面构建通量。通过不同的方法，基于动理学理论开发了几种动理学方案，例如基于无碰撞Boltzmann方程的动理学通量矢量分裂（KFVS）方案和基于Bhatnagar-Gross-Krook（BGK）方程的GKS，其中在数值通量的构造中考虑粒子碰撞。结果表明，GKS方法避免了在KFVS方案中为减少数值耗散而增加临时“碰撞”的带来的模糊性。在BGK型方案中，NS解决方案的气体动理学BGK方案已得到很好的发展，并已成功应用于从低速不可压缩到高超声速可压缩流动的连续流动模拟。

在本文中，我们提出了一种半隐式气体动理学方案（SIGKS），作为光滑流动的一种可替代的BGK型方案。所提出的方案的最显着特征在于，在单元界面处的通量的构造基于BGK方程的离散特征解，其来自拉格朗日相位。这种方法导致粒子碰撞和传输机制在一个时间步长内耦合在一起，这使得新方案表现出较低的数值耗散以及时间步长与粒子碰撞时间分离。此外，在方案中以隐式方式计算通量，与原始GKS方法相比，预期其具有改进的数值稳定性。进行一组数值模拟以验证当前方案作为NS方程可行的求解方案。

本文的其余部分组织如下。半隐式气体动理学方案在第2节中得出; SIGKS的数值通量误差分析和与原始GKS的比较在第3节中进行;第4节进行了数值试验，以验证新方案的性能，最后在第5节中得出了一些结论。

2.半隐式气体动理学方案

Boltzmann方程用单线态气体分布函数的演化方程表示多粒子动理学系统的行为。 其简化版之一是BGK模型：

其中是气体分布函数，是接近的平衡状态。 和都是空间，时间t，粒子速度xi;和内部变量eta;的函数。 粒子碰撞时间tau;与粘度和热传导系数有关。 均衡状态是麦克斯韦分布.

其中D是空间维度，K是内部自由度，rho;是密度，u是宏观速度，R是气体常数，T是气体温度。 分布函数f和保守变量W之间的关系是:

通量计算如下：

式中，=dxi;deta;是相空间中的体积元，deta;=hellip;,psi;由下式得出：

由于在粒子碰撞过程中质量、动量和能量是守恒的，所以f和g在时空的任何一点满足守恒约束：

为了开发有限体积方案，首先将计算域划分为一组控制体积。然后，将式（1）两边的psi;相乘，并在控制体Vi上的相空间和物理空间中对其积分，从tn到tn 1，由于粒子碰撞过程中保守变量的守恒，使得Vi中心保守变量的更新成为:

在这里

是跨越单元界面的宏观通量，h =Delta;t/ 2。中点规则用于对流项的时域积分。

根据方程式（6），更新平均守恒变量的关键因素是计算通量Fn 1/2，它可以由气体分布函数f（x，xi;，tn h）单独确定。 这里拉格朗日视角应用于f（x，xi;，tn h）的构造，我们将方程（1）在半个时间步长内沿特征线与单元界面处的端点（xb）进行积分，并用梯形规则来计算碰撞项：

通过它们在Xb的邻域泰勒展开式来近似f（Xbminus;xi;h，xi;，tn）和g（Xbminus;xi;h，xi;，tn），我们可以重写式（8）：

由于本方案的目标是求解耗散区域内NS的数值解，因此可以采用查普曼-恩斯科格展开法来近似分布函数。然后，可将两个近似值应用于等式（9）。首先，用fasymp;gminus;tau;（ ）的一阶Chapman-Enskog展开式来近似f，第二个近似式是asymp;，它与一阶Chapman-Enskog展开式一致，只包括流体力学变量rho;、u和T的一阶导数。到目前为止，分布函数在单元界面可以近似为：

其中g=，Ag=g。显然,a和A可以表示为：

式（10）中的隐式项g（, h）仅由守恒变量w（, h）确定，利用动量守恒，根据（10）可以明确地得到：

因此，尽管式（9）是隐式公式，但它可以精确地计算。此外，隐式可以提高格式的数值稳定性，这将在后面的数值部分显示出来。

将式（10）代入式（7），计算单元界面上的数值通量，并根据式（6）更新时刻的守恒变量。上述步骤可在下一个时间步长内重复。应该强调的是，分布的重构本质上是多维的，因此SIGKS可以自然地扩展到三维流动问题。此外，从上述偏差可以看出，与原来的GKS一样，所提出的气体动理学方案是更新宏观守恒流变量的有限体积法，微观分布函数仅用于计算宏观流量，在实际计算中没有出现。因此，本方案的内存和效率与原GKS和有限体积CFD方法相似。

最后，我们讨论了Prandtl数修正问题。由于在BGK模型中，所有的分子，不管速度如何，都有相同的粒子碰撞时间tau;，BGK方程给出了一个固定的Prandtl数，等于1。为了将Prandtl数修正到任何实际值，这里我们使用参考文献[4]中提出的热通量修正方法。

式中，是等式（4）中的能量通量，q是热通量

三、数值通量误差

我们现在分析所提出的SIGKS中的数值通量误差。为了简单起见，我们在下面的分析中假设=0。通过将动理学方程（1）沿其特征线积分，可以得到单元界面处的分布函数，其积分形式为：

其中= -（ -）。 然后，在一个时间步长中跨越单元i的界面的精确通量

因此,通过将一个在SIGKS方法中，由式（7）在一个时间步长内得到的数值通量与精确值比较。我们可以计算通量在计算中出现的误差。

为了清晰起见,我们首先改写由（10）重构的分布函数：

在这里，g=g 。另一方面,精确的分布函数可以通过调用一阶Chapman-Enskog展开和线性连续重构为：

基于式（17）、（18）, SIGKS中的数值通量在一个时间步长内的误差为：

对于连续流，这里,其误差为。

有趣的是, 与原始光滑流动气体动理学方案相比，本方案采用动理学方程的积分解来计算单元界面处的分布函数,例如：

对于光滑流动,上面的解可以由(18)近似。见参考文献[4,28]。此外,通过这种方案重构GKS是Navier-Stokes方程的显式Lax-Wendroff方案。 通过比较式(17)与(18),我们可以看到,主要的区别在于目前SIGKS和原始的GKS在于式(17)中的方括号中的后面三项成线性关系。由于项，这种差异使得目前SIGKS成为隐式方案，尽管它可以由首先获得守恒变量W（）精确地计算出来。额外的步骤将导致增加计算成本,但也可能提高数值稳定性。 对于Navier-Stokes方程，这种差异也将目前的SIGKS方案与经典Lax-Wendroff方案区分开来。

4、数值验证

在本节中,将模拟几个光滑的不可压缩和可压缩的例子以验证上节提出的方案,包括热 Couette 流、热 Poiseuille流、通过平板的层流、二维lid-driven腔流和激波结构。我们还通过模拟空腔流，就准确性、数值稳定性将目前的方案与原始的GKS相比较。

在所有的数值模拟中，碰撞时间tau;是由tau;=mu;/p决定的,在这里mu;是动力粘度，是压力; 时间步长是由CFL条件决定的,即:

eta;是CFL数，在这里设置其为0.5,除非另有说明,是最小网格间距,是最大的流速,是声速,gamma;是绝热系数;除非另有说明，否则气体被认为是单原子气体，gamma;=5/3; 中心差分法应用于近似守恒变量在空间上的梯度,从理论上来说，这将在二阶空间精度上取得很好的效果。

4.1 热Couette 流

热Couette流的分析解,用于计算SIGKS的精度。这种测试方法也在先前的研究中用来测试粘性和导热问题的数值方法。我们现在考虑的一个问题是：一种不可压缩的粘性流体在两个相距为H且无限大平板间流动，上板温度恒为，在水平方向上移动速度为U,下板温度恒为 (gt;)且静止。在模拟中,我们设置H=1,U=1,=1,=0,对于不可压缩极限,马赫数; 无滑动边界条件对上下板都适用,周期性边界条件用于流体两端。在粘度和热传导系数不变的假设下,分析温度和速度可以得到：

在这里y是较低板的高度,Pr是Prandtl数,Ec是Eckert数,定义为,为等压热容。

为了评估方案的精确度,测量了不同的网格分辨率温度的全局相对误差,全局相对误差定义为：

在这里T是通过目前的方案得到的温度,Te是分析结果。在仿真中,我们设置tau;=Delta;t，在垂直方向表示从1/10到1/40网格间距。如所示图1所示，对于Ec=20.,Pr=1.0和Pr=0.72,Ec=40，温度场相对误差的拟合曲线的斜率分别等于2.04和1.99,这证明新的方案在空间上具有二阶精度。在相同的初始条件下，网格间距为1/20,水平速度u分布沿垂直方向的模拟结果如图2所示,这与分析解吻合得很好。

图1 具有不同Ec和Pr的网格大小的数值误差

图2 水平速度的分布，沿垂直方向网格间距Delta;x=0.05

此外,我们还在的一个大范围内测试了当前的模型，代表了粘性耗散与热传导之间的比率。对于小值PrEc,温度在两板间沿垂直方向变化近乎为线性,这意味着与热传导相比，粘性散热的影响很弱。 随着的增加,粘性散热将成为最主要的,并且温度曲线偏离了线性分布。图3描述了当Ec=40时，在不同的Prandtl数下无因次温度分布,图4描述了当Pr=0.5，在不同的Eckert数下无因次温度分布,从中可以看出当前模型的数值结果与精确解高度一致，如图5所示,本方案可以准确地描述在2.0到100间的粘性和导热的流动。

图3 Ec=40时，在不同的Prandtl数下无因次温度分布

图4 Pr=0.5时，在不同的