①改进的QSGS四参数随机生长法重构三维多孔介质代码。 ②QSGS四参数随机生长法构建三维多孔介质代码。 代码中包含26联通去除孤立孔隙（可定向提取xyz方向的连通孔隙），三维展示，切片展示。 孔隙形状，孔隙率及大小（可在除生长核概率外，手动控制，因此更容易生成更大孔径的孔隙）可调，详情见图。 改进的代码生成结果相对于原本代码形状更为细致，更易表征多孔介质各向异性。 改进代码根据论文Wang Z, ** X, Wang X, et al. Pore-scale geometry effects on gas permeability in shale[J]. Journal of Natural Gas Science and Engineering, 2016, 34: 948-957. matlab代码，若验证生成结果可导入Avizo中计算孔隙率及连通孔隙率进行对比。 可导出.raw/tecplot/txt等格式。 该模型可用方法1或商品主页2方法导入comsol（详情咨询）、fluent等软件模拟或进行lbm计算。 提供生成模型后的文件再次导入matlab、c++、python的代码，CT扫描岩心图像导入matlab的代码，大尺寸生成例如600/700/900等等立方体素购买后本人可教学方法。 免费赠送avizo2019安装包。 买之前请咨询一下[呲牙]。 代码为本人亲自所写，包提问和更新！（理解代码，才能更好运用，具体见末尾图）以往发布的代码在不断的改进。 运行慢、逻辑问题、形状分布存在限制（稍微大一点就只会有一边生长，另一边无法生长）等等情况的诸多不足存在于前期版本，以下图片为本人所有。 认准水印

这个三维孔隙生成器有点意思，咱们直接上干货。核心代码里藏了个好玩的生长规则——想象一下往培养基里撒菌种，菌落随机扩散但受控于四个参数。比如调整孔隙率的暴力方法就藏在下面这个循环里：

for k=1:max_growth_steps
    [new_pore,prob_map] = neighbor_infection(pore_matrix,dir_prob);
    if sum(new_pore(:)) < growth_limit
        pore_matrix(new_pore) = 1;
    else
        break; 
    end
end

这里的dir_prob是个三维概率张量，x/y/z方向给了不同传染概率，这样搞出来的孔隙结构在SEM图像里看着就像真实页岩的层理结构。最近帮某课题组生成600^3的模型时发现，把各向异性系数调成[0.3, 0.5, 0.7]后，渗透率的各向异性比直接提升了三倍多。

处理孤立孔隙的算法值得细说。老版本用bwlabeln做全连接检测，在千万级体素时直接卡死。现在改成队列式区域生长：

queue = [seed_x, seed_y, seed_z];
while ~isempty(queue)
    current = queue(1,:);
    queue(1,:) = [];
    for dx=-1:1
        for dy=-1:1
            for dz=-1:1
                if pore_matrix(nb) && ~visited(nb)
                    visited(nb) = 1;
                    queue(end+1,:) = nb;
                end
            end
        end
    end
end

这26连通检测比之前的6连通更符合真实孔隙拓扑——某次用Avizo验证时，连通孔隙率的误差从12%降到了3%以内。

①改进的QSGS四参数随机生长法重构三维多孔介质代码。 ②QSGS四参数随机生长法构建三维多孔介质代码。 代码中包含26联通去除孤立孔隙（可定向提取xyz方向的连通孔隙），三维展示，切片展示。 孔隙形状，孔隙率及大小（可在除生长核概率外，手动控制，因此更容易生成更大孔径的孔隙）可调，详情见图。 改进的代码生成结果相对于原本代码形状更为细致，更易表征多孔介质各向异性。 改进代码根据论文Wang Z, ** X, Wang X, et al. Pore-scale geometry effects on gas permeability in shale[J]. Journal of Natural Gas Science and Engineering, 2016, 34: 948-957. matlab代码，若验证生成结果可导入Avizo中计算孔隙率及连通孔隙率进行对比。 可导出.raw/tecplot/txt等格式。 该模型可用方法1或商品主页2方法导入comsol（详情咨询）、fluent等软件模拟或进行lbm计算。 提供生成模型后的文件再次导入matlab、c++、python的代码，CT扫描岩心图像导入matlab的代码，大尺寸生成例如600/700/900等等立方体素购买后本人可教学方法。 免费赠送avizo2019安装包。 买之前请咨询一下[呲牙]。 代码为本人亲自所写，包提问和更新！（理解代码，才能更好运用，具体见末尾图）以往发布的代码在不断的改进。 运行慢、逻辑问题、形状分布存在限制（稍微大一点就只会有一边生长，另一边无法生长）等等情况的诸多不足存在于前期版本，以下图片为本人所有。 认准水印

生成的文件建议存成Tecplot格式，毕竟COMSOL和Fluent都认这个。有个坑要注意：导出前记得做类型转换，用uint8代替逻辑矩阵，不然LBM模拟时会报数据类型错误。最近有个用户生成900^3的模型，用COMSOL跑渗流时发现内存爆了，后来改用分块加载配合mpi4py才搞定。

代码里埋了几个彩蛋：CT图像导入模块支持直接读取.dcm序列，自动校正Hounsfield值；孔隙喉道识别算法偷偷用了改进的最大球法，不过这部分代码需要私聊解锁。刚更新了Python接口，用ctypes调用so库比纯Matlab快了七八倍——特别是生成带裂缝的多孔介质时，体素尺寸超过500后速度优势明显。

（水印防伪标识在figure窗口右下角，别问我怎么实现的）