Ecm with diffusion

examples/scripts/Compare_ECM_ECMD.py

@@ -0,0 +1,59 @@
+import pybamm
+
+model1 = pybamm.equivalent_circuit.Thevenin(
+    options={"diffusion element": "false"}
+)  # Model1: ECM Model definition
+model2 = pybamm.equivalent_circuit.Thevenin(
+    options={"diffusion element": "true"}
+)  # Model2: ECMD Model definition
+
+parameter_values1 = model1.default_parameter_values
+parameter_values2 = model2.default_parameter_values
+parameter_values2.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",
+        ),
+    ]
+)
+
+models = [
+    model1,
+    model2,
+]
+
+sims = []
+for model, parameter_values in zip(
+    [model1, model2], [parameter_values1, parameter_values2]
+):
+    sim = pybamm.Simulation(
+        model, experiment=experiment, parameter_values=parameter_values
+    )
+    sim.solve()  # Solve for specific time range (0 to 3600 seconds)
+    sims.append(sim)
+
+labels = [
+    "Thevenin (No Diffusion)",
+    "Thevenin (Diffusion, τ=580s)",  # Adjust labels as needed
+]
+pybamm.dynamic_plot(sims, time_unit="seconds", labels=labels)
+
+default_params = model1.default_parameter_values
+nominal_capacity = default_params["Nominal cell capacity [A.h]"]
+print(f"Nominal Cell Capacity (Ah): {nominal_capacity}")
+
+# Second plot to compare SoC and Surface SoC
+output_variables = ["SoC", "Surface SoC"]
+
+sim_diffusion = sims[1]  # Assuming model2 is at index 1 in sims list
+sim_diffusion.plot(output_variables=output_variables)

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/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:
@@ -59,7 +62,7 @@ class Thevenin(pybamm.BaseModel):
     """
 
     def __init__(
-        self, name="Thevenin Equivalent Circuit Model", options=None, build=True
+        self, name="Thevenin Equivalent Circuit Modelx", options=None, build=True
     ):
         super().__init__(name)
 
@@ -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,109 @@
+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)
+        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()