文末有完整代码分享链接
文件介绍
automain 为元胞自动机主函数
choosedirection 选择方向函数，主函数调用
judgedirection 判断位置函数，主函数调用
neighbor 求每个元胞的邻居函数，主函数调用
surfaceness 求表面粗糙度
porosity 求孔隙率
flux 求水通量
intercept 求盐截留率
matrixplot 可视化图形
关键代码
while num<=98%tempmatrix = simulation; %定义新的矩阵变量暂时保存当前画面for i=1:rmfor j=1:rnswitch simulation(i,j)case -1 %边界元胞，不作处理case 0 %膜孔元胞，不作处理case 1 %哌嗪元胞if ismember(simulation(i,j),grouppip{i,j})continue;elseif ismember(3,nextcell{i,j}) %若元胞上面邻居含有元胞3，则结合生成6direction=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)==4if swapi<=m1+m2continue;elseif swapi>m1+m2+min(size(find(simulation==1)),size(find(simulation==3)))/n+1 %哌嗪只能活跃于油相的表层,该层为期望膜厚度continue;endelsetemp=simulation(i,j); %交换元胞状态simulation(i,j)=simulation(swapi,swapj);simulation(swapi,swapj)=temp;endelseif 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)];endendendcase 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)==0simulation(swapi,swapj)=2;endendcase 3 %TMC元胞if ismember(simulation(i,j),groupTMC{i,j})continue;elseif ismember(1,nextcell{i,j}) %若元胞邻居含有元胞1，则结合生成6direction=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.7resultindex=choosedirection(direction,1,p_i); %1表示向上，依据概率选择最佳方向[swapi,swapj]=judgedirection(resultindex,i,j,rn);if (simulation(swapi,swapj)==4 && swapi>m1+m2) || simulation(swapi,swapj)==3temp=simulation(i,j);simulation(i,j)=simulation(swapi,swapj);simulation(swapi,swapj)=temp;endelseif 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)];endendendcase 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)>0if 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;endelseif 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;endelsecontinue;endif 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)==5continue;elseif swapi<=m1+m2 && simulation(swapi,swapj)==3continue;elseif swapi>=m1-1 && simulation(swapi,swapj)==2continue;elsesimulation(i,j)=simulation(swapi,swapj);simulation(swapi,swapj)=temp;end%矩阵元素聚类elseif i<m1direction=[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)==5continue;elseif swapi<=m1+m2 && simulation(swapi,swapj)==3continue;elseif swapi>=m1-1 && simulation(swapi,swapj)==2continue;elsesimulation(i,j)=simulation(swapi,swapj);simulation(swapi,swapj)=temp;endendendlist=neighbor(i,j,simulation);nextcell{i,j}=list;endend[row,col]=find(simulation==6); %PA的电荷效应len=length(row);for i=1:lenif row(i)>m1+m2active(row(i),col(i))=active(row(i),col(i))+1;elseif row(i)<m1active(row(i),col(i))=active(row(i),col(i))-1;endend效果展示
分享链接：
M00198-MATLAB界面聚合的元胞自动机模拟完整实现运行
