有限扩散集团凝聚模型第一讲： DLCA模型定义及MATLAB中的实现

有限扩散集团凝聚模型（DLCA）第一讲：

模型定义及MATLAB中的实现

作者：牟天蔚

 

话说有限扩散集团絮凝模型，它的外文名字叫做Diffusion Limited Cluster Aggregation（DLCA），听着名字好高大上，但这个名字说明不了它是干啥的。我们需要深入的讲解一下。

首先我给自动化专业的人讲解一个概念——凝聚，。在被污染的水中，污染物有很多种，有的是溶解在水中的（类似于盐水），有的是悬浮在水中的（类似水中的沙子），还有一些（重点）是即不溶于水的也不悬浮于水中的，这种东西叫做“胶体”，这种东西简直就是个坑，要说溶液吧，加热或者冷却一下东西就沉到水底了，悬浮物更不用说了，一张滤纸就直接搞定。恰恰是这个胶体，就是滚刀肉，你怎么弄都不行，而且这东西里面还有很多病毒，人一喝这水，那是说没就没啊。所以我们应该把这个胶体给干掉，怎么干掉，加药！加什么药，好多种，比如氯化铝就是一种。药怎么干掉它？就是把它给聚集到一块，然后变大，变大，变大，最后变成悬浮物了，最后一过滤就没啦！刚才说了，给它聚集到一块，这就是凝聚，用那群砖家狗的话说就是：胶体颗粒脱稳并形成微小聚集体的过程（说的叫神马，不就是粒子被吸到一块了么）。

其次，我要说一下有限扩散集团凝聚模型(DLCA)，这名字倒是够装X的了，你看他的英文名第一个字母D，扩散，扩散指的是啥呢，就是现在有胶体一个，这个胶体在水中可以来回的飞啊飞啊飞啊。字母L：L是有限的意思，就是飞啊飞啊飞啊，但是它不能乱飞，只能在我指定的空间里面飞翔。C指的是集团，就是说水里面不只一个胶体，有好多好多哦。A指的是凝聚，不解释了，上面已经说的很清楚了。用一张图解释如下：

 

之后，我们就要说一下DLCA模型的原理了，这部分有点难懂，所以我就按照一种大家容易理解的方式讲讲。最重要的我们要谨记一句箴言：走碰粘沉。意思是胶体在水中走动，然后碰上了，粘在一块，最后变成悬浮物沉降到水底。把走碰粘沉广义点说法，就是一个胶体游走，碰到另外一个后粘贴成一个粒子，这个粒子在继续游走，再碰上一个粒子，然后这两个又变一个。。。。。，最后把上面说的一个胶体游走变成许多胶体，你就想象吧，那么多东西来回乱撞，合体，再撞再合体。。。。。

最后，就到了我说的重点内容上面了，用MATLAB如何实现这个程序？第一步是输入一些参数，如初始的粒子个数，单个粒子的直径，这些粒子数减小到多少的时候停止程序，以及有限扩散的范围L（limited）。这就是说现在在一个长宽高为100的立方体里面，有5000个直径为1的粒子。

N=5000; %初始粒子数d0=1; %粒子的直径及粘附间距（中小球半径为0.5）nmin=2000; %结束时粒子集团总数L=100;  %设置立方体大小（L*L*L)

第二步，将这5000粒子的位置确定，我的确定方式是均匀分布，在这里我们先假设直径为0，然后把立方体划分为一个个体积为200的小立方体（共5000个），并让这5000个粒子放在分别放在这5000个立方体里面，并个每一个粒子编号1~5000。

xyz=[unifrnd(0,100,N,3),(1:N)']; %随机生成三维坐标（均匀分布）

第三步，每次粒子运动都是没有规律的，因此需要随机的改变位置，但是呢，实际中温度越高，粒子肯定是运动的越快，我们需要设定一个粒子运动的步长，也就是一步走多远，步长设置方法为：

S=r/N^a，公式中S为步长，r和a是常数，在程序里面，我取的是2和0.5.

para1=0.5;para2=2;step=para2/1^para1;

第四步就是产生随机数，然后把初始的粒子坐标加上，模拟粒子游走了。为了很好的模拟效果，我们知道sin(x)可以取到-1和1之间的任何值，所以说我们应该先产生一个随机数x让它在0~2pi之间(坐标系一圈，取到所有[-1,1]之间的值），带入sin()中，在乘上步长step，为了使得X,Y,Z轴走的长度不一致，那么我们应该多设置几个随机数让他们相乘，这样就可以模拟的更好，因此应该这么编写：

    theta=2*pi*rand(length(flag),1);gamma=pi*rand(length(flag),1);    for i=1:length(flag)        theta1(i)=theta(flag(i)); %flag指5000个粒子所在的位置是属于哪个团        gamma1(i)=gamma(flag(i));    end    xyz(:,1)=xyz(:,1)+step.*sin(gamma1).*cos(theta1);    xyz(:,2)=xyz(:,2)+step.*sin(gamma1).*sin(theta1);    xyz(:,3)=xyz(:,3)+step.*cos(gamma1);

第五步，就是开始不断的循环，直到到达结束时粒子集团总数为止。

第六步，统计絮凝的过程，在这里我不讲，下回书我会继续逐一的讲解。

得出的结果图是这个样子的

 

 

开始前（均匀分布）

结束后（2000团簇）

最后奉上代码，其一是非统计的代码，就是前五步；其二是统计的代码，就是包含第六步。

非统计的代码

%%%%%%%%%%%%%%%%%%% 作者：牟天蔚 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 欢迎访问本人博客：http://blog.sina.com.cn/u/1752049725 %%%%%%%%%%%%%%%%%%%% 版权所有，仅供学习交流，严禁用于商业用途，谢谢合作%%%%%%%%%clcclear%% 输入模块（初始化）N=5000; %初始粒子数d0=1; %粒子的直径及粘附间距（中小球半径为0.5）nmin=2000; %结束絮凝集团总数L=100;  %设置立方体大小（L*L*L)st=0; %运动次数xyz=[unifrnd(0,100,N,3),(1:N)']; %随机生成三维坐标（均匀分布）flag=(1:1:N)'; %生成初始团簇编号num=ones(size(flag,1),1);%初始化布朗运动步长para1=0.5;para2=2;step=para2/1^para1;%% 运算模块bh=unique(xyz(:,4));%团簇数目初始化d=[];  %最大回转半径记录gamma1=zeros(length(flag),1);theta1=zeros(length(flag),1);Bh=[]; %记录团簇数目变量tongji_lizishumu=[];while(length(bh)>=nmin)    theta=2*pi*rand(length(flag),1);    gamma=pi*rand(length(flag),1);    for i=1:length(flag)        tempxyz=xyz(flag==i,:);        if(~isempty(tempxyz))            x0=mean(tempxyz(:,1)); %重心计算            y0=mean(tempxyz(:,2)); %重心计算            z0=mean(tempxyz(:,3)); %重心计算r=sqrt((tempxyz(:,1)-x0).^2+(tempxyz(:,2)-y0).^2+(tempxyz(:,3)-z0).^2); %粒子到中心的距离        end        theta1(i)=theta(flag(i));        gamma1(i)=gamma(flag(i));    end    bh=unique(xyz(:,4));    Bh=[Bh;size(bh,1)];%记录团簇数目变量    xyz(:,1)=xyz(:,1)+step.*sin(gamma1).*cos(theta1);    xyz(:,2)=xyz(:,2)+step.*sin(gamma1).*sin(theta1);    xyz(:,3)=xyz(:,3)+step.*cos(gamma1);    for i=1:size(xyz,1)-1        for j=i+1:size(xyz,1)            if xyz(j,4)~=xyz(i,4)                d=sqrt((xyz(i,1)-xyz(j,1)).^2+(xyz(i,2)-xyz(j,2)).^2+(xyz(i,3)-xyz(j,3)).^2);                if d<=d0                    xyz(j,4)=xyz(i,4);                end            end        end    end    flag=xyz(:,4);    %找出最大的X_rmn个最大回转半径，并求平均值（X_rmn<=团簇数）    st=st+1 %步长加1end

 

统计后的代码

%%%%%%%%%%%%%%%%%%% 作者：牟天蔚 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 欢迎访问本人博客：http://blog.sina.com.cn/u/1752049725 %%%%%%%%%%%%%%%%%%%% 版权所有，仅供学习交流，严禁用于商业用途，谢谢合作 %%%%%clcclear%% 输入模块（初始化）N=5000; %初始粒子数d0=1; %粒子的直径及粘附间距（中小球半径为0.5）nmin=2000; %结束絮凝集团总数L=100;  %设置立方体大小（L*L*L)st=0; %运动次数Tuanchu_s=1;%单个团簇所含粒子数的变化过程%确定初始中心位置x_cen=L/2;y_cen=L/2;z_cen=L/2;xyz=[unifrnd(0,100,N,3),(1:N)']; %随机生成三维坐标（均匀分布）flag=(1:1:N)'; %生成初始团簇编号num=ones(size(flag,1),1);%重心初始位置x0=mean(xyz(:,1));y0=mean(xyz(:,2));z0=mean(xyz(:,3));r=sqrt((xyz(:,1)-x0).^2+(xyz(:,2)-y0).^2+(xyz(:,3)-z0).^2); %粒子到中心的距离Rm=max(r); %最大回转半径初始化Ra=mean(r); %平均回转半径Sr=1-((N*(0.5*d0)^3)/Rm^3); %团簇的孔隙率初始化%回转分形维数xx=log(sort(r));p1=polyfit(xx,log((1:N))',1);D=p1(1);%该时刻的粒度分型维数xx1=sort(r(r<=Rm));yy1=log((1:size(xx,1))');p2=polyfit(xx1,yy1,1);D2=p2(1);%初始化布朗运动步长para1=0.5;para2=2;step=para2/1^para1;X_rmn=30; %所含粒子数最多的个团簇重从大到小排列多少个，见99行%% 运算模块bh=unique(xyz(:,4));%团簇数目初始化d=[];  %最大回转半径记录gamma1=zeros(length(flag),1);theta1=zeros(length(flag),1);Bh=[]; %记录团簇数目变量single_Num=[];%单个粒子团簇数目初始化lizishu_double=[];lizishu_fourth=[];lizishu_sixth=[];lizishu_eighth=[];lizishu_tenth=[];xunzhaoRmn_R=[];xunzhaoRan_R=[];xuzhaoD_D=[];Sr_avg_S=[];xuzhaoD2_D=[];Tuanchu_s_tongji=[];xx1=[];yy1=[];Rmn=[]; %最大回转半径序列Rnum_lizi=[]; %各团簇粒子数序列Ran=[];D=[];DD=[];xxx1=[];yyy1=[];tongji_lizishumu=[];while(length(bh)>=nmin)    %num(i)=length(flag(flag==bh(i)));    %布朗运动    Rmn=[]; %最大回转半径序列    Rnum_lizi=[]; %各团簇粒子数序列    Ran=[];    D=[];    DD=[];    tongji_lizishumu=[];    theta=2*pi*rand(length(flag),1);    gamma=pi*rand(length(flag),1);    for i=1:length(flag)        xx1=[];        yy1=[];        tempxyz=xyz(flag==i,:);        if(~isempty(tempxyz))            x0=mean(tempxyz(:,1)); %重心计算            y0=mean(tempxyz(:,2)); %重心计算            z0=mean(tempxyz(:,3)); %重心计算            r=sqrt((tempxyz(:,1)-x0).^2+(tempxyz(:,2)-y0).^2+(tempxyz(:,3)-z0).^2); %粒子到中心的距离            Rm=max(r); %最大回转半径初始化            Rmn=[Rmn;Rm];%最大回转半径序列            Ra=mean(r); %平均回转半径            Ran=[Ran;Ra]; %平均回转半径序列            num_lizi=size(tempxyz,1); %粒子数初始化            Rnum_lizi=[Rnum_lizi;num_lizi]; %%各团簇粒子数序列            %该时刻的回转分型维数计算            tongji_lizishumu=[tongji_lizishumu;size(tempxyz,1)];            if(size(r,1)>2)                xx=log(sort(r));                p1=polyfit(xx,log(1:size(xx,1))',1);                D1=p1(1);                D=[D;D1];            elseif(size(r,1)==2)                 xx=log(sort(r));                 yy=log(1:size(xx,1))';                 D1=(yy(2)-yy(1))/(xx(2)-xx(1));                 D=[D;D1];            else                D1=0;                D=[D;D1];            end            %粒子分形维数计算            if max(tongji_lizishumu)==2 || max(tongji_lizishumu)==1                DD=[DD;0];            elseif max(tongji_lizishumu)==3                 tempxx1=log(2);                 tempyy1=log(size(find(tongji_lizishumu<=2),1));                 tempxx2=log(3);                 tempyy2=log(size(find(tongji_lizishumu<=3),1));                 D2=(tempyy2-tempyy1)/(tempxx2-tempxx1);                 DD=[DD;D2];                            else                for j=2:max(tongji_lizishumu)                    tempjisuan=j;                    tempxx1=log(j);                    xx1=[xx1;tempxx1];                    tempyy1=log(size(find(tongji_lizishumu<=j),1));                    yy1=[yy1;tempyy1];                end                p2=polyfit(xx1,yy1,1);                D2=p2(1);                DD=[DD;D2];                            end        end        theta1(i)=theta(flag(i));        gamma1(i)=gamma(flag(i));    end    xxx1=xx1;    yyy1=yy1;    bh=unique(xyz(:,4));    Bh=[Bh;size(bh,1)];%记录团簇数目变量    %单一粒子模拟    Tuanchu_s_num=size(find(xyz(:,4)==Tuanchu_s),1);    Tuanchu_s_tongji=[Tuanchu_s_tongji;Tuanchu_s_num];    xyz(:,1)=xyz(:,1)+step.*sin(gamma1).*cos(theta1);    xyz(:,2)=xyz(:,2)+step.*sin(gamma1).*sin(theta1);    xyz(:,3)=xyz(:,3)+step.*cos(gamma1);    for i=1:size(xyz,1)-1        for j=i+1:size(xyz,1)            if xyz(j,4)~=xyz(i,4)                d=sqrt((xyz(i,1)-xyz(j,1)).^2+(xyz(i,2)-xyz(j,2)).^2+(xyz(i,3)-xyz(j,3)).^2);                if d<=d0                    xyz(j,4)=xyz(i,4);                end            end        end    end    flag=xyz(:,4);    single_Num=[single_Num;size(find(Rnum_lizi==1),1)]; %单个粒子的数据集合    lizishu_double=[lizishu_double;size(find(Rnum_lizi>=2 & Rnum_lizi<4),1)];    lizishu_fourth=[lizishu_fourth;size(find(Rnum_lizi>=4 & Rnum_lizi<6),1)];    lizishu_sixth=[lizishu_sixth;size(find(Rnum_lizi>=6 & Rnum_lizi<8),1)];    lizishu_eighth=[lizishu_eighth;size(find(Rnum_lizi>=8 & Rnum_lizi<10),1)];    lizishu_tenth=[lizishu_tenth;size(find(Rnum_lizi>=10,1))];    %找出最大的X_rmn个最大回转半径，并求平均值（X_rmn<=团簇数）    tempRmn=[tongji_lizishumu Rmn];    tempRan=[tongji_lizishumu Ran];    xunzhaoRmn = sortrows(tempRmn,1);    xunzhaoRan = sortrows(tempRan,1);    xunzhaoRmn1=mean(xunzhaoRmn(end-X_rmn+1:end,2));    xunzhaoRan1=mean(xunzhaoRan(end-X_rmn+1:end,2));    xunzhaoRmn_R=[xunzhaoRmn_R;xunzhaoRmn1];    xunzhaoRan_R=[xunzhaoRan_R;xunzhaoRan1];    %平均回转分形维度    tempD=[tongji_lizishumu D];    xunzhaoD = sortrows(tempD,1);    xuzhaoD1=mean(xunzhaoD(end-X_rmn+1:end,2));    xuzhaoD_D=[xuzhaoD_D;xuzhaoD1];    %平均粒子分形维度    tempD2=[tongji_lizishumu DD];    xuzhaoD2=mean(find(~isnan(tempD2(:,2))));    xuzhaoD2_D=[xuzhaoD2_D;xuzhaoD2];    %平均孔隙率    Sr=1-((tongji_lizishumu.*(0.5*d0).^3)./(Rmn.^3)); %团簇的平均孔隙率    Sr_avg=mean(Sr(Sr>0));    Sr_avg_S=[Sr_avg_S;Sr_avg];    st=st+1 %步长加1end%% 制图模块%%%%粒子数目-回转半径%%%%%%%%%%%%%%%%%%%%%%%figurehold ontemp_plot=sortrows([Rnum_lizi Rmn],1);plot(temp_plot(:,1),temp_plot(:,2),'*');grid onxlabel('团簇中所含有的粒子数目');ylabel('团簇的最大回转半径');title('粒子数目随回转半径变化');hold off%%%%%团簇数目随絮凝时间的变化%%%%%%%%%%%%%%%figurehold onplot(1:st,Bh);grid onxlabel('絮凝时间');ylabel('团簇数目');title('团簇数目随絮凝时间的变化');hold off%%%%%仅含单个粒子的团簇数目随絮凝时间的变化%%%%%%%%%%%%%%%figurehold onplot(1:st,single_Num);grid onxlabel('絮凝时间');ylabel('仅含单个粒子的团簇数目');title('仅含单个粒子的团簇数目随絮凝时间的变化');hold off%%%%%不同团簇中粒子数目随絮凝时间的变化%%%%%%%%%%%%%%%%%%%%%%%%figurehold onplot(1:st,lizishu_double);plot(1:st,lizishu_fourth);plot(1:st,lizishu_sixth);plot(1:st,lizishu_eighth);plot(1:st,lizishu_tenth);legend('粒子数≥2','粒子数≥4','粒子数≥6','粒子数≥8','粒子数≥10');grid onxlabel('絮凝时间');ylabel('团簇的数目变化');title('不同团簇中粒子数目随絮凝时间的变化');hold off%%%%%粒子数目与团簇数的统计%%%%%%%%%%%%%%%%%%%%%%%%tongji_num=zeros(length(flag),1);tongji_num1=[];for i=1:length(flag)    tongji_num(i)=size(find(xyz(:,4)==i),1); %找出各团簇中的粒子数目endtongji_num=tongji_num(tongji_num~=0);for i=1:max(tongji_num)    tongji_num1(i,1)=size(find(tongji_num==i),1); %选出每一个团中有过少个粒子endtongji_num1=tongji_num1';figurehold ongrid onbar(1:max(tongji_num),tongji_num1);xlabel('团簇中所含有的粒子数');ylabel('团簇数目');title('粒子数目与团簇数的统计');hold off%%%%%絮凝过程中絮体长度特征的变化%%%%%%%%%%%%%%%%%%？figurehold onplot(1:st,xunzhaoRmn_R);plot(1:st,xunzhaoRan_R);grid onxlabel('絮凝时间');ylabel('团簇的尺寸');legend('最大回转半径','平均回转半径');title('絮凝过程中絮体长度特征的变化');hold off%%%%%回转分形维数随絮凝时间的变化%%%%%%%%%%%%%%%？figurehold oncccccc=size(find(xuzhaoD_D>=1000 | xuzhaoD_D==0),1);plot(1:st-cccccc,xuzhaoD_D(cccccc+1:end,:));grid onxlabel('絮凝时间');ylabel('平均回转分形维数');title('回转分形维数随絮凝时间的变化');hold off%%%%%絮凝过程中团簇的孔隙率变化%%%%%%%%%%%%%%%？figurehold onplot(1:st,Sr_avg_S);grid onxlabel('絮凝时间');ylabel('平均空隙率');title('絮凝过程中团簇的孔隙率变化');hold off%%%%%粒子分形维度%%%%%%%%%%%%%%%？figurehold onplot(1:st,xuzhaoD2_D);grid onxlabel('絮凝时间');ylabel('平均粒子分形维数');title('粒子分形维度');hold off%%%%%粒子分形维度计算图%%%%%%%%%%%%%%%？figurehold onplot(xx1,yy1);grid onxlabel('ln(n)');ylabel('ln(N(n)');title('粒子分形维度计算图');hold off%%%%%单个团簇所含粒子数的变化过程%%%%%%%%%%%%%%%？figurehold onplot(1:st,Tuanchu_s_tongji);grid onxlabel('絮凝时间');ylabel('团簇中所含有的粒子数');title('单个团簇所含粒子数的变化过程');hold off