mechanics for composite electrode

.gitignore

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

src/pybamm/models/full_battery_models/base_battery_model.py

@@ -594,12 +594,12 @@ 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"
-                )
+                # 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,9 @@ 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,7 @@ 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

@@ -60,7 +60,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,18 @@ 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 +84,21 @@ 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 +108,10 @@ 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

@@ -37,7 +37,7 @@ 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)

src/pybamm/parameters/lithium_ion_parameters.py

@@ -78,6 +78,7 @@ def _set_parameters(self):
         self.voltage_high_cut = self.elec.voltage_high_cut
         self.ocp_soc_0 = self.elec.ocp_soc_0
         self.ocp_soc_100 = self.elec.ocp_soc_100
+
 
         # Domain parameters
         for domain in self.domain_params.values():
@@ -268,20 +269,7 @@ 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"
@@ -306,14 +294,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):
         """
@@ -507,8 +487,19 @@ def _set_parameters(self):
         if self.options["particle shape"] == "spherical":
             self.a_typ = 3 * pybamm.xyz_average(self.epsilon_s) / self.R_typ
 
-        # Mechanical property
+        # Mechanical parameters
         self.nu = pybamm.Parameter(f"{pref}{Domain} electrode Poisson's ratio")
+        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")
 
         # Loss of active material parameters
         self.m_LAM = pybamm.Parameter(
@@ -774,6 +765,29 @@ def dxdU(self, U, T):
             dxdU += self.dxdU_j(U, T, i)
         return dxdU
 
+    def Omega(self, sto, T):
+        """Dimensional partial molar volume of Li in solid solution [m3.mol-1]"""
+        Domain = self.domain.capitalize()
+        inputs = {f"{self.phase_prefactor}{Domain} particle stoichiometry": sto, "Temperature [K]": T}
+        return pybamm.FunctionParameter(
+            f"{self.phase_prefactor}{Domain} electrode partial molar volume [m3.mol-1]", inputs
+        )
+    def E(self, sto, T):
+        """Dimensional Young's modulus"""
+        Domain = self.domain.capitalize()
+        inputs = {f"{self.phase_prefactor}{Domain} particle stoichiometry": sto, "Temperature [K]": T}
+        return pybamm.FunctionParameter(
+            f"{self.phase_prefactor}{Domain} electrode Young's modulus [Pa]", inputs
+        )
+    def k_cr(self, T):
+        """
+        Cracking rate for the electrode;
+        """
+        Domain = self.domain.capitalize()
+        return pybamm.FunctionParameter(
+            f"{self.phase_prefactor}{Domain} electrode cracking rate", {"Temperature [K]": T}
+        )
+
     def t_change(self, sto):
         """
         Volume change for the electrode; sto should be R-averaged
@@ -782,29 +796,8 @@ def t_change(self, sto):
         return pybamm.FunctionParameter(
             f"{self.phase_prefactor}{Domain} electrode volume change",
             {
-                f"{Domain} particle stoichiometry": sto,
+                f"{self.phase_prefactor}{Domain} electrode volume change",
             },
         )
 
-    def Omega(self, sto, T):
-        """Dimensional partial molar volume of Li in solid solution [m3.mol-1]"""
-        domain, Domain = self.domain_Domain
-        inputs = {
-            f"{self.phase_prefactor} particle stoichiometry": sto,
-            "Temperature [K]": T,
-        }
-        return pybamm.FunctionParameter(
-            f"{self.phase_prefactor}{Domain} electrode partial molar volume [m3.mol-1]",
-            inputs,
-        )
 
-    def E(self, sto, T):
-        """Dimensional Young's modulus"""
-        domain, Domain = self.domain_Domain
-        inputs = {
-            f"{self.phase_prefactor} particle stoichiometry": sto,
-            "Temperature [K]": T,
-        }
-        return pybamm.FunctionParameter(
-            f"{self.phase_prefactor}{Domain} electrode Young's modulus [Pa]", inputs
-        )