Plett OCP Hysteresis Model

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", "Plett", 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", "Plett"],
             "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 == "Plett":
+                    pybamm.citations.register("Wycisk2022")
+                    ocp_model = ocp_submodels.PlettOpenCircuitPotential
                 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 .plett_ocp import PlettOpenCircuitPotential

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

@@ -0,0 +1,150 @@
+#
+# from Wycisk 2022
+#
+import pybamm
+from . import BaseOpenCircuitPotential
+
+
+class PlettOpenCircuitPotential(BaseOpenCircuitPotential):
+    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
+
+        # get dQ, change in capacity with respect to stoichiometry
+        c_max = self.phase_param.c_max
+        epsilon_s_av = self.phase_param.epsilon_s_av
+        V_electrode = self.param.A_cc * self.domain_param.L
+        Li_max = c_max * V_electrode * epsilon_s_av
+        Q_max = Li_max * self.param.F / 3600
+        dQ = Q_max
+
+        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} particle size"])
+                h = pybamm.PrimaryBroadcast(h, [f"{domain} 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]
+
+            ocp_surf_eq = self.phase_param.U(sto_surf, T)
+            dUdT = self.phase_param.dUdT(sto_surf)
+
+            # Bulk OCP is from the average SOC and temperature
+            sto_bulk = variables[f"{Domain} electrode {phase_name}stoichiometry"]
+            T_bulk = pybamm.xyz_average(pybamm.size_average(T))
+            ocp_bulk_eq = self.phase_param.U(sto_bulk, T_bulk)
+
+            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
+
+            H = self.phase_param.H(sto_surf)
+            variables[f"{Domain} electrode {phase_name}equilibrium OCP [V]"] = (
+                ocp_surf_eq
+            )
+            variables[f"{Domain} electrode {phase_name}bulk equilibrium OCP [V]"] = (
+                ocp_bulk_eq
+            )
+            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
+            variables[
+                f"{Domain} electrode {phase_name}hysteresis state distribution"
+            ] = h
+            H_x_av = pybamm.x_average(H)
+            h_x_av = pybamm.x_average(h)
+            # check if psd
+            if domain_options["particle size"] == "distribution":
+                if f"{domain} electrode" in sto_surf.domains["primary"]:
+                    ocp_surf = ocp_surf_eq + H * h
+                elif 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
+                    else:
+                        ocp_surf = ocp_surf_eq + H * h
+                else:
+                    ocp_surf = ocp_surf_eq + H_x_av * h_x_av
+            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
+
+        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 = self.param.current_with_time
+        current = variables[f"{Domain} electrode interfacial current density [A.m-2]"]
+        # current = current.orphans[0]
+        # current = current.SecondaryBroadcast(current,f'{domain} electrode')
+        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.K
+        K_x = self.phase_param.K_x
+        h = variables[f"{Domain} electrode {phase_name}hysteresis state"]
+        # h_x_av = variables[f'X-averaged {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"]
+        # h_av = variables[f'{Domain} electrode {phase_name}hysteresis state distribution']
+        # self.initial_conditions[h_av] = pybamm.Scalar(0)
+        self.initial_conditions[h] = pybamm.Scalar(0)

pybamm/parameters/lithium_ion_parameters.py

@@ -535,6 +535,12 @@ 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.K = pybamm.Parameter(f"{pref}{Domain} particle hysteresis decay rate")
+        self.K_x = 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)
@@ -630,6 +636,17 @@ def U(self, sto, T, lithiation=None):
             out.print_name = r"U_\mathrm{p}(c^\mathrm{surf}_\mathrm{s,p}, T)"
         return out
 
+    def H(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}
+        h_ref = pybamm.FunctionParameter(
+            f"{self.phase_prefactor}{Domain} electrode OCP hysteresis [V]", inputs
+        )
+        return h_ref
+
     def dUdT(self, sto):
         """
         Dimensional entropic change of the open-circuit potential [V.K-1]
@@ -645,6 +662,8 @@ def dUdT(self, sto):
             inputs,
         )
 
+    # def dQdU(self,sto)
+
     def X_j(self, index):
         "Available host sites indexed by reaction j"
         domain = self.domain

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

@@ -68,6 +68,25 @@ def test_current_sigmoid_ocp(self):
         modeltest = tests.StandardModelTest(model, parameter_values=parameter_values)
         modeltest.test_all(skip_output_tests=True)
 
+    def test_plett_ocp(self):
+        options = {"open-circuit potential": ("Plett", "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 OCP [V]" "": parameter_values[
+                    "Negative electrode OCP [V]"
+                ],
+                "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_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_plett_ocp(self):
+        options = {"open-circuit potential": "Plett"}
+        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,13 @@ def test_well_posed_current_sigmoid_ocp_with_psd(self):
         }
         self.check_well_posedness(options)
 
+    def test_well_posed_plett_ocp_with_psd(self):
+        options = {
+            "open-circuit potential": "Plett",
+            "particle size": "distribution",
+        }
+        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_plett_ocp(self):
+        options = {
+            "open-circuit potential": "Plett",
+        }
+        model = pybamm.lithium_ion.MPM(options)
+        model.check_well_posedness()
+
 
 class TestMPMExternalCircuits(TestCase):
     def test_well_posed_voltage(self):