二维随机生长法-原创

测试案例

% demo1
% node1 = genQSGS2d([100 100],0.3,0.01,[0.2,0.1]);
% demo2
prob = [0.1 0.1 0.1 0.1 0.2 0.05 0.2 0.05]';
node2 = genQSGS2d([100 100],0.4,0.007,prob);

genQSGS2d.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 四参数随机生长算法
% 参考文献
% [1] Wang M ,  Wang J ,  Pan N , et al. Mesoscopic predictions of 
%     the effective thermal conductivity for microscale random porous media.[J]. 
%     Physical Review E Statistical Nonlinear & Soft Matter Physics, 2007, 75(3):036702.
% code @ 2021-12-21 All rights reserved.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [arrgrid] = genQSGS2d(box,Vs,cdd,D)
%GENQSGS2D 在给定区域QSGS
%   box -   盒子大小 参数形式 L或者[height,width]
%   Vs  -   固相体积分数
%   cdd -   随机生成种子概率
%   D   -   生长概率 参数形式为d或者[d1_4 d5_8]或者[d1 d2 d3 d4 d5 d6 d7 d8]
%    6   2   5
%      \ | /
%   3 —— 0 —— 1
%      / | \
%    7   4   8
% example
% node = genQSGS2d([100 100],0.3,0.01,[0.2,0.1]);
%% 输入参数设置
% 默认参数
if nargin < 4
    D = [0.2,0.1];
    if nargin <3
        cdd = 0.01;
        if nargin<2
            Vs = 0.5;
            if nargin<1
                box = [100,100];
            end
        end
    end
end
if length(box)==1
    len = box;wid = box;
elseif length(box) ==2
    len = box(1);wid = box(2);
else
    error('盒子大小参数输入有误！');
end

% 根据输入参数类型确定8个方向生长概率
if length(D)==1
    direcprob=D*ones(8,1); 
elseif length(D)==2
    direcprob = zeros(8,1);
    direcprob(1:4) = D(1);
    direcprob(5:8) = D(2);
elseif length(D)==8
    direcprob = D;
else
    error('生长概率参数D输入有误！');
end

if Vs<0||Vs>1
    error('固相体积分数Vs输入有误！');
end
if cdd<0
    error('初始种子概率cdd输入有误！')
end

%% 初始化参数
solid_need = round(Vs*len*wid); 
arrgrid = zeros(len,wid);
direc = [1,0,-1,0,1,-1,-1,1;0,1,0,-1,1,1,-1,-1]';

%% 初始生长核心
cdd_arr = rand(len,wid);
cdd_logic = (cdd_arr<cdd);

arrgrid(cdd_logic) = 1;
solid_num = sum(cdd_logic(:));  % 固体点数量
nsolid = solid_num;
solid_arr = find(cdd_logic==1); % 固体列表
cnt = 0;

flag = true;
%% 固相扩散生长
while solid_num <= solid_need && flag
    for id = 1:solid_num % 对每一个固相进行生长
        [sx,sy] = ind2sub(box,solid_arr(id));
        pos = [sx,sy] + direc;
        for ii=1:8
            if isinbox(pos(ii,:),box)
                if rand() < direcprob(ii) && arrgrid(pos(ii,1),pos(ii,2))==0
                    arrgrid(pos(ii,1),pos(ii,2))=1;
                    nsolid = nsolid +1;
                    solid_arr(nsolid) = sub2ind(box,pos(ii,1),pos(ii,2));               
                end
            end
        end
        if nsolid > solid_need %固体数量足够 退出for循环
            flag = false;
            break;
        end
    end
    solid_num = nsolid;
    volsolid = solid_num /len/wid;  % 体积分数
    cnt = cnt+1;
    
    if cnt > 1e5
        flag = false;
        disp('迭代次数过多，可能无法生成！');
    end
    if mod(cnt,2)==0
        fprintf('第%d次迭代，固相体积分数%f \n',cnt,volsolid);
    end
end

% 画灰度图
imagesc(arrgrid);colormap(gray);colorbar;

% 导出plt tecplot可读文件
if 0 == 1
    pltwrite('QSGS.plt',arrgrid);
end
end

function inbox = isinbox(ijk,box)
% 判断ijk是否在盒子内
    d = (ijk > 0) & ((box-ijk) >=0);
    inbox = (sum(d,2)==2);
end
function pltwrite(filename,data,hasHeader)
%PLTWRITE 输出单个矩阵到plt文件
if nargin==2
    hasHeader=1;
end
disp('开始输出plt文件');
fid = fopen(filename,'w+');
if length(size(data))==2
    [len,wid]=size(data); 
    if hasHeader
        fprintf(fid,'VARIABLES=\"x\",\"y\",\"issolid\"\n');
        fprintf(fid,'ZONE T=\"BOX\",I=%d,J=%d,F=POINT\n',len,wid);
    end
    for j=1:wid
        for i = 1:len
            fprintf(fid,'%d %d %d\n',i,j,data(i,j));
        end
    end
end
if length(size(data))==3
    [len,wid,hig]=size(data);
    if hasHeader
        fprintf(fid,'VARIABLES=\"x\",\"y\",\"z\",\"issolid\"\n');
        fprintf(fid,'ZONE T=\"BOX\",I=%d,J=%d,K=%d,F=POINT\n',len,wid,hig);
    end
    for k=1:hig
        for j=1:wid
            for i = 1:len
                fprintf(fid,'%d %d %d %d\n',i,j,k,data(i,j,k));
            end
        end
    end
end
fclose(fid);
disp('输出plt文件完毕');
end