Make particle mechanics work with composite electrodes (#4698)

* composite working

* composite swelling plus crack length works

* fix notebook

---------

Co-authored-by: Alexander Bills <[email protected]>
Co-authored-by: Robert Timms <[email protected]>

docs/source/examples/notebooks/models/loss_of_active_materials.ipynb

@@ -443,6 +443,8 @@
     "second = 0.1\n",
     "parameter_values.update(\n",
     "    {\n",
+    "        \"Primary: Negative electrode reference concentration for free of deformation [mol.m-3]\": 0.0,\n",
+    "        \"Secondary: Negative electrode reference concentration for free of deformation [mol.m-3]\": 0.0,\n",
     "        \"Primary: Negative electrode LAM constant proportional term [s-1]\": 1e-4 / 3600,\n",
     "        \"Secondary: Negative electrode LAM constant proportional term [s-1]\": 1e-4\n",
     "        / 3600\n",

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

@@ -326,9 +326,10 @@ 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} {phase_name}electrode roughness ratio"] - 1
             roughness_av = (
-                variables[f"X-averaged {domain} electrode roughness ratio"] - 1
+                variables[f"X-averaged {domain} {phase_name}electrode roughness ratio"]
+                - 1
             )
         else:
             roughness = 1

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

@@ -120,6 +120,7 @@ def _get_standard_concentration_variables(self, variables):
         """Update variables related to the SEI concentration."""
         domain, Domain = self.domain_Domain
         phase_param = self.phase_param
+        phase_name = phase_param.phase_name
         reaction_name = self.reaction_name
 
         # Set scales to one for the "no SEI" model so that they are not required
@@ -189,7 +190,7 @@ def _get_standard_concentration_variables(self, variables):
         # Concentration variables are handled slightly differently for SEI on cracks
         elif self.reaction == "SEI on cracks":
             L_cr = variables[f"{Domain} {reaction_name}thickness [m]"]
-            roughness = variables[f"{Domain} electrode roughness ratio"]
+            roughness = variables[f"{Domain} {phase_name}electrode roughness ratio"]
 
             n_SEI_cr = L_cr * L_to_n * (roughness - 1)  # SEI on cracks concentration
             if self.size_distribution:

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

@@ -30,8 +30,8 @@ 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_param.phase_name}particle crack length [m]": l_cr,
+            f"X-averaged {domain} {self.phase_param.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
@@ -127,21 +127,21 @@ def _get_mechanical_results(self, variables):
 
     def _get_standard_surface_variables(self, variables):
         domain, Domain = self.domain_Domain
-        phase_name = self.phase_name
+        phase_name = self.phase_param.phase_name
 
-        l_cr = variables[f"{Domain} particle crack length [m]"]
+        l_cr = variables[f"{Domain} {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} {phase_name}crack surface to volume ratio [m-1]": a_cr,
+            f"{Domain} {phase_name}electrode roughness ratio": roughness,
+            f"X-averaged {domain} {phase_name}electrode roughness ratio": roughness_xavg,
         }
         return variables

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

@@ -36,20 +36,20 @@ def __init__(self, param, domain, x_average, options, phase="primary"):
 
     def get_fundamental_variables(self):
         domain, Domain = self.domain_Domain
-
+        phase_name = self.phase_param.phase_name
         if self.x_average is True:
             l_cr_av = pybamm.Variable(
-                f"X-averaged {domain} particle crack length [m]",
+                f"X-averaged {domain} {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} {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)
@@ -58,22 +58,24 @@ def get_fundamental_variables(self):
 
     def get_coupled_variables(self, variables):
         domain, Domain = self.domain_Domain
-
+        phase_name = self.phase_param.phase_name
         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} {phase_name}particle surface tangential stress [Pa]"
+        ]
+        l_cr = variables[f"{Domain} {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} {phase_name}particle cracking rate [m.s-1]": dl_cr,
+                f"X-averaged {domain} {phase_name}particle cracking rate [m.s-1]": pybamm.x_average(
                     dl_cr
                 ),
             }
@@ -82,36 +84,44 @@ def get_coupled_variables(self, variables):
 
     def set_rhs(self, variables):
         domain, Domain = self.domain_Domain
-
+        phase_name = self.phase_param.phase_name
         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} {phase_name}particle crack length [m]"
+            ]
+            dl_cr = variables[
+                f"X-averaged {domain} {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} {phase_name}particle crack length [m]"]
+            dl_cr = variables[f"{Domain} {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
+        phase_name = self.phase_param.phase_name
+        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]"]
+            l_cr = variables[
+                f"X-averaged {domain} {phase_name}particle crack length [m]"
+            ]
         else:
-            l_cr = variables[f"{Domain} particle crack length [m]"]
+            l_cr = variables[f"{Domain} {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}
 
     def add_events_from(self, variables):
         domain, Domain = self.domain_Domain
-
+        phase_name = self.phase_param.phase_name
         if self.x_average is True:
-            l_cr = variables[f"X-averaged {domain} particle crack length [m]"]
+            l_cr = variables[
+                f"X-averaged {domain} {phase_name}particle crack length [m]"
+            ]
         else:
-            l_cr = variables[f"{Domain} particle crack length [m]"]
+            l_cr = variables[f"{Domain} {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} {phase_name} particle crack length larger than particle radius",
+                1 - pybamm.max(l_cr) / self.phase_param.R_typ,
             )
         )

src/pybamm/parameters/lead_acid_parameters.py

@@ -344,6 +344,7 @@ def __init__(self, phase, domain_param):
         self.domain = domain_param.domain
         self.main_param = domain_param.main_param
         self.geo = domain_param.geo.prim
+        self.phase_name = phase
 
     def _set_parameters(self):
         main = self.main_param

src/pybamm/parameters/lithium_ion_parameters.py

@@ -269,20 +269,6 @@ def _set_parameters(self):
         self.b_s = self.geo.b_s
         self.tau_s = self.geo.tau_s
 
-        # Mechanical parameters
-        self.c_0 = pybamm.Parameter(
-            f"{Domain} electrode reference concentration for free of deformation "
-            "[mol.m-3]"
-        )
-
-        self.l_cr_0 = pybamm.Parameter(f"{Domain} electrode initial crack length [m]")
-        self.w_cr = pybamm.Parameter(f"{Domain} electrode initial crack width [m]")
-        self.rho_cr = pybamm.Parameter(
-            f"{Domain} electrode number of cracks per unit area [m-2]"
-        )
-        self.b_cr = pybamm.Parameter(f"{Domain} electrode Paris' law constant b")
-        self.m_cr = pybamm.Parameter(f"{Domain} electrode Paris' law constant m")
-
         # Utilisation parameters
         self.u_init = pybamm.Parameter(
             f"Initial {domain} electrode interface utilisation"
@@ -307,15 +293,6 @@ def sigma(self, T):
             f"{Domain} electrode conductivity [S.m-1]", inputs
         )
 
-    def k_cr(self, T):
-        """
-        Cracking rate for the electrode;
-        """
-        Domain = self.domain.capitalize()
-        return pybamm.FunctionParameter(
-            f"{Domain} electrode cracking rate", {"Temperature [K]": T}
-        )
-
     def LAM_rate_current(self, i, T):
         """
         Dimensional rate of loss of active material as a function of applied current
@@ -516,6 +493,34 @@ def _set_parameters(self):
             f"{pref}{Domain} electrode reaction-driven LAM factor [m3.mol-1]"
         )
 
+        # Mechanical parameters
+        self.c_0 = pybamm.Parameter(
+            f"{pref}{Domain} electrode reference concentration for free of deformation "
+            "[mol.m-3]"
+        )
+
+        self.l_cr_0 = pybamm.Parameter(
+            f"{pref}{Domain} electrode initial crack length [m]"
+        )
+        self.w_cr = pybamm.Parameter(
+            f"{pref}{Domain} electrode initial crack width [m]"
+        )
+        self.rho_cr = pybamm.Parameter(
+            f"{pref}{Domain} electrode number of cracks per unit area [m-2]"
+        )
+        self.b_cr = pybamm.Parameter(f"{pref}{Domain} electrode Paris' law constant b")
+        self.m_cr = pybamm.Parameter(f"{pref}{Domain} electrode Paris' law constant m")
+
+    def k_cr(self, T):
+        """
+        Cracking rate for the electrode;
+        """
+        phase_prefactor = self.phase_prefactor
+        Domain = self.domain.capitalize()
+        return pybamm.FunctionParameter(
+            f"{phase_prefactor}{Domain} electrode cracking rate", {"Temperature [K]": T}
+        )
+
     def D(self, c_s, T, lithiation=None):
         """
         Dimensional diffusivity in particle. In the parameter sets this is defined as

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

@@ -462,7 +462,8 @@ def lico2_volume_change_Ai2020(sto):
                 "Secondary: Negative electrode partial molar volume [m3.mol-1]": 3.1e-06,
                 "Secondary: Negative electrode Young's modulus [Pa]": 15000000000.0,
                 "Secondary: Negative electrode Poisson's ratio": 0.3,
-                "Negative electrode reference concentration for free of deformation [mol.m-3]": 0.0,
+                "Primary: Negative electrode reference concentration for free of deformation [mol.m-3]": 0.0,
+                "Secondary: Negative electrode reference concentration for free of deformation [mol.m-3]": 0.0,
                 "Primary: Negative electrode volume change": graphite_volume_change_Ai2020,
                 "Secondary: Negative electrode volume change": graphite_volume_change_Ai2020,
                 "Positive electrode partial molar volume [m3.mol-1]": -7.28e-07,