高雷诺数流动下2D-MRT格子玻尔兹曼方法的稳定性研究：正交化策略与松弛参数

《Computers & Fluids》：Stabilizing the 2D-MRT Lattice Boltzmann Method for High-Reynolds-Number Flows: Orthogonalization Strategies and Relaxation Parameters

　　摘要多松弛时间（MRT）格子玻尔兹曼模型通常仅适用于高粘度流体，这限制了其在相对粗糙网格上模拟高雷诺数流动的应用。本研究探讨了二维MRT格子玻尔兹曼模型在模拟低粘度高雷诺数流动时的数值稳定性和精度。本文采用了格点声速为1/3的D2Q9格点。特别关注了过渡矩阵的构建方法，该方法通常是通过对若干选定矩进行格拉姆-施密特正交化得到的。通过对二维MRT格子玻尔兹曼模型进行线性稳定性分析，研究了正交化策略和松弛参数的影响。比较了三种正交化策略：不进行正交化、传统无权重正交化以及加权正交化。研究结果表明，对变换矩阵进行正交化对于低粘度条件下的稳定模拟至关重要。线性稳定性分析表明，声学模式主要受与体积粘度相关的se以及幽灵松弛率s?的影响，而剪切模式则受决定运动粘度的sν以及幽灵松弛率sq的影响。加权正交化显著扩大了s?的可行范围，而sq与sν之间的依赖关系是固有的，无法消除。基于这些研究结果，确定了可行的松弛参数范围，为稳定的D2Q9 MRT模拟提供了实用指导。通过一系列基准问题进一步验证了研究结果，结果表明MRT能够以接近2的收敛速率再现主要流动特征。此外，加权正交化使得在粘度低至10?5、雷诺数高达104的情况下仍能实现稳定模拟。这些发现凸显了采用加权正交化过渡矩阵的MRT方法在实际应用中作为模拟低粘度高雷诺数流动的一种高效且稳健的解决方案的潜力。

引言格子玻尔兹曼方法是一种用于模拟流体流动的数值方法，它源自格点气体自动机。该方法的提出是为了克服格点气体自动机的若干局限性。与格点气体自动机使用布尔变量表示粒子占据状态不同，格子玻尔兹曼方法采用连续的粒子分布函数，从而显著降低了格点气体自动机模拟中固有的统计噪声。此外，格子玻尔兹曼方法引入了松弛时间参数，使得调节传输系数更具灵活性，而这在格点气体自动机中是受限的。因此，格子玻尔兹曼方法能够以更高的数值稳定性和精度模拟更广泛的流体行为。格子玻尔兹曼方法可视为玻尔兹曼方程的离散化形式。除了空间和时间上的离散化外，连续的速度空间也被替换为一组有限的离散速度。分布函数通过碰撞和流动两个主要过程演化。为模拟碰撞项，通常采用巴特纳加尔-格罗斯-克鲁克近似，该方法假设分布函数会在特征松弛时间内向其平衡状态放松，这一方法有效地模拟了流体分子的碰撞过程。通过查普曼-恩斯科格分析，格子玻尔兹曼方法在流体力学极限下可以恢复到宏观的纳维-斯托克斯方程。BGK模型中的松弛时间与流体的剪切粘度直接相关。尽管取得了巨大成功，但BGK模型在模拟高速或低粘度流动时往往会出现不稳定现象。通常，格子玻尔兹曼模拟中的粘度必须设定在0.01或0.001左右，进一步降低粘度则会导致数值不稳定，这严重限制了BGK模型在高雷诺数流动中的应用。由于提高流速也会导致数值不稳定，唯一的替代方案就是细化计算网格，但这会大幅增加计算成本并降低模型的效率。因此，人们提出了多种策略来提高格子玻尔兹曼方法的稳定性和稳健性。规则化模型引入了更高阶的平衡分布函数以及规则化的碰撞算子，可以有效过滤掉分布函数中的高阶非平衡成分。熵格子玻尔兹曼方法也采用了新的平衡分布函数，因为通过赫米特展开得到的截断平衡分布违反了离散H定理。熵格子玻尔兹曼方法没有采用固定的松弛过程，而是将碰撞步骤表述为一个约束最小化问题。BGK模型的另一项重要改进是多松弛时间格子玻尔兹曼方法。该模型最初是为D2Q9格点提出的，后来又扩展到了D3Q15和D3Q19等三维格点。MRT模型的碰撞操作在矩空间中进行，在这个空间中，分布函数的不同矩——如质量、动量、能量通量和粘性应力——直接对应于宏观物理量。MRT框架的关键创新在于它能够为不同的矩分配不同的松弛率，从而实现对动量扩散和热传递等各种物理过程的更精细控制。尽管在速度空间和矩空间之间转换增加了复杂性，但MRT模型并没有显著增加计算时间。通过适当的代码优化，其计算成本与BGK模型相当，仅需适度增加内存，计算时间大约增加10–20%。在数值稳定性方面，MRT模型的表现远远优于BGK方案。例如，在液滴振荡模拟中，MRT模型可以在比BGK模型所能承受的低五倍的粘度条件下稳定运行。在某些具有挑战性的场景中，其稳定性甚至提升了一个数量级。基于原始矩MRT框架，还发展出了中心矩和累积量碰撞模型，以进一步提升稳定性。然而，MRT模型也存在一些缺点。其中一个问题是它引入了大量自由参数，其中许多参数的物理影响尚不完全清楚。在某些情况下，MRT模型的稳定性甚至低于BGK模型，这说明其松弛参数需要针对每次模拟进行仔细调整。这种调整过程可能耗时较长，且往往需要手动试错。虽然已经有一些努力试图自动化参数优化，但这些方法通常是为特定问题或模型设计的，尚未找到通用解决方案。另一条研究路线着眼于过渡矩阵的作用，该矩阵将分布函数从速度空间映射到矩空间。构建过渡矩阵通常包括两个步骤：首先，选择一组对应于守恒量（如密度和动量）以及非守恒的高阶矩的特征向量；其次，从这些特征向量中构建一组线性独立且相互正交的基向量。正交化的目的是为了分离各个矩的松弛过程，但其实际效果仍是一个未解之谜。一些研究表明，高阶矩的松弛仍可能受到低阶矩的影响。另有研究者发现，非正交过渡矩阵也能达到类似的精度，并且具有更高的计算效率。第三种方法是使用加权内积对特征向量进行正交化，其应用在相关文献中有报道。Fakhari等人利用这种方法解决了多相流问题，但模型中选择的粘度与其他MRT模型中的粘度处于同一数量级。Yang等人证明了加权正交矩阵能够确保加权L2稳定性，但这一结论仅限于静水条件。因此，正交化过程的必要性和作用仍在讨论中，尤其是在当前模型失效的情况下。同时，还需要进一步的理论和数值研究来明确正交化在MRT格子玻尔兹曼方法中的确切作用。为研究正交化过程在MRT模型中的作用，本文采用了线性稳定性分析方法，该方法由Lallemand等人引入到MRT模型中。线性稳定性分析用于研究系统在稳态条件附近对微小扰动的响应，该方法有助于系统地研究各种模型参数——如松弛率、宏观量以及过渡矩阵——如何影响模型的色散特性和耗散特性。与主要适用于长波区域且侧重于模型与宏观流体方程一致性的查普曼-恩斯科格展开不同，线性稳定性分析能够提供模型在更宽波数范围内的谱特性，因此它是评估格子玻尔兹曼方法稳健性的重要工具。本文的其余部分结构如下：第2节回顾MRT模型的基本概念；第3节对不同的MRT模型进行线性稳定性分析；第4节给出数值验证结果；第5节总结主要结论。

片段MRT格子玻尔兹曼方法为便于说明，本文仅考虑D2Q9格点。D2Q9格点由九个离散速度ci及相应的权重wi组成，其定义如下：ci={(0,0)i=0,(±1,0),(0,±1)i=1,?,4,(±1,±1)i=5,?,8，wi={4/9i=0,1/9i=1,?,4,1/36i=5,?,8。MRT格子玻尔兹曼方程表示为：f(x+ciΔt,t+Δt)=f(x,t)?M?1SM[f(x,t)?feq(x,t)]+M?1(I?S2)MG，其中fi表示沿ci方向的分布函数，f=[f0,f1,?,f8]T。S=diag(s0,s1,?,s8)为松弛矩阵，M为...

MRT格子玻尔兹曼方法的稳定性分析现有的关于LMRT的分析主要集中在Γuo上，而与Γno和Γwo相关的研究则相对较少。本节通过研究MRT模型的谱特性，阐明了正交化策略和松弛率对数值稳定性的影响。该分析揭示了MRT方案中的不稳定性机制，并为选择合适的松弛参数以模拟低粘度流动提供了指导。

数值结果基于线性系统响应的假设，得到了MRT的谱特性。然而在实际情况中，真实流体对振荡的反应本质上是非线性的。此外，流体中的扰动并非由单一单色波表示，而是由不同频率和波数的波的叠加构成的。为考虑这些复杂性，本节给出了使用MRT模型模拟低粘度流动的数值示例。

结论在本研究中，系统地分析了MRT模型在模拟低粘度高雷诺数流动时的数值稳定性和精度。首先通过线性稳定性分析研究了正交化策略和松弛参数的影响。分析结果表明，对变换矩阵M的基向量进行正交化对于低粘度流动至关重要。声学模式和剪切模式在...都会出现分裂、交换和混合现象。

作者贡献声明李顺阳：撰写——初稿、研究工作。李婉：数据整理。桂楠：监督、概念设计。杨行端：数据整理。蒋胜耀：资金获取。李顺阳|李婉|桂楠|杨行端|蒋胜耀中国北京清华大学核能与新能源技术研究院、先进核能技术协同创新中心、反应堆工程与安全重点实验室，邮编100084。