Composite porosity

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 - Added two more submodels (options) for the SEI: Lars von Kolzenberg (2020) model and Tunneling Limit model ([#4394](https://github.com/pybamm-team/PyBaMM/pull/4394))
+- Porosity change now works for composite electrode ([#4417](https://github.com/pybamm-team/PyBaMM/pull/4417))
 
 ## Breaking changes
 

src/pybamm/models/submodels/porosity/reaction_driven_porosity.py

@@ -25,36 +25,59 @@ def __init__(self, param, options, x_average):
     def get_coupled_variables(self, variables):
         eps_dict = {}
         for domain in self.options.whole_cell_domains:
-            if domain == "separator":
-                delta_eps_k = 0  # separator porosity does not change
-            else:
-                Domain = domain.split()[0].capitalize()
-                L_sei_k = variables[f"{Domain} total SEI thickness [m]"]
+            delta_eps_k = 0
+            if domain != "separator":  # separator porosity does not change
+                dom = domain.split()[0]
+                Domain = dom.capitalize()
                 if Domain == "Negative":
                     L_sei_0 = self.param.n.prim.L_sei_0
                 elif Domain == "Positive":
                     L_sei_0 = self.param.p.prim.L_sei_0
-                L_pl_k = variables[f"{Domain} lithium plating thickness [m]"]
-                L_dead_k = variables[f"{Domain} dead lithium thickness [m]"]
-                L_sei_cr_k = variables[f"{Domain} total SEI on cracks thickness [m]"]
                 roughness_k = variables[f"{Domain} electrode roughness ratio"]
+                SEI_option = getattr(self.options, dom)["SEI"]
+                phases_option = getattr(self.options, dom)["particle phases"]
+                phases = self.options.phases[dom]
+                for phase in phases:
+                    if phases_option == "1" and phase == "primary":
+                        # `domain` has one phase
+                        phase_name = ""
+                        pref = ""
+                    else:
+                        # `domain` has more than one phase
+                        phase_name = phase + " "
+                        pref = phase.capitalize() + ": "
+                    L_sei_k = variables[f"{Domain} total {phase_name}SEI thickness [m]"]
+                    if SEI_option == "none":
+                        L_sei_0 = pybamm.Scalar(0)
+                    else:
+                        L_sei_0 = pybamm.Parameter(f"{pref}Initial SEI thickness [m]")
+                    L_pl_k = variables[
+                        f"{Domain} {phase_name}lithium plating thickness [m]"
+                    ]
+                    L_dead_k = variables[
+                        f"{Domain} {phase_name}dead lithium thickness [m]"
+                    ]
+                    L_sei_cr_k = variables[
+                        f"{Domain} total {phase_name}SEI on cracks thickness [m]"
+                    ]
 
-                L_tot = (
-                    (L_sei_k - L_sei_0)
-                    + L_pl_k
-                    + L_dead_k
-                    + L_sei_cr_k * (roughness_k - 1)
-                )
+                    L_tot = (
+                        (L_sei_k - L_sei_0)
+                        + L_pl_k
+                        + L_dead_k
+                        + L_sei_cr_k * (roughness_k - 1)
+                    )
 
-                a_k = variables[
-                    f"{Domain} electrode surface area to volume ratio [m-1]"
-                ]
+                    a_k = variables[
+                        f"{Domain} electrode {phase_name}"
+                        "surface area to volume ratio [m-1]"
+                    ]
 
-                # This assumes a thin film so curvature effects are neglected.
-                # They could be included (e.g. 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_k = -a_k * L_tot
+                    # This assumes a thin film so curvature effects are neglected.
+                    # They could be included (e.g. 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_k += -a_k * L_tot
 
             domain_param = self.param.domain_params[domain.split()[0]]
             eps_k = domain_param.epsilon_init + delta_eps_k

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

@@ -356,6 +356,7 @@ def test_composite_graphite_silicon_sei(self):
             "particle phases": ("2", "1"),
             "open-circuit potential": (("single", "current sigmoid"), "single"),
             "SEI": "ec reaction limited",
+            "SEI porosity change": "true",
         }
         parameter_values = pybamm.ParameterValues("Chen2020_composite")
         name = "Negative electrode active material volume fraction"

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

@@ -581,6 +581,7 @@ def test_well_posed_composite_different_degradation(self):
         options = {
             "particle phases": ("2", "1"),
             "SEI": ("ec reaction limited", "none"),
+            "SEI porosity change": "true",
             "lithium plating": ("reversible", "none"),
             "open-circuit potential": (("current sigmoid", "single"), "single"),
         }
@@ -589,6 +590,7 @@ def test_well_posed_composite_different_degradation(self):
         options = {
             "particle phases": ("2", "1"),
             "SEI": (("ec reaction limited", "solvent-diffusion limited"), "none"),
+            "SEI porosity change": "true",
             "lithium plating": (("reversible", "irreversible"), "none"),
             "open-circuit potential": (("current sigmoid", "single"), "single"),
         }