Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Adds support for particle size distributions combined with particle mechanics. #4807

Merged
merged 9 commits into from
Feb 3, 2025
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 4 additions & 0 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
@@ -1,5 +1,9 @@
# [Unreleased](https://github.com/pybamm-team/PyBaMM/)

## Features

- Added support for particle size distributions combined with particle mechanics. ([#4807](https://github.com/pybamm-team/PyBaMM/pull/4807))

# [v25.1.1](https://github.com/pybamm-team/PyBaMM/tree/v25.1.1) - 2025-01-20

## Features
Expand Down
10 changes: 0 additions & 10 deletions src/pybamm/models/full_battery_models/base_battery_model.py
Original file line number Diff line number Diff line change
Expand Up @@ -545,21 +545,11 @@ def __init__(self, extra_options):
"'quadratic' and 'quartic' concentration profiles have not yet "
"been implemented for particle-size ditributions"
)
if options["particle mechanics"] != "none":
raise NotImplementedError(
"Particle mechanics submodels do not yet support particle-size"
" distributions."
)
if options["particle shape"] != "spherical":
raise NotImplementedError(
"Particle shape must be 'spherical' for particle-size distribution"
" submodels."
)
if options["stress-induced diffusion"] == "true":
raise NotImplementedError(
"stress-induced diffusion cannot yet be included in "
"particle-size distributions."
)
if options["thermal"] == "x-full":
raise NotImplementedError(
"X-full thermal submodels do not yet support particle-size"
Expand Down
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
#
# Class for varying active material volume fraction, driven by stress
# Class for varying active material volume fraction
#
import pybamm

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,10 @@ class BaseMechanics(pybamm.BaseSubModel):
"""

def __init__(self, param, domain, options, phase="primary"):
if options["particle size"] == "distribution":
self.size_distribution = True
else:
self.size_distribution = False
super().__init__(param, domain, options=options, phase=phase)

def _get_standard_variables(self, l_cr):
Expand All @@ -35,6 +39,71 @@ def _get_standard_variables(self, l_cr):
}
return variables

def _get_standard_size_distribution_variables(self, l_cr_dist):
domain, Domain = self.domain_Domain
l_cr_av_dist = pybamm.x_average(l_cr_dist)
variables = {
f"{Domain} {self.phase_param.phase_name}particle crack length distribution [m]": l_cr_dist,
f"X-averaged {domain} {self.phase_param.phase_name}particle crack length distribution [m]": l_cr_av_dist,
}
return variables

def _get_mechanical_size_distribution_results(self, variables):
domain, Domain = self.domain_Domain
phase_name = self.phase_param.phase_name
phase_param = self.phase_param
c_s_rav = variables[
f"R-averaged {domain} {phase_name}particle concentration distribution [mol.m-3]"
]
c_s_surf = variables[
f"{Domain} {phase_name}particle surface concentration distribution [mol.m-3]"
]
T = pybamm.PrimaryBroadcast(
variables[f"{Domain} electrode temperature [K]"],
[f"{domain} {phase_name}particle size"],
)
T = pybamm.PrimaryBroadcast(
T,
[f"{domain} {phase_name}particle"],
)

# use a tangential approximation for omega
c_0 = phase_param.c_0
R0 = phase_param.R
sto = variables[f"{Domain} {phase_name}particle concentration distribution"]
Omega = pybamm.r_average(phase_param.Omega(sto, T))

E0 = pybamm.r_average(phase_param.E(sto, T))
nu = phase_param.nu
# Ai2019 eq [10]
disp_surf = Omega * R0 / 3 * (c_s_rav - c_0)
# c0 reference concentration for no deformation
# stress evaluated at the surface of the particles
# Ai2019 eq [7] with r=R
stress_r_surf = pybamm.Scalar(0)
# Ai2019 eq [8] with r=R
# c_s_rav is already multiplied by 3/R^3 inside r_average
stress_t_surf = Omega * E0 * (c_s_rav - c_s_surf) / 3.0 / (1.0 - nu)

# Averages
stress_r_surf_av = pybamm.x_average(stress_r_surf)
stress_t_surf_av = pybamm.x_average(stress_t_surf)
disp_surf_av = pybamm.x_average(disp_surf)

variables.update(
{
f"{Domain} {phase_name}particle surface radial stress distribution [Pa]": stress_r_surf,
f"{Domain} {phase_name}particle surface tangential stress distribution [Pa]": stress_t_surf,
f"{Domain} {phase_name}particle surface displacement distribution [m]": disp_surf,
f"X-averaged {domain} {phase_name}particle surface "
"radial stress distribution [Pa]": stress_r_surf_av,
f"X-averaged {domain} {phase_name}particle surface "
"tangential stress distribution [Pa]": stress_t_surf_av,
f"X-averaged {domain} {phase_name}particle surface displacement distribution [m]": disp_surf_av,
}
)
return variables

def _get_mechanical_results(self, variables):
domain_param = self.domain_param
domain, Domain = self.domain_Domain
Expand All @@ -58,10 +127,11 @@ def _get_mechanical_results(self, variables):
]

# use a tangential approximation for omega
c_0 = phase_param.c_0
R0 = phase_param.R
sto = variables[f"{Domain} {phase_name}particle concentration"]
Omega = pybamm.r_average(phase_param.Omega(sto, T))
R0 = phase_param.R
c_0 = phase_param.c_0

E0 = pybamm.r_average(phase_param.E(sto, T))
nu = phase_param.nu
L0 = domain_param.L
Expand All @@ -79,7 +149,7 @@ def _get_mechanical_results(self, variables):
stress_r_surf = pybamm.Scalar(0)
# Ai2019 eq [8] with r=R
# c_s_rav is already multiplied by 3/R^3 inside r_average
stress_t_surf = Omega * E0 / 3.0 / (1.0 - nu) * (c_s_rav - c_s_surf)
stress_t_surf = Omega * E0 * (c_s_rav - c_s_surf) / 3.0 / (1.0 - nu)

# Averages
stress_r_surf_av = pybamm.x_average(stress_r_surf)
Expand Down
113 changes: 88 additions & 25 deletions src/pybamm/models/submodels/particle_mechanics/crack_propagation.py
Original file line number Diff line number Diff line change
Expand Up @@ -37,22 +37,50 @@ def __init__(self, param, domain, x_average, options, phase="primary"):
def get_fundamental_variables(self):
domain, Domain = self.domain_Domain
phase_name = self.phase_param.phase_name
if self.x_average is True:
l_cr_av = pybamm.Variable(
f"X-averaged {domain} {phase_name}particle crack length [m]",
domain="current collector",
scale=self.phase_param.l_cr_0,
)
if self.x_average:
if self.size_distribution:
l_cr_av_dist = pybamm.Variable(
f"X-averaged {domain} {phase_name}particle crack length distribution [m]",
domains={
"primary": f"{domain} particle size",
"secondary": "current collector",
},
scale=self.phase_param.l_cr_0,
)
l_cr_dist = pybamm.SecondaryBroadcast(
l_cr_av_dist, f"{domain} electrode"
)
l_cr_av = pybamm.size_average(l_cr_av_dist)
else:
l_cr_av = pybamm.Variable(
f"X-averaged {domain} {phase_name}particle crack length [m]",
domain="current collector",
scale=self.phase_param.l_cr_0,
)
l_cr = pybamm.PrimaryBroadcast(l_cr_av, f"{domain} electrode")
else:
l_cr = pybamm.Variable(
f"{Domain} {phase_name}particle crack length [m]",
domain=f"{domain} electrode",
auxiliary_domains={"secondary": "current collector"},
scale=self.phase_param.l_cr_0,
)
if self.size_distribution:
l_cr_dist = pybamm.Variable(
f"{Domain} {phase_name}particle crack length distribution [m]",
domains={
"primary": f"{domain} particle size",
"secondary": f"{domain} electrode",
"tertiary": "current collector",
},
scale=self.phase_param.l_cr_0,
)
l_cr = pybamm.size_average(l_cr_dist)
else:
l_cr = pybamm.Variable(
f"{Domain} {phase_name}particle crack length [m]",
domain=f"{domain} electrode",
auxiliary_domains={"secondary": "current collector"},
scale=self.phase_param.l_cr_0,
)

variables = self._get_standard_variables(l_cr)
if self.size_distribution:
variables.update(self._get_standard_size_distribution_variables(l_cr_dist))

return variables

Expand All @@ -61,14 +89,26 @@ def get_coupled_variables(self, variables):
phase_name = self.phase_param.phase_name
variables.update(self._get_standard_surface_variables(variables))
variables.update(self._get_mechanical_results(variables))
if self.size_distribution:
variables.update(self._get_mechanical_size_distribution_results(variables))
T = variables[f"{Domain} electrode temperature [K]"]
k_cr = self.phase_param.k_cr(T)
m_cr = self.phase_param.m_cr
b_cr = self.phase_param.b_cr
stress_t_surf = variables[
f"{Domain} {phase_name}particle surface tangential stress [Pa]"
]
l_cr = variables[f"{Domain} {phase_name}particle crack length [m]"]
if self.size_distribution:
stress_t_surf = variables[
f"{Domain} {phase_name}particle surface tangential stress distribution [Pa]"
]
else:
stress_t_surf = variables[
f"{Domain} {phase_name}particle surface tangential stress [Pa]"
]
if self.size_distribution:
l_cr = variables[
f"{Domain} {phase_name}particle crack length distribution [m]"
]
else:
l_cr = variables[f"{Domain} {phase_name}particle crack length [m]"]
# # compressive stress will not lead to crack propagation
dK_SIF = stress_t_surf * b_cr * pybamm.sqrt(np.pi * l_cr) * (stress_t_surf >= 0)
dl_cr = k_cr * (dK_SIF**m_cr) / 3600 # divide by 3600 to replace t0_cr
Expand All @@ -86,14 +126,24 @@ def set_rhs(self, variables):
domain, Domain = self.domain_Domain
phase_name = self.phase_param.phase_name
if self.x_average is True:
l_cr = variables[
f"X-averaged {domain} {phase_name}particle crack length [m]"
]
if self.size_distribution:
l_cr = variables[
f"X-averaged {domain} {phase_name}particle crack length distribution [m]"
]
else:
l_cr = variables[
f"X-averaged {domain} {phase_name}particle crack length [m]"
]
dl_cr = variables[
f"X-averaged {domain} {phase_name}particle cracking rate [m.s-1]"
]
else:
l_cr = variables[f"{Domain} {phase_name}particle crack length [m]"]
if self.size_distribution:
l_cr = variables[
f"{Domain} {phase_name}particle crack length distribution [m]"
]
else:
l_cr = variables[f"{Domain} {phase_name}particle crack length [m]"]
dl_cr = variables[f"{Domain} {phase_name}particle cracking rate [m.s-1]"]
self.rhs = {l_cr: dl_cr}

Expand All @@ -102,12 +152,25 @@ def set_initial_conditions(self, variables):
phase_name = self.phase_param.phase_name
l_cr_0 = self.phase_param.l_cr_0
if self.x_average is True:
l_cr = variables[
f"X-averaged {domain} {phase_name}particle crack length [m]"
]
if self.size_distribution:
l_cr = variables[
f"X-averaged {domain} {phase_name}particle crack length distribution [m]"
]
l_cr_0 = pybamm.PrimaryBroadcast(l_cr_0, f"{domain} particle size")
else:
l_cr = variables[
f"X-averaged {domain} {phase_name}particle crack length [m]"
]
else:
l_cr = variables[f"{Domain} {phase_name}particle crack length [m]"]
l_cr_0 = pybamm.PrimaryBroadcast(l_cr_0, f"{domain} electrode")
if self.size_distribution:
l_cr = variables[
f"{Domain} {phase_name}particle crack length distribution [m]"
]
l_cr_0 = pybamm.PrimaryBroadcast(l_cr_0, f"{domain} electrode")
l_cr_0 = pybamm.PrimaryBroadcast(l_cr_0, f"{domain} particle size")
else:
l_cr = variables[f"{Domain} {phase_name}particle crack length [m]"]
l_cr_0 = pybamm.PrimaryBroadcast(l_cr_0, f"{domain} electrode")
self.initial_conditions = {l_cr: l_cr_0}

def add_events_from(self, variables):
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -48,4 +48,6 @@ def get_fundamental_variables(self):
def get_coupled_variables(self, variables):
variables.update(self._get_standard_surface_variables(variables))
variables.update(self._get_mechanical_results(variables))
if self.size_distribution:
variables.update(self._get_mechanical_size_distribution_results(variables))
return variables
Original file line number Diff line number Diff line change
Expand Up @@ -525,6 +525,22 @@ def test_well_posed_psd(self):
options = {"particle size": "distribution", "surface form": "algebraic"}
self.check_well_posedness(options)

def test_well_posed_psd_swelling_and_cracking(self):
options = {
"particle size": "distribution",
"particle mechanics": "swelling and cracking",
"surface form": "algebraic",
}
self.check_well_posedness(options)

def test_well_posed_psd_swelling_only(self):
options = {
"particle size": "distribution",
"particle mechanics": "swelling only",
"surface form": "algebraic",
}
self.check_well_posedness(options)

def test_well_posed_transport_efficiency_Bruggeman(self):
options = {"transport efficiency": "Bruggeman"}
self.check_well_posedness(options)
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -78,21 +78,6 @@ def test_nonspherical_particle_not_implemented(self):
with pytest.raises(NotImplementedError):
pybamm.lithium_ion.MPM(options)

def test_loss_active_material_stress_negative_not_implemented(self):
options = {"loss of active material": ("stress-driven", "none")}
with pytest.raises(NotImplementedError):
pybamm.lithium_ion.MPM(options)

def test_loss_active_material_stress_positive_not_implemented(self):
options = {"loss of active material": ("none", "stress-driven")}
with pytest.raises(NotImplementedError):
pybamm.lithium_ion.MPM(options)

def test_loss_active_material_stress_both_not_implemented(self):
options = {"loss of active material": "stress-driven"}
with pytest.raises(NotImplementedError):
pybamm.lithium_ion.MPM(options)

def test_reversible_plating_with_porosity_not_implemented(self):
options = {
"lithium plating": "reversible",
Expand All @@ -101,11 +86,6 @@ def test_reversible_plating_with_porosity_not_implemented(self):
with pytest.raises(pybamm.OptionError, match="distributions"):
pybamm.lithium_ion.MPM(options)

def test_stress_induced_diffusion_not_implemented(self):
options = {"stress-induced diffusion": "true"}
with pytest.raises(NotImplementedError):
pybamm.lithium_ion.MPM(options)

def test_msmr(self):
options = {
"open-circuit potential": "MSMR",
Expand Down Expand Up @@ -186,25 +166,3 @@ def test_ec_reaction_limited_not_implemented(self):
}
with pytest.raises(NotImplementedError):
pybamm.lithium_ion.MPM(options)


class TestMPMWithMechanics:
def test_well_posed_negative_cracking_not_implemented(self):
options = {"particle mechanics": ("swelling and cracking", "none")}
with pytest.raises(NotImplementedError):
pybamm.lithium_ion.MPM(options)

def test_well_posed_positive_cracking_not_implemented(self):
options = {"particle mechanics": ("none", "swelling and cracking")}
with pytest.raises(NotImplementedError):
pybamm.lithium_ion.MPM(options)

def test_well_posed_both_cracking_not_implemented(self):
options = {"particle mechanics": "swelling and cracking"}
with pytest.raises(NotImplementedError):
pybamm.lithium_ion.MPM(options)

def test_well_posed_both_swelling_only_not_implemented(self):
options = {"particle mechanics": "swelling only"}
with pytest.raises(NotImplementedError):
pybamm.lithium_ion.MPM(options)