The Fluid Simulator That Doesn’t Solve the Fluid Equations

不解流体方程的流体模拟器 | Towards Data Science

物理学

不解流体方程的流体模拟器

实践教程：从统计力学和动力学理论到C++实现与超级计算机扩展

Ferran Alia

2026年7月25日

13分钟阅读

分享

使用Gemini生成的图片

上方图片展示了卡门涡街（在快速流动中任何钝体后形成的交替涡旋）。你可能在桥墩、飞机机翼和烟囱后面见过这种模式。这是流体力学中最具辨识度的现象之一。

我生成这张图片时，没有解过任何流体方程。

相反，这种格子玻尔兹曼方法（LBM）模拟将流体视为微观层面的真实形态：相互碰撞的粒子统计分布。宏观行为（涡旋、压力梯度、粘性阻尼）作为结果自然出现，而非人为输入。

LBM在计算物理中占据着独特的位置，介于分子动力学模拟和纳维-斯托克斯求解器之间，在介观尺度运行。在本文中，我将从基本原理推导LBM，用约200行C++代码实现，并在MareNostrum 5超级计算机上运行。完整源代码和简化版Python/NumPy实现可在GitHub上找到。

简要总结：LBM用离散速度格子上的简单更新规则取代难以处理的纳维-斯托克斯偏微分方程。在正确极限下可恢复纳维-斯托克斯方程，而非假设成立。这使LBM在理论上优雅、实践中快速且易于并行化。它在低马赫数、亚音速范围内表现最佳，在极低粘度下会变得数值不稳定。

为什么直观方法会很痛苦

流体模拟的标准路径是通过纳维-斯托克斯方程，这些宏观守恒定律支配着流体力学中几乎所有我们关心的现象。除了少数玩具几何结构外，它们在分析上难以处理，因此实践中需要在网格上离散化并数值求解：有限体积法、有限差分法或有限元法。

这些方法虽然有效，但对复杂几何结构存在隐藏成本。对于曲面边界（如圆柱体、车身、多孔介质），要么需要昂贵且难以自动化的贴体网格，要么需要在壁面附近进行专门的模板修正。每次几何变化都需要重新生成网格。这是叠加在物理问题上的工程问题。

LBM几乎完全规避了这个问题。原因在理解其起点后会变得清晰。

不同的起点：统计力学

LBM不是问“该点的流速场是什么？”，而是问“该点各个方向移动的粒子数量是多少？”

核心对象是单粒子分布函数f(x, v, t)，它告诉你在时间t，位置x附近以速度v移动的粒子预期数量。我们实际关心的宏观场只是f的矩：

$$ \rho(\mathbf{x},t) = m\int f\, d^3v \qquad \rho\,\mathbf{u}(\mathbf{x},t) = m\int \mathbf{v}\,f\, d^3v $$

在全局平衡状态（均匀密度、无整体流动）下，f呈现熟悉的麦克斯韦-玻尔兹曼形式：

$$ f^{eq}(\mathbf{v}) = \frac{\rho}{(2\pi RT)^{D/2}} \exp\!\left(-\frac{|\mathbf{v} – \mathbf{u}|^2}{2RT}\right) $$

远离平衡态时，$f$ 会呈现其他形状。支配其演化过程的方程是玻尔兹曼方程：

$$ \frac{\partial f}{\partial t} + \mathbf{v}\cdot\nabla_\mathbf{x} f = \Omega[f] $$

左边部分仅表示输运过程：粒子在空间中流动时会携带其速度分布。右边部分 $\Omega[f]$ 是碰撞算子，通过在不同速度区间重新分配粒子，将 $f$ 驱动回平衡态 $f^{\text{eq}}$。问题是 $\Omega[f]$ 是一个五维积分，需要对所有可能的碰撞伙伴和散射角度进行积分，包含了完整的微观碰撞物理过程。对于任何实际模拟来说，这完全无法处理。

圆柱流的速度场。作者提供

一行奇迹：BGK近似

1954年，Bhatnagar、Gross和Krook提出了一种巧妙的简化方法。他们没有详细建模每一次碰撞，而是捕捉碰撞的净效应：碰撞以与当前偏离平衡态程度成正比的速率，将 $f$ 驱向平衡态。

$$ \Omega[f] \approx -\frac{1}{\tau}\bigl(f – f^{eq}\bigr) $$

这就是 BGK 近似，它将一个难以处理的五维积分简化为标量乘法。参数 $\tau$ 是弛豫时间，表示系统在受到扰动后恢复到局部平衡态所需的平均时间。$\tau$ 值大意味着弛豫过程缓慢；$\tau$ 值小意味着弛豫过程迅速。

使 BGK 近似超越便捷技巧的关键在于 $\tau$ 的物理意义。通过 Chapman-Enskog 展开（一种系统性的微扰分析），可以证明当 Knudsen 数很小时（即每段平均自由路径发生大量碰撞），BGK 方程会精确地恢复出纳维-斯托克斯方程，其运动学粘度为：

$$ \nu = c_s^2\!\left(\tau – \frac{1}{2}\right) $$

注意这个关系式是在无量纲晶格单位下成立的；转换为物理单位需要单独进行长度和时间尺度校准。但在模拟中，$\tau$ 是唯一的物理调节参数：增大 $\tau$ 可获得粘稠流动；将 $\tau$ 推近 0.5 可模拟高雷诺数流动（不过低于约 0.55 的值可能导致数值不稳定）。

关键在于，这正是为什么 LBM（格子玻尔兹曼方法）不仅仅是纳维-斯托克斯方程的伪装。与纳维-斯托克斯方程的联系是通过推导得出的，而非假设：你正在求解一个动力学方程，该方程在适当极限下可被证明收敛于纳维-斯托克斯行为。这是统计力学中一个非常普遍的主题。

从连续到离散：D2Q9晶格

BGK方程仍然存在于连续速度空间中。为了模拟它，我们需要一组有限的速度。方法是对 $f^{\text{eq}}$ 在速度 $u$ 的幂级数上进行泰勒展开，该展开在 $u$ 远小于晶格声速（低马赫数假设）时有效：

$$ f_i^{eq} = \rho\, w_i\!\left[1 + 3(\mathbf{c}_i\cdot\mathbf{u}) + \frac{9}{2}(\mathbf{c}_i\cdot\mathbf{u})^2 – \frac{3}{2}|\mathbf{u}|^2\right] $$

连续速度 $c$ 被有限的离散速度 $c_i$ 取代，每个速度方向都携带一个权重 $w_i$，这些权重的选择是为了正确恢复速度矩。对于二维模拟，标准选择是 D2Q9：2维空间，9个速度方向（静止状态、4个轴向方向和4个对角方向）。

作者提供

权重由各向同性决定：静止粒子权重为 $4/9$，轴向方向权重为 $1/9$，对角方向权重为 $1/36$。晶格声速为 $c_s = 1/\sqrt{3}$。算法中的其他所有内容都由这九个数值推导得出。

算法：两个步骤，无限循环

掌握离散平衡后，完整的LBM算法在每个时间步仅需执行四个操作，通过循环直至收敛：

• 宏观量恢复：从当前分布函数f_i计算每个节点的ρ和u
• 碰撞：以1/τ速率将每个f_i松弛至f_i^eq（纯局部操作，无需邻居通信）
• 流动：将每个分布函数沿其传播方向传递至相邻节点
• 边界条件：在域边界和固体壁面重建缺失的分布函数

碰撞步骤是所有物理过程发生的核心环节。流动步骤仅涉及纯运动学的簿记操作。这两个步骤的交替作用产生了流体行为。

void collide() {
    const double omega = 1.0 / tau;
    #pragma omp parallel for collapse(2) schedule(static)
    for (int x = 0; x < Nx; x++) {
        for (int y = 0; y < Ny; y++) {
            if (solid[x][y]) continue;
            double feq[9];
            computeEquilibrium(rho[x][y], ux[x][y], uy[x][y], feq);
            for (int i = 0; i < 9; i++)
                f[x][y][i] -= omega * (f[x][y][i] - feq[i]);
        }
    }
}

流动实现对性能影响巨大。简单的推送模式（每个线程在(x,y)处将出站分布函数写入相邻节点）会导致伪共享：线程写入相邻内存区域，反复访问相同缓存行导致频繁失效。下述拉取（收集）模式解决了这个问题，因为每个线程仅从上游邻居读取数据，仅写入自己的节点。

void stream() {
    #pragma omp parallel for collapse(2) schedule(static)
    for (int x = 0; x < Nx; x++) {
        for (int y = 0; y < Ny; y++) {
            if (solid[x][y]) continue;
            for (int i = 0; i < 9; i++) {
                const int xsrc = x - EIX[i];
                const int ysrc = y - EIY[i];
                const bool oob = (xsrc < 0 || xsrc >= Nx ||
                                  ysrc < 0 || ysrc >= Ny);
                f_new[x][y][i] = (oob || solid[xsrc][ysrc])
                    ? f[x][y][OPP[i]]    // bounce-back from wall
                    : f[xsrc][ysrc][i];  // stream from upstream neighbour
            }
        }
    }
}

由于碰撞完全局部化且流动仅需最近邻读取，整个算法具有惊人的并行性：每个节点独立更新，无需全局求解、无需隐式系统、无需矩阵求逆。这种特性可无缝映射到多核CPU、GPU和分布式集群；基于CUDA的LBM实现通常能实现接近线性扩展，直至达到内存带宽限制。

边界条件：免费获得复杂几何处理

每次流动步骤后，固体壁面相邻的节点会缺失部分入射分布函数，因为这些分布函数本应从固体内部到达，而固体中没有流体可以流动。如何填补这些缺失的分布函数定义了边界条件。

反弹边界条件（Bounce-back BC）。作者：By author

反弹边界条件处理无滑移壁面。任何流入固体节点的分布函数都会沿相反方向反射：f_ī = f*_i。物理图像是粒子在壁面处反转速度，无滑移条件（边界速度为零）会自动出现，无需显式速度约束。物理壁面位于最后一个流体节点与第一个固体节点的中间位置，使该方案具有二阶空间精度。

优雅的结果是，任何障碍物形状都只是一个布尔标志。要添加一个圆柱体，只需将相关节点标记为实体。反弹规则会自动处理每个标记节点的边界物理，无需网格、无需模板修正，也无需在不同几何形状之间修改代码。

在入口处，我们采用Zou-He方案：施加固定速度uw，通过剩余已知分布函数和施加速度，一致地重建三个向右指向的缺失分布函数（f_1, f_5, f_8），从而在不独立施加的情况下恢复正确的入口密度。

在出口处，零梯度条件会复制倒数第二列的所有分布函数，使流体自由流出域而不反射虚假波向上游传播。

这些模拟中使用的边界条件。作者提供

结果：两种流动模式

模拟域是一个400×400的格子，其中心位于域长度四分之一处，略微垂直偏移以打破对称性，放置一个半径为40个格子单位的圆形障碍物。入口速度为uw = 0.1个格子单位，远低于cs ≈ 0.577（整个过程满足低马赫数假设）。两个τ值揭示了本质上不同的流动模式。

当τ = 0.75（Re ≈ 96）时，尾流稳定且对称。分布函数在所有位置都接近局部平衡态，因为碰撞足够剧烈，能在扰动发展之前将其抑制。流动是粘性且层流的。

当τ = 0.55（Re ≈ 480）时，情况完全改变。

τ = 0.55时的流动。作者提供

涡旋交替从上下表面脱落，形成卡门涡街。从统计力学的角度来看，这是关键见解：涡旋脱落是由分布函数f偏离麦克斯韦-玻尔兹曼平衡形状驱动的。较小的τ意味着碰撞使f向平衡态弛豫得更慢，允许非平衡部分f − f^eq增长到足够大，从而反馈到宏观动量场并维持振荡。完全处于平衡态的流体将不会出现这种结构。

布尔标志方法意味着改变几何形状无需任何成本。这里是相同的求解器，未经修改，应用于一排三个方形障碍物：

尾流相互耦合，涡街相互作用，丰富的多体流动结构涌现。代码唯一的变化是将不同节点标记为实体。

扩展到MareNostrum 5

在碰撞和流动物理循环中同时使用#pragma omp parallel for，求解器可以以最小的代价扩展到多个核心。关键选择是流动物理模式：拉取模式确保每个线程独占其输出内存，消除导致朴素推送实现早期停滞的虚假共享。

我们在巴塞罗那超级计算中心的MareNostrum 5全节点上进行了基准测试：2× Intel Xeon Platinum 8480+（Sapphire Rapids），跨两个NUMA域的112个核心。

实际与理想加速比对比。作者提供

三种性能区间清晰地显现出来。从1到32线程时，效率保持在75%以上，达到24.6倍的加速比，这真正实现了对内存密集型模板代码的有效扩展。当线程数增加到64时，效率下降至42%，因为内存控制器达到饱和状态（约30MB的分布函数数据的读写速度超过了总线的承载能力）。当线程数达到112时，加速比反而下降至23.5倍，因为跨物理插槽的内存访问导致一半的内存操作需要通过跨插槽互连，这在带宽已饱和的情况下进一步增加了延迟。对于400×400网格而言，单个插槽上32个核心是实际的最佳选择。

一个值得关注的数据：在完成流式拉取重构后，单线程基准性能从原来的113秒降至31.7秒。这种3.6倍的性能提升仅通过代码修改就实现，且是在未涉及任何并行处理的情况下，仅在单核上完成。内存访问模式对性能的影响远早于硬件层面的优化。

何时使用LBM（以及何时不适用）

当几何结构复杂或随时间变化、需要在大规模并行硬件上运行、或需要模拟介观尺度现象（如壁面附近需要经验修正的滑移速度）时，LBM是合适的工具。它也是多相流模拟的绝佳基础，这将是本系列文章下一篇文章的主题。

对于高马赫数流动（当马赫数超过0.3时低阶泰勒展开失效）、极高雷诺数（τ > 0.5的稳定性要求限制了粘度的最小值，除非使用更先进的碰撞算子）或需要从一开始就精确校准物理单位的情况，LBM则不是合适的工具。

下一步：Shan-Chen模型的多相流动

目前描述的求解器仅处理单组分、单相流体。许多最有趣的物理问题涉及两种不相溶的流体（如油和水、上升气泡、气流中的液滴动力学）。

Shan-Chen模型通过引入短程粒子间作用力，扩展了LBM以处理这些情况，从而产生具有相变的有效状态方程。界面张力和相分离是从碰撞统计中自然产生的，而不是作为边界条件强加的。在下一篇文章中，我将推导Shan-Chen模型，将其作为本求解器的适度扩展实现，并展示当两种流体争夺同一晶格时自发相分离的视觉效果。

结语

LBM最让我感到满意的，是它揭示了不同尺度之间的关系。纳维-斯托克斯方程似乎很基础：它们是流体力学方程，在我们能看见和触摸的尺度上书写。但它们实际上并不基础。它们是大量粒子碰撞的统计结果，在适当的极限下可以从动力学描述中恢复。

LBM使这种关系具体化并可计算。粘度不是你设定的参数，而是碰撞时间尺度；卡门涡街不是你编程的流体不稳定性，而是当碰撞无法跟上惯性时自然产生的现象。模拟中的每一个宏观现象，在其下的统计力学中都有直接的微观解释。

这种尺度间的桥梁，通过约200行C++代码具体实现，正是LBM的核心价值所在。

参考文献：

对于希望深入研究的读者，以下按技术难度递增顺序列出了权威参考文献：

Mohamad, A.A. (2019). 《格子玻尔兹曼方法：原理与工程应用（含代码示例）》Springer出版社。最易入门的读物，全书包含实例推导和代码演示。

Krüger, T. 等 (2017). 《格子玻尔兹曼方法：原理与实践》Springer出版社。综合性参考著作，涵盖推导过程、边界条件、多相模型及GPU实现等完整内容。

Bhatnagar, P.L., Gross, E.P. & Krook, M. (1954). 《气体碰撞过程的模型》《物理评论》94(3), 511–525。原始BGK论文。

Chapman, S. & Cowling, T.G. (1970). 《非均匀气体的数学理论》剑桥大学出版社。完整严谨的Chapman-Enskog展开理论解析。

本文的完整C++源代码、OpenMP基准测试脚本以及Python/NumPy参考实现均可在GitHub获取。

作者信息

查看Ferran Alia的全部文章

空气动力学

,

C++

深度解析

流体力学

仿真模拟

分享本文

• Facebook分享

• LinkedIn分享

• X平台分享

Towards Data Science是社区出版物。提交您的洞见以触达全球读者，并通过TDS作者支付计划获得收益。

请将href更新为实际投稿链接

为TDS撰写文章

✦ 结束CTA ✦