Merge pull request #1130 from pybamm-team/polynomial-concentration

Polynomial concentration

Author: rtimms

CHANGELOG.md

@@ -2,9 +2,10 @@
 
 ## Features
 
+-   Added particle submodel that use a polynomial approximation to the concentration within the electrode particles ([#1130](https://github.com/pybamm-team/PyBaMM/pull/1130))
 -   Added Modulo, Floor and Ceiling operators ([#1121](https://github.com/pybamm-team/PyBaMM/pull/1121))
 -   Added DFN model for a half cell ([#1121](https://github.com/pybamm-team/PyBaMM/pull/1121))
--   Automatically compute surface area per unit volume based on particle shape for li-ion models ([#1120])(https://github.com/pybamm-team/PyBaMM/pull/1120)
+-   Automatically compute surface area per unit volume based on particle shape for li-ion models ([#1120](https://github.com/pybamm-team/PyBaMM/pull/1120))
 -   Added "R-averaged particle concentration" variables ([#1118](https://github.com/pybamm-team/PyBaMM/pull/1118))
 -   Added support for sensitivity calculations to the casadi solver ([#1109](https://github.com/pybamm-team/PyBaMM/pull/1109))
 -   Added support for index 1 semi-explicit dae equations and sensitivity calculations to JAX BDF solver ([#1107](https://github.com/pybamm-team/PyBaMM/pull/1107))
@@ -22,9 +23,10 @@
 
 ## Breaking changes
 
--  The modules containing standard parameters are now classes so they can take options
-(e.g. `standard_parameters_lithium_ion` is now `LithiumIonParameters`) ([#1120])(https://github.com/pybamm-team/PyBaMM/pull/1120)
--  Renamed `quick_plot_vars` to `output_variables` in `Simulation` to be consistent with `QuickPlot`. Passing `quick_plot_vars` to `Simulation.plot()` has been deprecated and `output_variables` should be passed instead ([#1099](https://github.com/pybamm-team/PyBaMM/pull/1099))
+-   The "fast diffusion" particle option has been renamed "uniform profile" ([#1130](https://github.com/pybamm-team/PyBaMM/pull/1130)) 
+-   The modules containing standard parameters are now classes so they can take options
+(e.g. `standard_parameters_lithium_ion` is now `LithiumIonParameters`) ([#1120](https://github.com/pybamm-team/PyBaMM/pull/1120))
+-   Renamed `quick_plot_vars` to `output_variables` in `Simulation` to be consistent with `QuickPlot`. Passing `quick_plot_vars` to `Simulation.plot()` has been deprecated and `output_variables` should be passed instead ([#1099](https://github.com/pybamm-team/PyBaMM/pull/1099))
 
 
 # [v0.2.3](https://github.com/pybamm-team/PyBaMM/tree/v0.2.3) - 2020-07-01

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

@@ -7,5 +7,5 @@ Particle
   base_particle
   fickian_single_particle
   fickian_many_particles
-  fast_single_particle
-  fast_many_particles
+  polynomial_single_particle
+  polynomial_many_particles

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

@@ -0,0 +1,5 @@
+Polynomial Many Particles
+=========================
+
+.. autoclass:: pybamm.particle.PolynomialManyParticles
+    :members:

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

@@ -0,0 +1,5 @@
+Polynomial Single Particle
+===========================
+
+.. autoclass:: pybamm.particle.PolynomialSingleParticle
+    :members:

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

@@ -41,8 +41,8 @@ For more advanced usage, new sets of parameters, spatial methods and solvers can
 
 PyBaMM is built around an expression tree structure.
 
-- [The expression tree notebook](expression_tree/expression-tree.ipynb) explains how this works, from model creation to solution. 
-- [The broadcast notebook](expression_tree/broadcasts.ipynb) explains the different types of broadcast. 
+- [The expression tree notebook](expression_tree/expression-tree.ipynb) explains how this works, from model creation to solution.
+- [The broadcast notebook](expression_tree/broadcasts.ipynb) explains the different types of broadcast.
 
 The following notebooks are specific to different stages of the PyBaMM pipeline, such as choosing a model, spatial method, or solver.
 

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

@@ -0,0 +1,35 @@
+#
+# Compare models for diffusion within the electrode particles
+#
+import pybamm
+
+# load models
+models = [
+    pybamm.lithium_ion.DFN(
+        options={"particle": "Fickian diffusion"}, name="Fickian diffusion"
+    ),
+    pybamm.lithium_ion.DFN(
+        options={"particle": "uniform profile"}, name="uniform profile"
+    ),
+    pybamm.lithium_ion.DFN(
+        options={"particle": "quadratic profile"}, name="quadratic profile"
+    ),
+    pybamm.lithium_ion.DFN(
+        options={"particle": "quartic profile"}, name="quartic profile"
+    ),
+]
+
+# pick parameter values
+parameter_values = pybamm.ParameterValues(chemistry=pybamm.parameter_sets.Ecker2015)
+
+# create and solve simulations
+sims = []
+for model in models:
+    sim = pybamm.Simulation(model, parameter_values=parameter_values)
+    sim.solve([0, 3600])
+    sims.append(sim)
+    print("Particle model: {}".format(model.name))
+    print("Solve time: {}s".format(sim.solution.solve_time))
+
+# plot results
+pybamm.dynamic_plot(sims)

examples/notebooks/README.md

@@ -22,11 +22,11 @@
 model.submodels["positive electrode"] = pybamm.electrode.ohm.LeadingOrder(
     model.param, "Positive"
 )
-model.submodels["negative particle"] = pybamm.particle.FastSingleParticle(
-    model.param, "Negative"
+model.submodels["negative particle"] = pybamm.particle.PolynomialSingleParticle(
+    model.param, "Negative", "uniform profile"
 )
-model.submodels["positive particle"] = pybamm.particle.FastSingleParticle(
-    model.param, "Positive"
+model.submodels["positive particle"] = pybamm.particle.PolynomialSingleParticle(
+    model.param, "Positive", "uniform profile"
 )
 model.submodels["negative interface"] = pybamm.interface.InverseButlerVolmer(
     model.param, "Negative", "lithium-ion main"

examples/notebooks/compare-particle-diffusion-models.ipynb

@@ -193,3 +193,14 @@ eprint={2005.05127},
 archivePrefix={arXiv},
 primaryClass={physics.app-ph},
 }
+
+@article{subramanian2005,
+  title={Efficient macro-micro scale coupled modeling of batteries},
+  author={Subramanian, Venkat R and Diwakar, Vinten D and Tapriyal, Deepak},
+  journal={Journal of The Electrochemical Society},
+  volume={152},
+  number={10},
+  pages={A2002},
+  year={2005},
+  publisher={IOP Publishing}
+}

examples/notebooks/using-submodels.ipynb

@@ -46,7 +46,8 @@ class BaseBatteryModel(pybamm.BaseModel):
                 "potential pair" or "potential pair quite conductive".
             * "particle" : str, optional
                 Sets the submodel to use to describe behaviour within the particle.
-                Can be "Fickian diffusion" (default) or "fast diffusion".
+                Can be "Fickian diffusion" (default), "uniform profile",
+                "quadratic profile", or "quartic profile".
             * "particle shape" : str, optional
                 Sets the model shape of the electrode particles. This is used to
                 calculate the surface area per unit volume. Can be "spherical"
@@ -348,10 +349,21 @@ def options(self, extra_options):
                 raise pybamm.OptionError(
                     "cannot have transverse convection in 0D model"
                 )
-        if options["particle"] not in ["Fickian diffusion", "fast diffusion"]:
+        if options["particle"] not in [
+            "Fickian diffusion",
+            "fast diffusion",
+            "uniform profile",
+            "quadratic profile",
+            "quartic profile",
+        ]:
             raise pybamm.OptionError(
                 "particle model '{}' not recognised".format(options["particle"])
             )
+        if options["particle"] == "fast diffusion":
+            raise NotImplementedError(
+                "The 'fast diffusion' option has been renamed. "
+                "Use 'uniform profile' instead."
+            )
         if options["particle shape"] not in ["spherical", "user"]:
             raise pybamm.OptionError(
                 "particle shape '{}' not recognised".format(options["particle shape"])

examples/scripts/compare_SPM_diffusion_models.py

@@ -85,12 +85,20 @@ def set_particle_submodel(self):
             self.submodels["positive particle"] = pybamm.particle.FickianManyParticles(
                 self.param, "Positive"
             )
-        elif self.options["particle"] == "fast diffusion":
-            self.submodels["negative particle"] = pybamm.particle.FastManyParticles(
-                self.param, "Negative"
+        elif self.options["particle"] in [
+            "uniform profile",
+            "quadratic profile",
+            "quartic profile",
+        ]:
+            self.submodels[
+                "negative particle"
+            ] = pybamm.particle.PolynomialManyParticles(
+                self.param, "Negative", self.options["particle"]
             )
-            self.submodels["positive particle"] = pybamm.particle.FastManyParticles(
-                self.param, "Positive"
+            self.submodels[
+                "positive particle"
+            ] = pybamm.particle.PolynomialManyParticles(
+                self.param, "Positive", self.options["particle"]
             )
 
     def set_solid_submodel(self):

examples/scripts/compare-dae-solver.py examples/scripts/compare_dae_solver.py

@@ -104,12 +104,20 @@ def set_particle_submodel(self):
             self.submodels["positive particle"] = pybamm.particle.FickianSingleParticle(
                 self.param, "Positive"
             )
-        elif self.options["particle"] == "fast diffusion":
-            self.submodels["negative particle"] = pybamm.particle.FastSingleParticle(
-                self.param, "Negative"
+        elif self.options["particle"] in [
+            "uniform profile",
+            "quadratic profile",
+            "quartic profile",
+        ]:
+            self.submodels[
+                "negative particle"
+            ] = pybamm.particle.PolynomialSingleParticle(
+                self.param, "Negative", self.options["particle"]
             )
-            self.submodels["positive particle"] = pybamm.particle.FastSingleParticle(
-                self.param, "Positive"
+            self.submodels[
+                "positive particle"
+            ] = pybamm.particle.PolynomialSingleParticle(
+                self.param, "Positive", self.options["particle"]
             )
 
     def set_negative_electrode_submodel(self):

examples/scripts/compare_particle_models.py

@@ -106,12 +106,20 @@ def set_particle_submodel(self):
             self.submodels["positive particle"] = pybamm.particle.FickianSingleParticle(
                 self.param, "Positive"
             )
-        elif self.options["particle"] == "fast diffusion":
-            self.submodels["negative particle"] = pybamm.particle.FastSingleParticle(
-                self.param, "Negative"
+        elif self.options["particle"] in [
+            "uniform profile",
+            "quadratic profile",
+            "quartic profile",
+        ]:
+            self.submodels[
+                "negative particle"
+            ] = pybamm.particle.PolynomialSingleParticle(
+                self.param, "Negative", self.options["particle"]
             )
-            self.submodels["positive particle"] = pybamm.particle.FastSingleParticle(
-                self.param, "Positive"
+            self.submodels[
+                "positive particle"
+            ] = pybamm.particle.PolynomialSingleParticle(
+                self.param, "Positive", self.options["particle"]
             )
 
     def set_negative_electrode_submodel(self):

examples/scripts/create-model.py examples/scripts/create_model.py

@@ -134,6 +134,28 @@
     auxiliary_domains={"secondary": "current collector"},
     bounds=(0, 1),
 )
+c_s_n_rav = pybamm.Variable(
+    "R-averaged negative particle concentration",
+    domain="negative electrode",
+    auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, 1),
+)
+c_s_p_rav = pybamm.Variable(
+    "R-averaged positive particle concentration",
+    domain="positive electrode",
+    auxiliary_domains={"secondary": "current collector"},
+    bounds=(0, 1),
+)
+c_s_n_rxav = pybamm.Variable(
+    "R-X-averaged negative particle concentration",
+    domain="current collector",
+    bounds=(0, 1),
+)
+c_s_p_rxav = pybamm.Variable(
+    "R-X-averaged positive particle concentration",
+    domain="current collector",
+    bounds=(0, 1),
+)
 c_s_n_surf = pybamm.Variable(
     "Negative particle surface concentration",
     domain="negative electrode",
@@ -156,7 +178,25 @@
     domain="current collector",
     bounds=(0, 1),
 )
-
+# Average particle concentration gradient (for polynomial particle concentration
+# models). Note: we make the distinction here between the flux defined as
+# N = -D*dc/dr and the concentration gradient q = dc/dr
+q_s_n_rav = pybamm.Variable(
+    "R-averaged negative particle concentration gradient",
+    domain="negative electrode",
+    auxiliary_domains={"secondary": "current collector"},
+)
+q_s_p_rav = pybamm.Variable(
+    "R-averaged positive particle concentration gradient",
+    domain="positive electrode",
+    auxiliary_domains={"secondary": "current collector"},
+)
+q_s_n_rxav = pybamm.Variable(
+    "R-X-averaged negative particle concentration gradient", domain="current collector"
+)
+q_s_p_rxav = pybamm.Variable(
+    "R-X-averaged positive particle concentration gradient", domain="current collector"
+)
 
 # Porosity
 eps_n = pybamm.Variable(
@@ -181,13 +221,13 @@
 
 # Piecewise constant (for asymptotic models)
 eps_n_pc = pybamm.Variable(
-    "X-averaged negative electrode porosity", domain="current collector", bounds=(0, 1),
+    "X-averaged negative electrode porosity", domain="current collector", bounds=(0, 1)
 )
 eps_s_pc = pybamm.Variable(
     "X-averaged separator porosity", domain="current collector", bounds=(0, 1)
 )
 eps_p_pc = pybamm.Variable(
-    "X-averaged positive electrode porosity", domain="current collector", bounds=(0, 1),
+    "X-averaged positive electrode porosity", domain="current collector", bounds=(0, 1)
 )
 
 eps_piecewise_constant = pybamm.Concatenation(

examples/scripts/custom_model.py

@@ -1,5 +1,5 @@
 from .base_particle import BaseParticle
 from .fickian_many_particles import FickianManyParticles
 from .fickian_single_particle import FickianSingleParticle
-from .fast_many_particles import FastManyParticles
-from .fast_single_particle import FastSingleParticle
+from .polynomial_single_particle import PolynomialSingleParticle
+from .polynomial_many_particles import PolynomialManyParticles