Hello, I am trying to create SPMe model using the final equations from Marquis2019 with the help of pybamm.base_model() to get a better feel of creating models in pybamm. But I am unable to match the electrolyte dynamics with pybamm.lithium_ion.SPMe(). I have been trying to debug it for some time but to no avail. Would very much appreciate any additional input onto this. The code I used for this is as follows:
import pybamm
import numpy as np
import matplotlib.pyplot as plt
from Chen_2020 import *
model = pybamm.BaseModel()
####################Parameter and Variable definitions########################################
epsl_n = pybamm.Parameter('Porosity negative')
epsl_p = pybamm.Parameter('Porosity positive')
epsl_s = pybamm.Parameter('Porosity separator')
tor_n = pybamm.Parameter('Transport factor negative')
tor_p = pybamm.Parameter('Transport factor positive')
tor_s = pybamm.Parameter('Transport factor separator')
# De_n = pybamm.Parameter('Diffusivity negative [m2.s-1]')
De_p = pybamm.Parameter('Diffusivity positive [m2.s-1]')
De_s = pybamm.Parameter('Diffusivity separator [m2.s-1]')
L_n = pybamm.Parameter('Negative electrode thickness [m]')
L_p = pybamm.Parameter("Positive electrode thickness [m]")
L_s = pybamm.Parameter('Separator thickness [m]')
t = pybamm.Parameter('Transference number')
F= pybamm.Parameter('Faraday constant')
Iapp = pybamm.Parameter('Applied current [A]')
A_cell = pybamm.Parameter('Cell area [m2]')
ce_typ = pybamm.Parameter("Typical electrolyte concentration [mol.m-3]")
ce_n= pybamm.Variable("Negative", domain='negative electrode')
ce_p = pybamm.Variable("Positive", domain='positive electrode')
ce_s = pybamm.Variable("Separator", domain = 'separator')
###########Diffusivity functionss########################################################33
def De_n(ce_n):
Dn = electrolyte_diffusivity_Nyman2008(ce_n,298.15)
return Dn
def diffu_separator(ce_s):
Ds = electrolyte_diffusivity_Nyman2008(ce_s, 298.15)
return Ds
def diffu_positive(ce_p):
Dp = electrolyte_diffusivity_Nyman2008(ce_p, 298.15)
return Dp
diffu_pos = 8.794e-11* (pybamm.surf(ce_p)**2) - 3.972e-10* (pybamm.surf(ce_p) )+4.862e-10
diffu_neg = 8.794e-11* (pybamm.surf(ce_n)**2) - 3.972e-10* (pybamm.surf(ce_p) )+4.862e-10
diffu_sep = 8.794e-11* (pybamm.surf(ce_s)**2) - 3.972e-10* (pybamm.surf(ce_s) )+4.862e-10
#####################Definition of equations and model#########################################
Ne_n = -tor_n * diffu_neg * pybamm.grad(ce_n)
dcen_dt = -pybamm.div(Ne_n) / epsl_n + (1-t)*Iapp/(F*L_n*A_cell*epsl_n)
Ne_p = -tor_p *diffu_pos* pybamm.grad(ce_p)
dcep_dt = -pybamm.div(Ne_p)/ epsl_p + (t-1)*Iapp/(F*L_p*A_cell*epsl_p)
Ne_s = -tor_s*diffu_sep*pybamm.grad(ce_s)
dces_dt = -pybamm.div(Ne_s)/epsl_p
ce = pybamm.concatenation(ce_n, ce_s, ce_p)
dce_dt = pybamm.concatenation(dcen_dt, dces_dt, dcep_dt)
model.rhs={ce:dce_dt}
model.boundary_conditions= {
ce:{
'left':(pybamm.Scalar(0), 'Neumann'),
'right':(pybamm.Scalar(0),'Neumann')
}
}
model.initial_conditions = {ce: ce_typ}
###################Declaration of parameters and variables ################################
model.variables ={
'Negative': ce_n,
"Positive": ce_p,
"Separator": ce_s,
"Negative average":pybamm.XAverage(ce_n),
"Positive average": pybamm.XAverage(ce_p),
"Separator average": pybamm.XAverage(ce_s)
}
param = pybamm.ParameterValues(
{
'Porosity negative': 0.25,
'Porosity positive': 0.335,
'Porosity separator': 0.47,
'Transport factor negative': 0.25**1.5,
'Transport factor positive': 0.335 **1.5,
'Transport factor separator': 0.47 ** 1.5,
# 'Diffusivity negative [m2.s-1]': 8.794e-11 - 3.972e-10 + 4.862e-10,
'Diffusivity positive [m2.s-1]': 8.794e-11 - 3.972e-10 + 4.862e-10,
'Diffusivity separator [m2.s-1]': 8.794e-11 - 3.972e-10 + 4.862e-10,
'Negative electrode thickness [m]': 8.52e-5,
'Positive electrode thickness [m]': 7.56e-5,
'Separator thickness [m]': 1.2e-5,
'Transference number': 0.2594,
'Faraday constant': 96485.332,
'Cell area [m2]': 0.14,
"Typical electrolyte concentration [mol.m-3]" : 1000,
"Applied current [A]": -5
}
)
################Describing the geometry#########################
# intf1 = pybamm.Scalar(0.00008752)
# intf2= pybamm.Scalar(0.0000972)
# end = pybamm.Scalar(0.0001728)
xneg = pybamm.SpatialVariable('xneg', domain = {'primary': ['negative electrode']}, coord_sys='cartesian')
xpos = pybamm.SpatialVariable('xpos', domain = {'primary': ['positive electrode']}, coord_sys='cartesian')
xsep = pybamm.SpatialVariable('xsep', domain = {'primary': ['separator']}, coord_sys='cartesian')
geometry = {
'negative electrode': {xneg:{'min': pybamm.Scalar(0), 'max':pybamm.Scalar(1)}}, # 8.52e-5
'separator':{xsep:{'min':pybamm.Scalar(1), 'max':pybamm.Scalar(2)}}, # 9.72e-5
'positive electrode':{xpos:{'min':pybamm.Scalar(2), 'max':pybamm.Scalar(3)}} #17.28e-5
}
param.process_model(model)
param.process_geometry(geometry)
submesh_types ={
'negative electrode': pybamm.Uniform1DSubMesh,
'separator':pybamm.Uniform1DSubMesh,
'positive electrode': pybamm.Uniform1DSubMesh
}
var_pts = {xneg:20, xsep:20, xpos :20}
mesh = pybamm.Mesh(geometry, submesh_types, var_pts)
spatial_methods ={
'negative electrode':pybamm.FiniteVolume(),
'separator':pybamm.FiniteVolume(),
'positive electrode':pybamm.FiniteVolume()
}
disc = pybamm.Discretisation(mesh, spatial_methods)
disc.process_model(model)
solver= pybamm.ScipySolver()
t= np.linspace(0,3600)
solution = solver.solve(model,t)
ce_p = np.array(solution['Positive'].entries)
ce_s = np.array(solution['Separator'].entries)
ce_n = np.array(solution['Negative'].entries)
cp_avg = np.array(solution['Positive average'].entries)
cs_avg = np.array(solution['Separator average'].entries)
cn_avg= np.array(solution['Negative average'].entries)
plt.plot(t,cp_avg, color='blue')
plt.plot(t,cs_avg,color='green')
plt.plot(t,cn_avg, color='red')
# for i in range(0,20):
# plt.plot(t,ce_p[i], color='blue')
# plt.plot(t,ce_s[i],color='green')
# plt.plot(t,ce_n[i], color='red')
plt.show()
