LBM流体仿真

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

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

LBM Theory

Discretization

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

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

ρ (x, y) = \sum_{i = 1}^{9} f_{i} (x, y)

u (x, y) = \frac{1}{ρ (x, y)} \sum_{i = 1}^{9} f_{i} (x, y) v_{i}

其中每个格子在微观上9个方向的分布函数为 $f_{i}$ ，而速度 $v_{i}$ 则如下图所示。

Collide and Stream

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

f_{i}^{eq} = ρ w_{i} [1 + 3 (v_{i} \cdot u) + \frac{9}{2} (v_{i} \cdot u)^{2} - \frac{3}{2} u^{2}]

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

f_{i}^{new} = f_{i} + ω (f_{i}^{eq} - f_{i})

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

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

LBM Implementation

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

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

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

Lid-Driven Cavity Flow

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

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

R e = \frac{U L}{ν}

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

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

Kármán Vortex Street

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

R e = \frac{U D}{ν}

其中 $D$ 为圆柱直径。

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

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