Merge pull request pybamm-team#3439 from pybamm-team/rewrite-cooling-…

…terms simplify expression of cooling terms in x-lumped thermal models
agriyakhetarpal · Nov 7, 2023 · 33d2187 · 33d2187
2 parents d37a8f2 + ce8e56b
commit 33d2187
Show file tree

Hide file tree

Showing 4 changed files with 60 additions and 68 deletions.
diff --git a/docs/source/examples/notebooks/models/pouch-cell-model.ipynb b/docs/source/examples/notebooks/models/pouch-cell-model.ipynb
diff --git a/pybamm/models/submodels/thermal/lumped.py b/pybamm/models/submodels/thermal/lumped.py
@@ -56,10 +56,9 @@ def set_rhs(self, variables):
         # Newton cooling, accounting for surface area to volume ratio
         cell_surface_area = self.param.A_cooling
         cell_volume = self.param.V_cell
-        total_cooling_coefficient = (
-            -self.param.h_total * cell_surface_area / cell_volume
+        Q_cool_vol_av = (
+            -self.param.h_total * (T_vol_av - T_amb) * cell_surface_area / cell_volume
         )
-        Q_cool_vol_av = total_cooling_coefficient * (T_vol_av - T_amb)
 
         self.rhs = {
             T_vol_av: (Q_vol_av + Q_cool_vol_av) / self.param.rho_c_p_eff(T_vol_av)

diff --git a/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_1D_current_collectors.py b/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_1D_current_collectors.py
@@ -58,33 +58,29 @@ def set_rhs(self, variables):
         y = pybamm.standard_spatial_vars.y
         z = pybamm.standard_spatial_vars.z
 
-        # Account for surface area to volume ratio of pouch cell in surface and side
-        # cooling terms
-        cell_volume = self.param.L * self.param.L_y * self.param.L_z
-
+        # Calculate cooling, accounting for surface area to volume ratio of pouch cell
+        edge_area = self.param.L_z * self.param.L
         yz_surface_area = self.param.L_y * self.param.L_z
-        yz_surface_cooling_coefficient = (
+        cell_volume = self.param.L * self.param.L_y * self.param.L_z
+        Q_yz_surface = (
             -(self.param.n.h_cc(y, z) + self.param.p.h_cc(y, z))
+            * (T_av - T_amb)
             * yz_surface_area
             / cell_volume
         )
-
-        side_edge_area = self.param.L_z * self.param.L
-        side_edge_cooling_coefficient = (
+        Q_edge = (
             -(self.param.h_edge(0, z) + self.param.h_edge(self.param.L_y, z))
-            * side_edge_area
+            * (T_av - T_amb)
+            * edge_area
             / cell_volume
         )
-
-        total_cooling_coefficient = (
-            yz_surface_cooling_coefficient + side_edge_cooling_coefficient
-        )
+        Q_cool_total = Q_yz_surface + Q_edge
 
         self.rhs = {
             T_av: (
                 pybamm.div(self.param.lambda_eff(T_av) * pybamm.grad(T_av))
                 + Q_av
-                + total_cooling_coefficient * (T_av - T_amb)
+                + Q_cool_total
             )
             / self.param.rho_c_p_eff(T_av)
         }
@@ -94,7 +90,7 @@ def set_boundary_conditions(self, variables):
         T_amb = variables["Ambient temperature [K]"]
         T_av = variables["X-averaged cell temperature [K]"]
 
-        # find tab locations (top vs bottom)
+        # Find tab locations (top vs bottom)
         L_y = param.L_y
         L_z = param.L_z
         neg_tab_z = param.n.centre_z_tab
@@ -104,11 +100,10 @@ def set_boundary_conditions(self, variables):
         pos_tab_top_bool = pybamm.Equality(pos_tab_z, L_z)
         pos_tab_bottom_bool = pybamm.Equality(pos_tab_z, 0)
 
-        # calculate tab vs non-tab area on top and bottom
+        # Calculate tab vs non-tab area on top and bottom
         neg_tab_area = param.n.L_tab * param.n.L_cc
         pos_tab_area = param.p.L_tab * param.p.L_cc
         total_area = param.L * param.L_y
-
         non_tab_top_area = (
             total_area
             - neg_tab_area * neg_tab_top_bool
@@ -120,18 +115,22 @@ def set_boundary_conditions(self, variables):
             - pos_tab_area * pos_tab_bottom_bool
         )
 
-        # calculate effective cooling coefficients
+        # Calculate heat fluxes weighted by area
         # Note: can't do y-average of h_edge here since y isn't meshed. Evaluate at
         # midpoint.
-        top_cooling_coefficient = (
-            param.n.h_tab * neg_tab_area * neg_tab_top_bool
-            + param.p.h_tab * pos_tab_area * pos_tab_top_bool
-            + param.h_edge(L_y / 2, L_z) * non_tab_top_area
+        q_tab_n = -param.n.h_tab * (T_av - T_amb)
+        q_tab_p = -param.p.h_tab * (T_av - T_amb)
+        q_edge_top = -param.h_edge(L_y / 2, L_z) * (T_av - T_amb)
+        q_edge_bottom = -param.h_edge(L_y / 2, 0) * (T_av - T_amb)
+        q_top = (
+            q_tab_n * neg_tab_area * neg_tab_top_bool
+            + q_tab_p * pos_tab_area * pos_tab_top_bool
+            + q_edge_top * non_tab_top_area
         ) / total_area
-        bottom_cooling_coefficient = (
-            param.n.h_tab * neg_tab_area * neg_tab_bottom_bool
-            + param.p.h_tab * pos_tab_area * pos_tab_bottom_bool
-            + param.h_edge(L_y / 2, 0) * non_tab_bottom_area
+        q_bottom = (
+            q_tab_n * neg_tab_area * neg_tab_bottom_bool
+            + q_tab_p * pos_tab_area * pos_tab_bottom_bool
+            + q_edge_bottom * non_tab_bottom_area
         ) / total_area
 
         # just use left and right for clarity
@@ -141,21 +140,14 @@ def set_boundary_conditions(self, variables):
         self.boundary_conditions = {
             T_av: {
                 "left": (
-                    pybamm.boundary_value(
-                        bottom_cooling_coefficient * (T_av - T_amb),
-                        "left",
-                    )
-                    / pybamm.boundary_value(lambda_eff, "left"),
+                    pybamm.boundary_value(-q_bottom / lambda_eff, "left"),
                     "Neumann",
                 ),
                 "right": (
-                    pybamm.boundary_value(
-                        -top_cooling_coefficient * (T_av - T_amb), "right"
-                    )
-                    / pybamm.boundary_value(lambda_eff, "right"),
+                    pybamm.boundary_value(q_top / lambda_eff, "right"),
                     "Neumann",
                 ),
-            }
+            },
         }
 
     def set_initial_conditions(self, variables):

diff --git a/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_2D_current_collectors.py b/pybamm/models/submodels/thermal/pouch_cell/pouch_cell_2D_current_collectors.py
@@ -58,20 +58,22 @@ def set_rhs(self, variables):
         y = pybamm.standard_spatial_vars.y
         z = pybamm.standard_spatial_vars.z
 
+        # Calculate cooling
+        Q_yz_surface_W_per_m2 = -(self.param.n.h_cc(y, z) + self.param.p.h_cc(y, z)) * (
+            T_av - T_amb
+        )
+        Q_edge_W_per_m2 = -self.param.h_edge(y, z) * (T_av - T_amb)
+
         # Account for surface area to volume ratio of pouch cell in surface cooling
         # term
-        cell_volume = self.param.L * self.param.L_y * self.param.L_z
-
         yz_surface_area = self.param.L_y * self.param.L_z
-        yz_surface_cooling_coefficient = (
-            -(self.param.n.h_cc(y, z) + self.param.p.h_cc(y, z))
-            * yz_surface_area
-            / cell_volume
+        cell_volume = self.param.L * self.param.L_y * self.param.L_z
+        Q_yz_surface = pybamm.source(
+            Q_yz_surface_W_per_m2 * yz_surface_area / cell_volume, T_av
         )
-
         # Edge cooling appears as a boundary term, so no need to account for surface
         # area to volume ratio
-        edge_cooling_coefficient = -self.param.h_edge(y, z)
+        Q_edge = pybamm.source(Q_edge_W_per_m2, T_av, boundary=True)
 
         # Governing equations contain:
         #   - source term for y-z surface cooling
@@ -88,10 +90,8 @@ def set_rhs(self, variables):
             T_av: (
                 self.param.lambda_eff(T_av) * pybamm.laplacian(T_av)
                 + pybamm.source(Q_av, T_av)
-                + pybamm.source(yz_surface_cooling_coefficient * (T_av - T_amb), T_av)
-                + pybamm.source(
-                    edge_cooling_coefficient * (T_av - T_amb), T_av, boundary=True
-                )
+                + Q_yz_surface
+                + Q_edge
             )
             / self.param.rho_c_p_eff(T_av)
         }
@@ -102,24 +102,21 @@ def set_boundary_conditions(self, variables):
         y = pybamm.standard_spatial_vars.y
         z = pybamm.standard_spatial_vars.z
 
+        # Calculate heat fluxes
+        q_tab_n = -self.param.n.h_tab * (T_av - T_amb)
+        q_tab_p = -self.param.p.h_tab * (T_av - T_amb)
+        q_edge = -self.param.h_edge(y, z) * (T_av - T_amb)
+
         # Subtract the edge cooling from the tab portion so as to not double count
         # Note: tab cooling is also only applied on the current collector hence
-        # the (l_cn / l) and (l_cp / l) prefactors. We also still have edge cooling
+        # the (l_cn / l) and (l_cp / l) prefactors. We still have edge cooling
         # in the region: x in (0, 1)
-        h_tab_n_corrected = (self.param.n.L_cc / self.param.L) * (
-            self.param.n.h_tab - self.param.h_edge(y, z)
-        )
-        h_tab_p_corrected = (self.param.p.L_cc / self.param.L) * (
-            self.param.p.h_tab - self.param.h_edge(y, z)
-        )
-
-        negative_tab_bc = pybamm.boundary_value(
-            -h_tab_n_corrected * (T_av - T_amb) / self.param.n.lambda_cc(T_av),
+        negative_tab_bc = (self.param.n.L_cc / self.param.L) * pybamm.boundary_value(
+            (q_tab_n - q_edge) / self.param.n.lambda_cc(T_av),
             "negative tab",
         )
-        positive_tab_bc = pybamm.boundary_value(
-            -h_tab_p_corrected * (T_av - T_amb) / self.param.p.lambda_cc(T_av),
-            "positive tab",
+        positive_tab_bc = (self.param.p.L_cc / self.param.L) * pybamm.boundary_value(
+            (q_tab_p - q_edge) / self.param.p.lambda_cc(T_av), "positive tab"
         )
 
         self.boundary_conditions = {