Thermal coupling¶

Breathe Simulate uses a lumped thermal model by default. The cell self heats under load and exchanges heat with the ambient. Sometimes you already have a richer (for example 3D) thermal model of your own and want Breathe to run the electrochemistry under your cell temperature trajectory.

This is open loop coupling:

1. Your external thermal model produces a cell temperature trajectory T(t).
2. You feed T(t) into Breathe Simulate, which runs the electrochemistry under it and returns the heat generation Q(t) and voltage.
3. You feed Q(t) back into your thermal model to advance it.

In this mode Breathe's own thermal model is overridden, so the ambient temperature and heat transfer coefficient are unused.

This notebook covers:

1. Querying the design parameters of the Molicel P45B.
2. Running an isothermal simulation where the temperature is held constant.
3. Using a temperature timeseries object to couple an external thermal model.

In [17]:

from breathe_simulate import api_interface as api
from breathe_simulate import Cycler, TemperatureProfile
from breathe_simulate import enable_notebook_plotly

enable_notebook_plotly()

from breathe_simulate import api_interface as api from breathe_simulate import Cycler, TemperatureProfile from breathe_simulate import enable_notebook_plotly enable_notebook_plotly()

Query the design parameters¶

We can pull the design parameters for any base battery. Here we use the Ni65 (CATL-147) D/MP2, and also read its equilibrium KPIs to get the nominal capacity, which we use to set up the cycler in C rate.

In [18]:

base_params = api.get_design_parameters("Molicel P45B")
base_params

base_params = api.get_design_parameters("Molicel P45B") base_params

Out[18]:

{'NLiRatio': 1.192712619876931,
 'NPratio': 1.039106051327928,
 'Vmax_V': 4.200003037059671,
 'Vmin_V': 2.4999805201342147,
 'aluminiumThickness_um': 22.2,
 'anode': 'molicel_p45b_anode',
 'anodeArealLoading_mAhcm2': 3.8996207810294305,
 'anodeCathodeAreaRatio': 1.065290648694187,
 'anodeMassLoading_mgcm2': 8.44896611475494,
 'anodePorosity': 0.16663641853607114,
 'anodeTabResistance_uOhmAh': 0.0,
 'anodeThickness_um': 48.90000000000002,
 'cathode': 'molicel_p45b_cathode',
 'cathodeArealLoading_mAhcm2': 3.7528611983790308,
 'cathodeMassLoading_mgcm2': 14.188761420434503,
 'cathodePorosity': 0.14216134985845785,
 'cathodeTabResistance_uOhmAh': 0.0,
 'cathodeThickness_um': 37.48333333333334,
 'concentration_molL': 1.0,
 'copperThickness_um': 11.2,
 'electrodeArea_cm2': 1495.6199999999994,
 'electrolyte': 'molicel_p45b_electrolyte',
 'electrolyteBuffer_rel': 0.2,
 'format': 'molicel_p45b',
 'lamne': 0.0,
 'lampe': 0.0,
 'lli': 0.0,
 'separatorThickness_um': 16.8}

In [19]:

results = api.get_eqm_kpis("Molicel P45B")
baseline_capacity = results.capacity
print(f"Nominal capacity: {baseline_capacity:.3f} Ah")

results = api.get_eqm_kpis("Molicel P45B") baseline_capacity = results.capacity print(f"Nominal capacity: {baseline_capacity:.3f} Ah")

Nominal capacity: 4.500 Ah

An isothermal simulation¶

First, a baseline. We set up a 1C discharge and run it isothermally, so the cell temperature is held constant at the initial temperature. In this mode there is no self heating and no exchange with the ambient, so the ambient temperature and heat transfer coefficient have no effect.

In [20]:

cycler = Cycler(selected_unit="C", cell_capacity=baseline_capacity)
cycler_dict = cycler.cc_dch(
    -1.0, 2.5
)  # 1C discharge to 2.5 V (discharge current is negative)
print(cycler_dict)

cycler = Cycler(selected_unit="C", cell_capacity=baseline_capacity) cycler_dict = cycler.cc_dch( -1.0, 2.5 ) # 1C discharge to 2.5 V (discharge current is negative) print(cycler_dict)

{'cycle_type': 'CC_DCH', 'selected_unit': 'C', 'control_parameters': {'I_dch': -1.0, 'V_min': 2.5}}

In [21]:

isothermal_result = api.run_sim(
    base_battery="Molicel P45B",
    cycler=cycler_dict,
    initialSoC=1.0,
    initialTemperature_degC=25.0,
    isothermal=True,
)

isothermal_result = api.run_sim( base_battery="Molicel P45B", cycler=cycler_dict, initialSoC=1.0, initialTemperature_degC=25.0, isothermal=True, )

The cell temperature stays flat at 25 °C for the whole run:

In [22]:

isothermal_result.plot_dynamic_response("Cell temperature [°C]")

isothermal_result.plot_dynamic_response("Cell temperature [°C]")

In [23]:

isothermal_result.plot_voltage_response()

isothermal_result.plot_voltage_response()

Couple an external thermal model¶

Instead of holding the temperature constant (isothermal) or letting Breathe's lumped model evolve it, we can prescribe the cell temperature as a function of time. This is how you couple an external thermal model.

We describe the trajectory with a TemperatureProfile, a temperature timeseries object. It takes time breakpoints (in seconds) and the corresponding cell temperatures (in °C). In a real workflow these values come from your own thermal model. Here we use an illustrative profile that ramps up and cools back down.

In [24]:

profile = TemperatureProfile(
    time=[0, 600, 1200, 1800, 2400, 3000, 3600],  # seconds
    value=[25, 28, 33, 38, 36, 31, 29],  # °C, e.g. from your 3D thermal model
)
profile.plot()

profile = TemperatureProfile( time=[0, 600, 1200, 1800, 2400, 3000, 3600], # seconds value=[25, 28, 33, 38, 36, 31, 29], # °C, e.g. from your 3D thermal model ) profile.plot()

Now pass the profile to run_sim via cellTemperatureProfile. It works with any cycler. The protocol (here a 1C discharge) is unchanged, and only the thermal behaviour differs.

In [25]:

coupled_result = api.run_sim(
    base_battery="Molicel P45B",
    cycler=cycler_dict,
    initialSoC=1.0,
    cellTemperatureProfile=profile,
)

coupled_result = api.run_sim( base_battery="Molicel P45B", cycler=cycler_dict, initialSoC=1.0, cellTemperatureProfile=profile, )

The simulated cell temperature now follows the trajectory we prescribed, interpolated onto the solver time points:

In [26]:

coupled_result.plot_dynamic_response("Cell temperature [°C]")

coupled_result.plot_dynamic_response("Cell temperature [°C]")

And the heat generation is returned in the dynamic data. This is the quantity you feed back into your external thermal model to close the loop:

In [27]:

coupled_result.plot_dynamic_response("Heat generation total [W]")

coupled_result.plot_dynamic_response("Heat generation total [W]")

In [28]:

# The underlying timeseries are available as plain arrays to hand to your thermal model:
data = next(iter(coupled_result.dynamic_data.values()))
time_s = data["Time [s]"]
heat_W = data["Heat generation total [W]"]
print(f"{len(time_s)} samples, peak heat generation = {max(heat_W):.3f} W")

# The underlying timeseries are available as plain arrays to hand to your thermal model: data = next(iter(coupled_result.dynamic_data.values())) time_s = data["Time [s]"] heat_W = data["Heat generation total [W]"] print(f"{len(time_s)} samples, peak heat generation = {max(heat_W):.3f} W")

1000 samples, peak heat generation = 1.361 W

Closing the loop with an external thermal model¶

So far we ran a single prescribed temperature trajectory. In a real coupling you iterate. Your thermal model produces a temperature, Breathe Simulate returns the heat generation under it, that heat updates your thermal model, which produces a new temperature, and so on until the two agree.

Below is a small stand in for your external thermal model. Replace it with your own (for example 3D) model. Here it is a single node energy balance for a forced air cooled cell that turns the heat generation Breathe returns into a cell temperature.

In [29]:

import numpy as np

# A small stand in for your external thermal model. Replace it with your own.
# Single node energy balance for a forced air cooled cell:
#   thermal_mass * dT/dt = Q(t) - cooling_conductance * (T - T_ambient)
thermal_mass_J_per_K = 70.0  # roughly a 21700 cell (mass times specific heat)
cooling_conductance_W_per_K = 0.5  # forced air cooling to the surroundings
T_ambient_degC = 25.0


def external_thermal_model(time_s, heat_W):
    time_s = np.asarray(time_s, dtype=float)
    heat_W = np.asarray(heat_W, dtype=float)
    temperature = np.empty_like(time_s)
    temperature[0] = T_ambient_degC
    for i in range(len(time_s) - 1):
        dt = time_s[i + 1] - time_s[i]
        cooling = cooling_conductance_W_per_K * (temperature[i] - T_ambient_degC)
        temperature[i + 1] = (
            temperature[i] + dt * (heat_W[i] - cooling) / thermal_mass_J_per_K
        )
    return temperature

import numpy as np # A small stand in for your external thermal model. Replace it with your own. # Single node energy balance for a forced air cooled cell: # thermal_mass * dT/dt = Q(t) - cooling_conductance * (T - T_ambient) thermal_mass_J_per_K = 70.0 # roughly a 21700 cell (mass times specific heat) cooling_conductance_W_per_K = 0.5 # forced air cooling to the surroundings T_ambient_degC = 25.0 def external_thermal_model(time_s, heat_W): time_s = np.asarray(time_s, dtype=float) heat_W = np.asarray(heat_W, dtype=float) temperature = np.empty_like(time_s) temperature[0] = T_ambient_degC for i in range(len(time_s) - 1): dt = time_s[i + 1] - time_s[i] cooling = cooling_conductance_W_per_K * (temperature[i] - T_ambient_degC) temperature[i + 1] = ( temperature[i] + dt * (heat_W[i] - cooling) / thermal_mass_J_per_K ) return temperature

We start from a flat guess at the ambient temperature, then repeat: run Breathe Simulate under the current temperature, pass the returned heat generation to the thermal model, and build a new TemperatureProfile from the result. To keep the staggered exchange stable we blend each new temperature with the previous one (under relaxation), and we stop once the temperature stops moving.

In [30]:

profile = TemperatureProfile(time=[0.0, 3600.0], value=[T_ambient_degC, T_ambient_degC])

relaxation = 0.5  # blend new and previous temperature to keep the coupling stable
tolerance_degC = 0.05
max_iterations = 12
convergence = []

for iteration in range(max_iterations):
    result = api.run_sim(
        base_battery="Molicel P45B",
        cycler=cycler_dict,
        initialSoC=1.0,
        cellTemperatureProfile=profile,
    )
    data = next(iter(result.dynamic_data.values()))
    time_s = data["Time [s]"]
    heat_W = data["Heat generation total [W]"]

    previous = np.interp(time_s, profile.to_dict()["time"], profile.to_dict()["value"])
    modelled = external_thermal_model(time_s, heat_W)
    updated = (1.0 - relaxation) * previous + relaxation * modelled

    change_degC = float(np.max(np.abs(updated - previous)))
    convergence.append(change_degC)
    print(
        f"iteration {iteration}: max temperature change = {change_degC:.3f} degC, peak = {updated.max():.1f} degC"
    )

    profile = TemperatureProfile(time=time_s, value=updated.tolist())
    if change_degC < tolerance_degC:
        print("Converged: the electrochemical and thermal models agree.")
        break

profile = TemperatureProfile(time=[0.0, 3600.0], value=[T_ambient_degC, T_ambient_degC]) relaxation = 0.5 # blend new and previous temperature to keep the coupling stable tolerance_degC = 0.05 max_iterations = 12 convergence = [] for iteration in range(max_iterations): result = api.run_sim( base_battery="Molicel P45B", cycler=cycler_dict, initialSoC=1.0, cellTemperatureProfile=profile, ) data = next(iter(result.dynamic_data.values())) time_s = data["Time [s]"] heat_W = data["Heat generation total [W]"] previous = np.interp(time_s, profile.to_dict()["time"], profile.to_dict()["value"]) modelled = external_thermal_model(time_s, heat_W) updated = (1.0 - relaxation) * previous + relaxation * modelled change_degC = float(np.max(np.abs(updated - previous))) convergence.append(change_degC) print( f"iteration {iteration}: max temperature change = {change_degC:.3f} degC, peak = {updated.max():.1f} degC" ) profile = TemperatureProfile(time=time_s, value=updated.tolist()) if change_degC < tolerance_degC: print("Converged: the electrochemical and thermal models agree.") break

iteration 0: max temperature change = 1.121 degC, peak = 26.1 degC
iteration 1: max temperature change = 0.547 degC, peak = 26.7 degC
iteration 2: max temperature change = 0.267 degC, peak = 26.9 degC
iteration 3: max temperature change = 0.131 degC, peak = 27.1 degC
iteration 4: max temperature change = 0.064 degC, peak = 27.1 degC
iteration 5: max temperature change = 0.031 degC, peak = 27.2 degC
Converged: the electrochemical and thermal models agree.

The converged temperature and heat generation are now self consistent. We can look at the final coupled temperature and how quickly the loop converged.

In [31]:

profile.plot()

profile.plot()

In [32]:

import plotly.graph_objects as go

fig = go.Figure(go.Scatter(y=convergence, mode="lines+markers"))
fig.update_layout(
    title="Coupling convergence",
    xaxis_title="Iteration",
    yaxis_title="Max temperature change [°C]",
)
fig

import plotly.graph_objects as go fig = go.Figure(go.Scatter(y=convergence, mode="lines+markers")) fig.update_layout( title="Coupling convergence", xaxis_title="Iteration", yaxis_title="Max temperature change [°C]", ) fig

Notes¶

1. Ambient temperature and heat transfer coefficient are ignored when a cellTemperatureProfile is provided.

2. cellTemperatureProfile cannot be combined with isothermal=True or with batch sweeps (initial conditions given as a list).

3. A plain dict {"time": [...], "value": [...]} is also accepted, and you can build a profile straight from a DataFrame with TemperatureProfile.from_dataframe(df, time_column=..., temperature_column=...).

Version 8e4a326f