简介基于Lattice Boltzmann MethodLBM的顶盖驱动流C程序与何雅玲教授相关教材中的代码完全一致面向计算流体力学CFD初学者、LBM研究者以及数值模拟课程教学。通过该程序可在简单容器内模拟顶部边界移动引起的流体运动直观观察速度场、压力场随时间的演化从而深入理解顶盖驱动流这一经典验证算例以及LBM在处理复杂边界与并行计算方面的天然优势。资源包为RAR压缩格式内含1个cpp文件大小约1KB代码精炼适合直接阅读和运行便于对照理论推导与编程实现。已有351人浏览学习反映出其在LBM入门与教学场景中的实用价值。对刚接触CFD或LBM的读者而言这份程序不仅是一个可运行的数值实验更是连接教材理论、算法实现与流动现象分析的学习样本有助于快速掌握LBM建模思路与C实现要点。1. 顶盖驱动流是 LBM 的试金石题目里的 LDF 值得这样拆顶盖驱动流lid-driven cavity flow常缩写为 LDF大概是 LBM 学习者绕不开的第一个算例方腔四周封闭只有顶盖以恒定速度拖动流体腔内逐步形成主涡、角涡雷诺数一变涡结构跟着变。它像一块试金石边界处理差一点、参数换算错一点涡心位置就漂出基准解的可接受范围。这篇文章按国内外教材最常见的路径何雅玲老师的《格子Boltzmann方法的理论及应用》也沿用这条路线把方案讲清楚D2Q9 模型 BGK 碰撞 Zou-He/反弹边界从 Re 反推格子的松弛时间 tau跑通一个最小 Python 实现再和 Ghia 基准数据对比涡心与中心线剖面。适合刚把 LBM 公式看过一遍、想亲手复现基准解的人也适合想回过头把 tau、边界格式这些参数边界捋清楚的工程师。2. 从 Re 到 tau顶盖驱动流 LBM 的格子单位换算2.1 D2Q9 与 BGK 为什么是顶盖驱动流 LBM 的默认组合LBM 和传统 CFD 最直观的区别是不解压力泊松方程。密度是分布函数的零阶矩压力由状态方程 p rho * cs^2 直接给出压力场的地位从未知被求解降级成后处理读取。对顶盖驱动流这种封闭方腔压力本身不是考核重点速度场和涡结构才是所以 LBM 特别合适。碰撞算子层面BGK 单松弛格式因为代码量最少成为教学和预研的首选其演化方程是f_i(x c_i dt, t dt) f_i(x, t) - (f_i(x, t) - f_i^eq(x, t)) / tau这里的 i 取 0 到 8对应 D2Q9 的九个离散方向。零方向权重 4/9轴向四个方向权重 1/9对角线方向权重 1/36声速 cs sqrt(1/3) 是在 dt dx 1 的格子单位下算出来的。从 Chapman-Enskog 展开可以得到运动粘度和松弛时间的线性关系nu cs^2 * (tau - 0.5) * dt在格子单位下 dt 1cs^2 1/3于是 nu (tau - 0.5) / 3。这一条是后面所有参数换算的起点要确定格子粘性实际只需求出 tau要确 TAU先得把物理雷诺数翻译成格子单位。2.2 由雷诺数反推松弛时间 tau 的公式表顶盖驱动流的定义是 Re U_wall * L / nu。在格子单位里U_wall 就是顶盖速度 u_wall特征长度 L 取方腔两侧格点数之差 (N - 1)索引从 0 到 N-1两点之间的距离是 1。整理得到nu u_wall * (N - 1) / Retau 0.5 nu / cs^2 0.5 3 * u_wall * (N - 1) / Re每次换网格或换 Re都要重新用这一步反推 tau而不是随手把上一组参数复制过来。下面这组参数组合可以直接用作初始值也是我平时调试 LDF 的常用起点工况ReNu_wallnu格子单位tau低雷诺校验100650.050.0320.596常用基准10001290.050.00640.5192加密网格对照10002570.050.01280.5384高雷诺初探100002570.050.001280.50384提示tau 必须严格大于 0.5。表里最后一行的 tau 已经接近 0.503BGK 在高雷诺数下很容易出现棋盘式振荡这种情况通常需要换 MRT 或正则化 LBM这部分放到第 5 章展开。注意表里第二行和第三行的对比同一个 Re1000网格从 129 加密到 257 后格子粘性反而翻倍tau 离 0.5 更远了。原因在公式里——特征长度 L N-1 跟着网格走要保持 Re 不变nu 就得等比增大。刚接触 LBM 的人经常在这个位置犯迷糊网格越细不应该是粘性越小吗物理上是这样但前提是物理速度、物理长度和 dt/dx 的换算关系同步调整如果只换 N 而不重新定标时间步实际上模拟的是另一个物理算例。2.3 u_wall 的取值窗口Ma 数约束与何雅玲教材里的习惯做法u_wall 不是随便取的。顶盖驱动流 LBM 本质上是可压模型Ma u_wall / cs 太大时密度波动的非物理效应会污染速度场经验上限是 Ma 0.1~0.2。用 cs ≈ 0.577 折算u_wall 大致取 0.05 到 0.1 是安全的教学算例普遍落在 0.05。但 u_wall 又不能取得太小。固定 Re 和 N 之后u_wall 越小nu 越小tau 越靠近 0.5。tau 0.51 附近 BGK 的数值稳定性已经明显变差顶盖驱动流在 Re1000 时用 u_wall0.01、N129算出的 tau 是 0.50384跑几百步速度场就会开始抖。何雅玲老师的《格子Boltzmann方法的理论及应用》里顶盖驱动流算例的顶盖速度也基本落在这个量级原因不是教材保守而是这个窗口本身就是 LBM 的参数边界往下是稳定性下限往上是可压缩性上限。实际调试时我一般固定 u_wall0.05需要稳就往 0.1 方向提需要提高 Re 或减少可压缩影响就往 0.03 方向压。3. 用 Python Numpy 跑通最小顶盖驱动流 LBM3.1 流步前先给边界反弹与 Zou-He 的组合顺序LBM 循环里碰撞和流这两步本身没有条件分支边界条件全部体现在分布函数怎么改。顶盖驱动流 LBM 的最常见组合是底部和左右墙用反弹格式顶盖用 Zou-He 速度边界。反弹格式针对静止壁其思想是假设分布函数打到壁面后按原路弹回在代码里体现为三组方向互换左墙把从流体侧来的 f1、f5、f8 改成 f3、f6、f7右墙把 f3、f6、f7 改成 f1、f5、f8底墙把 f4、f7、f8 改成 f2、f5、f6。这三组赋值必须发生在流步之前因为赋值后这些分布函数要在接下来的流步里输送到相邻流体格点。顶盖是运动边界不能直接用反弹格式否则速度滑移没法处理。Zou-He 速度边界的做法是给定边界宏观速度后先由质量守恒关系反推出边界密度再假设未知方向分布的非平衡部分与对向已知分布互为相反数解出几个缺失方向。对顶盖而言法线指向 y 正方向未知方向是 f4、f7、f8已知方向 f0、f1、f2、f3、f5、f6。zou-he 给出的公式在 v0 时简化为f4 f2f7 f6 - 0.5*(f1-f3) - rhou_wall/6f8 f5 0.5(f1-f3) rho*u_wall/6。提示Zou-He 公式里的三项切向修正0.5*(f1-f3) 和 rho*u_wall/6不能省。删掉它们退化成反弹加速度修正的一阶格式中心线剖面的涡心偏差会明显变大。3.2 顶盖驱动流 LBM 的完整可运行代码下面的代码是 Numpy 实现N65、Re100 时几秒内可以跑完把 N 改成 129、Re 改成 1000 也能在普通笔记本上完成只是要多等几分钟。代码里注释标明了哪些行是 LBM 理论公式的工程落地点。import numpy as np N 65 Re 100.0 u_wall 0.05 cs2 1.0 / 3.0 tau 0.5 3.0 * u_wall * (N - 1) / Re max_step 8000 tol 1e-10 # D2Q9 九个方向的位移 X np.array([0, 1, 0, -1, 0, 1, -1, -1, 1]) Y np.array([0, 0, 1, 0, -1, 1, 1, -1, -1]) W np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36]) f np.tile(W, (N, N, 1)).astype(float) ux np.zeros((N, N)) uy np.zeros((N, N)) def feq(rho, ux, uy): cu X[None, None, :] * ux[..., None] Y[None, None, :] * uy[..., None] u2 ux * ux uy * uy return W[None, None, :] * rho[..., None] * ( 1.0 3.0 * cu 4.5 * cu * cu - 1.5 * u2[..., None]) for step in range(max_step): rho f.sum(axis-1) ux_new np.einsum(ijk,k-ij, f, X) / rho uy_new np.einsum(ijk,k-ij, f, Y) / rho # 残差控制取速度场的最大变化量 if step % 50 0: res np.max(np.abs(ux_new - ux)) np.max(np.abs(uy_new - uy)) if res tol: print(converged at step, step) break ux, uy ux_new, uy_new # 碰撞BGK 松弛到平衡态 f (feq(rho, ux, uy) - f) / tau # 静止壁反弹左墙、右墙、底墙 f[:, 0, 3] f[:, 0, 1] f[:, 0, 6] f[:, 0, 5] f[:, 0, 7] f[:, 0, 8] f[:, -1, 1] f[:, -1, 3] f[:, -1, 5] f[:, -1, 6] f[:, -1, 8] f[:, -1, 7] f[0, :, 2] f[0, :, 4] f[0, :, 5] f[0, :, 7] f[0, :, 6] f[0, :, 8] # 顶盖 Zou-He 速度边界v0、uu_wall rho_wall (f[-1, :, 0] f[-1, :, 1] f[-1, :, 3] 2.0 * (f[-1, :, 2] f[-1, :, 5] f[-1, :, 6])) f[-1, :, 4] f[-1, :, 2] f[-1, :, 7] (f[-1, :, 6] - 0.5 * (f[-1, :, 1] - f[-1, :, 3]) - rho_wall * u_wall / 6.0) f[-1, :, 8] (f[-1, :, 5] 0.5 * (f[-1, :, 1] - f[-1, :, 3]) rho_wall * u_wall / 6.0) # 流步np.roll 会周期卷绕四行边界已被重新赋值等效封闭方腔 f np.roll(f, shift(-X, -Y), axis(0, 1))速度残差用两次宏观速度的最大绝对值差来衡量比用分布函数差更直观。碰撞步发生在边界赋值之前边界节点的分布函数被碰撞更新后又马上覆盖所以边界节点的碰撞计算是无效功保留它纯粹是为了避免把内部节点单独挑出来代码更短。np.roll 的周期卷绕只影响四个边界行而边界行在这步之前已经被反弹和 Zou-He 完整赋值卷绕进来的分布不会再进入内部区域因此这种写法和显式循环的流步结果一致。3.3 边界方式怎么选反弹、Zou-He、非平衡外推的差别这个算例里我选了三面反弹 顶盖 Zou-He但它不是唯一方案。反弹格式实现最简单、对复杂几何最友好代价是静止壁精度只有一阶顶盖驱动流中底壁和侧壁的局部涡量误差会略微偏大。Zou-He 是速度边界二阶精度但有明显的方向性每个边界都要按法向换一套公式四个角点涉及多重边界的赋值冲突我的处理顺序是先反弹三面墙、最后用顶盖整行覆盖左右上角结果不会进入主流区。非平衡外推Guo 格式实现最对称先给边界节点赋已知宏观量 rho 和 u然后 f f_eq(rho, u) (f_相邻 - f_eq_相邻)把非平衡部分从相邻流体节点抄过来顶盖驱动流 LBM 用它也很常见何雅玲老师的教材里对这类外推思想有专门说明。要在三者之间切换基本只动边界处理这一段。非平衡外推不需要为每个方向写公式适合做参数扫描Zou-He 在顶盖上的基准校验结果最稳适合做主算例。三个版本的顶盖驱动流 LBM 我都跑过低 Re 下涡心坐标差别小于一个格点到 Re10000 才会拉开差距——到那个区间 BGK 本身也已经接近极限了。4. 用 Ghia 基准数据校验 LDF涡心与中心线速度剖面4.1 涡量与流函数的格子单位算法顶盖驱动流 LBM 跑完后的直接输出是速度场要跟文献基准对比需要从速度场里提取涡心和中心线剖面。涡量计算非常简单用中心差分即可dx dy 1Nf ux.shape[0] omega np.zeros_like(ux) omega[1:-1, 1:-1] (uy[1:-1, 2:] - uy[1:-1, :-2]) / 2.0 \ - (ux[2:, 1:-1] - ux[:-2, 1:-1]) / 2.0流函数 psi 用来定位涡心。psi 满足 laplace 型方程 ∇²psi -omega四边均为 0封闭方腔无穿透。用 SOR 迭代解这个方程的标准做法是边界值为 0内部点逐次更新 psi[i,j] 0.25*(四邻点之和 omega[i,j])迭代到残差小于 1e-9。主涡是顺时针涡psi 在涡心处取全局极小值用 np.argmin 找全局最小值位置就可以得到格点坐标。4.2 中心线速度剖面怎么跟基准数据对齐Ghia 基准数据是 1982 年发表在 Journal of Computational Physics 上的有限差分解是顶盖驱动流 LBM 最常用的对照表。提取剖面的代码只有四行j0 (N - 1) // 2 u_center ux[:, j0] # 垂直中线 x0.5 上的水平速度 v_center uy[j0, :] # 水平中线 y0.5 上的垂直速度把 u_center 按 y 坐标画线和 Ghia 表里 Re100 的 u 速度数据放在同一张图里v_center 按 x 坐标画线和 Ghia 的 v 速度数据对比。注意两个剖面的横轴坐标是 y/L 或 x/L 的无量纲位置在格子单位下就是索引除以 (N-1)。网格 65 时剖面是 65 个点直接和 Ghia 表里面的 33 个参考点对比网格 129 时插值对齐即可。主涡心的参考位置也是硬指标我整理了常用的三组基准值Re主涡心 x/L主涡心 y/L参考来源1000.61720.7344Ghia et al., 19824000.55470.6055Ghia et al., 198210000.53130.5625Ghia et al., 1982LBM 算出的涡心如果落在参考位置周围一个格点以内边界和参数就基本没问题偏差超过两三个格点优先怀疑 tau 算错、顶盖速度没有真正施加或者还没收敛就停了。4.3 收敛判据残差小于多少才能停顶盖驱动流是稳态问题但 LBM 迭代是显式时间推进密度波在方腔里的衰减速度决定了收敛步伐。只跑固定步数不可靠Re100 时 3000 步接近稳态Re1000 时 10000 步经常还有 1e-7 量级的缓慢漂移。用速度场的 L-infinity 残差判断时阈值取 1e-9 到 1e-10 比较稳阈值放宽到 1e-6涡心坐标看起来不变但中心线剖面在和 Ghia 数据对比时会差出一个马赫数量级的系统偏差这个偏差容易被误判成边界格式问题。判断是否真正收敛的辅助手段是监控主涡心坐标连续 500 步涡心不动比残差更可靠。残差的格子单位含义要搞清楚格子速度 1e-9 对应物理速度约 1e-9 * cs在剖面曲线上是看不见的所以这个阈值并不是小题大做。5. 三个顶盖驱动流 LBM 的实战技巧参数窗口、涡心追踪与角点噪声5.1 同一 Re 下 u_wall 的取舍tau 的公式把三个量绑在一起固定 Re 时u_wall 和 N 共同决定 tau。想用加大 N 来换更准的涡心又保持 u_wall0.05tau 反而变大可压缩误差没有下降网格只是让壁面边界层多几个点。想压制可压缩误差而把 u_wall 降到 0.02tau 又会滑向危险区。我的习惯是Re 1000 用 u_wall 0.05、tau 落在 0.52~0.60Re 超过 10000 时不再依赖 BGK改用 MRT 或者正则化碰撞算子后者在 tau 接近 0.51 时仍能维持稳定。用 N 做网格无关性验证时同时把 u_wall 按比例调整让 tau 保持在 0.52 到 0.60 窗口内而不是机械地固定 Re 和 u_wall 只改 N。5.2 用涡心坐标而不是残差判断稳态残差小不代表涡结构已经到达平衡顶盖驱动流在中等雷诺数下存在缓慢的二次涡迁移速度场残差 1e-8 时涡心每 1000 步可能还能移动小半个格点。我在后处理阶段的做法是保存多个快照分别计算流函数追踪主涡心psi, omega compute_stream_function(ux, uy) # 自己整理的 SOR 求解函数 i_c, j_c np.unravel_index(np.argmin(psi), psi.shape)把每 500 步的 (i_c, j_c) 序列列出来当连续两次快照的涡心坐标差值不超过 1 个格点时再认定收敛此时再去画中心线剖面和 Ghia 数据对比。这套逻辑也可以直接并进主循环做自适应停机只是每 500 步解一次泊松方程会拖慢纯 Python 代码我一般只在模拟结束后对保存的速度场做离线分析。5.3 角点噪声的处理顶盖两个上角点同时是移动壁和静止壁的交界角点的速度函数不连续Zou-He 赋值后局部密度会出现毫厘级的振荡。这个振荡顺着顶盖向下游传播时会被 BGK 耗散衰减距离大概几个格子主流区的主涡基本不受影响但在画涡量等值线时角点噪声会掩盖很小的角涡误导判断。处理办法是后处理时把最外一层网格剔除用内层区域画等值线如果想保角涡可以把顶盖角点的速度过渡写成线性 ramp让 u_wall 在边上两个格点内线性上升代价是主涡心会偏移大约零点几个格点对比基准解前要记得这一步是人为修改过的边界条件。这个取舍没有标准答案以你的对比目标为准要复现 Ghia 涡心就保留阶跃速度要观察角落涡结构就用 ramp 或直接避开最外层。最后再检查一遍 Re 定义里的 L 是 (N-1)一旦这里多写或少写 1整套参数都会偏。本文还有配套的精品资源点击获取 顶盖驱动流LBM模拟：参数换算、边界处理与Python实现
简介基于Lattice Boltzmann MethodLBM的顶盖驱动流C程序与何雅玲教授相关教材中的代码完全一致面向计算流体力学CFD初学者、LBM研究者以及数值模拟课程教学。通过该程序可在简单容器内模拟顶部边界移动引起的流体运动直观观察速度场、压力场随时间的演化从而深入理解顶盖驱动流这一经典验证算例以及LBM在处理复杂边界与并行计算方面的天然优势。资源包为RAR压缩格式内含1个cpp文件大小约1KB代码精炼适合直接阅读和运行便于对照理论推导与编程实现。已有351人浏览学习反映出其在LBM入门与教学场景中的实用价值。对刚接触CFD或LBM的读者而言这份程序不仅是一个可运行的数值实验更是连接教材理论、算法实现与流动现象分析的学习样本有助于快速掌握LBM建模思路与C实现要点。1. 顶盖驱动流是 LBM 的试金石题目里的 LDF 值得这样拆顶盖驱动流lid-driven cavity flow常缩写为 LDF大概是 LBM 学习者绕不开的第一个算例方腔四周封闭只有顶盖以恒定速度拖动流体腔内逐步形成主涡、角涡雷诺数一变涡结构跟着变。它像一块试金石边界处理差一点、参数换算错一点涡心位置就漂出基准解的可接受范围。这篇文章按国内外教材最常见的路径何雅玲老师的《格子Boltzmann方法的理论及应用》也沿用这条路线把方案讲清楚D2Q9 模型 BGK 碰撞 Zou-He/反弹边界从 Re 反推格子的松弛时间 tau跑通一个最小 Python 实现再和 Ghia 基准数据对比涡心与中心线剖面。适合刚把 LBM 公式看过一遍、想亲手复现基准解的人也适合想回过头把 tau、边界格式这些参数边界捋清楚的工程师。2. 从 Re 到 tau顶盖驱动流 LBM 的格子单位换算2.1 D2Q9 与 BGK 为什么是顶盖驱动流 LBM 的默认组合LBM 和传统 CFD 最直观的区别是不解压力泊松方程。密度是分布函数的零阶矩压力由状态方程 p rho * cs^2 直接给出压力场的地位从未知被求解降级成后处理读取。对顶盖驱动流这种封闭方腔压力本身不是考核重点速度场和涡结构才是所以 LBM 特别合适。碰撞算子层面BGK 单松弛格式因为代码量最少成为教学和预研的首选其演化方程是f_i(x c_i dt, t dt) f_i(x, t) - (f_i(x, t) - f_i^eq(x, t)) / tau这里的 i 取 0 到 8对应 D2Q9 的九个离散方向。零方向权重 4/9轴向四个方向权重 1/9对角线方向权重 1/36声速 cs sqrt(1/3) 是在 dt dx 1 的格子单位下算出来的。从 Chapman-Enskog 展开可以得到运动粘度和松弛时间的线性关系nu cs^2 * (tau - 0.5) * dt在格子单位下 dt 1cs^2 1/3于是 nu (tau - 0.5) / 3。这一条是后面所有参数换算的起点要确定格子粘性实际只需求出 tau要确 TAU先得把物理雷诺数翻译成格子单位。2.2 由雷诺数反推松弛时间 tau 的公式表顶盖驱动流的定义是 Re U_wall * L / nu。在格子单位里U_wall 就是顶盖速度 u_wall特征长度 L 取方腔两侧格点数之差 (N - 1)索引从 0 到 N-1两点之间的距离是 1。整理得到nu u_wall * (N - 1) / Retau 0.5 nu / cs^2 0.5 3 * u_wall * (N - 1) / Re每次换网格或换 Re都要重新用这一步反推 tau而不是随手把上一组参数复制过来。下面这组参数组合可以直接用作初始值也是我平时调试 LDF 的常用起点工况ReNu_wallnu格子单位tau低雷诺校验100650.050.0320.596常用基准10001290.050.00640.5192加密网格对照10002570.050.01280.5384高雷诺初探100002570.050.001280.50384提示tau 必须严格大于 0.5。表里最后一行的 tau 已经接近 0.503BGK 在高雷诺数下很容易出现棋盘式振荡这种情况通常需要换 MRT 或正则化 LBM这部分放到第 5 章展开。注意表里第二行和第三行的对比同一个 Re1000网格从 129 加密到 257 后格子粘性反而翻倍tau 离 0.5 更远了。原因在公式里——特征长度 L N-1 跟着网格走要保持 Re 不变nu 就得等比增大。刚接触 LBM 的人经常在这个位置犯迷糊网格越细不应该是粘性越小吗物理上是这样但前提是物理速度、物理长度和 dt/dx 的换算关系同步调整如果只换 N 而不重新定标时间步实际上模拟的是另一个物理算例。2.3 u_wall 的取值窗口Ma 数约束与何雅玲教材里的习惯做法u_wall 不是随便取的。顶盖驱动流 LBM 本质上是可压模型Ma u_wall / cs 太大时密度波动的非物理效应会污染速度场经验上限是 Ma 0.1~0.2。用 cs ≈ 0.577 折算u_wall 大致取 0.05 到 0.1 是安全的教学算例普遍落在 0.05。但 u_wall 又不能取得太小。固定 Re 和 N 之后u_wall 越小nu 越小tau 越靠近 0.5。tau 0.51 附近 BGK 的数值稳定性已经明显变差顶盖驱动流在 Re1000 时用 u_wall0.01、N129算出的 tau 是 0.50384跑几百步速度场就会开始抖。何雅玲老师的《格子Boltzmann方法的理论及应用》里顶盖驱动流算例的顶盖速度也基本落在这个量级原因不是教材保守而是这个窗口本身就是 LBM 的参数边界往下是稳定性下限往上是可压缩性上限。实际调试时我一般固定 u_wall0.05需要稳就往 0.1 方向提需要提高 Re 或减少可压缩影响就往 0.03 方向压。3. 用 Python Numpy 跑通最小顶盖驱动流 LBM3.1 流步前先给边界反弹与 Zou-He 的组合顺序LBM 循环里碰撞和流这两步本身没有条件分支边界条件全部体现在分布函数怎么改。顶盖驱动流 LBM 的最常见组合是底部和左右墙用反弹格式顶盖用 Zou-He 速度边界。反弹格式针对静止壁其思想是假设分布函数打到壁面后按原路弹回在代码里体现为三组方向互换左墙把从流体侧来的 f1、f5、f8 改成 f3、f6、f7右墙把 f3、f6、f7 改成 f1、f5、f8底墙把 f4、f7、f8 改成 f2、f5、f6。这三组赋值必须发生在流步之前因为赋值后这些分布函数要在接下来的流步里输送到相邻流体格点。顶盖是运动边界不能直接用反弹格式否则速度滑移没法处理。Zou-He 速度边界的做法是给定边界宏观速度后先由质量守恒关系反推出边界密度再假设未知方向分布的非平衡部分与对向已知分布互为相反数解出几个缺失方向。对顶盖而言法线指向 y 正方向未知方向是 f4、f7、f8已知方向 f0、f1、f2、f3、f5、f6。zou-he 给出的公式在 v0 时简化为f4 f2f7 f6 - 0.5*(f1-f3) - rhou_wall/6f8 f5 0.5(f1-f3) rho*u_wall/6。提示Zou-He 公式里的三项切向修正0.5*(f1-f3) 和 rho*u_wall/6不能省。删掉它们退化成反弹加速度修正的一阶格式中心线剖面的涡心偏差会明显变大。3.2 顶盖驱动流 LBM 的完整可运行代码下面的代码是 Numpy 实现N65、Re100 时几秒内可以跑完把 N 改成 129、Re 改成 1000 也能在普通笔记本上完成只是要多等几分钟。代码里注释标明了哪些行是 LBM 理论公式的工程落地点。import numpy as np N 65 Re 100.0 u_wall 0.05 cs2 1.0 / 3.0 tau 0.5 3.0 * u_wall * (N - 1) / Re max_step 8000 tol 1e-10 # D2Q9 九个方向的位移 X np.array([0, 1, 0, -1, 0, 1, -1, -1, 1]) Y np.array([0, 0, 1, 0, -1, 1, 1, -1, -1]) W np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36]) f np.tile(W, (N, N, 1)).astype(float) ux np.zeros((N, N)) uy np.zeros((N, N)) def feq(rho, ux, uy): cu X[None, None, :] * ux[..., None] Y[None, None, :] * uy[..., None] u2 ux * ux uy * uy return W[None, None, :] * rho[..., None] * ( 1.0 3.0 * cu 4.5 * cu * cu - 1.5 * u2[..., None]) for step in range(max_step): rho f.sum(axis-1) ux_new np.einsum(ijk,k-ij, f, X) / rho uy_new np.einsum(ijk,k-ij, f, Y) / rho # 残差控制取速度场的最大变化量 if step % 50 0: res np.max(np.abs(ux_new - ux)) np.max(np.abs(uy_new - uy)) if res tol: print(converged at step, step) break ux, uy ux_new, uy_new # 碰撞BGK 松弛到平衡态 f (feq(rho, ux, uy) - f) / tau # 静止壁反弹左墙、右墙、底墙 f[:, 0, 3] f[:, 0, 1] f[:, 0, 6] f[:, 0, 5] f[:, 0, 7] f[:, 0, 8] f[:, -1, 1] f[:, -1, 3] f[:, -1, 5] f[:, -1, 6] f[:, -1, 8] f[:, -1, 7] f[0, :, 2] f[0, :, 4] f[0, :, 5] f[0, :, 7] f[0, :, 6] f[0, :, 8] # 顶盖 Zou-He 速度边界v0、uu_wall rho_wall (f[-1, :, 0] f[-1, :, 1] f[-1, :, 3] 2.0 * (f[-1, :, 2] f[-1, :, 5] f[-1, :, 6])) f[-1, :, 4] f[-1, :, 2] f[-1, :, 7] (f[-1, :, 6] - 0.5 * (f[-1, :, 1] - f[-1, :, 3]) - rho_wall * u_wall / 6.0) f[-1, :, 8] (f[-1, :, 5] 0.5 * (f[-1, :, 1] - f[-1, :, 3]) rho_wall * u_wall / 6.0) # 流步np.roll 会周期卷绕四行边界已被重新赋值等效封闭方腔 f np.roll(f, shift(-X, -Y), axis(0, 1))速度残差用两次宏观速度的最大绝对值差来衡量比用分布函数差更直观。碰撞步发生在边界赋值之前边界节点的分布函数被碰撞更新后又马上覆盖所以边界节点的碰撞计算是无效功保留它纯粹是为了避免把内部节点单独挑出来代码更短。np.roll 的周期卷绕只影响四个边界行而边界行在这步之前已经被反弹和 Zou-He 完整赋值卷绕进来的分布不会再进入内部区域因此这种写法和显式循环的流步结果一致。3.3 边界方式怎么选反弹、Zou-He、非平衡外推的差别这个算例里我选了三面反弹 顶盖 Zou-He但它不是唯一方案。反弹格式实现最简单、对复杂几何最友好代价是静止壁精度只有一阶顶盖驱动流中底壁和侧壁的局部涡量误差会略微偏大。Zou-He 是速度边界二阶精度但有明显的方向性每个边界都要按法向换一套公式四个角点涉及多重边界的赋值冲突我的处理顺序是先反弹三面墙、最后用顶盖整行覆盖左右上角结果不会进入主流区。非平衡外推Guo 格式实现最对称先给边界节点赋已知宏观量 rho 和 u然后 f f_eq(rho, u) (f_相邻 - f_eq_相邻)把非平衡部分从相邻流体节点抄过来顶盖驱动流 LBM 用它也很常见何雅玲老师的教材里对这类外推思想有专门说明。要在三者之间切换基本只动边界处理这一段。非平衡外推不需要为每个方向写公式适合做参数扫描Zou-He 在顶盖上的基准校验结果最稳适合做主算例。三个版本的顶盖驱动流 LBM 我都跑过低 Re 下涡心坐标差别小于一个格点到 Re10000 才会拉开差距——到那个区间 BGK 本身也已经接近极限了。4. 用 Ghia 基准数据校验 LDF涡心与中心线速度剖面4.1 涡量与流函数的格子单位算法顶盖驱动流 LBM 跑完后的直接输出是速度场要跟文献基准对比需要从速度场里提取涡心和中心线剖面。涡量计算非常简单用中心差分即可dx dy 1Nf ux.shape[0] omega np.zeros_like(ux) omega[1:-1, 1:-1] (uy[1:-1, 2:] - uy[1:-1, :-2]) / 2.0 \ - (ux[2:, 1:-1] - ux[:-2, 1:-1]) / 2.0流函数 psi 用来定位涡心。psi 满足 laplace 型方程 ∇²psi -omega四边均为 0封闭方腔无穿透。用 SOR 迭代解这个方程的标准做法是边界值为 0内部点逐次更新 psi[i,j] 0.25*(四邻点之和 omega[i,j])迭代到残差小于 1e-9。主涡是顺时针涡psi 在涡心处取全局极小值用 np.argmin 找全局最小值位置就可以得到格点坐标。4.2 中心线速度剖面怎么跟基准数据对齐Ghia 基准数据是 1982 年发表在 Journal of Computational Physics 上的有限差分解是顶盖驱动流 LBM 最常用的对照表。提取剖面的代码只有四行j0 (N - 1) // 2 u_center ux[:, j0] # 垂直中线 x0.5 上的水平速度 v_center uy[j0, :] # 水平中线 y0.5 上的垂直速度把 u_center 按 y 坐标画线和 Ghia 表里 Re100 的 u 速度数据放在同一张图里v_center 按 x 坐标画线和 Ghia 的 v 速度数据对比。注意两个剖面的横轴坐标是 y/L 或 x/L 的无量纲位置在格子单位下就是索引除以 (N-1)。网格 65 时剖面是 65 个点直接和 Ghia 表里面的 33 个参考点对比网格 129 时插值对齐即可。主涡心的参考位置也是硬指标我整理了常用的三组基准值Re主涡心 x/L主涡心 y/L参考来源1000.61720.7344Ghia et al., 19824000.55470.6055Ghia et al., 198210000.53130.5625Ghia et al., 1982LBM 算出的涡心如果落在参考位置周围一个格点以内边界和参数就基本没问题偏差超过两三个格点优先怀疑 tau 算错、顶盖速度没有真正施加或者还没收敛就停了。4.3 收敛判据残差小于多少才能停顶盖驱动流是稳态问题但 LBM 迭代是显式时间推进密度波在方腔里的衰减速度决定了收敛步伐。只跑固定步数不可靠Re100 时 3000 步接近稳态Re1000 时 10000 步经常还有 1e-7 量级的缓慢漂移。用速度场的 L-infinity 残差判断时阈值取 1e-9 到 1e-10 比较稳阈值放宽到 1e-6涡心坐标看起来不变但中心线剖面在和 Ghia 数据对比时会差出一个马赫数量级的系统偏差这个偏差容易被误判成边界格式问题。判断是否真正收敛的辅助手段是监控主涡心坐标连续 500 步涡心不动比残差更可靠。残差的格子单位含义要搞清楚格子速度 1e-9 对应物理速度约 1e-9 * cs在剖面曲线上是看不见的所以这个阈值并不是小题大做。5. 三个顶盖驱动流 LBM 的实战技巧参数窗口、涡心追踪与角点噪声5.1 同一 Re 下 u_wall 的取舍tau 的公式把三个量绑在一起固定 Re 时u_wall 和 N 共同决定 tau。想用加大 N 来换更准的涡心又保持 u_wall0.05tau 反而变大可压缩误差没有下降网格只是让壁面边界层多几个点。想压制可压缩误差而把 u_wall 降到 0.02tau 又会滑向危险区。我的习惯是Re 1000 用 u_wall 0.05、tau 落在 0.52~0.60Re 超过 10000 时不再依赖 BGK改用 MRT 或者正则化碰撞算子后者在 tau 接近 0.51 时仍能维持稳定。用 N 做网格无关性验证时同时把 u_wall 按比例调整让 tau 保持在 0.52 到 0.60 窗口内而不是机械地固定 Re 和 u_wall 只改 N。5.2 用涡心坐标而不是残差判断稳态残差小不代表涡结构已经到达平衡顶盖驱动流在中等雷诺数下存在缓慢的二次涡迁移速度场残差 1e-8 时涡心每 1000 步可能还能移动小半个格点。我在后处理阶段的做法是保存多个快照分别计算流函数追踪主涡心psi, omega compute_stream_function(ux, uy) # 自己整理的 SOR 求解函数 i_c, j_c np.unravel_index(np.argmin(psi), psi.shape)把每 500 步的 (i_c, j_c) 序列列出来当连续两次快照的涡心坐标差值不超过 1 个格点时再认定收敛此时再去画中心线剖面和 Ghia 数据对比。这套逻辑也可以直接并进主循环做自适应停机只是每 500 步解一次泊松方程会拖慢纯 Python 代码我一般只在模拟结束后对保存的速度场做离线分析。5.3 角点噪声的处理顶盖两个上角点同时是移动壁和静止壁的交界角点的速度函数不连续Zou-He 赋值后局部密度会出现毫厘级的振荡。这个振荡顺着顶盖向下游传播时会被 BGK 耗散衰减距离大概几个格子主流区的主涡基本不受影响但在画涡量等值线时角点噪声会掩盖很小的角涡误导判断。处理办法是后处理时把最外一层网格剔除用内层区域画等值线如果想保角涡可以把顶盖角点的速度过渡写成线性 ramp让 u_wall 在边上两个格点内线性上升代价是主涡心会偏移大约零点几个格点对比基准解前要记得这一步是人为修改过的边界条件。这个取舍没有标准答案以你的对比目标为准要复现 Ghia 涡心就保留阶跃速度要观察角落涡结构就用 ramp 或直接避开最外层。最后再检查一遍 Re 定义里的 L 是 (N-1)一旦这里多写或少写 1整套参数都会偏。本文还有配套的精品资源点击获取