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.

A half cell model

In the previous notebook we saw how to compare full and reduced-order models. Both of these models were posed on a single domain: the negative electrode particle. Here we will see how to create a model which contains multiple domains.

We consider a problem posed on a half-cell geometry, which consists of a separator (\(-L_s<x<0\)) and a positive electrode (\(0<x<L_p\)). These two regions are considered “macroscale” domains. At each point in the positive electrode we treat a “microscale” problem to model diffusion of lithium within the positive particles (\(0<r<R_p\)). We will see how to create variables in each of the different domains so that the governing equations are properly coupled.

In the interest of simplicity we assume that the current in both the the solid and electrolyte is given by Ohm’s law, and ignore any concentration gradients in the electrolyte. The governing equations for charge conservation at the macroscale are then

\[\begin{split}i_e = -\kappa \nabla \phi_e, \quad \nabla i_e = a j, \quad -L_s<x<0, \\ i = -\sigma \nabla \phi, \quad \nabla \cdot i = -a j, \quad 0<x<L_p,\end{split}\]

where \(i\) and \(i_e\) are the current densities in the solid and electrolyte, respectively, \(\phi\) and \(\phi_e\) are the electric potentials in the solid and electrolyte, respectively, \(\sigma\) is the solid conductivity, \(\kappa\) is the ionic conductivity, \(a\) is the electrode surface area per unit volume and \(j\) the interfacial current density. The charge conservation equations are subject to the boundary conditions

\[\left.\phi_e\right\vert_{x=-L_s} = 0, \quad \left.i_e\right\vert_{x=L_p} = 0, \quad \left.i\right\vert_{x=0} = 0, \quad \left.i\right\vert_{x=L_p} = \frac{I_{\text{app}}}{A},\]

where \(I_{\text{app}}\) is the applied current and \(A\) is the electrode cross-sectional area. We then have an equation posed at each macroscopic point in the electrode (\(0<x<L_p\)) describing transport of lithium within the active material particles. That is,

\[\frac{\partial c}{\partial t} = \nabla \cdot (D \nabla c), \quad 0<r<R_p,\]

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, \(F\) Faraday’s constant, and \(c_0\) the initial concentration. For the interfacial current density we assume Butler-Volmer kinetics

\[\begin{split}j = \begin{cases} 0 &-L_s<x<0 \\ 2 j_0(c)\sinh\left(\frac{F}{2RT}(\phi-\phi_e-U(c))\right) &0<x<L_p \end{cases},\end{split}\]

where \(j_0\) is the exchange current density, and \(U\) is the open-circuit potential.

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

model = pybamm.BaseModel()

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

Let’s first define our model variables. We can define the electric potential in the positive electrode in the same way as we defined the concentration variables in the previous notebooks

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

The potential in the electrolyte spans two domains. To set this up we first define the electric potential in each of the domains

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

and then concatenate these two variables together to define a single variable that spans the separator and positive electrode

phi_e = pybamm.concatenation(phi_e_s, phi_e_p)

Note that in our formulation the separator will be on the left and the positive electrode will be on the right, so this is the order in which we concatenated the variables for the electrolyte potential in each domain.

The concentration in the electrode particles can vary both in \(r\) and \(x\), but diffusion only occurs in the \(r\) direction. In order to handle this situation we introduce the concept of “auxiliary domains”. These are domains in which quantities can vary, but spatial operators do not act. To set up our concentration variable we create a Variable which has domain “positive particle” and secondary domain “positive electrode”

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

Now spatial operators acting on c only act in the \(r\) direction (corresponding to the primary domain “positive particle”), but c can still depend on \(x\) (corresponding to the secondary domain “positive electrode”). For more details on the different domains (primary, secondary, etc.) see the broadcasts notebook.

Next we will define our parameters. As seen before, scalar parameters can be defined using the Parameter object

F = pybamm.Parameter("Faraday constant [C.mol-1]")
R = pybamm.Parameter("Molar gas constant [J.mol-1.K-1]")
T = pybamm.Parameter("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]")

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 [mol.m-3]")

Parameters can also have some functional dependence. We can define such parameters using the FunctionParameter object. We also need to specify the inputs (i.e. the variables on which the function depends). In our example both the exchange current density and the open-circuit potential will depend on the particle surface concentration.

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)

We also need to define the interfacial current, which is zero in the separator and given by the Butler-Volmer equation in the positive electrode

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)

Now we can write our governing equations, boundary and initial conditions. Note that we provide initial conditions for the algebraic equations. These are not really initial conditions, but are used as an initial guess for the solver.

# charge conservation equations
i = -sigma * pybamm.grad(phi)
i_e = -kappa * pybamm.grad(phi_e)
model.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.rhs = {c: dcdt}

# boundary conditions
model.boundary_conditions = {
    phi: {
        "left": (pybamm.Scalar(0), "Neumann"),
        "right": (-I_app / A / sigma, "Neumann"),
    },
    phi_e: {
        "left": (pybamm.Scalar(0), "Dirichlet"),
        "right": (pybamm.Scalar(0), "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.initial_conditions = {phi: U_init, phi_e: 0, c: c0}

Finally we can add any variables of interest to the model

model.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

As we have seen before, we can provide values for our parameters using the ParameterValues class. As well as providing scalar values, we also need to provide the functional form using by our FunctionParameter objects. Here we will define these functions locally, but we could provide the path to a function defined elsewhere or provide data (see the parameterization notebook).

from pybamm import tanh

# both functions will depend on the maximum concentration
c_max = pybamm.Parameter("Maximum concentration in positive electrode [mol.m-3]")


def exchange_current_density(c_surf):
    k = 6 * 10 ** (-7)  # reaction rate [(A/m2)(m3/mol)**1.5]
    c_e = 1000  # (constant) electrolyte concentration [mol.m-3]
    return k * c_e**0.5 * c_surf**0.5 * (c_max - c_surf) ** 0.5


def open_circuit_potential(c_surf):
    stretch = 1.062
    sto = stretch * c_surf / c_max

    u_eq = (
        2.16216
        + 0.07645 * tanh(30.834 - 54.4806 * sto)
        + 2.1581 * tanh(52.294 - 50.294 * sto)
        - 0.14169 * tanh(11.0923 - 19.8543 * sto)
        + 0.2051 * tanh(1.4684 - 5.4888 * sto)
        + 0.2531 * tanh((-sto + 0.56478) / 0.1316)
        - 0.02167 * tanh((sto - 0.525) / 0.006)
    )
    return u_eq

Now we can pass these functions, along with our scalar parameters, to ParameterValues

param = pybamm.ParameterValues(
    {
        "Surface area per unit volume [m-1]": 0.15e6,
        "Positive particle radius [m]": 10e-6,
        "Separator thickness [m]": 25e-6,
        "Positive electrode thickness [m]": 100e-6,
        "Electrode cross-sectional area [m2]": 2.8e-2,
        "Applied current [A]": 0.9,
        "Positive electrode conductivity [S.m-1]": 10,
        "Electrolyte conductivity [S.m-1]": 1,
        "Diffusion coefficient [m2.s-1]": 1e-13,
        "Faraday constant [C.mol-1]": 96485,
        "Initial concentration [mol.m-3]": 25370,
        "Molar gas constant [J.mol-1.K-1]": 8.314,
        "Temperature [K]": 298.15,
        "Maximum concentration in positive electrode [mol.m-3]": 51217,
        "Positive electrode exchange-current density [A.m-2]": exchange_current_density,
        "Positive electrode OCP [V]": open_circuit_potential,
    }
)

Next we define the geometry. In the same way that our variable for the concentration in the electrode particles had and “auxiliary domain”, our spatial variable \(r\) also has an auxiliary domain. This means that when the model in discretised there will be the correct number of particles included in the geometry - one for each point in \(x\).

r = pybamm.SpatialVariable(
    "r",
    domain=["positive particle"],
    auxiliary_domains={"secondary": "positive electrode"},
    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}},
}

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 = {
    "separator": pybamm.Uniform1DSubMesh,
    "positive electrode": pybamm.Uniform1DSubMesh,
    "positive particle": pybamm.Uniform1DSubMesh,
}
var_pts = {x_s: 10, x_p: 20, r: 30}
mesh = pybamm.Mesh(geometry, submesh_types, var_pts)

spatial_methods = {
    "separator": pybamm.FiniteVolume(),
    "positive electrode": pybamm.FiniteVolume(),
    "positive 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

We can now solve the model

# solve
solver = pybamm.CasadiSolver(root_tol=1e-3)
t_eval = np.linspace(0, 3600, 600)
solution = solver.solve(model, t_eval)

and plot the results. To make the plot we will use dynamic_plot which automatically creates a slider plot given a solution and a list of variables to plot. By nesting variables in the list we can plot two variables together on the same axes.

# plot
pybamm.dynamic_plot(
    solution,
    [
        "Positive electrode potential [V]",
        "Electrolyte potential [V]",
        "Positive electrode interfacial current density [A.m-2]",
        "Positive particle surface concentration [mol.m-3]",
        "Average positive particle surface concentration [mol.m-3]",
        ["Positive electrode OCP [V]", "Voltage [V]"],
    ],
)

2021-11-19 15:30:23,304 - [WARNING] processed_variable.get_spatial_scale(520): No length scale set for positive electrode. Using default of 1 [m].
2021-11-19 15:30:23,328 - [WARNING] processed_variable.get_spatial_scale(520): No length scale set for separator. Using default of 1 [m].
2021-11-19 15:30:23,367 - [WARNING] processed_variable.get_spatial_scale(520): No length scale set for positive electrode. Using default of 1 [m].
2021-11-19 15:30:23,395 - [WARNING] processed_variable.get_spatial_scale(520): No length scale set for positive electrode. Using default of 1 [m].

<pybamm.plotting.quick_plot.QuickPlot at 0x7ff07a82d130>

In the next notebook we consider a simple model for SEI growth, and show how to correctly pose the model in non-dimensional form and then create and solve it using pybamm.

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.