【LBM教程 9】LBM 稳定性和精度

关注同名微信公众号，获取更多信息

LBM通过离散化的分子运动来模拟流体的宏观行为，其理论基础扎根于统计物理学和分子动力学，而其数值计算的优势则在于天然适合并行计算和处理复杂边界问题。

然而，随着LBM在工业和科研领域中应用的不断深入，对其稳定性和数值精度的研究也越来越受到关注。

稳定性分析

稳定性在数值方法中至关重要，尤其对于LBM这种显式时间推进的方法更是如此。

在LBM中，稳定性主要受两个方面影响：碰撞过程中的松弛时间参数以及迁移过程中的离散效应。为了更直观地理解稳定性，常常采用Von Neumann稳定性分析方法，通过线性化分析，考察误差随时间演化的情况。

1、松弛时间与稳定性

在BGK碰撞模型中，松弛时间 $\tau$ 是调控系统向平衡状态恢复的主要参数。通常可以通过Chapman-Enskog展开将松弛时间与宏观流体的动力粘性系数联系起来，其关系为：

$$\nu =c_s^2\left(\tau -\frac{\Delta t}{2}\right)$$

为了确保流体粘性系数为正，必然要求：

$$\tau >\frac{\Delta t}{2}$$

这一不等式既是物理意义上的要求，也是数值稳定性的基本条件之一。

如果 $\tau$ 过小，则系统的数值误差和振荡将急剧增大，导致不稳定；而 $\tau$ 过大则会带来过度粘性化，使得精度下降。

2、Von Neumann稳定性分析

为了进行更细致的稳定性研究，我们将对LBM离散更新公式进行线性化处理。设扰动量为 $ \delta f_i(\mathbf{x}, t)$，令

$$f_i(\mathbf{x},t)=f_i^{(0)}+\delta f_i(\mathbf{x},t)$$

其中 $f_i^{(0)}$ 为平衡态分布函数。经过线性化处理后，更新公式可以写成：

$$\delta f_i(\mathbf{x}+{\mathbf{c}}_i\Delta t,t+\Delta t)=\left(1-\frac{1}{\tau }\right)\delta f_i(\mathbf{x},t)$$

对该关系进行傅里叶变换，假设扰动可以表示为平面波形式：

$$\delta f_i(\mathbf{x},t)={\hat{f}}_i\exp [i(\mathbf{k}\cdot \mathbf{x}-\omega t)]$$

代入上式后，我们可以得到放大因子 $G$：

$$G=\exp (-i\omega \Delta t)=\left(1-\frac{1}{\tau }\right)\exp (-i\mathbf{k}\cdot {\mathbf{c}}_i\Delta t)$$

对于稳定性来说，需要保证所有模式的幅值不超过1，即：

$$\left|G\right|\le 1$$

由此可以得到稳定性条件，这一条件在实际数值计算中常常需要结合具体的离散速度模型和边界条件综合考虑。需要指出的是，上述推导为理想化条件下的分析，实际应用中还需考虑非线性效应和边界处理带来的附加影响。

3、数值扩散与离散效应

在LBM中，数值稳定性不仅依赖于碰撞模型和松弛参数，还受到空间离散和时间离散引入的数值扩散（或耗散）效应的影响。

离散化引入的数值耗散在某种程度上可以缓解震荡问题，但过多的数值耗散则会使解的精度降低。通过对更新公式的高阶展开，可以分析出离散化误差主要是以 $\Delta t^2$ 和 $\Delta x^2$ 的形式出现，这要求在实际计算中合理选择网格大小和时间步长，确保这些误差项在数值解中处于可接受范围内。

在实际应用中，通常会通过调节参数或者采用高阶LBM模型来降低离散化误差，同时也可以利用多松弛时间模型（MRT）来改善稳定性和准确性。

精度分析及Chapman-Enskog展开

为了验证 LBM 能够正确恢复流体的宏观行为，我们需要通过 Chapman-Enskog 展开方法，从微观的离散演化方程出发，推导出宏观的 Navier-Stokes 方程。该过程不仅验证了 LBM 模型的正确性，也为数值精度的分析提供了理论基础。

1、Chapman-Enskog展开基本思想
Chapman-Enskog展开假设分布函数 $f_i$ 可以展开为多个阶次：

$$f_i=f_i^{(0)}+ϵf_i^{(1)}+{ϵ}^2{f}_i^{(2)}+\cdots$$

同时，对时间和空间的导数也进行多尺度展开：

$${\partial }_t=ϵ{\partial }_{t_1}+{ϵ}^2{\partial }_{t_2}+\cdots ,\nabla =ϵ{\nabla }_1$$

其中，$ϵ$ 为小参数，代表Knudsen数（无量纲自由程长度与特征长度的比值）。将离散更新公式进行Taylor展开，并将各阶项分离，可以得到一系列方程：

• 零阶：
  $$f_i^{(0)}=f_i^{eq}$$

• 一阶：
  $${\partial }_{t_1}{f}_i^{(0)}+{\mathbf{c}}_i\cdot {\nabla }_1{f}_i^{(0)}=-\frac{1}{\tau }{f}_i^{(1)}$$

• 二阶：
  $${\partial }_{t_2}{f}_i^{(0)}+{\partial }_{t_1}{f}_i^{(1)}+{\mathbf{c}}_i\cdot {\nabla }_1{f}_i^{(1)}+\frac{\Delta t}{2}{\left({\partial }_{t_1}+{\mathbf{c}}_i\cdot {\nabla }_1\right)}^2{f}_i^{(0)}=-\frac{1}{\tau }{f}_i^{(2)}$$

通过对上述各阶展开求和并对矩方法求和，可以得到连续介质的守恒方程和动量方程，从而证明在低马赫数和小Knudsen数的条件下，LBM能够恢复Navier-Stokes方程，其粘性系数与松弛时间之间的关系正如前文所述。

2、宏观方程的恢复

利用零阶和一阶展开，可以证明质量守恒方程：

$${\partial }_t\rho +\nabla \cdot (\rho \mathbf{u})=0$$

以及动量方程：

$${\partial }_t(\rho \mathbf{u})+\nabla \cdot (\rho \mathbf{uu})=-\nabla p+\nabla \cdot \left[\nu \left(\nabla \mathbf{u}+(\nabla \mathbf{u}{)}^T\right)\right]$$

其中压力 $p$ 与密度的关系可以由理想气体状态方程得到，对于 LBM 常取 $p=c_s^2\rho$。上述推导表明，LBM在合适的参数区间内，具有较高的数值精度，能够忠实恢复连续介质的基本动力学。

3、高阶修正项与精度讨论

在实际应用中，由于离散化的截断误差，高阶项会引入修正，这些修正项正是影响数值精度的重要来源。

例如，在一阶展开中，忽略高阶非线性项可能导致在高速流动或者边界效应明显的情形下误差增大。为了提高数值精度，研究者们提出了多种方法：

• 改进平衡分布函数的构造：在平衡态分布函数中加入更高阶的速度项展开，使得对流项的逼近更加精确；
• 采用高阶LBM模型：例如D3Q27等模型，通过增加离散速度数目来减少离散误差；
• 多松弛时间模型（MRT）：通过对不同矩进行不同松弛处理，分离剪切模量与体积模量的演化，进一步降低非物理耗散。

这些改进在理论上能够减小截断误差，但同时也会增加模型的复杂性，在实际编程实现中需要进行平衡取舍。

数值误差的来源与分析

LBM方法在数值计算中存在着多种误差来源，主要包括：

• 离散化误差：时间和空间的离散化会引入截断误差，其主要量级与 $\Delta t$ 和 $\Delta x$ 有关。
• 边界条件误差：在处理复杂边界时，常用的反弹边界（bounce-back）或修正边界条件会引入附加误差，特别是在处理曲面或斜边界时。
• 非线性累积误差：由于 LBM 更新过程中存在非线性碰撞过程，误差可能在多次迭代后累积，导致数值不稳定。
• 浮点运算误差：在实际计算机实现中，有限精度的浮点数运算也会对结果产生微小影响，尤其在长时间演化下可能积累成较大误差。

1、离散化误差与数值精度

在LBM中，离散化误差主要体现在两个方面：

• 时间离散误差：由于更新公式是显式时间推进方法，其局部截断误差与 $\Delta t$ 的高次幂有关。
• 空间离散误差：离散速度的选取导致在高阶矩的恢复上存在偏差，尤其在高雷诺数或大梯度区间，误差会较为明显。

对这些误差的分析主要依赖于Taylor展开法，通过对离散更新公式的展开，可以证明截断误差的主要项为 $O(\Delta t^2)$ 和 $O(\Delta x^2)$。例如，对于时间推进项：

$$f_i(\mathbf{x}+{\mathbf{c}}_i\Delta t,t+\Delta t)=f_i(\mathbf{x},t)+\Delta t\left({\partial }_t{f}_i+{\mathbf{c}}_i\cdot \nabla f_i\right)+\frac{\Delta t^2}{2}{\left({\partial }_t+{\mathbf{c}}_i\cdot \nabla \right)}^2{f}_i+\cdots$$

通过这种展开方法，可以定量分析离散化误差，并在模型设计时选择合适的网格和时间步长以确保总体误差控制在可接受范围内。

2、边界条件引起的误差

在LBM中，常见的边界条件包括：

• 反弹边界条件：简单易实现，但在某些情况下会引入一阶误差；
• Zou-He边界条件：适用于固定速度或压力边界，但在处理复杂几何边界时可能需要修正；
• 插值边界条件：通过对未知分布函数进行插值，可以提高边界精度，但计算量较大。

边界条件的处理直接影响到局部的稳定性和整体的数值精度。

数值模拟中，通常需要对不同边界条件进行对比分析，并通过数值实验确定最佳方案。

3、非线性累积误差与长期演化

在LBM数值计算中，由于碰撞过程的非线性，本身即使初始误差很小，在长时间迭代后也可能累积放大。

因此，在设计LBM数值模型时，不仅要关注局部截断误差，更要关注长期稳定性。常用的策略有：

• 自适应时间步长：在流动剧烈变化时缩小时间步长，从而降低误差累积；
• 误差校正算法：通过对局部误差进行修正（例如引入辅助场或修正项），实现全局精度的提升；
• 多尺度求解：将大尺度和小尺度误差分别处理，通过多重网格方法改善整体误差分布。

这些方法在实际应用中需要结合具体问题进行综合考虑，并通过数值实验验证其有效性。

实际应用讨论

在前述理论分析的基础上，我们讨论一些实际案例上关注的问题，以验证LBM在稳定性和精度方面的表现。实际数值实验通常会关注以下几点：

• 流场稳定性：在高雷诺数流动或者存在尖锐梯度的情况下，LBM数值解是否出现振荡或不稳定现象；
• 误差量化：通过与解析解或高精度数值解对比，量化LBM在不同网格和时间步长下的误差变化；
• 参数敏感性： 分析松弛时间、网格分辨率以及边界条件对整体误差和数值稳定性的影响。

例如，在模拟圆柱绕流问题时，通过改变松弛时间 $\tau$ 和网格分辨率，可以观察到：

• 当 $\tau$ 过小（接近 $\Delta t/2$）时，系统的数值振荡明显，且局部速度和压力场出现非物理解；
• 当网格分辨率较低时，边界附近的速度剖面存在较大误差，而采用高分辨率网格后，误差显著下降，但计算代价也随之增大。

在实际数值实验中，MRT模型在高雷诺数流动中能够更有效地抑制非物理振荡，从而获得更高的数值精度。同时，通过适当调整各矩的松弛系数，还可以进一步优化数值耗散特性，保证物理量的守恒性。

此外，对于边界条件的处理也是影响精度的重要因素。以反弹边界条件为例，虽然实现简单，但在某些几何复杂的情形下会引入一阶误差；因此，针对具体问题可能需要采用更高阶的边界修正方法，或者通过插值技术减小边界误差。

总结

展望未来，随着计算能力的不断提升和算法研究的不断深入，LBM 方法将继续在复杂流动模拟、湍流研究以及多物理耦合问题中发挥越来越重要的作用。如何在保证高精度的同时进一步提高数值稳定性，仍然是LBM研究的重要方向。未来的工作可能会聚焦于：

• 自适应网格与自适应时间步长： 根据局部误差自适应调整网格分辨率和时间步长，从而在计算效率和精度之间取得更好平衡；
• 高阶离散模型： 开发新型高阶离散速度模型，进一步降低截断误差，提高流场细节的捕捉能力；
• 耦合多物理场： 在热传导、电磁场、化学反应等多物理场问题中推广LBM方法，使其在复杂耦合问题中具有更广泛的适用性；
• 数据驱动与机器学习方法： 结合数据驱动方法和机器学习技术，对LBM中的参数进行优化调控，自动识别数值不稳定区域，并进行局部修正。

LBM方法以其独特的理论优势和数值实现的简便性，已在许多流体力学和多物理场问题中得到成功应用。对于稳定性和精度的深入研究，不仅有助于提升LBM的数值性能，也为新型数值方法的发展提供了重要参考。