diff --git a/CHANGELOG.md b/CHANGELOG.md index 0232023c2e..216453b669 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -36,6 +36,8 @@ - `RegPower` now evaluates the same operations in `evaluate` and CasADi, taking a literal exponent of 0.5 or 2 as square roots, so they agree bitwise; CasADi used a `hypot` form. With `delta=0`, `evaluate` and the Jacobian use `base ** exponent`, as CasADi did. ([#5955](https://github.com/pybamm-team/PyBaMM/pull/5955)) - CasADi linear and `pchip` interpolants on a uniform grid now find a point's segment by binary search. The `"exact"` lookup could pick the neighbouring segment at or one ulp beside a knot, giving the Jacobian that segment's slope. ([#5955](https://github.com/pybamm-team/PyBaMM/pull/5955)) - `QuickPlot`, and so `Simulation.plot`, now uses matplotlib's color cycle (`axes.prop_cycle`) for default colors when a style sheet or rcParams changes it; otherwise it keeps its red, blue, black, green, magenta and cyan colors. ([#5953](https://github.com/pybamm-team/PyBaMM/pull/5953)) +- Irreversible heating now includes the resistive loss across the SEI film when `"SEI film resistance"` is `"distributed"` or `"average"`, so it reaches `"Total heating [W]"` and the cell temperature. ([#5841](https://github.com/pybamm-team/PyBaMM/pull/5841)) +- Irreversible heating now includes the heat from SEI growth, SEI on cracks and lithium plating side reactions, so their electrical losses reach `"Total heating [W]"` and the cell temperature. ([#5875](https://github.com/pybamm-team/PyBaMM/pull/5875)) ## Optimizations diff --git a/packages/pybamm/src/pybamm/models/submodels/interface/base_interface.py b/packages/pybamm/src/pybamm/models/submodels/interface/base_interface.py index 6624b646f8..cbe7d2851c 100644 --- a/packages/pybamm/src/pybamm/models/submodels/interface/base_interface.py +++ b/packages/pybamm/src/pybamm/models/submodels/interface/base_interface.py @@ -391,6 +391,10 @@ def _get_standard_sei_film_overpotential_variables(self, eta_sei): } return variables + # Size average, as for eta_r in "_get_standard_overpotential_variables" + if eta_sei.domain and eta_sei.domain[0].endswith("particle size"): + eta_sei = pybamm.size_average(eta_sei) + # Average, and broadcast if necessary eta_sei_av = pybamm.x_average(eta_sei) if eta_sei.domain == []: diff --git a/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py b/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py index fa5e0b590f..2b7785547f 100644 --- a/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py +++ b/packages/pybamm/src/pybamm/models/submodels/thermal/base_thermal.py @@ -146,8 +146,13 @@ def _get_standard_coupled_variables(self, variables): f"Positive electrode {phase}volumetric interfacial current density [A.m-3]" ] eta_r_p = variables[f"Positive electrode {phase}reaction overpotential [V]"] - # Irreversible electrochemical heating - Q_rxn_p += a_j_p * eta_r_p + eta_sei_p = variables[ + f"Positive electrode {phase}SEI film overpotential [V]" + ] + # Irreversible heating. eta_r is the kinetic overpotential net of the + # film drop, so eta_r - eta_sei adds the resistive loss in the SEI film + Q_rxn_p += a_j_p * (eta_r_p - eta_sei_p) + Q_rxn_p += self._side_reaction_heating(variables, "positive", phase) # Reversible electrochemical heating dUdT_p = variables[f"Positive electrode {phase}entropic change [V.K-1]"] Q_rev_p += a_j_p * T_p * dUdT_p @@ -176,10 +181,21 @@ def _get_standard_coupled_variables(self, variables): ocp_options = (ocp_options,) * len(phase_names) if self.options.electrode_types["negative"] == "planar": - i_n = variables["Lithium metal total interfacial current density [A.m-2]"] + # The total current splits between lithium plating and SEI growth + j_pl = variables["Lithium metal plating current density [A.m-2]"] eta_r_n = variables["Lithium metal interface reaction overpotential [V]"] + eta_sei_n = variables["Negative electrode SEI film overpotential [V]"] + Q_rxn_n = j_pl * (eta_r_n - eta_sei_n) + if self.options.negative["SEI"] not in ["none", "constant"]: + j_sei = variables[ + "Negative electrode SEI interfacial current density [A.m-2]" + ] + delta_phi_n = variables[ + "Lithium metal interface surface potential difference [V]" + ] + Q_rxn_n += j_sei * (delta_phi_n - self.param.n.prim.U_sei) Q_rxn_n = pybamm.PrimaryBroadcast( - i_n * eta_r_n / self.param.n.L, + Q_rxn_n / self.param.n.L, ["negative electrode"], "current collector", ) @@ -199,8 +215,13 @@ def _get_standard_coupled_variables(self, variables): eta_r_n = variables[ f"Negative electrode {phase}reaction overpotential [V]" ] - # Irreversible electrochemical heating - Q_rxn_n += a_j_n * eta_r_n + eta_sei_n = variables[ + f"Negative electrode {phase}SEI film overpotential [V]" + ] + # Irreversible heating. eta_r is the kinetic overpotential net of the + # film drop, so eta_r - eta_sei adds the resistive loss in the SEI film + Q_rxn_n += a_j_n * (eta_r_n - eta_sei_n) + Q_rxn_n += self._side_reaction_heating(variables, "negative", phase) # Reversible electrochemical heating dUdT_n = variables[f"Negative electrode {phase}entropic change [V.K-1]"] Q_rev_n += a_j_n * T_n * dUdT_n @@ -383,6 +404,35 @@ def _current_collector_heating(self, variables): Q_s_cp = self.param.p.sigma_cc * pybamm.grad_squared(phi_s_cp) return Q_s_cn, Q_s_cp + def _side_reaction_heating(self, variables, domain, phase_name): + """ + Irreversible heating from SEI growth and lithium plating in a porous + electrode, as each side-reaction current times ``delta_phi - U`` for that + reaction, which includes any SEI film drop. + """ + Domain = domain.capitalize() + phase = phase_name.strip() or "primary" + domain_options = getattr(self.options, domain) + phase_options = getattr(domain_options, phase) + delta_phi = variables[f"{Domain} electrode surface potential difference [V]"] + prefix = f"{Domain} electrode {phase_name}" + suffix = "volumetric interfacial current density [A.m-3]" + + Q = 0 + if phase_options["SEI"] not in ["none", "constant"]: + U_sei = self.param.domain_params[domain].phase_params[phase].U_sei + reactions = ["SEI"] + if domain_options["SEI on cracks"] == "true": + reactions.append("SEI on cracks") + for reaction in reactions: + a_j_sei = variables[f"{prefix}{reaction} {suffix}"] + Q += a_j_sei * (delta_phi - U_sei) + if phase_options["lithium plating"] != "none": + # Lithium plating has zero open-circuit potential + a_j_pl = variables[f"{prefix}lithium plating {suffix}"] + Q += a_j_pl * delta_phi + return Q + def _heat_of_mixing(self, variables): """Compute heat of mixing source terms.""" diff --git a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py index 7f2609c4e8..01b3aa4454 100644 --- a/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py +++ b/packages/pybamm/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_thermal_models.py @@ -222,3 +222,83 @@ def R_contact(T): # cell heats up above the reference, so R (and hence the contact # overpotential magnitude) is larger than the constant case np.testing.assert_array_less(np.abs(dphi_const[1:]), np.abs(dphi_fn[1:])) + + @pytest.mark.parametrize( + "model, options, parameter_set, current", + [ + ( + pybamm.lithium_ion.SPMe, + {"SEI film resistance": "distributed"}, + "OKane2022", + 5, + ), + ( + pybamm.lithium_ion.SPMe, + {"SEI film resistance": "average"}, + "OKane2022", + 5, + ), + ( + pybamm.lithium_ion.DFN, + {"SEI film resistance": "distributed"}, + "OKane2022", + 5, + ), + ( + pybamm.lithium_ion.MPM, + {"SEI film resistance": "distributed"}, + "OKane2022", + 5, + ), + ( + pybamm.lithium_ion.SPMe, + {"working electrode": "positive", "SEI film resistance": "distributed"}, + "Xu2019", + 2.4e-3, + ), + ], + ) + def test_sei_film_heating(self, model, options, parameter_set, current): + # With the temperature held at ambient, the extra heat generated when the SEI + # film resistance is switched on must match the extra electrical energy lost + parameter_values = pybamm.ParameterValues(parameter_set) + # Xu2019 has no current collector or thermal parameters, so take them from + # Chen2020 + chen2020 = pybamm.ParameterValues("Chen2020") + parameter_values.update( + {k: chen2020[k] for k in chen2020 if k not in parameter_values}, + check_already_exists=False, + ) + if model is pybamm.lithium_ion.MPM: + parameter_values = pybamm.get_size_distribution_parameters(parameter_values) + parameter_values.update( + { + "Current function [A]": current, + "Initial SEI thickness [m]": 5e-8, + "Total heat transfer coefficient [W.m-2.K-1]": 1e6, + } + ) + t = np.linspace(0, 600, 61) + + voltage, heating = {}, {} + for film in ["none", options["SEI film resistance"]]: + sim = pybamm.Simulation( + model( + { + **options, + "SEI": "constant", + "SEI film resistance": film, + "thermal": "lumped", + } + ), + parameter_values=parameter_values, + ) + sol = sim.solve([0, 600]) + voltage[film] = sol["Voltage [V]"](t) + heating[film] = sol["Total heating [W]"](t) + + film = options["SEI film resistance"] + electrical_loss = np.trapezoid(current * (voltage["none"] - voltage[film]), t) + extra_heat = np.trapezoid(heating[film] - heating["none"], t) + assert electrical_loss > 0 + np.testing.assert_allclose(extra_heat, electrical_loss, rtol=1e-2) diff --git a/packages/pybamm/tests/integration/test_models/test_submodels/test_thermal/test_energy_balance.py b/packages/pybamm/tests/integration/test_models/test_submodels/test_thermal/test_energy_balance.py new file mode 100644 index 0000000000..6ef7fdcf3b --- /dev/null +++ b/packages/pybamm/tests/integration/test_models/test_submodels/test_thermal/test_energy_balance.py @@ -0,0 +1,123 @@ +import numpy as np +import pytest + +import pybamm + + +class TestEnergyBalance: + @pytest.mark.parametrize( + "options, parameter_set, current, initial_soc", + [ + ( + { + "SEI": "solvent-diffusion limited", + "SEI on cracks": "true", + "particle mechanics": "swelling and cracking", + "stress-induced diffusion": "true", + "lithium plating": "partially reversible", + }, + "OKane2022", + -10.0, + 0.1, + ), + ( + { + "particle phases": ("2", "1"), + "open-circuit potential": (("single", "current sigmoid"), "single"), + "SEI": "solvent-diffusion limited", + "lithium plating": "partially reversible", + "surface form": "algebraic", + }, + "Bonkile2024", + -10.0, + 0.1, + ), + ( + {"working electrode": "positive", "SEI": "reaction limited"}, + "Xu2019", + 2.4e-3, + None, + ), + ], + ) + def test_side_reaction_heat_matches_electrical_loss( + self, options, parameter_set, current, initial_soc + ): + # Ohmic plus irreversible heat must equal the power drawn by all reactions + # at their open-circuit potentials minus the terminal power + model = pybamm.lithium_ion.DFN( + { + "calculate heat source for isothermal models": "true", + "SEI film resistance": "distributed", + "total interfacial current density as a state": "true", + **options, + } + ) + param = pybamm.ParameterValues(parameter_set) + # Xu2019 has no SEI growth, current collector or thermal parameters + okane2022 = pybamm.ParameterValues("OKane2022") + param.update( + {k: okane2022[k] for k in okane2022 if k not in param}, + check_already_exists=False, + ) + param.update( + { + "Current function [A]": current, + "Ambient temperature [K]": 268.15, + "Initial temperature [K]": 268.15, + "Initial SEI thickness [m]": 50e-9, + } + ) + if parameter_set == "Xu2019": + # Grow SEI fast enough on the lithium metal for its heat to matter + param["SEI reaction exchange current density [A.m-2]"] = 1e-3 + sim = pybamm.Simulation(model, parameter_values=param) + solution = sim.solve([0, 1800], initial_soc=initial_soc) + t = solution.t + + area = param["Electrode width [m]"] * param["Electrode height [m]"] + stored_power = 0 + for domain in ["negative", "positive"]: + Domain = domain.capitalize() + if model.options.electrode_types[domain] == "planar": + # Lithium plating has zero open-circuit potential + j_sei = solution[ + "Negative electrode SEI interfacial current density [A.m-2]" + ] + stored_power += ( + j_sei(t) * param["SEI open-circuit potential [V]"] * area + ) + continue + x = sim.mesh[f"{domain} electrode"].nodes + phases = model.options.phases[domain] + for phase in phases: + phase_name = f"{phase} " if len(phases) > 1 else "" + prefix = f"{phase.capitalize()}: " if len(phases) > 1 else "" + U_sei = param[f"{prefix}SEI open-circuit potential [V]"] + + def field(name, prefix=f"{Domain} electrode {phase_name}", x=x): + return solution[f"{prefix}{name}"](t=t, x=x) + + power_density = field( + "volumetric interfacial current density [A.m-3]" + ) * field("open-circuit potential [V]") + for reaction in ["SEI ", "SEI on cracks "]: + power_density += ( + field( + f"{reaction}volumetric interfacial current density [A.m-3]" + ) + * U_sei + ) + # The default mesh is uniform within each electrode + thickness = param[f"{Domain} electrode thickness [m]"] + stored_power += power_density.mean(axis=0) * thickness * area + + heat = solution["Ohmic heating [W]"](t) + solution[ + "Irreversible electrochemical heating [W]" + ](t) + electrical_loss = -stored_power - solution["Current [A]"](t) * solution[ + "Voltage [V]" + ](t) + np.testing.assert_allclose( + np.trapezoid(heat, t), np.trapezoid(electrical_loss, t), rtol=1e-3 + ) diff --git a/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_half_cell_tests.py b/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_half_cell_tests.py index e8fb3d9eb6..66f7885a81 100644 --- a/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_half_cell_tests.py +++ b/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_half_cell_tests.py @@ -58,6 +58,15 @@ def test_well_posed_reaction_limited_sei(self): options = {"SEI": "reaction limited"} self.check_well_posedness(options) + def test_well_posed_reaction_limited_sei_lumped_thermal(self): + options = { + "SEI": "reaction limited", + "thermal": "lumped", + "SEI film resistance": "distributed", + "total interfacial current density as a state": "true", + } + self.check_well_posedness(options) + def test_well_posed_asymmetric_reaction_limited_sei(self): options = {"SEI": "reaction limited (asymmetric)"} self.check_well_posedness(options) diff --git a/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py b/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py index 965f5b8949..1a0cc7b534 100644 --- a/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py +++ b/packages/pybamm/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py @@ -372,6 +372,19 @@ def test_well_posed_sei_reaction_limited_on_cracks(self): } self.check_well_posedness(options) + def test_well_posed_side_reactions_lumped_thermal(self): + options = { + "SEI": "reaction limited", + "SEI on cracks": "true", + "particle mechanics": "swelling and cracking", + "lithium plating": "reversible", + "thermal": "lumped", + "SEI film resistance": "distributed", + "stress-induced diffusion": "true", + "total interfacial current density as a state": "true", + } + self.check_well_posedness(options) + def test_well_posed_sei_solvent_diffusion_limited_on_cracks(self): options = { "SEI": "solvent-diffusion limited",