Repository: pybamm-team/PyBaMM

Issue:
Empirical hysteresis 
We currently support an empirical “current sigmoid” model for OCV hysteresis, where users can give a lithiation and delithiation OCV and switch between them depending on the sign of the current. Including such a model has knock-on effects for other model parameters, which may also need to be depend on the sign of the current in order to accurately fit experimental data. In particular, we should add empirical hysteresis for exchange-current and diffusivity. 

Longer term we plan to add more sophisticated and self-consistent hysteresis models, but this is a good pragmatic solution for now. 

Discussion:

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 084a6f7db..17170ab05 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -2,9 +2,10 @@
 
 ## Breaking changes
 
-- PyBaMM now has optional dependencies that can be installed with the pattern `pip install pybamm[option]` e.g. `pybamm[plot]` ([#3044](https://github.com/pybamm-team/PyBaMM/pull/3044))
-- `pybamm_install_jax` is deprecated. It is now replaced with `pip install pybamm[jax]` ([#3163](https://github.com/pybamm-team/PyBaMM/pull/3163))
+- 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))
+- `pybamm_install_jax` is deprecated. It is now replaced with `pip install pybamm[jax]` ([#3163](https://github.com/pybamm-team/PyBaMM/pull/3163))
+- PyBaMM now has optional dependencies that can be installed with the pattern `pip install pybamm[option]` e.g. `pybamm[plot]` ([#3044](https://github.com/pybamm-team/PyBaMM/pull/3044))
 
 ## Bug fixes
 
diff --git a/examples/scripts/emperical_hysteresis.py b/examples/scripts/emperical_hysteresis.py
new file mode 100644
index 000000000..3f3fc7c64
--- /dev/null
+++ b/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]",
+    ]
+)
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index 4ea461113..7a3987410 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/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 @@ class BatteryModelOptions(pybamm.FuzzyDict):
                 "potential pair",
                 "potential pair quite conductive",
             ],
+            "diffusivity": ["single", "current sigmoid"],
             "dimensionality": [0, 1, 2],
             "electrolyte conductivity": [
                 "default",
@@ -199,6 +212,7 @@ class BatteryModelOptions(pybamm.FuzzyDict):
                 "composite",
                 "integrated",
             ],
+            "exchange-current density": ["single", "current sigmoid"],
             "hydrolysis": ["false", "true"],
             "intercalation kinetics": [
                 "symmetric Butler-Volmer",
@@ -543,6 +557,8 @@ class BatteryModelOptions(pybamm.FuzzyDict):
                         (
                             option
                             in [
+                                "diffusivity",
+                                "exchange-current density",
                                 "intercalation kinetics",
                                 "interface utilisation",
                                 "loss of active material",
diff --git a/pybamm/models/submodels/interface/base_interface.py b/pybamm/models/submodels/interface/base_interface.py
index bc22a639b..080c8f54a 100644
--- a/pybamm/models/submodels/interface/base_interface.py
+++ b/pybamm/models/submodels/interface/base_interface.py
@@ -64,6 +64,7 @@ class BaseInterface(pybamm.BaseSubModel):
         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 @@ class BaseInterface(pybamm.BaseSubModel):
                     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)
diff --git a/pybamm/models/submodels/particle/base_particle.py b/pybamm/models/submodels/particle/base_particle.py
index 3ef6ee4d7..94e9fb5c4 100644
--- a/pybamm/models/submodels/particle/base_particle.py
+++ b/pybamm/models/submodels/particle/base_particle.py
@@ -27,14 +27,27 @@ class BaseParticle(pybamm.BaseSubModel):
         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"
diff --git a/pybamm/models/submodels/particle/fickian_diffusion.py b/pybamm/models/submodels/particle/fickian_diffusion.py
index c85716373..7707284a3 100644
--- a/pybamm/models/submodels/particle/fickian_diffusion.py
+++ b/pybamm/models/submodels/particle/fickian_diffusion.py
@@ -197,7 +197,8 @@ class FickianDiffusion(BaseParticle):
                     "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(
diff --git a/pybamm/models/submodels/particle/polynomial_profile.py b/pybamm/models/submodels/particle/polynomial_profile.py
index 62fe38a8b..dda3a3066 100644
--- a/pybamm/models/submodels/particle/polynomial_profile.py
+++ b/pybamm/models/submodels/particle/polynomial_profile.py
@@ -194,7 +194,8 @@ class PolynomialProfile(BaseParticle):
             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"],
diff --git a/pybamm/models/submodels/particle/x_averaged_polynomial_profile.py b/pybamm/models/submodels/particle/x_averaged_polynomial_profile.py
index cef3c6ec0..45497c4c5 100644
--- a/pybamm/models/submodels/particle/x_averaged_polynomial_profile.py
+++ b/pybamm/models/submodels/particle/x_averaged_polynomial_profile.py
@@ -104,7 +104,8 @@ class XAveragedPolynomialProfile(PolynomialProfile):
         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 @@ class XAveragedPolynomialProfile(PolynomialProfile):
         # 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":
diff --git a/pybamm/parameters/lithium_ion_parameters.py b/pybamm/parameters/lithium_ion_parameters.py
index 6d23bc187..d489ac668 100644
--- a/pybamm/parameters/lithium_ion_parameters.py
+++ b/pybamm/parameters/lithium_ion_parameters.py
@@ -524,7 +524,7 @@ class ParticleLithiumIonParameters(BaseParameters):
         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 @@ class ParticleLithiumIonParameters(BaseParameters):
         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 @@ class ParticleLithiumIonParameters(BaseParameters):
             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 @@ class ParticleLithiumIonParameters(BaseParameters):
             "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,
         )


Test patch:
diff --git a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
index ee51f864c..b0dad88cc 100644
--- a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
+++ b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py
@@ -21,8 +21,10 @@ PRINT_OPTIONS_OUTPUT = """\
 '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'])
diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
index 99b2b1014..1911f2989 100644
--- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
+++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py
@@ -359,3 +359,11 @@ class BaseUnitTestLithiumIon:
     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)