英语原文共 13 页

几乎不可压缩流动的有限体积格子Boltzmann方法的气体动力学BGK格式

Weidong Li ， Wei Li

a中国北京市海淀区北京计算科学研究中心，邮编100094

b美国弗吉尼亚州诺福克市老多米尼加大学数学与统计系，邮编：23529

c清华大学深圳研究院，深圳518057

文章信息

文章历史：

收到日期：2017年6月17日

2017年12月12日修订

接受日期：2017年12月14日

2017年12月14日在线发布

关键词

有限体积点阵玻尔兹曼方法

气体动力学方案

bgk非结构网格

摘要

这项工作提出了有限体积格子Boltzmann方法（FVLBM）的气体动力学BGK方案。在本方法中，格子Boltzmann BGK方程的形式解析解用于确定单元界面上的通量。通过这种正式的分析解决方案，目前的FVLBM有两个显着特征。一个是格子Boltzmann方程（LBE）中碰撞项的影响包含在单元界面上的重建通量中，这减少了数值耗散，另一个是，不像FVLBM的传统通量方案，当前通量方案是空间和时间的函数，它为非定常流动一步提供时间和空间的二阶精度。此外由于碰撞不变，通过引入宏变量的预测步骤，平衡分布函数的隐含特征被消除，因此，本方案中的时间步骤仅受Courant-Friedrichs - Lewy（CFL）的限制状况。另外，讨论了边界条件的一些数值实现，并且为了验证本方案的准确性，还提供了几个二维不可压缩不稳定和稳定基准流的比较研究。清楚地验证了本方案的准确性和有效性。

介绍

格子Boltzmann方法已被建立为一个有效的不可压缩计算流体动力学（CFD）范例，在二维和三维空间中的均匀笛卡尔网格上工作。LBM 的均匀笛卡尔网格的限制极大地妨碍了其在工程应用中的应用。克服均匀笛卡尔网格约束的方法之一是用有限体积法求解格子Boltzmann 方程，该方法在复杂几何上具有很好的自适应性，并且在基于传统CFD的Navier-Stokes方程上取得了很大的成功。

有限体积格子Boltzmann方法（FVLBM）在过去二十年中见证了一些进展。Succi等人提出了第一个基于结构网格的有限体积LBM（FD-LBM）。Peng和Ubertini等人开发了基于非结构化网格的FV-LBM。在这些工作中，单元边界处的分布函数的通量直接通过它们在单元中心的值通过内插来计算，这被称为中心差分方案（CDS）。虽然CDS有低数值耗散，也缺乏鲁棒性。更糟糕的是，这些方案的时间步长不仅受到Courant-Friedrichs-Lewy（CFL）条件的限制，还受到放松时间的限制。为了提高稳定性， 提出了一些上风型有限体积格子Boltzmann方案[4-7]来获得单元边界分布函数的通量。这些方案提高了FVLBM的鲁棒性，并且格子Boltzmann方程中碰撞项的影响不包括在这些方案中。因此，这些方案具有与Navier-Stokes方程的动力学通量分裂方案[8]类似的行为，这是非常稳健的，但正如[9]中所指出的， 它们依赖于插值技术来减少其中的差异。尽可能多地插入左右状态以减少数值耗散。为了减少过度耗散，碰撞项应该包含在电池接口上的通量[9,10]，这是气体BGK方案[9]的关键成功，最近开发的离散统一气体动力学方案（DUGKS）[11]和统一气体动力学方案（UGKS）[12] 在连续统一极限。在目前的工作中，我们讨论了非结构网格上有限体积格子Boltzmann方法的气体动力学BGK方案的发展。提出的方案有两个优点。第一个是在格子Boltzmann BGK的形式解析解的基础上方程中，包括碰撞项的影响，这对于降低方案的大数值粘度是重要的; 第二是在一个时间步长中实现二阶时空精度，这对于非定常流动是有效的。
本文的其余部分安排如下。第2节详细讨论了非结构网格及气体动力学有限体积格子Boltzmann方法的公式，包括格子Boltzmann方程（LBE）的时空有限体积离散化，气体动力学BGK通量方案和数值实现的一些其他细节。第3节给出了目前气体动力学BGK的一些必要分析。然后，第4节提供了所提出的气体动力学有限体积格子Boltzmann方案 在以下流动情况下获得的数值结果：（a）稳定的Kovasznay流动;（b） 层流边界层在平板上流动;（c）不稳定的泰勒 - 格林涡流。最后，我们在第5节中总结了这篇论文。

气体动力学有限体积格式的数值公式

2.1基于非结构网格的时空有限体格子Boltzmann方法

图1.时空计算元素Omega;。

具有q离散ve-的离散速度Boltzmann方程位置{xi;alpha;|alpha; = 0, 1, . . ., (q minus; 1)}在d维空间x中可以给定为：

其中和alpha;是分布函数，并且分别对应于离散速度和q的碰撞项由具体的离散速度模型确定。碰撞项可以是BGK模型[13]或MRT模型[14,15]。

本文采用了BGK碰撞模型，即

式中和nu;是分别对应于离散速度、松弛时间、晶格声速和动力粘度的平衡分布函数。

通过介绍（D 1）维矢量中的梯度算子和其他欧氏空间-时间：

式中，可简写式（1）为：

在不失去一般性的情况下，除非另有规定，式（4）是暂时的一维和空间二维，即d=2。图1说明了体积△V的规则棱柱时空计算单元，其中，在时空元素的底部和顶部的三角形，分别有,区和,质心，同时分别是和平面。因为上面和下面的三角形是一致的，我们得出, △V = Si △t = Sirsquo;△t。图1所示的时空棱柱形单元V上等式（4）的积分得出：

将高斯-格林定理应用于式（5）可得出：

式中，a是时空体积的表面矢量△v。表面积分采用简单正交，碰撞项采用梯形法则，式（6）近似为：

其中

并且，在下面的内容中，除非指定了空间和时间坐标，并且是形心指数集，j是空间三角形的相邻空间三角形，即和共享一个公共边。此外，单元界面上的通量可以表示为

其中为时空元边表面中心的位置矢量，、、为单位法向量，中心的位置矢量和界面边缘i j的长度分别如图1所示，界面边缘i j为三角形i和三角形j共享的边界边。将式（8）代入式（7）给出

这是时间和空间的二阶精度。式（9）所定义的方案由于术语而隐含在碰撞项中，因此，当前FVBLM中的时间步仅由courant–friedrichs–lewy（cfl）条件为了使该方案更为明确，引入了宏变量的预测步骤。在预测步骤中定义

图2.tn超平面上的空间单元i。

把它代入式（9），我们得到

其右侧已知，因此确定。此外，由于碰撞不变量，宏观密度rho;和时间tn 1处的速度U可计算如下

然后计算出时刻的平衡分布函数。此外，，因此，可由式（10）中的计算得出。

2.2.单元界面上通量的气体动力学方案

为了减少过多的数值耗散，LBE中碰撞项的影响应包括在重构的通量对单元界面的影响[9,10]中，因此，用方程（1）中格点bgk方程的形式积分解来构造通量。图2所示超平面界面边缘ij上处格点bgk方程[9,12]的积分解为

式中 。

为了计算式（11）中的，首先应得到。利用式（13），我们需要重建形式积分解中出现的分布函数。对于二阶精度，平衡分布函数可确定为

其中，初始分布函数可表示为

替换方程。（14）和（15）在格点bgk方程的积分解中，即等式（13），我们得到

需要指出的是，公式（16）的推导没有连续假设，因此，公式（16）可用于模拟所有kn状态流，前提是流场中不存在不连续性，即无冲击[12]。在连续极限下，假设和chapman-enskog展开式给出是合理的。因此，等式（16）简化为

假设。可以看出，式（17）与时间的一阶泰勒展开式的结果具有相同的形式，这是完全一致的，因为如果，作为相空间中具有连续粒子速度的连续气体动力学bgk方案[9]，我们将替换为，或者如果单元界面上存在不连续性，就像在统一气体动力学中一样。方案[12]格点bgk方程的积分解将给出一个与一阶泰勒展开式完全不同的表达式。

从公式（17）我们得出

通过柯西-科瓦莱夫斯基程序[16]，时间导数可以交换为空间导数，即：

因此，式（18）变成

其中碰撞项的影响包含在通量中，通量是时间和空间的函数。此外，式（20）表明，目前有限体积达到了二阶时空精度，将式（20）代入式（9）得到

2.3.时刻界面边缘分布函数及其梯度

对于光滑流场，可以用线性插值法重建和时间的初始分布函数，以达到二阶精度，即：

其中和分别是从到和的距离和距离的比率。由于通常与图2中的不一致，的分布函数需要从的分布函数重建，也就是说，

图3.虚拟单元方法的说明。边界单元i和它的虚拟单元，它们的质心分别位于点和点。

由加权最小二乘法（5，6，17]得到的单元边缘，，在单元中心，处，从梯度和中计算出界面边缘iJ上的梯度。定义

然后我们可以计算如下[18]：

用右手边的第二项来避免奇偶解耦问题，用加权最小二乘法求出单元中心梯度和。

为了得到梯度，对目前的二阶方案合理地假设了。因此，式（20）变成

式中，为了保持守恒，平衡分布函数是从的宏观力矩得到的。

2.4.边界条件

在物理边界上，宏观Dirichlet和Neumann边界条件都是通过超平面上的重影单元法实现的，如图3所示。图中，空间单元是边界单元i对应边界边ab的虚拟单元，其质心位于。连接单元和i单元形心的线与边界AB相交于。

在边界边AB的宏观边界条件约束下，利用虚拟单元法可以得到上的宏观变量rho;和u。从而得到了平衡分布函数。通过线性插值，用等式（22）的类似方法计算处的平衡分布函数：

其中和分别是从到和的距离。处的非平衡分布函数可通过非平衡外推方案[6,7,19]计算，即：

式中，是处的非平衡分布函数。因此，我们在点有分布函数。

此外，的分布函数可以通过线性插值来计算：

一旦获得了的分布函数，就可以用加权最小二乘法计算单元i中的分布函数的梯度，如在内部单元中。

为了计算边界边界AB的通量，我们需要得到AB质心处的分布函数及其梯度，即。如式（23），我们有

其中，梯度可通过在图3所示的虚拟控制体上应用高斯格林定理来计算，即：

其中，分别为控制体的边界面向量，为控制体的体积，边界顶点a和b的分布函数，即和可通过反距离加权（idw）插值法计算[20]。此外，还合理地假设和使用来获得物理边界界面上的通量。

提出的气体动力学BGK方案的一些分析

1. U分量的速度等值线 （b）V分量的速度等值线

1. 压力轮廓 （d）流线

图4.Kovasznay流量在Re=40时的解析解。

1. 单元676 (b) 单元2706 (c) 单元10744 (d)单元43158

图5.Kovasznay流的网格

1. U分量的速度等值线 （b）V分量的速度等值线

1. 压力分布 （d）流线型

图6.用43158个近似均匀的三角形单元求解Re=40时的Kovasznay流

本节对本方案的性质进行了必要的讨论。

在连续体极限下，式（17）表明，目前的气体动力学方案成为连续体流的玻尔兹曼方程的一阶查普曼-恩斯科格展开解，即：

因此，本方案对不可压缩流动的数值模拟非常有效。更有趣的是，在高度非平衡、平滑的流动中，碰撞变得不重要，tau;gt;gt;△t和式（20）趋向于

这是玻尔兹曼方程在无碰撞状态下的解。此外，由于对流和碰撞的成分都包含在单元界面上，并且在单元中心保持隐含的碰撞项，因此本气体动力学BGK方案采用的时间步不受碰撞松弛时间tau;的控制，而是完全由CFL条件决定，即：

式中，sigma;为CFL数，i为对流光谱半径，定义如下：

因此，根据[10]中所提到的渐近保持（ap）方案的准则，本气体动力学方案可能具有渐近保持（ap）性质，而本气体动力学方案主要针对几乎不可压缩的流动设计。本方案AP的原因可能来自EQ起重作业人员。（17）至（20）。对大千牛流量模拟能力的论证将是我们今后的工作。

作为气体动力学方案之一，本方案与Dugks[11]具有许多共同的特点，另一种气体动力学方案，如在单元界面上使用时间的通量、通过等式（10）中给出的变换在单元中心使用隐式碰撞项等。因此，本方案与Dugks非常相似。然而，这

资料编号：[3723]