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.

Comparison of PyBaMM and COMSOL Discharge Curves#

In this notebook we compare the discharge curves obtained by solving the DFN model both in PyBaMM and COMSOL. Results are presented for a range of C-rates, and we see an excellent agreement between the two implementations. If you would like to compare internal variables please see the script compare_comsol_DFN which creates a slider plot comparing potentials and concentrations as functions of time and space for a given C-rate. For more information on the DFN model, see the DFN notebook.

First we need to import pybamm, and then change 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
import pickle
import matplotlib.pyplot as plt

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

We then create a dictionary of the C-rates we would like to solve for and compare. Note that the repository currently only contains COMSOL results for the C-rates listed below.

C_rates = {"01": 0.1, "05": 0.5, "1": 1, "2": 2, "3": 3}

We get the DFN model equations, geometry, and default parameters. Before processing the model, we adjust the electrode height and depth to be 1 m, and also adjust the electrode conductivities, to match the one-dimensional model we solved in COMSOL. The model is then processed using the default geometry and updated parameters. Finally, we create a mesh and discretise the model.

# load model and geometry
model = pybamm.lithium_ion.DFN()
geometry = model.default_geometry

# load parameters and process model and geometry
param = model.default_parameter_values
        "Electrode width [m]": 1,
        "Electrode height [m]": 1,
        "Negative electrode conductivity [S.m-1]": 126,
        "Positive electrode conductivity [S.m-1]": 16.6,
        "Current function [A]": "[input]",

# create mesh
var_pts = {"x_n": 31, "x_s": 11, "x_p": 31, "r_n": 11, "r_p": 11}
mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts)

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

We create the figure by looping over the dictionary of C-rates. In each step of the loop we load the COMSOL results from a .csv file and solve the DFN model in pybamm. The output variables are then processed, allowing us to plot the discharges curve computed using pybamm and COMSOL, and their absolute difference.

# create figure
fig, ax = plt.subplots(figsize=(15, 8))
discharge_curve = plt.subplot(211)
plt.xlim([0, 26])
plt.ylim([3.2, 3.9])
plt.xlabel(r"Discharge Capacity (Ah)")
plt.ylabel("Voltage (V)")
plt.title(r"Comsol $\cdots$ PyBaMM $-$")
voltage_difference_plot = plt.subplot(212)
plt.xlim([0, 26])
plt.xlabel(r"Discharge Capacity (Ah)")
plt.ylabel(r"$\vert V - V_{comsol} \vert$")
colors = iter(plt.cycler(color="bgrcmyk"))

# loop over C_rates dict to create plot
for key, C_rate in C_rates.items():
    # load the comsol results
    comsol_results_path = pybamm.get_parameters_filepath(
    comsol_variables = pickle.load(open(comsol_results_path, "rb"))
    comsol_time = comsol_variables["time"]
    comsol_voltage = comsol_variables["voltage"]

    # update current density
    current = 24 * C_rate

    # solve model at comsol times
    solver = pybamm.CasadiSolver(mode="fast")
    solution = solver.solve(
        model, comsol_time, inputs={"Current function [A]": current}
    time_in_seconds = solution["Time [s]"].entries
    # discharge capacity
    discharge_capacity = solution["Discharge capacity [A.h]"]
    discharge_capacity_sol = discharge_capacity(time_in_seconds)
    comsol_discharge_capacity = comsol_time * current / 3600

    # extract the voltage
    voltage = solution["Voltage [V]"]
    voltage_sol = voltage(time_in_seconds)

    # calculate the difference between the two solution methods
    end_index = min(len(time_in_seconds), len(comsol_time))
    voltage_difference = np.abs(voltage_sol[0:end_index] - comsol_voltage[0:end_index])

    # plot discharge curves and absolute voltage_difference
    color = next(colors)["color"]
        comsol_discharge_capacity, comsol_voltage, color=color, linestyle=":"
        label=f"{C_rate} C",
        discharge_capacity_sol[0:end_index], voltage_difference, color=color

    top=0.92, bottom=0.08, left=0.10, right=0.95, hspace=0.25, wspace=0.35
/var/folders/z4/5lmf5d5d23sc2gkhs__zfnfc0000gn/T/ipykernel_2839/4153409347.py:5: MatplotlibDeprecationWarning: Auto-removal of overlapping axes is deprecated since 3.6 and will be removed two minor releases later; explicitly call ax.remove() as needed.
  discharge_curve = plt.subplot(211)


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). Journal of Open Research Software, 9(1):14, 2021. doi:10.5334/jors.309.

[ ]: