Composite model for particle mechanics

.gitignore

@@ -138,3 +138,5 @@ results/
 
 # tests
 test_callback.log
+PublicPyBaMM_Composite/
+PyBaMM/

CHANGELOG.md

@@ -9,6 +9,7 @@
 - Added `BasicDFN` model for sodium-ion batteries ([#4451](https://github.com/pybamm-team/PyBaMM/pull/4451))
 - Added sensitivity calculation support for `pybamm.Simulation` and `pybamm.Experiment` ([#4415](https://github.com/pybamm-team/PyBaMM/pull/4415))
 - Added OpenMP parallelization to IDAKLU solver for lists of input parameters ([#4449](https://github.com/pybamm-team/PyBaMM/pull/4449))
+- Porosity change now works for composite electrode ([#4417](https://github.com/pybamm-team/PyBaMM/pull/4417))
 - Added phase-dependent particle options to LAM
   ([#4369](https://github.com/pybamm-team/PyBaMM/pull/4369))
 - Added a lithium ion equivalent circuit model with split open circuit voltages for each electrode (`SplitOCVR`). ([#4330](https://github.com/pybamm-team/PyBaMM/pull/4330))

src/pybamm/models/full_battery_models/base_battery_model.py

@@ -598,12 +598,13 @@ def __init__(self, extra_options):
         if isinstance(options["stress-induced diffusion"], str):
             if (
                 options["stress-induced diffusion"] == "true"
-                and options["particle mechanics"] == "none"
+                # and options["particle mechanics"] == "none"
             ):
-                raise pybamm.OptionError(
-                    "cannot have stress-induced diffusion without a particle "
-                    "mechanics model"
-                )
+                pass
+                # raise pybamm.OptionError(
+                #     "cannot have stress-induced diffusion without a particle "
+                #     "mechanics model"
+                # )
 
         if options["working electrode"] != "both":
             if options["thermal"] == "x-full":

src/pybamm/models/submodels/interface/base_interface.py

@@ -306,9 +306,14 @@ def _get_standard_volumetric_current_density_variables(self, variables):
         a_j_av = pybamm.x_average(a_j)
 
         if reaction_name == "SEI on cracks ":
-            roughness = variables[f"{Domain} electrode roughness ratio"] - 1
+            roughness = (
+                variables[f"{Domain} electrode {self.phase_name}roughness ratio"] - 1
+            )
             roughness_av = (
-                variables[f"X-averaged {domain} electrode roughness ratio"] - 1
+                variables[
+                    f"X-averaged {domain} electrode {self.phase_name}roughness ratio"
+                ]
+                - 1
             )
         else:
             roughness = 1

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

@@ -218,7 +218,9 @@ def _get_standard_concentration_variables(self, variables):
         elif self.reaction == "SEI on cracks":
             L_inner_cr = variables[f"{Domain} inner {reaction_name}thickness [m]"]
             L_outer_cr = variables[f"{Domain} outer {reaction_name}thickness [m]"]
-            roughness = variables[f"{Domain} electrode roughness ratio"]
+            roughness = variables[
+                f"{Domain} electrode {self.phase_name}roughness ratio"
+            ]
 
             n_inner_cr = L_inner_cr * L_to_n_inner * (roughness - 1)
             n_outer_cr = L_outer_cr * L_to_n_outer * (roughness - 1)

src/pybamm/models/submodels/particle/base_particle.py

@@ -29,7 +29,6 @@ def __init__(self, param, domain, options, phase="primary"):
 
     def _get_effective_diffusivity(self, c, T, current):
         domain, Domain = self.domain_Domain
-        domain_param = self.domain_param
         phase_param = self.phase_param
         domain_options = getattr(self.options, domain)
 
@@ -60,7 +59,7 @@ def _get_effective_diffusivity(self, c, T, current):
             E = pybamm.r_average(phase_param.E(sto, T))
             nu = phase_param.nu
             theta_M = Omega / (self.param.R * T) * (2 * Omega * E) / (9 * (1 - nu))
-            stress_factor = 1 + theta_M * (c - domain_param.c_0)
+            stress_factor = 1 + theta_M * (c - phase_param.c_0)
         else:
             stress_factor = 1
 

src/pybamm/models/submodels/particle_mechanics/base_mechanics.py

@@ -23,15 +23,15 @@ class BaseMechanics(pybamm.BaseSubModel):
 
     """
 
-    def __init__(self, param, domain, options, phase="primary"):
+    def __init__(self, param, domain, options, phase):
         super().__init__(param, domain, options=options, phase=phase)
 
     def _get_standard_variables(self, l_cr):
         domain, Domain = self.domain_Domain
         l_cr_av = pybamm.x_average(l_cr)
         variables = {
-            f"{Domain} particle crack length [m]": l_cr,
-            f"X-averaged {domain} particle crack length [m]": l_cr_av,
+            f"{Domain} {self.phase_name}particle crack length [m]": l_cr,
+            f"X-averaged {domain} {self.phase_name}particle crack length [m]": l_cr_av,
         }
         return variables
 
@@ -61,7 +61,7 @@ def _get_mechanical_results(self, variables):
         sto = variables[f"{Domain} {phase_name}particle concentration"]
         Omega = pybamm.r_average(phase_param.Omega(sto, T))
         R0 = phase_param.R
-        c_0 = domain_param.c_0
+        c_0 = phase_param.c_0
         E0 = pybamm.r_average(phase_param.E(sto, T))
         nu = phase_param.nu
         L0 = domain_param.L
@@ -129,19 +129,19 @@ def _get_standard_surface_variables(self, variables):
         domain, Domain = self.domain_Domain
         phase_name = self.phase_name
 
-        l_cr = variables[f"{Domain} particle crack length [m]"]
+        l_cr = variables[f"{Domain} {self.phase_name}particle crack length [m]"]
         a = variables[
             f"{Domain} electrode {phase_name}surface area to volume ratio [m-1]"
         ]
-        rho_cr = self.domain_param.rho_cr
-        w_cr = self.domain_param.w_cr
+        rho_cr = self.phase_param.rho_cr
+        w_cr = self.phase_param.w_cr
         roughness = 1 + 2 * l_cr * rho_cr * w_cr  # ratio of cracks to normal surface
         a_cr = (roughness - 1) * a  # crack surface area to volume ratio
 
         roughness_xavg = pybamm.x_average(roughness)
         variables = {
-            f"{Domain} crack surface to volume ratio [m-1]": a_cr,
-            f"{Domain} electrode roughness ratio": roughness,
-            f"X-averaged {domain} electrode roughness ratio": roughness_xavg,
+            f"{Domain} {self.phase_name}crack surface to volume ratio [m-1]": a_cr,
+            f"{Domain} electrode {self.phase_name}roughness ratio": roughness,
+            f"X-averaged {domain} electrode {self.phase_name}roughness ratio": roughness_xavg,
         }
         return variables

src/pybamm/models/submodels/particle_mechanics/crack_propagation.py

@@ -27,7 +27,7 @@ class CrackPropagation(BaseMechanics):
 
     """
 
-    def __init__(self, param, domain, x_average, options, phase="primary"):
+    def __init__(self, param, domain, x_average, options, phase):
         super().__init__(param, domain, options, phase)
         self.x_average = x_average
 
@@ -39,17 +39,17 @@ def get_fundamental_variables(self):
 
         if self.x_average is True:
             l_cr_av = pybamm.Variable(
-                f"X-averaged {domain} particle crack length [m]",
+                f"X-averaged {domain} {self.phase_name}particle crack length [m]",
                 domain="current collector",
-                scale=self.domain_param.l_cr_0,
+                scale=self.phase_param.l_cr_0,
             )
             l_cr = pybamm.PrimaryBroadcast(l_cr_av, f"{domain} electrode")
         else:
             l_cr = pybamm.Variable(
-                f"{Domain} particle crack length [m]",
+                f"{Domain} {self.phase_name}particle crack length [m]",
                 domain=f"{domain} electrode",
                 auxiliary_domains={"secondary": "current collector"},
-                scale=self.domain_param.l_cr_0,
+                scale=self.phase_param.l_cr_0,
             )
 
         variables = self._get_standard_variables(l_cr)
@@ -62,18 +62,20 @@ def get_coupled_variables(self, variables):
         variables.update(self._get_standard_surface_variables(variables))
         variables.update(self._get_mechanical_results(variables))
         T = variables[f"{Domain} electrode temperature [K]"]
-        k_cr = self.domain_param.k_cr(T)
-        m_cr = self.domain_param.m_cr
-        b_cr = self.domain_param.b_cr
-        stress_t_surf = variables[f"{Domain} particle surface tangential stress [Pa]"]
-        l_cr = variables[f"{Domain} particle crack length [m]"]
+        k_cr = self.phase_param.k_cr(T)
+        m_cr = self.phase_param.m_cr
+        b_cr = self.phase_param.b_cr
+        stress_t_surf = variables[
+            f"{Domain} {self.phase_name}particle surface tangential stress [Pa]"
+        ]
+        l_cr = variables[f"{Domain} {self.phase_name}particle crack length [m]"]
         # # compressive stress will not lead to crack propagation
         dK_SIF = stress_t_surf * b_cr * pybamm.sqrt(np.pi * l_cr) * (stress_t_surf >= 0)
         dl_cr = k_cr * (dK_SIF**m_cr) / 3600  # divide by 3600 to replace t0_cr
         variables.update(
             {
-                f"{Domain} particle cracking rate [m.s-1]": dl_cr,
-                f"X-averaged {domain} particle cracking rate [m.s-1]": pybamm.x_average(
+                f"{Domain} {self.phase_name}particle cracking rate [m.s-1]": dl_cr,
+                f"X-averaged {domain} {self.phase_name}particle cracking rate [m.s-1]": pybamm.x_average(
                     dl_cr
                 ),
             }
@@ -84,21 +86,29 @@ def set_rhs(self, variables):
         domain, Domain = self.domain_Domain
 
         if self.x_average is True:
-            l_cr = variables[f"X-averaged {domain} particle crack length [m]"]
-            dl_cr = variables[f"X-averaged {domain} particle cracking rate [m.s-1]"]
+            l_cr = variables[
+                f"X-averaged {domain} {self.phase_name}particle crack length [m]"
+            ]
+            dl_cr = variables[
+                f"X-averaged {domain} {self.phase_name}particle cracking rate [m.s-1]"
+            ]
         else:
-            l_cr = variables[f"{Domain} particle crack length [m]"]
-            dl_cr = variables[f"{Domain} particle cracking rate [m.s-1]"]
+            l_cr = variables[f"{Domain} {self.phase_name}particle crack length [m]"]
+            dl_cr = variables[
+                f"{Domain} {self.phase_name}particle cracking rate [m.s-1]"
+            ]
         self.rhs = {l_cr: dl_cr}
 
     def set_initial_conditions(self, variables):
         domain, Domain = self.domain_Domain
 
-        l_cr_0 = self.domain_param.l_cr_0
+        l_cr_0 = self.phase_param.l_cr_0
         if self.x_average is True:
             l_cr = variables[f"X-averaged {domain} particle crack length [m]"]
         else:
-            l_cr = variables[f"{Domain} particle crack length [m]"]
+            l_cr = variables[
+                f"X-averaged {domain} {self.phase_name}particle crack length [m]"
+            ]
             l_cr_0 = pybamm.PrimaryBroadcast(l_cr_0, f"{domain} electrode")
         self.initial_conditions = {l_cr: l_cr_0}
 
@@ -108,10 +118,12 @@ def add_events_from(self, variables):
         if self.x_average is True:
             l_cr = variables[f"X-averaged {domain} particle crack length [m]"]
         else:
-            l_cr = variables[f"{Domain} particle crack length [m]"]
+            l_cr = variables[
+                f"X-averaged {domain} {self.phase_name}particle crack length [m]"
+            ]
         self.events.append(
             pybamm.Event(
-                f"{domain} particle crack length larger than particle radius",
-                1 - pybamm.max(l_cr) / self.domain_param.prim.R_typ,
+                f"{domain} {self.phase_name}particle crack length larger than particle radius",
+                1 - pybamm.max(l_cr) / self.phase_param.R_typ,
             )
         )

src/pybamm/models/submodels/particle_mechanics/no_mechanics.py

@@ -22,7 +22,7 @@ class NoMechanics(BaseMechanics):
         Phase of the particle (default is "primary")
     """
 
-    def __init__(self, param, domain, options, phase="primary"):
+    def __init__(self, param, domain, options, phase):
         super().__init__(param, domain, options, phase)
 
     def get_fundamental_variables(self):
@@ -35,8 +35,8 @@ def get_fundamental_variables(self):
         variables = self._get_standard_variables(zero)
         variables.update(
             {
-                f"{Domain} particle cracking rate [m.s-1]": zero,
-                f"X-averaged {domain} particle cracking rate [m.s-1]": zero_av,
+                f"{Domain} {self.phase_name}particle cracking rate [m.s-1]": zero,
+                f"X-averaged {domain} {self.phase_name}particle cracking rate [m.s-1]": zero_av,
             }
         )
         return variables

src/pybamm/models/submodels/particle_mechanics/swelling_only.py

@@ -23,7 +23,7 @@ class SwellingOnly(BaseMechanics):
 
     """
 
-    def __init__(self, param, domain, options, phase="primary"):
+    def __init__(self, param, domain, options, phase):
         super().__init__(param, domain, options, phase)
 
         pybamm.citations.register("Ai2019")
@@ -39,8 +39,8 @@ def get_fundamental_variables(self):
         variables = self._get_standard_variables(zero)
         variables.update(
             {
-                f"{Domain} particle cracking rate [m.s-1]": zero,
-                f"X-averaged {domain} particle cracking rate [m.s-1]": zero_av,
+                f"{Domain} {self.phase_name}particle cracking rate [m.s-1]": zero,
+                f"X-averaged {domain} {self.phase_name}particle cracking rate [m.s-1]": zero_av,
             }
         )
         return variables

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

@@ -25,8 +25,9 @@ def __init__(self, param, options, x_average):
     def get_coupled_variables(self, variables):
         eps_dict = {}
         for domain in self.options.whole_cell_domains:
+            delta_eps_k = 0
             if domain == "separator":
-                delta_eps_k = 0  # separator porosity does not change
+                pass  # separator porosity does not change
             else:
                 Domain = domain.split()[0].capitalize()
                 L_sei_k = variables[f"{Domain} total SEI thickness [m]"]
@@ -37,7 +38,9 @@ def get_coupled_variables(self, variables):
                 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"]
+                roughness_k = variables[
+                    f"{Domain} electrode {self.phase_name}roughness ratio"
+                ]
 
                 L_tot = (
                     (L_sei_k - L_sei_0)
@@ -47,14 +50,15 @@ def get_coupled_variables(self, variables):
                 )
 
                 a_k = variables[
-                    f"{Domain} electrode surface area to volume ratio [m-1]"
+                    f"{Domain} electrode {self.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
+                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