Differential method from Axen et al. (2022) for OCP with hysteresis (…

…Issue pybamm-team#4808)

src/pybamm/CITATIONS.bib

@@ -47,6 +47,19 @@ @article{Andersson2019
   doi = {10.1007/s12532-018-0139-4},
 }
 
+@article{Axen2022,
+title = {Evaluation of hysteresis expressions in a lumped voltage prediction model of a NiMH battery system in stationary storage applications},
+journal = {Journal of Energy Storage},
+volume = {48},
+pages = {103985},
+year = {2022},
+issn = {2352-152X},
+doi = {https://doi.org/10.1016/j.est.2022.103985},
+url = {https://www.sciencedirect.com/science/article/pii/S2352152X22000329},
+author = {Jenny B\"{o}rjesson Ax\'{e}n and Henrik Ekstr\"{o}m and Erika Widenkvist Zetterstr\"{o}m and G\"{o}ran Lindbergh},
+keywords = {Battery modelling, OCV hysteresis, NiMH, Voltage Response Model, Physical model, Application verification},
+}
+
 @article{Baker2018,
   title={Multi-species, multi-reaction model for porous intercalation electrodes: Part I. Model formulation and a perturbation solution for low-scan-rate, linear-sweep voltammetry of a spinel lithium manganese oxide electrode},
   author={Baker, Daniel R and Verbrugge, Mark W},

src/pybamm/models/full_battery_models/base_battery_model.py

@@ -102,9 +102,10 @@ class BatteryModelOptions(pybamm.FuzzyDict):
                 reactions.
             * "open-circuit potential" : str
                 Sets the model for the open circuit potential. Can be "single"
-                (default), "current sigmoid", "Wycisk", or "MSMR". If "MSMR" then the "particle"
-                option must also be "MSMR". A 2-tuple can be provided for different
-                behaviour in negative and positive electrodes.
+                (default), "current sigmoid", "Wycisk", "Axen", or "MSMR".
+                If "MSMR" then the "particle" option must also be "MSMR".
+                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:
@@ -273,7 +274,13 @@ def __init__(self, extra_options):
                 "stress and reaction-driven",
             ],
             "number of MSMR reactions": ["none"],
-            "open-circuit potential": ["single", "current sigmoid", "MSMR", "Wycisk"],
+            "open-circuit potential": [
+                "single",
+                "current sigmoid",
+                "MSMR",
+                "Wycisk",
+                "Axen",
+            ],
             "operating mode": [
                 "current",
                 "voltage",

src/pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py

@@ -252,6 +252,9 @@ def set_open_circuit_potential_submodel(self):
                 elif ocp_option == "Wycisk":
                     pybamm.citations.register("Wycisk2022")
                     ocp_model = ocp_submodels.WyciskOpenCircuitPotential
+                elif ocp_option == "Axen":
+                    pybamm.citations.register("Axen2022")
+                    ocp_model = ocp_submodels.AxenOpenCircuitPotential
                 elif ocp_option == "MSMR":
                     ocp_model = ocp_submodels.MSMROpenCircuitPotential
                 self.submodels[f"{domain} {phase} open-circuit potential"] = ocp_model(

src/pybamm/models/submodels/interface/open_circuit_potential/__init__.py

@@ -3,5 +3,6 @@
 from .current_sigmoid_ocp import CurrentSigmoidOpenCircuitPotential
 from .msmr_ocp import MSMROpenCircuitPotential
 from .wycisk_ocp import WyciskOpenCircuitPotential
+from .axen_ocp import AxenOpenCircuitPotential
 
-__all__ = ['base_ocp', 'current_sigmoid_ocp', 'msmr_ocp', 'single_ocp', 'wycisk_ocp']
+__all__ = ['axen_ocp', 'base_ocp', 'current_sigmoid_ocp', 'msmr_ocp', 'single_ocp', 'wycisk_ocp']

src/pybamm/models/submodels/interface/open_circuit_potential/axen_ocp.py

@@ -0,0 +1,141 @@
+#
+# from Axen 2022
+#
+import pybamm
+from . import BaseOpenCircuitPotential
+
+
+class AxenOpenCircuitPotential(BaseOpenCircuitPotential):
+    """
+    Class for open-circuit potential with hysteresis based on the approach outlined by Axen et al. https://doi.org/10.1016/j.est.2022.103985.
+    This approach tracks the evolution of the hysteresis using an ODE system scaled by the volumetric interfacial current density and
+    a decay rate parameter.
+    """
+
+    def __init__(self, param, domain, reaction, options, phase="primary"):
+        super().__init__(param, domain, reaction, options=options, phase=phase)
+
+    def get_fundamental_variables(self):
+        domain, Domain = self.domain_Domain
+        phase_name = self.phase_name
+        h = pybamm.Variable(
+            f"{Domain} electrode {phase_name}hysteresis state",
+            domains={
+                "primary": f"{domain} electrode",
+                "secondary": "current collector",
+            },
+        )
+        return {
+            f"{Domain} electrode {phase_name}hysteresis state": h,
+        }
+
+    def get_coupled_variables(self, variables):
+        domain, Domain = self.domain_Domain
+        phase_name = self.phase_name
+
+        if self.reaction == "lithium-ion main":
+            T = variables[f"{Domain} electrode temperature [K]"]
+            h = variables[f"{Domain} electrode {phase_name}hysteresis state"]
+
+            # For "particle-size distribution" models, take distribution version
+            # of c_s_surf that depends on particle size.
+            domain_options = getattr(self.options, domain)
+            if domain_options["particle size"] == "distribution":
+                sto_surf = variables[
+                    f"{Domain} {phase_name}particle surface stoichiometry distribution"
+                ]
+                # If variable was broadcast, take only the orphan
+                if isinstance(sto_surf, pybamm.Broadcast) and isinstance(
+                    T, pybamm.Broadcast
+                ):
+                    sto_surf = sto_surf.orphans[0]
+                    T = T.orphans[0]
+                T = pybamm.PrimaryBroadcast(T, [f"{domain} {phase_name}particle size"])
+                h = pybamm.PrimaryBroadcast(h, [f"{domain} {phase_name}particle size"])
+            else:
+                sto_surf = variables[
+                    f"{Domain} {phase_name}particle surface stoichiometry"
+                ]
+                # If variable was broadcast, take only the orphan
+                if isinstance(sto_surf, pybamm.Broadcast) and isinstance(
+                    T, pybamm.Broadcast
+                ):
+                    sto_surf = sto_surf.orphans[0]
+                    T = T.orphans[0]
+
+            variables[
+                f"{Domain} electrode {phase_name}hysteresis state distribution"
+            ] = h
+
+            # Bulk OCP is from the average SOC and temperature
+            lith_ref = self.phase_param.U_hysteresis_branch(
+                sto_surf, T, branch="lithiation"
+            )
+            delith_ref = self.phase_param.U_hysteresis_branch(
+                sto_surf, T, branch="delithiation"
+            )
+            H = lith_ref - delith_ref
+            variables[f"{Domain} electrode {phase_name}OCP hysteresis [V]"] = H
+
+            delith_ref_x_av = pybamm.x_average(delith_ref)
+            H_x_av = pybamm.x_average(H)
+            h_x_av = pybamm.x_average(h)
+            variables[f"X-averaged {domain} electrode {phase_name}hysteresis state"] = (
+                h_x_av
+            )
+
+            # check if psd
+            if domain_options["particle size"] == "distribution":
+                # should always be true
+                if f"{domain} particle size" in sto_surf.domains["primary"]:
+                    # check if MPM Model
+                    if "current collector" in sto_surf.domains["secondary"]:
+                        ocp_surf = delith_ref_x_av + H_x_av * h_x_av
+                    # must be DFN with PSD model
+                    elif (
+                        f"{domain} electrode" in sto_surf.domains["secondary"]
+                        or f"{domain} {phase_name}particle size"
+                        in sto_surf.domains["primary"]
+                    ):
+                        ocp_surf = delith_ref + H * h
+            # must not be a psd
+            else:
+                ocp_surf = delith_ref + H * h
+
+            ocp_bulk = delith_ref_x_av + H_x_av * h_x_av
+
+            dUdT = self.phase_param.dUdT(sto_surf)
+
+        variables.update(self._get_standard_ocp_variables(ocp_surf, ocp_bulk, dUdT))
+        return variables
+
+    def set_rhs(self, variables):
+        domain, Domain = self.domain_Domain
+        phase_name = self.phase_name
+
+        i_vol = variables[
+            f"{Domain} electrode {phase_name}volumetric interfacial current density [A.m-3]"
+        ]
+        c_max = self.phase_param.c_max
+        epsl = self.phase_param.epsilon_s
+
+        K_lith = self.phase_param.hysteresis_decay_lithiation
+        K_delith = self.phase_param.hysteresis_decay_delithiation
+        h = variables[f"{Domain} electrode {phase_name}hysteresis state"]
+
+        i_vol_sign = pybamm.sign(i_vol)
+        S_lith = -i_vol * K_lith / (self.param.F * c_max * epsl)
+        S_delith = -i_vol * K_delith / (self.param.F * c_max * epsl)
+        f = 0.5 * (1 - i_vol_sign) + i_vol_sign * h
+        d = 0.5 * (1 - i_vol_sign) / (S_lith + 1e-10) + 0.5 * (1 + i_vol_sign) / (
+            S_delith + 1e-10
+        )
+        dhdt = f / d
+
+        self.rhs[h] = dhdt
+
+    def set_initial_conditions(self, variables):
+        domain, Domain = self.domain_Domain
+        phase_name = self.phase_name
+        h = variables[f"{Domain} electrode {phase_name}hysteresis state"]
+        self.initial_conditions[h] = self.phase_param.h_init

src/pybamm/parameters/lithium_ion_parameters.py

@@ -457,6 +457,12 @@ def _set_parameters(self):
         self.hysteresis_decay = pybamm.Parameter(
             f"{pref}{Domain} particle hysteresis decay rate"
         )
+        self.hysteresis_decay_lithiation = pybamm.Parameter(
+            f"{pref}{Domain} particle hysteresis decay rate for lithiation"
+        )
+        self.hysteresis_decay_delithiation = pybamm.Parameter(
+            f"{pref}{Domain} particle hysteresis decay rate for delithiation"
+        )
         self.hysteresis_switch = pybamm.Parameter(
             f"{pref}{Domain} particle hysteresis switching factor"
         )
@@ -642,6 +648,49 @@ def U(self, sto, T, lithiation=None):
             out.print_name = r"U_\mathrm{p}(c^\mathrm{surf}_\mathrm{s,p}, T)"
         return out
 
+    def U_hysteresis_branch(self, sto, T, branch="lithiation"):
+        """
+        Dimensional open-circuit potential [V] for lithiation or delithiation,
+        calculated as U(x,T) = U_ref(x) + dUdT(x) * (T - T_ref).
+        See the documentation for dUdT for more details.
+        """
+        # bound stoichiometry between tol and 1-tol. Adding 1/sto + 1/(sto-1) later
+        # will ensure that ocp goes to +- infinity if sto goes into that region
+        # anyway
+        Domain = self.domain.capitalize()
+        tol = pybamm.settings.tolerances["U__c_s"]
+        sto = pybamm.maximum(pybamm.minimum(sto, 1 - tol), tol)
+        inputs = {f"{self.phase_prefactor}{Domain} particle stoichiometry": sto}
+
+        # if "branch" takes a value different from "lithiation" or "delithiation",
+        # the function "u_ref" defaults to a hysteresis-free ocp
+        if branch == "lithiation":
+            u_ref = pybamm.FunctionParameter(
+                f"{self.phase_prefactor}{Domain} electrode lithiation OCP [V]", inputs
+            )
+        elif branch == "delithiation":
+            u_ref = pybamm.FunctionParameter(
+                f"{self.phase_prefactor}{Domain} electrode delithiation OCP [V]", inputs
+            )
+        else:
+            u_ref = pybamm.FunctionParameter(
+                f"{self.phase_prefactor}{Domain} electrode OCP [V]", inputs
+            )
+
+        dudt_func = self.dUdT(sto)
+        u_ref = u_ref + (T - self.main_param.T_ref) * dudt_func
+
+        # add a term to ensure that the OCP goes to infinity at 0 and -infinity at 1
+        # this will not affect the OCP for most values of sto
+        # see #1435
+        out = u_ref + 1e-6 * (1 / sto + 1 / (sto - 1))
+
+        if self.domain == "negative":
+            out.print_name = r"U_\mathrm{n}(c^\mathrm{surf}_\mathrm{s,n}, T)"
+        elif self.domain == "positive":
+            out.print_name = r"U_\mathrm{p}(c^\mathrm{surf}_\mathrm{s,p}, T)"
+        return out
+
     def dUdT(self, sto):
         """
         Dimensional entropic change of the open-circuit potential [V.K-1].