diff --git a/CHANGELOG.md b/CHANGELOG.md index 24c9be68e0..8322884eda 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Features +- Made composite electrode model compatible with particle size distribution ([#4687](https://github.com/pybamm-team/PyBaMM/pull/4687)) - Added `Symbol.post_order()` method to return an iterable that steps through the tree in post-order fashion. ([#4684](https://github.com/pybamm-team/PyBaMM/pull/4684)) - Added two more submodels (options) for the SEI: Lars von Kolzenberg (2020) model and Tunneling Limit model ([#4394](https://github.com/pybamm-team/PyBaMM/pull/4394)) diff --git a/docs/source/examples/notebooks/models/using-submodels.ipynb b/docs/source/examples/notebooks/models/using-submodels.ipynb index cce32fcefd..db889972f8 100644 --- a/docs/source/examples/notebooks/models/using-submodels.ipynb +++ b/docs/source/examples/notebooks/models/using-submodels.ipynb @@ -527,10 +527,10 @@ " model.param, \"positive\", model.options, cracks=True\n", ")\n", "model.submodels[\"Negative lithium plating\"] = pybamm.lithium_plating.NoPlating(\n", - " model.param, \"Negative\"\n", + " model.param, \"negative\"\n", ")\n", "model.submodels[\"Positive lithium plating\"] = pybamm.lithium_plating.NoPlating(\n", - " model.param, \"Positive\"\n", + " model.param, \"positive\"\n", ")" ] }, diff --git a/src/pybamm/expression_tree/averages.py b/src/pybamm/expression_tree/averages.py index 11538ea153..c95ccb4e5e 100644 --- a/src/pybamm/expression_tree/averages.py +++ b/src/pybamm/expression_tree/averages.py @@ -306,6 +306,10 @@ def xyz_average(symbol: pybamm.Symbol) -> pybamm.Symbol: return yz_average(x_average(symbol)) +def xyzs_average(symbol: pybamm.Symbol) -> pybamm.Symbol: + return xyz_average(size_average(symbol)) + + def r_average(symbol: pybamm.Symbol) -> pybamm.Symbol: """ Convenience function for creating an average in the r-direction. @@ -373,7 +377,15 @@ def size_average( # If symbol doesn't have a domain, or doesn't have "negative particle size" # or "positive particle size" as a domain, it's average value is itself if symbol.domain == [] or not any( - domain in [["negative particle size"], ["positive particle size"]] + domain + in [ + ["negative particle size"], + ["positive particle size"], + ["negative primary particle size"], + ["positive primary particle size"], + ["negative secondary particle size"], + ["positive secondary particle size"], + ] for domain in list(symbol.domains.values()) ): return symbol @@ -394,13 +406,30 @@ def size_average( else: if f_a_dist is None: geo = pybamm.geometric_parameters + name = "R" + if "negative" in symbol.domain[0]: + name += "_n" + elif "positive" in symbol.domain[0]: + name += "_p" + if "primary" in symbol.domain[0]: + name += "_prim" + elif "secondary" in symbol.domain[0]: + name += "_sec" R = pybamm.SpatialVariable( - "R", domains=symbol.domains, coord_sys="cartesian" + name, domains=symbol.domains, coord_sys="cartesian" ) - if ["negative particle size"] in symbol.domains.values(): + if ["negative particle size"] in symbol.domains.values() or [ + "negative primary particle size" + ] in symbol.domains.values(): f_a_dist = geo.n.prim.f_a_dist(R) - elif ["positive particle size"] in symbol.domains.values(): + elif ["negative secondary particle size"] in symbol.domains.values(): + f_a_dist = geo.n.sec.f_a_dist(R) + elif ["positive particle size"] in symbol.domains.values() or [ + "positive primary particle size" + ] in symbol.domains.values(): f_a_dist = geo.p.prim.f_a_dist(R) + elif ["positive secondary particle size"] in symbol.domains.values(): + f_a_dist = geo.p.sec.f_a_dist(R) return SizeAverage(symbol, f_a_dist) diff --git a/src/pybamm/geometry/battery_geometry.py b/src/pybamm/geometry/battery_geometry.py index 9be08ff619..3e1986507a 100644 --- a/src/pybamm/geometry/battery_geometry.py +++ b/src/pybamm/geometry/battery_geometry.py @@ -83,6 +83,18 @@ def battery_geometry( "negative particle size": {"R_n": {"min": R_min_n, "max": R_max_n}}, } ) + phases = int(options.negative["particle phases"]) + if phases >= 2: + geometry.update( + { + "negative primary particle size": { + "R_n_prim": {"min": R_min_n, "max": R_max_n} + }, + "negative secondary particle size": { + "R_n_sec": {"min": R_min_n, "max": R_max_n} + }, + } + ) if ( options is not None and options.positive["particle size"] == "distribution" @@ -95,6 +107,18 @@ def battery_geometry( "positive particle size": {"R_p": {"min": R_min_p, "max": R_max_p}}, } ) + phases = int(options.positive["particle phases"]) + if phases >= 2: + geometry.update( + { + "positive primary particle size": { + "R_p_prim": {"min": R_min_p, "max": R_max_p} + }, + "positive secondary particle size": { + "R_p_sec": {"min": R_min_p, "max": R_max_p} + }, + } + ) # Add current collector domains current_collector_dimension = options["dimensionality"] diff --git a/src/pybamm/geometry/standard_spatial_vars.py b/src/pybamm/geometry/standard_spatial_vars.py index 565908deb9..b6638694f7 100644 --- a/src/pybamm/geometry/standard_spatial_vars.py +++ b/src/pybamm/geometry/standard_spatial_vars.py @@ -98,6 +98,7 @@ }, coord_sys="cartesian", ) + R_p = pybamm.SpatialVariable( "R_p", domain=["positive particle size"], @@ -108,6 +109,44 @@ coord_sys="cartesian", ) +R_n_prim = pybamm.SpatialVariable( + "R_n_prim", + domain=["negative primary particle size"], + auxiliary_domains={ + "secondary": "negative electrode", + "tertiary": "current collector", + }, + coord_sys="cartesian", +) +R_p_prim = pybamm.SpatialVariable( + "R_p_prim", + domain=["positive primary particle size"], + auxiliary_domains={ + "secondary": "positive electrode", + "tertiary": "current collector", + }, + coord_sys="cartesian", +) + +R_n_sec = pybamm.SpatialVariable( + "R_n_sec", + domain=["negative secondary particle size"], + auxiliary_domains={ + "secondary": "negative electrode", + "tertiary": "current collector", + }, + coord_sys="cartesian", +) +R_p_sec = pybamm.SpatialVariable( + "R_p_sec", + domain=["positive secondary particle size"], + auxiliary_domains={ + "secondary": "positive electrode", + "tertiary": "current collector", + }, + coord_sys="cartesian", +) + # Domains at cell edges x_n_edge = pybamm.SpatialVariableEdge( "x_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 fd404cec75..e47a6c05c8 100644 --- a/src/pybamm/models/full_battery_models/base_battery_model.py +++ b/src/pybamm/models/full_battery_models/base_battery_model.py @@ -623,7 +623,6 @@ def __init__(self, extra_options): if options["particle phases"] not in ["1", ("1", "1")]: if not ( options["surface form"] != "false" - and options["particle size"] == "single" and options["particle"] == "Fickian diffusion" ): raise pybamm.OptionError( @@ -868,6 +867,10 @@ def default_var_pts(self): "z": 10, "R_n": 30, "R_p": 30, + "R_n_prim": 30, + "R_p_prim": 30, + "R_n_sec": 30, + "R_p_sec": 30, } # Reduce the default points for 2D current collectors if self.options["dimensionality"] == 2: @@ -888,6 +891,10 @@ def default_submesh_types(self): "positive secondary particle": pybamm.Uniform1DSubMesh, "negative particle size": pybamm.Uniform1DSubMesh, "positive particle size": pybamm.Uniform1DSubMesh, + "negative primary particle size": pybamm.Uniform1DSubMesh, + "positive primary particle size": pybamm.Uniform1DSubMesh, + "negative secondary particle size": pybamm.Uniform1DSubMesh, + "positive secondary particle size": pybamm.Uniform1DSubMesh, } if self.options["dimensionality"] == 0: base_submeshes["current collector"] = pybamm.SubMesh0D @@ -910,6 +917,10 @@ def default_spatial_methods(self): "positive secondary particle": pybamm.FiniteVolume(), "negative particle size": pybamm.FiniteVolume(), "positive particle size": pybamm.FiniteVolume(), + "negative primary particle size": pybamm.FiniteVolume(), + "positive primary particle size": pybamm.FiniteVolume(), + "negative secondary particle size": pybamm.FiniteVolume(), + "positive secondary particle size": pybamm.FiniteVolume(), } if self.options["dimensionality"] == 0: # 0D submesh - use base spatial method diff --git a/src/pybamm/models/submodels/active_material/base_active_material.py b/src/pybamm/models/submodels/active_material/base_active_material.py index ea0e826e09..4827c1f6d4 100644 --- a/src/pybamm/models/submodels/active_material/base_active_material.py +++ b/src/pybamm/models/submodels/active_material/base_active_material.py @@ -81,9 +81,19 @@ def _get_standard_active_material_variables(self, eps_solid): R = self.phase_param.R elif domain_options["particle size"] == "distribution": if self.domain == "negative": - R_ = pybamm.standard_spatial_vars.R_n + if self.phase_name == "primary ": + R_ = pybamm.standard_spatial_vars.R_n_prim + elif self.phase_name == "secondary ": + R_ = pybamm.standard_spatial_vars.R_n_sec + else: + R_ = pybamm.standard_spatial_vars.R_n elif self.domain == "positive": - R_ = pybamm.standard_spatial_vars.R_p + if self.phase_name == "primary ": + R_ = pybamm.standard_spatial_vars.R_p_prim + elif self.phase_name == "secondary ": + R_ = pybamm.standard_spatial_vars.R_p_sec + else: + R_ = pybamm.standard_spatial_vars.R_p R = pybamm.size_average(R_) R_av = pybamm.x_average(R) diff --git a/src/pybamm/models/submodels/interface/base_interface.py b/src/pybamm/models/submodels/interface/base_interface.py index 0ad08d5454..83f7c81b43 100644 --- a/src/pybamm/models/submodels/interface/base_interface.py +++ b/src/pybamm/models/submodels/interface/base_interface.py @@ -46,6 +46,8 @@ def __init__(self, param, domain, reaction, options, phase="primary"): self.reaction_name = self.phase_name + self.reaction_name self.reaction = reaction + domain_options = getattr(self.options, domain) + self.size_distribution = domain_options["particle size"] == "distribution" def _get_exchange_current_density(self, variables): """ @@ -93,8 +95,10 @@ def _get_exchange_current_density(self, variables): # "current collector" c_e = pybamm.PrimaryBroadcast(c_e, ["current collector"]) # broadcast c_e, T onto "particle size" - c_e = pybamm.PrimaryBroadcast(c_e, [f"{domain} particle size"]) - T = pybamm.PrimaryBroadcast(T, [f"{domain} particle size"]) + c_e = pybamm.PrimaryBroadcast( + c_e, [f"{domain} {phase_name}particle size"] + ) + T = pybamm.PrimaryBroadcast(T, [f"{domain} {phase_name}particle size"]) else: c_s_surf = variables[ @@ -215,13 +219,23 @@ def _get_standard_interfacial_current_variables(self, j): # Size average. For j variables that depend on particle size, see # "_get_standard_size_distribution_interfacial_current_variables" - if j.domain in [["negative particle size"], ["positive particle size"]]: + if j.domain in [ + ["negative particle size"], + ["positive particle size"], + ["negative primary particle size"], + ["positive primary particle size"], + ["negative secondary particle size"], + ["positive secondary particle size"], + ]: j = pybamm.size_average(j) # Average, and broadcast if necessary j_av = pybamm.x_average(j) if j.domain == []: j = pybamm.FullBroadcast(j, f"{domain} electrode", "current collector") - elif j.domain == ["current collector"]: + elif ( + j.domain == ["current collector"] + and getattr(self, "reaction_loc", None) != "interface" + ): j = pybamm.PrimaryBroadcast(j, f"{domain} electrode") variables = { @@ -230,6 +244,13 @@ def _get_standard_interfacial_current_variables(self, j): f"X-averaged {domain} electrode {reaction_name}" "interfacial current density [A.m-2]": j_av, } + if self.phase_name == reaction_name: + variables.update( + { + f"{Domain} electrode {reaction_name}interfacial current density [A.m-2]": j, + f"X-averaged {domain} electrode {reaction_name}interfacial current density [A.m-2]": j_av, + } + ) return variables @@ -298,7 +319,6 @@ def _get_standard_volumetric_current_density_variables(self, variables): a = variables[ f"{Domain} electrode {phase_name}surface area to volume ratio [m-1]" ] - j = variables[ f"{Domain} electrode {reaction_name}interfacial current density [A.m-2]" ] @@ -335,7 +355,14 @@ def _get_standard_overpotential_variables(self, eta_r): # Size average. For eta_r variables that depend on particle size, see # "_get_standard_size_distribution_overpotential_variables" - if eta_r.domain in [["negative particle size"], ["positive particle size"]]: + if eta_r.domain in [ + ["negative particle size"], + ["positive particle size"], + ["negative primary particle size"], + ["positive primary particle size"], + ["negative secondary particle size"], + ["positive secondary particle size"], + ]: eta_r = pybamm.size_average(eta_r) # X-average, and broadcast if necessary @@ -410,9 +437,11 @@ def _get_standard_size_distribution_interfacial_current_variables(self, j): if j.domains["secondary"] == [f"{domain} electrode"]: # x-average j_xav = pybamm.x_average(j) - else: + elif getattr(self, "reaction_loc", None) != "interface": j_xav = j j = pybamm.SecondaryBroadcast(j_xav, [f"{domain} electrode"]) + else: + j_xav = j variables = { f"{Domain} electrode {reaction_name}" diff --git a/src/pybamm/models/submodels/interface/kinetics/base_kinetics.py b/src/pybamm/models/submodels/interface/kinetics/base_kinetics.py index 7cfa83631b..c3b6ed5329 100644 --- a/src/pybamm/models/submodels/interface/kinetics/base_kinetics.py +++ b/src/pybamm/models/submodels/interface/kinetics/base_kinetics.py @@ -76,7 +76,9 @@ def get_coupled_variables(self, variables): self.reaction == "lithium-ion main" and domain_options["particle size"] == "distribution" ): - delta_phi = pybamm.PrimaryBroadcast(delta_phi, [f"{domain} particle size"]) + delta_phi = pybamm.PrimaryBroadcast( + delta_phi, [f"{domain} {phase_name}particle size"] + ) # Get exchange-current density. For MSMR models we calculate the exchange # current density for each reaction, then sum these to give a total exchange @@ -169,7 +171,7 @@ def get_coupled_variables(self, variables): elif j0.domain in ["current collector", ["current collector"]]: T = variables["X-averaged cell temperature [K]"] u = variables[f"X-averaged {domain} electrode interface utilisation"] - elif j0.domain == [f"{domain} particle size"]: + elif j0.domain == [f"{domain} {phase_name}particle size"]: if j0.domains["secondary"] != [f"{domain} electrode"]: T = variables["X-averaged cell temperature [K]"] u = variables[f"X-averaged {domain} electrode interface utilisation"] @@ -178,7 +180,7 @@ def get_coupled_variables(self, variables): u = variables[f"{Domain} electrode interface utilisation"] # Broadcast T onto "particle size" domain - T = pybamm.PrimaryBroadcast(T, [f"{domain} particle size"]) + T = pybamm.PrimaryBroadcast(T, [f"{domain} {phase_name}particle size"]) else: T = variables[f"{Domain} electrode temperature [K]"] u = variables[f"{Domain} electrode interface utilisation"] @@ -198,7 +200,9 @@ def get_coupled_variables(self, variables): else: j = self._get_kinetics(j0, ne, eta_r, T, u) - if j.domain == [f"{domain} particle size"]: + if j.domain == [f"{domain} particle size"] or j.domain == [ + f"{domain} {phase_name}particle size" + ]: # If j depends on particle size, get size-dependent "distribution" # variables first variables.update( diff --git a/src/pybamm/models/submodels/interface/lithium_plating/base_plating.py b/src/pybamm/models/submodels/interface/lithium_plating/base_plating.py index 0530476f9d..c34eca9e05 100644 --- a/src/pybamm/models/submodels/interface/lithium_plating/base_plating.py +++ b/src/pybamm/models/submodels/interface/lithium_plating/base_plating.py @@ -64,7 +64,9 @@ def _get_standard_concentration_variables(self, c_plated_Li, c_dead_Li): phase_name = self.phase_name phase_param = self.phase_param domain, Domain = self.domain_Domain - + if self.size_distribution is True: + c_plated_Li = pybamm.size_average(c_plated_Li) + c_dead_Li = pybamm.size_average(c_dead_Li) # Set scales to one for the "no plating" model so that they are not required # by parameter values in general if isinstance(self, pybamm.lithium_plating.NoPlating): @@ -137,3 +139,24 @@ def _get_standard_reaction_variables(self, j_stripping): } return variables + + def _get_standard_size_distribution_reaction_variables(self, j_stripping): + """ + A private function to obtain the standard variables which + can be derived from the lithum stripping interfacial reaction current + """ + domain, Domain = self.domain_Domain + phase_name = self.phase_name + j_stripping_sav = pybamm.size_average(j_stripping) + j_stripping_av = pybamm.x_average(j_stripping_sav) + + variables = { + f"{Domain} electrode {phase_name}lithium plating " + "interfacial current density distribution [A.m-2]": j_stripping, + f"{Domain} electrode {phase_name}lithium plating " + "interfacial current density [A.m-2]": j_stripping_sav, + f"X-averaged {domain} electrode {phase_name}lithium plating " + "interfacial current density [A.m-2]": j_stripping_av, + } + + return variables diff --git a/src/pybamm/models/submodels/interface/lithium_plating/no_plating.py b/src/pybamm/models/submodels/interface/lithium_plating/no_plating.py index ddc13cfb97..8efa664a92 100644 --- a/src/pybamm/models/submodels/interface/lithium_plating/no_plating.py +++ b/src/pybamm/models/submodels/interface/lithium_plating/no_plating.py @@ -20,12 +20,32 @@ def __init__(self, param, domain, options=None, phase="primary"): super().__init__(param, domain, options=options, phase=phase) def get_fundamental_variables(self): - zero = pybamm.FullBroadcast( - pybamm.Scalar(0), f"{self.domain} electrode", "current collector" - ) + phase_name = self.phase_name + if self.size_distribution is False: + zero = pybamm.FullBroadcast( + pybamm.Scalar(0), f"{self.domain} electrode", "current collector" + ) + else: + zero = pybamm.FullBroadcast( + pybamm.Scalar(0), + f"{self.domain} {phase_name}particle size", + { + "secondary": f"{self.domain} electrode", + "tertiary": "current collector", + }, + ) variables = self._get_standard_concentration_variables(zero, zero) + if self.size_distribution: + variables.update( + self._get_standard_size_distribution_overpotential_variables(zero) + ) + variables.update( + self._get_standard_size_distribution_reaction_variables(zero) + ) + else: + variables.update(self._get_standard_reaction_variables(zero)) variables.update(self._get_standard_overpotential_variables(zero)) - variables.update(self._get_standard_reaction_variables(zero)) + return variables def get_coupled_variables(self, variables): diff --git a/src/pybamm/models/submodels/interface/open_circuit_potential/base_ocp.py b/src/pybamm/models/submodels/interface/open_circuit_potential/base_ocp.py index 309eb343c0..e2c4a43551 100644 --- a/src/pybamm/models/submodels/interface/open_circuit_potential/base_ocp.py +++ b/src/pybamm/models/submodels/interface/open_circuit_potential/base_ocp.py @@ -49,7 +49,14 @@ def _get_standard_ocp_variables(self, ocp_surf, ocp_bulk, dUdT): reaction_name = self.reaction_name # Update size variables then size average. - if ocp_surf.domain in [["negative particle size"], ["positive particle size"]]: + if ocp_surf.domain in [ + ["negative particle size"], + ["positive particle size"], + ["negative primary particle size"], + ["positive primary particle size"], + ["negative secondary particle size"], + ["positive secondary particle size"], + ]: variables = self._get_standard_size_distribution_ocp_variables( ocp_surf, dUdT ) diff --git a/src/pybamm/models/submodels/interface/open_circuit_potential/current_sigmoid_ocp.py b/src/pybamm/models/submodels/interface/open_circuit_potential/current_sigmoid_ocp.py index fa55ba1079..8fc41d1711 100644 --- a/src/pybamm/models/submodels/interface/open_circuit_potential/current_sigmoid_ocp.py +++ b/src/pybamm/models/submodels/interface/open_circuit_potential/current_sigmoid_ocp.py @@ -36,7 +36,7 @@ def get_coupled_variables(self, variables): ): sto_surf = sto_surf.orphans[0] T = T.orphans[0] - T = pybamm.PrimaryBroadcast(T, [f"{domain} particle size"]) + T = pybamm.PrimaryBroadcast(T, [f"{domain} {phase_name}particle size"]) else: sto_surf = variables[ f"{Domain} {phase_name}particle surface stoichiometry" diff --git a/src/pybamm/models/submodels/interface/open_circuit_potential/single_ocp.py b/src/pybamm/models/submodels/interface/open_circuit_potential/single_ocp.py index 932e2e3eff..42a4b53ebb 100644 --- a/src/pybamm/models/submodels/interface/open_circuit_potential/single_ocp.py +++ b/src/pybamm/models/submodels/interface/open_circuit_potential/single_ocp.py @@ -25,7 +25,7 @@ def get_coupled_variables(self, variables): ): sto_surf = sto_surf.orphans[0] T = T.orphans[0] - T = pybamm.PrimaryBroadcast(T, [f"{domain} particle size"]) + T = pybamm.PrimaryBroadcast(T, [f"{domain} {phase_name}particle size"]) else: sto_surf = variables[ f"{Domain} {phase_name}particle surface stoichiometry" @@ -42,7 +42,7 @@ def get_coupled_variables(self, variables): # Bulk OCP is from the average SOC and temperature sto_bulk = variables[f"{Domain} electrode {phase_name}stoichiometry"] - T_bulk = pybamm.xyz_average(pybamm.size_average(T)) + T_bulk = pybamm.xyzs_average(T) ocp_bulk = self.phase_param.U(sto_bulk, T_bulk) elif self.reaction == "lithium metal plating": T = variables[f"{Domain} electrode temperature [K]"] diff --git a/src/pybamm/models/submodels/interface/sei/base_sei.py b/src/pybamm/models/submodels/interface/sei/base_sei.py index 8a4176d416..af5a929bb5 100644 --- a/src/pybamm/models/submodels/interface/sei/base_sei.py +++ b/src/pybamm/models/submodels/interface/sei/base_sei.py @@ -152,7 +152,11 @@ def _get_standard_concentration_variables(self, variables): L = variables[f"{Domain} {reaction_name}thickness [m]"] n_SEI = L * L_to_n - n_SEI_xav = pybamm.x_average(n_SEI) + if self.size_distribution: + n_SEI_sav = pybamm.size_average(n_SEI) + n_SEI_xav = pybamm.x_average(n_SEI_sav) + else: + n_SEI_xav = pybamm.x_average(n_SEI) n_SEI_av = pybamm.yz_average(n_SEI_xav) # Calculate change in SEI concentration with respect to initial state @@ -188,8 +192,13 @@ def _get_standard_concentration_variables(self, variables): roughness = variables[f"{Domain} electrode roughness ratio"] n_SEI_cr = L_cr * L_to_n * (roughness - 1) # SEI on cracks concentration - n_SEI_cr_xav = pybamm.x_average(n_SEI_cr) - n_SEI_cr_av = pybamm.yz_average(n_SEI_cr_xav) + if self.size_distribution: + n_SEI_cr_sav = pybamm.size_average(n_SEI_cr) + n_SEI_cr_xav = pybamm.x_average(n_SEI_cr_sav) + n_SEI_cr_av = pybamm.yz_average(n_SEI_cr_xav) + else: + n_SEI_cr_xav = pybamm.x_average(n_SEI_cr) + n_SEI_cr_av = pybamm.yz_average(n_SEI_cr_xav) # Calculate change in SEI cracks concentration # Initial state depends on roughness (to avoid division by zero) @@ -253,3 +262,36 @@ def _get_standard_reaction_variables(self, j_sei): ) return variables + + def _get_standard_reaction_distribution_variables(self, j_sei): + """ + A private function to obtain the standard variables which + can be derived from the SEI interfacial reaction current + + Parameters + ---------- + j_sei : :class:`pybamm.Symbol` + The SEI interfacial reaction current. + + Returns + ------- + variables : dict + The variables which can be derived from the SEI currents. + """ + domain, Domain = self.domain_Domain + j_sei = pybamm.size_average(j_sei) + variables = { + f"{Domain} electrode {self.reaction_name}" + "interfacial current density [A.m-2]": j_sei, + } + + if self.reaction_loc != "interface": + j_sei_av = pybamm.x_average(j_sei) + variables.update( + { + f"X-averaged {domain} electrode {self.reaction_name}" + "interfacial current density [A.m-2]": j_sei_av, + } + ) + + return variables diff --git a/src/pybamm/models/submodels/interface/sei/no_sei.py b/src/pybamm/models/submodels/interface/sei/no_sei.py index 27cea3b69b..10192a485f 100644 --- a/src/pybamm/models/submodels/interface/sei/no_sei.py +++ b/src/pybamm/models/submodels/interface/sei/no_sei.py @@ -32,12 +32,26 @@ def get_fundamental_variables(self): domain = self.domain.lower() if self.reaction_loc == "interface": zero = pybamm.PrimaryBroadcast(pybamm.Scalar(0), "current collector") + elif self.size_distribution: + zero = pybamm.FullBroadcast( + pybamm.Scalar(0), + f"{domain} {self.phase_name}particle size", + {"secondary": f"{domain} electrode", "tertiary": "current collector"}, + ) else: zero = pybamm.FullBroadcast( pybamm.Scalar(0), f"{domain} electrode", "current collector" ) variables = self._get_standard_thickness_variables(zero) - variables.update(self._get_standard_reaction_variables(zero)) + + if self.size_distribution: + variables.update(self._get_standard_reaction_distribution_variables(zero)) + variables.update( + self._get_standard_size_distribution_interfacial_current_variables(zero) + ) + else: + variables.update(self._get_standard_reaction_variables(zero)) + variables.update(self._get_standard_interfacial_current_variables(zero)) return variables def get_coupled_variables(self, variables): diff --git a/src/pybamm/models/submodels/interface/total_interfacial_current.py b/src/pybamm/models/submodels/interface/total_interfacial_current.py index 79e13e37f6..3df4ff757a 100644 --- a/src/pybamm/models/submodels/interface/total_interfacial_current.py +++ b/src/pybamm/models/submodels/interface/total_interfacial_current.py @@ -165,7 +165,8 @@ def _get_whole_cell_coupled_variables(self, variables): if domain == "separator": var_dict[domain] = zero_s else: - var_dict[domain] = variables[variable_template.format(domain + " ")] + var = variables[variable_template.format(domain + " ")] + var_dict[domain] = var var = pybamm.concatenation(*var_dict.values()) variables.update({variable_template.format(""): var}) diff --git a/src/pybamm/models/submodels/particle/base_particle.py b/src/pybamm/models/submodels/particle/base_particle.py index b774e58a0c..b5b8fb4499 100644 --- a/src/pybamm/models/submodels/particle/base_particle.py +++ b/src/pybamm/models/submodels/particle/base_particle.py @@ -185,6 +185,7 @@ def _get_distribution_variables(self, R): """ domain, Domain = self.domain_Domain phase_name = self.phase_name + Phase_prefactor = self.phase_param.phase_prefactor R_typ = self.phase_param.R_typ # [m] # Particle-size distribution (area-weighted) @@ -237,21 +238,21 @@ def _get_distribution_variables(self, R): variables = { f"{Domain} {phase_name}particle sizes": R / R_typ, - f"{Domain} {phase_name}particle sizes [m]": R, - f"{Domain} area-weighted {phase_name}particle-size" + f"{Phase_prefactor}{Domain} {phase_name}particle sizes [m]": R, + f"{Phase_prefactor}{Domain} area-weighted {phase_name}particle-size" " distribution [m-1]": f_a_dist, - f"{Domain} volume-weighted {phase_name}particle-size" + f"{Phase_prefactor}{Domain} volume-weighted {phase_name}particle-size" " distribution [m-1]": f_v_dist, - f"{Domain} number-based {phase_name}particle-size" + f"{Phase_prefactor}{Domain} number-based {phase_name}particle-size" " distribution [m-1]": f_num_dist, - f"{Domain} area-weighted mean particle radius [m]": R_a_mean, - f"{Domain} volume-weighted mean particle radius [m]": R_v_mean, - f"{Domain} number-based mean particle radius [m]": R_num_mean, - f"{Domain} area-weighted {phase_name}particle-size" + f"{Phase_prefactor}{Domain} area-weighted mean particle radius [m]": R_a_mean, + f"{Phase_prefactor}{Domain} volume-weighted mean particle radius [m]": R_v_mean, + f"{Phase_prefactor}{Domain} number-based mean particle radius [m]": R_num_mean, + f"{Phase_prefactor}{Domain} area-weighted {phase_name}particle-size" " standard deviation [m]": sd_a, - f"{Domain} volume-weighted {phase_name}particle-size" + f"{Phase_prefactor}{Domain} volume-weighted {phase_name}particle-size" " standard deviation [m]": sd_v, - f"{Domain} number-based {phase_name}particle-size" + f"{Phase_prefactor}{Domain} number-based {phase_name}particle-size" " standard deviation [m]": sd_num, # X-averaged sizes and distributions f"X-averaged {domain} {phase_name}particle sizes [m]": pybamm.x_average(R), diff --git a/src/pybamm/models/submodels/particle/fickian_diffusion.py b/src/pybamm/models/submodels/particle/fickian_diffusion.py index 634e2ce730..5ae5d3a4da 100644 --- a/src/pybamm/models/submodels/particle/fickian_diffusion.py +++ b/src/pybamm/models/submodels/particle/fickian_diffusion.py @@ -31,7 +31,7 @@ def __init__(self, param, domain, options, phase="primary", x_average=False): def get_fundamental_variables(self): domain, Domain = self.domain_Domain phase_name = self.phase_name - + Phase_prefactor = self.phase_param.phase_prefactor variables = {} if self.size_distribution is False: if self.x_average is False: @@ -81,7 +81,7 @@ def get_fundamental_variables(self): ) variables = self._get_distribution_variables(R) f_v_dist = variables[ - f"{Domain} volume-weighted {phase_name}" + f"{Phase_prefactor}{Domain} volume-weighted {phase_name}" "particle-size distribution [m-1]" ] else: @@ -130,6 +130,7 @@ def get_fundamental_variables(self): def get_coupled_variables(self, variables): domain, Domain = self.domain_Domain phase_name = self.phase_name + Phase_prefactor = self.phase_param.phase_prefactor if self.size_distribution is False: if self.x_average is False: @@ -219,7 +220,7 @@ def get_coupled_variables(self, variables): ) variables.update(self._get_standard_flux_distribution_variables(N_s)) # Size-averaged flux variables - R = variables[f"{Domain} {phase_name}particle sizes [m]"] + R = variables[f"{Phase_prefactor}{Domain} {phase_name}particle sizes [m]"] f_a_dist = self.phase_param.f_a_dist(R) D_eff = pybamm.Integral(f_a_dist * D_eff, R) N_s = pybamm.Integral(f_a_dist * N_s, R) diff --git a/src/pybamm/parameters/geometric_parameters.py b/src/pybamm/parameters/geometric_parameters.py index 0155bbe03c..2998e97318 100644 --- a/src/pybamm/parameters/geometric_parameters.py +++ b/src/pybamm/parameters/geometric_parameters.py @@ -165,6 +165,14 @@ def f_a_dist(self, R): Dimensional electrode area-weighted particle-size distribution """ Domain = self.domain.capitalize() + if len(R.domain) == 0: + self.phase_prefactor = "" + elif "primary" in R.domains["primary"][0]: + self.phase_prefactor = "Primary: " + elif "secondary" in R.domains["primary"][0]: + self.phase_prefactor = "Secondary: " + else: + self.phase_prefactor = "" inputs = {f"{self.phase_prefactor}{Domain} particle-size variable [m]": R} return pybamm.FunctionParameter( f"{self.phase_prefactor}{Domain} " diff --git a/src/pybamm/parameters/size_distribution_parameters.py b/src/pybamm/parameters/size_distribution_parameters.py index 60383fff50..6b7087b8c8 100644 --- a/src/pybamm/parameters/size_distribution_parameters.py +++ b/src/pybamm/parameters/size_distribution_parameters.py @@ -9,14 +9,31 @@ def get_size_distribution_parameters( param, R_n_av=None, + R_n_av_prim=None, + R_n_av_sec=None, R_p_av=None, + R_p_av_prim=None, + R_p_av_sec=None, sd_n=0.3, + sd_n_prim=0.3, + sd_n_sec=0.3, sd_p=0.3, + sd_p_prim=0.3, + sd_p_sec=0.3, R_min_n=None, + R_min_n_prim=None, + R_min_n_sec=None, R_min_p=None, + R_min_p_prim=None, + R_min_p_sec=None, R_max_n=None, + R_max_n_prim=None, + R_max_n_sec=None, R_max_p=None, + R_max_p_prim=None, + R_max_p_sec=None, working_electrode="both", + composite="neither", ): """ A convenience method to add standard area-weighted particle-size distribution @@ -62,54 +79,132 @@ def get_size_distribution_parameters( """ if working_electrode == "both": # Radii from given parameter set - R_n_typ = param["Negative particle radius [m]"] - + if composite == "negative" or composite == "both": + R_n_typ_prim = param["Primary: Negative particle radius [m]"] + R_n_typ_sec = param["Secondary: Negative particle radius [m]"] + # Set the mean particle radii + R_n_av_prim = R_n_av_prim or R_n_typ_prim + R_n_av_sec = R_n_av_sec or R_n_typ_sec + # Minimum radii + R_min_n_prim = R_min_n_prim or np.max([0, 1 - sd_n_prim * 5]) + R_min_n_sec = R_min_n_sec or np.max([0, 1 - sd_n_sec * 5]) + + # Max radii + R_max_n_prim = R_max_n_prim or (1 + sd_n_prim * 5) + R_max_n_sec = R_max_n_sec or (1 + sd_n_sec * 5) + + # Area-weighted particle-size distribution + def f_a_dist_n_prim(R): + return lognormal(R, R_n_av_prim, sd_n_prim * R_n_av_prim) + + def f_a_dist_n_sec(R): + return lognormal(R, R_n_av_sec, sd_n_sec * R_n_av_sec) + + param.update( + { + "Primary: Negative minimum particle radius [m]": R_min_n_prim + * R_n_av_prim, + "Primary: Negative maximum particle radius [m]": R_max_n_prim + * R_n_av_prim, + "Primary: Negative area-weighted " + + "particle-size distribution [m-1]": f_a_dist_n_prim, + "Secondary: Negative minimum particle radius [m]": R_min_n_sec + * R_n_av_sec, + "Secondary: Negative maximum particle radius [m]": R_max_n_sec + * R_n_av_sec, + "Secondary: Negative area-weighted " + + "particle-size distribution [m-1]": f_a_dist_n_sec, + }, + check_already_exists=False, + ) + else: + R_n_typ = param["Negative particle radius [m]"] + # Set the mean particle radii + R_n_av = R_n_av or R_n_typ + # Minimum radii + R_min_n = R_min_n or np.max([0, 1 - sd_n * 5]) + + # Max radii + R_max_n = R_max_n or (1 + sd_n * 5) + + # Area-weighted particle-size distribution + def f_a_dist_n(R): + return lognormal(R, R_n_av, sd_n * R_n_av) + + param.update( + { + "Negative minimum particle radius [m]": R_min_n * R_n_av, + "Negative maximum particle radius [m]": R_max_n * R_n_av, + "Negative area-weighted " + + "particle-size distribution [m-1]": f_a_dist_n, + }, + check_already_exists=False, + ) + + if composite == "positive" or composite == "both": + R_p_typ_prim = param["Primary: Positive particle radius [m]"] + R_p_typ_sec = param["Secondary: Positive particle radius [m]"] # Set the mean particle radii - R_n_av = R_n_av or R_n_typ - + R_p_av_prim = R_p_av_prim or R_p_typ_prim + R_p_av_sec = R_p_av_sec or R_p_typ_sec # Minimum radii - R_min_n = R_min_n or np.max([0, 1 - sd_n * 5]) + R_min_p_prim = R_min_p_prim or np.max([0, 1 - sd_p_prim * 5]) + R_min_p_sec = R_min_p_sec or np.max([0, 1 - sd_p_sec * 5]) # Max radii - R_max_n = R_max_n or (1 + sd_n * 5) + R_max_p_prim = R_max_p_prim or (1 + sd_p_prim * 5) + R_max_p_sec = R_max_p_sec or (1 + sd_p_sec * 5) # Area-weighted particle-size distribution - def f_a_dist_n(R): - return lognormal(R, R_n_av, sd_n * R_n_av) + def f_a_dist_p_prim(R): + return lognormal(R, R_p_av_prim, sd_p_prim * R_p_av_prim) + + def f_a_dist_p_sec(R): + return lognormal(R, R_p_av_sec, sd_p_sec * R_p_av_sec) param.update( { - "Negative minimum particle radius [m]": R_min_n * R_n_av, - "Negative maximum particle radius [m]": R_max_n * R_n_av, - "Negative area-weighted " - + "particle-size distribution [m-1]": f_a_dist_n, + "Primary: Positive minimum particle radius [m]": R_min_p_prim + * R_p_av_prim, + "Primary: Positive maximum particle radius [m]": R_max_p_prim + * R_p_av_prim, + "Primary: Positive area-weighted " + + "particle-size distribution [m-1]": f_a_dist_p_prim, + "Secondary: Positive minimum particle radius [m]": R_min_p_sec + * R_p_av_sec, + "Secondary: Positive maximum particle radius [m]": R_max_p_sec + * R_p_av_sec, + "Secondary: Positive area-weighted " + + "particle-size distribution [m-1]": f_a_dist_p_sec, }, check_already_exists=False, ) - # Radii from given parameter set - R_p_typ = param["Positive particle radius [m]"] + else: + # Radii from given parameter set + R_p_typ = param["Positive particle radius [m]"] - # Set the mean particle radii - R_p_av = R_p_av or R_p_typ + # Set the mean particle radii + R_p_av = R_p_av or R_p_typ - # Minimum radii - R_min_p = R_min_p or np.max([0, 1 - sd_p * 5]) + # Minimum radii + R_min_p = R_min_p or np.max([0, 1 - sd_p * 5]) - # Max radii - R_max_p = R_max_p or (1 + sd_p * 5) + # Max radii + R_max_p = R_max_p or (1 + sd_p * 5) - # Area-weighted particle-size distribution - def f_a_dist_p(R): - return lognormal(R, R_p_av, sd_p * R_p_av) + # Area-weighted particle-size distribution + def f_a_dist_p(R): + return lognormal(R, R_p_av, sd_p * R_p_av) - param.update( - { - "Positive minimum particle radius [m]": R_min_p * R_p_av, - "Positive maximum particle radius [m]": R_max_p * R_p_av, - "Positive area-weighted " + "particle-size distribution [m-1]": f_a_dist_p, - }, - check_already_exists=False, - ) + param.update( + { + "Positive minimum particle radius [m]": R_min_p * R_p_av, + "Positive maximum particle radius [m]": R_max_p * R_p_av, + "Positive area-weighted " + + "particle-size distribution [m-1]": f_a_dist_p, + }, + check_already_exists=False, + ) return param diff --git a/tests/integration/test_models/standard_output_tests.py b/tests/integration/test_models/standard_output_tests.py index ca5e27607d..4e6788724f 100644 --- a/tests/integration/test_models/standard_output_tests.py +++ b/tests/integration/test_models/standard_output_tests.py @@ -322,17 +322,17 @@ def __init__(self, model, param, disc, solution, operating_condition): # Take only the x-averaged of these for now, since variables cannot have # 4 domains yet self.c_s_n_dist = solution[ - "X-averaged negative particle concentration distribution" + f"X-averaged negative {self.phase_name_n}particle concentration distribution" ] self.c_s_p_dist = solution[ - "X-averaged positive particle concentration distribution" + f"X-averaged positive {self.phase_name_p}particle concentration distribution" ] self.c_s_n_surf_dist = solution[ - "Negative particle surface concentration distribution" + f"Negative {self.phase_name_n}particle surface concentration distribution" ] self.c_s_p_surf_dist = solution[ - "Positive particle surface concentration distribution" + f"Positive {self.phase_name_p}particle surface concentration distribution" ] def test_concentration_increase_decrease(self): diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py index 5ede55ff6d..7a728a744b 100644 --- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py +++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py @@ -52,6 +52,10 @@ def setup(self): "R_p": 3, "y": 5, "z": 5, + "R_n_prim": 3, + "R_n_sec": 3, + "R_p_prim": 3, + "R_p_sec": 3, } def test_basic_processing(self): @@ -62,6 +66,31 @@ def test_basic_processing(self): ) modeltest.test_all() + def test_composite(self): + options = { + "particle phases": ("2", "1"), + "open-circuit potential": (("single", "current sigmoid"), "single"), + "particle size": "distribution", + } + parameter_values = pybamm.ParameterValues("Chen2020_composite") + name = "Negative electrode active material volume fraction" + x = 0.1 + parameter_values.update( + {f"Primary: {name}": (1 - x) * 0.75, f"Secondary: {name}": x * 0.75} + ) + parameter_values = pybamm.get_size_distribution_parameters( + parameter_values, + composite="negative", + R_min_n_prim=0.9, + R_min_n_sec=0.9, + R_max_n_prim=1.1, + R_max_n_sec=1.1, + ) + # self.run_basic_processing_test(options, parameter_values=parameter_values) + model = pybamm.lithium_ion.DFN(options) + modeltest = tests.StandardModelTest(model, parameter_values=parameter_values) + modeltest.test_all() + def test_basic_processing_tuple(self): options = {"particle size": ("single", "distribution")} model = pybamm.lithium_ion.DFN(options) diff --git a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_newman_tobias.py b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_newman_tobias.py index fb091ef56d..c24a063f70 100644 --- a/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_newman_tobias.py +++ b/tests/integration/test_models/test_full_battery_models/test_lithium_ion/test_newman_tobias.py @@ -31,6 +31,9 @@ def test_composite_graphite_silicon_sei(self): def test_composite_reaction_driven_LAM(self): pass # skip this test + def test_particle_size_composite(self): + pass # skip this test + def test_composite_stress_driven_LAM(self): pass # skip this test diff --git a/tests/unit/test_citations.py b/tests/unit/test_citations.py index 9b1c9fa2d3..dccef9b51c 100644 --- a/tests/unit/test_citations.py +++ b/tests/unit/test_citations.py @@ -357,13 +357,13 @@ def test_sripad_2020(self): citations._reset() assert "Sripad2020" not in citations._papers_to_cite - pybamm.kinetics.Marcus(None, None, None, None, None) + pybamm.kinetics.Marcus(None, "negative", None, None, None) assert "Sripad2020" in citations._papers_to_cite assert "Sripad2020" in citations._citation_tags.keys() citations._reset() assert "Sripad2020" not in citations._papers_to_cite - pybamm.kinetics.MarcusHushChidsey(None, None, None, None, None) + pybamm.kinetics.MarcusHushChidsey(None, "negative", None, None, None) assert "Sripad2020" in citations._papers_to_cite assert "Sripad2020" in citations._citation_tags.keys() diff --git a/tests/unit/test_expression_tree/test_averages.py b/tests/unit/test_expression_tree/test_averages.py index 5f9736693e..9065b85ca9 100644 --- a/tests/unit/test_expression_tree/test_averages.py +++ b/tests/unit/test_expression_tree/test_averages.py @@ -201,7 +201,14 @@ def test_size_average(self): ) assert average_a == a - for domain in [["negative particle size"], ["positive particle size"]]: + for domain in [ + ["negative particle size"], + ["positive particle size"], + ["negative primary particle size"], + ["positive primary particle size"], + ["negative secondary particle size"], + ["positive secondary particle size"], + ]: a = pybamm.Symbol("a", domain=domain) R = pybamm.SpatialVariable("R", domain) av_a = pybamm.size_average(a) diff --git a/tests/unit/test_geometry/test_battery_geometry.py b/tests/unit/test_geometry/test_battery_geometry.py index c73ef89af1..30d417fcbd 100644 --- a/tests/unit/test_geometry/test_battery_geometry.py +++ b/tests/unit/test_geometry/test_battery_geometry.py @@ -77,6 +77,16 @@ def test_geometry(self): geometry["positive secondary particle"]["r_p_sec"]["max"] == geo.p.sec.R_typ ) + options = { + "particle phases": "2", + "particle size": "distribution", + } + geometry = pybamm.battery_geometry(options=options) + assert "negative primary particle size" in geometry + assert "positive primary particle size" in geometry + assert "negative secondary particle size" in geometry + assert "positive secondary particle size" in geometry + def test_geometry_error(self): with pytest.raises(pybamm.GeometryError, match="Invalid current"): pybamm.battery_geometry( diff --git a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py index be7ee861c4..554b65d880 100644 --- a/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py +++ b/tests/unit/test_models/test_full_battery_models/test_base_battery_model.py @@ -140,6 +140,10 @@ def test_default_var_pts(self): "z": 10, "R_n": 30, "R_p": 30, + "R_n_prim": 30, + "R_p_prim": 30, + "R_n_sec": 30, + "R_p_sec": 30, } model = pybamm.BaseBatteryModel({"dimensionality": 0}) assert var_pts == model.default_var_pts diff --git a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py index 973a0f348b..56f4c76f44 100644 --- a/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py +++ b/tests/unit/test_models/test_full_battery_models/test_lithium_ion/test_dfn.py @@ -27,6 +27,10 @@ def test_well_posed_size_distribution_tuple(self): options = {"particle size": ("single", "distribution")} self.check_well_posedness(options) + def test_well_posed_size_distribution_composite(self): + options = {"particle size": "distribution", "particle phases": "2"} + self.check_well_posedness(options) + def test_well_posed_current_sigmoid_ocp_with_psd(self): options = { "open-circuit potential": "current sigmoid", diff --git a/tests/unit/test_parameters/test_parameter_values.py b/tests/unit/test_parameters/test_parameter_values.py index cdcfb30ede..0e32f60476 100644 --- a/tests/unit/test_parameters/test_parameter_values.py +++ b/tests/unit/test_parameters/test_parameter_values.py @@ -893,7 +893,7 @@ def dist(R): var_av_proc = param.process_symbol(var_av) assert isinstance(var_av_proc, pybamm.SizeAverage) - R = pybamm.SpatialVariable("R", "negative particle size") + R = pybamm.SpatialVariable("R_n", "negative particle size") assert var_av_proc.f_a_dist == R**2 def test_process_not_constant(self): diff --git a/tests/unit/test_parameters/test_size_distribution_parameters.py b/tests/unit/test_parameters/test_size_distribution_parameters.py index 414b422055..a0f7847a32 100644 --- a/tests/unit/test_parameters/test_size_distribution_parameters.py +++ b/tests/unit/test_parameters/test_size_distribution_parameters.py @@ -40,3 +40,76 @@ def test_parameter_values(self): R_test = pybamm.Scalar(1.0) values.evaluate(param.n.prim.f_a_dist(R_test)) values.evaluate(param.p.prim.f_a_dist(R_test)) + + # check that the composite parameters works as well + params_composite = pybamm.ParameterValues("Chen2020_composite") + params_composite.update( + { + "Primary: Positive particle radius [m]": 1e-6, + "Secondary: Positive particle radius [m]": 1e-6, + "Negative particle radius [m]": 1e-6, + }, + check_already_exists=False, + ) + values_composite = pybamm.get_size_distribution_parameters( + params_composite, + composite="negative", + ) + assert "Primary: Negative maximum particle radius [m]" in values_composite + param = pybamm.LithiumIonParameters({"particle phases": ("2", "1")}) + R_test_n_prim = pybamm.Scalar(1e-6) + R_test_n_prim.domains = {"primary": ["negative primary particle size"]} + R_test_n_sec = pybamm.Scalar(1e-6) + R_test_n_sec.domains = {"primary": ["negative secondary particle size"]} + values_composite.evaluate(param.n.prim.f_a_dist(R_test_n_prim)) + values_composite.evaluate(param.n.sec.f_a_dist(R_test_n_sec)) + params_composite = pybamm.ParameterValues("Chen2020_composite") + params_composite.update( + { + "Primary: Positive particle radius [m]": 1e-6, + "Secondary: Positive particle radius [m]": 1e-6, + "Negative particle radius [m]": 1e-6, + }, + check_already_exists=False, + ) + values_composite = pybamm.get_size_distribution_parameters( + params_composite, + composite="positive", + ) + assert "Primary: Positive maximum particle radius [m]" in values_composite + param = pybamm.LithiumIonParameters({"particle phases": ("1", "2")}) + R_test_p_prim = pybamm.Scalar(1e-6) + R_test_p_prim.domains = {"primary": ["positive primary particle size"]} + R_test_p_sec = pybamm.Scalar(1e-6) + R_test_p_sec.domains = {"primary": ["positive secondary particle size"]} + values_composite.evaluate(param.p.prim.f_a_dist(R_test_p_prim)) + values_composite.evaluate(param.p.sec.f_a_dist(R_test_p_sec)) + params_composite = pybamm.ParameterValues("Chen2020_composite") + params_composite.update( + { + "Primary: Positive particle radius [m]": 1e-6, + "Secondary: Positive particle radius [m]": 1e-6, + "Negative particle radius [m]": 1e-6, + }, + check_already_exists=False, + ) + values_composite = pybamm.get_size_distribution_parameters( + params_composite, + composite="both", + ) + assert "Primary: Negative maximum particle radius [m]" in values_composite + assert "Primary: Positive maximum particle radius [m]" in values_composite + params_composite = pybamm.ParameterValues("Chen2020_composite") + params_composite.update( + { + "Primary: Positive particle radius [m]": 1e-6, + "Secondary: Positive particle radius [m]": 1e-6, + "Negative particle radius [m]": 1e-6, + }, + check_already_exists=False, + ) + values_composite = pybamm.get_size_distribution_parameters( + params_composite, + ) + assert "Primary: Negative maximum particle radius [m]" not in values_composite + assert "Primary: Positive maximum particle radius [m]" not in values_composite