from scipy.sparse import coo_matrix from scipy.sparse.linalg import spsolve from numpy import * from pylab import * # PyBand Porous Tafel (c) by Joshua Gallaway # PyBand Porous Tafel is licensed under a # Creative Commons Attribution-ShareAlike 3.0 Unported License. # This is set up to solve the nonlinear # boundary value problem for a porous # electrode in cartesian coordinates # with Tafel kinetics # N is number of species # NJ is number of spatial points # electronic current density i1 is c1 # ionic current density i2 is c2 # electronic potential phi1 is c3 # ionic potential phi2 is c4 N = 4 NJ = 100 # Set convergence tolerance and max number of iterations tol = 1e-12 itmax = 100 # Constants for the problem I = 0.1 # A/cm2 L = 1.0 # cm, enter as a float a = 23300 # 1/cm io = 2e-7 # A/cm2 n = 1 F = 96500 # C/mol R = 8.314 # J/mol-K T = 298 # K sigma = 20 # S/cm kappa = 0.06 # S/cm alpha = 0.5 beta = -(1-alpha)*n*F/R/T # Makes the mesh xx = linspace(0,L,NJ) h = L/(NJ-1) # =========== # INITIAL GUESSES # =========== def initguess(): cold = ones([N,NJ]) #cold[0,:]=0.05 #cold[1,:]=0.05 #cold[2,:]=-0.1 #cold[3,:]=0.0 cold[0,:]=5 cold[1,:]=5 cold[2,:]=5 cold[3,:]=5 # The initial guesses shouldn't matter return cold # =========== # FILLMAT # =========== def fillmat(cold): # first column is equation # second column is position # third column is species # initialize matrices sma = zeros([N,NJ,N]) smb = zeros([N,NJ,N]) smd = zeros([N,NJ,N]) smg = zeros([N,NJ]) # fill matrices smb[0,:,0] = 1 smb[0,:,1] = 1 smb[1,:,1] = -1 smb[2,:,2] = sigma smb[3,:,3] = kappa smd[1,:,2] = -beta*a*io*exp(beta*(cold[2,:]-cold[3,:])) smd[1,:,3] = beta*a*io*exp(beta*(cold[2,:]-cold[3,:])) smd[2,:,0] = 1 smd[3,:,1] = 1 smg[1,:] = (1-beta*cold[2,:]+beta*cold[3,:])*a*io*exp(beta*(cold[2,:]-cold[3,:])) # __________Boundary condition 1 smp = zeros([N,N]) sme = zeros([N,N]) smf = zeros([N,1]) sme[1,1] = 1 sme[3,3] = 1 smf[1] = I smp[2,2] = 1 # others smp[0,0] = 1 smp[0,1] = 1 # Insert (sme smp smf) into (smb smd smg) smb[:,0,:] = smp[:,:] smd[:,0,:] = sme[:,:] smg[:,0] = transpose(smf) # __________Boundary condition 2 smp = zeros([N,N]) sme = zeros([N,N]) smf = zeros([N,1]) sme[0,0] = 1 smf[0] = I # others #smp[1,1] = -1 #sme[1,2] = -beta*a*io*exp(beta*(cold[2,NJ-1]-cold[3,NJ-1])) #sme[1,3] = beta*a*io*exp(beta*(cold[2,NJ-1]-cold[3,NJ-1])) sme[1,1] = 1 smp[2,2] = sigma sme[2,0] = 1 #smf[1] = (1-beta*cold[2,NJ-1]+beta*cold[3,NJ-1])*a*io*exp(beta*(cold[2,NJ-1]-cold[3,NJ-1])) smf[1] = 0 smp[3,3] = kappa sme[3,1] = 1 smf[3] = 0 # Insert (sme smp smf) into (smb smd smg) smb[:,NJ-1,:] = smp[:,:] smd[:,NJ-1,:] = sme[:,:] smg[:,NJ-1] = transpose(smf) return sma, smb, smd, smg # =========== # ABDGXY # =========== def abdgxy(sma, smb, smd, smg): sma = transpose(sma, (0, 2, 1)) smb = transpose(smb, (0, 2, 1)) smd = transpose(smd, (0, 2, 1)) A = sma-h/2.0*smb B = -2.0*sma+h**2*smd D = sma+h/2.0*smb G = h**2*smg B[:,:,0] = h*smd[:,:,0]-1.5*smb[:,:,0] D[:,:,0] = 2.0*smb[:,:,0] G[:,0]=h*smg[:,0] X = -0.5*smb[:,:,0] A[:,:,NJ-1]=-2.0*smb[:,:,NJ-1] B[:,:,NJ-1]=h*smd[:,:,NJ-1]+1.5*smb[:,:,NJ-1] G[:,NJ-1]=h*smg[:,NJ-1] Y=0.5*smb[:,:,NJ-1] ABD = concatenate((A, B, D), axis=1) BC1 = concatenate((B[:,:,0] , D[:,:,0] , X), axis=1) BC2 = concatenate((Y , A[:,:,NJ-1] , B[:,:,NJ-1]), axis=1) ABD[:,:,0] = BC1 ABD[:,:,NJ-1] = BC2 return ABD, G # =========== # BAND # =========== def band(ABD, G): # If you have Newman's BAND code compiled in fortran you could # invoke that here to solve the matrix problem. # Here we use SciPy's matrix problem solver, which # is a bit slower. BMrow = reshape(arange(1,N*NJ+1), (NJ,N)) BMrow = BMrow[:, :, newaxis] BMrow = transpose(BMrow, (1, 2, 0)) BMrow = BMrow[:,[0 for i in range(3*N)],:] a = arange(1,3*N+1) a = a[newaxis,:] a = repeat(a,N,0) a = a[:,:,newaxis] a = repeat(a,NJ,2) b = arange(0,(N)*(NJ-3)+N,N) b = hstack((b[0], b, b[len(b)-1])) b = b[newaxis,newaxis,:] b = repeat(b,N,0) b = repeat(b,3*N,1) BMcol = a + b BMcol = BMcol - 1 BMrow = BMrow - 1 BMrow = ravel(BMrow) BMcol = ravel(BMcol) ABD = ravel(ABD) BigMat = coo_matrix((ABD, (BMrow, BMcol)), shape=(N*NJ, N*NJ)).tocsc() BigG = transpose(G) BigG = ravel(BigG) # print BigMat.todense() delc = spsolve(BigMat, BigG) return delc # =========== # MAIN PROGRAM # =========== cold = initguess() it = 1 did = False for it in range(1,itmax): if it in (1,5,10,15,20,25,30,35,40): # You could plot selected iterations this way # to watch the problem converge pass sma, smb, smd, smg = fillmat(cold) ABD, G = abdgxy(sma, smb, smd, smg) delc = band(ABD, G) delc = delc.reshape((NJ, N)) delc = transpose(delc) error = cold - delc maxerror = amax(abs(error)) print it, ', ', maxerror cold = delc if maxerror < tol: did = True print 'Converged in ' +str(it)+ ' steps.' break if not did: print 'The program did not converge.' # =========== # PLOT RESULTS # =========== fig, (ax1,ax2) = plt.subplots(2,1,figsize=(6, 6)) fig.subplots_adjust(left=0.15) params = {'mathtext.default': 'regular' } rcParams.update(params) xs = xx c1 = cold[0,:] c2 = cold[1,:] c3 = cold[2,:] c4 = cold[3,:] ax1.plot(xs,c1,'-o', color='blue', markeredgecolor='k') ax1.plot(xs,c2,'-s', color='green', markeredgecolor='k') ax1.legend(['$i_{1}$','$i_{2}$'],loc='center right') ax1.set_ylabel('$I\/(A/cm^{2})$') ax2.plot(xs,c3,'-o', color='blue', markeredgecolor='k') ax2.plot(xs,c4,'-s', color='green', markeredgecolor='k' ) ax2.legend(['$\phi_{1}$','$\phi_{2}$'],loc='center right') #xlabel('x / cm') ax2.set_xlabel('$x\/(cm)$') ax2.set_ylabel('$\phi\/(V)$') ax1.set_xlim(0, L) ax2.set_xlim(0, L) ax1.tick_params(axis='both', direction='in', bottom=True, top=True, left=True, right=True) ax2.tick_params(axis='both', direction='in', bottom=True, top=True, left=True, right=True) show() saver = [xs, c1, c2, c3, c4] saver = zip(*saver) savetxt('porous1.csv', saver, delimiter=",") savefig('porous1.png', dpi=300) #close('all')