Ecm with diffusion

CHANGELOG.md

@@ -1,14 +1,15 @@
 # [Unreleased](https://github.com/pybamm-team/PyBaMM/)
 
-- Replaced rounded faraday constant with its exact value in `bpx.py` for better comparison between different tools
-
 ## Features
 
 - Added additional user-configurable options to the (`IDAKLUSolver`) and adjusted the default values to improve performance. ([#4282](https://github.com/pybamm-team/PyBaMM/pull/4282))
+- Added the diffusion element to be used in the Thevenin model. ([#4254](https://github.com/pybamm-team/PyBaMM/pull/4254))
 
 ## Optimizations
 
 - Improved performance and reliability of DAE consistent initialization. ([#4301](https://github.com/pybamm-team/PyBaMM/pull/4301))
+- Replaced rounded Faraday constant with its exact value in `bpx.py` for better comparison between different tools. ([#4290](https://github.com/pybamm-team/PyBaMM/pull/4290))
+
 ## Bug Fixes
 
 - Fixed bug where IDAKLU solver failed when `output variables` were specified and an event triggered. ([#4300](https://github.com/pybamm-team/PyBaMM/pull/4300))

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 @@ class Thevenin(pybamm.BaseModel):
                 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 __init__(
     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 @@ def set_rc_submodels(self):
             )
             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 @@ def set_submodels(self, build):
         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 @@ def default_quick_plot_variables(self):
                 "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)
+        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)
+        pybamm.citations.register("Fan2022")
+        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()

tests/integration/test_models/test_full_battery_models/test_equivalent_circuit/test_thevenin.py

@@ -7,3 +7,15 @@ def test_basic_processing(self):
         model = pybamm.equivalent_circuit.Thevenin()
         modeltest = tests.StandardModelTest(model)
         modeltest.test_all()
+
+    def test_diffusion(self):
+        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
+        )
+        modeltest = tests.StandardModelTest(model, parameter_values=parameter_values)
+        modeltest.test_all()

tests/unit/test_citations.py

@@ -345,6 +345,18 @@ def test_msmr(self):
         self.assertIn("Verbrugge2017", citations._papers_to_cite)
         self.assertIn("Verbrugge2017", citations._citation_tags.keys())
 
+    def test_thevenin(self):
+        citations = pybamm.citations
+
+        citations._reset()
+        pybamm.equivalent_circuit.Thevenin()
+        self.assertNotIn("Fan2022", citations._papers_to_cite)
+        self.assertNotIn("Fan2022", citations._citation_tags.keys())
+
+        pybamm.equivalent_circuit.Thevenin(options={"diffusion element": "true"})
+        self.assertIn("Fan2022", citations._papers_to_cite)
+        self.assertIn("Fan2022", citations._citation_tags.keys())
+
     def test_parameter_citations(self):
         citations = pybamm.citations
 

tests/unit/test_models/test_full_battery_models/test_equivalent_circuit/test_thevenin.py

@@ -11,6 +11,28 @@ def test_standard_model(self):
         model = pybamm.equivalent_circuit.Thevenin()
         model.check_well_posedness()
 
+    def test_default_properties(self):
+        model = pybamm.equivalent_circuit.Thevenin()
+        x = model.variables["x ECMD"]
+
+        # test var_pts
+        self.assertEqual(model.default_var_pts, {x: 20})
+
+        # test geometry
+        self.assertEqual(
+            model.default_geometry, {"ECMD particle": {x: {"min": 0, "max": 1}}}
+        )
+
+        # test spatial methods
+        self.assertIsInstance(
+            model.default_spatial_methods["ECMD particle"], pybamm.FiniteVolume
+        )
+
+        # test submesh types
+        self.assertEqual(
+            model.default_submesh_types, {"ECMD particle": pybamm.Uniform1DSubMesh}
+        )
+
     def test_changing_number_of_rcs(self):
         options = {"number of rc elements": 0}
         model = pybamm.equivalent_circuit.Thevenin(options=options)
@@ -33,6 +55,11 @@ def test_changing_number_of_rcs(self):
             model = pybamm.equivalent_circuit.Thevenin(options=options)
             model.check_well_posedness()
 
+    def test_diffusion_element(self):
+        options = {"diffusion element": "true"}
+        model = pybamm.equivalent_circuit.Thevenin(options=options)
+        model.check_well_posedness(post_discretisation=True)
+
     def test_calculate_discharge_energy(self):
         options = {"calculate discharge energy": "true"}
         model = pybamm.equivalent_circuit.Thevenin(options=options)