尧图网站开发 
简介基于D2Q9格子玻尔兹曼方法LBM的多孔介质渗流Matlab仿真源码面向流体力学与渗流物理方向的新手和科研人员用来解决多孔介质中流场模拟的建模与程序实现问题。压缩包内含5个文件包含3个.m源文件、1个txt说明和1个docx算法文档整体仅16KB体量虽小但功能完整目前已有1794人浏览学习属于经过亲测校正、可正常运行的全套项目代码。运行porous.m与SpeedDataOutput等脚本可计算多孔介质中的速度场数据并以Flash格式输出结果帮助直观理解LBM的D2Q9离散模型、孔隙边界处理以及流场演化规律相关应用或科研中可基于这些m文件快速修改参数扩展不同孔隙结构模拟。txt文件简述项目运行流程docx文档补充无约束条件下Prim算法的Matlab实现为类似数值计算提供参考整体实用度高适合边看代码边调试在入门LBM与多孔介质模拟时节省大量编程时间。1. 多孔介质流动模拟为什么值得用 D2Q9 模型做渗流、油藏、岩土或过滤材料仿真的人大多遇到过传统 CFD 在复杂孔隙结构里难收敛的问题。Navier-Stokes 方程在孔隙尺度上需要非常细的网格而且要反复处理速度-压力耦合算一个三维岩心切片往往要跑几小时甚至几天。D2Q9 模型用格子玻尔兹曼方法LBM换了个思路不直接解宏观方程而是让粒子分布函数在离散格子上碰撞和迁移宏观速度与压力从分布函数的矩里恢复出来。多孔介质只需在格子上标出固体节点边界条件天然简单复杂几何不再成为网格生成的噩梦。这套 Matlab 实现把孔隙尺度流动、达西定律验证、速度场输出一条线串通适合两类人一类是想把 LBM 底层编程完全吃透的研究生另一类是工程师想快速评估多孔结构渗透率但不想立刻陷入 OpenFOAM 或 FIUENT 的庞大体系。资源包里porous.m是主程序SpeedDataOutput.m和SpeedDataOutputFlash.m负责把瞬时速度场和稳定后的结果导成矩阵文件。下面会直接对着代码拆解并给出可复现的运行路径和参数标定方法。2. D2Q9 格子玻尔兹曼方法的多孔介质建模原理2.1 九速离散模型与演化方程D2Q9 表示二维九速度方向是 LBM 最经典的速度集合。九个方向权重 w_i 分别是中心方向 4/9轴向方向 1/9对角方向 1/36。每个格子存 9 个分布函数 f_i每一步执行碰撞和迁移f_i(x e_i dt, t dt) f_i(x, t) - (f_i(x, t) - f_i^eq(x, t)) / tau其中 tau 是无量纲松弛时间和流体运动粘度 nu 的关系是nu (tau - 0.5) * cs^2 * dtcs1/sqrt(3) 是格子声速。宏观密度求和各方向分布函数速度则是动量求和除以密度。这套方程通过 Chapman-Enskog 展开可以还原不可压 Navier-Stokes 方程所以 LBM 是微观看是粒子、宏观看是流体的方法不是纯粹的元胞自动机。2.2 多孔介质在格子上的表示方式多孔介质建模的常见做法是用一个布尔矩阵标记固体格点1 表示孔喉流体可以经过0 表示固体骨架。生成方式是随机概率加连通性控制也可以从 CT 图像二值化直接导入。在碰撞-迁移过程中固体格点不参与流体演化遇到固体边界的分布函数执行反弹格式。这个方式的优点是无需在固体表面构造曲线边界精度对于渗透率定性研究完全足够。比随机生成更重要的是保证入口和出口处不出现全封闭的截面积否则流动无法启动。我一般会先生成随机多孔矩阵再用图像形态学检查连通性确保从入口到出口至少存在一条连续流体路径。代码里的porous.m采用随机筛孔方式参数包括多孔介质尺寸、孔隙率 phi、固体障碍块大小以及边界驱动方式。2.3 参数无量纲化与 tau 的选取LBM 里所有物理量都用格子单位你要把物理参数转换到格子空间。关键是无量纲压力梯度或者驱动力。常见做法是在入口施加恒定外力 g_x模拟压力梯度出口为自由流动。松弛时间 tau 必须保持在 0.5 到 1.0 之间接近 0.5 时粘度很小容易数值不稳定超过 1.5 则耗散过大模拟速度慢且精度下降。% 多孔介质流动参数设置示例 clear; clc; nx 200; ny 100; % 计算域尺寸 phi 0.8; % 孔隙率 0.8 tau 0.8; % 松弛时间nu(tau-0.5)/3 gx 1e-5; % 外加驱动力相当于压力梯度 porous rand(ny, nx) phi; % 随机生成流体区域1为流体这里的rand生成均匀分布随机数 phi让每个格点以 phi 概率成为流体。注意这个写法没有排除孤立固体区域但孔隙率较高时通常不影响整体趋势。tau 取 0.8 时格子粘度约 0.1是稳定性和扩散性的折中。驱动力gx如果过大会导致高马赫数误差一般使最大流速不超过 0.1 格子单位。3. Matlab 源码结构与核心实现3.1 项目文件功能拆解包内包含porous.m、SpeedDataOutput.m、SpeedDataOutputFlash.m、说明.txt四个主要文件。porous.m负责初始化、主循环和宏观量统计SpeedDataOutput.m在模拟稳定后输出速度分量 u、v 和压力 p 到 .mat 文件SpeedDataOutputFlash.m则是逐时间步输出瞬时场用于观察流动发展过程。这是 LBM 项目里常见的数据输出分层稳定前用 Flash 模式监控稳定后一次性输出最终场。3.2 分布函数数组与碰撞过程% D2Q9 速度方向定义 cx [0, 1, 0, -1, 0, 1, -1, -1, 1]; cy [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]; % 分布函数初始化赋予平衡态密度1速度0 f zeros(ny, nx, 9); for i 1:9 f(:, :, i) w(i); end这里 cx、cy 是离散速度分量第 1 个方向是静止粒子2~5 是轴向6~9 是对角方向。初始化直接设置成平衡态可以避免早期高频振荡。注意zeros(ny, nx, 9)的顺序Y 在前 X 在后是为了后面索引写成porous(y,x)更自然但在做reshape时要保持维度一致。主循环里的碰撞步将分布函数向平衡态松弛% 碰撞计算宏观量再松弛 rho sum(f, 3); ux (sum(f(:,:,[2 5 6 9]), 3) sum(f(:,:,[4 7 8]), 3)); % 上面这行是示意真正实现要按概率加权求速度实际更严谨的写法是先从分布函数算密度和动量密度再做feq计算。下面这段是porous.m中碰撞核的完整结构for step 1:maxstep % 宏观量 rho sum(f, 3); ux (f(:,:,2) f(:,:,6) f(:,:,9) - f(:,:,4) - f(:,:,7) - f(:,:,8)) ./ rho; uy (f(:,:,3) f(:,:,6) f(:,:,7) - f(:,:,5) - f(:,:,8) - f(:,:,9)) ./ rho; % 平衡态分布 for i 1:9 um cx(i)*ux cy(i)*uy; feq(:, :, i) w(i) .* rho .* (1 3*um 4.5*um.^2 - 1.5*(ux.^2 uy.^2)); end % 碰撞 f f - (f - feq) / tau; % 迁移周期边界用 circshift反弹固体用手动交换 for i 1:9 f(:, :, i) circshift(f(:, :, i), [cx(i), cy(i)]); end endcircshift实现周期性迁移但这对多孔介质并不完全合适左右边界应该是压力驱动而不是周期。原代码的解决方式是入口加恒定密度、出口外推靠外力实现所以你在具体工程里不要照搬这段 circshift要看porous.m里的边界处理。3.3 固体格点的反弹格式多孔介质固体格点在迁移前需要标记否则粒子会穿进骨架。反弹格式最简单的实现迁移前记录那些指向固体格点的分布函数迁移后把这些值沿反方向写回。% 标记固体 solid ~porous; % 迁移后处理以 direction 2向右为例 % 如果右侧邻居是固体则分布函数反弹回来相当于反射 for i 1:9 xshift cx(i); yshift cy(i); fx mod(x xshift - 1, nx) 1; % 周期索引 fy mod(y yshift - 1, ny) 1; % 对固体目标做反弹这里的代码要结合胞元索引操作 end实际项目中为速度考虑反弹通常写成对每个固体格点把从流体格点迁移过来的值收集然后把对应相反方向的分量写回相邻流体格点。也可以用半反弹边界公式f_opposite(x_f, tdt) f_i(x_f, tdt) - 2 * w_i * rho_wall * (e_i dot u_wall)当壁面速度为 0 时简化成直接反转。porous.m为了通用性用的是标准反弹加固液界面修正保证无滑移边界条件。4. 模拟运行、达西定律验证与结果输出4.1 运行流程与收敛判断启动模拟前先运行porous.m它会调用参数初始化、生成多孔介质、执行主循环并在稳定后调用输出函数。判断是否达到稳态有一个实用标准监测出入口质量流量当相邻 2000 步的流量变化小于 1e-8 时视为收敛。下面这段代码可以加在主循环末尾% 记录入口处平均速度 if mod(step, 100) 0 ux_in(step/100) mean(mean(ux(1:ny, 2))); if length(ux_in) 20 diff abs(ux_in(end) - ux_in(end-5)); if diff 1e-8 fprintf(收敛于第 %d 步\n, step); break; end end end这里ux_in数组长度按 100 步为单位递增ux(1:ny, 2)取入口列所有 Y 格点的水平速度均值。注意不能只比相邻两步因为低速流动的随机涨落可能让你误判比较间隔 500 步以上的差值更可靠。4.2 渗透率计算与达西定律流场稳定后达西定律给出表观速度Q -K * A * dP / (mu * L)其中 Q 是体积流量A 是入口截面积L 是长度dP 是进出口压差。在 LBM 格子单位中压差由密度差换算得来dP (rho_in - rho_out) * cs^2。利用模拟得到的流量可反算渗透率 K% 计算渗透率格子单位 u_inlet mean(mean(ux(1:ny, 3))); % 入口区平均速度 Q u_inlet * ny; % 单位厚度流量 dP (rho_in - rho_out) * (1/3); % 压力差 K_LBM Q * nu * nx / dP; % 格子渗透率 K_dimensionless K_LBM / (nx * ny); % 无量纲化便于对比nu就是(tau-0.5)/3。这个参数建议写在输出文件里方便你换算到物理单位。如果你要对比文献中的无量纲渗透率通常把它对孔隙平均直径的平方做归一化。SpeedDataOutput.m就负责把 ux、uy、rho、K 值统一写入.mat文件后续用 MATLAB 的plot或者导入 ParaView 都能可视化。4.3 输出文件格式及后处理技巧SpeedDataOutput.m保存的变量一般包括ux_final、uy_final、rho_final、porous_matrix、parameters结构体。参数结构体中记录了 nx、ny、tau、phi、gx 和渗透率。这种保存方式比把数据写成文本快得多而且便于跨脚本二次分析。如果你需要用 Python 进行机器学习后续处理可以在 MATLAB 里执行save(flow_result.mat, ux_final, uy_final, rho_final, parameters, -v7.3);然后用 Pythonscipy.io.loadmat加载。压降和速度场之间的对应关系通常是达西流验证的核心。下面这张表给出了不同孔隙率下典型参数组合可以用于快速试算孔隙率 phi松弛时间 tau驱动力 gx收敛步数参考稳定性表现0.90.71e-515000稳定0.80.81e-525000稳定0.70.95e-640000较慢但稳定0.61.02e-660000高阻低速0.51.11e-6慢接近渗透极限孔隙率越低有效通道越少要让粒子绕行需要更小驱动力否则流速过快导致数值失稳。tau 大于 1 时耗散大收敛步数明显增加所以不建议孔隙率低于 0.5 时使用单一松弛时间模型而应切换到多松弛模型。5. 从单相流到多孔介质渗透率预测的调试与扩展5.1 调试必看的三个信号遇到模拟发散第一反应不是降低驱动力而是先看三个信号。第一个是最大速度是否超过 0.2LBM 的可压缩误差和速度的三次方相关速度超过 0.3 时压力场会严重畸变。第二个是负密度是否出现rho在碰撞后若出现负值说明 tau 太接近 0.5 或者初始扰动太大。第三个是入口流量是否持续推进多孔介质入口如果被固体堵住ux_in会一直为零此时需要重新生成孔隙矩阵。我调试porous.m时常用的方法是把SpeedDataOutputFlash.m里输出时间间隔调小观察流动前沿的推进形态。如果前沿像锯齿一样杂乱大部分原因是反弹格式索引写错了固体格点的分布函数没有完全反转。更稳妥的办法是先在一个 50x50 的完全自由空间里跑通 LBM再加入多孔介质骨架这样能隔离算法错误和物理参数错误。5.2 扩展非均匀孔隙率与渗透率张量计算原资源是均匀随机多孔介质实际岩心样本往往有层理和裂缝。扩展方法是把porous矩阵从随机布尔矩阵换成按深度渐变概率生成的场% 生成垂直方向孔隙率渐变的例子 for j 1:ny phi_j phi_min (phi_max - phi_min) * j / ny; porous(j, :) rand(1, nx) phi_j; end这种做法能观察层状介质中的优势流通道配合末端出口分段收集流量可以拟合出渗透率的各向异性。计算渗透率张量需要分别施加 X 方向压力梯度和 Y 方向压力梯度得到两个速度场组成 2x2 矩阵。由于原代码只包含单方向驱动力你需要把驱动力同时加到 Y 方向再跑一次分别记录 Qx 和 Qy才能构建完整张量。5.3 与 FIUENT 参数设定的相互校准工程上常用 FIUENT 做多孔介质宏观模拟LBM 结果可以作为微观参数输入。FIUENT 里多孔介质需要设定粘性阻力系数和惯性阻力系数这些系数的来源原本要依靠实验或经验公式。现在你可以从 LBM 模拟得到达西渗透率 K再换算成 FIUENT 的粘性阻力系数 1/α% 在 LBM 结果基础上换算 FIUENT 参数 alpha K_LBM; % 格子渗透率 physical_K alpha * dx_unit^2; % 转化为物理渗透率dx_unit 为格子尺寸 visc_resistance 1 / physical_K; % 粘性阻力系数 1/α其中dx_unit是每个格子对应的物理长度比如每个格子代表 1 微米那么物理渗透率单位就是 μm²。把这个值乘以 1e-12 换成 m² 后再填到 FIUENT 的多孔介质面板中。要注意 LBM 的渗透率基于达西流动只适用于低雷诺数如果 FIUENT 里流速较高还需要补充 Forchheimer 修正项这时 LBM 需要额外做不同压差下的多组模拟来拟合惯性系数。5.4 零速区块的抑制与输出正确性校验porous.m常被忽略的一个坑是随机生成的孔隙率虽然宏观达标但会存在封闭孔隙导致流体绕过这些盲区实际流速比达西定律预测的低。要验证结果是否合理可以用相对渗透率归一法将 LBM 测得的渗透率除以同尺寸完全自由空间的渗透率该理论值约等于 nx²/12比值应小于孔隙率且大于孔隙率的平方量级。如果比值异常小检查连通性。另外输出速度场后用 MATLAB 画流线时要屏蔽固体格点% 绘制流线固体区域标为灰 ux_plot ux_final; ux_plot(solid) NaN; uy_plot uy_final; uy_plot(solid) NaN; figure; contourf(rho_final, 20); hold on; streamline(ux_plot, uy_plot);把固体格点设为NaN后绘图函数会自动跳过非法值流线不会穿过骨架。Flash 输出时记得保存时增加datetime时间戳否则多个瞬时场数据会互相覆盖这在小步长调试时很容易发生。以上就是这套 D2Q9 多孔介质 LBM 模拟从模型到工程交付的完整落地方案。本文还有配套的精品资源点击获取 
简介基于D2Q9格子玻尔兹曼方法LBM的多孔介质渗流Matlab仿真源码面向流体力学与渗流物理方向的新手和科研人员用来解决多孔介质中流场模拟的建模与程序实现问题。压缩包内含5个文件包含3个.m源文件、1个txt说明和1个docx算法文档整体仅16KB体量虽小但功能完整目前已有1794人浏览学习属于经过亲测校正、可正常运行的全套项目代码。运行porous.m与SpeedDataOutput等脚本可计算多孔介质中的速度场数据并以Flash格式输出结果帮助直观理解LBM的D2Q9离散模型、孔隙边界处理以及流场演化规律相关应用或科研中可基于这些m文件快速修改参数扩展不同孔隙结构模拟。txt文件简述项目运行流程docx文档补充无约束条件下Prim算法的Matlab实现为类似数值计算提供参考整体实用度高适合边看代码边调试在入门LBM与多孔介质模拟时节省大量编程时间。1. 多孔介质流动模拟为什么值得用 D2Q9 模型做渗流、油藏、岩土或过滤材料仿真的人大多遇到过传统 CFD 在复杂孔隙结构里难收敛的问题。Navier-Stokes 方程在孔隙尺度上需要非常细的网格而且要反复处理速度-压力耦合算一个三维岩心切片往往要跑几小时甚至几天。D2Q9 模型用格子玻尔兹曼方法LBM换了个思路不直接解宏观方程而是让粒子分布函数在离散格子上碰撞和迁移宏观速度与压力从分布函数的矩里恢复出来。多孔介质只需在格子上标出固体节点边界条件天然简单复杂几何不再成为网格生成的噩梦。这套 Matlab 实现把孔隙尺度流动、达西定律验证、速度场输出一条线串通适合两类人一类是想把 LBM 底层编程完全吃透的研究生另一类是工程师想快速评估多孔结构渗透率但不想立刻陷入 OpenFOAM 或 FIUENT 的庞大体系。资源包里porous.m是主程序SpeedDataOutput.m和SpeedDataOutputFlash.m负责把瞬时速度场和稳定后的结果导成矩阵文件。下面会直接对着代码拆解并给出可复现的运行路径和参数标定方法。2. D2Q9 格子玻尔兹曼方法的多孔介质建模原理2.1 九速离散模型与演化方程D2Q9 表示二维九速度方向是 LBM 最经典的速度集合。九个方向权重 w_i 分别是中心方向 4/9轴向方向 1/9对角方向 1/36。每个格子存 9 个分布函数 f_i每一步执行碰撞和迁移f_i(x e_i dt, t dt) f_i(x, t) - (f_i(x, t) - f_i^eq(x, t)) / tau其中 tau 是无量纲松弛时间和流体运动粘度 nu 的关系是nu (tau - 0.5) * cs^2 * dtcs1/sqrt(3) 是格子声速。宏观密度求和各方向分布函数速度则是动量求和除以密度。这套方程通过 Chapman-Enskog 展开可以还原不可压 Navier-Stokes 方程所以 LBM 是微观看是粒子、宏观看是流体的方法不是纯粹的元胞自动机。2.2 多孔介质在格子上的表示方式多孔介质建模的常见做法是用一个布尔矩阵标记固体格点1 表示孔喉流体可以经过0 表示固体骨架。生成方式是随机概率加连通性控制也可以从 CT 图像二值化直接导入。在碰撞-迁移过程中固体格点不参与流体演化遇到固体边界的分布函数执行反弹格式。这个方式的优点是无需在固体表面构造曲线边界精度对于渗透率定性研究完全足够。比随机生成更重要的是保证入口和出口处不出现全封闭的截面积否则流动无法启动。我一般会先生成随机多孔矩阵再用图像形态学检查连通性确保从入口到出口至少存在一条连续流体路径。代码里的porous.m采用随机筛孔方式参数包括多孔介质尺寸、孔隙率 phi、固体障碍块大小以及边界驱动方式。2.3 参数无量纲化与 tau 的选取LBM 里所有物理量都用格子单位你要把物理参数转换到格子空间。关键是无量纲压力梯度或者驱动力。常见做法是在入口施加恒定外力 g_x模拟压力梯度出口为自由流动。松弛时间 tau 必须保持在 0.5 到 1.0 之间接近 0.5 时粘度很小容易数值不稳定超过 1.5 则耗散过大模拟速度慢且精度下降。% 多孔介质流动参数设置示例 clear; clc; nx 200; ny 100; % 计算域尺寸 phi 0.8; % 孔隙率 0.8 tau 0.8; % 松弛时间nu(tau-0.5)/3 gx 1e-5; % 外加驱动力相当于压力梯度 porous rand(ny, nx) phi; % 随机生成流体区域1为流体这里的rand生成均匀分布随机数 phi让每个格点以 phi 概率成为流体。注意这个写法没有排除孤立固体区域但孔隙率较高时通常不影响整体趋势。tau 取 0.8 时格子粘度约 0.1是稳定性和扩散性的折中。驱动力gx如果过大会导致高马赫数误差一般使最大流速不超过 0.1 格子单位。3. Matlab 源码结构与核心实现3.1 项目文件功能拆解包内包含porous.m、SpeedDataOutput.m、SpeedDataOutputFlash.m、说明.txt四个主要文件。porous.m负责初始化、主循环和宏观量统计SpeedDataOutput.m在模拟稳定后输出速度分量 u、v 和压力 p 到 .mat 文件SpeedDataOutputFlash.m则是逐时间步输出瞬时场用于观察流动发展过程。这是 LBM 项目里常见的数据输出分层稳定前用 Flash 模式监控稳定后一次性输出最终场。3.2 分布函数数组与碰撞过程% D2Q9 速度方向定义 cx [0, 1, 0, -1, 0, 1, -1, -1, 1]; cy [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]; % 分布函数初始化赋予平衡态密度1速度0 f zeros(ny, nx, 9); for i 1:9 f(:, :, i) w(i); end这里 cx、cy 是离散速度分量第 1 个方向是静止粒子2~5 是轴向6~9 是对角方向。初始化直接设置成平衡态可以避免早期高频振荡。注意zeros(ny, nx, 9)的顺序Y 在前 X 在后是为了后面索引写成porous(y,x)更自然但在做reshape时要保持维度一致。主循环里的碰撞步将分布函数向平衡态松弛% 碰撞计算宏观量再松弛 rho sum(f, 3); ux (sum(f(:,:,[2 5 6 9]), 3) sum(f(:,:,[4 7 8]), 3)); % 上面这行是示意真正实现要按概率加权求速度实际更严谨的写法是先从分布函数算密度和动量密度再做feq计算。下面这段是porous.m中碰撞核的完整结构for step 1:maxstep % 宏观量 rho sum(f, 3); ux (f(:,:,2) f(:,:,6) f(:,:,9) - f(:,:,4) - f(:,:,7) - f(:,:,8)) ./ rho; uy (f(:,:,3) f(:,:,6) f(:,:,7) - f(:,:,5) - f(:,:,8) - f(:,:,9)) ./ rho; % 平衡态分布 for i 1:9 um cx(i)*ux cy(i)*uy; feq(:, :, i) w(i) .* rho .* (1 3*um 4.5*um.^2 - 1.5*(ux.^2 uy.^2)); end % 碰撞 f f - (f - feq) / tau; % 迁移周期边界用 circshift反弹固体用手动交换 for i 1:9 f(:, :, i) circshift(f(:, :, i), [cx(i), cy(i)]); end endcircshift实现周期性迁移但这对多孔介质并不完全合适左右边界应该是压力驱动而不是周期。原代码的解决方式是入口加恒定密度、出口外推靠外力实现所以你在具体工程里不要照搬这段 circshift要看porous.m里的边界处理。3.3 固体格点的反弹格式多孔介质固体格点在迁移前需要标记否则粒子会穿进骨架。反弹格式最简单的实现迁移前记录那些指向固体格点的分布函数迁移后把这些值沿反方向写回。% 标记固体 solid ~porous; % 迁移后处理以 direction 2向右为例 % 如果右侧邻居是固体则分布函数反弹回来相当于反射 for i 1:9 xshift cx(i); yshift cy(i); fx mod(x xshift - 1, nx) 1; % 周期索引 fy mod(y yshift - 1, ny) 1; % 对固体目标做反弹这里的代码要结合胞元索引操作 end实际项目中为速度考虑反弹通常写成对每个固体格点把从流体格点迁移过来的值收集然后把对应相反方向的分量写回相邻流体格点。也可以用半反弹边界公式f_opposite(x_f, tdt) f_i(x_f, tdt) - 2 * w_i * rho_wall * (e_i dot u_wall)当壁面速度为 0 时简化成直接反转。porous.m为了通用性用的是标准反弹加固液界面修正保证无滑移边界条件。4. 模拟运行、达西定律验证与结果输出4.1 运行流程与收敛判断启动模拟前先运行porous.m它会调用参数初始化、生成多孔介质、执行主循环并在稳定后调用输出函数。判断是否达到稳态有一个实用标准监测出入口质量流量当相邻 2000 步的流量变化小于 1e-8 时视为收敛。下面这段代码可以加在主循环末尾% 记录入口处平均速度 if mod(step, 100) 0 ux_in(step/100) mean(mean(ux(1:ny, 2))); if length(ux_in) 20 diff abs(ux_in(end) - ux_in(end-5)); if diff 1e-8 fprintf(收敛于第 %d 步\n, step); break; end end end这里ux_in数组长度按 100 步为单位递增ux(1:ny, 2)取入口列所有 Y 格点的水平速度均值。注意不能只比相邻两步因为低速流动的随机涨落可能让你误判比较间隔 500 步以上的差值更可靠。4.2 渗透率计算与达西定律流场稳定后达西定律给出表观速度Q -K * A * dP / (mu * L)其中 Q 是体积流量A 是入口截面积L 是长度dP 是进出口压差。在 LBM 格子单位中压差由密度差换算得来dP (rho_in - rho_out) * cs^2。利用模拟得到的流量可反算渗透率 K% 计算渗透率格子单位 u_inlet mean(mean(ux(1:ny, 3))); % 入口区平均速度 Q u_inlet * ny; % 单位厚度流量 dP (rho_in - rho_out) * (1/3); % 压力差 K_LBM Q * nu * nx / dP; % 格子渗透率 K_dimensionless K_LBM / (nx * ny); % 无量纲化便于对比nu就是(tau-0.5)/3。这个参数建议写在输出文件里方便你换算到物理单位。如果你要对比文献中的无量纲渗透率通常把它对孔隙平均直径的平方做归一化。SpeedDataOutput.m就负责把 ux、uy、rho、K 值统一写入.mat文件后续用 MATLAB 的plot或者导入 ParaView 都能可视化。4.3 输出文件格式及后处理技巧SpeedDataOutput.m保存的变量一般包括ux_final、uy_final、rho_final、porous_matrix、parameters结构体。参数结构体中记录了 nx、ny、tau、phi、gx 和渗透率。这种保存方式比把数据写成文本快得多而且便于跨脚本二次分析。如果你需要用 Python 进行机器学习后续处理可以在 MATLAB 里执行save(flow_result.mat, ux_final, uy_final, rho_final, parameters, -v7.3);然后用 Pythonscipy.io.loadmat加载。压降和速度场之间的对应关系通常是达西流验证的核心。下面这张表给出了不同孔隙率下典型参数组合可以用于快速试算孔隙率 phi松弛时间 tau驱动力 gx收敛步数参考稳定性表现0.90.71e-515000稳定0.80.81e-525000稳定0.70.95e-640000较慢但稳定0.61.02e-660000高阻低速0.51.11e-6慢接近渗透极限孔隙率越低有效通道越少要让粒子绕行需要更小驱动力否则流速过快导致数值失稳。tau 大于 1 时耗散大收敛步数明显增加所以不建议孔隙率低于 0.5 时使用单一松弛时间模型而应切换到多松弛模型。5. 从单相流到多孔介质渗透率预测的调试与扩展5.1 调试必看的三个信号遇到模拟发散第一反应不是降低驱动力而是先看三个信号。第一个是最大速度是否超过 0.2LBM 的可压缩误差和速度的三次方相关速度超过 0.3 时压力场会严重畸变。第二个是负密度是否出现rho在碰撞后若出现负值说明 tau 太接近 0.5 或者初始扰动太大。第三个是入口流量是否持续推进多孔介质入口如果被固体堵住ux_in会一直为零此时需要重新生成孔隙矩阵。我调试porous.m时常用的方法是把SpeedDataOutputFlash.m里输出时间间隔调小观察流动前沿的推进形态。如果前沿像锯齿一样杂乱大部分原因是反弹格式索引写错了固体格点的分布函数没有完全反转。更稳妥的办法是先在一个 50x50 的完全自由空间里跑通 LBM再加入多孔介质骨架这样能隔离算法错误和物理参数错误。5.2 扩展非均匀孔隙率与渗透率张量计算原资源是均匀随机多孔介质实际岩心样本往往有层理和裂缝。扩展方法是把porous矩阵从随机布尔矩阵换成按深度渐变概率生成的场% 生成垂直方向孔隙率渐变的例子 for j 1:ny phi_j phi_min (phi_max - phi_min) * j / ny; porous(j, :) rand(1, nx) phi_j; end这种做法能观察层状介质中的优势流通道配合末端出口分段收集流量可以拟合出渗透率的各向异性。计算渗透率张量需要分别施加 X 方向压力梯度和 Y 方向压力梯度得到两个速度场组成 2x2 矩阵。由于原代码只包含单方向驱动力你需要把驱动力同时加到 Y 方向再跑一次分别记录 Qx 和 Qy才能构建完整张量。5.3 与 FIUENT 参数设定的相互校准工程上常用 FIUENT 做多孔介质宏观模拟LBM 结果可以作为微观参数输入。FIUENT 里多孔介质需要设定粘性阻力系数和惯性阻力系数这些系数的来源原本要依靠实验或经验公式。现在你可以从 LBM 模拟得到达西渗透率 K再换算成 FIUENT 的粘性阻力系数 1/α% 在 LBM 结果基础上换算 FIUENT 参数 alpha K_LBM; % 格子渗透率 physical_K alpha * dx_unit^2; % 转化为物理渗透率dx_unit 为格子尺寸 visc_resistance 1 / physical_K; % 粘性阻力系数 1/α其中dx_unit是每个格子对应的物理长度比如每个格子代表 1 微米那么物理渗透率单位就是 μm²。把这个值乘以 1e-12 换成 m² 后再填到 FIUENT 的多孔介质面板中。要注意 LBM 的渗透率基于达西流动只适用于低雷诺数如果 FIUENT 里流速较高还需要补充 Forchheimer 修正项这时 LBM 需要额外做不同压差下的多组模拟来拟合惯性系数。5.4 零速区块的抑制与输出正确性校验porous.m常被忽略的一个坑是随机生成的孔隙率虽然宏观达标但会存在封闭孔隙导致流体绕过这些盲区实际流速比达西定律预测的低。要验证结果是否合理可以用相对渗透率归一法将 LBM 测得的渗透率除以同尺寸完全自由空间的渗透率该理论值约等于 nx²/12比值应小于孔隙率且大于孔隙率的平方量级。如果比值异常小检查连通性。另外输出速度场后用 MATLAB 画流线时要屏蔽固体格点% 绘制流线固体区域标为灰 ux_plot ux_final; ux_plot(solid) NaN; uy_plot uy_final; uy_plot(solid) NaN; figure; contourf(rho_final, 20); hold on; streamline(ux_plot, uy_plot);把固体格点设为NaN后绘图函数会自动跳过非法值流线不会穿过骨架。Flash 输出时记得保存时增加datetime时间戳否则多个瞬时场数据会互相覆盖这在小步长调试时很容易发生。以上就是这套 D2Q9 多孔介质 LBM 模拟从模型到工程交付的完整落地方案。本文还有配套的精品资源点击获取