Repository: pybamm-team/PyBaMM
Pull request #731 (2019-11-13 21:57:25)

Issue 699 remove autograd
Accidentally merged 579 into 699 instead of master but now can't push straight to master ...



Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index fb7475c74..e8b6f28f3 100644
--- a/CHANGELOG.md
+++ b/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))
diff --git a/docs/source/models/submodels/index.rst b/docs/source/models/submodels/index.rst
index ff2cc824d..5377dcd86 100644
--- a/docs/source/models/submodels/index.rst
+++ b/docs/source/models/submodels/index.rst
@@ -14,3 +14,4 @@ Submodels
     particle/index
     porosity/index
     thermal/index
+    tortuosity/index
diff --git a/docs/source/models/submodels/tortuosity/base_tortuosity.rst b/docs/source/models/submodels/tortuosity/base_tortuosity.rst
new file mode 100644
index 000000000..0d3a2261f
--- /dev/null
+++ b/docs/source/models/submodels/tortuosity/base_tortuosity.rst
@@ -0,0 +1,6 @@
+Base Model
+===========
+
+.. autoclass:: pybamm.tortuosity.BaseModel
+    :members:
+
diff --git a/docs/source/models/submodels/tortuosity/bruggeman_tortuosity.rst b/docs/source/models/submodels/tortuosity/bruggeman_tortuosity.rst
new file mode 100644
index 000000000..0e7297377
--- /dev/null
+++ b/docs/source/models/submodels/tortuosity/bruggeman_tortuosity.rst
@@ -0,0 +1,8 @@
+Bruggeman Model 
+===============
+
+.. autoclass:: pybamm.tortuosity.Bruggeman
+    :members:
+
+
+
diff --git a/docs/source/models/submodels/tortuosity/index.rst b/docs/source/models/submodels/tortuosity/index.rst
new file mode 100644
index 000000000..595a9bb30
--- /dev/null
+++ b/docs/source/models/submodels/tortuosity/index.rst
@@ -0,0 +1,9 @@
+Tortuosity
+==========
+
+.. toctree::
+  :maxdepth: 1
+
+  base_tortuosity
+  bruggeman_tortuosity
+
diff --git a/input/parameters/lead-acid/anodes/lead_Sulzer2019/parameters.csv b/input/parameters/lead-acid/anodes/lead_Sulzer2019/parameters.csv
index 82b3a8125..56fc70558 100644
--- a/input/parameters/lead-acid/anodes/lead_Sulzer2019/parameters.csv
+++ b/input/parameters/lead-acid/anodes/lead_Sulzer2019/parameters.csv
@@ -12,7 +12,8 @@ Negative electrode open-circuit potential [V],[function]lead_ocp_Bode1977,,
 ,,,
 # Microstructure,,,
 Negative electrode surface area density [m-1],2300000,,
-Negative electrode Bruggeman coefficient,1.5,,
+Negative electrode Bruggeman coefficient (electrolyte),1.5,,
+Negative electrode Bruggeman coefficient (electrode),1.5,,
 Negative electrode morphological parameter,0.6,srinivasan2003mathematical,
 Negative electrode capacity [C.m-3],3473000000,,
 ,,,
diff --git a/input/parameters/lead-acid/cathodes/lead_dioxide_Sulzer2019/parameters.csv b/input/parameters/lead-acid/cathodes/lead_dioxide_Sulzer2019/parameters.csv
index ca1c45498..cd38f78b4 100644
--- a/input/parameters/lead-acid/cathodes/lead_dioxide_Sulzer2019/parameters.csv
+++ b/input/parameters/lead-acid/cathodes/lead_dioxide_Sulzer2019/parameters.csv
@@ -12,7 +12,8 @@ Positive electrode open-circuit potential [V],[function]lead_dioxide_ocp_Bode197
 ,,,
 # Microstructure,,,
 Positive electrode surface area density [m-1],23000000,,
-Positive electrode Bruggeman coefficient,1.5,,
+Positive electrode Bruggeman coefficient (electrolyte),1.5,,
+Positive electrode Bruggeman coefficient (electrode),1.5,,
 Positive electrode morphological parameter,0.6,srinivasan2003mathematical,
 Positive electrode capacity [C.m-3],2745000000,,
 ,,,
diff --git a/input/parameters/lead-acid/separators/agm_Sulzer2019/parameters.csv b/input/parameters/lead-acid/separators/agm_Sulzer2019/parameters.csv
index 8625bf67d..539bddd46 100644
--- a/input/parameters/lead-acid/separators/agm_Sulzer2019/parameters.csv
+++ b/input/parameters/lead-acid/separators/agm_Sulzer2019/parameters.csv
@@ -2,4 +2,5 @@ Name [units],Value,Reference,Notes
 # Empty rows and rows starting with ‘#’ will be ignored,,,
 ,,,
 Maximum porosity of separator,0.92,,
-Separator Bruggeman coefficient,1.5,,
+Separator Bruggeman coefficient (electrolyte),1.5,,
+Separator Bruggeman coefficient (electrode),1.5,,
diff --git a/input/parameters/lithium-ion/anodes/graphite_mcmb2528_Marquis2019/parameters.csv b/input/parameters/lithium-ion/anodes/graphite_mcmb2528_Marquis2019/parameters.csv
index 28a1733cb..3d1469d95 100644
--- a/input/parameters/lithium-ion/anodes/graphite_mcmb2528_Marquis2019/parameters.csv
+++ b/input/parameters/lithium-ion/anodes/graphite_mcmb2528_Marquis2019/parameters.csv
@@ -8,10 +8,12 @@ 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,
+Negative electrode Bruggeman coefficient (electrolyte),1.5,Scott Moura FastDFN,
+Negative electrode Bruggeman coefficient (electrode),1.5,Scott Moura FastDFN,
 ,,,
 # Interfacial reactions,,,
 Negative electrode cation signed stoichiometry,-1,,
diff --git a/input/parameters/lithium-ion/cathodes/lico2_Marquis2019/parameters.csv b/input/parameters/lithium-ion/cathodes/lico2_Marquis2019/parameters.csv
index 99f2b3840..d1fabddca 100644
--- a/input/parameters/lithium-ion/cathodes/lico2_Marquis2019/parameters.csv
+++ b/input/parameters/lithium-ion/cathodes/lico2_Marquis2019/parameters.csv
@@ -8,10 +8,12 @@ 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,
+Positive electrode Bruggeman coefficient (electrolyte),1.5,Scott Moura FastDFN,
+Positive electrode Bruggeman coefficient (electrode),1.5,Scott Moura FastDFN,
 ,,,
 # Interfacial reactions,,,
 Positive electrode cation signed stoichiometry,-1,,
diff --git a/input/parameters/lithium-ion/separators/separator_Marquis2019/parameters.csv b/input/parameters/lithium-ion/separators/separator_Marquis2019/parameters.csv
index c07313986..d1af7d5bd 100644
--- a/input/parameters/lithium-ion/separators/separator_Marquis2019/parameters.csv
+++ b/input/parameters/lithium-ion/separators/separator_Marquis2019/parameters.csv
@@ -2,7 +2,8 @@ Name [units],Value,Reference,Notes
 # Empty rows and rows starting with ‘#’ will be ignored,,,
 ,,,
 Separator porosity,1,Scott Moura FastDFN,
-Separator Bruggeman coefficient,1.5,Scott Moura FastDFN,
+Separator Bruggeman coefficient (electrolyte),1.5,Scott Moura FastDFN,
+Separator Bruggeman coefficient (electrode),1.5,Scott Moura FastDFN,
 Separator density [kg.m-3],397,,
 Separator specific heat capacity [J.kg-1.K-1],700,,
 Separator thermal conductivity [W.m-1.K-1],0.16,,
diff --git a/pybamm/__init__.py b/pybamm/__init__.py
index ced6de35e..4e3d07864 100644
--- a/pybamm/__init__.py
+++ b/pybamm/__init__.py
@@ -164,6 +164,7 @@ from .models.submodels import (
     particle,
     porosity,
     thermal,
+    tortuosity,
 )
 
 #
diff --git a/pybamm/models/full_battery_models/base_battery_model.py b/pybamm/models/full_battery_models/base_battery_model.py
index a50592bca..90a381447 100644
--- a/pybamm/models/full_battery_models/base_battery_model.py
+++ b/pybamm/models/full_battery_models/base_battery_model.py
@@ -415,7 +415,9 @@ class BaseBatteryModel(pybamm.BaseModel):
         # 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 @@ class BaseBatteryModel(pybamm.BaseModel):
                         )
                         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 @@ class BaseBatteryModel(pybamm.BaseModel):
 
         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":
diff --git a/pybamm/models/full_battery_models/lead_acid/full.py b/pybamm/models/full_battery_models/lead_acid/full.py
index 056a3efc3..a93339cb4 100644
--- a/pybamm/models/full_battery_models/lead_acid/full.py
+++ b/pybamm/models/full_battery_models/lead_acid/full.py
@@ -37,6 +37,7 @@ class Full(BaseModel):
         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()
diff --git a/pybamm/models/full_battery_models/lead_acid/higher_order.py b/pybamm/models/full_battery_models/lead_acid/higher_order.py
index a03afdc51..8151cfe03 100644
--- a/pybamm/models/full_battery_models/lead_acid/higher_order.py
+++ b/pybamm/models/full_battery_models/lead_acid/higher_order.py
@@ -49,6 +49,7 @@ class BaseHigherOrderModel(BaseModel):
         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()
 
diff --git a/pybamm/models/full_battery_models/lead_acid/loqs.py b/pybamm/models/full_battery_models/lead_acid/loqs.py
index cf6bea3ed..14f6bf591 100644
--- a/pybamm/models/full_battery_models/lead_acid/loqs.py
+++ b/pybamm/models/full_battery_models/lead_acid/loqs.py
@@ -37,6 +37,7 @@ class LOQS(BaseModel):
         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 @@ class LOQS(BaseModel):
             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:
diff --git a/pybamm/models/full_battery_models/lithium_ion/dfn.py b/pybamm/models/full_battery_models/lithium_ion/dfn.py
index 626822962..b3acfb85a 100644
--- a/pybamm/models/full_battery_models/lithium_ion/dfn.py
+++ b/pybamm/models/full_battery_models/lithium_ion/dfn.py
@@ -35,6 +35,7 @@ class DFN(BaseModel):
 
         self.set_reactions()
         self.set_porosity_submodel()
+        self.set_tortuosity_submodels()
         self.set_convection_submodel()
         self.set_interfacial_submodel()
         self.set_particle_submodel()
diff --git a/pybamm/models/full_battery_models/lithium_ion/spm.py b/pybamm/models/full_battery_models/lithium_ion/spm.py
index d11f38abe..fb1e40538 100644
--- a/pybamm/models/full_battery_models/lithium_ion/spm.py
+++ b/pybamm/models/full_battery_models/lithium_ion/spm.py
@@ -33,6 +33,7 @@ class SPM(BaseModel):
         super().__init__(options, name)
 
         self.set_porosity_submodel()
+        self.set_tortuosity_submodels()
         self.set_convection_submodel()
         self.set_interfacial_submodel()
         self.set_particle_submodel()
diff --git a/pybamm/models/full_battery_models/lithium_ion/spme.py b/pybamm/models/full_battery_models/lithium_ion/spme.py
index 974b60759..13084f91a 100644
--- a/pybamm/models/full_battery_models/lithium_ion/spme.py
+++ b/pybamm/models/full_battery_models/lithium_ion/spme.py
@@ -37,6 +37,7 @@ class SPMe(BaseModel):
 
         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 @@ class SPMe(BaseModel):
 
         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)
diff --git a/pybamm/models/submodels/base_submodel.py b/pybamm/models/submodels/base_submodel.py
index 26f41a2af..e91d5d676 100644
--- a/pybamm/models/submodels/base_submodel.py
+++ b/pybamm/models/submodels/base_submodel.py
@@ -68,13 +68,25 @@ class BaseSubModel:
 
     @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):
diff --git a/pybamm/models/submodels/electrode/ohm/base_ohm.py b/pybamm/models/submodels/electrode/ohm/base_ohm.py
index b2bb673e5..f9174f894 100644
--- a/pybamm/models/submodels/electrode/ohm/base_ohm.py
+++ b/pybamm/models/submodels/electrode/ohm/base_ohm.py
@@ -25,21 +25,19 @@ class BaseModel(BaseElectrode):
 
     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}
diff --git a/pybamm/models/submodels/electrode/ohm/composite_ohm.py b/pybamm/models/submodels/electrode/ohm/composite_ohm.py
index 696179856..201cbcc8d 100644
--- a/pybamm/models/submodels/electrode/ohm/composite_ohm.py
+++ b/pybamm/models/submodels/electrode/ohm/composite_ohm.py
@@ -34,13 +34,13 @@ class Composite(BaseModel):
         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 @@ class Composite(BaseModel):
             ]
             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 @@ class Composite(BaseModel):
     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 @@ class Composite(BaseModel):
 
         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}
diff --git a/pybamm/models/submodels/electrode/ohm/full_ohm.py b/pybamm/models/submodels/electrode/ohm/full_ohm.py
index fac16021b..c2e64a00f 100644
--- a/pybamm/models/submodels/electrode/ohm/full_ohm.py
+++ b/pybamm/models/submodels/electrode/ohm/full_ohm.py
@@ -36,16 +36,14 @@ class Full(BaseModel):
     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 @@ class Full(BaseModel):
 
         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 @@ class Full(BaseModel):
 
         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",
diff --git a/pybamm/models/submodels/electrode/ohm/surface_form_ohm.py b/pybamm/models/submodels/electrode/ohm/surface_form_ohm.py
index f78453e49..143a1ee87 100644
--- a/pybamm/models/submodels/electrode/ohm/surface_form_ohm.py
+++ b/pybamm/models/submodels/electrode/ohm/surface_form_ohm.py
@@ -31,13 +31,13 @@ class SurfaceForm(BaseModel):
         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 @@ class SurfaceForm(BaseModel):
             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(
diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/base_higher_order_stefan_maxwell_conductivity.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/base_higher_order_stefan_maxwell_conductivity.py
index b6c5bf756..88f97a86a 100644
--- a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/base_higher_order_stefan_maxwell_conductivity.py
+++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/base_higher_order_stefan_maxwell_conductivity.py
@@ -44,9 +44,9 @@ class BaseHigherOrder(BaseModel):
             "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 @@ class BaseHigherOrder(BaseModel):
         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"]:
diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/full_stefan_maxwell_conductivity.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/full_stefan_maxwell_conductivity.py
index 89272892c..2f2857451 100644
--- a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/full_stefan_maxwell_conductivity.py
+++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/full_stefan_maxwell_conductivity.py
@@ -34,11 +34,11 @@ class Full(BaseModel):
     def get_coupled_variables(self, variables):
         param = self.param
         T = variables["Cell temperature"]
-        eps = variables["Porosity"]
+        tor = variables["Electrolyte tortuosity"]
         c_e = variables["Electrolyte concentration"]
         phi_e = variables["Electrolyte potential"]
 
-        i_e = (param.kappa_e(c_e, T) * (eps ** param.b) * param.gamma_e / param.C_e) * (
+        i_e = (param.kappa_e(c_e, T) * tor * param.gamma_e / param.C_e) * (
             param.chi(c_e) * (1 + param.Theta * T) * pybamm.grad(c_e) / c_e
             - pybamm.grad(phi_e)
         )
diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/full_surface_form_stefan_maxwell_conductivity.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/full_surface_form_stefan_maxwell_conductivity.py
index 996f9fdfe..4cde01bab 100644
--- a/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/full_surface_form_stefan_maxwell_conductivity.py
+++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/conductivity/surface_potential_form/full_surface_form_stefan_maxwell_conductivity.py
@@ -130,21 +130,18 @@ class BaseModel(BaseStefanMaxwellConductivity):
 
     def _get_conductivities(self, variables):
         param = self.param
-        eps = variables[self.domain + " electrode porosity"]
+        tor_e = variables[self.domain + " electrolyte tortuosity"]
+        tor_s = variables[self.domain + " electrode tortuosity"]
         c_e = variables[self.domain + " electrolyte concentration"]
         T = variables[self.domain + " electrode temperature"]
         if self.domain == "Negative":
             sigma = param.sigma_n
-            b = param.b_n
         elif self.domain == "Positive":
             sigma = param.sigma_p
-            b = param.b_p
-        sigma_eff = sigma * (1 - eps) ** b
-        conductivity = (
-            param.kappa_e(c_e, T)
-            * (eps ** b)
-            / (param.C_e / param.gamma_e + param.kappa_e(c_e, T) / sigma_eff)
-        )
+
+        kappa_eff = param.kappa_e(c_e, T) * tor_e
+        sigma_eff = sigma * tor_s
+        conductivity = kappa_eff / (param.C_e / param.gamma_e + kappa_eff / sigma_eff)
 
         return conductivity, sigma_eff
 
@@ -190,11 +187,11 @@ class BaseModel(BaseStefanMaxwellConductivity):
         i_boundary_cc = variables["Current collector current density"]
         c_e_s = variables["Separator electrolyte concentration"]
         phi_e_n = variables["Negative electrolyte potential"]
-        eps_s = variables["Separator porosity"]
+        tor_s = variables["Separator porosity"]
         T = variables["Separator temperature"]
 
         chi_e_s = param.chi(c_e_s)
-        kappa_s_eff = param.kappa_e(c_e_s, T) * (eps_s ** param.b_s)
+        kappa_s_eff = param.kappa_e(c_e_s, T) * tor_s
 
         phi_e_s = pybamm.PrimaryBroadcast(
             pybamm.boundary_value(phi_e_n, "right"), "separator"
@@ -226,14 +223,14 @@ class BaseModel(BaseStefanMaxwellConductivity):
 
         x_n = pybamm.standard_spatial_vars.x_n
         x_p = pybamm.standard_spatial_vars.x_p
-        eps = variables[self.domain + " electrode porosity"]
+        tor = variables[self.domain + " electrode tortuosity"]
         i_boundary_cc = variables["Current collector current density"]
         i_e = variables[self.domain + " electrolyte current density"]
 
         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.IndefiniteIntegral(i_s / conductivity, x_n)
 
         elif self.domain == "Positive":
@@ -241,7 +238,7 @@ class BaseModel(BaseStefanMaxwellConductivity):
             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
diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/composite_stefan_maxwell_diffusion.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/composite_stefan_maxwell_diffusion.py
index dc7cd379d..574a65501 100644
--- a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/composite_stefan_maxwell_diffusion.py
+++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/composite_stefan_maxwell_diffusion.py
@@ -29,7 +29,7 @@ class Composite(Full):
 
     def get_coupled_variables(self, variables):
 
-        eps_0 = variables["Leading-order porosity"]
+        tor_0 = variables["Leading-order electrolyte tortuosity"]
         c_e_0_av = variables["Leading-order x-averaged electrolyte concentration"]
         c_e = variables["Electrolyte concentration"]
         # i_e = variables["Electrolyte current density"]
@@ -40,7 +40,7 @@ class Composite(Full):
 
         whole_cell = ["negative electrode", "separator", "positive electrode"]
         N_e_diffusion = (
-            -(eps_0 ** param.b)
+            -tor_0
             * pybamm.PrimaryBroadcast(param.D_e(c_e_0_av, T_0), whole_cell)
             * pybamm.grad(c_e)
         )
diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/first_order_stefan_maxwell_diffusion.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/first_order_stefan_maxwell_diffusion.py
index f0914b417..a4c773cd0 100644
--- a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/first_order_stefan_maxwell_diffusion.py
+++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/first_order_stefan_maxwell_diffusion.py
@@ -40,6 +40,9 @@ class FirstOrder(BaseModel):
         eps_n_0 = variables["Leading-order x-averaged negative electrode porosity"]
         eps_s_0 = variables["Leading-order x-averaged separator porosity"]
         eps_p_0 = variables["Leading-order x-averaged positive electrode porosity"]
+        tor_n_0 = variables["Leading-order x-averaged negative electrolyte tortuosity"]
+        tor_s_0 = variables["Leading-order x-averaged separator tortuosity"]
+        tor_p_0 = variables["Leading-order x-averaged positive electrolyte tortuosity"]
         deps_n_0_dt = variables[
             "Leading-order x-averaged negative electrode porosity change"
         ]
@@ -70,9 +73,9 @@ class FirstOrder(BaseModel):
         )
 
         # Diffusivities
-        D_e_n = (eps_n_0 ** param.b_n) * param.D_e(c_e_0, T_0)
-        D_e_s = (eps_s_0 ** param.b_s) * param.D_e(c_e_0, T_0)
-        D_e_p = (eps_p_0 ** param.b_p) * param.D_e(c_e_0, T_0)
+        D_e_n = tor_n_0 * param.D_e(c_e_0, T_0)
+        D_e_s = tor_s_0 * param.D_e(c_e_0, T_0)
+        D_e_p = tor_p_0 * param.D_e(c_e_0, T_0)
 
         # Fluxes
         N_e_n_1 = -pybamm.outer(rhs_n, x_n)
diff --git a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/full_stefan_maxwell_diffusion.py b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/full_stefan_maxwell_diffusion.py
index 907627cfd..e007a5e6f 100644
--- a/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/full_stefan_maxwell_diffusion.py
+++ b/pybamm/models/submodels/electrolyte/stefan_maxwell/diffusion/full_stefan_maxwell_diffusion.py
@@ -31,7 +31,7 @@ class Full(BaseModel):
 
     def get_coupled_variables(self, variables):
 
-        eps = variables["Porosity"]
+        tor = variables["Electrolyte tortuosity"]
         c_e = variables["Electrolyte concentration"]
         # i_e = variables["Electrolyte current density"]
         v_box = variables["Volume-averaged velocity"]
@@ -39,7 +39,7 @@ class Full(BaseModel):
 
         param = self.param
 
-        N_e_diffusion = -(eps ** param.b) * param.D_e(c_e, T) * pybamm.grad(c_e)
+        N_e_diffusion = -tor * param.D_e(c_e, T) * pybamm.grad(c_e)
         # N_e_migration = (param.C_e * param.t_plus) / param.gamma_e * i_e
         # N_e_convection = c_e * v_box
 
diff --git a/pybamm/models/submodels/interface/diffusion_limited/full_diffusion_limited.py b/pybamm/models/submodels/interface/diffusion_limited/full_diffusion_limited.py
index e085ded46..a9f4bfc9f 100644
--- a/pybamm/models/submodels/interface/diffusion_limited/full_diffusion_limited.py
+++ b/pybamm/models/submodels/interface/diffusion_limited/full_diffusion_limited.py
@@ -26,10 +26,10 @@ class FullDiffusionLimited(BaseModel):
     def _get_diffusion_limited_current_density(self, variables):
         param = self.param
         if self.domain == "Negative":
-            eps_s = variables["Separator porosity"]
+            tor_s = variables["Separator tortuosity"]
             c_ox_s = variables["Separator oxygen concentration"]
             N_ox_neg_sep_interface = (
-                -pybamm.boundary_value(eps_s ** param.b_s, "left")
+                -pybamm.boundary_value(tor_s, "left")
                 * param.curlyD_ox
                 * pybamm.BoundaryGradient(c_ox_s, "left")
             )
diff --git a/pybamm/models/submodels/oxygen_diffusion/composite_oxygen_diffusion.py b/pybamm/models/submodels/oxygen_diffusion/composite_oxygen_diffusion.py
index a9f1c8411..41c1f54d2 100644
--- a/pybamm/models/submodels/oxygen_diffusion/composite_oxygen_diffusion.py
+++ b/pybamm/models/submodels/oxygen_diffusion/composite_oxygen_diffusion.py
@@ -32,21 +32,16 @@ class Composite(Full):
 
     def get_coupled_variables(self, variables):
 
-        eps_0 = separator_and_positive_only(variables["Leading-order porosity"])
+        tor_0 = separator_and_positive_only(
+            variables["Leading-order electrolyte tortuosity"]
+        )
         c_ox = variables["Separator and positive electrode oxygen concentration"]
         # TODO: allow charge and convection?
         # v_box = pybamm.Scalar(0)
 
         param = self.param
 
-        b = pybamm.Concatenation(
-            pybamm.FullBroadcast(param.b_s, ["separator"], "current collector"),
-            pybamm.FullBroadcast(
-                param.b_p, ["positive electrode"], "current collector"
-            ),
-        )
-
-        N_ox_diffusion = -(eps_0 ** b) * param.curlyD_ox * pybamm.grad(c_ox)
+        N_ox_diffusion = -tor_0 * param.curlyD_ox * pybamm.grad(c_ox)
 
         N_ox = N_ox_diffusion  # + c_ox * v_box
         # Flux in the negative electrode is zero
diff --git a/pybamm/models/submodels/oxygen_diffusion/first_order_oxygen_diffusion.py b/pybamm/models/submodels/oxygen_diffusion/first_order_oxygen_diffusion.py
index f3d5d1799..9422c2d1d 100644
--- a/pybamm/models/submodels/oxygen_diffusion/first_order_oxygen_diffusion.py
+++ b/pybamm/models/submodels/oxygen_diffusion/first_order_oxygen_diffusion.py
@@ -37,12 +37,14 @@ class FirstOrder(BaseModel):
         x_p = pybamm.standard_spatial_vars.x_p
 
         # Unpack
-        eps_s_0_av = variables["Leading-order x-averaged separator porosity"]
-        eps_p_0_av = variables["Leading-order x-averaged positive electrode porosity"]
+        tor_s_0_av = variables["Leading-order x-averaged separator tortuosity"]
+        tor_p_0_av = variables[
+            "Leading-order x-averaged positive electrolyte tortuosity"
+        ]
 
         # Diffusivities
-        D_ox_s = (eps_s_0_av ** param.b_s) * param.curlyD_ox
-        D_ox_p = (eps_p_0_av ** param.b_p) * param.curlyD_ox
+        D_ox_s = tor_s_0_av * param.curlyD_ox
+        D_ox_p = tor_p_0_av * param.curlyD_ox
 
         # Reactions
         sj_ox_p = sum(
diff --git a/pybamm/models/submodels/oxygen_diffusion/full_oxygen_diffusion.py b/pybamm/models/submodels/oxygen_diffusion/full_oxygen_diffusion.py
index 77c1dd908..515153a5c 100644
--- a/pybamm/models/submodels/oxygen_diffusion/full_oxygen_diffusion.py
+++ b/pybamm/models/submodels/oxygen_diffusion/full_oxygen_diffusion.py
@@ -67,21 +67,14 @@ class Full(BaseModel):
 
     def get_coupled_variables(self, variables):
 
-        eps = separator_and_positive_only(variables["Porosity"])
+        tor = separator_and_positive_only(variables["Electrolyte tortuosity"])
         c_ox = variables["Separator and positive electrode oxygen concentration"]
         # TODO: allow charge and convection?
         v_box = pybamm.Scalar(0)
 
         param = self.param
 
-        b = pybamm.Concatenation(
-            pybamm.FullBroadcast(param.b_s, ["separator"], "current collector"),
-            pybamm.FullBroadcast(
-                param.b_p, ["positive electrode"], "current collector"
-            ),
-        )
-
-        N_ox_diffusion = -(eps ** b) * param.curlyD_ox * pybamm.grad(c_ox)
+        N_ox_diffusion = -tor * param.curlyD_ox * pybamm.grad(c_ox)
 
         N_ox = N_ox_diffusion + c_ox * v_box
         # Flux in the negative electrode is zero
diff --git a/pybamm/models/submodels/porosity/base_porosity.py b/pybamm/models/submodels/porosity/base_porosity.py
index 158adb1d7..72012656a 100644
--- a/pybamm/models/submodels/porosity/base_porosity.py
+++ b/pybamm/models/submodels/porosity/base_porosity.py
@@ -33,21 +33,30 @@ class BaseModel(pybamm.BaseSubModel):
             "X-averaged positive electrode porosity": pybamm.x_average(eps_p),
         }
 
+        # activate material volume fractions
+        eps_solid_n = 1 - eps_n - self.param.epsilon_inactive_n
+        eps_solid_s = 1 - eps_s - self.param.epsilon_inactive_s
+        eps_solid_p = 1 - eps_p - self.param.epsilon_inactive_p
+        eps_solid = pybamm.Concatenation(eps_solid_n, eps_solid_s, eps_solid_p)
+
+        am = "active material volume fraction"
+        variables.update(
+            {
+                am.capitalize(): eps_solid,
+                "Negative electrode " + am: eps_solid_n,
+                "Separator " + am: eps_solid,
+                "Positive electrode " + am: eps_solid_p,
+                "X-averaged negative electrode " + am: pybamm.x_average(eps_solid_n),
+                "X-averaged separator " + am: pybamm.x_average(eps_solid),
+                "X-averaged positive electrode " + am: pybamm.x_average(eps_solid_p),
+            }
+        )
+
         if set_leading_order is True:
-            variables.update(
-                {
-                    "Leading-order negative electrode porosity": eps_n,
-                    "Leading-order separator porosity": eps_s,
-                    "Leading-order positive electrode porosity": eps_p,
-                    "Leading-order x-averaged "
-                    + "negative electrode porosity": pybamm.x_average(eps_n),
-                    "Leading-order x-averaged separator porosity": pybamm.x_average(
-                        eps_s
-                    ),
-                    "Leading-order x-averaged "
-                    + "positive electrode porosity": pybamm.x_average(eps_p),
-                }
-            )
+            leading_order_variables = {
+                "Leading-order " + name.lower(): var for name, var in variables.items()
+            }
+            variables.update(leading_order_variables)
 
         return variables
 
@@ -83,3 +92,12 @@ class BaseModel(pybamm.BaseSubModel):
             )
 
         return variables
+
+    def set_events(self, variables):
+        eps_n = variables["Negative electrode porosity"]
+        eps_p = variables["Positive electrode porosity"]
+        self.events["Zero negative electrode porosity cut-off"] = pybamm.min(eps_n)
+        self.events["Max negative electrode porosity cut-off"] = pybamm.max(eps_n) - 1
+        self.events["Zero positive electrode porosity cut-off"] = pybamm.min(eps_p)
+        self.events["Max positive electrode porosity cut-off"] = pybamm.max(eps_p) - 1
+
diff --git a/pybamm/models/submodels/porosity/full_reaction_driven_porosity.py b/pybamm/models/submodels/porosity/full_reaction_driven_porosity.py
index 3f65460ce..7f911abed 100644
--- a/pybamm/models/submodels/porosity/full_reaction_driven_porosity.py
+++ b/pybamm/models/submodels/porosity/full_reaction_driven_porosity.py
@@ -52,4 +52,4 @@ class Full(BaseModel):
 
     def set_initial_conditions(self, variables):
         eps = variables["Porosity"]
-        self.initial_conditions = {eps: self.param.eps_init}
+        self.initial_conditions = {eps: self.param.epsilon_init}
diff --git a/pybamm/models/submodels/porosity/leading_reaction_driven_porosity.py b/pybamm/models/submodels/porosity/leading_reaction_driven_porosity.py
index 137746edf..e2eb57e69 100644
--- a/pybamm/models/submodels/porosity/leading_reaction_driven_porosity.py
+++ b/pybamm/models/submodels/porosity/leading_reaction_driven_porosity.py
@@ -60,6 +60,6 @@ class LeadingOrder(BaseModel):
         eps_s_av = variables["X-averaged separator porosity"]
         eps_p_av = variables["X-averaged positive electrode porosity"]
 
-        self.initial_conditions = {eps_n_av: self.param.eps_n_init}
-        self.initial_conditions.update({eps_s_av: self.param.eps_s_init})
-        self.initial_conditions.update({eps_p_av: self.param.eps_p_init})
+        self.initial_conditions = {eps_n_av: self.param.epsilon_n_init}
+        self.initial_conditions.update({eps_s_av: self.param.epsilon_s_init})
+        self.initial_conditions.update({eps_p_av: self.param.epsilon_p_init})
diff --git a/pybamm/models/submodels/tortuosity/__init__.py b/pybamm/models/submodels/tortuosity/__init__.py
new file mode 100644
index 000000000..cf0411c9d
--- /dev/null
+++ b/pybamm/models/submodels/tortuosity/__init__.py
@@ -0,0 +1,2 @@
+from .base_tortuosity import BaseModel
+from .bruggeman_tortuosity import Bruggeman
diff --git a/pybamm/models/submodels/tortuosity/base_tortuosity.py b/pybamm/models/submodels/tortuosity/base_tortuosity.py
new file mode 100644
index 000000000..7c67705d7
--- /dev/null
+++ b/pybamm/models/submodels/tortuosity/base_tortuosity.py
@@ -0,0 +1,52 @@
+#
+# Base class for tortuosity
+#
+import pybamm
+
+
+class BaseModel(pybamm.BaseSubModel):
+    """Base class for tortuosity
+
+    Parameters
+    ----------
+    param : parameter class
+        The parameters to use for this submodel
+    phase : str
+        The material for the model ('electrolyte' or 'electrode').
+
+    **Extends:** :class:`pybamm.BaseSubModel`
+    """
+
+    def __init__(self, param, phase):
+        super().__init__(param)
+        self.phase = phase
+
+    def _get_standard_tortuosity_variables(self, tor, set_leading_order=False):
+        tor_n, tor_s, tor_p = tor.orphans
+
+        variables = {
+            self.phase + " tortuosity": tor,
+            "Negative " + self.phase.lower() + " tortuosity": tor_n,
+            "Positive " + self.phase.lower() + " tortuosity": tor_p,
+            "X-averaged negative "
+            + self.phase.lower()
+            + " tortuosity": pybamm.x_average(tor_n),
+            "X-averaged positive "
+            + self.phase.lower()
+            + " tortuosity": pybamm.x_average(tor_p),
+        }
+        if self.phase == "Electrolyte":
+            variables.update(
+                {
+                    "Separator tortuosity": tor_s,
+                    "X-averaged separator tortuosity": pybamm.x_average(tor_s),
+                }
+            )
+
+        if set_leading_order is True:
+            leading_order_variables = {
+                "Leading-order " + name.lower(): var for name, var in variables.items()
+            }
+            variables.update(leading_order_variables)
+
+        return variables
diff --git a/pybamm/models/submodels/tortuosity/bruggeman_tortuosity.py b/pybamm/models/submodels/tortuosity/bruggeman_tortuosity.py
new file mode 100644
index 000000000..0ccc5a8a4
--- /dev/null
+++ b/pybamm/models/submodels/tortuosity/bruggeman_tortuosity.py
@@ -0,0 +1,37 @@
+#
+# Class for Bruggeman tortuosity
+#
+import pybamm
+
+from .base_tortuosity import BaseModel
+
+
+class Bruggeman(BaseModel):
+    """Submodel for Bruggeman tortuosity
+
+    **Extends:** :class:`pybamm.tortuosity.BaseModel`
+    """
+
+    def __init__(self, param, phase, set_leading_order=False):
+        super().__init__(param, phase)
+        self.set_leading_order = set_leading_order
+
+    def get_coupled_variables(self, variables):
+        param = self.param
+
+        if self.phase == "Electrolyte":
+            eps_n, eps_s, eps_p = variables["Porosity"].orphans
+            tor = pybamm.Concatenation(
+                eps_n ** param.b_e_n, eps_s ** param.b_e_s, eps_p ** param.b_e_p
+            )
+        elif self.phase == "Electrode":
+            eps_n, eps_s, eps_p = variables["Active material volume fraction"].orphans
+            tor = pybamm.Concatenation(
+                eps_n ** param.b_s_n, eps_s ** param.b_s_s, eps_p ** param.b_s_p
+            )
+
+        variables.update(
+            self._get_standard_tortuosity_variables(tor, self.set_leading_order)
+        )
+
+        return variables
diff --git a/pybamm/parameters/geometric_parameters.py b/pybamm/parameters/geometric_parameters.py
index 3b70db486..d10d4a026 100644
--- a/pybamm/parameters/geometric_parameters.py
+++ b/pybamm/parameters/geometric_parameters.py
@@ -36,14 +36,12 @@ a_n_dim = pybamm.Parameter("Negative electrode surface area density [m-1]")
 a_p_dim = pybamm.Parameter("Positive electrode surface area density [m-1]")
 R_n = pybamm.Parameter("Negative particle radius [m]")
 R_p = pybamm.Parameter("Positive particle radius [m]")
-b_n = pybamm.Parameter("Negative electrode Bruggeman coefficient")
-b_s = pybamm.Parameter("Separator Bruggeman coefficient")
-b_p = pybamm.Parameter("Positive electrode Bruggeman coefficient")
-b = pybamm.Concatenation(
-    pybamm.FullBroadcast(b_n, ["negative electrode"], "current collector"),
-    pybamm.FullBroadcast(b_s, ["separator"], "current collector"),
-    pybamm.FullBroadcast(b_p, ["positive electrode"], "current collector"),
-)
+b_e_n = pybamm.Parameter("Negative electrode Bruggeman coefficient (electrolyte)")
+b_e_s = pybamm.Parameter("Separator Bruggeman coefficient (electrolyte)")
+b_e_p = pybamm.Parameter("Positive electrode Bruggeman coefficient (electrolyte)")
+b_s_n = pybamm.Parameter("Negative electrode Bruggeman coefficient (electrode)")
+b_s_s = pybamm.Parameter("Separator Bruggeman coefficient (electrode)")
+b_s_p = pybamm.Parameter("Positive electrode Bruggeman coefficient (electrode)")
 
 # --------------------------------------------------------------------------------------
 "Dimensionless Parameters"
diff --git a/pybamm/parameters/standard_parameters_lead_acid.py b/pybamm/parameters/standard_parameters_lead_acid.py
index 496d5fb85..ef92f84dc 100644
--- a/pybamm/parameters/standard_parameters_lead_acid.py
+++ b/pybamm/parameters/standard_parameters_lead_acid.py
@@ -72,12 +72,18 @@ sigma_cp_dimensional = sigma_p_dim
 # Microstructure
 a_n_dim = pybamm.geometric_parameters.a_n_dim
 a_p_dim = pybamm.geometric_parameters.a_p_dim
-b_n = pybamm.geometric_parameters.b_n
-b_s = pybamm.geometric_parameters.b_s
-b_p = pybamm.geometric_parameters.b_p
-b = pybamm.geometric_parameters.b
+b_e_n = pybamm.geometric_parameters.b_e_n
+b_e_s = pybamm.geometric_parameters.b_e_s
+b_e_p = pybamm.geometric_parameters.b_e_p
+b_s_n = pybamm.geometric_parameters.b_s_n
+b_s_s = pybamm.geometric_parameters.b_s_s
+b_s_p = pybamm.geometric_parameters.b_s_p
 xi_n = pybamm.Parameter("Negative electrode morphological parameter")
 xi_p = pybamm.Parameter("Positive electrode morphological parameter")
+# no binder
+epsilon_inactive_n = pybamm.Scalar(0)
+epsilon_inactive_s = pybamm.Scalar(0)
+epsilon_inactive_p = pybamm.Scalar(0)
 
 # Electrochemical reactions
 # Main
@@ -393,19 +399,19 @@ voltage_high_cut = (
 # Electrolyte volumetric capacity
 Q_e_max = (l_n * eps_n_max + l_s * eps_s_max + l_p * eps_p_max) / (s_p - s_n)
 Q_e_max_dimensional = Q_e_max * c_e_typ * F
-capacity = Q_e_max_dimensional * 8 * A_cs * L_x
+capacity = Q_e_max_dimensional * n_electrodes_parallel * A_cs * L_x
 
 # Initial conditions
 q_init = pybamm.Parameter("Initial State of Charge")
 c_e_init = q_init
 c_ox_init = c_ox_init_dim / c_ox_typ
-eps_n_init = eps_n_max - beta_surf_n * Q_e_max / l_n * (1 - q_init)
-eps_s_init = eps_s_max
-eps_p_init = eps_p_max + beta_surf_p * Q_e_max / l_p * (1 - q_init)
-eps_init = pybamm.Concatenation(
-    pybamm.FullBroadcast(eps_n_init, ["negative electrode"], "current collector"),
-    pybamm.FullBroadcast(eps_s_init, ["separator"], "current collector"),
-    pybamm.FullBroadcast(eps_p_init, ["positive electrode"], "current collector"),
+epsilon_n_init = eps_n_max - beta_surf_n * Q_e_max / l_n * (1 - q_init)
+epsilon_s_init = eps_s_max
+epsilon_p_init = eps_p_max + beta_surf_p * Q_e_max / l_p * (1 - q_init)
+epsilon_init = pybamm.Concatenation(
+    pybamm.FullBroadcast(epsilon_n_init, ["negative electrode"], "current collector"),
+    pybamm.FullBroadcast(epsilon_s_init, ["separator"], "current collector"),
+    pybamm.FullBroadcast(epsilon_p_init, ["positive electrode"], "current collector"),
 )
 curlyU_n_init = Q_e_max * (1.2 - q_init) / (Q_n_max * l_n)
 curlyU_p_init = Q_e_max * (1.2 - q_init) / (Q_p_max * l_p)
diff --git a/pybamm/parameters/standard_parameters_lithium_ion.py b/pybamm/parameters/standard_parameters_lithium_ion.py
index 73a9f4b79..62fbcf646 100644
--- a/pybamm/parameters/standard_parameters_lithium_ion.py
+++ b/pybamm/parameters/standard_parameters_lithium_ion.py
@@ -82,10 +82,12 @@ a_k_dim = pybamm.Concatenation(
 )
 R_n = pybamm.geometric_parameters.R_n
 R_p = pybamm.geometric_parameters.R_p
-b_n = pybamm.geometric_parameters.b_n
-b_s = pybamm.geometric_parameters.b_s
-b_p = pybamm.geometric_parameters.b_p
-b = pybamm.geometric_parameters.b
+b_e_n = pybamm.geometric_parameters.b_e_n
+b_e_s = pybamm.geometric_parameters.b_e_s
+b_e_p = pybamm.geometric_parameters.b_e_p
+b_s_n = pybamm.geometric_parameters.b_s_n
+b_s_s = pybamm.geometric_parameters.b_s_s
+b_s_p = pybamm.geometric_parameters.b_s_p
 
 # Electrochemical reactions
 ne_n = pybamm.Parameter("Negative electrode electrons in reaction")
@@ -277,6 +279,11 @@ epsilon = pybamm.Concatenation(
     pybamm.FullBroadcast(epsilon_s, ["separator"], "current collector"),
     pybamm.FullBroadcast(epsilon_p, ["positive electrode"], "current collector"),
 )
+epsilon_s_n = pybamm.Parameter("Negative electrode active material volume fraction")
+epsilon_s_p = pybamm.Parameter("Positive electrode active material volume fraction")
+epsilon_inactive_n = 1 - epsilon_n - epsilon_s_n
+epsilon_inactive_s = 1 - epsilon_s
+epsilon_inactive_p = 1 - epsilon_p - epsilon_s_p
 a_n = a_n_dim * R_n
 a_p = a_p_dim * R_p
 


Test patch:
diff --git a/tests/unit/test_expression_tree/test_symbol.py b/tests/unit/test_expression_tree/test_symbol.py
index 840185b10..55956c551 100644
--- a/tests/unit/test_expression_tree/test_symbol.py
+++ b/tests/unit/test_expression_tree/test_symbol.py
@@ -272,6 +272,7 @@ class TestSymbol(unittest.TestCase):
 
         variables = {
             "Porosity": param.epsilon,
+            "Electrolyte tortuosity": param.epsilon ** 1.5,
             "Porosity change": deps_dt,
             "Volume-averaged velocity": v_box,
             "Negative electrode interfacial current density": one_n,
diff --git a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_composite_stefan_maxwell_conductivity.py b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_composite_stefan_maxwell_conductivity.py
index c3e55d659..bfe620ec2 100644
--- a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_composite_stefan_maxwell_conductivity.py
+++ b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_composite_stefan_maxwell_conductivity.py
@@ -22,6 +22,9 @@ class TestComposite(unittest.TestCase):
             "Leading-order x-averaged negative electrode porosity": a,
             "Leading-order x-averaged separator porosity": a,
             "Leading-order x-averaged positive electrode porosity": a,
+            "Leading-order x-averaged negative electrolyte tortuosity": a,
+            "Leading-order x-averaged separator tortuosity": a,
+            "Leading-order x-averaged positive electrolyte tortuosity": a,
             "X-averaged cell temperature": a,
         }
         submodel = pybamm.electrolyte.stefan_maxwell.conductivity.Composite(param)
diff --git a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_full_stefan_maxwell_conductivity.py b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_full_stefan_maxwell_conductivity.py
index 2aecbf88a..80870cbc6 100644
--- a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_full_stefan_maxwell_conductivity.py
+++ b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_full_stefan_maxwell_conductivity.py
@@ -12,7 +12,7 @@ class TestFull(unittest.TestCase):
         param = pybamm.standard_parameters_lithium_ion
         a = pybamm.Scalar(0)
         variables = {
-            "Porosity": a,
+            "Electrolyte tortuosity": a,
             "Electrolyte concentration": a,
             "Negative electrode interfacial current density": pybamm.FullBroadcast(
                 a, "negative electrode", "current collector"
diff --git a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_surface_form/test_full_surface_form_stefan_maxwell_conductivity.py b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_surface_form/test_full_surface_form_stefan_maxwell_conductivity.py
index bc4bc47d7..530ee7206 100644
--- a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_surface_form/test_full_surface_form_stefan_maxwell_conductivity.py
+++ b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_conductivity/test_surface_form/test_full_surface_form_stefan_maxwell_conductivity.py
@@ -23,6 +23,8 @@ class TestFull(unittest.TestCase):
                 a, "current collector"
             ),
             "Negative electrode porosity": a_n,
+            "Negative electrolyte tortuosity": a_n,
+            "Negative electrode tortuosity": a_n,
             "Negative electrolyte concentration": a_n,
             "Negative electrode interfacial current density": a_n,
             "Electrolyte potential": pybamm.Concatenation(a_n, a_s, a_p),
@@ -55,6 +57,8 @@ class TestFull(unittest.TestCase):
             "Separator electrolyte potential": a_s,
             "Separator electrolyte current density": a_s,
             "Positive electrode porosity": a_p,
+            "Positive electrolyte tortuosity": a_p,
+            "Positive electrode tortuosity": a_p,
             "Positive electrolyte concentration": a_p,
             "Positive electrode interfacial current density": a_p,
             "Positive electrode temperature": a_p,
diff --git a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_diffusion/test_full_stefan_maxwell_diffusion.py b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_diffusion/test_full_stefan_maxwell_diffusion.py
index a45ec88bb..87b8de8d9 100644
--- a/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_diffusion/test_full_stefan_maxwell_diffusion.py
+++ b/tests/unit/test_models/test_submodels/test_electrolyte/test_stefan_maxwell/test_diffusion/test_full_stefan_maxwell_diffusion.py
@@ -13,6 +13,7 @@ class TestFull(unittest.TestCase):
         a = pybamm.Scalar(0)
         variables = {
             "Porosity": a,
+            "Electrolyte tortuosity": a,
             "Porosity change": a,
             "Volume-averaged velocity": a,
             "Electrolyte concentration": a,
diff --git a/tests/unit/test_models/test_submodels/test_oxygen_diffusion/test_full_oxygen_diffusion.py b/tests/unit/test_models/test_submodels/test_oxygen_diffusion/test_full_oxygen_diffusion.py
index 5618f0200..cf11d742a 100644
--- a/tests/unit/test_models/test_submodels/test_oxygen_diffusion/test_full_oxygen_diffusion.py
+++ b/tests/unit/test_models/test_submodels/test_oxygen_diffusion/test_full_oxygen_diffusion.py
@@ -17,6 +17,11 @@ class TestFull(unittest.TestCase):
                 pybamm.FullBroadcast(a, "separator", "current collector"),
                 pybamm.FullBroadcast(a, "positive electrode", "current collector"),
             ),
+            "Electrolyte tortuosity": pybamm.Concatenation(
+                pybamm.FullBroadcast(a, "negative electrode", "current collector"),
+                pybamm.FullBroadcast(a, "separator", "current collector"),
+                pybamm.FullBroadcast(a, "positive electrode", "current collector"),
+            ),
             "Porosity change": pybamm.Concatenation(
                 pybamm.FullBroadcast(a, "negative electrode", "current collector"),
                 pybamm.FullBroadcast(a, "separator", "current collector"),
diff --git a/tests/unit/test_models/test_submodels/test_porosity/test_base_porosity.py b/tests/unit/test_models/test_submodels/test_porosity/test_base_porosity.py
index eba6e0d35..4025bde77 100644
--- a/tests/unit/test_models/test_submodels/test_porosity/test_base_porosity.py
+++ b/tests/unit/test_models/test_submodels/test_porosity/test_base_porosity.py
@@ -9,8 +9,10 @@ import unittest
 
 class TestBaseModel(unittest.TestCase):
     def test_public_functions(self):
+        a = pybamm.Scalar(0)
+        variables = {"Negative electrode porosity": a, "Positive electrode porosity": a}
         submodel = pybamm.porosity.BaseModel(None)
-        std_tests = tests.StandardSubModelTests(submodel)
+        std_tests = tests.StandardSubModelTests(submodel, variables)
         std_tests.test_all()
 
 
diff --git a/tests/unit/test_models/test_submodels/test_tortuosity/__init__.py b/tests/unit/test_models/test_submodels/test_tortuosity/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/tests/unit/test_models/test_submodels/test_tortuosity/test_base_tortuosity.py b/tests/unit/test_models/test_submodels/test_tortuosity/test_base_tortuosity.py
new file mode 100644
index 000000000..6ab387911
--- /dev/null
+++ b/tests/unit/test_models/test_submodels/test_tortuosity/test_base_tortuosity.py
@@ -0,0 +1,28 @@
+#
+# Test base tortuosity submodel
+#
+
+import pybamm
+import tests
+import unittest
+
+
+class TestBaseModel(unittest.TestCase):
+    def test_public_functions(self):
+        submodel = pybamm.tortuosity.BaseModel(None, "Electrode")
+        std_tests = tests.StandardSubModelTests(submodel)
+        std_tests.test_all()
+
+        submodel = pybamm.tortuosity.BaseModel(None, "Electrolyte")
+        std_tests = tests.StandardSubModelTests(submodel)
+        std_tests.test_all()
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    pybamm.settings.debug_mode = True
+    unittest.main()
diff --git a/tests/unit/test_models/test_submodels/test_tortuosity/test_bruggeman_tortuosity.py b/tests/unit/test_models/test_submodels/test_tortuosity/test_bruggeman_tortuosity.py
new file mode 100644
index 000000000..eed98a060
--- /dev/null
+++ b/tests/unit/test_models/test_submodels/test_tortuosity/test_bruggeman_tortuosity.py
@@ -0,0 +1,35 @@
+#
+# Test Bruggeman tortuosity submodel
+#
+
+import pybamm
+import tests
+import unittest
+
+
+class TestBruggeman(unittest.TestCase):
+    def test_public_functions(self):
+        param = pybamm.standard_parameters_lithium_ion
+        a = pybamm.Concatenation(
+            pybamm.FullBroadcast(0, ["negative electrode"], "current collector"),
+            pybamm.FullBroadcast(0, ["separator"], "current collector"),
+            pybamm.FullBroadcast(0, ["positive electrode"], "current collector"),
+        )
+        variables = {"Porosity": a, "Active material volume fraction": a}
+        submodel = pybamm.tortuosity.Bruggeman(param, "Electrolyte")
+        std_tests = tests.StandardSubModelTests(submodel, variables)
+        std_tests.test_all()
+
+        submodel = pybamm.tortuosity.Bruggeman(param, "Electrode")
+        std_tests = tests.StandardSubModelTests(submodel, variables)
+        std_tests.test_all()
+
+
+if __name__ == "__main__":
+    print("Add -v for more debug output")
+    import sys
+
+    if "-v" in sys.argv:
+        debug = True
+    pybamm.settings.debug_mode = True
+    unittest.main()
diff --git a/tests/unit/test_parameters/test_standard_parameters_lead_acid.py b/tests/unit/test_parameters/test_standard_parameters_lead_acid.py
index e5f7a6193..0521776d5 100644
--- a/tests/unit/test_parameters/test_standard_parameters_lead_acid.py
+++ b/tests/unit/test_parameters/test_standard_parameters_lead_acid.py
@@ -157,9 +157,15 @@ class TestStandardParametersLeadAcid(unittest.TestCase):
         # Test that relevant parameters have changed as expected
         self.assertLess(param_eval_update["q_init"], param_eval["q_init"])
         self.assertLess(param_eval_update["c_e_init"], param_eval["c_e_init"])
-        self.assertLess(param_eval_update["eps_n_init"], param_eval["eps_n_init"])
-        self.assertEqual(param_eval_update["eps_s_init"], param_eval["eps_s_init"])
-        self.assertLess(param_eval_update["eps_p_init"], param_eval["eps_p_init"])
+        self.assertLess(
+            param_eval_update["epsilon_n_init"], param_eval["epsilon_n_init"]
+        )
+        self.assertEqual(
+            param_eval_update["epsilon_s_init"], param_eval["epsilon_s_init"]
+        )
+        self.assertLess(
+            param_eval_update["epsilon_p_init"], param_eval["epsilon_p_init"]
+        )
         self.assertGreater(
             param_eval_update["curlyU_n_init"], param_eval["curlyU_n_init"]
         )