Merge pull request #627 from pybamm-team/issue-623-soc_enhancement

Issue 623 soc enhancement

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')

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

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)

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

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

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),

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)