Merge pull request #1617 from pybamm-team/issue-1505-porosity

Issue 1505 porosity

CHANGELOG.md

@@ -36,6 +36,7 @@ example notebook ([#1602](https://github.com/pybamm-team/PyBaMM/pull/1602))
 
 ## Bug fixes
 
+-   Porosity variation for SEI and plating models is calculated from the film thickness rather than from a separate ODE ([#1617](https://github.com/pybamm-team/PyBaMM/pull/1617))
 -   Fixed a bug where the order of the indexing for the entries of variables discretised using FEM was incorrect ([#1556](https://github.com/pybamm-team/PyBaMM/pull/1556))
 -   Fix broken module import for spyder when running a script twice ([#1555](https://github.com/pybamm-team/PyBaMM/pull/1555))
 -   Fixed ElectrodeSOH model for multi-dimensional simulations ([#1548](https://github.com/pybamm-team/PyBaMM/pull/1548))

docs/source/models/submodels/porosity/index.rst

@@ -7,4 +7,5 @@ Porosity
   base_porosity
   constant_porosity
   reaction_driven_porosity
+  reaction_driven_porosity_ode
 

docs/source/models/submodels/porosity/reaction_driven_porosity_ode.rst

@@ -0,0 +1,8 @@
+Reaction-driven Model as an ODE
+===============================
+
+.. autoclass:: pybamm.porosity.ReactionDrivenODE
+    :members:
+
+
+

examples/scripts/calendar_ageing.py

@@ -35,6 +35,7 @@
 
     sim.solve(t_eval=t_eval, solver=solver)
     sims.append(sim)
+
 pb.dynamic_plot(
     sims,
     [

examples/scripts/cycling_ageing.py

@@ -30,24 +30,28 @@
             "Hold at 4.2 V until C/20",
             "Rest for 30 minutes",
             "Discharge at C/3 until 2.8 V",
+            "Rest for 30 minutes",
         ),
         (
             "Charge at 1 C until 4.2 V",
             "Hold at 4.2 V until C/20",
             "Rest for 30 minutes",
             "Discharge at 1 C until 2.8 V",
+            "Rest for 30 minutes",
         ),
         (
             "Charge at 1 C until 4.2 V",
             "Hold at 4.2 V until C/20",
             "Rest for 30 minutes",
             "Discharge at 2 C until 2.8 V",
+            "Rest for 30 minutes",
         ),
         (
             "Charge at 1 C until 4.2 V",
             "Hold at 4.2 V until C/20",
             "Rest for 30 minutes",
-            "Discharge at 3 C until 2.8 V",
+            "Discharge at 3 C until 2.8 V (10 second period)",
+            "Rest for 30 minutes",
         ),
     ]
 )

pybamm/models/full_battery_models/lead_acid/full.py

@@ -54,7 +54,7 @@ def __init__(self, options=None, name="Full model", build=True):
         pybamm.citations.register("Sulzer2019physical")
 
     def set_porosity_submodel(self):
-        self.submodels["porosity"] = pybamm.porosity.ReactionDriven(
+        self.submodels["porosity"] = pybamm.porosity.ReactionDrivenODE(
             self.param, self.options, False
         )
 

pybamm/models/full_battery_models/lead_acid/higher_order.py

@@ -198,7 +198,7 @@ def set_full_porosity_submodel(self):
         Update porosity submodel, now that we have the spatially heterogeneous
         interfacial current densities
         """
-        self.submodels["full porosity"] = pybamm.porosity.ReactionDriven(
+        self.submodels["full porosity"] = pybamm.porosity.ReactionDrivenODE(
             self.param, self.options, False
         )
 
@@ -271,7 +271,7 @@ def set_full_porosity_submodel(self):
         Update porosity submodel, now that we have the spatially heterogeneous
         interfacial current densities
         """
-        self.submodels["full porosity"] = pybamm.porosity.ReactionDriven(
+        self.submodels["full porosity"] = pybamm.porosity.ReactionDrivenODE(
             self.param, self.options, False
         )
 

pybamm/models/full_battery_models/lead_acid/loqs.py

@@ -95,7 +95,7 @@ def set_current_collector_submodel(self):
 
     def set_porosity_submodel(self):
 
-        self.submodels["leading-order porosity"] = pybamm.porosity.ReactionDriven(
+        self.submodels["leading-order porosity"] = pybamm.porosity.ReactionDrivenODE(
             self.param, self.options, True
         )
 

pybamm/models/submodels/interface/lithium_plating/base_plating.py

@@ -70,6 +70,7 @@ def _get_standard_concentration_variables(self, c_plated_Li):
             f"X-averaged {domain} lithium plating concentration": c_plated_Li_av,
             f"X-averaged {domain} lithium plating concentration"
             " [mol.m-3]": c_plated_Li_av * c_scale,
+            f"{Domain} lithium plating thickness": L_plated_Li,
             f"{Domain} lithium plating thickness [m]": L_plated_Li * L_scale,
             f"X-averaged {domain} lithium plating thickness [m]": L_plated_Li_av
             * L_scale,

pybamm/models/submodels/interface/sei/base_sei.py

@@ -55,6 +55,8 @@ def _get_standard_thickness_variables(self, L_inner, L_outer):
 
         L_inner_av = pybamm.x_average(L_inner)
         L_outer_av = pybamm.x_average(L_outer)
+        L_tot = L_inner + L_outer
+        L_tot_av = L_inner_av + L_outer_av
 
         variables = {
             "Inner " + domain + " SEI thickness": L_inner,
@@ -65,6 +67,10 @@ def _get_standard_thickness_variables(self, L_inner, L_outer):
             "Outer " + domain + " SEI thickness [m]": L_outer * L_scale,
             "X-averaged outer " + domain + " SEI thickness": L_outer_av,
             "X-averaged outer " + domain + " SEI thickness [m]": L_outer_av * L_scale,
+            self.domain + " electrode SEI thickness": L_tot,
+            self.domain + " electrode SEI thickness [m]": L_tot * L_scale,
+            "X-averaged " + domain + " SEI thickness": L_tot_av,
+            "X-averaged " + domain + " SEI thickness [m]": L_tot_av * L_scale,
         }
 
         # Get variables related to the total thickness

pybamm/models/submodels/porosity/__init__.py

@@ -1,3 +1,4 @@
 from .base_porosity import BaseModel
 from .constant_porosity import Constant
 from .reaction_driven_porosity import ReactionDriven
+from .reaction_driven_porosity_ode import ReactionDrivenODE

pybamm/models/submodels/porosity/reaction_driven_porosity.py

@@ -1,12 +1,12 @@
 #
-# Class for reaction driven porosity changes
+# Class for reaction driven porosity changes as a multiple of SEI/plating thicknesses
 #
 import pybamm
 from .base_porosity import BaseModel
 
 
 class ReactionDriven(BaseModel):
-    """Reaction-driven porosity changes
+    """Reaction-driven porosity changes as a multiple of SEI/plating thicknesses
 
     Parameters
     ----------
@@ -24,87 +24,28 @@ def __init__(self, param, options, x_average):
         super().__init__(param, options)
         self.x_average = x_average
 
-    def get_fundamental_variables(self):
-        if self.x_average is True:
-            eps_n_pc = pybamm.standard_variables.eps_n_pc
-            eps_s_pc = pybamm.standard_variables.eps_s_pc
-            eps_p_pc = pybamm.standard_variables.eps_p_pc
-
-            eps_n = pybamm.PrimaryBroadcast(eps_n_pc, "negative electrode")
-            eps_s = pybamm.PrimaryBroadcast(eps_s_pc, "separator")
-            eps_p = pybamm.PrimaryBroadcast(eps_p_pc, "positive electrode")
-        else:
-            eps_n = pybamm.standard_variables.eps_n
-            eps_s = pybamm.standard_variables.eps_s
-            eps_p = pybamm.standard_variables.eps_p
-        variables = self._get_standard_porosity_variables(eps_n, eps_s, eps_p)
-
-        return variables
-
     def get_coupled_variables(self, variables):
+        L_sei_n = variables["Total negative electrode SEI thickness [m]"]
+        L_sei_0 = self.param.L_inner_0_dim + self.param.L_outer_0_dim
+        L_pl_n = variables["Negative electrode lithium plating thickness [m]"]
 
-        if self.x_average is True:
-            j_n = variables["X-averaged negative electrode interfacial current density"]
-            j_p = variables["X-averaged positive electrode interfacial current density"]
-            j_sei_n = variables[
-                "X-averaged negative electrode SEI interfacial current density"
-            ]
-            j_plating = variables[
-                "X-averaged negative electrode lithium plating "
-                "interfacial current density"
-            ]
-            deps_s_dt = pybamm.PrimaryBroadcast(0, "current collector")
-        else:
-            j_n = variables["Negative electrode interfacial current density"]
-            j_p = variables["Positive electrode interfacial current density"]
-            j_sei_n = variables["Negative electrode SEI interfacial current density"]
-            j_plating = variables[
-                "Negative electrode lithium plating interfacial current density"
-            ]
-            deps_s_dt = pybamm.FullBroadcast(
-                0, "separator", auxiliary_domains={"secondary": "current collector"}
-            )
+        L_tot = (L_sei_n - L_sei_0) + L_pl_n
 
-        deps_n_dt = -self.param.beta_surf_n * j_n
-        if self.options["SEI porosity change"] == "true":
-            beta_sei_n = self.param.beta_sei_n
-            deps_n_dt += beta_sei_n * j_sei_n
-        if self.options["lithium plating porosity change"] == "true":
-            beta_plating = self.param.beta_plating
-            deps_n_dt += beta_plating * j_plating
+        a_n = variables["Negative electrode surface area to volume ratio [m-1]"]
 
-        deps_p_dt = -self.param.beta_surf_p * j_p
+        # This assumes a thin film so curvature effects are neglected. They could be
+        # included (i.e. for a sphere it is a_n * (L_tot + L_tot ** 2 / R_n + L_tot **
+        # 3 / (3 * R_n ** 2))) but it is not clear if it is relevant or not.
+        delta_eps_n = -a_n * L_tot
 
-        variables.update(
-            self._get_standard_porosity_change_variables(
-                deps_n_dt, deps_s_dt, deps_p_dt
-            )
+        eps_n = self.param.epsilon_n_init + delta_eps_n
+        eps_s = pybamm.FullBroadcast(
+            self.param.epsilon_s_init, "separator", "current collector"
         )
+        # no SEI or plating in the positive electrode
+        eps_p = pybamm.FullBroadcast(
+            self.param.epsilon_p_init, "positive electrode", "current collector"
+        )
+        variables = self._get_standard_porosity_variables(eps_n, eps_s, eps_p)
 
         return variables
-
-    def set_rhs(self, variables):
-        if self.x_average is True:
-            for domain in ["negative electrode", "separator", "positive electrode"]:
-                eps_av = variables["X-averaged " + domain + " porosity"]
-                deps_dt_av = variables["X-averaged " + domain + " porosity change"]
-                self.rhs.update({eps_av: deps_dt_av})
-        else:
-            eps = variables["Porosity"]
-            deps_dt = variables["Porosity change"]
-            self.rhs = {eps: deps_dt}
-
-    def set_initial_conditions(self, variables):
-        if self.x_average is True:
-            eps_n_av = variables["X-averaged negative electrode porosity"]
-            eps_s_av = variables["X-averaged separator porosity"]
-            eps_p_av = variables["X-averaged positive electrode porosity"]
-
-            self.initial_conditions = {
-                eps_n_av: self.param.epsilon_n_init,
-                eps_s_av: self.param.epsilon_s_init,
-                eps_p_av: self.param.epsilon_p_init,
-            }
-        else:
-            eps = variables["Porosity"]
-            self.initial_conditions = {eps: self.param.epsilon_init}

pybamm/models/submodels/porosity/reaction_driven_porosity_ode.py

@@ -0,0 +1,92 @@
+#
+# Class for reaction driven porosity changes as an ODE
+#
+import pybamm
+from .base_porosity import BaseModel
+
+
+class ReactionDrivenODE(BaseModel):
+    """Reaction-driven porosity changes as an ODE
+
+    Parameters
+    ----------
+    param : parameter class
+        The parameters to use for this submodel
+    options : dict
+        Options dictionary passed from the full model
+    x_average : bool
+        Whether to use x-averaged variables (SPM, SPMe, etc) or full variables (DFN)
+
+    **Extends:** :class:`pybamm.porosity.BaseModel`
+    """
+
+    def __init__(self, param, options, x_average):
+        super().__init__(param, options)
+        self.x_average = x_average
+
+    def get_fundamental_variables(self):
+        if self.x_average is True:
+            eps_n_pc = pybamm.standard_variables.eps_n_pc
+            eps_s_pc = pybamm.standard_variables.eps_s_pc
+            eps_p_pc = pybamm.standard_variables.eps_p_pc
+
+            eps_n = pybamm.PrimaryBroadcast(eps_n_pc, "negative electrode")
+            eps_s = pybamm.PrimaryBroadcast(eps_s_pc, "separator")
+            eps_p = pybamm.PrimaryBroadcast(eps_p_pc, "positive electrode")
+        else:
+            eps_n = pybamm.standard_variables.eps_n
+            eps_s = pybamm.standard_variables.eps_s
+            eps_p = pybamm.standard_variables.eps_p
+        variables = self._get_standard_porosity_variables(eps_n, eps_s, eps_p)
+
+        return variables
+
+    def get_coupled_variables(self, variables):
+
+        if self.x_average is True:
+            j_n = variables["X-averaged negative electrode interfacial current density"]
+            j_p = variables["X-averaged positive electrode interfacial current density"]
+            deps_s_dt = pybamm.PrimaryBroadcast(0, "current collector")
+        else:
+            j_n = variables["Negative electrode interfacial current density"]
+            j_p = variables["Positive electrode interfacial current density"]
+            deps_s_dt = pybamm.FullBroadcast(
+                0, "separator", auxiliary_domains={"secondary": "current collector"}
+            )
+
+        deps_n_dt = -self.param.beta_surf_n * j_n
+        deps_p_dt = -self.param.beta_surf_p * j_p
+
+        variables.update(
+            self._get_standard_porosity_change_variables(
+                deps_n_dt, deps_s_dt, deps_p_dt
+            )
+        )
+
+        return variables
+
+    def set_rhs(self, variables):
+        if self.x_average is True:
+            for domain in ["negative electrode", "separator", "positive electrode"]:
+                eps_av = variables["X-averaged " + domain + " porosity"]
+                deps_dt_av = variables["X-averaged " + domain + " porosity change"]
+                self.rhs.update({eps_av: deps_dt_av})
+        else:
+            eps = variables["Porosity"]
+            deps_dt = variables["Porosity change"]
+            self.rhs = {eps: deps_dt}
+
+    def set_initial_conditions(self, variables):
+        if self.x_average is True:
+            eps_n_av = variables["X-averaged negative electrode porosity"]
+            eps_s_av = variables["X-averaged separator porosity"]
+            eps_p_av = variables["X-averaged positive electrode porosity"]
+
+            self.initial_conditions = {
+                eps_n_av: self.param.epsilon_n_init,
+                eps_s_av: self.param.epsilon_s_init,
+                eps_p_av: self.param.epsilon_p_init,
+            }
+        else:
+            eps = variables["Porosity"]
+            self.initial_conditions = {eps: self.param.epsilon_init}

tests/unit/test_models/test_submodels/test_porosity/test_full_reaction_driven_porosity.py

@@ -12,13 +12,12 @@ def test_public_functions(self):
         # param = pybamm.LeadAcidParameters()
         param = pybamm.LithiumIonParameters()
         a_n = pybamm.PrimaryBroadcast(pybamm.Scalar(0), ["negative electrode"])
-        a_p = pybamm.PrimaryBroadcast(pybamm.Scalar(0), ["positive electrode"])
 
         variables = {
-            "Negative electrode interfacial current density": a_n,
-            "Negative electrode SEI interfacial current density": a_n,
-            "Positive electrode interfacial current density": a_p,
-            "Negative electrode lithium plating interfacial current density": a_n,
+            "Total negative electrode SEI thickness [m]": a_n,
+            "Negative electrode lithium plating thickness [m]": a_n,
+            "Negative electrode surface area to volume ratio [m-1]": a_n,
+
         }
         options = {
             "SEI": "ec reaction limited",

tests/unit/test_models/test_submodels/test_porosity/test_leading_reaction_driven_porosity.py

@@ -13,11 +13,9 @@ def test_public_functions(self):
         param = pybamm.LithiumIonParameters()
         a = pybamm.PrimaryBroadcast(pybamm.Scalar(0), "current collector")
         variables = {
-            "X-averaged negative electrode interfacial current density": a,
-            "X-averaged negative electrode SEI interfacial current density": a,
-            "X-averaged positive electrode interfacial current density": a,
-            "X-averaged negative electrode lithium plating "
-            "interfacial current density": a,
+            "Total negative electrode SEI thickness [m]": a,
+            "Negative electrode lithium plating thickness [m]": a,
+            "Negative electrode surface area to volume ratio [m-1]": a,
         }
         options = {
             "SEI": "ec reaction limited",