递归正则化格子玻尔兹曼方法：强制实现、边界处理及其在大涡模拟中的应用

《Computers & Fluids》：Recursive regularized lattice Boltzmann method: Forcing implementation, boundary treatment and application in large eddy simulation

　　李玉倩|齐一鸣|沈杰|陈松泽|牛中华|彭成|王连平
中国南方科技大学，广东省湍流研究与应用重点实验室，复杂流动与软物质研究中心，力学与航空航天工程系，深圳，518055

摘要
众所周知，格子玻尔兹曼方法（LBM）的递归正则化（RR）方案在模拟高雷诺数流动时表现出更好的数值稳定性。然而，该方案中的外体力处理和边界处理仅被简要讨论过。此外，离散粒子速度集及其展开阶数对方案整体性能（计算效率、数值稳定性和准确性）的影响尚未完全理解。在本研究中，推导了带有强迫项的连续玻尔兹曼方程的非平衡分布的递归关系，并讨论了不同离散速度模型和Hermite展开阶数的截断误差。然后，提供了玻尔兹曼方程的离散化方法和实现细节。此外，还提出了一种正则化边界处理方法。通过几个基准案例来说明RR-LBM的优势。首先，使用Taylor–Green涡流案例表明，由于去除了O(Ma3)阶误差，其准确性得到了提高。其次，使用盖驱动腔流案例展示了我们提出的边界处理的改进数值稳定性以及在方案中包含三阶非平衡分布的重要性。第三，通过层流和湍流通道流案例验证了强迫项的实现。研究发现，包含27个粒子速度的三阶Hermite展开（即27V-33）可以提供与原始RR-LBM相当的结果。粗网格模拟表明，27V-33方案比19速度模型（19V-33）更稳定。最后，将RR方案扩展到大涡模拟中，结果表明所提出的27V-33模型具有合理的准确性。

引言
在过去三十年中，格子玻尔兹曼方法（LBM）已成为Navier–Stokes（N–S）方程的另一种求解方法，并因其简单的算法结构、固有的并行效率以及在处理复杂边界条件方面的灵活性而被广泛用于模拟复杂流动[1]、[2]、[3]。LBM包括两个主要步骤：碰撞步骤，描述粒子向其平衡状态的松弛；以及流动步骤，将分布函数在空间和时间上传播。在宏观矩空间中，LBM实际上表示了一组矩的演化方程，碰撞步骤可以被视为不同矩的松弛过程。标准LBM使用Bhatnagar–Gross–Krook（BGK）模型或单松弛时间（SRT）模型[4]、[5]、[6]，其中所有非平衡矩在碰撞项中以单一松弛时间松弛到其平衡对应值。尽管BGK模型实现简单且计算效率高，但它存在数值不稳定性问题，尤其是在高雷诺数流动中。普遍认为LBM的不稳定性源于碰撞核，这是由于高阶幽灵矩的非物理松弛[7]。数值不稳定性的另一个来源是边界条件的不一致处理[8]。

为了提高LBM的数值稳定性，提出了几种碰撞模型。多松弛时间（MRT）模型已被证明可以提高LBM的数值稳定性[7]、[9]、[10]。在MRT模型中，不同矩以不同的速率松弛，从而更好地控制各种动能模式的数值耗散并提高数值稳定性。然而，MRT模型由于分布函数和矩空间之间的转换而引入了额外的复杂性。同时，寻找最优松弛参数需要额外的努力，例如使用线性稳定性理论。为了进一步提高MRT方法的数值稳定性，提出了几种方法，其中碰撞在与移动框架相关的矩空间中进行，例如Dubois等人的中心矩LBM[11]和Geier等人的累积LBM[12]。这两种方法都显示出比MRT方法更好的稳定性。此外，通过适当选择松弛时间，累积LBM可以实现四阶精度[13]、[14]。虽然累积LBM通过提高伽利略不变性和减少矩之间的非物理耦合来增强稳定性，但它需要在碰撞步骤中对分布函数、中心矩和累积量进行额外的转换[12]，从而增加了算法复杂性和计算成本。

另一种提高LBM数值稳定性的方法是正则化方法，该方法通过正则化步骤约束高阶矩。现代正则化LBM由Latt和Chopard提出[15]，其中仅保留与粘性剪切应力相关的二阶矩，所有三阶及更高阶矩被排除。值得注意的是，Ladd之前也采用了类似的策略[16]、[17]、[18]，他认为在弱可压缩流动的碰撞算子中应仅保留二阶张量（粘性剪切应力），而不会影响N–S方程的高阶非流体动力学模式可以忽略。然而，在他的工作中，应力张量的对称部分和不对称部分是分别松弛的。这种方案实际上是一种多松弛时间（MRT）方法。Zhang等人[19]使用Shan等人提出的基于Hermite展开的框架将正则化方法扩展到三阶[20]。他们认为，与非平衡分布类似，非平衡分布也应该在指定的Hermite展开阶数（例如三阶）处截断。这种正则化方法称为投影正则化，因为非平衡系数是直接从分布函数计算得出的，将分布函数投影到Hermite空间。正则化步骤以非常低的计算成本增强了稳定性，因为只需要简单的矩计算和相应的投影。

2015年，Malaspinas[21]提出了一种递归正则化（RR）碰撞模型。这种方法引入了非平衡系数的递归关系，即可以从低阶非平衡系数递归获得三阶及更高阶的非平衡系数。递归关系是使用Chapman–Enskog分析和低Knudsen数近似推导出来的。与BGK或MRT方法相比，RR模型表现出更高的稳定性和准确性。然而，递归关系是在没有考虑外体力情况下推导出来的，尽管平衡和非平衡分布函数都在Hermite空间中展开到了六阶，但这种截断的合理性以及展开阶数的讨论并未明确给出。后来，Jacob等人[22]使用D3Q19（三维，19个粒子速度）模型开发了一种混合递归正则化（HRR）模型，其中平衡和非平衡分布函数都在Hermite空间中展开到了三阶。在他们的工作中，三阶Hermite系数被正交化以确保恢复平衡分布的三阶矩。与原始RR模型[21]不同，后者是从非平衡分布的二阶矩重构应力张量的。在HRR中，偏差应力张量是通过非平衡分布的二阶矩和有限差分方法组合得到的。引入了一个权重因子来控制数值耗散，从而提高了数值稳定性。随后，Feng等人考虑了HRR中的外部强迫，并成功将其应用于气象流动[23]、[24]、[25]，但没有提供带有强迫项的递归关系的详细推导，也没有充分讨论实现细节。

在数值方案中，边界条件对于获得准确和稳定的结果非常重要。Latt和Chopard[26]比较了正则化LBM的不同边界处理方法，发现仅在低雷诺数下重建未知分布的方法更准确，而在高雷诺数下替换所有分布的方法更稳定。Verschaeve和Müller[27]将边界条件从平面壁扩展到曲面无滑移壁，用于正则化LBM，其中使用了有限差分来重建非平衡分布。然而，这些方法仅限于二维情况，需要进一步验证三维问题。关于递归正则化方法，Fend等人[28]为可压缩流动的HRR中的曲面边界提出了一种基于有限差的方案，这需要空间插值，从而增加了计算成本。通常，为了效率，需要一种局部简单且稳定的边界处理方法。

模拟高雷诺数流动具有挑战性，因为需要高分辨率来捕捉所有湍流尺度。对于壁面边界流动，由于固体壁附近的精细结构以及获得统计收敛结果所需的较长平均时间，计算成本更高。大涡模拟（LES）为非稳态湍流提供了一种有用且相对高效的模拟方法。在LES中，对Navier–Stokes方程应用空间滤波器，并引入涡粘性来表示从已解析尺度到未解析尺度的有效能量传递。未解析尺度的贡献由具有特定涡粘性模型的亚格子尺度（SGS）应力项表示。Hou等人[29]首次尝试将LES结合到LBM中，他们在LBM中采用了Smagorinsky涡粘性模型来模拟二维盖驱动腔流。通常，Smagorinsky模型的耗散相对较大，需要壁函数或动态模型来校正边界节点处的速度和应力，这需要额外的努力和计算时间[30]、[31]。作为替代方案，LBM中考虑了不需要额外壁处理的壁适应型亚格子模型，包括壁适应型局部涡粘性（WALE）模型[32]和Vreman模型[33]。由于稳定性问题，LBM中的LES通常使用MRT模型。由于递归正则化格子玻尔兹曼方法（RR-LBM）可以提高LBM的稳定性，Wang等人[34]将RR-LBM与LES集成，并比较了MRT-LBM和RR-LBM的结果。他们发现，尽管MRT-LBM的耗散性较低，但它会在湍流中产生虚假的非物理噪声，表明RR-LBM具有更好的数值稳定性。然而，他们的工作中使用了Smagorinsky模型和扩散界面浸没边界方法。使用其他涡粘性模型和边界处理方法检查RR-LBM的大涡模拟数值结果仍需进一步研究。

如上所述，RR-LBM已被证明可以提高数值稳定性。尽管如此，仍有几个关键问题需要澄清，包括带有强迫项的非平衡系数递归关系的严格推导、平衡和非平衡分布的Hermite展开的截断阶数，以及边界条件的实现。在本文中，我们将首先推导带有外体力的一般情况的递归关系，并讨论平衡和非平衡分布的展开阶数对RR-LBM整体性能的影响。然后，提供玻尔兹曼方程的具体离散化和实现细节，包括提出的局部边界处理方法。使用几个基准案例验证了该方案，包括Taylor–Green涡流、盖驱动腔流、层流通道流和湍流通道流。最后，检查了使用RR-LBM的LES结果。

本文的其余部分组织如下：第2节介绍方法论，重点介绍带有强迫项的递归关系的推导、粒子速度和Hermite展开截断阶数对恢复的宏观方程的影响以及边界处理。第3节报告了使用四个基准案例的数值结果，并提出了降阶RR-LBM，随后将LES结合到RR-LBM中以模拟湍流通道流。第4节提出了主要结论。

**节选**
带有强迫项的玻尔兹曼-BGK方程可以写为：
?f/?t + ξ_j?f/?x_j + bj?f/?ξ_j = ?1/τ[f ? feq]
其中，右侧采用Bhatnagar-Gross-Krook（BGK）近似作为碰撞核。分布函数 f(x,ξ,t) 在空间位置 x 和时间 t 处定义，表示具有速度 ξ 的粒子的概率密度。bj 是单位质量上的外力的 j 分量，其中 j=1,2,3 表示三维空间中的三个方向。feq 是局部平衡分布函数。

在本节中，我们首先使用周期性泰勒-格林涡流模型来检验不同速度模型以及平衡和非平衡分布的不同展开阶数的准确性，该模型忽略了固体壁边界处理的影响。接下来，在盖板驱动的腔体流中验证了所提出的边界处理方法，并通过层流通道流评估了外体力作用的实现情况。然后，我们提出了一个递归正则化方案。

为了提高格子玻尔兹曼方法在模拟三维高雷诺数壁面湍流时的数值稳定性，我们将最近提出的递归正则化（RR）方案扩展到包含强迫项的一般情况，并提出了一种递归正则化的边界处理方法。此外，还讨论了使用不同离散粒子速度和赫米特展开阶数时的数值准确性和稳定性。

作者贡献声明：
李玉倩：撰写 – 审稿与编辑、撰写 – 原稿、可视化、验证、方法论、研究、形式分析、数据管理。
齐一鸣：撰写 – 审稿与编辑。
沈杰：撰写 – 审稿与编辑、项目管理、方法论、研究、形式分析、概念化。
陈松泽：撰写 – 审稿与编辑、概念化。
王连平：撰写 – 审稿与编辑、监督、资源管理、项目管理、研究、资金支持。

利益冲突声明：
作者声明他们没有已知的可能会影响本文研究工作的财务利益或个人关系。

致谢：
本研究得到了国家自然科学基金（NSFC 编号 T2250710183、U2241269、91852205）、NSFC “非线性力学中的多尺度问题”优秀研究群体计划（编号 12588201）、广东省湍流研究与应用重点实验室（2023B1212060001）以及粤港澳大湾区数据驱动流体力学与工程应用联合实验室（2020B1212030001）的支持。计算资源由...提供。