企业官网建设流程全解析

简介：格子玻尔兹曼方法（LBM）在多孔介质流体模拟中的应用是计算流体力学的重要方向，这份基于LBM_P的MATLAB代码面向煤炭、石油、岩土等领域科研人员与工程师，适用于渗流分析、污染物迁移及流固耦合等场景。包内共1个m文件，压缩包约3KB，代码结构紧凑，便于直接读取并调整边界条件、渗透率等参数，借此重现多孔介质中的速度场与压力分布，也可作为拓展LBM模型的起点。目前已有1018人学习下载，适合具备一定流体力学与MATLAB基础、希望快速上手LBM_P算法的中高级使用者。通过研读代码，可更直观理解LBM_P如何处理复杂孔隙结构与固液界面，并围绕瓦斯排放预测、注水采油优化、地下水污染扩散评估等实际工程问题开展数值实验。

1. 多孔介质流动为什么找上格子玻尔兹曼

地下水的迁移、燃料电池气体扩散层里的氧气输运、催化剂载体里的反应物扩散，最终都要归结到孔隙尺度的流动。这类几何随机、孔道曲折、边界复杂的算例，用宏观 CFD 先吃一个贴体网格的生成，往往比解方程本身还费劲。格子玻尔兹曼方法（Lattice Boltzmann Method，LBM）恰好把这一步绕开了，它用规则格点加一个布尔数组记录固体位置，流体在格点上碰撞、迁移，遇到固体格点直接反弹。LBM_P 里这个 P 指的就是 Porous：同一套 D2Q9 内核，把几何换成多孔骨架，就能算渗透率、迂曲度和内部流场。这套路径适合做渗流力学、岩心分析、燃料电池和过滤材料的人，新手按格点逻辑也能把结果跑出来。

2. 格子玻尔兹曼的孔道视角：从 BGK 碰撞到多孔骨架边界

2.1 从分布函数到宏观流速：D2Q9 与 BGK 的演化

多孔介质里的流动一般雷诺数很小，很多场景甚至落在蠕动流区间，宏观上仍由纳维-斯托克斯方程支配。但 LBM 不直接解 NS 方程组，而是追踪粒子分布函数 (f_i)，在二维最常用 D2Q9 离散速度模型。演化过程可以压缩成三步：碰撞、迁移、外力修正。BGK 单松弛碰撞算子的演化方程写出来就是

[ f_i(x+c_i \Delta t, t+\Delta t) = f_i(x,t) - \frac{1}{\tau}\left(f_i - f_i^{eq}\right) + F_i ]

其中 (\tau) 是无量纲松弛时间，与运动粘度 (\nu) 的关系在 D2Q9 里是 (\nu = (\tau-0.5)/3)。宏观密度和速度通过分布函数的零阶矩和一阶矩恢复：

[ \rho = \sum_i f_i, \quad u = \frac{1}{\rho}\sum_i f_i c_i ]

D2Q9 的九个离散速度方向是（0,0）、(±1,0)、(0,±1)、(±1,±1)，对应的权重分别为 4/9、1/9、1/36。这个权重表格是整套代码的起点。

平衡分布函数是 LBM 把介观量拉回宏观流场的桥梁。下面这段代码把平衡分布写成可复用的函数，之后不管几何怎么换，这层不用动。

def equilibrium_d2q9(rho, ux, uy, w, ex, ey): # 计算 D2Q9 的平衡分布函数 feq = np.zeros((w.size,) + rho.shape) u2 = ux * ux + uy * uy for i in range(w.size): eu = ex[i] * ux + ey[i] * uy feq[i] = w[i] * rho * (1.0 + 3.0 * eu + 4.5 * eu * eu - 1.5 * u2) return feq

这段代码的核心逻辑是让平衡分布随局部密度和速度重新分配权重。初始时若流场速度为零，rho 取 1.0，那么九个方向的分布正好就是权重本身；速度起来后，eu 的高阶项会补偿动量的方向性。参数里 tau 通常取 0.55~1.0，因为 (\tau) 太靠近 0.5 时粘度近似为零，数值上会出现高频振荡；取太大则流动进入强耗散区，压力梯度需要放很大才能推动流体，导致低速区域的舍入误差放大。

2.2 生成多孔几何的三种路径：规则阵列、随机圆盘、数字岩心

多孔介质 LBM 的“几何”说到底就是一张布尔掩膜，固体格点记 True，流体格点记 False。换结构不用改求解器，只换掩膜，这是 LBM 在多孔问题上比贴体网格方案省事得多的根本原因。常见做法有三条路径。

规则阵列适合先检验代码：一排圆柱排成正方形或菱形，Darcy 渗透率在文献里有对照值，跑完能确认边界格式没有系统性偏差。随机圆盘是开发阶段性价比最高的选择，生成快且孔隙率和孔喉分布都接近天然松散介质。数字岩心最接近真实，但需要先做图像滤波、连通域分析和分辨率判断；如果喉道直径只有四五个格点，LBM 的反弹边界误差会盖过真实信号，这时候再好的 CT 数据也救不回来。

2.3 多孔介质边界的反弹细节：有效壁面位置与驱动方式

LBM 处理固体边界的最简单方案是反弹边界：迁移时分布函数进入固体格点后，沿来路反向弹回。这种格式的实现几乎不增加计算量，但代价是有效壁面位置不在流体格点中心，也不在固体格点中心，而是在两者之间的中点。对孔隙介质来说，这意味着一个名义上 8 格宽的喉道，实际流动截面可能只有 7 格，渗透率偏差可能达到百分之十几。孔隙越窄，这个偏差越致命。所以多孔介质 LBM 的一个前置要求是喉道宽度至少要 10 个格点，经验上取 15~20 格更稳妥。

驱动方式同样影响结果。入口出口做成压力边界时，Zou-He 格式的数值反射容易在孔隙内部激起微弱的伪振荡，收敛后平均速度也许看不出问题，但局部流场会带着非物理条纹。更稳的做法是在整个计算域施加恒定的体积力，等效于一个均匀压力梯度，配合周期边界实现无限长多孔介质段。体积力驱动的另一个好处是压力梯度的数值大小直接由外力强度给出，后续由达西定律反算渗透率时少一道对进出口压差的统计误差。

3. 组装一个可以算渗透率的多孔介质 LBM 流程

3.1 用随机圆盘生成带目标孔隙率的多孔骨架

生成多孔几何时，最常遇到的问题是孔隙率看起来对，实际流动通道却被少量“死胡同”占掉一大块。以下代码生成非重叠随机圆盘骨架，并用取模操作让圆盘跨边界时自动周期复制，保证后面周期边界可用。

import numpy as np nx, ny = 256, 256 target_phi = 0.65 radius = 4 rng = np.random.default_rng(2024) solid = np.zeros((nx, ny), dtype=bool) needed = int((1.0 - target_phi) * nx * ny) while solid.sum() < needed: cx, cy = rng.integers(0, nx), rng.integers(0, ny) for dx in range(-radius, radius + 1): for dy in range(-radius, radius + 1): if dx * dx + dy * dy <= radius * radius: # 取模实现周期复制，避免边界处圆盘被截断 x, y = (cx + dx) % nx, (cy + dy) % ny solid[x, y] = True porosity = 1.0 - solid.mean() print("actual porosity:", porosity)

这段循环每次投一个圆盘，每次增加约 (\pi r^2) 个固体格点，所以实际孔隙率会略低于目标值。若要精确到小数点后两位，可以在投盘前预计算本次新增格点数，超出 needed 时跳过一次投放。这里的 radius 是圆盘半径的格点数量，取 4 时喉道宽度约为 8 格，只适合快速演示；正式计算建议把圆盘半径放大到 8~10 格，同时把计算域扩大到 512 以上，否则孔隙率统计和渗透率结果都不可信。

3.2 碰撞、迁移、反弹三段主循环的实现

有了掩膜后，主循环是标准的碰撞-迁移-反弹三步。碰撞前先恢复宏观量，碰撞时只处理流体格点；迁移用 np.roll 实现周期位移；反弹在固体格点内交换反向分布函数。

# D2Q9 基础数组 w = np.array([4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36]) ex = np.array([0, 1, -1, 0, 0, 1, -1, 1, -1]) ey = np.array([0, 0, 0, 1, -1, 1, -1, -1, 1]) opp = np.array([0, 2, 1, 4, 3, 6, 5, 8, 7]) tau = 0.8 Fx = 1e-5 # 体积力，等效压力梯度 fluid = ~solid f = np.zeros((9, nx, ny)) for i in range(9): f[i] = w[i] # 初始化为均匀平衡分布，rho=1，u=0 for step in range(100000): # 恢复宏观量 rho = f.sum(axis=0) ux = (f[1] - f[2] + f[5] - f[6] + f[7] - f[8]) / rho uy = (f[3] - f[4] + f[5] + f[6] - f[7] - f[8]) / rho feq = equilibrium_d2q9(rho, ux, uy, w, ex, ey) # 碰撞：只对流体格点 f[:, fluid] = f[:, fluid] - (f[:, fluid] - feq[:, fluid]) / tau # 体积力驱动，采用外力项加权近似 for i in range(9): f[i, fluid] += w[i] * Fx * ex[i] # 迁移：周期边界 for i in range(9): f[i] = np.roll(np.roll(f[i], ex[i], axis=0), ey[i], axis=1) # 反弹：固体格点内交换反向分布 for i in range(1, 9): if i < opp[i]: tmp = f[i, solid].copy() f[i, solid] = f[opp[i], solid] f[opp[i], solid] = tmp if step % 1000 == 0: res = np.abs(ux[fluid] - ux_prev[fluid]).mean() if res < 1e-9: break ux_prev = ux.copy()

这里的碰撞项把分布函数按松弛时间拉向平衡态，(\tau=0.8) 对应运动粘度 (0.1) 格子单位，粘性耗散适中。体积力项加在迁移之前，外力会沿 x 方向驱动流体；注意这是教学用的简化外力实现，严格做法需要在平衡速度里加入 ( \tau F/\rho ) 的修正项，孔隙率极低或速度梯度很大时两种写法会有可观测差异。迁移后分布函数会流入固体格点，反弹段将固体格点内的分布函数成对反向交换，等效于把进入固体的粒子弹回流体，有效无滑移边界落在流体与固体格点之间。最后用流体格点的速度变化残差做稳态判据，这个阈值要参考入口体积力和孔隙率来调，外力越小收敛越慢，但残差阈值不能跟着放松，否则渗透率会被未稳定的慢速流动带偏。

3.3 多孔介质 LBM 的关键参数表与量纲约定

把格子单位换算到物理单位是新手最容易翻车的地方。LBM 计算里长度、时间和质量各有一个独立尺度，渗透率换算公式为 ( k_{phys} = k_{lb} \cdot (\Delta x_{phys}/\Delta t_{phys})^2 )，其他量都由这三者推出来。下表给出我一般会用的初始化参数区间。

孔隙雷诺数用孔隙平均速度乘上孔喉特征长度再除以运动粘度，在格子单位里可以直接算。如果 Re 大于 1，惯性项不可忽略，Darcy 定律的线性关系不再成立，这时候再套渗透率公式会得到一个随驱动力变化的人造渗透率。Knudsen 数以分子自由程与孔喉直径之比估计，格点单位下很难直接给出，但可以记住一条经验：tau 离 0.5 越远，非连续效应越小；tau 超过 1.2 后反而是强耗散掩盖了物理粘度。

4. 从多孔介质 LBM 结果反推渗透率：达西定律与五个稳定性坑

4.1 达西定律的渗透率提取与单位换算

稳态后的流场可以直接套达西定律算渗透率。对体积力驱动的周期结构，宏观压力梯度就是 (-\Delta p/L = Fx)，体积平均速度取流体格点上的速度平均值，公式是

[ k_{lb} = \frac{\nu_{lb} \cdot \langle u_x \rangle_{fluid}}{F_x} ]

代码实现很简短：

# 稳态后的渗透率（格子单位） nu_lb = (tau - 0.5) / 3.0 u_avg = ux[fluid].mean() # 孔隙内体积平均速度 k_lb = nu_lb * u_avg / Fx # Fx 为体积力，等效压力梯度 # 换算到物理单位 dx_phys = 1e-6 # 每个格点对应 1 微米 dt_phys = 1e-8 # 时间步，由粘度匹配确定 k_phys = k_lb * (dx_phys / dt_phys) ** 2 print("permeability (m^2):", k_phys)

注意 ux[fluid].mean() 是流体格点速度的算术平均，它没有乘孔隙率，因为取平均的集合已经是流体域。若要当场验证孔隙率是否真的参与了流动，可以对照经验公式 Kozeny-Carman：( k \approx d_p^2 \phi^3 / (180(1-\phi)^2) )，其中 (d_p) 是固体圆盘直径的物理尺寸。两者能对上量级就说明主循环没有系统性错误；偏差超过两三倍时，先检查喉道宽度是否足够，再检查反弹边界是否把有效截面缩了。单位换算里，(\Delta t_{phys}) 不是随便定的，常见做法是先定 (\Delta x_{phys})，再根据希望匹配的物理运动粘度反推时间步，使得格子粘度和物理粘度满足同一个无量纲数。

4.2 多孔介质 LBM 常见的五个坑与对应排查手段

跑多孔介质 LBM 时，渗透率结果不对的原因通常很集中。下面五类问题我基本每次都先排查一遍。

最后一条在多孔介质里尤其隐蔽。格子玻尔兹曼的反弹边界有效位置在流体与固体格点之间，等效于把固体圆盘的半径扩大 0.5 格；宏观 CFD 的壁面严格落在几何边界上。小孔隙下这 0.5 格误差会显著改变渗透率。处理方法是在统计几何参数时把固体边界外扩 0.5 格作为实际流动边界，或者用更精细的插值反弹边界格式来降低壁面位置误差。对多数工程判断来说，先保证喉道分辨率够、再用两种尺度验证，比花大力气改边界格式更划算。

5. 进阶技巧：用迂曲度复核多孔介质 LBM 的稳态结果

5.1 迂曲度的格点算法

渗透率是标量，能描述整体阻力，却无法暴露流场内部的细节错误。迂曲度（tortuosity）描述的是流体实际绕行路径与宏观直线距离的比值，对骨架产生的流动异常敏感。稳态流场中，迂曲度可以用速度矢量的分布直接估计：

[ T = \frac{\langle |\mathbf{u}| \rangle_{fluid}}{\langle |u_x| \rangle_{fluid}} ]

分子是速度模长的体积平均，分母是轴向速度绝对值的体积平均。代码实现如下：

uxf = ux[fluid] uyf = uy[fluid] speed = np.sqrt(uxf**2 + uyf**2) tort = speed.mean() / np.abs(uxf).mean() print("tortuosity:", tort)

这段计算的核心是区分平均方向与平均路径长度。常见错误是用 (\sqrt{\langle u_x^2\rangle}) 做分母，这会因对称流动中 (u_y) 的均值几乎为零，结果把横向脉动当成迂曲度的一部分，算出来的 T 系统性偏大。正确做法是保留绝对值。对随机圆盘骨架，孔隙率 0.65 时迂曲度通常在 1.2~1.5 之间，孔隙率越低，T 越大；如果算出来 T 小于 1.05，多半是流场没有充分发展或体积力已经超过达西区间，流线被拉直了。

5.2 与 Fluent、COMSOL 的对照验证方法

等渗透率和迂曲度都从 LBM 里出来，建议再和宏观流体仿真做一次交叉验证。Fluent 和 COMSOL 都能导入同一张二值掩膜生成几何，但这里有一个常被忽略的细节：对照时不能只比渗透率。渗透率是一个积分量，不同的局部流场分布可以凑出同一个值，所以先把速度场导出来，再算一次同样的迂曲度和沿主流方向的流线分布。COMSOL 里可以直接用速度分量构造自定义表达式，Fluent 则在后处理里用 Custom Field Function 定义绝对速度和轴向速度，两者计算迂曲度的口径要保持一致。LBM 速度场与 COMSOL 或 Fluent 的偏差维持在 5% 以内，基本可以判定反弹边界和多孔骨架的设置没有问题；偏差集中在某个局部区域时，优先检查该区域的喉道分辨率和二值化阈值，而不是怀疑求解器本身。

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