简介：这份MATLAB脚本资源面向需要构建多孔介质几何模型的研究生、科研人员及工程仿真从业者，用于在二维平面内随机生成圆形孔隙结构，可作为渗流分析、材料孔隙率建模或图像处理实验的前处理工具。压缩包内仅含1个m文件，体积约987B，属于轻量级源码，直接运行即可得到随机分布的圆孔图形，便于按需调整孔径、数量与分布范围。已有4152人学习下载，说明该思路在相关课程与课题中具有较高的参考价值。读者可获得一套可复用的随机圆孔生成逻辑，理解如何借助MATLAB的随机数函数与绘图命令控制孔隙位置和尺寸，并在此基础上扩展为多孔介质网格、孔隙率统计或有限元前处理流程，适合作为入门级脚本范例快速嵌入自己的项目。

1. 随机孔隙生成：从一张SEM图到可计算的多孔介质模型

做材料、岩土、电池极片或者过滤膜方向的开发者，大概率都遇到过同一个卡点：手里有扫描电镜图，能看出材料内部是典型的多孔结构，但要把这种“乱中有序”的孔隙形态搬进仿真软件里做渗流、扩散或力学计算，却不知道从哪下手。手工画几个圆孔太假，直接拿真实图像二值化又常常连不成可计算的网格。这个用 MATLAB 做随机孔隙生成的项目，解决的正是这一步——它不依赖任何商业建模软件，纯代码生成可控孔隙率、可控孔径分布、可导出为通用网格格式的三维多孔介质模型。适合正在做多孔介质仿真、需要批量构造随机样本、或者想把孔隙结构参数化研究的研究生和工程开发。下面我按“它到底生成什么 → 怎么跑起来 → 怎么改参数 → 怎么避坑”的顺序，把这份源码拆开讲。

2. 项目结构与技术选型：为什么是 MATLAB 而不是 Python

2.1 目录组织与文件职责

拿到源码包后，先别急着运行，花两分钟把目录结构看清楚，后面改代码会顺很多。这个项目的组织方式很典型，属于“一个主入口 + 若干功能函数 + 输出目录”的结构，没有复杂的包管理，也没有面向对象的重型封装，对新手比较友好。

random_porous_media/ ├── main.m % 主入口，所有参数在这里集中配置 ├── generate_pores.m % 核心：随机孔隙生成算法 ├── pore_size_dist.m % 孔径分布控制（正态/均匀/对数正态） ├── build_voxel.m % 体素化：把几何孔隙转成三维矩阵 ├── export_mesh.m % 导出为 STL / VTK / 文本网格 ├── utils/ │ ├── overlap_check.m % 孔隙重叠检测 │ └── boundary_clip.m % 边界裁剪，防止孔隙超出立方体 ├── data/ │ └── sample_sem.png % 可选的参考图，用于对照孔径范围 └── output/ ├── model.vtk % 默认导出的体素/网格文件 └── pores_info.mat % 每个孔隙的球心、半径、体积记录

这个结构里，main.m是唯一需要你手动改参数的地方，其余函数尽量保持“纯函数”风格——给定输入就返回确定输出，不依赖全局变量。这一点对二次开发很关键：你想换一种孔径分布，只需要改pore_size_dist.m，不用动主流程。

2.2 为什么选 MATLAB 做这件事

很多人第一反应是“Python 有 trimesh、有 scipy，为什么不用 Python”。在这个具体场景下，MATLAB 有三个实际优势。第一，矩阵操作是它的原生能力，而多孔介质体素化本质上就是三维矩阵的布尔运算，A(idx) = 1这种写法比 NumPy 更直观，调试时可以直接在变量编辑器里看三维切片。第二，isosurface、patch、stlwrite这些函数开箱即用，不需要额外装一堆几何库。第三，很多材料实验室的服务器上 MATLAB 是标配，而 Python 环境反而要自己配。当然，如果你团队全是 Python 栈，核心算法逻辑完全可以平移，只是导出部分要换成trimesh或pyvista。

提示：源码里用到的stlwrite在旧版本 MATLAB 中可能不存在，R2018b 之后才内置。如果你的版本较早，需要单独下载同名函数放到 utils 目录。

2.3 核心数据流：从参数到网格

整个项目的数据流可以概括为四步：参数配置 → 随机孔隙生成 → 体素化 → 网格导出。每一步的输出都是下一步的输入，中间结果会保存到output目录，方便你中断后复查。理解这条链路，后面看代码就不会迷路。

% main.m 中的核心调用链（简化版） params = struct(); params.L = 100; % 立方体边长，单位微米 params.target_porosity = 0.3;% 目标孔隙率 params.r_mean = 8; % 平均孔径半径 params.r_std = 2; % 孔径标准差 params.dist_type = 'normal'; % 分布类型 [pores, porosity] = generate_pores(params); % 生成孔隙列表 voxel = build_voxel(pores, params.L); % 体素化 export_mesh(voxel, 'output/model.vtk'); % 导出

这段代码里，params把所有可调参数集中在一起，这是我最推荐的做法——不要把这些数字散落在各个函数里。generate_pores返回的pores是一个结构体数组，每个元素包含center、radius、volume三个字段。build_voxel把它转成L×L×L的逻辑矩阵，1表示孔隙，0表示实体。最后export_mesh负责格式转换。参数说明：L决定模型物理尺寸，target_porosity是你要逼近的孔隙率，r_mean和r_std控制孔径分布形态。

3. 核心算法拆解：随机孔隙是怎么“长”出来的

3.1 随机球体生成与重叠检测

随机孔隙最朴素的实现就是“往立方体里随机扔球”，但直接扔会出问题：球和球之间可能重叠，重叠后实际孔隙率会低于目标值，而且孔径分布也会失真。所以核心逻辑是“生成一个、检测一次、通过才保留”。

function [pores, porosity] = generate_pores(params) L = params.L; target = params.target_porosity; pores = struct('center', {}, 'radius', {}, 'volume', {}); total_vol = 0; max_iter = 10000; % 防止死循环 iter = 0; while total_vol / L^3 < target && iter < max_iter iter = iter + 1; % 1. 按指定分布生成半径 r = pore_size_dist(params); % 2. 随机生成球心，留出边界余量 c = rand(1,3) * (L - 2*r) + r; % 3. 重叠检测 if overlap_check(c, r, pores) continue; % 重叠则丢弃，重新生成 end % 4. 接受该孔隙 new_pore = struct('center', c, 'radius', r, ... 'volume', 4/3*pi*r^3); pores(end+1) = new_pore; total_vol = total_vol + new_pore.volume; end porosity = total_vol / L^3; end

逐段说明：max_iter是必须的，否则当目标孔隙率设得过高（比如 0.6 以上）时，球体几乎无法再塞进去，程序会卡死。rand(1,3) * (L - 2*r) + r这一行保证球心到六个面的距离都大于半径，避免孔隙被边界切掉。overlap_check是性能瓶颈，它要遍历已有孔隙列表，判断球心距离是否小于两半径之和。如果孔隙数量上千，这里会明显变慢，常见优化是用空间网格分桶，只检查邻近格子里的孔隙。

注意：overlap_check默认是“完全不允许重叠”。如果你做的是连通孔隙研究，可以改成允许部分重叠，但要在参数里加一个overlap_tol控制重叠比例，否则孔隙率会失控。

3.2 孔径分布控制：正态、均匀还是对数正态

真实多孔材料的孔径往往不是单一尺寸，而是有一个分布。项目里pore_size_dist.m提供了三种模式，切换只需要改params.dist_type。

function r = pore_size_dist(params) switch params.dist_type case 'normal' r = params.r_mean + params.r_std * randn(); r = max(r, 0.1); % 防止出现负半径 case 'uniform' r = params.r_min + (params.r_max - params.r_min) * rand(); case 'lognormal' mu = log(params.r_mean^2 / sqrt(params.r_std^2 + params.r_mean^2)); sigma = sqrt(log(1 + params.r_std^2 / params.r_mean^2)); r = lognrnd(mu, sigma); otherwise error('未知的分布类型: %s', params.dist_type); end end

正态分布适合孔径集中在某个均值附近的情况，比如某些烧结金属。均匀分布适合孔径范围明确但无集中趋势的场景。对数正态分布最贴近真实岩土和生物材料，因为小孔数量多、大孔数量少。参数说明：r_mean是均值，r_std是标准差，r_min和r_max只在均匀分布下生效。max(r, 0.1)这一行是血泪经验——正态分布理论上可能生成负值，不加保护会导致后续体积计算出错。

3.3 体素化：把几何球体转成三维矩阵

有了孔隙列表，下一步是把它变成仿真软件能读的体素矩阵。体素化的本质是遍历每个体素中心，判断它是否落在任意一个孔隙球内。

function voxel = build_voxel(pores, L) resolution = 1; % 每个体素的边长，单位微米 n = round(L / resolution); voxel = false(n, n, n); [X, Y, Z] = ndgrid(1:n, 1:n, 1:n); for k = 1:length(pores) c = pores(k).center / resolution; r = pores(k).radius / resolution; % 利用球体方程快速筛选 dist2 = (X - c(1)).^2 + (Y - c(2)).^2 + (Z - c(3)).^2; voxel(dist2 <= r^2) = true; end end

这段代码里，resolution决定网格精度，设得越小模型越精细，但内存和计算量按三次方增长。ndgrid生成三个三维坐标矩阵，然后用球体方程dist2 <= r^2做布尔判断。这里有个性能技巧：不要对每个体素循环判断，而是对每个孔隙一次性计算整个网格的距离场，利用 MATLAB 的向量化操作。参数说明：resolution = 1表示一个体素对应一微米，如果你的模型尺度是毫米级，建议改成 0.5 或更小，否则孔隙细节会丢失。

4. 运行与参数调优：从默认配置到自己的模型

4.1 第一次运行：检查环境与路径

拿到源码后，第一步不是直接跑main.m，而是确认三件事：MATLAB 版本、当前工作目录、输出目录是否存在。很多“跑不起来”的问题其实都是路径问题。

% 在命令行执行，检查环境 ver('matlab') % 查看版本 pwd % 确认当前目录 addpath(genpath('utils')); % 把 utils 加入搜索路径 if ~exist('output', 'dir') mkdir('output'); % 自动创建输出目录 end

addpath(genpath('utils'))这一行很关键，因为overlap_check和boundary_clip在子目录里，不加路径会报“未定义函数”。mkdir的判断也是常见做法，避免因为目录不存在导致导出失败。做完这三步，再运行main.m，正常情况下命令行会输出实际孔隙率和生成耗时。

4.2 参数怎么改：孔隙率、孔径、分辨率

参数调优是这个项目最常用的操作，下面这张表把关键参数、作用和推荐范围列清楚，改的时候直接对照。

改参数时有个经验：先把target_porosity设低一点（比如 0.2），确认流程跑通，再逐步提高。因为孔隙率越高，重叠检测的失败率越高，生成时间会非线性增长。如果发现程序跑了几分钟还没结束，先看命令行输出的迭代次数，如果iter接近max_iter，说明目标孔隙率设得太高，或者孔径太大导致空间不够。

4.3 导出格式与下游对接

生成完体素矩阵后，导出格式决定了你能对接哪些仿真软件。项目默认导出 VTK，因为它是体素数据最通用的格式，ParaView、Avizo、COMSOL 都能读。

function export_mesh(voxel, filename) [n1, n2, n3] = size(voxel); [X, Y, Z] = meshgrid(1:n2, 1:n1, 1:n3); % 提取等值面，生成三角网格 fv = isosurface(X, Y, Z, double(voxel), 0.5); % 导出为 STL stlwrite(filename, fv); % 同时保存体素矩阵，方便后续处理 save('output/voxel.mat', 'voxel'); end

isosurface把体素矩阵转成三角面片，阈值 0.5 表示“孔隙和实体的分界”。stlwrite导出 STL，适合做 3D 打印或导入有限元软件。同时保存voxel.mat是个好习惯，因为体素数据可以重新导出成其他格式，不用重新生成。参数说明：filename建议带路径，比如'output/model.vtk'，否则会写到当前目录。如果你的仿真软件需要体素格式而不是表面网格，直接用voxel.mat即可。

5. 避坑与排查：那些让我重跑过好几次的问题

5.1 孔隙率对不上：目标 0.3 实际只有 0.22

现象：设置target_porosity = 0.3，但程序输出的porosity只有 0.22 左右。原因通常是重叠检测太严格，很多球被丢弃后没有补足。解决方法是把overlap_check里的判断条件放宽，允许少量重叠，或者在循环结束后检查实际孔隙率，如果低于目标值，自动降低半径重新生成一批小孔填充。我一般会在generate_pores末尾加一个补偿循环，专门用小半径孔隙填补剩余空间。

5.2 程序卡死：迭代次数跑满也不结束

现象：命令行一直不返回，iter涨到max_iter才停。原因一般是目标孔隙率设得过高，或者r_mean太大导致空间很快被占满。解决方法是先降低target_porosity到 0.2 测试，如果正常，再逐步提高。另一个办法是动态调整半径——当剩余空间变小时，自动减小新生成孔隙的半径，这样能逼近更高的孔隙率。

5.3 导出文件打不开：STL 面片法向错误

现象：STL 导入 ParaView 或 3D 打印切片软件后显示异常，面片缺失或法向朝内。原因是isosurface生成的面片法向取决于体素矩阵的数值方向，如果voxel是逻辑矩阵，double(voxel)后孔隙为 1、实体为 0，等值面 0.5 提取的是孔隙边界，法向可能朝内。解决方法是在导出前用flip翻转体素矩阵，或者用patch的FaceNormals手动校正。我一般会在export_mesh里加一个voxel = ~voxel的选项，根据下游软件的要求切换。

5.4 内存爆掉：分辨率设太小导致矩阵过大

现象：resolution = 0.2时，L = 100的模型会产生500×500×500的逻辑矩阵，虽然逻辑矩阵只占 1 字节每元素，但ndgrid生成的三个双精度坐标矩阵会瞬间吃掉几十 GB 内存。解决方法是不要用ndgrid生成全局坐标矩阵，而是对每个孔隙用局部窗口计算。具体做法是只计算孔隙包围盒范围内的体素，这样内存占用与孔隙数量成正比，而不是与网格总数成正比。

5.5 孔径分布失真：小孔太多或大孔太少

现象：用对数正态分布时，生成的孔径直方图和理论曲线对不上。原因通常是lognrnd的参数换算写错，或者r_mean和r_std的物理含义混淆。解决方法是先用histfit画一下生成半径的分布，和理论 PDF 对比。如果偏差大，检查mu和sigma的换算公式——mu是半径对数的均值，不是半径均值，这两个差一个指数关系，很容易搞混。

6. 二次开发：把随机孔隙生成接进自己的仿真流程

6.1 批量生成与参数扫描

单次生成一个模型只能做定性观察，真正做研究需要批量生成不同孔隙率、不同孔径的样本，然后统计渗流或力学性能。最直接的做法是把main.m改成一个函数，接受参数结构体，返回孔隙率和导出路径，然后在外层用parfor并行扫描。

function run_batch() porosities = 0.15:0.05:0.45; r_means = [5, 8, 12]; results = []; for p = porosities for r = r_means params = struct('L', 100, 'target_porosity', p, ... 'r_mean', r, 'r_std', r*0.2, ... 'dist_type', 'lognormal'); [~, actual_porosity] = generate_pores(params); fname = sprintf('output/model_p%.2f_r%d.vtk', p, r); export_mesh(build_voxel( ... generate_pores(params), params.L), fname); results = [results; p, r, actual_porosity]; end end writematrix(results, 'output/batch_results.csv'); end

这段代码的关键是把generate_pores调用两次的问题——第一次拿孔隙列表，第二次拿实际孔隙率。更好的写法是让generate_pores返回两个输出，避免重复计算。parfor替换for可以并行加速，但要注意每个 worker 的随机数种子要独立设置，否则生成的模型会重复。参数说明：porosities和r_means是你要扫描的变量，results记录实际孔隙率用于后续校验。

6.2 从独立孔隙到连通孔隙

默认生成的孔隙是互相独立的球体，但真实多孔材料往往是连通的。改造思路有两种：一是允许球体重叠，重叠部分自然形成连通通道；二是在球体之间额外生成圆柱形喉道。第一种改动小，只需要把overlap_check的返回值改成“允许部分重叠”，并加一个重叠比例上限。第二种更真实，但需要额外写一个generate_throats函数，在邻近孔隙对之间生成圆柱。

我一般会先用第一种快速验证连通性对性能的影响，如果效果明显，再投入时间做第二种。连通性判断可以用bwlabeln对体素矩阵做连通域标记，如果只有一个连通域，说明孔隙完全连通。

6.3 验证生成结果是否可信

生成完模型不能直接扔进仿真，要先做几个基本校验。第一，实际孔隙率和目标值的偏差是否在 5% 以内。第二，孔径分布的直方图是否和设定分布吻合。第三，用isosurface可视化一下，看有没有明显的孤立孔隙或异常大孔。第四，如果做渗流，可以先在体素矩阵上用简单的格子玻尔兹曼或达西定律估算渗透率，和文献值对比。

% 快速校验：连通域数量 CC = bwconncomp(voxel, 6); % 6 邻域连通 num_pores = CC.NumObjects; fprintf('连通域数量: %d\n', num_pores); if num_pores > 1 warning('存在孤立孔隙，检查重叠检测参数'); end

bwconncomp用 6 邻域判断连通，如果连通域数量大于 1，说明有孤立孔隙。这时候要么放宽重叠条件，要么手动删除小连通域。这个校验我每次生成新参数后都会跑一遍，比肉眼看图靠谱得多。

6.4 一个让我改了三版的教训

最早我图省事，把resolution和L设成固定值，结果换了一个材料体系后，孔径尺度变了，生成的模型要么全是孔要么全是实体。后来我强制自己每次改参数前先算一下L / r_mean的比值，确保它在 8 到 20 之间——低于 8 边界效应太强，高于 20 孔隙太小看不清。从那以后我每次生成新模型都强制走一遍“孔隙率校验 + 连通域检查 + 孔径直方图”这三步，再也没出现过跑完仿真才发现模型不对的情况。希望这个习惯也能帮到你。

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