Merge pull request #2529 from pybamm-team/porosity-times-concentration

Porosity times concentration

CHANGELOG.md

@@ -2,15 +2,18 @@
 
 ## Features
 
+- Added variables "Loss of lithium due to loss of active material in negative/positive electrode [mol]". These should be included in the calculation of "total lithium in system" to make sure that lithium is truly conserved. ([#2529](https://github.com/pybamm-team/PyBaMM/pull/2529))
 - `initial_soc` can now be a string "x V", in which case the simulation is initialized to start from that voltage ([#2508](https://github.com/pybamm-team/PyBaMM/pull/2508))
 - The `ElectrodeSOH` solver can now calculate electrode balance based on a target "cell capacity" (requires cell capacity "Q" as input), as well as the default "cyclable cell capacity" (requires cyclable lithium capacity "Q_Li" as input). Use the keyword argument `known_value` to control which is used. ([#2508](https://github.com/pybamm-team/PyBaMM/pull/2508))
 
 ## Bug fixes
 
+- Fixed "constant concentration" electrolyte model so that "porosity times concentration" is conserved when porosity changes ([#2529](https://github.com/pybamm-team/PyBaMM/pull/2529))
 - Fix installation on `Google Colab` (`pybtex` and `Colab` issue) ([#2526](https://github.com/pybamm-team/PyBaMM/pull/2526))
 
 ## Breaking changes
 
+- Renamed "Negative/Positive electrode SOC" to "Negative/Positive electrode stoichiometry" to avoid confusion with cell SOC ([#2529](https://github.com/pybamm-team/PyBaMM/pull/2529))
 - Removed external variables and submodels. InputParameter should now be used in all cases ([#2502](https://github.com/pybamm-team/PyBaMM/pull/2502))
 - Trying to use a solver to solve multiple models results in a RuntimeError exception ([#2481](https://github.com/pybamm-team/PyBaMM/pull/2481))
 - Inputs for the `ElectrodeSOH` solver are now (i) "Q_Li", the total cyclable capacity of lithium in the electrodes (previously "n_Li", the total number of moles, n_Li = 3600/F \* Q_Li) (ii) "Q_n", the capacity of the negative electrode (previously "C_n"), and "Q_p", the capacity of the positive electrode (previously "C_p") ([#2508](https://github.com/pybamm-team/PyBaMM/pull/2508))

examples/notebooks/models/electrode-state-of-health.ipynb

@@ -88,8 +88,8 @@
                 "spm_sol.plot([\n",
                 "    \"Terminal voltage [V]\", \n",
                 "    \"Current [A]\", \n",
-                "    \"Negative electrode SOC\",\n",
-                "    \"Positive electrode SOC\",\n",
+                "    \"Negative electrode stoichiometry\",\n",
+                "    \"Positive electrode stoichiometry\",\n",
                 "])"
             ]
         },
@@ -290,8 +290,8 @@
             ],
             "source": [
                 "t = spm_sol[\"Time [h]\"].data\n",
-                "x_spm = spm_sol[\"Negative electrode SOC\"].data\n",
-                "y_spm = spm_sol[\"Positive electrode SOC\"].data\n",
+                "x_spm = spm_sol[\"Negative electrode stoichiometry\"].data\n",
+                "y_spm = spm_sol[\"Positive electrode stoichiometry\"].data\n",
                 "\n",
                 "x_0 = esoh_sol[\"x_0\"].data * np.ones_like(t)\n",
                 "y_0 = esoh_sol[\"y_0\"].data * np.ones_like(t)\n",
@@ -635,4 +635,4 @@
     },
     "nbformat": 4,
     "nbformat_minor": 4
-}
+}

pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py

@@ -178,17 +178,19 @@ def set_degradation_variables(self):
 
         # Lithium lost to side reactions
         # Different way of measuring LLI but should give same value
-        LLI_sei = self.variables["Loss of lithium to SEI [mol]"]
-        LLI_reactions = LLI_sei
+        n_Li_lost_sei = self.variables["Loss of lithium to SEI [mol]"]
+        n_Li_lost_reactions = n_Li_lost_sei
         if "negative electrode" in domains:
-            LLI_sei_cracks = self.variables["Loss of lithium to SEI on cracks [mol]"]
-            LLI_pl = self.variables["Loss of lithium to lithium plating [mol]"]
-            LLI_reactions += LLI_sei_cracks + LLI_pl
+            n_Li_lost_sei_cracks = self.variables[
+                "Loss of lithium to SEI on cracks [mol]"
+            ]
+            n_Li_lost_pl = self.variables["Loss of lithium to lithium plating [mol]"]
+            n_Li_lost_reactions += n_Li_lost_sei_cracks + n_Li_lost_pl
 
         self.variables.update(
             {
-                "Total lithium lost to side reactions [mol]": LLI_reactions,
-                "Total capacity lost to side reactions [A.h]": LLI_reactions
+                "Total lithium lost to side reactions [mol]": n_Li_lost_reactions,
+                "Total capacity lost to side reactions [A.h]": n_Li_lost_reactions
                 * param.F
                 / 3600,
             }

pybamm/models/submodels/active_material/constant_active_material.py

@@ -35,4 +35,11 @@ def get_fundamental_variables(self):
             self._get_standard_active_material_change_variables(deps_solid_dt)
         )
 
+        variables.update(
+            {
+                "Loss of lithium due to loss of active material "
+                f"in {domain} electrode [mol]": pybamm.Scalar(0)
+            }
+        )
+
         return variables

pybamm/models/submodels/active_material/loss_active_material.py

@@ -53,6 +53,16 @@ def get_fundamental_variables(self):
                 auxiliary_domains={"secondary": "current collector"},
             )
         variables = self._get_standard_active_material_variables(eps_solid)
+        lli_due_to_lam = pybamm.Variable(
+            "Loss of lithium due to loss of active material "
+            f"in {domain} electrode [mol]"
+        )
+        variables.update(
+            {
+                "Loss of lithium due to loss of active material "
+                f"in {domain} electrode [mol]": lli_due_to_lam
+            }
+        )
         return variables
 
     def get_coupled_variables(self, variables):
@@ -133,7 +143,23 @@ def set_rhs(self, variables):
                 f"{Domain} electrode active material volume fraction change"
             ]
 
-        self.rhs = {eps_solid: deps_solid_dt}
+        # Loss of lithium due to loss of active material
+        # See eq 37 in "Sulzer, Valentin, et al. "Accelerated battery lifetime
+        # simulations using adaptive inter-cycle extrapolation algorithm."
+        # Journal of The Electrochemical Society 168.12 (2021): 120531.
+        lli_due_to_lam = variables[
+            "Loss of lithium due to loss of active material "
+            f"in {domain} electrode [mol]"
+        ]
+        # Multiply by mol.m-3 * m3 to get mol
+        c_s_av = variables[f"Average {domain} particle concentration [mol.m-3]"]
+        V = self.domain_param.L * self.param.A_cc
+
+        self.rhs = {
+            # minus sign because eps_solid is decreasing and LLI measures positive
+            lli_due_to_lam: -c_s_av * V * pybamm.x_average(deps_solid_dt),
+            eps_solid: deps_solid_dt,
+        }
 
     def set_initial_conditions(self, variables):
         domain, Domain = self.domain_Domain
@@ -148,3 +174,9 @@ def set_initial_conditions(self, variables):
         else:
             eps_solid = variables[f"{Domain} electrode active material volume fraction"]
             self.initial_conditions = {eps_solid: eps_solid_init}
+
+        lli_due_to_lam = variables[
+            "Loss of lithium due to loss of active material "
+            f"in {domain} electrode [mol]"
+        ]
+        self.initial_conditions[lli_due_to_lam] = pybamm.Scalar(0)

pybamm/models/submodels/electrolyte_diffusion/constant_concentration.py

@@ -23,12 +23,13 @@ def __init__(self, param, options=None):
         super().__init__(param, options)
 
     def get_fundamental_variables(self):
-        c_e_dict = {
-            domain: pybamm.FullBroadcast(1, domain, "current collector")
+        eps_c_e_dict = {
+            domain: self.param.domain_params[domain.split()[0]].epsilon_init * 1
             for domain in self.options.whole_cell_domains
         }
-        variables = self._get_standard_concentration_variables(c_e_dict)
-
+        variables = self._get_standard_porosity_times_concentration_variables(
+            eps_c_e_dict
+        )
         N_e = pybamm.FullBroadcastToEdges(
             0,
             [domain for domain in self.options.whole_cell_domains],
@@ -40,15 +41,21 @@ def get_fundamental_variables(self):
         return variables
 
     def get_coupled_variables(self, variables):
-        eps_c_e_dict = {}
+        c_e_dict = {}
         for domain in self.options.whole_cell_domains:
             Domain = domain.capitalize()
             eps_k = variables[f"{Domain} porosity"]
-            c_e_k = variables[f"{Domain.split()[0]} electrolyte concentration"]
-            eps_c_e_dict[domain] = eps_k * c_e_k
-        variables.update(
-            self._get_standard_porosity_times_concentration_variables(eps_c_e_dict)
+            eps_c_e_k = variables[f"{Domain} porosity times concentration"]
+            c_e_k = eps_c_e_k / eps_k
+            c_e_dict[domain] = c_e_k
+
+        variables["Electrolyte concentration concatenation"] = pybamm.concatenation(
+            *c_e_dict.values()
         )
+        variables.update(self._get_standard_domain_concentration_variables(c_e_dict))
+
+        c_e = variables["Porosity times concentration"] / variables["Porosity"]
+        variables.update(self._get_standard_whole_cell_concentration_variables(c_e))
 
         return variables
 

pybamm/models/submodels/particle/base_particle.py

@@ -139,7 +139,7 @@ def _get_total_concentration_variables(self, variables):
 
         variables.update(
             {
-                f"{Domain} electrode {phase_name}SOC": c_s_vol_av,
+                f"{Domain} electrode {phase_name}stoichiometry": c_s_vol_av,
                 f"{Domain} electrode {phase_name}volume-averaged "
                 "concentration": c_s_vol_av,
                 f"{Domain} electrode {phase_name}volume-averaged "

pybamm/solvers/processed_variable.py

@@ -6,6 +6,7 @@
 import numpy as np
 import pybamm
 import scipy.interpolate as interp
+from scipy.integrate import cumulative_trapezoid
 
 
 class ProcessedVariable(object):
@@ -61,6 +62,7 @@ def __init__(
 
         # Set timescale
         self.timescale = solution.timescale_eval
+        self.t_pts_nondim = solution.t
         self.t_pts = solution.t * self.timescale
 
         # Store length scales
@@ -114,30 +116,24 @@ def initialise_0D(self):
         # initialise empty array of the correct size
         entries = np.empty(len(self.t_pts))
         idx = 0
-        last_t = 0
+
+        entries = np.empty(len(self.t_pts))
+        idx = 0
         # Evaluate the base_variable index-by-index
         for ts, ys, inputs, base_var_casadi in zip(
             self.all_ts, self.all_ys, self.all_inputs_casadi, self.base_variables_casadi
         ):
             for inner_idx, t in enumerate(ts):
                 t = ts[inner_idx]
                 y = ys[:, inner_idx]
-                if self.cumtrapz_ic is not None:
-                    if idx == 0:
-                        new_val = t * base_var_casadi(t, y, inputs).full()[0, 0]
-                        entries[idx] = self.cumtrapz_ic + (
-                            t * base_var_casadi(t, y, inputs).full()[0, 0]
-                        )
-                    else:
-                        new_val = (t - last_t) * (
-                            base_var_casadi(t, y, inputs).full()[0, 0]
-                        )
-                        entries[idx] = new_val + entries[idx - 1]
-                else:
-                    entries[idx] = base_var_casadi(t, y, inputs).full()[0, 0]
+                entries[idx] = float(base_var_casadi(t, y, inputs))
 
                 idx += 1
-                last_t = t
+
+        if self.cumtrapz_ic is not None:
+            entries = cumulative_trapezoid(
+                entries, self.t_pts_nondim, initial=float(self.cumtrapz_ic)
+            )
 
         # set up interpolation
         if len(self.t_pts) == 1:

pybamm/util.py

@@ -57,6 +57,13 @@ def __getitem__(self, key):
         try:
             return super().__getitem__(key)
         except KeyError:
+            if key in ["Negative electrode SOC", "Positive electrode SOC"]:
+                domain = key.split(" ")[0]
+                raise KeyError(
+                    f"Variable '{domain} electrode SOC' has been renamed to "
+                    f"'{domain} electrode stoichiometry' to avoid confusion "
+                    "with cell SOC"
+                )
             best_matches = self.get_best_matches(key)
             raise KeyError(f"'{key}' not found. Best matches are {best_matches}")
 

tests/integration/test_models/standard_output_tests.py

@@ -293,8 +293,13 @@ def __init__(self, model, param, disc, solution, operating_condition):
         self.N_s_n = solution[f"Negative {self.phase_name_n}particle flux"]
         self.N_s_p = solution[f"Positive {self.phase_name_p}particle flux"]
 
-        self.c_SEI_tot = solution["Loss of lithium to SEI [mol]"]
-        self.c_Li_tot = solution["Loss of lithium to lithium plating [mol]"]
+        self.n_Li_side = solution["Total lithium lost to side reactions [mol]"]
+        self.n_Li_LAM_n = solution[
+            "Loss of lithium due to loss of active material in negative electrode [mol]"
+        ]
+        self.n_Li_LAM_p = solution[
+            "Loss of lithium due to loss of active material in positive electrode [mol]"
+        ]
 
         if model.options["particle size"] == "distribution":
             # These concentration variables are only present for distribution models.
@@ -404,8 +409,9 @@ def test_conservation(self):
         c_s_tot = (
             self.c_s_n_tot(self.solution.t)
             + self.c_s_p_tot(self.solution.t)
-            + self.c_SEI_tot(self.solution.t)
-            + self.c_Li_tot(self.solution.t)
+            + self.n_Li_side(self.solution.t)
+            + self.n_Li_LAM_n(self.solution.t)
+            + self.n_Li_LAM_p(self.solution.t)
         )
         diff = (c_s_tot[1:] - c_s_tot[:-1]) / c_s_tot[:-1]
         if self.model.options["particle"] == "quartic profile":
@@ -800,20 +806,32 @@ def __init__(self, model, param, disc, solution, operating_condition):
         self.LLI = solution["Loss of lithium inventory [%]"]
         self.n_Li_lost = solution["Total lithium lost [mol]"]
         self.n_Li_lost_rxn = solution["Total lithium lost to side reactions [mol]"]
+        self.n_Li_lost_LAM_n = solution[
+            "Loss of lithium due to loss of active material in negative electrode [mol]"
+        ]
+        self.n_Li_lost_LAM_p = solution[
+            "Loss of lithium due to loss of active material in positive electrode [mol]"
+        ]
 
     def test_degradation_modes(self):
         """Test degradation modes are between 0 and 100%"""
         np.testing.assert_array_less(-3e-3, self.LLI(self.t))
         np.testing.assert_array_less(-1e-13, self.LAM_ne(self.t))
         np.testing.assert_array_less(-1e-13, self.LAM_pe(self.t))
+        np.testing.assert_array_less(-1e-13, self.n_Li_lost_LAM_n(self.t))
+        np.testing.assert_array_less(-1e-13, self.n_Li_lost_LAM_p(self.t))
         np.testing.assert_array_less(self.LLI(self.t), 100)
         np.testing.assert_array_less(self.LAM_ne(self.t), 100)
         np.testing.assert_array_less(self.LAM_pe(self.t), 100)
 
     def test_lithium_lost(self):
         """Test the two ways of measuring lithium lost give the same value"""
         np.testing.assert_array_almost_equal(
-            self.n_Li_lost(self.t), self.n_Li_lost_rxn(self.t), decimal=3
+            self.n_Li_lost(self.t),
+            self.n_Li_lost_rxn(self.t)
+            + self.n_Li_lost_LAM_n(self.t)
+            + self.n_Li_lost_LAM_p(self.t),
+            decimal=5,
         )
 
     def test_all(self):

tests/integration/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py

@@ -209,12 +209,12 @@ def test_sei_asymmetric_ec_reaction_limited(self):
         )
         self.run_basic_processing_test(options, parameter_values=parameter_values)
 
-    def test_loss_active_material_stress_negative(self):
+    def test_loss_active_material_stress_positive(self):
         options = {"loss of active material": ("none", "stress-driven")}
         parameter_values = pybamm.ParameterValues("Ai2020")
         self.run_basic_processing_test(options, parameter_values=parameter_values)
 
-    def test_loss_active_material_stress_positive(self):
+    def test_loss_active_material_stress_negative(self):
         options = {"loss of active material": ("stress-driven", "none")}
         parameter_values = pybamm.ParameterValues("Ai2020")
         self.run_basic_processing_test(options, parameter_values=parameter_values)

tests/unit/test_util.py

@@ -69,6 +69,9 @@ def test_fuzzy_dict(self):
         with self.assertRaisesRegex(KeyError, "'test3' not found. Best matches are "):
             d.__getitem__("test3")
 
+        with self.assertRaisesRegex(KeyError, "stoichiometry"):
+            d.__getitem__("Negative electrode SOC")
+
     def test_get_parameters_filepath(self):
         tempfile_obj = tempfile.NamedTemporaryFile("w", dir=".")
         self.assertTrue(