Merge pull request #2672 from pybamm-team/#2669-particle-tuple

Allow particle size to be a tuple and fix MPM half-cell

Author: rtimms

CHANGELOG.md

@@ -1,5 +1,9 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
+## Features
+
+- The "particle size" option can now be a tuple to allow different behaviour in each electrode([#2672](https://github.com/pybamm-team/PyBaMM/pull/2672)).
+
 # [v23.1](https://github.com/pybamm-team/PyBaMM/tree/v23.1) - 2023-01-31
 
 ## Features

examples/notebooks/models/using-submodels.ipynb

@@ -44,36 +44,56 @@ def battery_geometry(
         "separator": {"x_s": {"min": l_n, "max": l_n_l_s}},
         "positive electrode": {"x_p": {"min": l_n_l_s, "max": 1}},
     }
+
     # Add particle domains
     if include_particles is True:
         zero_one = {"min": 0, "max": 1}
-        geometry.update(
-            {
-                "negative particle": {"r_n": zero_one},
-                "positive particle": {"r_p": zero_one},
-            }
-        )
         for domain in ["negative", "positive"]:
-            phases = int(getattr(options, domain)["particle phases"])
-            if phases >= 2:
+            if options.electrode_types[domain] == "porous":
                 geometry.update(
                     {
-                        f"{domain} primary particle": {"r_n_prim": zero_one},
-                        f"{domain} secondary particle": {"r_n_sec": zero_one},
+                        f"{domain} particle": {f"r_{domain[0]}": zero_one},
                     }
                 )
+                phases = int(getattr(options, domain)["particle phases"])
+                if phases >= 2:
+                    geometry.update(
+                        {
+                            f"{domain} primary particle": {
+                                f"r_{domain[0]}_prim": zero_one
+                            },
+                            f"{domain} secondary particle": {
+                                f"r_{domain[0]}_sec": zero_one
+                            },
+                        }
+                    )
+
     # Add particle size domains
-    if options is not None and options["particle size"] == "distribution":
+    if (
+        options is not None
+        and options.negative["particle size"] == "distribution"
+        and options.electrode_types["negative"] == "porous"
+    ):
         R_min_n = geo.n.prim.R_min
-        R_min_p = geo.p.prim.R_min
         R_max_n = geo.n.prim.R_max
-        R_max_p = geo.p.prim.R_max
         geometry.update(
             {
                 "negative particle size": {"R_n": {"min": R_min_n, "max": R_max_n}},
+            }
+        )
+    if (
+        options is not None
+        and options.positive["particle size"] == "distribution"
+        and options.electrode_types["positive"] == "porous"
+    ):
+        R_min_p = geo.p.prim.R_min
+        R_max_p = geo.p.prim.R_max
+        geometry.update(
+            {
                 "positive particle size": {"R_p": {"min": R_min_p, "max": R_max_p}},
             }
         )
+
     # Add current collector domains
     current_collector_dimension = options["dimensionality"]
     if form_factor == "pouch":

pybamm/geometry/battery_geometry.py

@@ -73,11 +73,6 @@ class BatteryModelOptions(pybamm.FuzzyDict):
                 "stress-driven", "reaction-driven", or "stress and reaction-driven".
                 A 2-tuple can be provided for different behaviour in negative and
                 positive electrodes.
-            * "particle phases": str
-                Number of phases present in the electrode. A 2-tuple can be provided for
-                different behaviour in negative and positive electrodes.
-                For example, set to ("2", "1") for a negative electrode with 2 phases,
-                e.g. graphite and silicon.
             * "operating mode" : str
                 Sets the operating mode for the model. This determines how the current
                 is set. Can be:
@@ -97,7 +92,19 @@ class BatteryModelOptions(pybamm.FuzzyDict):
             * "particle" : str
                 Sets the submodel to use to describe behaviour within the particle.
                 Can be "Fickian diffusion" (default), "uniform profile",
-                "quadratic profile", or "quartic profile".
+                "quadratic profile", or "quartic profile". A 2-tuple can be provided 
+                for different behaviour in negative and positive electrodes.
+            * "particle mechanics" : str
+                Sets the model to account for mechanical effects such as particle
+                swelling and cracking. Can be "none" (default), "swelling only",
+                or "swelling and cracking".
+                A 2-tuple can be provided for different behaviour in negative and
+                positive electrodes.
+            * "particle phases": str
+                Number of phases present in the electrode. A 2-tuple can be provided for
+                different behaviour in negative and positive electrodes.
+                For example, set to ("2", "1") for a negative electrode with 2 phases,
+                e.g. graphite and silicon.
             * "particle shape" : str
                 Sets the model shape of the electrode particles. This is used to
                 calculate the surface area to volume ratio. Can be "spherical"
@@ -106,12 +113,6 @@ class BatteryModelOptions(pybamm.FuzzyDict):
                 Sets the model to include a single active particle size or a
                 distribution of sizes at any macroscale location. Can be "single"
                 (default) or "distribution". Option applies to both electrodes.
-            * "particle mechanics" : str
-                Sets the model to account for mechanical effects such as particle
-                swelling and cracking. Can be "none" (default), "swelling only",
-                or "swelling and cracking".
-                A 2-tuple can be provided for different behaviour in negative and
-                positive electrodes.
             * "SEI" : str
                 Set the SEI submodel to be used. Options are:
 
@@ -528,7 +529,7 @@ def __init__(self, extra_options):
             ):
                 raise pybamm.OptionError(
                     "If there are multiple particle phases: 'surface form' cannot be "
-                    "'false', 'particle size' must be 'false', 'particle' must be "
+                    "'false', 'particle size' must be 'single', 'particle' must be "
                     "'Fickian diffusion'. Also the following must "
                     "be 'none': 'particle mechanics', "
                     "'loss of active material', 'lithium plating'"
@@ -554,9 +555,10 @@ def __init__(self, extra_options):
                                 "interface utilisation",
                                 "loss of active material",
                                 "open circuit potential",
-                                "particle mechanics",
                                 "particle",
+                                "particle mechanics",
                                 "particle phases",
+                                "particle size",
                                 "stress-induced diffusion",
                             ]
                             and isinstance(value, tuple)
@@ -845,12 +847,14 @@ def options(self, extra_options):
                 raise pybamm.OptionError(
                     "x-average side reactions cannot be 'false' for SPM models"
                 )
-        if isinstance(self, pybamm.lithium_ion.SPM) and not isinstance(
-            self, pybamm.lithium_ion.MPM
-        ):
-            if options["particle size"] == "distribution":
+        if isinstance(self, pybamm.lithium_ion.SPM):
+            if (
+                "distribution" in options["particle size"]
+                and options["surface form"] == "false"
+            ):
                 raise pybamm.OptionError(
-                    "'particle size' should be 'single' for SPM and SPMe models"
+                    "surface form must be 'algebraic' or 'differential' if "
+                    " 'particle size' contains a 'distribution'"
                 )
         if isinstance(self, pybamm.lead_acid.BaseModel):
             if options["thermal"] != "isothermal" and options["dimensionality"] != 0:

pybamm/models/full_battery_models/base_battery_model.py

@@ -54,5 +54,7 @@ def __init__(self, options=None, name="Many-Particle Model", build=True):
     @property
     def default_parameter_values(self):
         default_params = super().default_parameter_values
-        default_params = pybamm.get_size_distribution_parameters(default_params)
+        default_params = pybamm.get_size_distribution_parameters(
+            default_params, electrode=self.options["working electrode"]
+        )
         return default_params

pybamm/models/full_battery_models/lithium_ion/mpm.py

@@ -83,10 +83,11 @@ def _get_standard_active_material_variables(self, eps_solid):
             # R_n, R_p. For a size distribution, calculate the area-weighted
             # mean using the distribution instead. Then the surface area is
             # calculated the same way
-            if self.options["particle size"] == "single":
+            domain_options = getattr(self.options, domain)
+            if domain_options["particle size"] == "single":
                 R = self.phase_param.R
                 R_dim = self.phase_param.R_dimensional
-            elif self.options["particle size"] == "distribution":
+            elif domain_options["particle size"] == "distribution":
                 if self.domain == "negative":
                     R_ = pybamm.standard_spatial_vars.R_n
                 elif self.domain == "positive":

pybamm/models/submodels/active_material/base_active_material.py

@@ -73,7 +73,8 @@ def _get_exchange_current_density(self, variables):
         if self.reaction == "lithium-ion main":
             # For "particle-size distribution" submodels, take distribution version
             # of c_s_surf that depends on particle size.
-            if self.options["particle size"] == "distribution":
+            domain_options = getattr(self.options, domain)
+            if domain_options["particle size"] == "distribution":
                 c_s_surf = variables[
                     f"{Domain} {phase_name}particle surface concentration distribution"
                 ]

pybamm/models/submodels/interface/base_interface.py

@@ -72,16 +72,20 @@ def get_coupled_variables(self, variables):
                 delta_phi = delta_phi.orphans[0]
         # For "particle-size distribution" models, delta_phi must then be
         # broadcast to "particle size" domain
+        domain_options = getattr(self.options, domain)
         if (
             self.reaction == "lithium-ion main"
-            and self.options["particle size"] == "distribution"
+            and domain_options["particle size"] == "distribution"
         ):
             delta_phi = pybamm.PrimaryBroadcast(delta_phi, [f"{domain} particle size"])
 
         # Get exchange-current density
         j0 = self._get_exchange_current_density(variables)
         # Get open-circuit potential variables and reaction overpotential
-        if self.options["particle size"] == "distribution":
+        if (
+            domain_options["particle size"] == "distribution"
+            and self.options.electrode_types[domain] == "porous"
+        ):
             ocp = variables[
                 f"{Domain} electrode {reaction_name}open circuit potential distribution"
             ]

pybamm/models/submodels/interface/kinetics/base_kinetics.py

@@ -12,13 +12,26 @@ def get_coupled_variables(self, variables):
         m_lith = pybamm.sigmoid(current, 0, k)  # for lithation (current < 0)
         m_delith = 1 - m_lith  # for delithiation (current > 0)
 
-        Domain = self.domain.capitalize()
+        domain, Domain = self.domain_Domain
         phase_name = self.phase_name
 
         if self.reaction == "lithium-ion main":
             T = variables[f"{Domain} electrode temperature"]
-            # Particle size distribution is not yet implemented
-            if self.options["particle size"] != "distribution":
+            # For "particle-size distribution" models, take distribution version
+            # of c_s_surf that depends on particle size.
+            domain_options = getattr(self.options, domain)
+            if domain_options["particle size"] == "distribution":
+                c_s_surf = variables[
+                    f"{Domain} {phase_name}particle surface concentration distribution"
+                ]
+                # If variable was broadcast, take only the orphan
+                if isinstance(c_s_surf, pybamm.Broadcast) and isinstance(
+                    T, pybamm.Broadcast
+                ):
+                    c_s_surf = c_s_surf.orphans[0]
+                    T = T.orphans[0]
+                T = pybamm.PrimaryBroadcast(T, [f"{domain} particle size"])
+            else:
                 c_s_surf = variables[
                     f"{Domain} {phase_name}particle surface concentration"
                 ]

pybamm/models/submodels/interface/open_circuit_potential/current_sigmoid_ocp.py

@@ -14,7 +14,8 @@ def get_coupled_variables(self, variables):
             T = variables[f"{Domain} electrode temperature"]
             # For "particle-size distribution" models, take distribution version
             # of c_s_surf that depends on particle size.
-            if self.options["particle size"] == "distribution":
+            domain_options = getattr(self.options, domain)
+            if domain_options["particle size"] == "distribution":
                 c_s_surf = variables[
                     f"{Domain} {phase_name}particle surface concentration distribution"
                 ]