Merge branch 'pybamm-team:develop' into issue#2028

CHANGELOG.md

@@ -1,5 +1,9 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
+## Features
+
+- Added support for particle size distributions combined with particle mechanics. ([#4807](https://github.com/pybamm-team/PyBaMM/pull/4807))
+
 # [v25.1.1](https://github.com/pybamm-team/PyBaMM/tree/v25.1.1) - 2025-01-20
 
 ## Features

src/pybamm/models/full_battery_models/base_battery_model.py

@@ -545,21 +545,11 @@ def __init__(self, extra_options):
                     "'quadratic' and 'quartic' concentration profiles have not yet "
                     "been implemented for particle-size ditributions"
                 )
-            if options["particle mechanics"] != "none":
-                raise NotImplementedError(
-                    "Particle mechanics submodels do not yet support particle-size"
-                    " distributions."
-                )
             if options["particle shape"] != "spherical":
                 raise NotImplementedError(
                     "Particle shape must be 'spherical' for particle-size distribution"
                     " submodels."
                 )
-            if options["stress-induced diffusion"] == "true":
-                raise NotImplementedError(
-                    "stress-induced diffusion cannot yet be included in "
-                    "particle-size distributions."
-                )
             if options["thermal"] == "x-full":
                 raise NotImplementedError(
                     "X-full thermal submodels do not yet support particle-size"

src/pybamm/models/submodels/active_material/loss_active_material.py

@@ -1,5 +1,5 @@
 #
-# Class for varying active material volume fraction, driven by stress
+# Class for varying active material volume fraction
 #
 import pybamm
 

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

@@ -24,6 +24,10 @@ class BaseMechanics(pybamm.BaseSubModel):
     """
 
     def __init__(self, param, domain, options, phase="primary"):
+        if options["particle size"] == "distribution":
+            self.size_distribution = True
+        else:
+            self.size_distribution = False
         super().__init__(param, domain, options=options, phase=phase)
 
     def _get_standard_variables(self, l_cr):
@@ -35,6 +39,71 @@ def _get_standard_variables(self, l_cr):
         }
         return variables
 
+    def _get_standard_size_distribution_variables(self, l_cr_dist):
+        domain, Domain = self.domain_Domain
+        l_cr_av_dist = pybamm.x_average(l_cr_dist)
+        variables = {
+            f"{Domain} {self.phase_param.phase_name}particle crack length distribution [m]": l_cr_dist,
+            f"X-averaged {domain} {self.phase_param.phase_name}particle crack length distribution [m]": l_cr_av_dist,
+        }
+        return variables
+
+    def _get_mechanical_size_distribution_results(self, variables):
+        domain, Domain = self.domain_Domain
+        phase_name = self.phase_param.phase_name
+        phase_param = self.phase_param
+        c_s_rav = variables[
+            f"R-averaged {domain} {phase_name}particle concentration distribution [mol.m-3]"
+        ]
+        c_s_surf = variables[
+            f"{Domain} {phase_name}particle surface concentration distribution [mol.m-3]"
+        ]
+        T = pybamm.PrimaryBroadcast(
+            variables[f"{Domain} electrode temperature [K]"],
+            [f"{domain} {phase_name}particle size"],
+        )
+        T = pybamm.PrimaryBroadcast(
+            T,
+            [f"{domain} {phase_name}particle"],
+        )
+
+        # use a tangential approximation for omega
+        c_0 = phase_param.c_0
+        R0 = phase_param.R
+        sto = variables[f"{Domain} {phase_name}particle concentration distribution"]
+        Omega = pybamm.r_average(phase_param.Omega(sto, T))
+
+        E0 = pybamm.r_average(phase_param.E(sto, T))
+        nu = phase_param.nu
+        # Ai2019 eq [10]
+        disp_surf = Omega * R0 / 3 * (c_s_rav - c_0)
+        # c0 reference concentration for no deformation
+        # stress evaluated at the surface of the particles
+        # Ai2019 eq [7] with r=R
+        stress_r_surf = pybamm.Scalar(0)
+        # Ai2019 eq [8] with r=R
+        # c_s_rav is already multiplied by 3/R^3 inside r_average
+        stress_t_surf = Omega * E0 * (c_s_rav - c_s_surf) / 3.0 / (1.0 - nu)
+
+        # Averages
+        stress_r_surf_av = pybamm.x_average(stress_r_surf)
+        stress_t_surf_av = pybamm.x_average(stress_t_surf)
+        disp_surf_av = pybamm.x_average(disp_surf)
+
+        variables.update(
+            {
+                f"{Domain} {phase_name}particle surface radial stress distribution [Pa]": stress_r_surf,
+                f"{Domain} {phase_name}particle surface tangential stress distribution [Pa]": stress_t_surf,
+                f"{Domain} {phase_name}particle surface displacement distribution [m]": disp_surf,
+                f"X-averaged {domain} {phase_name}particle surface "
+                "radial stress distribution [Pa]": stress_r_surf_av,
+                f"X-averaged {domain} {phase_name}particle surface "
+                "tangential stress distribution [Pa]": stress_t_surf_av,
+                f"X-averaged {domain} {phase_name}particle surface displacement distribution [m]": disp_surf_av,
+            }
+        )
+        return variables
+
     def _get_mechanical_results(self, variables):
         domain_param = self.domain_param
         domain, Domain = self.domain_Domain
@@ -58,10 +127,11 @@ def _get_mechanical_results(self, variables):
         ]
 
         # use a tangential approximation for omega
+        c_0 = phase_param.c_0
+        R0 = phase_param.R
         sto = variables[f"{Domain} {phase_name}particle concentration"]
         Omega = pybamm.r_average(phase_param.Omega(sto, T))
-        R0 = phase_param.R
-        c_0 = phase_param.c_0
+
         E0 = pybamm.r_average(phase_param.E(sto, T))
         nu = phase_param.nu
         L0 = domain_param.L
@@ -79,7 +149,7 @@ def _get_mechanical_results(self, variables):
         stress_r_surf = pybamm.Scalar(0)
         # Ai2019 eq [8] with r=R
         # c_s_rav is already multiplied by 3/R^3 inside r_average
-        stress_t_surf = Omega * E0 / 3.0 / (1.0 - nu) * (c_s_rav - c_s_surf)
+        stress_t_surf = Omega * E0 * (c_s_rav - c_s_surf) / 3.0 / (1.0 - nu)
 
         # Averages
         stress_r_surf_av = pybamm.x_average(stress_r_surf)

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

@@ -37,22 +37,50 @@ 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} {phase_name}particle crack length [m]",
-                domain="current collector",
-                scale=self.phase_param.l_cr_0,
-            )
+        if self.x_average:
+            if self.size_distribution:
+                l_cr_av_dist = pybamm.Variable(
+                    f"X-averaged {domain} {phase_name}particle crack length distribution [m]",
+                    domains={
+                        "primary": f"{domain} particle size",
+                        "secondary": "current collector",
+                    },
+                    scale=self.phase_param.l_cr_0,
+                )
+                l_cr_dist = pybamm.SecondaryBroadcast(
+                    l_cr_av_dist, f"{domain} electrode"
+                )
+                l_cr_av = pybamm.size_average(l_cr_av_dist)
+            else:
+                l_cr_av = pybamm.Variable(
+                    f"X-averaged {domain} {phase_name}particle crack length [m]",
+                    domain="current collector",
+                    scale=self.phase_param.l_cr_0,
+                )
             l_cr = pybamm.PrimaryBroadcast(l_cr_av, f"{domain} electrode")
         else:
-            l_cr = pybamm.Variable(
-                f"{Domain} {phase_name}particle crack length [m]",
-                domain=f"{domain} electrode",
-                auxiliary_domains={"secondary": "current collector"},
-                scale=self.phase_param.l_cr_0,
-            )
+            if self.size_distribution:
+                l_cr_dist = pybamm.Variable(
+                    f"{Domain} {phase_name}particle crack length distribution [m]",
+                    domains={
+                        "primary": f"{domain} particle size",
+                        "secondary": f"{domain} electrode",
+                        "tertiary": "current collector",
+                    },
+                    scale=self.phase_param.l_cr_0,
+                )
+                l_cr = pybamm.size_average(l_cr_dist)
+            else:
+                l_cr = pybamm.Variable(
+                    f"{Domain} {phase_name}particle crack length [m]",
+                    domain=f"{domain} electrode",
+                    auxiliary_domains={"secondary": "current collector"},
+                    scale=self.phase_param.l_cr_0,
+                )
 
         variables = self._get_standard_variables(l_cr)
+        if self.size_distribution:
+            variables.update(self._get_standard_size_distribution_variables(l_cr_dist))
 
         return variables
 
@@ -61,14 +89,26 @@ def get_coupled_variables(self, variables):
         phase_name = self.phase_param.phase_name
         variables.update(self._get_standard_surface_variables(variables))
         variables.update(self._get_mechanical_results(variables))
+        if self.size_distribution:
+            variables.update(self._get_mechanical_size_distribution_results(variables))
         T = variables[f"{Domain} electrode temperature [K]"]
         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]"]
+        if self.size_distribution:
+            stress_t_surf = variables[
+                f"{Domain} {phase_name}particle surface tangential stress distribution [Pa]"
+            ]
+        else:
+            stress_t_surf = variables[
+                f"{Domain} {phase_name}particle surface tangential stress [Pa]"
+            ]
+        if self.size_distribution:
+            l_cr = variables[
+                f"{Domain} {phase_name}particle crack length distribution [m]"
+            ]
+        else:
+            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
@@ -86,14 +126,24 @@ 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} {phase_name}particle crack length [m]"
-            ]
+            if self.size_distribution:
+                l_cr = variables[
+                    f"X-averaged {domain} {phase_name}particle crack length distribution [m]"
+                ]
+            else:
+                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} {phase_name}particle crack length [m]"]
+            if self.size_distribution:
+                l_cr = variables[
+                    f"{Domain} {phase_name}particle crack length distribution [m]"
+                ]
+            else:
+                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}
 
@@ -102,12 +152,25 @@ def set_initial_conditions(self, variables):
         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} {phase_name}particle crack length [m]"
-            ]
+            if self.size_distribution:
+                l_cr = variables[
+                    f"X-averaged {domain} {phase_name}particle crack length distribution [m]"
+                ]
+                l_cr_0 = pybamm.PrimaryBroadcast(l_cr_0, f"{domain} particle size")
+            else:
+                l_cr = variables[
+                    f"X-averaged {domain} {phase_name}particle crack length [m]"
+                ]
         else:
-            l_cr = variables[f"{Domain} {phase_name}particle crack length [m]"]
-            l_cr_0 = pybamm.PrimaryBroadcast(l_cr_0, f"{domain} electrode")
+            if self.size_distribution:
+                l_cr = variables[
+                    f"{Domain} {phase_name}particle crack length distribution [m]"
+                ]
+                l_cr_0 = pybamm.PrimaryBroadcast(l_cr_0, f"{domain} electrode")
+                l_cr_0 = pybamm.PrimaryBroadcast(l_cr_0, f"{domain} particle size")
+            else:
+                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):

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

@@ -48,4 +48,6 @@ def get_fundamental_variables(self):
     def get_coupled_variables(self, variables):
         variables.update(self._get_standard_surface_variables(variables))
         variables.update(self._get_mechanical_results(variables))
+        if self.size_distribution:
+            variables.update(self._get_mechanical_size_distribution_results(variables))
         return variables

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

@@ -525,6 +525,22 @@ def test_well_posed_psd(self):
         options = {"particle size": "distribution", "surface form": "algebraic"}
         self.check_well_posedness(options)
 
+    def test_well_posed_psd_swelling_and_cracking(self):
+        options = {
+            "particle size": "distribution",
+            "particle mechanics": "swelling and cracking",
+            "surface form": "algebraic",
+        }
+        self.check_well_posedness(options)
+
+    def test_well_posed_psd_swelling_only(self):
+        options = {
+            "particle size": "distribution",
+            "particle mechanics": "swelling only",
+            "surface form": "algebraic",
+        }
+        self.check_well_posedness(options)
+
     def test_well_posed_transport_efficiency_Bruggeman(self):
         options = {"transport efficiency": "Bruggeman"}
         self.check_well_posedness(options)

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

@@ -78,21 +78,6 @@ def test_nonspherical_particle_not_implemented(self):
         with pytest.raises(NotImplementedError):
             pybamm.lithium_ion.MPM(options)
 
-    def test_loss_active_material_stress_negative_not_implemented(self):
-        options = {"loss of active material": ("stress-driven", "none")}
-        with pytest.raises(NotImplementedError):
-            pybamm.lithium_ion.MPM(options)
-
-    def test_loss_active_material_stress_positive_not_implemented(self):
-        options = {"loss of active material": ("none", "stress-driven")}
-        with pytest.raises(NotImplementedError):
-            pybamm.lithium_ion.MPM(options)
-
-    def test_loss_active_material_stress_both_not_implemented(self):
-        options = {"loss of active material": "stress-driven"}
-        with pytest.raises(NotImplementedError):
-            pybamm.lithium_ion.MPM(options)
-
     def test_reversible_plating_with_porosity_not_implemented(self):
         options = {
             "lithium plating": "reversible",
@@ -101,11 +86,6 @@ def test_reversible_plating_with_porosity_not_implemented(self):
         with pytest.raises(pybamm.OptionError, match="distributions"):
             pybamm.lithium_ion.MPM(options)
 
-    def test_stress_induced_diffusion_not_implemented(self):
-        options = {"stress-induced diffusion": "true"}
-        with pytest.raises(NotImplementedError):
-            pybamm.lithium_ion.MPM(options)
-
     def test_msmr(self):
         options = {
             "open-circuit potential": "MSMR",
@@ -186,25 +166,3 @@ def test_ec_reaction_limited_not_implemented(self):
         }
         with pytest.raises(NotImplementedError):
             pybamm.lithium_ion.MPM(options)
-
-
-class TestMPMWithMechanics:
-    def test_well_posed_negative_cracking_not_implemented(self):
-        options = {"particle mechanics": ("swelling and cracking", "none")}
-        with pytest.raises(NotImplementedError):
-            pybamm.lithium_ion.MPM(options)
-
-    def test_well_posed_positive_cracking_not_implemented(self):
-        options = {"particle mechanics": ("none", "swelling and cracking")}
-        with pytest.raises(NotImplementedError):
-            pybamm.lithium_ion.MPM(options)
-
-    def test_well_posed_both_cracking_not_implemented(self):
-        options = {"particle mechanics": "swelling and cracking"}
-        with pytest.raises(NotImplementedError):
-            pybamm.lithium_ion.MPM(options)
-
-    def test_well_posed_both_swelling_only_not_implemented(self):
-        options = {"particle mechanics": "swelling only"}
-        with pytest.raises(NotImplementedError):
-            pybamm.lithium_ion.MPM(options)