不可压缩流动的高精度非标准格子玻尔兹曼方法

马 超，吴 杰

(南京航空航天大学 航空学院, 南京 210016)

0 引 言

在过去的三十年间，LBM 已经被发展到诸多领域，除了作为可压、不可压Navier-Stokes （N-S）方程的替代工具，还可以用作其他偏微分方程的解法器。尤其是在不可压流动模拟方面，LBM 被广泛认为是一种简单高效的数值工具。标准LBM 的碰撞迁移过程使得编程简单、易于并行、数值耗散低。然而，这样的求解方式也存在一些不足。例如，物理空间和离散速度空间的耦合导致标准LBM 必须使用均匀笛卡尔网格，对复杂外形的几何适应性差，时间步长和空间步长必须一致，即Courant-Friedrichs-Lewy（CFL）数为1，标准的单松弛BGK（Bhatnagar-Gross-Krook）格式在高雷诺数情况下数值稳定性差。此外，尽管LBM 的数值耗散低，然而它始终还是只有二阶空间精度。应用LBM 进行湍流模拟往往需要大量的网格，计算效率低。提高标准LBM 的精度是较为困难的，目前主要的策略是对速度空间高阶离散或者对平衡态函数高阶展开。这些不足给LBM 应用到高速、复杂外形的工程问题造成了一定困难。

格子玻尔兹曼方程（LBE）可以看作是离散速度玻尔兹曼方程（DVBE）的一种差分格式。已证明，物理对称性和格子对称性可以分开处理[1]，物理空间的离散可以独立于速度空间的离散，这样就可以将离散速度和计算网格进行解耦。在这一思想指导下建立起来的方法一般称为非标准格子玻尔兹曼方法（OLBM）。OLBM 主要分为两大类，一类是在欧拉框架下采用传统的数值离散方法，比如利用有限差分[2]、有限体积[3]、有限元法[4]直接求解DVBE（以下称为欧拉法）；另一类是在半拉格朗日（semi-Lagrangian）框架下借助插值方法实现粒子分布函数的演化（以下称为半拉格朗日法）[5-7]。早期的OLBM 大多只有二阶精度，并且相比标准LBM 会引入更大的数值耗散。近十多年来，高精度格式开始引起学者们的广泛关注。这类格式拥有较小的数值耗散，在多尺度流动的模拟上具有明显优势，例如湍流、气动噪声。和低精度的格式相比，高精度格式在取得相同误差条件下需要的网格要少得多，因而提高了计算效率。利用高精度格式求解N-S 方程数值模拟不可压流的算法一般较为复杂，主要因为N-S 方程的对流项是非线性，黏性项包含二阶偏导数，以及压强和速度的耦合。最常用的投影法或者压强修正法都需要求解压强泊松方程，该方程不容易收敛。非线性的对流项对于通过插值基函数在网格单元内构造数值解高阶多项式的高精度方法来说，还可能会产生混淆误差。而DVBE 可以看成一组线性偏微分方程，它的曲线坐标形式相比N-S 方程也要简洁得多，所以用高精度格式求解DVBE 要容易得多。事实上，相关工作早在21 世纪初就有展开。Shi 等[8]采用间断伽辽金谱元法求解DVBE；Düster 等[9]也发展了一种间断伽辽金格子玻尔兹曼法。然而，由于这两种方法缺少网格收敛性研究，它们的收敛阶能否取得高阶并不清楚。后来，Hejranfar 和Ezzatneshan[10-11]提出了一种紧致有限差分格子玻尔兹曼法（CFDLBM）并做了收敛性分析，验证了该方法能够取得四阶精度。由于中心型的紧致有限差分格式没有数值耗散，为了确保计算的稳定，必须引入过滤器。他们将CFDLBM 和标准LBM做了比较，当取得相似的误差时，CFDLBM 的效率更高。为了避免过滤操作、简化计算过程，后续又有相关工作采用带有数值耗散的高阶差分格式。如利用五阶WENO 格式[12]以及五阶迎风型紧致有限差分格式[13]求解DVBE。有限差分对于复杂外形的几何适应性较差，为了改进这一问题，Li[14]将高精度谱差分方法（SD）和DVBE 结合，提出了SDLBM。通过坐标变换形函数，SDLBM 可以应用非结构网格。Li 对SDLBM 做了详细的收敛性分析，证明了在被压缩性误差影响之前，SDLBM 可以获得高阶精度，并且算法简单、易于并行。

以上这些方法的时间离散都采用显式格式。当雷诺数较高时，DVBE 碰撞项的刚性将严重限制时间步长，大大降低计算效率。目前应对这种问题的方法主要有两种，一种是将碰撞项隐式处理，直接求解DVBE（以下称为直接法），常用的方法有Guo 和Zhao[15]提出的基于半隐格式的显式方法和隐式-显式Runge-Kutta 方法（IMEX）等[16]；另一种是将计算分为两步（以下称为分步法），首先执行碰撞步，这和标准LBM 一样，然后在迁移步采用高精度格式求解一个纯对流方程。通过这种方法学者们也提出了不同的高精度OLBM，如Min 和Lee[17]提出的谱元间断伽辽金格子玻尔兹曼法以及Zadehgol 等[18]提出的节点间断伽辽金格子玻尔兹曼法。然而，鲜有文献比较分步法和碰撞项隐式处理的直接法之间的数值表现。

另一种构造高精度OLBM 的思路是在半拉格朗日框架下采用高阶插值[19]。由于不需要求解对流方程，没有CFL 条件的限制，这种方法能在大时间步长下保持稳定，并且每个时间步空间算子（插值）只演化一次，因此这种方法计算效率较高。然而事实上，半拉格朗日法误差和时间步长是成正比的。要想获得高精度，还是需要取较小的时间步长。此外，半拉格朗日法的离散误差相比分步法要更大[20]。因此，这几种高精度OLBM 的优劣依然值得进一步讨论。

本文将最近较为流行的通量重构（FR）格式引入OLBM，发展了一种高阶通量重构格子玻尔兹曼方法（FRLBM）。FR 格式为多种单元内多自由度的高阶格式，如节点间断有限元（NDG）和SD，提供了一个统一的框架[21]，且算法更加简单，计算效率更高[22]。在FRLBM 的框架下，本文比较了直接法和分步法的精度以及稳定性，并且对比了基于半隐格式的显式方法和一阶IMEX、二阶IMEX 以及三阶IMEX 对不同雷诺数下CFL 数取值的影响。研究结果旨在进一步加深对高精度OLBM 的认识，为发展更加高效高精度的OLBM 提供参考。

1 FRLBM 计算方法

1.1 控制方程

直接法求解的方程为式（1），而分步法则先采用θ方法将式（1）在时间区间[t,t+ Δt]沿特征线方向积分，得到[17]：

本文采用高精度FR 格式对式（1）和式（11）进行空间离散。

1.2 空间离散

本文仅讨论四边形单元下FR 格式，考虑以下的二维守恒律：

式中，K表示定义该网格单元的节点数量，M为坐标转换函数。根据坐标变换关系，守恒方程（12）在标准单元下的表达形式为：

为了在单元内构造高阶的连续通量插值多项式，需要定义两类点集，即求解点和通量点。求解点本文应用Gauss-Legendre 点，通量点为每个方向上的边界点。三阶FR 方法的标准单元如图1 所示。

图1 三阶FR 方法标准单元求解点（圆形）和通量点（正方形）分布Fig. 1 Arrangement of the solution points (circles) and the interface flux points (squares) in a standard element for the third-order FR method

1.3 时间离散

如果采用显式时间离散格式求解方程（1），每个单元i的时间步长受限于CFL 条件以及碰撞松弛时间，如下式所示：

1.4 边界条件

本文数值算例中涉及的无滑移壁面及远场边界处理，采用非平衡态外推的方法[27]，将边界上的分布函数写成平衡态和非平衡态之和的形式：

2 数值算例

2.1 Taylor-Green 涡

Taylor-Green 涡问题被广泛用来研究非定常不可压流动算法的收敛性，这个问题的解析解为：

2.2 顶盖驱动方腔流

表1 直接法和分步法不同网格下的误差及收敛精度Table 1 Error and convergence order of direct method and two-step method for different grid resolution

图2 直接法和分步法不同CFL 数下的误差和计算时间Fig. 2 Errors and computational time at different CFL numbers by direct method and two-step method

为了更好地捕捉方腔角落处的涡，采用非均匀网格：

其中，a=tanhβ ， β=1.5。当Re= 3 200，采用的网格单元数为30×30；当Re= 5 000 和7 500，采用的网格单元数 为40×40；当Re= 10 000，采 用 的 网 格 单 元 数 为50×50。40×40 非均匀网格单元的分布如图3 所示。图4 为p4FRLBM 计算得到的四个不同雷诺数下的流线。从图中可以看出，Re= 3 200 时右下角出现小涡，Re= 5 000 时左下角出现小涡。随着雷诺数的增加，底角三级涡越来越明显。由于FR 格式的高阶精度和高分辨率特性以及非均匀网格的使用，这些流动细节被FRLBM 清楚地捕获。图5 定量比较了方腔中线的无量纲速度，其中左侧纵坐标和右侧纵坐标分别为x方向和y方向无量纲速度。本文结果和相关文献[10,12,29]吻合良好，尽管本文采用了较少的网格单元。此外，本文还比较了4 种时间离散方式在4 个Re下收敛所能取到的最大CFL 数，如图6 所示。可以发现，基于半隐格式的显式方法和一阶IMEX 所能取到的最大CFL 数相同，为0.25。并且，随着雷诺数的增大，CFL 数降低。而二阶IMEX 和三阶IMEX 即使在雷诺数较大时，CFL 数依然可以取为1 附近，且三阶IMEX 最大CFL 数略大于二阶IMEX。考虑到三阶IMEX 要增加一步空间算子，二阶IMEX 是较好的时间离散方式。

图3 顶盖驱动方腔流40×40 非均匀网格单元的分布Fig. 3 Non-uniform mesh distribution with 40×40 for lid-driven cavity flow

图4 顶盖驱动方腔流在Re = 3 200、5 000、7 500 和10 000 时的流线图Fig. 4 Streamlines for lid-driven cavity flow at Re = 3 200, 5 000, 7 500 and 10 000

图5 顶盖驱动方腔流在Re = 3 200、5 000、7 500 和10 000 时中线上的速度型比较Fig. 5 Comparison of centerline velocity profiles for lid-driven cavity flow at Re = 3 200, 5 000, 7 500 and 10 000

图6 顶盖驱动方腔流在Re = 3 200、5 000、7 500 和10 000 时四种时间离散格式的最大CFL 数比较Fig. 6 Comparison of maximum CFL numbers of the four time discretization methods for lid-driven cavity flow at Re = 3 200,5 000, 7 500 and 10 000

2.3 圆柱绕流

图7 圆柱绕流计算域及贴体网格单元的分布Fig. 7 Computational domain and body-fitted mesh distribution for flow over a circular cylinder

圆柱绕流是验证算法模拟复杂外形绕流的经典算例之一，本文计算了瞬态雷诺数Re= 200 的情况。采用沿径向延展的非均匀O 型贴体网格，网格单元为1 次曲线单元，即式（13）中的K= 8。径向延展的变化由下式确定：为CL的一半。此外，将本文计算得到的CD、CL，以及斯特劳哈尔数St和相关文献[7,11,30,31]做了定量比较，如表2 所示，可以发现吻合良好。斯特劳哈尔数的定义为St=D/u0Tp， 其中D为圆柱直径，u0为来流速度，Tp为CL随时间变化的峰值周期。因此，本文的方法模拟复杂外形绕流是可靠的。

图8 圆柱绕流Re = 200 时的瞬时流线及涡量等值线Fig. 8 The instantaneous streamlines and vorticity contours for flow over a circular cylinder at Re = 200

图9 圆柱绕流Re = 200 时升阻力随时间的周期性变化Fig. 9 Periodic variations of drag and lift coefficients for flow over a circular cylinder at Re = 200

表2 圆柱绕流Re = 200 时结果和文献比较Table 2 Comparison of the results of flow over a circular cylinder at Re = 200 with the available results.

3 结 论

OLBM 可以改善标准LBM 的几何适应性，本文回顾了不可压流动高精度OLBM 的发展过程。相比基于N-S 方程的高精度不可压流算法，基于DVBE 的高精度算法要容易实现得多。本文基于高精度通量重构格式发展了一种通量重构格子玻尔兹曼方法，并在此框架下比较了直接求解DVBE 和先碰撞、再求解纯对流方程这两种方法的数值表现。研究发现，在小时间步长下，这两种方法误差接近，采用高精度格式都能获得高阶的收敛阶。然而，当时间步长增大时，分步法的误差出现明显上升，而直接法误差几乎不变。数值测试显示，在取得相似误差时直接法能取更大的时间步长，运算周期较短、计算效率更高。此外，采用二阶隐式-显式Runge-kutta 推进的直接法，稳定性优于分步法。接着，通过顶盖驱动方腔流算例，比较了4 种碰撞项隐式处理的时间推进方法在不同雷诺数下所能取到的最大CFL 数。数值结果表明，二阶IMEX 是较好的时间离散方法。然而，尽管IMEX 能使时间步长的选取摆脱松弛时间的限制，但由于对流项依然显式处理，在模拟壁湍流问题时较小的边界层网格厚度同样会严重限制时间步长。因此，发展低存储、全隐式的时间推进方法是OLBM 的重要发展方向之一，目前已有学者开展了一些工作[32-33]。本文采用的碰撞模型是单松弛，为了增强稳定性，也可以采用更先进的碰撞模型，比如多松弛等。此外，高精度OLBM 还可以直接拓展到湍流、多相流、热流等应用中[34-35]，目前相关工作正在开展中。