Plett OCP Hysteresis Model

CHANGELOG.md

@@ -14,6 +14,7 @@
 - Renamed "electrode diffusivity" to "particle diffusivity" as a non-breaking change with a deprecation warning ([#3624](https://github.com/pybamm-team/PyBaMM/pull/3624))
 - Add support for BPX version 0.4.0 which allows for blended electrodes and user-defined parameters in BPX([#3414](https://github.com/pybamm-team/PyBaMM/pull/3414))
 - Added `by_submodel` feature in `print_parameter_info` method to allow users to print parameters and types of submodels in a tabular and readable format ([#3628](https://github.com/pybamm-team/PyBaMM/pull/3628))
+- Added `WyciskOpenCircuitPotential` for differential capacity hysteresis state open-circuit potential submodel ([#3593](https://github.com/pybamm-team/PyBaMM/pull/3593))
 
 ## Bug Fixes
 

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

@@ -7,3 +7,4 @@ Open-circuit potential models
   current_sigmoid_ocp
   single_ocp
   msmr_ocp
+  wycisk_ocp

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

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

docs/source/examples/index.rst

@@ -50,6 +50,7 @@ The notebooks are organised into subfolders, and can be viewed in the galleries
     notebooks/models/compare-particle-diffusion-models.ipynb
     notebooks/models/composite_particle.ipynb
     notebooks/models/coupled-degradation.ipynb
+    notebooks/models/differential-capacity-hysteresis-state.ipynb
     notebooks/models/DFN-with-particle-size-distributions.ipynb
     notebooks/models/DFN.ipynb
     notebooks/models/electrode-state-of-health.ipynb

pybamm/CITATIONS.bib

@@ -691,3 +691,16 @@ @article{landesfeind2019temperature
   year={2019},
   publisher={The Electrochemical Society}
 }
+
+@article{Wycisk2022,
+title = {Modified Plett-model for modeling voltage hysteresis in lithium-ion cells},
+journal = {Journal of Energy Storage},
+volume = {52},
+pages = {105016},
+year = {2022},
+issn = {2352-152X},
+doi = {https://doi.org/10.1016/j.est.2022.105016},
+url = {https://www.sciencedirect.com/science/article/pii/S2352152X22010192},
+author = {Dominik Wycisk and Marc Oldenburger and Marc Gerry Stoye and Toni Mrkonjic and Arnulf Latz},
+keywords = {Lithium-ion battery, Voltage hysteresis, Plett-model, Silicon–graphite anode},
+}

pybamm/models/full_battery_models/base_battery_model.py

@@ -102,7 +102,7 @@ class BatteryModelOptions(pybamm.FuzzyDict):
                 reactions.
             * "open-circuit potential" : str
                 Sets the model for the open circuit potential. Can be "single"
-                (default), "current sigmoid", or "MSMR". If "MSMR" then the "particle"
+                (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.
             * "operating mode" : str
@@ -263,7 +263,7 @@ def __init__(self, extra_options):
                 "stress and reaction-driven",
             ],
             "number of MSMR reactions": ["none"],
-            "open-circuit potential": ["single", "current sigmoid", "MSMR"],
+            "open-circuit potential": ["single", "current sigmoid", "MSMR", "Wycisk"],
             "operating mode": [
                 "current",
                 "voltage",

pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py

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

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

@@ -2,3 +2,4 @@
 from .single_ocp import SingleOpenCircuitPotential
 from .current_sigmoid_ocp import CurrentSigmoidOpenCircuitPotential
 from .msmr_ocp import MSMROpenCircuitPotential
+from .wycisk_ocp import WyciskOpenCircuitPotential

pybamm/models/submodels/interface/open_circuit_potential/wycisk_ocp.py

@@ -0,0 +1,160 @@
+#
+# from Wycisk 2022
+#
+import pybamm
+from . import BaseOpenCircuitPotential
+
+
+class WyciskOpenCircuitPotential(BaseOpenCircuitPotential):
+    """
+    Class for open-circuit potential with hysteresis based on the approach outlined by Wycisk :footcite:t:'Wycisk2022'.
+    This approach employs a differential capacity hysteresis state variable. The decay and switching of the hysteresis state
+    is tunable via two additional parameters. The hysteresis state is updated based on the current and the differential capacity.
+    """
+
+    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
+        phase = self.phase
+
+        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
+            sto_bulk = variables[f"{Domain} electrode {phase_name}stoichiometry"]
+            c_scale = self.phase_param.c_max
+            variables[f"Total lithium in {phase} phase in {domain} electrode [mol]"] = (
+                sto_bulk * c_scale
+            )  # c_s_vol * L * A
+
+            ocp_surf_eq = self.phase_param.U(sto_surf, T)
+            variables[f"{Domain} electrode {phase_name}equilibrium OCP [V]"] = (
+                ocp_surf_eq
+            )
+
+            T_bulk = pybamm.xyz_average(pybamm.size_average(T))
+            ocp_bulk_eq = self.phase_param.U(sto_bulk, T_bulk)
+            variables[f"{Domain} electrode {phase_name}bulk equilibrium OCP [V]"] = (
+                ocp_bulk_eq
+            )
+
+            H = self.phase_param.hysteresis(sto_surf)
+            variables[f"{Domain} electrode {phase_name}OCP hysteresis [V]"] = H
+
+            dU = self.phase_param.U(sto_surf, T_bulk).diff(sto_surf)
+            dQ = self.phase_param.Q(sto_surf).diff(sto_surf)
+            dQdU = dQ / dU
+            variables[
+                f"{Domain} electrode {phase_name}differential capacity [A.s.V-1]"
+            ] = dQdU
+
+            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 = ocp_surf_eq + 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 = ocp_surf_eq + H * h
+            # must not be a psd
+            else:
+                ocp_surf = ocp_surf_eq + H * h
+
+            H_s_av = pybamm.size_average(H_x_av)
+            h_s_av = pybamm.size_average(h_x_av)
+
+            ocp_bulk = ocp_bulk_eq + H_s_av * h_s_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
+
+        current = variables[
+            f"{Domain} electrode {phase_name}interfacial current density [A.m-2]"
+        ]
+        # check if composite or not
+        if phase_name != "":
+            Q_cell = variables[f"{Domain} electrode {phase_name}phase capacity [A.h]"]
+        else:
+            Q_cell = variables[f"{Domain} electrode capacity [A.h]"]
+
+        dQdU = variables[
+            f"{Domain} electrode {phase_name}differential capacity [A.s.V-1]"
+        ]
+        dQdU = dQdU.orphans[0]
+        K = self.phase_param.hysteresis_decay
+        K_x = self.phase_param.hysteresis_switch
+        h = variables[f"{Domain} electrode {phase_name}hysteresis state"]
+
+        dhdt = (
+            K * (current / (Q_cell * (dQdU**K_x))) * (1 - pybamm.sign(current) * h)
+        )  #! current is backwards for a halfcell
+        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] = pybamm.Scalar(0)

pybamm/parameters/lithium_ion_parameters.py

@@ -512,6 +512,14 @@ def _set_parameters(self):
         eps_c_init_av = pybamm.xyz_average(
             self.epsilon_s * pybamm.r_average(self.c_init)
         )
+        # if self.options['open-circuit potential'] == 'Plett':
+        self.hysteresis_decay = pybamm.Parameter(
+            f"{pref}{Domain} particle hysteresis decay rate"
+        )
+        self.hysteresis_switch = pybamm.Parameter(
+            f"{pref}{Domain} particle hysteresis switching factor"
+        )
+        self.h_init = pybamm.Scalar(0)
 
         if self.options["open-circuit potential"] != "MSMR":
             self.U_init = self.U(self.sto_init_av, main.T_init)
@@ -641,6 +649,33 @@ def U(self, sto, T, lithiation=None):
             out.print_name = r"U_\mathrm{p}(c^\mathrm{surf}_\mathrm{s,p}, T)"
         return out
 
+    def hysteresis(self, sto):
+        """Dimensional hysteresis [V]"""
+        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}
+        lith_ref = pybamm.FunctionParameter(
+            f"{self.phase_prefactor}{Domain} electrode lithiation OCP [V]", inputs
+        )
+        delith_ref = pybamm.FunctionParameter(
+            f"{self.phase_prefactor}{Domain} electrode delithiation OCP [V]", inputs
+        )
+        h_ref = abs(delith_ref - lith_ref) / 2
+        # h_ref = pybamm.FunctionParameter(
+        #     f"{self.phase_prefactor}{Domain} electrode OCP hysteresis [V]", inputs
+        # )
+        return h_ref
+
+    def Q(self, sto):
+        """Capacity change as a function of stoichiometry"""
+        c_max = self.c_max
+        epsilon_s_av = self.epsilon_s_av
+        V_electrode = self.main_param.A_cc * self.domain_param.L
+        Li_max = c_max * V_electrode * epsilon_s_av
+        Q_max = Li_max * self.main_param.F / 3600
+        return Q_max * sto
+
     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

@@ -68,6 +68,28 @@ def test_current_sigmoid_ocp(self):
         modeltest = tests.StandardModelTest(model, parameter_values=parameter_values)
         modeltest.test_all(skip_output_tests=True)
 
+    def test_wycisk_ocp(self):
+        options = {"open-circuit potential": ("Wycisk", "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": 1,
+                "Negative particle hysteresis switching factor": 1,
+                # "Negative electrode OCP hysteresis [V]": lambda sto: 1,
+            },
+            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

@@ -33,7 +33,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'])
+'open-circuit potential': 'single' (possible: ['single', 'current sigmoid', 'MSMR', 'Wycisk'])
 'operating mode': 'current' (possible: ['current', 'voltage', 'power', 'differential power', 'explicit power', 'resistance', 'differential resistance', 'explicit resistance', 'CCCV'])
 'particle': 'Fickian diffusion' (possible: ['Fickian diffusion', 'fast 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

@@ -406,6 +406,10 @@ def test_well_posed_current_sigmoid_ocp(self):
         options = {"open-circuit potential": "current sigmoid"}
         self.check_well_posedness(options)
 
+    def test_well_posed_wycisk_ocp(self):
+        options = {"open-circuit potential": "Wycisk"}
+        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

@@ -35,6 +35,20 @@ def test_well_posed_current_sigmoid_ocp_with_psd(self):
         }
         self.check_well_posedness(options)
 
+    def test_well_posed_wycisk_ocp_with_psd(self):
+        options = {
+            "open-circuit potential": "Wycisk",
+            "particle size": "distribution",
+        }
+        self.check_well_posedness(options)
+
+    def test_well_posed_wycisk_ocp_with_composite(self):
+        options = {
+            "open-circuit potential": (("Wycisk", "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

@@ -116,6 +116,13 @@ def test_msmr(self):
         model = pybamm.lithium_ion.MPM(options)
         model.check_well_posedness()
 
+    def test_wycisk_ocp(self):
+        options = {
+            "open-circuit potential": "Wycisk",
+        }
+        model = pybamm.lithium_ion.MPM(options)
+        model.check_well_posedness()
+
 
 class TestMPMExternalCircuits(TestCase):
     def test_well_posed_voltage(self):