简介面向计算流体力学初学者与格子玻尔兹曼方法研究者这份资源提供基于格子玻尔兹曼方法的顶盖驱动流C模拟程序对应何雅玲教授著作中的标准代码实现。顶盖驱动流是计算流体动力学中检验数值算法的经典问题格子玻尔兹曼方法以离散玻尔兹曼方程为根基能较自然地处理复杂边界与并行计算因此被广泛用于微流动、多相流等场景。程序通过顶盖移动驱动方腔流体运动可清晰呈现速度场、压力场等流动特性既适合用于验证算法与教学演示也可作为二次开发的基础框架。压缩包为rar格式仅含1个cpp源文件整体大小约1KB代码量精简但结构完整便于逐行阅读和修改有助于将理论落实到实际编程。该资源已有351人学习如果希望快速理解顶盖驱动流这一经典案例或需要一份可直接运行、对照教材学习的参考代码这份资料能提供直观支持。同时简洁的单一文件结构也降低了环境配置成本适合作为课程设计或科研入门的起步模板。1. LDF与顶盖驱动流LBM一个被低估却决定精度的核心数组在格子Boltzmann方法LBM里LDF通常指lattice distribution function也就是格子分布函数。顶盖驱动流作为CFD最经典的基准算例几乎每个学LBM的人都跑过但真正把LDF从初始化、碰撞、迁移到边界覆盖这一整条链路理清的人并不多。何雅玲等作者在相关专著中把顶盖驱动流当作标准验证算例不是因为网格简单而是因为四个边界的运动状态完全不同上盖以恒定速度拖动其余三壁静止四个角点还存在速度间断。LDF在这些位置上的行为会直接暴露索引写错、边界方向反了、碰撞迁移顺序不合适等问题。下面从BGK方程和D2Q9模型入手把LDF在顶盖驱动流LBM里的组织方式、边界修正和结果验证说透。适合已经跑通过均匀流LBM、正准备处理复杂边界的工程师也适合想搞清非平衡外推和反弹边界差异的读者。2. 顶盖驱动流LBM的LDF更新方程与D2Q9数组排布2.1 BGK方程中LDF的三步演进顶盖驱动流采用不可压BGK近似LDF的更新方程可以写成f_i(x c_i δt, t δt) - f_i(x, t) - (f_i - f_i^eq) / τ这里的 f_i(x,t) 就是第 i 个离散速度方向上的LDFτ 是无量纲松弛时间c_i 是D2Q9离散速度。顶盖驱动流没有外力项碰撞、迁移、边界修正三步循环即可。宏观密度和宏观速度由LDF的零阶矩和一阶矩恢复ρ Σ f_iρu Σ f_i c_iLDF的核心价值在于它是碰撞和迁移的真正载体。边界条件、压力修正、动量交换都要通过LDF来写而不是直接操作速度。顶盖驱动流之所以适合练手是因为它能把LDF在四个方向上的行为都覆盖到。2.1.1 D2Q9速度索引如何对应物理方向先固定一组常用的D2Q9索引顺序cx np.array([0, 1, 0, -1, 0, 1, -1, -1, 1]) cy np.array([0, 0, 1, 0, -1, 1, 1, -1, -1]) omega np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36])索引0静止粒子1到4东、北、西、南四个正方向5到8东北、西北、西南、东南四个对角方向。omega[i]是对应的权重系数。LBM中几乎所有边界处理都会直接用到这些索引所以最好把它作为常量写死不要在迁移函数里临时计算。2.2 LDF在顶盖驱动流中的轴顺序和内存布局常见做法是让数组形状为f[ny, nx, 9]第一维是 y第二维是 x第三维是离散速度索引。这样写有几个好处宏观量计算可以直接对axis2求和迁移时用切片操作也容易和物理方向对齐。对象推荐写法说明离散分布函数f[ny, nx, 9]最后一位是速度索引宏观密度rho[ny, nx]对LDF求和宏观速度ux[ny, nx]、uy[ny, nx]矩恢复得到计算域x 0 为左壁y 0 为底壁顶盖位于 y ny-1初始化时常见做法是直接填充平衡态分布def init_field(nx, ny, rho01.0): f np.zeros((ny, nx, 9), dtypenp.float64) rho np.full((ny, nx), rho0, dtypenp.float64) ux np.zeros((ny, nx), dtypenp.float64) uy np.zeros((ny, nx), dtypenp.float64) for i in range(9): f[..., i] omega[i] * rho0 return f, rho, ux, uy这里直接用omega[i] * rho0作为初始LDF是因为静止流场的平衡态分布就是这个形式。如果初值偏离平衡态太远前几百步会激发出明显的高频压力波顶盖驱动流这类低速算例虽然不会崩但收敛速度会变慢。2.3 碰撞和迁移的顺序选择推荐的做法是先碰撞后迁移从当前LDF计算宏观量用BGK碰撞算子更新LDF把碰撞后的LDF沿离散速度方向迁移对边界上缺失的LDF做覆盖。def collide(f, tau): rho f.sum(axis2) rho[rho 1e-12] 1e-12 ux (f[..., 1] - f[..., 3] f[..., 5] - f[..., 6] - f[..., 7] f[..., 8]) / rho uy (f[..., 2] - f[..., 4] f[..., 5] f[..., 6] - f[..., 7] - f[..., 8]) / rho feq np.zeros_like(f) for i in range(9): cu 3.0 * (cx[i] * ux cy[i] * uy) feq[..., i] omega[i] * rho * ( 1.0 cu 0.5 * cu * cu - 1.5 * (ux * ux uy * uy) ) f - (f - feq) / tau return rho, ux, uyrho先做一个小幅下限截断是为了避免除以零。顶盖驱动流的密度场不会出现真空但这个保护对初学者调试其他算例时很有用。2.3.1 迁移函数不要用np.roll直接回绕很多示例代码用np.roll做迁移但np.roll会把模型另一边的LDF卷进边界顶盖驱动流是封闭方腔必须逐方向切片迁移def stream(f): f_new np.empty_like(f) f_new[:, :, 0] f[:, :, 0] f_new[:, 1:, 1] f[:, :-1, 1] f_new[1:, :, 2] f[:-1, :, 2] f_new[:, :-1, 3] f[:, 1:, 3] f_new[:-1, :, 4] f[1:, :, 4] f_new[1:, 1:, 5] f[:-1, :-1, 5] f_new[1:, :-1, 6] f[:-1, 1:, 6] f_new[:-1, :-1, 7] f[1:, 1:, 7] f_new[:-1, 1:, 8] f[1:, :-1, 8] return f_new以方向1东向为例f[..., 1][x]的新值来自f[..., 1][x-1]因此用f_new[:, 1:, 1] f[:, :-1, 1]。边界节点 x0 的方向1不会在这一步被填上它要留给后面的边界函数处理。先想清楚哪个方向来自哪个邻居再写切片否则一跑读取基准数据时很快会发现问题。3. 用何雅玲专著中常用的非平衡外推给顶盖驱动流LBM做边界3.1 顶盖驱动流边界条件如何落到LDF上顶盖驱动流的上壁以速度 U_lid 向右运动左、右、底壁速度为零。在LBM中宏观速度边界需要转化为“哪些LDF分量是未知的”。用标准反弹只能得到一阶精度而且角点处容易产生不对称工程代码常用非平衡外推。非平衡外推的核心思路是边界上未知的LDF等于该位置的平衡态部分加上内部邻近格点上LDF的非平衡部分。写成常见形式f_i(x_b) f_i^eq(ρ_b, u_b) f_i(x_f) - f_i^eq(ρ_f, u_f)其中 x_b 是边界格点x_f 是沿内法线方向距离一个格距的流体格点。这个格式不像反弹边界那样需要事先猜一个虚拟格点对弯曲边界和压力驱动边界也更容易扩展也是何雅玲等作者编写的LBM教材中重点介绍的方向。3.2 四个边界上到底哪几个方向需要覆盖先明确边界条件中常见的坑左壁 x0 上的未知LDF不是朝左的而是朝右的。原因是迁移方程中f_i(x_b) 由 x_b - c_i 处的LDF迁移而来只有源点落在计算域外这个分量才是未知的。左壁外侧的粒子要进入计算域必须具有 cx 0 的速度因此未知方向是 1、5、8。边界未知方向索引内部邻居位置左壁 x01, 5, 8x1右壁 xnx-13, 6, 7xnx-2底壁 y02, 5, 6y1顶盖 yny-14, 7, 8yny-2如果把左壁的方向写成了3、6、7顶盖驱动流跑出来镜像效应会特别明显。3.3 顶盖驱动流的非平衡外推实现下面这段代码实现了四个壁面的非平衡外推重点在顶盖处要强加 U_liddef apply_bc(f, U_lid): ny, nx f.shape[:2] rho f.sum(axis2) ux (f[..., 1] - f[..., 3] f[..., 5] - f[..., 6] - f[..., 7] f[..., 8]) / rho uy (f[..., 2] - f[..., 4] f[..., 5] f[..., 6] - f[..., 7] - f[..., 8]) / rho feq equilibrium(rho, ux, uy) # 左壁x0未知方向 1、5、8内部邻居 x1 for i in (1, 5, 8): f[:, 0, i] feq[:, 1, i] f[:, 1, i] - feq[:, 1, i] # 右壁xnx-1未知方向 3、6、7内部邻居 xnx-2 for i in (3, 6, 7): f[:, -1, i] feq[:, -2, i] f[:, -2, i] - feq[:, -2, i] # 底壁y0未知方向 2、5、6内部邻居 y1 for i in (2, 5, 6): f[0, :, i] feq[1, :, i] f[1, :, i] - feq[1, :, i] # 顶盖yny-1未知方向 4、7、8注意给定 uxU_lid ux_top ux.copy() ux_top[-1, :] U_lid uy_top uy.copy() uy_top[-1, :] 0.0 feq_top equilibrium(rho, ux_top, uy_top) for i in (4, 7, 8): f[-1, :, i] feq_top[-1, :, i] f[-2, :, i] - feq_top[-2, :, i]代码中的equilibrium函数就是上一章的平衡态分布函数。这里要注意左右底壁的平衡态是由当前流场ux、uy计算的但顶盖处必须把ux_top[-1, :]强制改成 U_lid否则边界会不断把顶盖速度“抹”回零。3.3.1 角点重叠时的处理顺序四个角的LDF同时属于两个边界非平衡外推的公式在两个壁面上可能会覆盖同一个索引。左下角方向5既被左壁的外推覆盖又被底壁的外推覆盖。常见做法是先处理左右壁再处理底和顶盖让后面的代码覆盖前面。这种处理在顶盖驱动流中是标准且稳定的。如果发现角落出现小范围振荡可以专写一条角点分支分别计算角点密度和速度而不是依赖边界覆盖顺序。4. 顶盖驱动流LBM主循环LDF初始化、碰撞迁移与基准数据对比4.1 完整可运行的最小实现把前面的初始化、碰撞、迁移和边界组合起来就是一个可跑的顶盖驱动流LBM程序。下面代码用NumPy实现网格128×128Re400import numpy as np def equilibrium(rho, ux, uy): feq np.zeros(rho.shape (9,), dtypenp.float64) for i in range(9): cu 3.0 * (cx[i] * ux cy[i] * uy) feq[..., i] omega[i] * rho * ( 1.0 cu 0.5 * cu * cu - 1.5 * (ux * ux uy * uy) ) return feq nx ny 128 U_lid 0.1 Re 400.0 tau 0.5 3.0 * U_lid * nx / Re print(tau , tau) f, rho, ux, uy init_field(nx, ny) for step in range(20000): rho, ux, uy collide(f, tau) f stream(f) apply_bc(f, U_lid) if step % 1000 0 and step 0: # 用中心线速度的最大变化来判断是否收敛 if step % 5000 0: print(step, step, rho_std, rho.std()) u_center ux[:, nx // 2].copy() np.save(centerline_u.npy, u_center)这里的tau是通过Re定义的。顶盖驱动流的雷诺数表达式为Re U_lid × N / ν其中 ν (τ - 0.5) / 3N是顶盖长度对应的格点数。因此τ 0.5 3 × U_lid × N / ReReNU_lidtau说明1001000.10.8适合先跑通4001280.10.596经典验证点10001280.10.538需要更多迭代步10002560.050.538降马赫数更稳定4.2 参加基准数据时看什么顶盖驱动流最常对比的数据是两条中心线上的速度剖面x 0.5 中心线处 ux 沿 y 的分布以及 y 0.5 中线处 uy 沿 x 的分布。Ghia et al. 的经典数据是标准参照。跑完以后把u_center归一化为ux / U_lid再画在同一坐标系里通常能直接看出问题。import matplotlib.pyplot as plt y np.arange(ny) / (ny - 1) plt.plot(u_center / U_lid, y, labelLBM) plt.gca().invert_yaxis() plt.legend() plt.savefig(centerline_velocity.png, dpi150)注意把 y 轴从顶盖y1到底板y0invert_yaxis()这样和文献图的坐标系一致。4.3 各种输运参数对结果的影响如果中心和文献对不齐优先检查三个参数。第一是迭代步数。低Re 100时大约5000步就能得到可接受的剖面Re1000时建议跑30000步以上。第二是 U_lid 的取值。U_lid 太大会让马赫数 Ma U_lid / c_s 超过0.2出现可压缩误差一般控制在0.1以下。第三是网格分辨率。Re1000时128×128只能看趋势256×256的剖面会更接近基准数据但代价是迭代步数和内存同时上升。4.4 运行时报错 nan 的定位顺序LBM代码跑出nan绝大多数情况下不是方程错了而是LDF出现了负值。碰撞后出现负值通常发生在 tau 太接近0.5时边界外推时密度rho出现小于0的值也会引起。排查时可以每隔几百步打印f.min()if step % 500 0: print(step, f.min(), rho.min(), rho.max())一旦f.min()开始小于0先降低 U_lid再把 tau 提高到0.6附近看看。如果仍然负值检查迁移方向索引是否和离散速度数组对齐尤其是对角线方向5到8。5. 顶盖驱动流LBM跑完后用LDF重建的这三个检查最容易被跳过5.1 只用流线图不能发现边界索引错误流线图看起来很漂亮时中心线速度剖面可能误差超过20%。建议每次跑完都直接提取ux[:, nx // 2]把它和Ghia数据画在一起。如果剖面在顶盖附近没有回到 U_lid说明顶盖边界强加速度失败如果在底板附近出现符号反转说明底壁的未知方向索引可能写反了。5.2 用总质量和LDF负值概率检查边界泄漏非平衡外推虽然不保证严格质量守恒但顶盖驱动流运行稳定后全场平均密度应该接近初始值。可以把这段检查放在最后mass_mean f.sum(axis2).mean() neg_ratio (f 0).mean() print(mean rho:, mass_mean, negative LDF ratio:, neg_ratio)如果mass_mean漂移超过0.1%重点看四个角点的LDF覆盖顺序如果neg_ratio大于0说明tau太低或 U_lid 太高。5.3 定位涡心时用流函数而不是压力等值线顶盖驱动流的涡心位置是宏观速度接近零的位置直接从速度场找极值容易受到数值噪声干扰。常见做法是从底边界向上积分流函数psi np.cumsum(ux, axis0) vortex np.unravel_index(np.argmin(psi), psi.shape) vortex_y vortex[0] / (ny - 1) vortex_x vortex[1] / (nx - 1) print(vortex center:, vortex_x, vortex_y)对于顺时针主涡流函数最小值对应的就是涡心。这个值可以和Ghia数据中的涡心坐标对照。Re100时涡心大约在 x≈0.62、y≈0.74附近Re400时会稍微向中心下移并偏向右侧。若涡心明显偏到右上角优先检查顶盖附近的边界外推是不是没有区分已知和未知方向若涡心沿对角线震荡通常说明角点的边界覆盖顺序不一致。本文还有配套的精品资源点击获取