简介：这份资源围绕格子Boltzmann方法（LBM）的流动数值模拟展开，面向流体力学、计算物理方向的学习者与研究人员，尤其适合希望用C++动手实现LBM、理解其算法流程的中高级读者。压缩包为tgz格式，体积约1.79MB，内含基于OpenLatticeBoltzmann项目0.7r1版本的源码与示例，涵盖二维与三维多种典型流动场景，如管道流动、绕流、水槽波浪生成等，可用于流场分析、湍流与传热研究的入门实践。资源以C++实现，兼顾计算效率与模块化设计，便于参数配置、几何导入与结果分析。目前已有804人学习下载，读者可借助其中的核心算法、数据结构与示例工程，快速上手LBM建模，并在此基础上进行二次开发，扩展至多相流或流固耦合等特定研究需求，是理解并应用该方法的一份实用参考。

1. 从零手写 LBM 流动模拟：一份能跑通的 C++ 源码拆解

顶盖驱动方腔流是 CFD 圈子里公认的“Hello World”，但真拿格子玻尔兹曼方法（LBM）从零写一遍，很多人卡在第一个时间步就翻车——密度场炸了、速度场发散、边界条件写反。我手里这份基于 LBM 的流动模拟 C++ 源码包，就是冲着这个痛点来的：它不依赖任何第三方 CFD 库，纯 C++ 实现 D2Q9 模型，包含碰撞、迁移、宏观量提取、边界处理四个核心模块，编译完直接出速度场和流线数据。适合两类人：一是想搞懂 LBM 底层数据结构的 C++ 开发者，二是需要快速搭一个可改可调的流动模拟骨架做算法验证的研究生。源码结构干净，没有过度封装，每个函数对应一个物理步骤，改参数不用翻三层继承。

2. LBM 核心原理与 D2Q9 数据结构：为什么用分布函数而不是解 N-S 方程

2.1 从玻尔兹曼方程到格子离散：LBM 到底在算什么

传统 CFD 解的是宏观连续方程，LBM 走的是另一条路：它模拟流体微观粒子的分布函数演化。你可以把每个格点想象成一个小房间，里面站着 9 个方向的“信使”，每个信使手里拿着一个数值，代表该方向上的粒子分布密度。每个时间步，信使先在自己房间里按碰撞规则交换手里的数值，然后整体往各自方向挪一格到隔壁房间。宏观的密度和速度，就是这 9 个数值的加权和。

这套逻辑的数学根基是玻尔兹曼输运方程，经过速度空间离散（D2Q9 就是 2 维 9 速度）、空间离散（格子）、时间离散后，得到格子玻尔兹曼方程：

f_i(x + e_i·Δt, t + Δt) = f_i(x, t) + Ω_i(f)

其中 f_i 是第 i 个方向的分布函数，e_i 是离散速度向量，Ω_i 是碰撞算子。最常见的碰撞模型是 BGK 近似，把碰撞简化成向平衡态松弛：

Ω_i = -(f_i - f_i^eq) / τ

τ 是松弛时间，直接决定运动黏度 ν = c_s²(τ - 0.5)Δt。c_s 是格子声速，D2Q9 里等于 1/√3。这个公式是整个模拟的命门——τ 必须大于 0.5，否则黏度为负，模拟必炸。很多新手第一次跑出 NaN，九成是 τ 设到了 0.5 以下。

D2Q9 的 9 个速度方向不是随便定的，它要满足质量、动量、能量守恒的矩条件。具体权重和速度向量如下表：

平衡态分布函数 f_i^eq 的公式是：

f_i^eq = w_i · ρ · [1 + (e_i·u)/c_s² + (e_i·u)²/(2c_s⁴) - u²/(2c_s²)]

这个公式看着复杂，代码里就是几行乘加。真正容易写错的是速度向量和权重的对应关系——索引 5 到 8 的对角方向，e_ix 和 e_iy 的符号组合必须和权重表严格对齐，错一个符号，流场就会出现诡异的各向异性。

2.2 源码里的数据结构：一维数组存二维分布

打开源码的lbm.h，核心数据结构就两个数组：

// lbm.h 核心数据结构节选 const int NX = 256; // x 方向格点数 const int NY = 256; // y 方向格点数 const int Q = 9; // D2Q9 模型速度方向数 // 分布函数存储：f[i][y][x] 用一维数组展平，方便缓存命中 // 索引计算：idx = i * NX * NY + y * NX + x double* f; // 当前时间步分布函数 double* f_new; // 迁移后的分布函数 double* rho; // 宏观密度场 double* ux; // 宏观速度 x 分量 double* uy; // 宏观速度 y 分量

这里有个设计选择值得说：为什么用一维数组而不是double f[Q][NY][NX]？因为 C++ 里多维数组的内存布局虽然是连续的，但编译器对指针别名和边界检查的处理会影响向量化。展平成一维后，手动算索引，编译器更容易生成 SIMD 指令。代价是代码可读性下降，但 LBM 每个时间步要遍历所有格点，性能收益是实打实的。

速度向量和权重通常定义成静态常量数组：

// 速度向量：ex[i], ey[i] 对应第 i 个方向 const double ex[Q] = {0, 1, 0, -1, 0, 1, -1, -1, 1}; const double ey[Q] = {0, 0, 1, 0, -1, 1, 1, -1, -1}; // 权重 const double w[Q] = {4.0/9, 1.0/9, 1.0/9, 1.0/9, 1.0/9, 1.0/36, 1.0/36, 1.0/36, 1.0/36};

注意 ex 和 ey 的排列顺序：索引 1 是 (1,0)，索引 2 是 (0,1)，索引 3 是 (-1,0)，索引 4 是 (0,-1)，然后 5 到 8 是对角。这个顺序不是唯一的，但一旦定了，后面所有循环都必须按这个顺序来。我见过有人把对角方向的顺序换了，结果流场出现棋盘格振荡，查了一整天。

2.3 碰撞与迁移：两个循环的顺序不能反

LBM 的一个时间步分两步：碰撞（collision）和迁移（streaming）。碰撞在本地格点完成，只涉及当前格点的 9 个分布函数；迁移把碰撞后的分布函数沿各自速度方向送到邻居格点。

// 碰撞步骤：计算平衡态并松弛 for (int y = 0; y < NY; y++) { for (int x = 0; x < NX; x++) { int idx = y * NX + x; double rho_local = 0.0, ux_local = 0.0, uy_local = 0.0; // 先算宏观量 for (int i = 0; i < Q; i++) { double fi = f[i * NX * NY + idx]; rho_local += fi; ux_local += fi * ex[i]; uy_local += fi * ey[i]; } ux_local /= rho_local; uy_local /= rho_local; // 再算平衡态并碰撞 for (int i = 0; i < Q; i++) { double eu = ex[i] * ux_local + ey[i] * uy_local; double feq = w[i] * rho_local * (1.0 + eu / (cs * cs) + (eu * eu) / (2.0 * cs * cs * cs * cs) - (ux_local * ux_local + uy_local * uy_local) / (2.0 * cs * cs)); int idf = i * NX * NY + idx; f_new[idf] = f[idf] - (f[idf] - feq) / tau; } } }

碰撞循环里，宏观量的计算必须在平衡态之前，因为 feq 依赖 ρ 和 u。这里有个性能细节：f[i * NX * NY + idx]这种索引方式，当 i 变化时跨步很大，缓存不友好。更优的做法是把分布函数按格点连续存储，即f[idx * Q + i]，这样内层循环访问连续内存。源码里用的是前者，因为更直观，但如果你要跑大网格，建议改成后者，实测有 20% 到 30% 的性能差距。

迁移步骤紧接着碰撞：

// 迁移步骤：把 f_new 沿速度方向送到邻居 for (int y = 0; y < NY; y++) { for (int x = 0; x < NX; x++) { for (int i = 0; i < Q; i++) { int x_src = x - (int)ex[i]; // 来源格点 int y_src = y - (int)ey[i]; // 周期性边界：越界则绕回 if (x_src < 0) x_src += NX; if (x_src >= NX) x_src -= NX; if (y_src < 0) y_src += NY; if (y_src >= NY) y_src -= NY; f[i * NX * NY + y * NX + x] = f_new[i * NX * NY + y_src * NX + x_src]; } } }

迁移的方向容易搞反：x_src = x - ex[i]表示“从上游拿”，而不是“往下游送”。两种写法在周期性边界下等价，但在有壁面的情况下，边界处理逻辑完全不同。源码用的是“拉”的方式，配合边界条件更自然。

提示：碰撞和迁移必须用两个独立的数组，不能原地更新。如果碰撞后直接写回 f，迁移时就会用到已经被覆盖的数据，结果完全错误。这是 LBM 实现里最经典的坑之一。

3. 编译、运行与参数调优：从源码到速度场云图

3.1 编译环境与依赖：纯 C++ 标准库，零外部依赖

这份源码的一个好处是没有任何第三方依赖，不需要 MPI、不需要 OpenMP、不需要 CUDA，只要有 C++11 以上的编译器就能跑。我一般用 g++ 直接编译：

# 编译：开启 O3 优化和自动向量化 g++ -O3 -march=native -std=c++11 -o lbm_sim main.cpp lbm.cpp -lm # 运行：默认参数是 256x256 网格，10000 步 ./lbm_sim

-O3是必须的，LBM 的计算密集度很高，不开优化跑 256×256 网格 10000 步大概要几分钟，开了之后能压到几十秒。-march=native让编译器针对当前 CPU 生成 AVX 指令，对碰撞循环里的乘加运算提升明显。如果你在 ARM 机器上跑，把-march=native换成-mcpu=native。

源码的main.cpp里参数是硬编码的，但都集中在文件开头：

// main.cpp 可调参数区 const int MAX_ITER = 10000; // 总时间步数 const double tau = 0.8; // 松弛时间，决定黏度 const double U_lid = 0.1; // 顶盖驱动速度（格子单位） const int OUTPUT_INTERVAL = 1000; // 每隔多少步输出一次场数据

tau = 0.8对应运动黏度 ν = (0.8 - 0.5)/3 ≈ 0.1。这个值不算小，流场会比较“稠”，适合看稳态涡结构。如果你想看更接近湍流的瞬态，把 tau 降到 0.55 到 0.6 之间，但要注意稳定性——tau 越接近 0.5，马赫数稍微大一点就会发散。

3.2 顶盖驱动方腔的边界条件实现

顶盖驱动方腔的边界条件分三类：顶盖是运动壁面，左右和底部是静止壁面，四个角是奇点。源码里用的是“反弹格式”处理静止壁面，顶盖用“动量注入”方式。

// 边界处理：在迁移之后调用 void applyBoundaryConditions(double* f, double* f_new) { // 顶盖：y = NY-1，速度 U_lid 沿 x 正方向 int y = NY - 1; for (int x = 0; x < NX; x++) { int idx = y * NX + x; double rho_wall = 0.0; // 先算壁面密度 for (int i = 0; i < Q; i++) rho_wall += f[i * NX * NY + idx]; // 对指向壁面内部的分布函数，用平衡态重构 // 这里以方向 4（0,-1）和方向 7、8 为例 double ux_wall = U_lid, uy_wall = 0.0; for (int i : {4, 7, 8}) { double eu = ex[i] * ux_wall + ey[i] * uy_wall; double feq = w[i] * rho_wall * (1.0 + eu / (cs * cs) + (eu * eu) / (2.0 * cs * cs * cs * cs) - (ux_wall * ux_wall + uy_wall * uy_wall) / (2.0 * cs * cs)); f[i * NX * NY + idx] = feq; } } // 左右和底部：标准反弹格式 // 以左壁 x=0 为例，方向 1（1,0）的分布函数来自方向 3（-1,0） for (int y = 0; y < NY; y++) { int idx = y * NX + 0; f[1 * NX * NY + idx] = f[3 * NX * NY + idx]; } // 右壁、底部类似，略 }

顶盖的处理是新手最容易翻车的地方。方向 4、7、8 是“指向壁面内部”的方向，它们的分布函数不能从邻居迁移得到，必须用壁面速度重构。重构时用的密度是壁面格点的局部密度，不是全局平均密度——用错会导致质量不守恒，长时间跑下来总质量漂移。

反弹格式的原理很简单：粒子撞到静止壁面后原路返回。所以左壁 x=0 处，方向 1（向右）的分布函数等于方向 3（向左）的分布函数。但要注意，这个赋值发生在迁移之后，因为迁移已经把方向 3 的粒子送到了壁面格点。

3.3 输出格式与可视化：用 Python 画速度场

源码每 1000 步输出一次ux、uy、rho三个二进制文件，格式是简单的 double 数组，按行优先排列。用 Python 读进来画图：

import numpy as np import matplotlib.pyplot as plt NX, NY = 256, 256 # 读取二进制场数据 ux = np.fromfile('ux_10000.bin', dtype=np.float64).reshape(NY, NX) uy = np.fromfile('uy_10000.bin', dtype=np.float64).reshape(NY, NX) rho = np.fromfile('rho_10000.bin', dtype=np.float64).reshape(NY, NX) # 画速度幅值云图 speed = np.sqrt(ux**2 + uy**2) plt.figure(figsize=(8, 6)) plt.imshow(speed, origin='lower', cmap='jet') plt.colorbar(label='Velocity magnitude') plt.title('Lid-driven cavity flow, tau=0.8, 10000 steps') plt.xlabel('X') plt.ylabel('Y') plt.savefig('velocity_field.png', dpi=150) plt.show() # 画流线 plt.figure(figsize=(8, 6)) plt.streamplot(np.arange(NX), np.arange(NY), ux, uy, density=2, color=speed, cmap='jet') plt.title('Streamlines') plt.savefig('streamlines.png', dpi=150) plt.show()

reshape(NY, NX)的顺序不能反，因为 C++ 输出时是 y 外层循环、x 内层循环，对应行优先。如果画出来的图上下颠倒，把origin='lower'改成origin='upper'就行。

判断模拟是否收敛，看两个指标：一是总动能随时间的变化，稳态后应该趋于常数；二是涡心位置，顶盖驱动方腔在 Re=1000 左右时主涡心大约在 (0.53, 0.57) 附近（归一化坐标）。如果涡心跑到角落或者出现多个涡，说明边界条件或者 tau 有问题。

注意：源码输出的二进制文件没有头信息，纯 double 数组。如果你改了 NX 或 NY，Python 读取时的 reshape 参数必须同步改，否则数据会错位。

4. 避坑与排查：LBM 模拟发散的五个血泪教训

4.1 现象：跑几百步后密度场出现 NaN

原因：松弛时间 tau 太接近 0.5，或者顶盖速度 U_lid 太大导致马赫数超限。LBM 的弱可压缩性要求马赫数 Ma = U_lid / c_s 小于 0.3 左右，c_s = 1/√3 ≈ 0.577，所以 U_lid 不应超过 0.17。源码默认 0.1 是安全的，但有人为了“加快收敛”调到 0.3，基本几十步就炸。

解决：把 U_lid 降回 0.1 以下，同时检查 tau 是否大于 0.55。如果必须用高雷诺数，应该增加网格分辨率而不是提高速度。

4.2 现象：流场出现棋盘格振荡，速度场一正一负交替

原因：碰撞和迁移的循环顺序写反了，或者迁移时用了原地更新。LBM 的碰撞必须在本地完成，迁移必须从旧数组读、往新数组写。如果碰撞后直接覆盖了 f，迁移时读到的就是已经碰撞过的数据，物理上完全错误。

解决：严格用两个数组，碰撞写 f_new，迁移从 f_new 读到 f。源码里f和f_new的角色在每个时间步末尾交换，不要手动去改。

4.3 现象：总质量随时间缓慢漂移，跑久了密度整体偏移

原因：边界条件的密度重构用了全局平均密度而不是局部密度。顶盖壁面的平衡态重构需要该格点的局部密度，如果用全局 ρ0 = 1.0，壁面附近的质量就不守恒。

解决：在边界处理函数里，先遍历该格点的 9 个分布函数求和得到局部密度，再用这个密度算 feq。源码里rho_wall就是这么来的。

4.4 现象：编译通过但运行时报段错误

原因：数组索引越界。最常见的是迁移步骤里x_src或y_src算出负数后没有正确加 NX/NY，或者加了之后仍然越界。另一个可能是f和f_new分配的内存大小是Q * NX * NY，但索引时用了i * NX * NY + y * NX + x，如果 i 到了 9 就溢出了。

解决：在迁移循环里加断言assert(x_src >= 0 && x_src < NX)，编译时加-g -fsanitize=address跑一遍，能直接定位越界位置。

4.5 现象：模拟结果和文献对不上，涡心位置偏差大

原因：网格分辨率不够。顶盖驱动方腔在 Re=1000 时，256×256 网格勉强够用，但如果 Re 更高或者要看更精细的涡结构，需要 512×512 甚至 1024×1024。另一个原因是收敛判据太松，10000 步可能还没到稳态。

解决：先跑一个粗网格看趋势，确认边界条件没问题后，再加密网格。收敛判据用速度场的变化率，比如max|u_new - u_old| < 1e-6，而不是固定步数。

5. 进阶玩法：从方腔流到多孔介质，改三处就能跑

源码的骨架是通用的，换成其他流动场景只需要改三个地方：边界条件函数、初始场、以及输出逻辑。我拿多孔介质渗流举个例子——这是 LBM 的经典应用场景，比方腔流更贴近实际工程。

多孔介质的做法是在流场里随机放置一批固体格点，粒子撞到固体格点就反弹。不需要改碰撞和迁移，只需要在边界处理里加一个判断：

// 多孔介质：固体格点反弹 // 先用一个 bool 数组标记固体 bool* is_solid = new bool[NX * NY]; // 随机初始化，孔隙率 0.7 for (int i = 0; i < NX * NY; i++) is_solid[i] = (rand() / (double)RAND_MAX) > 0.7; // 在迁移之后，对固体格点做反弹 for (int y = 0; y < NY; y++) { for (int x = 0; x < NX; x++) { int idx = y * NX + x; if (!is_solid[idx]) continue; // 固体格点：所有方向的分布函数反弹 // 方向 i 的分布函数来自相反方向 i_op int i_op[Q] = {0, 3, 4, 1, 2, 7, 8, 5, 6}; for (int i = 0; i < Q; i++) { f[i * NX * NY + idx] = f_new[i_op[i] * NX * NY + idx]; } } }

i_op数组是方向索引的相反映射：0 对 0，1 对 3，2 对 4，5 对 7，6 对 8。这个映射表必须和速度向量的定义严格一致，错一个就漏质量。

压力驱动用周期性边界加一个体积力实现。在碰撞步骤里给速度加一个增量：

// 在碰撞循环里，算完宏观速度后加体积力 ux_local += force_x * tau; // force_x 是 x 方向体积力 uy_local += force_y * tau;

体积力的大小控制雷诺数，force_x = 1e-5对应很慢的渗流，1e-3就快很多。跑够步数后，统计平均速度，用达西定律反算渗透率：K = ν * u_avg / force_x。这个值和孔隙率、固体分布有关，可以拿来验证代码是否正确。

验证方法很简单：把孔隙率设成 1.0（没有固体），结果应该退化成方腔流或者周期性槽道流。如果对不上，说明反弹逻辑或者体积力实现有问题。我一般会先跑这个退化测试，确认无误后再加固体。

从那以后我每次改 LBM 代码，都强制走一遍“退化测试—粗网格—细网格”三步，再急也不跳过。希望帮到你。

本文还有配套的精品资源，点击获取