#834 correcting errors on reaction limited model

examples/scripts/sei_growth.py

@@ -0,0 +1,13 @@
+import pybamm as pb
+
+pb.set_logging_level("DEBUG")
+pb.settings.debug_mode = True
+
+options = {"sei": "reaction limited"}
+model = pb.lithium_ion.DFN(options)
+
+sim = pb.Simulation(model)
+
+solver = pb.CasadiSolver(mode="safe")
+sim.solve(solver=solver)
+sim.plot(["Terminal voltage [V]", "Total SEI thickness [m]"])

input/parameters/lithium-ion/seis/Example/README.md

@@ -0,0 +1,9 @@
+# SEI parameters
+
+Some example parameters for SEI growth from the papers:
+
+> Ramadass, P., Haran, B., Gomadam, P. M., White, R., & Popov, B. N. (2004). Development of first principles capacity fade model for Li-ion cells. Journal of the Electrochemical Society, 151(2), A196-A203.
+> Ploehn, H. J., Ramadass, P., & White, R. E. (2004). Solvent diffusion model for aging of lithium-ion battery cells. Journal of The Electrochemical Society, 151(3), A456-A462.
+> Single, F., Latz, A., & Horstmann, B. (2018). Identifying the mechanism of continued growth of the solid–electrolyte interphase. ChemSusChem, 11(12), 1950-1955.
+
+Note: this parameter set does not claim to be representative of the true parameter values. Instead these are parameter values that were used to fit SEI models to observed experimental data in the referenced papers.

input/parameters/lithium-ion/seis/Example/parameters.csv

@@ -0,0 +1,19 @@
+Name [units],Value,Reference,Notes
+# Empty rows and rows starting with ‘#’ will be ignored,,,
+,,,
+# SEI properties,,,
+Inner SEI reaction proportion,0.5,,
+Inner SEI partial molar volume [m3.mol-1],95.86E-18, use same value,
+Outer SEI partial molar volume [m3.mol-1],95.86E-18, use same value,
+SEI reaction exchange current density [A.m-2],1.5E-6, Ramadass paper,
+SEI resistance per unit thickness [Ohm.m-1],1,Ramadass paper,
+Outer SEI solvent diffusivity [m2.s-1],2.5E-22,Single paper,
+Bulk solvent concentration [mol.m-3],2.636E3, Ploehn paper,
+Ratio of inner and outer SEI exchange current densities,1, Assume same,
+Inner SEI open-circuit potential [V],0.1,,
+Outer SEI open-circuit potential [V],0.8,,
+Inner SEI electron conducticity [S.m-1],8.95E-14, Single paper,
+Inner SEI lithium interstitial diffusivity [m2.s-1],1E-15, Single paper,
+Lithium interstitial reference concentration [mol.m-3],15, Single paper,
+Initial inner SEI thickness [m], 7.5E-9, 1/2 of initial thickness in Single paper,
+Initial outer SEI thickness [m], 7.5E-9, 1/2 of initial thickness in Single paper,

pybamm/__init__.py

@@ -174,6 +174,7 @@ def version(formatted=False):
     porosity,
     thermal,
     tortuosity,
+    sei,
 )
 
 #
@@ -202,6 +203,7 @@ def version(formatted=False):
 from .parameters import electrical_parameters
 from .parameters import thermal_parameters
 from .parameters import standard_parameters_lithium_ion, standard_parameters_lead_acid
+from .parameters import sei_parameters
 from .parameters.print_parameters import print_parameters, print_evaluated_parameters
 from .parameters import parameter_sets
 

pybamm/models/full_battery_models/base_battery_model.py

@@ -60,6 +60,11 @@ class BaseBatteryModel(pybamm.BaseModel):
                 variable for instead of solving in PyBaMM. The entries of the lists
                 are strings that correspond to the submodel names in the keys
                 of `self.submodels`.
+            * "sei" : str
+                Set the sei submodel to be used. If None, no SEI model is use. Options
+                are "reaction limited", "solvent diffusion limited",
+                "electron migration limited", or "lithium interstitial
+                diffusion limited".
 
 
     **Extends:** :class:`pybamm.BaseModel`
@@ -155,6 +160,7 @@ def options(self, extra_options):
             "thermal": "isothermal",
             "thermal current collector": False,
             "external submodels": [],
+            "sei": None,
         }
         options = default_options
         # any extra options overwrite the default options

pybamm/models/full_battery_models/lithium_ion/dfn.py

@@ -44,6 +44,7 @@ def __init__(self, options=None, name="Doyle-Fuller-Newman model", build=True):
         self.set_electrolyte_submodel()
         self.set_thermal_submodel()
         self.set_current_collector_submodel()
+        self.set_sei_submodel()
 
         if build:
             self.build_model()
@@ -102,6 +103,10 @@ def set_electrolyte_submodel(self):
             self.param, self.reactions
         )
 
+    def set_sei_submodel(self):
+        if self.options["sei"] == "reaction limited":
+            self.submodels["sei"] = pybamm.sei.ReactionLimited(self.param)
+
     @property
     def default_geometry(self):
         dimensionality = self.options["dimensionality"]

pybamm/models/submodels/sei/__init__.py

@@ -0,0 +1,2 @@
+from .base_sei import BaseModel
+from .reaction_limited import ReactionLimited

pybamm/models/submodels/sei/base_sei.py

@@ -37,37 +37,95 @@ def _get_standard_thickness_variables(self, L_inner, L_outer):
         variables : dict
             The variables which can be derived from the SEI thicknesses.
         """
+        sp = pybamm.sei_parameters
 
         L_inner_av = pybamm.x_average(L_inner)
         L_outer_av = pybamm.x_average(L_outer)
 
+        L_sei = L_inner + L_outer
+        L_sei_av = pybamm.x_average(L_sei)
+
         n_inner = L_inner  # inner SEI concentration
         n_outer = L_outer  # outer SEI concentration
         n_inner_av = pybamm.x_average(L_inner)
         n_outer_av = pybamm.x_average(L_outer)
 
-        n_SEI = n_inner + n_outer  # SEI concentration
+        n_SEI = n_inner + n_outer / sp.v_bar  # SEI concentration
         n_SEI_av = pybamm.x_average(n_SEI)
 
-        Q_sei = n_SEI_av * self.param.L_n * self.param.L_y * self.L_z
+        Q_sei = n_SEI_av * self.param.L_n * self.param.L_y * self.param.L_z
+
+        L_scale = sp.L_sei_0_dim
+        n_scale = sp.L_sei_0 * sp.a_n / sp.V_bar_inner_dimensional
+        n_outer_scale = sp.L_sei_0 * sp.a_n / sp.V_bar_outer_dimensional
 
         variables = {
             "Inner SEI thickness": L_inner,
-            "Inner SEI thickness [m]": L_inner,
+            "Inner SEI thickness [m]": L_inner * L_scale,
             "X-averaged inner SEI thickness": L_inner_av,
-            "X-averaged inner SEI thickness [m]": L_inner_av,
+            "X-averaged inner SEI thickness [m]": L_inner_av * L_scale,
             "Outer SEI thickness": L_outer,
-            "Outer SEI thickness [m]": L_outer,
+            "Outer SEI thickness [m]": L_outer * L_scale,
             "X-averaged outer SEI thickness": L_outer_av,
-            "X-averaged outer SEI thickness [m]": L_outer_av,
-            "Inner SEI concentration [mol.m-3]": n_inner,
-            "X-averaged inner SEI concentration [mol.m-3]": n_inner_av,
-            "Outer SEI concentration [mol.m-3]": n_outer,
-            "X-averaged outer SEI concentration [mol.m-3]": n_outer_av,
-            "SEI concentration [mol.m-3]": n_SEI,
-            "X-averaged SEI concentration [mol.m-3]": n_SEI_av,
-            "Loss of lithium to SEI [mols]": Q_sei,
+            "X-averaged outer SEI thickness [m]": L_outer_av * L_scale,
+            "Total SEI thickness": L_sei,
+            "Total SEI thickness [m]": L_sei * L_scale,
+            "X-averaged total SEI thickness": L_sei_av,
+            "X-averaged total SEI thickness [m]": L_sei_av * L_scale,
+            "Inner SEI concentration [mol.m-3]": n_inner * n_scale,
+            "X-averaged inner SEI concentration [mol.m-3]": n_inner_av * n_scale,
+            "Outer SEI concentration [mol.m-3]": n_outer * n_outer_scale,
+            "X-averaged outer SEI concentration [mol.m-3]": n_outer_av * n_outer_scale,
+            "SEI concentration [mol.m-3]": n_SEI * n_scale,
+            "X-averaged SEI concentration [mol.m-3]": n_SEI_av * n_scale,
+            "Loss of lithium to SEI [mols]": Q_sei * n_scale,
         }
 
         return variables
 
+    def _get_standard_reaction_variables(self, j_inner, j_outer):
+        """
+            A private function to obtain the standard variables which
+            can be derived from the SEI interfacial reaction current
+
+            Parameters
+            ----------
+            j_inner : :class:`pybamm.Symbol`
+                The inner SEI interfacial reaction current.
+            j_outer : :class:`pybamm.Symbol`
+                The outer SEI interfacial reaction current.
+
+            Returns
+            -------
+            variables : dict
+                The variables which can be derived from the SEI thicknesses.
+        """
+
+        j_i_av = pybamm.x_average(j_inner)
+        j_o_av = pybamm.x_average(j_outer)
+
+        j_sei = j_inner + j_outer
+        j_sei_av = pybamm.x_average(j_sei)
+
+        sp = pybamm.sei_parameters
+        j_scale = sp.F * sp.L_sei_0 / sp.V_bar_inner_dimensional / sp.tau_discharge
+
+        variables = {
+            "Inner SEI reaction interfacial current density": j_inner,
+            "Inner SEI reaction interfacial current density [A.m-2]": j_inner * j_scale,
+            "X-averaged inner SEI reaction interfacial current density": j_i_av,
+            "X-averaged inner SEI reaction interfacial current density [A.m-2]": j_i_av
+            * j_scale,
+            "Outer SEI reaction interfacial current density": j_outer,
+            "Outer SEI reaction interfacial current density [A.m-2]": j_outer * j_scale,
+            "X-averaged outer SEI reaction interfacial current density": j_o_av,
+            "X-averaged outer SEI reaction interfacial current density [A.m-2]": j_o_av
+            * j_scale,
+            "SEI reaction interfacial current density": j_sei,
+            "SEI reaction interfacial current density [A.m-2]": j_sei * j_scale,
+            "X-averaged SEI reaction interfacial current density": j_sei_av,
+            "X-averaged SEI reaction interfacial current density [A.m-2]": j_sei_av
+            * j_scale,
+        }
+
+        return variables

pybamm/models/submodels/sei/reaction_limited.py

@@ -26,32 +26,46 @@ def get_fundamental_variables(self):
         L_inner = pybamm.standard_variables.L_inner
         L_outer = pybamm.standard_variables.L_outer
 
-        variables = self.get_standard_thickness_variables(L_inner, L_outer)
+        variables = self._get_standard_thickness_variables(L_inner, L_outer)
 
         return variables
 
     def get_coupled_variables(self, variables):
-        c_s = variables[self.domain + " particle concentration"]
-        T_k = pybamm.PrimaryBroadcast(
-            variables[self.domain + " electrode temperature"],
-            [self.domain.lower() + " particle"],
-        )
+        phi_s_n = variables["Negative electrode potential"]
+        phi_e_n = variables["Negative electrolyte potential"]
+        j_n = variables["Negative electrode interfacial current density"]
+        L_sei = variables["Total SEI thickness"]
 
-        if self.domain == "Negative":
-            N_s = -self.param.D_n(c_s, T_k) * pybamm.grad(c_s)
-        elif self.domain == "Positive":
-            N_s = -self.param.D_p(c_s, T_k) * pybamm.grad(c_s)
+        C_sei = pybamm.sei_parameters.C_sei_reaction
+        R_sei = pybamm.sei_parameters.R_sei
+        alpha = 0.5
+        # alpha = pybamm.sei_parameters.alpha
 
-        variables.update(self._get_standard_flux_variables(N_s, N_s))
+        # need to revise for thermal case
+        j_sei = (1 / C_sei) * pybamm.exp(-(phi_s_n - phi_e_n - j_n * L_sei * R_sei))
 
-        if self.domain == "Negative":
-            x = pybamm.standard_spatial_vars.x_n
-            R = pybamm.FunctionParameter("Negative particle distribution in x", x)
-            variables.update({"Negative particle distribution in x": R})
+        j_inner = alpha * j_sei
+        j_outer = (1 - alpha) * j_sei
 
-        elif self.domain == "Positive":
-            x = pybamm.standard_spatial_vars.x_p
-            R = pybamm.FunctionParameter("Positive particle distribution in x", x)
-            variables.update({"Positive particle distribution in x": R})
+        variables.update(self._get_standard_reaction_variables(j_inner, j_outer))
 
         return variables
+
+    def set_rhs(self, variables):
+        L_inner = variables["Inner SEI thickness"]
+        L_outer = variables["Outer SEI thickness"]
+        j_inner = variables["Inner SEI reaction interfacial current density"]
+        j_outer = variables["Outer SEI reaction interfacial current density"]
+
+        v_bar = pybamm.sei_parameters.v_bar
+
+        self.rhs = {L_inner: j_inner, L_outer: v_bar * j_outer}
+
+    def set_initial_conditions(self, variables):
+        L_inner = variables["Inner SEI thickness"]
+        L_outer = variables["Outer SEI thickness"]
+
+        L_inner_0 = pybamm.sei_parameters.L_inner_0
+        L_outer_0 = pybamm.sei_parameters.L_outer_0
+
+        self.initial_conditions = {L_inner: L_inner_0, L_outer: L_outer_0}

pybamm/parameters/parameter_sets.py

@@ -37,6 +37,7 @@
     "cathode": "lico2_Marquis2019",
     "electrolyte": "lipf6_Marquis2019",
     "experiment": "1C_discharge_from_full_Marquis2019",
+    "sei": "Example",
 }
 
 NCA_Kim2011 = {
@@ -47,6 +48,7 @@
     "cathode": "nca_Kim2011",
     "electrolyte": "lipf6_Kim2011",
     "experiment": "1C_discharge_from_full_Kim2011",
+    "sei": "Example",
 }
 
 #

pybamm/parameters/parameter_values.py

@@ -109,14 +109,21 @@ def update_from_chemistry(self, chemistry):
         # Create path to file
         path = os.path.join(pybamm.root_dir(), "input", "parameters", base_chemistry)
         # Load each component name
-        for component_group in [
+
+        component_groups = [
             "cell",
             "anode",
             "cathode",
             "separator",
             "electrolyte",
             "experiment",
-        ]:
+        ]
+
+        # add sei parameters if provided
+        if "sei" in chemistry:
+            component_groups += ["sei"]
+
+        for component_group in component_groups:
             # Make sure component is provided
             try:
                 component = chemistry[component_group]

pybamm/parameters/sei_parameters.py

@@ -39,8 +39,32 @@ def alpha(m_ratio, phi_s, phi_e, U_inner, U_outer):
     "Lithium interstitial reference concentration [mol.m-3]"
 )
 
+L_inner_0_dim = pybamm.Parameter("Initial inner SEI thickness [m]")
+L_outer_0_dim = pybamm.Parameter("Initial outer SEI thickness [m]")
+
+L_sei_0_dim = L_inner_0_dim + L_outer_0_dim
+
 # --------------------------------------------------------------------------------------
 # Dimensionless parameters
 
-# write as C_SEI_reaction etc ...
+U_n_ref = pybamm.standard_parameters_lithium_ion.U_n_ref
+F = pybamm.standard_parameters_lithium_ion.F
+R = pybamm.standard_parameters_lithium_ion.R
+tau_discharge = pybamm.standard_parameters_lithium_ion.tau_discharge
+T_ref = pybamm.standard_parameters_lithium_ion.T_ref
+
+a_n = pybamm.standard_parameters_lithium_ion.a_n_dim
+L_x = pybamm.standard_parameters_lithium_ion.L_x
+
+I_typ = pybamm.electrical_parameters.I_typ
+
+C_sei_reaction = (
+    F * L_sei_0_dim / (m_sei_dimensional * tau_discharge * V_bar_inner_dimensional)
+) * pybamm.exp(-(F * U_n_ref / (2 * R * T_ref)))
+
+R_sei = F * I_typ * R_sei_dimensional / (a_n * L_x)
+
+v_bar = V_bar_outer_dimensional / V_bar_inner_dimensional
 
+L_inner_0 = L_inner_0_dim / L_sei_0_dim
+L_outer_0 = L_outer_0_dim / L_sei_0_dim