Tip

An interactive online version of this notebook is available, which can be accessed via

---

Alternatively, you may download this notebook and run it offline.

Attention

You are viewing this notebook on the latest version of the documentation, where these notebooks may not be compatible with the stable release of PyBaMM since they can contain features that are not yet released. We recommend viewing these notebooks from the stable version of the documentation. To install the latest version of PyBaMM that is compatible with the latest notebooks, build PyBaMM from source.

Creating a Simple Model for SEI Growth

Before adding a new model, please read the contribution guidelines

In this notebook, we will run through the steps involved in creating a new model within pybamm. We will then solve and plot the outputs of the model. We have chosen to implement a very simple model of SEI growth. We first give a brief derivation of the model and discuss how to nondimensionalise the model so that we can show the full process of model conception to solution within a single notebook.

Note: if you run the entire notebook and then try to evaluate the earlier cells, you will likely receive an error. This is because the state of objects is mutated as it is passed through various stages of processing. In this case, we recommend that you restart the Kernel and then evaluate cells in turn through the notebook.

A Simple Model of Solid Electrolyte Interphase (SEI) Growth

The SEI is a porous layer that forms on the surfaces of negative electrode particles from the products of electrochemical reactions which consume lithium and electrolyte solvents. In the first few cycles of use, a lithium-ion battery loses a large amount of capacity; this is generally attributed to lithium being consumed to produce SEI. However, after a few cycles, the rate of capacity loss slows at a rate often (but not always) reported to scale with the square root of time. SEI growth is therefore often considered to be limited in some way by a diffusion process.

Dimensional Model

In our simple SEI model, we consider a one-dimensional SEI which extends from the surface of a planar negative electrode at \(x=0\) until \(x=L\), where \(L\) is the thickness of the SEI. Since the SEI is porous, there is some electrolyte within the region \(x\in[0, L]\) and therefore some concentration of solvent, \(c\). Within the porous SEI, the solvent is transported via a diffusion process according to:

\[\begin{split}\frac{\partial c}{\partial t} = - \frac{\partial N}{\partial x}, \quad N = - D(c) \frac{\partial c}{\partial x} \label{1}\\\end{split}\]

where \(t\) is the time, \(N\) is the solvent flux, and \(D(c)\) is the effective solvent diffusivity (a function of the solvent concentration).

On the electrode-SEI surface (\(x=0\)) the solvent is consumed by the SEI growth reaction, \(R\). We assume that diffusion of solvent in the bulk electrolyte (\(x>L\)) is fast so that on the SEI-electrolyte surface (\(x=L\)) the concentration of solvent is fixed at the value \(c_{\infty}\). Therefore, the boundary conditions are

\[N|_{x=0} = - R, \quad c|_{x=L} = c_{\infty},\]

We also assume that the concentration of solvent within the SEI is initially uniform and equal to the bulk electrolyte solvent concentration, so that the initial condition is

\[c|_{t=0} = c_{\infty}\]

Since the SEI is growing, we require an additional equation for the SEI thickness. The thickness of the SEI grows at a rate proportional to the SEI growth reaction \(R\), where the constant of proportionality is the partial molar volume of the reaction products, \(\hat{V}\). We also assume that the SEI is initially of thickness \(L_0\). Therefore, we have

\[\frac{d L}{d t} = \hat{V} R, \quad L|_{t=0} = L_0\]

Finally, we assume for the sake of simplicity that the SEI growth reaction is irreversible and that the potential difference across the SEI is constant. The reaction is also assumed to be proportional to the concentration of solvent at the electrode-SEI surface (\(x=0\)). Therefore, the reaction flux is given by

\[R = k c|_{x=0}\]

where \(k\) is the reaction rate constant (which is in general dependent upon the potential difference across the SEI).

Fixing the moving boundary

The model above is a moving boundary problem as it is defined in \(x\in[0, L]\). In order to solve it, we need to fix the boundary by introducing the scaling

\[x = L \xi.\]

Then, applying the chain rule we have

\[\frac{\partial}{\partial x} \rightarrow \frac{1}{L} \frac{\partial}{\partial \xi}, \quad \text{ and } \quad \frac{\partial}{\partial t} \rightarrow \frac{\partial}{\partial t} - \frac{L'(t)}{L(t)} \xi \frac{\partial}{\partial \xi}.\]

Then, (1) can be rewritten as

\[\frac{\partial c}{\partial t} = \frac{\hat{V} R}{L} \xi \frac{\partial c}{\partial \xi} + \frac{1}{L^2} \frac{\partial}{\partial \xi} \left( D(c) \frac{\partial c}{\partial \xi}\right)\]

Entering the Model into PyBaMM

As always, we begin by importing PyBaMM and changing our working directory to the root of the pybamm/ folder.

%pip install "pybamm[plot,cite]" -q    # install PyBaMM if it is not installed
import pybamm
import numpy as np
import os

os.chdir(pybamm.__path__[0] + "/..")

[notice] A new release of pip is available: 23.0.1 -> 23.2
[notice] To update, run: pip install --upgrade pip
Note: you may need to restart the kernel to use updated packages.

A model is defined in six steps: 1. Initialise model 2. Define parameters and variables 3. State governing equations 4. State boundary conditions 5. State initial conditions 6. State output variables

We shall proceed through each step to enter our simple SEI growth model.

Initialise model

We first initialise the model using the BaseModel class. This sets up the required structure for our model.

model = pybamm.BaseModel()

Define parameters and variables

In our SEI model, we have two dimensionless parameters, \(k\) and \(\hat{V}\), and one dimensionless function \(D(c)\), which are all given in terms of the dimensional parameters, see (5). In pybamm, inputs are dimensional, so we first state all the dimensional parameters. We then define the dimensionless parameters, which are expressed an non-dimensional groupings of dimensional parameters. To define the dimensional parameters, we use the Parameter object to create parameter symbols. Parameters which are functions are defined using FunctionParameter object and should be defined within a python function as shown.

# dimensional parameters
k = pybamm.Parameter("Reaction rate constant [m.s-1]")
L_0 = pybamm.Parameter("Initial thickness [m]")
V_hat = pybamm.Parameter("Partial molar volume [m3.mol-1]")
c_inf = pybamm.Parameter("Bulk electrolyte solvent concentration [mol.m-3]")


def D(cc):
    return pybamm.FunctionParameter(
        "Diffusivity [m2.s-1]", {"Solvent concentration [mol.m-3]": cc}
    )

We now define the dimensionless variables in our model. Since these are the variables we solve for directly, we do not need to write them in terms of the dimensional variables. We simply use SpatialVariable and Variable to create the required symbols:

xi = pybamm.SpatialVariable("xi", domain="SEI layer", coord_sys="cartesian")
c = pybamm.Variable("Solvent concentration [mol.m-3]", domain="SEI layer")
L = pybamm.Variable("SEI thickness [m]")

State governing equations

We can now use the symbols we have created for our parameters and variables to write out our governing equations. Note that before we use the reaction flux and solvent flux, we must derive new symbols for them from the defined parameter and variable symbols. Each governing equation must also be stated in the explicit form d/dt = rhs since pybamm only stores the right hand side (rhs) and assumes that the left hand side is the time derivative. The governing equations are then simply

# SEI reaction flux
R = k * pybamm.BoundaryValue(c, "left")

# solvent concentration equation
N = -1 / L * D(c) * pybamm.grad(c)
dcdt = (V_hat * R) / L * pybamm.inner(xi, pybamm.grad(c)) - 1 / L * pybamm.div(N)

# SEI thickness equation
dLdt = V_hat * R

Once we have stated the equations, we can add them to the model.rhs dictionary. This is a dictionary whose keys are the variables being solved for, and whose values correspond right hand sides of the governing equations for each variable.

model.rhs = {c: dcdt, L: dLdt}

State boundary conditions

We only have boundary conditions on the solvent concentration equation. We must state where a condition is Neumann (on the gradient) or Dirichlet (on the variable itself).

The boundary condition on the electrode-SEI (x=0) boundary is:

\[N|_{\xi=0} = - R, \quad N|_{\xi=0} = - D(c|_{\xi=0} )\left.\frac{\partial c}{\partial \xi}\right|_{\xi=0}\]

which is a Neumann condition. To implement this boundary condition in pybamm, we must first rearrange the equation so that the gradient of the concentration, \(\nabla c|_{x=0}\), is the subject. Therefore we have

\[\left.\frac{\partial c}{\partial \xi}\right|_{\xi=0} = \frac{R L}{D(c|_{\xi=0} )}\]

which we enter into pybamm as

D_left = pybamm.BoundaryValue(
    D(c), "left"
)  # pybamm requires BoundaryValue(D(c)) and not D(BoundaryValue(c))
grad_c_left = R * L / D_left

On the SEI-electrolyte boundary (\(\xi=1\)), we have the boundary condition

\[c|_{\xi=1} = c_∞\]

which is a Dirichlet condition and is just entered as

c_right = c_inf

We now load these boundary conditions into the model.boundary_conditions dictionary in the following way, being careful to state the type of boundary condition:

model.boundary_conditions = {
    c: {"left": (grad_c_left, "Neumann"), "right": (c_right, "Dirichlet")}
}

State initial conditions

There are two initial conditions in our model:

\[c|_{t=0} = c_∞, \quad L|_{t=0} = L_0\]

which are simply written in pybamm as

c_init = c_inf
L_init = L_0

and then included into the model.initial_conditions dictionary:

model.initial_conditions = {c: c_init, L: L_init}

State output variables

We already have everything required in model for the model to be used and solved, but we have not yet stated what we actually want to output from the model. PyBaMM allows users to output any combination of symbols as an output variable therefore allowing the user the flexibility to output important quantities without further tedious postprocessing steps.

Some useful outputs for this simple model are: - the SEI thickness - the SEI growth rate - the solvent concentration

These are added to the model by adding entries to the model.variables dictionary

model.variables = {
    "SEI thickness [m]": L,
    "SEI growth rate [m]": dLdt,
    "Solvent concentration [mol.m-3]": c,
}

The model is now fully defined and ready to be used. If you plan on reusing the model several times, you can additionally set model defaults which may include: a default geometry to run the model on, a default set of parameter values, a default solver, etc.

Using the Model

The model will now behave in the same way as any of the inbuilt PyBaMM models. However, to demonstrate that the model works we display the steps involved in solving the model but we will not go into details within this notebook.

# define geometry
geometry = pybamm.Geometry(
    {"SEI layer": {xi: {"min": pybamm.Scalar(0), "max": pybamm.Scalar(1)}}}
)


def Diffusivity(cc):
    return cc * 10 ** (-12)


# parameter values (not physically based, for example only!)
param = pybamm.ParameterValues(
    {
        "Reaction rate constant [m.s-1]": 1e-6,
        "Initial thickness [m]": 1e-6,
        "Partial molar volume [m3.mol-1]": 10,
        "Bulk electrolyte solvent concentration [mol.m-3]": 1,
        "Diffusivity [m2.s-1]": Diffusivity,
    }
)

# process model and geometry
param.process_model(model)
param.process_geometry(geometry)

# mesh and discretise
submesh_types = {"SEI layer": pybamm.Uniform1DSubMesh}
var_pts = {xi: 100}
mesh = pybamm.Mesh(geometry, submesh_types, var_pts)

spatial_methods = {"SEI layer": pybamm.FiniteVolume()}
disc = pybamm.Discretisation(mesh, spatial_methods)
disc.process_model(model)

<pybamm.models.base_model.BaseModel at 0x7f3a8005b490>

# solve
solver = pybamm.ScipySolver()
t = [0, 100]  # solve for 100s
solution = solver.solve(model, t)

# post-process output variables
L_out = solution["SEI thickness [m]"]
c_out = solution["Solvent concentration [mol.m-3]"]

Using these outputs, we can now plot the SEI thickness as a function of time and also the solvent concentration profile within the SEI. We use a slider to plot the concentration profile at different times. Note that, even though our model is written in nondimensional form, the processed variables are functions of dimensional space and time (in SI units).

import matplotlib.pyplot as plt

# plot SEI thickness in microns as a function of t in microseconds
# and concentration in mol/m3 as a function of x in microns
L_0_eval = param.evaluate(L_0)
xi = np.linspace(0, 1, 100)  # dimensionless space


def plot(t):
    _, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5))
    ax1.plot(solution.t, L_out(solution.t) * 1e6)
    ax1.plot(t, L_out(t) * 1e6, "r.")
    ax1.set_ylabel(r"SEI thickness [$\mu$m]")
    ax1.set_xlabel(r"t [s]")

    ax2.plot(xi * L_out(t) * 1e6, c_out(t, xi))
    ax2.set_ylim(0, 1.1)
    ax2.set_xlim(0, L_out(solution.t[-1]) * 1e6)
    ax2.set_ylabel("Solvent concentration [mol.m-3]")
    ax2.set_xlabel(r"x [$\mu$m]")

    plt.tight_layout()
    plt.show()


import ipywidgets as widgets

widgets.interact(
    plot, t=widgets.FloatSlider(min=0, max=solution.t[-1], step=0.1, value=0)
);

Formally adding your model

The purpose of this notebook has been to go through the steps involved in getting a simple model working within PyBaMM. However, if you plan on reusing your model and want greater flexibility then we recommend that you create a new class for your model. We have set out instructions on how to do this in the “Adding a Model” tutorial in the documentation.

References

The relevant papers for this notebook are:

pybamm.print_citations()

[1] Joel A. E. Andersson, Joris Gillis, Greg Horn, James B. Rawlings, and Moritz Diehl. CasADi – A software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 11(1):1–36, 2019. doi:10.1007/s12532-018-0139-4.
[2] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, and others. Array programming with NumPy. Nature, 585(7825):357–362, 2020. doi:10.1038/s41586-020-2649-2.
[3] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). Journal of Open Research Software, 9(1):14, 2021. doi:10.5334/jors.309.
[4] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, and others. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods, 17(3):261–272, 2020. doi:10.1038/s41592-019-0686-2.