简介基于Matlab开发的多孔介质孔隙网络建模软件包面向石油工程、环境工程、生物工程及材料科学等领域的研究者与工程师用于构建虚拟孔隙网络模型模拟并预测多孔材料内流体的流动、传输及反应特性同时降低对专业编程技能的要求。资源共398个文件压缩包约10.11MB包含54个m功能脚本、196个vtk可视化数据、72个dat数据文件、12张jpg图片以及html文档、MAT数据等覆盖图像处理、孔隙网络生成、物理特性分析与流体传输模拟等完整模块。已有83人学习/浏览。用户可获得完整的Matlab源码、示例数据集与说明文档内含Berea砂岩等多组孔隙网络数据便于直接运行验证支持脚本修改和GUI交互可灵活开展参数化研究在多孔介质结构—性能关系分析与工程预测中具有实用价值。1. 多孔介质建模为什么偏偏是孔隙网络模型多孔介质建模这活儿圈里人应该都不陌生。岩石、混凝土、催化剂载体、燃料电池交换膜、土壤团粒这些东西从微米到百米尺度都充满孔洞和裂缝直接做全尺寸微观模拟算力上根本不现实。搞数值模拟的人都知道可选路径无非三条连续介质模型把孔隙率渗透率当平滑场、直接孔隙尺度模拟扫描图像后做CFD或格子Boltzmann以及孔隙网络模型PNMPore Network Modeling。PNM的思路特别“工程化”把复杂的孔隙空间抽象成一组“节点喉道”的网络节点代表孔隙体存储流体的位置喉道代表连接通道流体经过的瓶颈。这种降维处理丢掉了很多几何细节留下的却是最影响宏观流动的两样东西——连通性和喉道尺寸分布。1956年Fatt第一次用电阻网络类比多孔介质流动之后PNM一路发展到现在已经是从数字岩心预测渗透率、相对渗透率、毛细压力的主力手段之一。我今天要聊的这个Matlab软件包就是把整套PNM流程封装好的研究工具。它解决了什么问题最直接的一点你不用从零开始写网格生成、拉普拉斯矩阵组装、两相流侵位逻辑这些底层代码拿到手改改参数就能跑通一个从“孔隙结构”到“渗透率预测”的完整链路。适合谁用搞岩心分析的研究生、做燃料电池/电解槽膜电极结构优化的工程师、研究土壤水盐运移的水文方向同学以及所有想快速给多孔材料做“虚拟实验”的人。Matlab本身生态成熟矩阵运算、稀疏矩阵求解、可视化一套全通拿来写网络模型属于“杀鸡用牛刀但确实好用”的典型场景。2. 软件包的整体架构与建模流程拆解2.1 模块划分从几何到流动的四层结构这个软件包在结构上很清晰大致可以分成四个模块网络生成模块、几何参数赋值模块、单相流求解模块、多相流求解模块。有些版本还带了电导率、扩散和传热模拟的接口这跟OpenPNM这类Python工具的架构逻辑是同源的。网络生成模块干的是“造骨架”的活。输入可以是规则格网上的随机半径分布也可以是从CT扫描图像提取的真实孔隙-喉道拓扑。如果是后者常见的提取路线有两条最大球法把孔隙空间分解为最大内切球的集合和中轴线法基于距离变换寻找孔隙空间骨架。这个软件包走的是标准路线从二值化图像出发做距离变换用分水岭或最大球聚类识别孔隙体再根据接触关系构建喉道连接。几何赋值模块就是把识别结果翻译成流动计算需要的东西孔隙半径、喉道半径、喉道长度、形状因子shape factor。这部分有个关键点——不可约的几何假设。实际岩石喉道截面绝不是正圆而是各种三角形的蚀刻形态。软件包的默认处理方式是采用三角形/方形截面假设计算形状因子G A / P²并用这个因子修正传导率公式。别小看这一步同样的半径圆形截面和三角形截面的Hagen-Poiseuille修正系数可以差到30%以上直接决定渗透率预测的偏差方向。单相流求解模块解决的是绝对渗透率问题本质上是组装一个压力网络并求解线性系统。每个喉道的流量用Hagen-Poiseuille定律Q g · ΔP / (μ · L)其中g是几何传导率圆形截面情况下g πr⁴/8。对规则圆形喉道半径四次方这个关系意味着喉道半径的微小误差会被放大成渗透率的巨大误差这也是为什么几何提取的精度直接决定预测可信度。多相流求解模块就更复杂了。常用的有两种工作模式侵位渗流invasion percolation和稳态两相流。侵位渗流模拟的是缓慢排水过程湿相如水被非湿相如油或气逐喉道替换替入顺序由入口毛细压力决定Pc 2σcosθ / r这是Yong-Laplace方程的喉道版本。软件包用事件驱动的方式逐步推进每次寻找当前最小突破压力的喉道侵位之后重新更新邻接关系直到非湿相到达出口。稳态两相流则需要在每个网络节点上同时求解两相的质量守恒方程配合相对渗透率滞后模型计算量大得多但能给出饱和度-相对渗透率曲线。2.2 为什么用Matlab写这套东西这个问题我从工程师角度说一下。Python生态里OpenPNM确实很强但Matlab的优势在于矩阵操作跟PNM的数学结构是天然匹配的。PNM的流场求解本质是在稀疏邻接矩阵上做线性代数Matlab的稀疏矩阵、PCG迭代求解器、反斜杠运算符\都是几十年工业验证过的写起来比你手动调scipy.sparse还要顺手。可视化也是加分项scatter3画孔隙球体plot画喉道圆柱一句代码出图做论文插图比自己处理3D渲染省太多时间。另外Matlab的App Designer和脚本一体化很适合教学和快速验证。你可以在一个脚本里改网络尺寸、改孔径分布、改边界条件然后立刻看到渗透率变化趋势。对于做参数敏感性分析这种场景这种交互性远比编译型语言舒服。2.3 网络拓扑的三种获取方式实际操作中网络拓扑怎么来直接决定了后面所有计算的意义。软件包一般支持三种构建方式我分别说下它们的使用场景和坑。第一种是规则格网随机扰动。你在20×20×20网格上每个格点放一个球体孔隙按高斯或对数正态分布随机赋半径相邻节点之间用圆柱喉道连接。这种网络不会对应任何真实样品但适合做方法验证和敏感性分析因为你可以精确控制孔径分布、配位数平均连接数等参数用来检验模拟程序是否正常。第二种是随机球体堆积。用DEM或者随机填充算法生成一堆不重叠的球体球间空隙作为孔隙空间然后提取网络。这是模拟颗粒材料砂土、陶瓷烧结体的标准做法。坑在于提出来以后孔隙体和喉道的定义方式不同结果差异很大你需要统一用的是Delaunay三角化划分Voronoi胞元的方法还是最大球法。第三种是数字岩心图像。这就回到了CT/Micro-CT扫描。二值化、去噪、距离变换、分水岭分割整个链条中有大量参数要调。软件包的优势在这里体现它提供了一整套形态学后处理工具可以把图像直接转成网络数据存成结构体数组struct array然后供后续计算调用。缺点是算得慢一张1000³体素的图像提取一次可能要跑几个小时而且内存吃紧。3. 核心求解细节与参数选择3.1 渗透率计算到底怎么算的这部分是软件包的核心也是大家最容易误解的地方。绝对渗透率的计算流程说起来很简单在x方向施加一个压力差ΔP求解整个网络的压力分布统计通过截面的总流量Q然后反推Darcy定律K Q · μ · L / (A · ΔP)这里A是整个网络的横截面积L是网络长度μ是流体粘度。但真正落地的时候有几个细节要特别注意。压力分布的求解需要组装一个节点压力线性方程组。对每个孔隙体i所有连接的喉道j流入的流量总和等于零不可压缩稳态条件Σ_j g_{ij}(P_i - P_j) 0对所有内部节点成立。入口和出口节点设为固定压力Dirichlet边界。剩下的就是组装一个N×N的稀疏矩阵然后用A\b或者PCG求解。这里我强烈建议用Matlab的sparse函数组装矩阵别用全矩阵——20万节点的网络全矩阵直接内存爆炸稀疏矩阵只要几十MB。求完压力后出口边界上所有喉道流量之和就是总流量Q。注意这时候要用穿过出口面的喉道流量求和不是所有连接到出口节点的喉道都算很容易漏算或重复计算。配位数平均每个节点连接的喉道数对渗透率影响极大。规则六方网格配位数是6钻石网格是4随机网络通常落在3到5之间。配位数越高连通性越好渗透率呈非线性上升。所以如果你发现模拟渗透率比实验低很多先检查配位数是不是偏低了。3.2 传导率模型的选择传导率是喉道几何和流体性质的桥梁也是各种PNM实现里差别最大的地方。这个软件包提供了至少三种传导率选项。第一种是圆管公式g πr⁴ / (8μ)最简单但只适用于圆形截面且长度等于喉道长度的理想情况。第二种是方形/三角形截面修正公式。实际喉道形状不是圆流动阻力更大。对三角形截面引入形状因子G后传导率用修正的Poiseuille公式计算常数系数不同。这个版本更贴近真实岩心代价是需要额外知道喉道形状因子。第三种是考虑孔隙体阻力的等效传导率。真实喉道的两端有孔隙球体流体在孔隙和喉道之间发生额外的收缩/扩张流动损失。软件包一般会算一个串联等效传导率1/g_eff 1/g_throat 1/(2·g_pore)其中g_pore是孔隙体的贡献。这一步对高孔隙比孔隙半径远大于喉道半径的网络特别重要能避免系统性地高估渗透率。我实际测下来对基于砂岩数字岩心提取的网络用圆管公式算出的渗透率通常偏高50%到一倍用修正模型后和实验值误差能控制到20%以内。所以做真实样品研究时偷懒选圆管公式是要付出代价的。3.3 两相流模拟参数谁是湿相、接触角怎么定两相流模拟的参数设置比单相流容易翻车。最关键的是三件事湿相定义、接触角和入口毛细压力范围。湿相的定义取决于岩石的表面润湿性。砂岩天然亲水水是湿相但含油岩心经过老化处理后可能转为混合润湿这时候湿相的定义要跟着变。软件包允许用户指定每个喉道的接触角可以做到空间非均匀润湿性。这个参数对相对渗透率曲线的影响非常大——同样的孔隙结构亲水系统和亲油系统的相对渗透率交叉点可能差出10%的饱和度。入口毛细压力范围的设置决定了模拟过程的阶段划分。如果你用侵位渗流模拟排水初始压力要设在刚好超过最大喉道对应的入口压力因为最大喉道最先被侵入P_entry 2σcosθ / r_max别小看这一步入口压力设太高了会直接跳过很多中间体状态导致饱和度跳跃式变化曲线不平滑。我的经验是初始压力设在最大喉道入口压力的0.8倍然后按指数或线性步长逐步提高。每个压力步下需要让弛豫过程完全停止即没有新的喉道被侵入再记录状态。还有一个坑毛细压力滞后。排驱后如果做吸入imbibition由于陷捕trapping效应曲线不回原路。软件包默认支持这个滞后模拟但需要确认你设置的是排水还是吸入路径。很多新手在这里困惑明明同一个网络两个方向测出的曲线完全不同这是正常物理现象不是程序bug。4. 实操从零跑通一个渗透率预测算例4.1 数据准备与软件包调用这一步我以最常见的规则网络为例完整走一遍流程。% 定义网络尺寸与孔径分布参数 nx 15; ny 15; nz 15; % 网络节点数 mean_r 15e-6; % 平均孔隙半径 15微米 std_r 4e-6; % 标准差 mean_rt 6e-6; % 平均喉道半径 std_rt 1.5e-6; porosity_target 0.25; % 目标孔隙度 % 调用网络生成函数软件包自带 network generate_network(grid, [nx ny nz], ... pore_radius, [mean_r std_r], throat_radius, [mean_rt std_rt], ... shape, circle, target_porosity, porosity_target);生成函数返回的network结构体里包含pore_coords节点坐标、pore_radius、throat_conns喉道两端节点索引、throat_radius和throat_length这些字段。生成完毕后先画出来看一眼确认没有孤立节点和异常大喉道再往下走。4.2 单相流求解与渗透率提取接下来调用求解模块% 设置流体参数 mu 1e-3; % 水粘度 1 mPa·s dp 100; % 施加压差 100 Pa direction 1; % 沿x方向流动 % 计算绝对渗透率 [K, P, Q] compute_absolute_permeability(network, mu, dp, direction); fprintf(渗透率 %.4f D\n, K*1.01325e12); % 单位转换为达西内部逻辑我之前说过了组装稀疏矩阵、求解压力场、统计流量、反推Darcy渗透率。这里提醒一句单位制务必统一最好全程用SI单位米、帕、秒。换算成达西只用最后一步。我第一次跑的时候就是用微米、毫帕混着来结果渗透率出来差了10的6次方量级排查了半天才发现是单位问题。求解完之后可以用plot_pressure_field(network, P)看压力分布压力应该沿x方向平滑递减。为了确认结果是可信的应该做两个验证一是把网格密度翻倍看渗透率是否基本不变判断是否达到REV——代表性单元体积二是取不同随机种子跑多个实现看渗透率的统计方差。如果网格从10³到20³渗透率还在剧烈变化说明你的网络尺寸太小没有代表性要继续加大。4.3 两相流模拟与相对渗透率曲线两相流模拟调用方式类似但要注意设置更多的参数% 设置两相流参数 sigma 0.03; % 界面张力 30 mN/m theta 30 * pi/180; % 接触角 30度水湿 pc_max 5000; % 最大毛细压力 Pa n_steps 30; % 压力步数 % 侵位渗流排水模拟 [Sw, Krw, Kro, pc] invasion_percolation(network, ... sigma, sigma, theta, theta, pc_max, pc_max, n_steps, n_steps); % 绘制相对渗透率曲线 figure; plot(Sw, Krw, b-o, Sw, Kro, r-s); xlabel(S_w); ylabel(Kr); legend(Krw (water), Kro (oil)); grid on;排水过程从完全饱和水Sw1开始非湿相油逐步侵位。你得到的相对渗透率曲线应该呈现典型特征水的相对渗透率随Sw下降而下降油的相对渗透率随Sw下降而上升两者交叉点通常在Sw0.5到0.6之间。这里有个教训侵位渗流默认两端连通条件但如果网络在某个方向不连通存在死端孔隙你会看到Sw降到某个值后死活不再下降那是因为剩余的水被完全困住了。这其实是真实物理但如果你需要对比实验数据务必确认你的网络连通性跟实际样品一致。死端孔隙比例太高会导致模拟的残余油饱和度异常偏高。5. 常见问题与排查技巧实录5.1 矩阵奇异或不收敛这是最频繁踩的坑。症状是求解压场时Matlab报警“Matrix is singular”或者PCG迭代不收敛。绝大多数情况下原因是网络中存在孤立节点也就是某个孔隙体没有任何喉道连接到主网络上。孤立节点的压力方程没法确定矩阵就奇异了。排查方法很简单用conncomp求连通分量看是否有大于1的连通块% 构建邻接矩阵 A sparse([throat_conns(:,1); throat_conns(:,2)], ... [throat_conns(:,2); throat_conns(:,1)], 1, n_nodes, n_nodes); % 检查连通分量 bins conncomp(graph(A)); fprintf(连通分量数量: %d\n, max(bins));如果发现不连通要么重新生成网络时设置最小配位数约束要么在提取图像时调整孔隙合并阈值保证网络中不存在完全孤立的孔隙体。对规则网络来说生成时加一个“至少连接一个邻居”的约束就够了。5.2 渗透率对网格尺寸过度敏感你可能会遇到这个情况网格从10³加密到20³渗透率涨了30%以上。这在规则网格上往往说明边界效应太大。网络尺寸太小边界上的节点贡献了过多比例而边界节点的连通性通常比内部节点差导致渗透率被低估。解法是增大网络尺寸并同时看多个方向的渗透率来稳定评估。有时候你会发现x和y方向渗透率明明是各向同性的但算出来偏差10%以上这就是尺寸不够的表现。5.3 内存不足与求解速度慢20万节点的网络稀疏矩阵求解本身只要几秒但如果网络生成阶段用全矩阵去存邻接关系可能早就崩了。我的建议是全流程优先用稀疏存储。Matlab的sparse非常成熟神经网络、流体网络、图论问题这些场景都是它的主场。如果网络特别大超过50万节点反斜杠\可能不够快换成pcg带不完全Cholesky预条件子L ichol(A); [x, flag] pcg(A, b, 1e-8, 500, L, L);实测下来50万节点问题用PCG大概能比直接求解快三到五倍内存占用也小一圈。5.4 接触角和界面张力对结果影响巨大最后提醒一个容易被忽略的点接触角和界面张力不是随便给的必须查参考文献或实验数据。我之前做过一个燃料电池气体扩散层的模拟只是把接触角从120度改成110度相对渗透率曲线整体移动了5个饱和度点。你以为你在做孔隙结构优化实际上你只是在调润湿性参数。这提醒我们任何参数更改前都要做好敏感性分析确定哪些参数是决定性的。6. 一些个人体会和后续扩展思路用这套软件包跑了大概小半年我把自己的体会沉淀一下。第一件是PNM的预测精度天花板不在于求解器而在于几何提取质量。数字岩心分割参数一变渗透率差出一倍很正常。所以真正花时间的地方不该是调求解器选项而是把图像分割、孔隙识别做到位。第二件是Matlab写PNM的优点是调试方便脚本一跑中间变量全在workspace里随便点开看结构体字段比在Python里反复print要舒服很多。后续想扩展的方向我认为至少有三个值得动手。一是把传热模块串进去做热流耦合的孔隙尺度模拟这在燃料电池和地热领域都有需求。二是加入反应流模块模拟溶解-沉淀过程导致孔隙结构演化对酸化压裂和CO₂矿化封存场景很有价值。三是把提取和模拟封装成App或函数库方便实验室其他人用而不需要每个人去啃脚本。有一点想特别叮嘱这类工具代码哪怕软件包再完善也一定要自己通读核心求解函数的每一行理解矩阵是怎么组装的、边界条件是怎么施加的。因为孔隙网络模型的灵活性太高任何标准代码都无法覆盖所有场景的边界情况最终改代码的还是你自己。最后分享一个干活时的小技巧每次跑模拟前先构造一个只有几个节点的微型网络比如3×3×3手动验证一遍结果能算出来再上大网络。这个习惯帮我避免了很多低级的组装修bug省下的调试时间远超验证那几分钟的成本。本文还有配套的精品资源点击获取