diff --git a/CHANGELOG.md b/CHANGELOG.md index f683a96db..b75e88531 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -2,6 +2,7 @@ ## Features +- [#974](https://github.com/pybop-team/PyBOP/pull/974) - Adds voltage components to each grouped model as well as asymmetric and multiphase Butler-Volmer kinetics. - [#969](https://github.com/pybop-team/PyBOP/pull/969) - Updates synthetic data and adds example script for thermal parameterisation. - [#965](https://github.com/pybop-team/PyBOP/pull/965) - Adds synthetic data and example scripts for OCV parameterisation. - [#963](https://github.com/pybop-team/PyBOP/pull/963) - Adds an example for generating synthetic data from a specification and exporting it to a PyProBE-compatible parquet file. diff --git a/examples/notebooks/battery_parameterisation/lgm50_pulse_validation.ipynb b/examples/notebooks/battery_parameterisation/lgm50_pulse_validation.ipynb index 9509c47cf..b091a7316 100644 --- a/examples/notebooks/battery_parameterisation/lgm50_pulse_validation.ipynb +++ b/examples/notebooks/battery_parameterisation/lgm50_pulse_validation.ipynb @@ -223,7 +223,7 @@ "id": "12", "metadata": {}, "source": [ - "We can construct the model and parameter values for a two-RC circuit model with initial values as listed below. Note, the \"Initial SOC\" is shifted slightly to better match the zero degree data." + "We can construct the model and parameter values for a two-RC circuit model with initial values as listed below. Note, the \"Initial SoC\" is shifted slightly to better match the zero degree data." ] }, { diff --git a/examples/scripts/battery_parameterisation/gitt_fitting.py b/examples/scripts/battery_parameterisation/gitt_fitting.py index 6abe016f2..8c067cec6 100644 --- a/examples/scripts/battery_parameterisation/gitt_fitting.py +++ b/examples/scripts/battery_parameterisation/gitt_fitting.py @@ -47,7 +47,7 @@ ) # Group the parameters -grouped_parameter_values = pybop.lithium_ion.SPDiffusion.create_grouped_parameters( +grouped_parameter_values = pybop.li_half_cell.SPDiffusion.create_grouped_parameters( parameter_values ) @@ -60,11 +60,13 @@ gitt_parameter_data = gitt_fit() # Plot the functional parameters -pybop.plot.dataset(gitt_parameter_data, signal=["Particle diffusion time scale [s]"]) +pybop.plot.dataset( + gitt_parameter_data, signal=["Positive particle diffusion time scale [s]"] +) pybop.plot.dataset(gitt_parameter_data, signal=["Series resistance [Ohm]"]) # Run the identified model -identified_model = pybop.lithium_ion.SPDiffusion(build=True) +identified_model = pybop.li_half_cell.SPDiffusion(build=True) grouped_parameter_values.update(gitt_fit.best_inputs) grouped_parameter_values["Current function [A]"] = pybamm.Interpolant( dataset["Time [s]"], dataset["Current [A]"], pybamm.t @@ -76,7 +78,7 @@ # Return to the original model and update the diffusivity value diffusivity = np.mean( parameter_values["Positive particle radius [m]"] ** 2 - / gitt_parameter_data["Particle diffusion time scale [s]"] + / gitt_parameter_data["Positive particle diffusion time scale [s]"] ) parameter_values.update({"Positive particle diffusivity [m2.s-1]": diffusivity}) diff --git a/examples/scripts/battery_parameterisation/gitt_pulse.py b/examples/scripts/battery_parameterisation/gitt_pulse.py index 1ce651e42..be48cd90d 100644 --- a/examples/scripts/battery_parameterisation/gitt_pulse.py +++ b/examples/scripts/battery_parameterisation/gitt_pulse.py @@ -30,7 +30,7 @@ ) # Group the parameters -grouped_parameter_values = pybop.lithium_ion.SPDiffusion.create_grouped_parameters( +grouped_parameter_values = pybop.li_half_cell.SPDiffusion.create_grouped_parameters( parameter_values ) diff --git a/examples/scripts/comparison_examples/gitt_models.py b/examples/scripts/comparison_examples/gitt_models.py index d5df9d907..913385a85 100644 --- a/examples/scripts/comparison_examples/gitt_models.py +++ b/examples/scripts/comparison_examples/gitt_models.py @@ -2,6 +2,7 @@ import pybamm import pybop +from pybop.models.li_half_cell import SPDiffusion, WeppnerHuggins # Define model and parameter values model_options = {"working electrode": "positive"} @@ -30,13 +31,13 @@ } ) -for model in [pybop.lithium_ion.WeppnerHuggins(), pybop.lithium_ion.SPDiffusion()]: +for model in [WeppnerHuggins(), SPDiffusion()]: # GITT target parameter diffusion_parameter = pybop.Parameter(pybop.Gaussian(5000, 1000)) - if isinstance(model, pybop.lithium_ion.WeppnerHuggins): + if isinstance(model, WeppnerHuggins): # Group parameter values - grouped_parameter_values = ( - pybop.lithium_ion.WeppnerHuggins.create_grouped_parameters(parameter_values) + grouped_parameter_values = WeppnerHuggins.create_grouped_parameters( + parameter_values ) # We can fit only the duration of the pulse @@ -49,7 +50,7 @@ dataset["Discharge capacity [A.h]"][-1] - dataset["Discharge capacity [A.h]"][0] ) - * (grouped_parameter_values["Theoretical electrode capacity [A.s]"] / 3600) + * grouped_parameter_values["Theoretical electrode capacity [A.h]"] ) grouped_parameter_values.update( { @@ -61,7 +62,7 @@ # Fitting parameters grouped_parameter_values.update( { - "Particle diffusion time scale [s]": diffusion_parameter, + "Positive particle diffusion time scale [s]": diffusion_parameter, "Reference voltage [V]": pybop.Parameter( initial_value=grouped_parameter_values["Reference voltage [V]"], ), @@ -70,14 +71,14 @@ else: # Group parameter values - grouped_parameter_values = ( - pybop.lithium_ion.SPDiffusion.create_grouped_parameters(parameter_values) + grouped_parameter_values = SPDiffusion.create_grouped_parameters( + parameter_values ) # Fitting parameters grouped_parameter_values.update( { - "Particle diffusion time scale [s]": diffusion_parameter, + "Positive particle diffusion time scale [s]": diffusion_parameter, "Series resistance [Ohm]": pybop.Parameter( initial_value=grouped_parameter_values["Series resistance [Ohm]"], ), @@ -87,7 +88,7 @@ # Build the problem gitt_dataset = ( dataset.get_subset(pulse_index) - if isinstance(model, pybop.lithium_ion.WeppnerHuggins) + if isinstance(model, WeppnerHuggins) else dataset ) simulator = pybop.pybamm.Simulator( @@ -97,13 +98,15 @@ problem = pybop.Problem(simulator, cost) # Build the optimisation problem - optim = pybop.SciPyMinimize(problem) + options = pybop.SciPyMinimizeOptions(method="Nelder-Mead") + optim = pybop.SciPyMinimize(problem, options) # Run the optimisation problem result = optim.run() print(result) print( - "Diffusion time [s]:", result.best_inputs["Particle diffusion time scale [s]"] + "Diffusion time [s]:", + result.best_inputs["Positive particle diffusion time scale [s]"], ) # Plot the timeseries output diff --git a/examples/scripts/dfn_parameterisation/3_cell_balancing.py b/examples/scripts/dfn_parameterisation/3_cell_balancing.py index c67fd8407..acd05271c 100644 --- a/examples/scripts/dfn_parameterisation/3_cell_balancing.py +++ b/examples/scripts/dfn_parameterisation/3_cell_balancing.py @@ -354,7 +354,7 @@ def solve_batch(self, inputs, calculate_sensitivities: bool = False): "Maximum negative stoichiometry": x_100, "Minimum positive stoichiometry": y_100, "Maximum positive stoichiometry": y_0, - "Measured cell capacity [A.s]": Q_soc / CE * 3600, + "Measured cell capacity [A.h]": Q_soc / CE, "Coulombic efficiency": CE, "Negative electrode pOCP [V]": negative_ocp_function, "Positive electrode pOCP [V]": positive_ocp_function, diff --git a/pybop/__init__.py b/pybop/__init__.py index 228d1e7f2..8d78e365f 100644 --- a/pybop/__init__.py +++ b/pybop/__init__.py @@ -56,6 +56,7 @@ # Model classes # from .models import lithium_ion +from .models import li_half_cell from .models._exponential_decay import ExponentialDecayModel from .models.lithium_ion.utils import Interpolant, InverseOCV diff --git a/pybop/applications/gitt_methods.py b/pybop/applications/gitt_methods.py index 02151c4d1..ca6fa3499 100644 --- a/pybop/applications/gitt_methods.py +++ b/pybop/applications/gitt_methods.py @@ -36,7 +36,9 @@ def __init__( ): self.parameter_values = parameter_values.copy() self.parameters = { - "Particle diffusion time scale [s]": pybop.Parameter(bounds=[0, np.inf]), + "Positive particle diffusion time scale [s]": pybop.Parameter( + bounds=[0, np.inf] + ), "Series resistance [Ohm]": pybop.Parameter(bounds=[0, np.inf]), } self.cost = cost or pybop.RootMeanSquaredError @@ -44,7 +46,7 @@ def __init__( self.optimiser_options = optimiser_options or self.optimiser.default_options() # Create model - self.model = pybop.lithium_ion.SPDiffusion() + self.model = pybop.li_half_cell.SPDiffusion() self.problem = None def __call__( @@ -110,7 +112,9 @@ def __init__( self.optimiser_options = optimiser_options or self.optimiser.default_options() # Set up OCV root-finding function - self.inverse_ocp = pybop.InverseOCV(parameter_values["Electrode OCP [V]"]) + self.inverse_ocp = pybop.InverseOCV( + parameter_values["Positive electrode OCP [V]"] + ) # Initialise single pulse fitter self.pulse_fit = GITTPulseFit( @@ -148,7 +152,9 @@ def __call__(self) -> pybop.Dataset: # Log the result self.pulses.append(pulse_result) diffusion_time.append( - pulse_result.best_inputs["Particle diffusion time scale [s]"] + pulse_result.best_inputs[ + "Positive particle diffusion time scale [s]" + ] ) series_resistance.append( pulse_result.best_inputs["Series resistance [Ohm]"] @@ -169,14 +175,16 @@ def __call__(self) -> pybop.Dataset: self.parameter_data = pybop.Dataset( { "Stoichiometry": np.asarray(stoichiometry), - "Particle diffusion time scale [s]": np.asarray(diffusion_time), + "Positive particle diffusion time scale [s]": np.asarray( + diffusion_time + ), "Series resistance [Ohm]": np.asarray(series_resistance), cost_name: np.asarray(best_cost), } if len(stoichiometry) > 1 and stoichiometry[-1] > stoichiometry[0] else { "Stoichiometry": np.flipud(np.asarray(stoichiometry)), - "Particle diffusion time scale [s]": np.flipud( + "Positive particle diffusion time scale [s]": np.flipud( np.asarray(diffusion_time) ), "Series resistance [Ohm]": np.flipud(np.asarray(series_resistance)), @@ -187,8 +195,8 @@ def __call__(self) -> pybop.Dataset: # Compute mean values self.best_inputs = { - "Particle diffusion time scale [s]": np.mean( - self.parameter_data["Particle diffusion time scale [s]"] + "Positive particle diffusion time scale [s]": np.mean( + self.parameter_data["Positive particle diffusion time scale [s]"] ), "Series resistance [Ohm]": np.mean( self.parameter_data["Series resistance [Ohm]"] diff --git a/pybop/models/alternative_functions.py b/pybop/models/alternative_functions.py new file mode 100644 index 000000000..991e93d82 --- /dev/null +++ b/pybop/models/alternative_functions.py @@ -0,0 +1,39 @@ +import pybamm + +""" Alternative functions written as classes to allow pickling. """ + + +class FunctionalDiffusionTime: + def __init__(self, r2_scale, D, c_scale): + self.r2_scale = r2_scale + self.D = D + self.c_scale = c_scale + + def __call__(self, sto, T): + r2_scale, D, c_scale = self.r2_scale, self.D, self.c_scale + return r2_scale / D(sto * c_scale, T) + + +class AsymmetricButlerVolmer: + def __init__(self, alpha): + self.alpha = alpha # cathodic transfer coefficient + + def __call__(self, sto_surf, sto_e, eta_RT_F): + alpha = self.alpha + j0 = sto_surf**alpha * (sto_e * (1 - sto_surf)) ** (1 - alpha) + return j0 * (pybamm.exp((1 - alpha) * eta_RT_F) - pybamm.exp(-alpha * eta_RT_F)) + + +class MultiphaseButlerVolmer: + def __init__(self, alpha, omega): + self.alpha = alpha + self.omega = omega + + def __call__(self, sto_surf, sto_e, eta_RT_F): + alpha, omega = self.alpha, self.omega + j0 = ( + sto_surf ** (alpha * omega) + * (1 - sto_surf) ** ((1 - alpha) * omega) + * sto_e ** (1 - alpha) + ) + return j0 * (pybamm.exp((1 - alpha) * eta_RT_F) - pybamm.exp(-alpha * eta_RT_F)) diff --git a/pybop/models/li_half_cell/__init__.py b/pybop/models/li_half_cell/__init__.py new file mode 100644 index 000000000..5ecba98cd --- /dev/null +++ b/pybop/models/li_half_cell/__init__.py @@ -0,0 +1,5 @@ +# +# Import lithium half cell models +# +from .sp_diffusion import SPDiffusion +from .weppner_huggins import WeppnerHuggins diff --git a/pybop/models/li_half_cell/base_model.py b/pybop/models/li_half_cell/base_model.py new file mode 100644 index 000000000..d75de6089 --- /dev/null +++ b/pybop/models/li_half_cell/base_model.py @@ -0,0 +1,112 @@ +import pybamm +from pybamm import FunctionParameter, Parameter +from pybamm import lithium_ion as pybamm_lithium_ion + +from pybop.models.lithium_ion.utils import InverseOCV + + +class BaseHalfCellModel(pybamm_lithium_ion.BaseModel): + """ + A base model for PyBOP's grouped-parameter lithium-ion half-cell models. + + Parameters + ---------- + name : str, optional + The name of the model. + **model_kwargs : optional + Valid PyBaMM model option keys and their values, for example: + options : dict, optional + A dictionary of options to customise the behaviour of the PyBaMM model. + build : bool, optional + If True, the model is built upon creation (default: False). + """ + + def __init__(self, name="Base Model", **model_kwargs): + super().__init__(name=name, **model_kwargs) + + def build_model(self): + """ + Build model variables and equations + Credit: PyBaMM + """ + self._build_model() + + self._built = True + pybamm.logger.info(f"Finish building {self.name}") + + @staticmethod + def set_initial_state( + initial_value, + parameter_values, + direction=None, + param=None, + inplace=True, + options=None, + inputs=None, + tol=1e-6, + ): + """ + Set the value of the initial state of charge. + + Parameters + ---------- + initial_value : float + Target initial value. + If float, interpreted as stoichiometry, must be between 0 and 1. + If string e.g. "4 V", interpreted as voltage, must be between V_min and V_max. + parameter_values : :class:`pybamm.ParameterValues` + Parameters and their corresponding values. + param : :class:`pybamm.LithiumIonParameters`, optional + The symbolic parameter set to use for the simulation. + If not provided, the default parameter set will be used. + inplace: bool, optional + If True, replace the parameters values in place. Otherwise, return a new set of + parameter values. Default is True. + options : dict-like, optional + A dictionary of options to be passed to the model, see + :class:`pybamm.BatteryModelOptions`. + inputs : dict, optional + A dictionary of input parameters to pass to the model when solving. + tol : float, optional + The tolerance for the solver used to compute the initial stoichiometries. + A lower value results in higher precision but may increase computation time. + Default is 1e-6. + """ + parameter_values = parameter_values if inplace else parameter_values.copy() + + if isinstance(initial_value, str) and initial_value.endswith("V"): + V_init = float(initial_value[:-1]) + V_min = parameter_values.evaluate( + Parameter("Lower voltage cut-off [V]"), inputs=inputs + ) + V_max = parameter_values.evaluate( + Parameter("Upper voltage cut-off [V]"), inputs=inputs + ) + + if not V_min - tol <= V_init <= V_max + tol: + raise ValueError( + f"Initial voltage {V_init}V is outside the voltage limits ({V_min}, {V_max})." + ) + + def ocv_function(sto_p): + U_p = FunctionParameter( + "Positive electrode OCP [V]", + {"Positive particle stoichiometry": sto_p}, + ) + return parameter_values.evaluate(U_p, inputs=inputs).squeeze() + + inverse_ocv = InverseOCV(ocv_function) + sto_p = inverse_ocv(V_init) + + elif isinstance(initial_value, int | float): + sto_p = initial_value + + else: + raise ValueError("Initial value must be a float or a string ending in 'V'.") + + if not 0 <= sto_p <= 1: + raise ValueError("Initial stoichiometry should be between 0 and 1.") + + parameter_values["Initial stoichiometry"] = sto_p + + return parameter_values diff --git a/pybop/models/li_half_cell/sp_diffusion.py b/pybop/models/li_half_cell/sp_diffusion.py new file mode 100644 index 000000000..366144e4c --- /dev/null +++ b/pybop/models/li_half_cell/sp_diffusion.py @@ -0,0 +1,267 @@ +import pybamm +from pybamm import ( + Event, + FunctionParameter, + Parameter, + ParameterValues, + PrimaryBroadcast, + Scalar, + SpatialVariable, + Variable, +) +from pybamm import t as pybamm_t + +from pybop.models.li_half_cell.base_model import BaseHalfCellModel + + +class SPDiffusion(BaseHalfCellModel): + """ + Diffusion model for a single, spherical particle representing a half-cell for GITT. + + Note: the working electrode is the positive electrode. + + Parameters + ---------- + name : str, optional + The name of the model. + **model_kwargs : optional + Valid PyBaMM model option keys and their values, for example: + options : dict, optional + A dictionary of options to customise the behaviour of the PyBaMM model. + build : bool, optional + If True, the model is built upon creation (default: False). + """ + + def __init__(self, name="Single Particle Diffusion Model", **model_kwargs): + super().__init__(name=name, **model_kwargs) + + ###################### + # Variables + ###################### + # Variables that depend on time only are created without a domain + Q = Variable("Discharge capacity [A.h]") + Qt = Variable("Throughput capacity [A.h]") + + # Variables that vary spatially are created with a domain + sto_p = Variable("Positive particle stoichiometry", domain="positive particle") + + # Surf takes the surface value of a variable, i.e. its boundary value on the + # right side. This is also accessible via `boundary_value(x, "right")`, with + # "left" providing the boundary value of the left side + sto_p_surf = pybamm.surf(sto_p) + + # Events specify points at which a solution should terminate + tol = pybamm.settings.tolerances["U__c_s"] + self.events += [ + Event( + "Minimum positive particle surface stoichiometry", + pybamm.min(sto_p_surf) - tol, + ), + Event( + "Maximum positive particle surface stoichiometry", + (1 - tol) - pybamm.max(sto_p_surf), + ), + ] + + ###################### + # Parameters + ###################### + # Parameters are purely symbolic at this stage, and will be set by the + # `ParameterValues` class when the model is processed. + + # Grouped parameters + Q_th_p = Parameter("Theoretical electrode capacity [A.h]") * 3600 + + sto_p_init = Parameter("Initial stoichiometry") + + ###################### + # Input current (positive on discharge) + ###################### + I = self.param.current_with_time + + ###################### + # State of Charge + ###################### + # The `rhs` dictionary contains differential equations, with the key being the + # variable in the d/dt + self.rhs[Q] = I / 3600 + self.rhs[Qt] = abs(I) / 3600 + # Initial conditions must be provided for the ODEs + self.initial_conditions[Q] = Scalar(0) + self.initial_conditions[Qt] = Scalar(0) + + ###################### + # Diffusion within the particle + ###################### + # The div and grad operators will be converted to the appropriate matrix + # multiplication at the discretisation stage + self.rhs[sto_p] = pybamm.div(pybamm.grad(sto_p) / self.tau_d(sto_p)) + + # Boundary conditions must be provided for equations with spatial derivatives + j_p = -I / (3 * Q_th_p) + self.boundary_conditions[sto_p] = { + "left": (Scalar(0), "Neumann"), + "right": (-self.tau_d(sto_p_surf) * j_p, "Neumann"), + } + + self.initial_conditions[sto_p] = sto_p_init + + ###################### + # Cell voltage + ###################### + U = self.U(sto_p_surf) + V = U - self.R0(sto_p_surf) * I + + # Save the initial OCV + self.param.ocv_init = self.U(sto_p_init) + + # Events specify points at which a solution should terminate + self.events += [ + Event("Minimum voltage [V]", V - self.param.voltage_low_cut), + Event("Maximum voltage [V]", self.param.voltage_high_cut - V), + ] + + ###################### + # (Some) variables + ###################### + self.variables = { + "Positive particle stoichiometry": sto_p, + "Positive particle surface stoichiometry": PrimaryBroadcast( + sto_p_surf, "positive particle" + ), + "Time [s]": pybamm_t, + "Current [A]": I, + "Current variable [A]": I, # for compatibility with pybamm.Experiment + "Discharge capacity [A.h]": Q, + "Throughput capacity [A.h]": Qt, + "Voltage [V]": V, + "Voltage expression [V]": V, # for compatibility with "voltage as a state" + "Open-circuit voltage [V]": U, + } + + def U(self, sto): + """ + Dimensional open-circuit potential [V]. + Credit: PyBaMM + """ + inputs = {"Positive particle surface stoichiometry": sto} + out = FunctionParameter("Positive electrode OCP [V]", inputs) + + out.print_name = r"U_\mathrm{p}(c^\mathrm{surf}_\mathrm{s,p})" + return out + + def tau_d(self, sto): + """ + Dimensional solid-state diffusion time scale [s]. + """ + inputs = {"Positive particle surface stoichiometry": sto} + return FunctionParameter("Positive particle diffusion time scale [s]", inputs) + + def R0(self, sto): + """ + Series resistance [Ohm]. + """ + inputs = {"Positive particle surface stoichiometry": sto} + return FunctionParameter("Series resistance [Ohm]", inputs) + + @property + def default_parameter_values(self) -> ParameterValues: + param = ParameterValues("Xu2019") + return self.create_grouped_parameters(param) + + @property + def default_quick_plot_variables(self): + return [ + "Positive particle stoichiometry", + "Positive particle surface stoichiometry", + "Current [A]", + {"Open-circuit voltage [V]", "Voltage [V]"}, + ] + + @property + def default_var_pts(self): + r_p = SpatialVariable( + "r_p", domain=["positive particle"], coord_sys="spherical polar" + ) + return {r_p: 20} + + @property + def default_geometry(self): + return {"positive particle": {"r_p": {"min": 0, "max": 1}}} + + @property + def default_submesh_types(self): + return {"positive particle": pybamm.Uniform1DSubMesh} + + @property + def default_spatial_methods(self): + return {"positive particle": pybamm.FiniteVolume()} + + @staticmethod + def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterValues: + """ + Create a parameter set for the Single Particle Diffusion Model from a + PyBaMM lithium-ion ParameterValues object. + + Note: the working electrode is the positive electrode. + + Parameters + ---------- + parameter_values : pybamm.ParameterValues + Parameters and their corresponding values. + + Returns + ------- + parameter_values : pybamm.ParameterValues + A new set of parameters and their values. + """ + param = parameter_values + + # Unpack physical parameters + F = pybamm.constants.F.value + T = param["Ambient temperature [K]"] + alpha_p = param["Positive electrode active material volume fraction"] + c_max_p = param["Maximum concentration in positive electrode [mol.m-3]"] + L_p = param["Positive electrode thickness [m]"] + R_p = param["Positive particle radius [m]"] + D_p = param["Positive particle diffusivity [m2.s-1]"] + + # Get reference exchange current density [A.m-2] + ce0 = param["Initial concentration in electrolyte [mol.m-3]"] + j0_p = param.evaluate( + param["Positive electrode exchange-current density [A.m-2]"]( + ce0, c_max_p / 2, c_max_p, T + ) + ) + + # Compute the cell area + A = param["Electrode height [m]"] * param["Electrode width [m]"] + + # Compute the initial stoichiometry + sto_p_init = ( + param["Initial concentration in positive electrode [mol.m-3]"] / c_max_p + ) + + # Grouped parameters + Q_th_p = F * alpha_p * c_max_p * L_p * A / 3600 + tau_d_p = R_p**2 / D_p + + # Estimate the series resistance, neglecting conductivity losses + RT_F = pybamm.constants.R.value * T / F + Rct_typ = (RT_F * R_p) / (3 * alpha_p * L_p * A * j0_p) + R0 = Rct_typ + param["Contact resistance [Ohm]"] + + parameter_dictionary = { + "Nominal cell capacity [A.h]": param["Nominal cell capacity [A.h]"], + "Current function [A]": param["Current function [A]"], + "Initial stoichiometry": sto_p_init, + "Lower voltage cut-off [V]": param["Lower voltage cut-off [V]"], + "Upper voltage cut-off [V]": param["Upper voltage cut-off [V]"], + "Positive electrode OCP [V]": param["Positive electrode OCP [V]"], + "Theoretical electrode capacity [A.h]": Q_th_p, + "Positive particle diffusion time scale [s]": tau_d_p, + "Series resistance [Ohm]": R0, + } + parameter_values = ParameterValues(values=parameter_dictionary) + parameter_values._set_initial_state = SPDiffusion.set_initial_state # noqa: SLF001 + return parameter_values diff --git a/pybop/models/lithium_ion/weppner_huggins.py b/pybop/models/li_half_cell/weppner_huggins.py similarity index 83% rename from pybop/models/lithium_ion/weppner_huggins.py rename to pybop/models/li_half_cell/weppner_huggins.py index 37b752a42..78f007f2a 100644 --- a/pybop/models/lithium_ion/weppner_huggins.py +++ b/pybop/models/li_half_cell/weppner_huggins.py @@ -3,10 +3,10 @@ from pybamm import DummySolver, Parameter, ParameterValues from pybamm import t as pybamm_t -from pybop.models.lithium_ion.base_model import BaseGroupedModel +from pybop.models.li_half_cell.base_model import BaseHalfCellModel -class WeppnerHuggins(BaseGroupedModel): +class WeppnerHuggins(BaseHalfCellModel): """ Represents the Weppner & Huggins model to fit diffusion coefficients to GITT data. @@ -48,12 +48,12 @@ def __init__(self, name="Weppner & Huggins model", **model_kwargs): # Parameters are purely symbolic at this stage, and will be set by the # `ParameterValues` class when the model is processed. - Q_th = Parameter("Theoretical electrode capacity [A.s]") + Q_th_p = Parameter("Theoretical electrode capacity [A.h]") * 3600 U = Parameter("Reference voltage [V]") U_prime = Parameter("Derivative of the OCP wrt stoichiometry [V]") - tau_d = Parameter("Particle diffusion time scale [s]") + tau_d_p = Parameter("Positive particle diffusion time scale [s]") ###################### # Input current (positive on discharge) @@ -64,9 +64,9 @@ def __init__(self, name="Weppner & Huggins model", **model_kwargs): # Governing equations ###################### # Surface stoichiometry - sto_surf = 2 * I / (3 * Q_th) * (pybamm_t * tau_d / np.pi) ** 0.5 + sto_p_surf = 2 * I / (3 * Q_th_p) * (pybamm_t * tau_d_p / np.pi) ** 0.5 # Linearised voltage - V = U + U_prime * sto_surf + V = U + U_prime * sto_p_surf ###################### # (Some) variables @@ -136,25 +136,25 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal # Unpack physical parameters F = pybamm.constants.F.value - alpha = param["Positive electrode active material volume fraction"] - c_max = param["Maximum concentration in positive electrode [mol.m-3]"] - L = param["Positive electrode thickness [m]"] - R = param["Positive particle radius [m]"] - D = param["Positive particle diffusivity [m2.s-1]"] + alpha_p = param["Positive electrode active material volume fraction"] + c_max_p = param["Maximum concentration in positive electrode [mol.m-3]"] + L_p = param["Positive electrode thickness [m]"] + R_p = param["Positive particle radius [m]"] + D_p = param["Positive particle diffusivity [m2.s-1]"] # Compute the cell area A = param["Electrode height [m]"] * param["Electrode width [m]"] # Grouped parameters - Q_th = F * alpha * c_max * L * A - tau_d = R**2 / D + Q_th_p = F * alpha_p * c_max_p * L_p * A / 3600 + tau_d_p = R_p**2 / D_p parameter_dictionary = { "Current function [A]": param["Current function [A]"], "Reference voltage [V]": 4, "Derivative of the OCP wrt stoichiometry [V]": -1, - "Theoretical electrode capacity [A.s]": Q_th, - "Particle diffusion time scale [s]": tau_d, + "Theoretical electrode capacity [A.h]": Q_th_p, + "Positive particle diffusion time scale [s]": tau_d_p, } parameter_values = ParameterValues(values=parameter_dictionary) parameter_values._set_initial_state = WeppnerHuggins.set_initial_state # noqa: SLF001 diff --git a/pybop/models/lithium_ion/__init__.py b/pybop/models/lithium_ion/__init__.py index 7dcbbc032..1745f9b42 100644 --- a/pybop/models/lithium_ion/__init__.py +++ b/pybop/models/lithium_ion/__init__.py @@ -1,8 +1,6 @@ # # Import lithium ion models # -from .sp_diffusion import SPDiffusion from .grouped_spm import GroupedSPM from .grouped_spme import GroupedSPMe -from .weppner_huggins import WeppnerHuggins from .cell_temperature import CellTemperature diff --git a/pybop/models/lithium_ion/base_model.py b/pybop/models/lithium_ion/base_model.py index 2f4733060..9844b29aa 100644 --- a/pybop/models/lithium_ion/base_model.py +++ b/pybop/models/lithium_ion/base_model.py @@ -2,6 +2,10 @@ from pybamm import FunctionParameter, Parameter from pybamm import lithium_ion as pybamm_lithium_ion +from pybop.models.alternative_functions import ( + AsymmetricButlerVolmer, + MultiphaseButlerVolmer, +) from pybop.models.lithium_ion.utils import InverseOCV @@ -52,7 +56,7 @@ def set_initial_state( ---------- initial_value : float Target initial value. - If float, interpreted as SOC, must be between 0 and 1. + If float, interpreted as SoC, must be between 0 and 1. If string e.g. "4 V", interpreted as voltage, must be between V_min and V_max. parameter_values : :class:`pybamm.ParameterValues` Parameters and their corresponding values. @@ -77,10 +81,10 @@ def set_initial_state( if isinstance(initial_value, str) and initial_value.endswith("V"): V_init = float(initial_value[:-1]) V_min = parameter_values.evaluate( - pybamm.Parameter("Lower voltage cut-off [V]"), inputs=inputs + Parameter("Lower voltage cut-off [V]"), inputs=inputs ) V_max = parameter_values.evaluate( - pybamm.Parameter("Upper voltage cut-off [V]"), inputs=inputs + Parameter("Upper voltage cut-off [V]"), inputs=inputs ) if not V_min - tol <= V_init <= V_max + tol: @@ -104,8 +108,7 @@ def ocv_function(soc): "Negative electrode OCP [V]", {"Negative particle stoichiometry": sto_n}, ) - - return parameter_values.evaluate(U_p - U_n, inputs=inputs) + return parameter_values.evaluate(U_p - U_n, inputs=inputs).squeeze() inverse_ocv = InverseOCV(ocv_function) soc = inverse_ocv(V_init) @@ -117,8 +120,39 @@ def ocv_function(soc): raise ValueError("Initial value must be a float or a string ending in 'V'.") if not 0 <= soc <= 1: - raise ValueError("Initial SOC should be between 0 and 1.") + raise ValueError("Initial SoC should be between 0 and 1.") parameter_values["Initial SoC"] = soc return parameter_values + + @staticmethod + def symmetric_butler_volmer(sto_surf, sto_e, eta_RT_F): + """ + Dimensionless Butler-Volmer exchange rate. + """ + j0 = (sto_surf * sto_e * (1 - sto_surf)) ** 0.5 + return 2 * j0 * pybamm.sinh(0.5 * eta_RT_F) + + @staticmethod + def get_asymmetric_butler_volmer(domain: str): + """ + Get the asymmetric Butler-Volmer exchange rate function for the "positive" or + "negative" domain. + """ + Domain = domain.capitalize() + alpha = pybamm.Parameter(f"{Domain} electrode charge transfer coefficient") + + return AsymmetricButlerVolmer(alpha) + + @staticmethod + def get_multiphase_butler_volmer(domain: str): + """ + Get the multiphase Butler-Volmer exchange rate function for the "positive" or + "negative" domain. + """ + Domain = domain.capitalize() + alpha = pybamm.Parameter(f"{Domain} electrode charge transfer coefficient") + omega = pybamm.Parameter(f"{Domain} electrode charge transfer ideality factor") + + return MultiphaseButlerVolmer(alpha, omega) diff --git a/pybop/models/lithium_ion/cell_temperature.py b/pybop/models/lithium_ion/cell_temperature.py index 3fbb32f31..3b5c40132 100644 --- a/pybop/models/lithium_ion/cell_temperature.py +++ b/pybop/models/lithium_ion/cell_temperature.py @@ -52,7 +52,7 @@ def __init__(self, name="Cell Temperature Model", **model_kwargs): # Parameters are purely symbolic at this stage, and will be set by the # `ParameterValues` class when the model is processed. - Q_meas = Parameter("Measured cell capacity [A.s]") + Q_meas = Parameter("Measured cell capacity [A.h]") * 3600 soc_init = Parameter("Initial SoC") x_0 = Parameter("Minimum negative stoichiometry") @@ -145,21 +145,12 @@ def __init__(self, name="Cell Temperature Model", **model_kwargs): def U(self, sto, domain): """ - Dimensional open-circuit potential [V], calculated as U(x) = U_ref(x). + Dimensional open-circuit potential [V]. Credit: PyBaMM """ - # bound stoichiometry between tol and 1-tol. Adding 1/sto + 1/(sto-1) later - # will ensure that ocp goes to +- infinity if sto goes into that region - # anyway Domain = domain.capitalize() - tol = pybamm.settings.tolerances["U__c_s"] - sto = pybamm.maximum(pybamm.minimum(sto, 1 - tol), tol) inputs = {f"{Domain} particle surface stoichiometry": sto} - u_ref = FunctionParameter(f"{Domain} electrode OCP [V]", inputs) - - # add a term to ensure that the OCP goes to infinity at 0 and -infinity at 1 - # this will not affect the OCP for most values of sto - out = u_ref + 1e-6 * (1 / sto + 1 / (sto - 1)) + out = FunctionParameter(f"{Domain} electrode OCP [V]", inputs) if domain == "negative": out.print_name = r"U_\mathrm{n}(c^\mathrm{surf}_\mathrm{s,n})" @@ -251,8 +242,8 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal soc_init = (sto_p_init - y_0) / (y_100 - y_0) # Compute the capacity within the stoichiometry limits - Q_th_p = F * alpha_p * c_max_p * L_p * A - Q_th_n = F * alpha_n * c_max_n * L_n * A + Q_th_p = F * alpha_p * c_max_p * L_p * A / 3600 + Q_th_n = F * alpha_n * c_max_n * L_n * A / 3600 Q_meas_p = (y_0 - y_100) * Q_th_p Q_meas_n = (x_100 - x_0) * Q_th_n if abs(Q_meas_n / Q_meas_p - 1) > 1e-6: @@ -279,7 +270,7 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal "Upper voltage cut-off [V]": param["Upper voltage cut-off [V]"], "Positive electrode OCP [V]": param["Positive electrode OCP [V]"], "Negative electrode OCP [V]": param["Negative electrode OCP [V]"], - "Measured cell capacity [A.s]": Q_meas, + "Measured cell capacity [A.h]": Q_meas, "Negative electrode OCP entropic change [V.K-1]": param[ "Negative electrode OCP entropic change [V.K-1]" ], diff --git a/pybop/models/lithium_ion/grouped_spm.py b/pybop/models/lithium_ion/grouped_spm.py index 0f12ebfa8..251854db0 100644 --- a/pybop/models/lithium_ion/grouped_spm.py +++ b/pybop/models/lithium_ion/grouped_spm.py @@ -1,4 +1,3 @@ -import numpy as np import pybamm from pybamm import ( Event, @@ -14,6 +13,7 @@ get_min_max_stoichiometries, ) +from pybop.models.alternative_functions import FunctionalDiffusionTime from pybop.models.lithium_ion.base_model import BaseGroupedModel @@ -62,6 +62,9 @@ def __init__(self, name="Grouped Single Particle Model", **model_kwargs): Q = Variable("Discharge capacity [A.h]") Qt = Variable("Throughput capacity [A.h]") + v_s_n = Variable("Negative particle surface voltage variable [V]") + v_s_p = Variable("Positive particle surface voltage variable [V]") + # Variables that vary spatially are created with a domain sto_n = Variable( "Negative particle stoichiometry", @@ -79,22 +82,23 @@ def __init__(self, name="Grouped Single Particle Model", **model_kwargs): sto_p_surf = pybamm.surf(sto_p) # Events specify points at which a solution should terminate + tol = pybamm.settings.tolerances["U__c_s"] self.events += [ Event( "Minimum negative particle surface stoichiometry", - pybamm.min(sto_n_surf) - 0.01, + pybamm.min(sto_n_surf) - tol, ), Event( "Maximum negative particle surface stoichiometry", - (1 - 0.01) - pybamm.max(sto_n_surf), + (1 - tol) - pybamm.max(sto_n_surf), ), Event( "Minimum positive particle surface stoichiometry", - pybamm.min(sto_p_surf) - 0.01, + pybamm.min(sto_p_surf) - tol, ), Event( "Maximum positive particle surface stoichiometry", - (1 - 0.01) - pybamm.max(sto_p_surf), + (1 - tol) - pybamm.max(sto_p_surf), ), ] @@ -116,11 +120,8 @@ def __init__(self, name="Grouped Single Particle Model", **model_kwargs): y_0 = Parameter("Maximum positive stoichiometry") # Grouped parameters - Q_th_p = Parameter("Measured cell capacity [A.s]") / (y_0 - y_100) - Q_th_n = Parameter("Measured cell capacity [A.s]") / (x_100 - x_0) - - tau_d_p = Parameter("Positive particle diffusion time scale [s]") - tau_d_n = Parameter("Negative particle diffusion time scale [s]") + Q_th_p = Parameter("Measured cell capacity [A.h]") * 3600 / (y_0 - y_100) + Q_th_n = Parameter("Measured cell capacity [A.h]") * 3600 / (x_100 - x_0) tau_ct_p = Parameter("Positive electrode charge transfer time scale [s]") tau_ct_n = Parameter("Negative electrode charge transfer time scale [s]") @@ -159,66 +160,59 @@ def __init__(self, name="Grouped Single Particle Model", **model_kwargs): ###################### # Primary broadcasts are used to broadcast scalar quantities across a domain # into a vector of the right shape, for multiplying with other vectors - alpha = 0.5 # cathodic transfer coefficient - j0_n = sto_n_surf**alpha * (1 - sto_n_surf) ** (1 - alpha) / tau_ct_n - j0_p = sto_p_surf**alpha * (1 - sto_p_surf) ** (1 - alpha) / tau_ct_p - if not include_double_layer: - # Assuming alpha = 0.5 - j_n = PrimaryBroadcast(I / (3 * Q_th_n), "negative electrode") - j_p = PrimaryBroadcast(-I / (3 * Q_th_p), "positive electrode") - eta_n = 2 * RT_F * pybamm.arcsinh(j_n / (2 * j0_n)) - eta_p = 2 * RT_F * pybamm.arcsinh(j_p / (2 * j0_p)) - v_s_n = pybamm.x_average(eta_n + U_n) - v_s_p = pybamm.x_average(eta_p + U_p) + + # Overpotentials + eta_n = PrimaryBroadcast(v_s_n - U_n, "negative electrode") + eta_p = PrimaryBroadcast(v_s_p - U_p, "positive electrode") + + # Exchange rates + j_n = self.j(sto_n_surf, 1.0, eta_n / RT_F, "negative") / tau_ct_n + j_p = self.j(sto_p_surf, 1.0, eta_p / RT_F, "positive") / tau_ct_p ###################### # Double layer ###################### if include_double_layer: - # Additional variables - v_s_n = Variable("Negative particle surface voltage variable [V]") - v_s_p = Variable("Positive particle surface voltage variable [V]") - # Additional parameters C_p = Parameter("Positive electrode capacitance [F]") C_n = Parameter("Negative electrode capacitance [F]") - # Overpotentials - eta_n = PrimaryBroadcast(v_s_n - U_n, "negative electrode") - eta_p = PrimaryBroadcast(v_s_p - U_p, "positive electrode") - - # Exchange current - j_n = j0_n * ( - pybamm.exp((1 - alpha) * eta_n / RT_F) - - pybamm.exp(-alpha * eta_n / RT_F) - ) - j_p = j0_p * ( - pybamm.exp((1 - alpha) * eta_p / RT_F) - - pybamm.exp(-alpha * eta_p / RT_F) - ) - # Electrode surface potentials - self.rhs[v_s_n] = 1 / C_n * (I - 3 * Q_th_n * pybamm.x_average(j_n)) - self.rhs[v_s_p] = 1 / C_p * (-I - 3 * Q_th_p * pybamm.x_average(j_p)) - self.initial_conditions[v_s_n] = U_n_init - self.initial_conditions[v_s_p] = U_p_init + self.rhs[v_s_n] = (I - 3 * Q_th_n * pybamm.x_average(j_n)) / C_n + self.rhs[v_s_p] = (-I - 3 * Q_th_p * pybamm.x_average(j_p)) / C_p + else: + self.algebraic[v_s_n] = I - 3 * Q_th_n * pybamm.x_average(j_n) + self.algebraic[v_s_p] = -I - 3 * Q_th_p * pybamm.x_average(j_p) + + self.initial_conditions[v_s_n] = U_n_init + self.initial_conditions[v_s_p] = U_p_init ###################### # Particles ###################### # The div and grad operators will be converted to the appropriate matrix # multiplication at the discretisation stage - self.rhs[sto_n] = pybamm.div(pybamm.grad(sto_n) / tau_d_n) - self.rhs[sto_p] = pybamm.div(pybamm.grad(sto_p) / tau_d_p) + self.rhs[sto_n] = pybamm.div( + pybamm.grad(sto_n) / self.tau_d(sto_n, T, "negative") + ) + self.rhs[sto_p] = pybamm.div( + pybamm.grad(sto_p) / self.tau_d(sto_p, T, "positive") + ) # Boundary conditions must be provided for equations with spatial derivatives self.boundary_conditions[sto_n] = { "left": (Scalar(0), "Neumann"), - "right": (-tau_d_n * pybamm.x_average(j_n), "Neumann"), + "right": ( + -self.tau_d(sto_n_surf, T, "negative") * pybamm.x_average(j_n), + "Neumann", + ), } self.boundary_conditions[sto_p] = { "left": (Scalar(0), "Neumann"), - "right": (-tau_d_p * pybamm.x_average(j_p), "Neumann"), + "right": ( + -self.tau_d(sto_p_surf, T, "positive") * pybamm.x_average(j_p), + "Neumann", + ), } self.initial_conditions[sto_n] = sto_n_init @@ -238,6 +232,52 @@ def __init__(self, name="Grouped Single Particle Model", **model_kwargs): Event("Maximum voltage [V]", self.param.voltage_high_cut - V), ] + ###################### + # Voltage components + ###################### + # Include the following variables to enable plotting via PyBaMM's plot_voltage_components + ocp_n_bulk = self.U( + pybamm.x_average(Q_th_n * pybamm.r_average(sto_n)) + / pybamm.x_average(Q_th_n), + "negative", + ) + ocp_p_bulk = self.U( + pybamm.x_average(Q_th_p * pybamm.r_average(sto_p)) + / pybamm.x_average(Q_th_p), + "positive", + ) + voltage_components = { + "Battery voltage [V]": V, + "Battery open-circuit voltage [V]": ocp_p_bulk - ocp_n_bulk, + "Battery particle concentration overpotential [V]": ( + (pybamm.x_average(self.U(sto_p_surf, "positive")) - ocp_p_bulk) + - (pybamm.x_average(self.U(sto_n_surf, "negative")) - ocp_n_bulk) + ), + "X-averaged battery reaction overpotential [V]": pybamm.x_average(eta_p) + - pybamm.x_average(eta_n), + "X-averaged battery concentration overpotential [V]": Scalar(0), + "X-averaged battery electrolyte ohmic losses [V]": Scalar(0), + "X-averaged battery solid phase ohmic losses [V]": Scalar(0), + "Contact overpotential [V]": R0 * I, # includes Ohmic losses in this model + # and split by electrode + "Negative electrode bulk open-circuit potential [V]": ocp_n_bulk, + "Positive electrode bulk open-circuit potential [V]": ocp_p_bulk, + "Negative particle concentration overpotential [V]": pybamm.x_average( + self.U(sto_n_surf, "negative") + ) + - ocp_n_bulk, + "Positive particle concentration overpotential [V]": pybamm.x_average( + self.U(sto_p_surf, "positive") + ) + - ocp_p_bulk, + "X-averaged negative electrode reaction overpotential [V]" + "": pybamm.x_average(eta_n), + "X-averaged positive electrode reaction overpotential [V]" + "": pybamm.x_average(eta_p), + "X-averaged battery negative solid phase ohmic losses [V]": Scalar(0), + "X-averaged battery positive solid phase ohmic losses [V]": Scalar(0), + } + ###################### # (Some) variables ###################### @@ -266,33 +306,25 @@ def __init__(self, name="Grouped Single Particle Model", **model_kwargs): + eta_p - pybamm.boundary_value(eta_p, "right"), "Time [s]": pybamm_t, + "Time [h]": pybamm_t / 3600, "Current [A]": I, "Current variable [A]": I, # for compatibility with pybamm.Experiment "Discharge capacity [A.h]": Q, "Throughput capacity [A.h]": Qt, "Voltage [V]": V, "Voltage expression [V]": V, # for compatibility with "voltage as a state" - "Battery voltage [V]": V, "Open-circuit voltage [V]": U_p - U_n, + **voltage_components, } def U(self, sto, domain): """ - Dimensional open-circuit potential [V], calculated as U(x) = U_ref(x). + Dimensional open-circuit potential [V]. Credit: PyBaMM """ - # bound stoichiometry between tol and 1-tol. Adding 1/sto + 1/(sto-1) later - # will ensure that ocp goes to +- infinity if sto goes into that region - # anyway Domain = domain.capitalize() - tol = pybamm.settings.tolerances["U__c_s"] - sto = pybamm.maximum(pybamm.minimum(sto, 1 - tol), tol) inputs = {f"{Domain} particle surface stoichiometry": sto} - u_ref = FunctionParameter(f"{Domain} electrode OCP [V]", inputs) - - # add a term to ensure that the OCP goes to infinity at 0 and -infinity at 1 - # this will not affect the OCP for most values of sto - out = u_ref + 1e-6 * (1 / sto + 1 / (sto - 1)) + out = FunctionParameter(f"{Domain} electrode OCP [V]", inputs) if domain == "negative": out.print_name = r"U_\mathrm{n}(c^\mathrm{surf}_\mathrm{s,n})" @@ -300,6 +332,28 @@ def U(self, sto, domain): out.print_name = r"U_\mathrm{p}(c^\mathrm{surf}_\mathrm{s,p})" return out + def tau_d(self, sto, T, domain): + """ + Dimensional solid-state diffusion time scale [s]. + """ + Domain = domain.capitalize() + inputs = {f"{Domain} particle surface stoichiometry": sto, "Temperature [K]": T} + return FunctionParameter(f"{Domain} particle diffusion time scale [s]", inputs) + + def j(self, sto_surf, sto_e, eta_RT_F, domain): + """ + Dimensionless exchange rate. + """ + Domain = domain.capitalize() + inputs = { + f"{Domain} particle surface stoichiometry": sto_surf, + f"{Domain} electrode electrolyte stoichiometry": sto_e, + f"{Domain} electrode dimensionless overpotential": eta_RT_F, + } + return FunctionParameter( + f"{Domain} electrode dimensionless exchange rate", inputs + ) + @property def default_parameter_values(self) -> ParameterValues: param = ParameterValues("Chen2020") @@ -313,8 +367,8 @@ def default_parameter_values(self) -> ParameterValues: @property def default_quick_plot_variables(self): return [ - "Negative particle surface stoichiometry", - "Positive particle surface stoichiometry", + "Negative particle stoichiometry", + "Positive particle stoichiometry", "Current [A]", { "Negative electrode potential [V]", @@ -414,6 +468,7 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal # Unpack physical parameters F = pybamm.constants.F.value + T = param["Ambient temperature [K]"] alpha_p = param["Positive electrode active material volume fraction"] alpha_n = param["Negative electrode active material volume fraction"] c_max_p = param["Maximum concentration in positive electrode [mol.m-3]"] @@ -430,8 +485,6 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal b_n = param["Negative electrode Bruggeman coefficient (electrolyte)"] Cdl_p = param["Positive electrode double-layer capacity [F.m-2]"] Cdl_n = param["Negative electrode double-layer capacity [F.m-2]"] - m_p = 3.42e-6 # (A/m2)(m3/mol)**1.5 - m_n = 6.48e-7 # (A/m2)(m3/mol)**1.5 sigma_p = ( param["Positive electrode conductivity [S.m-1]"] * alpha_p ** param["Positive electrode Bruggeman coefficient (electrode)"] @@ -446,7 +499,19 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal L_s = param["Separator thickness [m]"] epsilon_sep = param["Separator porosity"] b_sep = param["Separator Bruggeman coefficient (electrolyte)"] - sigma_e = param["Electrolyte conductivity [S.m-1]"] # (ce0, T) + kappa_e = param["Electrolyte conductivity [S.m-1]"] # (ce0, T) + + # Get reference exchange current density [A.m-2] + j0_p = param.evaluate( + param["Positive electrode exchange-current density [A.m-2]"]( + ce0, c_max_p / 2, c_max_p, T + ) + ) + j0_n = param.evaluate( + param["Negative electrode exchange-current density [A.m-2]"]( + ce0, c_max_n / 2, c_max_n, T + ) + ) # Compute the cell area and thickness A = param["Electrode height [m]"] * param["Electrode width [m]"] @@ -457,7 +522,7 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal L_p / (3 * epsilon_p**b_p) + L_s / (epsilon_sep**b_sep) + L_n / (3 * epsilon_n**b_n) - ) / (sigma_e * A) + ) / (kappa_e * A) Rs = (L_p / sigma_p + L_n / sigma_n) / (3 * A) R0 = Re + Rs + param["Contact resistance [Ohm]"] @@ -469,8 +534,8 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal soc_init = (sto_p_init - y_0) / (y_100 - y_0) # Compute the capacity within the stoichiometry limits - Q_th_p = F * alpha_p * c_max_p * L_p * A - Q_th_n = F * alpha_n * c_max_n * L_n * A + Q_th_p = F * alpha_p * c_max_p * L_p * A / 3600 + Q_th_n = F * alpha_n * c_max_n * L_n * A / 3600 Q_meas_p = (y_0 - y_100) * Q_th_p Q_meas_n = (x_100 - x_0) * Q_th_n if abs(Q_meas_n / Q_meas_p - 1) > 1e-6: @@ -481,11 +546,18 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal # Grouped parameters Q_meas = (Q_meas_n + Q_meas_p) / 2 - tau_d_p = R_p**2 / D_p - tau_d_n = R_n**2 / D_n + try: + tau_d_p = R_p**2 / D_p + except TypeError: + tau_d_p = FunctionalDiffusionTime(R_p**2, D_p, c_max_p) + + try: + tau_d_n = R_n**2 / D_n + except TypeError: + tau_d_n = FunctionalDiffusionTime(R_n**2, D_n, c_max_n) - tau_ct_p = F * R_p / (m_p * np.sqrt(ce0)) - tau_ct_n = F * R_n / (m_n * np.sqrt(ce0)) + tau_ct_p = c_max_p * F * R_p / (2 * j0_p) + tau_ct_n = c_max_n * F * R_n / (2 * j0_n) C_p = 3 * alpha_p * Cdl_p * L_p * A / R_p C_n = 3 * alpha_n * Cdl_n * L_n * A / R_n @@ -496,8 +568,8 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal parameter_dictionary = { "Nominal cell capacity [A.h]": param["Nominal cell capacity [A.h]"], "Current function [A]": param["Current function [A]"], - "Ambient temperature [K]": param["Ambient temperature [K]"], - "Initial temperature [K]": param["Ambient temperature [K]"], + "Ambient temperature [K]": T, + "Initial temperature [K]": T, "Initial SoC": soc_init, "Minimum negative stoichiometry": x_0, "Maximum negative stoichiometry": x_100, @@ -507,9 +579,11 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal "Upper voltage cut-off [V]": param["Upper voltage cut-off [V]"], "Positive electrode OCP [V]": param["Positive electrode OCP [V]"], "Negative electrode OCP [V]": param["Negative electrode OCP [V]"], - "Measured cell capacity [A.s]": Q_meas, + "Measured cell capacity [A.h]": Q_meas, "Positive particle diffusion time scale [s]": tau_d_p, "Negative particle diffusion time scale [s]": tau_d_n, + "Positive electrode dimensionless exchange rate": GroupedSPM.symmetric_butler_volmer, + "Negative electrode dimensionless exchange rate": GroupedSPM.symmetric_butler_volmer, "Positive electrode charge transfer time scale [s]": tau_ct_p, "Negative electrode charge transfer time scale [s]": tau_ct_n, "Positive electrode capacitance [F]": C_p, diff --git a/pybop/models/lithium_ion/grouped_spme.py b/pybop/models/lithium_ion/grouped_spme.py index d3713eab0..e9c027110 100644 --- a/pybop/models/lithium_ion/grouped_spme.py +++ b/pybop/models/lithium_ion/grouped_spme.py @@ -1,4 +1,3 @@ -import numpy as np import pybamm from pybamm import ( Event, @@ -15,6 +14,7 @@ get_min_max_stoichiometries, ) +from pybop.models.alternative_functions import FunctionalDiffusionTime from pybop.models.lithium_ion.base_model import BaseGroupedModel @@ -57,7 +57,7 @@ def __init__( doi = {10.1149/1945-7111/add41b} } """ - ) + ) # Note that the electrode electrolyte timescales have been replaced by relative transport efficiencies ###################### # Variables @@ -66,6 +66,9 @@ def __init__( Q = Variable("Discharge capacity [A.h]") Qt = Variable("Throughput capacity [A.h]") + v_s_n = Variable("Negative particle surface voltage variable [V]") + v_s_p = Variable("Positive particle surface voltage variable [V]") + # Variables that vary spatially are created with a domain sto_n = Variable( "Negative particle stoichiometry", @@ -99,22 +102,36 @@ def __init__( sto_p_surf = pybamm.surf(sto_p) # Events specify points at which a solution should terminate + tol = pybamm.settings.tolerances["U__c_s"] self.events += [ Event( "Minimum negative particle surface stoichiometry", - pybamm.min(sto_n_surf) - 0.01, + pybamm.min(sto_n_surf) - tol, ), Event( "Maximum negative particle surface stoichiometry", - (1 - 0.01) - pybamm.max(sto_n_surf), + (1 - tol) - pybamm.max(sto_n_surf), ), Event( "Minimum positive particle surface stoichiometry", - pybamm.min(sto_p_surf) - 0.01, + pybamm.min(sto_p_surf) - tol, ), Event( "Maximum positive particle surface stoichiometry", - (1 - 0.01) - pybamm.max(sto_p_surf), + (1 - tol) - pybamm.max(sto_p_surf), + ), + # model does not capture electrolyte depletion, use the DFN instead + Event( + "Minimum negative electrode electrolyte stoichiometry", + pybamm.min(sto_e_n) - 0, + ), + Event( + "Minimum separator electrolyte stoichiometry", + pybamm.min(sto_e_sep) - 0, + ), + Event( + "Minimum positive electrode electrolyte stoichiometry", + pybamm.min(sto_e_p) - 0, ), ] @@ -136,12 +153,9 @@ def __init__( y_0 = Parameter("Maximum positive stoichiometry") # Grouped parameters - Q_th_p = Parameter("Measured cell capacity [A.s]") / (y_0 - y_100) - Q_th_n = Parameter("Measured cell capacity [A.s]") / (x_100 - x_0) - Q_e = Parameter("Reference electrolyte capacity [A.s]") - - tau_d_p = Parameter("Positive particle diffusion time scale [s]") - tau_d_n = Parameter("Negative particle diffusion time scale [s]") + Q_th_p = Parameter("Measured cell capacity [A.h]") * 3600 / (y_0 - y_100) + Q_th_n = Parameter("Measured cell capacity [A.h]") * 3600 / (x_100 - x_0) + Q_e = Parameter("Reference electrolyte capacity [A.h]") * 3600 tau_ct_p = Parameter("Positive electrode charge transfer time scale [s]") tau_ct_n = Parameter("Negative electrode charge transfer time scale [s]") @@ -156,9 +170,9 @@ def __init__( zeta_n = Parameter("Negative electrode relative porosity") zeta_p = Parameter("Positive electrode relative porosity") - tau_e_n = Parameter("Negative electrode electrolyte diffusion time scale [s]") - tau_e_sep = Parameter("Separator electrolyte diffusion time scale [s]") - tau_e_p = Parameter("Positive electrode electrolyte diffusion time scale [s]") + tau_e = Parameter("Electrolyte diffusion time scale [s]") + beta_n = Parameter("Negative electrode relative transport efficiency") + beta_p = Parameter("Positive electrode relative transport efficiency") ###################### # Input current (positive on discharge) @@ -197,74 +211,63 @@ def __init__( ###################### # Primary broadcasts are used to broadcast scalar quantities across a domain # into a vector of the right shape, for multiplying with other vectors - alpha = 0.5 # cathodic transfer coefficient - j0_n = ( - sto_n_surf**alpha * (sto_e_n * (1 - sto_n_surf)) ** (1 - alpha) / tau_ct_n + + # Overpotentials + eta_n = (v_s_n - U_n) + (2 * RT_F * (1 - t_plus)) * ( + pybamm.x_average(pybamm.log(sto_e_n)) - pybamm.log(sto_e_n) ) - j0_p = ( - sto_p_surf**alpha * (sto_e_p * (1 - sto_p_surf)) ** (1 - alpha) / tau_ct_p + eta_p = (v_s_p - U_p) + (2 * RT_F * (1 - t_plus)) * ( + pybamm.x_average(pybamm.log(sto_e_p)) - pybamm.log(sto_e_p) ) - if not include_double_layer: - # Assuming alpha = 0.5 - j_n = PrimaryBroadcast(I / (3 * Q_th_n), "negative electrode") - j_p = PrimaryBroadcast(-I / (3 * Q_th_p), "positive electrode") - eta_n = 2 * RT_F * pybamm.arcsinh(j_n / (2 * j0_n)) - eta_p = 2 * RT_F * pybamm.arcsinh(j_p / (2 * j0_p)) - v_s_n = pybamm.x_average(eta_n + U_n) - v_s_p = pybamm.x_average(eta_p + U_p) + + # Exchange rates + j_n = self.j(sto_n_surf, sto_e_n, eta_n / RT_F, "negative") / tau_ct_n + j_p = self.j(sto_p_surf, sto_e_p, eta_p / RT_F, "positive") / tau_ct_p ###################### # Double layer ###################### if include_double_layer: - # Additional variables - v_s_n = Variable("Negative particle surface voltage variable [V]") - v_s_p = Variable("Positive particle surface voltage variable [V]") - # Additional parameters C_p = Parameter("Positive electrode capacitance [F]") C_n = Parameter("Negative electrode capacitance [F]") - # Overpotentials - eta_n = (v_s_n - U_n) + (2 * RT_F * (1 - t_plus)) * ( - pybamm.x_average(pybamm.log(sto_e_n)) - pybamm.log(sto_e_n) - ) - eta_p = (v_s_p - U_p) + (2 * RT_F * (1 - t_plus)) * ( - pybamm.x_average(pybamm.log(sto_e_p)) - pybamm.log(sto_e_p) - ) - - # Exchange current - j_n = j0_n * ( - pybamm.exp((1 - alpha) * eta_n / RT_F) - - pybamm.exp(-alpha * eta_n / RT_F) - ) - j_p = j0_p * ( - pybamm.exp((1 - alpha) * eta_p / RT_F) - - pybamm.exp(-alpha * eta_p / RT_F) - ) - # Electrode surface potentials - self.rhs[v_s_n] = 1 / C_n * (I - 3 * Q_th_n * pybamm.x_average(j_n)) - self.rhs[v_s_p] = 1 / C_p * (-I - 3 * Q_th_p * pybamm.x_average(j_p)) - self.initial_conditions[v_s_n] = U_n_init - self.initial_conditions[v_s_p] = U_p_init + self.rhs[v_s_n] = (I - 3 * Q_th_n * pybamm.x_average(j_n)) / C_n + self.rhs[v_s_p] = (-I - 3 * Q_th_p * pybamm.x_average(j_p)) / C_p + else: + self.algebraic[v_s_n] = I - 3 * Q_th_n * pybamm.x_average(j_n) + self.algebraic[v_s_p] = -I - 3 * Q_th_p * pybamm.x_average(j_p) + + self.initial_conditions[v_s_n] = U_n_init + self.initial_conditions[v_s_p] = U_p_init ###################### # Particles ###################### # The div and grad operators will be converted to the appropriate matrix # multiplication at the discretisation stage - self.rhs[sto_n] = pybamm.div(pybamm.grad(sto_n) / tau_d_n) - self.rhs[sto_p] = pybamm.div(pybamm.grad(sto_p) / tau_d_p) + self.rhs[sto_n] = pybamm.div( + pybamm.grad(sto_n) / self.tau_d(sto_n, T, "negative") + ) + self.rhs[sto_p] = pybamm.div( + pybamm.grad(sto_p) / self.tau_d(sto_p, T, "positive") + ) # Boundary conditions must be provided for equations with spatial derivatives self.boundary_conditions[sto_n] = { "left": (Scalar(0), "Neumann"), - "right": (-tau_d_n * pybamm.x_average(j_n), "Neumann"), + "right": ( + -self.tau_d(sto_n_surf, T, "negative") * pybamm.x_average(j_n), + "Neumann", + ), } self.boundary_conditions[sto_p] = { "left": (Scalar(0), "Neumann"), - "right": (-tau_d_p * pybamm.x_average(j_p), "Neumann"), + "right": ( + -self.tau_d(sto_p_surf, T, "positive") * pybamm.x_average(j_p), + "Neumann", + ), } self.initial_conditions[sto_n] = sto_n_init @@ -274,45 +277,38 @@ def __init__( # Electrolyte ###################### self.rhs[sto_e_n] = ( - pybamm.div(pybamm.grad(sto_e_n) / tau_e_n - (t_plus * I / Q_e) * x_n / l_n) + pybamm.div( + pybamm.grad(sto_e_n) * beta_n / tau_e - (t_plus * I / Q_e) * x_n / l_n + ) + (3 / Q_e) * Q_th_n * j_n / l_n ) / zeta_n self.rhs[sto_e_sep] = pybamm.div( - pybamm.grad(sto_e_sep) / tau_e_sep - t_plus * I / Q_e + pybamm.grad(sto_e_sep) / tau_e - t_plus * I / Q_e ) self.rhs[sto_e_p] = ( pybamm.div( - pybamm.grad(sto_e_p) / tau_e_p - (t_plus * I / Q_e) * (1 - x_p) / l_p + pybamm.grad(sto_e_p) * beta_p / tau_e + - (t_plus * I / Q_e) * (1 - x_p) / l_p ) + (3 / Q_e) * Q_th_p * j_p / l_p ) / zeta_p self.boundary_conditions[sto_e_n] = { "left": (Scalar(0), "Neumann"), - "right": ( - tau_e_n * pybamm.boundary_gradient(sto_e_sep, "left") / tau_e_sep, - "Neumann", - ), + "right": (pybamm.boundary_gradient(sto_e_sep, "left") / beta_n, "Neumann"), } self.boundary_conditions[sto_e_sep] = { "left": (pybamm.boundary_value(sto_e_n, "right"), "Dirichlet"), "right": (pybamm.boundary_value(sto_e_p, "left"), "Dirichlet"), } self.boundary_conditions[sto_e_p] = { - "left": ( - tau_e_p * pybamm.boundary_gradient(sto_e_sep, "right") / tau_e_sep, - "Neumann", - ), + "left": (pybamm.boundary_gradient(sto_e_sep, "right") / beta_p, "Neumann"), "right": (Scalar(0), "Neumann"), } - self.initial_conditions[sto_e_n] = PrimaryBroadcast( - Scalar(1), "negative electrode" - ) - self.initial_conditions[sto_e_sep] = PrimaryBroadcast(Scalar(1), "separator") - self.initial_conditions[sto_e_p] = PrimaryBroadcast( - Scalar(1), "positive electrode" - ) + self.initial_conditions[sto_e_n] = Scalar(1) + self.initial_conditions[sto_e_sep] = Scalar(1) + self.initial_conditions[sto_e_p] = Scalar(1) ###################### # Cell voltage @@ -328,6 +324,52 @@ def __init__( Event("Maximum voltage [V]", self.param.voltage_high_cut - V), ] + ###################### + # Voltage components + ###################### + # Include the following variables to enable plotting via PyBaMM's plot_voltage_components + ocp_n_bulk = self.U( + pybamm.x_average(Q_th_n * pybamm.r_average(sto_n)) + / pybamm.x_average(Q_th_n), + "negative", + ) + ocp_p_bulk = self.U( + pybamm.x_average(Q_th_p * pybamm.r_average(sto_p)) + / pybamm.x_average(Q_th_p), + "positive", + ) + voltage_components = { + "Battery voltage [V]": V, + "Battery open-circuit voltage [V]": ocp_p_bulk - ocp_n_bulk, + "Battery particle concentration overpotential [V]": ( + (pybamm.x_average(self.U(sto_p_surf, "positive")) - ocp_p_bulk) + - (pybamm.x_average(self.U(sto_n_surf, "negative")) - ocp_n_bulk) + ), + "X-averaged battery reaction overpotential [V]": pybamm.x_average(eta_p) + - pybamm.x_average(eta_n), + "X-averaged battery concentration overpotential [V]": eta_e, + "X-averaged battery electrolyte ohmic losses [V]": Scalar(0), + "X-averaged battery solid phase ohmic losses [V]": Scalar(0), + "Contact overpotential [V]": R0 * I, # includes Ohmic losses in this model + # and split by electrode + "Negative electrode bulk open-circuit potential [V]": ocp_n_bulk, + "Positive electrode bulk open-circuit potential [V]": ocp_p_bulk, + "Negative particle concentration overpotential [V]": pybamm.x_average( + self.U(sto_n_surf, "negative") + ) + - ocp_n_bulk, + "Positive particle concentration overpotential [V]": pybamm.x_average( + self.U(sto_p_surf, "positive") + ) + - ocp_p_bulk, + "X-averaged negative electrode reaction overpotential [V]" + "": pybamm.x_average(eta_n), + "X-averaged positive electrode reaction overpotential [V]" + "": pybamm.x_average(eta_p), + "X-averaged battery negative solid phase ohmic losses [V]": Scalar(0), + "X-averaged battery positive solid phase ohmic losses [V]": Scalar(0), + } + ###################### # (Some) variables ###################### @@ -368,33 +410,25 @@ def __init__( - pybamm.log(pybamm.concatenation(sto_e_n, sto_e_sep, sto_e_p)) ), "Time [s]": pybamm_t, + "Time [h]": pybamm_t / 3600, "Current [A]": I, "Current variable [A]": I, # for compatibility with pybamm.Experiment "Discharge capacity [A.h]": Q, "Throughput capacity [A.h]": Qt, "Voltage [V]": V, "Voltage expression [V]": V, # for compatibility with "voltage as a state" - "Battery voltage [V]": V, "Open-circuit voltage [V]": U_p - U_n, + **voltage_components, } def U(self, sto, domain): """ - Dimensional open-circuit potential [V], calculated as U(x) = U_ref(x). + Dimensional open-circuit potential [V]. Credit: PyBaMM """ - # bound stoichiometry between tol and 1-tol. Adding 1/sto + 1/(sto-1) later - # will ensure that ocp goes to +- infinity if sto goes into that region - # anyway Domain = domain.capitalize() - tol = pybamm.settings.tolerances["U__c_s"] - sto = pybamm.maximum(pybamm.minimum(sto, 1 - tol), tol) inputs = {f"{Domain} particle surface stoichiometry": sto} - u_ref = FunctionParameter(f"{Domain} electrode OCP [V]", inputs) - - # add a term to ensure that the OCP goes to infinity at 0 and -infinity at 1 - # this will not affect the OCP for most values of sto - out = u_ref + 1e-6 * (1 / sto + 1 / (sto - 1)) + out = FunctionParameter(f"{Domain} electrode OCP [V]", inputs) if domain == "negative": out.print_name = r"U_\mathrm{n}(c^\mathrm{surf}_\mathrm{s,n})" @@ -402,6 +436,28 @@ def U(self, sto, domain): out.print_name = r"U_\mathrm{p}(c^\mathrm{surf}_\mathrm{s,p})" return out + def tau_d(self, sto, T, domain): + """ + Dimensional solid-state diffusion time scale [s]. + """ + Domain = domain.capitalize() + inputs = {f"{Domain} particle surface stoichiometry": sto, "Temperature [K]": T} + return FunctionParameter(f"{Domain} particle diffusion time scale [s]", inputs) + + def j(self, sto_surf, sto_e, eta_RT_F, domain): + """ + Dimensionless exchange rate. + """ + Domain = domain.capitalize() + inputs = { + f"{Domain} particle surface stoichiometry": sto_surf, + f"{Domain} electrode electrolyte stoichiometry": sto_e, + f"{Domain} electrode dimensionless overpotential": eta_RT_F, + } + return FunctionParameter( + f"{Domain} electrode dimensionless exchange rate", inputs + ) + @property def default_parameter_values(self) -> ParameterValues: param = ParameterValues("Chen2020") @@ -418,9 +474,9 @@ def default_parameter_values(self) -> ParameterValues: @property def default_quick_plot_variables(self): return [ - "Negative particle surface stoichiometry", + "Negative particle stoichiometry", "Electrolyte stoichiometry", - "Positive particle surface stoichiometry", + "Positive particle stoichiometry", "Current [A]", { "Negative electrode potential [V]", @@ -521,6 +577,7 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal # Unpack physical parameters F = pybamm.constants.F.value + T = param["Ambient temperature [K]"] alpha_p = param["Positive electrode active material volume fraction"] alpha_n = param["Negative electrode active material volume fraction"] c_max_p = param["Maximum concentration in positive electrode [mol.m-3]"] @@ -537,8 +594,6 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal b_n = param["Negative electrode Bruggeman coefficient (electrolyte)"] Cdl_p = param["Positive electrode double-layer capacity [F.m-2]"] Cdl_n = param["Negative electrode double-layer capacity [F.m-2]"] - m_p = 3.42e-6 # (A/m2)(m3/mol)**1.5 - m_n = 6.48e-7 # (A/m2)(m3/mol)**1.5 sigma_p = ( param["Positive electrode conductivity [S.m-1]"] * alpha_p ** param["Positive electrode Bruggeman coefficient (electrode)"] @@ -555,7 +610,19 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal epsilon_sep = param["Separator porosity"] b_sep = param["Separator Bruggeman coefficient (electrolyte)"] t_plus = param["Cation transference number"] - sigma_e = param["Electrolyte conductivity [S.m-1]"] # (ce0, T) + kappa_e = param["Electrolyte conductivity [S.m-1]"] # (ce0, T) + + # Get reference exchange current density [A.m-2] + j0_p = param.evaluate( + param["Positive electrode exchange-current density [A.m-2]"]( + ce0, c_max_p / 2, c_max_p, T + ) + ) + j0_n = param.evaluate( + param["Negative electrode exchange-current density [A.m-2]"]( + ce0, c_max_n / 2, c_max_n, T + ) + ) # Compute the cell area and thickness A = param["Electrode height [m]"] * param["Electrode width [m]"] @@ -566,7 +633,7 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal L_p / (3 * epsilon_p**b_p) + L_s / (epsilon_sep**b_sep) + L_n / (3 * epsilon_n**b_n) - ) / (sigma_e * A) + ) / (kappa_e * A) Rs = (L_p / sigma_p + L_n / sigma_n) / (3 * A) R0 = Re + Rs + param["Contact resistance [Ohm]"] @@ -578,8 +645,8 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal soc_init = (sto_p_init - y_0) / (y_100 - y_0) # Compute the capacity within the stoichiometry limits - Q_th_p = F * alpha_p * c_max_p * L_p * A - Q_th_n = F * alpha_n * c_max_n * L_n * A + Q_th_p = F * alpha_p * c_max_p * L_p * A / 3600 + Q_th_n = F * alpha_n * c_max_n * L_n * A / 3600 Q_meas_p = (y_0 - y_100) * Q_th_p Q_meas_n = (x_100 - x_0) * Q_th_n if abs(Q_meas_n / Q_meas_p - 1) > 1e-6: @@ -589,20 +656,27 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal # Grouped parameters Q_meas = (Q_meas_n + Q_meas_p) / 2 - Q_e = F * epsilon_sep * ce0 * L * A + Q_e = F * epsilon_sep * ce0 * L * A / 3600 zeta_p = epsilon_p / epsilon_sep zeta_n = epsilon_n / epsilon_sep - tau_d_p = R_p**2 / D_p - tau_d_n = R_n**2 / D_n + try: + tau_d_p = R_p**2 / D_p + except TypeError: + tau_d_p = FunctionalDiffusionTime(R_p**2, D_p, c_max_p) + + try: + tau_d_n = R_n**2 / D_n + except TypeError: + tau_d_n = FunctionalDiffusionTime(R_n**2, D_n, c_max_n) - tau_e_p = epsilon_sep * L**2 / (epsilon_p**b_p * De) - tau_e_n = epsilon_sep * L**2 / (epsilon_n**b_n * De) - tau_e_sep = epsilon_sep * L**2 / (epsilon_sep**b_sep * De) + tau_e = epsilon_sep * L**2 / (epsilon_sep**b_sep * De) + beta_p = epsilon_p**b_p / epsilon_sep**b_sep + beta_n = epsilon_n**b_n / epsilon_sep**b_sep - tau_ct_p = F * R_p / (m_p * np.sqrt(ce0)) - tau_ct_n = F * R_n / (m_n * np.sqrt(ce0)) + tau_ct_p = c_max_p * F * R_p / (2 * j0_p) + tau_ct_n = c_max_n * F * R_n / (2 * j0_n) C_p = 3 * alpha_p * Cdl_p * L_p * A / R_p C_n = 3 * alpha_n * Cdl_n * L_n * A / R_n @@ -613,8 +687,8 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal parameter_dictionary = { "Nominal cell capacity [A.h]": param["Nominal cell capacity [A.h]"], "Current function [A]": param["Current function [A]"], - "Ambient temperature [K]": param["Ambient temperature [K]"], - "Initial temperature [K]": param["Ambient temperature [K]"], + "Ambient temperature [K]": T, + "Initial temperature [K]": T, "Initial SoC": soc_init, "Minimum negative stoichiometry": x_0, "Maximum negative stoichiometry": x_100, @@ -624,15 +698,17 @@ def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterVal "Upper voltage cut-off [V]": param["Upper voltage cut-off [V]"], "Positive electrode OCP [V]": param["Positive electrode OCP [V]"], "Negative electrode OCP [V]": param["Negative electrode OCP [V]"], - "Measured cell capacity [A.s]": Q_meas, - "Reference electrolyte capacity [A.s]": Q_e, + "Measured cell capacity [A.h]": Q_meas, + "Reference electrolyte capacity [A.h]": Q_e, "Positive electrode relative porosity": zeta_p, "Negative electrode relative porosity": zeta_n, "Positive particle diffusion time scale [s]": tau_d_p, "Negative particle diffusion time scale [s]": tau_d_n, - "Positive electrode electrolyte diffusion time scale [s]": tau_e_p, - "Negative electrode electrolyte diffusion time scale [s]": tau_e_n, - "Separator electrolyte diffusion time scale [s]": tau_e_sep, + "Electrolyte diffusion time scale [s]": tau_e, + "Positive electrode relative transport efficiency": beta_p, + "Negative electrode relative transport efficiency": beta_n, + "Positive electrode dimensionless exchange rate": GroupedSPMe.symmetric_butler_volmer, + "Negative electrode dimensionless exchange rate": GroupedSPMe.symmetric_butler_volmer, "Positive electrode charge transfer time scale [s]": tau_ct_p, "Negative electrode charge transfer time scale [s]": tau_ct_n, "Positive electrode capacitance [F]": C_p, diff --git a/pybop/models/lithium_ion/sp_diffusion.py b/pybop/models/lithium_ion/sp_diffusion.py deleted file mode 100644 index 6798aec17..000000000 --- a/pybop/models/lithium_ion/sp_diffusion.py +++ /dev/null @@ -1,311 +0,0 @@ -import pybamm -from pybamm import ( - Event, - FunctionParameter, - Parameter, - ParameterValues, - PrimaryBroadcast, - Scalar, - SpatialVariable, - Variable, -) -from pybamm import t as pybamm_t - -from pybop.models.lithium_ion.base_model import BaseGroupedModel -from pybop.models.lithium_ion.utils import InverseOCV - - -class SPDiffusion(BaseGroupedModel): - """ - Diffusion model for a single, spherical particle representing a half-cell for GITT. - - Note: the working electrode is the positive electrode. - - Parameters - ---------- - name : str, optional - The name of the model. - **model_kwargs : optional - Valid PyBaMM model option keys and their values, for example: - options : dict, optional - A dictionary of options to customise the behaviour of the PyBaMM model. - build : bool, optional - If True, the model is built upon creation (default: False). - """ - - def __init__(self, name="Single Particle Diffusion Model", **model_kwargs): - super().__init__(name=name, **model_kwargs) - - ###################### - # Variables - ###################### - # Variables that depend on time only are created without a domain - Q = Variable("Discharge capacity [A.h]") - Qt = Variable("Throughput capacity [A.h]") - - # Variables that vary spatially are created with a domain - sto = Variable("Particle stoichiometry", domain="particle") - sto_surf = pybamm.surf(sto) - - # Events specify points at which a solution should terminate - self.events += [ - Event( - "Minimum particle surface stoichiometry", - pybamm.min(sto_surf) - 0.01, - ), - Event( - "Maximum particle surface stoichiometry", - (1 - 0.01) - pybamm.max(sto_surf), - ), - ] - - ###################### - # Parameters - ###################### - # Parameters are purely symbolic at this stage, and will be set by the - # `ParameterValues` class when the model is processed. - - Q_th = Parameter("Theoretical electrode capacity [A.s]") - - sto_init = Parameter("Initial stoichiometry") - - ###################### - # Input current (positive on discharge) - ###################### - I = self.param.current_with_time - - ###################### - # State of Charge - ###################### - # The `rhs` dictionary contains differential equations, with the key being the - # variable in the d/dt - self.rhs[Q] = I / 3600 - self.rhs[Qt] = abs(I) / 3600 - # Initial conditions must be provided for the ODEs - self.initial_conditions[Q] = Scalar(0) - self.initial_conditions[Qt] = Scalar(0) - - ###################### - # Diffusion within the particle - ###################### - # The div and grad operators will be converted to the appropriate matrix - # multiplication at the discretisation stage - self.rhs[sto] = pybamm.div(pybamm.grad(sto) / self.tau_d(sto)) - - # Boundary conditions must be provided for equations with spatial derivatives - j = -I / (3 * Q_th) - self.boundary_conditions[sto] = { - "left": (Scalar(0), "Neumann"), - "right": (-self.tau_d(sto_surf) * j, "Neumann"), - } - - self.initial_conditions[sto] = sto_init - - ###################### - # Cell voltage - ###################### - U = self.U(sto_surf) - V = U - self.R0(sto_surf) * I - - # Save the initial OCV - self.param.ocv_init = self.U(sto_init) - - ###################### - # (Some) variables - ###################### - self.variables = { - "Particle stoichiometry": sto, - "Particle surface stoichiometry": PrimaryBroadcast(sto_surf, "particle"), - "Time [s]": pybamm_t, - "Current [A]": I, - "Current variable [A]": I, # for compatibility with pybamm.Experiment - "Discharge capacity [A.h]": Q, - "Throughput capacity [A.h]": Qt, - "Voltage [V]": V, - "Voltage expression [V]": V, # for compatibility with "voltage as a state" - "Open-circuit voltage [V]": U, - } - - def U(self, sto): - """ - Dimensional open-circuit potential [V], calculated as U(x) = U_ref(x). - Credit: PyBaMM - """ - # bound stoichiometry between tol and 1-tol. Adding 1/sto + 1/(sto-1) later - # will ensure that ocp goes to +- infinity if sto goes into that region - # anyway - tol = pybamm.settings.tolerances["U__c_s"] - sto = pybamm.maximum(pybamm.minimum(sto, 1 - tol), tol) - inputs = {"Particle surface stoichiometry": sto} - u_ref = FunctionParameter("Electrode OCP [V]", inputs) - - # add a term to ensure that the OCP goes to infinity at 0 and -infinity at 1 - # this will not affect the OCP for most values of sto - out = u_ref + 1e-6 * (1 / sto + 1 / (sto - 1)) - - out.print_name = r"U(c^\mathrm{surf}_\mathrm{s})" - return out - - def tau_d(self, sto): - """ - Diffusion time scale [s] for lithium in the particles. - """ - inputs = {"Particle surface stoichiometry": sto} - return FunctionParameter("Particle diffusion time scale [s]", inputs) - - def R0(self, sto): - """ - Series resistance [Ohm]. - """ - inputs = {"Particle surface stoichiometry": sto} - return FunctionParameter("Series resistance [Ohm]", inputs) - - @property - def default_parameter_values(self) -> ParameterValues: - param = ParameterValues("Xu2019") - return self.create_grouped_parameters(param) - - @property - def default_quick_plot_variables(self): - return [ - "Particle stoichiometry", - "Particle surface stoichiometry", - "Current [A]", - {"Open-circuit voltage [V]", "Voltage [V]"}, - ] - - @property - def default_var_pts(self): - r = SpatialVariable("r", domain=["particle"], coord_sys="spherical polar") - return {r: 20} - - @property - def default_geometry(self): - r = SpatialVariable("r", domain=["particle"], coord_sys="spherical polar") - return {"particle": {r: {"min": 0, "max": 1}}} - - @property - def default_submesh_types(self): - return {"particle": pybamm.Uniform1DSubMesh} - - @property - def default_spatial_methods(self): - return {"particle": pybamm.FiniteVolume()} - - @staticmethod - def create_grouped_parameters(parameter_values: ParameterValues) -> ParameterValues: - """ - Create a parameter set for the Single Particle Diffusion Model from a - PyBaMM lithium-ion ParameterValues object. - - Note: the working electrode is the positive electrode. - - Parameters - ---------- - parameter_values : pybamm.ParameterValues - Parameters and their corresponding values. - - Returns - ------- - parameter_values : pybamm.ParameterValues - A new set of parameters and their values. - """ - param = parameter_values - - # Unpack physical parameters - F = pybamm.constants.F.value - alpha = param["Positive electrode active material volume fraction"] - c_max = param["Maximum concentration in positive electrode [mol.m-3]"] - L = param["Positive electrode thickness [m]"] - R = param["Positive particle radius [m]"] - D = param["Positive particle diffusivity [m2.s-1]"] - sto_init = ( - param["Initial concentration in positive electrode [mol.m-3]"] / c_max - ) - ocp = param["Positive electrode OCP [V]"] - - # Compute the cell area - A = param["Electrode height [m]"] * param["Electrode width [m]"] - - # Grouped parameters - Q_th = F * alpha * c_max * L * A - tau_d = R**2 / D - - parameter_dictionary = { - "Nominal cell capacity [A.h]": param["Nominal cell capacity [A.h]"], - "Current function [A]": param["Current function [A]"], - "Initial stoichiometry": sto_init, - "Electrode OCP [V]": ocp, - "Theoretical electrode capacity [A.s]": Q_th, - "Particle diffusion time scale [s]": tau_d, - "Series resistance [Ohm]": 1, - } - parameter_values = ParameterValues(values=parameter_dictionary) - parameter_values._set_initial_state = SPDiffusion.set_initial_state # noqa: SLF001 - return parameter_values - - @staticmethod - def set_initial_state( - initial_value, - parameter_values, - direction=None, - param=None, - inplace=True, - options=None, - inputs=None, - tol=1e-6, - ): - """ - Set the value of the initial state of charge. - - Parameters - ---------- - initial_value : float - Target initial value. - If float, interpreted as SOC, must be between 0 and 1. - If string e.g. "4 V", interpreted as voltage, must be between V_min and V_max. - parameter_values : :class:`pybamm.ParameterValues` - Parameters and their corresponding values. - param : :class:`pybamm.LithiumIonParameters`, optional - The symbolic parameter set to use for the simulation. - If not provided, the default parameter set will be used. - inplace: bool, optional - If True, replace the parameters values in place. Otherwise, return a new set of - parameter values. Default is True. - options : dict-like, optional - A dictionary of options to be passed to the model, see - :class:`pybamm.BatteryModelOptions`. - inputs : dict, optional - A dictionary of input parameters to pass to the model when solving. - tol : float, optional - The tolerance for the solver used to compute the initial stoichiometries. - A lower value results in higher precision but may increase computation time. - Default is 1e-6. - """ - parameter_values = parameter_values if inplace else parameter_values.copy() - - if isinstance(initial_value, str) and initial_value.endswith("V"): - V_init = float(initial_value[:-1]) - - def ocv_function(sto): - U = FunctionParameter( - "Electrode OCP [V]", {"Particle stoichiometry": sto} - ) - - return parameter_values.evaluate(U, inputs=inputs) - - inverse_ocv = InverseOCV(ocv_function) - sto = inverse_ocv(V_init) - - elif isinstance(initial_value, int | float): - sto = initial_value - - else: - raise ValueError("Initial value must be a float or a string ending in 'V'.") - - if not 0 <= sto <= 1: - raise ValueError("Initial stoichiometry should be between 0 and 1.") - - parameter_values["Initial stoichiometry"] = sto - - return parameter_values diff --git a/pybop/models/lithium_ion/utils.py b/pybop/models/lithium_ion/utils.py index 6448306b0..e5402429c 100644 --- a/pybop/models/lithium_ion/utils.py +++ b/pybop/models/lithium_ion/utils.py @@ -133,7 +133,7 @@ def solve_batch(self, inputs, calculate_sensitivities: bool = False): for x in inputs: diff = np.abs(ocv_function(x["Root"]) - self.ocv_value) sol = Solution() - sol.set_solution_variable("Difference", data=np.asarray([diff])) + sol.set_solution_variable("Difference", data=np.atleast_1d(diff)) solutions.append(sol) return solutions diff --git a/tests/integration/models/test_gitt_models.py b/tests/integration/models/test_gitt_models.py index 66f5d5f96..7c690c5a8 100644 --- a/tests/integration/models/test_gitt_models.py +++ b/tests/integration/models/test_gitt_models.py @@ -14,8 +14,8 @@ # Parameter configurations DIFFUSION_PARAMS = [ - ("Theoretical electrode capacity [A.s]", 10), - ("Particle diffusion time scale [s]", 2000), + ("Theoretical electrode capacity [A.h]", 0.003), + ("Positive particle diffusion time scale [s]", 2000), ] @@ -28,8 +28,8 @@ class TestGITTModels: @pytest.fixture( params=[ - pybop.lithium_ion.WeppnerHuggins(), - pybop.lithium_ion.SPDiffusion(), + pybop.li_half_cell.WeppnerHuggins(), + pybop.li_half_cell.SPDiffusion(), ], scope="module", ) diff --git a/tests/integration/test_applications.py b/tests/integration/test_applications.py index 248305a41..f3173c9d6 100644 --- a/tests/integration/test_applications.py +++ b/tests/integration/test_applications.py @@ -48,13 +48,14 @@ def charge_dataset(self, parameter_values): ) def test_interpolant(self, parameter_values, discharge_dataset): - parameter_values = pybop.lithium_ion.SPDiffusion.create_grouped_parameters( + parameter_values = pybop.li_half_cell.SPDiffusion.create_grouped_parameters( parameter_values ) - parameter_values["Electrode OCP [V]"] = pybop.Interpolant( + parameter_values["Positive electrode OCP [V]"] = pybop.Interpolant( discharge_dataset["Stoichiometry"], discharge_dataset["Voltage [V]"] ) - model = pybop.lithium_ion.SPDiffusion(build=True) + parameter_values.set_initial_state(0.9) + model = pybop.li_half_cell.SPDiffusion(build=True) t_eval = np.linspace(0, 10, 100) solution = pybamm.Simulation(model, parameter_values=parameter_values).solve( t_eval=t_eval, t_interp=t_eval @@ -186,25 +187,25 @@ def pulse_data(self, half_cell_model, half_cell_parameter_values): def test_gitt_pulse_fit( self, half_cell_model, half_cell_parameter_values, pulse_data ): - parameter_values = pybop.lithium_ion.SPDiffusion.create_grouped_parameters( + parameter_values = pybop.li_half_cell.SPDiffusion.create_grouped_parameters( half_cell_parameter_values ) - diffusion_time = parameter_values["Particle diffusion time scale [s]"] + diffusion_time = parameter_values["Positive particle diffusion time scale [s]"] gitt_fit = pybop.GITTPulseFit(parameter_values=parameter_values) gitt_result = gitt_fit(gitt_pulse=pulse_data) np.testing.assert_allclose( - gitt_result.best_inputs["Particle diffusion time scale [s]"], + gitt_result.best_inputs["Positive particle diffusion time scale [s]"], diffusion_time, rtol=5e-2, ) def test_gitt_fit(self, half_cell_model, half_cell_parameter_values, pulse_data): - parameter_values = pybop.lithium_ion.SPDiffusion.create_grouped_parameters( + parameter_values = pybop.li_half_cell.SPDiffusion.create_grouped_parameters( half_cell_parameter_values ) - diffusion_time = parameter_values["Particle diffusion time scale [s]"] + diffusion_time = parameter_values["Positive particle diffusion time scale [s]"] with pytest.raises( ValueError, match="The initial current in the pulse dataset must be zero." @@ -224,7 +225,7 @@ def test_gitt_fit(self, half_cell_model, half_cell_parameter_values, pulse_data) gitt_parameter_data = gitt_fit() np.testing.assert_allclose( - gitt_parameter_data["Particle diffusion time scale [s]"], + gitt_parameter_data["Positive particle diffusion time scale [s]"], np.asarray([diffusion_time]), rtol=5e-2, ) diff --git a/tests/unit/test_models.py b/tests/unit/test_models.py index f3bac60b1..28db1fae7 100644 --- a/tests/unit/test_models.py +++ b/tests/unit/test_models.py @@ -1,9 +1,13 @@ import numpy as np import pybamm import pytest +from matplotlib.axes._axes import Axes +from matplotlib.figure import Figure import pybop +GROUPED_MODEL = pybop.lithium_ion.GroupedSPM | pybop.lithium_ion.GroupedSPMe + class TestModels: """ @@ -16,13 +20,14 @@ class TestModels: params=[ pybop.ExponentialDecayModel(), pybop.lithium_ion.CellTemperature(), - pybop.lithium_ion.WeppnerHuggins(), - pybop.lithium_ion.SPDiffusion(), pybop.lithium_ion.GroupedSPM(), pybop.lithium_ion.GroupedSPM(options={"surface form": "differential"}), pybop.lithium_ion.GroupedSPMe(), pybop.lithium_ion.GroupedSPMe(options={"surface form": "differential"}), + pybop.li_half_cell.WeppnerHuggins(), + pybop.li_half_cell.SPDiffusion(), ], + ids=lambda val: f"{type(val).__name__}", scope="module", ) def model(self, request): @@ -47,11 +52,24 @@ def test_model_simulation(self, model): fig = solution.plot() assert isinstance(fig, pybamm.QuickPlot) + if isinstance(model, GROUPED_MODEL): + for split in [False, True]: + fig, ax = solution.plot_voltage_components(split_by_electrode=split) + assert isinstance(fig, Figure) + assert isinstance(ax, Axes) + + if not isinstance(model, pybop.ExponentialDecayModel): + experiment = pybamm.Experiment(["Discharge at 1C for 1 minute"]) + solution = pybamm.Simulation(model, experiment=experiment).solve( + calc_esoh=False + ) + assert solution["Time [s]"].data[-1] == 60 + def test_set_initial_state(self, model): if isinstance(model, pybop.ExponentialDecayModel): pass # Only testing the battery models for now - elif isinstance(model, pybop.lithium_ion.WeppnerHuggins): + elif isinstance(model, pybop.li_half_cell.WeppnerHuggins): param = model.default_parameter_values with pytest.raises( ValueError, @@ -60,7 +78,7 @@ def test_set_initial_state(self, model): param.set_initial_state(0.5) else: - if isinstance(model, pybop.lithium_ion.SPDiffusion): + if isinstance(model, pybop.li_half_cell.SPDiffusion): initial_state = "Initial stoichiometry" else: initial_state = "Initial SoC" @@ -69,7 +87,7 @@ def test_set_initial_state(self, model): param.set_initial_state(0.5) assert param[initial_state] == 0.5 - param.set_initial_state("2.8 V") + param.set_initial_state("3.8 V") assert 0 <= param[initial_state] <= 1 with pytest.raises( @@ -81,12 +99,33 @@ def test_set_initial_state(self, model): with pytest.raises(ValueError, match="should be between 0 and 1."): param.set_initial_state(-1) - if not isinstance(model, pybop.lithium_ion.SPDiffusion): + if not isinstance(model, pybop.li_half_cell.SPDiffusion): with pytest.raises( ValueError, match=r"V is outside the voltage limits" ): param.set_initial_state("-1 V") + def test_alternative_functions(self, model): + if isinstance(model, GROUPED_MODEL): + parameter_values = model.default_parameter_values + + # Use alternative Butler-Volmer kinetics + parameter_values["Negative electrode charge transfer coefficient"] = 0.6 + parameter_values["Negative electrode dimensionless exchange rate"] = ( + model.get_asymmetric_butler_volmer("negative") + ) + parameter_values["Positive electrode charge transfer coefficient"] = 0.6 + parameter_values["Positive electrode charge transfer ideality factor"] = 0.6 + parameter_values["Positive electrode dimensionless exchange rate"] = ( + model.get_multiphase_butler_volmer("positive") + ) + + t_eval = np.linspace(0, 10, 11) + solution = pybamm.Simulation( + model, parameter_values=parameter_values + ).solve(t_eval=t_eval, t_interp=t_eval) + np.testing.assert_allclose(solution["Time [s]"].data, t_eval) + class TestModelUtils: """