Open AccessArticle

Thermophysical Characterization of Paraffin Wax Based on Mass-Accommodation Methods Applied to a Cylindrical Thermal Energy-Storage Unit

Valter Silva-Nava

Ernesto M. Hernández-Cooper

Jesús Enrique Chong-Quero

and

José A. Otero

Tecnologico de Monterrey, Escuela de Ingeniería y Ciencias, Carr. al Lago de Guadalupe Km. 3.5, Atizapán de Zaragoza 52926, Estado de Mexico, Mexico

Author to whom correspondence should be addressed.

Molecules 2022, 27(4), 1189; https://doi.org/10.3390/molecules27041189

Submission received: 8 January 2022 / Revised: 3 February 2022 / Accepted: 5 February 2022 / Published: 10 February 2022

(This article belongs to the Special Issue Phase Change Materials 2.0)

Download

Browse Figures

Versions Notes

Abstract

Two mass-accommodation methods are proposed to describe the melting of paraffin wax used as a phase-change material in a centrally heated annular region. The two methods are presented as models where volume changes produced during the phase transition are incorporated through total mass conservation. The mass of the phase-change material is imposed as a constant, which brings an additional equation of motion. Volume changes in a cylindrical unit are pictured in two different ways. On the one hand, volume changes in the radial direction are proposed through an equation of motion where the outer radius of the cylindrical unit is promoted as a dynamical variable of motion. On the other hand, volume changes along the axial symmetry axis of the cylindrical unit are proposed through an equation of motion, where the excess volume of liquid constitutes the dynamical variable. The energy–mass balance at the liquid–solid interface is obtained according to each method of conceiving volume changes. The resulting energy–mass balance at the interface constitutes an equation of motion for the radius of the region delimited by the liquid–solid interface. Subtle differences are found between the equations of motion for the interface. The differences are consistent with mass conservation and local mass balance at the interface. Stationary states for volume changes and the radius of the region delimited by the liquid–solid interface are obtained for each mass-accommodation method. We show that the relationship between these steady states is proportional to the relationship between liquid and solid densities when the system is close to the high melting regime. Experimental tests are performed in a vertical annular region occupied by a paraffin wax. The boundary conditions used in the experimental tests produce a thin liquid layer during a melting process. The experimental results are used to characterize the phase-change material through the proposed models in this work. Finally, the thermodynamic properties of the paraffin wax are estimated by minimizing the quadratic error between the temperature readings within the phase-change material and the temperature field predicted by the proposed model.

Keywords:

mass-accommodation methods; phase-change materials; cylindrical thermal energy-storage unit

1. Introduction

The thermal energy density of systems based on latent heat-storage units can be increased by using the latent heat of materials as an additional form of energy storage. The energy density of thermal energy-storage (TES) units that are only based on sensible heat is significantly lower than energy-density values achieved on latent heat thermal energy-storage (LHTES) units. The subject of solar energy harvesting in concentrating solar power plants for thermoelectric generation [1,2,3] and for domestic water heating systems [4] has become an appealing subject in research studies that are focused on the potential applications of these systems. The subject of heat storage presents one of several alternative applications that are aimed at reducing fossil fuel consumption, which has been recently representing a concerning issue due to the worldwide increase in greenhouse gas emissions. Experimental studies have been performed on several types of LHTES units to enhance heat-transfer rates and heat-storage properties. The transfer rate of thermal energy between a heat-transfer fluid (HTF) and the phase-change material (PCM) used to store energy represents a crucial parameter in TES units. The intermittence of solar energy produced by the solar irradiance oscillations during a 24 h period presents a challenging problem on backup systems based on LHTES units, due to the low thermal conductivity of PCMs used in these types of applications.

Domestic heat water applications based on PCMs have also been extensively studied. Paraffin and salt hydrates have been widely used in these kinds of applications due to their low cost and relatively high energy-storage capacity in narrow operating temperature ranges [5]. Paraffin has been used as a PCM during melting and solidification experiments in a tilted annular region, where liquid water acting as a HTF is flowing through an inner cavity [5]. The outer surface of the wall was tilted to improve the rate of heat transfer during a melting (charging) and solidification (discharge) process. Thermal energy-storage systems that combine heat-exchange strategies through copper rods and graphite particles to increase the heat-transfer rates in a cylindrical unit, have also been studied [6]. Additionally, 2D models have been used to analyse the thermal performance of truncated conical systems and cylindrical heat-storage units with fins to enhance heat transfer rates [7]. The relationship between melting and solidification times and the temperature of the HTF has been experimentally determined in a cylindrical unit, where paraffin was used as the PCM [8]. The authors determined the thermal behaviour of the PCM through temperature measurements inside the PCM. Charging times have also been experimentally obtained in a shell-and-tube LHTES unit with different ratios of the tube–shell radius [9]. The authors measured the time-dependent temperature field within the PCM, and determined the ratio of the tube–shell radius with lower melting times and higher energy densities. Numerical predictions have been validated through experimental estimations of the liquid–solid front dynamics in cylindrical systems with electrical heating by a central rod [10]. The authors determined the PCMs melting fractions and temperature variations for different electrical power values. Melting and solidification experiments on three different paraffin types were carried out in tilted cylindrical units to determine the effects of the HTF temperature and flow rate, on the thermal performance of the PCM [11].

The interest in achieving higher energy densities has led to the usage of LHTES systems. Operating temperature ranges and materials are selected according to the type of application. The thermodynamics of PCMs required to analyse LHTES devices demands more sophisticated mathematical models and numerical methods to describe the dynamics of the phase transition. Finite volume element methods have been used to describe the freezing of supercooled liquid water [12]. The authors did not consider heat transfer through natural convection and volume changes upon freezing of liquid water. Heat exchange from natural convection has been taken into account during the phase-change process in confined systems; however, density changes induced by pressure increments during the freezing of liquid water are not considered by assuming incompressible phases [13,14]. The effects of natural convection have been considered for the prediction of the melting fraction, which was experimentally estimated from temperature field measurements in cylindrical units [10,15]. The authors, however, did not consider the volume changes of the system during the phase transition, since equal densities in the liquid and solid phase were assumed. Effects on the energy stored, charging times, and melting fractions produced by considering volume changes during phase transitions at constant pressure have been addressed on planar configurations [16,17,18,19]. The authors did not take into account the effects of supercooling (superheating) and natural convection during the solidification (melting) process, despite considering high temperature gradients. Mass-accommodation methods have also been used in planar cavities where volume changes were incorporated during the freezing of liquid water [13,14].

Volume changes have been incorporated by imposing total mass as a constant of the motion, where the length of the system was promoted as a dynamical variable. The equation of motion for the system’s length guarantees that mass is not created or destroyed during the phase-change process [16]. The consistency of the obtained solutions through the additional equation of motion was verified through the behaviour of the PCM during melting and solidification in adiabatic systems. Volume changes during phase transitions in cylindrical configurations can also be taken into account through total mass conservation. Other authors have introduced an additional equation of motion for the outer radius in cylindrical and spherical geometries [20]. The motion of the external radius is governed through an adiabatic boundary condition, which would create and destroy mass during the phase change process. Two different mass-accommodation methods for cylindrical configurations are proposed in this work. One of the proposed models introduces an additional equation of motion for the outer radius of the cylindrical unit and incorporates volume changes in the radial direction. A second mass-accommodation method takes into account volume changes along the axial direction of the cylindrical unit. Melting of solid produces an excess volume of liquid that may scatter throughout the top of the storage unit, or it may be pictured as being frequently extracted from the system. Mass conservation applied along the axial direction introduces an additional equation of motion for the excess volume of liquid, according to the second method used to accommodate mass in the system. Exact steady-state solutions for the external radius and for the radius of the region delimited by the liquid–solid interface are obtained in each case. Additionally, an experimental setup is designed to analyse the temperature field within paraffin wax used as a PCM, and placed in a vertical annular region with a rigid outer wall. Finally, the second method was applied to estimate the thermodynamic parameters of the paraffin through a least-square minimization procedure.

2. Mass-Accommodation Methods

The system under consideration consists of an aluminium container of length L that constitutes an annular region with internal radius

r_{0}

and external radius R. The region is occupied by a PCM and heat is transferred in a direction perpendicular to the axial symmetry axis z. Thermal energy is transported through a HTF that flows through a copper tube of radius

r_{0}

along the symmetry axis z. The temperature at the inner wall (copper–PCM interface) of the cylindrical unit is kept constant by the HTF, and it is homogeneously distributed along the inner wall in contact with the PCM. The temperature at the outer wall (aluminium–PCM interface) is approximately constant in time and homogeneously distributed. The temperature

T_{H}

at the copper–PCM interface is higher than the melting temperature

T_{m}

of the PCM, and the temperature

T_{C}

at the aluminium–PCM interface is below

T_{m}

. Thermal energy is transferred radially outwards and during the charging process, a portion of solid PCM is transformed into liquid. The homogeneous distribution of the temperature at each boundary surface guarantees that heat transfer along the z direction is negligible. The phase transition takes place at constant pressure and the temperature at the liquid–solid interface is constant and equal to the liquid–solid saturation temperature at thermodynamic equilibrium

T_{m}

. Therefore, superheating of solid phase during the melting process is not considered.

The temperature dependence of the thermodynamic variables is not considered. Buoyancy effects in the liquid phase are not incorporated since the temperature dependence of the liquid density is not being considered in this work. On the one hand, thermal expansion of liquid produced by temperature gradients, along with the orientation of a homogeneously heated surface relative to the gravitational field, can give rise to natural convection induced through buoyancy phenomena within the liquid [21]. On the other hand, it has been found experimentally that temperature changes in the liquid and solid phases are dominated by conduction when small liquid fractions are formed in annular regions where the heated surface is concentric to the outer surface and oriented along the axial direction [10]. The experimental results obtained in this work are used to determine the thermodynamic properties of paraffin wax at low melting fractions and when the system is in the conductive regime. Volume changes produced during the melting process due to the density difference between liquid and solid were incorporated by using two kinds of mass-accommodation methods.

2.1. Mass Accommodation through Radial Changes

The energy–mass balance equation at the liquid–solid interface depends on the method considered to incorporate volume changes. During melting or solidification, the system can expand or shrink in the radial direction and volume changes can be incorporated through total mass conservation. Through this method, one of the boundaries becomes a dynamic variable of motion as previously described in rectangular systems [16,22] and in confined PCMs [13,14], where incompressible phases were assumed. The energy–mass balance at the interface during a melting process in a cylindrical unit can be obtained by estimating the mass of liquid generated between t and

t + Δ t

as follows:

\frac{Δ M_{ℓ}}{Δ t} = \frac{M_{ℓ} (t + Δ t) - M_{ℓ} (t)}{Δ t},

(1)

where the mass of liquid at any instant in time t is:

M_{ℓ} (t) = π ρ_{ℓ} (\bar{r} {(t)}^{2} - r_{0}^{2}) L,

(2)

where

ρ_{ℓ}

is the density of the liquid phase and

\bar{r} (t)

is the radius of the boundary constituted by the liquid–solid interface at time t. The net amount of thermal energy absorbed by the solid at time t during a melting process can be obtained through the time derivative of the previous equation, as follows:

2 π L Δ h_{m} ρ_{ℓ} \bar{r} (t) \frac{d \bar{r} (t)}{d t} = ϕ (\bar{r} (t), t),

(3)

where

Δ h_{m}

is the enthalpy of formation or latent heat of fusion and

ϕ (r (t), t)

is the net thermal flux through the liquid–solid interface, given by:

ϕ (\bar{r} (t), t) = 2 π L \bar{r} (t) (- k_{ℓ} \frac{\partial T_{ℓ} (r, t)}{\partial r} |_{r = \bar{r} (t)} + k_{s} \frac{\partial T_{s} (r, t)}{\partial r} |_{r = \bar{r} (t)}),

(4)

where

k_{ℓ}

(

k_{s}

) is the thermal conductivity of the liquid(solid) phase and

T_{ℓ} (r, t)

(

T_{s} (r, t)

) is the temperature distribution in the liquid(solid) phase. Combining the last two equations, the energy–mass balance at

\bar{r} (t)

is given by

ρ_{ℓ} Δ h_{m} \frac{d \bar{r} (t)}{d t} = - k_{ℓ} \frac{\partial T_{ℓ} (r, t)}{\partial r} |_{r = \bar{r} (t)} + k_{s} \frac{\partial T_{s} (r, t)}{\partial r} |_{r = \bar{r} (t)} .

(5)

Additionally, the phase-change process can be described by estimating the mass of melted solid between t and

t + Δ t

, where

Δ M_{s} = M_{s} (t) - M_{s} (t + Δ t)

. Using the mass of solid phase located between

\bar{r} (t)

and

R (t)

, an equivalent energy–mass balance at the interface can be obtained through the time derivative of

M_{s} (t)

as follows:

2 π ρ_{s} L Δ h_{m} (R (t) \frac{d R (t)}{d t} - \bar{r} (t) \frac{d \bar{r} (t)}{d t}) = ϕ (r (t), t),

(6)

where

R (t)

is the outer radius that becomes a variable of the motion to preserve mass during the melting process. Substituting the net thermal flux at the interface given by Equation (2), an equivalent energy–mass balance at the interface can be found as follows

ρ_{s} Δ h_{m} (\frac{R (t)}{\bar{r} (t)} \frac{d R (t)}{d t} - \frac{d \bar{r} (t)}{d t}) = - k_{ℓ} \frac{\partial T_{ℓ} (r, t)}{\partial r} |_{r = \bar{r} (t)} + k_{s} \frac{\partial T_{s} (r, t)}{\partial r} |_{r = \bar{r} (t)} .

(7)

Equations (5) or (7) must be coupled to the differential equation obtained by imposing total mass as a constant of the motion. The outer radius

R (t)

becomes a dynamical variable to incorporate volume changes during the melting process. Imposing mass conservation through the time derivative of total mass, the additional equation of motion for

R (t)

is given by

ρ_{ℓ} \bar{r} (t) \frac{d \bar{r} (t)}{d t} + ρ_{s} (R (t) \frac{d R (t)}{d t} - \bar{r} (t) \frac{d \bar{r} (t)}{d t}) = 0 .

(8)

The time evolution of

\bar{r} (t)

and

R (t)

may be described through Equations (5) and (8), or equivalently, Equations (7) and (8), where the temperature distribution at each phase is found through the local energy balance given by:

ρ_{i} C_{i} \frac{\partial T_{i} (r, t)}{\partial t} = \frac{k_{i}}{r} \frac{\partial}{\partial r} (r \frac{\partial T_{i} (r, t)}{\partial r}),

(9)

where

C_{i}

is the specific heat capacity of phase i. Additionally, it is straightforward to show that Equations (5) or (7), can also be applied when the system is subjected to boundary conditions that produce solidification of liquid phase [16].

2.2. Mass Accommodation through Axial Growth

The second method that can be used to incorporate volume changes during a melting process consists of estimating the excess volume of liquid that grows beyond the top surface and when the outer surface at

r = R

is constant in time. The height of the liquid column increases beyond the top surface during a melting process, given that for most PCMs the liquid density is lower than the density of the solid phase. Mass accommodation must incorporate the height of the liquid column

Δ z (t)

as a dynamical variable instead of the outer radius R, as shown in Figure 1. The variable

Δ z (t)

, is related to the excess volume of liquid, as illustrated in Figure 1.

The current method promotes

Δ z (t)

instead of the outer radius R as the dynamical variable of motion. The energy–mass balance at the liquid–solid interface will show subtle changes in comparison with Equation (5) when considering volume displacements along the axial symmetry axis z. In this scenario, the mass of liquid

M_{ℓ} (t)

at some time t has a slightly more complicated form than

M_{ℓ} (t)

shown through Equation (1), and is given by:

M_{ℓ} (t) = π ρ_{ℓ} (\bar{r} {(t)}^{2} - r_{0}^{2}) (L + Δ z (t)),

(10)

where

Δ z (t)

represents the height of the excess volume of liquid, as illustrated in Figure 1. Alternatively, the expression for

M_{s} (t)

adopts a more simpler form since the outer radius R is constant; then,

M_{s} (t)

in this case is given by

M_{s} (t) = π ρ_{s} L (R^{2} - \bar{r} {(t)}^{2}) .

(11)

The last equation only incorporates

\bar{r} (t)

as a dynamical variable. In this case, the energy–mass balance at the interface adopts a simpler form by considering the amount of solid transformed into liquid, instead of using the mass of liquid given by Equation (10). The mass of melted solid can be found as

\frac{Δ M_{s}}{Δ t} = \frac{M_{s} (t) - M_{s} (t + Δ t)}{Δ t} .

(12)

The rate of melted solid is therefore, equal to

- d M_{s} (t) / d t

in the limit when

Δ t \to 0

. The rate of thermal energy absorbed at the liquid–solid interface can be obtained as:

2 π L Δ h_{m} ρ_{s} \bar{r} (t) \frac{d \bar{r} (t)}{d t} = ϕ (\bar{r} (t), t),

(13)

where

ϕ (\bar{r} (t), t)

is the net thermal flux at the interface given by Equation (4). The energy–mass balance equation is obtained by substituting the expression for

ϕ (\bar{r} (t), t)

as follows:

ρ_{s} Δ h_{m} \frac{d \bar{r} (t)}{d t} = - k_{ℓ} \frac{\partial T_{ℓ} (r, t)}{\partial r} |_{r = \bar{r} (t)} + k_{s} \frac{\partial T_{s} (r, t)}{\partial r} |_{r = \bar{r} (t)} .

(14)

which is almost identical to Equation (5); however, the solid density

ρ_{s}

must be used instead of

ρ_{ℓ}

when considering volume changes along the axial direction. The density that must appear in the energy-mass balance equation at the interface results from the manner in which volume changes are being incorporated and is consistent with mass conservation.

The excess volume of liquid, which is related to the variable

Δ z (t)

, can be estimated by imposing total mass conservation through the time derivative of the total mass

M (t) = M_{ℓ} (t) + M_{s} (t)

, as follows:

\bar{r} (t) \frac{d \bar{r} (t)}{d t} δ z (t) + \frac{1}{2} (\bar{r} {(t)}^{2} - r_{0}^{2}) \frac{d δ z (t)}{d t} - \bar{r} (t) \frac{d \bar{r} (t)}{d t} (\frac{ρ_{s}}{ρ_{ℓ}} - 1) = 0,

(15)

where

δ z (t) = Δ z (t) / L

represents the proportion of excess liquid and the volume of this liquid can be estimated as

Δ V_{ℓ} (t) = π L δ z (t) (\bar{r} {(t)}^{2} - r_{0}^{2})

. Equations (14) and (15) can be solved for the dynamic variables

\bar{r} (t)

and

δ z (t)

. On the one hand, if the excess volume of liquid is being frequently removed, local heat balance at each phase may be applied through Equation (9), neglecting the effects of heat transfer produced by the liquid scattered through the top surface. On the other hand, if

δ z (t) ≪ 1

during the entire charging process, heat transferred by the scattered liquid is negligible.

The phase transition can also be pictured by considering the rate of liquid mass formed during the melting process. An equivalent energy-mass balance equation at the interface can be found by considering the time derivative of the liquid mass

M_{ℓ} (t)

given by Equation (10) and using Equation (4) as follows

ρ_{ℓ} Δ h_{m} [\frac{d \bar{r} (t)}{d t} (1 + δ z (t)) + \frac{1}{2} (1 + {(\frac{r_{0}}{\bar{r} (t)})}^{2}) \frac{d δ z (t)}{d t}] = - k_{ℓ} \frac{\partial T_{ℓ} (r, t)}{\partial r} |_{r = \bar{r} (t)} + k_{s} \frac{\partial T_{s} (r, t)}{\partial r} |_{r = \bar{r} (t)} .

(16)

Conceiving the phase change in this way, the dynamical variables

\bar{r} (t)

and

δ z (t)

are coupled through Equation (16) and the time evolution of these variables must be obtained through the simultaneous solution of Equations (15) and (16). The system of equations

\bar{r} (t)

and

δ z (t)

and shown through Equations (15) and (16) is equivalent to the system given by Equations (14) and (15). Numerical solutions are easier to implement by picturing the phase-change process through Equations (14) and (15). In this situation,

δ z (t)

is not present in Equation (14). The problem can be solved only through Equation (14), independently of the equation of motion for

δ z (t)

, as if liquid mass was being destroyed in the process. Additionally, the amount of mass destroyed or that must be constantly removed can be obtained from the solution to Equation (15) and the value of

\bar{r} (t)

2.3. Steady-State Regime

The steady-state solutions in both cases were used to verify the consistency of the numerical solutions. The steady-state solutions were used to validate the numerical method considered in this work. Although, the mathematical model where the outer radius

R (t)

is promoted as a dynamical variable was not used for the analysis of the experimental results, it was also solved in this section for comparison with the solutions obtained through the second mass-accommodation method.

Asymptotic time limits can be found by using the steady-state solutions to the heat equation in each phase when the system is subjected to isothermal boundary conditions. The general solution to Equation (9) in the steady-state regime is the classical logarithmic function given by:

T_{i}^{(s s)} (r) = A_{i} ln (r) + B_{i},

(17)

where

T_{i}^{(s s)}

represents the temperature profile at phase i in the steady state (ss), and

A_{i}

B_{i}

are the corresponding constants of integration. The constants shown in the last equation can be found through the boundary conditions at

r = r_{0}

r = {\bar{r}}_{s s}

and

r = R_{s s}

, where

{\bar{r}}_{s s}

and

R_{s s}

represent the steady-state values of the region delimited by the liquid–solid interface and outer radius, respectively. The system is subjected to the homogeneous isothermal boundary conditions given by:

\begin{matrix} T (r_{0}, t) & = & T_{H}, \\ T (\bar{r} (t), t) & = & T_{m}, \\ T (R (t), t) & = & T_{C} . \end{matrix}

(18)

Applying the boundary conditions shown through Equation (18), the steady-state solution for the temperature profiles in each phase can be obtained in a straightforward manner, as follows:

\begin{matrix} T_{ℓ}^{(s s)} (r) & = & - \frac{Δ T_{H}}{ln (\frac{{\bar{r}}_{s s}}{r_{0}})} ln (\frac{r}{r_{0}}) + T_{H}, \\ T_{s}^{(s s)} (r) & = & - \frac{Δ T_{C}}{ln (\frac{R_{s s}}{{\bar{r}}_{s s}})} ln (\frac{r}{R_{s s}}) + T_{C}, \end{matrix}

(19)

where

Δ T_{H} = T_{H} - T_{m}

and

Δ T_{C} = T_{m} - T_{C}

. The value of

{\bar{r}}_{s}

can be obtained through the solution of Equation (5) in the steady state, as follows:

k_{s} \frac{d T_{s}^{(s s)} (r)}{d r} |_{r = {\bar{r}}_{s s}} = k_{ℓ} \frac{d T_{ℓ}^{(s s)} (r)}{d r} |_{r = {\bar{r}}_{s s}} .

(20)

The temperature profiles in the steady state given by Equation (19) can be substituted in the last equation to obtain an expression for

{\bar{r}}_{s s}

in terms of

R_{s s}

, which after some algebra is found, as follows

{\bar{r}}_{s s} = r_{0} {(\frac{R_{s s}}{r_{0}})}^{γ} with γ = \frac{k_{ℓ} Δ T_{H}}{k_{ℓ} Δ T_{H} + k_{s} Δ T_{C}} .

(21)

Total mass conservation can be used to obtain an additional condition for

{\bar{r}}_{s s}

and

R_{s s}

. Mass conservation in the steady state is given by:

ρ_{ℓ} ({\bar{r}}_{s s}^{2} - r_{0}^{2}) + ρ_{s} (R_{s s}^{2} - {\bar{r}}_{s s}^{2}) = ρ_{ℓ} (\bar{r} {(0)}^{2} - r_{0}^{2}) + ρ_{s} (R {(0)}^{2} - \bar{r} {(0)}^{2}),

(22)

where

\bar{r} (0)

is the initial radius of the region delimited by the liquid–solid interface and

R (0)

represents the initial value of the outer radius. The right hand side of Equation (22) is the initial mass per unit length of the cylindrical unit, which should be equal to its steady-state value. Equations (21) and (22) constitute a system of nonlinear equations for

{\bar{r}}_{s s}

and

R_{s s}

The solution for

{\bar{r}}_{s s}

that corresponds to the second mass-accommodation method can be found through Equation (21), when the outer radius is constant and equal to its initial value

R_{s s} = R (0)

. However, in this method, mass conservation is pictured as illustrated in Figure 1 and Equation (22) must be changed in order to accommodate mass along the axial direction, as follows:

ρ_{ℓ} ({\bar{r}}_{s s}^{2} - r_{0}^{2}) (1 + δ z_{s s}) + ρ_{s} (R {(0)}^{2} - {\bar{r}}_{s s}^{2}) = ρ_{ℓ} (\bar{r} {(0)}^{2} - r_{0}^{2}) + ρ_{s} (R {(0)}^{2} - \bar{r} {(0)}^{2}),

(23)

The steady-state value for

δ z_{s s}

can be obtained from Equation (23), by substituting the solution of

{\bar{r}}_{s s}

estimated from Equation (21) when

R_{s s} = R (0)

. The proportion of excess liquid

δ z_{s s}

in the steady state is then given by:

δ z_{s s} = \frac{(ρ_{s} / ρ_{ℓ} - 1) ({(R (0) / r_{0})}^{2 γ} - {(\bar{r} (0) / r_{0})}^{2})}{{(R (0) / r_{0})}^{2 γ} - 1},

(24)

where

γ

is defined through Equation (21). Finally, the proportion of scattered liquid mass

Δ M_{ℓ}^{(s s)}

compared to the total mass of liquid in the steady state

M_{ℓ}^{(s s)}

and defined as

δ^{(s s)} m_{ℓ} = Δ M_{ℓ}^{(s s)} / M_{ℓ}^{(s s)}

, is given by:

δ^{(s s)} m_{ℓ} = \frac{δ z_{s s}}{1 + δ z_{s s}},

(25)

where

δ z_{s s}

is given by Equation (24).

The consistency of the numerical solutions can be verified through comparison of the asymptotic time values obtained from the finite difference method, with the steady-state solutions given through Equations (21), (22), (24) and (25). Several numerical examples were performed by probing different values of

γ

. The relation between the steady-state values for

{\bar{r}}_{s s}

according to each mass-accommodation method will be obtained for increasing values of

γ

A special case can be found in the limit

k_{ℓ} ≫ k_{s}

. The steady-state value for

{\bar{r}}_{s s}

becomes equal to

R_{s s}

when

k_{ℓ} ≫ k_{s}

, and all the solid mass is melted during the process. In this limit, where

γ \to 1

and high melting fractions are expected, the steady-state value for

R_{s s}

can be obtained from total mass conservation described through Equation (22) as follows

R_{s s}^{2} = (1 - ρ_{s} / ρ_{ℓ}) \bar{r} {(0)}^{2} + \frac{ρ_{s}}{ρ_{ℓ}} R^{2} (0) .

(26)

Additionally, when the initial radius of the region delimited by the liquid–solid interface

\bar{r} (0) ≪ R (0)

, the steady-state value for

R_{s s}

at high melting fractions is approximately given by

R_{s s} \approx \sqrt{\frac{ρ_{s}}{ρ_{ℓ}}} R (0) .

(27)

Equation (27) predicts the maximum possible value for

R_{s s} = {\bar{r}}_{s s}

when

k_{ℓ} ≫ k_{s}

. The outer radius of the second mass-accommodation method is constant, and in this limit

{\bar{r}}_{s s} = R (0)

. The relation between the steady-state values of the radii

{\bar{r}}_{s s}

according to each mass-accommodation method in this limit is given by:

\frac{{\bar{r}}_{s s}^{(1)}}{{\bar{r}}_{s s}^{(2)}} = \sqrt{\frac{ρ_{s}}{ρ_{ℓ}}},

(28)

where

{\bar{r}}_{s s}^{(1)}

({\bar{r}}_{s s}^{(2)})

is the radius of the region delimited by the liquid–solid interface according to the first (second) method. Equation (28) predicts the maximum possible value for the relation between

{\bar{r}}_{s s}^{(1)}

and

{\bar{r}}_{s s}^{(2)}

according to the mass-accommodation methods discussed in this work. Additionally, this relation is also valid for the steady-state values of the outer radius, since

{\bar{r}}_{s s}^{(1)} = R_{s s}

and

{\bar{r}}_{s s}^{(2)} = R (0)

when

k_{ℓ} ≫ k_{s}

Equations (27) and (28) that correspond to high melting fractions are only valid when the thermal expansion of the liquid phase is not considered or is negligible within the temperature operating range of the cylindrical unit. Equations (27) and (28) were obtained by considering liquid PCMs that lie in the conductive regime, since natural convection induced through buoyancy effects is not considered in this work. Additionally, corrections to volume changes that arise from the thermal expansion of the liquid and solid are expected. Consequently, Equations (27) and (28) are only valid when

γ ≪ 1

, where

γ

is the dimensionless parameter defined through Equation (21), and some care must be taken when applying these equations to PCMs with high values of

α

or thermal units with high temperature gradients.

Numerical Examples

Figure 2 shows the behaviour in the steady-state regime of

{\bar{r}}_{s s}^{(1)} / {\bar{r}}_{s s}^{(2)}

and

R_{s s} / R (0)

for increasing values of

k_{ℓ}

. The thermodynamic properties of the PCM that belong to paraffin wax are shown in Table 1 [23]. The internal radius of the cylindrical unit is

r_{0} = 0.00635 m

. The initial values of the region delimited by the liquid–solid interface and outer radius are

\bar{r} (0) = 0.01 m

and

R (0) = 0.108 m

, respectively. Initial temperature profiles are constant and equal to

T_{ℓ} (r, 0) = 343.15 K

in the liquid domain and

T_{s} (r, 0) = 290.15 K

in the solid domain. Isothermal boundary conditions are applied to the system, where the copper–PCM interface at

r = r_{0}

is constant and equal to

T_{H} = 343.15 K

, and the outer surface at the aluminium–PCM interface

r = R (t)

is kept at a constant temperature value of

T_{C} = 290.15 K

The implicit finite difference method (FDM) [24] with a second-order approximation on the space derivatives was used to solve the mathematical model described through Equations (5), (8) and (9). Additionally, the FDM was also used to solve the model proposed through Equations (9), (14) and (15). Second-order approximations to the spacial derivatives and a first-order approximation to the time derivative were used to discretize the local energy–mass balance given by Equation (9) as follows:

(\frac{β_{i}}{2 r^{n}} - λ_{i}) T_{i}^{n - 1, j} + (2 λ_{i} + 1) T_{i}^{n, j} + (- \frac{β_{i}}{2 r^{n}} - λ_{i}) T_{i}^{n + 1, j} = T_{i}^{n, j - 1},

(29)

with

\begin{matrix} T_{i}^{n, j} = T (r^{n}, t^{j}), T_{i}^{n, j - 1} = T (r^{n}, t^{j - 1}), \end{matrix}

(30)

\begin{matrix} T_{i}^{n + 1, j} = T (r^{n + 1}, t^{j}), T_{i}^{n - 1, j} = T (r^{n - 1}, t^{j}), \end{matrix}

(31)

\begin{matrix} λ_{i} = \frac{k_{i} Δ t}{ρ_{i} C_{i} Δ r_{i}^{2}}, β_{i} = \frac{k_{i} Δ t}{ρ_{i} C_{i} Δ r_{i}} \end{matrix}

(32)

where

i = ℓ (s)

is used for liquid (solid) phase, n represents the nth node, j the jth time level and

Δ r_{i}

the node separation in each phase. The dimensionless parameter at each phase

λ_{i}

defined in the last equation depends on the thermodynamic properties of each phase, the node separation

Δ r_{i}

and the time step

Δ t

used during the simulations. Forward and backward second-order approximations to the space derivatives that appear in Equations (5) and (14) were used as follows:

\begin{matrix} {\bar{r}}^{j + 1} = {\bar{r}}^{j} - \frac{Δ t k_{ℓ}}{2 ρ_{i} Δ r_{ℓ} Δ h_{m}} (3 T_{ℓ}^{n, j} - 4 T_{ℓ}^{n - 1, j} + T_{ℓ}^{n - 2, j}) + \\ \frac{Δ t k_{s}}{2 ρ_{i} Δ r_{s} Δ h_{m}} (- 3 T_{s}^{n, j} + 4 T_{s}^{n + 1, j} - T_{s}^{n + 2, j}), \end{matrix}

(33)

where

ρ_{i} = ρ_{ℓ}

ρ_{i} = ρ_{s}

represents the density that is used in Equations (8) or (15), respectively. A first-order approximation to the time derivative of

\bar{r}

was used with an explicit scheme to estimate the radius of the region delimited by the liquid–solid interface at the next time level

{\bar{r}}^{j + 1}

Explicit schemes with a first-order approximation on the time derivatives were used to determine the outer radius

R^{j + 1}

or the excess liquid

δ z^{j + 1}

at the next time level according to each mass-accommodation method. In this approximation, the outer radius

R^{j + 1}

can be obtained from the discretized form of Equation (8) and using the value of

{\bar{r}}^{j + 1}

given by the previous equation with

ρ_{i} = ρ_{ℓ}

as follows

(\frac{ρ_{ℓ}}{ρ_{s}} - 1) {\bar{r}}^{j} ({\bar{r}}^{j + 1} - {\bar{r}}^{j}) + R^{j} (R^{j + 1} - R^{j}) = 0 .

(34)

The amount of excess liquid

δ z^{j + 1}

when using the second mass-accommodation method, can be obtained from the discretized form of Equation (15) and using

ρ_{i} = ρ_{s}

when estimating

{\bar{r}}^{j + 1}

, as follows

(1 - \frac{ρ_{s}}{ρ_{ℓ}} + δ z^{j}) {\bar{r}}^{j} ({\bar{r}}^{j + 1} - {\bar{r}}^{j}) + \frac{1}{2} ({({\bar{r}}^{j})}^{2} - r_{0}^{2}) (δ z^{j + 1} - δ z^{j}) = 0 .

(35)

Finally, a total number of 50 nodes in the liquid layer and 100 nodes in the solid phase were used during all simulations performed in this work. The time step on each of the numerical simulations performed is

Δ t = 0.001 s

The numerical results shown in Figure 2 correspond to the time asymptotic values for

\bar{r} (t)

and

R (t)

. The numerical results are compared with the steady-state values for

{\bar{r}}_{s s}

and

R_{s s}

given by the solution to the nonlinear system of Equations (21) and (22), and also with the steady-state value obtained through Equation (22) when the outer radius is constant and equal to

R (0)

Time asymptotic values were obtained by solving each model through the FDM and for several values of

k_{ℓ}

in the range

k_{ℓ} / k_{s} = [0.83, 4166.7]

. Figure 2 shows the logarithmic behaviour of

{\bar{r}}_{\lim}^{(1, FDM)} / {\bar{r}}_{\lim}^{(2, FDM)}

and

R_{\lim}^{(1, FDM)} / R (0)

, which is consistent with the result shown through Equation (28).

3. Experimental Setup

The characterization of the PCM was performed through the experimental setup shown in Figure 3a. Additionally, the setup is shown schematically in Figure 3b. The components used to perform the experimental tests and the thermocouple distribution on the whole experimental array are shown schematically in Figure 3b. The experimental array consists of the following main components:

1: A 10 L cylindrical container which constitutes a vertical annular region that is used to store paraffin wax. Four arrays of thermocouples distributed in concentric circles were placed inside the PCM in its solid phase. The inner radius of the annular region is formed by a $0.5 in$ copper tube that is placed along the axial symmetry axis of the cylindrical unit. Liquid water was gradually heated and circulated through the copper tube for thermal energy transfer at the copper–PCM interface.
2: A data acquisition system for temperature processing and data collection through the thermocouple array.
3: A system designed to control the liquid-water temperature and mass flow was developed.

Temperature sensing at the cylindrical TES unit consists of 22 K-type thermocouples as shown in Figure 4. Thermocouples used to collect data within the PCM were distributed in four sets of concentric circles with radii

r_{1} = 4.1 cm

r_{2} = 5.8 cm

r_{3} = 7.5 cm

and

r_{4} = 9.2 cm

. Additional thermocouples were used to measure the surrounding air temperature or ambient temperature, the copper–PCM interface temperature, the aluminium–PCM interface temperature, the inlet/outlet temperature of the liquid water, the temperature of the water at three equally spaced positions along the direction of the HTF mass flow, and the air temperature in direct contact with the aluminum shell.

The data acquisition system was developed through a data collection board (National Instruments 6062E DAQ PCMCIA), a module for thermocouple signal conditioning (NI SCXI-1102B), a 32-channel isothermal terminal block (NI SCXI-1303) and a rack for instrument housing (NI SCXI-1000). The National Instruments modules are shown in Figure 5a and a laptop with a code developed in LabView for monitoring and storing data is shown in Figure 5b. The HTF (water) mass flow and temperature was controlled through a Lauda™ ECO E 10 S Heating Thermostatic Bath. Temperature data were collected with the data acquisition system just described. Constant mass flow of liquid water which acted as the HTF was fixed at

10 liters / \min

. The HTF was gradually heated through the Lauda thermostatic bath and maximum temperature values were fixed at

70^{\circ} C

. Finally, experimental tests were carried out during a four hour period and temperature data was collected during the entire duration of the tests through the data acquisition system previously described.

4. Results and Discussion

The second mass-accommodation method, described in Section 2, was applied to analyse the experimental results obtained from a melting (charging) process of paraffin wax used as the PCM. The thermodynamic variables of the PCM were estimated by assuming that heat transfer within the paraffin is dominated by conduction and through the model presented in Section 2. The PCM is stored inside an annular region where water constitutes the HTF, and circulates through an inner copper tube with a radius of

r_{0} = 6.35 mm

concentric to an aluminium surface with an outer radius of

R = 10.8 cm

. Thermocouples for temperature measurements within the PCM domain were placed at a height equal to

h = L / 2

, where

L = 10 cm

represents the total height of the heat storage unit. Thermal energy is transferred to the PCM through liquid water that circulates within the inner tube. The water acting as the HTF absorbs thermal energy from a heat bath with a thermostat fixed to

70^{\circ} C

. The outer radius is in contact with the surrounding air at ambient temperature. Liquid water is gradually heated and circulated through the inner tube from an initial temperature of

17.5^{\circ}

C until the water temperature reaches a steady-state value of

70 \pm 1^{\circ}

C. Temperature data collection within the paraffin started from the instant in which the temperature at

r = r_{0}

reached the melting temperature

T_{m}

of the PCM.

The temperature sensing was performed through a distribution of thermocouples in four concentric circles at fixed radii of:

r_{1} = 4.1 cm

r_{2} = 5.8 cm

r_{3} = 7.5 cm

and

r_{4} = 9.2 cm

. Four thermocouples were placed at equal angular separations in each concentric circle, as shown in Figure 4. Therefore, the total number of thermocouples used to measure the temperature profile within the paraffin wax was 16. Each set of four sensors along the radial direction were slightly tilted, as shown in Figure 4, to minimize errors in temperature measurements due to the thermal energy absorbed by the nearest thermocouples. The temperature at each radius

r_{i}

was estimated through the temperature average obtained from the four thermocouples distributed along each concentric circle. Additionally, one thermocouple was placed at the copper–PCM interface to estimate the temperature at

r = r_{0}

. Finally, four thermocouples were placed at the aluminium–PCM interface to determine experimental values of the temperature at

r = R

Figure 6 shows the average temperature as a function of time, obtained from the thermocouple readings at the copper–PCM and aluminium–PCM interface. Temperature values were registered approximately every second from the instant in which the copper–PCM interface reaches an average value equal to

T_{m} = 55^{\circ} C

. The temperature at the copper–PCM interface was obtained from the experimental values shown in Figure 6 and defined as a piecewise function for the numerical simulations. Two sections were obtained and each section was approximated through a polynomial fit with the highest correlation, as illustrated in Figure 6. Additionally, the temperature at the aluminium–PCM interface for the numerical simulations was obtained through a polynomial fit with the highest correlation.

The boundary condition at the copper–PCM interface used to solve the model described in Section 2 and shown in Figure 6 is given by:

T_{ℓ} (r_{0}, t) = \{\begin{matrix} 0.027 t + 24.37, & 1191.7 \leq t \leq 1637.4 \\ 3.34 \times 10^{- 12} t^{3} - 1.01 \times 10^{- 7} t^{2} + 10.06 \times 10^{- 4} t + 65.86, & 1638.3 \leq t \leq 14402, \end{matrix}

(36)

where the temperature ranges obtained in the time domain

1191.7 \leq t \leq 1637.4

and

1638.3 \leq t \leq 14402

are

[56.55, 68.58]^{\circ}

C and

[67.25, 69.38]^{\circ}

C, respectively.

The root-mean-squared error (

rmse

) of each polynomial fit was obtained as follows:

rmse = \sqrt{\frac{1}{N} \sum_{i = 1}^{N} {(T_{i}^{(\exp)} - T_{i}^{(fit)})}^{2}},

(37)

where N is the total number of observations or temperature readings and

T^{(\exp)}

(

T^{(fit)}

) represents the experimental (fitted) temperature values. The

rmse

obtained from the linear and cubic functions shown in Equation (36) is

rmse = 0.1804^{^{\circ}} C

and

rmse = 0.1147^{\circ}

C, respectively. Additionally, the correlation coefficient

r_{c}^{2}

was determined from several polynomial fits, and the function with the highest value of

r_{c}^{2}

is shown in Equation (36). The correlation coefficient was obtained through the following relation:

r_{c}^{2} = 1 - \frac{\sum_{i = 1}^{N} {(T_{i}^{(\exp)} - T_{i}^{(fit)})}^{2}}{\sum_{i = 1}^{N} {(T_{i}^{(\exp)} - T_{avg}^{(\exp)})}^{2}},

(38)

where

T_{avg}^{(\exp)}

(

T_{avg}^{(fit)}

) is the average temperature of the experimental (fitted) data. The highest correlation coefficient found corresponds to a linear and cubic function where

r_{c}^{2} = 0.9973

and

r_{c}^{2} = 0.9495

, respectively.

Finally, the boundary condition at the aluminium–PCM interface was determined from the experimental data shown in Figure 6 as follows:

T_{s} (R (t), t) = - 7.56 \times 10^{- 9} t^{2} + 6.13 \times 10^{- 4} t + 16.78 for 1191.7 \leq t \leq 14402,

(39)

where the temperature range in the time domain given by the above equation is

[17.50, 24.04]^{\circ}

C as illustrated in Figure 6. The polynomial shown in the last equation corresponds to the function with the highest correlation of

r_{c}^{2} = 0.9985

and a root-mean-squared error of

rmse = 0.0731^{\circ}

The temperature at each value of

r_{i}

considered in the experimental setup shown in Figure 4 was obtained through the temperature-averaged readings, registered by the set of thermocouples distributed along each concentric circle. The average temperature at each radius was registered at equally spaced time intervals, as shown in Figure 7.

Liquid and solid samples of paraffin wax were prepared to estimate the liquid and solid densities. The density of liquid paraffin was determined by pouring a fixed volume of PCM on a beaker previously placed on a scale. Four samples of liquid PCM with different volumes at

70^{\circ}

C were used to estimate the density of the liquid phase, and an average density of

ρ_{ℓ} = 735.25 kg / m^{3}

was estimated. Similarly, four solid cylindrical samples of PCM with different volumes at

20^{\circ}

C were used to determine the density of the solid phase. The volume and mass of each solid cylinder was measured and an average density of solid PCM

ρ_{s} = 849.39 kg / m^{3}

was estimated. The melting temperature

T_{m} = 55^{\circ}

C was determined through the liquid–solid coexistence of the paraffin wax close to thermodynamic equilibrium.

Thermal conductivities and specific heat capacities were estimated through the second mass-accommodation method described in Section 2 and using the nonhomogeneous isothermal boundary conditions given by Equations (36) and (39). The temperature dependence of the thermodynamic variables in the operating temperature range of the cylindrical unit of

[17, 70]^{\circ}

C was not considered. The latent heat of fusion can be approximated as

Δ h_{m} = (C_{ℓ} - C_{s}) T_{m}

assuming that

C_{ℓ}

and

C_{s}

are close to their saturation values at

T_{m} = 55^{\circ} C

. The mass-accommodation method described through Equations (14) and (15) for the dynamical variables

\bar{r} (t)

and

δ z (t)

, and the local energy balance at each phase given by Equation (9) was implemented through the implicit FDM previously described. The initial radius of the region delimited by the liquid–solid interface

\bar{r} (0)

was very close to the radius of the inner tube, since data collection started very close to the melting temperature of the PCM. A value of

\bar{r} (0) = 6.36 mm

was used and a temperature of

55^{\circ}

C was established at each node within the initial liquid layer. The HTF was circulated at ambient temperature and gradually heated before data collection, until the copper–PCM interface reached the melting temperature of the PCM. During this previous stage, the PCM was found in its solid state and a logarithmic temperature profile was registered by the thermocouples when data acquisition started. The logarithmic temperature distribution obtained from the gradual heating of the solid PCM was used as the initial temperature profile in the solid phase. Low melting fractions are expected in the range of thermal conductivities considered during the analysis of experimental data and according to the models introduced in this work.

Figure 7 shows the results obtained for the temperature at each concentric circle of radius

r_{i}

. The time evolution of the average temperature value at each radial coordinate

r_{i}

is shown in Figure 7. The FDM solutions to the model previously described are also shown in solid lines. The numerical solutions were obtained for several possible values of

k_{ℓ}

k_{s}

C_{ℓ}

and

C_{s}

. The numerical results were compared with the experimental temperature values shown in Figure 7, and a quadratic error function was defined to find the set of thermodynamical parameters that best reproduce the experimental data. The quadratic error used was defined as follows:

E_{T} (k_{ℓ}, k_{s}, C_{ℓ}, C_{s}) = \sum_{j = 1}^{n} {(T_{i}^{(\exp)} (t_{j}) - T_{i}^{(F D M)} (t_{j}))}^{2},

(40)

where

T_{i}^{(\exp)} (t_{j})

represents the average temperature registered by the ith thermocouple and

T_{i}^{(F D M)} (t_{j})

is the temperature obtained through the numerical solution of the second mass-accommodation method at each time value

t_{j}

shown in Figure 7 and at each thermocouple radial position

r_{i}

. The error function is evaluated for each particular set of thermodynamic parameters

k_{ℓ}, k_{s}, C_{ℓ}, C_{s}

, where a total of

10^{4}

sets of thermodynamic parameters was investigated in a range of possible values to find the least square error. The result is illustrated in Figure 7, where the average temperatures registered at each sensor are shown in red circles. The error bars correspond to the standard deviation obtained from the temperature readings registered by each of the four sensors distributed along a particular concentric circle of radius

r_{i}

. Figure 7 also shows the numerical result obtained with the set of thermodynamic parameters that minimize the error defined through Equation (40). The numerical results shown in Figure 7 correspond to the solutions of the second mass-accommodation method described previously and defined through Equations (9), (14) and (15).

Collected temperature values

T_{i}^{(\exp)}

through the thermocouple data acquisition system are shown in Figure 7 at the radial coordinates

r_{1} = 4.1 cm

r_{2} = 5.8 cm

r_{3} = 7.5 cm

and

r_{4} = 9.2 cm

. The solutions obtained with the FDM and according to the second mass-accommodation method

T_{i}^{(2, FDM)}

constitute the numerical solutions with the set of thermodynamic parameters shown in Table 2 that minimize the error given by Equation (39). According to the definition of the latent heat as the difference between liquid and solid enthalpies at the saturation temperature, and according to the specific heat capacities shown in Table 2, the latent heat of fusion of the paraffin wax estimated in this work is

Δ h_{m} = 206.73 kJ / kg

According to the second mass-accommodation model and the type of boundary conditions applied during the experimental tests, small values of

\bar{r} (t) ≪ R (0)

are expected since

{\bar{r}}_{s s} ≪ R (0)

when

k_{ℓ} ≪ k_{s}

. Close to this limit, when

γ ≪ 1

, and according to Equation (21), the radial displacement of the liquid–solid interface is very small compared to the outer radius as shown in Figure 8. The time evolution of the fraction of melted solid

f_{s} = (M_{s} (0) - M_{s} (t)) / M_{s} (0)

was obtained during the time domain shown in Figure 8. The fraction of melted solid was estimated through the set of thermodynamic variables that minimizes the quadratic error given by Equation (40) and shown in Table 2. The results illustrated in Figure 8 show that heat transfer within the liquid phase lies in the conductive regime where a small thickness of liquid PCM layer is formed, and very small fractions of melted solid

f_{s} \approx 10^{- 3}

are observed.

5. Conclusions

In this work, two mass-accommodation methods that consider volume changes during PCM melting in a vertical annular region were proposed. Mass conservation introduces an additional equation of motion that depends on the way in which mass is distributed. On the one hand, the outer radius was promoted as a dynamical variable to accommodate mass in the radial direction. On the other hand, the excess volume of liquid becomes a dynamical variable that is used to accommodate mass in the axial direction. We also found that the local energy–mass balance at the liquid–solid interface depends on the particular mass-accommodation method being used. Exact analytical expressions were found for the steady-state values in each case. Additionally, the relationship between the dynamical variables in the steady state and for high melting rates is proportional to the relationship between liquid and solid densities. The numerical solutions were verified through the steady-state solutions found in this work and for increasing values of

k_{ℓ}

. The thermodynamic properties of a paraffin were estimated through a minimization process. Experimental measurements of the temperature field within the PCM were used to characterize the behaviour of the PCM at low melting rates. A series of numerical simulations with different thermodynamic parameter values were performed to find the set of thermodynamic variables that minimizes the error between the experimental and numerical results. Finally, we confirm that conduction constitutes the main heat-transfer mechanism in cylindrical units with low melting rates.

Author Contributions

Conceptualization and discussion of mathematical models was conducted by V.S.-N., J.A.O. and E.M.H.-C.; numerical simulations were developed and obtained by V.S.-N. and J.A.O.; experimental setup and tests designed and performed by V.S.-N., J.E.C.-Q. and J.A.O.; the manuscript was written and revised by V.S.-N., J.A.O., E.M.H.-C. and J.E.C.-Q. and figures were elaborated by J.A.O. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations and symbols are used in this manuscript:

CSP	Concentrating Solar Power
PCM	Phase-change material
TES	Thermal energy storage
LHTES	Latent heat thermal energy storage
HTF	Heat-transfer fluid
FDM	Finite difference method
rmse	Root-mean-squared error
$k_{ℓ}$	Thermal conductivity of the liquid
$k_{s}$	Thermal conductivity of the solid
$C_{ℓ}$	Specific heat capacity of the liquid
$C_{s}$	Specific heat capacity of the solid
$ρ_{ℓ}$	Liquid density
$ρ_{s}$	Solid density
$Δ h_{m}$	Latent heat of fusion
$T_{m}$	Melting temperature
L	Height of cylindrical unit
R	Outer radius
$r_{0}$	Copper tube radius
$\bar{r}$	Radius of region delimited by the liquid–solid interface
$R_{s s}$	Steady state of outer radius
${\bar{r}}_{s s}$	Steady-state value of $\bar{r}$
$T_{ℓ} (r_{0}, t)$	Temperature at the copper–PCM interface
$T_{s} (R, t)$	Temperature at the aluminium–PCM interface
$γ$	Dimensionless exponent
$r_{c}$	Correlation coefficient
$f_{s}$	Fraction of melted solid

References

Gil, A.; Medrano, M.; Martorell, I.; Lázaro, A.; Dolado, P.; Zalba, B.; Cabeza, L.F. State of the art on high temperature thermal energy storage for power generation. Part 1-Concepts, materials and modellization. Renew. Sustain. Energy Rev. 2010, 14, 31–55. [Google Scholar] [CrossRef]
Mathur, A.; Kasetty, R.; Oxley, J.; Mendez, J.; Nithyanandam, K. Using encapsulated phase change salts for concentrated solar power plant. Energy Procedia 2014, 49, 908–915. [Google Scholar] [CrossRef] [Green Version]
Liu, M.; Saman, W.; Bruno, F. Review on storage materials and thermal performance enhancement techniques for high temperature phase change thermal storage systems. Renew. Sustain. Energy Rev. 2012, 16, 2118–2132. [Google Scholar] [CrossRef]
Cabeza, L.F.; Ibanez, M.; Sole, C.; Roca, J.; Nogues, M. Experimentation with a water tank including a PCM module. Sol. Energy Mater. Sol. Cells 2006, 90, 1273–1282. [Google Scholar] [CrossRef]
Akgün, M.; Aydın, O.; Kaygusuz, K. Experimental study on melting/solidification characteristics of a paraffin as PCM. Energy Convers. Manag. 2007, 48, 669–678. [Google Scholar] [CrossRef]
Mao, Q.; Li, Y. Experimental and numerical investigation on enhancing heat transfer performance of a phase change thermal storage tank. J. Energy Storage 2020, 31, 101725. [Google Scholar] [CrossRef]
Mao, Q.; Li, Y.; Li, G.; Badiei, A. Study on the influence of tank structure and fin configuration on heat transfer performance of phase change thermal storage system. Energy 2021, 235, 121382. [Google Scholar] [CrossRef]
Avci, M.; Yazici, M.Y. Experimental study of thermal energy storage characteristics of a paraffin in a horizontal tube-in-shell storage unit. Energy Convers. Manag. 2013, 73, 271–277. [Google Scholar] [CrossRef]
Seddegh, S.; Wang, X.; Joybari, M.M.; Haghighat, F. Investigation of the effect of geometric and operating parameters on thermal behavior of vertical shell-and-tube latent heat energy storage systems. Energy 2017, 137, 69–82. [Google Scholar] [CrossRef] [Green Version]
Chen, Y.; Peng, C.; Guo, Y. Experimental and numerical study on melting process of paraffin in a vertical annular cylinder. Therm. Sci. 2019, 23, 525–535. [Google Scholar] [CrossRef] [Green Version]
Akgün, M.; Aydın, O.; Kaygusuz, K. Thermal energy storage performance of paraffin in a novel tube-in-shell system. Appl. Therm. Eng. 2008, 28, 405–413. [Google Scholar] [CrossRef]
Kaiser, J.; Adami, S.; Akhatov, I.; Adams, N. A semi-implicit conservative sharp-interface method for liquid-solid phase transition. Int. J. Heat Mass Transf. 2020, 155, 119800. [Google Scholar] [CrossRef]
Dallaire, J.; Gosselin, L. Various ways to take into account density change in solid-liquid phase change models: Formulation and consequences. Int. J. Heat Mass Transf. 2016, 103, 672. [Google Scholar] [CrossRef]
Dallaire, J.; Gosselin, L. Numerical modeling of solid-liquid phase change in a closed 2D cavity with density change, elastic wall and natural convection. Int. J. Heat Mass Transf. 2017, 114, 903. [Google Scholar] [CrossRef]
Wang, S.; Guo, Y.; Peng, C.; Wang, W. Experimental study of paraffin melting in cylindrical cavity with central electric heating rod. J. Environ. Manag. 2019, 237, 264–271. [Google Scholar] [CrossRef] [PubMed]
Santiago, R.D.; Hernández, E.M.; Otero, J.A. Constant mass model for the liquid-solid phase transition on a one-dimensional Stefan problem: Transient and steady state regimes. Int. J. Therm. Sci. 2017, 118, 40. [Google Scholar] [CrossRef]
Santiago Acosta, R.D.; Otero, J.A.; Hernández Cooper, E.M.; Pérez-Álvarez, R. Thermal expansion effects on the one-dimensional liquid-solid phase transition in high temperature phase change materials. AIP Adv. 2019, 9, 025125. [Google Scholar] [CrossRef] [Green Version]
Rodriguez-Aleman, S.; Hernandez-Cooper, E.M.; Otero, J.A. Consequences of total thermal balance during melting and solidification of high temperature phase change materials. Therm. Sci. Eng. Prog. 2020, 20, 100750. [Google Scholar] [CrossRef]
Rodríguez-Alemán, S.; Hernández-Cooper, E.M.; Pérez-Álvarez, R.; Otero, J.A. Effects of Total Thermal Balance on the Thermal Energy Absorbed or Released by a High-Temperature Phase Change Material. Molecules 2021, 26, 365. [Google Scholar] [CrossRef] [PubMed]
Caldwell, J.; Chiu, C. Numerical solution to two-phase Stefan problems by the heat balance integral method. WIT Trans. Eng. Sci. 1994, 5, 10. [Google Scholar]
Madruga, S.; Curbelo, J. Effect of the inclination angle on the transient melting dynamics and heat transfer of a phase change material. Phys. Fluids 2021, 33, 055110. [Google Scholar] [CrossRef]
Hetmaniok, E.; Slota, D.; Zielonka, A. Solution of the Direct Alloy Solidification Problem including the Phenomenon of Material Shrinkage. Therm. Sci. 2017, 21, 105. [Google Scholar] [CrossRef] [Green Version]
Haji-Sheikh, A.; Eftekhar, J.; Lou, D. Some thermophysical properties of paraffin wax as a thermal storage medium. In Proceedings of the 3rd Joint Thermophysics, Fluids, Plasma and Heat Transfer Conference, St. Louis, MO, USA, 7–11 June 1982; p. 846. [Google Scholar]
Hernández, E.; Castillo, F.; Santiago, R.; Martínez, R.; Otero, J.; Oseguera, J.; Jiménez, A. Numerical and semi-analytical solutions of the compound layer growth kinetics in cylindrical surfaces during plasma nitriding of pure iron. In Journal of Physics: Conference Series; IOP Publishing: Bristol, UK, 2017; Volume 792, p. 012056. [Google Scholar]

Figure 1. Schematic representation of liquid volume growth in the axial direction. The height of the column that represents the excess liquid at any time t is

Δ z (t)

. The volume of this liquid at some time t is

Δ V_{ℓ} (t) = (\bar{r} {(t)}^{2} - r_{0}^{2}) (L + Δ z (t))

and represents the liquid that will scatter throughout the top surface of the cylinder or the volume of liquid that must be removed from the cylindrical unit.

Figure 1. Schematic representation of liquid volume growth in the axial direction. The height of the column that represents the excess liquid at any time t is

Δ z (t)

. The volume of this liquid at some time t is

Δ V_{ℓ} (t) = (\bar{r} {(t)}^{2} - r_{0}^{2}) (L + Δ z (t))

and represents the liquid that will scatter throughout the top surface of the cylinder or the volume of liquid that must be removed from the cylindrical unit.

Figure 2. Logarithmic relation between the two mass-accommodation methods previously discussed

ln ({\bar{r}}^{(1)} / {\bar{r}}^{(2)})

and

ln (R^{(1)} / R^{(2)})

, obtained from the exact steady-state solutions given through Equations (21) and (22), and the asymptotic time values estimated through the FDM. Asterisk and cross symbols are used to represent the relation obtained through the exact steady-state values. Empty circles and squares represent the relation between the asymptotic time limits according to the numerical solutions for each mass-accommodation method. The dashed line corresponds to the predicted relation for high melting fractions according to Equation (28).

Figure 2. Logarithmic relation between the two mass-accommodation methods previously discussed

ln ({\bar{r}}^{(1)} / {\bar{r}}^{(2)})

and

ln (R^{(1)} / R^{(2)})

Figure 3. (a) Experimental setup with the cylindrical unit and thermocouple array. (b) Schematic representation of the experimental setup with the following components: 1. Lauda Thermostatic Bath, 2. Thermal energy-storage unit, 3. SCXI-1000 National Instruments module for thermocouple signal conditioning, 4. Laptop for data processing, A. 1 K-type thermocouple for heat bath temperature sensing, B. 1 K-type thermocouple for HTF inlet temperature sensing, C. 22 K-type thermocouple array for temperature sensing at the copper–PCM interface, aluminium–PCM interface and PCM temperature field estimation and D. 1 K-type thermocouple for HTF outlet temperature sensing. E. 1 K-type thermocouple for ambient temperature sensing.

Figure 4. Cylindrical unit with thermocouple array for temperature sensing.

Figure 5. Data acquisition and processing system. (a) Data acquisition system. (b) Laptop for temperature data processing.

Figure 6. Nonhomogeneous isothermal boundary conditions. (a) Symbols represent time-dependent temperature values registered by the thermocouple located at the copper–PCM interface and the solid line corresponds to the best fit with the highest correlation. (b) Time-dependent temperature readings at the aluminium–PCM interface. Symbols represent the temperature values registered by the thermocouples located at

r = 0.108 m

and the solid line is the best fit with the highest correlation.

r = 0.108 m

and the solid line is the best fit with the highest correlation.

Figure 7. Time evolution of the temperature at each thermocouple radial position according to the the experimental and numerical results. Experimental values of the average temperature at each radial coordinate (a)

r_{1} = 4.1 cm

, (b)

r_{2} = 5.8 cm

, (c)

r_{3} = 7.5 cm

and (d)

r_{4} = 9.2 cm

, respectively are shown in red circles and the result obtained through the FDM is shown in solid lines. The numerical result was obtained through the solution of the model described by Equations (9), (14) and (15) and through the set of thermodynamic variables

k_{ℓ}, k_{s}, C_{ℓ} and C_{s}

with the lowest quadratic error shown by Equation (40).

r_{1} = 4.1 cm

, (b)

r_{2} = 5.8 cm

, (c)

r_{3} = 7.5 cm

and (d)

r_{4} = 9.2 cm

k_{ℓ}, k_{s}, C_{ℓ} and C_{s}

with the lowest quadratic error shown by Equation (40).

Figure 8. (a) Numerical solution to the liquid–solid interface motion. The results illustrate the relation between the liquid’s thickness and the outer radius of the thermal unit, according to the numerical solution of the proposed model with the set of thermodynamic variables shown in Table 2. (b) Time evolution of

f_{s}

obtained from the numerical solutions to the proposed model and with the set of parameters shown in Table 2.

f_{s}

obtained from the numerical solutions to the proposed model and with the set of parameters shown in Table 2.

Table 1. Thermodynamic properties and liquid–solid saturation values of paraffin wax according to the authors of Ref. [23].

$T_{m}$	$Δ h_{m}$	$k_{ℓ}$ ( $k_{s}$ )	$C_{ℓ}$ ( $C_{s}$ )	$ρ_{ℓ}$ ( $ρ_{s}$ )
$K$	( $kJ / kg$ )	( $W / m \cdot K$ )	( $kJ / kg \cdot K$ )	( $kg / m^{3}$ )
317	$266.0$	$0.24$ ( $0.24$ )	$2.95$ ( $2.51$ )	760 (818)

Table 2. Thermodynamic properties and liquid–solid saturation values of paraffin wax estimated through the experimental measurements and minimization of the quadratic error defined in this work.

$T_{m}$	$Δ h_{m}$	$k_{ℓ}$ ( $k_{s}$ )	$C_{ℓ}$ ( $C_{s}$ )	$ρ_{ℓ}$ ( $ρ_{s}$ )
$K$	( $kJ / kg$ )	( $W / m \cdot K$ )	( $kJ / kg \cdot K$ )	( $kg / m^{3}$ )
328	$206.73$	$0.05$ ( $0.317$ )	$3.20$ ( $2.57$ )	760(818)

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Silva-Nava, V.; Hernández-Cooper, E.M.; Chong-Quero, J.E.; Otero, J.A. Thermophysical Characterization of Paraffin Wax Based on Mass-Accommodation Methods Applied to a Cylindrical Thermal Energy-Storage Unit. Molecules 2022, 27, 1189. https://doi.org/10.3390/molecules27041189

AMA Style

Silva-Nava V, Hernández-Cooper EM, Chong-Quero JE, Otero JA. Thermophysical Characterization of Paraffin Wax Based on Mass-Accommodation Methods Applied to a Cylindrical Thermal Energy-Storage Unit. Molecules. 2022; 27(4):1189. https://doi.org/10.3390/molecules27041189

Chicago/Turabian Style

Silva-Nava, Valter, Ernesto M. Hernández-Cooper, Jesús Enrique Chong-Quero, and José A. Otero. 2022. "Thermophysical Characterization of Paraffin Wax Based on Mass-Accommodation Methods Applied to a Cylindrical Thermal Energy-Storage Unit" Molecules 27, no. 4: 1189. https://doi.org/10.3390/molecules27041189

APA Style

Silva-Nava, V., Hernández-Cooper, E. M., Chong-Quero, J. E., & Otero, J. A. (2022). Thermophysical Characterization of Paraffin Wax Based on Mass-Accommodation Methods Applied to a Cylindrical Thermal Energy-Storage Unit. Molecules, 27(4), 1189. https://doi.org/10.3390/molecules27041189

Article Menu

Thermophysical Characterization of Paraffin Wax Based on Mass-Accommodation Methods Applied to a Cylindrical Thermal Energy-Storage Unit

Abstract

1. Introduction

2. Mass-Accommodation Methods

2.1. Mass Accommodation through Radial Changes

2.2. Mass Accommodation through Axial Growth

2.3. Steady-State Regime

Numerical Examples

3. Experimental Setup

4. Results and Discussion

5. Conclusions

Author Contributions

Funding

Data Availability Statement

Conflicts of Interest

Abbreviations

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI