Tip

An interactive online version of this notebook is available, which can be accessed via # A step towards the Single Particle Model#

In the previous notebook we saw how to solve a PDE model in pybamm. Now it is time to solve a real-life battery problem! We consider the problem of spherical diffusion in the negative electrode particle within the single particle model. That is,

$\frac{\partial c}{\partial t} = \nabla \cdot (D \nabla c),$

with the following boundary and initial conditions:

$\left.\frac{\partial c}{\partial r}\right\vert_{r=0} = 0, \quad \left.\frac{\partial c}{\partial r}\right\vert_{r=R} = -\frac{j}{FD}, \quad \left.c\right\vert_{t=0} = c_0,$

where $$c$$ is the concentration, $$r$$ the radial coordinate, $$t$$ time, $$R$$ the particle radius, $$D$$ the diffusion coefficient, $$j$$ the interfacial current density, $$F$$ Faraday’s constant, and $$c_0$$ the initial concentration.

In this example we use the following parameters:

Symbol

Units

Value

$$R$$

m

$$10 \times 10^{-6}$$

$$D$$

m$${^2}$$ s$$^{-1}$$

$$3.9 \times 10^{-14}$$

$$j$$

A m$$^{-2}$$

$$1.4$$

$$F$$

C mol$$^{-1}$$

$$96485$$

$$c_0$$

mol m$$^{-3}$$

$$2.5 \times 10^{4}$$

## Setting up the model#

As before, we begin by importing the PyBaMM library into this notebook, along with any other packages we require, and start with an empty pybamm.BaseModel

:

%pip install pybamm[plot,cite] -q    # install PyBaMM if it is not installed
import pybamm
import numpy as np
import matplotlib.pyplot as plt

model = pybamm.BaseModel()

Note: you may need to restart the kernel to use updated packages.


We then define all of the model variables and parameters. Parameters are created using the pybamm.Parameter class and are given informative names (with units). Later, we will provide parameter values and the Parameter objects will be turned into numerical values. For more information please see the parameter values notebook.

:

R = pybamm.Parameter("Particle radius [m]")
D = pybamm.Parameter("Diffusion coefficient [m2.s-1]")
j = pybamm.Parameter("Interfacial current density [A.m-2]")
c0 = pybamm.Parameter("Initial concentration [mol.m-3]")

c = pybamm.Variable("Concentration [mol.m-3]", domain="negative particle")


Now we define our model equations, boundary and initial conditions, as in the previous example.

:

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

# boundary conditions
lbc = pybamm.Scalar(0)
rbc = -j / F / D
model.boundary_conditions = {c: {"left": (lbc, "Neumann"), "right": (rbc, "Neumann")}}

# initial conditions
model.initial_conditions = {c: c0}


Finally, we add any variables of interest to the dictionary model.variables

:

model.variables = {
"Concentration [mol.m-3]": c,
"Surface concentration [mol.m-3]": pybamm.surf(c),
"Flux [mol.m-2.s-1]": N,
}


## Using the model#

In order to discretise and solve the model we need to provide values for all of the parameters. This is done via the pybamm.ParameterValues class, which accepts a dictionary of parameter names and values

:

param = pybamm.ParameterValues(
{
"Diffusion coefficient [m2.s-1]": 3.9e-14,
"Interfacial current density [A.m-2]": 1.4,
"Initial concentration [mol.m-3]": 2.5e4,
}
)


Here all of the parameters are simply scalars, but they can also be functions or read in from data (see parameter values notebook).

As in the previous example, we define the particle geometry. Note that in this example the definition of the geometry contains a parameter, the particle radius $$R$$

:

r = pybamm.SpatialVariable("r", domain=["negative particle"], coord_sys="spherical polar")
geometry = {"negative particle": {r: {"min": pybamm.Scalar(0), "max": R}}}


Both the model and geometry can now be processed by the parameter class. This replaces the parameters with the values

:

param.process_model(model)
param.process_geometry(geometry)


We can now set up our mesh, choose a spatial method, and discretise our model

:

submesh_types = {"negative particle": pybamm.Uniform1DSubMesh}
var_pts = {r: 20}
mesh = pybamm.Mesh(geometry, submesh_types, var_pts)

spatial_methods = {"negative particle": pybamm.FiniteVolume()}
disc = pybamm.Discretisation(mesh, spatial_methods)
disc.process_model(model);


The model is now discretised and ready to be solved.

### Solving the model#

As is the previous example, we choose a solver and times at which we want the solution returned.

:

# solve
solver = pybamm.ScipySolver()
t = np.linspace(0, 3600, 600)
solution = solver.solve(model, t)

# post-process, so that the solution can be called at any time t or space r
# (using interpolation)
c = solution["Concentration [mol.m-3]"]
c_surf = solution["Surface concentration [mol.m-3]"]

# plot
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13, 4))

ax1.plot(solution.t, c_surf(solution.t))
ax1.set_xlabel("Time [s]")
ax1.set_ylabel("Surface concentration [mol.m-3]")

r = mesh["negative particle"].nodes # radial position
time = 1000  # time in seconds
ax2.plot(r * 1e6, c(t=time, r=r), label="t={}[s]".format(time))

2021-11-19 15:29:22,931 - [WARNING] processed_variable.get_spatial_scale(520): No length scale set for negative particle. Using default of 1 [m]. 