简介：这份资源提供了一套基于 Matlab 实现的格子玻尔兹曼方法（LBM）流体仿真代码，专门用于模拟流体在多孔介质中的流动过程，面向计算机、电子信息工程、数学等专业的大学生及科研人员，可服务于课程设计、期末大作业与毕业设计等场景。压缩包共包含 11 个文件，以 7 个 m 脚本文件为核心，辅以 3 张 png 结果图与 1 份 pdf 说明文档，整体约 3.45MB，结构紧凑、便于直接运行与二次开发。代码采用参数化编程思路，关键参数可灵活调整，注释详尽、逻辑清晰，覆盖泊肃叶流、方腔流及多孔介质流动等多个算例，并附赠案例数据，用户无需从零建模即可快速获得仿真结果。目前已有 82 人学习关注，适合希望深入理解 LBM 碰撞与迁移机制、掌握数值仿真技巧并完成相关课题的读者参考使用。

1. 从一张砂砾切片说起：LBM 为什么适合啃多孔介质这块硬骨头

多孔介质里的流动，最反直觉的一点是：你明明知道达西定律能算平均流量，却很难说清流体在孔隙拐角处到底怎么拐、哪里会滞留、哪里会优先冲出一条通道。岩石、土壤、电池电极、催化载体、纤维滤芯，这些结构的共同点是固体骨架随机、孔隙尺度小、界面面积大。用传统 CFD 贴体网格去追每一个孔隙，前处理就能把人耗死。格子玻尔兹曼方法（LBM）换了个思路：不直接解 Navier-Stokes，而是让粒子分布函数在规则格子上碰撞、迁移，边界条件用反弹格式处理复杂固体壁面，天然适合多孔介质这种几何。Matlab 虽然不是 LBM 的主流生产工具，但它的矩阵化写法、可视化能力和调试便利性，让它成为理解算法、验证小规模算例、快速出图的合适载体。这篇笔记面向想用 Matlab 把 LBM 跑通、并真正模拟出多孔介质渗流的人：从 D2Q9 模型怎么摆、多孔介质几何怎么生成、边界怎么设，到代码结构、参数标定和踩坑记录，一步步落到能复现的程度。

2. D2Q9 与多孔介质：把物理模型拆成可写的矩阵

2.1 为什么选 D2Q9 而不是 D3Q19 起步

LBM 的维度与速度方向数用 DnQm 表示，n 是空间维度，m 是离散速度方向数。二维问题最常用 D2Q9，九个方向：一个静止方向、四个轴向、四个对角。选它起步的理由很实际：每个格点的分布函数只有 9 个分量，内存和计算量可控，Matlab 里用三维数组f(Ny, Nx, 9)就能装下；边界反弹、周期、压力出口这些格式在 D2Q9 上资料最全，出错时容易对照。三维 D3Q19 当然更贴近真实多孔介质，但格点数一上去，Matlab 的循环版本会慢到让人怀疑人生，除非做向量化或 mex 加速。我的建议是：先用 D2Q9 把碰撞迁移、边界、多孔介质生成、渗透率统计这条链路走通，再考虑升维。

D2Q9 的九个速度向量和权重是固定的：

权重之和为 1，这是保证宏观密度能从分布函数零阶矩还原的前提。速度向量的排列顺序会影响后续迁移步的索引写法，建议一开始就固定成上表顺序，别中途换。

2.2 多孔介质几何：随机圆颗粒与四参数控制

多孔介质在 LBM 里通常表示为一个布尔矩阵solid(Ny, Nx)，true 表示固体，false 表示流体。生成方式有很多：随机圆颗粒、随机矩形、分形、CT 图像二值化。最容易复现且参数直观的是随机圆颗粒法，用四个参数控制：孔隙率目标phi、颗粒半径范围[rmin, rmax]、颗粒数量上限Nmax、最小间距gap。核心逻辑是反复随机撒圆，若新圆与已有圆重叠或间距不足就丢弃，直到达到目标孔隙率或尝试次数上限。

function solid = gen_porous(Ny, Nx, phi, rmin, rmax, Nmax, gap) % 生成随机圆颗粒多孔介质，返回逻辑矩阵 solid solid = false(Ny, Nx); [xx, yy] = meshgrid(1:Nx, 1:Ny); target_solid = (1 - phi) * Ny * Nx; % 目标固体格点数 count = 0; tries = 0; while count < target_solid && tries < Nmax tries = tries + 1; r = rmin + (rmax - rmin) * rand(); cx = r + 1 + (Nx - 2*r - 2) * rand(); cy = r + 1 + (Ny - 2*r - 2) * rand(); mask = (xx - cx).^2 + (yy - cy).^2 <= (r + gap)^2; if any(solid(mask)) % 与已有颗粒冲突 continue; end solid(mask) = true; count = sum(solid(:)); end end

这段代码里gap是颗粒间最小流体间隙，太小会导致孔隙连通性差、渗透率极低，太大则孔隙率控制不准。Nmax是尝试次数上限，防止高固体含量时死循环。实际调参时，phi取 0.5 到 0.8 比较常见，rmin建议不小于 3 个格子，否则边界反弹的离散误差会很明显。生成后一定要用连通域检查流体区域是否贯通，否则算出来的渗透率没有意义。

2.3 从分布函数到宏观量：密度、速度、压力的还原

每个时间步，先由分布函数求宏观密度和速度：

rho = sum(f, 3); % 密度 = 零阶矩 ux = sum(f .* reshape(cx, 1, 1, 9), 3) ./ rho; uy = sum(f .* reshape(cy, 1, 1, 9), 3) ./ rho;

压力用状态方程p = rho * cs^2，其中cs^2 = 1/3（格子单位）。这里有个容易翻车的点：Matlab 的sum(f, 3)对f(Ny, Nx, 9)是按第三维求和，结果维度是Ny×Nx，没问题；但如果你把f存成(9, Ny, Nx)，求和维度就要改成 1，速度点乘也要相应调整。两种存储方式都有人用，关键是全篇统一，别在迁移步里混用。

3. 碰撞与迁移：Matlab 里怎么写才不慢又不乱

3.1 BGK 碰撞算子与松弛时间的关系

LBM 的演化分两步：碰撞和迁移。碰撞最常用 BGK 近似：

function f = collide(f, rho, ux, uy, tau, w, cx, cy) % BGK 碰撞：f = f - (f - feq)/tau feq = zeros(size(f)); for k = 1:9 cu = cx(k)*ux + cy(k)*uy; feq(:,:,k) = w(k) * rho .* (1 + 3*cu + 4.5*cu.^2 - 1.5*(ux.^2 + uy.^2)); end f = f - (f - feq) / tau; end

tau是松弛时间，和运动黏度nu = (tau - 0.5) / 3（格子单位）。tau必须大于 0.5，否则黏度为负，计算直接发散。工程上tau取 0.6 到 1.0 比较稳，太接近 0.5 会数值不稳定，太大则黏性过强、雷诺数上不去。平衡态分布feq的展开式里系数 3、4.5、1.5 来自 D2Q9 的声速平方cs^2 = 1/3，写错一个系数，宏观方程就还原不出 Navier-Stokes。

3.2 迁移步的索引陷阱与向量化写法

迁移就是每个格点的分布函数沿对应方向搬到邻居。最直观的循环写法：

for k = 1:9 f(:, :, k) = circshift(f(:, :, k), [cy(k), cx(k)]); end

circshift默认周期边界，对多孔介质内部没问题，但进出口需要额外处理。注意circshift的第二个参数是[行偏移, 列偏移]，对应[y, x]，而速度向量cx是 x 方向、cy是 y 方向，顺序别写反，否则流场会整体转 90 度，这种错误肉眼看不出来，但渗透率会完全不对。向量化写法比逐点循环快一个数量级，Matlab 里优先用circshift。

3.3 多孔介质边界：反弹格式的三种写法

固体壁面用反弹格式：流体粒子撞到固体，速度反向。标准半步反弹（half-way bounce-back）在迁移后执行：

for k = 1:9 kk = opp(k); % 相反方向索引 tmp = f(:, :, k); tmp(solid) = f(solid + (kk-1)*Ny*Nx); % 简化示意，实际用逻辑索引 f(:, :, k) = tmp; end

实际写的时候用逻辑索引更清晰：先保存碰撞后的f，再对固体格点把f(:,:,k)替换成邻居的f(:,:,opp(k))。opp是相反方向映射表，比如 1↔3、2↔4、5↔7、6↔8。反弹格式的精度和边界位置有关，半步反弹把壁面放在格点中间，比整步反弹更准。多孔介质里固体格点很多，反弹步是性能热点，尽量向量化，别用双层循环。

3.4 进出口边界：压力驱动与速度驱动的选择

模拟多孔介质渗流，通常用压力驱动：左端固定密度rho_in，右端固定密度rho_out，密度差产生流动。实现方式是在进出口各加一层虚拟格点，每步把平衡态分布赋给虚拟层：

% 左端入口，给定 rho_in，速度由外推或零梯度估计 feq_in = compute_feq(rho_in, ux_in, uy_in, w, cx, cy); f(:, 1, :) = feq_in; % 右端出口，给定 rho_out feq_out = compute_feq(rho_out, ux_out, uy_out, w, cx, cy); f(:, end, :) = feq_out;

入口速度ux_in一般设为零或由上一时间步外推，避免入口反射。密度差drho = rho_in - rho_out控制驱动力，太大（比如超过 0.1）会导致可压缩性误差，太小则流动太弱、统计噪声大。常见做法是rho_in = 1.0 + drho/2，rho_out = 1.0 - drho/2，drho取 0.001 到 0.01。

4. 把代码跑起来：主循环、参数标定与渗透率统计

4.1 主循环骨架与收敛判据

主循环按「宏观量 → 碰撞 → 迁移 → 边界 → 收敛检查」顺序：

for step = 1:maxStep rho = sum(f, 3); ux = sum(f .* reshape(cx,1,1,9), 3) ./ rho; uy = sum(f .* reshape(cy,1,1,9), 3) ./ rho; f = collide(f, rho, ux, uy, tau, w, cx, cy); f = stream(f, cx, cy); f = apply_bounce(f, solid, opp); f = apply_inlet_outlet(f, rho_in, rho_out, w, cx, cy); if mod(step, 100) == 0 res = norm(ux(:) - ux_old(:)) / norm(ux(:)); if res < 1e-6, break; end ux_old = ux; end end

收敛判据用速度场相对变化，阈值 1e-6 到 1e-8。多孔介质里速度场局部波动大，建议每 100 步检查一次，别每步都算范数，否则统计开销比计算还大。

4.2 格子单位与物理单位的换算

LBM 算出来的是格子单位，要对应到真实物理量，需要选特征长度、特征速度、特征密度。常见做法：取颗粒直径或孔隙特征尺寸为特征长度L0，取真实渗流速度或声速为特征速度U0，则时间尺度T0 = L0/U0，黏度换算nu_phys = nu_lattice * L0 * U0。雷诺数Re = U0 * L0 / nu_phys必须和格子雷诺数一致。这一步最容易出错的是把格子声速当成物理声速，记住cs = 1/sqrt(3)只是格子单位下的数值，不是真实声速。

4.3 渗透率计算：达西定律反推与统计稳态

渗透率K由达西定律反推：K = -nu * <u> / (dp/dx)，其中<u>是流体区域平均速度，dp/dx是压力梯度。压力用p = rho/3，进出口密度差对应压力差。统计时只对流体格点求平均，固体格点排除：

fluid_mask = ~solid; u_avg = mean(ux(fluid_mask)); dpdx = (rho_out - rho_in) / 3 / Nx; K = -nu * u_avg / dpdx;

注意符号：如果rho_in > rho_out，压力梯度为负，流速为正，K为正。统计要在收敛后进行，最好再平均若干时间步，减少波动。多孔介质渗透率对几何随机性敏感，同一组参数生成多个样本取平均，才能得到有代表性的值。

5. 避坑与排查：那些让结果悄悄错掉的细节

5.1 现象：算了几千步速度场还是零

原因通常是驱动力没加上，或者进出口边界被反弹步覆盖。检查顺序：先看rho_in和rho_out是否真的不同，再看进出口赋值是否在反弹之后执行。如果反弹步把进出口虚拟层也当成固体处理，密度差就传不进去。解决：进出口层不参与反弹，或者在反弹后用平衡态强制覆盖。

5.2 现象：计算几步后密度爆炸到 NaN

原因多半是tau太接近 0.5，或者初始密度设得太大。tau = 0.5时黏度为零，数值不稳定。解决：tau至少取 0.6，初始密度全场设 1.0，进出口密度差先从小值（0.001）试起。如果还炸，检查平衡态展开式系数是否写错，尤其是cu.^2的系数 4.5 和ux.^2+uy.^2的系数 1.5。

5.3 现象：渗透率比文献值高一个数量级

原因通常是固体边界反弹没生效，或者固体矩阵和分布函数索引错位。检查solid矩阵的true位置是否和f的格点一一对应，反弹步是否真的替换了固体格点的分布函数。另一个常见原因是孔隙率统计时把边界虚拟层也算进流体，导致平均速度偏高。解决：统计时排除进出口虚拟层，只统计内部流体格点。

5.4 现象：流场出现棋盘格振荡

原因通常是碰撞和迁移的耦合顺序不对，或者circshift方向写反导致信息传播错乱。标准顺序是先碰撞后迁移，迁移用circshift时偏移量是[cy(k), cx(k)]。如果顺序反了，或者偏移量写成[cx(k), cy(k)]，流场会出现非物理振荡。解决：固定顺序，用一个小算例（比如单圆绕流）验证流线是否对称。

5.5 现象：Matlab 跑得比预期慢很多

原因通常是碰撞和反弹用了逐点循环，或者每步都做可视化。Matlab 的循环开销大，碰撞和反弹尽量向量化。可视化每几百步做一次，别每步imagesc。如果还是慢，考虑把核心循环写成 mex，或者用gpuArray加速，但先确认算法正确再优化。

6. 进阶技巧：用连通域和迂曲度验证你的多孔介质

跑通基本流程后，真正决定结果可信度的是多孔介质几何的物理合理性。两个我常用的验证指标：流体连通域和迂曲度。连通域用bwlabel检查流体区域是否贯通，如果最大连通域占比远小于 1，说明孔隙被固体隔断，渗透率会偏低甚至为零。迂曲度用流线路径长度除以直线距离估计，反映孔隙通道的弯曲程度，迂曲度越大，相同孔隙率下渗透率越低。

fluid = ~solid; [L, num] = bwlabel(fluid, 4); % 4 连通 sizes = histcounts(L, 1:num+1); max_frac = max(sizes) / sum(fluid(:)); if max_frac < 0.95 warning('流体连通域占比仅 %.2f，孔隙可能不贯通', max_frac); end

bwlabel的第二个参数 4 表示四连通，多孔介质里用四连通更严格，能发现对角缝隙造成的假连通。如果最大连通域占比低于 0.95，建议重新生成几何或调整gap。另一个习惯是：每次改参数后，先跑一个低分辨率算例看流线是否合理，再上高分辨率。我踩过最深的坑是几何生成后没检查连通性，算了一晚上收敛曲线很漂亮，结果渗透率比预期低两个数量级，回头一看流体区域被一圈固体围成了孤岛。希望这些能帮到你。

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