简介这套基于Matlab开发的多孔介质孔隙网络建模软件包集成了从图像处理到流体传输模拟的完整工作流面向从事油气储层、地下水渗流、材料科学等研究的科研人员与工程师可用于构建虚拟孔隙网络并预测渗透率、孔隙率等物性参数。资源共398个文件压缩包约10.11MB主要包含Matlab脚本、VTK可视化数据、DAT格式数据文件等其中脚本负责建模与计算VTK/DAT文件提供孔隙网络结构与示例数据支撑。目前已有86人学习使用适合希望借助现成代码快速开展PNM模拟的Matlab用户。包内涵盖图像分割与骨架化、孔隙网络提取、物性计算、多相流模拟及参数化研究等核心模块并附有Berea砂岩等经典样例数据用户可结合脚本与数据直接运行、修改参数或嵌入自己的研究中有效降低多孔介质建模的门槛。1. 从体素到预测值Matlab孔隙网络建模到底能省多少事拿到一块数字岩心或者CT扫描的多孔电极样品很多人的第一反应是拿格子玻尔兹曼直接算渗透率结果一个几百体素边长的小模型就够工作站跑一整夜。基于Matlab的多孔介质孔隙网络建模软件包走的是另一条路先把复杂的孔隙空间抽象成由球孔隙和杆喉道组成的网络再把流动模拟降到几百个节点上求解同一份数据几分钟内就能给出孔隙度、孔径分布、配位数、渗透率和毛管压力曲线。它解决的是“采样容易、计算难”的典型痛点适合岩石物理、电化学储能材料、陶瓷膜等方向的研究生和一线工程师前提是你接受“统计等价”这个精度档位。2. 先把结构看清楚孔隙网络模型的核心假设与软件包模块2.1 为什么要做降维从数字岩心直接模拟说起数字岩心技术普及之后问题从“采不到数据”变成了“算不动模型”。一个1000×1000×1000体素的二值图直接做流体模拟内存和算力都很紧张。我曾在某实验室的工作站上跑过一次格子玻尔兹曼模拟600个体素边长的小模型细分网格跑了两天换成同一份数据抽出的孔隙网络模型两分钟出渗透率误差控制在可以接受的范围内。这个时间差就是孔隙网络模型存在的意义。孔隙网络模型的降维思路非常朴素把连续的孔隙空间压缩成两类几何体——孔隙体相对宽敞的位置抽象成球和喉道连接孔隙的最窄通道抽象成圆柱或锥形杆。流动阻力由喉道半径和长度决定储集能力由孔隙体积决定这样就绕开了在复杂边界上求解纳维-斯托克斯方程的麻烦。软件包内所有模块都围绕这个假设展开所以用之前必须想清楚一件事你的样品是不是“球杆模型”能描述的。随机分布的多孔材料没问题裂缝主导或者强非均质样品误差会明显放大这部分判断在后续参数标定里还要反复提到。2.2 模块划分与一次完整的处理流程解压软件包之后典型的目录结构是下面这样root/ netgen/ % 网络生成骨架化、最大球提取、孔喉匹配 netstats/ % 网络统计孔径、喉径、配位数、孔喉比 simulate/ % 流动模拟单相渗透率、两相入侵渗流 visualize/ % 可视化与导出球棒图、毛管压力曲线 example/ % 示例数据与一键运行脚本 run_pipeline.m % 主入口脚本先解释每个模块的职责netgen负责把体素数据变成网络结构是整个流程中最耗时也最看参数的一部分netstats对生成的网络做统计输出孔径分布、配位数这些中间量用来检查网络是否合理simulate在网络上做流动模拟单相出渗透率两相出毛管压力-饱和度曲线和相对渗透率visualize把网络画成球棒图把曲线导出成CSV。通常我拿到软件包会先跑一遍example目录下的数据。example里自带一个模拟生成的均质样品的二值体素数据体素边长在文档里写的是1微米。直接运行run_pipeline.m会依次调用四个模块并输出三张图网络图、孔径分布直方图、毛管压力曲线。这一步的目的是确认环境没问题再换成自己的数据。值得注意的一点是主脚本里的默认参数针对的都是微米级体素如果你手里是纳米级CT数据单位换算问题在第5章会专门展开。2.3 输入数据约定与最容易被忽略的细节软件包要求输入是一个三维逻辑数组true代表孔隙false代表固体颗粒。这看起来很简单但我在实际使用中见过不少翻车案例有人直接把uint8的灰度体素数据丢进去结果孔隙度算出来百分之九十九因为整个数组都是非零整数。正确的做法是先做阈值二值化再交给软件包。下面是加载数据并做基本检查的示意代码% 加载体素数据假定文件里是一个三维uint8数组 raw load(voxel_data.mat).vol; % 三维数组单位是体素 % 用OTSU阈值做二值化避免手动选阈值 bw raw graythresh(raw); % 转成逻辑类型true代表孔隙 bw logical(bw); % 统计基本参数先确认数据没读反 porosity nnz(bw) / numel(bw); fprintf(体素尺寸: %d x %d x %d\n, size(bw)); fprintf(孔隙度: %.3f\n, porosity); % 软件包要求的输入结构 data.voxel_size 1.0; % 体素边长单位微米必须根据实际CT分辨率修改 data.porosity_map bw;代码里的voxel_size是最容易被忽视的参数。假设你的CT分辨率是每个体素0.5微米这里写成1.0后面算出来的渗透率会直接偏大一个量级。另一个细节是孔隙度的检验先用自己的脚本统计一下孔隙度再和样品的历史数据对比如果对不上说明二值化阈值有问题或者数据里混入了噪声体素这时候不要继续往下跑先解决数据问题。3. 从体素到网络骨架化、最大球提取与参数标定3.1 骨架化保留拓扑、砍掉毛刺网络生成的第一步是提取孔隙空间的骨架。骨架是孔隙空间的中轴线保留了完整的拓扑结构但把体积压缩成了线状结构。Matlab内置的bwskel函数可以直接用但直接作用在大三维数组上会很慢而且会产生大量短枝——这些短枝是噪声体素造成的假分支会影响后续的孔喉匹配。软件包的netgen模块默认的做法是两步走先做距离变换再在距离场上找局部极大值作为候选孔隙中心最后用骨架线把这些候选中心连接起来。距离变换的结果就是每个孔隙体素到最近固体边界的距离这个距离值本质上定义了能从该体素塞进去的最大内切球半径是后续所有算法的地基。下面是一段距离变换与候选中心提取的示意代码% 三维距离变换每个孔隙体素到最近固体边界距离 D bwdist(~bw, euclidean); % 只保留距离大于阈值的体素过滤表面噪声 min_r 2; % 最小半径阈值单位体素过小会保留大量噪声球 D_f D; D_f(D_f min_r) 0; % 找局部极大值作为候选孔隙中心 % imregionalmax 会返回局部极大位置这里把它和D_f叠加 cand imregionalmax(D_f); % 提取候选中心坐标与半径 [idx, idy, idz] ind2sub(size(cand), find(cand)); radii D_f(cand); fprintf(候选孔隙中心数量: %d\n, numel(radii));代码里min_r是关键参数默认值是2个体素。如果样品里有大量微孔这个值可以降到1如果样品颗粒较粗、噪声多建议提到3或4。min_r设得太小的后果是网络里塞满了半径只有一两个体素的假孔隙喉道数量虚高配位数统计失真设得太大又会丢掉真实的小孔。我在处理某碳酸盐岩样品时min_r从2改到4网络孔隙度直接掉了8个百分点孔径分布曲线也变了形所以这个参数不要拍脑袋定要对着孔径分布直方图反复调。3.2 最大球算法把孔隙空间装配成球与杆有了候选孔隙中心下一步是最大球装配。这个算法的思路可以这样理解以每个候选中心为球心向四周膨胀直到碰到固体边界得到一个球然后遍历所有球去掉被更大球完全覆盖的球剩下的就是“局部最大球”代表孔隙体两个相邻最大球的交集区域就是喉道。这套思路在油气行业的数字岩心分析里用了很多年比起直接对骨架做拓扑分析它更贴近“孔隙-喉道”的地质语义。核心装配参数有三个我一般会放到一个结构体里统一管理参数名默认值作用调整建议min_radius2过滤小于该半径的孔隙球单位体素微孔样品下调到1粗颗粒样品上调merge_overlap0.6两个球重叠体积比例超过该值时合并为一个球0.5~0.7之间过小导致孔隙数量虚增boundary_crop3距离边界N层体素内的球直接删除消除边界假孔数据边界效应明显时加大到5min_radius上一节已经说了。merge_overlap这个参数比较微妙重叠比例设成0.6意味着相邻两个大小相近的球只要有60%的体积重合就会合并成一个设成0.9的话原本属于同一个大孔隙的小球会被拆成好几个导致配位数偏高。boundary_crop是用来对付边界效应的——三维体素数据在切割边缘处难免有部分孔隙被截断这些截断处会产生又大又扁的假孔隙球直接裁掉N层是最省事的处理方式。装配完成的网络会输出一个net结构体里面记录了每个孔隙球的中心坐标和半径每条喉道的两端连接哪个球、喉道半径和长度。到了这一步可以先用netstats统计一下配位数和孔喉比正常的多孔材料配位数在3到6之间孔喉比在2到8之间。如果你的数据落在这个区间外面先别急着往下跑模拟回头查参数。% 运行网络生成 net netgen_extract(bw, min_radius, 2, ... merge_overlap, 0.6, ... boundary_crop, 3); % 检查网络统计量是否落在合理区间 stats netstats_summary(net); fprintf(孔隙数量: %d\n, stats.num_pores); fprintf(喉道数量: %d\n, stats.num_throats); fprintf(平均配位数: %.2f\n, stats.mean_coordination); fprintf(孔喉比: %.2f\n, stats.pore_throat_ratio);3.3 参数标定不要让一套参数跑遍所有样品孔隙网络建模里有个说法叫“参数玄学”指的是同样的代码、同样的流程换个样品参数就得重调。我刚用这套软件包的时候拿均质样品标定好的参数直接套在一块裂缝发育的样品上结果渗透率比实验值大了将近20倍原因就是min_radius和merge_overlap这套组合在裂缝区域表现完全不同。我现在的做法是每换一类样品先取一小块代表性区域做参数扫描。具体操作是固定其他参数把min_radius从1扫到5每跑一组就统计网络孔隙度和原始孔隙度的偏差、孔径分布形状选偏差最小且分布形态最接近实验压汞数据的组合。这套做法需要多跑几次但比盲目套参数可靠得多。标定完成后把参数写回脚本文件保存下次同类样品直接复用不要每次都从头调。参数标定的过程是在跑偏与修正之间来回拉锯软件包只是工具真正决定网络质量的是你对样品的理解和参数的取舍。4. 把网络变成预测值单相渗透率、两相毛管压力与可视化4.1 单相流动模拟把网络当成电阻网络网络建好之后单相渗透率的计算可以理解成求解一个大型电阻网络。每个喉道是电阻孔隙节点是电压测量点驱动力是进出口压差。喉道的传导率按圆管层流公式计算这是Hagen-Poiseuille定律的直接应用。如果喉道是圆柱形传导率跟半径的四次方成正比意味着喉道半径的微小误差会被急剧放大这也是第5章要重点排查渗透率偏差的原因。软件包的simulate模块会把网络自动组装成稀疏线性方程组然后用共轭梯度法求解节点压力再反算总流量和等效渗透率。使用者需要提供的参数很少入口压力、出口压力、流体黏度。下面是调用示意% 单相流参数设置 params.p_in 101325; % 入口压力单位Pa params.p_out 0; % 出口压力单位Pa params.viscosity 1e-3; % 流体黏度单位Pa.s模拟水 params.length 50; % 模型沿流动方向长度单位um params.area 2500; % 横截面积单位um^2 % 运行单相模拟 result simulate_single_phase(net, params); % 输出的绝对渗透率单位是um^21um^2 1.01325 D K_md result.perm * 1.01325e12; fprintf(绝对渗透率: %.3f mD\n, K_md);代码里params.length和params.area必须和你的数据区域一致。很多人在这里忘记了软件包内部所有长度单位都是体素net里记录的半径和长度都是“多少个体素”而不是物理长度。如果体素边长是0.5微米喉道半径实际是r_voxel乘以0.5。软件包在simulate_single_phase内部会读取data.voxel_size做换算但你必须在netgen_extract之前就把voxel_size设对否则网络生成阶段就已经错了。4.2 两相入侵渗流从毛管压力曲线到相对渗透率两相流模拟用的是入侵渗流算法思路是非润湿相从入口面逐步侵入网络每次只侵入当前毛管压力下能够进入的最小喉道。喉道能否被侵入由Young-Laplace方程决定毛管压力等于两倍界面张力乘以接触角余弦再除以喉道半径。喉道越细需要的毛管压力越大所以随着压差升高非润湿相逐步占领越来越小的孔隙空间。入侵渗流的输出是毛管压力-饱和度曲线Pc-Sw曲线和相对渗透率曲线。Pc-Sw曲线是评价储层岩石最常用的曲线之一也是和实验压汞数据对照的主要依据。下面是两相模拟的核心调用% 两相参数设置水为润湿相油为非润湿相 params2.gamma 0.03; % 界面张力单位N/m params2.theta 30; % 接触角单位度水湿系统常见取值 params2.wetting water; % 润湿相标识 params2.nonwetting oil; % 运行入侵渗流 pc_sw simulate_two_phase(net, params2); % 导出毛管压力曲线横坐标是润湿相饱和度 csvwrite(pc_sw_curve.csv, [pc_sw.saturation, pc_sw.capillary_pressure]); % 绘制Pc-Sw曲线 semilogy(pc_sw.saturation, pc_sw.capillary_pressure, o-); xlabel(S_w (润湿相饱和度)); ylabel(P_c (Pa)); grid on;接触角theta的选择直接影响曲线位置。水湿系统常见的接触角在0到50度之间油湿系统可能到140度如果接触角给错毛管压力曲线整体平移跟实验数据对不上。gamma界面张力也需要根据实际流体体系设定模拟油水系统用0.03左右没问题模拟气水系统要改成0.07左右。4.3 可视化与导出不只为了好看可视化模块的价值被很多人低估了。网络球棒图能直观告诉你孔隙和喉道在空间上怎么分布哪里有大孔聚集哪里是瓶颈通道。毛管压力曲线导成CSV后可以跟实验室压汞曲线直接叠加对比。下面是一个球棒图绘制的代码示例% 绘制孔隙网络三维视图 figure(Color, w); for i 1:numel(net.throat_radius) % 每条喉道画一条线连接两端孔隙 p1 net.pore_center(net.throat_from(i), :); p2 net.pore_center(net.throat_to(i), :); line([p1(1) p2(1)], [p1(2) p2(2)], [p1(3) p2(3)], ... LineWidth, 1.5, Color, [0.6 0.6 0.6]); end hold on; % 孔隙球用半径映射颜色 scatter3(net.pore_center(:,1), net.pore_center(:,2), net.pore_center(:,3), ... net.pore_radius*30, net.pore_radius, filled); colorbar; xlabel(Voxel); ylabel(Voxel); zlabel(Voxel); title(提取出的孔隙网络);这段代码把喉道画成线条孔隙画成带半径颜色映射的球体。颜色深浅反映孔径大小一眼就能看出网络里有没有异常的大孔或断点。我一般会在跑模拟前先画一张这样的图相当于给网络做一次目检如果看到某个区域全是孤立球、几乎没有连线说明那个区域的喉道在提取阶段被误删了要回头检查参数而不是急着跑模拟。5. 避坑五处容易翻车的细节与排查手段5.1 渗透率比实验值高一个量级喉道半径的“半径错觉”最早用这套软件包跑某砂岩样品时算出来的渗透率比实验室气测值高了将近8倍。排查了很久问题出在喉道半径的取值上。软件包的netgen模块在提取喉道时记录的是两个孔隙球之间的最大内切球半径而不是喉道最窄处的实际收缩半径。真实喉道往往是锥形的最窄处比最大内切球半径小很多。解决方法是给喉道半径乘一个收缩系数常见做法是取0.7到0.85之间具体值要根据孔径分布和实验渗透率联调。从那以后我每跑一个新样品都会先算一次单相渗透率跟实验值对比偏差超过50%就先查收缩系数。5.2 两相模拟出现非润湿相“短路”两相入侵渗流跑出来的Pc-Sw曲线在低饱和度段出现一条水平直线这说明非润湿相在极低压差下直接穿透了整个网络。检查网络结构后发现有个别孔隙球靠一条半径异常大的喉道直接连通了入口和出口形成了“短路通道”。原因是最小球过滤参数min_radius设太小一条由噪声体素构成的假大孔被保留了。解决办法是提高min_radius或者在网络生成后做一次死端剪除——把只连接一个孔隙的喉道全部删掉这些死端对流动没有贡献但会干扰入侵路径判断。软件包的netgen_extract有剪除开关默认开启我是在调试时手动关掉了才出的问题。5.3 模型一大Matlab直接内存耗尽处理500立方体素的样品时Matlab报内存不足崩在距离变换那一步。原因是bwdist在三维数组上会产生一个double类型的距离场单个维度double占8字节500的三次方就是1GB再加上后续的中间变量很容易撑爆内存。解决思路有两个一是把数据切成多个有重叠的块分块处理再合并网络代价是边界效应会多一点二是把数据降采样到250体素左右先跑通流程确认网络质量没问题再上全分辨率。实际项目里我常用后者先用粗网格验证参数趋势再用细网格出最终结果。5.4 网络孔隙度和原始孔隙度对不上网络生成后netstats输出的孔隙度比原始二值图算出的孔隙度低了15%这属于比较严重的偏差。逐个分析影响因素后发现主要原因是boundary_crop3裁掉了边界层的孔隙以及min_radius3把大量半径只有1到2个体素的微孔全部过滤掉了。微孔占比高的样品这两个参数的影响会特别明显。解决办法是把boundary_crop降到1或者0同时把min_radius设成2如果噪声还严重就先用中值滤波对二值图做一次平滑而不是直接提高半径阈值。网络孔隙度和原始孔隙度的偏差控制在3个百分点以内比较合理超过这个数就要怀疑参数设置把真实孔隙过滤掉了。5.5 曲线毛刺和跳变单位换算的锅有一次跑相对渗透率曲线结果曲线在某个饱和度点上突然跳变看起来很不自然。排查到最后发现是输入的体素边长前后不一致生成网络时voxel_size设置的1.0微米而实际数据是0.5微米导致喉道半径的物理值被高估了一倍毛管压力和渗透率全部偏小。毛刺和跳变本质上是网络里某些喉道的半径出现了离群值单位不统一是这类离群值最常见的来源。我每次换数据都会先检查单位换算把voxel_size、长度、面积这些参数写在一个配置结构体里统一管理不散落在各处硬编码。单位问题属于那种“看起来像算法问题、实际是数据问题”的典型排查优先级应该排在最高。6. 进阶技巧批量扫描标定参数并交叉验证网络质量6.1 批量扫描参数选出最优组合单个样品调参效率太低我后来写了一个批量扫描脚本循环改变min_radius和merge_overlap自动统计每组参数下的网络孔隙度偏差和配位数输出一个对比表。下面是扫描脚本的骨架min_r_list 1:5; overlap_list 0.4:0.1:0.8; results []; for mr min_r_list for ov overlap_list net netgen_extract(bw, min_radius, mr, merge_overlap, ov); stats netstats_summary(net); poro_err abs(stats.porosity - porosity_ground_truth) / porosity_ground_truth; results(end1, :) [mr, ov, stats.num_pores, stats.mean_coordination, poro_err]; end end % 按孔隙度误差升序排列打印前10组参数 results sortrows(results, 5); disp(array2table(results(1:10, :), ... VariableNames, {min_radius,overlap,num_pores,coord,poro_err}));6.2 用经验公式交叉验证网络质量批量扫描选出的参数组合还要过一道交叉验证的关。我习惯把多个不同孔隙度的样品放在一起跑把软件包预测的渗透率跟Kozeny-Carman经验公式估算值做对比。如果网络参数合理两组数据在对数坐标下应该呈平行趋势如果偏差过大说明网络在孔隙连接性上系统性失真。这个验证方法成本低也不需要额外实验数据适合在没有实验值时快速判断网络建模质量。从那以后我每换一类样品都会强制走一遍这套流程先画网络球棒图目检再对渗透率验证参数最后批量扫描确认稳定性。希望这个流程也能帮你在自己的数据上少踩几个坑。本文还有配套的精品资源点击获取
