企业官网建设流程全解析

简介：针对多孔介质建模与数值模拟中的几何建模需求，一个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 控制连通性与生成三维变体

很多时候我们不仅关心孔隙率，还关心孔道是否连通。比如做渗流模拟时，必须保证孔隙从入口连到出口，否则流体根本流不过去。

检查连通性用bwconncomp：

CC = 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. 常见问题速查表

把我在使用过程中遇到的问题整理成了一张表，很多都是网上不好搜到的细节：

还有一些经验和建议要补充：第一，先固定参数的物理意义。L是实验样品的特征尺寸，Rmin、Rmax来自CT或扫描电镜统计的孔径分布，随机场方法里的sigma对应材料的相关长度。不要把参数当成随便调的游戏，每个参数最好都跟实际材料对得上，这样仿真的结果才有工程参考价值。

第二，建议把生成和统计功能分开封装。生成算法输出(cx, cy, r)或二值mask，统计孔隙率、连通性的代码单独写。这样换模型时不用大量改动，也是我做这类项目一贯的模块化习惯。

第三，大模型生成要善于利用稀疏思想。比如把网格点按概率采样，随机抽取一部分点来判断是否在孔内，用抽样比例换算孔隙率，这样精度虽然略降但速度快很多。对初步筛选结构候选方案来说，比精确计算网格快得多。

最后分享一点我个人的体会。用MATLAB做随机孔隙建模，最舒服的一点是它把随机数生成、图像滤波、形态学操作、可视化这些工具都集成了，不用在几个软件之间来回倒腾。我自己实际做项目时，通常先用这套MATLAB脚本快速生成几十种孔隙结构的候选样本，统计孔隙率和连通性把明显不可行的筛掉，再精挑两三组导入到COMSOL里跑详细的多物理场模拟。这个流程已经帮我在两个材料研究项目里省下了大量“生成结构—发现不合适—重新建模”的时间。

如果你后续做的东西跟渗透率相关，建议把二维圆孔模型升级成三维随机场模型后再导出做模拟，因为二维和三维的渗流路径差别很大，二维里看起来连通的系统放到三维里可能完全封闭。还有个值得尝试的扩展方向是把随机场生成和分形维数结合，用分形布朗运动替代高斯随机场，能模拟出更接近天然岩石的粗糙孔隙表面，感兴趣的话可以从MATLAB的wfbm函数开始入手。

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