li metal model runs with EC SEI but solver fails

examples/scripts/DFN_half_cell.py

@@ -5,15 +5,15 @@
 import pybamm
 import numpy as np
 
-pybamm.set_logging_level("INFO")
+pybamm.set_logging_level("DEBUG")
 
 # load model
-options = {"working electrode": "positive"}
+options = {"working electrode": "positive", "SEI": "ec reaction limited"}
 model1 = pybamm.lithium_ion.DFN(options=options)
-model2 = pybamm.lithium_ion.BasicDFNHalfCell(options=options)
+# model2 = pybamm.lithium_ion.BasicDFNHalfCell(options=options)
 
 sols = []
-for model in [model1, model2]:
+for model in [model1]:
     # create geometry
     geometry = model.default_geometry
 

pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py

@@ -119,52 +119,47 @@ def set_degradation_variables(self):
         )
 
     def set_sei_submodel(self):
+        if self.half_cell:
+            reaction_loc = "interface"
+        elif self.x_average:
+            reaction_loc = "x-average"
+        else:
+            reaction_loc = "full electrode"
 
-        # SEI
         if self.options["SEI"] == "none":
             self.submodels["sei"] = pybamm.sei.NoSEI(self.param, self.options)
-
-        if self.options["SEI"] == "constant":
+        elif self.options["SEI"] == "constant":
             self.submodels["sei"] = pybamm.sei.ConstantSEI(self.param)
-
         elif self.options["SEI"] == "reaction limited":
             self.submodels["sei"] = pybamm.sei.ReactionLimited(
-                self.param, self.x_average
+                self.param, reaction_loc, self.options
             )
-
         elif self.options["SEI"] == "solvent-diffusion limited":
             self.submodels["sei"] = pybamm.sei.SolventDiffusionLimited(
-                self.param, self.x_average
+                self.param, reaction_loc, self.options
             )
-
         elif self.options["SEI"] == "electron-migration limited":
             self.submodels["sei"] = pybamm.sei.ElectronMigrationLimited(
-                self.param, self.x_average
+                self.param, reaction_loc, self.options
             )
-
         elif self.options["SEI"] == "interstitial-diffusion limited":
             self.submodels["sei"] = pybamm.sei.InterstitialDiffusionLimited(
-                self.param, self.x_average
+                self.param, reaction_loc, self.options
             )
-
         elif self.options["SEI"] == "ec reaction limited":
             self.submodels["sei"] = pybamm.sei.EcReactionLimited(
-                self.param, self.x_average
+                self.param, reaction_loc, self.options
             )
 
     def set_lithium_plating_submodel(self):
-
-        # negative electrode
         if self.options["lithium plating"] == "none":
             self.submodels["lithium plating"] = pybamm.lithium_plating.NoPlating(
                 self.param, self.options
             )
-
         elif self.options["lithium plating"] == "reversible":
             self.submodels[
                 "lithium plating"
             ] = pybamm.lithium_plating.ReversiblePlating(self.param, self.x_average)
-
         elif self.options["lithium plating"] == "irreversible":
             self.submodels[
                 "lithium plating"

pybamm/models/full_battery_models/lithium_ion/dfn.py

@@ -90,16 +90,23 @@ def set_interfacial_submodel(self):
         # Set the counter-electrode model for the half-cell model
         # The negative electrode model will be ignored
         if self.half_cell:
-            self.submodels[
-                "counter electrode interface"
-            ] = pybamm.interface.InverseButlerVolmer(
-                self.param, "Negative", "lithium metal plating", self.options
-            )  # assuming symmetric reaction for now so we can take the inverse
-            self.submodels[
-                "counter electrode interface current"
-            ] = pybamm.interface.CurrentForInverseButlerVolmerLithiumMetal(
-                self.param, "Negative", "lithium metal plating", self.options
-            )
+            if self.options["surface form"] == "false":
+                self.submodels[
+                    "counter electrode interface"
+                ] = pybamm.interface.InverseButlerVolmer(
+                    self.param, "Negative", "lithium metal plating", self.options
+                )  # assuming symmetric reaction for now so we can take the inverse
+                self.submodels[
+                    "counter electrode interface current"
+                ] = pybamm.interface.CurrentForInverseButlerVolmerLithiumMetal(
+                    self.param, "Negative", "lithium metal plating", self.options
+                )
+            else:
+                self.submodels[
+                    "counter electrode interface"
+                ] = pybamm.interface.ButlerVolmer(
+                    self.param, "Negative", "lithium metal plating", self.options
+                )
 
     def set_particle_submodel(self):
 
@@ -162,9 +169,16 @@ def set_solid_submodel(self):
         # Set the counter-electrode model for the half-cell model
         # The negative electrode model will be ignored
         if self.half_cell:
-            self.submodels[
-                "counter electrode potential"
-            ] = pybamm.electrode.ohm.LithiumMetalExplicit(self.param, self.options)
+            if self.options["SEI"] in ["none", "constant"]:
+                self.submodels[
+                    "counter electrode potential"
+                ] = pybamm.electrode.ohm.LithiumMetalExplicit(self.param, self.options)
+            else:
+                self.submodels[
+                    "counter electrode potential"
+                ] = pybamm.electrode.ohm.LithiumMetalSurfaceForm(
+                    self.param, self.options
+                )
 
     def set_electrolyte_submodel(self):
 

pybamm/models/standard_variables.py

@@ -310,6 +310,15 @@ def __init__(self):
             domain=["negative electrode"],
             auxiliary_domains={"secondary": "current collector"},
         )
+        # For SEI reaction at the li metal/separator interface in a li metal model
+        self.L_inner_interface = pybamm.Variable(
+            "Inner SEI thickness",
+            domain=["current collector"],
+        )
+        self.L_outer_interface = pybamm.Variable(
+            "Outer SEI thickness",
+            domain=["current collector"],
+        )
 
     def __setattr__(self, name, value):
         value.print_name = name

pybamm/models/submodels/electrode/ohm/__init__.py

@@ -3,4 +3,4 @@
 from .full_ohm import Full
 from .leading_ohm import LeadingOrder
 from .surface_form_ohm import SurfaceForm
-from .li_metal_explicit import LithiumMetalExplicit
+from .li_metal import LithiumMetalExplicit, LithiumMetalSurfaceForm

pybamm/models/submodels/electrode/ohm/base_ohm.py

@@ -25,6 +25,8 @@ def __init__(self, param, domain, options=None, set_positive_potential=True):
         super().__init__(param, domain, options, set_positive_potential)
 
     def set_boundary_conditions(self, variables):
+        if self.half_cell:
+            return
 
         if self.domain == "Negative":
             phi_s_cn = variables["Negative current collector potential"]

pybamm/models/submodels/electrode/ohm/li_metal.py

@@ -0,0 +1,136 @@
+#
+# Subodels for a lithium metal electrode
+#
+import pybamm
+from .base_ohm import BaseModel
+
+
+class LithiumMetalSurfaceForm(BaseModel):
+    """Explicit model for potential drop across a lithium metal electrode.
+
+    Parameters
+    ----------
+    param : parameter class
+        The parameters to use for this submodel
+    options : dict, optional
+        A dictionary of options to be passed to the model.
+
+    **Extends:** :class:`pybamm.electrode.ohm.BaseModel`
+    """
+
+    def __init__(self, param, options=None):
+        super().__init__(param, "Negative", options=options)
+
+    def get_fundamental_variables(self):
+        ocp_ref = self.param.U_n_ref
+        pot_scale = self.param.potential_scale
+
+        delta_phi = pybamm.Variable(
+            "Lithium metal interface surface potential difference",
+            domain="current collector",
+        )
+        variables = {
+            "Lithium metal interface surface potential difference": delta_phi,
+            "Lithium metal interface surface potential difference [V]": ocp_ref
+            + delta_phi * pot_scale,
+        }
+
+        return variables
+
+    def get_coupled_variables(self, variables):
+        param = self.param
+
+        i_boundary_cc = variables["Current collector current density"]
+        T_n = variables["Negative current collector temperature"]
+        l_n = param.l_n
+        delta_phi_s = i_boundary_cc * l_n / param.sigma_n(T_n)
+        delta_phi_s_dim = param.potential_scale * delta_phi_s
+
+        phi_s_cn = variables["Negative current collector potential"]
+        delta_phi = variables["Lithium metal interface surface potential difference"]
+
+        # Potentials at the anode/separator interface
+        phi_s = phi_s_cn - delta_phi_s
+        phi_e = phi_s - delta_phi
+
+        variables.update(
+            {
+                "Negative electrode potential drop": delta_phi_s,
+                "Negative electrode potential drop [V]": delta_phi_s_dim,
+                "X-averaged negative electrode ohmic losses": delta_phi_s / 2,
+                "X-averaged negative electrode ohmic losses [V]": delta_phi_s_dim / 2,
+                "Lithium metal interface electrode potential": phi_s,
+                "Lithium metal interface electrolyte potential": phi_e,
+            }
+        )
+        return variables
+
+    def set_initial_conditions(self, variables):
+        delta_phi = variables["Lithium metal interface surface potential difference"]
+        delta_phi_init = self.param.U_n_init
+
+        self.initial_conditions = {delta_phi: delta_phi_init}
+
+    def set_algebraic(self, variables):
+        sum_j = variables[
+            "Sum of x-averaged "
+            + self.domain.lower()
+            + " electrode interfacial current densities"
+        ]
+
+        sum_j_av = variables[
+            "X-averaged "
+            + self.domain.lower()
+            + " electrode total interfacial current density"
+        ]
+        delta_phi = variables["Lithium metal interface surface potential difference"]
+
+        self.algebraic[delta_phi] = sum_j_av - sum_j
+
+
+class LithiumMetalExplicit(BaseModel):
+    """Explicit model for potential drop across a lithium metal electrode.
+
+    Parameters
+    ----------
+    param : parameteslackr class
+        The parameters to use for this submodel
+    options : dict, optional
+        A dictionary of options to be passed to the model.
+
+    **Extends:** :class:`pybamm.electrode.ohm.BaseModel`
+    """
+
+    def __init__(self, param, options=None):
+        super().__init__(param, "Negative", options=options)
+
+    def get_coupled_variables(self, variables):
+        param = self.param
+
+        i_boundary_cc = variables["Current collector current density"]
+        T_n = variables["Negative current collector temperature"]
+        l_n = param.l_n
+        delta_phi_s = i_boundary_cc * l_n / param.sigma_n(T_n)
+        delta_phi_s_dim = param.potential_scale * delta_phi_s
+
+        phi_s_cn = variables["Negative current collector potential"]
+        delta_phi = variables["Lithium metal interface surface potential difference"]
+
+        phi_s = phi_s_cn - delta_phi_s
+        phi_e = phi_s - delta_phi
+
+        variables.update(
+            {
+                "Negative electrode potential drop": delta_phi_s,
+                "Negative electrode potential drop [V]": delta_phi_s_dim,
+                "X-averaged negative electrode ohmic losses": delta_phi_s / 2,
+                "X-averaged negative electrode ohmic losses [V]": delta_phi_s_dim / 2,
+                "Lithium metal interface electrode potential": phi_s,
+                "Lithium metal interface electrolyte potential": phi_e,
+            }
+        )
+
+        return variables
+
+    def set_boundary_conditions(self, variables):
+        pass

pybamm/models/submodels/electrolyte_conductivity/base_electrolyte_conductivity.py

@@ -379,12 +379,7 @@ def set_boundary_conditions(self, variables):
         phi_e = variables["Electrolyte potential"]
 
         if self.half_cell:
-            phi_s_cn = variables["Negative current collector potential"]
-            delta_phi_s = variables["Negative electrode potential drop"]
-            delta_phi = variables["Negative electrode surface potential difference"]
-
-            phi_s = phi_s_cn - delta_phi_s
-            phi_e_ref = phi_s - delta_phi
+            phi_e_ref = variables["Lithium metal interface electrolyte potential"]
             lbc = (phi_e_ref, "Dirichlet")
         else:
             lbc = (pybamm.Scalar(0), "Neumann")