GPU Gems 2 · Chapter 47. Flow Simulation with Complex Boundaries

直接解 Navier–Stokes 往往很贵，而图形应用通常更需要“足够可信、可以实时更新”的流场。LBM 把连续流体换成规则格点：每个格点保存几组沿离散方向移动的 packet，局部规则反复更新，就能在宏观上产生密度和速度场。

本章的难点不只是流体公式，还包括障碍物和动态边界。复杂网格必须变成 lattice 能理解的离散节点；如果每帧在 CPU 上 voxelize 再把边界传回 GPU，数据传输会吞掉模拟收益。因此边界生成、LBM 更新、粒子 advect 和渲染都尽量留在 GPU。

1. 从 packet distributions 到 D2Q9

Lattice Boltzmann Method 把空间划成 Cartesian lattice。每个节点保存多个分布 f_qi，离散时间步中先让它们在局部碰撞，再沿对应速度方向传播到邻居。宏观 density 是所有分布的和，velocity 是按方向加权后的动量。

D2Q9 有 9 个 packet distributions：中心项、上/下/左/右四个轴向项，以及四个对角项。3D 版本常用 D3Q19，每个节点有 19 个方向，内存主要就消耗在这些分布数组上。

rho = sum(q = 0 .. Q - 1, f[q])
u   = sum(q = 0 .. Q - 1, f[q] * e[q]) / rho

D2Q9 和 D3Q19 不是两种渲染材质，而是状态布局契约：每个方向都有自己的邻居偏移、权重和反向方向。选择模型时要同时考虑目标维度、边界精度、纹理数量和每步需要更新的带宽。

2. equilibrium、collision 和 streaming

collision 在每个节点局部完成。简化地说，先由 rho 和 u 计算 f_eq，再用 f_post = f + ω(f_eq − f) 松弛到平衡。ω 影响粘性和稳定性，过激的参数可能放大数值误差。

streaming 把 collision 后的 packet 按方向搬运。中心分布留在原节点，其他方向写入对应邻居；GPU 实现通常从 source textures 读取上一时间步，在 target textures 中形成下一步。

一个时间步的逻辑顺序是：读分布纹理 → 生成 density/velocity → 生成 equilibrium distributions → collision 与 streaming → 应用边界和 outflow → 交换分布纹理。每个阶段可以由 fragment program 对应纹理 texel 执行，3D 数据则拆成切片堆或 tiled 2D texture。

3. GPU 数据布局和内存路径

将所有相同速度向量的分布分到一组数组，能让每一张 2D texture 保持原 lattice 的布局。D2Q9 就是 9 张方向纹理或经过通道 packing 的若干纹理；每次更新时同一 lattice 坐标从各方向取样，计算一个节点的宏观量和新分布。

这里的性能边界通常是带宽，而不是单条公式。需要减少重复读取、合理复用 density/velocity 的中间纹理，并用双缓冲隔离时间步。原章的 GPU 实现还把计算、边界处理、粒子位置更新和渲染串在同一显存路径中，避免每步搬运大型状态数组。

bindTextures(sourceDistributions);
renderQuad(equilibriumAndCollision);
renderQuad(streamingAndBoundary);
swap(sourceDistributions, targetDistributions);

这段伪代码只表达资源关系：具体 shader 可以合并阶段，也可以为了边界分支单独处理。验证时应先对照 CPU LBM 的每个方向分布，再比较 density、velocity 和粒子轨迹，不能只看最终画面是否像流体。

4. 复杂边界：从连续几何到 lattice 节点

voxelization 把任意模型变成 Cartesian grid volume。静态障碍可以预计算；移动或变形物体则需要每次几何改变都重新生成边界节点。除了位置，还需要 wall velocity 和用于边界平面计算的系数。

简单的 slicing 会沿一个方向逐层移动 clip slab，render pass 数等于 slice 数。原章从三个正交方向进行 depth peeling，得到覆盖的体素集合，允许少量重复，因为重复不会损害体素边界的准确性。GPU 当前 fragment 无法 scatter 到任意输出位置，所以先把体素属性放入 vertex array，再通过 vertex program 将每个体素点放回正确位置。

5. boundary condition：bounce-back 与动态壁面

boundary condition 只对内部节点的 LBM 方程还不够。边界 link 两侧的两个分布方向相反且共线；简单的 bounce-back 把撞到固体的分布反射回去，曲壁和移动壁面则要加入交点位置与 u_w 修正。

对于动态边界，voxel array 还携带 voxel position、wall velocity 和 polygon plane coefficients。渲染边界顶点数组时，可以通过不同颜色通道和顶点平移，让一个 boundary node 收到它需要的相反方向分布。边界处理必须和 source/target 分布的时间步严格对齐，否则几何虽然更新了，packet 仍会按照旧墙面反射。

6. 可视化：让速度场可读

LBM 状态本身不容易直接观察，因此从入口注入带颜色的 particles，按当前 velocity field 更新位置。粒子越过出口或停在零速度区域时回收到入口；总粒子数保持固定，便于使用 vertex array 渲染。

先预测：把圆形障碍物换成条形障碍物，哪一侧更可能产生长尾涡流？把显示粒子关闭后，速度格点仍然应该变化，因为粒子只是可视化层，不是 LBM 的状态变量。

实验用小型 D2Q9 网格实际执行 equilibrium、collision、streaming 和边界反射；它刻意牺牲规模换取可读性。真实 3D D3Q19 系统的内存和带宽更大，但数据流、双缓冲和动态边界的契约相同。

小结

• LBM 用离散方向分布表示流体，D2Q9 和 D3Q19 决定每个节点的状态规模。
• collision 让分布靠近平衡，streaming 把 packet 沿 lattice link 传播到邻居。
• GPU 用方向纹理和 source/target 双缓冲把每个时间步写成规则数据流。
• 动态几何必须重新 voxelize，边界节点要携带位置、壁面速度和几何属性。
• 粒子 advect 是观察速度场的可视化层，不能替代分布状态或边界条件。