I. Cơ sở mô phỏng lỗ rỗng ngẫu nhiên hai chiều

1. Thuật toán cốt lõi: Phương pháp xếp chồng hình cầu ngẫu nhiên

Thông qua việc sinh các hình cầu có vị trí và bán kính ngẫu nhiên, đảm bảo không chồng lấn, mô phỏng cấu trúc lỗ rỗng. Các tham số quan trọng gồm:

• Độ rỗng (Porosity): tỷ lệ thể tích lỗ rỗng
• Phân bố bán kính: phân bố chuẩn hoặc phân bố đều
• Khoảng cách an toàn: ngăn các hình cầu dính nhau khi chia lưới

Khung mã MATLAB:

function pores = generate2DPores(targetPorosity, maxRadius, domainLength)
    % targetPorosity: độ rỗng mục tiêu | maxRadius: bán kính lớn nhất
    % domainLength: cạnh miền | pores: [x, y, r]
    pores = zeros(0, 3);
    domainArea = domainLength^2;
    poreArea = 0;
    maxAttempts = 1e5;
    attempts = 0;

    while poreArea / domainArea < targetPorosity && attempts < maxAttempts
        attempts = attempts + 1;
        x = domainLength * rand();
        y = domainLength * rand();
        r = max(0.1, abs(normrnd(maxRadius / 2, maxRadius / 4)));

        if ~hasCollision2D(pores, x, y, r)
            pores(end + 1, :) = [x, y, r];
            poreArea = poreArea + pi * r^2;
        end
    end
end

function collision = hasCollision2D(pores, x, y, r)
    collision = false;
    for k = 1:size(pores, 1)
        dx = pores(k, 1) - x;
        dy = pores(k, 2) - y;
        if hypot(dx, dy) < (pores(k, 3) + r) * 1.05
            collision = true;
            return;
        end
    end
end

II. Mở rộng mô phỏng lỗ rỗng ngẫu nhiên ba chiều

1. Hướng cải tiến thuật toán

• Kiểm tra va chạm ba chiều: cần xét 26 hướng lân cận
• Tối ưu tính toán thể tích: dùng cây bát phân (octree) để tăng tốc truy vấn không gian
• Kiểm soát phân bố bán kính: hỗ trợ phân bố chuẩn cắt cụt (truncated normal)

Triển khai MATLAB:

function spheres = generate3DPores(targetPorosity, meanRadius, stdRadius, domainLength)
    % targetPorosity: độ rỗng mục tiêu | meanRadius: bán kính trung bình
    % stdRadius: độ lệch chuẩn | domainLength: cạnh miền
    spheres = zeros(0, 4); % [x, y, z, r]
    totalVolume = domainLength^3;
    poreVolume = 0;
    maxSpheres = 1e6;

    while poreVolume / totalVolume < targetPorosity && size(spheres, 1) < maxSpheres
        r = max(0.1, abs(normrnd(meanRadius, stdRadius)));
        pos = domainLength * rand(1, 3);

        if ~hasCollision3D(spheres, pos, r)
            spheres(end + 1, :) = [pos, r];
            poreVolume = poreVolume + (4 / 3) * pi * r^3;
        end
    end
end

function collision = hasCollision3D(spheres, pos, r)
    collision = false;
    for k = 1:size(spheres, 1)
        dist = norm(pos - spheres(k, 1:3));
        if dist < (spheres(k, 4) + r) * 1.1
            collision = true;
            return;
        end
    end
end

III. Các kỹ thuật tối ưu chính

1. Tăng tốc phân vùng không gian (thuật toán octree)

classdef OctreeNode
    properties
        bounds      % [xmin ymin zmin; xmax ymax zmax]
        children = []
        sphereIdx = []
    end

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

        function idx = query(obj, pos, r)
            idx = [];
            if isempty(obj.children)
                idx = obj.sphereIdx;
                return;
            end

            for k = 1:numel(obj.children)
                childBounds = obj.children(k).bounds;
                if any(pos + r < childBounds(1, :)) || ...
                   any(pos - r > childBounds(2, :))
                    continue;
                end
                idx = [idx; obj.children(k).query(pos, r)];
            end
        end
    end
end

2. Phương pháp trực quan hóa

% Trực quan hóa ba chiều (dùng vol3d)
vox = generateSphereVoxel(spheres, gridSize, domainLength);
vol3d('CData', vox, 'Parent', gca);
axis equal; colormap(gray); shading interp;

function vox = generateSphereVoxel(spheres, gridSize, domainLength)
    % Sinh dữ liệu voxel cho các hình cầu
    xAxis = linspace(0, domainLength, gridSize);
    [X, Y, Z] = ndgrid(xAxis, xAxis, xAxis);
    vox = zeros(gridSize, gridSize, gridSize);

    for k = 1:size(spheres, 1)
        center = spheres(k, 1:3);
        r = spheres(k, 4);
        mask = (X - center(1)).^2 + (Y - center(2)).^2 + (Z - center(3)).^2 < r^2;
        vox(mask) = 1;
    end
end

IV. Kiểm soát tham số và kiểm chứng

1. Bảng tham số chính

2. Phương pháp kiểm chứng

% Kiểm chứng độ rỗng
truePorosity = 0.65;
simPorosity = sum((4/3)*pi*spheres(:,4).^3) / domainLength^3;
porosityError = abs(truePorosity - simPorosity);

% Kiểm chứng tính liên thông
connComp = bwconncomp(vox);
if connComp.NumRegions > 1
    error('Môi trường xốp tồn tại vùng cô lập!');
end

V. Ví dụ ứng dụng kỹ thuật

1. Mô phỏng thấm trong đá

% Sinh môi trường xốp
spheres = generate3DPores(0.3, 0.2, 0.05, 100);

% Xuất sang mô hình COMSOL
comsol.model.geom('geom1').create('box1', 'Box');
comsol.model.geom('geom1').feature('box1').set('size', [100, 100, 100]);
for k = 1:size(spheres, 1)
    comsol.model.geom('geom1').feature('sphere', 'create', ...
        'pos', spheres(k, 1:3), 'radius', spheres(k, 4));
end
comsol.model.geom('geom1').feature('subtract').set('input', {'box1', 'sphere'});

2. Mô hình hóa tầng chứa dầu khí

• Thiết lập tham số: độ rỗng 0.2-0.4, phân bố bán kính lệch (phản ánh môi trường trầm tích)
• Hậu xử lý: tính tensor độ thấm

% Tính độ thấm (công thức Kozeny-Carman)
k = (phi^3 * d^2) / (180 * (1 - phi)^2); % d là đường kính lỗ đặc trưng

VI. Giải pháp cho các vấn đề thường gặp

1. Hiệu suất tính toán thấp

   • Sử dụng phân vùng không gian octree (tăng tốc > 50 lần)
   • Tính toán song song (thay vòng lặp for bằng parfor)

2. Lỗ rỗng không liên thông

   • Thêm thuật toán sửa tính liên thông:

   function spheres = ensureConnectivity(spheres, domainLength)
       % Dùng BFS để nối các vùng cô lập
       visited = false(size(spheres, 1), 1);
       queue = 1;
       visited(queue) = true;
   
       while ~isempty(queue)
           current = queue(1);
           queue(1) = [];
           neighbors = findNeighbors(spheres, current);
           for idx = neighbors
               if ~visited(idx)
                   visited(idx) = true;
                   queue(end + 1) = idx; %#ok<AGROW>
               end
           end
       end
   
       spheres(visited == false, :) = [];
   end
   

3. Tích tụ ở biên

   • Dùng phương pháp phản xạ gương:

   function pos = reflectBoundary(pos, r, domainLength)
       for d = 1:3
           if pos(d) - r < 0
               pos(d) = 2*r - pos(d);
           elseif pos(d) + r > domainLength
               pos(d) = 2*domainLength - 2*r - pos(d);
           end
       end
   end
   

VII. Hướng ứng dụng mở rộng

1. Ghép nối đa trường vật lý
   • Mô phỏng ghép nhiệt - dòng - rắn (cần xuất sang mô hình COMSOL)
   • Tiến hóa lỗ rỗng do phản ứng hóa học (thêm biến trường pha)
2. Hỗ trợ học máy
   • Dùng GAN sinh cấu trúc lỗ rỗng phức tạp
   • Trích xuất đặc trưng lỗ rỗng dựa trên CNN
3. Kiểm chứng thực nghiệm
   • Kiểm chứng vi cấu trúc bằng in 3D
   • Đối chiếu với chụp cắt lớp X-quang