#1748 add base fickian model

pybamm/models/full_battery_models/base_battery_model.py

@@ -124,6 +124,11 @@ class BatteryModelOptions(pybamm.FuzzyDict):
             * "SEI porosity change" : str
                 Whether to include porosity change due to SEI formation, can be "false"
                 (default) or "true".
+            * "stress induced diffusion" : str
+                Whether to includes stress induced diffusion, can be "false" or "true".
+                The default is "false" if "particle mechanics" is "none" and "true"
+                otherwise. A 2-tuple can be provided for different behaviour in negative
+                and positive electrodes.
             * "surface form" : str
                 Whether to use the surface formulation of the problem. Can be "false"
                 (default), "differential" or "algebraic".
@@ -186,6 +191,7 @@ def __init__(self, extra_options):
             ],
             "SEI film resistance": ["none", "distributed", "average"],
             "SEI porosity change": ["true", "false"],
+            "stress induced diffusion": ["true", "false"],
             "surface form": ["false", "differential", "algebraic"],
             "thermal": ["isothermal", "lumped", "x-lumped", "x-full"],
             "total interfacial current density as a state": ["true", "false"],
@@ -238,17 +244,40 @@ def __init__(self, extra_options):
         # The "SEI film resistance" option will still be overridden by extra_options if
         # provided
 
-        # Change the default for swelling based on which LAM option is
+        # Change the default for particle mechanics based on which LAM option is
         # provided
         # return "none" if option not given
         lam_option = extra_options.get("loss of active material", "none")
         if "stress-driven" in lam_option:
             default_options["particle mechanics"] = "swelling only"
         else:
             default_options["particle mechanics"] = "none"
-        # The "SEI film resistance" option will still be overridden by extra_options if
+        # The "particle mechanics" option will still be overridden by extra_options if
         # provided
 
+        # Change the default for stress induced diffusion based on which particle
+        # mechanics option is provided
+        mechanics_options = extra_options.get("particle mechanics", "none")
+        if mechanics_options == "none":
+            default_options["stress induced diffusion"] = "false"
+        else:
+            if isinstance(self.options["particle mechanics"], str):
+                mech_left = self.options["particle mechanics"]
+                mech_right = self.options["particle mechanics"]
+            else:
+                mech_left, mech_right = self.options["particle mechanics"]
+            if mech_left == "none":
+                default_left = "false"
+            else:
+                default_left = "true"
+            if mech_right == "none":
+                default_right = "false"
+            else:
+                default_right = "true"
+            default_options["stress induced diffusion"] = (default_left, default_right)
+        # The "stress induced diffusion" option will still be overridden by
+        # extra_options if provided
+
         options = pybamm.FuzzyDict(default_options)
         # any extra options overwrite the default options
         for name, opt in extra_options.items():
@@ -283,10 +312,6 @@ def __init__(self, extra_options):
 
         # Options not yet compatible with particle-size distributions
         if options["particle size"] == "distribution":
-            if options["SEI"] != "none":
-                raise NotImplementedError(
-                    "SEI submodels do not yet support particle-size distributions."
-                )
             if options["lithium plating"] != "none":
                 raise NotImplementedError(
                     "Lithium plating submodels do not yet support particle-size "
@@ -299,9 +324,18 @@ def __init__(self, extra_options):
                 )
             if options["particle shape"] != "spherical":
                 raise NotImplementedError(
-                    "Particle shape must be 'spherical' for particle-size distributions"
+                    "Particle shape must be 'spherical' for particle-size distribution"
                     " submodels."
                 )
+            if options["SEI"] != "none":
+                raise NotImplementedError(
+                    "SEI submodels do not yet support particle-size distributions."
+                )
+            if options["stress induced diffusion"] == "true":
+                raise NotImplementedError(
+                    "Stress induced diffusion cannot yet be included in "
+                    "particle-size distributions."
+                )
             if options["thermal"] == "x-full":
                 raise NotImplementedError(
                     "X-full thermal submodels do not yet support particle-size"

pybamm/models/standard_variables.py

@@ -175,12 +175,12 @@ def __init__(self):
             auxiliary_domains={"secondary": "current collector"},
             bounds=(0, 1),
         )
-        self.c_s_n_rxav = pybamm.Variable(
+        self.c_s_n_av = pybamm.Variable(
             "Average negative particle concentration",
             domain="current collector",
             bounds=(0, 1),
         )
-        self.c_s_p_rxav = pybamm.Variable(
+        self.c_s_p_av = pybamm.Variable(
             "Average positive particle concentration",
             domain="current collector",
             bounds=(0, 1),
@@ -220,11 +220,11 @@ def __init__(self):
             domain="positive electrode",
             auxiliary_domains={"secondary": "current collector"},
         )
-        self.q_s_n_rxav = pybamm.Variable(
+        self.q_s_n_av = pybamm.Variable(
             "Average negative particle concentration gradient",
             domain="current collector",
         )
-        self.q_s_p_rxav = pybamm.Variable(
+        self.q_s_p_av = pybamm.Variable(
             "Average positive particle concentration gradient",
             domain="current collector",
         )

pybamm/models/submodels/particle/base_fickian.py

@@ -0,0 +1,63 @@
+#
+# Base class for particles with Fickian diffusion
+#
+import pybamm
+from base_particle import BaseParticle
+
+
+class BaseFickian(BaseParticle):
+    """
+    Base class for molar conservation in particles with Fickian diffusion.
+
+    Parameters
+    ----------
+    param : parameter class
+        The parameters to use for this submodel
+    domain : str
+        The domain of the model either 'Negative' or 'Positive'
+    options: dict
+        A dictionary of options to be passed to the model.
+        See :class:`pybamm.BaseBatteryModel`
+
+    **Extends:** :class:`pybamm.particle.BaseParticle`
+    """
+
+    def __init__(self, param, domain, options):
+        super().__init__(param, domain, options)
+
+        if self.options["stress induced diffusion"] == "true":
+            pybamm.citations.register("Ai2019")
+            pybamm.citations.register("Deshpande2012")
+
+    def _get_effective_diffusivity(self, c, T):
+        if self.domain == "Negative":
+            D = self.param.D_n(c, T)
+        elif self.domain == "Positive":
+            D = self.param.D_p(c, T)
+
+        # Account for stress induced diffusion
+        if self.options["stress induced diffusion"] == "true":
+            if self.domain == "Negative":
+                theta = self.param.theta_n
+                c_0 = self.param.c_0_n
+            elif self.domain == "Positive":
+                theta = self.param.theta_p
+                c_0 = self.param.c_0_p
+
+            D_eff = D * (1 + theta * (c - c_0) / (1 + self.param.Theta * T))
+        else:
+            D_eff = D
+
+        return D_eff
+
+    def _get_standard_flux_variables(self, N_s, N_s_xav, D_eff):
+        variables = {
+            self.domain + " particle flux": N_s,
+            "X-averaged " + self.domain.lower() + " particle flux": N_s_xav,
+            self.domain + " effective diffusivity": D_eff,
+            "X-averaged "
+            + self.domain.lower()
+            + " effective diffusivity": pybamm.x_average(D_eff),
+        }
+
+        return variables

pybamm/models/submodels/particle/base_particle.py

@@ -14,12 +14,15 @@ class BaseParticle(pybamm.BaseSubModel):
         The parameters to use for this submodel
     domain : str
         The domain of the model either 'Negative' or 'Positive'
+    options: dict
+        A dictionary of options to be passed to the model.
+        See :class:`pybamm.BaseBatteryModel`
 
     **Extends:** :class:`pybamm.BaseSubModel`
     """
 
-    def __init__(self, param, domain):
-        super().__init__(param, domain)
+    def __init__(self, param, domain, options=None):
+        super().__init__(param, domain, options=options)
 
     def _get_standard_concentration_variables(
         self, c_s, c_s_xav=None, c_s_rav=None, c_s_av=None, c_s_surf=None
@@ -148,4 +151,3 @@ def _get_standard_flux_variables(self, N_s, N_s_xav):
         }
 
         return variables
-

pybamm/models/submodels/particle/no_distribution/fickian_diffusion.py

@@ -2,10 +2,10 @@
 # Class for particles with Fickian diffusion and x-dependence
 #
 import pybamm
-from ..base_particle import BaseParticle
+from ..base_fickian import BaseFickian
 
 
-class FickianDiffusion(BaseParticle):
+class FickianDiffusion(BaseFickian):
     """
     Class for molar conservation in particles, employing Fick's law, and allowing
     variation in the electrode domain. I.e., the concentration varies with r
@@ -17,12 +17,15 @@ class FickianDiffusion(BaseParticle):
         The parameters to use for this submodel
     domain : str
         The domain of the model either 'Negative' or 'Positive'
+    options: dict
+        A dictionary of options to be passed to the model.
+        See :class:`pybamm.BaseBatteryModel`
 
     **Extends:** :class:`pybamm.particle.BaseParticle`
     """
 
-    def __init__(self, param, domain):
-        super().__init__(param, domain)
+    def __init__(self, param, domain, options):
+        super().__init__(param, domain, options)
 
     def get_fundamental_variables(self):
         if self.domain == "Negative":
@@ -42,13 +45,10 @@ def get_coupled_variables(self, variables):
             [self.domain.lower() + " particle"],
         )
 
-        if self.domain == "Negative":
-            N_s = -self.param.D_n(c_s, T) * pybamm.grad(c_s)
-
-        elif self.domain == "Positive":
-            N_s = -self.param.D_p(c_s, T) * pybamm.grad(c_s)
+        D_eff = self._get_effective_diffusivity(c_s, T)
+        N_s = -D_eff * pybamm.grad(c_s)
 
-        variables.update(self._get_standard_flux_variables(N_s, N_s))
+        variables.update(self._get_standard_flux_variables(N_s, N_s, D_eff))
         variables.update(self._get_total_concentration_variables(variables))
 
         return variables
@@ -67,21 +67,12 @@ def set_rhs(self, variables):
     def set_boundary_conditions(self, variables):
 
         c_s = variables[self.domain + " particle concentration"]
-        T = pybamm.PrimaryBroadcast(
-            variables[self.domain + " electrode temperature"],
-            c_s.domain,
-        )
+        D_eff = variables[self.domain + " effective diffusivity"]
         j = variables[self.domain + " electrode interfacial current density"]
         R = variables[self.domain + " particle radius"]
 
         if self.domain == "Negative":
-            rbc = (
-                -self.param.C_n
-                * j
-                * R
-                / self.param.a_R_n
-                / pybamm.surf(self.param.D_n(c_s, T))
-            )
+            rbc = -self.param.C_n * j * R / self.param.a_R_n / pybamm.surf(D_eff)
 
         elif self.domain == "Positive":
             rbc = (
@@ -90,7 +81,7 @@ def set_boundary_conditions(self, variables):
                 * R
                 / self.param.a_R_p
                 / self.param.gamma_p
-                / pybamm.surf(self.param.D_p(c_s, T))
+                / pybamm.surf(D_eff)
             )
 
         self.boundary_conditions = {