program pitting_corrosion implicit none character(len=20) filename double precision, parameter:: total_t = 4. ! time for corrosion [s] double precision, parameter:: dx = 0.1e-5 ! spacing of computational grids [m] double precision, parameter:: dy = 0.1e-5 double precision, parameter:: dz = 0.1e-5 double precision, parameter:: dxdx = dx*dx double precision, parameter:: dydy = dy*dy double precision, parameter:: dzdz = dz*dz double precision, parameter:: dxdy = dx*dy double precision, parameter:: dxdz = dx*dz double precision, parameter:: dydz = dy*dz double precision, parameter:: dxdydz = dx*dy*dz double precision, parameter:: dt = 1.0e-7 ! time increment [s] double precision, parameter:: square = (0.2e-3)**2. ! size of collector [m2] double precision, parameter:: RRR = 5.*dx ! radius of hole [m] integer, parameter:: nend = nint(total_t/dt) ! total number of time steps integer, parameter:: num_grain = 1+1 ! number of grains (=1) + electrolyte integer, parameter:: file_num = 100 ! the number of vtk file for output integer, parameter:: nout = nend/file_num ! data output interval integer, parameter:: ini_loop = 2000000 ! loop for initial diffusion integer, parameter:: nd = 2**5 ! number of finite difference grids along x direction integer, parameter:: ndm = nd-1 integer, parameter:: nd2 = nd/2 integer, parameter:: nd3 = nd*3./4. integer, parameter:: nd3m = nd3-1 integer, parameter:: mlax = nd integer, parameter:: mlay = nd integer, parameter:: mlaz = nd integer, parameter:: loop_max=10000 ! maximum number of iteration for Gauss-saidel method integer, dimension(6,0:nd,0:nd,0:nd):: no ! grain number double precision, parameter:: RR = 8.31446 ! gas constant [J/K/mol] double precision, parameter:: T = 273.0 + 60.0 ! temperature [K] double precision, parameter:: RT = RR*T double precision, parameter:: F = 96485.3 ! Faradays constant [C/mol] double precision, parameter:: pi = 3.141592 double precision, dimension(mlax,mlay,mlaz):: p,pp ! phase-field variable at time t and t+dt double precision, dimension(mlax,mlay,mlaz):: c_H,cc_H ! H ion concentration double precision, dimension(mlax,mlay,mlaz):: c_Fe,cc_Fe ! Fe ion concentration double precision, dimension(mlax,mlay,mlaz):: c_OH,cc_OH ! OH ion concentration double precision, dimension(mlax,mlay,mlaz):: c_Cl,cc_Cl ! Cl ion concentration double precision, dimension(mlax,mlay,mlaz):: v,vv,pre_v ! potential double precision, dimension(mlax,mlay,mlaz):: D_Fe,D_H ! diffusion coefficient double precision, dimension(mlax,mlay,mlaz):: D_OH,D_Cl ! diffusion coefficient double precision, dimension(mlax,mlay,mlaz):: sgm_eff ! conductivity double precision, dimension(mlax,mlay,mlaz):: dv double precision, dimension(mlax,mlay,mlaz):: I_Fe,I_H ! current double precision, dimension(mlax,mlay,mlaz):: rou ! charge density double precision, dimension(mlax,mlay,mlaz):: Jd,Jm ! diffusion term and intercaltion flux term double precision, dimension(mlax,mlay,mlaz):: dpi ! for Allen-Cahn double precision, dimension(mlax,mlay,mlaz):: acc2 double precision, parameter:: mob_i = 5.e-7 ! phase-field mobility for interfacial part double precision, parameter:: mob_s = 5.e-5 double precision, parameter:: del = 3.*dx double precision, parameter:: aaa = sqrt(3.*del*0.5/2.2)! gradient energy coefficient double precision, parameter:: www = 3.*2.2/del ! potential height double precision, parameter:: zzz_H = 1.0 ! valence of H ion double precision, parameter:: zzz_Fe = 2.0 ! valence of Fe ion double precision, parameter:: zzz_OH = -1.0 ! valence of OH ion double precision, parameter:: zzz_Cl = -1.0 ! valence of Cl ion double precision, parameter:: zzz = 2.0 ! valence of the ion double precision, parameter:: alph = 0.5 ! permeability coefficient of Butler-Volmer kinetics double precision, parameter:: beta = 0.5 ! permeability coefficient of Butler-Volmer kinetics double precision, parameter:: azF = alph*zzz*F double precision, parameter:: bzF = beta*zzz*F double precision, parameter:: Ds_Fe = 1.52e-9 ! Fe ion diffusion coefficient in solution [m2/s] double precision, parameter:: Ds_H = 15.e-9 ! H ion diffusion coefficient in solution [m2/s] double precision, parameter:: Ds_OH = 9.e-9 ! OH ion diffusion coefficient in solution [m2/s] double precision, parameter:: Ds_Cl = 4.e-9 ! Cl ion diffusion coefficient in solution [m2/s] double precision, parameter:: De = 0. ! Ions diffusion coefficient in electrode [m2/s] double precision, parameter:: Vm = 7.09e-6 ! molar volume [m3/mol] double precision, parameter:: cs = 1./Vm ! site density [mol/m3] double precision, parameter:: c0 = 1.e+3 ! Fe2+ bulk concentration double precision, parameter:: pH = 2.5 ! pH double precision, parameter:: sgm_e = 1.e+7 ! conductivity of electrode [S/m] double precision, parameter:: sgm_s = 100. ! conductivity of electrolyte [S/m] double precision, parameter:: Eeq_Fe = -0.54 ! equilibrium potential [V] double precision, parameter:: EEE = -0.42 ! double precision, parameter:: conv_v = 1.e-10 ! convergent calculation double precision, parameter:: SOR = 1.2 ! parameter for SOR method double precision, parameter:: KKK1 = 1.e+6 ! equilibrium constant double precision, parameter:: KKK3 = 1.e+6 ! equilibrium constant integer :: nnn integer:: l,m,n,ll,mm,nn integer:: lp,lm,mp,np,nm,loop,num,number double precision:: r,r0,ppp1,ppp2,ppp3,ppp4,activ1,activ2,Ia_Fe,phi0,radi double precision:: ips,abc,cFe_total,cH_total,xxx,yyy,pall,mob_r double precision:: multi,dpi_all,x0,y0,Ia,ic,acc double precision:: Ecorr,Icorr,max_v,eta,vlct,prev_p,i_star double precision:: x1,x2,x3,x4,x5,x6,xx double precision:: y1,y2,y3,y4,y5,y6,yy double precision:: z1,z2,z3,z4,z5,z6,zz double precision:: s1,s2,s3,s4,s5,s6,ss double precision:: t1,t2,t3,t4,t5,t6,tt double precision:: u1,u2,u3,u4,u5,u6,uu double precision:: v11,v12,v13,v14,v15,v16,v17 double precision:: v21,v22,v23,v24 double precision:: ini_cH,ini_cOH,ini_cCl,ini_c,DDD integer:: ii,jj,i,j,ie,iw,jn,js integer:: nstep,iout,kk,time integer:: h,z,zf,zb double precision, dimension(0:nd,0:nd,0:nd):: n_x, n_y, n_z !***********************************initial setting**************************************** p = 0. c_Fe = 0. c_H = 0. c_OH = 0. c_Cl = 0. v = 0. dv = 0. pre_v = 0. I_Fe = 0. I_H = 0. rou = 0. Jd = 0. Jm = 0. dpi = 0. vv = EEE ! initial current density Ia = 0. Ic = 0. if(pH.eq.2.5 .and. EEE.lt.-0.41 .and. EEE.gt.-0.43)then Ia = 110. / 100. mob_r = 2.e+6 endif if(Ia.eq.0)then write(*,*) 'Current denisty I is zero' stop else write(*,*) '=============================================================' write(*,*) 'Simulate corr. for',real(nend*dt),'sec.' write(*,*) 'pH =',real(pH),'and EEE =',real(EEE) write(*,*) '---> Ia is',real(Ia),'A/m2' endif ! initial setting of phase-field variable r0 = nd3*dx do n=1,mlaz do m=1,mlay do l=1,mlax if(l.le.nd3) p(l,m,n) = 1. if(l.gt.nd3) p(l,m,n) = 0. r = sqrt(dble((m-nd2)*(m-nd2)*dydy+(n-nd2)*(n-nd2)*dzdz)) if(r.le.RRR) then r = real(l*dx) -r0 p(l,m,n) = 0.5*(1.-tanh(sqrt(2.*www)/(2.*aaa)*r)) endif enddo enddo enddo write(*,*) '=============================================================' ! initial concentration-field call inidiff(nd,p,c_Cl,c_OH,c_H,cc_Cl,cc_OH,cc_H,D_Cl,D_OH,D_H,ini_loop,dxdx,dydy,dzdz,dt,pH,Ds_Fe,Ds_H,Ds_OH,Ds_Cl,De) cc_H = c_H cc_Fe = c_Fe cc_OH = c_OH cc_Cl = c_Cl pp = p !---- setting of parameter "multi" eta=EEE-Eeq_Fe dpi_all=0.0 number=0 do n=1,mlaz do m=1,mlay do l=1,mlax lp=l+1 if(l.eq.mlax) lp=mlax lm=l-1 if(l.eq.1) lm=1 mp=m+1 if(m.eq.mlay) mp=mlay mm=m-1 if(m.eq.1) mm=1 np=n+1 if(n.eq.mlaz) np=mlaz nm=n-1 if(n.eq.1) nm=1 ppp1 = 2.*www*p(l,m,n)*(p(l,m,n)-1.)*(2.*p(l,m,n)-1.) !! double-well potential term ppp2 = aaa*((p(lp,m,n)-2.*p(l,m,n)+p(lm,m,n))/dxdx & +(p(l,mp,n)-2.*p(l,m,n)+p(l,mm,n))/dydy & +(p(l,m,np)-2.*p(l,m,n)+p(l,m,nm))/dzdz) !! gradient energy term activ1 = KKK1*c_OH(l,m,n)/c0 activ2 = c_Fe(l,m,n)/c0/KKK3/c_H(l,m,n)/c0 ppp3 = 6.*p(l,m,n)*(1.-p(l,m,n))*(activ1*exp((azF*eta)/RT)-activ2*exp(-(bzF*eta)/RT)) !! potential term if(c_H(l,m,n).le.1.e-10) ppp3=0 dpi(l,m,n)= -mob_i*( ppp1 - ppp2 ) - mob_r*ppp3 !! non-linear phase-field model dpi_all=dpi_all+dpi(l,m,n) enddo enddo enddo multi=Ia/F/zzz_Fe*square/(dabs(dpi_all)*cs*dxdydz) if(multi.le.0)then write(*,*) 'multi is error!!','multi is',multi stop endif !*********************** output of initial data *********************** write(filename,'(a,a)') 'ini_data','.dat' open(58,file=filename,position='append') write(58,*) 'time length is',real(nend*dt),'sec.' write(58,*) 'Total step is',real(nend/1e+6),'x 10^6' write(58,*) 'cs=',real(cs),'mol/m3' write(58,*) 'square=',real(square),'mm2' write(58,*) 'EEE=',real(EEE),'V' write(58,*) 'mob_i=',real(mob_i*multi) write(58,*) 'mob_r=',real(mob_r*multi) close(58) iout = 0 nstep = 0 acc2 = 0. call fpc(iout,mlax,mlay,mlaz,p,pp,dpi,multi,dt,acc2) call feionc(iout,mlax,mlay,mlaz,c_Fe,cc_Fe) call Hionc(iout,mlax,mlay,mlaz,c_H,cc_H,pp) call OHionc(iout,mlax,mlay,mlaz,c_OH,cc_OH) call Clionc(iout,mlax,mlay,mlaz,c_Cl,cc_Cl) call vvv(iout,mlax,mlay,mlaz,vv) !====================================== MAIN LOOP START ========================================== write(*,*) '----------------------Main loop start----------------------' write(*,'((a),1i4)') ' ***file***',iout do nstep=1,nend if(mod(nstep,nout/100).eq.0 .and. mod(nstep,nout).eq.0)then write(*,'(1i10,(a))') nstep,' steps' write(*,'((a),1i4)') ' ***file***',iout+1 endif if(mod(nstep,nout/100).eq.0 .and. mod(nstep,nout).ne.0)then write(*,'(1i10,(a))') nstep,' steps' endif !-------- solving phase-field equation -------- eta=EEE-Eeq_Fe dpi_all=0.0 number=0 do n=1,mlaz do m=1,mlay do l=1,mlax r = sqrt(dble((m-nd2)*(m-nd2)*dydy+(n-nd2)*(n-nd2)*dzdz)) lp=l+1 if(l.eq.mlax) lp=mlax if(r.gt.RRR .and. l.eq.nd3) lp=nd3 lm=l-1 if(l.eq.1) lm=1 if(r.gt.RRR .and. l.eq.nd3+1) lm=nd3+1 mp=m+1 if(m.eq.mlay) mp=mlay mm=m-1 if(m.eq.1) mm=1 np=n+1 if(n.eq.mlaz) np=mlaz nm=n-1 if(n.eq.1) nm=1 ppp1 = 2.*www*p(l,m,n)*(p(l,m,n)-1.)*(2.*p(l,m,n)-1.) ! double-well potential term ppp2 = aaa*((p(lp,m,n)-2.*p(l,m,n)+p(lm,m,n))/dxdx & +(p(l,mp,n)-2.*p(l,m,n)+p(l,mm,n))/dydy & +(p(l,m,np)-2.*p(l,m,n)+p(l,m,nm))/dzdz) ! gradient energy term activ1 = KKK1*c_OH(l,m,n)/c0 if(c_H(l,m,n).eq.0) c_H(l,m,n)=1.e-20 activ2 = c_Fe(l,m,n)/c0/KKK3/c_H(l,m,n)/c0 ppp3 = 6.*p(l,m,n)*(1.-p(l,m,n))*(activ1*exp((azF*eta)/RT)-activ2*exp(-(bzF*eta)/RT)) ! potential term dpi(l,m,n)= -mob_i*( ppp1 - ppp2 ) - mob_r*ppp3 !! non-linear phase-field model if(l.gt.nd3m .and. r.gt.RRR) dpi(l,m,n)=0. dpi_all=dpi_all+dpi(l,m,n) enddo enddo enddo acc = 0. acc2 = 0. do n=1,mlaz do m=1,mlay do l=1,mlax acc = 3.-LOG10(c_H(l,m,n)) ! value of pH if(acc.ge.2.5)then acc = 2.5 acc2(l,m,n) = 1. else if(acc.gt.-2.4 .and. acc.le.-1.0)then acc2(l,m,n) = (-60792*acc*acc*acc*acc -238835*acc*acc*acc -227116*acc*acc -92012*acc -4712.2) / (Ia*100.) else if(acc.gt.-1.0 .and. acc.le.0.0)then acc2(l,m,n) = (-34543*acc*acc*acc -6705.3*acc*acc -7036.4*acc +3949.4) / (Ia*100.) else if(acc.gt.0.0 .and. acc.le.1.0)then acc2(l,m,n) = (1281.9*acc*acc -3982.9*acc +3974.9) / (Ia*100.) else if(acc.gt.1.0 .and. acc.le.2.5)then acc2(l,m,n) = (471.25*acc*acc -2384.3*acc +3145.1) / (Ia*100.) endif if(acc2(l,m,n).lt.0)then write(*,*) 'acc2 error' stop else if(acc2(l,m,n).lt.1)then acc2(l,m,n) = 1. endif pp(l,m,n) = p(l,m,n) + multi*acc2(l,m,n)*dpi(l,m,n)*dt ! update of phase-field variable ! pre-calculation for poisson equation abc = 3.*pp(l,m,n)*pp(l,m,n)-2*pp(l,m,n)*pp(l,m,n)*pp(l,m,n) sgm_eff(l,m,n)=sgm_e*abc+sgm_s*(1.-abc) rou(l,m,n)=-F*zzz_Fe*cs*multi*acc2(l,m,n)*dpi(l,m,n) ! pre-calculation for cahn-hilliard equation D_Fe(l,m,n) = De * abc + Ds_Fe * (1.-abc) D_H (l,m,n) = De * abc + Ds_H * (1.-abc) D_OH(l,m,n) = De * abc + Ds_OH * (1.-abc) D_Cl(l,m,n) = De * abc + Ds_Cl * (1.-abc) enddo enddo enddo ! output of result if(mod(nstep,nout).eq.0)then iout=iout+1 call fpc(iout,mlax,mlay,mlaz,p,pp,dpi,multi,dt,acc2) endif ! solving Poisson equation do loop=1,1000000 max_v=0.0 do n=1,mlaz do m=1,mlay do l=1,mlax lp=l+1 if(l.eq.mlax) lp=mlax lm=l-1 if(l.eq.1) lm=1 mp=m+1 if(m.eq.mlay) mp=mlay mm=m-1 if(m.eq.1) mm=1 np=n+1 if(n.eq.mlaz) np=mlaz nm=n-1 if(n.eq.1) nm=1 if(p(l,m,n).ge.0.99) vv(l,m,n)=EEE pre_v(l,m,n)=vv(l,m,n) if(p(l,m,n).ge.0.99) goto 3939 if(l.ge.nd3m .and. p(lm,m,n).ge.0.9)then vv(l,m,n)=vv(l,m,n)+SOR*((vv(lp,m,n)+vv(l,m,n)+vv(l,mp,n)+vv(l,mm,n) & +vv(l,m,np)+vv(l,m,nm)+rou(l,m,n)*dxdx/sgm_s)/6.-vv(l,m,n)) else vv(l,m,n)=vv(l,m,n)+SOR*((vv(lp,m,n)+vv(lm,m,n)+vv(l,mp,n)+vv(l,mm,n) & +vv(l,m,np)+vv(l,m,nm)+rou(l,m,n)*dxdx/sgm_s)/6.-vv(l,m,n)) endif 3939 continue dv(l,m,n)=(dabs(pre_v(l,m,n)-vv(l,m,n))) if(max_v.lt.dv(l,m,n))then max_v=dv(l,m,n) endif enddo enddo enddo if(max_v.lt.conv_v) goto 1582 enddo write(*,*) 'Loop is more than 1000000!!!' stop 1582 continue ! output of simulation result if(mod(nstep,nout).eq.0)then call vvv(iout,mlax,mlay,mlaz,vv) endif ! solving Cahn-Hilliard equation !============= Fe2+ ion ============= Ia_Fe=0. do n=1,mlaz do m=1,mlay do l=1,mlax lp=l+1 if(l.eq.mlax) lp=mlax lm=l-1 if(l.eq.1) lm=1 mp=m+1 if(m.eq.mlay) mp=mlay mm=m-1 if(m.eq.1) mm=1 np=n+1 if(n.eq.mlaz) np=mlaz nm=n-1 if(n.eq.1) nm=1 x1= ((D_Fe(lp,m,n)+D_Fe(l ,m,n))*(c_Fe(lp,m,n)-c_Fe(l ,m,n)))/(2.*dxdx) x2=0.-((D_Fe(l ,m,n)+D_Fe(lm,m,n))*(c_Fe(l ,m,n)-c_Fe(lm,m,n)))/(2.*dxdx) xx=x1+x2 y1= ((D_Fe(l,mp,n)+D_Fe(l,m ,n))*(c_Fe(l,mp,n)-c_Fe(l,m ,n)))/(2.*dydy) y2=0.-((D_Fe(l,m ,n)+D_Fe(l,mm,n))*(c_Fe(l,m ,n)-c_Fe(l,mm,n)))/(2.*dydy) yy=y1+y2 z1= ((D_Fe(l,m,np)+D_Fe(l,m,n ))*(c_Fe(l,m,np)-c_Fe(l,m,n )))/(2.*dzdz) z2=0.-((D_Fe(l,m,n )+D_Fe(l,m,nm))*(c_Fe(l,m,n )-c_Fe(l,m,nm)))/(2.*dzdz) zz=z1+z2 Jd(l,m,n)=xx+yy+zz s1=(D_Fe(lp,m,n)+D_Fe(l ,m,n))*(c_Fe(lp,m,n)+c_Fe(l ,m,n))*(vv(lp,m ,n)-vv(l ,m,n))/(4.*dxdx) s2=(D_Fe(l ,m,n)+D_Fe(lm,m,n))*(c_Fe(l ,m,n)+c_Fe(lm,m,n))*(vv(l ,m ,n)-vv(lm,m,n))/(4.*dxdx) ss=s1-s2 t1=(D_Fe(l,mp,n)+D_Fe(l,m ,n))*(c_Fe(l,mp,n)+c_Fe(l,m ,n))*(vv(l,mp,n)-vv(l,m ,n))/(4.*dydy) t2=(D_Fe(l,m ,n)+D_Fe(l,mm,n))*(c_Fe(l,m ,n)+c_Fe(l,mm,n))*(vv(l,m ,n)-vv(l,mm,n))/(4.*dydy) tt=t1-t2 u1=(D_Fe(l,m,np)+D_Fe(l,m,n ))*(c_Fe(l,m,np)+c_Fe(l,m,n ))*(vv(l,m,np)-vv(l,m,n ))/(4.*dzdz) u2=(D_Fe(l,m,n )+D_Fe(l,m,nm))*(c_Fe(l,m,n )+c_Fe(l,m,nm))*(vv(l,m,n )-vv(l,m,nm))/(4.*dzdz) uu=u1-u2 Jm(l,m,n)=(ss+tt+uu)*zzz_Fe*F/RT I_Fe(l,m,n)=0. I_Fe(l,m,n)=dabs(cs*multi*acc2(l,m,n)*dpi(l,m,n)) cc_Fe(l,m,n) = c_Fe(l,m,n) + (Jd(l,m,n) + Jm(l,m,n) + I_Fe(l,m,n))*dt Ia_Fe=Ia_Fe+I_Fe(l,m,n) !============= H+ ion============= x1= ((D_H(lp,m,n)+D_H(l ,m,n))*(c_H(lp,m,n)-c_H(l ,m,n)))/(2.*dxdx) x2=0.-((D_H(l ,m,n)+D_H(lm,m,n))*(c_H(l ,m,n)-c_H(lm,m,n)))/(2.*dxdx) xx=x1+x2 y1= ((D_H(l,mp,n)+D_H(l,m ,n))*(c_H(l,mp,n)-c_H(l,m ,n)))/(2.*dydy) y2=0.-((D_H(l,m ,n)+D_H(l,mm,n))*(c_H(l,m ,n)-c_H(l,mm,n)))/(2.*dydy) yy=y1+y2 z1= ((D_H(l,m,np)+D_H(l,m,n ))*(c_H(l,m,np)-c_H(l,m,n )))/(2.*dzdz) z2=0.-((D_H(l,m,n )+D_H(l,m,nm))*(c_H(l,m,n )-c_H(l,m,nm)))/(2.*dzdz) zz=z1+z2 Jd(l,m,n)=xx+yy+zz s1=(D_H(lp,m,n)+D_H(l ,m,n))*(c_H(lp,m,n)+c_H(l ,m,n))*(vv(lp,m ,n)-vv(l ,m,n))/(4.*dxdx) s2=(D_H(l ,m,n)+D_H(lm,m,n))*(c_H(l ,m,n)+c_H(lm,m,n))*(vv(l ,m ,n)-vv(lm,m,n))/(4.*dxdx) ss=s1-s2 t1=(D_H(l,mp,n)+D_H(l,m ,n))*(c_H(l,mp,n)+c_H(l,m ,n))*(vv(l,mp,n)-vv(l,m ,n))/(4.*dydy) t2=(D_H(l,m ,n)+D_H(l,mm,n))*(c_H(l,m ,n)+c_H(l,mm,n))*(vv(l,m ,n)-vv(l,mm,n))/(4.*dydy) tt=t1-t2 u1=(D_H(l,m,np)+D_H(l,m,n ))*(c_H(l,m,np)+c_H(l,m,n ))*(vv(l,m,np)-vv(l,m,n ))/(4.*dzdz) u2=(D_H(l,m,n )+D_H(l,m,nm))*(c_H(l,m,n )+c_H(l,m,nm))*(vv(l,m,n )-vv(l,m,nm))/(4.*dzdz) uu=u1-u2 Jm(l,m,n)=(ss+tt+uu)*zzz_H*F/RT I_H(l,m,n)=0. I_H(l,m,n) = I_Fe(l,m,n)*2. cc_H(l,m,n) = c_H(l,m,n) + (Jd(l,m,n) + Jm(l,m,n) + I_H(l,m,n))*dt !=============== OH- ion =============== x1= ((D_OH(lp,m,n)+D_OH(l ,m,n))*(c_OH(lp,m,n)-c_OH(l ,m,n)))/(2.*dxdx) x2=0.-((D_OH(l ,m,n)+D_OH(lm,m,n))*(c_OH(l ,m,n)-c_OH(lm,m,n)))/(2.*dxdx) xx=x1+x2 y1= ((D_OH(l,mp,n)+D_OH(l,m ,n))*(c_OH(l,mp,n)-c_OH(l,m ,n)))/(2.*dydy) y2=0.-((D_OH(l,m ,n)+D_OH(l,mm,n))*(c_OH(l,m ,n)-c_OH(l,mm,n)))/(2.*dydy) yy=y1+y2 z1= ((D_OH(l,m,np)+D_OH(l,m,n ))*(c_OH(l,m,np)-c_OH(l,m,n )))/(2.*dzdz) z2=0.-((D_OH(l,m,n )+D_OH(l,m,nm))*(c_OH(l,m,n )-c_OH(l,m,nm)))/(2.*dzdz) zz=z1+z2 Jd(l,m,n)=xx+yy+zz s1=(D_OH(lp,m,n)+D_OH(l ,m,n))*(c_OH(lp,m,n)+c_OH(l ,m,n))*(vv(lp,m,n)-vv(l ,m,n))/(4.*dxdx) s2=(D_OH(l ,m,n)+D_OH(lm,m,n))*(c_OH(l ,m,n)+c_OH(lm,m,n))*(vv(l ,m,n)-vv(lm,m,n))/(4.*dxdx) ss=s1-s2 t1=(D_OH(l,mp,n)+D_OH(l,m ,n))*(c_OH(l,mp,n)+c_OH(l,m ,n))*(vv(l,mp,n)-vv(l,m ,n))/(4.*dydy) t2=(D_OH(l,m ,n)+D_OH(l,mm,n))*(c_OH(l,m ,n)+c_OH(l,mm,n))*(vv(l,m ,n)-vv(l,mm,n))/(4.*dydy) tt=t1-t2 u1=(D_OH(l,m,np)+D_OH(l,m,n ))*(c_OH(l,m,np)+c_OH(l,m,n ))*(vv(l,m,np)-vv(l,m,n ))/(4.*dzdz) u2=(D_OH(l,m,n )+D_OH(l,m,nm))*(c_OH(l,m,n )+c_OH(l,m,nm))*(vv(l,m,n )-vv(l,m,nm))/(4.*dzdz) uu=u1-u2 Jm(l,m,n)=(ss+tt+uu)*zzz_OH*F/RT cc_OH(l,m,n) = c_OH(l,m,n) + (Jd(l,m,n) + Jm(l,m,n))*dt if(l.eq.mlax) cc_OH(l,m,n) = c_OH(l,m,n) !=============== Cl- ion =============== x1= ((D_Cl(lp,m,n)+D_Cl(l ,m,n))*(c_Cl(lp,m,n)-c_Cl(l ,m,n)))/(2.*dxdx) x2=0.-((D_Cl(l ,m,n)+D_Cl(lm,m,n))*(c_Cl(l ,m,n)-c_Cl(lm,m,n)))/(2.*dxdx) xx=x1+x2 y1= ((D_Cl(l,mp,n)+D_Cl(l,m ,n))*(c_Cl(l,mp,n)-c_Cl(l,m ,n)))/(2.*dydy) y2=0.-((D_Cl(l,m ,n)+D_Cl(l,mm,n))*(c_Cl(l,m ,n)-c_Cl(l,mm,n)))/(2.*dydy) yy=y1+y2 z1= ((D_Cl(l,m,np)+D_Cl(l,m,n ))*(c_Cl(l,m,np)-c_Cl(l,m,n )))/(2.*dzdz) z2=0.-((D_Cl(l,m,n )+D_Cl(l,m,nm))*(c_Cl(l,m,n )-c_Cl(l,m,nm)))/(2.*dzdz) zz=z1+z2 Jd(l,m,n)=xx+yy+zz s1=(D_Cl(lp,m,n)+D_Cl(l ,m,n))*(c_Cl(lp,m,n)+c_Cl(l ,m,n))*(vv(lp,m,n)-vv(l ,m,n))/(4.*dxdx) s2=(D_Cl(l ,m,n)+D_Cl(lm,m,n))*(c_Cl(l ,m,n)+c_Cl(lm,m,n))*(vv(l ,m,n)-vv(lm,m,n))/(4.*dxdx) ss=s1-s2 t1=(D_Cl(l,mp,n)+D_Cl(l,m ,n))*(c_Cl(l,mp,n)+c_Cl(l,m ,n))*(vv(l,mp,n)-vv(l,m ,n))/(4.*dydy) t2=(D_Cl(l,m ,n)+D_Cl(l,mm,n))*(c_Cl(l,m ,n)+c_Cl(l,mm,n))*(vv(l,m ,n)-vv(l,mm,n))/(4.*dydy) tt=t1-t2 u1=(D_Cl(l,m,np)+D_Cl(l,m,n ))*(c_Cl(l,m,np)+c_Cl(l,m,n ))*(vv(l,m,np)-vv(l,m,n ))/(4.*dzdz) u2=(D_Cl(l,m,n )+D_Cl(l,m,nm))*(c_Cl(l,m,n )+c_Cl(l,m,nm))*(vv(l,m,n )-vv(l,m,nm))/(4.*dzdz) uu=u1-u2 Jm(l,m,n)=(ss+tt+uu)*zzz_Cl*F/RT cc_Cl(l,m,n) = c_Cl(l,m,n) + (Jd(l,m,n) + Jm(l,m,n))*dt if(l.eq.mlax) cc_Cl(l,m,n) = c_Cl(l,m,n) enddo enddo enddo !-------- output of results ------- if(mod(nstep,nout).eq.0)then call feionc(iout,mlax,mlay,mlaz,c_Fe,cc_Fe) call Hionc(iout,mlax,mlay,mlaz,c_H,cc_H,pp) call OHionc(iout,mlax,mlay,mlaz,c_OH,cc_OH) call Clionc(iout,mlax,mlay,mlaz,c_Cl,cc_Cl) endif time = nstep * dt ! updating data do n=1,mlaz do m=1,mlay do l=1,mlax p(l,m,n) = pp(l,m,n) c_Fe(l,m,n)= cc_Fe(l,m,n) c_H(l,m,n) = cc_H(l,m,n) c_OH(l,m,n)= cc_OH(l,m,n) c_Cl(l,m,n)= cc_Cl(l,m,n) enddo enddo enddo enddo ! for main loop end program pitting_corrosion !======================================= MAIN LOOP END ============================================ !=============== SUBROUTINES ================ subroutine fpc(iout,mlax,mlay,mlaz,p,pp,dpi,multi,dt,acc2) implicit none character(len=12) filename double precision, dimension(mlax,mlay,mlaz):: p,pp,dpi,acc2 integer:: iout,mlax,mlay,mlaz,l,m,n double precision :: multi,dt write(filename,'(a,i3.3,a)') 'p',iout,'.vtk' open(1,file=filename) write(1,'(a)') '# vtk DataFile Version 3.0' write(1,'(a)') 'output.vtk' write(1,'(a)') 'ASCII' write(1,'(a)') 'DATASET STRUCTURED_POINTS' write(1,'(a,3i5)') 'DIMENSIONS', mlax,mlay,mlaz write(1,'(a,3f4.1)') 'ORIGIN ', 0.0, 0.0, 0.0 write(1,'(a,3i2)') 'ASPECT_RATIO', 1, 1, 1 write(1,'(a,1i11)') 'POINT_DATA', mlax*mlay*mlaz write(1,'(a)') 'SCALARS volume_scalars double' write(1,'(a)') 'LOOKUP_TABLE default' do n=1,mlaz do m=1,mlay do l=1,mlax write(1,*) pp(l,m,n) enddo enddo enddo close(1) return end !======== Fe2+ ion =========== subroutine feionc(iout,mlax,mlay,mlaz,c_Fe,cc_Fe) implicit none character(len=12) filename double precision, dimension(mlax,mlay,mlaz):: c_Fe,cc_Fe integer:: iout,mlax,mlay,mlaz,l,m,n write(filename,'(a,i3.3,a)') 'c_Fe',iout,'.vtk' open(2,file=filename) write(2,'(a)') '# vtk DataFile Version 3.0' write(2,'(a)') 'output.vtk' write(2,'(a)') 'ASCII' write(2,'(a)') 'DATASET STRUCTURED_POINTS' write(2,'(a,3i5)') 'DIMENSIONS', mlax,mlay,mlaz write(2,'(a,3f4.1)') 'ORIGIN ', 0.0, 0.0, 0.0 write(2,'(a,3i2)') 'ASPECT_RATIO', 1, 1, 1 write(2,'(a,1i11)') 'POINT_DATA', mlax*mlay*mlaz write(2,'(a)') 'SCALARS volume_scalars double' write(2,'(a)') 'LOOKUP_TABLE default' do n=1,mlaz do m=1,mlay do l=1,mlax write(2,*) cc_Fe(l,m,n) enddo enddo enddo close(2) return end !======== H+ ion and pH =========== subroutine Hionc(iout,mlax,mlay,mlaz,c_H,cc_H,pp) implicit none character(len=12) filename double precision, dimension(mlax,mlay,mlaz):: c_H,cc_H,pp integer:: iout,mlax,mlay,mlaz,l,m,n double precision :: ooo write(filename,'(a,i3.3,a)') 'c_H',iout,'.vtk' open(2,file=filename) write(2,'(a)') '# vtk DataFile Version 3.0' write(2,'(a)') 'output.vtk' write(2,'(a)') 'ASCII' write(2,'(a)') 'DATASET STRUCTURED_POINTS' write(2,'(a,3i5)') 'DIMENSIONS', mlax,mlay,mlaz write(2,'(a,3f4.1)') 'ORIGIN ', 0.0, 0.0, 0.0 write(2,'(a,3i2)') 'ASPECT_RATIO', 1, 1, 1 write(2,'(a,1i11)') 'POINT_DATA', mlax*mlay*mlaz write(2,'(a)') 'SCALARS volume_scalars double' write(2,'(a)') 'LOOKUP_TABLE default' do n=1,mlaz do m=1,mlay do l=1,mlax write(2,*) cc_H(l,m,n) enddo enddo enddo close(2) write(filename,'(a,i3.3,a)') 'pH',iout,'.vtk' open(2,file=filename) write(2,'(a)') '# vtk DataFile Version 3.0' write(2,'(a)') 'output.vtk' write(2,'(a)') 'ASCII' write(2,'(a)') 'DATASET STRUCTURED_POINTS' write(2,'(a,3i5)') 'DIMENSIONS', mlax,mlay,mlaz write(2,'(a,3f4.1)') 'ORIGIN ', 0.0, 0.0, 0.0 write(2,'(a,3i2)') 'ASPECT_RATIO', 1, 1, 1 write(2,'(a,1i11)') 'POINT_DATA', mlax*mlay*mlaz write(2,'(a)') 'SCALARS volume_scalars double' write(2,'(a)') 'LOOKUP_TABLE default' do n=1,mlaz do m=1,mlay do l=1,mlax ooo = 0. ooo = 3.-LOG10(cc_H(l,m,n)) if(pp(l,m,n).ge.0.9) ooo = 100. write(2,*) ooo enddo enddo enddo close(2) return end !========== OH- ion ============== subroutine OHionc(iout,mlax,mlay,mlaz,c_OH,cc_OH) implicit none character(len=12) filename double precision, dimension(mlax,mlay,mlaz):: c_OH,cc_OH integer:: iout,mlax,mlay,mlaz,l,m,n write(filename,'(a,i3.3,a)') 'c_OH',iout,'.vtk' open(2,file=filename) write(2,'(a)') '# vtk DataFile Version 3.0' write(2,'(a)') 'output.vtk' write(2,'(a)') 'ASCII' write(2,'(a)') 'DATASET STRUCTURED_POINTS' write(2,'(a,3i5)') 'DIMENSIONS', mlax,mlay,mlaz write(2,'(a,3f4.1)') 'ORIGIN ', 0.0, 0.0, 0.0 write(2,'(a,3i2)') 'ASPECT_RATIO', 1, 1, 1 write(2,'(a,1i11)') 'POINT_DATA', mlax*mlay*mlaz write(2,'(a)') 'SCALARS volume_scalars double' write(2,'(a)') 'LOOKUP_TABLE default' do n=1,mlaz do m=1,mlay do l=1,mlax write(2,*) cc_OH(l,m,n) enddo enddo enddo close(2) return end !============= Cl- ion =============== subroutine Clionc(iout,mlax,mlay,mlaz,c_Cl,cc_Cl) implicit none character(len=12) filename double precision, dimension(mlax,mlay,mlaz):: c_Cl,cc_Cl integer:: iout,mlax,mlay,mlaz,l,m,n write(filename,'(a,i3.3,a)') 'c_Cl',iout,'.vtk' open(2,file=filename) write(2,'(a)') '# vtk DataFile Version 3.0' write(2,'(a)') 'output.vtk' write(2,'(a)') 'ASCII' write(2,'(a)') 'DATASET STRUCTURED_POINTS' write(2,'(a,3i5)') 'DIMENSIONS', mlax,mlay,mlaz write(2,'(a,3f4.1)') 'ORIGIN ', 0.0, 0.0, 0.0 write(2,'(a,3i2)') 'ASPECT_RATIO', 1, 1, 1 write(2,'(a,1i11)') 'POINT_DATA', mlax*mlay*mlaz write(2,'(a)') 'SCALARS volume_scalars double' write(2,'(a)') 'LOOKUP_TABLE default' do n=1,mlaz do m=1,mlay do l=1,mlax write(2,*) cc_Cl(l,m,n) enddo enddo enddo close(2) return end !============ potential =================== subroutine vvv(iout,mlax,mlay,mlaz,vv) implicit none character(len=12) filename double precision, dimension(mlax,mlay,mlaz):: vv integer:: iout,mlax,mlay,mlaz,l,m,n write(filename,'(a,i3.3,a)') 'vv',iout,'.vtk' open(2,file=filename) write(2,'(a)') '# vtk DataFile Version 3.0' write(2,'(a)') 'output.vtk' write(2,'(a)') 'ASCII' write(2,'(a)') 'DATASET STRUCTURED_POINTS' write(2,'(a,3i5)') 'DIMENSIONS', mlax,mlay,mlaz write(2,'(a,3f4.1)') 'ORIGIN ', 0.0, 0.0, 0.0 write(2,'(a,3i2)') 'ASPECT_RATIO', 1, 1, 1 write(2,'(a,1i11)') 'POINT_DATA', mlax*mlay*mlaz write(2,'(a)') 'SCALARS volume_scalars double' write(2,'(a)') 'LOOKUP_TABLE default' do n=1,mlaz do m=1,mlay do l=1,mlax write(2,*) vv(l,m,n) enddo enddo enddo close(2) return end !============ initial distribution of ions =========== subroutine inidiff(nd,p,c_Cl,c_OH,c_H,cc_Cl,cc_OH,cc_H,D_Cl,D_OH,D_H,ini_loop,dxdx,dydy,dzdz,dt,pH,Ds_Fe,Ds_H,Ds_OH,Ds_Cl,De) implicit none character(len=12) filename integer:: l,m,n,lp,lm,mp,mm,nd,np,nm,loop,ini_loop double precision :: dxdx,dydy,dzdz,x1,x2,xx,y1,y2,yy,z1,z2,zz,dt,pH,abc,Ds_Fe,Ds_H,Ds_OH,Ds_Cl,De double precision, dimension(nd,nd,nd):: p,c_Cl,c_OH,c_H,cc_Cl,cc_OH,cc_H,D_Cl,D_OH,D_H,D_Fe write(*,*) '----------------Initial diffusion start!!----------------' do n=1,nd do m=1,nd do l=1,nd c_H(l,m,n) = (1.-p(l,m,n))*(10.**( 0.-pH))*1000. ![mol/L] -> [mol/m3] c_OH(l,m,n) = (1.-p(l,m,n))*(10.**(-13.+pH))*1000. ![mol/L] -> [mol/m3] c_Cl(l,m,n) = (1.-p(l,m,n))*1000. ![mol/L] -> [mol/m3] enddo enddo enddo D_Fe = 0. D_H = 0. D_OH = 0. D_Cl = 0. do n=1,nd do m=1,nd do l=1,nd abc = 3.*p(l,m,n)*p(l,m,n)-2*p(l,m,n)*p(l,m,n)*p(l,m,n) D_Fe(l,m,n) = De * abc + Ds_H * (1.-abc) D_H (l,m,n) = De * abc + Ds_H * (1.-abc) D_OH(l,m,n) = De * abc + Ds_H * (1.-abc) D_Cl(l,m,n) = De * abc + Ds_H * (1.-abc) enddo enddo enddo do loop = 1, ini_loop do n=1,nd do m=1,nd do l=1,nd lp=l+1 if(l.eq.nd) lp=nd lm=l-1 if(l.eq.1) lm=1 mp=m+1 if(m.eq.nd) mp=nd mm=m-1 if(m.eq.1) mm=1 np=n+1 if(n.eq.nd) np=nd nm=n-1 if(n.eq.1) nm=1 !============= H ============= x1= ((D_H(lp,m,n)+D_H(l ,m,n))*(c_H(lp,m,n)-c_H(l ,m,n)))/(2.*dxdx) x2=0.-((D_H(l ,m,n)+D_H(lm,m,n))*(c_H(l ,m,n)-c_H(lm,m,n)))/(2.*dxdx) xx=x1+x2 y1= ((D_H(l,mp,n)+D_H(l,m ,n))*(c_H(l,mp,n)-c_H(l,m ,n)))/(2.*dydy) y2=0.-((D_H(l,m ,n)+D_H(l,mm,n))*(c_H(l,m ,n)-c_H(l,mm,n)))/(2.*dydy) yy=y1+y2 z1= ((D_H(l,m,np)+D_H(l,m,n ))*(c_H(l,m,np)-c_H(l,m,n )))/(2.*dzdz) z2=0.-((D_H(l,m,n )+D_H(l,m,nm))*(c_H(l,m,n )-c_H(l,m,nm)))/(2.*dzdz) zz=z1+z2 cc_H(l,m,n) = c_H(l,m,n) + (xx+yy+zz)*dt if(l.eq.nd) cc_H(l,m,n)=c_H(l,m,n) !=============== OH ion =============== x1= ((D_OH(lp,m,n)+D_OH(l ,m,n))*(c_OH(lp,m,n)-c_OH(l ,m,n)))/(2.*dxdx) x2=0.-((D_OH(l ,m,n)+D_OH(lm,m,n))*(c_OH(l ,m,n)-c_OH(lm,m,n)))/(2.*dxdx) xx=x1+x2 y1= ((D_OH(l,mp,n)+D_OH(l,m ,n))*(c_OH(l,mp,n)-c_OH(l,m ,n)))/(2.*dydy) y2=0.-((D_OH(l,m ,n)+D_OH(l,mm,n))*(c_OH(l,m ,n)-c_OH(l,mm,n)))/(2.*dydy) yy=y1+y2 z1= ((D_OH(l,m,np)+D_OH(l,m,n ))*(c_OH(l,m,np)-c_OH(l,m,n )))/(2.*dzdz) z2=0.-((D_OH(l,m,n )+D_OH(l,m,nm))*(c_OH(l,m,n )-c_OH(l,m,nm)))/(2.*dzdz) zz=z1+z2 cc_OH(l,m,n) = c_OH(l,m,n) + (xx+yy+zz)*dt if(l.eq.nd) cc_OH(l,m,n) = c_OH(l,m,n) !=============== Cl ion =============== x1= ((D_Cl(lp,m,n)+D_Cl(l ,m,n))*(c_Cl(lp,m,n)-c_Cl(l ,m,n)))/(2.*dxdx) x2=0.-((D_Cl(l ,m,n)+D_Cl(lm,m,n))*(c_Cl(l ,m,n)-c_Cl(lm,m,n)))/(2.*dxdx) xx=x1+x2 y1= ((D_Cl(l,mp,n)+D_Cl(l,m ,n))*(c_Cl(l,mp,n)-c_Cl(l,m ,n)))/(2.*dydy) y2=0.-((D_Cl(l,m ,n)+D_Cl(l,mm,n))*(c_Cl(l,m ,n)-c_Cl(l,mm,n)))/(2.*dydy) yy=y1+y2 z1= ((D_Cl(l,m,np)+D_Cl(l,m,n ))*(c_Cl(l,m,np)-c_Cl(l,m,n )))/(2.*dzdz) z2=0.-((D_Cl(l,m,n )+D_Cl(l,m,nm))*(c_Cl(l,m,n )-c_Cl(l,m,nm)))/(2.*dzdz) zz=z1+z2 cc_Cl(l,m,n) = c_Cl(l,m,n) + (xx+yy+zz)*dt if(l.eq.nd) cc_Cl(l,m,n) = c_Cl(l,m,n) enddo enddo enddo c_H = cc_H c_OH= cc_OH c_Cl= cc_Cl enddo !for loop write(*,*) 'Initial diffusion is finished' return end !=========================================================