LBM-P多孔介质渗流模拟：从D2Q9演化方程到渗透率验证

简介一份基于LBM_P的格子玻尔兹曼方法MATLAB实现面向需要模拟多孔介质中流体流动的科研人员与工程技术人员。该方法通过离散化微观粒子碰撞与迁移过程在保留流体动力学特征的同时简化计算LBM_P进一步引入渗透率和固液界面处理可准确模拟孔隙、裂缝和喉道中的渗流、污染物迁移及热交换等复杂问题适用于煤炭开采、油藏工程和地质环境评估等场景。压缩包内包含1个m文件仅3KB代码结构紧凑、无多余依赖便于直接阅读、调试与二次开发。已有1018人学习下载适合具备流体力学基础并希望快速上手LBM_P算法的中高级用户。借助该代码读者可掌握多孔介质流动模拟的核心流程灵活调整边界条件与渗透率等参数并针对具体工程问题进行数值实验与结果分析。1. 从多孔介质的细观尺度说起LBM-P 为什么比网格法更顺手多孔介质里的流体模拟标题里出现的“LBM_P”通常不是某个商业软件代号而是一套把格子玻尔兹曼方法LBM和孔隙率参数 P 绑定的建模思路。网格法软件 fluent 或 comsol 做这类问题时先要把多孔骨架几何重新切网格孔隙率一变就要重做前处理LBM 天然把多孔骨架离散在同一个标量场上孔隙率 P 只是格点上的 0/1 标记改标记就是改结构。这就是为什么在岩心渗流、滤芯、催化床、电池电极这类低雷诺数复杂通道问题里LBM-P 经常被列在参考方案的前排。本文把这套方法从演化方程讲到可运行的 Python 代码再落到渗透率验证适合想绕过商业软件前处理流程、自己控制几何和参数的工程师。2. 格子玻尔兹曼 D2Q9 演化方程把孔隙率 P 写进碰撞与反弹2.1 离散速度与权重D2Q9 为什么够用格子玻尔兹曼的核心不是直接解 Navier-Stokes而是在每个格点上维护一组离散速度方向的分布函数通过“碰撞-迁移”两步迭代逼近宏观流动。D2Q9 表示二维 9 个速度方向是二维多孔介质研究的默认起点。9 个方向恰好能满足低马赫数下密度、动量、动量通量的短矩匹配太少会丢各向异性太多对二维收益不大。速度与权重的对应关系是整套代码的地基参数如下。方向索引速度 (cx, cy)权重 w0(0, 0)4/91(1, 0)1/92(0, 1)1/93(-1, 0)1/94(0, -1)1/95(1, 1)1/366(-1, 1)1/367(-1, -1)1/368(1, -1)1/36权重是高斯积分的二维求积节点它们在后面的平衡态分布计算中直接参与。说句实在的这组权重不需要背代码里写成数组就行但要保证索引和速度严格对应否则迁移方向乱掉密度场会在几万步内被无中生有的“能量”冲散。2.2 BGK 碰撞与力项孔隙率 P 的两种进法最常见的 LBM 碰撞模型是 BGK 单松弛近似演化方程写作f_i(x c_i dt, t dt) f_i(x, t) - omega * (f_i - f_i_eq) F_iomega 1 / tautau 是松弛时间决定运动粘度。f_i_eq 是平衡态分布由宏观速度和密度算出满足流动的连续性方程和动量方程。公式略代码里能直接展开。孔隙率 P 在这个框架里有两种完全不同的进法。第一种是显式几何法用一个二维布尔数组标记固体格点流体格点上跑正常碰撞迁移固体格点不参与流体演化用反弹边界施加无滑移。此时孔隙率是几何标记隐含出来的不单独出现在方程里。第二种是体积平均法把整个区域当多孔介质在平衡态里乘上局部孔隙率 P或在力项里加达西阻力。前者适合能画出骨架细节的小尺度模型后者适合大尺度工程模型。多孔介质 LBM 代码中 80% 的算例走的都是显式几何法标题里这个“LBM_P”可以理解为“带孔隙率标记的 LBM”。2.3 反弹边界没有 bounce-back 的多孔 LBM 是失效的骨架边界如果处理错计算出的渗透率会差一个数量级。显式几何法里最常用的是半程反弹流体粒子碰到固体格点速度反向等于在固体表面构造了无滑移条件。实现上就是在迁移之后把固体格点上的分布函数按相反方向交换即方向 1 和 3 交换、2 和 4 交换、5 与 7、6 与 8 两对也分别交换。这个操作用一行 NumPy 索引就能完成但它必须发生在迁移之后。顺序错了结果相当于在固体边界上开了个半格宽度的缝隙渗透率会被高估而且不容易靠肉眼察觉。真正的多孔介质模拟里骨架半径通常只有 5~10 个格点边界误差跨越半个流道足以让整条渗透率曲线失效。2.4 稳定性边界tau 与雷诺数要一起看BGK 的松弛时间 tau 不能太小。tau 越接近 0.5格子粘度越低流动越容易进入不稳定振荡tau 太大又会让流动过度耗散速度场被抹平。常见做法是把 tau 放在 0.6~1.0 之间对应的格子粘度 nu (tau - 0.5) / 3 约在 0.033~0.167。同时要盯住孔隙内的实际雷诺数Re u_pore * d_p / nu。多孔介质中达西速度 q 和孔隙平均流速 u_pore 相差一个孔隙率u_pore q / eps。Re 超过 10 后惯性效应明显达西定律不再适用再用线性渗透率公式算出来的值会随驱动压力变化。LBM-P 适合的算区是 Re 10这个条件在岩石渗流、滤芯过滤场景下基本都满足。3. 用 Python 写一个最小 LBM-P 求解器核心循环与可复现参数3.1 核心循环碰撞、迁移、反弹三步走下面这段代码可以在几分钟内跑出一个二维多孔介质的稳态渗流场。为了让读者能独立复现方向表写全注释按物理步分解。import numpy as np # D2Q9: (cx, cy) 顺序与权重表保持严格一致 C [(0, 0), (1, 0), (0, 1), (-1, 0), (0, -1), (1, 1), (-1, 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]) # 反弹时配对的方向索引 OPP [0, 3, 4, 1, 2, 7, 8, 5, 6] def lbm_step(f, omega, fx, solid): rho f.sum(axis2) # 宏观速度y 方向正方向取向下 ux (f[:,:,1] f[:,:,5] f[:,:,8] - f[:,:,3] - f[:,:,6] - f[:,:,7]) / np.maximum(rho, 1e-12) uy (f[:,:,2] f[:,:,5] f[:,:,6] - f[:,:,4] - f[:,:,7] - f[:,:,8]) / np.maximum(rho, 1e-12) f_new np.empty_like(f) for i, (cx, cy) in enumerate(C): # BGK 碰撞 cu 3.0 * (cx * ux cy * uy) feq W[i] * rho * (1.0 cu 0.5 * cu*cu - 1.5 * (ux*ux uy*uy)) # 外力项只沿 x 方向等效 dP/dx -fx force 3.0 * W[i] * fx * cx f_i f[:,:,i] - omega * (f[:,:,i] - feq) force # 迁移从上游位置取旧值 f_new[:,:,i] np.roll(np.roll(f_i, cx, axis0), cy, axis1) # 固体格点反弹实现无滑移边界 f_new[solid] f_new[solid][:, OPP] return f_new碰撞里force 3.0 * W[i] * fx * cx这一项是线性力对应压力梯度。它不等于直接给每个格点加一个均匀速度而是通过动量注入让流场自己发展出压力驱动下的速度分布。稳态时注入的动量被骨架上的固体边界耗散掉平均速度不再上升此时流场就是达西流状态。迁移部分的np.roll参数值得解释。方向 1 的 (cx, cy) (1, 0)表示一个小时步后格点 (x, y) 上的 f1 来自上游格点 (x-1, y)。np.roll(f_i, 1, axis0)正是把数组沿 x 轴右移一位逻辑上等价于取上游值。反向方向则用负数偏移。这里用周期迁移配套做法是把上下两行格点也标记为固体后面说几何时再展开。3.2 格子单位与物理单位参数映射表LBM 的麻烦在于所有量都是格子单位和 SI 单位之间是缩放关系。参数映射不做调出来的 tau 和 fx 无法对应任何真实流速。物理量格子单位表示换算说明长度1 格 dx算渗透率后按 dx 换算如 1 格 1 μm时间步1 步 dt稳态达到的步数换算到秒密度rho ≈ 1.0不可压缩流入口参考密度取 1运动粘度nu (tau - 0.5) / 3tau 取 0.8 时 nu 0.1压力梯度dp/dx fxfx 是每步每格注入的动量渗透率K q * nu / fx结果单位是格²再乘 dx² 到 m²一个实用的经验是先设格子粘度 nu再根据目标物理粘度反推 tau网格尺寸定了之后格子长度 dx 自然决定速度量级。fx 取 1e-5 到 1e-6 比较稳。如果 fx 过大局部速度超过 0.1 格子/步压缩性误差会显现压力场出现波纹。3.3 跑通主流程初始化、残差监控、稳态判定上层驱动代码需要完成初始化、残差监控和固体区域设置。初始时全流场密度设为 1分布函数直接赋平衡态速度取 0。残差可以监控速度场两次采样之间的最大变化降到 1e-6 以下视为稳态。def run_simulation(nx, ny, solid, tau0.8, fx1e-5, steps5000): omega 1.0 / tau f np.zeros((nx, ny, 9)) rho_init 1.0 for i, (cx, cy) in enumerate(C): f[:,:,i] W[i] * rho_init # 零速度平衡态 for step in range(steps): f lbm_step(f, omega, fx, solid) if step % 500 0: # 监控全局平均速度是否进入平稳段 rho f.sum(axis2) ux (f[:,:,1] f[:,:,5] f[:,:,8] - f[:,:,3] - f[:,:,6] - f[:,:,7]) / np.maximum(rho, 1e-12) umax np.abs(ux[~solid]).max() print(fstep {step}, u_max {umax:.6f}) return f需要注意solid数组必须用布尔类型并且在计算宏观速度时要排除固体格点。固体内分布函数在反弹后仍保留数值但这些值没有物理意义统计速度时要用[~solid]掩码过滤。判断稳态不要只看某个格点要看全场最大速度的变化率和全局密度漂移后者能暴露隐藏的质量源或泄漏边界。4. 反向验证随机圆盘堆积的渗透率计算与 Kozeny-Carman 对照4.1 用数组刻多孔骨架六方排列圆盘网格法里多孔几何要生成 STL 再划分网格LBM-P 里几何就是布尔数组。这里生成六方排列的圆盘骨架接触点可控方便算孔隙率。def make_hex_porous(nx, ny, R, gap2): solid np.zeros((nx, ny), dtypebool) row 0 y0 R 2 # 避开最外圈固体壁面 while y0 R ny - 2: if row % 2 0: xs np.arange(R 2, nx - R, 2 * R gap) else: xs np.arange(2 * R gap, nx - R, 2 * R gap) R for x in xs: yy, xx np.ogrid[:ny, :nx] solid[(xx - x)**2 (yy - y0)**2 R**2] True y0 int(np.sqrt(3) * (R gap / 2)) row 1 # 上下两行强制设固体配合周期迁移形成壁面 solid[0, :] True solid[-1, :] True eps 1.0 - solid.mean() return solid, eps圆盘半径 R 是影响计算量的关键参数。R 5 时每个圆大概占 78 个格点孔隙通道只有 3~4 个格点宽反弹边界误差占比偏高R 10 以上更稳。gap控制圆盘间距gap2时可减少接触点处的奇异性也让迁移过程中邻格更平滑。接触点若完全相切反弹会出现“尖角效应”局部速度偏高渗透率偏大。4.2 渗透率怎么从速度场里剥出来稳态速度场中包含了两类速度孔隙内平均速度 u_pore 和达西速度 q。计算渗透率必须用达西速度即总流量除以总面积而不是孔内平均速度。两者相差一个孔隙率倍数。# 假设 f 已通过 run_simulation 达到稳态 rho f.sum(axis2) ux (f[:,:,1] f[:,:,5] f[:,:,8] - f[:,:,3] - f[:,:,6] - f[:,:,7]) / np.maximum(rho, 1e-12) fluid ~solid u_pore np.mean(ux[fluid]) # 孔隙空间平均流速 eps fluid.mean() # 孔隙率 q u_pore * eps # 达西速度 nu (0.8 - 0.5) / 3 # 与 tau 对应 # 达西定律q (K / nu) * fx K_lbm q * nu / fx # Kozeny-Carman 经验公式对照 dp 2 * R # 特征粒径 K_kc eps**3 * dp**2 / (150 * (1 - eps)**2) print(fLBM K {K_lbm:.4f}, KC K {K_kc:.4f}, ratio {K_lbm / K_kc:.2f})这段代码最后算出的K_lbm单位是格²。换算思路若 1 格 1 μm则 1 格² 1 μm² 1e-12 m²。多孔介质渗透率常用毫达西表示1 mD 约等于 1e-15 m²所以 K_lbm 0.05 时对应约 50 mD量级接近粉砂岩。换算看起来简单但漏乘 dx² 是最常见的低级错误。4.3 和理论公式的偏差范围怎么判断算对了Kozeny-Carman 公式本身是经验关系系数 150 对圆盘规则排列并不是绝对准确合理的比值落在 0.8~1.5 之间。如果比值偏差超出一个数量级优先检查三件事孔隙率是否和预期一致、孔隙通道分辨率是否足够、fx 是否落在线性渗流区。更细的验证做法是扫一组孔隙率。固定 R改变 gap得到孔隙率从 0.4 到 0.7 的多个算例分别算 K_lbm和 Kozeny-Carman 画在同一张图上。趋势一致基本说明边界处理和力项正确。若只在某个孔隙率偏差多半是骨架接触点导致局部网格过渡区问题增大 R 重新测。5. LBM-P 的调试顺序与收敛判定从 tau 到 REV5.1 先看三个体检指标拿到算例结果后不要急着算渗透率先做三个快速体检。第一全局密度总和漂移是否稳定理想情况每 100 步变化小于 1e-6若持续上涨说明边界泄漏。第二速度场中固体内部是否出现明显非零速度出现则反弹顺序错或 solid 掩码没生效。第三沿流向的截面压力分布应该接近线性弯曲毛刺说明通道内有未处理的高频振荡此时调小 fx 或增大 tau。5.2 REV 收敛性算一次不算数多孔介质模拟最终要回答“这个结构的渗透率是多少”单次算例容易受随机几何和边界影响。正确做法是算代表性体积单元REV固定孔隙率和粒径比例把计算域放大 1.5 倍和 2 倍分别重算 K_lbm。若 K_lbm 变化小于 5%认为当前网格规模和边界条件组合达到了代表性若 K 翻倍说明骨架稀疏单个样本没有代表性需要扩大区域或多取几个随机种子取平均。R 5 的圆盘在 200×100 网格下往往不够 REVR 8~10、网格 500×200 时才稳定。这一步少则占用额外一倍的机时但不做的话后面所有对比都可能建立在偶然几何上。5.3 用达西区和 Forchheimer 区判定收尾线性渗透率只对低 Re 成立而同一个多孔结构在不同驱动压力下会呈现不同的流动状态。最终报告里给出渗透率时建议额外输出一个 Re 值Re_pore q * dp / nu if Re_pore 5: print(Darcy regime, linear K is valid) else: print(进入 Forchheimer 区需补充二次阻力项)Re_pore 超过 10 后速度与压差不再是线性关系再用单一渗透率描述会误导后续工程换算。把 tau、fx、REV 三个参数的验证顺序固化下来以后换几何、换流体都不用重新猜参数。本文还有配套的精品资源点击获取