diff --git a/examples/scripts/SPMe_SOC.py b/examples/scripts/SPMe_SOC.py new file mode 100644 index 0000000000..411eb50f67 --- /dev/null +++ b/examples/scripts/SPMe_SOC.py @@ -0,0 +1,131 @@ +# +# Example to show the state of charge of a battery using the SPMe model +# Initial conditions are specified to start each electrode in 1/2 charged +# state. A charge and discharge are performed with current chosen to be +# 1C rate when electrode dimensions are euqal. +# Coulomb counting is performed to calculate the capacity of the +# battery within the operating voltage limits and maximum particle concs. +# The anode thickenss is varied to highlight the importance of electrode +# sizing to enable full lithium utilization +# +import pybamm +import numpy as np +import matplotlib.pyplot as plt + +plt.close("all") +pybamm.set_logging_level(30) + +factor = 6.38 +capacities = [] +specific_capacities = [] +l_p = 1e-4 +thicknesses = np.linspace(1.0, 2.5, 11) * l_p +for l_n in thicknesses: + e_ratio = np.around(l_n / l_p, 3) + # Dimensions + h = 0.137 + w = 0.207 / factor + A = h * w + l_s = 2.5e-5 + l1d = (l_n + l_p + l_s) + vol = h * w * l1d + vol_cm3 = vol * 1e6 + tot_cap = 0.0 + tot_time = 0.0 + fig, axes = plt.subplots(1, 2, sharey=True) + I_mag = 1.01 / factor + print('*' * 30) + for enum, I_app in enumerate([-1.0, 1.0]): + I_app *= I_mag + # load model + model = pybamm.lithium_ion.SPMe() + # create geometry + geometry = model.default_geometry + # load parameter values and process model and geometry + param = model.default_parameter_values + param.update( + {"Electrode height [m]": h, + "Electrode width [m]": w, + "Negative electrode thickness [m]": l_n, + "Positive electrode thickness [m]": l_p, + "Separator thickness [m]": l_s, + "Lower voltage cut-off [V]": 2.8, + "Upper voltage cut-off [V]": 4.7, + "Maximum concentration in negative electrode [mol.m-3]": 25000, + "Maximum concentration in positive electrode [mol.m-3]": 50000, + "Initial concentration in negative electrode [mol.m-3]": 12500, + "Initial concentration in positive electrode [mol.m-3]": 25000, + "Negative electrode surface area density [m-1]": 180000.0, + "Positive electrode surface area density [m-1]": 150000.0, + "Typical current [A]": I_app, + } + ) + param.process_model(model) + param.process_geometry(geometry) + s_var = pybamm.standard_spatial_vars + var_pts = {s_var.x_n: 5, s_var.x_s: 5, s_var.x_p: 5, + s_var.r_n: 5, s_var.r_p: 10} + # set mesh + mesh = pybamm.Mesh(geometry, model.default_submesh_types, var_pts) + # discretise model + disc = pybamm.Discretisation(mesh, model.default_spatial_methods) + disc.process_model(model) + # solve model + t_eval = np.linspace(0, 0.2, 100) + sol = model.default_solver.solve(model, t_eval) + var = "Positive electrode average extent of lithiation" + xpext = pybamm.ProcessedVariable(model.variables[var], + sol.t, sol.y, mesh=mesh) + var = "Negative electrode average extent of lithiation" + xnext = pybamm.ProcessedVariable(model.variables[var], + sol.t, sol.y, mesh=mesh) + var = "X-averaged positive particle surface concentration" + xpsurf = pybamm.ProcessedVariable(model.variables[var], + sol.t, sol.y, mesh=mesh) + var = "X-averaged negative particle surface concentration" + xnsurf = pybamm.ProcessedVariable(model.variables[var], + sol.t, sol.y, mesh=mesh) + time = pybamm.ProcessedVariable(model.variables["Time [h]"], + sol.t, sol.y, mesh=mesh) + # Coulomb counting + time_hours = time(sol.t) + dc_time = np.around(time_hours[-1], 3) + # Capacity mAh + cap = np.absolute(I_app * 1000 * dc_time) + cap_time = np.absolute(I_app * 1000 * time_hours) + axes[enum].plot(cap_time, + xnext(sol.t), 'r-', label='Average Neg') + axes[enum].plot(cap_time, + xpext(sol.t), 'b-', label='Average Pos') + axes[enum].plot(cap_time, + xnsurf(sol.t), 'r--', label='Surface Neg') + axes[enum].plot(cap_time, + xpsurf(sol.t), 'b--', label='Surface Pos') + axes[enum].set_xlabel('Capacity [mAh]') + handles, labels = axes[enum].get_legend_handles_labels() + axes[enum].legend(handles, labels) + if I_app < 0.0: + axes[enum].set_ylabel('Extent of Lithiation, Elecrode Ratio: ' + + str(e_ratio)) + axes[enum].title.set_text('Charge') + else: + axes[enum].title.set_text('Discharge') + print('Applied Current', I_app, 'A', 'Time', + dc_time, 'hrs', 'Capacity', cap, 'mAh') + tot_cap += cap + tot_time += dc_time + + print('Anode : Cathode thickness', e_ratio) + print('Total Charge/Discharge Time', tot_time, 'hrs') + print('Total Capacity', np.around(tot_cap, 3), 'mAh') + specific_cap = np.around(tot_cap, 3) / vol_cm3 + print('Total Capacity', specific_cap, 'mAh.cm-3') + capacities.append(tot_cap) + specific_capacities.append(specific_cap) + +fig, (ax1, ax2) = plt.subplots(2, 1, sharex=True) +ax1.plot(thicknesses / l_p, capacities) +ax2.plot(thicknesses / l_p, specific_capacities) +ax1.set_ylabel('Capacity [mAh]') +ax2.set_ylabel('Specific Capacity [mAh.cm-3]') +ax2.set_xlabel('Anode : Cathode thickness') diff --git a/pybamm/__init__.py b/pybamm/__init__.py index 6db639ad3a..c4bb58b099 100644 --- a/pybamm/__init__.py +++ b/pybamm/__init__.py @@ -109,6 +109,7 @@ def version(formatted=False): z_average, yz_average, boundary_value, + r_average ) from .expression_tree.functions import * from .expression_tree.parameter import Parameter, FunctionParameter diff --git a/pybamm/expression_tree/unary_operators.py b/pybamm/expression_tree/unary_operators.py index a24ef0c4df..15194dd432 100644 --- a/pybamm/expression_tree/unary_operators.py +++ b/pybamm/expression_tree/unary_operators.py @@ -1036,3 +1036,30 @@ def boundary_value(symbol, side): # Otherwise, calculate boundary value else: return BoundaryValue(symbol, side) + + +def r_average(symbol): + """convenience function for creating an average in the r-direction + + Parameters + ---------- + symbol : :class:`pybamm.Symbol` + The function to be averaged + + Returns + ------- + :class:`Symbol` + the new averaged symbol + """ + # If symbol doesn't have a particle domain, its r-averaged value is itself + if symbol.domain not in [["positive particle"], ["negative particle"]]: + new_symbol = symbol.new_copy() + new_symbol.parent = None + return new_symbol + # If symbol is a Broadcast, its average value is its child + elif isinstance(symbol, pybamm.Broadcast): + return symbol.orphans[0] + else: + r = pybamm.SpatialVariable("r", symbol.domain) + v = pybamm.Broadcast(pybamm.Scalar(1), symbol.domain) + return Integral(symbol, r) / Integral(v, r) diff --git a/pybamm/models/submodels/particle/base_particle.py b/pybamm/models/submodels/particle/base_particle.py index e454405b55..ab3e5c7f0f 100644 --- a/pybamm/models/submodels/particle/base_particle.py +++ b/pybamm/models/submodels/particle/base_particle.py @@ -26,12 +26,16 @@ def _get_standard_concentration_variables(self, c_s, c_s_xav): c_s_surf = pybamm.surf(c_s, set_domain=True) c_s_surf_av = pybamm.x_average(c_s_surf) + geo_param = pybamm.geometric_parameters if self.domain == "Negative": c_scale = self.param.c_n_max + active_volume = geo_param.a_n_dim * geo_param.R_n / 3 elif self.domain == "Positive": c_scale = self.param.c_p_max - + active_volume = geo_param.a_p_dim * geo_param.R_p / 3 + c_s_r_av = pybamm.r_average(c_s_xav) + c_s_r_av_vol = active_volume * c_s_r_av variables = { self.domain + " particle concentration": c_s, self.domain + " particle concentration [mol.m-3]": c_s * c_scale, @@ -48,6 +52,12 @@ def _get_standard_concentration_variables(self, c_s, c_s_xav): "X-averaged " + self.domain.lower() + " particle surface concentration [mol.m-3]": c_scale * c_s_surf_av, + self.domain + " electrode active volume fraction": active_volume, + self.domain + + " electrode volume-averaged concentration": c_s_r_av_vol, + self.domain + " electrode " + + "volume-averaged concentration [mol.m-3]": c_s_r_av_vol * c_scale, + self.domain + " electrode average extent of lithiation": c_s_r_av } return variables diff --git a/pybamm/models/submodels/particle/fickian/fickian_many_particles.py b/pybamm/models/submodels/particle/fickian/fickian_many_particles.py index 91ce9a325a..159a0a710f 100644 --- a/pybamm/models/submodels/particle/fickian/fickian_many_particles.py +++ b/pybamm/models/submodels/particle/fickian/fickian_many_particles.py @@ -30,6 +30,8 @@ def get_fundamental_variables(self): elif self.domain == "Positive": c_s = pybamm.standard_variables.c_s_p + # TODO: implement c_s_xav for Fickian many particles (tricky because this + # requires averaging a secondary domain) variables = self._get_standard_concentration_variables(c_s, c_s) return variables diff --git a/pybamm/solvers/base_solver.py b/pybamm/solvers/base_solver.py index 49ce4515f6..63dd82ad6d 100644 --- a/pybamm/solvers/base_solver.py +++ b/pybamm/solvers/base_solver.py @@ -73,7 +73,7 @@ def solve(self, model, t_eval): solution.total_time = timer.time() - start_time solution.set_up_time = set_up_time - pybamm.logger.info("Finish solving {} ({})".format(model.name, termination)) + pybamm.logger.warning("Finish solving {} ({})".format(model.name, termination)) pybamm.logger.info( "Set-up time: {}, Solve time: {}, Total time: {}".format( timer.format(solution.set_up_time), diff --git a/tests/unit/test_expression_tree/test_unary_operators.py b/tests/unit/test_expression_tree/test_unary_operators.py index 2e4f08cbf6..38649e5ee3 100644 --- a/tests/unit/test_expression_tree/test_unary_operators.py +++ b/tests/unit/test_expression_tree/test_unary_operators.py @@ -275,6 +275,31 @@ def test_average(self): with self.assertRaises(pybamm.DomainError): pybamm.x_average(a) + def test_r_average(self): + a = pybamm.Scalar(1) + average_a = pybamm.r_average(a) + self.assertEqual(average_a.id, a.id) + + average_broad_a = pybamm.r_average(pybamm.Broadcast(a, ["negative particle"])) + self.assertEqual(average_broad_a.evaluate(), np.array([1])) + + for domain in [ + ["negative particle"], + ["positive particle"] + ]: + a = pybamm.Symbol("a", domain=domain) + r = pybamm.SpatialVariable("r", domain) + av_a = pybamm.r_average(a) + self.assertIsInstance(av_a, pybamm.Division) + self.assertIsInstance(av_a.children[0], pybamm.Integral) + self.assertEqual(av_a.children[0].integration_variable[0].domain, r.domain) + # electrode domains go to current collector when averaged + self.assertEqual(av_a.domain, []) + + a = pybamm.Symbol("a", domain="bad domain") + with self.assertRaises(pybamm.DomainError): + pybamm.x_average(a) + def test_yz_average(self): a = pybamm.Scalar(1) z_average_a = pybamm.z_average(a)