Ecm with diffusion

CHANGELOG.md

@@ -1,5 +1,9 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
+## Features
+
+- Added the diffusion element to be used in the Thevenin model. ([#4254](https://github.com/pybamm-team/PyBaMM/pull/4254))
+
 # [v24.5rc0](https://github.com/pybamm-team/PyBaMM/tree/v24.5rc0) - 2024-05-01
 
 ## Features

examples/scripts/run_ecmd.py

@@ -0,0 +1,28 @@
+import pybamm
+
+pybamm.set_logging_level("INFO")
+
+model = pybamm.equivalent_circuit.Thevenin(options={"diffusion element": "true"})
+parameter_values = model.default_parameter_values
+
+parameter_values.update(
+    {"Diffusion time constant [s]": 580}, check_already_exists=False
+)
+
+experiment = pybamm.Experiment(
+    [
+        (
+            "Discharge at C/10 for 10 hours or until 3.3 V",
+            "Rest for 30 minutes",
+            "Rest for 2 hours",
+            "Charge at 100 A until 4.1 V",
+            "Hold at 4.1 V until 5 A",
+            "Rest for 30 minutes",
+            "Rest for 1 hour",
+        ),
+    ]
+)
+
+sim = pybamm.Simulation(model, experiment=experiment, parameter_values=parameter_values)
+sim.solve()
+sim.plot()

pybamm/CITATIONS.bib

@@ -794,3 +794,13 @@ @article{Wycisk2022
   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},
 }
+
+@article{Fan2022,
+  title={Data-driven identification of lithium-ion batteries: A nonlinear equivalent circuit model with diffusion dynamics},
+  author={Fan, Chuanxin and O’Regan, Kieran and Li, Liuying and Higgins, Matthew D and Kendrick, Emma and Widanage, Widanalage D},
+  journal={Applied Energy},
+  volume={321},
+  pages={119336},
+  year={2022},
+  publisher={Elsevier}
+}

pybamm/models/full_battery_models/equivalent_circuit/thevenin.py

@@ -30,6 +30,9 @@
                 throughput capacity in addition to discharge capacity. Must be one of
                 "true" or "false". "false" is the default, since calculating discharge
                 energy can be computationally expensive for simple models like SPM.
+            * "diffusion element" : str
+                Whether to include the diffusion element to the model. Must be one of
+                "true" or "false". "false" is the default.
             * "operating mode" : str
                 Sets the operating mode for the model. This determines how the current
                 is set. Can be:
@@ -73,6 +76,7 @@
     def set_options(self, extra_options=None):
         possible_options = {
             "calculate discharge energy": ["false", "true"],
+            "diffusion element": ["false", "true"],
             "operating mode": OperatingModes("current"),
             "number of rc elements": NaturalNumberOption(1),
         }
@@ -165,6 +169,18 @@
             )
             self.element_counter += 1
 
+    def set_diffusion_submodel(self):
+        if self.options["diffusion element"] == "false":
+            self.submodels["Diffusion"] = (
+                pybamm.equivalent_circuit_elements.NoDiffusion(self.param, self.options)
+            )
+        elif self.options["diffusion element"] == "true":
+            self.submodels["Diffusion"] = (
+                pybamm.equivalent_circuit_elements.DiffusionElement(
+                    self.param, self.options
+                )
+            )
+
     def set_thermal_submodel(self):
         self.submodels["Thermal"] = pybamm.equivalent_circuit_elements.ThermalSubModel(
             self.param, self.options
@@ -180,6 +196,7 @@
         self.set_ocv_submodel()
         self.set_resistor_submodel()
         self.set_rc_submodels()
+        self.set_diffusion_submodel()
         self.set_thermal_submodel()
         self.set_voltage_submodel()
 
@@ -227,3 +244,25 @@
                 "Irreversible heat generation [W]",
             ],
         ]
+
+    @property
+    def default_var_pts(self):
+        x = pybamm.SpatialVariable(
+            "x ECMD", domain=["ECMD particle"], coord_sys="Cartesian"
+        )
+        return {x: 20}
+
+    @property
+    def default_geometry(self):
+        x = pybamm.SpatialVariable(
+            "x ECMD", domain=["ECMD particle"], coord_sys="Cartesian"
+        )
+        return {"ECMD particle": {x: {"min": 0, "max": 1}}}
+
+    @property
+    def default_submesh_types(self):
+        return {"ECMD particle": pybamm.Uniform1DSubMesh}
+
+    @property
+    def default_spatial_methods(self):
+        return {"ECMD particle": pybamm.FiniteVolume()}

pybamm/models/submodels/equivalent_circuit_elements/__init__.py

@@ -3,6 +3,7 @@
 from .rc_element import RCElement
 from .thermal import ThermalSubModel
 from .voltage_model import VoltageModel
+from .diffusion_element import NoDiffusion, DiffusionElement
 
 __all__ = ['ocv_element', 'rc_element', 'resistor_element', 'thermal',
            'voltage_model']

pybamm/models/submodels/equivalent_circuit_elements/diffusion_element.py

@@ -0,0 +1,110 @@
+import pybamm
+
+
+class NoDiffusion(pybamm.BaseSubModel):
+    """
+    Without Diffusion element for
+    equivalent circuits.
+
+    Parameters
+    ----------
+    param : parameter class
+        The parameters to use for this submodel
+    options : dict, optional
+        A dictionary of options to be passed to the model.
+    """
+
+    def __init__(self, param, options=None):
+        super().__init__(param)
+        pybamm.citations.register("Fan2022")
+        self.model_options = options
+
+    def get_coupled_variables(self, variables):
+        z = pybamm.PrimaryBroadcast(variables["SoC"], "ECMD particle")
+        x = pybamm.SpatialVariable(
+            "x ECMD", domain=["ECMD particle"], coord_sys="Cartesian"
+        )
+        z_surf = pybamm.surf(z)
+        eta_diffusion = pybamm.Scalar(0)
+
+        variables.update(
+            {
+                "Distributed SoC": z,
+                "x ECMD": x,
+                "Diffusion overpotential [V]": eta_diffusion,
+                "Surface SoC": z_surf,
+            }
+        )
+
+        return variables
+
+
+class DiffusionElement(pybamm.BaseSubModel):
+    """
+    With Diffusion element for
+    equivalent circuits.
+
+    Parameters
+    ----------
+    param : parameter class
+        The parameters to use for this submodel
+    options : dict, optional
+        A dictionary of options to be passed to the model.
+    """
+
+    def __init__(self, param, options=None):
+        super().__init__(param)
+        self.model_options = options
+
+    def get_fundamental_variables(self):
+        z = pybamm.Variable("Distributed SoC", domain="ECMD particle")
+        x = pybamm.SpatialVariable(
+            "x ECMD", domain=["ECMD particle"], coord_sys="Cartesian"
+        )
+        variables = {
+            "Distributed SoC": z,
+            "x ECMD": x,
+        }
+        return variables
+
+    def get_coupled_variables(self, variables):
+        z = variables["Distributed SoC"]
+        soc = variables["SoC"]
+        z_surf = pybamm.surf(z)
+        eta_diffusion = -(self.param.ocv(z_surf) - self.param.ocv(soc))
+
+        variables.update(
+            {
+                "Diffusion overpotential [V]": eta_diffusion,
+                "Surface SoC": z_surf,
+            }
+        )
+
+        return variables
+
+    def set_rhs(self, variables):
+        cell_capacity = self.param.cell_capacity
+        current = variables["Current [A]"]
+        z = variables["Distributed SoC"]
+
+        # governing equations
+        dzdt = pybamm.div(pybamm.grad(z)) / self.param.tau_D
+        self.rhs = {z: dzdt}
+
+        # boundary conditions
+        lbc = pybamm.Scalar(0)
+        rbc = -self.param.tau_D * current / (cell_capacity * 3600)
+        self.boundary_conditions = {
+            z: {"left": (lbc, "Neumann"), "right": (rbc, "Neumann")}
+        }
+
+    def set_initial_conditions(self, variables):
+        z = variables["Distributed SoC"]
+        self.initial_conditions = {z: self.param.initial_soc}
+
+    def set_events(self, variables):
+        z_surf = variables["Surface SoC"]
+        self.events += [
+            pybamm.Event("Minimum surface SoC", z_surf),
+            pybamm.Event("Maximum surface SoC", 1 - z_surf),
+        ]

pybamm/models/submodels/equivalent_circuit_elements/ocv_element.py

@@ -56,7 +56,7 @@ def set_initial_conditions(self, variables):
 
     def set_events(self, variables):
         soc = variables["SoC"]
-        self.events = [
+        self.events += [
             pybamm.Event("Minimum SoC", soc),
             pybamm.Event("Maximum SoC", 1 - soc),
         ]

pybamm/models/submodels/equivalent_circuit_elements/voltage_model.py

@@ -30,7 +30,9 @@ def get_coupled_variables(self, variables):
         for i in range(number_of_elements):
             overpotential += variables[f"Element-{i} overpotential [V]"]
 
-        voltage = ocv + overpotential
+        diffusion_overpotential = variables["Diffusion overpotential [V]"]
+
+        voltage = ocv + overpotential + diffusion_overpotential
 
         # Power and Resistance
         current = variables["Current [A]"]

pybamm/parameters/ecm_parameters.py

@@ -4,6 +4,7 @@
 class EcmParameters:
     def __init__(self):
         self.cell_capacity = pybamm.Parameter("Cell capacity [A.h]")
+        self.tau_D = pybamm.Parameter("Diffusion time constant [s]")
 
         self._set_current_parameters()
         self._set_voltage_parameters()