LBM流体仿真

Created in February 09, 2025

2025   ·   CG   Simulation

基于LBM方法实现二维流体仿真。

格子Boltzmann方法(Lattice Boltzmann Method, LBM)是一种基于介观(mesoscopic)模拟尺度的计算流体力学(Computational Fluid Dynamics, CFD)方法。和传统CFD方法相比，LBM具有介于微观分子动力学模型和宏观连续模型的介观模型特点，因此具备流体相互作用描述简单、复杂边界易于设置、易于并行计算等优势。在这篇博客中，我们将从零开始一步步构建一个LBM二维仿真器，并利用NVIDIA Warp库实现流体实时仿真。

LBM Theory

Discretization

LBM是基于网格的计算方法。在宏观层面，我们首先需要将流场空间离散成网格，每个网格记录下流体在该位置上的密度\(\rho\)以及速度\(\mathbf{u}\)。这些宏观量是我们最终希望从模拟中获得的结果，它们描述了流体的整体行为。

而在微观层面上，每个格子还包含了一组粒子分布函数(distribution function)，它们描述了粒子在各个微观方向上流动的概率。以二维LBM最常用的D2Q9模型为例，流体在每个格子上的密度与流速可以由9个微观方向上的分布函数来描述:

\[\rho(x, y) = \sum_{i=1}^9 f_i (x, y)\] \[\mathbf{u} (x, y) = \frac{1}{\rho (x, y)} \sum_{i=1}^9 f_i (x, y) \mathbf{v}_i\]

其中每个格子在微观上9个方向的分布函数为\(f_i\)，而速度\(\mathbf{v}_i\)则如下图所示。

Collide and Stream

LBM在进行流体仿真时，核心操作包括两个步骤：碰撞(collide)和流动(stream)。这里”碰撞”是指在每个格子上的粒子实现均衡(equilibrium)。在D2Q9模型中，均衡分布的计算公式为：

\[f^{\text{eq}}_i = \rho w_i \bigg[ 1 + 3 (\mathbf{v}_i \cdot \mathbf{u}) + \frac{9}{2} (\mathbf{v}_i \cdot \mathbf{u})^2 - \frac{3}{2} \mathbf{u}^2 \bigg]\]

其中，\(w_i\)为每个微观方向上速度的权重。均衡分布计算公式的推导过程可以参考wiki，这里不多做赘述。计算出均衡分布后就更新每个微观方向上的分布函数

\[f^{\text{new}}_i = f_i + \omega (f^{\text{eq}}_i - f_i)\]

上式中\(\omega = \frac{1}{\tau}\)为弛豫参数，\(\tau\)为弛豫时间，它反映了流体的黏性。可以看到，LBM中的碰撞步骤实际上是对原始分布函数\(f_i\)以及均衡分布\(f^{\text{eq}}_i\)进行线性插值。

计算出新的分布函数\(f^{\text{new}}_i\)后，粒子将沿微观速度方向传播到相邻的网格点，这一过程称为流动(stream)。通过交替进行碰撞和流动步骤，我们就能实现流体仿真。

LBM Implementation

LBM进行流体仿真的基本步骤可以概括为：

  1. 初始化流场
  2. 循环迭代:
    1. 计算碰撞
    2. 计算迁移
    3. 更新流场

除了上述基本步骤外，进行流体仿真时还需要计算相关的物理量并施加边界条件等操作。这些步骤需要根据具体的仿真问题进行处理，相关细节将在后续的实例中进一步说明。整体代码可以参考git仓库。

Lid-Driven Cavity Flow

顶盖驱动流(lid-driven cavity flow)是流体动力学中常用的基准问题，用于验证计算方法。我们可以设想一个方形舱体，舱内充满流体，并在顶部施加一个恒定的水平方向速度。此时，流体将在水平流速的驱动下流动，直到达到稳定状态。

顶盖驱动流的雷诺数计算公式为：

\[Re = \frac{U L}{\nu}\]

其中流体的特征速度为\(U\)，舱体尺寸为\(L\)，流体粘度为\(\nu\)。在仿真过程中，通常会给定流体的雷诺数、特征速度和舱体尺寸，进而计算出粘度\(\nu\)进行仿真。而舱体的边界条件为除了顶部边界外其它边界的流体流速均为0。

运行CavityFlow.py可以得到仿真结果如下。

Kármán Vortex Street

卡门涡街(Kármán vortex street)是自然界中常见的流体现象，通常出现在流体经过钝体时，会产生周期性的漩涡脱落。要模拟卡门涡街，最直观的方法是求解圆柱绕流问题。我们可以在矩形流场中放置一个圆柱，当流体流经圆柱时，可能会形成卡门涡街现象。实际上，卡门涡街的产生与流体的雷诺数密切相关，只有当雷诺数足够大时，才会观察到这一现象。对于圆柱绕流问题，雷诺数的计算公式为：

\[Re = \frac{U D}{\nu}\]

其中\(D\)为圆柱直径。

和顶盖驱动流相比，圆柱绕流的边界条件设置要稍微复杂一些。对于矩形流场的边界，我们需要设置上下边界上的流速为0表示无滑移边界；在左边入口边界上设置流场的来流速度为固定值\(U\)；而右侧出口边界则设置为Neumann边界条件，即流体速度等于左侧对应格子的速度(速度梯度为0)，表示出口处流动无扰动。此外还需要设置圆柱表面上的边界条件，当流体接触到圆柱表面时会速度会变为0，表示无滑移边界。

运行KarmanVortexStreet.py可以得到仿真结果如下。

MRT

Reference