matlab-程序-晶粒长大

2.Ln=200;%格点边长3.L=zeros(Ln);%格点矩阵4.Q=120;%总取向数5.step_num=500;%MC总步数6.interval_save_jpg=20;%图形存储间隔7.interval_stastics=2;%晶粒平均参数和相对密度统计间隔8.stastics_data=zeros(step_num/interval_stastics,5);%存储每interval_stastics次MCS后的平均晶粒尺寸和相对密度，存储格式为(MCS,graincount,averagearea,averagediameter,relativedensity)9.10.11.%烧结模拟过程参数赋值12.T=1;%温度参数13.J1=1;%晶界能14.15.16.%初始结构的格点赋值17.rand_l=randperm(Ln^2);18.fori=1:Ln^2*V_pore19.ifrem(rand_l(i),Ln)==020.x=Ln;21.y=fix(rand_l(i)/Ln);22.else23.x=rem(rand_l(i),Ln);24.y=fix(rand_l(i)/Ln)+1;25.end26.L(x,y)=-1;%标识孔洞区域27.end28.fori=(Ln^2*V_pore+1):Ln^229.ifrem(rand_l(i),Ln)==030.x=Ln;31.y=fix(rand_l(i)/Ln);32.else33.x=rem(rand_l(i),Ln);34.y=fix(rand_l(i)/Ln)+1;35.end36.rand_Q=randperm(Q);37.L(x,y)=rand_Q(1);%标识晶粒区域38.end39.40.temp_L=zeros(Ln+2);%标识边界区域标示为0，便于后续处理41.temp_L(2:Ln+1,2:Ln+1)=L;42.L=temp_L;43.Ln=Ln+2;%此时L边长Ln=Ln+244.45.s=[-1-146.-1047.-1148.0-149.0150.1-151.1052.11];%便于随即选取所选格点周围相邻的一个格点53.54.55.%开始CAS模拟56.forstep=1:step_num57.58.step%显示MCS进程59.60.rand_l=randperm(Ln^2);61.62.fori=1:Ln^2%随机选取格点63.64.%ifrem(i,1000)==0%显示选取格点进程65.%i66.%end67.68.ifrem(rand_l(i),Ln)==0%把rand_l(i)转换成实际坐标L(x,y)69.x=Ln;70.y=fix(rand_l(i)/Ln);71.else72.x=rem(rand_l(i),Ln);73.y=fix(rand_l(i)/Ln)+1;74.end75.76.ifL(x,y)~=0%如果格点不在边界区域则开始模拟77.78.%---------------------如果选取点为晶粒格点---------------------79.ifL(x,y)~=-1%如果选取点为晶粒格点80.81.ss=s+[xy;xy;xy;xy;xy;xy;xy;xy];%存储所选晶粒格点L(x,y)周围格点坐标82.83.diff_grain_count=0;%计数所选晶粒格点L(x,y)周围与之不同的晶粒格点数84.same_grain_count=0;%计数所选晶粒格点L(x,y)周围与之相同的晶粒格点数85.diff_ss=zeros(8,3);%存储所选晶粒格点L(x,y)周围与之不同的晶粒格点坐标及其取向度值86.forii=1:887.ifL(ss(ii,1),ss(ii,2))~=088.ifL(ss(ii,1),ss(ii,2))~=L(x,y)&L(ss(ii,1),ss(ii,2))~=-189.diff_grain_count=diff_grain_count+1;90.diff_ss(diff_grain_count,1)=ss(ii,1);91.diff_ss(diff_grain_count,2)=ss(ii,2);92.diff_ss(diff_grain_count,3)=L(ss(ii,1),ss(ii,2));93.end94.ifL(ss(ii,1),ss(ii,2))==L(x,y)95.same_grain_count=same_grain_count+1;96.end97.ifL(ss(ii,1),ss(ii,2))==-198.pore_count=pore_count+1;99.end100.end101.end102.total_grain_count=diff_grain_count+same_grain_count;103.104.ifdiff_grain_count~=0%如果所选晶粒格点L(x,y)周围有与之不同的晶粒格点105.106.BG_energy=J1*diff_grain_count;%BG_energy为格点所处晶界能107.108.diff_ss_1=diff_ss(1:diff_grain_count,3);%diff_ss_1存储被选格点周围取向度值与之不同的格点取向度值109.diff_ss_2=unique(diff_ss_1);%去除diff_ss_1中的重复元素，并存储到diff_ss_2110.rand_ll=randperm(length(diff_ss_2));111.temp_Q=diff_ss_2(rand_ll(1));112.change_BG_energy=J1*(total_grain_count-length(find(diff_ss_1==temp_Q)));113.114.ifchange_BG_energy=BG_energy115.L(x,y)=temp_Q;116.end117.118.ifchange_BG_energyBG_energy119.set_probability=rand();120.ifexp(-(change_BG_energy-BG_energy)/T)=set_probability121.L(x,y)=temp_Q;122.end123.end124.125.end%对应于如果所选晶粒格点L(x,y)周围有与之不同的晶粒格点126.127.end%对应于如果选取点为晶粒格点128.%---------------------如果选取点为晶粒格点---------------------129.130.end%对应于如果格点不在边界区域则开始模拟131.132.end%对应于随机选取格点133.134.%========================================================后处理过程========================================================135.136.%后处理1开始---------------每interval_save_jpg次MCS后存储图形矩阵---------------%137.ifrem(step,interval_save_jpg)==0138.figure1=figure('visible','off','PaperPosition',[3.0679.2814.8111.1],'PaperSize',[20.9829.68]);139.axes1=axes('Layer','top','YDir','reverse','Parent',figure1);140.axis(axes1,[0.5Ln-2+0.50.5Ln-2+0.5]);141.image1=image('CData',L(2:Ln-2+1,2:Ln-2+1),'CDataMapping','scaled','XData',[1Ln],'YData',[1Ln],'Parent',axes1);142.ifV_pore==0143.jpg_name=strcat(num2str(step),'_','Ln=',num2str(Ln-2),'_','Q=',num2str(Q),'_','T=',num2str(T),'_','J1=',num2str(J1),'.jpg');144.else145.jpg_name=strcat(num2str(step),'_','Ln=',num2str(Ln-2),'_','Q=',num2str(Q),'_','Vpore=',num2str(V_pore),'_','T=',num2str(T),'_','J1=',num2str(J1),'_','J2=',num2str(J2),'.jpg');146.end147.saveas(image1,jpg_name,'jpg');148.end149.%后处理1结束---------------每interval_save_jpg次MCS后存储图形矩阵---------------%150.151.152.%后处理2开始---------------每interval_stastics次MCS后统计平均晶粒尺寸，并存入stastics_data中---------------%153.ifrem(step,interval_stastics)==0154.155.stastics_data(step/interval_stastics,1)=step;%存储此次step156.157.temp_L=L(2:Ln-1,2:Ln-1);%去掉标示为0的边界158.Q_exist=unique(temp_L);%Q_exist存储L矩阵中仍存在的晶粒取向度Q159.ifQ_exist(1)==-1160.Q_exist=Q_exist(2:length(Q_exist));161.end162.Q_length=length(Q_exist);%Q_length为L矩阵中仍存在的晶粒取向度个数（不包括孔洞）163.164.forqq=1:Q_length%只统计具有在L矩阵中仍存在的取向度的晶粒的个数165.temp_L=L;166.temp_L(temp_L~=Q_exist(qq))=0;167.temp_L=bwlabel(temp_L,8);168.now_grainnum=max(max(temp_L));%返回目前取向度为Q_exist(qq)的晶粒数目169.stastics_data(step/interval_stastics,2)=stastics_data(step/interval_stastics,2)+now_grainnum;%累加晶粒个数170.end171.172.stastics_data(step/interval_stastics,3)=(Ln-2)^2*(1-V_pore)/stastics_data(step/interval_stastics,2);%统计晶粒平均面积173.174.stastics_data(step/interval_stastics,4)=sqrt(stastics_data(step/interval_stastics,3));%统计晶粒平均直径175.176.end177.%后处理2结束---------------每interval_stastics次MCS后统计平均晶粒尺寸，并存入stastics_data中---------------%178.179.180.%========================================================后处理过程========================================================18