Meaning of the thermal outputs and implementation

[alexgUCSD]

Hi all,

I would like to compute the heat generation and SOC from models in PyBaMM (probably SPM and DFN) to compare to ECM outputs.

Just for some context, I have a working ECM model implemented in Python based on this paper "Battery Performance Modeling on Maxwell X-57" from Chin et al. (http://openmdao.org/pubs/chin_battery_performance_x57_2019.pdf). This model outputs the SOC as well as the overall heat generation.

The idea is to use a physics-based model such as the one implemented in PyBaMM instead of the ECM. In that case, for comparison, we would like to get the same outputs as in the ECM from the models in PyBaMM, i.e., overall heat generation and SOC. The cell temperature would also be great. The code I am using at the moment is below:

import pybamm
import numpy as np
import matplotlib.pyplot as plt

# inputs
options = {"thermal": "lumped"}
DFN_21700 = pybamm.lithium_ion.DFN(options=options)
geometry = pybamm.battery_geometry(form_factor="cylindrical")
chemistry=pybamm.parameter_sets.Mohtat2020
parameters = pybamm.ParameterValues(chemistry=chemistry)
parameters["Initial concentration in positive electrode [mol.m-3]"] = 38000.

# Experiment setup
experiment = pybamm.Experiment(
    ["Discharge at 1C for 2 minutes (1 second period)",
    "Rest for 20 minutes (1 second period)"]*10
)

# Solving the model
var_pts = {"x_n": 30, "x_s": 30, "x_p": 30, "r_n": 10, "r_p": 10}
pybamm.set_logging_level("NOTICE")
sim_DFN = pybamm.Simulation(DFN_21700, experiment=experiment, parameter_values=parameters, var_pts=var_pts)
sim_DFN.solve(initial_soc=1)

# Outputs
solution_DFN = sim_DFN.solution
t_DFN = solution_DFN["Time [s]"].entries

SOC_DFN_negative = solution_DFN["Negative electrode SOC"].entries
SOC_DFN_positive = solution_DFN["Positive electrode SOC"].entries

Total_heating_DFN = solution_DFN["Total heating [W.m-3]"].entries
Ohmic_heating_DFN = solution_DFN["Ohmic heating [W.m-3]"].entries
Cell_Temp_DFN = solution_DFN["Cell temperature [K]"].entries

# Plots
fig, axs = plt.subplots(1, 2)

axs[0].plot(t_DFN, 100*SOC_DFN_negative)
axs[0].set_title("Negative SOC as a function of time")
axs[0].set(xlabel="Time [in s]")
axs[0].set(ylabel="Negative SOC [in %]")

axs[1].plot(t_DFN, 100*SOC_DFN_positive)
axs[1].set_title("Positive SOC as a function of time")
axs[1].set(xlabel="Time [in s]")
axs[1].set(ylabel="Positive SOC [in %]")
plt.show()

Cell_volume = parameters["Cell volume [m3]"]
Total_heating_DFN = np.average(Total_heating_DFN, axis=0)*Cell_volume
Ohmic_heating_DFN = np.average(Ohmic_heating_DFN, axis=0)*Cell_volume
Cell_Temp_DFN = np.average(Cell_Temp_DFN, axis=0)

fig, axs = plt.subplots(1, 3)
axs[0].plot(t_DFN, Total_heating_DFN)
axs[0].set_title("Total heating as a function of time")
axs[0].set(xlabel="Time [in s]")
axs[0].set(ylabel="Total heating [W]")

axs[1].plot(t_DFN, Ohmic_heating_DFN)
axs[1].set_title("Joule heating as a function of time")
axs[1].set(xlabel="Time [in s]")
axs[1].set(ylabel="Joule heating [W]")

axs[2].plot(t_DFN, Cell_Temp_DFN)
axs[2].set_title("Average cell temperature as a function of time")
axs[2].set(xlabel="Time [in s]")
axs[2].set(ylabel="Cell temperature [K]")
plt.show()

As I am not very much experienced with battery modeling I am not too sure of how to interpret the "Negative electrode SOC" and "Positive electrode SOC". In addition, for the "Total heating" and "Cell temperature" is it a spatial distribution? How do I compute the overall heat generation of the cell at a given time-step? Should I sum the different entries and multiply by the volume?

Any explanation, help or comment regarding the code or the understanding of the model is more than welcome. I am very new to the topic so I am sorry if some questions seem basic.

Thank you a lot in advance for your help!
Alex

[TomTranter]

Hi @alexgUCSD so your question isn't basic and thanks for asking. The SoC of a battery is a little bit subjective because typically it will be defined between some voltage limits for safety and performance reasons. PyBaMM simulates the battery as a combination of the two electrodes and the voltage is determined by the combination of the states of each of them. In your code above you are pulling out variables defined for each electrode which are actually the average amount of lithium in the particles of each electrode divided by the absolute maximum that those particles could take so this is independent of the overall cell voltage limits you set. Separately the experiment class allows you to set the cell SoC and this does take into account voltage limits and will adjust the initial concentrations in each electrode based on the physical parameters of each electrode. Behind the scenes the simulation class is actually running a little experiment to do the cell balancing to work out what the concs should be and you can do this yourself to get a better idea of what's going on. There is a notebook related to simulating long experiments and examining state of health which includes the code you need to run explicitly here. Notice that the cell SoC rarely maps to the electrode SoC 1:1 as this would involve completely removing and filling up the lithium in the particles and would result in very low and high voltages. On the thermal question the code you have already looks good to me. Total heating should include the ohmic component. You can do model.variables.search('heating') to see all the components

[TomTranter]

I modified your code to calculate cell_SoC

import pybamm
import numpy as np
import matplotlib.pyplot as plt

# inputs
options = {"thermal": "lumped"}
DFN_21700 = pybamm.lithium_ion.DFN(options=options)
geometry = pybamm.battery_geometry(form_factor="cylindrical")
chemistry=pybamm.parameter_sets.Mohtat2020
parameters = pybamm.ParameterValues(chemistry=chemistry)
parameters["Initial concentration in positive electrode [mol.m-3]"] = 38000.

# SoC
# Calculate stoichiometries at 100% SOC
esoh_model = pybamm.lithium_ion.ElectrodeSOH()
esoh_sim = pybamm.Simulation(esoh_model, parameter_values=parameters)
param = DFN_21700.param

Vmin = parameters['Lower voltage cut-off [V]']
Vmax = parameters['Upper voltage cut-off [V]']
Cn = parameters.evaluate(param.n.cap_init)
Cp = parameters.evaluate(param.p.cap_init)
n_Li_init = parameters.evaluate(param.n_Li_particles_init)

esoh_sol = esoh_sim.solve(
    [0],
    inputs={"V_min": Vmin, "V_max": Vmax, "C_n": Cn, "C_p": Cp, "n_Li": n_Li_init},
)
print(f"Initial negative electrode SOC: {esoh_sol['x_100'].data[0]:.3f}")
print(f"Initial positive electrode SOC: {esoh_sol['y_100'].data[0]:.3f}")

# Update parameter values with initial conditions
c_n_max = parameters.evaluate(param.n.c_max)
c_p_max = parameters.evaluate(param.p.c_max)
parameters.update(
    {
        "Initial concentration in negative electrode [mol.m-3]": esoh_sol["x_100"].data[0] * c_n_max,
        "Initial concentration in positive electrode [mol.m-3]": esoh_sol["y_100"].data[0] * c_p_max,
    }
)


# Experiment setup
experiment = pybamm.Experiment(
    ["Discharge at 1C for 2 minutes (1 second period)",
    "Rest for 20 minutes (1 second period)"]*10
)

# Solving the model
var_pts = {"x_n": 30, "x_s": 30, "x_p": 30, "r_n": 10, "r_p": 10}
pybamm.set_logging_level("NOTICE")
sim_DFN = pybamm.Simulation(DFN_21700, experiment=experiment, parameter_values=parameters, var_pts=var_pts)
sim_DFN.solve(initial_soc=1)

# Outputs
solution_DFN = sim_DFN.solution
t_DFN = solution_DFN["Time [s]"].entries

SOC_DFN_negative = solution_DFN["Negative electrode SOC"].entries
SOC_DFN_positive = solution_DFN["Positive electrode SOC"].entries


x_100 = esoh_sol['x_100'].data[0]
x_0 = esoh_sol['x_0'].data[0]
x = solution_DFN['Average negative particle concentration'].entries
cell_SoC_x = 100*(x-x_0)/(x_100-x_0)

y_100 = esoh_sol['y_100'].data[0]
y_0 = esoh_sol['y_0'].data[0]
y = solution_DFN['Average positive particle concentration'].entries
cell_SoC_y = 100*(y-y_0)/(y_100-y_0)

Total_heating_DFN = solution_DFN["Total heating [W.m-3]"].entries
Ohmic_heating_DFN = solution_DFN["Ohmic heating [W.m-3]"].entries
Cell_Temp_DFN = solution_DFN["Cell temperature [K]"].entries

print('SoC from x and y should match:', np.allclose(cell_SoC_x, cell_SoC_y))

# Plots
fig, axs = plt.subplots(1, 2)

axs[0].plot(t_DFN, cell_SoC_x)
axs[0].set_title("Negative SOC as a function of time")
axs[0].set(xlabel="Time [in s]")
axs[0].set(ylabel="Negative SOC [in %]")

axs[1].plot(t_DFN, cell_SoC_y)
axs[1].set_title("Positive SOC as a function of time")
axs[1].set(xlabel="Time [in s]")
axs[1].set(ylabel="Positive SOC [in %]")
plt.show()

Cell_volume = parameters["Cell volume [m3]"]
Total_heating_DFN = np.average(Total_heating_DFN, axis=0)*Cell_volume
Ohmic_heating_DFN = np.average(Ohmic_heating_DFN, axis=0)*Cell_volume
Cell_Temp_DFN = np.average(Cell_Temp_DFN, axis=0)

fig, axs = plt.subplots(1, 3)
axs[0].plot(t_DFN, Total_heating_DFN)
axs[0].set_title("Total heating as a function of time")
axs[0].set(xlabel="Time [in s]")
axs[0].set(ylabel="Total heating [W]")

axs[1].plot(t_DFN, Ohmic_heating_DFN)
axs[1].set_title("Joule heating as a function of time")
axs[1].set(xlabel="Time [in s]")
axs[1].set(ylabel="Joule heating [W]")

axs[2].plot(t_DFN, Cell_Temp_DFN)
axs[2].set_title("Average cell temperature as a function of time")
axs[2].set(xlabel="Time [in s]")
axs[2].set(ylabel="Cell temperature [K]")
plt.show()

[TomTranter]

@tinosulzer @rtimms could this have been done in any more succinct way?

[rtimms]

This works with 23.4.1

import pybamm
import numpy as np
import matplotlib.pyplot as plt

# Set up model and parameters
options = {"thermal": "lumped"}
DFN_21700 = pybamm.lithium_ion.DFN(options=options)
param = DFN_21700.param
parameters = pybamm.ParameterValues("Mohtat2020")
parameters["Initial concentration in positive electrode [mol.m-3]"] = 38000.0

# Calculate stoichiometries at 100% SOC
x_0, x_100, y_100, y_0 = pybamm.lithium_ion.get_min_max_stoichiometries(
    parameters, param
)

# Update parameter values with initial conditions
c_n_max = parameters["Maximum concentration in negative electrode [mol.m-3]"]
c_p_max = parameters["Maximum concentration in positive electrode [mol.m-3]"]
parameters.update(
    {
        "Initial concentration in negative electrode [mol.m-3]": x_100 * c_n_max,
        "Initial concentration in positive electrode [mol.m-3]": y_100 * c_p_max,
    }
)

# Experiment setup
experiment = pybamm.Experiment(
    [
        "Discharge at 1C for 2 minutes (1 second period)",
        "Rest for 20 minutes (1 second period)",
    ]
    * 10
)

# Solving the model
var_pts = {"x_n": 30, "x_s": 30, "x_p": 30, "r_n": 10, "r_p": 10}
pybamm.set_logging_level("NOTICE")
sim_DFN = pybamm.Simulation(
    DFN_21700, experiment=experiment, parameter_values=parameters, var_pts=var_pts
)
sim_DFN.solve()

# Outputs
solution_DFN = sim_DFN.solution
t_DFN = solution_DFN["Time [s]"].entries

x = solution_DFN["Average negative particle concentration"].entries
cell_SoC_x = 100 * (x - x_0) / (x_100 - x_0)

y = solution_DFN["Average positive particle concentration"].entries
cell_SoC_y = 100 * (y - y_0) / (y_100 - y_0)

Total_heating_DFN = solution_DFN["Total heating [W.m-3]"].entries
Ohmic_heating_DFN = solution_DFN["Ohmic heating [W.m-3]"].entries
Cell_Temp_DFN = solution_DFN["Cell temperature [K]"].entries

print("SoC from x and y should match:", np.allclose(cell_SoC_x, cell_SoC_y))

# Plots
fig, axs = plt.subplots(1, 2)

axs[0].plot(t_DFN, cell_SoC_x)
axs[0].set_title("Negative SOC as a function of time")
axs[0].set(xlabel="Time [in s]")
axs[0].set(ylabel="Negative SOC [in %]")

axs[1].plot(t_DFN, cell_SoC_y)
axs[1].set_title("Positive SOC as a function of time")
axs[1].set(xlabel="Time [in s]")
axs[1].set(ylabel="Positive SOC [in %]")

Cell_volume = parameters["Cell volume [m3]"]
Total_heating_DFN = np.average(Total_heating_DFN, axis=0) * Cell_volume
Ohmic_heating_DFN = np.average(Ohmic_heating_DFN, axis=0) * Cell_volume
Cell_Temp_DFN = np.average(Cell_Temp_DFN, axis=0)

fig, axs = plt.subplots(1, 3)
axs[0].plot(t_DFN, Total_heating_DFN)
axs[0].set_title("Total heating as a function of time")
axs[0].set(xlabel="Time [in s]")
axs[0].set(ylabel="Total heating [W]")

axs[1].plot(t_DFN, Ohmic_heating_DFN)
axs[1].set_title("Joule heating as a function of time")
axs[1].set(xlabel="Time [in s]")
axs[1].set(ylabel="Joule heating [W]")

axs[2].plot(t_DFN, Cell_Temp_DFN)
axs[2].set_title("Average cell temperature as a function of time")
axs[2].set(xlabel="Time [in s]")
axs[2].set(ylabel="Cell temperature [K]")
plt.show()

[Alistair-Douglas]

I just want to say the help and commitment you guys go to creating this community is fantastic.

To make an effort to update a reply almost a year later is just so good.

I'm using PyBamm for my dissertation to evaluate the value of using thermal management to increase the profit of FCAS contracts (UK) and this piece of software is making it possible for someone who is a general novice to both Python and battery chemistry.

Thank you so much to the whole team.

[alexgUCSD]

Hi @TomTranter,
First of all, thank you so much for the detailed answer. I think I now understand the SoC part. However, when I am trying to run the suggested code, I have the following error Cn = parameters.evaluate(param.n.cap_init) AttributeError: 'LithiumIonParameters' object has no attribute 'n'. It does seem that your code should be working based on the notebook you shared so I am not sure why it is not running.
In addition, for the thermal part, I am just hesitating with this line of code Total_heating_DFN = np.average(Total_heating_DFN, axis=0)*Cell_volume. It seems that the heat is computed at each discretization point and I would like the global quantity at each time-step in Watt. Hence, the average and the volume. Is that the appropriate way to obtain this value or is there another option that I am not seeing?
Again, thanks a lot for your answer!

[TomTranter]

what version of pybamm are you running? the n.cap_init is new syntax I think. For the heating we already provide averaged heating outputs. If you search heating in the variables you should see the ones you need

[alexgUCSD]

Hi,
I was running pybamm 22.3 and I just updated to 22.6 and it is now running! I have also found the variables I was looking for.
Thank you for all the help!
Alex