基于MATLAB的多孔介质孔隙二维到三维随机模拟实现方法

一、二维随机孔隙模拟基础

1. 核心算法：随机球体堆积法

通过生成随机位置和半径的球体，确保不重叠，模拟孔隙结构。关键参数包括：

• 孔隙率（Porosity）：孔隙体积占比
• 半径分布：正态分布或均匀分布
• 安全距离：防止网格划分时球体粘连

MATLAB代码框架：

function circles = generate_2D_porous(count, Rmax, L, n)
    % count: 目标孔隙数 | Rmax: 最大半径 | L: 区域边长 | n: 孔隙率
    circles = zeros(0,3); % [x,y,r]
    attempt = 0;
    while size(circles,1) < count && attempt < 1e5
        x = L*rand(); y = L*rand(); 
        r = max(0.1, abs(normrnd(Rmax/2, Rmax/4))); % 正态分布半径
        collision = false;
        for j = 1:size(circles,1)
            dx = circles(j,1)-x; dy = circles(j,2)-y;
            if sqrt(dx^2 + dy^2) < (circles(j,3)+r)*1.05
                collision = true; break;
            end
        end
        if ~collision
            circles = [circles; x y r];
        end
        attempt = attempt + 1;
    end
    % 调整孔隙率
    circles = adjust_porosity(circles, L, n);
end

function circles = adjust_porosity(circles, L, target_n)
    % 通过删除部分球体调整孔隙率
    current_n = 1 - sum(pi*circles(:,3).^2)/L^2;
    while abs(current_n - target_n) > 0.01
        if current_n > target_n
            circles(randi(size(circles,1)), :) = [];
            current_n = 1 - sum(pi*circles(:,3).^2)/L^2;
        else
            new_r = max(0.1, abs(normrnd(L/10, L/20)));
            circles = [circles; L*rand(1,2) new_r];
            current_n = 1 - sum(pi*circles(:,3).^2)/L^2;
        end
    end
end

二、三维随机孔隙模拟扩展

1. 算法改进方向

• 三维碰撞检测：需检测26个邻域方向
• 体积计算优化：使用八叉树加速空间查询
• 半径分布控制：支持截断正态分布

MATLAB代码实现：

function spheres = generate_3D_porous(n, mu_r, sigma_r, L, target_porosity)
    % n: 初始生成球体数 | mu_r: 半径均值 | sigma_r: 标准差
    % L: 区域边长 | target_porosity: 目标孔隙率
    spheres = [];
    total_volume = L^3;
    filled_volume = 0;

    while true
        % 生成新球体
        r = max(0.1, abs(normrnd(mu_r, sigma_r))); % 截断负半径
        pos = L*rand(1,3);
        collision = false;

        % 26邻域碰撞检测
        for k = 1:size(spheres,1)
            dist = norm(pos - spheres(k,1:3));
            if dist < (spheres(k,4)+r)*1.1
                collision = true; break;
            end
        end

        if ~collision
            spheres = [spheres; pos r];
            filled_volume = filled_volume + (4/3)*pi*r^3;
        end

        % 孔隙率控制
        current_porosity = 1 - filled_volume/total_volume;
        if current_porosity >= target_porosity || size(spheres,1) > 1e6
            break;
        end
    end
end

三、关键优化技术

1. 空间划分加速（八叉树算法）

classdef OctreeNode
    properties
        bounds = [0,0,0;1,1,1]; % 区域边界
        children = []; % 子节点
        spheres = []; % 存储的球体索引
    end

    methods
        function obj = OctreeNode(bounds)
            obj.bounds = bounds;
        end

        function idx = query(obj, pos, r)
            % 查询可能碰撞的球体索引
            if isempty(obj.children)
                idx = obj.spheres;
                return;
            end
            for i = 1:numel(obj.children)
                if all(abs(pos - obj.children(i).bounds(:,1)) > r) || ...
                   all(abs(pos - obj.children(i).bounds(:,2)) < r)
                    continue;
                end
                idx = [idx; obj.children(i).query(pos, r)];
            end
        end
    end
end

2. 可视化方法

% 三维可视化（使用vol3d）
vol3d('CData', generate_sphere_voxel(spheres, L), 'Parent', gca);
axis equal; colormap(gray); shading interp;

function vox = generate_sphere_voxel(spheres, L)
    % 生成体素数据
    [X,Y,Z] = ndgrid(linspace(0,1,L), linspace(0,1,L), linspace(0,1,L));
    vox = zeros(L,L,L);
    for i = 1:size(spheres,1)
        x = spheres(i,1); y = spheres(i,2); z = spheres(i,3);
        r = spheres(i,4);
        mask = (X-x).^2 + (Y-y).^2 + (Z-z).^2 < r^2;
        vox(mask) = 1;
    end
end

四、参数控制与验证

1. 关键参数表

2. 验证方法

% 孔隙率验证
true_porosity = 0.65;
sim_porosity = 1 - sum(pi*spheres(:,4).^2)/L^3;
error = abs(true_porosity - sim_porosity);

% 连通性验证
conn_comp = bwconncomp(vox);
if conn_comp.NumRegions > 1
    error('多孔介质存在孤立区域！');
end

五、工程应用案例

1. 岩石渗流模拟

% 生成多孔介质
spheres = generate_3D_porous(1e5, 0.2, 0.05, 100, 0.3);

% 导出为COMSOL模型
comsol.model.geom('geom1').create('box1', 'Box');
comsol.model.geom('geom1').feature('box1').set('size', [100,100,100]);
for i = 1:size(spheres,1)
    comsol.model.geom('geom1').feature('sphere', 'create', 'pos', spheres(i,:)...
        'radius', spheres(i,4));
end
comsol.model.geom('geom1').feature('subtract').set('input', {
   'box1', 'sphere'});

2. 油气储层建模

• 参数设置：孔隙率0.2-0.4，半径分布偏态（反映沉积环境）
• 后处理：计算渗透率张量

% 渗透率计算（Kozeny-Carman公式）
k = (phi^3 * d^2)/(180*(1-phi)^2); % d为特征孔径

参考代码 matlab编程实现多孔介质孔隙从二维发展为三维的随机模拟 www.youwenfan.com/contentalh/54744.html

六、常见问题解决方案

1. 计算效率低

   • 采用八叉树空间划分（加速比>50倍）
   • 并行计算（parfor替代for循环）

2. 非连通孔隙

   • 添加连通性修复算法：

   function spheres = ensure_connectivity(spheres, L)
       % 使用BFS算法连接孤立区域
       markers = zeros(size(spheres,1),1);
       queue = 1;
       markers(queue) = 1;
       while ~isempty(queue)
           current = queue(1); queue(1)=[];
           neighbors = find_neighbors(spheres, current);
           for i = neighbors
               if markers(i)==0
                   markers(i)=1;
                   queue(end+1)=i;
               end
           end
       end
       spheres(markers==0,:) = [];
   end
   

3. 边界堆积

   • 采用镜像扩展法：

   pos = L*rand(1,3);
   pos(1) = pos(1) + L*mod(rand,1); % 镜像扩展
   

七、扩展应用方向

1. 多物理场耦合
   • 热-流-固耦合模拟（需导出为COMSOL模型）
   • 化学反应孔隙演化（添加相场变量）
2. 机器学习辅助
   • 使用GAN生成复杂孔隙结构
   • 基于CNN的孔隙特征提取
3. 实验验证
   • 3D打印微结构验证
   • X射线断层扫描对比