diff --git a/docs/source/examples/notebooks/models/loss_of_active_materials.ipynb b/docs/source/examples/notebooks/models/loss_of_active_materials.ipynb index 54ace50fc6..249ae509a6 100644 --- a/docs/source/examples/notebooks/models/loss_of_active_materials.ipynb +++ b/docs/source/examples/notebooks/models/loss_of_active_materials.ipynb @@ -443,6 +443,8 @@ "second = 0.1\n", "parameter_values.update(\n", " {\n", + " \"Primary: Negative electrode reference concentration for free of deformation [mol.m-3]\": 0.0,\n", + " \"Secondary: Negative electrode reference concentration for free of deformation [mol.m-3]\": 0.0,\n", " \"Primary: Negative electrode LAM constant proportional term [s-1]\": 1e-4 / 3600,\n", " \"Secondary: Negative electrode LAM constant proportional term [s-1]\": 1e-4\n", " / 3600\n", diff --git a/src/pybamm/models/full_battery_models/base_battery_model.py b/src/pybamm/models/full_battery_models/base_battery_model.py index e47a6c05c8..83d62df278 100644 --- a/src/pybamm/models/full_battery_models/base_battery_model.py +++ b/src/pybamm/models/full_battery_models/base_battery_model.py @@ -525,14 +525,14 @@ def __init__(self, extra_options): # Options not yet compatible with particle-size distributions if options["particle size"] == "distribution": - if options["heat of mixing"] != "false": + if options["lithium plating porosity change"] != "false": raise NotImplementedError( - "Heat of mixing submodels do not yet support particle-size " - "distributions." + "Lithium plating porosity change not yet supported for particle-size" + " distributions." ) - if options["lithium plating"] != "none": + if options["heat of mixing"] != "false": raise NotImplementedError( - "Lithium plating submodels do not yet support particle-size " + "Heat of mixing submodels do not yet support particle-size " "distributions." ) if options["particle"] in ["quadratic profile", "quartic profile"]: diff --git a/src/pybamm/models/submodels/interface/base_interface.py b/src/pybamm/models/submodels/interface/base_interface.py index 83f7c81b43..229aa5fdf2 100644 --- a/src/pybamm/models/submodels/interface/base_interface.py +++ b/src/pybamm/models/submodels/interface/base_interface.py @@ -326,9 +326,10 @@ def _get_standard_volumetric_current_density_variables(self, variables): a_j_av = pybamm.x_average(a_j) if reaction_name == "SEI on cracks ": - roughness = variables[f"{Domain} electrode roughness ratio"] - 1 + roughness = variables[f"{Domain} {phase_name}electrode roughness ratio"] - 1 roughness_av = ( - variables[f"X-averaged {domain} electrode roughness ratio"] - 1 + variables[f"X-averaged {domain} {phase_name}electrode roughness ratio"] + - 1 ) else: roughness = 1 diff --git a/src/pybamm/models/submodels/interface/sei/base_sei.py b/src/pybamm/models/submodels/interface/sei/base_sei.py index af5a929bb5..d801e25e61 100644 --- a/src/pybamm/models/submodels/interface/sei/base_sei.py +++ b/src/pybamm/models/submodels/interface/sei/base_sei.py @@ -120,6 +120,7 @@ def _get_standard_concentration_variables(self, variables): """Update variables related to the SEI concentration.""" domain, Domain = self.domain_Domain phase_param = self.phase_param + phase_name = phase_param.phase_name reaction_name = self.reaction_name # Set scales to one for the "no SEI" model so that they are not required @@ -189,7 +190,7 @@ def _get_standard_concentration_variables(self, variables): # Concentration variables are handled slightly differently for SEI on cracks elif self.reaction == "SEI on cracks": L_cr = variables[f"{Domain} {reaction_name}thickness [m]"] - roughness = variables[f"{Domain} electrode roughness ratio"] + roughness = variables[f"{Domain} {phase_name}electrode roughness ratio"] n_SEI_cr = L_cr * L_to_n * (roughness - 1) # SEI on cracks concentration if self.size_distribution: diff --git a/src/pybamm/models/submodels/particle/base_particle.py b/src/pybamm/models/submodels/particle/base_particle.py index b5b8fb4499..e47039b973 100644 --- a/src/pybamm/models/submodels/particle/base_particle.py +++ b/src/pybamm/models/submodels/particle/base_particle.py @@ -29,7 +29,6 @@ def __init__(self, param, domain, options, phase="primary"): def _get_effective_diffusivity(self, c, T, current): domain, Domain = self.domain_Domain - domain_param = self.domain_param phase_param = self.phase_param domain_options = getattr(self.options, domain) @@ -60,7 +59,7 @@ def _get_effective_diffusivity(self, c, T, current): E = pybamm.r_average(phase_param.E(sto, T)) nu = phase_param.nu theta_M = Omega / (self.param.R * T) * (2 * Omega * E) / (9 * (1 - nu)) - stress_factor = 1 + theta_M * (c - domain_param.c_0) + stress_factor = 1 + theta_M * (c - phase_param.c_0) else: stress_factor = 1 diff --git a/src/pybamm/models/submodels/particle_mechanics/base_mechanics.py b/src/pybamm/models/submodels/particle_mechanics/base_mechanics.py index 1301722da0..c41c12fe1b 100644 --- a/src/pybamm/models/submodels/particle_mechanics/base_mechanics.py +++ b/src/pybamm/models/submodels/particle_mechanics/base_mechanics.py @@ -30,8 +30,8 @@ def _get_standard_variables(self, l_cr): domain, Domain = self.domain_Domain l_cr_av = pybamm.x_average(l_cr) variables = { - f"{Domain} particle crack length [m]": l_cr, - f"X-averaged {domain} particle crack length [m]": l_cr_av, + f"{Domain} {self.phase_param.phase_name}particle crack length [m]": l_cr, + f"X-averaged {domain} {self.phase_param.phase_name}particle crack length [m]": l_cr_av, } return variables @@ -61,7 +61,7 @@ def _get_mechanical_results(self, variables): sto = variables[f"{Domain} {phase_name}particle concentration"] Omega = pybamm.r_average(phase_param.Omega(sto, T)) R0 = phase_param.R - c_0 = domain_param.c_0 + c_0 = phase_param.c_0 E0 = pybamm.r_average(phase_param.E(sto, T)) nu = phase_param.nu L0 = domain_param.L @@ -127,21 +127,21 @@ def _get_mechanical_results(self, variables): def _get_standard_surface_variables(self, variables): domain, Domain = self.domain_Domain - phase_name = self.phase_name + phase_name = self.phase_param.phase_name - l_cr = variables[f"{Domain} particle crack length [m]"] + l_cr = variables[f"{Domain} {phase_name}particle crack length [m]"] a = variables[ f"{Domain} electrode {phase_name}surface area to volume ratio [m-1]" ] - rho_cr = self.domain_param.rho_cr - w_cr = self.domain_param.w_cr + rho_cr = self.phase_param.rho_cr + w_cr = self.phase_param.w_cr roughness = 1 + 2 * l_cr * rho_cr * w_cr # ratio of cracks to normal surface a_cr = (roughness - 1) * a # crack surface area to volume ratio roughness_xavg = pybamm.x_average(roughness) variables = { - f"{Domain} crack surface to volume ratio [m-1]": a_cr, - f"{Domain} electrode roughness ratio": roughness, - f"X-averaged {domain} electrode roughness ratio": roughness_xavg, + f"{Domain} {phase_name}crack surface to volume ratio [m-1]": a_cr, + f"{Domain} {phase_name}electrode roughness ratio": roughness, + f"X-averaged {domain} {phase_name}electrode roughness ratio": roughness_xavg, } return variables diff --git a/src/pybamm/models/submodels/particle_mechanics/crack_propagation.py b/src/pybamm/models/submodels/particle_mechanics/crack_propagation.py index b51f1d1ebd..5a84f870d8 100644 --- a/src/pybamm/models/submodels/particle_mechanics/crack_propagation.py +++ b/src/pybamm/models/submodels/particle_mechanics/crack_propagation.py @@ -36,20 +36,20 @@ 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} particle crack length [m]", + f"X-averaged {domain} {phase_name}particle crack length [m]", domain="current collector", - scale=self.domain_param.l_cr_0, + scale=self.phase_param.l_cr_0, ) l_cr = pybamm.PrimaryBroadcast(l_cr_av, f"{domain} electrode") else: l_cr = pybamm.Variable( - f"{Domain} particle crack length [m]", + f"{Domain} {phase_name}particle crack length [m]", domain=f"{domain} electrode", auxiliary_domains={"secondary": "current collector"}, - scale=self.domain_param.l_cr_0, + scale=self.phase_param.l_cr_0, ) variables = self._get_standard_variables(l_cr) @@ -58,22 +58,24 @@ def get_fundamental_variables(self): def get_coupled_variables(self, variables): domain, Domain = self.domain_Domain - + phase_name = self.phase_param.phase_name variables.update(self._get_standard_surface_variables(variables)) variables.update(self._get_mechanical_results(variables)) T = variables[f"{Domain} electrode temperature [K]"] - k_cr = self.domain_param.k_cr(T) - m_cr = self.domain_param.m_cr - b_cr = self.domain_param.b_cr - stress_t_surf = variables[f"{Domain} particle surface tangential stress [Pa]"] - l_cr = variables[f"{Domain} particle crack length [m]"] + 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]"] # # 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 variables.update( { - f"{Domain} particle cracking rate [m.s-1]": dl_cr, - f"X-averaged {domain} particle cracking rate [m.s-1]": pybamm.x_average( + f"{Domain} {phase_name}particle cracking rate [m.s-1]": dl_cr, + f"X-averaged {domain} {phase_name}particle cracking rate [m.s-1]": pybamm.x_average( dl_cr ), } @@ -82,36 +84,44 @@ def get_coupled_variables(self, variables): 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} particle crack length [m]"] - dl_cr = variables[f"X-averaged {domain} particle cracking rate [m.s-1]"] + 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} particle crack length [m]"] - dl_cr = variables[f"{Domain} particle cracking rate [m.s-1]"] + 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} def set_initial_conditions(self, variables): domain, Domain = self.domain_Domain - - l_cr_0 = self.domain_param.l_cr_0 + 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} particle crack length [m]"] + l_cr = variables[ + f"X-averaged {domain} {phase_name}particle crack length [m]" + ] else: - l_cr = variables[f"{Domain} particle crack length [m]"] + 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): 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} particle crack length [m]"] + l_cr = variables[ + f"X-averaged {domain} {phase_name}particle crack length [m]" + ] else: - l_cr = variables[f"{Domain} particle crack length [m]"] + l_cr = variables[f"{Domain} {phase_name}particle crack length [m]"] self.events.append( pybamm.Event( - f"{domain} particle crack length larger than particle radius", - 1 - pybamm.max(l_cr) / self.domain_param.prim.R_typ, + f"{domain} {phase_name} particle crack length larger than particle radius", + 1 - pybamm.max(l_cr) / self.phase_param.R_typ, ) ) diff --git a/src/pybamm/parameters/lead_acid_parameters.py b/src/pybamm/parameters/lead_acid_parameters.py index 7da6a19410..7beef48142 100644 --- a/src/pybamm/parameters/lead_acid_parameters.py +++ b/src/pybamm/parameters/lead_acid_parameters.py @@ -344,6 +344,7 @@ def __init__(self, phase, domain_param): self.domain = domain_param.domain self.main_param = domain_param.main_param self.geo = domain_param.geo.prim + self.phase_name = phase def _set_parameters(self): main = self.main_param diff --git a/src/pybamm/parameters/lithium_ion_parameters.py b/src/pybamm/parameters/lithium_ion_parameters.py index 1ea20aeb3e..1ccd21975a 100644 --- a/src/pybamm/parameters/lithium_ion_parameters.py +++ b/src/pybamm/parameters/lithium_ion_parameters.py @@ -269,20 +269,6 @@ def _set_parameters(self): self.b_s = self.geo.b_s self.tau_s = self.geo.tau_s - # Mechanical parameters - self.c_0 = pybamm.Parameter( - f"{Domain} electrode reference concentration for free of deformation " - "[mol.m-3]" - ) - - self.l_cr_0 = pybamm.Parameter(f"{Domain} electrode initial crack length [m]") - self.w_cr = pybamm.Parameter(f"{Domain} electrode initial crack width [m]") - self.rho_cr = pybamm.Parameter( - f"{Domain} electrode number of cracks per unit area [m-2]" - ) - self.b_cr = pybamm.Parameter(f"{Domain} electrode Paris' law constant b") - self.m_cr = pybamm.Parameter(f"{Domain} electrode Paris' law constant m") - # Utilisation parameters self.u_init = pybamm.Parameter( f"Initial {domain} electrode interface utilisation" @@ -307,15 +293,6 @@ def sigma(self, T): f"{Domain} electrode conductivity [S.m-1]", inputs ) - def k_cr(self, T): - """ - Cracking rate for the electrode; - """ - Domain = self.domain.capitalize() - return pybamm.FunctionParameter( - f"{Domain} electrode cracking rate", {"Temperature [K]": T} - ) - def LAM_rate_current(self, i, T): """ Dimensional rate of loss of active material as a function of applied current @@ -516,6 +493,34 @@ def _set_parameters(self): f"{pref}{Domain} electrode reaction-driven LAM factor [m3.mol-1]" ) + # Mechanical parameters + self.c_0 = pybamm.Parameter( + f"{pref}{Domain} electrode reference concentration for free of deformation " + "[mol.m-3]" + ) + + self.l_cr_0 = pybamm.Parameter( + f"{pref}{Domain} electrode initial crack length [m]" + ) + self.w_cr = pybamm.Parameter( + f"{pref}{Domain} electrode initial crack width [m]" + ) + self.rho_cr = pybamm.Parameter( + f"{pref}{Domain} electrode number of cracks per unit area [m-2]" + ) + self.b_cr = pybamm.Parameter(f"{pref}{Domain} electrode Paris' law constant b") + self.m_cr = pybamm.Parameter(f"{pref}{Domain} electrode Paris' law constant m") + + def k_cr(self, T): + """ + Cracking rate for the electrode; + """ + phase_prefactor = self.phase_prefactor + Domain = self.domain.capitalize() + return pybamm.FunctionParameter( + f"{phase_prefactor}{Domain} electrode cracking rate", {"Temperature [K]": T} + ) + def D(self, c_s, T, lithiation=None): """ Dimensional diffusivity in particle. In the parameter sets this is defined as diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py index 768872e454..2456a89b60 100644 --- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py +++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/base_lithium_ion_tests.py @@ -462,7 +462,8 @@ def lico2_volume_change_Ai2020(sto): "Secondary: Negative electrode partial molar volume [m3.mol-1]": 3.1e-06, "Secondary: Negative electrode Young's modulus [Pa]": 15000000000.0, "Secondary: Negative electrode Poisson's ratio": 0.3, - "Negative electrode reference concentration for free of deformation [mol.m-3]": 0.0, + "Primary: Negative electrode reference concentration for free of deformation [mol.m-3]": 0.0, + "Secondary: Negative electrode reference concentration for free of deformation [mol.m-3]": 0.0, "Primary: Negative electrode volume change": graphite_volume_change_Ai2020, "Secondary: Negative electrode volume change": graphite_volume_change_Ai2020, "Positive electrode partial molar volume [m3.mol-1]": -7.28e-07, diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py index ad02212840..5f693264d3 100644 --- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py +++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_mpm.py @@ -123,6 +123,18 @@ def test_wycisk_ocp(self): model = pybamm.lithium_ion.MPM(options) model.check_well_posedness() + def test_mpm_with_lithium_plating(self): + options = { + "lithium plating": "irreversible", + } + model = pybamm.lithium_ion.MPM(options) + model.check_well_posedness() + options = { + "lithium plating": "reversible", + } + model = pybamm.lithium_ion.MPM(options) + model.check_well_posedness() + class TestMPMExternalCircuits: def test_well_posed_voltage(self): @@ -196,15 +208,3 @@ def test_well_posed_both_swelling_only_not_implemented(self): options = {"particle mechanics": "swelling only"} with pytest.raises(NotImplementedError): pybamm.lithium_ion.MPM(options) - - -class TestMPMWithPlating: - def test_well_posed_reversible_plating_not_implemented(self): - options = {"lithium plating": "reversible"} - with pytest.raises(NotImplementedError): - pybamm.lithium_ion.MPM(options) - - def test_well_posed_irreversible_plating_not_implemented(self): - options = {"lithium plating": "irreversible"} - with pytest.raises(NotImplementedError): - pybamm.lithium_ion.MPM(options)