Phase-Change Process Inside a Small-Radii Cylinder Subjected to Cyclic Convective Boundary Conditions: A Numerical Study

Phase-change materials (PCMs) are widely used in energy storage applications because of their high latent heat [1–3]. The temperature and density of PCMs remain nearly constant during the phase change [4], which makes them suitable for a variety of thermal management applications including heat exchangers [1,5], electronic cooling [6,7], building applications [8], solar energy storage [9], and spacecraft thermal control [10,11]. In these applications, PCMs are often encapsulated in a variety of enclosures (e.g., cylinders, spheres, annular cavities, etc.) to reduce their reactivity with the environment, control the volume change, and provide a supporting structure [12].

Phase-change processes in the encapsulated geometries are inherently complex and transient, since the propagation of the solid–liquid interfaces depend on the rate and direction of heat transfer. In larger geometries, the flow induced by the buoyancy can be important. Apart from experimental and analytical approaches, different numerical methods based on enthalpy formulation models, effective heat capacity models, and temperature-transforming models have been used in the literature to simulate the complex phase-change process [13–19]. Analytical solutions are obtained by tracking the two-phase interface and solving each phase independently while coupling the heat transfer at the interface. However, the enthalpy formulation [20,21] or the effective heat capacity [22] models utilize a single differential equation that governs both phases, and hence there is no need to track the interface. In the enthalpy formulation model, the total enthalpy of the PCM is decomposed into latent and sensible heat to obtain a differential equation that includes the melt fraction. On the other hand, in the effective capacity model, a large heat capacity is used in the phase-change temperature range to account for the latent heat effects [17,22,23].

Many experimental, analytical, and numerical studies on the melting and solidification of the PCM in encapsulated geometries have been reported in the literature by imposing constant temperature or heat flux at the encapsulation boundary and for a variety of encapsulation geometries such as spheres [24–26], rectangular containers [27–30], and horizontal [31–33] and vertical cylinders [34–43]. These studies provide a good understanding of the heat transfer and thermal storage characteristics of the encapsulated PCMs (EPCMs) under fixed temperature or heat flux boundary conditions.

In many practical applications such as PCM heat exchanger design [1], buildings exposed to solar radiation with embedded PCMs for energy management, and waste heat thermal storage, the encapsulated PCM is exposed to cyclic thermal conditions. Due to the transient nature of the melting front and the low thermal conductivity associated with the PCMs, the heat transfer rate could vary significantly under various cycling frequencies. Moreover, multiple phase bands inside the encapsulated domain may form [44], which are not studied extensively in the literature. Among the limited number of studies that considered a cyclic boundary condition, Ho and Chu [45] imposed a sinusoidal surface temperature to study a natural convection-dominated melting process. Casano and Piva [46] performed numerical and experimental studies of the melting and freezing process in a plane PCM slab under cyclic temperature boundary conditions and studied the effects of Fourier and Stefan numbers on the energy storage capacity. They observed that the mean value of the energy storage in the system decreases at higher values of Fourier and Stefan numbers. Mazzeo et al. [47] solved a Stefan problem exposed to a periodic boundary condition analytically by considering a single melting interface. More recently, Kant et al. [48] studied the melting behavior of PCMs under sinusoidally varying heat flux conditions and reported that the time of melting reduces under a variable heat flux when compared to an equivalent constant heat flux.

This paper focuses on the thermal analysis of PCM encapsulated in a long cylindrical polymer tube of a small diameter (larger aspect ratio) under a cyclic convective boundary condition. The effects of natural convection are neglected because of the small tube diameters [49]. The motivation for the research stems from the need to accurately capture the transient movement of the phase-change interface and heat transfer rate to and from an encapsulated PCM exposed to cyclic heating and cooling [1] over a broad range of heating and cooling times that feature multiple phase-change interfaces. The paper employs a numerical approach and characterizes the effect of PCM thermophysical properties, tube radius, ambient air temperature, convective heat transfer coefficient, and the duration of heating and cooling on the phase-change behavior and thermal performance of the encapsulated PCM. Results are presented using dimensionless parameters (i.e., Fourier, Stefan, and Biot numbers) to provide a general framework for utilizing PCM in applications where the cyclic heating and cooling is required or beneficial. For example, Karmakar et al. [1] proposed an EPCM heat exchanger designed to cool the water flow exiting a steam condenser in a steam power plant. In the EPCM heat exchanger, warm water from the condenser flows over an array of EPCM tubes and melts the PCM. This is followed by a period where ambient air flows over the tubes to refreeze the PCM and make it available for the next period of melting under the warm water flow. Thus, the cylindrical EPCM tubes were periodically subjected to heating and cooling. In addition to cooling in power plants [1], there is a potential application of similar EPCM heat exchangers in automotive and building heating and cooling applications, data center cooling, and a variety of other industrial applications.

where h is the specific enthalpy, and T is the temperature.

with the convergence threshold of $ϵ=10−6$⁠. Accuracy of the numerical simulations depends on the convergence threshold, especially for steep enthalpy curves in the melting temperature range as the error in temperature propagates into the error in enthalpy with a factor equal to the rate of change of enthalpy with temperature (i.e., $eh=eT×dh/dT$⁠). Hence, it is imperative to use lower threshold levels for the steep melting curves. To start the subiteration process, the temperature at the previous time-step estimates the initial temperature of $T0t+Δt$⁠, which results in the explicit estimation of heat transfer to each cell. It is important to note that the above procedure requires the knowledge of a functional relationship between the enthalpy and the temperature of the PCM material.

where $Ts=Tm−0.01$ and $Tl=Tm+0.01$⁠. Numerical simulations involving Eqs. (11) and (16) are conducted for a one-dimensional axisymmetric cylindrical domain of radius $1×10−3$ m, where a uniformly distributed grid of 500 cells and a time-step of $Δt=1×10−5$ s are employed. Figure 1 shows the analytical solution and numerical predictions of the location of the solid–liquid interface as a function of time. On average, the predictions deviate by 0.7% from the analytical solution.

Numerical simulations are performed for a cylindrical domain of radius of 5 m for which Eqs. (11) and (19) are solved. The numerical domain is discretized uniformly in the radial direction with 4000 cells, and a time-step of $Δt=1s$ is utilized. Figure 2 shows the analytical and predicted liquid and solid fronts as a function of time. Current numerical predictions agree well with the analytical solution with an average deviation of 2.5% which validates the accuracy of the numerical approach.

This section presents the ability of our solver to predict the cyclic melting and freezing process of a phase-change material encapsulated in a small diameter vertical cylinder. The experimental test setup by Kapilow et al. [49] was designed to evaluate the melting and freezing characteristics of a commercial PCM (RT35HC) [61] encapsulated inside a cylindrical high-density polyethylene (HDPE) tube subjected to convective boundary condition. They reported the measured PCM volume phase fraction and the overall heat transfer coefficient of the entire system over time. Figure 3 shows the schematic of the encapsulated PCM cross section where the air flows from left to right. The inner and outer tube diameters are $2.78×10−3$ m and $3.18×10−3$ m, respectively. Based on the outer tube diameter as the length scale, the Reynolds number of the air flowing over the tube is 1797. The reported Rayleigh number of the experimental tests is very low (50–1200) due to the small tube diameter and thus suggests insignificant effects of natural convection [49]. Hence, the natural convection effects are also neglected in the current study.

where a = 0.45775, b = 34.86900, c = 1.213776, and d = 0.454482. The fitted $hlat, hsen$⁠, and $htot$ are shown with lines in Fig. 4. It should be noted that the reported cooling and heating enthalpy curve in the manufacturer's data sheet slightly differs [61], and hence the averaged enthalpy curve is considered here. Utilizing the reported enthalpy curve of the RT35HC eliminates any assumptions regarding the enthalpy function, melting temperature range, and potentially leads to more accurate numerical simulations.

where $T$ is the half cycle. In Kapilow's experiment [49], $T$ is large enough that allows PCM to fully melt and freeze, and the experiment is only conducted over one cycle of melting and freezing. Equation (21), however, provides a general cyclic square-wave function which alternates the air temperature between $Tcold$ and $Thot$⁠.

where $Rair$ is the thermal resistance of the air, R_e is the thermal resistance of the HDPE tube, r_i and r_o are the inner and outer radius of the HDPE tube, respectively, and h_air is the heat transfer coefficient around the outer periphery of the tube.

To validate the numerical solver with the experimental data, we perform simulations at a low frequency to allow full melting and solidification. The average liquid phase fraction and the predicted overall heat transfer are compared with the measurements of Kapilow et al. [49]. Figure 5 shows the average liquid melt fraction γ of the PCM material with a cyclic convective boundary condition with $T=400 s$⁠. To ensure a grid-independent solution, we perform numerical simulations with three different grids of 10 × 15, 20 × 30, and 40 × 60 grids in radial and azimuthal directions. The dependency of the solution on the time-step is also evaluated by conducting the simulations with the highest grid resolution with $Δt=$ 5, 1, and 0.1 s, which corresponds to $Δt/tmelting$ of $3.37×10−2, 6.76×10−3$⁠, and $6.76×10−4$⁠, respectively (Fig. 6). The computed average liquid fraction over time is grid-independent for all evaluated grids (Fig. 5) and does not change for $Δt$ less than 1 s. Numerical predictions agree well with the measurements during the melting period (⁠$200 s<t<400 s$⁠), but the rate of freezing (⁠$400 s<t<600 s$⁠) is slightly overpredicted. As the numerical procedures are validated with the analytical solutions, this discrepancy may be associated with:

1. The difference between the RT35HC enthalpy curve during melting and freezing [61].

2. Uncertainty in the thermophysical characteristics of RT35HC (see different values reported in Refs. [61–63]).

3. Variation in the air velocity along the axis of the vertical cylinder in the experiments leads to a lower overall heat transfer coefficient.

This emphasizes that attention needs to be made in characterizing the thermophysical properties of the PCM materials.

Figure 7 shows the heat transfer rate predictions as a function of time during the melting and freezing process. The integration of the rate of change of the latent enthalpy over the numerical domain determines the latent portion of the total heat transfer rate. Predictions of the latent portion of the heat transfer rate agree well with the measurements [49]. The PCM temperature rises from the initial temperature after the air temperature increases at $t=200 s$⁠. The phase change only begins after the PCM temperature in the outer portion of the PCM reaches the melting temperature range, and hence the $qlat$ increases sharply for $t>200 s$ and reaches a maximum value which corresponds to the lowest possible resistance of the PCM to the heat transfer. During the melting process, the melt front progresses inward, and hence PCM exerts an extra resistance to the heat transfer, which is associated with the decrease in the latent heat transfer rate as observed in Ref. [49]. The freezing process behaves similarly with the only difference that the freezing rate is faster than the melting due to the higher thermal conductivity of the PCM in its solid phase and more substantial temperature difference with the ambient (see Table 3).

Simulations are performed for the same range of Ste and Bi number and for $0.1≤ζ≤1$ and continue until the PCM shows a cyclic behavior. Due to the uniform heat transfer coefficient assumption around the circular tube, all simulations are conducted using a one-dimensional grid with 40 grid points uniformly distributed from the center to the outer radius of the PCM using $Δt̃=0.01$⁠. The selected $Δt̃$ provides $Δt̃/T̃ref<5.1×10−3$ for all simulations.

with a mean and maximum deviation of 2.3% and 3.9% from the predictions, respectively. Figure 9 compares the predicted $T̃ref$ using numerical simulations with the ones calculated using the proposed correlation Eq. (35) at selected Biot numbers. The applicable range of the above correlation is $0.01≤Ste≤1$ and $0.1≤Bi≤10$⁠. Although the effect of the Ste and Bi on the $T̃ref$ is clear, it might be also necessary to understand the effect of each individual physical quantity (including the thermophysical properties of the PCM, and tube radius) on the actual time required to complete the heating or cooling process $Tref$⁠. The Appendix provides more insights on this matter. However, to simplify the interpretation of the Ste and Bi on the PCM behavior, in the context of the following analysis and without loss of generality, we consider constant thermophysical properties and a fixed tube radius. Hence, $Tref∝T̃ref$ and the variation of the Ste and Bi numbers are limited to the changes in the values of temperature difference (⁠$Thot$ − $Tcold$⁠) and the overall heat transfer coefficient U, respectively.

To understand the effect of Ste, Fig. 10 shows the temporal variation of the melt fraction and the dimensionless heat transfer rate of the system as a function of Fourier number $t̃$ for Ste values of 0.01, 0.1, and 1 with a fixed Biot number of one. Under a fixed Bi number, the corresponding reference Fourier number $T̃ref$ for these three cases are 155, 20, and 6, respectively. This process consists of three phases: an initial sensible heat transfer where the value of γ remains unchanged, followed by a combined latent and sensible heat transfer where the value of γ changes from 0 to 1, and then again a pure sensible heat transfer toward the end with a constant value of γ. Varying Ste only affects the second phase of heat transfer as when the value of γ changes over time, and thus the initial and final sensible heat transfer phases (i.e., $∂γ/∂t̃=0$⁠) are unaffected by Ste (see Eq. (27)). By definition, the Ste number represents the ratio of sensible heat to the latent heat. Consequently, at very low values of Ste that corresponds to a small temperature difference, the PCM melts at a very low heat transfer rate (see $q̃$ for Ste = 0.01 in Fig. 10), and thus the duration of PCM melting is the highest. Figure 11 shows this behavior where the surface temperature takes the longest (i.e., larger Fourier number) to reach the ambient condition. On the other hand, when the Ste number increases to 1, the rate of heat transfer in the PCM system is the highest (see $q̃$ for Ste = 1 in Fig. 10) due to the highest temperature difference, and thus the duration of melting turns out to be the least. The plot of the surface temperature at Ste = 1 in Fig. 11 also confirms this behavior as it reaches the ambient temperature at the fastest rate.

Figure 12 shows the effect of the Bi number on the temporal variation of melt fraction and the heat transfer rate of the system. Under a fixed Ste number, the corresponding reference Fourier number for these three cases are $T̃ref$ = 141, 20, and 7, respectively. Similar to the previous scenario, the entire process of heat transfer in the system consists of two purely sensible heat transfer portions and a combined latent and sensible heat transfer. The rise in the Bi number corresponds to an increase in the values of the convective heat transfer coefficient, and as such, the reference Fourier number decreases due to the increased heat transfer from the ambient to the PCM system. At a small value of the Bi number, the convective resistance dominates the resistance due to conduction, and hence the heat transfer rate is not significantly affected by the variation of the conductive resistance inside the tube during the phase-change process and remains almost constant (see Fig. 12)). The heat transfer rate is also the lowest at this Bi number due to the small heat transfer coefficient. The surface temperature for Bi = 0.1 (shown in Fig. 13) remains close to the PCM melt temperature over most of the period. When the Bi number rises to 1, the heat transfer rate decreases almost linearly with time as the convective resistance is no longer dominant, and hence the heat transfer rate is affected by the variations in the conductive thermal resistance inside the tube. Due to the higher Bi number, the melting rate also increases. Finally, at Bi = 10, the conductive resistance becomes more dominant such that the heat transfer rate in the system shows the strongest dependence on the conductive resistance of the PCM during the phase-change process, and hence, decreases very rapidly with the progression of the phase-change interface toward the center of the tube. The rate of melting is the highest under this condition as the value of the external heat transfer coefficient is the highest. The surface temperature at Bi = 10 approaches the ambient temperature very quickly compared to the previous conditions.

The discussion below presents the effect of the cyclic exposure of the PCM cylinder to hot and cold ambient temperature under various fractions of $T̃ref$⁠. For a fixed Ste = 0.1 and Bi = 1, Fig. 14 shows the variation of melt fraction as a function of Fourier number for ζ = 0.1, 0.5, and 1. This plot corresponds to the Fourier number interval of $2T̃ref$ after each of the cases has reached a cyclic state; ζ = 0.1 has 10 cycles, ζ = 0.5 has 2 cycles, and ζ = 1 has 1 cycle of alternating hot and cold ambient conditions. At ζ = 1, the heating and cooling process completes in each half period. However, reducing the value of ζ to 0.5 prohibits the completion of the heating and cooling process, and the melt fraction approximately varies between 0.1 and 0.9. For the lowest value of ζ, the melt fraction varies between 0.4 and 0.6 as the heating and cooling periods shorten. It is important to note that at the cyclic equilibrium condition, the melt fraction swings around 0.5 as both the heat transfer coefficient, and the difference between the melting temperature and the ambient temperature are identical for these heating and cooling processes. It is clear that having different values for these parameters during heating and melting could alter the cyclic equilibrium behavior of the PCM.

Figure 15 shows the normalized heat transfer rate distribution for the same Fourier number intervals as in Fig. 14. The time-averaged heat transfer rate for ζ = 0.1 is around 0.95, whereas it decreases to values of 0.81 and 0.57 when ζ is increased to 0.5 and 1, respectively. Thus, under the current conditions, exposure of the PCM to alternate heating or cooling for a shorter duration than $T̃ref$ is more effective in enhancing heat transfer. This interesting behavior can be explained in terms of the surface temperature at $r̃$ = 1 for various ζ values as shown in Fig. 16. At ζ = 0.1, the surface temperature oscillates closely around the PCM phase-change temperature (θ = 0), and therefore the rate of heat transfer is greater due to the higher temperature difference between the PCM outer surface and the ambient conditions. When the value of ζ increases to 0.5, the temporal variation of the PCM surface temperature moves away from the PCM melt temperature toward the ambient temperature, which decreases the effective temperature difference and the time-averaged heat transfer rate. At the highest value of ζ = 1, the PCM surface temperature reaches the ambient temperature, which further reduces the heat transfer rate, and thus the average heat transfer rate is lower.

As discussed earlier, the cyclic variation of the ambient temperature with a period lower than $T̃ref$ could potentially generate multiple phase-change interfaces. The outer portion of PCM is always exposed to this variation in the boundary condition, and hence the melting or freezing fronts are generated at this outer location. Here, we define the phase-change interface as any location in which the sign of θ changes (since θ = 0 represents the melting temperature). For significantly longer periods, a phase-change front initiates at $r̃=1$ and progresses toward the core of the tube (⁠$r̃=0$⁠) such that eventually, the whole PCM material is either liquid or solid. Figure 17 shows this behavior. When ζ is reduced to 0.5, a new phase-change front will be initiated at $r̃=1$ before the former front reaches the core (see Fig. 18). Hence, the case of $ζ=0.5$ shows a maximum of two phase-change fronts. The predictions at the ζ value of 0.1 (Fig. 19) shows a maximum of three simultaneous phase-change fronts. It is interesting to note that the speed at which these phase fronts move toward the core is also affected by the value of ζ and increases for a lower values of ζ.

In this study, we numerically investigate the melting or solidification behavior of phase-change materials encapsulated in a small-radii cylinder. We further study the case of the encapsulated PCM undergoing a cyclic heating and cooling process through a convective boundary condition with variable ambient temperature (square-wave) at the cylinder surface. A new correlation is proposed to estimate the reference Fourier number that corresponds to complete melting or freezing as a function of Ste and Bi. Higher Bi and Ste accelerates the freezing/melting process which is completed at smaller Fourier numbers. Studying the effects of Bi and Ste numbers on the rate of heat transfer and surface temperature helps us understand the impact of these parameters on the reference Fourier number.

The effects of the time period of cyclic heating and cooling on the time-averaged heat transfer rate are studied. Shorter periods of cyclic heating enhance the time-averaged heat transfer rate due to lower surface temperature variations. The tube consists of multiple phase-change fronts simultaneously; for lower ζ values (⁠$ζ=T̃/T̃ref$⁠), the corresponding phase-change fronts move faster toward the center of the tube.

The proposed correlation and the effect of Bi and Ste numbers on overall heat transfer and melting time can be used for applications that utilize small-radii cylindrical encapsulated PCM exposed to cyclic convective boundary conditions [1]. The significant changes observed in the thermal behavior of the encapsulated PCM under different cyclic conditions indicate the importance of conducting such studies for PCMs encapsulated in other geometries such as spheres or rectangular enclosures.

This work was supported by a grant from the Advanced Research Project Agency-Energy (DE-AR0000572) under the Advanced Research in Dry cooling (ARID) program with the Electric Power Research Institute (EPRI) as the prime contractor. This financial support is gratefully acknowledged.

• Advanced Research Project Agency-Energy (DE-AR0000572; Funder ID: 10.13039/100000015).

c =

specific heat (J kg⁻¹ K⁻¹)

D =

cylinder diameter (m)

e =

error

h =

specific enthalpy (J kg⁻¹)

H =

enthalpy (J)

k =

thermal conductivity (W m⁻¹ K⁻¹)

L =

specific enthalpy of fusion (J kg⁻¹)

m =

mass (kg)

q =

heat transfer rate (W)

Q =

heat transfer per unit mass (J kg⁻¹)

$Q$ =

heat (J)

r =

radial distance from cylinder center (m)

R =

thermal resistance (W⁻¹ m² K)

$R$ =

melt front location (m)

t =

time (s)

T =

temperature (K)

$T$ =

time period of the half cycle of melting/freezing (s)

U =

heat transfer coefficient (W⁻¹ m² K)

V =

volume (m³)

$V$ =

cell volume (m³)

$Δt$ =

time-step (s)

Greek Symbols

α =

thermal diffusivity (m² s⁻¹)

γ =

liquid volume fraction

ϵ =

convergence criteria

ζ =

fraction of the dimensionless half cycle time period of melting/freezing

θ =

dimensionless temperature

λ =

empirical constant

ν =

kinematic viscosity (m² s⁻¹)

ρ =

density (kg m⁻³)

Superscripts and Subscripts

cold =

lower limit of the ambient temperature

e =

HDPE tube

hot =

upper limit of the ambient temperature

i =

inner

l =

liquid

lat =

latent heat

m =

melting

n =

subiteration number

o =

outer

ref =

reference

s =

solid

su =

eutectic

$∞$ =

ambient

Dimensionless Groups

Bi =

Biot number

$q̃$ =

dimensionless heat transfer rate

$r̃$ =

dimensionless radial distance from cylinder center

Ste =

Stefan number

$t̃$ =

Fourier number, tα/L²

$T̃$ =

dimensionless time period of the half cycle of melting/freezing (s)

Hence, for a constant Stefan number, Eq. (35) indicates that the increase in Biot number monotonically decreases the reference Fourier number $T̃ref$⁠. The increase in the Biot number can be due to higher values of h or R or a lower value of k. Equation (A2), however, indicates that the reference time $Tref$ increases for larger length scales or lower values of the thermal conductivity, but decreases for higher heat transfer coefficients.

At a constant Biot number, $T̃ref$ monotonically decreases with Stefan number (see Eq. (35)) which can be due to higher $ΔT$ or thermal capacity, or lower latent heat of fusion. However, although lowering L or increasing $ΔT$ still decreases $Tref$⁠, an increase in the thermal capacity (i.e., higher Stefan number) actually increases the reference time. Furthermore, Eqs. (A2) and (A4) show that for a fixed Bi and Ste number, increasing the density of the PCM increases the time required to complete the heating or cooling process.