building half cell model by using DFN model with option: "working electrode" as "negative"

[ChuljungKim89]

Hi all,

I am working with my master thesis at the moment and not familiar with python.

But I am trying to build a half cell model with pybamm, namely with DFN model. With option for working electrode "negative", I get an keyerror as follows:
KeyError: "'Maximum concentration in negative electrode [mol.m-3]' not found. Best matches are ['Maximum concentration in positive electrode [mol.m-3]', 'Initial concentration in positive electrode [mol.m-3]', 'Initial concentration in electrolyte [mol.m-3]']"

but when I use option for working electrode as "positive", as in example note introduced, I don't get any error but it is about cathod.
I saw in the "BasicDFNHalfCell" code file, that for anode I just have to put parameters for anode, instead of cathod. Is it so simple or do I have to define all parameters for anode by myself by using "pybamm.Pameter("Parameter name")" and "pybamm.ParameterValues({"Parameter name": "value"}"?

---

And I have another question about geometry.
As I mentioned above, I am trying to use coin cell geometry and I have tried to build my model based on the example notes in PyBaMM repository. But the thing is, that I am not sure whether my geometry is well implemented or not. Could anyone have a look at the code and tell me, if it is well defined and implemented or not? If not, what would be the problem?
(At this point, I have no idea what "Electrode cross-sectional area [m2]" meant. Could anyone explain, what it is about and if it is also for coin cell relevant? When I erase this parameter value, I get an error by solving problem)

---

Below is my script for anode half cell model for coin cell geometry.
When I use model with "pybamm.BaseModel( )" with this code, I have no problem solving, but if I try to use model "pybamm.lithium_ion.DFN( )" with option for working electrode "positive" and parameters for graphite anode half cell, I have also another error "RecursionError: maximum recursion depth exceeded".

Thanks for your help ahead.

import pybamm
import numpy as np
import matplotlib.pyplot as plt
pybamm.set_logging_level("INFO")

halfcell = {"working electrode": "positive"}

#model_gr = pybamm.BaseModel()
model_gr = pybamm.lithium_ion.DFN(options = halfcell)

define electric potential in the working electrode.

working electrode is discribed as "cathode" but

for anode as working electrode, you can just put

paramters for anode to cathode parameter values

phi = pybamm.Variable("Positive electrode potential [V]", domain="positive electrode")

phi_e_s = pybamm.Variable("Separator electrolyte potential [V]", domain="separator")
phi_e_p = pybamm.Variable("Positive electrolyte potential [V]", domain="positive electrode")

phi_e = pybamm.concatenation(phi_e_s, phi_e_p)

defining concentration variable: lithium concentration in electrode particle

c = pybamm.Variable(
"Positive particle concentration [mol.m-3]",
domain="positive particle",
auxiliary_domains={
"secondary": "positive electrode",
},
)

some constants

#from scipy import constants
F = pybamm.Parameter("Faraday constant [C.mol-1]")
R = pybamm.Parameter("Molar gas constant [J.mol-1.K-1]")
T = pybamm.Parameter("Ambient temperature [K]")
T0 = pybamm.Parameter("Initial temperature [K]")

a = pybamm.Parameter("Surface area per unit volume [m-1]")
R_p = pybamm.Parameter("Positive particle radius [m]")
L_s = pybamm.Parameter("Separator thickness [m]")
L_p = pybamm.Parameter("Positive electrode thickness [m]")
A = pybamm.Parameter("Electrode cross-sectional area [m2]")
R_curr = pybamm.Parameter("Current collector radius [m]")

sigma = pybamm.Parameter("Positive electrode conductivity [S.m-1]")
kappa = pybamm.Parameter("Electrolyte conductivity [S.m-1]")
D = pybamm.Parameter("Diffusion coefficient [m2.s-1]")

I_app = -pybamm.Parameter("Applied current [A]")
c0 = pybamm.Parameter("Initial concentration in positive electrode [mol.m-3]")
c_max = pybamm.Parameter("Maximum concentration in positive electrode [mol.m-3]")
capacity = pybamm.Parameter("Nominal cell capacity [A.h]")

function parameters

c_surf = pybamm.surf(c) # get the surface concentration
inputs = {"Positive particle surface concentration [mol.m-3]": c_surf}
j0 = pybamm.FunctionParameter("Positive electrode exchange-current density [A.m-2]", inputs)
U = pybamm.FunctionParameter("Positive electrode OCP [V]", inputs)

basic kinetics and variables for the equations

j_s = pybamm.PrimaryBroadcast(0, "separator")
j_p = 2 * j0 * pybamm.sinh((F / 2 / R / T) * (phi - phi_e_p - U))
j = pybamm.concatenation(j_s, j_p)

governing equations

charge conservation equations

i = -sigma * pybamm.grad(phi)
i_e = -kappa * pybamm.grad(phi_e)

model_gr.algebraic = {
phi: pybamm.div(i) + a * j_p,
phi_e: pybamm.div(i_e) - a * j,
}

particle equations (mass conservation)

N = -D * pybamm.grad(c) # flux
dcdt = -pybamm.div(N)
model_gr.rhs = {c: dcdt}

boundary conditions

model_gr.boundary_conditions = {
phi: {"left": (pybamm.Scalar(1), "Neumann"), "right": (-I_app / A / sigma, "Neumann")},
phi_e: {"left": (pybamm.Scalar(0), "Dirichlet"), "right": (pybamm.Scalar(1), "Neumann")},
c: {"left": (pybamm.Scalar(0), "Neumann"), "right": (-j_p / F / D, "Neumann")}
}

initial conditions

inputs = {"Initial concentration [mol.m-3]": c0}
U_init = pybamm.FunctionParameter("Positive electrode OCP [V]", inputs)
model_gr.initial_conditions = {
phi: U_init,
phi_e: 0,
c: c0
}

add variables I want to know

model_gr.variables = {
"Positive electrode potential [V]": phi,
"Electrolyte potential [V]": phi_e,
"Positive particle concentration [mol.m-3]": c,
"Positive particle surface concentration [mol.m-3]": c_surf,
"Average positive particle surface concentration [mol.m-3]": pybamm.x_average(c_surf),
"Positive electrode interfacial current density [A.m-2]": j_p,
"Positive electrode OCP [V]":pybamm.boundary_value(U, "right"),
"Voltage [V]": pybamm.boundary_value(phi, "right"),
}

######## Using the model #########

both functions will depend on the maximum concentration

from pybamm import exp, tanh

def exchange_current_density_gr(c_surf):
k = 6.48 * 10 **(-7) # reaction rate [(A/m2)(m3/mol)1.5] LG50 Chen2020
c_e = 1000 # (constant) electrolyte concentration [mol.m-3]
# Arrhenius effect is neglected, because of small cell and low C-rate
return k * c_e 0.5 * c_surf ** 0.5 * (c_max - c_surf) ** 0.5

Below is for graphite OCV

def open_circuit_potential_gr(c_surf):
stretch = 1.062
sto = stretch *c_surf / c_max
print("sto = " ,sto)

u_eq = (
    1.9793 * exp(-39.3631 * sto)
    + 0.2482
    - 0.0909 * tanh(29.8538 * (sto - 0.1234))
    - 0.04478 * tanh(14.9159 * (sto - 0.2769))
    - 0.0205 * tanh(30.4444 * (sto - 0.6103))
)

return u_eq

def applied_current_density(C_rates):
return -(C_rates * capacity)

C_rates = 0.001 # [1/h]
#C_rates = [0.005, 0.1, 0.2] # [1/h]

determining charging time for each C-rate

'''
charging_times = []
h_evals = []
for i in range(len(C_rates)):
charging_time = (1 / C_rates[i]) # in [s]
h_evals.append(np.linspace(0, (charging_time / charging_time * 100), 600))
charging_times.append(np.linspace(0, charging_time, 600)*3600)
print("evaluation time for C-rate="+str(C_rates[i]), ": ", charging_time, ' hours')
'''
charging_times = (1/C_rates) # in [s]
h_evals = np.linspace(0, (charging_times/charging_times)*100, 600)

print("evaluation time check: ", charging_times)
print("SoC check: ", h_evals)

I_app = applied_current_density(C_rates) # reference current value for C/10

param_gr = pybamm.ParameterValues(
{
"Surface area per unit volume [m-1]": 0.37e6,
"Positive particle radius [m]": 10e-6,
"Separator thickness [m]": 260e-6,
"Positive electrode thickness [m]": 45.5e-6,
"Electrode cross-sectional area [m2]": 2.8e-2,
"Applied current [A]": I_app,
"Positive electrode conductivity [S.m-1]": 13,
"Electrolyte conductivity [S.m-1]": 1e-2,
"Diffusion coefficient [m2.s-1]": 3.9e-14,
"Faraday constant [C.mol-1]": 96485,
"Initial concentration in positive electrode [mol.m-3]": 650,
"Molar gas constant [J.mol-1.K-1]": 8.314,
"Ambient temperature [K]": 298.15,
"Maximum concentration in positive electrode [mol.m-3]": 51217,
"Positive electrode exchange-current density [A.m-2]": exchange_current_density_gr,
"Positive electrode OCP [V]": open_circuit_potential_gr,
"Current collector radius [m]": 0.8e-3,
"Nominal cell capacity [A.h]": 97e-3,
"Initial temperature [K]": T0
}
)

Geometry definition

r = pybamm.SpatialVariable(
"r",
domain=["positive particle"],
auxiliary_domains={
"secondary": "positive electrode"
},
coord_sys="spherical polar")

r_cc = pybamm.SpatialVariable(
"r_cc",
domain=["current collector"],
coord_sys="spherical polar")

x_s = pybamm.SpatialVariable("x_s", domain=["separator"], coord_sys="cartesian")
x_p = pybamm.SpatialVariable("x_p", domain=["positive electrode"], coord_sys="cartesian")

geometry = {
"separator": {x_s: {"min": -L_s, "max": 0}},
"positive electrode": {x_p: {"min": 0, "max": L_p}},
"positive particle": {r: {"min": 0, "max": R_p}},
"current collector": {r_cc: {"min": 0, "max": R_curr}}
}

param_gr.process_model(model_gr)
param_gr.process_geometry(geometry)

submesh_types = {
"separator": pybamm.Uniform1DSubMesh,
"positive electrode": pybamm.Uniform1DSubMesh,
"positive particle": pybamm.Uniform1DSubMesh,
"current collector": pybamm.Uniform1DSubMesh,
}
var_pts = {x_s: 10, x_p: 20, r: 30, r_cc: 10}
mesh = pybamm.Mesh(geometry, submesh_types, var_pts)

spatial_methods = {
"separator": pybamm.FiniteVolume(),
"positive electrode": pybamm.FiniteVolume(),
"positive particle": pybamm.FiniteVolume(),
"current collector": pybamm.FiniteVolume(),
}
disc = pybamm.Discretisation(mesh, spatial_methods)
disc.process_model(model_gr);

solve

solver = pybamm.CasadiSolver(mode='safe', root_tol=1e-3)

solutions = solver.solve(model_gr, t_eval=charging_times)

#plot

V_gr = solutions["Positive electrode OCP [V]"].entries
plt.plot(h_evals, V_gr, label = "$C_\mathrm{rate}=$"+str(C_rates))
plt.xlabel("State of Charge [%]")
plt.ylabel("OCV [V]")
plt.title("Graphite anode half-cell")
plt.legend()

[valentinsulzer]

See #3198