Merge pull request #1529 from tobykirk/add-mpm

Many-Particle Models

CHANGELOG.md

@@ -3,7 +3,10 @@
 ## Features
 
 -   Added temperature dependence on electrode electronic conductivity ([#1570](https://github.com/pybamm-team/PyBaMM/pull/1570))
--   Added fitted expressions for OCPs for the Chen2020 parameter set ([#1526](https://github.com/pybamm-team/PyBaMM/pull/1526))
+-   Added a new lithium-ion model `MPM` or Many-Particle Model, with a distribution of particle sizes in each electrode. ([#1529](https://github.com/pybamm-team/PyBaMM/pull/1529))
+-   Added 2 new submodels for lithium transport in a size distribution of electrode particles: Fickian diffusion (`FickianSingleSizeDistribution`) and uniform concentration profile (`FastSingleSizeDistribution`). ([#1529](https://github.com/pybamm-team/PyBaMM/pull/1529))
+-   Added a "particle size" domain to the default lithium-ion geometry, including plotting capabilities (`QuickPlot`) and processing of variables (`ProcessedVariable`). ([#1529](https://github.com/pybamm-team/PyBaMM/pull/1529))
+-   Added fitted expressions for OCPs for the Chen2020 parameter set ([#1526](https://github.com/pybamm-team/PyBaMM/pull/1497))
 -   Added `initial_soc` argument to `Simualtion.solve` for specifying the initial SOC when solving a model ([#1512](https://github.com/pybamm-team/PyBaMM/pull/1512))
 -   Added `print_name` to some symbols ([#1495](https://github.com/pybamm-team/PyBaMM/pull/1495), [#1497](https://github.com/pybamm-team/PyBaMM/pull/1497))
 -   Added Base Parameters class and SymPy in dependencies ([#1495](https://github.com/pybamm-team/PyBaMM/pull/1495))

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

examples/notebooks/using-submodels.ipynb

@@ -403,7 +403,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "In the Single Particle Model, the overpotential can be obtianed by inverting the Butler-Volmer relation, so we choose the `InverseButlerVolmer` submodel for the interface, with the \"main\" lithium-ion reaction. Because of how the current is implemented, we also need to separately specify the `CurrentForInverseButlerVolmer` submodel"
+    "In the Single Particle Model, the overpotential can be obtianed by inverting the Butler-Volmer relation, so we choose the `InverseButlerVolmer` submodel for the interface, with the \"main\" lithium-ion reaction (and default lithium ion options). Because of how the current is implemented, we also need to separately specify the `CurrentForInverseButlerVolmer` submodel"
    ]
   },
   {
@@ -414,10 +414,14 @@
    "source": [
     "model.submodels[\n",
     "    \"negative interface\"\n",
-    "] = pybamm.interface.InverseButlerVolmer(model.param, \"Negative\", \"lithium-ion main\")\n",
+    "] = pybamm.interface.InverseButlerVolmer(\n",
+    "    model.param, \"Negative\", \"lithium-ion main\", options=model.options\n",
+    ")\n",
     "model.submodels[\n",
     "    \"positive interface\"\n",
-    "] = pybamm.interface.InverseButlerVolmer(model.param, \"Positive\", \"lithium-ion main\")\n",
+    "] = pybamm.interface.InverseButlerVolmer(\n",
+    "    model.param, \"Positive\", \"lithium-ion main\", options=model.options\n",
+    ")\n",
     "model.submodels[\n",
     "    \"negative interface current\"\n",
     "] = pybamm.interface.CurrentForInverseButlerVolmer(\n",
@@ -540,23 +544,14 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
-   "language": "python",
-   "name": "python3"
+   "display_name": "Python 2.7.17 64-bit",
+   "name": "python375jvsc74a57bd0fd69f43f58546b570e94fd7eba7b65e6bcc7a5bbc4eab0408017d18902915d69"
   },
   "language_info": {
-   "codemirror_mode": {
-    "name": "ipython",
-    "version": 3
-   },
-   "file_extension": ".py",
-   "mimetype": "text/x-python",
    "name": "python",
-   "nbconvert_exporter": "python",
-   "pygments_lexer": "ipython3",
-   "version": "3.7.5"
+   "version": ""
   }
  },
  "nbformat": 4,
  "nbformat_minor": 2
-}
+}

examples/scripts/custom_model.py

@@ -36,10 +36,10 @@
     model.param, "Positive", "uniform profile"
 )
 model.submodels["negative interface"] = pybamm.interface.InverseButlerVolmer(
-    model.param, "Negative", "lithium-ion main"
+    model.param, "Negative", "lithium-ion main", options=model.options
 )
 model.submodels["positive interface"] = pybamm.interface.InverseButlerVolmer(
-    model.param, "Positive", "lithium-ion main"
+    model.param, "Positive", "lithium-ion main", options=model.options
 )
 model.submodels[
     "negative interface current"

pybamm/CITATIONS.txt

@@ -411,6 +411,31 @@
   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

@@ -176,6 +176,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

@@ -104,29 +104,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")
@@ -211,14 +231,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

@@ -1272,10 +1272,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:
@@ -1421,6 +1427,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`