文末有完整代码分享链接
文件介绍
automain 为元胞自动机主函数
choosedirection 选择方向函数，主函数调用
judgedirection 判断位置函数，主函数调用
neighbor 求每个元胞的邻居函数，主函数调用
surfaceness 求表面粗糙度
porosity 求孔隙率
flux 求水通量
intercept 求盐截留率
matrixplot 可视化图形
关键代码
while num<=98
%tempmatrix = simulation; %定义新的矩阵变量暂时保存当前画面
for i=1:rm
for j=1:rn
switch simulation(i,j)
case -1 %边界元胞，不作处理
case 0 %膜孔元胞，不作处理
case 1 %哌嗪元胞
if ismember(simulation(i,j),grouppip{i,j})
continue;
else
if ismember(3,nextcell{i,j}) %若元胞上面邻居含有元胞3，则结合生成6
direction=find(nextcell{i,j}==3);
randindex=randperm(length(direction));
resultindex=direction(randindex(1)); %假设水分子运动为布朗运动
[swapi,swapj]=judgedirection(resultindex,i,j,rn);
simulation(i,j)=6; %哌嗪位被PA替代
simulation(swapi,swapj)=4; %TMC位被油相溶剂替代
active(i,j)=2; %假设pa只反应两次
elseif ismember(0,nextcell{i,j}) || ismember(2,nextcell{i,j}) || ismember(4,nextcell{i,j})
direction=union(find(nextcell{i,j}==0),union(find(nextcell{i,j}==2),find(nextcell{i,j}==4))); %寻找四个方向中含有元胞0,1,2,4的方向
p_max=1;p_min=0.5; %概率变化范围
if i<=m1 %若哌嗪在水相
p_i=p_max-(i-1)*(p_max-p_min)/m1;
resultindex=choosedirection(direction,2,p_i);
%resultindex=2;
elseif i>m1+m2 %若哌嗪在油相
p_i=p_min+(i-m1-m2)*(p_max-p_min)/m3;
resultindex=choosedirection(direction,1,p_i);
else %若哌嗪在界面中，等概率扩散
%p_i=p_max-(i-1)*(p_max-p_min)/m1;
resultindex=choosedirection(direction,2,1);
%resultindex=2;
end
%由于界面聚合，水相哌嗪在界面处移动概率最小，两边概率最大，概率随着行数逐渐变化,依据概率选择最佳方向
[swapi,swapj]=judgedirection(resultindex,i,j,rn); %获得选定方向元胞的行列号
if simulation(swapi,swapj)==0 %如果扩散到膜孔,膜孔元胞置为哌嗪元胞，哌嗪元胞变成水分子元胞
simulation(swapi,swapj)=1;
simulation(i,j)=2;
elseif simulation(swapi,swapj)==4
if swapi<=m1+m2
continue;
elseif swapi>m1+m2+min(size(find(simulation==1)),size(find(simulation==3)))/n+1 %哌嗪只能活跃于油相的表层,该层为期望膜厚度
continue;
end
else
temp=simulation(i,j); %交换元胞状态
simulation(i,j)=simulation(swapi,swapj);
simulation(swapi,swapj)=temp;
end
elseif ismember(1,nextcell{i,j}) %若邻居有自身，则依玻尔兹曼因子可能形成团簇
if rand>=w_pip && i>2 && i<rm-1 && j>2 && j<rn-1 %弹性碰撞
direction=find(nextcell{i,j}==1);
randindex=randperm(length(direction));
resultindex=direction(randindex(1)); %假设水分子运动为布朗运动
[swapi,swapj]=judgedirection(resultindex,i,j,rn);
temp=simulation(2*i-swapi,2*j-swapj); %交换元胞状态
simulation(2*i-swapi,2*j-swapj)=simulation(i,j);
simulation(i,j)=temp;
else %形成团簇
grouppip{i,j}=[simulation(i,j),simulation(swapi,swapj)];
end
end
end
case 2 %水分子元胞
if ismember(0,nextcell{i,j}) %若水分子扩散到膜孔
direction=find(nextcell{i,j}==0);
randindex=randperm(length(direction));
resultindex=direction(randindex(1)); %假设水分子运动为布朗运动
[swapi,swapj]=judgedirection(resultindex,i,j,rn);
if simulation(swapi,swapj)==0
simulation(swapi,swapj)=2;
end
end
case 3 %TMC元胞
if ismember(simulation(i,j),groupTMC{i,j})
continue;
else
if ismember(1,nextcell{i,j}) %若元胞邻居含有元胞1，则结合生成6
direction=find(nextcell{i,j}==1);
randindex=randperm(length(direction));
resultindex=direction(randindex(1));
[swapi,swapj]=judgedirection(resultindex,i,j,rn);
simulation(i,j)=4;
simulation(swapi,swapj)=6; %哌嗪位被PA替代
active(swapi,swapj)=2;
elseif ismember(4,nextcell{i,j})
direction=find(nextcell{i,j}==4);
p_max=1;p_min=0.5;
p_i=p_min+(i-m1-m2)*(p_max-p_min)/m3;%由于界面聚合，油相哌嗪要往界面去的概率更大，假设概率为0.7
resultindex=choosedirection(direction,1,p_i); %1表示向上，依据概率选择最佳方向
[swapi,swapj]=judgedirection(resultindex,i,j,rn);
if (simulation(swapi,swapj)==4 && swapi>m1+m2) || simulation(swapi,swapj)==3
temp=simulation(i,j);
simulation(i,j)=simulation(swapi,swapj);
simulation(swapi,swapj)=temp;
end
elseif ismember(3,nextcell{i,j})
if rand>=w_TMC && i>2 && i<rm-1 && j>2 && j<rn-1 %弹性碰撞
direction=find(nextcell{i,j}==3);
randindex=randperm(length(direction));
resultindex=direction(randindex(1)); %假设水分子运动为布朗运动
[swapi,swapj]=judgedirection(resultindex,i,j,rn);
temp=simulation(2*i-swapi,2*j-swapj); %交换元胞状态
simulation(2*i-swapi,2*j-swapj)=simulation(i,j);
simulation(i,j)=temp;
else %形成团簇
groupTMC{i,j}=[simulation(i,j),simulation(swapi,swapj)];
end
end
end
case 4 %HEX油相溶剂元胞
case 5 %基膜元胞
case 6 %PA元胞
p1=0.4286;p2=0.5714; %经计算，p1表示6和1反应再生成6的概率，p2表示6和3反应再生成6的概率
if active(i,j)>0
if ismember(1,nextcell{i,j}) && ismember(3,nextcell{i,j}) %如果周围仍然有单体
flag=(p1*rand>p2*rand);
if flag==1 %和1发生反应
direction=find(nextcell{i,j}==1);
randindex=randperm(length(direction));
resultindex=randindex(1);
[swapi,swapj]=judgedirection(resultindex,i,j,rn);
simulation(swapi,swapj)=6;
active(swapi,swapj)=1;
active(i,j)=active(i,j)-1;
else %和3发生反应
direction=find(nextcell{i,j}==3);
resultindex=direction(1);
[swapi,swapj]=judgedirection(resultindex,i,j,rn);
simulation(swapi,swapj)=6;
active(swapi,swapj)=1;
active(i,j)=active(i,j)-1;
end
elseif ismember(1,nextcell{i,j}) || ismember(3,nextcell{i,j}) %若6既和1又和3相邻
direction=union(find(nextcell{i,j}==1),find(nextcell{i,j}==3));
resultindex=direction(1);
[swapi,swapj]=judgedirection(resultindex,i,j,rn);
simulation(swapi,swapj)=6;
active(swapi,swapj)=1;
active(i,j)=active(i,j)-1;
end
else
continue;
end
if i>m1+m2 %若生成的PA处于油相中，则应该以一定的规律附着在界面
direction=[1 2 3 4];
p_max=1;p_min=0.5;
p_i=p_min+(i-m1-m2)*(p_max-p_min)/m3;
resultindex=choosedirection(direction,1,p_i);
[swapi,swapj]=judgedirection(resultindex,i,j,rn);
temp=simulation(i,j);
if simulation(swapi,swapj) ==0 || simulation(swapi,swapj)==5
continue;
elseif swapi<=m1+m2 && simulation(swapi,swapj)==3
continue;
elseif swapi>=m1-1 && simulation(swapi,swapj)==2
continue;
else
simulation(i,j)=simulation(swapi,swapj);
simulation(swapi,swapj)=temp;
end
%矩阵元素聚类
elseif i<m1
direction=[1 2 3 4];
p_max=1;p_min=0.5;
p_i=p_min-(i-m1)*(p_max-p_min)/m1;
resultindex=choosedirection(direction,2,p_i);
[swapi,swapj]=judgedirection(resultindex,i,j,rn);
temp=simulation(i,j);
if simulation(swapi,swapj) ==0 || simulation(swapi,swapj)==5
continue;
elseif swapi<=m1+m2 && simulation(swapi,swapj)==3
continue;
elseif swapi>=m1-1 && simulation(swapi,swapj)==2
continue;
else
simulation(i,j)=simulation(swapi,swapj);
simulation(swapi,swapj)=temp;
end
end
end
list=neighbor(i,j,simulation);
nextcell{i,j}=list;
end
end
[row,col]=find(simulation==6); %PA的电荷效应
len=length(row);
for i=1:len
if row(i)>m1+m2
active(row(i),col(i))=active(row(i),col(i))+1;
elseif row(i)<m1
active(row(i),col(i))=active(row(i),col(i))-1;
end
end效果展示
分享链接：