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

…Issue pybamm-team#4808)

Test functions have been added.
Minor modifications have been added to the module axen_ocp to pass all the tests.

docs/source/api/models/submodels/interface/open_circuit_potential/axen_ocp.rst

@@ -0,0 +1,5 @@
+Axen Open Circuit Potential
+=============================
+
+.. autoclass:: pybamm.open_circuit_potential.AxenOpenCircuitPotential
+    :members:

docs/source/api/models/submodels/interface/open_circuit_potential/index.rst

@@ -8,3 +8,4 @@ Open-circuit potential models
   single_ocp
   msmr_ocp
   wycisk_ocp
+  axen_ocp

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

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

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

@@ -12,9 +12,6 @@ class AxenOpenCircuitPotential(BaseOpenCircuitPotential):
     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
@@ -68,12 +65,8 @@ def get_coupled_variables(self, variables):
             ] = 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"
-            )
+            lith_ref = self.phase_param.U(sto_surf, T, "lithiation")
+            delith_ref = self.phase_param.U(sto_surf, T, "delithiation")
             H = lith_ref - delith_ref
             variables[f"{Domain} electrode {phase_name}OCP hysteresis [V]"] = H
 
@@ -102,7 +95,11 @@ def get_coupled_variables(self, variables):
             else:
                 ocp_surf = delith_ref + H * h
 
-            ocp_bulk = delith_ref_x_av + H_x_av * h_x_av
+            delith_ref_s_av = pybamm.size_average(delith_ref_x_av)
+            H_s_av = pybamm.size_average(H_x_av)
+            h_s_av = pybamm.size_average(h_x_av)
+
+            ocp_bulk = delith_ref_s_av + H_s_av * h_s_av
 
             dUdT = self.phase_param.dUdT(sto_surf)
 
@@ -127,10 +124,12 @@ def set_rhs(self, variables):
         S_delith = -i_vol * K_delith / (self.param.F * c_max * epsl)
 
         i_vol_sign = pybamm.sign(i_vol)
-        f = 0.5 * (1 - i_vol_sign) + i_vol_sign * h
-        d = 0.5 * (1 - i_vol_sign) * S_lith + 0.5 * (1 + i_vol_sign) * S_delith
+        signed_h = 0.5 * (1 - i_vol_sign) + i_vol_sign * h
+        rate_coefficient = (
+            0.5 * (1 - i_vol_sign) * S_lith + 0.5 * (1 + i_vol_sign) * S_delith
+        )
 
-        dhdt = f * d
+        dhdt = rate_coefficient * signed_h
 
         self.rhs[h] = dhdt
 

src/pybamm/parameters/lithium_ion_parameters.py

@@ -648,49 +648,6 @@ 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].

tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py

@@ -90,6 +90,30 @@ def test_wycisk_ocp(self):
         modeltest = tests.StandardModelTest(model, parameter_values=parameter_values)
         modeltest.test_all(skip_output_tests=True)
 
+    def test_axen_ocp(self):
+        options = {"open-circuit potential": ("Axen", "single")}
+        model = pybamm.lithium_ion.MPM(options)
+        parameter_values = pybamm.ParameterValues("Chen2020")
+        parameter_values = pybamm.get_size_distribution_parameters(parameter_values)
+        parameter_values.update(
+            {
+                "Negative electrode lithiation OCP [V]": lambda sto: parameter_values[
+                    "Negative electrode OCP [V]"
+                ](sto)
+                - 0.1,
+                "Negative electrode delithiation OCP [V]": lambda sto: parameter_values[
+                    "Negative electrode OCP [V]"
+                ](sto)
+                + 0.1,
+                "Negative particle hysteresis decay rate for lithiation": 15,
+                "Negative particle hysteresis decay rate for delithiation": 15,
+                "Initial hysteresis state in negative electrode": 0.5,
+            },
+            check_already_exists=False,
+        )
+        modeltest = tests.StandardModelTest(model, parameter_values=parameter_values)
+        modeltest.test_all(skip_output_tests=True)
+
     def test_voltage_control(self):
         options = {"operating mode": "voltage"}
         model = pybamm.lithium_ion.MPM(options)

tests/unit/test_models/test_full_battery_models/test_base_battery_model.py

@@ -34,7 +34,7 @@
 'lithium plating porosity change': 'false' (possible: ['false', 'true'])
 'loss of active material': 'stress-driven' (possible: ['none', 'stress-driven', 'reaction-driven', 'current-driven', 'stress and reaction-driven'])
 'number of MSMR reactions': 'none' (possible: ['none'])
-'open-circuit potential': 'single' (possible: ['single', 'current sigmoid', 'MSMR', 'Wycisk'])
+'open-circuit potential': 'single' (possible: ['single', 'current sigmoid', 'MSMR', 'Wycisk', 'Axen'])
 'operating mode': 'current' (possible: ['current', 'voltage', 'power', 'differential power', 'explicit power', 'resistance', 'differential resistance', 'explicit resistance', 'CCCV'])
 'particle': 'Fickian diffusion' (possible: ['Fickian diffusion', 'uniform profile', 'quadratic profile', 'quartic profile', 'MSMR'])
 'particle mechanics': 'swelling only' (possible: ['none', 'swelling only', 'swelling and cracking'])

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

@@ -503,6 +503,10 @@ def test_well_posed_wycisk_ocp(self):
         options = {"open-circuit potential": "Wycisk"}
         self.check_well_posedness(options)
 
+    def test_well_posed_axen_ocp(self):
+        options = {"open-circuit potential": "Axen"}
+        self.check_well_posedness(options)
+
     def test_well_posed_msmr(self):
         options = {
             "open-circuit potential": "MSMR",

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

@@ -52,6 +52,20 @@ def test_well_posed_wycisk_ocp_with_composite(self):
         }
         self.check_well_posedness(options)
 
+    def test_well_posed_axen_ocp_with_psd(self):
+        options = {
+            "open-circuit potential": "Axen",
+            "particle size": "distribution",
+        }
+        self.check_well_posedness(options)
+
+    def test_well_posed_axen_ocp_with_composite(self):
+        options = {
+            "open-circuit potential": (("Axen", "single"), "single"),
+            "particle phases": ("2", "1"),
+        }
+        self.check_well_posedness(options)
+
     def test_well_posed_external_circuit_explicit_power(self):
         options = {"operating mode": "explicit power"}
         self.check_well_posedness(options)

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

@@ -123,6 +123,13 @@ def test_wycisk_ocp(self):
         model = pybamm.lithium_ion.MPM(options)
         model.check_well_posedness()
 
+    def test_axen_ocp(self):
+        options = {
+            "open-circuit potential": "Axen",
+        }
+        model = pybamm.lithium_ion.MPM(options)
+        model.check_well_posedness()
+
     def test_mpm_with_lithium_plating(self):
         options = {
             "lithium plating": "irreversible",