Many-Particle Models

docs/source/expression_tree/unary_operator.rst

@@ -81,6 +81,8 @@ Unary Operators
 
 .. autofunction:: pybamm.r_average
 
+.. autofunction:: pybamm.size_average
+
 .. autofunction:: pybamm.z_average
 
 .. autofunction:: pybamm.yz_average

docs/source/models/lithium_ion/index.rst

@@ -6,6 +6,7 @@ Lithium-ion Models
   base_lithium_ion_model
   spm
   spme
+  mpm
   dfn
   newman_tobias
   yang2017

docs/source/models/lithium_ion/mpm.rst

@@ -0,0 +1,5 @@
+Many Particle Model (MPM)
+===========================
+
+.. autoclass:: pybamm.lithium_ion.MPM
+    :members:

docs/source/models/submodels/particle/index.rst

@@ -9,3 +9,4 @@ Particle
   fickian_many_particles
   polynomial_single_particle
   polynomial_many_particles
+  size_distribution/index

docs/source/models/submodels/particle/size_distribution/base_distribution.rst

@@ -0,0 +1,6 @@
+Particle Size Distribution Base Model
+=====================================
+
+.. autoclass:: pybamm.particle.BaseSizeDistribution
+    :members:
+

docs/source/models/submodels/particle/size_distribution/fast_single_distribution.rst

@@ -0,0 +1,5 @@
+Fast Single Size Distribution
+=============================
+
+.. autoclass:: pybamm.particle.FastSingleSizeDistribution
+    :members:

docs/source/models/submodels/particle/size_distribution/fickian_single_distribution.rst

@@ -0,0 +1,6 @@
+Fickian Single Size Distribution
+================================
+
+.. autoclass:: pybamm.particle.FickianSingleSizeDistribution
+    :members:
+

docs/source/models/submodels/particle/size_distribution/index.rst

@@ -0,0 +1,9 @@
+Particle Size Distribution
+==========================
+
+.. toctree::
+  :maxdepth: 1
+
+  base_distribution
+  fickian_single_distribution
+  fast_single_distribution

pybamm/CITATIONS.txt

@@ -405,6 +405,31 @@ doi={10.1149/2.0661810jes}
   publisher={Elsevier}
 }
 
+@misc{Kirk2020,
+  author = {Kirk, Toby L. and Evans, Jack and Please, Colin P. and Chapman, S. Jonathan},
+  title         = {Modelling electrode heterogeneity in lithium-ion batteries: unimodal and bimodal particle-size distributions},
+  year          = {2020},
+  howpublished  = {arXiv:2006.12208},
+  eprint        = {2006.12208},
+  url           = {https://arxiv.org/abs/2006.12208},
+  archiveprefix = {arXiv},
+  primaryclass  = {physics.app-ph},
+}
+
+@article{Kirk2021,
+	author = {Toby L. Kirk and Colin P. Please and S. Jon Chapman},
+	title = {Physical Modelling of the Slow Voltage Relaxation Phenomenon in Lithium-Ion Batteries},
+	journal = {Journal of The Electrochemical Society},
+	year = 2021,
+	month = {jun},
+	publisher = {The Electrochemical Society},
+	volume = {168},
+	number = {6},
+	pages = {060554},
+	doi = {10.1149/1945-7111/ac0bf7},
+	url = {https://doi.org/10.1149/1945-7111/ac0bf7},
+}
+
 @article{Xu2019,
   title={Evolution of Dead Lithium Growth in Lithium Metal Batteries: Experimentally Validated Model of the Apparent Capacity Loss},
   author={Xu, Shanshan and Chen, Kuan-Hung and Dasgupta, Neil P and Siegel, Jason B and Stefanopoulou, Anna G},

pybamm/__init__.py

@@ -173,6 +173,7 @@ def version(formatted=False):
 from .parameters.thermal_parameters import thermal_parameters, ThermalParameters
 from .parameters.lithium_ion_parameters import LithiumIonParameters
 from .parameters.lead_acid_parameters import LeadAcidParameters
+from .parameters.size_distribution_parameters import *
 from .parameters import parameter_sets
 
 #

pybamm/expression_tree/broadcasts.py

@@ -100,29 +100,49 @@ def check_and_set_domains(
         self, child, broadcast_type, broadcast_domain, broadcast_auxiliary_domains
     ):
         """See :meth:`Broadcast.check_and_set_domains`"""
-        # Can only do primary broadcast from current collector to electrode or particle
-        # or from electrode to particle. Note current collector to particle *is* allowed
+        # Can only do primary broadcast from current collector to electrode,
+        # particle-size or particle or from electrode to particle-size or particle.
+        # Note e.g. current collector to particle *is* allowed
         if child.domain == []:
             pass
         elif child.domain == ["current collector"] and broadcast_domain[0] not in [
             "negative electrode",
             "separator",
             "positive electrode",
+            "negative particle size",
+            "positive particle size",
             "negative particle",
             "positive particle",
         ]:
             raise pybamm.DomainError(
                 """Primary broadcast from current collector domain must be to electrode
-                or separator or particle domains"""
+                or separator or particle or particle size domains"""
+            )
+        elif (
+            child.domain[0]
+            in [
+                "negative electrode",
+                "separator",
+                "positive electrode",
+            ]
+            and broadcast_domain[0] not in [
+                "negative particle",
+                "positive particle",
+                "negative particle size",
+                "positive particle size",
+            ]
+        ):
+            raise pybamm.DomainError(
+                """Primary broadcast from electrode or separator must be to particle
+                or particle size domains"""
             )
         elif child.domain[0] in [
-            "negative electrode",
-            "separator",
-            "positive electrode",
+            "negative particle size",
+            "positive particle size",
         ] and broadcast_domain[0] not in ["negative particle", "positive particle"]:
             raise pybamm.DomainError(
-                """Primary broadcast from electrode or separator must be to particle
-                domains"""
+                """Primary broadcast from particle size domain must be to particle
+                domain"""
             )
         elif child.domain[0] in ["negative particle", "positive particle"]:
             raise pybamm.DomainError("Cannot do primary broadcast from particle domain")
@@ -207,14 +227,29 @@ def check_and_set_domains(
         if child.domain[0] in [
             "negative particle",
             "positive particle",
+        ] and broadcast_domain[0] not in [
+            "negative particle size",
+            "positive particle size",
+            "negative electrode",
+            "separator",
+            "positive electrode",
+        ]:
+            raise pybamm.DomainError(
+                """Secondary broadcast from particle domain must be to particle-size,
+                electrode or separator domains"""
+            )
+        if child.domain[0] in [
+            "negative particle size",
+            "positive particle size",
         ] and broadcast_domain[0] not in [
             "negative electrode",
             "separator",
             "positive electrode",
+            "current collector"
         ]:
             raise pybamm.DomainError(
-                """Secondary broadcast from particle domain must be to electrode or
-                separator domains"""
+                """Secondary broadcast from particle size domain must be to
+                electrode or separator or current collector domains"""
             )
         elif child.domain[0] in [
             "negative electrode",

pybamm/expression_tree/unary_operators.py

@@ -1258,10 +1258,16 @@ def x_average(symbol):
         # Even if domain is "negative electrode", "separator", or
         # "positive electrode", and we know l, we still compute it as Integral(1, x)
         # as this will be easier to identify for simplifications later on
-        if symbol.domain == ["negative particle"]:
+        if (
+            symbol.domain == ["negative particle"] or
+            symbol.domain == ["negative particle size"]
+        ):
             x = pybamm.standard_spatial_vars.x_n
             l = geo.l_n
-        elif symbol.domain == ["positive particle"]:
+        elif (
+            symbol.domain == ["positive particle"] or
+            symbol.domain == ["positive particle size"]
+        ):
             x = pybamm.standard_spatial_vars.x_p
             l = geo.l_p
         else:
@@ -1407,6 +1413,67 @@ def r_average(symbol):
         return Integral(symbol, r) / Integral(v, r)
 
 
+def size_average(symbol):
+    """convenience function for averaging over particle size R using the area-weighted
+    particle-size distribution.
+
+    Parameters
+    ----------
+    symbol : :class:`pybamm.Symbol`
+        The function to be averaged
+    Returns
+    -------
+    :class:`Symbol`
+        the new averaged symbol
+    """
+    # Can't take average if the symbol evaluates on edges
+    if symbol.evaluates_on_edges("primary"):
+        raise ValueError(
+            """Can't take the size-average of a symbol that evaluates on edges"""
+        )
+
+    # 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"]]
+        for domain in list(symbol.domains.values())
+    ):
+        new_symbol = symbol.new_copy()
+        new_symbol.parent = None
+        return new_symbol
+
+    # If symbol is a primary broadcast to "particle size", take the orphan
+    elif isinstance(symbol, pybamm.PrimaryBroadcast) and symbol.domain in [
+        ["negative particle size"], ["positive particle size"]
+    ]:
+        return symbol.orphans[0]
+    # If symbol is a secondary broadcast to "particle size" from "particle",
+    # take the orphan
+    elif (
+        isinstance(symbol, pybamm.SecondaryBroadcast) and
+        symbol.domains["secondary"] in [
+            ["negative particle size"], ["positive particle size"]
+        ]
+    ):
+        return symbol.orphans[0]
+    # Otherwise, perform the integration in R
+    else:
+        R = pybamm.SpatialVariable(
+            "R",
+            domain=symbol.domain,
+            auxiliary_domains=symbol.auxiliary_domains,
+            coord_sys="cartesian",
+        )
+        geo = pybamm.geometric_parameters
+        if ["negative particle size"] in list(symbol.domains.values()):
+            f_a_dist = geo.f_a_dist_n(R)
+        elif ["positive particle size"] in list(symbol.domains.values()):
+            f_a_dist = geo.f_a_dist_p(R)
+
+        # take average using Integral and distribution f_a_dist
+        return Integral(f_a_dist * symbol, R) / Integral(f_a_dist, R)
+
+
 def boundary_value(symbol, side):
     """
     convenience function for creating a :class:`pybamm.BoundaryValue`

pybamm/geometry/battery_geometry.py

@@ -4,14 +4,21 @@
 import pybamm
 
 
-def battery_geometry(include_particles=True, current_collector_dimension=0):
+def battery_geometry(
+    include_particles=True,
+    options=None,
+    current_collector_dimension=0,
+):
     """
     A convenience function to create battery geometries.
 
     Parameters
     ----------
     include_particles : bool
         Whether to include particle domains
+    options : dict
+        Dictionary of model options. Necessary for "particle-size geometry",
+        relevant for lithium-ion chemistries.
     current_collector_dimensions : int, default
         The dimensions of the current collector. Should be 0 (default), 1 or 2
 
@@ -39,6 +46,25 @@ def battery_geometry(include_particles=True, current_collector_dimension=0):
                 "positive particle": {var.r_p: {"min": 0, "max": 1}},
             }
         )
+    # Add particle size domains
+    if (
+        options is not None and
+        options["particle size"] == "distribution"
+    ):
+        R_min_n = geo.R_min_n
+        R_min_p = geo.R_min_p
+        R_max_n = geo.R_max_n
+        R_max_p = geo.R_max_p
+        geometry.update(
+            {
+                "negative particle size": {
+                    var.R_n: {"min": R_min_n, "max": R_max_n}
+                },
+                "positive particle size": {
+                    var.R_p: {"min": R_min_p, "max": R_max_p}
+                },
+            }
+        )
 
     if current_collector_dimension == 0:
         geometry["current collector"] = {var.z: {"position": 1}}

pybamm/geometry/standard_spatial_vars.py

@@ -50,6 +50,25 @@
     coord_sys="spherical polar",
 )
 
+R_n = pybamm.SpatialVariable(
+    "R_n",
+    domain=["negative particle size"],
+    auxiliary_domains={
+        "secondary": "negative electrode",
+        "tertiary": "current collector",
+    },
+    coord_sys="cartesian",
+)
+R_p = pybamm.SpatialVariable(
+    "R_p",
+    domain=["positive particle size"],
+    auxiliary_domains={
+        "secondary": "positive electrode",
+        "tertiary": "current collector",
+    },
+    coord_sys="cartesian",
+)
+
 # Domains at cell edges
 x_n_edge = pybamm.SpatialVariableEdge(
     "x_n",
@@ -101,3 +120,22 @@
     },
     coord_sys="spherical polar",
 )
+
+R_n_edge = pybamm.SpatialVariableEdge(
+    "R_n",
+    domain=["negative particle size"],
+    auxiliary_domains={
+        "secondary": "negative electrode",
+        "tertiary": "current collector",
+    },
+    coord_sys="cartesian",
+)
+R_p_edge = pybamm.SpatialVariableEdge(
+    "R_p",
+    domain=["positive particle size"],
+    auxiliary_domains={
+        "secondary": "positive electrode",
+        "tertiary": "current collector",
+    },
+    coord_sys="cartesian",
+)