Matlab LBM多孔介质流动仿真：D2Q9模型与边界条件实现

简介：这份资源提供了一套基于 Matlab 的 LBM 流体仿真代码，用于模拟流体流经多孔介质的流动过程，面向计算机、电子信息工程、数学等专业的大学生及科研人员，可服务于课程设计、期末大作业与毕业设计等场景。资源包共 11 个文件，包含 7 个 m 脚本、3 张 png 结果图与 1 份 pdf 说明文档，整体约 3.45MB，脚本覆盖泊肃叶流、方腔流及多孔介质流动等多个算例，图片与文档便于对照验证仿真结果。代码采用参数化编程，参数修改方便，注释详尽，编程思路清晰，附赠案例数据可直接运行，无需从零搭建模型。目前已有 82 人学习，适合希望快速上手 LBM 数值仿真、理解碰撞与迁移迭代流程并完成相关课题实践的读者参考使用。

1. 从一份 Matlab LBM 代码包说起：多孔介质流动仿真到底能跑出什么

如果你做过多孔介质渗流、电池电极孔隙尺度流动或者土壤渗透率估算，大概率绕不开一个尴尬：商业 CFD 软件网格剖分到孔隙级别就崩，算一次要等半天，参数还调不动。这份「使用 Matlab 代码实现 LBM 流体仿真以模拟流经多孔介质的流动」的资源包，走的是另一条路——用格子玻尔兹曼方法（Lattice Boltzmann Method，LBM）在 Matlab 里把多孔介质流动从头搭出来。LBM 的妙处在于它不直接解 Navier-Stokes 方程，而是模拟流体粒子在格点上的碰撞和迁移，边界处理天然适合复杂孔隙几何，代码量比传统有限体积法小一个量级。这份资源适合三类人：想入门 LBM 但被 C++ 劝退的流体方向研究生、需要在 Matlab 里快速验证多孔介质渗透率的研究者、以及想拿一份能跑通的代码当二次开发底座的工程师。下面我按「代码结构 → 核心算法 → 边界处理 → 避坑 → 进阶验证」的顺序，把这份包拆开讲透。

2. 代码包结构与 LBM 核心算法：D2Q9 模型怎么在 Matlab 里落地

2.1 先看清目录：哪些文件是入口，哪些是工具函数

拿到一个 Matlab 仿真包，最怕的就是打开一堆 .m 文件不知道从哪跑。常见做法是先看根目录有没有main.m或run_simulation.m这类入口脚本，再看有没有README或注释头。这份 LBM 多孔介质包的结构一般会包含这几类文件：

入口脚本里通常先定义格子数nx、ny，松弛时间tau，然后循环调用碰撞和迁移。如果你打开发现没有明确入口，找那个包含for t = 1:nt时间循环的文件，那就是主脚本。

2.2 D2Q9 模型：九个方向的速度和权重怎么设

LBM 的核心是分布函数 $f_i(x,t)$，表示在位置 $x$、时间 $t$、沿第 $i$ 个方向运动的粒子数。D2Q9 是二维九方向模型，九个离散速度向量和对应权重是固定的：

% D2Q9 离散速度向量 (c) 和权重 (w) % 方向编号: 1-9，其中 1 是静止，2-5 是轴向，6-9 是对角 ex = [0, 1, 0, -1, 0, 1, -1, -1, 1]; % x 方向分量 ey = [0, 0, 1, 0, -1, 1, 1, -1, -1]; % y 方向分量 w = [4/9, 1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36]; % 权重 % 平衡态分布函数计算 function feq = equilibrium(rho, ux, uy, ex, ey, w) % rho: 宏观密度, ux/uy: 宏观速度分量 % 返回 1x9 的平衡态分布 cu = ex * ux + ey * uy; % 点积 u2 = ux^2 + uy^2; feq = rho .* w .* (1 + 3*cu + 4.5*cu.^2 - 1.5*u2); end

这段代码里ex、ey定义了粒子可能的运动方向，w是各方向的权重，静止方向权重最大（4/9），对角方向最小（1/36）。平衡态函数equilibrium是 LBM 的基石，它把宏观密度和速度映射到九个方向的分布上。参数rho和ux、uy在每个格点都要算，通常从当前分布函数求和得到：rho = sum(f)，ux = sum(f.*ex)/rho，uy = sum(f.*ey)/rho。

2.3 碰撞与迁移：一个时间步里到底发生了什么

LBM 的每个时间步分两步：碰撞（collision）和迁移（streaming）。碰撞让分布函数向平衡态靠拢，迁移让分布函数沿各自方向跑到相邻格点。BGK 碰撞算子是最常用的：

% BGK 碰撞 + 迁移，单时间步 function [f, rho, ux, uy] = lbm_step(f, ex, ey, w, tau, nx, ny) % f: 当前分布函数，尺寸 [nx, ny, 9] % tau: 松弛时间，决定粘度 nu = (tau - 0.5)/3 % 1. 计算宏观量 rho = sum(f, 3); ux = sum(f .* reshape(ex, 1, 1, 9), 3) ./ rho; uy = sum(f .* reshape(ey, 1, 1, 9), 3) ./ rho; % 2. 碰撞：向平衡态松弛 feq = zeros(nx, ny, 9); for i = 1:9 cu = ex(i)*ux + ey(i)*uy; feq(:,:,i) = rho .* w(i) .* (1 + 3*cu + 4.5*cu.^2 - 1.5*(ux.^2+uy.^2)); end f = f - (1/tau) * (f - feq); % BGK 碰撞 % 3. 迁移：沿各方向搬运 for i = 1:9 f(:,:,i) = circshift(f(:,:,i), [ex(i), ey(i)]); end end

碰撞步里f = f - (1/tau)*(f - feq)是 BGK 近似，tau越接近 0.5 粘度越小，但低于 0.5 会数值不稳定。迁移步用circshift实现周期性边界，这是最省事的写法，但多孔介质仿真通常需要更复杂的边界。注意circshift的第二个参数是[ex(i), ey(i)]，对应 x 和 y 方向的位移，正负号决定搬运方向。这段代码跑起来后，每个时间步更新一次f，宏观量rho、ux、uy可以实时可视化。

3. 多孔介质几何生成与边界条件：让流体在孔隙里真正流起来

3.1 多孔介质几何：随机圆球、CT 切片还是规则阵列

多孔介质的几何决定了流动路径。这份资源包里常见的几何生成方式有三种：随机放置圆形颗粒、读取 CT 扫描切片二值化、规则圆柱阵列。随机圆球法最容易上手：

% 生成随机圆形颗粒的多孔介质几何 function solid = generate_porous_media(nx, ny, num_particles, radius_range) % solid: 逻辑矩阵，true 表示固体，false 表示流体 solid = false(nx, ny); for k = 1:num_particles r = radius_range(1) + rand*(radius_range(2)-radius_range(1)); cx = rand * nx; cy = rand * ny; [X, Y] = meshgrid(1:nx, 1:ny); dist = sqrt((X-cx).^2 + (Y-cy).^2); solid(dist < r) = true; end % 确保入口出口不被堵死 solid(1:2, :) = false; solid(end-1:end, :) = false; end

num_particles控制颗粒数量，radius_range控制粒径范围。生成后solid为 true 的格点是固体，流体不能穿过。最后两行强制清空入口出口两行，防止颗粒堵住流动通道。孔隙率可以通过1 - sum(solid(:))/(nx*ny)算出来，一般多孔介质孔隙率在 0.3 到 0.6 之间。如果孔隙率太低，流体几乎流不动，需要减少颗粒数或缩小粒径。

3.2 边界条件：入口速度、出口压力、壁面反弹

多孔介质 LBM 的边界条件比普通管道流复杂。入口通常用速度边界（Zou-He 或 Guo 格式），出口用压力边界，固体壁面用反弹格式（bounce-back）。反弹格式最简单：流体粒子撞到固体就原路返回。

% 固体壁面反弹边界（半步反弹） function f = bounce_back(f, solid, ex, ey) % solid: 逻辑矩阵，true 为固体 % 对每个固体格点，把相反方向的分布函数交换 [nx, ny] = size(solid); for i = 1:9 % 找相反方向 opp = find(ex == -ex(i) & ey == -ey(i)); if isempty(opp), continue; end % 固体格点处，把 i 方向的分布换成 opp 方向 temp = f(:,:,i); f(:,:,i) = f(:,:,opp); f(:,:,opp) = temp; end end

半步反弹的物理含义是：粒子撞到壁面后，沿相反方向弹回，动量守恒。这段代码里opp找的是与i方向相反的编号，然后交换两个方向的分布函数。注意反弹边界要放在迁移之后、碰撞之前，顺序错了结果会发散。入口速度边界常用 Zou-He 格式，根据目标速度反推入口处的分布函数，出口则固定密度（压力）让流体自然流出。

3.3 主循环：把初始化、碰撞、迁移、边界串起来

把前面几块拼起来，主循环大概长这样：

% 主循环框架 nx = 200; ny = 100; % 格子数 tau = 0.8; % 松弛时间 nt = 5000; % 时间步数 u_in = 0.05; % 入口速度（格子单位） % 初始化 f = zeros(nx, ny, 9); rho0 = 1.0; for i = 1:9 f(:,:,i) = rho0 * w(i); % 均匀平衡态初始化 end solid = generate_porous_media(nx, ny, 80, [3, 8]); % 时间循环 for t = 1:nt % 1. 碰撞 [f, rho, ux, uy] = collision(f, ex, ey, w, tau); % 2. 迁移 f = streaming(f, ex, ey); % 3. 边界条件 f = bounce_back(f, solid, ex, ey); f = inlet_velocity(f, u_in, rho0, ex, ey, w); f = outlet_pressure(f, rho0, ex, ey, w); % 4. 每 500 步可视化 if mod(t, 500) == 0 imagesc(sqrt(ux.^2 + uy.^2)'); axis equal; colorbar; title(['Step ', num2str(t)]); drawnow; end end

u_in = 0.05是格子单位速度，实际物理速度要乘以特征速度。tau = 0.8对应粘度nu = (0.8-0.5)/3 ≈ 0.1，这个值比较稳。时间步nt一般要跑到流场稳定，可以通过监测出口流量或总动能判断。可视化用imagesc画速度幅值，axis equal保证格子比例正确。如果跑几百步后速度场出现棋盘振荡，多半是tau太小或边界处理有问题。

4. 避坑与排查：LBM 多孔介质仿真最容易翻车的五个地方

4.1 仿真跑着跑着就发散，密度出现负值

现象：时间步到几百或几千步时，rho出现负数或 NaN，速度场爆炸。原因：tau太接近 0.5，或者入口速度u_in太大，马赫数超过 0.3 后可压缩性误差急剧放大。解决：把tau调到 0.7 以上，u_in降到 0.1 以下。LBM 是弱可压缩求解器，格子单位速度超过 0.3 基本不可信。如果必须高雷诺数，考虑用 MRF 碰撞模型或加密网格。

4.2 多孔介质里流体完全不动，速度全是零

现象：跑了几千步，速度场几乎为零，压力差也建立不起来。原因：入口出口被固体颗粒堵死，或者反弹边界把入口也当成壁面处理了。解决：检查solid矩阵的入口出口两行是否被清空。反弹边界只对固体格点生效，入口出口格点要排除。可以在bounce_back里加判断if solid(i,j)再交换。

4.3 渗透率算出来和文献差一个量级

现象：用达西定律反算渗透率，结果和参考文献差 10 倍以上。原因：格子单位没换算到物理单位，或者孔隙率、粒径统计方式不一致。解决：LBM 的渗透率在格子单位下是K_lattice = nu * u_avg / grad_p，换算到物理单位要乘特征长度的平方。另外确认孔隙率定义是面积比还是体积比，二维仿真通常用面积比。粒径分布也要和文献对齐，随机圆球的粒径范围不同，渗透率差几倍很正常。

4.4 边界处出现非物理的密度波

现象：入口或出口附近密度波动很大，速度场出现回流。原因：Zou-He 速度边界和压力边界不匹配，或者出口太靠近多孔介质区域。解决：在入口出口和多孔介质之间留一段缓冲区域，一般 10 到 20 个格子。出口压力边界用外推格式比固定密度更稳。如果还不行，把入口速度再降低，或者改用 Guo 格式的力项驱动流动。

4.5 Matlab 循环太慢，跑一次要几个小时

现象：for t = 1:nt循环里每一步都要几秒，nt = 10000根本跑不完。原因：Matlab 的逐元素循环效率低，尤其是circshift和for i = 1:9这种嵌套。解决：把九个方向的分布函数向量化，用reshape和permute替代循环。迁移步可以用矩阵索引一次性完成，不用circshift。另外把nx、ny控制在 200 以内，超过 500 建议转 C++ 或 GPU。Matlab 的parfor对 LBM 帮助有限，因为每个时间步有依赖。

5. 进阶验证：用达西定律反算渗透率，判断结果能不能信

跑完仿真只是第一步，关键是要验证结果对不对。多孔介质 LBM 最直接的验证方法是达西定律：在低雷诺数下，流量和压力梯度成正比，比例系数就是渗透率。具体操作是跑几组不同压力梯度，拟合流量-压力梯度曲线，看线性度好不好。

% 达西定律验证：多组压力梯度算渗透率 grad_p_list = [1e-5, 2e-5, 5e-5, 1e-4]; % 不同压力梯度 u_avg_list = zeros(size(grad_p_list)); for k = 1:length(grad_p_list) % 用对应压力梯度跑仿真，取稳定后平均速度 u_avg_list(k) = run_simulation(grad_p_list(k)); end % 线性拟合 p = polyfit(grad_p_list, u_avg_list, 1); K_lattice = nu * p(1); % 格子单位渗透率 fprintf('格子单位渗透率: %.6f\n', K_lattice); % 换算到物理单位（假设特征长度 L_ref） L_ref = 1e-3; % 1 mm K_physical = K_lattice * L_ref^2; fprintf('物理渗透率: %.3e m^2\n', K_physical);

这段代码里grad_p_list是人为施加的压力梯度，run_simulation返回稳定后的平均速度。polyfit拟合斜率，乘以粘度nu得到格子单位渗透率。换算到物理单位时，L_ref是格子对应的物理长度，比如一个格子代表 1 毫米，那L_ref = 1e-3米。如果拟合的 R² 低于 0.95，说明雷诺数太高，达西定律不适用，需要降低压力梯度重跑。

另一个验证角度是看速度场是否合理：孔隙内部速度应该比固体附近高，回流区应该出现在颗粒背面，整体流线应该绕过固体而不是穿过。如果速度场看起来像随机噪声，多半是边界条件或碰撞参数有问题。我一般会先跑一个单圆柱绕流的算例，确认代码能复现卡门涡街，再上多孔介质。从那以后我每次拿到新的 LBM 代码，都强制走一遍「单圆柱验证 → 孔隙率扫描 → 达西拟合」三步，少一步都不敢信结果。希望这份拆解能帮到你，代码包里的参数和边界按自己的几何改一改就能跑起来。

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