LBM多孔介质流动MATLAB仿真：D2Q9与D3Q19实现与渗透率标定

简介：面向LBM（格子玻尔兹曼方法）初学者与相关专业课程设计人群，内容是国外经典的D2Q9与D3Q19模型Matlab实现。代码采用参数化编程，清晰标注注释，可方便修改参数并直接运行，适合用于本科毕设、课程作业中多孔介质流动等流体仿真场景。压缩包共10个文件，以5个m主程序文件为核心，辅以3张结果示意图便于对照，另含参考rar数据与asv自动保存文件，整体大小仅78KB，小巧易用。已有151人学习，代码结构紧凑，涵盖二维与三维LBM建模、边界处理及数据输出等完整流程，可帮助读者从零理解经典模型的核心实现。作者为资深算法工程师，包内附带可直接运行的案例数据，遇到问题可通过私信获取支持。

1. LBM D2Q9与D3Q19的matlab仿真包解决什么问题

解压“国外经典LBM D2Q9 D3Q19模型的matlab代码.zip”，最先看到的是Porous_Medium_Flows.m、Porous_Medium_Flows_3D.m、LBMval.m和几个save开头的函数。文件名已经交代了这套代码的边界：二维用D2Q9速度模型做多孔介质流动，三维用D3Q19做同一件事。LBM（格子玻尔兹曼方法）把流体离散成粒子分布函数在一个格子上的碰撞与迁移，不需要像传统CFD那样驱动压力泊松方程迭代，特别适合孔隙结构复杂、边界不规则的场合，这也是它在多孔介质、微流控、渗流计算里比有限体积法更吃香的原因。

这套代码的价值在于它不是孤立的算法演示，而是一套能直接跑出压力场、速度场和VTK后处理文件的教学框架。压缩包里附了3.png、4.png等运行结果图，LBMval.m承担二维BGK的迭代核心，savevtk.m和savevtkvector.m负责把结果导出成ParaView能读的格式，整套数据流是闭合的。适合正在做多孔介质流动模拟、渗透率标定、LBM课程设计或毕业设计的人，也适合想快速读懂BGK落地细节的工程师。如果MATLAB环境没配好，按matlab下载安装教程装好2021a即可直接运行，代码本身兼容2014以上版本。

2. D2Q9离散速度模型与LBM碰撞迁移的matlab实现

2.1 九速格子与平衡态分布函数

LBM的基本变量不是速度场，而是分布函数 f_i(x,t)，下标 i 表示离散速度方向。D2Q9的“9”指二维平面上取9个方向：静止、4个正交方向、4个对角方向。每个方向有各自的权重 w_i，权重决定了宏观流动恢复Navier-Stokes方程时的误差阶数。

宏观密度 ρ 和宏观速度 u 由分布函数的零阶矩和一阶矩恢复：

ρ = Σ f_i，ρu = Σ f_i · c_i

BGK模型假设粒子分布函数以特征时间 τ 松弛到平衡态 f_i^eq，平衡态分布写成：

f_i^eq = ρ·w_i·(1 + 3(c_i·u) + 4.5(c_i·u)² - 1.5(u·u))

这里的系数3、4.5、1.5来自声速 c_s = 1/√3 的格子单位设定，不是经验调参，而是从Chapman-Enskog展开中恢复不可压缩Navier-Stokes方程的必要条件。多孔介质流动的Re数通常远小于1，BGK的误差阶能覆盖这一范围，这也是教学代码普遍选择BGK而不是MRT的原因。

2.2 碰撞-迁移主循环：LBMval.m的函数封装

我复现LBMval.m的典型写法，把9个方向的速度向量、权重和反向索引都放在函数开头，避免每个时间步重复构造：

function [rho, ux, uy, f] = LBMval(nx, ny, tau, u0, obstacle, nsteps) % D2Q9单松弛BGK求解器 % nx, ny: 网格尺寸; tau: 松弛时间; u0: 初始特征速度; obstacle: 固相掩模 cx = [0 1 0 -1 0 1 -1 -1 1]; cy = [0 0 1 0 -1 1 1 -1 -1]; opp = [0 3 4 1 2 7 8 5 6]; w = [4/9 1/9 1/9 1/9 1/9 1/36 1/36 1/36 1/36]; rho = ones(nx, ny); ux = u0 * ones(nx, ny); uy = zeros(nx, ny); f = zeros(nx, ny, 9); % 按平衡态初始化分布函数 for k = 1:9 cu = 3 * (cx(k)*ux + cy(k)*uy); f(:,:,k) = rho .* w(k) .* (1 + cu + 0.5*cu.^2 - 1.5*(ux.^2 + uy.^2)); end for t = 1:nsteps % 宏观量更新：零阶矩和一阶矩 rho = sum(f, 3); ux = sum(bsxfun(@times, f, reshape(cx,1,1,9)), 3) ./ rho; uy = sum(bsxfun(@times, f, reshape(cy,1,1,9)), 3) ./ rho; % BGK碰撞：向平衡态松弛 for k = 1:9 cu = 3 * (cx(k)*ux + cy(k)*uy); feq = rho .* w(k) .* (1 + cu + 0.5*cu.^2 - 1.5*(ux.^2 + uy.^2)); f(:,:,k) = f(:,:,k) - (f(:,:,k) - feq) / tau; end % 迁移：circshift实现周期性移位 for k = 1:9 f(:,:,k) = circshift(f(:,:,k), [cy(k), cx(k)]); end % 反弹边界：对固相格点做反向修正 if nargin > 4 f_old = f; for k = 1:9 f(:,:,k) = f(:,:,k) .* (1-obstacle) + f_old(:,:,opp(k)) .* obstacle; end end end end

这段代码严格按“宏观量更新→碰撞→迁移→边界”四个阶段组织。bsxfun在这里的作用是把9×1的方向向量广播到每层分布函数对应的空间维度，2016a之后的MATLAB可以直接用隐式扩展，但bsxfun能保证在旧版本上不报错，符合这个压缩包标注的2014/2019a/2021a兼容范围。

circshift(f(:,:,k), [cy(k), cx(k)])这一句是对应关系最容易写错的地方。MATLAB数组第一维是行方向，对应物理上的y轴；第二维是列方向，对应物理上的x轴。所以当方向矢量为(cx=1, cy=0)时，迁移要沿列方向右移一格，circshift第二参数[0, 1]；而cy=1表示沿行方向移动，写成[-1, 0]或[1, 0]取决于坐标系方向约定。源码里如果跑出来的流场整体反向，优先检查这里，这是所有LBM移植新手都会踩的坑。

2.3 参数化编程：tau、u0和网格尺寸怎么改

LBMval.m把物理参数全部提到函数入口，这是典型的参数化编程思路。运行Porous_Medium_Flows.m时只需要调整以下几组：

• nx、ny：网格分辨率。想复现3.png的流场需要保持与原图一致的网格数，我自己会把nx取120、ny取60，多孔介质区域放在中间，进出口留出充分发展段。
• tau：松弛时间，决定粘度。tau=0.6对应格子粘度nu=(tau-0.5)/3≈0.033，这个值在低速流动下兼顾稳定性和收敛速度。
• u0：初始特征速度。u0=0.05是安全值，超过0.1容易在障碍物锐角处出现局部负密度。
• nsteps：迭代步数。判断收敛最直接的方式是监控rho场的变化量，而不是固定迭代步数。

这套代码的注释密度较高，变量命名基本能对上公式符号，适合直接改造。唯一的性能隐患是碰撞段的9次循环，网格小于200×200时无感，超过500×500建议把循环展开成向量运算，否则单步耗时会有明显增加。

3. 反弹边界与多孔介质障碍物建模：Porous_Medium_Flows.m拆解

3.1 文件构成与障碍物掩模生成

压缩包里同时出现Porous_Medium_Flows.m、Porous_Medium_Flows.asv和Porous_Medium_Flows.rar，这三份本质是同一套代码。.m是主脚本，.asv是MATLAB编辑器在非正常关闭时自动保存的备份版本，内容和.m几乎一致，可以直接忽略；.rar则是打包拷贝，用于网盘分发。冷启动只需要打开Porous_Medium_Flows.m。

多孔介质在LBM里不需要像有限元那样贴体网格划分，只需要把固体格点标记为1、流体格点标记为0，生成一张布尔掩模即可。掩模的分辨率直接决定了孔隙结构的精度，这也是LBM处理复杂多孔介质的最大优势——孔隙形状信息只存在于掩模里，与求解器完全解耦。

构造障碍物的常见做法有三种，按计算量排序：

Porous_Medium_Flows.m里大概率用的是前两种，因为随机圆盘生成简单，且方便控制孔隙率。下面是一段我常用的随机圆盘掩模生成代码：

% 随机圆盘多孔介质，保证圆盘不越过边界 rng(42); % 固定随机种子，保证结果可复现 numDisc = 60; obstacle = false(nx, ny); for d = 1:numDisc rr = 2 + 5 * rand(1); % 圆盘半径2~7 cx0 = rr + 1 + (nx - 2*rr - 2) * rand(1); % 圆心x留出边界 cy0 = rr + 1 + (ny - 2*rr - 2) * rand(1); % 圆心y留出边界 [X, Y] = meshgrid(1:nx, 1:ny); obstacle = obstacle | ((X - cx0).^2 + (Y - cy0).^2 <= rr^2); end

rng(42)保证每次运行生成的孔隙分布完全一致，便于对比不同tau下的流场差异，这是实验类代码必须养成的习惯。圆心坐标减去半径再留1格余量，可以保证圆盘不会贴死边界，否则压差驱动的主通道可能被截断。遍历生成后用|做逻辑或合并，得到的obstacle矩阵直接传给LBMval.m。

3.2 反弹边界算子的向量化实现

多孔介质流动仿真中，流体与固体格点之间的相互作用靠边界条件实现。LBM最常用的“标准反弹格式”（halfway bounce-back）假设边界位于固体格点与流体格点之间的中点，分布函数沿原路反弹。它的离散表达是：

f_i(x, t+1) = f_opp(i)(x, t)

其中opp(i)是i的反方向，对应关系在LBMval.m里已经预定义。这个格式虽然只有一阶精度，但它天然保证质量守恒，对模拟渗透率这类全局宏观量，稳定性和守恒性比单点精度更重要。

在实现上，我一般把迁移和反弹合并成一步处理，避免额外循环：

% 标准反弹格式：与迁移同步完成 fluid = ~obstacle; % 流体掩模 f_old = f; % 保留迁移前分布 for k = 1:9 f(:,:,k) = circshift(f(:,:,k), [cy(k), cx(k)]); end for k = 1:9 f(:,:,k) = f(:,:,k) .* fluid + f_old(:,:,opp(k)) .* obstacle; end

这段代码的执行顺序是：先完成全部格点的周期性迁移，再把固体格点上的分布函数替换为迁移前的反向分布。用fluid和obstacle两个逻辑矩阵做掩模，比用find索引固体格点逐一赋值快几个数量级，因为LBM在MATLAB中的性能瓶颈就是循环中的数组索引开销。注意这里必须用f_old，因为f(:,:,k)在这一轮已经被覆盖，如果直接用f做反向取值会把迁移后的值反弹回去，流场会立刻发散。

3.3 主程序运行流程与收敛判据

Porous_Medium_Flows.m作为驱动脚本，严格控制了求解器的输入输出。一个可直接运行的版本大概长这样：

nx = 120; ny = 60; tau = 0.6; u0 = 0.05; nsteps = 2000; obstacle = generate_porous_medium(nx, ny); % 上文随机圆盘函数 [rho, ux, uy] = LBMval(nx, ny, tau, u0, obstacle, nsteps);

跑完一次大约需要几十秒，取决于网格规模和MATLAB版本。收敛判据建议写在驱动脚本里而不是求解器内部，方便观察中间状态：

% 每隔100步检测一次rho场的最大变化量 err = max(abs(rho(:) - rho_old(:))); if err < 1e-9 break; end

rho是压力在LBM中的投影（实际压力p = rho·c_s²），所以rho场的收敛等价于压力场收敛。多孔介质流动从初始均匀场到达定常态，通常需要几百到几千步，收敛越慢的地方往往是圆盘间隙最小的喉道区域，这些区域速度大、压力梯度大，对松弛时间tau的扰动最敏感。

4. 从D2Q9扩展到D3Q19：三维模型实现与内存约束

4.1 十九速方向表与权重

把多孔介质流动从二维推广到三维，速度模型从D2Q9换成D3Q19。D3Q19包含1个静止方向、6个轴向方向、12个面对角方向，比D3Q27少8个立体对角方向，是教学和工程中最常见的三维LBM选择——少8个方向意味着每个格点少存储8个浮点数，在三维网格下这个差距会放大成数GB的内存差异，而精度损失对低Re多孔介质流动可以忽略。

c = [0 1 -1 0 0 0 0 1 -1 1 -1 1 -1 1 -1 0 0 0 0; 0 0 0 1 -1 0 0 1 1 -1 -1 0 0 0 0 1 -1 1 -1; 0 0 0 0 0 1 -1 0 0 0 0 1 1 -1 -1 1 1 -1 -1]; w = [1/3, repmat(1/18,1,6), repmat(1/36,1,12)];

w的前三个分组对应静止、轴向、面对角三类权重。与D2Q9的权重体系相对照：D2Q9是4/9、1/9、1/36，D3Q19是1/3、1/18、1/36，系数差异来自三维空间下各方向角度的立体角占比。迁移的方向索引同样要按c矩阵的维度顺序写，否则三维方向表中任何一位错位都会让流场整体紊乱。

4.2 三维主循环的改动要点

把LBMval.m升维到三维，修改集中在四个地方：

• 分布函数从f(nx, ny, 9)变成f(nx, ny, nz, 19)。
• circshift的位移参数从[cy, cx]变成[cz, cy, cx]，对应MATLAB三维数组的行、列、页三个维度。
• 宏观量更新从第3维求和改为第4维求和。
• 反弹边界的掩模从二维矩阵变成三维逻辑数组。

核心迭代代码只摘一处：

% 三维迁移：cz对应第一维，cy对应第二维，cx对应第三维 for k = 1:19 f(:,:,:,k) = circshift(f(:,:,:,k), [c(3,k), c(2,k), c(1,k)]); end % 宏观量统计 rho = sum(f, 4); ux = sum(bsxfun(@times, f, reshape(c(1,:),1,1,1,19)), 4) ./ rho; uy = sum(bsxfun(@times, f, reshape(c(2,:),1,1,1,19)), 4) ./ rho; uz = sum(bsxfun(@times, f, reshape(c(3,:),1,1,1,19)), 4) ./ rho;

第4维求和的写法比二维多了一个维度，但逻辑完全一样：零阶矩得到rho，一阶矩除以rho得到速度分量。代码里reshape把方向分量扩成四维，第四维长度19对应19个方向，bsxfun会把速度分量广播到每个格点。

4.3 三维网格规模与内存边界

D3Q19的三维模型最大的瓶颈不是计算量而是内存。分布函数有19个分量，每个分量都是三维数组，double精度下每个格点要占152字节。做一组快速估算：

Porous_Medium_Flows_3D.m里如果直接在循环里写zeros(nx, ny, nz, 19)，默认double类型，256³网格还没开始迭代就已经占满一半内存。改成single是白捡一倍容量：

rho = single(rho); f = single(f);

代价是计算精度降低到约7位有效数字，对渗透率标定这种宏观量完全够用。注意混合使用的时候要保证所有数组都转成single，否则MATLAB会做一次隐式类型转换，吃掉性能红利。三维多孔介质算例建议从64³起步，先跑通圆柱阵列再换随机结构，直接上256³的随机孔隙结构大概率会碰到内存不足或超长时间迭代。

5. LBM数值稳定性、参数区间与典型排错

5.1 松弛时间与格子粘度的换算约束

LBM参数不能随便给，tau的取值直接决定了数值格式的稳定性上限。格子单位下的运动粘度为：

nu = (tau - 0.5) / 3

tau必须大于0.5，等于0.5时粘度为零，数值上表现出无粘特征，迭代必然发散。实际使用中tau的常见安全区间是0.55到0.8，对应nu在0.016到0.1之间。tau越接近0.5，粘度越小，流动越容易在障碍物锐角处出现振荡；tau过大则流动过于耗散，孔隙喉道处的速度梯度被抹平，渗透率被低估。

% 由物理粘度换算格子tau的常见做法 % nu_p: 物理运动粘度, U_p: 特征速度, N: 特征网格数 nu_lb = nu_p * N / (L_p * U_p); % 格子粘度 tau_lb = 3 * nu_lb + 0.5;

这里有一个隐藏约束：tau换算出来如果小于0.5，意味着网格分辨率不够跑这个粘度的流动，只能加网格或降低特征速度，而不是强行把tau设为0.5附近。我见过很多初学者在这个环节直接把tau硬改成0.6，结果就是Re数完全对不上，和实验数据没法比较。

5.2 特征速度上限与低马赫数约束

LBM恢复Navier-Stokes方程依赖小马赫数展开，特征速度u0在格子单位下必须远小于声速c_s=0.577。工程经验是把特征速度控制在0.05～0.1之间，超过0.2时数值误差以马赫数平方的速度增长，密度场会出现明显的压缩伪影。

% 检查当前速度场是否越界 umax = max(sqrt(ux.^2 + uy.^2), [], 'all'); if umax > 0.2 warning('特征速度偏高，建议降低u0或增加网格数'); end

注意这里的“速度”是格子速度，不是物理速度。多孔介质流动中，平均速度很低，但在孔隙喉道处会局部加速到平均速度的几倍，所以即使用u0=0.05启动，喉道处的局部速度也可能超过0.15。这也是为什么LBM多孔介质模拟里，网格分辨率不能太低——太粗的网格会把喉道速度放大到失稳边界。

5.3 失败现象对照表与诊断流程

最隐蔽的问题是“障碍物堵死主通道”。随机圆盘生成时如果圆心和半径分布不当，可能在中间形成一条横向贯穿的固体带，流体只能走上下绕行，宏观渗透率骤降几个数量级。用matlab图像处理工具箱里的bwconncomp可以直接验证入口到出口的连通性：

cc = bwconncomp(fluid); % 统计包含第2列的连通域是否同时包含第nx-1列 entryMask = any(fluid(:,1,:), [1 3]); exitMask = any(fluid(:,nx,:), [1 3]); labels = labelmatrix(cc); valid = false; for comp = 1:cc.NumObjects if entryMask(comp) && exitMask(comp) valid = true; break; end end if ~valid error('入口到出口不存在流体通路，请重新生成多孔介质'); end

bwconncomp把连通的流体区域打上同一组标签，然后检查同一连通域是否同时覆盖入口和出口区域。这个检查放在模拟之前执行，能省下好几个小时的无效计算。模拟过程中如果再出现NaN，就逐步把u0从0.05往下调，同时观察是哪个位置的密度先变成负值——那个位置通常是固体尖角，对应着数值震荡最强的区域。

6. savevtk导出流场与Darcy渗透率标定

6.1 savevtk输出VTK格式供ParaView可视化

MATLAB自带的绘图函数适合展示二维结果，三维流场还是交给ParaView。压缩包里的savevtk.m和savevtkvector.m负责把LBM结果写成VTK格式。VTK的STRUCTURED_POINTS格式结构很简单，标量场和矢量场各写一次头信息，然后把数据按行主序展开：

function savevtk(filename, field) % 标量场ASCII输出，field为二维或三维数组 fid = fopen(filename, 'w'); fprintf(fid, '# vtk DataFile Version 2.0\n'); fprintf(fid, 'LBM output\n'); fprintf(fid, 'ASCII\n'); fprintf(fid, 'DATASET STRUCTURED_POINTS\n'); fprintf(fid, 'DIMENSIONS %d %d %d\n', size(field,1), size(field,2), max([size(field,3),1])); fprintf(fid, 'ORIGIN 0 0 0\n'); fprintf(fid, 'SPACING 1 1 1\n'); fprintf(fid, 'POINT_DATA %d\n', numel(field)); fprintf(fid, 'SCALARS rho float\n'); fprintf(fid, 'LOOKUP_TABLE default\n'); fprintf(fid, '%g\n', field(:)); fclose(fid); end

保存后在ParaView里打开，先选择“Solid Color”并切换到rho字段，就能看到压力场分布。想看速度场时用savevtkvector.m输出uvw三个分量，然后在ParaView的Filters里加一个Glyph，以速度矢量作为方向源，箭头密度调到合适值即可直接观察涡流和绕流结构。

6.2 用Darcy定律标定渗透率

多孔介质流动模拟完成后，最该做的一件事是用Darcy定律验证结果，这比肉眼看流线图可靠得多。Darcy定律在格子单位下写成：

U = -(K/ν)·(dP/dx)

其中U是截面平均流速，K是渗透率，ν是格子粘度，dP/dx是压力梯度。稳态收敛后，从rho场取进出口两侧的平均压力差：

% 取左右两侧各1/4区域的平均压强，减少局部扰动 nx = size(rho,1); left = mean(rho(1:round(nx/4), :), 'all'); right = mean(rho(nx-round(nx/4):nx, :), 'all'); dP = right - left; % 格子单位压差 Lx = nx; % 流动方向长度 Ubar = mean(sqrt(ux.^2 + uy.^2), 'all'); % 体积平均流速 nu = (tau - 0.5) / 3; K_num = Ubar * nu * Lx / dP;

K_num就是数值渗透率。把它和半经验公式做对比，一种常用的参照是Kozeny-Carman模型：

K_ck = d_p² · ε³ / (180 · (1-ε)²)

其中d_p是颗粒等效直径，ε是孔隙率。随机圆盘模型中d_p取圆盘直径，ε用流体格点数除以总格点数计算。两者数量级一致、偏差在2倍以内，说明模型、边界条件、收敛判据都没有原则性错误；偏差超过一个数量级，优先检查是主通道被堵死还是收敛判据设置太松。渗透率对边界条件极其敏感，反弹格式的具体写法、进出口是否加了缓冲区，都会让K_num产生可观测的差异，这也是LBM结果对比实验时最常被审稿人质疑的地方。用这组脚本把K_num和K_ck同时打出来，标定工作就完成了一半。

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