c CV simulation - quasi-reversible Redox in solution - macro electrode c D_A /= D_B c........1.........2.........3.........4.........5.........6.........7.. c23456789|123456789|123456789|123456789|123456789|123456789|123456789|12 implicit real*8 (a-h,o-z) parameter(npa=100000) dimension CX_O(npa), CX_R(npa), hx(npa) dimension alpha_O(npa), beta_O(npa), gamma_O(npa) dimension alpha_R(npa), beta_R(npa), gamma_R(npa) dimension alpha_mod_O(npa), gamma_mod_R(npa) c --------------------------------------------------------------- pi = acos(-1.0d0) write (6,*) 'pi =', pi c ---------!! input section !! change these parameters --------------- Ei = 0.3d0 ! initial potential / V Ev = -0.3d0 ! vertex potential / V vCV = 100.0d0 ! scan rate / mV s-1 Delta_E = 1.0d0 ! potential step / mV h0 = 1.0d-8 ! mesh distance closest to electrode / m omega_x = 1.01d0 ! expansion coefficient E0 = 0.0d0 ! formal potential / V alpha = 0.5d0 ! transfer coefficient BV_k0 = 1.0d-4 ! standard kinetic constant / m s-1 c_O_inf = 1.0d-3 ! initial Ox concentration / mol dm-3 D_O = 1.0d-9 ! Diffusion coefficient of Ox / m2 s-1 D_R = 1.0d-9 ! Diffusion coefficient of Red / m2 s-1 eps = 1.5d-3 ! electrode radius / m Temp = 25.0d0+273.15d0 ! temperature / K open (2,file="output.csv") c --- unit conversion vCV = vCV*1.0d-3 ! mV s-1 -> V s-1 Delta_E = Delta_E*1.0d-3 ! mV -> V c_O_inf = c_O_inf*1.0d3 ! mol dm-3 -> mol m-3 c--- physical constant & parameter Fara = 9.648533212d4 ! Faraday constant / C mol-1 Rgas = 8.314462618d0 ! gas constant / J mol-1 K-1 F_RT = Fara/(Rgas*Temp) c - Dimensionless -------------------------------------------- d_ls_O = 1.0d0 ! dimensionless diffusion coefficient of Ox d_ls_R = D_R/D_O ! dimensionless diffusion coefficient of Red h0_ls = h0/eps ! dimensionless mesh distance closest to electrode BV_K0_ls = BV_k0*eps/D_O ! dimensionless standard kinetic constant theta_i = F_RT*(Ei-E0) theta_v = F_RT*(Ev-E0) Delta_theta = F_RT*Delta_E sigma = eps**2.0d0/D_O*F_RT*vCV ! dimensionless scan rate write (6,*) 'd_ls_R =', d_ls_R write (6,*) 'h0_ls, BV_K0_ls =', h0_ls, BV_K0_ls write (6,*) 'theta_i, theta_v =', theta_i, theta_v write (6,*) 'Delta_theta, sigma =', Delta_theta, sigma c - Calculate Delta_T, m_T --------------------------------- Delta_T = Delta_theta/sigma T_max = 2.0d0*(abs(theta_i - theta_v))/sigma ! whole step number mT = int(T_max/Delta_T) T_switch = abs(theta_i - theta_v)/sigma mT_v = int(T_switch/Delta_T) ! step number at vertex potential write (6,*) 'Delta_T, T_max =', Delta_T, T_max write (6,*) 'mT_v =', mT_v c - Calculate n_X ---------------------------------------- X_max = 6.0d0*sqrt(max(1.0d0,d_ls_R)*(T_max)) icount = 1 hx(1) = h0_ls hx_max = hx(1) do i=2,npa ! *1 icount = i hx(i) = h0_ls*omega_x**real(i-1,kind(0d0)) hx_max = hx_max + hx(i) if (hx_max > X_max) exit end do nX = icount ! *2 if (nX > npa) then stop 1550 end if write (6,*) 'X_max, nX =', X_max, nX c ------------------------------------------------------------ c --- begin simulation --------------------------------------- c ------------------------------------------------------------ write (2,*) 'potential /V', ',', 'current density / mA cm-2' ! *3 CX_O = 1.0d0 ! *4 CX_R = 0.0d0 c --- CV start --------------------------------------------- do k=1, mT+1 ! *5 if (k == 1) then ! *6 theta = theta_i else if (2 <= k .and. k <= mT_v+1) then theta = theta - Delta_theta else theta = theta + Delta_theta end if c - Calculate Thomas coefficients@*7 alpha_O = 0.0d0 beta_O = 0.0d0 gamma_O = 0.0d0 alpha_R = 0.0d0 beta_R = 0.0d0 gamma_R = 0.0d0 c --- i = 1 alpha_O(1) = -h0_ls*exp((1-alpha)*theta)*BV_K0_ls beta_O(1) = 1.0d0 + h0_ls*exp(-alpha*theta)*BV_K0_ls gamma_O(1) = -1.0d0 alpha_R(1) = -h0_ls*exp(-alpha*theta)*BV_K0_ls/d_ls_R beta_R(1) = 1.0d0 + & h0_ls*exp((1-alpha)*theta)*BV_K0_ls/d_ls_R gamma_R(1) = -1.0d0 c --- i = 2 to n-1 do i=2, nX-1 alpha_O(i) = -2.0d0*Delta_T/(hx(i-1)*(hx(i)+hx(i-1))) gamma_O(i) = -2.0d0*Delta_T/(hx(i)*(hx(i)+hx(i-1))) beta_O(i) = -alpha_O(i) - gamma_O(i) + 1.0d0 alpha_R(i) = -2.0d0*Delta_T*d_ls_R & /(hx(i-1)*(hx(i)+hx(i-1))) gamma_R(i) = -2.0d0*Delta_T*d_ls_R & /(hx(i)*(hx(i)+hx(i-1))) beta_R(i) = -alpha_R(i) - gamma_R(i) + 1.0d0 end do c --- i = n alpha_O (nx) = 0.0d0 beta_O(nx) = 1.0d0 gamma_O(nx) = 0.0d0 alpha_R(nx) = 0.0d0 beta_R(nx) = 1.0d0 gamma_R(nx) = 0.0d0 c --- modified alpha, modified gamma *8 alpha_mod_O = 0.0d0 gamma_mod_R = 0.0d0 alpha_mod_O(nX) = 0.0d0 do i = nX-1, 1, -1 alpha_mod_O(i) = alpha_O(i) & /(beta_O(i) - gamma_O(i)*alpha_mod_O(i+1)) end do gamma_mod_R(1) = gamma_R(1) & /(beta_R(1) - alpha_R(1)*alpha_mod_O(1)) do i = 2, nX-1 gamma_mod_R(i) = gamma_R(i) & /(beta_R(i) - alpha_R(i)*gamma_mod_R(i-1)) end do c --- delta_mod, CX calc *9 CX_O(1) = 0.0d0 CX_R(1) = 0.0d0 c --- delta_O_mod *10 CX_O(nX) = 1.0d0 do i= nX-1, 1, -1 CX_O(i) = (CX_O(i) - gamma_O(i)*CX_O(i+1)) & /(beta_O(i) - gamma_O(i)*alpha_mod_O(i+1)) end do c --- delta_R_mod CX_R(1) = (CX_R(1) - alpha_R(1)*CX_O(1)) & /(beta_R(1) - alpha_R(1)*alpha_mod_O(1)) do i = 2, nX CX_R(i) = (CX_R(i) - alpha_R(i)*CX_R(i-1)) & /(beta_R(i) - alpha_R(i)*gamma_mod_R(i-1)) end do c --- back substitution *11 CX_R(nX) = 0.0d0 do i = nX-1, 1, -1 CX_R(i) = CX_R(i) - gamma_mod_R(i)*CX_R(i+1) end do CX_O(1) = CX_O(1) - alpha_mod_O(1)*CX_R(1) do i = 2, nX CX_O(i) = CX_O(i) - alpha_mod_O(i)*CX_O(i-1) end do if (CX_O(nX) /= 1.0d0) then write(6,*) 'check h0 or omega' stop 908 end if c ------ output current ----------------------------------------- cur_ls = -(CX_O(2) - CX_O(1))/h0_ls ! *12 potential = theta/F_RT+E0 ! *13 cur = cur_ls*(Fara*c_O_inf*D_O)/eps*1.0d3/1.0d4 ! mA cm-2 write(2,*) potential, ',', cur c --- CV finish --------------------------------------------- end do ! *5 close(2) stop end