Issue 579 volume fractions

CHANGELOG.md

@@ -2,6 +2,8 @@
 
 ## Features
 
+-   Separated active and inactive material volume fractions ([#726](https://github.com/pybamm-team/PyBaMM/pull/726))
+-   Added submodels for tortuosity ([#726](https://github.com/pybamm-team/PyBaMM/pull/726))
 -   Simplified the interface for setting current functions ([#723](https://github.com/pybamm-team/PyBaMM/pull/723))
 -   Added Heaviside operator ([#723](https://github.com/pybamm-team/PyBaMM/pull/723))
 -   Added Simulation class ([#693](https://github.com/pybamm-team/PyBaMM/pull/693))

docs/source/models/submodels/index.rst

@@ -14,3 +14,4 @@ Submodels
     particle/index
     porosity/index
     thermal/index
+    tortuosity/index

docs/source/models/submodels/tortuosity/base_tortuosity.rst

@@ -0,0 +1,6 @@
+Base Model
+===========
+
+.. autoclass:: pybamm.tortuosity.BaseModel
+    :members:
+

docs/source/models/submodels/tortuosity/bruggeman_tortuosity.rst

@@ -0,0 +1,8 @@
+Bruggeman Model 
+===============
+
+.. autoclass:: pybamm.tortuosity.Bruggeman
+    :members:
+
+
+

docs/source/models/submodels/tortuosity/index.rst

@@ -0,0 +1,9 @@
+Tortuosity
+==========
+
+.. toctree::
+  :maxdepth: 1
+
+  base_tortuosity
+  bruggeman_tortuosity
+

input/parameters/lithium-ion/anodes/graphite_mcmb2528_Marquis2019/parameters.csv

@@ -8,7 +8,8 @@ Negative electrode diffusivity [m2.s-1],[function]graphite_mcmb2528_diffusivity_
 Negative electrode OCP [V],[function]graphite_mcmb2528_ocp_Dualfoil1998,
 ,,,
 # Microstructure,,,
-Negative electrode porosity,0.3,Scott Moura FastDFN,
+Negative electrode porosity,0.3,Scott Moura FastDFN,electrolyte volume fraction
+Negative electrode active material volume fraction,0.7,,assuming zero binder volume fraction
 Negative particle radius [m],1E-05,Scott Moura FastDFN,
 Negative electrode surface area density [m-1],180000,Scott Moura FastDFN,
 Negative electrode Bruggeman coefficient,1.5,Scott Moura FastDFN,

input/parameters/lithium-ion/cathodes/lico2_Marquis2019/parameters.csv

@@ -8,7 +8,8 @@ Positive electrode diffusivity [m2.s-1],[function]lico2_diffusivity_Dualfoil1998
 Positive electrode OCP [V],[function]lico2_ocp_Dualfoil1998,
 ,,,
 # Microstructure,,,
-Positive electrode porosity,0.3,Scott Moura FastDFN,
+Positive electrode porosity,0.3,Scott Moura FastDFN,electrolyte volume fraction
+Positive electrode active material volume fraction,0.7,,assuming zero binder volume fraction
 Positive particle radius [m],1E-05,Scott Moura FastDFN,
 Positive electrode surface area density [m-1],150000,Scott Moura FastDFN,
 Positive electrode Bruggeman coefficient,1.5,Scott Moura FastDFN,

pybamm/__init__.py

@@ -164,6 +164,7 @@ def version(formatted=False):
     particle,
     porosity,
     thermal,
+    tortuosity,
 )
 
 #

pybamm/models/full_battery_models/base_battery_model.py

@@ -415,7 +415,9 @@ def build_model(self):
         # return to it later and try again. If setting coupled variables fails and
         # there are no more submodels to try, raise an error.
         submodels = list(self.submodels.keys())
+        count = 0
         while len(submodels) > 0:
+            count += 1
             for submodel_name, submodel in self.submodels.items():
                 if submodel_name in submodels:
                     pybamm.logger.debug(
@@ -429,7 +431,7 @@ def build_model(self):
                         )
                         submodels.remove(submodel_name)
                     except KeyError as key:
-                        if len(submodels) == 1:
+                        if len(submodels) == 1 or count == 100:
                             # no more submodels to try
                             raise pybamm.ModelError(
                                 """Submodel "{}" requires the variable {}, but it cannot be found.
@@ -491,6 +493,14 @@ def build_model(self):
 
         self._built = True
 
+    def set_tortuosity_submodels(self):
+        self.submodels["electrolyte tortuosity"] = pybamm.tortuosity.Bruggeman(
+            self.param, "Electrolyte"
+        )
+        self.submodels["electrode tortuosity"] = pybamm.tortuosity.Bruggeman(
+            self.param, "Electrode"
+        )
+
     def set_thermal_submodel(self):
 
         if self.options["thermal"] == "isothermal":

pybamm/models/full_battery_models/lead_acid/full.py

@@ -37,6 +37,7 @@ def __init__(self, options=None, name="Full model", build=True):
         self.set_reactions()
         self.set_interfacial_submodel()
         self.set_porosity_submodel()
+        self.set_tortuosity_submodels()
         self.set_convection_submodel()
         self.set_electrolyte_submodel()
         self.set_solid_submodel()

pybamm/models/full_battery_models/lead_acid/higher_order.py

@@ -49,6 +49,7 @@ def __init__(self, options=None, name="Composite model", build=True):
         self.set_full_interface_submodel()
         self.set_full_convection_submodel()
         self.set_full_porosity_submodel()
+        self.set_tortuosity_submodels()
         self.set_thermal_submodel()
         self.set_current_collector_submodel()
 

pybamm/models/full_battery_models/lead_acid/loqs.py

@@ -37,6 +37,7 @@ def __init__(self, options=None, name="LOQS model", build=True):
         self.set_interfacial_submodel()
         self.set_convection_submodel()
         self.set_porosity_submodel()
+        self.set_tortuosity_submodels()
         self.set_negative_electrode_submodel()
         self.set_electrolyte_submodel()
         self.set_positive_electrode_submodel()
@@ -70,6 +71,14 @@ def set_porosity_submodel(self):
             self.param
         )
 
+    def set_tortuosity_submodels(self):
+        self.submodels[
+            "leading-order electrolyte tortuosity"
+        ] = pybamm.tortuosity.Bruggeman(self.param, "Electrolyte")
+        self.submodels[
+            "leading-order electrode tortuosity"
+        ] = pybamm.tortuosity.Bruggeman(self.param, "Electrode")
+
     def set_convection_submodel(self):
 
         if self.options["convection"] is False:

pybamm/models/full_battery_models/lithium_ion/dfn.py

@@ -35,6 +35,7 @@ def __init__(self, options=None, name="Doyle-Fuller-Newman model", build=True):
 
         self.set_reactions()
         self.set_porosity_submodel()
+        self.set_tortuosity_submodels()
         self.set_convection_submodel()
         self.set_interfacial_submodel()
         self.set_particle_submodel()

pybamm/models/full_battery_models/lithium_ion/spm.py

@@ -33,6 +33,7 @@ def __init__(self, options=None, name="Single Particle Model", build=True):
         super().__init__(options, name)
 
         self.set_porosity_submodel()
+        self.set_tortuosity_submodels()
         self.set_convection_submodel()
         self.set_interfacial_submodel()
         self.set_particle_submodel()

pybamm/models/full_battery_models/lithium_ion/spme.py

@@ -37,6 +37,7 @@ def __init__(
 
         self.set_reactions()
         self.set_porosity_submodel()
+        self.set_tortuosity_submodels()
         self.set_convection_submodel()
         self.set_interfacial_submodel()
         self.set_particle_submodel()
@@ -53,6 +54,14 @@ def set_porosity_submodel(self):
 
         self.submodels["porosity"] = pybamm.porosity.Constant(self.param)
 
+    def set_tortuosity_submodels(self):
+        self.submodels["electrolyte tortuosity"] = pybamm.tortuosity.Bruggeman(
+            self.param, "Electrolyte", True
+        )
+        self.submodels["electrode tortuosity"] = pybamm.tortuosity.Bruggeman(
+            self.param, "Electrode", True
+        )
+
     def set_convection_submodel(self):
 
         self.submodels["convection"] = pybamm.convection.NoConvection(self.param)

pybamm/models/submodels/base_submodel.py

@@ -68,13 +68,25 @@ def domain(self):
 
     @domain.setter
     def domain(self, domain):
-        if domain in ["Negative", "Separator", "Positive"]:
+        ok_domain_list = [
+            "Negative",
+            "Separator",
+            "Positive",
+            "Negative electrode",
+            "Negative electrolyte",
+            "Separator electrolyte",
+            "Positive electrode",
+            "Positive electrolyte",
+        ]
+        if domain in ok_domain_list:
             self._domain = domain
         elif domain is None:
             pass
         else:
             raise pybamm.DomainError(
-                "Domain must be either 'Negative' or 'Positive' not {}".format(domain)
+                "Domain '{}' not recognised (must be one of {})".format(
+                    domain, ok_domain_list
+                )
             )
 
     def set_domain_for_broadcast(self):

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

@@ -25,21 +25,19 @@ def __init__(self, param, domain, reactions=None, set_positive_potential=True):
 
     def set_boundary_conditions(self, variables):
 
-        phi_s = variables[self.domain + " electrode potential"]
-        eps = variables[self.domain + " electrode porosity"]
-        i_boundary_cc = variables["Current collector current density"]
-        phi_s_cn = variables["Negative current collector potential"]
-
         if self.domain == "Negative":
+            phi_s_cn = variables["Negative current collector potential"]
             lbc = (phi_s_cn, "Dirichlet")
             rbc = (pybamm.Scalar(0), "Neumann")
 
         elif self.domain == "Positive":
             lbc = (pybamm.Scalar(0), "Neumann")
-            sigma_eff = self.param.sigma_p * (1 - eps) ** self.param.b_p
+            i_boundary_cc = variables["Current collector current density"]
+            sigma_eff = self.param.sigma_p * variables["Positive electrode tortuosity"]
             rbc = (
                 i_boundary_cc / pybamm.boundary_value(-sigma_eff, "right"),
                 "Neumann",
             )
 
+        phi_s = variables[self.domain + " electrode potential"]
         self.boundary_conditions[phi_s] = {"left": lbc, "right": rbc}

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

@@ -34,13 +34,13 @@ def get_coupled_variables(self, variables):
         x_n = pybamm.standard_spatial_vars.x_n
         x_p = pybamm.standard_spatial_vars.x_p
 
-        eps_0 = variables[
-            "Leading-order x-averaged " + self.domain.lower() + " electrode porosity"
+        tor_0 = variables[
+            "Leading-order x-averaged " + self.domain.lower() + " electrode tortuosity"
         ]
         phi_s_cn = variables["Negative current collector potential"]
 
         if self._domain == "Negative":
-            sigma_eff_0 = self.param.sigma_n * (1 - eps_0) ** self.param.b_n
+            sigma_eff_0 = self.param.sigma_n * tor_0
             phi_s = pybamm.PrimaryBroadcast(
                 phi_s_cn, "negative electrode"
             ) + pybamm.outer(
@@ -54,7 +54,7 @@ def get_coupled_variables(self, variables):
             ]
             phi_e_p_av = variables["X-averaged positive electrolyte potential"]
 
-            sigma_eff_0 = self.param.sigma_p * (1 - eps_0) ** self.param.b_p
+            sigma_eff_0 = self.param.sigma_p * tor_0
 
             const = (
                 delta_phi_p_av
@@ -80,8 +80,8 @@ def get_coupled_variables(self, variables):
     def set_boundary_conditions(self, variables):
 
         phi_s = variables[self.domain + " electrode potential"]
-        eps_0 = variables[
-            "Leading-order x-averaged " + self.domain.lower() + " electrode porosity"
+        tor_0 = variables[
+            "Leading-order x-averaged " + self.domain.lower() + " electrode tortuosity"
         ]
         i_boundary_cc_0 = variables["Leading-order current collector current density"]
 
@@ -91,7 +91,7 @@ def set_boundary_conditions(self, variables):
 
         elif self.domain == "Positive":
             lbc = (pybamm.Scalar(0), "Neumann")
-            sigma_eff_0 = self.param.sigma_p * (1 - eps_0) ** self.param.b_p
+            sigma_eff_0 = self.param.sigma_p * tor_0
             rbc = (-i_boundary_cc_0 / sigma_eff_0, "Neumann")
 
         self.boundary_conditions[phi_s] = {"left": lbc, "right": rbc}

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

@@ -36,16 +36,14 @@ def get_fundamental_variables(self):
     def get_coupled_variables(self, variables):
 
         phi_s = variables[self.domain + " electrode potential"]
-        eps = variables[self.domain + " electrode porosity"]
+        tor = variables[self.domain + " electrode tortuosity"]
 
         if self.domain == "Negative":
             sigma = self.param.sigma_n
-            b = self.param.b_n
         elif self.domain == "Positive":
             sigma = self.param.sigma_p
-            b = self.param.b_p
 
-        sigma_eff = sigma * (1 - eps) ** b
+        sigma_eff = sigma * tor
         i_s = -sigma_eff * pybamm.grad(phi_s)
 
         variables.update({self.domain + " electrode effective conductivity": sigma_eff})
@@ -72,7 +70,7 @@ def set_boundary_conditions(self, variables):
 
         phi_s = variables[self.domain + " electrode potential"]
         phi_s_cn = variables["Negative current collector potential"]
-        eps = variables[self.domain + " electrode porosity"]
+        tor = variables[self.domain + " electrode tortuosity"]
         i_boundary_cc = variables["Current collector current density"]
 
         if self.domain == "Negative":
@@ -81,7 +79,7 @@ def set_boundary_conditions(self, variables):
 
         elif self.domain == "Positive":
             lbc = (pybamm.Scalar(0), "Neumann")
-            sigma_eff = self.param.sigma_p * (1 - eps) ** self.param.b_p
+            sigma_eff = self.param.sigma_p * tor
             rbc = (
                 i_boundary_cc / pybamm.boundary_value(-sigma_eff, "right"),
                 "Neumann",

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

@@ -31,13 +31,13 @@ def get_coupled_variables(self, variables):
         x_p = pybamm.standard_spatial_vars.x_p
         i_boundary_cc = variables["Current collector current density"]
         i_e = variables[self.domain + " electrolyte current density"]
-        eps = variables[self.domain + " electrode porosity"]
+        tor = variables[self.domain + " electrode tortuosity"]
         phi_s_cn = variables["Negative current collector potential"]
 
         i_s = pybamm.PrimaryBroadcast(i_boundary_cc, self.domain_for_broadcast) - i_e
 
         if self.domain == "Negative":
-            conductivity = param.sigma_n * (1 - eps) ** param.b_n
+            conductivity = param.sigma_n * tor
             phi_s = pybamm.PrimaryBroadcast(
                 phi_s_cn, "negative electrode"
             ) - pybamm.IndefiniteIntegral(i_s / conductivity, x_n)
@@ -47,7 +47,7 @@ def get_coupled_variables(self, variables):
             phi_e_s = variables["Separator electrolyte potential"]
             delta_phi_p = variables["Positive electrode surface potential difference"]
 
-            conductivity = param.sigma_p * (1 - eps) ** param.b_p
+            conductivity = param.sigma_p * tor
             phi_s = -pybamm.IndefiniteIntegral(
                 i_s / conductivity, x_p
             ) + pybamm.PrimaryBroadcast(

pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/base_higher_order_stefan_maxwell_conductivity.py

@@ -44,9 +44,9 @@ def get_coupled_variables(self, variables):
             "X-averaged negative electrode surface potential difference"
         ]
         phi_s_n_av = variables["X-averaged negative electrode potential"]
-        eps_n_av = variables["Leading-order x-averaged negative electrode porosity"]
-        eps_s_av = variables["Leading-order x-averaged separator porosity"]
-        eps_p_av = variables["Leading-order x-averaged positive electrode porosity"]
+        tor_n_av = variables["Leading-order x-averaged negative electrolyte tortuosity"]
+        tor_s_av = variables["Leading-order x-averaged separator tortuosity"]
+        tor_p_av = variables["Leading-order x-averaged positive electrolyte tortuosity"]
 
         T_av = variables["X-averaged cell temperature"]
         T_av_n = pybamm.PrimaryBroadcast(T_av, "negative electrode")
@@ -63,9 +63,9 @@ def get_coupled_variables(self, variables):
         x_p = pybamm.standard_spatial_vars.x_p
 
         # bulk conductivities
-        kappa_n_av = param.kappa_e(c_e_av, T_av) * eps_n_av ** param.b_n
-        kappa_s_av = param.kappa_e(c_e_av, T_av) * eps_s_av ** param.b_s
-        kappa_p_av = param.kappa_e(c_e_av, T_av) * eps_p_av ** param.b_p
+        kappa_n_av = param.kappa_e(c_e_av, T_av) * tor_n_av
+        kappa_s_av = param.kappa_e(c_e_av, T_av) * tor_s_av
+        kappa_p_av = param.kappa_e(c_e_av, T_av) * tor_p_av
 
         chi_av = param.chi(c_e_av)
         if chi_av.domain == ["current collector"]: