简介：这份资源提供了一套基于 Matlab 实现的格子玻尔兹曼方法（LBM）流体仿真代码，专门用于模拟流体流经多孔介质的流动过程，面向计算机、电子信息工程、数学等专业的大学生及科研人员，可服务于课程设计、期末大作业与毕业设计等场景。压缩包共 11 个文件，包含 7 个 m 脚本、3 张 png 结果图与 1 份 pdf 说明文档，整体约 3.45MB，脚本覆盖泊肃叶流、方腔流及多孔介质流动等多个算例，图片与文档便于对照验证仿真结果。代码采用参数化编程，参数修改方便，注释详尽、思路清晰，附赠案例数据可直接运行，读者能借此快速理解 LBM 的碰撞与迁移迭代流程，掌握多孔介质边界处理与数值实验方法，并在此基础上开展算法改进与工程分析。目前已有 82 人学习。

1. 从一张砂芯截面的流速云图说起：Matlab 做 LBM 多孔介质仿真到底能干什么

如果你手头有一块电池电极、一段岩心或者一张金属泡沫的 CT 切片，想知道流体从一端压进去之后内部每一条孔隙通道的流速分布，商用 CFD 前处理往往先把你劝退：几何清理、体网格划分、孔隙尺度网格量爆炸，一套流程走下来一天就没了。格子玻尔兹曼方法（LBM）在这类问题上有个天然优势——它本来就是定义在规则格子上的介观方法，多孔骨架只需要标记成固体节点，不需要贴体网格，用 Matlab 就能把整个求解器写清楚。这篇笔记讲的就是怎么用 Matlab 从零搭一个能跑通多孔介质渗流的 LBM 求解器，包括 D2Q9 模型、反弹边界、孔隙率与渗透率的提取，以及那些让结果发散的参数坑。适合做电池、岩土、多孔材料、微流控方向、手上有 Matlab 又不想被商用软件绑住的研究生和工程师。整套代码量不大，核心求解器两百行以内，但每一步的物理含义和数值边界都得交代清楚，否则跑出来的速度场就是一张好看的废图。

2. LBM 求解多孔介质渗流的物理底座与 D2Q9 离散

2.1 为什么多孔介质流动适合用 LBM 而不是传统 CFD

传统有限体积法处理多孔介质有两条路：一是把多孔区当成黑箱，用 Darcy 定律加经验阻力源项，好处是快，坏处是拿不到孔隙内部的速度场，孔隙尺度效应全靠经验系数兜；二是做孔隙尺度贴体网格，能拿到真实流场，但多孔骨架的比表面积大、喉道窄，网格量随孔隙率下降急剧上升，而且骨架与流体交界面的网格质量很难保证。LBM 走的是第三条路：把流体离散成在规则格子上迁移和碰撞的分布函数，固体骨架就是一组不参与迁移的格子节点，流体撞上去按反弹规则原路返回。这样几何处理退化成一张 0/1 掩膜图，孔隙率、连通性、各向异性全部由掩膜决定，前处理成本几乎为零。

从物理上讲，LBM 求解的是介观 Boltzmann 方程的离散速度版本，在低马赫数下通过 Chapman-Enskog 展开可以恢复到不可压 Navier-Stokes 方程。多孔介质渗流通常处于低雷诺数、低马赫数区间，正好落在 LBM 的舒适区。更关键的是，LBM 的局部性极强，碰撞和迁移只依赖自身和最近邻节点，天然适合用矩阵化操作在 Matlab 里向量化实现，不需要稀疏线性求解器，也不会有传统压力泊松方程那种全局迭代收敛问题。代价是时间步进显式，稳定域受松弛时间限制，这个后面避坑章节会重点讲。

2.2 D2Q9 模型、平衡态分布与宏观量提取

D2Q9 是二维九速度模型，九个离散速度方向覆盖了自身静止、四个轴向、四个对角。权重和速度向量是固定的，写死在代码里即可。平衡态分布函数由密度和速度决定，碰撞用 BGK 单松弛近似，迁移就是把每个方向的分布函数沿对应速度搬到邻居节点。宏观密度是九个方向分布函数之和，宏观速度是分布函数加权速度之和除以密度。

% D2Q9 离散速度与权重（行向量，便于向量化） ex = [0, 1, 0, -1, 0, 1, -1, -1, 1]; ey = [0, 0, 1, 0, -1, 1, 1, -1, -1]; w = [4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36]; % 平衡态分布函数，rho 与 ux,uy 为同尺寸矩阵 function feq = equilibrium(rho, ux, uy, ex, ey, w) feq = zeros([9, size(rho)]); for k = 1:9 eu = ex(k)*ux + ey(k)*uy; feq(k,:,:) = w(k) * rho .* (1 + 3*eu + 4.5*eu.^2 - 1.5*(ux.^2 + uy.^2)); end end

这段代码里ex、ey、w是 D2Q9 的标准配置，顺序不能乱，因为迁移步骤会按这个顺序索引邻居。平衡态公式里的系数 3 和 4.5 来自声速的平方为 1/3，这是格子单位下的约定，不要改成物理声速。rho和ux、uy是二维矩阵，feq是三维数组，第一维是方向。用三维数组而不是元胞数组，是为了后面迁移时能直接用 circshift 做向量化，速度差一个数量级。

2.3 多孔介质掩膜生成与反弹边界

多孔介质在 LBM 里就是一张逻辑掩膜，solid = true的节点不参与碰撞迁移，流体节点撞到固体时分布函数沿原方向反弹。生成掩膜有三种常见做法：规则圆柱阵列、随机圆颗粒堆积、从 CT 图像二值化导入。规则阵列适合验证，随机堆积适合看统计规律，CT 导入适合真实材料。下面给一个随机圆颗粒堆积的掩膜生成，控制孔隙率在 0.6 左右。

% 生成随机圆颗粒多孔介质掩膜 Nx = 200; Ny = 100; % 网格尺寸 target_porosity = 0.6; % 目标孔隙率 solid = false(Ny, Nx); rng(1); % 固定随机种子，保证可复现 while nnz(~solid)/numel(solid) > target_porosity cx = randi(Nx); cy = randi(Ny); r = randi([3, 6]); % 颗粒半径范围 [X, Y] = meshgrid(1:Nx, 1:Ny); circle = (X-cx).^2 + (Y-cy).^2 <= r^2; solid = solid | circle; end porosity = nnz(~solid)/numel(solid); fprintf('实际孔隙率 = %.4f\n', porosity);

rng(1)是后悔药，不固定种子每次跑出来的孔隙结构都不一样，渗透率对比就失去意义。颗粒半径范围决定比表面积，半径越小比表面积越大，相同孔隙率下渗透率越低。while循环的终止条件是孔隙率降到目标值以下，实际孔隙率会略低于目标值，因为最后一个颗粒可能压过头。如果需要精确孔隙率，改成先随机撒点再按面积排序逐个添加，或者用形态学腐蚀微调。掩膜生成后建议用imshow或imagesc看一眼连通性，孤立孔隙对渗透率没有贡献，但会拉低孔隙率统计值。

3. 用 Matlab 把 LBM 主循环跑起来：迁移、碰撞与收敛判据

3.1 初始化、时间步进与向量化迁移

主循环的结构是固定的：初始化平衡态，每个时间步先迁移再碰撞，然后施加边界条件，最后提取宏观量并判断收敛。迁移用 circshift 实现，注意 circshift 的方向和离散速度符号要对应，这是最容易翻车的地方。碰撞就是 BGK 松弛，松弛时间由粘度决定。

% 初始化 rho = ones(Ny, Nx); ux = zeros(Ny, Nx); uy = zeros(Ny, Nx); f = equilibrium(rho, ux, uy, ex, ey, w); tau = 0.8; % 松弛时间，决定粘度 omega = 1/tau; % 松弛频率 maxIter = 20000; tol = 1e-6; for iter = 1:maxIter % 迁移：沿每个离散速度方向搬移 fprev = f; for k = 1:9 f(k,:,:) = circshift(squeeze(fprev(k,:,:)), [ey(k), ex(k)]); end % 碰撞：BGK 松弛 rho = squeeze(sum(f, 1)); ux = squeeze(sum(f .* reshape(ex,9,1,1), 1)) ./ rho; uy = squeeze(sum(f .* reshape(ey,9,1,1), 1)) ./ rho; feq = equilibrium(rho, ux, uy, ex, ey, w); f = f - omega * (f - feq); % 反弹边界：固体节点上的分布函数反向 for k = 1:9 kk = 10 - k; % 反向方向索引 fk = squeeze(f(k,:,:)); fk(solid) = squeeze(fprev(kk,:,:))(solid); f(k,:,:) = reshape(fk, 1, Ny, Nx); end % 收敛判据 if mod(iter, 100) == 0 err = max(abs(ux(:) - ux_old(:))); if err < tol, break; end ux_old = ux; end end

迁移步骤里circshift的第二个参数是[ey(k), ex(k)]，先 y 后 x，因为 Matlab 矩阵第一维是行即 y 方向。这个顺序写反了流场会朝错误方向跑，而且不会报错，属于典型的玄学 bug。反弹边界里kk = 10 - k是 D2Q9 的反向索引映射，前提是速度向量按标准顺序排列，如果自己改了顺序这个映射就失效。反弹实现用的是上一时间步的分布函数，所以要先存fprev。收敛判据用速度场的最大变化量，每 100 步检查一次，阈值 1e-6 对格子单位来说比较严格，实际调试可以放宽到 1e-5。

3.2 压力驱动与速度驱动的入口出口边界

多孔介质渗流通常用压力驱动，入口出口给定密度差，或者用速度驱动给定入口速度。压力驱动更接近 Darcy 实验条件，实现也简单：入口列密度设为rho_in，出口列密度设为rho_out，中间线性插值作为初始猜测。速度驱动需要在入口施加平衡态分布并修正密度，实现稍复杂，但收敛更快。

% 压力驱动边界：入口高密度，出口低密度 rho_in = 1.0 + 0.001; rho_out = 1.0 - 0.001; % 每个时间步在碰撞后施加 rho(:,1) = rho_in; rho(:,end) = rho_out; ux(:,1) = 0; uy(:,1) = 0; ux(:,end) = 0; uy(:,end) = 0; f(:,:,1) = equilibrium(rho(:,1), ux(:,1), uy(:,1), ex, ey, w); f(:,:,end) = equilibrium(rho(:,end), ux(:,end), uy(:,end), ex, ey, w);

密度差 0.002 对应很小的压力梯度，保证马赫数远低于 0.1，压缩性误差可以忽略。密度差太大会导致入口附近出现可压缩效应，速度场偏离 Darcy 线性关系。入口出口的速度设为零是简化处理，更严格的做法是用零梯度外推，但对低雷诺数渗流影响不大。上下边界用周期性边界还是固壁边界取决于物理问题，如果模拟的是岩心夹持器内的流动，上下应该是固壁无滑移，用反弹即可；如果模拟的是周期性单胞，上下用周期性边界，用 circshift 在 y 方向首尾相接。

3.3 渗透率提取与 Darcy 定律验证

跑收敛之后，提取整个域的平均速度，结合压力梯度和流体粘度，用 Darcy 定律反算渗透率。这是验证求解器正确性的关键一步，也是多孔介质仿真的核心输出。

% 提取渗透率 u_avg = mean(ux(~solid)); % 流体节点平均 x 速度 mu = (tau - 0.5) / 3; % 格子单位粘度 dpdx = (rho_in - rho_out) / Nx / 3; % 压力梯度，声速平方 1/3 K = mu * u_avg / dpdx; % Darcy 渗透率，格子单位 fprintf('平均速度 = %.6e, 粘度 = %.6e, 渗透率 = %.6e\n', u_avg, mu, K);

粘度公式(tau-0.5)/3来自 BGK 的 Chapman-Enskog 展开，tau 必须大于 0.5，否则粘度为负，直接发散。压力梯度里除以 3 是因为格子声速平方为 1/3，密度差乘以声速平方才是压力差。渗透率是格子单位，要换算成物理单位需要确定特征长度和特征时间，通常用颗粒直径作为特征长度，用粘度和密度确定特征时间。验证方法是改变压力梯度看平均速度是否线性变化，如果线性度好说明处于 Darcy 区，如果偏离说明进入了惯性区，需要改用 Forchheimer 修正。

4. 多孔介质 LBM 仿真的避坑与排查清单

4.1 松弛时间接近 0.5 导致发散

现象：迭代几百步后速度场出现棋盘振荡，残差突然增大到 NaN。原因：tau 太接近 0.5，粘度趋近于零，BGK 碰撞的稳定性急剧下降，格子尺度上的数值噪声被放大。解决：把 tau 保持在 0.6 以上，最好在 0.7 到 1.0 之间。如果物理问题要求低粘度，改用多松弛时间 MRT 模型，或者用 entropic LBM，但代码复杂度会上升。我一般先在 tau=0.8 下跑通，确认流场合理后再调低 tau 看稳定边界。

4.2 孔隙率统计正确但渗透率为零

现象：掩膜孔隙率显示 0.6，但平均速度几乎为零，流场只在入口附近有扰动。原因：随机颗粒堆积时出现了贯穿整个 y 方向的固体墙，把流场在 x 方向截断了，流体无法从入口流到出口。解决：生成掩膜后做一次连通性检查，用bwlabel标记流体连通域，只保留与入口出口相连的连通域，其余孤立孔隙填成固体。更稳妥的做法是生成掩膜时在 y 方向留出若干条畅通通道，或者用更小的颗粒半径降低堵塞概率。

4.3 circshift 方向写反导致流场镜像

现象：入口在左侧，但速度场显示流体从右侧流入，或者流场上下颠倒。原因：circshift 的位移参数顺序是[行位移, 列位移]，行对应 y，列对应 x，而离散速度向量里 ex 是 x 方向，ey 是 y 方向，写代码时容易把[ey(k), ex(k)]写成[ex(k), ey(k)]。解决：在迁移循环里加一句断言，或者先用一个已知的均匀流场测试迁移步骤，均匀流场迁移后应该保持不变。这个 bug 不报错，只能靠物理直觉发现，属于血泪经验。

4.4 反弹边界用当前步分布函数导致质量不守恒

现象：总质量随时间缓慢漂移，长时间跑下来密度整体上升或下降。原因：反弹边界应该用碰撞前的分布函数，如果用了碰撞后的，固体节点上的分布函数被松弛修改过，反弹回去的质量不守恒。解决：在碰撞之前保存fprev，反弹时从fprev取反向分布函数。代码里fprev = f的位置必须在迁移之前，且反弹循环里用的是fprev(kk,:,:)而不是f(kk,:,:)。

4.5 收敛判据只看速度不看密度

现象：速度场已经稳定，但密度场还在缓慢变化，继续跑下去速度场又变了。原因：低马赫数下密度变化很小，速度收敛快于密度收敛，只监控速度会过早停止迭代。解决：同时监控密度和速度的最大变化量，取两者较大值作为收敛指标。另外压力驱动问题里总质量会随密度差调整，密度场需要更长时间才能达到稳态，建议把收敛阈值放宽到 1e-5 并增加检查间隔到 500 步。

5. 从单相渗流到相对渗透率曲线：一个可复用的进阶技巧

单相渗流跑通之后，最自然的进阶是两相流，比如水驱油或者气体扩散层里的气液两相。两相 LBM 常见的有 Shan-Chen 伪势模型和颜色梯度模型，前者实现简单但界面宽度大，后者界面锐利但代码复杂。如果目标只是相对渗透率曲线，有个取巧的办法：用单相求解器分别计算不同饱和度下的有效渗透率，饱和度通过随机填充固体节点来模拟，虽然物理上不严格，但趋势和实验数据对得上，适合快速筛选。具体做法是先生成基础掩膜，再按目标饱和度随机把一部分流体节点标记为第二相固体，跑单相求解器，得到该饱和度下的有效渗透率，归一化后就是相对渗透率。

% 用单相求解器扫饱和度，近似相对渗透率曲线 saturations = 0.3:0.1:0.9; K_eff = zeros(size(saturations)); for i = 1:length(saturations) solid_tmp = solid; fluid_idx = find(~solid); n_block = round((1 - saturations(i)) * numel(fluid_idx)); block_idx = fluid_idx(randperm(numel(fluid_idx), n_block)); solid_tmp(block_idx) = true; % 调用单相求解器，返回该掩膜下的渗透率 K_eff(i) = run_single_phase(solid_tmp, tau, rho_in, rho_out); end K_rel = K_eff / K_eff(end); % 以最大饱和度为基准归一化 plot(saturations, K_rel, 'o-'); xlabel('饱和度'); ylabel('相对渗透率');

这个技巧的关键在于randperm的随机种子要和基础掩膜一致，否则不同饱和度下的孔隙结构不可比。run_single_phase就是前面主循环封装成的函数，输入掩膜和边界条件，输出渗透率。扫一条曲线大概需要跑十次单相求解器，每次几千步，Matlab 在 200x100 网格上单次大约几十秒，整条曲线几分钟能出来。如果要做完整的两相驱替模拟，还是得上 Shan-Chen 或者相场模型，但那个代码量和调试成本是另一个量级。我自己的习惯是先用这个近似方法快速看趋势，确认有信号之后再投入两相求解器的开发。另外提醒一句，Matlab 的向量化虽然方便，但三维数组在循环里反复 squeeze 和 reshape 会有内存拷贝开销，网格超过 500x500 之后建议改用单精度或者分块处理，否则内存和速度都会成为瓶颈。希望帮到你。

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