Issue 1778 - heat of mixing (pybamm-team#2837)

* fix * fix exception * added a Newfile * pseudo * pybamm-team#1778 add heat of mixing option * pybamm-team#1778 fixed indentation * pybamm-team#1778 fixed test * pybamm-team#1778 fixed test * pybamm-team#1788 added heat of mixing tests * First commit * pybamm-team#1778 fixed test * First iteration * pybamm-team#1788 try to fix diff of U * first push from mac * Fixed full broadcast * Fixed domain in the integral * pybamm-team#1778 fix parameter values example * Added children[0].children[0] * Added half cell * style: pre-commit fixes * pybamm-team#1788 revert some unrelated changes * change ocv hack * Revert "change ocv hack" This reverts commit 3ad0c75. * pybamm-team#1778 fix heat of mixing * pybamm-team#1839 fix thermal submodels to take x_average * pybamm-team#1778 add example heat of mixing * ruff * style: pre-commit fixes * pybamm-team#1778 fix lead acid models * fix SPM with the right broadcast of temperature * pybamm-team#1778 refactor heat of mixing * style: pre-commit fixes * Add Richardson2021 citation for heat of mixing * Fixed heat of mixing sign * Rewritten heat of mixing example, added comparison plot to compare with the paper * Fixed ruff formatting errors in the heat of mixing example script * pybamm-team#1778 fix x-full * pybamm-team#1778 fix tests * pybamm-team#1778 improve coverage * style: pre-commit fixes * pybamm-team#1778 Valentin's suggestion to remove the dUdsto function and compute the derivative in the model directly instead * pybamm-team#1778 fix failing test * pybamm-team#1778 add heat of mixing to CHANGELOG --------- Co-authored-by: Alec Bills <48105066+abillscmu@users.noreply.github.com> Co-authored-by: smitasahu2 <smitasahuiitd@gmail.com> Co-authored-by: Afgr1087 <andres.galvis@fem.unicamp.br> Co-authored-by: Ferran Brosa Planella <Ferran.Brosa-Planella@warwick.ac.uk> Co-authored-by: Ivan Korotkin <i.korotkin@soton.ac.uk> Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
js1tr3 · Aug 12, 2024 · ffeaf85 · ffeaf85
1 parent b7f20b4
commit ffeaf85
Show file tree

Hide file tree

Showing 14 changed files with 326 additions and 17 deletions.
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -12,7 +12,8 @@
 - Add support for BPX version 0.4.0 which allows for blended electrodes and user-defined parameters in BPX([#3414](https://github.com/pybamm-team/PyBaMM/pull/3414))
 - Added `by_submodel` feature in `print_parameter_info` method to allow users to print parameters and types of submodels in a tabular and readable format ([#3628](https://github.com/pybamm-team/PyBaMM/pull/3628))
 - Added `WyciskOpenCircuitPotential` for differential capacity hysteresis state open-circuit potential submodel ([#3593](https://github.com/pybamm-team/PyBaMM/pull/3593))
-- Added `time` as an option for `Experiment.termination`. Now allows solving up to a user-specified time while also allowing different cycles and steps in an experiment to be handled normally. ([#4073](https://github.com/pybamm-team/PyBaMM/pull/4073))
+- Transport efficiency submodel has new options from the literature relating to different tortuosity factor models and also a new option called "tortuosity factor" for specifying the value or function directly as parameters ([#3437](https://github.com/pybamm-team/PyBaMM/pull/3437))
+- Heat of mixing source term can now be included into thermal models ([#2837](https://github.com/pybamm-team/PyBaMM/pull/2837))
 
 ## Bug Fixes
 

diff --git a/examples/scripts/compare_lithium_ion_heat_of_mixing.py b/examples/scripts/compare_lithium_ion_heat_of_mixing.py
@@ -0,0 +1,112 @@
+#
+# Compare SPMe model with and without heat of mixing
+#
+import pybamm
+import numpy as np
+import matplotlib.pyplot as plt
+
+pybamm.set_logging_level("INFO")
+
+# load models
+models = [
+    pybamm.lithium_ion.SPMe(
+        {"heat of mixing": "true", "thermal": "lumped"}, name="SPMe with heat of mixing"
+    ),
+    pybamm.lithium_ion.SPMe({"thermal": "lumped"}, name="SPMe without heat of mixing"),
+]
+
+# set parametrisation (Ecker et al., 2015)
+parameter_values = pybamm.ParameterValues("Ecker2015")
+
+# set mesh
+# (increased number of points in particles to avoid oscillations for high C-rates)
+var_pts = {"x_n": 16, "x_s": 8, "x_p": 16, "r_n": 128, "r_p": 128}
+
+# set the constant current discharge C-rate
+C_rate = 5
+
+# set the simulation time and increase the number of points
+# for better integration in time
+t_eval = np.linspace(0, 720, 360)
+
+# set up an extra plot with the heat of mixing vs time in each electrode and
+# the integrated heat of mixing vs time in each electrode to compare with
+# Figure 6(a) from Richardson et al. (2021)
+fig, axs = plt.subplots(2, len(models), figsize=(12, 7))
+
+# extract some of the parameters
+L_n = parameter_values["Negative electrode thickness [m]"]
+L_p = parameter_values["Positive electrode thickness [m]"]
+A = parameter_values["Electrode width [m]"] * parameter_values["Electrode height [m]"]
+
+# create and run simulations
+sims = []
+for m, model in enumerate(models):
+    sim = pybamm.Simulation(
+        model, parameter_values=parameter_values, var_pts=var_pts, C_rate=C_rate
+    )
+    sim.solve(t_eval)
+    sims.append(sim)
+
+    # extract solution for plotting
+    sol = sim.solution
+
+    # extract variables from the solution
+    time = sol["Time [h]"].entries
+    Q_mix = sol["Heat of mixing [W.m-3]"].entries
+
+    # heat of mixing in negative and positive electrodes multiplied by the electrode
+    # width, represents the integral of heat of mixing term across each of the
+    # electrodes (W.m-2)
+    Q_mix_n = Q_mix[0, :] * L_n
+    Q_mix_p = Q_mix[-1, :] * L_p
+
+    # heat of mixing integrals (J.m-2)
+    Q_mix_n_int = 0.0
+    Q_mix_p_int = 0.0
+
+    # data for plotting
+    Q_mix_n_plt = []
+    Q_mix_p_plt = []
+
+    # performs integration in time
+    for i, t in enumerate(time[1:]):
+        dt = (t - time[i]) * 3600  # seconds
+        Q_mix_n_avg = (Q_mix_n[i] + Q_mix_n[i + 1]) * 0.5
+        Q_mix_p_avg = (Q_mix_p[i] + Q_mix_p[i + 1]) * 0.5
+        # convert J to kJ and divide the integral by the electrode area A to compare
+        # with Figure 6(a) from Richardson et al. (2021)
+        Q_mix_n_int += Q_mix_n_avg * dt / 1000 / A
+        Q_mix_p_int += Q_mix_p_avg * dt / 1000 / A
+        Q_mix_n_plt.append(Q_mix_n_int)
+        Q_mix_p_plt.append(Q_mix_p_int)
+
+    # plots heat of mixing in each electrode vs time in minutes
+    axs[0, m].plot(time * 60, Q_mix_n, ls="-", label="Negative electrode")
+    axs[0, m].plot(time * 60, Q_mix_p, ls="--", label="Positive electrode")
+    axs[0, m].set_title(f"{model.name}")
+    axs[0, m].set_xlabel("Time [min]")
+    axs[0, m].set_ylabel("Heat of mixing [W.m-2]")
+    axs[0, m].grid(True)
+    axs[0, m].legend()
+
+    # plots integrated heat of mixing in each electrode vs time in minutes
+    axs[1, m].plot(time[1:] * 60, Q_mix_n_plt, ls="-", label="Negative electrode")
+    axs[1, m].plot(time[1:] * 60, Q_mix_p_plt, ls="--", label="Positive electrode")
+    axs[1, m].set_xlabel("Time [min]")
+    axs[1, m].set_ylabel("Integrated heat of mixing [kJ.m-2]")
+    axs[1, m].grid(True)
+    axs[1, m].legend()
+
+# plot
+pybamm.dynamic_plot(
+    sims,
+    output_variables=[
+        "X-averaged cell temperature [K]",
+        "X-averaged heat of mixing [W.m-3]",
+        "X-averaged total heating [W.m-3]",
+        "Heat of mixing [W.m-3]",
+        "Voltage [V]",
+        "Current [A]",
+    ],
+)
diff --git a/pybamm/CITATIONS.bib b/pybamm/CITATIONS.bib
@@ -418,6 +418,17 @@ @article{Richardson2020
   doi = {10.1016/j.electacta.2020.135862},
 }
 
+@article{Richardson2021,
+  title = {Heat Generation and a Conservation Law for Chemical Energy in Li-ion Batteries},
+  author = {Richardson, Giles and Korotkin, Ivan},
+  journal = {Electrochimica Acta},
+  volume = {392},
+  pages = {138909},
+  year = {2021},
+  publisher = {Elsevier},
+  doi = {10.1016/j.electacta.2021.138909},
+}
+
 @article{Sripad2020,
   title={Kinetics of lithium electrodeposition and stripping},
   author={Sripad, Shashank and Korff, Daniel and DeCaluwe, Steven C and Viswanathan, Venkatasubramanian},

diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
@@ -238,6 +238,7 @@ def __init__(self, extra_options):
                 "integrated",
             ],
             "exchange-current density": ["single", "current sigmoid"],
+            "heat of mixing": ["false", "true"],
             "hydrolysis": ["false", "true"],
             "intercalation kinetics": [
                 "symmetric Butler-Volmer",
@@ -511,6 +512,11 @@ def __init__(self, extra_options):
 
         # Options not yet compatible with particle-size distributions
         if options["particle size"] == "distribution":
+            if options["heat of mixing"] != "false":
+                raise NotImplementedError(
+                    "Heat of mixing submodels do not yet support particle-size "
+                    "distributions."
+                )
             if options["lithium plating"] != "none":
                 raise NotImplementedError(
                     "Lithium plating submodels do not yet support particle-size "
@@ -1243,7 +1249,10 @@ def set_thermal_submodel(self):
             if self.options["dimensionality"] == 0:
                 thermal_submodel = pybamm.thermal.pouch_cell.OneDimensionalX
 
-        self.submodels["thermal"] = thermal_submodel(self.param, self.options)
+        x_average = getattr(self, "x_average", False)
+        self.submodels["thermal"] = thermal_submodel(
+            self.param, self.options, x_average
+        )
 
     def set_current_collector_submodel(self):
         if self.options["current collector"] in ["uniform"]:

diff --git a/pybamm/models/submodels/thermal/base_thermal.py b/pybamm/models/submodels/thermal/base_thermal.py
@@ -2,6 +2,7 @@
 # Base class for thermal effects
 #
 import pybamm
+import numpy as np
 
 
 class BaseThermal(pybamm.BaseSubModel):
@@ -16,8 +17,12 @@ class BaseThermal(pybamm.BaseSubModel):
         A dictionary of options to be passed to the model.
     """
 
-    def __init__(self, param, options=None):
+    def __init__(self, param, options=None, x_average=False):
         super().__init__(param, options=options)
+        self.x_average = x_average
+
+        if self.options["heat of mixing"] == "true":
+            pybamm.citations.register("Richardson2021")
 
     def _get_standard_fundamental_variables(self, T_dict):
         """
@@ -175,15 +180,20 @@ def _get_standard_coupled_variables(self, variables):
             Q_rev_n, pybamm.FullBroadcast(0, "separator", "current collector"), Q_rev_p
         )
 
+        # Heat of mixing
+        Q_mix_s_n, Q_mix_s_s, Q_mix_s_p = self._heat_of_mixing(variables)
+        Q_mix = pybamm.concatenation(Q_mix_s_n, Q_mix_s_s, Q_mix_s_p)
+
         # Total heating
-        Q = Q_ohm + Q_rxn + Q_rev
+        Q = Q_ohm + Q_rxn + Q_rev + Q_mix
 
         # Compute the X-average over the entire cell, including current collectors
         # Note: this can still be a function of y and z for higher-dimensional pouch
         # cell models
         Q_ohm_av = self._x_average(Q_ohm, Q_ohm_s_cn, Q_ohm_s_cp)
         Q_rxn_av = self._x_average(Q_rxn, 0, 0)
         Q_rev_av = self._x_average(Q_rev, 0, 0)
+        Q_mix_av = self._x_average(Q_mix, 0, 0)
         Q_av = self._x_average(Q, Q_ohm_s_cn, Q_ohm_s_cp)
 
         # Compute the integrated heat source per unit simulated electrode-pair area
@@ -192,11 +202,14 @@ def _get_standard_coupled_variables(self, variables):
         Q_ohm_Wm2 = Q_ohm_av * param.L
         Q_rxn_Wm2 = Q_rxn_av * param.L
         Q_rev_Wm2 = Q_rev_av * param.L
+        Q_mix_Wm2 = Q_mix_av * param.L
         Q_Wm2 = Q_av * param.L
+
         # Now average over the electrode height and width
         Q_ohm_Wm2_av = self._yz_average(Q_ohm_Wm2)
         Q_rxn_Wm2_av = self._yz_average(Q_rxn_Wm2)
         Q_rev_Wm2_av = self._yz_average(Q_rev_Wm2)
+        Q_mix_Wm2_av = self._yz_average(Q_mix_Wm2)
         Q_Wm2_av = self._yz_average(Q_Wm2)
 
         # Compute total heat source terms (in W) over the *entire cell volume*, not
@@ -208,6 +221,7 @@ def _get_standard_coupled_variables(self, variables):
         Q_ohm_W = Q_ohm_Wm2_av * n_elec * A
         Q_rxn_W = Q_rxn_Wm2_av * n_elec * A
         Q_rev_W = Q_rev_Wm2_av * n_elec * A
+        Q_mix_W = Q_mix_Wm2_av * n_elec * A
         Q_W = Q_Wm2_av * n_elec * A
 
         # Compute volume-averaged heat source terms over the *entire cell volume*, not
@@ -216,6 +230,7 @@ def _get_standard_coupled_variables(self, variables):
         Q_ohm_vol_av = Q_ohm_W / V
         Q_rxn_vol_av = Q_rxn_W / V
         Q_rev_vol_av = Q_rev_W / V
+        Q_mix_vol_av = Q_mix_W / V
         Q_vol_av = Q_W / V
 
         # Effective heat capacity
@@ -228,6 +243,7 @@ def _get_standard_coupled_variables(self, variables):
                 "Ohmic heating [W.m-3]": Q_ohm,
                 "X-averaged Ohmic heating [W.m-3]": Q_ohm_av,
                 "Volume-averaged Ohmic heating [W.m-3]": Q_ohm_vol_av,
+                "Volume-averaged heat of mixing [W.m-3]": Q_mix_vol_av,
                 "Ohmic heating per unit electrode-pair area [W.m-2]": Q_ohm_Wm2,
                 "Ohmic heating [W]": Q_ohm_W,
                 # Irreversible
@@ -244,6 +260,12 @@ def _get_standard_coupled_variables(self, variables):
                 "Volume-averaged reversible heating [W.m-3]": Q_rev_vol_av,
                 "Reversible heating per unit electrode-pair area " "[W.m-2]": Q_rev_Wm2,
                 "Reversible heating [W]": Q_rev_W,
+                # Mixing
+                "Heat of mixing [W.m-3]": Q_mix,
+                "X-averaged heat of mixing [W.m-3]": Q_mix_av,
+                "Volume-averaged heating of mixing [W.m-3]": Q_mix_vol_av,
+                "Heat of mixing per unit electrode-pair area " "[W.m-2]": Q_mix_Wm2,
+                "Heat of mixing [W]": Q_mix_W,
                 # Total
                 "Total heating [W.m-3]": Q,
                 "X-averaged total heating [W.m-3]": Q_av,
@@ -290,6 +312,74 @@ 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 _heat_of_mixing(self, variables):
+        """Compute heat of mixing source terms."""
+        param = self.param
+
+        if self.options["heat of mixing"] == "true":
+            F = pybamm.constants.F.value
+            pi = np.pi
+
+            # Compute heat of mixing in negative electrode
+            if self.options.electrode_types["negative"] == "planar":
+                Q_mix_s_n = pybamm.FullBroadcast(
+                    0, ["negative electrode"], "current collector"
+                )
+            else:
+                a_n = variables["Negative electrode surface area to volume ratio [m-1]"]
+                R_n = variables["Negative particle radius [m]"]
+                N_n = a_n / (4 * pi * R_n**2)
+                if self.x_average:
+                    c_n = variables[
+                        "X-averaged negative particle concentration [mol.m-3]"
+                    ]
+                    T_n = variables["X-averaged negative electrode temperature [K]"]
+                else:
+                    c_n = variables["Negative particle concentration [mol.m-3]"]
+                    T_n = variables["Negative electrode temperature [K]"]
+                T_n_part = pybamm.PrimaryBroadcast(T_n, ["negative particle"])
+                dc_n_dr2 = pybamm.inner(pybamm.grad(c_n), pybamm.grad(c_n))
+                D_n = param.n.prim.D(c_n, T_n_part)
+                dUeq_n = param.n.prim.U(c_n / param.n.prim.c_max, T_n_part).diff(c_n)
+                integrand_r_n = D_n * dc_n_dr2 * dUeq_n
+                integration_variable_r_n = [
+                    pybamm.SpatialVariable("r", domain=integrand_r_n.domain)
+                ]
+                integral_r_n = pybamm.Integral(integrand_r_n, integration_variable_r_n)
+                Q_mix_s_n = -F * N_n * integral_r_n
+
+            # Compute heat of mixing in positive electrode
+            a_p = variables["Positive electrode surface area to volume ratio [m-1]"]
+            R_p = variables["Positive particle radius [m]"]
+            N_p = a_p / (4 * pi * R_p**2)
+            if self.x_average:
+                c_p = variables["X-averaged positive particle concentration [mol.m-3]"]
+                T_p = variables["X-averaged positive electrode temperature [K]"]
+            else:
+                c_p = variables["Positive particle concentration [mol.m-3]"]
+                T_p = variables["Positive electrode temperature [K]"]
+            T_p_part = pybamm.PrimaryBroadcast(T_p, ["positive particle"])
+            dc_p_dr2 = pybamm.inner(pybamm.grad(c_p), pybamm.grad(c_p))
+            D_p = param.p.prim.D(c_p, T_p_part)
+            dUeq_p = param.p.prim.U(c_p / param.p.prim.c_max, T_p_part).diff(c_p)
+            integrand_r_p = D_p * dc_p_dr2 * dUeq_p
+            integration_variable_r_p = [
+                pybamm.SpatialVariable("r", domain=integrand_r_p.domain)
+            ]
+            integral_r_p = pybamm.Integral(integrand_r_p, integration_variable_r_p)
+            Q_mix_s_p = -F * N_p * integral_r_p
+            Q_mix_s_s = pybamm.FullBroadcast(0, ["separator"], "current collector")
+        else:
+            Q_mix_s_n = pybamm.FullBroadcast(
+                0, ["negative electrode"], "current collector"
+            )
+            Q_mix_s_p = pybamm.FullBroadcast(
+                0, ["positive electrode"], "current collector"
+            )
+            Q_mix_s_s = pybamm.FullBroadcast(0, ["separator"], "current collector")
+
+        return Q_mix_s_n, Q_mix_s_s, Q_mix_s_p
+
     def _x_average(self, var, var_cn, var_cp):
         """
         Computes the X-average over the whole cell (including current collectors)

diff --git a/pybamm/models/submodels/thermal/isothermal.py b/pybamm/models/submodels/thermal/isothermal.py
@@ -18,8 +18,8 @@ class Isothermal(BaseThermal):
         A dictionary of options to be passed to the model.
     """
 
-    def __init__(self, param, options=None):
-        super().__init__(param, options=options)
+    def __init__(self, param, options=None, x_average=False):
+        super().__init__(param, options=options, x_average=x_average)
 
     def get_fundamental_variables(self):
         # Set the x-averaged temperature to the ambient temperature, which can be

diff --git a/pybamm/models/submodels/thermal/lumped.py b/pybamm/models/submodels/thermal/lumped.py
@@ -20,8 +20,8 @@ class Lumped(BaseThermal):
 
     """
 
-    def __init__(self, param, options=None):
-        super().__init__(param, options=options)
+    def __init__(self, param, options=None, x_average=False):
+        super().__init__(param, options=options, x_average=x_average)
         pybamm.citations.register("Timms2021")
 
     def get_fundamental_variables(self):

diff --git a/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_1D_current_collectors.py b/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_1D_current_collectors.py
@@ -22,8 +22,8 @@ class CurrentCollector1D(BaseThermal):
 
     """
 
-    def __init__(self, param, options=None):
-        super().__init__(param, options=options)
+    def __init__(self, param, options=None, x_average=True):
+        super().__init__(param, options=options, x_average=x_average)
         pybamm.citations.register("Timms2021")
 
     def get_fundamental_variables(self):

diff --git a/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_2D_current_collectors.py b/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_2D_current_collectors.py
@@ -22,8 +22,8 @@ class CurrentCollector2D(BaseThermal):
 
     """
 
-    def __init__(self, param, options=None):
-        super().__init__(param, options=options)
+    def __init__(self, param, options=None, x_average=True):
+        super().__init__(param, options=options, x_average=x_average)
         pybamm.citations.register("Timms2021")
 
     def get_fundamental_variables(self):

diff --git a/pybamm/models/submodels/thermal/pouch_cell/x_full.py b/pybamm/models/submodels/thermal/pouch_cell/x_full.py
@@ -24,8 +24,8 @@ class OneDimensionalX(BaseThermal):
 
     """
 
-    def __init__(self, param, options=None):
-        super().__init__(param, options=options)
+    def __init__(self, param, options=None, x_average=False):
+        super().__init__(param, options=options, x_average=x_average)
         pybamm.citations.register("Timms2021")
 
     def get_fundamental_variables(self):