Journal of Energy Storage 81 (2024) 110482

Contents lists available at ScienceDirect

Journal of Energy Storage

journal homepage: www.elsevier.com/locate/est

Research Papers

Numerical investigation on thermal performance of phase change materials

embedded in functionally graded metal foam
Tianyu Si a, Wei Cui a, Ting Ma a, Lin Lu b, Qiuwang Wang a, *
a
Key Laboratory of Thermo-Fluid Science and Engineering, MOE, Xi'an Jiaotong University, Xi'an, Shaanxi 710049, PR China
b
Department of Building Environment and Energy Engineering, The Hong Kong Polytechnic University, 999077, Hong Kong

A R T I C L E I N F O A B S T R A C T

Keywords: Latent heat thermal energy storage (LHTES) technology based on phase change materials (PCMs) embedded in
Phase change material metal foam has recently become a popular method to store thermal energy and release it to other processes with
Thermal energy storage thermal energy demand, which has been applied in solar energy, electronic devices, waste heat recovery, etc. In
Metal foam
this study, we proposed a novel design of metal foam with porosity changing over the horizontal and vertical
Linear porosity
coordinates. The thermal resistance model method was applied to explore the heat transfer mechanism and linear
porosity was proved to be able to improve the convection intensity in different stages of melting. Also, metal
foam with horizontal linear porosity or vertical linear porosity can both shorten the total melting time and in­
crease the average power density, while composite linear porosity can combine the enhancing effect of these two
structures and achieve an optimal enhancing effect. The optimal linear porosity gradient can enhance the
average power density by 15.8 % compared to those of uniform porosity. Furthermore, the heat storage per­
formance of linear porosity with different enclosure aspect ratios was also discussed. The linear porosity
structure was proved to be able to balance the nonuniformity of melting caused by natural convection.

thermal management and energy storage at a relatively stable

1. Introduction temperature.
However, most PCMs suffer from disadvantage of low thermal con­
Since the issues of energy shortage and environmental problems are ductivity [11], generally below 0.5 W⋅m− 2⋅K− 1, which may result in a
getting increasingly prominent and severe, the development of new low charging rate and undesirable temperature rise in LHTES systems.
energy resources and energy conservation is becoming important. The heat transfer efficiency of LHTES system can be enhanced by adding
Thermal energy storage (TES) system has recently become a popular media with high thermal conductivity such as metal fins [12], nano-
method to store surplus energy and offer thermal management capable particles [13,14], expanded graphite [15,16], and metal foams
of absorbing and releasing heat [1]. High efficient TES systems have [17–19], etc., into PCMs or using encapsulated PCMs [20]. Among these
been investigated to deal with the discrepancy between energy gener­ techniques, metal foam, as a common porous media, is a valuable
ation and consumption [2,3] in broad applications such as building enhancer due to its desirable properties of large specific surface area and
energy efficiency [4], solar energy storage [5], thermal management of good heat conducting capacity [18]. Via experimental and numerical
electric devices [6], industrial waste heat recovery [7], etc. investigations [21,22], composite PCMs (CPCMs), fabricated by filling
TES technologies can be divided into sensible heat thermal energy metal foam with PCMs, have been approved to possess both better heat
storage (SHTES), latent heat thermal energy storage (LHTES), and transfer rate and high energy storage density. Furthermore, it should be
thermochemical heat storage (TCHS) [8]. Among these techniques, the noted that the system parameters of the metal foams, especially the
LHTES systems based on PCMs appear to be more efficient for their high porosity, have a crucial influence on the heat transfer behavior of CPCMs
thermal storage density and stability in operating temperature [9]. With [23]. It is vital to find out the appropriate porosity for the engineering
high latent heat values, PCMs can achieve higher energy density design of LHTES systems. Xiao et al. [18] experimentally tested the
compared with SHTES materials [10]. Thus, the LHTES systems tend to effective thermal conductivities of the CPCMs prepared by paraffin and
be smaller than the SHTES systems. Furthermore, with the stable different metal foams. The thermal conductivities of the CPCMs were
working temperature of PCMs, the LHTES system can also be utilized for enhanced by over thirty times larger than that of pure paraffin. And the

* Corresponding author.
E-mail address: wangqw@mail.xjtu.edu.cn (Q. Wang).

https://doi.org/10.1016/j.est.2024.110482
Received 21 October 2023; Received in revised form 31 December 2023; Accepted 3 January 2024
Available online 11 January 2024
2352-152X/© 2024 Elsevier Ltd. All rights reserved.
T. Si et al. Journal of Energy Storage 81 (2024) 110482

Nomenclature σ dimensionless foam ligament radius

μ viscosity (kg⋅s− 1⋅m− 1)
Amush mushy zone constant (kg⋅m− 3⋅s) γ thermal expansion coefficient (K− 1)
cp specific heat capacity (J⋅kg− 1⋅K− 1) β liquid fraction
C inertia coefficient δ small constant
df fiber diameter (m) χ tortuosity coefficient
dp pore diameter (m) ω pore density
e geometric variable
E energy storage density (J⋅m− 3) Subscripts
Fo Fourier number a ambient
h heat transfer coefficient (W⋅m− 2⋅K− 1) ave. average
hsf interstitial heat transfer coefficient (W⋅m− 2⋅K− 1) cond conduction
k thermal conductivity (W⋅m− 1⋅K− 1) conv convection
K permeability (m2) eff effective
L latent heat (J⋅kg− 1) f fluid region
Lx enclosure length (m) h heating side
Ly enclosure height (m) m melted
Nu Nusselt number max maximum
p pressure (Pa) mf metal foam
P power density (W) PCMs phase change materials
Pr Prandtl number tot total
Q average heat flux (W⋅m− 2) w wall
R thermal resistance (K⋅W− 1) x in the horizontal direction
Re Reynolds number y in the vertical direction
t time (s) Abbreviations
T temperature (K) LHTES latent heat thermal energy storage
Tm melting temperature (K) PCMs phase change materials
Ts solidification temperature (K) TES thermal energy storage
u, v horizontal and linear velocity components (m⋅s− 1) SHTES sensible heat thermal energy storage
vh flowing velocity TCHS thermochemical heat storage
U velocity vector (m⋅s− 1) CPCMs composite phase change materials
V volume (m3) HTF heat transfer fluid
x, y Cartesian coordinates (m) PPI pores per inch
Greek symbols UP uniform porosity
α thermal diffusivity (m2⋅s− 1) HLP horizontal linear porosity
λ linear porosity gradient (m− 1) VLP vertical linear porosity
ε porosity CLP composite linear porosity
ρ density (kg⋅m− 3)

enhancement ratio of thermal conductivity is enlarged with the decrease novel layouts of metal foams employed in CPCMs.
of the porosity. Li et al. [24] investigated the effects of porosity on the Recently, metal foam with non-uniform structures of porosity has
wall temperature and the temperature nonuniformity inside the CPCMs. been studied. The non-uniform arrangement of porosity indicates the
Results show that the addition of copper foam in PCMs brings about both purpose of improving the heat transfer performance in the local area,
the enhancement of effective thermal conductivity and the suppression thus achieving an overall improvement in heat transfer performance.
of natural convection of liquid PCMs. Tao et al. [25] studied the effects Normally, the non-uniform arrangement types can be divided into par­
of metal foams porosity and pore density on the melting rate, heat tial filling [25,26], gradient porosity [29,30], and linear porosity
storage capacity, and heat storage density of PCMs. It is found that with [31,32]. Li et al. [26] proposed a nonuniform structure with a region of
the increase of porosity, CPCMs melting rate is improved and the ther­ lower porosity than elsewhere. This structure provides higher effective
mal energy storage density keeps stable, but the power density is greatly thermal conductivity of CPCMs and accelerates the development of the
reduced. Fteiti et al. [27] used a pseudo-random generator to create melting interface in the later stage of melting, which helps to achieve
structures of random porosity in composite metal foam-phase change high power density and energy storage density. Wang et al. [33]
materials, which revealed that random porous medium could provide a experimentally explored the thermal performances of LHTES units with
lower rate of PCM melting for a similar average porosity. Zhang et al. uniform porosity and horizontal gradient porosity of metal foam. The
[28] investigated the effect of porosity and pore density on metling average porosity of metal foam is kept unchanged to avoid the reduction
behavior of CPCMs. The thermal performance of the CPCMs is greatly in heat storage density. It is found that the application of gradient
affected by porosity than pore density, and the thermal storage rate porosity can generate a reduction of melting time by 37.6 % and
increased with porosity decreased. Based on the research mentioned improve the temperature uniformity at the same time compared with the
above, on the one hand, by simply reducing the porosity of metal foams, one with uniform metal foam. Zheng et al. [29] proposed a structure
the enhancement ratio of effective thermal conductivity can be with horizontal and vertical gradient porosity. The horizontal gradient
improved. On the other hand, the problems of temperature nonunifor­ porosity can decrease the thermal resistance of the unit while the ver­
mity and reduction of heat storage density still exist. To further improve tical gradient porosity effectively fix the problem of melting nonuni­
the thermal performance of LHTES systems, it is necessary to look for formity brought by natural convection. CPCMs with gradient porosity in

2
T. Si et al. Journal of Energy Storage 81 (2024) 110482

both horizontal and vertical directions possessed the shortest melting aspect ratio (Lx/Ly) are discussed by changing Lx (Lx = 10, 30, 50, 70, 90
time, which is 12.07 % shorter than uniform structure. mm). Initially, the temperature of the CPCM is set to be equal to the
Nevertheless, to the best knowledge of authors, the effect of linear ambient temperature (Ta = 25 ◦ C). During simulation, the HTF (Th =
porosity on the heat transfer mechanism of LHTES system was rarely 50 ◦ C) flows through the heating wall at a constant flowing velocity (vh
considered and the thermal resistance model method was rarely applied. = 1.6 m/s). The thermophysical properties of octadecane, copper, and
Furthermore, there is a lack of design strategy for horizontal and vertical water are listed in Table 1.
linear porosity for engineering applications. Therefore, in this study, a The copper foam is designed to have a pore size of 10 PPI (pores per
LHTES system filled with metal foam with horizontal and vertical linear inch) and a changing porosity (ε) based on the horizontal and vertical
porosity under the effect of heating wall is proposed. To evaluate the coordinates. We propose the following equation to describe the layout of
ability of the LHTES system to store and release surplus energy reliably the porosity:
and efficiently, numerical simulations were performed to find out the /
x − Lx /2 y − Ly 2
impact of linear porosity on the melting behavior, heat transfer mech­ ε(x, y) = εave + λx + λy (1)
anism and thermal storage performance of LHTES systems. From an Lx Ly
engineering perspective, the impact mechanism of enclosure aspect ratio
where λx is the porosity gradient in the horizontal direction. λy is the
on natural convection and porosity design strategy is also analyzed.
porosity gradient in the vertical direction. εave is the average porosity,
taken as 0.8. Four different basic layouts of copper foam are proposed,
2. Physical and mathematical model

2.1. Physical model

Table 1
Thermophysical properties of n-octadecane [11], copper, and water.
A schematic of the LHTES system filled with CPCM is depicted in
Fig. 1 (a). A 2D model of a rectangular enclosure with a heating wall is Thermophysical properties n-Octadecane Copper Water

depicted on the left side. A kind of paraffin wax, n-Octadecane (C18H38) ρ/kg⋅m− 3
775 8978 988.1
which has a melting temperature of 301.4 K [11] is selected as PCM and k/W⋅m− 1⋅K− 1 0.2 387.6 0.65
cp/J⋅kg− 1⋅K− 1 2400 381 4178
copper foam is composited to enhance the heat conduction. Water is
μ/kg⋅s− 1⋅m− 1 0.003 0.00055
selected as the HTF flowing through the heating surface. The enclosure L/J⋅kg− 1 245,000
length Lx and enclosure height Ly remain constant (Lx × Ly = 60 × 120 Tm/K 301.4
mm2) in Sections 4.1 and 4.2. In Section 4.3, the cases with different γ/K− 1 0.001

Fig. 1. Schematics of (a) the LHTES system filled with CPCMs, (b) metal foam with uniform structure (λx = λy = 0), (c) metal foam with horizontal linear porosity (λx
> 0, λy = 0), (d) metal foam with vertical linear porosity (λx = 0, λy > 0), (e) metal foam with composite linear porosity (λx > 0, λy > 0).

3
T. Si et al. Journal of Energy Storage 81 (2024) 110482

including uniform porosity (UP, λx = λy = 0), horizontal linear porosity [ ] ∂T ( )

( ) ( ) ( ) ∂T ∂T
(HLP, λx > 0, λy = 0), vertical linear porosity (VLP, λx = 0, λy > 0) and (1 − ε) ρcp mf + ε ρcp PCMs + ε ρcp PCMs u + v
composite linear porosity (CLP, λx > 0, λy > 0), which are depicted in ∂t ∂x ∂y
(5)
Fig. 1 (b), (c), (d) and (e) respectively. The parameter settings of all the 2 ∂β
= keff ∇ T − ερL
cases are listed in Table 2, while εmax is the maximum porosity which is ∂t
taken when λx = Lx, λy = Ly. Considering the limitation of manufacturing The fourth term on the right side in Eq. (3) is the factor of resistance
and geometry, εmax is set to be no >0.96 in all cases. considered in the mushy zone. Amush is mushy zone constant, taken as
105. δ is a constant term that is set to be 10− 3 [35]. The liquid fraction of
2.2. Mathematical model PCMs (β) equals 0 when the PCMs are unmelted and equals 1 when the
PCMs are completely melted. In the mushy region, β lies between 0 and
To depict the transient process of the melting process of the CPCMs, 1. The value of β can be calculated by Eq. (4).
several governing equations are developed and the following assump­ The permeability (K), inertial coefficient (C), and interfacial surface
tions are accepted to simplify the mathematical model. area (Asf) mentioned above are influenced by metal foam parameters
(1) The liquid flow is considered incompressible and laminar; like porosity (ε), pore density (ω), pore diameter (dp), and fiber diameter
(2) The buoyancy-driven flow is subject to the Boussinesq (df), which can be determined by the following equations [36].
approximation. ( )
(3) Thermophysical properties of copper foam and PCMs remain 4π 1 − 1
χ = 2 + 2cos + cos (2ε − 1) (7)
constant. 3 3
(4) The volume change during the melting process is negligible.
To simulate the melting process of the PCMs, the enthalpy porosity K=
ε2 dk2
(8)
method is utilized. and the thermal equilibrium model is used to simu­ 36χ (χ − 1)
late the heat transfer between metal foam and PCMs. The governing ( )−
equations and boundary conditions are as follows [21]:
1.63
1 df
C = 0.00212 √̅̅̅̅(1 − ε)− 0.132
(9)
Continuity equation: K dp

∂ρ ( →)
(2) χ
∂t
+ ρ ∇⋅U = 0 dk = dp (10)
3− χ
Momentum equation:
3
22.4 × 10−
(
ρ ∂u ρ ∂(uu) ∂(uv)
) dp = (11)
+ + = ω
ε ∂t ε2 ∂x ∂y √̅̅̅̅̅̅̅̅̅̅̅
( ) 1− ε
−
∂p μ 2
+ ∇ u−
μ ρC (1 − β)2
+ √̅̅̅̅ |u| u − Amush 3 u df = 1.18 dp (12)
∂x ε K 3π
K β +δ
( ) (3)
ρ ∂v ρ ∂(vu) ∂(vv) where dp is the pore diameter, df is the fiber diameter and χ is the tor­
+ + =
ε ∂t ε2 ∂x ∂y tuosity coefficient. The effective thermal conductivity (keff) is used in the
( ) energy equations to represent the thermal conductivity of PCMs and
∂p μ 2 μ ρC (1 − β)2
− + ∇ v− + √̅̅̅̅ |v| v − Amush 3 v + ρgγ(T − Tm ) metal foam, which is calculated by the following equations developed by
∂x ε K K β +δ
Boomsma and Poulikakos [37]:
⎧
⎪ β=1 T > Tm 1
⎪
⎪
⎨ keff = √̅̅̅ (14)
T − Ts 2(RA + RB + RC + RD )
β= Ts < T < Tm (4)
⎪
⎪ Tm − Ts
⎪
⎩ 4σ
β=0 T < Ts RA = (15)
(2e2 + πσ(1 − e) )kmf + (4 − 2e2 − πσ(1 − e) )kPCMs
Energy equations:
(e − 2σ )2
RB = (16)
(e − 2σ )e2 kmf + (2e − 4σ − (e − 2σ )e2 )kPCMs

( √̅̅̅ )2
2 − 2e
RC = ( √̅̅
̅ ) ( √̅̅̅ ( √̅̅̅ ) ) (17)
2πσ 2 1 − 2 2 e kmf + 2 2 − 2e − πσ2 1 − 2 2 e kPCMs
Table 2
Parameter settings of the linear porosity for each case.
2e
RD = (18)
Case λx λy εmax εave e2 kmf + (4 − e2 )kPCMs
UP 1–1 0 0 0.80 0.8 √̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
HLP 2–1 0.08 0 0.84 0.8 √√̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅̅
√ ( √̅̅̅ )̅
2–2 0.16 0 0.88 √ 2 2 − (5/8)e3 2 − 2ε
2–3 0.24 0 0.92 σ=√ ( √̅̅̅ ) , e = 0.339 (19)
2–4 0.32 0 0.96
π 3− 4 2e − e
VLP 3–1 0 0.08 0.84 0.8
3–2 0 0.16 0.88 where kmf and kPCMs is the thermal conductivity of metal foam and PCMs
3–3 0 0.24 0.92
respectively. σ is the dimensionless foam ligament radius.
3–4 0 0.32 0.96
CLP 4–1 0.08 0.24 0.96 0.8
The initial condition and boundary conditions are set as follows:
4–2 0.16 0.16 0.96 Initial condition:
4–3 0.24 0.08 0.96
4–4 0.08 0.16 0.92 T(x, y) = Ta , t = 0, 0 ≤ x ≤ Lx , 0 ≤ y ≤ Ly (20)
4–5 0.16 0.08 0.92
4–6 0.08 0.08 0.88
Boundary conditions:

4
T. Si et al. Journal of Energy Storage 81 (2024) 110482

⎧
⎪ ∂T porosity. Nine temperature probes were located in the units and one of
⎪ = hh (Th − T) 0 ≤ y ≤ Ly , x = 0
⎪
⎪
⎪
⎪ ∂n them was selected to verify the simulation results. The temperature
⎪
⎪
⎪
⎪ ∂T curves of the specified point in the case of gradient porosity and uniform
⎪
⎪
⎨
∂n
=0 0 ≤ y ≤ Ly , x = Lx porosity are depicted in Figs. 3 (a) and (b) respectively. It can be found
(21) that the simulation results are reasonably agreeable. Therefore, the
⎪ ∂T
⎪
⎪
⎪
⎪
⎪ ∂n
=0 0 ≤ x ≤ Lx , y = Ly simulation model applied in this paper with nonuniform metal foams
⎪
⎪
⎪
⎪ and the HTF boundary conditions can be verified. More detailed dis­
⎪ ∂T
⎪
⎩ =0 0 ≤ x ≤ Lx , y = 0 cussion was presented in our previous paper [34], which included the
∂n
comparisons with experiment results.
kw ρ vh Ly cp,w μw
hh = 0.664Prw 1/2 Rew 1/2 , Rew = w , Prw = (22) 4. Results and discussion
Ly μw kw

In this paper, the melting process and thermal storage performance

3. Numerical procedure and model validation
of CPCMs in different working conditions are comprehensively studied.
The melting interfaces and liquid fraction curves are characterized, and
3.1. Numerical procedure
thermal storage performance is discussed under different layouts and
working conditions. The results and discussion will include three sec­
To governing equations raised in the last section are solved using
tions: (1) analysis of the melting behavior and thermal resistance of
ANSYS FLUENT. The governing equations are discretized by the finite
CPCMs with linear porosity in different directions; (2) effect of linear
volume method (FVM) with a double precision solver. Second Order
porosity structures on the thermal storage performance of CPCMs; (3)
Implicit Method for Pressure Linked Equations (SIMPLE) algorithm is
comparison of the thermal storage performance of CPCMs with different
adopted for the pressure-velocity coupling. A user-defined function
aspect ratios and analysis of the intrinsic connections between the
(UDF) is utilized to calculate all the parameters of metal foam changing
porous structures and melting nonuniformity.
depending on the porosity. Residual convergence values for continuity,
momentum, and energy equations are set to be 10− 6, 10− 6, and 10− 10
4.1. Analysis of melting behavior and thermal resistance of CPCMs with
respectively.
linear porosity
Different grid sizes are decided and tested to find out the appropriate
grid size, and so as the time step. Fig. 2 (a) shows the temperature curve
In this section, the effect of linear porosity on the melting behavior
of a specified point (x = 0.1 m, y = 0.6 m) in the case of 80 × 160, 100 ×
and thermal resistance of CPCMs is analyzed, the melting processes of
200, and 160 × 320 grid sizes with the model of Lx × Ly = 60 × 120
four typical layouts, including UP (λx = λy = 0), HLP (λx = 0.16, λy = 0),
mm2. Fig. 2 (b) shows the validation of time step independence with
VLP (λx = 0, λy = 0.16), CLP (λx = 0.16, λy = 0.16) are discussed. To
time steps of 0.05, 0.1, and 0.2 s. The maximum difference of liquid
further discuss and compare the melting rates of CPCMs, the dimen­
fraction is <2 %. Hence, 100 × 200 gird size and a time step of 0.1 s are
sionless parameters including average melting fraction (βm) and Fourier
chosen for the present study. For the cases with different Lx (Lx = 10, 30,
number (Fo) are defined. It is given as:
50, 70, 90 mm) in Section 4.3, meshes with higher grid density is applied
∫
to ensure the calculation accuracy, the grid sizes are 10 × 120, 30 × 120, 1 αt kPCM
βm = βdVPCM ; Fo = 2 ; α = (23)
50 × 120, 70 × 120, 90 × 120, respectively. VPCM L ρcp
VPCM

3.2. Model validation where VPCM is the total volume of PCMs, and α is the thermal diffusivity
of PCMs.
The main character of the melting in this model is melting nonuni­ To further explain the heat transfer mechanisms of CPCMs, the
formity caused by natural convection and metal foams with linear thermal resistance model method is adopted. Fig. 4 (a) shows the
porosity with the heating and cooling of HTF. Therefore, the model of schematic of melting process of CPCMs. The liquid PCMs absorb heat
PCMs composite with nonuniform metal foams and the HTF boundary from the heating wall and flow upward, release the heat at the melting
conditions should be verified. The present work compares the simula­ interface, and flow down to the bottom. The upper melting interface
tion results with the experiments carried out by Wang et al. [33]. This moves faster than the lower part due to the natural convection, leading
study porposed a type of CPCM in heat storage units composed with to a tilted melting interface. The heating wall temperature (Tw(y)) and
gradient porosity metal foam and compared it to the one with uniform horizontal distance between solid liquid interface and heatingwall (Lx,

Fig. 2. (a) Grid independence validation (b) Time step independence validation.

5
T. Si et al. Journal of Energy Storage 81 (2024) 110482

Fig. 3. Comparison of the temperature evolution of specified point between simulation and experiments carried by Wang et al. [30] (a) gradient porosity and (b)
uniform porosity.

Fig. 4. Schematic of (a) melting process (b) thermal resistance network in CPCMs.

f(y)) are adopted. Fig. 4 (b) shows the schematic of thermal resistance When Fo = 0.010 and 0.020, with more PCMs melting, natural con­
network, including convective thermal resistance of heating fluid (Rh, vection tends to affect the heat transfer process. Nuf continues to rise
conv) and convective thermal resistance of melted fluid region (Rf,conv). during this time period, indicating that the heat transfer mechanism is
To quantitatively analyze the thermal resistance and heat transfer dominated by heat convection gradually.
mechanism, the conductive thermal resistance of melted fluid region (Rf, The influence of linear porosity can also be seen in Fig. 5. In the case
cond) and average Nusselt number (Nuf ) are adopted. The above pa­ of HLP and CLP, the horizontal linear porosity can accelerate the melting
rameters are defined as follows: rate without changing the shape of the interface. It can be inferred that
linear porosity can improve the thermal performance of the LHTES
1
Rh,conv = (24) system mainly by changing the effective thermal conductivity of CPCMs.
hh Ly
From Fig. 5 (b), the horizontal linear porosity results in a higher Nuf . In
Rf,cond the early stage of melting, horizontal linear porosity accelerate the
Nuf = (25) melting process and the convection intensity becomes higher. In the late
Rf,conv
stage of melting, the horizontal linear porosity results in a lower porosity
∫Ly near the solid-liquid interface, which brings a faster flow rate near the
Rf,cond =
1
, Rf,conv =
1
=
Tw − Tm
(26) interface and higher Nuf .
hw Ly Qin From Fig. 5 (a), it is worth noting that in the case of VLP and CLP,
L∫
x,f (y)
0 dx
with the influence of vertical linear porosity, the effective thermal
1
dy
keff (x, y)
0 conductivity of the lower part of the enclosure is increased while that of
To investigate the melting behavior and thermal resistance under the the upper part is reduced, so that the melting uniformity can be signif­
different structures of linear porosity, Fig. 5 (a) compares the melting icantly improved by vertical linear porosity, making the melting inter­
interface when Fo = 0.005, 0.010, and 0.020 under the four different face almost parallel to the vertical wall in the later stage. From Fig. 5 (b),
layouts, and Fig. 5 (b) shows the evolution of Nuf of CPCMs. When Fo = the vertical linear porosity significantly results in a higher Nuf when
0.005, the melting interface develops nearly parallel to the heating wall 0.02 < Fo < 0.03. The vertical linear porosity makes the interface par­
allel to the vertical wall with larger area, which results in a higher
as shown in Fig. 5 (a). Nuf is close to 1 in Fig. 5 (b), indicating that heat
convection intensity in this stage. When the CPCMs melt completely (Fo
conduction plays a central role in the heat transfer inside the CPCMs.

6
T. Si et al. Journal of Energy Storage 81 (2024) 110482

⎛
∫ ∫
1 ⎝
E= LβρPCM dVPCM + cp,PCM (TPCM − Ta )ρPCM dVPCM
mPCM
VPCM VPCM
⎞ (27)
∫
⎟
+ cp,mf (T − Ta )ρmf dVmf ⎠
Vmf

E
P= (28)
t

where the VPCM and Vmf indicate the volume of PCM and metal foam
respectively. The first term on the right side of Eq. (27) indicates the
latent heat of PCMs while the other two terms indicate the sensible heat
of PCM and metal foam respectively.

4.2.1. Effects of horizontal linear porosity

Fig. 6 (a) depicts the solid-liquid interface of HLP in the case of λx =
0, 0.16, and 0.32 at different times of melting (Fo = 0.01 and 0.02). As
shown in Fig. 6 (a), it is worth noting that the solid-liquid interface is
getting tilted with the development of the melting process rather than
being parallel to the heating wall. It can be inferred that the layout of
HLP has little impact on this phenomenon of nonuniform melting caused
x y by natural convection. In the preliminary stage of melting (Fo = 0.01),
x y
with the increase of λx the solid-liquid interface moves faster while the
average liquid fraction of CPCMs is getting higher. The increase of the
x y
HLP brings about the increase of the effective thermal conductivity
x y
adjacent to the heating wall, hence the melting rate would be acceler­
ated. Later in the melting process (Fo = 0.02), the effective thermal
conductivity away from the heating wall is reduced by the layout of HLP,
hence, the melting rate appears to decrease at the later stage of melting
Nu

under the impact of HLP. The curve of the average liquid fraction in the
case of λx = 0, 0.08, 0.16, 0.24, and 0.32 while λy = 0 is depicted in Fig. 6
(b). The average liquid fraction also rises faster in the early stage and
slower in the late stage with the increase of λx, as shown in Fig. 6 (b).
To evaluate the heat storage capacity in different cases, the evolution
of the average heat flux (Qin) of the heating wall is shown in Fig. 7. In all
cases Qin reaches the maximum value in the beginning and decreases
quickly (Fo < 0.01). Later, the decline of Qin is getting slowed down
Fo gradually in the mid-stage of melting (0.01 < Fo < 0.02). In the last stage
of melting (0.02 < Fo < 0.03), Qin experiences a dramatic decline again
Fig. 5. (a) Melting interface of CPCMs with four typical layouts of porosity at due to the challenging for melting of remaining solid region. It is worth
Fo = 0.005, 0.010, and 0.020 (b) Evolution of Nuf of CPCMs with four typical noting that a turning point exists in Fig. 7 when Fo = 0.015. Before the
layouts of porosity. turning point, increasing λx results in a higher Qin while after the turning
point the opposite trend is achieved. This trend indicates that the
> 003), the Nuf begins to decrease due to the drop of temperature dif­ effective thermal conductivity improved by HLP has impact on the heat
ference in CPCMs. storage performance of LHTES systems.
It can be inferred that linear porosity can improve the thermal per­ To quantificationally evaluate the energy storage efficiency and en­
formance of the LHTES system mainly by affecting the effective thermal ergy storage capacity with the layout of HLP, the power density, and Fo
conductivity of CPCMs and the intensity of natural convection in when βm reaches 1 are defined as average power density (Pave) and Fotot
different regions of CPCMs during different melting stage. Therefore, it respectively. Fig. 8 shows Fotot by the left Y-axis and Pave by the right Y-
is necessary to explore the design method of linear porosity to optimize axis under λx = 0, 0.08, 0.16, 0.24, and 0.32. As the heat transfer per­
the thermal storage performance of LHTES systems, which will be dis­ formance is transformed by HLP, the total thermal resistance would
cussed in Sections 4.2 and 4.3. change depending on λx. From Fig. 8, with the increase of λx, the Fotot
first decreases and then increases while Pave first increases and then
4.2. Effect of linear porosity on thermal storage performance of LHTES decreases. In the case of λx = 0.16, the optimal Fotot is obtained, which
system decreases by 3.6 % compared with UP. The optimal Pave is obtained
when λx = 0.24, which increases by 6.1 % compared with UP. Consid­
From the view of engineering, energy storage density and power ering the strengthening effect in the early stage and weakening effect in
density are two primary targets to evaluate the thermal storage perfor­ the late stage overall, there exists an optimal value of λx to achieve the
mance of the LHTES system. In this section, energy storage density (E) best energy storage efficiency of LHTES systems.
and power density (P) are defined to indicate the amount of energy
stored and the energy storage rate per unit mass, which includes sensible 4.2.2. Effects of vertical linear porosity
heat and latent heat. The definitions of energy storage density and During the melting process, nonuniform melting appears in the case
power density are shown in the following equations: of UP and HLP, which causes an undesirable decrease in the melting
rate. Considering the metal foam with VLP can improve the melting

7
T. Si et al. Journal of Energy Storage 81 (2024) 110482

Fig. 6. (a) Solid-liquid interface (Fo = 0.01 and 0.02) and (b) evolution of average liquid fraction (βm) of CPCMs in the case of HLP under different λx.

Fig. 7. Evolution of the average heat flux (Qin) in the case of HLP under different λx.

uniformity, it may also improve the thermal storage performance of process (Fo = 0.02), the melting uniformity of PCMs is greatly improved
CPCMs. Fig. 9 (a) shows the solid-liquid interface of VLP in the case of λy in the case of VLP. Furthermore, from Fig. 9 (b) it can be found that VLP
= 0, 0.16, and 0.32 at different times of melting (Fo = 0.01 and 0.02). begins to enhance the melting rate in the later stage of melting.
The curve of the average liquid fraction in the case of λy = 0, 0.08, 0.16, The evolution of Qin in the cases of VLP is also depicted in Fig. 10.
0.24, and 0.32 while λx = 0 is depicted in Fig. 9 (b). Totally, the curves of Qin show a similar trend as the cases of HLP. In the
In the preliminary stage of melting (Fo = 0.01), the melting interface early stage (Fo < 0.01), Qin begins with a high value and drops rapidly.
keeps almost parallel to the heating surface. The VLP enhances the In the middle stage (0.02 < Fo < 0.03), the curves flatten out. In the last
thermal conductivity in the lower part of the enclosure and reduces that stage (0.02 < Fo < 0.03), the curve drop significantly again. As we can
in the upper part. With the increase of λx, the PCMs in the lower part see from Fig. 10, there also exists a turning point when Fo = 0.02 in these
tend to melt faster than the upper part, causing a change in the shape of VLP cases. Before the turning point, the Qin curves of all cases are almost
the melting interface as shown in Fig. 9 (a). However, the curves in Fig. 9 coincided with each other, while after the turning point the curves start
(b) almost overlap in the preliminary stage of melting, proving that VLP to look different. When λy = 0.24, the curve is at the highest position.
has little effect on the early-stage melting rate. Later in the melting From Fig. 9 (a) and (b) it can be inferred that the turning point appears

8
T. Si et al. Journal of Energy Storage 81 (2024) 110482

Fig. 8. Total Fourier number (Fotot) and average power density (Pave) under different λx.

Fig. 9. (a) Solid-liquid interface (Fo = 0.01 and 0.02) and (b) evolution of average liquid fraction (βm) of CPCMs in the case of VLP under different λy.

when the melting interface reaches the right wall, which is also the increases the melting speed in the later stage. However, Pave in the case
moment that the remain solid region starts to harm the melting rate. This of λy = 0.32 is lower than that of λy = 0.24, indicating that too high λy
trend indicates that the distribution of effective thermal conductivity may have a negative effect on energy storage. When λy = 0.24 the
depending on VLP can improve the heat storage performance of LHTES melting interface appears great uniformity when the melting process is
systems by improving the melting nonuniformity. almost complete, so continuing to increase the λy would make a remain
In the case of VLP, the average power density and total energy solid region in the upper right corner, which does harm to the melting
storage density are also discussed. Fig. 11 shows the total Fourier rate in the late stage. From Fig. 11, in the case of λy = 0.24, the optimal
number by the left Y-axis and average power density by the right Y-axis Fotot and Pave are obtained, which decreases by 14.6 % and increases by
under λy = 0, 0.08, 0.16, 0.24, and 0.32. The total thermal resistance 9.9 % respectively compared with UP. In conclusion, there exists an
cannot be reduced in VLP, which is different from the case of HLP. optimal value of λy to reduce the negative effect of melting nonunifor­
Nonetheless, the average power density Pave first increases and then mity and achieve a better energy storage efficiency of LHTES systems.
decreases with the increase of λy. This trend is mainly because the VLP
enhances the heat transfer of the lower part and weakens that of the 4.2.3. Effects of composite linear porosity
upper part. It dramatically improves the melting uniformity and The thermal enhancing effects of HLP and VLP are discussed in detail

9
T. Si et al. Journal of Energy Storage 81 (2024) 110482

Fig. 10. Evolution of the average heat flux (Qin) in the case of VLP under different λy.

Fig. 11. Total Fourier number (Fotot) and average power density (Pave) under different λy.

in the last two sections and both HLP and VLP can accelerate the average 500 s. The evolution of the average liquid fraction of UP, HLP, VLP, and
power density. Hence, it can be inferred that CLP is a promising way to CLP is depicted in Fig. 12 (b).
improve the thermal storage performance of the LHTES system. In this In the early stage of melting (Fo = 0.01), it is already analyzed that
section, the CLP (λx = λy = 0.16) is compared with the HLP (λx = 0.16, λy the HLP can accelerate the melting rate and VLP can change the incli­
= 0) and VLP (λx = 0, λy = 0.16). Fig. 12 (a) shows the solid-liquid in­ nation of the melting interface. Compared with VLP, the case of CLP has
terfaces of UP, HLP, VLP, and CLP at different melting times of 100 and the same inclination of melting interface and the interface of CLP moves

10
T. Si et al. Journal of Energy Storage 81 (2024) 110482

Fig. 12. (a) Solid-liquid interface of CPCM (Fo = 0.01 and 0.02) and (b) evolution of average liquid fraction (βm) of CPCMs in the case of UP, HLP, VLP, and CLP.

faster. Compared with HLP, the case of CLP which has a changed melting 4.2.2. Compared with ULP, CLP has almost the same Qin in the early
interface with PCMs in the lower part melts faster. From the average stage (Fo < 0.15) and higher Qin in the late stage (Fo > 0.15). Compared
liquid fraction curves shown in Fig. 12 (b), The curve of CLP is almost with VLP, CLP has higher Qin in the early stage and lower Qin in the late
coincided with that of HLP in the early stage, which is higher than that of stage. Compared with UP, CLP has higher Qin in most of the time while
UP and VLP. In the later stage of melting (Fo = 0.02), the solid-liquid lower Qin just for a short time before the melting process finishes. The
interface of CLP still moves faster than that of VLP, with the same distribution of effective thermal conductivity depending on CLP can
shape. Compared with HLP, the melting uniformity of PCM is greatly obtain the improvement features of HLP and VLP at the same time.
improved in the case of CLP. The liquid fraction curves in Fig. 12 (b) In Fig. 14, the total Fourier number and average power density in the
indicate that the CLP achieves the highest average liquid fraction of all case of UP, HLP, VLP, and CLP are also depicted. In the case of CLP, the
four cases, with optimal energy storage efficiency both in the early stage linear porosity exists both in the horizontal and vertical directions and it
and late stage. possesses the ability to enhance the effective thermal conductivity in the
The evolution of Qin in the cases of UP, HLP, VLP, and CLP is also early stage and improve the melting uniformity in the later stage.
depicted in Fig. 13. The Qin curve of CLP possesses the characteristics of Furthermore, it can reduce the remaining solid PCM in the steady stage
both HLP and VLP, which have been discussed in Section 4.2.1 and to the maximum extent. So it is reasonable that CLP can achieve better

Fig. 13. Evolution of the average heat flux (Qin) in the cases of UP, HLP, VLP, and CLP.

11
T. Si et al. Journal of Energy Storage 81 (2024) 110482

Fig. 14. Total Fourier number (Fotot) and average power density (Pave) in the case of UP, HLP, VLP, and CLP.

Fotot and Pave. So for UP, HLP, VLP, and CLP, the Fotot values are 0.0287, a series of linear porosity structures with different aspect ratios (Lx/Ly =
0.0276, 0.0260, and 0.0244 while Pave values are 1440.4, 1516.5, 1/12, 3/12, 5/12, 7/12, 9/12) are simulated and analyzed.
1545.8, and 1668.1 W, respectively. The Fotot and Pave of CLP are 15.0 % As shown in Fig. 15, the x coordinate represents λx while the y co­
lower and 15.8 % higher than those of UP. ordinate represents λy. The coordinate is limited in the triangle space
because the upper limit of εmax is set to be 0.96. Pave is chosen to be the
evaluation indicator in contour maps. In Fig. 15 (a), the enclosure length
4.3. Analysis of the design strategy of composite linear porosity with Lx = 0.05 m. The Pave first increases and then decreases both in the di­
different aspect ratios rection of the x and y axis and the optimal value of λx and λy can be
found. The split line that indicates λy = 0.16 is marked in Fig. 15 (a).
The effects of linear porosity on the thermal performance of LHTES Around the split line, 6 typical points are selected, which respectively
systems are discussed in the last section. However, from an engineering correspond to the 6 contours in Fig. 15 (b). Point 1 and 2 are above the
perspective, the dimension and aspect ratio of the enclosure would be split line, point 3 and 4 are on the split line, and point 5 and point 6 are
changed to adapt to specific applications. Also, the porosity value has its below the split line. Fig. 15 (b) shows the solid-liquid interfaces of cases
natural limitation which should be considered in the design of λx and λy. correspond to the 6 typical points at the late stage of melting (Fo = 0.14).
In this section, to explore a design strategy of composite linear porosity,

Fig. 15. (a) Contour maps of average power density (Pave) with Lx/Ly = 5/12. (b) Comparison of solid-liquid interfaces in different cases when Fo = 0.14.

12
T. Si et al. Journal of Energy Storage 81 (2024) 110482

In the cases of point 1 and 2 (λy = 0.24), the remaining solid region is on split line moves upward and the natural convection intensity becomes
the top of the enclosure for the overoptimization of the VLP. In the cases stronger. It can be inferred that the higher aspect ratio causes stronger
of point 3 and 4 (λy = 0.16), the solid-liquid interface is approximately natural convection in the enclosure, thus a higher λy is needed to balance
parallel to the right wall, indicating a good opitimization effect. In the the melting nonuniformity in vertical direction. When Lx/Ly is relatively
cases of point 5 and 6 (λy = 0.08), with the underoptimization of the high (Lx/Ly = 7/12, 9/12), the split lines disappear in the contour. With
VLP, the remaining solid region is still on the right bottom corner of the the development of the melting process, the natural convection intensity
enclosure for the strong natural convection. Generally, it can be inferred become so strong that the VLP structure cannot balance the melting
that in the cases of the points close to the split line, the opitimization nonuniformity even λy is at its maximum. It is worth noting that the
effect of VLP can balance the influence of natural convection, leading to contour maps of Lx/Ly = 7/12 and 9/12 are similar in shape, while these
a higher Pave value of the split line. two cases have similar natural convection intensity. This trend indicates
Furthermore, since the split line possesses the optimal value of Pave, that with the increase of Lx/Ly, the natural convection intensity first
it is critical to explore the intrinsic connection between the split line, increases and then becomes stable, meanwhile the shape of the contour
natural convection intensity and aspect ratios. Fig. 16 (a) demonstrates map of Pave also stabilizes. The opitimization effect of linear porosity is
the contour maps of Pave with different aspect ratios (Lx/Ly = 1/12, 3/ intrinsically related to the natural convection intensity.
12, 5/12, 7/12, 9/12). The flow velocity contours of these 5 enclosure Compare the 5 contour maps of Pave, it can be seen that the optimal λx
sizes at the end of melting in the case of UP are shown in Fig. 16 (b), is achieved in the similar position, which is independent from Lx. This
which indicate the intensity of natural convection in the enclosure. trend confirms that the aspect ratio and natural convection intensity
When Lx/Ly is relatively low (Lx/Ly = 1/12, 3/12, 5/12), the split line does not have much impact on the opitimization effect of HLP.
can be easily found in the contour maps. With the increase of Lx/Ly, the Comparing the simulation data of different aspect ratios, the maximum

Fig. 16. (a) Contour maps of average power density (Pave) with different aspect ratios (Lx/Ly = 1/12, 3/12, 5/12, 7/12, 9/12). (b) Comparison of flow velocity
contours with different aspect ratios (Lx/Ly = 1/12, 3/12, 5/12, 7/12, 9/12) in the case of UP at the end of melting.

13
T. Si et al. Journal of Energy Storage 81 (2024) 110482

Pave is 3349.5, 2466.7, 1899.1, 1521.7, 1331.6 W when Lx/Ly = 1/12, 3/ Science Foundation of China (Grant No. 51721004).
12, 5/12, 7/12, 9/12, respectively. Obviously the Pave is getting lower
with the increase of Lx, thus a lower aspect ratio is recommended in References
practical applications.
[1] T. Yang, W.P. King, N. Miljkovic, Phase change material-based thermal energy
storage, Cell Rep. Phys. Sci. 2 (2021) 100540.
5. Conclusions [2] F. Agyenim, N. Hewitt, P. Eames, M. Smyth, A review of materials, heat transfer
and phase change problem formulation for latent heat thermal energy storage
systems (LHTESS), Renew. Sust. Energ. Rev. 14 (2010) 615–628.
To explore the heat transfer mechanism and improve the thermal
[3] E. Oró, A. de Gracia, A. Castell, M. Farid, L. Cabeza, Review on phase change
energy storage performance of the LHTES system, an analysis CPCMs materials (PCMs) for cold thermal energy storage applications, Appl. Energy 99
with two-dimensional linear porosity metal foam is proposed. The effect (2012) 513–533.
of linear porosity structure, including UP, HLP, VLP and CLP on melting [4] Q. Wang, C. Zhao, Parametric investigations of using a PCM curtain for energy
efficient buildings, Energy Buildings 94 (2015) 33–42.
process and thermal storage performance are analyzed respectively. The [5] J. Mahdi, R. Pal Singh, H. Taqi Al-Najjar, E. Nsofor, Efficient thermal management
results support the following conclusions. of the photovoltaic/phase change material system with innovative exterior metal-
(1) The effect of linear porosity on melting behavior and thermal foam layer, Sol. Energy 216 (2021) 411–427.
[6] J. Weng, X. Yang, G. Zhang, D. Ouyang, M. Chen, J. Wang, Optimization of the
resistance network of CPCMs are analyzed. Average Nusselt number detailed factors in a phase-change-material module for battery thermal
(Nuf ) is adopted to study the changing rules of thermal conduction management, Int. J. Heat Mass Transf. 138 (2019) 126–134.
resistance and convection resistance. The results show that both HLP [7] D. Li, J. Wang, Y. Ding, H. Yao, Y. Huang, Dynamic thermal management for
industrial waste heat recovery based on phase change material thermal storage,
and VLP can enhance the thermal convection intensity and increase Nuf Appl. Energy 236 (2019) 1168–1182.
in different stages of melting. [8] Y. Tao, Y. He, A review of phase change material and performance enhancement
method for latent heat storage system, Renew. Sust. Energ. Rev. 93 (2018)
(2) When the aspect ratio is constant, the enhancement ratio of
245–259.
average power density (Pave) is first getting higher and then lower with [9] J. Mahdi, E. Nsofor, Solidification enhancement in a triplex-tube latent heat energy
the increase of λx in the case of HLP. The optimal Pave is obtained when storage system using nanoparticles-metal foam combination, Energy 126 (2017)
λx = 0.24, which increases by 6.1 % compared with UP. In the case of 501–512.
[10] B. Zalba, J. Marı́n, L. Cabeza, H. Mehling, Review on thermal energy storage with
VLP, Pave is first getting higher and then lower with the increase of λy. phase change: materials, heat transfer analysis and applications, Appl. Therm. Eng.
The optimal Pave is increased by 9.9 % when λy = 0.16 compared with 23 (2003) 251–283.
UP. In the case of CLP, the maximum Pave is achieved, which is 15.0 % [11] W. Su, J. Darkwa, G. Kokogiannakis, Review of solid–liquid phase change materials
and their encapsulation technologies, Renew. Sust. Energ. Rev. 48 (2015) 373–391.
higher than those of UP. [12] X. Yang, J. Yu, T. Xiao, Z. Hu, Y. He, Design and operating evaluation of a finned
(3) By analyzing the contour maps of Pave with different enclosure shell-and-tube thermal energy storage unit filled with metal foam, Appl. Energy
aspect ratios, the split lines which indicates the optimal λy can be found 261 (2020) 114385.
[13] J. Mahdi, E. Nsofor, Melting enhancement in triplex-tube latent heat energy
in the contour maps. With the increase of Lx/Ly, the natural convection storage system using nanoparticles-metal foam combination, Appl. Energy 191
intensity first increases and then becomes stable, while the split line (2017) 22–34.
moves upward and eventually disappears. However, the enclosure [14] W. Cui, X. Li, X. Li, L. Lu, T. Ma, Q. Wang, Combined effects of nanoparticles and
ultrasonic field on thermal energy storage performance of phase change materials
aspect ratio does not have much impact on the opitimization effect of
with metal foam, Appl. Energy 309 (2022) 118465.
HLP. For the general Pave gets lower with the increase of Lx/Ly, a lower [15] C. Zhao, Z. Wu, Heat transfer enhancement of high temperature thermal energy
aspect ratio is recommended in practical applications. storage using metal foams and expanded graphite, Sol. Energy Mat. Sol. C. 95
(2011) 636–643.
This paper investigated the heat transfer mechanism and thermal
[16] W. Chen, X. Liang, S. Wang, Y. Ding, X. Gao, Z. Zhang, SiO2 hydrophilic
performance of CPCMs with two-dimensional linear porosity metal modification of expanded graphite to fabricate form-stable ternary nitrate
foam. It is perspective to explore more non-uniform structures or composite room temperature phase change material for thermal energy storage,
complicated pore structures in CPCMs in the future work. Efficient Chem. Eng. J. 413 (2021) 127549.
[17] X. Li, W. Cui, T. Simon, T. Ma, T. Cui, Q. Wang, Pore-scale analysis on selection of
optimization methods like topology optimization can be applied to composite phase change materials for photovoltaic thermal management, Appl.
further improve the heat transfer efficiency and store more energy. Energy 302 (2021) 117558.
[18] X. Xiao, P. Zhang, M. Li, Effective thermal conductivity of open-cell metal foams
impregnated with pure paraffin for latent heat storage, Int. J. Therm. Sci. 81
CRediT authorship contribution statement (2014) 94–105.
[19] Z. Wang, H. Zhang, B. Dou, G. Zhang, W. Wu, X. Zhou, Effect of copper metal foam
Tianyu Si: Writing – original draft, Methodology, Data curation, proportion on heat transfer enhancement in the melting process of phase change
materials, Appl. Therm. Eng. 201 (2022) 117778.
Conceptualization. Wei Cui: Writing – review & editing, Methodology. [20] A. Sari, C. Alkan, C. Bilgin, Micro/nano encapsulation of some paraffin eutectic
Ting Ma: Writing – review & editing. Lin Lu: Writing – review & editing. mixtures with poly(methyl methacrylate) shell: preparation, characterization and
Qiuwang Wang: Writing – review & editing, Supervision, Funding latent heat thermal energy storage properties, Appl. Energy 136 (2014) 217–227.
[21] Y. Tian, C. Zhao, A numerical investigation of heat transfer in phase change
acquisition. materials (PCMs) embedded in porous metals, Energy 36 (2011) 5539–5546.
[22] R. Ferfera, B. Madani, Thermal characterization of a heat exchanger equipped with
a combined material of phase change material and metallic foams, Int. J. Heat Mass
Declaration of competing interest Transf. 148 (2020) 119162.
[23] K. Lafdi, O. Mesalhy, S. Shaikh, Experimental study on the influence of foam
porosity and pore size on the melting of phase change materials, J. Appl. Phys. 102
The authors declare that they have no known competing financial (2007) 083549.
interests or personal relationships that could have appeared to influence [24] W. Li, Z. Qu, Y. He, W. Tao, Experimental and numerical studies on melting phase
the work reported in this paper. change heat transfer in open-cell metallic foams filled with paraffin, Appl. Therm.
Eng. 37 (2012) 1–9.
[25] Y. Tao, Y. You, Y. He, Lattice Boltzmann simulation on phase change heat transfer
Data availability in metal foams/paraffin composite phase change material, Appl. Therm. Eng. 93
(2016) 476–485.
[26] X. Li, J. Duan, T. Simon, T. Ma, T. Cui, Q. Wang, Nonuniform metal foam design
Data will be made available on request.
and pore-scale analysis of a tilted composite phase change material system for
photovoltaics thermal management, Appl. Energy 298 (2021) 117203.
Acknowledgement [27] M. Fteiti, M. Ghalambaz, M. Sheremet, et al., The impact of random porosity
distribution on the composite metal foam-phase change heat transfer for thermal
energy storage, J. Energy Storage 60 (2023) 106586.
This work was supported by the State Key Program of National [28] D. Zhang, C. Zhu, B. Huang, X. Duan, L. Gong, M. Xu, Thermal behaviors and
Natural Science Foundation of China (Grant No. 52130609), and the performance of phase change materials embedded in sparse porous skeleton
structure for thermal energy storage, J. Energy Storage 62 (2023) 106849.
Foundation for Innovative Research Groups of the National Natural

14
T. Si et al. Journal of Energy Storage 81 (2024) 110482

[29] Z. Zheng, C. Yang, Y. Xu, X. Cai, Effect of metal foam with two-dimensional [33] Z. Wang, J. Wu, D. Lei, H. Liu, J. Li, Z. Wu, Experimental study on latent thermal
porosity gradient on melting behavior in a rectangular cavity, Renew. Energy 172 energy storage system with gradient porosity copper foam for mid-temperature
(2021) 802–815. solar energy application, Appl. Energy 261 (2020) 114472.
[30] H. Ahmadi, N. Variji, A. Kaabinejadian, M. Moghimi, M. Siavashi, Optimal design [34] W. Cui, T. Si, X. Li, X. Li, L. Lu, T. Ma, Q. Wang, Heat transfer analysis of phase
and sensitivity analysis of energy storage for concentrated solar power plants using change material composited with metal foam-fin hybrid structure in inclination
phase change material by gradient metal foams, J. Energy Storage 35 (2021) container by numerical simulation and artificial neural network, Energy Rep. 8
102233. (2022) 10203–10218.
[31] J. Kong, H. Zuo, K. Zeng, Y. Lu, H. Xu, X. Zhang, H. Yang, H. Chen, Investigation of [35] H. Shmueli, G. Ziskind, R. Letan, Melting in a vertical cylindrical tube: numerical
optimization methods for metal foam with two-dimensional porosity gradient in investigation and comparison with experiments, Int. J. Heat Mass Transf. 53 (2010)
shell-and-tube latent heat storage, J. Energy Storage 63 (2023) 107004. 4082–4091.
[32] J. Yang, L. Yang, C. Xu, X. Du, Numerical analysis on thermal behavior of [36] J. Fourie, J. Du Plessis, Pressure drop modelling in cellular metallic foams, Chem.
solid–liquid phase change within copper foam with varying porosity, Int. J. Heat Eng. Sci. 57 (2002) 2781–2789.
Mass Transf. 84 (2015) 1008–1018. [37] K. Boomsma, K. Poulikakos, On the effective thermal conductivity of a three-
dimensionally structured fluid-saturated metal foam, Int. J. Heat Mass Transf. 44
(2001) 827–836.

15