In the context of global energy transition, the integration of renewable sources like wind and solar power into the grid is paramount. However, these sources exhibit intermittent and fluctuating output, necessitating robust energy storage solutions to ensure grid stability, peak shaving, and frequency regulation. The battery energy storage system (BESS) has emerged as a key technology for this purpose, leveraging high-energy-density lithium-ion batteries. During charge and discharge cycles, these batteries generate significant heat due to irreversible and reversible electrochemical reactions. Inadequate thermal management can lead to elevated temperatures and substantial temperature gradients within the battery pack, accelerating degradation, reducing lifespan, and in extreme cases, triggering thermal runaway—a critical safety hazard. Studies indicate that a temperature difference of 10–15 K within a battery pack can degrade lifespan by 30–50%. Therefore, effective thermal management is essential to maintain the maximum temperature difference below 5 K, ensuring the longevity, safety, and economic viability of the battery energy storage system.
Among various cooling methods, liquid cooling is widely recognized for its high efficiency and compactness in managing thermal loads. This study focuses on optimizing the liquid cooling performance for a battery cluster within a large-scale battery energy storage system. We develop a comprehensive thermo-fluid simulation model based on actual dimensions to investigate the impact of outlet positioning and cold plate arrangement on temperature uniformity and heat dissipation. The goal is to enhance the thermal performance of the battery energy storage system, thereby supporting its reliable operation in stationary applications like grid-scale energy storage.
The battery cluster under investigation is a modular assembly typical of industrial battery energy storage systems. It comprises 25 battery modules, each containing 16 lithium iron phosphate (LiFePO4) cells with a nominal capacity of 120 Ah. The overall cluster dimensions are 1.52 m in length, 0.75 m in width, and 2.36 m in height. Each cell is mounted on a liquid cold plate measuring 545 mm × 367 mm × 8 mm. The cooling system features a manifold with inlet and outlet pipes; coolant flows from the inlet, branches into parallel streams feeding the cold plates, and converges at the outlet. Initially, the cold plates are positioned at the bottom of the cells. We explore design modifications, including relocating the outlet and reorienting the cold plates to the cell sides, to improve thermal uniformity. This optimization is crucial for the battery energy storage system to operate efficiently under high-load conditions.
To simulate the thermal behavior, we establish a detailed computational fluid dynamics (CFD) model. The battery cells are treated as homogeneous anisotropic solids with internal heat generation. Key thermal properties—density, specific heat capacity, and directional thermal conductivities—are derived from the constituent materials using mixture rules. For instance, the effective product of density and specific heat is calculated as:
$$ \rho c_p = \frac{\sum_i \rho_i c_{p,i} V_i}{\sum_i V_i} $$
where \( \rho \) is density, \( c_p \) is specific heat capacity, and \( V_i \) is the volume of each component (e.g., electrodes, separator, casing). The thermal conductivities in the longitudinal, transverse, and through-plane directions are computed using series-parallel resistance models. For components in parallel, the effective conductivity is:
$$ \lambda_1 = \frac{\sum_i L_i \lambda_i}{\sum_i L_i} $$
and for series arrangements:
$$ \lambda_2 = \frac{\sum_i L_i}{\sum_i L_i / \lambda_i} $$
Here, \( \lambda \) denotes thermal conductivity and \( L \) represents thickness. These calculations yield the anisotropic properties essential for accurate thermal modeling in the battery energy storage system.
The heat generation rate within the battery cells during discharge is modeled using the simplified Bernardi equation:
$$ q = q_i + q_r = \frac{1}{V} \left( I^2 R – I T \frac{dU_{OCV}}{dT} \right) $$
where \( q \) is the volumetric heat generation rate (W/m³), \( I \) is the current (A), \( R \) is the internal resistance (Ω), \( T \) is temperature (K), and \( dU_{OCV}/dT \) is the temperature-entropy coefficient. The internal resistance \( R \) and coefficient \( dU_{OCV}/dT \) vary with the state of charge (SOC). We performed Hybrid Pulse Power Characterization (HPPC) tests to obtain these parameters across SOC ranges, as summarized in the following table for typical discharge rates:
| SOC (%) | Internal Resistance, R (mΩ) | dUOCV/dT (mV/K) |
|---|---|---|
| 100 | 0.85 | -0.12 |
| 80 | 0.92 | -0.10 |
| 60 | 1.05 | -0.08 |
| 40 | 1.20 | -0.06 |
| 20 | 1.40 | -0.04 |
| 0 | 1.65 | -0.02 |
For simulation, we assume an average discharge rate of 0.5C, resulting in a constant heat generation rate of approximately 3 W per cell. The key physical parameters of the battery cell are consolidated below:
| Parameter | Value |
|---|---|
| Cell Model | LFP 120 Ah |
| Dimensions (mm) | 48 × 173 × 170 |
| Nominal Voltage | 3.2 V |
| Mass | 2.86 kg |
| Density, ρ | 2026 kg/m³ |
| Specific Heat, cp | 1020 J/(kg·K) |
| In-plane Thermal Conductivity | 2.0 W/(m·K) |
| Cross-plane Thermal Conductivity | 1.25 W/(m·K) |
| Top/Bottom Thermal Conductivity | 4.0 W/(m·K) |
The coolant is a water-glycol mixture with properties approximated as constant: density \( \rho_f = 1050 \, \text{kg/m}^3 \), specific heat \( c_{p,f} = 3500 \, \text{J/(kg·K)} \), thermal conductivity \( \lambda_f = 0.45 \, \text{W/(m·K)} \), and dynamic viscosity \( \mu = 0.0015 \, \text{Pa·s} \). The flow is turbulent, and we employ the RNG k-epsilon model to capture turbulence effects. The governing equations for fluid flow and heat transfer are solved using the finite volume method in ANSYS Fluent. The continuity, momentum, and energy equations for the fluid domain are:
$$ \frac{\partial \rho_f}{\partial t} + \nabla \cdot (\rho_f \vec{v}_f) = 0 $$
$$ \frac{\partial (\rho_f \vec{v}_f)}{\partial t} + \nabla \cdot (\rho_f \vec{v}_f \vec{v}_f) = -\nabla p + \nabla \cdot (\mu \nabla \vec{v}_f) + \rho_f \vec{g} $$
$$ \frac{\partial}{\partial t} (\rho c_p T) + \nabla \cdot (\rho c_p \vec{v} T) = \nabla \cdot (\lambda \nabla T) $$
For the solid battery domain, the energy equation includes the heat source term:
$$ \frac{\partial}{\partial t} (\rho_b c_{p,b} T) = \nabla \cdot (\lambda_b \nabla T) + q $$
And for the aluminum cold plates and pipes:
$$ \frac{\partial}{\partial t} (\rho_{al} c_{p,al} T) = \nabla \cdot (\lambda_{al} \nabla T) $$
The turbulence kinetic energy \( k \) and dissipation rate \( \epsilon \) are modeled as:
$$ \frac{\partial}{\partial t} (\rho_f k) + \nabla (\rho_f k \vec{v}_f) = \nabla \left( \frac{\mu}{\sigma_k} \nabla k \right) + G_k – \rho_f \epsilon $$
$$ \frac{\partial}{\partial t} (\rho_f \epsilon) + \nabla (\rho_f \epsilon \vec{v}_f) = \nabla \left( \frac{\mu}{\sigma_\epsilon} \nabla \epsilon \right) + \frac{\epsilon}{k} (C_{1\epsilon} G_k – C_{2\epsilon} \rho_f \epsilon) $$
where \( G_k = \mu (\nabla \vec{v}_f + (\nabla \vec{v}_f)^T) \cdot \nabla \vec{v}_f \), and constants are \( C_{1\epsilon} = 1.42 \), \( C_{2\epsilon} = 1.68 \), \( \sigma_k = 0.7194 \), \( \sigma_\epsilon = 0.7194 \). The SIMPLE algorithm is used for pressure-velocity coupling, with second-order upwind schemes for momentum, energy, and turbulence quantities. Boundary conditions include a mass flow inlet of 9 L/min (0.15 kg/s) at 293 K, pressure outlet at atmospheric pressure, and initial/environmental temperature of 303 K. A polyhedral mesh is generated, and grid independence is verified; a mesh with 4.8 million cells yields results within 1.2% deviation of finer meshes, balancing accuracy and computational cost for this battery energy storage system model.
We first investigate the effect of outlet location on coolant distribution and thermal performance. Four configurations are analyzed, as illustrated schematically (though not shown here per guidelines). In all cases, the inlet is fixed on the left side of the supply manifold. The outlets are positioned: (a) bottom-left, (b) bottom-right, (c) top-left, and (d) top-right. The goal is to achieve balanced flow distribution among the parallel cold plates. The temperature distribution and maximum temperature difference across the battery cluster are evaluated after a 4000-second discharge simulation. The results are summarized in the table below:
| Configuration | Outlet Position | Max Battery ΔT (K) | Flow Distribution Uniformity |
|---|---|---|---|
| A | Bottom-left | 16.0 | Poor |
| B | Bottom-right | 15.5 | Poor |
| C | Top-left | 13.2 | Moderate |
| D | Top-right | 13.0 | Good |
Configuration D, with the outlet on the top-right opposite the inlet, yields the most uniform flow because the flow paths to each cold plate have nearly equal lengths—this is a principle of hydraulic balancing. The maximum temperature difference reduces from 16 K to 13 K, an improvement of 18.8%. The temperature non-uniformity primarily arises from uneven coolant supply; modules closer to the outlet experience lower flow rates in configurations A and B due to longer path resistance. To quantify flow distribution, we define a uniformity index \( \eta \) based on the standard deviation of flow rates among cold plates:
$$ \eta = 1 – \frac{\sigma_{\dot{m}}}{\bar{\dot{m}}} $$
where \( \sigma_{\dot{m}} \) is the standard deviation and \( \bar{\dot{m}} \) is the mean flow rate. The indices for the four configurations are:
| Configuration | Uniformity Index, η |
|---|---|
| A | 0.65 |
| B | 0.68 |
| C | 0.78 |
| D | 0.92 |
Higher \( \eta \) indicates better balance. Configuration D’s design minimizes flow maldistribution, which is critical for maintaining temperature homogeneity in a large battery energy storage system. The evolution of the average temperature difference between the hottest and coldest modules over time further demonstrates this; for configuration D, the difference stabilizes around 4.6 K after 2000 seconds, whereas for A and B, it grows steadily to over 10 K. This underscores the importance of hydraulic design in liquid-cooled battery energy storage systems.
Next, we explore the impact of cold plate orientation. The baseline design places cold plates beneath the cells (bottom cooling). We propose a novel arrangement where cold plates are positioned vertically between the cells, contacting the larger side surfaces (side cooling). This increases the contact area from approximately 0.083 m² per cell (bottom) to 0.294 m² per cell (side), and reduces the conduction thickness through the cell to the cooled surface. The thermal resistance network can be modeled as a series of resistances. For bottom cooling, the dominant resistance is through the cell thickness in the vertical direction. For side cooling, the resistance is through the cell width, which is smaller due to larger area and shorter path. The overall thermal resistance \( R_{th} \) from the heat source to coolant is approximated as:
$$ R_{th} = \frac{L}{\lambda A} $$
where \( L \) is conduction length, \( \lambda \) is thermal conductivity, and \( A \) is contact area. For side cooling, \( L \) is reduced and \( A \) increased, significantly lowering \( R_{th} \). We simulate both arrangements using the optimized outlet configuration D. The results after 4000 seconds of discharge are compelling:
| Cold Plate Arrangement | Max Battery ΔT (K) | Average Cell Temperature (K) | Cooling Efficiency Gain |
|---|---|---|---|
| Bottom Cooling | 13.0 | 308.5 | Baseline |
| Side Cooling | 3.7 | 306.2 | 71.5% ΔT reduction |
The side cooling configuration reduces the maximum temperature difference from 13 K to 3.7 K, a 71.5% improvement, well within the 5 K requirement for battery energy storage systems. The temperature distribution becomes remarkably uniform. The cooling efficiency \( \epsilon_c \) can be defined as:
$$ \epsilon_c = \frac{q_{dissipated}}{q_{generated}} $$
where \( q_{dissipated} \) is the heat removed by the coolant and \( q_{generated} \) is the total heat generated. For side cooling, \( \epsilon_c \) approaches 0.98, compared to 0.85 for bottom cooling, indicating superior heat extraction. This is attributed to enhanced convective heat transfer due to larger wetted area and reduced conductive resistance. The temperature evolution over time shows that for side cooling, the maximum ΔT peaks early and then declines, stabilizing below 4 K, whereas for bottom cooling, ΔT monotonically increases. This behavior aligns with the thermal time constant \( \tau \), given by:
$$ \tau = \frac{\rho c_p V}{h A} $$
where \( h \) is the convective heat transfer coefficient, \( V \) is cell volume, and \( A \) is cooled area. For side cooling, \( A \) is larger, reducing \( \tau \) and allowing faster thermal equilibrium. Additionally, the pressure drop across the cooling system is analyzed. Side cooling may increase flow resistance due to more complex channels, but our simulation shows only a marginal rise from 12 kPa to 15 kPa, which is manageable for pump selection in a battery energy storage system.
To generalize the findings, we derive dimensionless correlations for Nusselt number \( Nu \) and friction factor \( f \) for the cold plate designs. For turbulent flow in rectangular channels, the Dittus-Boelter equation is adapted:
$$ Nu = 0.023 Re^{0.8} Pr^{0.4} $$
where \( Re \) is Reynolds number and \( Pr \) is Prandtl number. The experimental correlation for our specific channel geometry yields:
$$ Nu = 0.021 Re^{0.8} Pr^{0.33} \left( \frac{A_{contact}}{A_{total}} \right)^{0.15} $$
The friction factor is given by:
$$ f = 0.184 Re^{-0.2} + \Delta f_{geometry} $$
These correlations help in scaling the design for other battery energy storage system configurations. Moreover, we evaluate the energy efficiency of the cooling system itself. The pump power \( P_{pump} \) is:
$$ P_{pump} = \frac{\dot{m} \Delta p}{\rho \eta_{pump}} $$
where \( \dot{m} \) is mass flow rate, \( \Delta p \) is pressure drop, and \( \eta_{pump} \) is pump efficiency. For our system, \( P_{pump} \) is about 5 W, negligible compared to the megawatt-scale energy throughput of the battery energy storage system, affirming the practicality of liquid cooling.
The optimization also considers transient operational scenarios, such as varying discharge rates and ambient temperatures. We simulate a duty cycle with 1C discharge for 1800 seconds, followed by 0.25C charge for 3600 seconds. The side cooling configuration maintains ΔT below 5 K throughout, whereas bottom cooling exceeds 10 K during high discharge. This robustness is vital for battery energy storage systems facing dynamic grid demands. Furthermore, we examine the impact of coolant temperature. Lowering inlet temperature from 293 K to 283 K reduces average temperature but can increase temperature gradients if flow distribution is poor. With side cooling, the gradient remains low, demonstrating its effectiveness across operating conditions.
In conclusion, this study demonstrates significant thermal performance improvements in a battery energy storage system through liquid cooling optimization. By redesigning the outlet location to achieve hydraulic balance, the maximum temperature difference is reduced by 18.8%. More profoundly, reorienting cold plates from the bottom to the sides of cells enhances contact area and reduces thermal resistance, slashing the maximum temperature difference by 71.5% to 3.7 K, within the recommended 5 K threshold. These modifications ensure better temperature uniformity, which prolongs battery life, enhances safety, and boosts the overall efficiency of the battery energy storage system. The methodologies and correlations presented here provide a framework for designing advanced thermal management systems for large-scale battery energy storage systems, supporting the integration of renewable energy and grid stability. Future work could explore hybrid cooling techniques, advanced coolant formulations, and real-time adaptive control to further optimize the battery energy storage system performance under diverse environmental and load conditions.