case_study_3.zip_phase field_phase field model_phase-field

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % PHASE-FIELD FINITE-DIFFRENCE CODE FOR SINTERING % % (OPTIMIZED FOR MATLAB/OCTAVE) % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %== get intial wall time: time0=clock(); format long; %-- Simulation cell parameters: Nx = 100; Ny = 100; NxNy= Nx*Ny; dx = 0.5; dy = 0.5; %--- Time integration parameters: nstep = 5000; nprint = 50; dtime= 1.0e-4; %--- Material specific Parameters: npart = 2; % coefm = 5.0; coefk = 2.0; coefl = 5.0; % dvol = 0.040; dvap = 0.002; dsur = 16.0; dgrb = 1.6; %---- %prepare microstructure: %---- iflag =1; [npart,etas,con] = micro_sint_pre(Nx,Ny,npart,iflag); %-- initialize eta: for i=1:Nx for j=1:Ny eta(i,j) = 0.0; end end for istep =1:nstep %-- evolve concentration: iflag =1; for i=1:Nx for j=1:Ny jp=j+1; jm=j-1; ip=i+1; im=i-1; jp=j+1; jm=j-1; ip=i+1; im=i-1; if(im == 0) im=Nx; end if(ip == (Nx+1)) ip=1; end if(jm == 0) jm = Ny; end if(jp == (Ny+1)) jp=1; end hne=con(ip,j); hnw=con(im,j); hns=con(i,jm); hnn=con(i,jp); hnc=con(i,j); lap_con(i,j) =(hnw + hne + hns + hnn -4.0*hnc)/(dx*dy); %--- derivative of free energy: [dfdcon,dfdeta]=free_energ_sint_v1(i,j,con,eta,etas,npart,iflag); dummy(i,j) = dfdcon - 0.5*coefm*lap_con(i,j); end %for j end %for i %-- for i=1:Nx for j=1:Ny jp=j+1; jm=j-1; ip=i+1; im=i-1; jp=j+1; jm=j-1; ip=i+1; im=i-1; if(im == 0) im=Nx; end if(ip == (Nx+1)) ip=1; end if(jm == 0) jm = Ny; end if(jp == (Ny+1)) jp=1; end hne=dummy(ip,j); hnw=dummy(im,j); hns=dummy(i,jm); hnn=dummy(i,jp); hnc=dummy(i,j); lap_dummy(i,j) =(hnw + hne + hns + hnn -4.0*hnc)/(dx*dy); %-- Mobility: phi = con(i,j)^3 *(10.0-15.0*con(i,j) + 6.0*con(i,j)^2); sum =0.0; for ipart =1:npart for jpart =1:npart if(ipart ~= jpart) sum = sum + etas(i,j,ipart)*etas(i,j,jpart); end end end mobil = dvol*phi + dvap*(1.0-phi) + dsur*con(i,j)*(1.0-con(i,j)) + 2.0*dgrb*sum; %-- time integration: con(i,j) = con(i,j) + dtime*mobil*lap_dummy(i,j); %-- for small deviations: if(con(i,j) >= 0.9999); con(i,j)= 0.9999; end if(con(i,j) < 0.00001); con(i,j) = 0.00001; end end %j end %i %-- evolve etas: iflag =2; for ipart = 1:npart for i=1:Nx for j=1:Ny eta(i,j) = etas(i,j,ipart); end end for i=1:Nx for j=1:Ny jp=j+1; jm=j-1; ip=i+1; im=i-1; jp=j+1; jm=j-1; ip=i+1; im=i-1; if(im == 0) im=Nx; end if(ip == (Nx+1)) ip=1; end if(jm == 0) jm = Ny; end if(jp == (Ny+1)) jp=1; end hne=eta(ip,j); hnw=eta(im,j); hns=eta(i,jm); hnn=eta(i,jp); hnc=eta(i,j); lap_eta(i,j) =(hnw + hne + hns + hnn -4.0*hnc)/(dx*dy); %--- derivative of free energy: [dfdcon,dfdeta]=free_energ_sint_v1(i,j,con,eta,etas,npart,iflag); %-- Time integration eta(i,j) = eta(i,j) - dtime * coefl*(dfdeta - 0.5 *coefk*lap_eta(i,j)); %-- for small deviations: if(eta(i,j) >= 0.9999); eta(i,j) = 0.9999; end if(eta(i,j) < 0.0001); eta(i,j) = 0.0001; end end %j end %i %-- for i=1:Nx for j=1:Ny etas(i,j,ipart) = eta (i,j); end end end %ipart %---- print results if((mod(istep,nprint) == 0) || (istep == 1) ) fprintf('done step: %5d\n',istep); %fname1 =sprintf('time_%d.out',istep); %out1 = fopen(fname1,'w'); %for i=1:Nx %for j=1:Ny %sum = 0.0; %for ipart=1:npart %sum = sum +etas(i,j,ipart)^2; %end %fprintf(out1,'%5d %5d %14.6e %14.6e\n',i,j,con(i,j),sum); %end %end %fclose(out1); %--- write vtk file: phi2 =zeros(Nx,Ny); for i=1:Nx for j=1:Ny for ipart=1:npart phi2(i,j) = phi2(i,j) +etas(i,j,ipart)^2; end end end write_vtk_grid_values(Nx,Ny,dx,dy,istep,con,phi2); end %end if end %istep %--- calculate compute time: compute_time = etime(clock(),time0); fprintf('Compute Time: %10d\n',compute_time);