Porosity change with particle size distributions

CHANGELOG.md

@@ -2,6 +2,7 @@
 
 ## Features
 
+- Porosity change now works for particle size distributions ([#5417](https://github.com/pybamm-team/PyBaMM/pull/5417))
 - Improved the performance of processed variables by replacing `casadi.vertcat` input stacking with numpy vectors. ([#5413](https://github.com/pybamm-team/PyBaMM/pull/5413))
 - Preserve custom variables and events in built-in model to_config ([#5411](https://github.com/pybamm-team/PyBaMM/pull/5411))
 - Allow out of bounds initial state of charge to enable initialising a simulation at a voltage outside the voltage limits. ([#5386](https://github.com/pybamm-team/PyBaMM/pull/5386))
@@ -18,12 +19,14 @@ as initial conditions. ([#5311](https://github.com/pybamm-team/PyBaMM/pull/5311)
 
 ## Bug fixes
 
+- SEI concentration is always a size-averaged quantity if particle size distribution present ([#5417](https://github.com/pybamm-team/PyBaMM/pull/5417))
 - Fixed a bug in the exchange current density calculation for MSMR models. ([#5404](https://github.com/pybamm-team/PyBaMM/pull/5404))
 - Fixed a bug where when converting `ExpressionFunctionParameter` to source code, `Interpolant` objects were being reduced to just their input variable names (e.g., sto) instead of preserving the full constructor call with data arrays. ([#5393](https://github.com/pybamm-team/PyBaMM/pull/5393))
 - Fixed a bug when using 2 phases with particle size distribution the bounds of the particle size distribution were always taken from the primary phase. ([#5415](https://github.com/pybamm-team/PyBaMM/pull/5415))
 
 ## Breaking changes
 
+- Porosity change no longer excludes initial SEI thickness, i.e. the initial SEI now takes up space. ([#5417](https://github.com/pybamm-team/PyBaMM/pull/5417))
 - The mass matrix inverse is no longer computed during discretisation. Solvers instead use sparse linear solves. ([#5391](https://github.com/pybamm-team/PyBaMM/pull/5391))
 - Dropped JAX support on macOS with Intel (x86_64) processors. JAX dropped macOS Intel wheels in version 0.5.0, and the minimum JAX version has been bumped to >=0.7.0 for Python 3.14 compatibility. macOS users require Apple Silicon (M-series) for JAX features. ([#5374](https://github.com/pybamm-team/PyBaMM/pull/5374))
 - Added a small regularisation term to the exchange current density which slightly modifies the functional form of the kinetics as stoichiometry approaches 0 or 1. ([#5371](https://github.com/pybamm-team/PyBaMM/pull/5371))

src/pybamm/models/full_battery_models/base_battery_model.py

@@ -631,16 +631,6 @@ def __init__(self, extra_options):
 
         # Options not yet compatible with particle-size distributions
         if options["particle size"] == "distribution":
-            if options["lithium plating porosity change"] != "false":
-                raise pybamm.OptionError(
-                    "Lithium plating porosity change not yet supported for particle-size"
-                    " distributions."
-                )
-            if options["SEI porosity change"] == "true":
-                raise NotImplementedError(
-                    "SEI porosity change submodels do not yet support particle-size "
-                    "distributions."
-                )
             if options["heat of mixing"] != "false":
                 raise NotImplementedError(
                     "Heat of mixing submodels do not yet support particle-size "

src/pybamm/models/submodels/interface/sei/base_sei.py

@@ -84,11 +84,9 @@ def _get_standard_concentration_variables(self, c_sei):
             }
         else:
             # c_sei is a bulk quantity [mol.m-3]
-            if self.size_distribution:
-                c_sei_sav = pybamm.size_average(c_sei)
-                c_sei_xav = pybamm.x_average(c_sei_sav)
-            else:
-                c_sei_xav = pybamm.x_average(c_sei)
+            if self.size_distribution is True:
+                c_sei = pybamm.size_average(c_sei)
+            c_sei_xav = pybamm.x_average(c_sei)
             c_sei_av = pybamm.yz_average(c_sei_xav)
             variables = {
                 f"{Domain} {self.reaction_name}concentration [mol.m-3]": c_sei,

src/pybamm/models/submodels/porosity/reaction_driven_porosity.py

@@ -31,6 +31,7 @@ def get_coupled_variables(self, variables):
                 dom = domain.split()[0]
                 Domain = dom.capitalize()
                 SEI_option = getattr(self.options, dom)["SEI"]
+                plating_option = getattr(self.options, dom)["lithium plating"]
                 phases_option = getattr(self.options, dom)["particle phases"]
                 phases = self.options.phases[dom]
                 for phase in phases:
@@ -42,40 +43,34 @@ def get_coupled_variables(self, variables):
                         # `domain` has more than one phase
                         phase_name = phase + " "
                         pref = phase.capitalize() + ": "
-                    a_k = variables[
-                        f"{Domain} electrode {phase_name}"
-                        "surface area to volume ratio [m-1]"
-                    ]
                     if SEI_option == "none":
-                        L_sei_0 = pybamm.Scalar(0)
+                        V_bar_sei = pybamm.Scalar(0)
+                    else:
+                        V_bar_sei = pybamm.Parameter(
+                            f"{pref}SEI partial molar volume [m3.mol-1]"
+                        )
+                    if plating_option == "none":
+                        V_bar_Li = pybamm.Scalar(0)
                     else:
-                        L_sei_0 = pybamm.Parameter(f"{pref}Initial SEI thickness [m]")
-                    L_sei_k = variables[f"{Domain} {phase_name}SEI thickness [m]"]
-                    L_pl_k = variables[
-                        f"{Domain} {phase_name}lithium plating thickness [m]"
+                        V_bar_Li = pybamm.Parameter(
+                            "Lithium metal partial molar volume [m3.mol-1]"
+                        )
+                    c_sei_k = variables[
+                        f"{Domain} {phase_name}SEI concentration [mol.m-3]"
                     ]
-                    L_dead_k = variables[
-                        f"{Domain} {phase_name}dead lithium thickness [m]"
+                    c_sei_cr_k = variables[
+                        f"{Domain} {phase_name}SEI on cracks concentration [mol.m-3]"
                     ]
-                    L_sei_cr_k = variables[
-                        f"{Domain} {phase_name}SEI on cracks thickness [m]"
+                    c_sei_tot_k = c_sei_k + c_sei_cr_k
+                    c_pl_k = variables[
+                        f"{Domain} {phase_name}lithium plating concentration [mol.m-3]"
                     ]
-                    roughness_k = variables[
-                        f"{Domain} {phase_name}electrode roughness ratio"
+                    c_dead_k = variables[
+                        f"{Domain} {phase_name}dead lithium concentration [mol.m-3]"
                     ]
+                    c_pl_tot_k = c_pl_k + c_dead_k
 
-                    L_tot = (
-                        (L_sei_k - L_sei_0)
-                        + L_pl_k
-                        + L_dead_k
-                        + L_sei_cr_k * (roughness_k - 1)
-                    )
-
-                    # This assumes a thin film so curvature effects are neglected.
-                    # They could be included (e.g. for a sphere it is
-                    # a_n * (L_tot + L_tot ** 2 / R_n + L_tot ** # 3 / (3 * R_n ** 2)))
-                    # but it is not clear if it is relevant or not.
-                    delta_eps_k += -a_k * L_tot
+                    delta_eps_k += V_bar_sei * c_sei_tot_k + V_bar_Li * c_pl_tot_k
 
             domain_param = self.param.domain_params[domain.split()[0]]
             eps_k = domain_param.epsilon_init + delta_eps_k

tests/unit/test_models/test_full_battery_models/test_base_battery_model.py

@@ -339,13 +339,6 @@ def test_options(self):
                     "plating porosity change"
                 }
             )
-        with pytest.raises(pybamm.OptionError, match=r"distributions"):
-            pybamm.BaseBatteryModel(
-                {
-                    "particle size": "distribution",
-                    "lithium plating porosity change": "true",
-                }
-            )
 
         # contact resistance
         with pytest.raises(pybamm.OptionError, match=r"contact resistance"):

tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py

@@ -79,14 +79,6 @@ def test_nonspherical_particle_not_implemented(self):
         with pytest.raises(NotImplementedError):
             pybamm.lithium_ion.MPM(options)
 
-    def test_reversible_plating_with_porosity_not_implemented(self):
-        options = {
-            "lithium plating": "reversible",
-            "lithium plating porosity change": "true",
-        }
-        with pytest.raises(pybamm.OptionError, match=r"distributions"):
-            pybamm.lithium_ion.MPM(options)
-
     def test_msmr(self):
         options = {
             "open-circuit potential": "MSMR",
@@ -166,11 +158,3 @@ def test_interstitial_diffusion_limited(self):
         options = {"SEI": "interstitial-diffusion limited"}
         model = pybamm.lithium_ion.MPM(options)
         model.check_well_posedness()
-
-    def test_ec_reaction_limited_not_implemented(self):
-        options = {
-            "SEI": "ec reaction limited",
-            "SEI porosity change": "true",
-        }
-        with pytest.raises(NotImplementedError):
-            pybamm.lithium_ion.MPM(options)