Merge pull request #1913 from pybamm-team/surface-form-half-cell

Surface form half cell

CHANGELOG.md

@@ -2,8 +2,10 @@
 
 ## Features
 
+-   Half-cell models can now be run with "surface form" ([#1913](https://github.com/pybamm-team/PyBaMM/pull/1913))
+-   Added option for different kinetics on anode and cathode ([#1913](https://github.com/pybamm-team/PyBaMM/pull/1913))
+-   Allow `pybamm.Solution.save_data()` to return a string if filename is None, and added json to_format option ([#1909](https://github.com/pybamm-team/PyBaMM/pull/1909))
 -   Added an option to force install compatible versions of jax and jaxlib if already installed using CLI ([#1881](https://github.com/pybamm-team/PyBaMM/pull/1881))
--   Allow pybamm.Solution.save_data() to return a string if filename is None, and added json to_format option ([#1909](https://github.com/pybamm-team/PyBaMM/pull/1909)
 
 ## Optimizations
 

pybamm/input/parameters/lithium_ion/experiments/1C_discharge_from_full_Xu2019/parameters.csv

@@ -13,6 +13,6 @@ Lower voltage cut-off [V],3.5,,
 Upper voltage cut-off [V],4.2,,
 ,,,
 # Initial conditions
-Initial concentration in positive electrode [mol.m-3],4630,48230*0.096,
+Initial concentration in positive electrode [mol.m-3],4631,48230*0.096 adjusted so that initial voltage is <4.2V,
 Initial concentration in electrolyte [mol.m-3],1000,,
 Initial temperature [K],298.15,,

pybamm/models/full_battery_models/base_battery_model.py

@@ -52,6 +52,8 @@ class BatteryModelOptions(pybamm.FuzzyDict):
                 Model for intercalation kinetics. Can be "symmetric Butler-Volmer"
                 (default), "asymmetric Butler-Volmer", "linear", "Marcus", or
                 "Marcus-Hush-Chidsey" (which uses the asymptotic form from Zeng 2014).
+                A 2-tuple can be provided for different behaviour in negative and
+                positive electrodes.
             * "interface utilisation": str
                 Can be "full" (default), "constant", or "current-driven".
             * "lithium plating" : str
@@ -392,6 +394,7 @@ def __init__(self, extra_options):
                         (
                             option
                             in [
+                                "intercalation kinetics",
                                 "interface utilisation",
                                 "loss of active material",
                                 "particle mechanics",
@@ -852,17 +855,17 @@ def summary_variables(self, value):
     def set_summary_variables(self):
         self._summary_variables = []
 
-    @property
-    def intercalation_kinetics(self):
-        if self.options["intercalation kinetics"] == "symmetric Butler-Volmer":
+    def get_intercalation_kinetics(self, domain):
+        options = getattr(self.options, domain.lower())
+        if options["intercalation kinetics"] == "symmetric Butler-Volmer":
             return pybamm.kinetics.SymmetricButlerVolmer
-        elif self.options["intercalation kinetics"] == "asymmetric Butler-Volmer":
+        elif options["intercalation kinetics"] == "asymmetric Butler-Volmer":
             return pybamm.kinetics.AsymmetricButlerVolmer
-        elif self.options["intercalation kinetics"] == "linear":
+        elif options["intercalation kinetics"] == "linear":
             return pybamm.kinetics.Linear
-        elif self.options["intercalation kinetics"] == "Marcus":
+        elif options["intercalation kinetics"] == "Marcus":
             return pybamm.kinetics.Marcus
-        elif self.options["intercalation kinetics"] == "Marcus-Hush-Chidsey":
+        elif options["intercalation kinetics"] == "Marcus-Hush-Chidsey":
             return pybamm.kinetics.MarcusHushChidsey
 
     @property

pybamm/models/full_battery_models/lead_acid/full.py

@@ -81,12 +81,11 @@ def set_convection_submodel(self):
             ] = pybamm.convection.through_cell.Full(self.param)
 
     def set_intercalation_kinetics_submodel(self):
-        self.submodels["negative interface"] = self.intercalation_kinetics(
-            self.param, "Negative", "lead-acid main", self.options
-        )
-        self.submodels["positive interface"] = self.intercalation_kinetics(
-            self.param, "Positive", "lead-acid main", self.options
-        )
+        for domain in ["Negative", "Positive"]:
+            intercalation_kinetics = self.get_intercalation_kinetics(domain)
+            self.submodels[domain.lower() + " interface"] = intercalation_kinetics(
+                self.param, domain, "lead-acid main", self.options
+            )
 
     def set_solid_submodel(self):
         if self.options["surface form"] == "false":

pybamm/models/full_battery_models/lithium_ion/base_lithium_ion_model.py

@@ -294,7 +294,8 @@ def set_li_metal_counter_electrode_submodels(self):
             self.submodels[
                 "counter electrode potential"
             ] = pybamm.electrode.ohm.LithiumMetalSurfaceForm(self.param, self.options)
-            self.submodels["counter electrode interface"] = self.intercalation_kinetics(
+            neg_intercalation_kinetics = self.get_intercalation_kinetics("Negative")
+            self.submodels["counter electrode interface"] = neg_intercalation_kinetics(
                 self.param, "Negative", "lithium metal plating", self.options
             )
 

pybamm/models/full_battery_models/lithium_ion/dfn.py

@@ -71,13 +71,11 @@ def set_convection_submodel(self):
         ] = pybamm.convection.through_cell.NoConvection(self.param, self.options)
 
     def set_intercalation_kinetics_submodel(self):
-
-        self.submodels["negative interface"] = self.intercalation_kinetics(
-            self.param, "Negative", "lithium-ion main", self.options
-        )
-        self.submodels["positive interface"] = self.intercalation_kinetics(
-            self.param, "Positive", "lithium-ion main", self.options
-        )
+        for domain in ["Negative", "Positive"]:
+            intercalation_kinetics = self.get_intercalation_kinetics(domain)
+            self.submodels[domain.lower() + " interface"] = intercalation_kinetics(
+                self.param, domain, "lithium-ion main", self.options
+            )
 
     def set_particle_submodel(self):
         for domain in ["Negative", "Positive"]:
@@ -122,8 +120,12 @@ def set_solid_submodel(self):
             submod_n = pybamm.electrode.ohm.Full(self.param, "Negative", self.options)
             submod_p = pybamm.electrode.ohm.Full(self.param, "Positive", self.options)
         else:
-            submod_n = pybamm.electrode.ohm.SurfaceForm(self.param, "Negative")
-            submod_p = pybamm.electrode.ohm.SurfaceForm(self.param, "Positive")
+            submod_n = pybamm.electrode.ohm.SurfaceForm(
+                self.param, "Negative", self.options
+            )
+            submod_p = pybamm.electrode.ohm.SurfaceForm(
+                self.param, "Positive", self.options
+            )
 
         self.submodels["negative electrode potential"] = submod_n
         self.submodels["positive electrode potential"] = submod_p

pybamm/models/full_battery_models/lithium_ion/spm.py

@@ -98,13 +98,11 @@ def set_intercalation_kinetics_submodel(self):
                 self.param, "Positive", "lithium-ion main", self.options
             )
         else:
-            self.submodels["negative interface"] = self.intercalation_kinetics(
-                self.param, "Negative", "lithium-ion main", self.options
-            )
-
-            self.submodels["positive interface"] = self.intercalation_kinetics(
-                self.param, "Positive", "lithium-ion main", self.options
-            )
+            for domain in ["Negative", "Positive"]:
+                intercalation_kinetics = self.get_intercalation_kinetics(domain)
+                self.submodels[domain.lower() + " interface"] = intercalation_kinetics(
+                    self.param, domain, "lithium-ion main", self.options
+                )
 
     def set_particle_submodel(self):
         for domain in ["Negative", "Positive"]:

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

@@ -29,7 +29,8 @@ def _get_li_metal_interface_variables(self, delta_phi_s, phi_s, phi_e):
 
 
 class LithiumMetalSurfaceForm(LithiumMetalBaseModel):
-    """Explicit model for potential drop across a lithium metal electrode.
+    """Model for potential drop across a lithium metal electrode, with a
+    differential or algebraic equation for the surface potential difference
 
     Parameters
     ----------
@@ -83,15 +84,33 @@ def set_initial_conditions(self, variables):
 
         self.initial_conditions = {delta_phi: delta_phi_init}
 
+    def set_rhs(self, variables):
+        if self.options["surface form"] == "differential":
+            j_pl = variables["Lithium metal plating current density"]
+            j_sei = variables["SEI interfacial current density"]
+            sum_j = j_pl + j_sei
+
+            i_cc = variables["Current collector current density"]
+            delta_phi = variables[
+                "Lithium metal interface surface potential difference"
+            ]
+
+            C_dl = self.param.C_dl_n
+
+            self.rhs[delta_phi] = 1 / C_dl * (i_cc - sum_j)
+
     def set_algebraic(self, variables):
-        j_pl = variables["Lithium metal plating current density"]
-        j_sei = variables["SEI interfacial current density"]
-        sum_j = j_pl + j_sei
+        if self.options["surface form"] != "differential":  # also catches "false"
+            j_pl = variables["Lithium metal plating current density"]
+            j_sei = variables["SEI interfacial current density"]
+            sum_j = j_pl + j_sei
 
-        i_cc = variables["Current collector current density"]
-        delta_phi = variables["Lithium metal interface surface potential difference"]
+            i_cc = variables["Current collector current density"]
+            delta_phi = variables[
+                "Lithium metal interface surface potential difference"
+            ]
 
-        self.algebraic[delta_phi] = i_cc - sum_j
+            self.algebraic[delta_phi] = i_cc - sum_j
 
 
 class LithiumMetalExplicit(LithiumMetalBaseModel):

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

@@ -16,13 +16,15 @@ class SurfaceForm(BaseModel):
         The parameters to use for this submodel
     domain : str
         Either 'Negative' or 'Positive'
+    options : dict, optional
+        A dictionary of options to be passed to the model.
 
 
     **Extends:** :class:`pybamm.electrode.ohm.BaseModel`
     """
 
-    def __init__(self, param, domain):
-        super().__init__(param, domain)
+    def __init__(self, param, domain, options=None):
+        super().__init__(param, domain, options=options)
 
     def get_coupled_variables(self, variables):
 
@@ -56,9 +58,8 @@ def get_coupled_variables(self, variables):
         variables.update(self._get_standard_current_variables(i_s))
 
         if (
-            "Negative electrode current density" in variables
-            and "Positive electrode current density" in variables
-        ):
+            self.half_cell or "Negative electrode current density" in variables
+        ) and "Positive electrode current density" in variables:
             variables.update(self._get_standard_whole_cell_variables(variables))
 
         return variables

pybamm/models/submodels/electrolyte_conductivity/surface_potential_form/full_surface_form_conductivity.py

@@ -71,14 +71,18 @@ def get_coupled_variables(self, variables):
 
             i_boundary_cc = variables["Current collector current density"]
             c_e_s = variables["Separator electrolyte concentration"]
-            phi_e_n = variables["Negative electrolyte potential"]
+            if self.half_cell:
+                phi_e_n_s = variables["Lithium metal interface electrolyte potential"]
+            else:
+                phi_e_n = variables["Negative electrolyte potential"]
+                phi_e_n_s = pybamm.boundary_value(phi_e_n, "right")
             tor_s = variables["Separator porosity"]
             T = variables["Separator temperature"]
 
             chi_e_s = param.chi(c_e_s, T)
             kappa_s_eff = param.kappa_e(c_e_s, T) * tor_s
 
-            phi_e = pybamm.boundary_value(phi_e_n, "right") + pybamm.IndefiniteIntegral(
+            phi_e = phi_e_n_s + pybamm.IndefiniteIntegral(
                 (1 + param.Theta * T) * chi_e_s / c_e_s * pybamm.grad(c_e_s)
                 - param.C_e * i_boundary_cc / kappa_s_eff,
                 x_s,

pybamm/models/submodels/interface/interface_utilisation/base_utilisation.py

@@ -22,12 +22,21 @@ class BaseModel(pybamm.BaseSubModel):
     def __init__(self, param, domain, options):
         super().__init__(param, domain, options=options)
 
-    def _get_standard_interface_utilisation_variables(self, u):
-        u_av = pybamm.x_average(u)
+    def _get_standard_interface_utilisation_variables(self, u_var):
+        u = pybamm.maximum(u_var, 1e-8)
+        u_var_av = pybamm.x_average(u_var)
+        u_av = pybamm.maximum(u_var_av, 1e-8)
         if self.half_cell and self.domain == "Negative":
-            variables = {"Lithium metal interface utilisation": u_av}
+            variables = {
+                "Lithium metal interface utilisation variable": u_var_av,
+                "Lithium metal interface utilisation": u_av,
+            }
         else:
             variables = {
+                self.domain + " electrode interface utilisation variable": u_var,
+                "X-averaged "
+                + self.domain.lower()
+                + " electrode interface utilisation variable": u_var_av,
                 self.domain + " electrode interface utilisation": u,
                 "X-averaged "
                 + self.domain.lower()

pybamm/models/submodels/interface/interface_utilisation/current_driven_utilisation.py

@@ -32,31 +32,31 @@ def get_fundamental_variables(self):
         if self.reaction_loc == "full electrode":
             if self.domain == "Negative":
                 u = pybamm.Variable(
-                    "Negative electrode interface utilisation",
+                    "Negative electrode interface utilisation variable",
                     domain="negative electrode",
                     auxiliary_domains={"secondary": "current collector"},
                 )
             else:
                 u = pybamm.Variable(
-                    "Positive electrode interface utilisation",
+                    "Positive electrode interface utilisation variable",
                     domain="positive electrode",
                     auxiliary_domains={"secondary": "current collector"},
                 )
         elif self.reaction_loc == "x-average":
             if self.domain == "Negative":
                 u_xav = pybamm.Variable(
-                    "X-averaged negative electrode interface utilisation",
+                    "X-averaged negative electrode interface utilisation variable",
                     domain="current collector",
                 )
             else:
                 u_xav = pybamm.Variable(
-                    "X-averaged positive electrode interface utilisation",
+                    "X-averaged positive electrode interface utilisation variable",
                     domain="current collector",
                 )
             u = pybamm.PrimaryBroadcast(u_xav, self.domain.lower() + " electrode")
         else:
             u = pybamm.Variable(
-                "Lithium metal interface utilisation",
+                "Lithium metal interface utilisation variable",
                 domain="current collector",
             )
         variables = self._get_standard_interface_utilisation_variables(u)
@@ -65,12 +65,14 @@ def get_fundamental_variables(self):
 
     def set_rhs(self, variables):
         if self.reaction_loc == "full electrode":
-            u = variables[self.domain + " electrode interface utilisation"]
+            u = variables[self.domain + " electrode interface utilisation variable"]
             a = variables[self.domain + " electrode surface area to volume ratio"]
             j = variables[self.domain + " electrode interfacial current density"]
         elif self.reaction_loc == "x-average":
             u = variables[
-                "X-averaged " + self.domain.lower() + " electrode interface utilisation"
+                "X-averaged "
+                + self.domain.lower()
+                + " electrode interface utilisation variable"
             ]
             a = variables[
                 "X-averaged "
@@ -83,7 +85,7 @@ def set_rhs(self, variables):
                 + " electrode interfacial current density"
             ]
         else:
-            u = variables["Lithium metal interface utilisation"]
+            u = variables["Lithium metal interface utilisation variable"]
             a = 1
             j = variables["Lithium metal total interfacial current density"]
 
@@ -92,17 +94,19 @@ def set_rhs(self, variables):
         else:
             beta = self.param.beta_utilisation_p
 
-        self.rhs = {u: beta * a * j}
+        self.rhs = {u: beta * a * u * j}
 
     def set_initial_conditions(self, variables):
         if self.reaction_loc == "full electrode":
-            u = variables[self.domain + " electrode interface utilisation"]
+            u = variables[self.domain + " electrode interface utilisation variable"]
         elif self.reaction_loc == "x-average":
             u = variables[
-                "X-averaged " + self.domain.lower() + " electrode interface utilisation"
+                "X-averaged "
+                + self.domain.lower()
+                + " electrode interface utilisation variable"
             ]
         else:
-            u = variables["Lithium metal interface utilisation"]
+            u = variables["Lithium metal interface utilisation variable"]
 
         if self.domain == "Negative":
             u_init = self.param.u_n_init

pybamm/solvers/casadi_solver.py

@@ -433,7 +433,12 @@ def integer_bisect():
 
         y0 = coarse_solution.y[:, event_idx_lower]
         dense_step_sol = self._run_integrator(
-            model, y0, inputs_dict, inputs, t_window_event_dense, use_grid=use_grid,
+            model,
+            y0,
+            inputs_dict,
+            inputs,
+            t_window_event_dense,
+            use_grid=use_grid,
             extract_sensitivities_in_solution=False,
         )
 
@@ -649,7 +654,7 @@ def _run_integrator(
             Setting to True or False will override this behaviour, forcing the
             sensitivities to be extracted or not (it is up to the caller to determine if
             the sensitivities are in fact present)
-       """
+        """
 
         pybamm.logger.debug("Running CasADi integrator")
 

tests/integration/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_half_cell_tests.py

@@ -96,3 +96,11 @@ def test_current_driven_utilisation(self):
             check_already_exists=False,
         )
         self.run_basic_processing_test(options, parameter_values=parameter_values)
+
+    def test_surface_form_differential(self):
+        options = {"surface form": "differential"}
+        self.run_basic_processing_test(options)
+
+    def test_surface_form_algebraic(self):
+        options = {"surface form": "algebraic"}
+        self.run_basic_processing_test(options)