简介针对多孔介质建模与数值模拟中的几何建模需求一个MATLAB脚本可在平面区域内随机生成圆孔孔隙替代手工构造随机分布的繁琐流程。脚本设计简洁适合材料科学、岩石力学或图像模拟方向的初学者快速上手也可作为子程序嵌入有限元、离散元等更大规模的前处理流程便于在真实算例中复用。压缩包整体仅987B包含1个.m文件核心逻辑集中在单个函数文件中无额外依赖阅读门槛低适合二次修改。已有4112人学习下载在随机结构生成话题中具有不错的参考价值。使用时可自定义圆孔数量、尺寸范围与生成域得到不同分布的孔隙布局为后续网格划分、物质输运或力学分析提供可重复的几何模型基底也有助于对比随机种子、边界条件等因素对生成效果的影响便于学习者理解随机几何生成中的关键控制变量。1. 为什么需要随机孔隙模型先说一个直接的问题为什么要费劲去生成随机孔隙在材料科学、岩石力学、土木工程、过滤分离、增材制造这些领域孔隙结构直接决定了材料的宏观性能。比如混凝土的耐久性和孔隙率强相关岩石的渗透率由孔道连通性主导电池电极和燃料电池的多孔层影响着离子传输效率甚至连食品冻干、制药片剂崩解这种事也跟内部孔隙分布脱不开关系。做这类研究时最理想的办法是拿真实样品去做CT扫描获得三维孔隙结构但扫描成本不低而且样品一多就麻烦。还有一种办法是直接在数值软件里建立理想化模型比如规则排列的圆柱孔、球孔但这类太理想的结构和真实材料差异很大。折中方案就是用MATLAB生成随机孔隙模型——控制孔径分布、孔隙率、空间分布规律得到一批具有统计特征的随机结构再导入COMSOL、Abaqus、Fluent或者自己写有限差分程序去做渗流、传热、应力分析。这篇内容适合谁刚接触多孔介质建模的研究生、需要批量生成结构样本做参数扫描的工程师以及想把孔隙生成这个环节快速集成到现有仿真流程里的朋友。我会从最基础的二维圆孔生成讲起扩展到三维球形孔再讲一种更贴近真实多孔材料的不规则孔道生成方法最后聊透孔隙率控制、常见坑和性能优化。代码都是可以直接复制跑的环境是MATLAB R2024a但向前兼容到R2021a应该没有太大问题。2. 二维随机圆孔生成先把最直观的模型跑通2.1 模型定义与参数约定二维随机圆孔是最容易上手的方案也是很多人做多孔介质模拟的第一步。模型可以这样定义在一个L x L的正方形区域里随机放置若干圆形孔洞每个圆孔的半径在[Rmin, Rmax]间随机取值目标是让孔洞总面积占整个区域面积的比例达到预设值——这个比例就是孔隙率porosity。这里的核心矛盾是如果完全随机放置圆与圆之间会大量重叠实际孔隙率会低于理论值。不处理重叠问题的随机孔放进仿真软件里会得到“两个孔连成一个大孔”的假连通结果渗透率偏大根本不可用。所以生成算法的关键不是随机而是“随机但不重叠”。为了解决这个问题我常用一种思路逐次随机生成圆孔每生成一个就检测跟已有圆孔是否重叠如果重叠就丢弃重来。2.2 边界膨胀与重叠判定有两个细节需要提前说明边界问题和重叠判定方式。对于边界圆孔不能只满足圆心在区域内还必须保证整个圆都在区域内也就是圆心坐标要在[r, L-r]范围内。如果不处理这一点靠边界的孔会被截断孔隙率计算也会失真。有些场景刻意需要边界截断孔那另当别论但作为通用模型我建议先保证圆完整落在区域内。重叠判定的公式很简单两个圆的圆心距离小于两圆半径之和即重叠。放在代码里就是d2 (xi - xj)^2 (yi - yj)^2; if d2 (ri rj)^2 % 重叠 end注意这里用距离平方和半径直接比较避免开平方运算能省一点算力。当圆孔数量达到几千时这种细节会明显影响运行速度。2.3 代码实现下面这个函数是经过简化和参数化的版本可以直接保存为genRandomCircles.m使用function [cx, cy, r] genRandomCircles(L, Rmin, Rmax, targetPorosity, maxAttempts) % 在 LxL 区域内生成不重叠的随机圆孔 % 返回圆心坐标 cx, cy 和半径 r % % L : 区域边长 % Rmin, Rmax : 半径范围 % targetPorosity: 目标孔隙率0~1 % maxAttempts : 每个圆的最大尝试次数防止死循环 % 预估最大圆孔数量按最小半径全部填满时作为上限 nMax ceil(4 * targetPorosity * L^2 / (pi * Rmin^2)); cx zeros(nMax, 1); cy zeros(nMax, 1); r zeros(nMax, 1); n 0; % 预生成半径池并按从大到小排序 % 先放大的圆后续用小的填补空隙能更接近目标孔隙率 rPool Rmin (Rmax - Rmin) * rand(nMax, 1); [rPool, idx] sort(rPool, descend); for k 1:nMax placed false; for attempt 1:maxAttempts % 在扣除边界余量后的区域内随机取圆心 rc rPool(k); xc (L - 2*rc) * rand() rc; yc (L - 2*rc) * rand() rc; % 与已有圆做重叠检测 overlap false; for i 1:n if (cx(i) - xc)^2 (cy(i) - yc)^2 (r(i) rc)^2 overlap true; break; end end if ~overlap n n 1; cx(n) xc; cy(n) yc; r(n) rc; placed true; break; end end if ~placed % 尝试多次仍放不下说明区域接近饱和提前终止 break; end end % 截取有效部分 cx cx(1:n); cy cy(1:n); r r(1:n); % 输出实际的孔隙率 totalArea sum(pi * r.^2); fprintf(目标孔隙率: %.4f, 实际孔隙率: %.4f\n, targetPorosity, totalArea / L^2); end调用方式也非常简单rng(2024); % 固定随机种子保证结果可复现 [cx, cy, r] genRandomCircles(100, 2, 6, 0.3, 200);跑完之后用以下代码快速可视化figure; hold on; rectangle(Position, [0 0 100 100], LineWidth, 2); viscircles([cx cy], r, Color, k, LineWidth, 0.5); axis equal; axis([0 100 0 100]); grid on; title(2D Random Pores);输出效果就是一块100 x 100的区域内散布着大小不一的黑色圆孔互不重叠。我跑了多次当目标孔隙率在30%左右时实际孔隙率一般能控制在28%~32%之间误差主要来自孔间空隙的“浪费”——圆与圆之间始终存在三角形或四边形的小间隙这些间隙无法被圆孔填满这是这种算法天生的限制。2.4 为什么用栅格化计算孔隙率你可能注意到函数里打印的实际孔隙率是用圆面积累加计算的也就是sum(pi * r.^2) / L^2这个值代表“孔洞在几何上的面积占比”的理论上限。但如果把圆孔放到有限元网格里网格单元要么被孔完全覆盖要么完全保留或者部分覆盖这时候的等效孔隙率会略低于几何值。因此在实际工程中我更推荐用栅格化方法来统计孔隙率在高分辨率网格上判断每个网格点是否落在任意圆内然后统计落在圆内的网格点数占总网格数的比例。代码如下gridN 500; [Xg, Yg] meshgrid(linspace(0, L, gridN), linspace(0, L, gridN)); mask false(gridN, gridN); for i 1:length(cx) mask mask | ((Xg - cx(i)).^2 (Yg - cy(i)).^2 r(i)^2); end porosity sum(mask(:)) / gridN^2;这样得到的是“像素级”的孔隙率和有限元网格分辨率相关更接近实际计算中用到的值。有一点可以留个印象网格分辨率越高计算出的孔隙率越接近几何真实值但耗时也越长一般取gridN 500左右足够。3. 三维球形孔隙扩展3.1 从2D到3D的改动二维模型跑通之后扩展到三维球形孔隙几乎是水到渠成的事核心改动只有三处圆心坐标从二维变成三维半径判定从圆变球可视化更复杂一点。其他逻辑完全一致。函数可以写成function [cx, cy, cz, r] genRandomSpheres(L, Rmin, Rmax, targetPorosity, maxAttempts) nMax ceil(6 * targetPorosity * L^3 / (4/3 * pi * Rmin^3)); cx zeros(nMax, 1); cy zeros(nMax, 1); cz zeros(nMax, 1); r zeros(nMax, 1); n 0; rPool Rmin (Rmax - Rmin) * rand(nMax, 1); [rPool, idx] sort(rPool, descend); for k 1:nMax rc rPool(k); placed false; for attempt 1:maxAttempts xc (L - 2*rc) * rand() rc; yc (L - 2*rc) * rand() rc; zc (L - 2*rc) * rand() rc; overlap false; for i 1:n d2 (cx(i) - xc)^2 (cy(i) - yc)^2 (cz(i) - zc)^2; if d2 (r(i) rc)^2 overlap true; break; end end if ~overlap n n 1; cx(n) xc; cy(n) yc; cz(n) zc; r(n) rc; placed true; break; end end if ~placed break; end end cx cx(1:n); cy cy(1:n); cz cz(1:n); r r(1:n); totalVolume sum(4/3 * pi * r.^3); fprintf(目标孔隙率: %.4f, 实际孔隙率: %.4f\n, targetPorosity, totalVolume / L^3); end注意nMax的计算公式变了三维球的体积是4/3 * pi * r^3如果全部用最小半径的球填满目标孔隙率需要的数量上限为targetPorosity * L^3 / (4/3 * pi * Rmin^3)这里保守起见乘了6给后续重叠检测失败留下的余量空间。3.2 三维孔隙率计算三维孔隙率仍用栅格化但这里有一个需要留意的内存问题纯三维网格meshgrid在gridN 300时就会产生300^3 2700万个网格点内存占用轻松超过200MB再大就可能报内存不足。我的经验是这样处理首先把gridN控制在150~200占用内存较小。其次尽量不用显式三维布尔矩阵反复做或运算而是每个球单独生成一个子区域的掩码再合并。下面是一个折中版本gridN 150; step L / (gridN - 1); [xg, yg, zg] ndgrid(linspace(0, L, gridN), linspace(0, L, gridN), linspace(0, L, gridN)); mask3 false(gridN, gridN, gridN); for i 1:length(cx) dist3 (xg - cx(i)).^2 (yg - cy(i)).^2 (zg - cz(i)).^2; mask3 mask3 | (dist3 r(i)^2); end porosity3D sum(mask3(:)) / numel(mask3);球数量在几百的时候逐球加掩码的耗时基本可以接受。如果球数量到了几千甚至上万建议改用分块处理把三维区域切成若干子块逐块统计落在孔内的点数最后汇总。这个思路后面在性能优化部分再说。3.3 可视化与导出三维可视化最简单的方案是画球表面figure; hold on; for i 1:length(cx) [sx, sy, sz] sphere(30); surf(sx * r(i) cx(i), sy * r(i) cy(i), sz * r(i) cz(i), ... FaceColor, interp, EdgeColor, none, FaceAlpha, 0.6); end axis equal; light; material dull;孔径如果只有几十个且球较大这种画法很直观。但球数量一多surf对象会很卡建议每10个球显示一个或者用scatter3只显示孔心位置做位置分布检查。如果后续要导入有限元软件一般不是直接用曲面网格而是把三维孔隙结构体素化输出。可以把mask3保存为.mat文件或者写成二值化图像序列for k 1:gridN imwrite(uint8(mask3(:,:,k)) * 255, sprintf(slice_%03d.tif, k)); end这样每张图像就是一个切片的孔隙分布很多仿真软件可以直接读取图像序列重建几何。4. 不规则连续孔道基于随机场阈值分割4.1 圆孔模型的局限与随机场思路圆孔模型有一个明显的问题真实多孔材料的孔隙往往不是孤立的圆球而是蜿蜒曲折、相互连通的通道。比如砂岩孔隙、海绵状结构、陶瓷膜过滤层里面孔的形状极其不规则孔与孔之间经常有细小的喉道连接。这种连通性对渗透率、扩散系数有决定性影响靠圆孔堆叠模拟不出来。这时候就该引入随机场阈值分割的方法。核心思想是先在每个网格点上生成一个随机数构成一个随机场然后对这个场做平滑处理最后取一个阈值——大于阈值的部分是骨架小于阈值的部分是孔隙。听起来简单但它背后有一个很直观的物理对应孔隙的形成本身就类似一个随机过程通过空间相关长度的控制可以让孔道聚集形成连通的网络而不是均匀撒点。4.2 基于高斯随机场的实现MATLAB里实现随机场阈值分割非常方便。下面是一个完整的例子rng(10); % 生成高斯白噪声场 field randn(200, 200); % 高斯滤波做平滑控制孔道尺度 field imgaussfilt(field, 3); % 归一化到0~1 field (field - min(field(:))) / (max(field(:)) - min(field(:))); % 设定目标孔隙率比如35% targetPorosity 0.35; threshold prctile(field(:), targetPorosity * 100); por field threshold; figure; imshow(por); title(sprintf(Random Field Pores, porosity%.2f, mean(por(:))));这段代码中imgaussfilt(field, 3)中的3是高斯滤波的标准差单位是像素。sigma越大平滑程度越高孔道整体尺度越大、结构越连续sigma越小孔道越碎、越孤立。我拿不同sigma试过sigma 1孔隙呈现很多细小且不连通的小孔看起来像筛子sigma 3~5开始出现比较连续的不规则孔道类似海绵结构sigma 8以上孔隙变成少数几大块连续区域更像裂缝网络。这个参数的物理意义是“材料的特征相关长度”——可以理解为孔隙在空间中相互关联的尺度对应到真实材料里就是孔径和孔间距的量级。没有CT扫描数据做标定时可以先根据想要的孔径大小反推sigma大致等于目标特征孔径的三分之一到一半。4.3 控制连通性与生成三维变体很多时候我们不仅关心孔隙率还关心孔道是否连通。比如做渗流模拟时必须保证孔隙从入口连到出口否则流体根本流不过去。检查连通性用bwconncompCC bwconncomp(por); numComponents CC.NumObjects; fprintf(连通分量数量: %d\n, numComponents);如果连通分量数量很大说明孔隙结构很碎主流方向上的渗流路径很少。想提升连通性可以增大sigma或者对孔隙做形态学膨胀操作把窄喉道连接起来% 对孔隙做一次半径为1像素的膨胀增加连通性 por_dilated imdilate(por, strel(disk, 1));三维情况完全同理把二维矩阵换成三维数组即可field3 randn(150, 150, 150); field3 imgaussfilt3(field3, 3); field3 (field3 - min(field3(:))) / (max(field3(:)) - min(field3(:))); threshold3 prctile(field3(:), targetPorosity * 100); por3 field3 threshold3;这里用到了imgaussfilt3是三维高斯滤波函数直接处理三维数组不需要额外写卷积。三维数组比较大时建议分块生成随机场再拼接但要注意分块边界处滤波会出边界效应一般留一部分重叠区然后裁剪掉。5. 孔隙率精确控制与参数调优5.1 二分法逼近目标孔隙率阈值分割法中我先用的prctile一步到位原理是根据经验分布的分位数直接找到阈值让目标孔隙率精确匹配。但这个方法有个前提随机场经过平滑和归一化之后像素值的分布已知并且单调prctile才可靠。实际中imgausgsfilt处理后的场分布接近正态分布但在边界区域会有些偏差直接用分位数一般没问题分布不对称时会略有偏差。如果遇到比较挑剔的情况比如要求孔隙率精确到千分位可以换成二分法搜阈值lo 0; hi 1; for iter 1:20 t (lo hi) / 2; p mean(field(:) t); if p targetPorosity hi t; else lo t; end end threshold (lo hi) / 2; por field threshold;二分法迭代20次精度已经远超需求而且不依赖场的分布形式对任何单调分布都有效。这个技巧我用到很多其他场景比如调整灰度图像的阈值以得到指定面积的二值区域统一套路。5.2 圆孔生成中孔隙率偏低的原因与对策随机投放圆孔得到的实际孔隙率总低于目标值前面提过核心原因是间隙浪费。圆与圆之间的空隙区域面积总和在孔隙率较高时会变得非常明显。对于这个情况有三个优化方向第一调整半径分布。半径越均匀堆积效率越高孔隙率误差越小半径差异过大时小圆会被大圆“挤”出去导致最终孔数偏少。把Rmin/Rmax的比值控制在0.3以上差别会好很多。第二允许边界截断。如果研究对象本身就是宏观区域内的随机孔洞边界截断是合理假设这样边界处的孔能正常放置实际孔隙率更接近目标。但要清楚导入有限元软件后边界孔会形成开口需要根据工况判断是否保留了这些开口。第三增加预留空间。把目标孔隙率在生成时临时调高几个百分点比如目标是0.3时按0.32去生成生成后再用栅格化重新统计一次不合格就重新生成。这个方法听着不优雅但在工程上很实用尤其是批量生成样本时非常省事。5.3 随机种子、性能与可复现性做研究发论文时评审会要求“数据可重现”。MATLAB的随机数生成器默认每次启动时状态不同所以同一段代码两次运行结果不一样。想复现结果就一条命令rng(2024);这里的2024换成任意整数都可以固定后每次运行生成的孔隙结构完全一致。另外建议在脚本开头固定rng后面所有记录结果的文件名里带上rng的数字比如porosity_2024.mat这样翻旧账的时候能精确对回去。性能方面二维重叠检测在圆孔数量达到几千时双重循环会变慢。两个优化技巧第一在循环内先检查径向距离的粗筛条件比如abs(xi - xj) (ri rj)时直接跳过计算因为两个圆的x方向都分不开y方向根本不用算第二如果孔数量特别大考虑把区域划分成网格索引只检查周围几个格子里的圆这是一个典型的空间索引加速思路MATLAB 中可以用rangesearch函数实现我试过在2000个圆时速度能快10倍以上。6. 常见问题速查表把我在使用过程中遇到的问题整理成了一张表很多都是网上不好搜到的细节问题现象原因解决方案生成结果孔隙率远低于目标值目标0.35实际只有0.22圆孔间隙浪费半径差异太大缩小半径范围或调高生成目标值做预补偿程序运行很长时间不结束圆孔填到一定数量后大量尝试都失败区域已接近饱和maxAttempts设太大设置合理的maxAttempts比如200次失败就break整体用while循环并设最大孔数三维mask计算内存不足报“内存不足”错误gridN太大三维网格爆炸gridN取150以下改用分块遍历用单精度或逻辑数组代替double随机场生成的孔隙太碎、不连通bwconncomp显示上千个连通分量sigma过小相关性不足增大imgaussfilt的sigma对孔隙做形态学膨胀两次运行结果不一样论文图无法复现没固定随机种子开头加rng(固定数字)可视化卡顿绘制大量球或圆时太慢图形对象太多用viscircles时每隔几个画一个或用散射图替代三维用AlphaShape合并后再显示prctile阈值得到的孔隙率和目标有偏差目标0.3实际0.285随机场分布不均匀或边界效应改用二分法阈值搜索确保精确匹配还有一些经验和建议要补充第一先固定参数的物理意义。L是实验样品的特征尺寸Rmin、Rmax来自CT或扫描电镜统计的孔径分布随机场方法里的sigma对应材料的相关长度。不要把参数当成随便调的游戏每个参数最好都跟实际材料对得上这样仿真的结果才有工程参考价值。第二建议把生成和统计功能分开封装。生成算法输出(cx, cy, r)或二值mask统计孔隙率、连通性的代码单独写。这样换模型时不用大量改动也是我做这类项目一贯的模块化习惯。第三大模型生成要善于利用稀疏思想。比如把网格点按概率采样随机抽取一部分点来判断是否在孔内用抽样比例换算孔隙率这样精度虽然略降但速度快很多。对初步筛选结构候选方案来说比精确计算网格快得多。最后分享一点我个人的体会。用MATLAB做随机孔隙建模最舒服的一点是它把随机数生成、图像滤波、形态学操作、可视化这些工具都集成了不用在几个软件之间来回倒腾。我自己实际做项目时通常先用这套MATLAB脚本快速生成几十种孔隙结构的候选样本统计孔隙率和连通性把明显不可行的筛掉再精挑两三组导入到COMSOL里跑详细的多物理场模拟。这个流程已经帮我在两个材料研究项目里省下了大量“生成结构—发现不合适—重新建模”的时间。如果你后续做的东西跟渗透率相关建议把二维圆孔模型升级成三维随机场模型后再导出做模拟因为二维和三维的渗流路径差别很大二维里看起来连通的系统放到三维里可能完全封闭。还有个值得尝试的扩展方向是把随机场生成和分形维数结合用分形布朗运动替代高斯随机场能模拟出更接近天然岩石的粗糙孔隙表面感兴趣的话可以从MATLAB的wfbm函数开始入手。本文还有配套的精品资源点击获取