Issue 1781 utilization

CHANGELOG.md

@@ -2,7 +2,8 @@
 
 ## Features
 
--   The name of a parameter set can be passed to `ParameterValues` as a string, e.g. `ParameterValues("Chen2020")`
+-   Added submodels for interface utilisation ([#1821](https://github.com/pybamm-team/PyBaMM/pull/1821))
+-   The name of a parameter set can be passed to `ParameterValues` as a string, e.g. `ParameterValues("Chen2020")` ([#1822](https://github.com/pybamm-team/PyBaMM/pull/1822))
 -   Reformatted SEI growth models into a single submodel with conditionals ([#1808](https://github.com/pybamm-team/PyBaMM/pull/1808))
 -   Stress-induced diffusion is now a separate model option instead of being automatically included when using the particle mechanics submodels ([#1797](https://github.com/pybamm-team/PyBaMM/pull/1797))
 -   `Experiment`s with drive cycles can be solved ([#1793](https://github.com/pybamm-team/PyBaMM/pull/1793))
@@ -19,7 +20,7 @@
 
 ## Breaking changes
 
--   The `chemistry` keyword argument in `ParameterValues` has been deprecated. Use `ParameterValues(chem)` instead of `ParameterValues(chemistry=chem)`
+-   The `chemistry` keyword argument in `ParameterValues` has been deprecated. Use `ParameterValues(chem)` instead of `ParameterValues(chemistry=chem)` ([#1822](https://github.com/pybamm-team/PyBaMM/pull/1822))
 -   Raise error when trying to convert an `Interpolant` with the "pchip" interpolator to CasADI ([#1791](https://github.com/pybamm-team/PyBaMM/pull/1791))
 -   Raise error if `Concatenation` is used directly with `Variable` objects (`concatenation` should be used instead) ([#1789](https://github.com/pybamm-team/PyBaMM/pull/1789))
 -   Made jax, jaxlib and the PyBaMM JaxSolver optional ([#1767](https://github.com/pybamm-team/PyBaMM/pull/1767))

examples/scripts/compare_interface_utilisation.py

@@ -0,0 +1,57 @@
+#
+# Example showing how to change the interface utilisation
+#
+
+import pybamm
+import numpy as np
+
+pybamm.set_logging_level("INFO")
+
+models = [
+    pybamm.lithium_ion.DFN({"interface utilisation": "full"}, name="full utilisation"),
+    pybamm.lithium_ion.DFN(
+        {"interface utilisation": "constant"}, name="constant utilisation"
+    ),
+    pybamm.lithium_ion.DFN(
+        {"interface utilisation": "current-driven"}, name="current-driven utilisation"
+    ),
+]
+
+param = models[0].default_parameter_values
+# add the user supplied parameters
+param.update(
+    {
+        "Initial negative electrode interface utilisation": 0.9,
+        "Initial positive electrode interface utilisation": 0.8,
+        "Negative electrode current-driven interface utilisation factor "
+        "[m3.mol-1]": -4e-5,
+        "Positive electrode current-driven interface utilisation factor "
+        "[m3.mol-1]": 4e-5,
+    },
+    check_already_exists=False,
+)
+
+# set up and solve simulations
+solutions = []
+t_eval = np.linspace(0, 3600, 100)
+
+for model in models:
+
+    sim = pybamm.Simulation(model, parameter_values=param)
+    solution = sim.solve(t_eval)
+    solutions.append(solution)
+
+# plot solutions
+pybamm.dynamic_plot(
+    solutions,
+    [
+        "Negative particle surface concentration [mol.m-3]",
+        "Positive particle surface concentration [mol.m-3]",
+        "Negative electrode interfacial current density [A.m-2]",
+        "Positive electrode interfacial current density [A.m-2]",
+        "Negative electrode interface utilisation",
+        "Positive electrode interface utilisation",
+        "Electrolyte potential [V]",
+        "Terminal voltage [V]",
+    ],
+)

examples/scripts/compare_lithium_ion_half_cell.py

@@ -3,7 +3,7 @@
 #
 import pybamm
 
-pybamm.set_logging_level("INFO")
+pybamm.set_logging_level("DEBUG")
 
 # load models
 models = [
@@ -20,4 +20,6 @@
     sims.append(sim)
 
 # plot
-pybamm.dynamic_plot(sims)
+pybamm.dynamic_plot(
+    sims, ["Lithium metal interface utilisation", "Terminal voltage [V]"]
+)

examples/scripts/custom_model.py

@@ -43,6 +43,12 @@
 model.submodels["positive interface"] = pybamm.interface.InverseButlerVolmer(
     model.param, "Positive", "lithium-ion main", options=model.options
 )
+model.submodels["negative interface utilisation"] = pybamm.interface_utilisation.Full(
+    model.param, "Negative", model.options
+)
+model.submodels["positive interface utilisation"] = pybamm.interface_utilisation.Full(
+    model.param, "Positive", model.options
+)
 model.submodels[
     "negative interface current"
 ] = pybamm.interface.CurrentForInverseButlerVolmer(

pybamm/__init__.py

@@ -149,6 +149,7 @@ def version(formatted=False):
 )
 from .models.submodels.interface import sei
 from .models.submodels.interface import lithium_plating
+from .models.submodels.interface import interface_utilisation
 
 #
 # Geometry

pybamm/models/full_battery_models/base_battery_model.py

@@ -47,7 +47,9 @@ class BatteryModelOptions(pybamm.FuzzyDict):
                 Whether to include hydrolysis in the model. Only implemented for
                 lead-acid models. Can be "false" (default) or "true". If "true", then
                 "surface form" cannot be 'false'.
-            * "lithium plating" : str, optional
+            * "interface utilisation": str
+                Can be "full" (default), "constant", or "current-driven".
+            * "lithium plating" : str
                 Sets the model for lithium plating. Can be "none" (default),
                 "reversible" or "irreversible".
             * "loss of active material" : str
@@ -69,9 +71,7 @@ class BatteryModelOptions(pybamm.FuzzyDict):
             * "particle shape" : str
                 Sets the model shape of the electrode particles. This is used to
                 calculate the surface area to volume ratio. Can be "spherical"
-                (default), "user" or "no particles". For the "user" option the surface
-                area per unit volume can be passed as a parameter, and is therefore not
-                necessarily consistent with the particle shape.
+                (default), or "no particles".
             * "particle size" : str
                 Sets the model to include a single active particle size or a
                 distribution of sizes at any macroscale location. Can be "single"
@@ -159,9 +159,10 @@ def __init__(self, extra_options):
                 "composite",
                 "integrated",
             ],
-            "hydrolysis": ["true", "false"],
+            "hydrolysis": ["false", "true"],
+            "interface utilisation": ["full", "constant", "current-driven"],
             "lithium plating": ["none", "reversible", "irreversible"],
-            "lithium plating porosity change": ["true", "false"],
+            "lithium plating porosity change": ["false", "true"],
             "loss of active material": ["none", "stress-driven", "reaction-driven"],
             "operating mode": ["current", "voltage", "power", "CCCV"],
             "particle": [
@@ -172,7 +173,7 @@ def __init__(self, extra_options):
                 "quartic profile",
             ],
             "particle mechanics": ["none", "swelling only", "swelling and cracking"],
-            "particle shape": ["spherical", "user", "no particles"],
+            "particle shape": ["spherical", "no particles"],
             "particle size": ["single", "distribution"],
             "SEI": [
                 "none",
@@ -184,37 +185,18 @@ def __init__(self, extra_options):
                 "ec reaction limited",
             ],
             "SEI film resistance": ["none", "distributed", "average"],
-            "SEI porosity change": ["true", "false"],
-            "stress-induced diffusion": ["true", "false"],
+            "SEI porosity change": ["false", "true"],
+            "stress-induced diffusion": ["false", "true"],
             "surface form": ["false", "differential", "algebraic"],
             "thermal": ["isothermal", "lumped", "x-lumped", "x-full"],
-            "total interfacial current density as a state": ["true", "false"],
+            "total interfacial current density as a state": ["false", "true"],
             "working electrode": ["both", "negative", "positive"],
         }
 
         default_options = {
-            "cell geometry": "none",
-            "convection": "none",
-            "current collector": "uniform",
-            "dimensionality": 0,
-            "electrolyte conductivity": "default",
-            "external submodels": [],
-            "hydrolysis": "false",
-            "lithium plating": "none",
-            "lithium plating porosity change": "false",
-            "loss of active material": "none",
-            "operating mode": "current",
-            "particle": "Fickian diffusion",
-            "particle mechanics": "none",
-            "particle shape": "spherical",
-            "particle size": "single",
-            "SEI": "none",
-            "SEI porosity change": "false",
-            "surface form": "false",
-            "thermal": "isothermal",
-            "total interfacial current density as a state": "false",
-            "working electrode": "both",
+            name: options[0] for name, options in self.possible_options.items()
         }
+        default_options["external submodels"] = []
 
         # Change the default for cell geometry based on which thermal option is provided
         extra_options = extra_options or {}
@@ -381,6 +363,7 @@ def __init__(self, extra_options):
                         (
                             option
                             in [
+                                "interface utilisation",
                                 "loss of active material",
                                 "particle mechanics",
                                 "particle",
@@ -872,6 +855,42 @@ def set_current_collector_submodel(self):
                 submodel = pybamm.current_collector.PotentialPair2plus1D(self.param)
         self.submodels["current collector"] = submodel
 
+    def set_interface_utilisation_submodel(self):
+        # this option can either be a string (both sides the same) or a 2-tuple
+        # to indicate different options in negative and positive electrodes
+        if isinstance(self.options["interface utilisation"], str):
+            util_left = self.options["interface utilisation"]
+            util_right = self.options["interface utilisation"]
+        else:
+            util_left, util_right = self.options["interface utilisation"]
+
+        if self.half_cell:
+            domains = [[util_left, "Counter"], [util_right, "Positive"]]
+        else:
+            domains = [[util_left, "Negative"], [util_right, "Positive"]]
+        for util, domain in domains:
+            name = domain.lower() + " interface utilisation"
+            if domain == "Counter":
+                domain = "Negative"
+            if util == "full":
+                self.submodels[name] = pybamm.interface_utilisation.Full(
+                    self.param, domain, self.options
+                )
+            elif util == "constant":
+                self.submodels[name] = pybamm.interface_utilisation.Constant(
+                    self.param, domain, self.options
+                )
+            elif util == "current-driven":
+                if self.half_cell and domain == "Negative":
+                    reaction_loc = "interface"
+                elif self.x_average:
+                    reaction_loc = "x-average"
+                else:
+                    reaction_loc = "full electrode"
+                self.submodels[name] = pybamm.interface_utilisation.CurrentDriven(
+                    self.param, domain, self.options, reaction_loc
+                )
+
     def set_voltage_variables(self):
 
         ocp_n = self.variables["Negative electrode open circuit potential"]

pybamm/models/full_battery_models/lead_acid/full.py

@@ -36,6 +36,7 @@ def __init__(self, options=None, name="Full model", build=True):
 
         self.set_external_circuit_submodel()
         self.set_interfacial_submodel()
+        self.set_interface_utilisation_submodel()
         self.set_porosity_submodel()
         self.set_active_material_submodel()
         self.set_tortuosity_submodels()

pybamm/models/full_battery_models/lead_acid/higher_order.py

@@ -36,6 +36,7 @@ def __init__(self, options=None, name="Composite model", build=True):
 
         self.set_external_circuit_submodel()
         self.set_leading_order_model()
+        self.set_interface_utilisation_submodel()
         # Electrolyte submodel to get first-order concentrations
         self.set_electrolyte_diffusion_submodel()
         self.set_other_species_diffusion_submodels()

pybamm/models/full_battery_models/lead_acid/loqs.py

@@ -35,6 +35,7 @@ def __init__(self, options=None, name="LOQS model", build=True):
 
         self.set_external_circuit_submodel()
         self.set_interfacial_submodel()
+        self.set_interface_utilisation_submodel()
         self.set_convection_submodel()
         self.set_porosity_submodel()
         self.set_active_material_submodel()

pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py

@@ -136,7 +136,10 @@ def set_degradation_variables(self):
         # Lithium lost to side reactions
         # Different way of measuring LLI but should give same value
         LLI_sei = self.variables["Loss of lithium to SEI [mol]"]
-        LLI_pl = self.variables["Loss of lithium to lithium plating [mol]"]
+        if self.half_cell:
+            LLI_pl = pybamm.Scalar(0)
+        else:
+            LLI_pl = self.variables["Loss of lithium to lithium plating [mol]"]
 
         LLI_reactions = LLI_sei + LLI_pl
         self.variables.update(
@@ -168,9 +171,7 @@ def set_summary_variables(self):
             "Total lithium lost from particles [mol]",
             "Total lithium lost from electrolyte [mol]",
             "Loss of lithium to SEI [mol]",
-            "Loss of lithium to lithium plating [mol]",
             "Loss of capacity to SEI [A.h]",
-            "Loss of capacity to lithium plating [A.h]",
             "Total lithium lost to side reactions [mol]",
             "Total capacity lost to side reactions [A.h]",
             # Resistance
@@ -182,6 +183,8 @@ def set_summary_variables(self):
                 "Negative electrode capacity [A.h]",
                 "Loss of active material in negative electrode [%]",
                 "Total lithium in negative electrode [mol]",
+                "Loss of lithium to lithium plating [mol]",
+                "Loss of capacity to lithium plating [A.h]",
             ]
 
         self.summary_variables = summary_variables
@@ -321,5 +324,7 @@ def set_li_metal_counter_electrode_submodels(self):
         # Models added specifically for the counter electrode have been labelled with
         # "counter electrode" so as not to be caught by this check
         self.submodels = {
-            k: v for k, v in self.submodels.items() if not k.startswith("negative")
+            k: v
+            for k, v in self.submodels.items()
+            if not (k.startswith("negative") or k == "lithium plating")
         }

pybamm/models/full_battery_models/lithium_ion/dfn.py

@@ -37,6 +37,7 @@ def __init__(self, options=None, name="Doyle-Fuller-Newman model", build=True):
 
         self.set_external_circuit_submodel()
         self.set_porosity_submodel()
+        self.set_interface_utilisation_submodel()
         self.set_crack_submodel()
         self.set_active_material_submodel()
         self.set_tortuosity_submodels()

pybamm/models/full_battery_models/lithium_ion/spm.py

@@ -37,6 +37,7 @@ def __init__(self, options=None, name="Single Particle Model", build=True):
 
         self.set_external_circuit_submodel()
         self.set_porosity_submodel()
+        self.set_interface_utilisation_submodel()
         self.set_crack_submodel()
         self.set_active_material_submodel()
         self.set_tortuosity_submodels()

pybamm/models/full_battery_models/lithium_ion/spme.py

@@ -40,6 +40,7 @@ def __init__(
 
         self.set_external_circuit_submodel()
         self.set_porosity_submodel()
+        self.set_interface_utilisation_submodel()
         self.set_crack_submodel()
         self.set_active_material_submodel()
         self.set_tortuosity_submodels()