Thermal management of cylindrical lithium-ion batteries using internal and external phase change material placement

A

= area [m²];

Bi

= Biot number [/];

C

= capacity [Ah];

c

= specific heat [$J/(kg·K)$];

c_e

= concentration in the electrolyte phase [mol/m³];

c_s

= concentration in the solid phase [mol/m³];

D

= diffusion coefficient [m²/s];

d

= thickness [m];

E

= energy [J/mol];

F

= Faraday’s constant [C/mol];

$f±$

= mean molar activity coefficient [/];

g

= gravity acceleration [m/s²];

H_l

= latent heat of fusion [J/kg];

h

= convective heat transfer coefficient [$J/(m2·K)$];

I

= improvement [%];

i

= Interfacial current density [A/m²];

i₀

= Exchange current density [A/m²];

k

= thermal conductivity [$W/(m·K)$];

k₀

= kinetic rate constant [m/s];

L

= characteristic length [m];

l

= height [m];

Nu

= Nusselt number [/];

Pr

= Prandtl number [/];

$Q˙$

= heat generation [W/m³];

R

= universal gas constant [$J/(mol·K)$];

R_s

= solid phase particle radius [m];

r

= radial coordinate [m];

Ra

= Raileigh number [/];

S_h

= latent heat term [W/m³];

T

= temperature [K];

t

= time [s];

$tLi+$

= electrolyte ionic current fraction by lithium ions [/];

U

= voltage [V];

V

= volume [m³];

x

= axial coordinate [m]; and

z

= axial coordinate [m].

α

= thermal diffusivity [m²/s];

α_a

= transfer coefficient [/];

β

= thermal expansion coefficient [1/K];

ε

= volume fraction [/];

η

= overpotential [V];

θ

= angular coordinate [rad];

ν

= kinematic viscosity [m²/s];

ρ

= density [kg/m³];

σ

= electrical conductivity [S/m];

φ

= electric potential [V];

χ

= correction factor [/]; and

ψ

= liquid phase fraction [/].

a

= anodic

act

= activation;

amb

= ambient;

ap

= applied;

av

= average;

c

= cathodic

cc

= current collector;

co

= core;

n

= negative;

e

= electrolyte;

ele

= electrode;

eq

= equilibrium;

f

= film;

gen

= generation;

jou

= Joule;

l

= liquid;

lim

= limit;

man

= mandrel;

max

= maximum;

n

= negative;

nom

= nominal;

op

= operational;

opt

= optimal;

p

= positive;

ref

= reference;

RL

= representative layer;

s

= solid phase;

sav

= saving;

sep

= separator;

sol

= solidification; and

sur

= surface.

eff

= effective.

BMS

= battery management system;

LCO

= lithium cobalt oxide;

LIB

= lithium ion battery;

P2D

= pseudo-two-dimensional;

PCM

= phase change material;

RMS

= root mean square;

SoC

= state of charge;

SPM

= single particle model; and

TR

= thermal runaway.

The rapid advancement of electrification technologies has placed lithium-ion batteries (LIBs) at the center of modern energy solutions, particularly in the fields of electric mobility (Veza et al., 2024) and renewable energy storage (Yang et al., 2018). Despite their high energy density, efficiency and cycle life, LIBs face critical challenges related to thermal management. Excessive heat generation during high charge/discharge cycles can lead to capacity fade (Atalay et al., 2020), performance loss and, in extreme cases, thermal runaway (TR) (Feng et al., 2019; Mallick and Gayen, 2023), posing significant safety risks.

To prevent thermal problems due to battery overheating, different multi-physical models are proposed to predict heat generation within the battery, so that the appropriate thermal management system for the battery pack can be designed. Models must be multi-physical, as their thermal generation depends on electrical or chemical parameters (Jordan et al., 2024). The models can be multi-dimensional, to correctly describe the behavior of the cell, or compact and fast, that can be integrated into battery management systems (BMSs) (Li et al., 2020). Multi-dimensional battery models are primarily classified as electro-thermal or electrochemical–thermal. Electro-thermal models estimate heat generation using equivalent electrical circuit parameters, such as Thevenin-based models, which represent battery behavior through an equivalent electrical circuit (Nikolian et al., 2016). In contrast, electrochemical–thermal models use governing equations to analyze the lithium-ion concentration within the battery, enabling a multi-scale analysis: from the microscale of electrode solid particles to the macroscale of the entire cell. Both types of models share a common approach to heat generation calculation, which is based on the formulation proposed by (Bernardi et al., 1985). This method accounts for a reversible heat component, which arises due to entropy variations within the system. This entropy-related heat contribution depends on the state of charge (SoC) and the entropy coefficient, and it influences the overall thermal behavior of the battery during charge and discharge cycles. Notably, the entropy coefficient can assume positive or negative values depending on the SoC and whether the battery is charging or discharging. Consequently, the sign of the reversible heat term also changes, leading to either a reduction or an increase in the total heat generation. One of the most widely adopted electrochemical–thermal models is the pseudo-two-dimensional (P2D) model (Doyle et al., 1993; Fuller et al., 1994), which provides a detailed representation of electrode kinetics and transport phenomena, or the single particle model (SPM) (Zhang et al., 2000) which is a model simplification, assuming that the concentration of lithium ions in the liquid phase is uniform throughout the battery, and the electric potential of the solid phase is uniform across the electrode. Different P2D model simplifications are proposed to improve simulation time (Jokar et al., 2016). For compact and fast models, the approach is primarily based on electrical parameters and are an optimal solution for temperature prediction in a real-time analysis. In this case, the lumped model (Catalano et al., 2023) can be used to reduce computational time. Also, innovative batteries such as sodium-sulfur batteries are studied through a lumped thermal model (Csemány, 2025).

To mitigate temperature issues, effective thermal management strategies are essential to maintain optimal operating temperatures, avoid the hotspot that can lead to thermal runaway (Afzal et al., 2022) and enhance both battery lifespan and efficiency (Lin et al., 2021). Various approaches have been proposed, including air (Chen et al., 2024a), and liquid cooling systems (Wu et al., 2024), heat pipes (Su et al., 2025) and advanced thermal interface materials (Lee et al., 2023). An emerging approach are dielectric fluids thermal management, Menale et al. (2019) analyzed two different dielectric liquids, Clearco-50cSt and Galden HT135 and compared them to air, proving that Galden HT135 is an efficient liquid for thermal management, while Williams et al. (2024) studied the phase change process under pool boiling conditions for a dielectric immersion cooling with Novec 7000. Also, passive and hybrid thermal management solutions (Zhao et al., 2020) are gaining increasing attention for their ability to balance thermal performance with space and energy efficiency, making them particularly attractive for electric vehicle and stationary storage applications. Finally, one of the most studied passive solutions is with phase change materials (PCM) (Pilali et al., 2025), leveraging their ability to absorb and store heat through latent heat during phase transition. Unlike conventional cooling methods, phase change materials offer intrinsic temperature stabilization during their solid-to-liquid transition. As the battery generates heat during operation, this thermal energy is absorbed by the PCM in the form of latent heat. While absorbing heat, the PCM remains at its nearly-constant phase change temperature, which limits the heat accumulation in the battery itself and reduces its temperature rise. Recent studies have investigated the enhancement of PCM thermal conductivity via the addition of nanoparticles, giving rise to nano-enhanced PCM. For instance, a numerical–experimental study on 18,650 Li-ion cells was conducted using paraffin wax combined with different nanoparticles such as CuO, Al₂O₃ and TiO₂ at varying concentrations to improve heat dissipation (Vyas et al., 2024). Also with this application, different hybrid solutions are proposed (Khan et al., 2023), for example with fins (Zhang et al., 2023), metal foam integration (Buonomo et al., 2018) or honeycomb structures (Sutheesh et al., 2024; Chen et al., 2024b). For instance, a hybrid BTMS combining PCM with a liquid-cooled plate arranged in a honeycomb configuration has been proposed for cylindrical cells, showing improved temperature uniformity and dissipation effectiveness through 3D simulations (Zhuang et al., 2021). However, these approaches typically require external structural elements, increasing the system’s size and complexity. In contrast, this study introduces a novel passive cooling configuration by integrating PCM internally inside the battery’s hollow mandrel, without increasing the external footprint. This internal PCM can be combined with a standard external PCM layer to create a hybrid layout that maximizes battery thermal dissipation. This is an opportunity to investigate such configurations using a coupled electrochemical–thermal model. This approach is of particular interest for compact systems, such as portable electronics and electric vehicles, where space is a critical constraint.

To address the associated challenges, a multiscale coupled electrochemical–thermal model of a lithium cobalt oxide (LCO) 18650 cylindrical battery is developed to investigate the thermal behavior under different passive thermal management configurations. The analysis focuses on three PCM arrangements: internal (within the mandrel), external (surrounding the cell) and a combined internal/external layout, assessing their effectiveness in containing temperature rise while minimizing the impact on the battery pack volume. The PCM analyzed is an inorganic compound, the sodium thiosulphate pentahydrate, the electrochemical model is based on the pseudo-two-dimensional approach, while the thermal model is treated as two-dimensional, leveraging symmetry along the angular direction. These strategies offer significant advantages in heat dissipation compared to a battery without PCM, with the internal PCM approach also providing space optimization benefits, allowing for the development of more compact battery packs. In addition, a parametric study on PCM external thickness is conducted to determine the optimal PCM layer required to maintain the battery’s operating temperature below two critical thresholds: a battery temperature limit of 50°C (Rodrigues et al., 2017; Jiang et al., 2020), as recommended by manufacturers to prevent accelerated aging and thermal runaway, and an operational temperature limit of 45°C, set to provide a safety margin and avoid approaching the maximum threshold too closely. This conservative limit also accounts for the progressive aging of the cell, which may lead to reduced thermal performance over time and an increased electrical resistance – an aspect not foreseen by the thermal management system – thus justifying the degree of thermal oversizing to ensure long-term reliability. This analysis is performed for both the external-only and internal/external PCM configurations, highlighting the potential space-saving benefits of a combined PCM approach.

The main research gaps and objectives of this study are summarized as follows:

• Gap 1: Most existing works on PCM-enhanced thermal management focus on external-only or hybrid solutions (e.g. fins, metal foams), which increase the overall battery footprint.

• Gap 2: The spatial efficiency of internal PCM configurations in cylindrical cells remains underexplored in comparison to more conventional external PCM approaches.

• Objective 1: To develop a fully coupled electrochemical–thermal model that captures both cell dynamics and PCM phase change effects with internal/external placement.

• Objective 2: To compare the thermal performance of internal, external and combined PCM configurations under identical volume and thermal boundary conditions.

• Objective 3: To assess the potential space savings in percentage terms and thermal benefits of internal PCM via a parametric analysis on PCM thickness.

This approach aims to support the design of compact and thermally stable battery systems, particularly for high energy density applications such as electric mobility and portable electronics.

In this section, the coupled electrochemical–thermal model of the battery and the thermal model for the phase change material are presented. The electrochemical–thermal model integrates a pseudo-two-dimensional electrochemical framework to simulate the internal electrochemical reactions, with a two-dimensional axisymmetric thermal model to predict heat generation $Q˙$ and evaluate temperature distribution $T$ within the battery. The interaction between these models enables a comprehensive analysis of the battery’s thermal response under various operating conditions. The coupling method between the models is also analyzed, showing how certain electrochemical parameters, such as the electrolyte diffusion coefficient $De$⁠, the electrolyte conductivity $σe$⁠, the mean molar activity coefficient $f±$ and the potential equilibrium $Ueq$⁠, are temperature-dependent and affect the calculation of heat generation $Q˙$ within the battery core, leading to a non-linearity of the problem. Meanwhile, the thermal behavior of the phase change material is described by an enthalpy-based model that accounts for latent heat absorption and release during phase transitions. A commercial finite element analysis solver, COMSOL Multiphysics, was used for this analysis.

A lithium-ion battery is composed of multiple repeating layers of current collectors, electrodes and separators, which are tightly wound into a cylindrical structure known as a “jelly roll”. The layers are typically wrapped around a support structure called mandrel, to provide mechanical stability. In this design, both the positive and negative current collectors are coated of their respective electrode, forming a compact, multi-layered configuration that enhances energy density. During charge and discharge, lithium ions migrate between the two electrodes. The model uses a pseudo-two-dimensional (P2D) approach (Doyle et al., 1993; Fuller et al., 1994), incorporating two interconnected spatial dimensions, as illustrated in Figure 1 during a discharge. The first dimension, $x$⁠, represents the cell thickness and captures charge and mass transport in the electrodes and electrolyte, along with the electrochemical reactions at the electrode-electrolyte interface. The second dimension corresponds to the particle radius $r$ and is used to solve the lithium diffusion equations within the electrode’s active material in spherical coordinates. This coupled approach enables a detailed representation of both macroscopic and microscopic transport phenomena within the battery. To enhance computational efficiency, the following assumptions are introduced:

• Current distribution is considered one-dimensional along the x-direction.

• The current collectors are assumed to have sufficiently high electrical conductivity $σ$⁠, ensuring a homogeneous potential distribution along the x-direction.

• The active material particles in the electrodes are considered to have a uniform size distribution with the same particle radius $Rs$⁠.

• Side reactions and gas generation during electrochemical processes are neglected.

Due to the second assumption, the electrochemical model does not include the current collectors, and no electrochemical reactions are present within them. As a result, the model does not account for the ohmic drop along with the current collectors. However, these components will be incorporated into the thermal model. The electrochemical computational domain includes only a representative layer, consisting of a single layer of each component. This simplification is justified by the assumption that all layers exhibit identical behavior. In the figure below, $dele,p$⁠, $dsep$ and $dele,n$ represent the positive electrode, separator and negative electrode thickness, while $dRL=dele,p+dsep+dele,n$ is the total thickness of the representative layer.

During the discharging process, lithium ions migrate from the negative electrode to the positive electrode, passing through the electrolyte. This process reversed during charge. The concentration of lithium ions in the solid phase $cs$⁠, within both electrodes is governed by the mass conservation equation. It is expressed as a function of the radial coordinate $r$ within a spherical particle and time $t$⁠:

where $Ds$ is the diffusion coefficient within the solid phase. As a boundary condition for the spherical particle, two Neumann conditions are used at the particle center [equation (2)] and at the particle surface [equation (3)]. In the first one, a symmetry condition is imposed at the particle center, ensuring that there is no net lithium flux at $r = 0$⁠. Since the diffusion process is radially symmetric, there is no preferential direction for lithium movement at the core, leading to a zero-concentration gradient. The second boundary condition links the lithium flux at the particle surface to the electrochemical reaction occurring at the electrode–electrolyte interface, ensuring that lithium flux at the surface matches the rate at which lithium is inserted into or extracted from the particle due to the applied current:

where $i$ represents the interfacial current density and $F$ is the Faraday’s constant.

As initial condition, a homogeneous concentration in the solid particle $cs,0$ is imposed. In Figure 2, a schematic representation of the boundary conditions for the lithium-ion concentration within the solid phase is shown.

The concentration of lithium ions in the electrolyte $ce$⁠, varies as a function of both cell thickness $x$ and time $t$⁠. It is described through equation (4) that represents the conservation of lithium ions in the electrolyte phase:

It consists of the following key terms:

• Time-dependent term$εe∂ce∂t$⁠: representing the accumulation or depletion of lithium ions in the electrolyte over time, where $εe$ is the electrolyte volume fraction, which varies across the battery components, accounting that only a portion of the battery’s volume is occupied by the electrolyte.

• Diffusive term$∂∂x(Deeff∂ce∂x)$⁠: describing how lithium ions diffuse through the electrolyte, driven by concentration gradients. where $Deeff$ is the effective diffusion coefficient that accounts for the presence of a porous structure composed of the electrodes and the separator, which affects ion movement. For this reason, the effective diffusion coefficient $Deeff$ is given by the following equation:

where$De$ is the intrinsic diffusion coefficient of lithium ions in the electrolyte and $p$ is the Bruggeman porosity exponent, assumed to be 1.5 for both electrodes and the separator. This empirical factor accounts for the tortuous pathways that lithium ions follow in a porous medium, reducing the effective diffusivity compared to a bulk electrolyte. This formulation ensures that the transport properties of lithium ions are correctly adjusted based on the microstructure of the battery materials:

• Source term$1-tLi+AsiF$⁠: representing the effect of electrochemical reactions at the electrode-electrolyte interface, where $tLi+$ is defined as the fraction of the total ionic current in the electrolyte that is carried by lithium ions, assuming a uniform electrolyte composition and $As$ is the specific interfacial area, defined in equation (6), that scales the reaction rate by representing the available active surface for lithium-ion exchange at the electrode-electrolyte interface:

where $εs$ represents the volume fraction occupied by the solid phase within the porous electrode.

At both ends of the electrolyte domain, $x = 0$ and $x=dRL$⁠, there is a Neumann’s boundary condition assuming no flux of lithium ions:

These boundary conditions indicate that lithium ions are not lost or gained at the boundaries, ensuring conservation of mass. Physically, this implies that the electrodes act as barriers, preventing lithium-ion diffusion beyond the defined electrolyte region. In the region between the electrodes and the separator, Li-ion concentration is assumed continuous across the surfaces. As an initial condition, a homogeneous concentration in the electrolyte $ce,0$ is settled.

The potential within the solid phase $φs$ obeys the charge conservation equation through Ohm’s law, which can be expressed as in equation (8):

where $σseff$ represents the effective electrical conductivity of the solid phase, as it depends on the porous structure and it is corrected through the Bruggeman coefficient $p$ equal to 1.5 and $σseff=σsεsp$⁠, where $σs$ is the intrinsic electrical conductivity of the active material. The boundary conditions governing this equation are:

• At the interface between the positive electrode and the positive current collector, $x=0$⁠, all the applied current is transported through the solid particles that coat the current collector:

where $iap$ is the applied current density per unit area of the positive electrode plate:

• At the interface between the negative electrode and the negative current collector, $x=dRL$⁠, a Dirichlet boundary condition with a zero-potential for a grounded reference condition:

• At the interfaces between the electrodes and the separator (i.e. $x=dele,p$⁠, $x=dRL-dele,n$⁠), charge transport occurs exclusively through the electrolyte, leading to a Neumann boundary condition, with a zero-flux condition at the electrode–separator interfaces. This occurs because, at these locations, the charge is no longer transported through the solid phase but rather through the electrolyte:

As initial conditions, the solid-phase potential $φs$ is set equal to the equilibrium potential $Ueq$ throughout the positive electrode thickness $dele,p$⁠, while a reference solid-phase zero-potential $φs$ is imposed across the negative electrode thickness $dele,n$⁠.

The potential in the electrolyte phase $φe$ is governed by the following equation:

This equation consists of three main contributions:

1. The conductive term$∂∂xσeeff∂φe∂x$⁠: Describing ionic charge transport by conduction within the electrolyte, driven by the electric potential gradient $∂φe∂x$⁠. The effective ionic conductivity of the electrolyte $σeeff$ is affected from porosity and tortuosity effects, using the Bruggeman correlation $σeeff=σeεep$ with the Bruggeman coefficient $p$ equal to 1.5.

2. Diffusive and electro migrative term$∂∂x2RTσeeffF1-tLi+1+dln f±dln ce∂ln ce∂x$⁠: This term accounts for the effects of ion diffusion (due to concentration gradients) and electromigration (due to Coulombic force acting on ions in an electric field) in the electrolyte. $2RTσeeff/F$ represents the thermodynamic contribution to diffusion, dependent on temperature $T$⁠, where $R$ is the universal gas constant. $f±$ is the mean molar activity coefficient which corrects for deviations from ideal behavior in electrolyte solutions due to ion interactions.

3. Electrochemical term$Asi$⁠: due to the charge transfer reaction at the electrode-electrolyte interface.

The same boundary condition is applied at the electrode–current collector interfaces (i.e. $x = 0$⁠, $x=dRL$⁠) assuming that all current $i$ is transported through the solid phase. This leads to a of Neumann condition for the electrolyte, as shown in equation (13). Potential continuity is assumed across the electrodes and the separator. In summary, in Figure 3 all the boundary conditions are represented for the lithium-ion concentration within the electrolyte, and for the solid and electrolyte potential, without representing continuity conditions:

As initial conditions, the electrolyte-phase potential $φe$ is set to the equilibrium potential $Ueq$ across the positive electrode thickness $dele,p$ and a reference electrolyte-phase zero-potential $φe$ is assigned throughout the negative electrode thickness $dele,n$⁠. In the separator region, the electrolyte-phase potential $φe$ is assumed to transition linearly between the values in the adjacent electrodes.

The electrochemical reaction kinetics at the interface between the solid electrode and the electrolyte follows the Butler–Volmer equation, describing the current density $i$ as in the following equation:

where $i0$ represents the exchange current density, defined in equation (15), $αa$ and $αc$ denote the anodic and cathodic transfer coefficients, respectively, representing the directionality and balance of the overall reaction and $η$ is the overpotential, which quantifies the deviation between the actual potential $φs-φe$ and the equilibrium voltage $Ueq$⁠, defined in equation (16):

where $k0$ is the kinetic rate constant, $ce,ref$ is the reference electrolyte concentration, equal to the initial electrolyte concentration $ce,0$$cs,max$ is the maximum concentration in the solid phase and $cs,R$ indicates the solid phase concentration at the surface of the active particle, with $r=Rs$⁠. The equilibrium voltage $Ueq$ is a temperature-dependent value, approximated using a first-order Taylor series expansion:

Here, $Ueq,ref$ represents the reference equilibrium voltage, while $Tref$ denotes the corresponding reference temperature $∂Ueq∂T$⁠.

The computational domain of the battery thermal model is composed of a mandrel, a core and a case. The thermal behavior of the battery is modeled using cylindrical coordinates $(r,θ,z)$⁠. Due to the inherent symmetry of the cell, of the boundary conditions and of the thermal conductivity $k$ along the angular direction, the coordinate $θ$ is neglected, allowing the model to be simplified to a two-dimensional axisymmetric framework. This reduction, which considers only radial $r$ and axial $z$ coordinates, provides computational efficiency while maintaining a high degree of accuracy compared to full three-dimensional simulations, as done in most thermal models coupled with P2D models (Alkhedher et al., 2024). The dimensional reduction is justified by the uniform internal heat generation $Q˙$within the cylindrical cell core and a uniform thermal resistance along the angular direction $θ$⁠. As a result, the external surface temperature $T$ remains nearly uniform along the circumference, preventing localized hotspots that are commonly observed in pouch cells (Grandjean et al., 2017; Dileep et al., 2024). This symmetry simplifies the analysis while accurately capturing the dominant heat transfer mechanisms in the radial and axial directions. To simplify the computational domain and reduce computation time, the battery core is treated as a single homogeneous component that encapsulates all internal materials. The thermal conductivity $k$ is considered anisotropic (Huang et al., 2020), differing along the axial $z$ and radial $r$ directions. In particular, the axial thermal conductivity $kz$ is computed through equivalent thermal network theory of layers in parallel, assuming an axial heat flux, while the radial conductivity $kr$ follows a series layer’s model, as detailed in Taheri et al. (2013) and shown in equations (18) and (19). The theory of equivalent thermal networks can be simplified by applying the formulations for flat sheets instead of cylinders, as the battery’s radius of curvature (on the order of millimeters, starting from the mandrel radius $rman$⁠) is significantly larger than the thickness of its layers (in the order of micrometers). The density $ρ$ and specific heat $c$ are determined as volume- and mass-averaged values, respectively, as presented in equations (20) and (21). The battery reference length for the core thermal properties calculation, is the representative layer thickness $dRL$⁠, considering also the current collectors’ thicknesses $dcc,p$ and $dcc,n$⁠:

Heat transfer within the battery follows the energy conservation equation:

To model heat generation inside the core, the formulation proposed by Bernardi et al. (1985), Rao and Newman (1997) and Gu and Wang, (2000) is used. The total heat source in the one-dimensional domain $Q˙1D$⁠, is calculated along the $x$ direction in the electrochemical model domain thickness, expressed in the following equations:

These terms account for contributions from all internal components, including the electrodes, separator and current collectors, detailed in the following bullets:

• Active polarization heat$Q˙act$⁠: Under open-circuit voltage conditions $Ueq$⁠, the system remains in local equilibrium, with lithium-ion potentials in the electrolyte and electrode material matching. However, once current flows, restoring equilibrium requires energy dissipation, resulting in activation polarization heat $Q˙act$⁠.

• Ohmic heat$Q˙ohm$⁠: It arises from electronic and ionic resistances in different phases and is described by Joule’s law, where the first term represents electronic ohmic losses in the solid phase, the second term accounts for ionic losses in the electrolyte and the third term describes ionic migration heating effects.

• Reversible heat$Q˙rev$⁠: Based on the difference in energy levels between reactants and products, the electrochemical reaction either absorbs or releases heat to maintain thermodynamic equilibrium. The phenomenon is associated with entropy variations within the system, which can lead to either exothermic or endothermic reactions, depending on the entropy change. These reactions, in turn, influence the overall heat generation $Q˙1D$⁠, either mitigating or exacerbating it. The effect is strongly dependent on the current density $i$⁠, which changes sign during charge and discharge cycles, and on the SoC through the entropy coefficient $∂Ueq/∂T$ sign, that represents the temperature dependence of the equilibrium potential $Ueq$ and is directly associated with entropy changes.

The heat generation calculated using equation (23) corresponds to the electrochemical model’s computational domain, which does not account for the thicknesses of the current collectors $dcc,p$ and $dcc,n$⁠. Consequently, the computed heat generation does not represent the correct volumetric heat generation. To incorporate the contribution of these components, the heat generation obtained from the P2D model, $Q˙1D$ must be scaled by a correction factor. This factor $χRL$ is defined as the ratio between the representative layer thickness, $dRL$ and the total thickness, including the current collectors, given by the following equation:

Finally, the mandrel and the can are not heat-generating components.

Because ohmic drops along the current collectors are omitted as an assumption, the potential error introduced by this simplification is assessed. The omission of current collector ohmic losses can slightly underestimate the core heat generation. Since the heat generation in the current collectors is exclusively due to ohmic effects, it can be approximated as $Q˙ohm,cc,j=i2σcc·dcc,idRL+dcc,p+dcc,n$ for a current collector $j$⁠, positive or negative. Considering a current density $i = 50 A/m2$ and the relevant thicknesses $d$ as reported in the Sections 3.1–3.2, respectively, and using electrical conductivities $σcc,j$ equal to $37.8 MS/m$ and $59.6 MS/m$ for aluminum (positive) and copper (negative) current collectors respectively (Goutam et al., 2017), the resulting heat generation is $Q˙ohm,cc,p = 4.21×10-6 W/m3$ and $Q˙ohm,cc,n = 1.87×10-6 W/m3$⁠, These values would add to the total heat generation $Q˙$ predicted by the coupled electrochemical–thermal model, but are negligible in magnitude. In fact, they are several orders of magnitude lower than the heat generation obtained in our setup, which is on the order of $104 W/m3$⁠. Therefore, the impact of this simplification is negligible and does not significantly affect the thermal predictions.

A Robin boundary condition is imposed along the external surface of the battery (i.e. the case), characterized by a fixed convective heat transfer coefficient $h$ and an ambient temperature $Tamb$⁠, while temperature continuity is ensured across the different regions. As the initial condition, the entire computational domain is uniformly set to the ambient temperature $Tamb$⁠. A schematic representation of the computational domain, the components and of the applied boundary condition is depicted in Figure 4:

Since the thermal model depends on several electrochemical parameters, and while key properties like the electrolyte diffusion coefficient $De$⁠, the electrolyte conductivity $σe$⁠, the mean molar activity coefficient $f±$ and the potential equilibrium $Ueq$ exhibit temperature dependence, coupling the thermal and electrochemical models is essential for accuracy. To achieve this, the average core surface temperature $Tco,av$ computed from the thermal model at each time step, is used to update the electrochemical model, as described in equations (30) and (31). $De$⁠, $σe$ and $f±$ temperature-dependency are calculated through Arrhenius equation [equation (31)], while the potential equilibrium $Ueq$ is calculated through equation (17):

where $Aco$ is the core area in the 2D domain, $p$ is the temperature-dependent parameter and $Eact,p$ and represents the activation energy associated with the parameter $p$⁠, governing its sensitivity to temperature variation. The updated electrochemical parameters are then fed back into the thermal model to determine heat generation within the core, introducing additional non-linearity into the system, as depicted in the flowchart in Figure 5. This iterative coupling ensures real-time updates of temperature-dependent parameters, maintaining consistency between the thermal and electrochemical models.

When simulating phase change materials, three distinct regions can be identified: a solid phase, a liquid phase and an intermediate mushy region where both phases coexist. In the numerical model, the heat transfer within the liquid phase is assumed to occur primarily through conduction, with natural convection effects being disregarded. This assumption is justified by the small thickness of the PCM layer in the given application, where buoyancy-driven flow can be considered negligible. To validate this simplification, a scaling analysis was performed based on the Rayleigh number:

where, $g$ gravity acceleration, $β$ is the thermal expansion coefficient, using a standard value equal to 10⁻³ for molten salts that will be used in this work, $ΔT$ is the temperature difference, $L$ is the characteristic length, $ν$ is the kinematic viscosity and $α=k/ρc$ is the thermal diffusivity. Two characteristic Rayleigh numbers were calculated using different length and temperature scales, representative of each geometry:

• Internal PCM: The characteristic length was taken as the mandrel radius $rman$⁠, as the heat flow direction is along the radius, and a temperature difference of $ΔT=1°C$ was assumed, based on the near-isothermal behavior of the small mandrel region (from center to outer PCM interface).

• External PCM layer: The characteristic length was taken as the actual PCM thickness $dPCM$used in the external configuration (as discussed later, in Section 3.2). For the temperature gradient, a maximum temperature difference $ΔT=30°C$ was considered, corresponding to a worst-case battery surface temperature versus ambient, with the maximum battery operating temperature equal to 50°C (Rodrigues et al., 2017; Jiang et al., 2020) and the ambient temperature equal to $Tamb=20°C$⁠.

For the internal PCM configuration, the Rayleigh number obtained is equal to 86, while for the external configuration, the Rayleigh number is even smaller. Given that these values are significantly low to the typical critical values (≈1,000) for considering natural convection phenomenon and considering the radial orientation of the heat flow (orthogonal to gravity), the simplification of neglecting natural convection is justified.

In addition, the following assumptions are considered:

• The PCM is treated as a homogeneous material with isotropic characteristics.

• Volume variations during phase transition are negligible, and the formation of air gaps due to expansion or contraction is disregarded.

In the absence of natural convection, the governing energy conservation equation for the PCM, formulated using an enthalpy-based method to account for latent heat effects (Fragnito et al., 2022), reduces to:

where the term $Sh$ is used to model the release or absorption of latent heat during melting or solidification, expressed as:

where $Hl$ is the latent heat of fusion, and $ψ$ is a dimensionless variable that represents the fraction of liquid phase, calculated as a function of temperature $T$⁠:

$Vl$ and $VPCM$ are the liquid fraction volume and the entire PCM volume, respectively, $Tsol$ is the solidification temperature limit and $Tl$ is the melting temperature limit. Continuity in the first derivative is maintained at the interfaces between the fully solid phase and the phase transition region, as well as between the fully liquid phase and the phase transition region. This prevents abrupt variations in the temperature $T$ and avoids discontinuities.

The material properties are defined as temperature-dependent functions, expressed in terms of the liquid fraction $ψ(T)$⁠, allowing for a continuous and thermodynamically consistent transition between solid and liquid phases. This approach is necessary because thermophysical properties such as density $ρ$⁠, thermal conductivity $k$ and specific heat $c$ undergo significant variation during phase change and a smooth interpolation based on $ψ(T)$⁠, ensures realistic representation of the transient melting behavior:

In this section a detailed description is provided of the electrical, chemical and thermal properties of the analyzed battery, as well as the thermal characteristics of the phase change material used in the study. In addition, the setup used for validation with the results from (O’Regan et al., 2022; Cianciullo et al., 2022) is outlined and are discussed the different configurations designed to assess the benefits of PCM integration – internally within a hollow mandrel, externally around the battery and in a combined internal–external arrangement. Finally, a study on the optimal PCM thickness in the external and internal–external configuration is performed to minimize the overall battery pack volume while ensuring that the temperature remains within safe operating limits.

The battery analyzed is a LCO 18650 cylindrical battery with a capacity $C$ of 3.5 Ah and a nominal voltage $Unom = 4.2 V$⁠. The dimensions that are necessary for both the electrochemical and the thermal model are depicted in Figure 6(a) and (b), respectively, while the battery electrochemical–thermal properties, and the phase change material thermal properties are detailed in Tables 1–3. In particular, the phase change material used is the sodium thiosulfate pentahydrate (STP), a molten salt with key thermal parameters detailed in references (Moynihan, 1966; Hadjieva et al., 2000; Landini et al., 2019). The selected PCM was chosen for its favorable thermal properties for battery applications. It exhibits a melting temperature range (40–48°C), which is close to the maximum limit $Tlim = 50°C$ recommended by manufacturers (Jiang et al., 2020), high latent heat (210 kJ/kg), chemical stability and non-flammability. These features make it a suitable candidate for compact passive thermal protection systems in lithium-ion cells. As the battery temperature reaches the phase change region of the PCM, the latent heat $Hl$ acts as a thermal buffer by maintaining the PCM temperature nearly constant. During this period, the absorbed latent heat temporarily stores the energy dissipated by the battery, thereby slowing the rise in cell temperature $∂T/∂t$ and stabilizing the thermal behavior of the system. This extended thermal regulation provides a longer response time to implement protective measures, reducing cell degradation and ensuring safety for individuals in proximity to the battery. Molten salts were chosen over paraffins because paraffins are highly flammable and may pose safety hazards near batteries.

The thermal properties of each core battery component are specified, taken from Cianciullo et al. (2022), and these have been used in the calculation of the core’s overall thermal properties. The electrochemical parameters are derived from (Melcher et al., 2016), where the electrolyte is LiPF₆ in 1:1 EC:DEC, with the electrolyte phase diffusion coefficient $De$⁠, phase conductivity $φe$⁠, mean molar activity coefficient $f±$ and the entropy coefficient $∂Ueq∂T$ of both positive and negative electrodes have been evaluated through the COMSOL Multiphysics material library of the battery components, listed in Table 1. These properties have been used both for the numerical comparison with (Cianciullo et al., 2022) and for the simulations presented in this study (Table 1 to 3).

To ensure consistency and reliability of the proposed coupled electrochemical–thermal model, a validation was performed against the experimental results provided by (O’Regan et al., 2022). The simulation replicates the same operating conditions adopted in their study: constant 1C discharge (where 1C represents the applied current required to fully discharge the battery within one hour), ambient temperature $Tamb = 25°C$ (considered also as initial temperature) and a convective heat transfer coefficient $h = 15 W/m2K$⁠. Following the full discharge phase, the model also simulates the subsequent rest period up to 10,000 s, allowing for a direct comparison of the cell’s thermal relaxation behavior. The electrochemical and thermal properties of the battery used in this validation case are set equal to those reported in (O’Regan et al., 2022), to ensure full compatibility with the experimental reference. The comparison focuses on the evolution of the average surface temperature $Tsur,avg$ as the experimental studies do not specify the exact location of the thermocouple on the cell surface. Therefore, the simulated value is taken as the average temperature along the cylindrical surface of the cell. Given the nearly isothermal behavior of the surface – supported by the relatively high axial thermal conductivity $kz$ – this approximation introduces minimal error and remains representative of the measured experimental quantity. Then, numerical results for the average battery temperature are compared with those reported in (Cianciullo et al., 2022), which validated their model through the thermal runaway ignition time $tTR$ – the time required for a battery to reach the initiation of thermal runaway – with a previously validated study using the P2D model (Melcher et al., 2016). A charge/discharge cycle is simulated, using battery characteristics reported in Section 3.1, lasting 500 s per full cycle, for a time period $t$ = 7,000 s, beginning with the charging phase at a current density $i$ of ±80A/m², as depicted in Figure 7. The analysis considers two different convective heat transfer coefficients, $h$⁠, set to 4 W/m²K and 6 W/m²K and an ambient temperature $Tamb$ of 20°C, corresponding to a natural convection with air. Assuming an electrode surface area of 0.1 m², as determined in (Saw et al., 2013), the current density at 1C is 35A/m². In this study, the comparison is conducted at 2.3C.

The different configurations analyzed are a function of the PCM position: internally within a hollow mandrel, externally around the battery can and together in an internal/external configuration, as shown in Figure 8. The volume of the internal PCM is equal to that of the mandrel, under the assumption that the enclosing container prevents any PCM leakage into the battery, thereby avoiding direct contact. The containers for both internal and external PCM are not explicitly modeled, as their thickness – on the order of a tenth of a millimeter – make their thermal influence negligible. The comparison between internal and external PCM configurations is conducted for an equivalent PCM volume, leading to an external PCM layer thickness $d$ of 0.22 mm. This assumption allows for an assessment of the potential advantages of an internal PCM relative to an external one, acknowledging that such a thickness is minimal in the context of a battery pack. However, the integration of an internal PCM would not contribute to an increase in the outer casing thickness. When scaled across multiple cells within an industrial battery pack, this distinction could be significant in terms of overall footprint. For fair comparison across configurations (internal, external, internal/external), volume constraints were applied where necessary, ensuring that improvements in thermal performance are not merely due to additional PCM volume, which would naturally lead to better thermal management, but rather to more efficient spatial placement. At the interface between the battery surface and the PCM domain (internal, external or combined), temperature continuity is ensured. This condition accounts for heat transfer between the cell and the adjacent PCM material, effectively coupling the thermal behavior of both regions. Then, for the external PCM, the same Robin boundary condition is applied as for the reference case. The analysis focuses on the average temperature $T$ plot by each configuration, the time required for the battery to reach the temperature limit $Tlim$ of 50°C and on the volume averaged PCM melting fraction $ψ$⁠, aiming to determine which configuration offers superior thermal management while optimizing battery efficiency and spatial constraints for potential battery pack applications. To compare the time to reach $Tlim$⁠, a time improvement percentage is defined as in equation (37):

where $tstop$ is the time when the maximum temperature $Tmax$ reaches the temperature limit $Tlim$ for the configuration under analysis, while $tstop,ref$ represents the corresponding time $tstop$ for the compared reference configuration. In this study, the reference configuration represents the single battery without PCM. The studies have been done with the same pulse discharge/charge cycle from the comparison with the work of (Cianciullo et al., 2022) for a time $t$ of 7,000 s, but a different current density $i$ of ±50A/m² (i.e. 1.5C), that correspond to a standard current for industrial applications. The corresponding boundary conditions used are an ambient temperature $Tamb$ of 20°C, and two different heat transfer coefficients $h$ = 8 W/m²K and 10 W/m²K, corresponding to two different natural convection conditions in air.

In addition, a parametric analysis that varies the external PCM thickness $tPCM$⁠, was conducted to determine the optimal thickness $topt$ required to prevent the battery from exceeding the temperature limit $Tlim$ = 50°C and the operational temperature limit $Tlim,op$⁠, set at 45°C. This threshold was set as a safety measure that battery customers could apply for designing their battery packs, to prevent the battery from reaching the temperature limit $Tlim$⁠, to mitigate the risk of potential damage and ensure safer operating conditions. The simulations have been performed with an external ambient temperature $Tamb$ = 20°C, and a heat transfer coefficient $h$ = 10 W/m²K, for a time $t$ of 7,000 s, making also a grid consistence analysis to assure the validity of the results. The test was performed on the most challenging configuration, with a 2 mm PCM external layer, in the internal and external configuration, making it the most sensitive to discretization errors. Based on the outcome of this test, grid independence is assumed to hold for all other configurations as well, including those with the mandrel and with reduced PCM thicknesses.

Since the simulation time is extended, with this heat transfer coefficient $h$⁠, the configuration goes to a thermal equilibrium after a certain period, and for this reason a trend of the maximum temperature $Tmax$ (i.e. the final temperature, at $t = 7,000 s$⁠) as a function of external PCM thickness $tPCM$ in both external and internal/external configuration is presented. The objective is to assess the spatial benefits of incorporating an internal PCM, optimizing the battery pack design by reducing its footprint, quantifying a volume saving percentage, expressed in equation (38). The parametric analysis focuses on the variation of PCM thickness, as this geometric parameter directly impacts the spatial footprint of the battery unit. This study aims to determine the minimum external PCM thickness in the internal/external configuration needed to keep the cell temperature within safe limits, and to compare it with the external-only configuration yielding the same maximum temperature $Tmax$⁠:

where $rref2$ is the total radius (battery+PCM) of the reference system to be compared.

The validation of the model with the results of O’Regan et al. (2022) and Cianciullo et al. (2022) are presented in Figures 9–10 terms of the average surface temperature $Tsur,avg$ and of the maximum battery temperature $Tmax$⁠, respectively, for different convective heat transfer coefficients $h$⁠. For the experimental results of O’Regan et al., the temperature trend is consistent during both the discharge and relaxation phases, as shown in Figure 9. The only discrepancy is that the experimental battery temperature drops below 25°C, indicated as the ambient temperature $Tamb$⁠, whereas in the simulation it does not, since 25°C is assumed as the minimum possible temperature. For the results of Cianciullo et al. (2022), the temperature trends indicate significantly high values in both configurations (Figure 10), suggesting that the thermal runaway mechanism is likely to be triggered. The similarity of the temperature curves $T$ confirms that the model accurately represents the typical thermal behavior of a battery. The observed fluctuating temperature trend $T$ results from the alternating charge and discharge phases, which cause variations in current density $i$ and, consequently, in reversible heat generation $Q˙rev$⁠. Since the entropy coefficient $∂Ueq/∂T$ is predominantly positive, the system experiences a lower temperature rise $∂T/∂t$ during discharging, as the entropy-driven heat contribution partially offsets the heat generated during the process As a result, the total heat generation is partially offset, leading to a reduced temperature increase $∂T/∂t$ during discharge. To better quantify the model’s accuracy, both the root mean square (RMS) error and the maximum error were evaluated, as shown in Table 4. For the case by O’Regan et al. (2022) a RMS error of 1.01°C and a maximum error of 2.50°C were observed. Although experimental uncertainty values are not explicitly reported in the cited study, a typical thermocouple accuracy of ± 1°C is assumed, consistent with standard practices of thermocouples. For the comparison against (Cianciullo et al., 2022), the RMS error was 1.15°C for $h = 4 W/m2K$ and 2.02°C for $h = 6 W/m2K$⁠, with corresponding maximum errors of 3.93°C and 4.01°C. Although these error values may appear relatively high, they occur over very wide temperature ranges (up to 223°C), which naturally amplifies the absolute deviation. This trend highlights how the error magnitude increases with the thermal excursion yet remains acceptable relative to the simulated and measured temperature scales. Remaining discrepancies may result from several factors, including uncertainties in material properties and modeling assumptions. Key thermophysical and electrochemical parameters – such as thermal conductivity $k$⁠, specific heat capacity $c$⁠, density $ρ$ and electrical conductivity $σ$⁠, – may vary with temperature $T$⁠, yet are typically treated as constant in numerical models. In addition, electrochemical properties such as the entropy coefficient $∂Ueq/∂T$ and ionic diffusivity $σe$ are sensitive to the SoC, aging effects and electrode microstructure, which may deviate from literature values or experimental results. These deviations, combined with the simplification of a uniform convective boundary conditions may contribute to the residual error observed between simulation and experiment. Nevertheless, considering all these factors and the overall agreement with multiple data sets, the model can be considered validated for the intended thermal analysis applications (Table 4).

To verify mesh independence, a sensitivity analysis was performed for the internal/external PCM configuration, using an external PCM thickness $d$ of 2 mm The monitored variable was the temperature at the mid-height along the battery axis, in the region originally occupied by the mandrel, now filled with PCM (⁠$Tman,max$⁠), which corresponds to the thermal hotspot in this configuration. The temperature was evaluated at the end of the charge/discharge cycle, when the accumulated heat is maximal. As shown in Figure 11, mesh refinement leads to a monotonic decrease in $Tman,max$⁠. The goal of the mesh refinement study was to achieve convergence of $Tman,max$up to the second decimal place. As shown in Figure 11, the difference between the two finest meshes (30,260 and 62,697 elements) is zero at this precision, indicating no variation in the second decimal digit. Since further refinement does not affect the result at the target precision, the mesh with 30,260 elements is selected. This choice balances numerical accuracy with computational efficiency and confirms that the solution is mesh-independent at the chosen precision. This mesh is therefore used for all subsequent simulations.

The maximum temperature evolution over time is analyzed for the three configurations under investigation: internal, external and combined internal/external PCM, along with a reference case featuring a single battery configuration without PCM. The results are presented for two convective heat transfer coefficients $h = 8 W/m2K$ and $h = 10 W/m2K$ in Figures 12(a) and 13(a). As expected, the reference configuration exhibits the highest temperatures, as it does not benefit from the latent heat absorption associated with the PCM phase transition. Conversely, the internal/external configuration yields the best thermal performance, as it effectively uses twice the PCM volume compared to the other two configurations. The internal and external PCM configurations produce similar temperature trends, with the internal PCM reaching the temperature limit $Tlim$ slightly earlier. The superior performance of the internal/external PCM configuration can be explained by a more favorable redistribution of latent heat buffering. While the internal PCM offers high thermal contact with the core, it lacks direct exposure to the ambient, leading to slower heat rejection. On the contrary, the external PCM benefits from enhanced convective cooling due to its position but responds later to temperature spikes generated at the battery core. By combining the two, the system achieves a more balanced thermal gradient, which slows down the propagation of heat from the center and prolongs the phase change process. This effect is particularly visible in the time-dependent melting fraction curve, in Figures 12(b) and 13(b). The rate at which the PCM transitions varies $∂ψ/∂t$ depends on the applied charge or discharge, as these influence heat generation $Q˙$ through the reversible heat effects $Q˙rev$⁠, which are positive during charge and negative during discharge, thereby reducing or increasing the overall heat generation $Q˙$⁠. As soon as the phase transition begins, approximately at 1,000 s, the temperature rise $∂T/∂t$ significantly decreases compared to the case without PCM. This effect is due to the latent heat absorption $Hl$⁠, which creates a temperature plateau as the system approaches the limit temperature region $Tlim$⁠. In all cases analyzed, the temperature limit $Tlim$ is exceeded. Although the absolute reduction in peak temperature is moderate under the given electrochemical and thermal conditions, the PCM demonstrates its effectiveness by significantly extending the safe operating time of the battery system. This behavior is expected in systems without strong active cooling, where the PCM acts primarily as a thermal buffer. Through latent heat absorption, it delays the onset of critical temperatures rather than drastically reducing the peak temperature. Therefore, the thermal benefit of PCM is more meaningfully evaluated in terms of time gained before reaching the limit temperature $Tlim$⁠, rather than the magnitude of temperature drop itself. This is clearly reflected in the time improvement $It$ to reach the limit temperature $tstop$⁠. The external PCM, compared to the internal configurations, performs marginally better in terms of thermal management efficiency. As reported in Table 5, for $h = 8 W/m2K$⁠, the time improvement $It$ in reaching the limit temperature is 32.8% for the internal PCM, 48.2% for the external PCM and 78.8% for the internal/external PCM configuration, compared to the reference case. Similarly, Table 6 shows that for $h = 10 W/m2K$⁠, the time improvements $It$ increase to 40.5%, 58.6% and 120.5%, respectively. This trend confirms that, in all configurations, a higher convective heat transfer coefficient delays the attainment of the limit temperature $Tlim$ is reached later for higher convective heat transfer coefficients $h$ across all configurations. This is because increasing$h$ enhances heat dissipation to the surroundings, particularly in configurations where the PCM is placed on the outer side, which in turn slows down the phase transition rate $∂ψ/∂t$⁠. As a result, the PCM remains in the phase change regime for a longer period, thereby providing more effective thermal buffering. Consequently, the thermal performance is improved for $h = 10 W/m2K$ as the prolonged phase transition extends the PCM’s temperature-regulating capacity. In summary, these results confirm the higher effectiveness of the internal/external arrangement, which benefits both from increased PCM volume and synergistic heat dissipation dynamics. In addition, the duration of the phase transition is longer in the internal–external PCM configuration than in the other cases, as it benefits from a greater PCM volume. However, if the objective is space optimization for battery packs, the internal PCM configuration provides comparable thermal performance while improving spatial integration, obtaining a volume saving $Vs$ of 4.7% compared to the PCM external configuration.