有限扩散集团凝聚模型(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/Na,公式中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 % %%%%%%%%%%%%%%%%%%% 版权所有,仅供学习交流,严禁用于商业用途,谢谢合作%%%%%%%%% clc clear %% 输入模块(初始化) 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 %步长加1 end
统计后的代码
%%%%%%%%%%%%%%%%%%% 作者:牟天蔚 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%% 欢迎访问本人博客:http://blog.sina.com.cn/u/1752049725 % %%%%%%%%%%%%%%%%%%% 版权所有,仅供学习交流,严禁用于商业用途,谢谢合作 %%%%% clc clear %% 输入模块(初始化) 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 %步长加1 end %% 制图模块 %%%%粒子数目-回转半径%%%%%%%%%%%%%%%%%%%%%%% figure hold on temp_plot=sortrows([Rnum_lizi Rmn],1); plot(temp_plot(:,1),temp_plot(:,2),'*'); grid on xlabel('团簇中所含有的粒子数目'); ylabel('团簇的最大回转半径'); title('粒子数目随回转半径变化'); hold off %%%%%团簇数目随絮凝时间的变化%%%%%%%%%%%%%%% figure hold on plot(1:st,Bh); grid on xlabel('絮凝时间'); ylabel('团簇数目'); title('团簇数目随絮凝时间的变化'); hold off %%%%%仅含单个粒子的团簇数目随絮凝时间的变化%%%%%%%%%%%%%%% figure hold on plot(1:st,single_Num); grid on xlabel('絮凝时间'); ylabel('仅含单个粒子的团簇数目'); title('仅含单个粒子的团簇数目随絮凝时间的变化'); hold off %%%%%不同团簇中粒子数目随絮凝时间的变化%%%%%%%%%%%%%%%%%%%%%%%% figure hold on plot(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 on xlabel('絮凝时间'); 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); %找出各团簇中的粒子数目 end tongji_num=tongji_num(tongji_num~=0); for i=1:max(tongji_num) tongji_num1(i,1)=size(find(tongji_num==i),1); %选出每一个团中有过少个粒子 end tongji_num1=tongji_num1'; figure hold on grid on bar(1:max(tongji_num),tongji_num1); xlabel('团簇中所含有的粒子数'); ylabel('团簇数目'); title('粒子数目与团簇数的统计'); hold off %%%%%絮凝过程中絮体长度特征的变化%%%%%%%%%%%%%%%%%%? figure hold on plot(1:st,xunzhaoRmn_R); plot(1:st,xunzhaoRan_R); grid on xlabel('絮凝时间'); ylabel('团簇的尺寸'); legend('最大回转半径','平均回转半径'); title('絮凝过程中絮体长度特征的变化'); hold off %%%%%回转分形维数随絮凝时间的变化%%%%%%%%%%%%%%%? figure hold on cccccc=size(find(xuzhaoD_D>=1000 | xuzhaoD_D==0),1); plot(1:st-cccccc,xuzhaoD_D(cccccc+1:end,:)); grid on xlabel('絮凝时间'); ylabel('平均回转分形维数'); title('回转分形维数随絮凝时间的变化'); hold off %%%%%絮凝过程中团簇的孔隙率变化%%%%%%%%%%%%%%%? figure hold on plot(1:st,Sr_avg_S); grid on xlabel('絮凝时间'); ylabel('平均空隙率'); title('絮凝过程中团簇的孔隙率变化'); hold off %%%%%粒子分形维度%%%%%%%%%%%%%%%? figure hold on plot(1:st,xuzhaoD2_D); grid on xlabel('絮凝时间'); ylabel('平均粒子分形维数'); title('粒子分形维度'); hold off %%%%%粒子分形维度计算图%%%%%%%%%%%%%%%? figure hold on plot(xx1,yy1); grid on xlabel('ln(n)'); ylabel('ln(N(n)'); title('粒子分形维度计算图'); hold off %%%%%单个团簇所含粒子数的变化过程%%%%%%%%%%%%%%%? figure hold on plot(1:st,Tuanchu_s_tongji); grid on xlabel('絮凝时间'); ylabel('团簇中所含有的粒子数'); title('单个团簇所含粒子数的变化过程'); hold off