Make particle size distribution work with composite electrode.

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+- Made composite electrode model compatible with particle size distribution ([#4687](https://github.com/pybamm-team/PyBaMM/pull/4687))
 - Added `Symbol.post_order()` method to return an iterable that steps through the tree in post-order fashion. ([#4684](https://github.com/pybamm-team/PyBaMM/pull/4684))
 - Added two more submodels (options) for the SEI: Lars von Kolzenberg (2020) model and Tunneling Limit model ([#4394](https://github.com/pybamm-team/PyBaMM/pull/4394))
 

docs/source/examples/notebooks/models/using-submodels.ipynb

@@ -527,10 +527,10 @@
     "    model.param, \"positive\", model.options, cracks=True\n",
     ")\n",
     "model.submodels[\"Negative lithium plating\"] = pybamm.lithium_plating.NoPlating(\n",
-    "    model.param, \"Negative\"\n",
+    "    model.param, \"negative\"\n",
     ")\n",
     "model.submodels[\"Positive lithium plating\"] = pybamm.lithium_plating.NoPlating(\n",
-    "    model.param, \"Positive\"\n",
+    "    model.param, \"positive\"\n",
     ")"
    ]
   },

src/pybamm/expression_tree/averages.py

@@ -306,6 +306,10 @@
     return yz_average(x_average(symbol))
 
 
+def xyzs_average(symbol: pybamm.Symbol) -> pybamm.Symbol:
+    return xyz_average(size_average(symbol))
+
+
 def r_average(symbol: pybamm.Symbol) -> pybamm.Symbol:
     """
     Convenience function for creating an average in the r-direction.
@@ -373,7 +377,15 @@
     # If symbol doesn't have a domain, or doesn't have "negative particle size"
     #  or "positive particle size" as a domain, it's average value is itself
     if symbol.domain == [] or not any(
-        domain in [["negative particle size"], ["positive particle size"]]
+        domain
+        in [
+            ["negative particle size"],
+            ["positive particle size"],
+            ["negative primary particle size"],
+            ["positive primary particle size"],
+            ["negative secondary particle size"],
+            ["positive secondary particle size"],
+        ]
         for domain in list(symbol.domains.values())
     ):
         return symbol
@@ -394,13 +406,44 @@
     else:
         if f_a_dist is None:
             geo = pybamm.geometric_parameters
-            R = pybamm.SpatialVariable(
-                "R", domains=symbol.domains, coord_sys="cartesian"
-            )
-            if ["negative particle size"] in symbol.domains.values():
+            if "negative" in symbol.domain[0]:
+                if "primary" in symbol.domain[0]:
+                    R = pybamm.SpatialVariable(
+                        "R_n_prim", domains=symbol.domains, coord_sys="cartesian"
+                    )
+                elif "secondary" in symbol.domain[0]:
+                    R = pybamm.SpatialVariable(
+                        "R_n_sec", domains=symbol.domains, coord_sys="cartesian"
+                    )
+                else:
+                    R = pybamm.SpatialVariable(
+                        "R_n", domains=symbol.domains, coord_sys="cartesian"
+                    )
+            else:
+                if "primary" in symbol.domain[0]:
+                    R = pybamm.SpatialVariable(
+                        "R_p_prim", domains=symbol.domains, coord_sys="cartesian"
+                    )
+                elif "secondary" in symbol.domain[0]:
+                    R = pybamm.SpatialVariable(
+                        "R_p_sec", domains=symbol.domains, coord_sys="cartesian"
+                    )
+                else:
+                    R = pybamm.SpatialVariable(
+                        "R_p", domains=symbol.domains, coord_sys="cartesian"
+                    )
+            if ["negative particle size"] in symbol.domains.values() or [
+                "negative primary particle size"
+            ] in symbol.domains.values():
                 f_a_dist = geo.n.prim.f_a_dist(R)
-            elif ["positive particle size"] in symbol.domains.values():
+            elif ["negative secondary particle size"] in symbol.domains.values():
+                f_a_dist = geo.n.sec.f_a_dist(R)
+            elif ["positive particle size"] in symbol.domains.values() or [
+                "positive primary particle size"
+            ] in symbol.domains.values():
                 f_a_dist = geo.p.prim.f_a_dist(R)
+            elif ["positive secondary particle size"] in symbol.domains.values():
+                f_a_dist = geo.p.sec.f_a_dist(R)
         return SizeAverage(symbol, f_a_dist)
 
 

src/pybamm/geometry/battery_geometry.py

@@ -83,6 +83,18 @@
                 "negative particle size": {"R_n": {"min": R_min_n, "max": R_max_n}},
             }
         )
+        phases = int(options.negative["particle phases"])
+        if phases >= 2:
+            geometry.update(
+                {
+                    "negative primary particle size": {
+                        "R_n_prim": {"min": R_min_n, "max": R_max_n}
+                    },
+                    "negative secondary particle size": {
+                        "R_n_sec": {"min": R_min_n, "max": R_max_n}
+                    },
+                }
+            )
     if (
         options is not None
         and options.positive["particle size"] == "distribution"
@@ -95,6 +107,18 @@
                 "positive particle size": {"R_p": {"min": R_min_p, "max": R_max_p}},
             }
         )
+        phases = int(options.positive["particle phases"])
+        if phases >= 2:
+            geometry.update(
+                {
+                    "positive primary particle size": {
+                        "R_p_prim": {"min": R_min_p, "max": R_max_p}
+                    },
+                    "positive secondary particle size": {
+                        "R_p_sec": {"min": R_min_p, "max": R_max_p}
+                    },
+                }
+            )
 
     # Add current collector domains
     current_collector_dimension = options["dimensionality"]

src/pybamm/geometry/standard_spatial_vars.py

@@ -98,6 +98,7 @@
     },
     coord_sys="cartesian",
 )
+
 R_p = pybamm.SpatialVariable(
     "R_p",
     domain=["positive particle size"],
@@ -108,6 +109,44 @@
     coord_sys="cartesian",
 )
 
+R_n_prim = pybamm.SpatialVariable(
+    "R_n_prim",
+    domain=["negative primary particle size"],
+    auxiliary_domains={
+        "secondary": "negative electrode",
+        "tertiary": "current collector",
+    },
+    coord_sys="cartesian",
+)
+R_p_prim = pybamm.SpatialVariable(
+    "R_p_prim",
+    domain=["positive primary particle size"],
+    auxiliary_domains={
+        "secondary": "positive electrode",
+        "tertiary": "current collector",
+    },
+    coord_sys="cartesian",
+)
+
+R_n_sec = pybamm.SpatialVariable(
+    "R_n_sec",
+    domain=["negative secondary particle size"],
+    auxiliary_domains={
+        "secondary": "negative electrode",
+        "tertiary": "current collector",
+    },
+    coord_sys="cartesian",
+)
+R_p_sec = pybamm.SpatialVariable(
+    "R_p_sec",
+    domain=["positive secondary particle size"],
+    auxiliary_domains={
+        "secondary": "positive electrode",
+        "tertiary": "current collector",
+    },
+    coord_sys="cartesian",
+)
+
 # Domains at cell edges
 x_n_edge = pybamm.SpatialVariableEdge(
     "x_n",

src/pybamm/models/full_battery_models/base_battery_model.py

@@ -623,7 +623,6 @@ def __init__(self, extra_options):
         if options["particle phases"] not in ["1", ("1", "1")]:
             if not (
                 options["surface form"] != "false"
-                and options["particle size"] == "single"
                 and options["particle"] == "Fickian diffusion"
             ):
                 raise pybamm.OptionError(
@@ -868,6 +867,10 @@ def default_var_pts(self):
             "z": 10,
             "R_n": 30,
             "R_p": 30,
+            "R_n_prim": 30,
+            "R_p_prim": 30,
+            "R_n_sec": 30,
+            "R_p_sec": 30,
         }
         # Reduce the default points for 2D current collectors
         if self.options["dimensionality"] == 2:
@@ -888,6 +891,10 @@ def default_submesh_types(self):
             "positive secondary particle": pybamm.Uniform1DSubMesh,
             "negative particle size": pybamm.Uniform1DSubMesh,
             "positive particle size": pybamm.Uniform1DSubMesh,
+            "negative primary particle size": pybamm.Uniform1DSubMesh,
+            "positive primary particle size": pybamm.Uniform1DSubMesh,
+            "negative secondary particle size": pybamm.Uniform1DSubMesh,
+            "positive secondary particle size": pybamm.Uniform1DSubMesh,
         }
         if self.options["dimensionality"] == 0:
             base_submeshes["current collector"] = pybamm.SubMesh0D
@@ -910,6 +917,10 @@ def default_spatial_methods(self):
             "positive secondary particle": pybamm.FiniteVolume(),
             "negative particle size": pybamm.FiniteVolume(),
             "positive particle size": pybamm.FiniteVolume(),
+            "negative primary particle size": pybamm.FiniteVolume(),
+            "positive primary particle size": pybamm.FiniteVolume(),
+            "negative secondary particle size": pybamm.FiniteVolume(),
+            "positive secondary particle size": pybamm.FiniteVolume(),
         }
         if self.options["dimensionality"] == 0:
             # 0D submesh - use base spatial method

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

@@ -81,9 +81,19 @@
                 R = self.phase_param.R
             elif domain_options["particle size"] == "distribution":
                 if self.domain == "negative":
-                    R_ = pybamm.standard_spatial_vars.R_n
+                    if self.phase_name == "primary ":
+                        R_ = pybamm.standard_spatial_vars.R_n_prim
+                    elif self.phase_name == "secondary ":
+                        R_ = pybamm.standard_spatial_vars.R_n_sec
+                    else:
+                        R_ = pybamm.standard_spatial_vars.R_n
                 elif self.domain == "positive":
-                    R_ = pybamm.standard_spatial_vars.R_p
+                    if self.phase_name == "primary ":
+                        R_ = pybamm.standard_spatial_vars.R_p_prim
+                    elif self.phase_name == "secondary ":
+                        R_ = pybamm.standard_spatial_vars.R_p_sec
+                    else:
+                        R_ = pybamm.standard_spatial_vars.R_p
                 R = pybamm.size_average(R_)
 
             R_av = pybamm.x_average(R)

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

@@ -46,6 +46,8 @@
             self.reaction_name = self.phase_name + self.reaction_name
 
         self.reaction = reaction
+        domain_options = getattr(self.options, domain)
+        self.size_distribution = domain_options["particle size"] == "distribution"
 
     def _get_exchange_current_density(self, variables):
         """
@@ -93,8 +95,10 @@
                     # "current collector"
                     c_e = pybamm.PrimaryBroadcast(c_e, ["current collector"])
                 # broadcast c_e, T onto "particle size"
-                c_e = pybamm.PrimaryBroadcast(c_e, [f"{domain} particle size"])
-                T = pybamm.PrimaryBroadcast(T, [f"{domain} particle size"])
+                c_e = pybamm.PrimaryBroadcast(
+                    c_e, [f"{domain} {phase_name}particle size"]
+                )
+                T = pybamm.PrimaryBroadcast(T, [f"{domain} {phase_name}particle size"])
 
             else:
                 c_s_surf = variables[
@@ -215,13 +219,23 @@
 
         # Size average. For j variables that depend on particle size, see
         # "_get_standard_size_distribution_interfacial_current_variables"
-        if j.domain in [["negative particle size"], ["positive particle size"]]:
+        if j.domain in [
+            ["negative particle size"],
+            ["positive particle size"],
+            ["negative primary particle size"],
+            ["positive primary particle size"],
+            ["negative secondary particle size"],
+            ["positive secondary particle size"],
+        ]:
             j = pybamm.size_average(j)
         # Average, and broadcast if necessary
         j_av = pybamm.x_average(j)
         if j.domain == []:
             j = pybamm.FullBroadcast(j, f"{domain} electrode", "current collector")
-        elif j.domain == ["current collector"]:
+        elif (
+            j.domain == ["current collector"]
+            and getattr(self, "reaction_loc", None) != "interface"
+        ):
             j = pybamm.PrimaryBroadcast(j, f"{domain} electrode")
 
         variables = {
@@ -230,6 +244,13 @@
             f"X-averaged {domain} electrode {reaction_name}"
             "interfacial current density [A.m-2]": j_av,
         }
+        if self.phase_name == reaction_name:
+            variables.update(
+                {
+                    f"{Domain} electrode {reaction_name}interfacial current density [A.m-2]": j,
+                    f"X-averaged {domain} electrode {reaction_name}interfacial current density [A.m-2]": j_av,
+                }
+            )
 
         return variables
 
@@ -298,10 +319,18 @@
             a = variables[
                 f"{Domain} electrode {phase_name}surface area to volume ratio [m-1]"
             ]
-
-        j = variables[
-            f"{Domain} electrode {reaction_name}interfacial current density [A.m-2]"
-        ]
+        if "plating" in reaction_name:
+            j = variables[
+                f"{Domain} electrode {reaction_name}interfacial current density [A.m-2]"
+            ]
+        elif phase_name in reaction_name:
+            j = variables[
+                f"{Domain} electrode {reaction_name}interfacial current density [A.m-2]"
+            ]
+        else:
+            j = variables[
+                f"{Domain} electrode {phase_name}{reaction_name}interfacial current density [A.m-2]"
+            ]
         a_j = a * j
         a_j_av = pybamm.x_average(a_j)
 
@@ -410,9 +439,11 @@
         if j.domains["secondary"] == [f"{domain} electrode"]:
             # x-average
             j_xav = pybamm.x_average(j)
-        else:
+        elif getattr(self, "reaction_loc", None) != "interface":
             j_xav = j
             j = pybamm.SecondaryBroadcast(j_xav, [f"{domain} electrode"])
+        else:
+            j_xav = j
 
         variables = {
             f"{Domain} electrode {reaction_name}"

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

@@ -76,7 +76,9 @@ def get_coupled_variables(self, variables):
             self.reaction == "lithium-ion main"
             and domain_options["particle size"] == "distribution"
         ):
-            delta_phi = pybamm.PrimaryBroadcast(delta_phi, [f"{domain} particle size"])
+            delta_phi = pybamm.PrimaryBroadcast(
+                delta_phi, [f"{domain} {phase_name}particle size"]
+            )
 
         # Get exchange-current density. For MSMR models we calculate the exchange
         # current density for each reaction, then sum these to give a total exchange
@@ -169,7 +171,7 @@ def get_coupled_variables(self, variables):
         elif j0.domain in ["current collector", ["current collector"]]:
             T = variables["X-averaged cell temperature [K]"]
             u = variables[f"X-averaged {domain} electrode interface utilisation"]
-        elif j0.domain == [f"{domain} particle size"]:
+        elif j0.domain == [f"{domain} {phase_name}particle size"]:
             if j0.domains["secondary"] != [f"{domain} electrode"]:
                 T = variables["X-averaged cell temperature [K]"]
                 u = variables[f"X-averaged {domain} electrode interface utilisation"]
@@ -178,7 +180,7 @@ def get_coupled_variables(self, variables):
                 u = variables[f"{Domain} electrode interface utilisation"]
 
             # Broadcast T onto "particle size" domain
-            T = pybamm.PrimaryBroadcast(T, [f"{domain} particle size"])
+            T = pybamm.PrimaryBroadcast(T, [f"{domain} {phase_name}particle size"])
         else:
             T = variables[f"{Domain} electrode temperature [K]"]
             u = variables[f"{Domain} electrode interface utilisation"]
@@ -198,7 +200,9 @@ def get_coupled_variables(self, variables):
         else:
             j = self._get_kinetics(j0, ne, eta_r, T, u)
 
-        if j.domain == [f"{domain} particle size"]:
+        if j.domain == [f"{domain} particle size"] or j.domain == [
+            f"{domain} {phase_name}particle size"
+        ]:
             # If j depends on particle size, get size-dependent "distribution"
             # variables first
             variables.update(