With the global increase in energy demand and the urgency of carbon neutrality, lithium-ion energy storage cells have become a cornerstone of renewable energy systems. However, these energy storage cells generate significant heat during charge and discharge cycles, which can degrade performance and pose safety risks. Liquid cooling is widely adopted for its high efficiency and cost-effectiveness. In this study, we propose a novel stamped flow channel structure for liquid cooling of an energy storage cell pack. We establish a flow-thermal-solid coupled numerical model using CFD to analyze the internal flow field and overall temperature distribution. Our goal is to understand key factors affecting heat dissipation and flow resistance, and to optimize the channel geometry for improved temperature uniformity across the energy storage cells.
We focus on a battery pack consisting of 52 prismatic lithium iron phosphate (LFP) energy storage cells arranged in a 4×13 configuration. The pack sits on a liquid cold plate composed of a base plate, a channel plate (with a stamped serpentine channel), and inlet/outlet nozzles. Thermal interface material with a conductivity of 0.8 W/(m·K) is applied between the cells and the base plate to reduce contact resistance, while foam pads (0.022 W/(m·K)) separate adjacent cells to minimize lateral heat transfer. The channel dimensions and material properties are summarized in Table 1 and Table 2.
| Component | Dimensions (mm) |
|---|---|
| Base plate | Length 1145 × Width 800 × Thickness 1.2 |
| Channel plate | Length 1145 × Width 800 × Thickness 1.0 |
| Energy storage cell (LFP 314 Ah) | Width 173.8±0.8, Thickness 71.7±0.8, Height 207.2±0.8 |
| Nozzle inner diameter | 10 |
| Material | Density (kg/m³) | Viscosity (kg/(m·s)) | Specific heat (J/(kg·K)) | Thermal conductivity (W/(m·K)) |
|---|---|---|---|---|
| Coolant (50% ethylene glycol solution) | 1074.194 | 4.216e-3 | 3273 | 0.382 |
| Aluminum (base & channel plates) | 2719 | – | 871 | 202.4 |
| LFP energy storage cell (anisotropic) | 2284 | – | 1060 | kx=9.04, ky=11, kz=3.56 |
The flow channel, shown schematically in Figure (see inserted image), is 4.5 mm deep. It adopts a one-inlet to three-parallel-branch configuration, with two of those branches further splitting into sub-branches, and finally all branches merge into two main outlets before exiting. The entire path follows a serpentine pattern to enhance heat transfer area. We discretize the three-dimensional model using an unstructured polyhedral mesh with local refinement and prism layers (3 layers, first layer height 0.2 mm) at the channel wall to capture near-wall flow. The total mesh count after grid independence study is about 4 million cells (Table 3).
| Mesh ID | Max cell size (mm) | Min cell size (mm) | Number of cells | Outlet total pressure (Pa) | Deviation (%) | Avg top temperature of pack (°C) | Deviation (%) |
|---|---|---|---|---|---|---|---|
| #1 | 17 | 1.5 | 2,817,518 | 705.07 | – | 34.17 | – |
| #2 | 16 | 1.2 | 4,055,112 | 701.46 | 0.51 | 34.51 | 0.99 |
| #3 | 15 | 1 | 5,547,631 | 700.81 | 0.09 | 34.77 | 0.75 |
We employ the Reynolds-averaged Navier-Stokes (RANS) equations with the SST k-ω turbulence model. The governing equations for mass, momentum, and energy are:
$$ \frac{\partial \rho}{\partial t} + \frac{\partial (\rho u_i)}{\partial x_i} = 0 $$
$$ \frac{\partial (\rho u_i)}{\partial t} + \frac{\partial (\rho u_i u_j)}{\partial x_j} = -\frac{\partial p}{\partial x_i} + \frac{\partial}{\partial x_j} \left[ \mu \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) – \rho \overline{u_i’ u_j’} \right] $$
$$ \frac{\partial (\rho e)}{\partial t} + \frac{\partial (\rho u_i h)}{\partial x_i} = \frac{\partial}{\partial x_i} \left( k_{\text{eff}} \frac{\partial T}{\partial x_i} \right) $$
where $\rho$ is density, $p$ is pressure, $\mu$ is dynamic viscosity, $u_i$ are mean velocity components, $-\rho \overline{u_i’ u_j’}$ are Reynolds stresses, $e = h + u^2/2$ is total energy, and $k_{\text{eff}}$ is effective thermal conductivity. The SST k-ω model adds transport equations for turbulent kinetic energy $k$ and specific dissipation rate $\omega$:
$$ \frac{\partial (\rho k)}{\partial t} + \frac{\partial (\rho u_j k)}{\partial x_j} = \frac{\partial}{\partial x_j} \left( \Gamma_k \frac{\partial k}{\partial x_j} \right) + G_k – Y_k $$
$$ \frac{\partial (\rho \omega)}{\partial t} + \frac{\partial (\rho u_j \omega)}{\partial x_j} = \frac{\partial}{\partial x_j} \left( \Gamma_\omega \frac{\partial \omega}{\partial x_j} \right) + G_\omega – Y_\omega + D_\omega $$
We validate our CFD model by comparing simulated surface temperatures of a single LFP energy storage cell at 1C, 2C, and 3C discharge rates against published experimental data (Li, 2013). Excellent agreement is obtained with errors below 2%, confirming the model’s reliability. For the full pack simulation, we set a constant heat generation rate of 14 W per energy storage cell at 0.5C discharge, an initial temperature of 35°C, and a coolant inlet mass flow rate of 5 L/min at 18°C. The outlet is a pressure outlet at 0 Pa. All external surfaces are adiabatic, and radiation is neglected. The SIMPLEC algorithm and second-order upwind schemes are used.
Our results for the baseline case (5 L/min) show a pressure drop of 9.828 kPa across the channel. The flow velocity distribution is fairly uniform, with no significant recirculation zones. The coolant temperature rises from 18°C at the inlet to 20.49°C at the outlet (ΔT=2.49°C). The temperature distribution on the top faces of the energy storage cell pack is shown in Figure (inserted image). The highest temperature (35.17°C) appears in the fourth row of cells, while the lowest (33.77°C) appears in the first row, giving a maximum temperature difference of 1.40°C. This non-uniformity is because the first and third rows are directly above the three inlet branches where the coolant is coolest, whereas the fourth row sits above the warmest outlet branches.
We investigate the effect of inlet mass flow rate (5, 8, 10, 12, 15 L/min) on thermal performance and flow resistance. Table 4 summarizes the results.
| Flow rate (L/min) | Avg top temperature (°C) | Max temperature difference (°C) | Pressure drop (kPa) | Thermal resistance (°C/kW) |
|---|---|---|---|---|
| 5 | 34.00 | 1.40 | 9.828 | 23.58 |
| 8 | 33.10 | 1.22 | 20.215 | 22.87 |
| 10 | 32.74 | 1.30 | 30.112 | 22.45 |
| 12 | 32.50 | 1.35 | 42.347 | 22.12 |
| 15 | 32.26 | 1.39 | 59.441 | 21.79 |
As the flow rate increases, the average temperature decreases, but the rate of improvement slows: from 0.30°C/(L/min) initially to 0.08°C/(L/min) at the highest flow. Meanwhile, the pressure drop increases dramatically (from 9.8 to 59.4 kPa), while the thermal resistance decreases only modestly. The maximum temperature difference exhibits a minimum of 1.22°C at 8 L/min before rising again, indicating that there is an optimal flow rate for temperature uniformity. Clearly, simply increasing the flow rate is not an efficient way to improve cooling; a balance must be struck between pumping power and thermal performance.
To further improve the temperature uniformity, we modify the channel geometry by merging the four outlet branches into two branches and adding convex hull structures (dimples) in each outlet branch. These convex hulls are intended to enhance flow disturbance and thus convective heat transfer in the hottest region. We test three configurations: 3 pairs, 4 pairs, and 10 pairs of convex hulls per branch. The swirling strength method is used to quantify the vortex intensity in the flow. The swirling strength $\lambda_{ci}$ is derived from the decomposition of the velocity gradient tensor:
$$ \nabla \mathbf{u} = \text{symmetric part} + \text{antisymmetric part} $$
The imaginary part of the complex eigenvalue pair gives the local swirling strength. Figure (inserted image) shows the vortex structures for each configuration at $\lambda_{ci}=0.01$. Clearly, the presence of convex hulls generates strong secondary flows and vortices, breaking the thermal boundary layer and enhancing mixing.
We perform thermal simulations for each convex hull configuration at 5 L/min inlet flow. The maximum temperature difference across the energy storage cell pack top faces is shown in Table 5, along with the average swirling strength in the outlet branches.
| Configuration | Number of convex hull pairs | Average swirling strength (s⁻¹) | Max temperature difference on pack top (°C) |
|---|---|---|---|
| Original | 0 | 2.1 | 1.40 |
| Option 1 | 3 | 5.8 | 1.25 |
| Option 2 | 4 | 6.9 | 1.18 |
| Option 3 | 10 | 8.4 | 1.12 |
As the number of convex hull pairs increases, the average swirling strength increases and the maximum temperature difference decreases, confirming that enhanced disturbance improves temperature uniformity. However, the rate of improvement diminishes after 4 pairs, suggesting that further increases yield marginal benefits. Additionally, the pressure drop rises slightly with more convex hulls (from 9.83 kPa for the original to 10.21 kPa for 10 pairs), still acceptable. Therefore, from a manufacturing perspective, 4 pairs of convex hulls provide a good balance between performance and complexity.
In summary, we have developed a CFD model for a liquid-cooled energy storage cell pack with a novel stamped channel. Our key findings are:
- At 5 L/min, the flow resistance is 9.828 kPa and the maximum temperature difference across the energy storage cells is 1.40°C, showing good uniformity.
- Increasing the flow rate reduces average temperature but greatly increases pressure drop, and the improvement in temperature uniformity peaks at around 8 L/min.
- Adding convex hull structures in the outlet branches enhances flow disturbance and improves temperature uniformity, reducing the maximum temperature difference to 1.12°C with 10 pairs.
- A moderate number of convex hulls (4 pairs) is recommended for practical applications to achieve a substantial improvement without excessive complexity.
This study provides valuable insights into the thermal management design of lithium-ion energy storage cell packs, emphasizing the importance of channel geometry optimization and flow rate control to achieve both low thermal gradients and manageable pumping power. Future work will explore more complex geometries, such as ribbed channels or pin fins, and experimental validation of the optimized design.