#3161 add emperical hysteresis for exchange-current and diffusivity

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+- Added option to use an empirical hysteresis model for the diffusivity and exchange-current density ([#3194](https://github.com/pybamm-team/PyBaMM/pull/3194))
 - Double-layer capacity can now be provided as a function of temperature ([#3174](https://github.com/pybamm-team/PyBaMM/pull/3174))
 
 ## Bug fixes

examples/scripts/emperical_hysteresis.py

@@ -0,0 +1,151 @@
+#
+# Empirical hysteresis modelling
+#
+import pybamm
+from numpy import array
+
+pybamm.set_logging_level("INFO")
+
+# Load model
+options = {
+    "open-circuit potential": ("current sigmoid", "single"),
+    "exchange-current density": ("current sigmoid", "single"),
+    "diffusivity": ("current sigmoid", "single"),
+}
+model = pybamm.lithium_ion.SPMe(options)
+
+# Load parameter values and add (de)lithiation OCPs, exchange-current densities,
+# and diffusion coefficients for the negative electrode. Note: these are only intended
+# to be illustrative, and are not based on any particular data.
+parameter_values = pybamm.ParameterValues("Chen2020")
+
+
+def ocp_lithiation(sto):
+    R = 8.3145
+    T = 298.15
+    F = 96485
+    p = array(
+        [
+            5.03704200e02,
+            -8.82372514e-06,
+            4.56628658e02,
+            5.61890927e02,
+            3.87812964e00,
+            -4.02931829e00,
+            1.52206158e03,
+            -1.59630048e01,
+            -2.67563363e-01,
+            5.20566396e-02,
+            -6.09073875e00,
+            2.65427032e-01,
+            1.46011356e-02,
+            -1.64753973e01,
+            7.39882291e-01,
+        ]
+    )
+    u_eq = (
+        p[0] * pybamm.exp(-p[1] * sto)
+        + p[2]
+        + p[3] * pybamm.tanh(p[4] * (sto - p[5]))
+        + p[6] * pybamm.tanh(p[7] * (sto - p[8]))
+        + p[9] * pybamm.tanh(p[10] * (sto - p[11]))
+        + p[12] * pybamm.tanh(p[13] * (sto - p[14]))
+        - 0.5 * R * T / F * pybamm.log(sto)
+    )
+    return u_eq
+
+
+def ocp_delithiation(sto):
+    R = 8.3145
+    T = 298.15
+    F = 96485
+    p = array(
+        [
+            3.91897019e01,
+            2.55972942e00,
+            -3.12343803e01,
+            2.75190487e00,
+            9.23453652e00,
+            2.15796783e-02,
+            -8.81134487e01,
+            -2.12131911e00,
+            3.08859341e-01,
+            9.90149054e01,
+            -2.58589571e00,
+            4.21646546e-01,
+            3.91986770e01,
+            2.93997264e00,
+            4.68431323e-01,
+        ]
+    )
+    u_eq = (
+        p[0] * pybamm.exp(-p[1] * sto)
+        + p[2]
+        + p[3] * pybamm.tanh(p[4] * (sto - p[5]))
+        + p[6] * pybamm.tanh(p[7] * (sto - p[8]))
+        + p[9] * pybamm.tanh(p[10] * (sto - p[11]))
+        + p[12] * pybamm.tanh(p[13] * (sto - p[14]))
+        - 0.5 * R * T / F * pybamm.log(sto)
+    )
+    return u_eq
+
+
+def ocp_average(sto):
+    return (ocp_lithiation(sto) + ocp_delithiation(sto)) / 2
+
+
+def exchange_current_density_lithiation(c_e, c_s_surf, c_s_max, T):
+    m_ref = 9e-7  # (A/m2)(m3/mol)**1.5
+    return m_ref * c_e**0.5 * c_s_surf**0.5 * (c_s_max - c_s_surf) ** 0.5
+
+
+def exchange_current_density_delithiation(c_e, c_s_surf, c_s_max, T):
+    m_ref = 6e-7  # (A/m2)(m3/mol)**1.5
+    return m_ref * c_e**0.5 * c_s_surf**0.5 * (c_s_max - c_s_surf) ** 0.5
+
+
+def exchange_current_density_average(sto):
+    return (
+        exchange_current_density_lithiation(sto)
+        + exchange_current_density_delithiation(sto)
+    ) / 2
+
+
+parameter_values.update(
+    {
+        "Negative electrode OCP [V]": ocp_average,
+        "Negative electrode lithiation OCP [V]": ocp_lithiation,
+        "Negative electrode delithiation OCP [V]": ocp_delithiation,
+        "Negative electrode exchange-current density [A.m-2]"
+        "": exchange_current_density_average,
+        "Negative electrode lithiation exchange-current density [A.m-2]"
+        "": exchange_current_density_lithiation,
+        "Negative electrode delithiation exchange-current density [A.m-2]"
+        "": exchange_current_density_delithiation,
+        "Negative electrode diffusivity [m2.s-1]": 3.3e-14,
+        "Negative electrode lithiation diffusivity [m2.s-1]": 4e-14,
+        "Negative electrode delithiation diffusivity [m2.s-1]": 2.6e-14,
+    },
+    check_already_exists=False,
+)
+
+# Create experiment and run simulation
+experiment = pybamm.Experiment(
+    ["Discharge at 1 C until 2.5 V", "Charge at 1 C until 4.2 V"]
+)
+sim = pybamm.Simulation(model, experiment=experiment, parameter_values=parameter_values)
+sim.solve()
+
+# Plot
+sim.plot(
+    [
+        "X-averaged negative particle surface concentration [mol.m-3]",
+        "X-averaged positive particle surface concentration [mol.m-3]",
+        "X-averaged negative electrode exchange current density [A.m-2]",
+        "X-averaged positive electrode exchange current density [A.m-2]",
+        "X-averaged negative particle concentration [mol.m-3]",
+        "X-averaged positive particle concentration [mol.m-3]",
+        "Current [A]",
+        "Voltage [V]",
+    ]
+)

pybamm/models/full_battery_models/base_battery_model.py

@@ -43,12 +43,20 @@ class BatteryModelOptions(pybamm.FuzzyDict):
             * "current collector" : str
                 Sets the current collector model to use. Can be "uniform" (default),
                 "potential pair" or "potential pair quite conductive".
+            * "diffusivity" : str
+                Sets the model for the diffusivity. Can be "single" 
+                (default) or "current sigmoid". A 2-tuple can be provided for different 
+                behaviour in negative and positive electrodes.                    
             * "dimensionality" : int
                 Sets the dimension of the current collector problem. Can be 0
                 (default), 1 or 2.
             * "electrolyte conductivity" : str
                 Can be "default" (default), "full", "leading order", "composite" or
                 "integrated".
+            * "exchange-current density" : str
+                Sets the model for the exchange-current density. Can be "single" 
+                (default) or "current sigmoid". A 2-tuple can be provided for different 
+                behaviour in negative and positive electrodes.                       
             * "hydrolysis" : str
                 Whether to include hydrolysis in the model. Only implemented for
                 lead-acid models. Can be "false" (default) or "true". If "true", then
@@ -72,6 +80,10 @@ class BatteryModelOptions(pybamm.FuzzyDict):
                 "stress-driven", "reaction-driven", or "stress and reaction-driven".
                 A 2-tuple can be provided for different behaviour in negative and
                 positive electrodes.
+            * "open-circuit potential" : str
+                Sets the model for the open circuit potential. Can be "single" 
+                (default) or "current sigmoid". A 2-tuple can be provided for different 
+                behaviour in negative and positive electrodes.                
             * "operating mode" : str
                 Sets the operating mode for the model. This determines how the current
                 is set. Can be:
@@ -191,6 +203,7 @@ def __init__(self, extra_options):
                 "potential pair",
                 "potential pair quite conductive",
             ],
+            "diffusivity": ["single", "current sigmoid"],
             "dimensionality": [0, 1, 2],
             "electrolyte conductivity": [
                 "default",
@@ -199,6 +212,7 @@ def __init__(self, extra_options):
                 "composite",
                 "integrated",
             ],
+            "exchange-current density": ["single", "current sigmoid"],
             "hydrolysis": ["false", "true"],
             "intercalation kinetics": [
                 "symmetric Butler-Volmer",
@@ -543,6 +557,8 @@ def __init__(self, extra_options):
                         (
                             option
                             in [
+                                "diffusivity",
+                                "exchange-current density",
                                 "intercalation kinetics",
                                 "interface utilisation",
                                 "loss of active material",

pybamm/models/submodels/interface/base_interface.py

@@ -64,6 +64,7 @@ def _get_exchange_current_density(self, variables):
         phase_param = self.phase_param
         domain, Domain = self.domain_Domain
         phase_name = self.phase_name
+        domain_options = getattr(self.options, domain)
 
         c_e = variables[f"{Domain} electrolyte concentration [mol.m-3]"]
         T = variables[f"{Domain} electrode temperature [K]"]
@@ -108,8 +109,23 @@ def _get_exchange_current_density(self, variables):
                     c_s_surf = c_s_surf.orphans[0]
                     c_e = c_e.orphans[0]
                     T = T.orphans[0]
-
-            j0 = phase_param.j0(c_e, c_s_surf, T)
+            # Get main reaction exchange-current density (may have empirical hysteresis)
+            if domain_options["exchange-current density"] == "single":
+                j0 = phase_param.j0(c_e, c_s_surf, T)
+            elif domain_options["exchange-current density"] == "current sigmoid":
+                current = variables["Total current density [A.m-2]"]
+                k = 100
+                if Domain == "Positive":
+                    lithiation_current = current
+                elif Domain == "Negative":
+                    lithiation_current = -current
+                m_lith = pybamm.sigmoid(
+                    0, lithiation_current, k
+                )  # lithiation_current > 0
+                m_delith = 1 - m_lith  # lithiation_current < 0
+                j0_lith = phase_param.j0(c_e, c_s_surf, T, "lithiation")
+                j0_delith = phase_param.j0(c_e, c_s_surf, T, "delithiation")
+                j0 = m_lith * j0_lith + m_delith * j0_delith
 
         elif self.reaction == "lithium metal plating":
             # compute T on the surface of the anode (interface with separator)

pybamm/models/submodels/particle/base_particle.py

@@ -27,14 +27,27 @@ def __init__(self, param, domain, options, phase="primary"):
         domain_options = getattr(self.options, domain)
         self.size_distribution = domain_options["particle size"] == "distribution"
 
-    def _get_effective_diffusivity(self, c, T):
+    def _get_effective_diffusivity(self, c, T, current):
         param = self.param
-        domain = self.domain
+        domain, Domain = self.domain_Domain
         domain_param = self.domain_param
         phase_param = self.phase_param
+        domain_options = getattr(self.options, domain)
 
-        # Get diffusivity
-        D = phase_param.D(c, T)
+        # Get diffusivity (may have empirical hysteresis)
+        if domain_options["diffusivity"] == "single":
+            D = phase_param.D(c, T)
+        elif domain_options["diffusivity"] == "current sigmoid":
+            k = 100
+            if Domain == "Positive":
+                lithiation_current = current
+            elif Domain == "Negative":
+                lithiation_current = -current
+            m_lith = pybamm.sigmoid(0, lithiation_current, k)  # lithiation_current > 0
+            m_delith = 1 - m_lith  # lithiation_current < 0
+            D_lith = phase_param.D(c, T, "lithiation")
+            D_delith = phase_param.D(c, T, "delithiation")
+            D = m_lith * D_lith + m_delith * D_delith
 
         # Account for stress-induced difftusion by defining a multiplicative
         # "stress factor"

pybamm/models/submodels/particle/fickian_diffusion.py

@@ -197,7 +197,8 @@ def get_coupled_variables(self, variables):
                     "current density distribution [A.m-2]"
                 ]
 
-        D_eff = self._get_effective_diffusivity(c_s, T)
+        current = variables["Total current density [A.m-2]"]
+        D_eff = self._get_effective_diffusivity(c_s, T, current)
         N_s = -D_eff * pybamm.grad(c_s)
 
         variables.update(

pybamm/models/submodels/particle/polynomial_profile.py

@@ -194,7 +194,8 @@ def get_coupled_variables(self, variables):
             T = pybamm.PrimaryBroadcast(
                 variables[f"{Domain} electrode temperature [K]"], [f"{domain} particle"]
             )
-            D_eff = self._get_effective_diffusivity(c_s, T)
+            current = variables["Total current density [A.m-2]"]
+            D_eff = self._get_effective_diffusivity(c_s, T, current)
             r = pybamm.SpatialVariable(
                 f"r_{domain[0]}",
                 domain=[f"{domain} particle"],

pybamm/models/submodels/particle/x_averaged_polynomial_profile.py

@@ -104,7 +104,8 @@ def get_coupled_variables(self, variables):
         R = variables[f"X-averaged {domain} particle radius [m]"]
 
         if self.name != "uniform profile":
-            D_eff_av = self._get_effective_diffusivity(c_s_av, T_av)
+            current = variables["Total current density [A.m-2]"]
+            D_eff_av = self._get_effective_diffusivity(c_s_av, T_av, current)
             i_boundary_cc = variables["Current collector current density [A.m-2]"]
             a_av = variables[
                 f"X-averaged {domain} electrode surface area to volume ratio [m-1]"
@@ -183,7 +184,8 @@ def get_coupled_variables(self, variables):
         # Set flux based on polynomial order
         if self.name != "uniform profile":
             T_xav = pybamm.PrimaryBroadcast(T_av, [f"{domain} particle"])
-            D_eff_xav = self._get_effective_diffusivity(c_s_xav, T_xav)
+            current = variables["Total current density [A.m-2]"]
+            D_eff_xav = self._get_effective_diffusivity(c_s_xav, T_xav, current)
             D_eff = pybamm.SecondaryBroadcast(D_eff_xav, [f"{domain} electrode"])
             variables.update(self._get_standard_diffusivity_variables(D_eff))
         if self.name == "uniform profile":

pybamm/parameters/lithium_ion_parameters.py

@@ -524,7 +524,7 @@ def _set_parameters(self):
         if main.options["particle shape"] == "spherical":
             self.a_typ = 3 * pybamm.xyz_average(self.epsilon_s) / self.R_typ
 
-    def D(self, c_s, T):
+    def D(self, c_s, T, lithiation=None):
         """
         Dimensional diffusivity in particle. In the parameter sets this is defined as
         a function of stoichiometry (dimensionless), but in the models we use it as a
@@ -535,16 +535,21 @@ def D(self, c_s, T):
         sto = c_s / self.c_max
         tol = pybamm.settings.tolerances["D__c_s"]
         sto = pybamm.maximum(pybamm.minimum(sto, 1 - tol), tol)
+        if lithiation is None:
+            lithiation = ""
+        else:
+            lithiation = lithiation + " "
         inputs = {
             f"{self.phase_prefactor}{Domain} particle stoichiometry": sto,
             "Temperature [K]": T,
         }
         return pybamm.FunctionParameter(
-            f"{self.phase_prefactor}{Domain} electrode diffusivity [m2.s-1]",
+            f"{self.phase_prefactor}{Domain} electrode {lithiation}"
+            "diffusivity [m2.s-1]",
             inputs,
         )
 
-    def j0(self, c_e, c_s_surf, T):
+    def j0(self, c_e, c_s_surf, T, lithiation=None):
         """Dimensional exchange-current density [A.m-2]"""
         tol = pybamm.settings.tolerances["j0__c_e"]
         c_e = pybamm.maximum(c_e, tol)
@@ -553,6 +558,10 @@ def j0(self, c_e, c_s_surf, T):
             pybamm.minimum(c_s_surf, (1 - tol) * self.c_max), tol * self.c_max
         )
         domain, Domain = self.domain_Domain
+        if lithiation is None:
+            lithiation = ""
+        else:
+            lithiation = lithiation + " "
         inputs = {
             "Electrolyte concentration [mol.m-3]": c_e,
             f"{Domain} particle surface concentration [mol.m-3]": c_s_surf,
@@ -561,7 +570,7 @@ def j0(self, c_e, c_s_surf, T):
             "Temperature [K]": T,
         }
         return pybamm.FunctionParameter(
-            f"{self.phase_prefactor}{Domain} electrode "
+            f"{self.phase_prefactor}{Domain} electrode {lithiation}"
             "exchange-current density [A.m-2]",
             inputs,
         )

tests/unit/test_models/test_full_battery_models/test_base_battery_model.py

@@ -21,8 +21,10 @@
 'contact resistance': 'false' (possible: ['false', 'true'])
 'convection': 'none' (possible: ['none', 'uniform transverse', 'full transverse'])
 'current collector': 'uniform' (possible: ['uniform', 'potential pair', 'potential pair quite conductive'])
+'diffusivity': 'single' (possible: ['single', 'current sigmoid'])
 'dimensionality': 0 (possible: [0, 1, 2])
 'electrolyte conductivity': 'default' (possible: ['default', 'full', 'leading order', 'composite', 'integrated'])
+'exchange-current density': 'single' (possible: ['single', 'current sigmoid'])
 'hydrolysis': 'false' (possible: ['false', 'true'])
 'intercalation kinetics': 'symmetric Butler-Volmer' (possible: ['symmetric Butler-Volmer', 'asymmetric Butler-Volmer', 'linear', 'Marcus', 'Marcus-Hush-Chidsey'])
 'interface utilisation': 'full' (possible: ['full', 'constant', 'current-driven'])

tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py

@@ -359,3 +359,11 @@ def test_well_posed_particle_phases_sei(self):
     def test_well_posed_current_sigmoid_ocp(self):
         options = {"open-circuit potential": "current sigmoid"}
         self.check_well_posedness(options)
+
+    def test_well_posed_current_sigmoid_exchange_current(self):
+        options = {"exchange-current density": "current sigmoid"}
+        self.check_well_posedness(options)
+
+    def test_well_posed_current_sigmoid_diffusivity(self):
+        options = {"diffusivity": "current sigmoid"}
+        self.check_well_posedness(options)