有限扩散集团凝聚模型(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