有限擴散集團凝聚模型第一講： 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/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