Skip to main content

Libraries

PyBaMM Model Devel...

PDE models in... [pybamm]

Libraries

PyBaMM Model Devel...

Single Partic... [pybamm]

Single Particle Model
PyBaMM documentation by the PyBaMM Team

PyBaMM documentation by the PyBaMM Team

unknown

Single Particle Model

The Single Particle Model

Now we will extend our PDE model to the full single particle model. The single particle model is a system of PDEs that describes the behaviour of a lithium-ion battery electrode particle.

The Single Particle Model state equations

The SPM model only needs to solve for the concentration of lithium ions in the positive and negative electrodes, cnc_n and cpc_p. The diffusion of lithium ions in each electrode particle cic_i is given by:
∂ci∂t=∇⋅(Di∇ci)\frac{\partial c_i}{\partial t} = \nabla \cdot (D_i \nabla c_i)
subject to the following boundary and initial conditions:
∂ci∂r∣r=0=0,∂c∂r∣r=Ri=−ji,c∣t=0=ci0\left.\frac{\partial c_i}{\partial r}\right\vert_{r=0} = 0, \quad \left.\frac{\partial c}{\partial r}\right\vert_{r=R_i} = -j_i, \quad \left.c\right\vert_{t=0} = c^0_i
where cic_i is the concentration of lithium ions in the positive (i=ni=n) or negative (i=pi=p) electrode, DiD_i is the diffusion coefficient, jij_i is the interfacial current density, and ci0c^0_i is the concentration at the particle surface.
The primary addition here in comparison with our previous PDE model is the fluxes of lithium ions in the positive and negative electrodes jij_i. These are dependent on the applied current II:
jn=IanδnFA,jp=−IapδpFA,\begin{align*} j_n &= \frac{I}{a_n \delta_n F \mathcal{A}}, \qquad j_p &= \frac{-I}{a_p \delta_p F \mathcal{A}}, \end{align*}
where ai=3ϵi/Ria_i = 3 \epsilon_i / R_i is the specific surface area of the electrode, ϵi\epsilon_i is the volume fraction of active material, δi\delta_i is the thickness of the electrode, FF is the Faraday constant, and A\mathcal{A} is the electrode surface area.

Function Parameters in PyBaMM

The applied current I(t)I(t) is a parameter to the model, but unlike the other parameters we have seen so far, it is a function of time. We can define this using pybamm.FunctionParameter:
import pybamm
I = pybamm.FunctionParameter("Current function [A]", {"Time [s]": pybamm.t})
This will allow us to define the applied current as a function of time when we create our pybamm.ParameterValues object.
def current(t):
    return 1
param = pybamm.ParameterValues(
    {
        "Current function [A]": current
    }
)

Domains in PyBaMM

In the SPM model we have two different spatial domains, one for the positive electrode and one for the negative electrode. Remember that when we define a variable in PyBaMM, we can specify the domain using the domain keyword argument:
c_n = pybamm.Variable("Negative particle concentration [mol.m-3]", domain="negative particle")

SPM governing equations

Copy your PDE model from the previous challenge to a new file, and modify it to include the state equations for the concentration of lithium ions in both the positive and negative electrodes. Note that now you will have two variables, cnc_n and cpc_p, that must be defined on separate spatial domains, and a number of new parameters. Define the applied current II as a input parameter that is a function of time using pybamm.FunctionParameter.

Output variables for the Single Particle Model

Now that we have defined the equations to solve, we turn to the output variables that we need to calculate from the state variables cnc_n and cpc_p. The terminal voltage of the battery is given by:
V=Up(xps)−Un(xns)+ηp−ηnV = U_p(x_p^s) - U_n(x_n^s) + \eta_p - \eta_n
where UiU_i is the open circuit potential (OCP) of the electrode, xis=ci(r=Ri)/cimaxx_i^s = c_i(r=R_i) / c_i^{max} is the surface stoichiometry, and ηi\eta_i is the overpotential.
Assuming Butler-Volmer kinetics and αi=0.5\alpha_i = 0.5, the overpotential is given by:
ηi=2RTFsinh⁡−1(jiF2i0,i)\eta_i = \frac{2RT}{F} \sinh^{-1} \left( \frac{j_i F}{2i_{0,i}} \right)
where the exchange current density i0,ii_{0,i} is given by:
i0,i=kiFceci(r=Ri)cimax−ci(r=Ri)i_{0,i} = k_i F \sqrt{c_e} \sqrt{c_i(r=R_i)} \sqrt{c_i^{max} - c_i(r=R_i)}
where cec_e is the concentration of lithium ions in the electrolyte, and kik_i is the reaction rate constant.

Surface stoichiometry

Lets revisit the definition of the surface stoichiometry xisx_i^s, which is
xis=ci(r=Ri)cimaxx_i^s = \frac{c_i(r=R_i)}{c_i^{max}}
cimaxc_i^{max} is the maximum concentration of lithium ions in the electrode, and is a parameter of the model. However, ci(r=Ri)c_i(r=R_i) is the concentration of lithium ions at the surface of the electrode particle. How can we express this in PyBaMM, given that we only have the concentration cic_i defined on the whole domain?
To get the surface concentration, we can use the pybamm.boundary_value or pybamm.surf functions. The pybamm.boundary_value function returns the value of a variable at the boundary of a domain (either "left" or "right"), and the pybamm.surf function returns the value of a variable at the right boundary of a domain. For example, to get the surface concentration of cnc_n we can use:
c_n_surf = pybamm.surf(c_n)
or
c_n_surf = pybamm.boundary_value(c_n, "right")

Open Circuit Potentials (OCPs)

The OCPs UiU_i are functions of the surface stoichiometries xisx_i^s, and we can define them using pybamm.FunctionParameter in a similar way to the applied current II. For example, to define the OCP of the positive electrode as a function of the surface stoichiometry xpsx_p^s:
U_p = pybamm.FunctionParameter("Positive electrode OCP [V]", {"stoichiometry": "x_p_s"})

PyBaMM's built-in functions

PyBaMM has a number of built-in functions that can be used in expressions. For the SPM model, you will need to use the pybamm.sqrt and pybamm.arcsinh functions:
four = pybamm.Scalar(4)
two = pybamm.sqrt(four)
arcsinh_four = pybamm.arcsinh(four)
You can see a list of all the functions available in PyBaMM in the documentation.

SPM output variables

Define the output variables for the SPM model. You will need to define the overpotentials ηn\eta_n and ηp\eta_p, the exchange current densities i0,ni_{0,n} and i0,pi_{0,p}, and the terminal voltage VV. You will also need to define the OCPs UnU_n and UpU_p as functions of the surface stoichiometries xnsx_n^s and xpsx_p^s. You can use pybamm.FunctionParameter to define the OCPs as functions of the surface stoichiometries.
Define the following output variables for the model
  • Terminal voltage VV
  • Surface concentration in negative particle cnsc_n^s

Discretising and solving the Single Particle Model

Now that we have defined the SPM model, we can discretise and solve it using PyBaMM. We can use the same meshing and discretisation methods as in the previous section, but we now need to discretise the model on two different spatial domains.

Discretise and solve the SPM model

Discretise and solve the SPM model using the same methods as in the previous section. The following parameter values object copies the parameters from the PyBaMM Chen2020 model, feel free to use this to define the parameters for the SPM model.
param = pybamm.ParameterValues("Chen2020")
# PyBaMM parameters provide the exchange current density directly, rather than the reaction rate, so define here
param.update({
    "Positive electrode reaction rate [m.s-1]": 1e-3,
    "Negative electrode reaction rate [m.s-1]": 1e-3,
}, check_already_exists=False)