An interactive online version of this notebook is available, which can be accessed via Open this notebook in Google Colab

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

Doyle-Fuller-Newman Model (DFN)#

Model Equations#

The DFN comprises equations for charge and mass conservation in the solid the solid and electrolyte, and also prescribes behaviour for the electrochemical reactions occurring on the interface between the solid an electrolyte. For more information please see [2] or other standard texts.

Below we summarise the model, with all parameters give in the table at the end of this notebook. Here we use a roman subscript \(\text{k} \in \text{n, s, p}\) is used to denote the regions negative electrode, separator, and positive electrode, respectively.

The model equations for the DFN read:

Charge conservation:#

\[\begin{split}\frac{\partial i_{\text{e,k}}}{\partial x} = \begin{cases} a_{\text{k}}j_{\text{k}}, \quad &\text{k} = \text{n, p}\\ 0, \qquad &\text{k} = \text{s} \end{cases},\end{split}\]
\[\begin{split}i_{\text{e,k}} = \epsilon_{\text{k}}^{\text{b}} \kappa_{\text{e}}(c_{\text{e,k}}) \left( - \frac{\partial \phi_{\text{e,k}}}{\partial x} + 2(1-t^+)\frac{RT}{F}\frac{\partial}{\partial x}\left(\log(c_{\text{e,k}})\right)\right), \text{k} \in \text{n, s, p}, \\ I-i_{\text{e,k}} = - \sigma_{\text{k}} \frac{\partial \phi_{\text{e,k}}}{\partial x}, \text{k} \in \text{n, s, p}.\end{split}\]

Mass conservation:#

\[\begin{split}\epsilon_{\text{k}} \frac{\partial c_{\text{e,k}}}{\partial t} = -\frac{\partial N_{\text{e,k}}}{\partial x} + \frac{1}{F}\frac{\partial i_{\text{e,k}}}{\partial x}, \text{k} \in \text{n, s, p},\\ N_{\text{e,k}} = -\epsilon_{\text{k}}^{\text{b}} D_{\text{e}}(c_{\text{e,k}}) \frac{\partial c_{\text{e,k}}}{\partial x} + \frac{t^+}{F} i_{\text{e,k}}, \\ \text{k} \in \text{n, s, p}, \\ \frac{\partial c_{\text{s,k}}}{\partial t} = -\frac{1}{r_{\text{k}}^2} \frac{\partial}{\partial r_{\text{k}}} \left(r_{\text{k}}^2 N_{\text{s,k}}\right), \\ \text{k} \in \text{n, p},\\ N_{\text{s,k}} = -D_{\text{s,k}}(c_{\text{s,k}}) \frac{\partial c_{\text{s,k}}}{\partial r_{\text{k}}}, \\ \text{k} \in \text{n, p}.\end{split}\]

Electrochemical reactions:#

\[\begin{split}j_{\text{k}} = 2 j_{\text{0,k}} \sinh\left(\frac{ F\eta_{\text{k}}}{2RT} \right), \\ \text{k} \in \text{n, p}, \\ j_{\text{0,k}} = c_{\text{s,k}}^{1/2} (1-c_{\text{s,k}})^{1/2}c_{\text{e,k}}^{1/2}\big|_{r_{\text{k}}=1}, \\ \text{k} \in \text{n, p}, \\ \eta_{\text{k}} = \phi_{\text{s,k}} - \phi_{\text{e,k}} - U_{\text{k}}(c_{\text{s,k}}\big|_{r_{\text{k}}=1}), \\ \text{k} \in \text{n, p}.\end{split}\]

These are to be solved subject to the following boundary conditions:


\[\begin{split}i_{\text{e,n}}\big|_{x=0} = 0, \quad i_{\text{e,p}}\big|_{x=L}=0, \\ \phi_{\text{e,n}}\big|_{x=L_{\text{n}}} = \phi_{\text{e,s}}\big|_{x=L_{\text{n}}}, \quad i_{\text{e,n}}\big|_{x=L_{\text{n}}} = i_{\text{e,s}}\big\vert_{x=L_{\text{n}}} = I, \\ \phi_{\text{e,s}}\big|_{x=L-L_{\text{p}}} = \phi_{\text{e,p}}\big|_{x=L-L_{\text{p}}}, \quad i_{\text{e,s}}\big|_{x=L-L_{\text{p}}} = i_{\text{e,p}}\big|_{x=L-L_{\text{p}}} = I.\end{split}\]

Concentration in the electrolyte:#

\[\begin{split}N_{\text{e,n}}\big|_{x=0} = 0, \quad N_{\text{e,p}}\big|_{x=L}=0,\\ c_{\text{e,n}}\big|_{x=L_{\text{n}}} = c_{\text{e,s}}|_{x=L_{\text{n}}}, \quad N_{\text{e,n}}\big|_{x=L_{\text{n}}}=N_{\text{e,s}}\big|_{x=L_{\text{n}}}, \\ c_{\text{e,s}}|_{x=L-L_{\text{p}}}=c_{\text{e,p}}|_{x=L-L_{\text{p}}}, \quad N_{\text{e,s}}\big|_{x=L-L_{\text{p}}}=N_{\text{e,p}}\big|_{x=L-L_{\text{p}}}.\end{split}\]

Concentration in the electrode active material:#

\[N_{\text{s,k}}\big|_{r_{\text{k}}=0} = 0, \quad \text{k} \in \text{n, p}, \quad \ \ N_{\text{s,k}}\big|_{r_{\text{k}}=R_{\text{k}}} = \frac{j_{\text{k}}}{F}, \quad \text{k} \in \text{n, p}.\]

Reference potential:#

\[\phi_{\text{s,cn}} = 0, \quad \boldsymbol{x} \in \partial \Omega_{\text{tab,n}}.\]

And the initial conditions:#

\[\begin{split}c_{\text{s,k}}(x,r,0) = c_{\text{s,k,0}}, \quad \text{k} \in \text{n, p},\\ c_{\text{e,k}}(x,0) = c_{\text{e,0}}, \quad \text{k} \in \text{n, s, p}.\end{split}\]

Example solving DFN using PyBaMM#

Below we show how to solve the DFN model, using the default geometry, mesh, parameters, discretisation and solver provided with PyBaMM. For a more detailed example, see the notebook on the SPM.

In order to show off all the different points at which the process of setting up and solving a model in PyBaMM can be customised we explicitly handle the stages of choosing a geometry, setting parameters, discretising the model and solving the model. However, it is often simpler in practice to use the Simulation class, which handles many of the stages automatically, as shown here.

First we need to import pybamm, along with numpy which we will use in this notebook.

%pip install "pybamm[plot,cite]" -q    # install PyBaMM if it is not installed
import pybamm
import numpy as np
WARNING: You are using pip version 21.0.1; however, version 21.1.2 is available.
You should consider upgrading via the '/home/user/Documents/PyBaMM/env/bin/python3.8 -m pip install --upgrade pip' command.
Note: you may need to restart the kernel to use updated packages.

We then load the DFN model and default geometry, and process them both using the default parameters.

# load model
model = pybamm.lithium_ion.DFN()

# create geometry
geometry = model.default_geometry

# load parameter values and process model and geometry
param = model.default_parameter_values

The next step is to set the mesh and discretise the model. Again, we choose the default settings.

# set mesh
mesh = pybamm.Mesh(geometry, model.default_submesh_types, model.default_var_pts)

# discretise model
disc = pybamm.Discretisation(mesh, model.default_spatial_methods)

The model is now ready to be solved. We select the default DAE solver for the DFN. Note that in order to successfully solve the system of DAEs we are required to give consistent initial conditions. This is handled automatically by PyBaMM during the solve operation.

# solve model
solver = model.default_solver
t_eval = np.linspace(0, 3600, 300)  # time in seconds
solution = solver.solve(model, t_eval)

To get a quick overview of the model outputs we can use the QuickPlot class, which plots a common set of useful outputs. The method Quickplot.dynamic_plot makes a slider widget.

quick_plot = pybamm.QuickPlot(
    solution, ["Positive electrode interfacial current density [A.m-2]"]


The relevant papers for this notebook are:

[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] Marc Doyle, Thomas F. Fuller, and John Newman. Modeling of galvanostatic charge and discharge of the lithium/polymer/insertion cell. Journal of the Electrochemical society, 140(6):1526–1533, 1993. doi:10.1149/1.2221597.
[3] 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.
[4] Scott G. Marquis, Valentin Sulzer, Robert Timms, Colin P. Please, and S. Jon Chapman. An asymptotic derivation of a single particle model with electrolyte. Journal of The Electrochemical Society, 166(15):A3693–A3706, 2019. doi:10.1149/2.0341915jes.
[5] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). ECSarXiv. February, 2020. doi:10.1149/osf.io/67ckj.