Correlations for Melting Time and Melt Fraction for Bottom Charged Tube in Tube Phase Change Material Heat Exchanger


 Multiple factors govern the Thermo-hydraulic behaviour of Latent heat storage devices. The correlation among these factors varies from case to case. In this work, a concentric tube in tube latent heat storage system is numerically modelled for the bottom charging case. Fixed grid enthalpy porosity approach is adopted to account for phase change. The numerical model’s independence is achieved by testing mesh size, time step, and maximum iterations per time step. The computational approach is validated against the experimental data. Non-dimensional parameters viz Rayleigh Number (3.04x105 to 65.75 x105), Stefan Number (0.2 to 1), Reynolds Number (600 to 3000), and L/D ratio (2 to 15) are varied in the respective ranges mentioned in parenthesis. Stefan number is found to have a major influence on the Melt Fraction and Melting time, compared to Rayleigh Number and Reynolds Number. Correlations are presented for quantifying the melt fraction and dimensionless melting time.


Introduction
Energy demand has increased the need for diligent management of energy. Proportionately, Thermal Energy Storage (TES) has attracted a lot of attention in recent decades. TES has the potential to maximize the benefits of the energy sector by storage and intermittent use. Among TES systems, latent heat storage is unique in their ability to provide more heat storage compared to sensible heat storage.
Various kinds [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17] of latent heat storage devices have been adopted to conserve energy. Multiple parameters influence the design of these latent heat storage devices. It would be handy for designers if these parameters are correlated through mathematical expressions. Efforts made by multiple researchers in this regard are presented hereunder. These correlations form two major divisions. The first is nondimensional melting time, i.e., Fourier Number (Fo) equated with one or more dependent parameters. The second one is Melt Fraction (MF), varying from zero to one, expressed in terms of multiple non-dimensional parameters. In this section, Fo correlations are discussed first and followed by MF correlations.
Riley [1] conducted a study on the solidification of metals in cylindrical and spherical systems for Stefan Numbers ranging from 4 to 20. This study has proposed the following mathematical correlation (Eq.1) for solidification time. In continuation to the above work, Hasan [4] employed stearic acid as PCM in a tube in tube system for 0.2≤ St ≤ 0.357 and proposed a similar correlation. cm = 0.14 +

0.25
The above correlations for the dimensionless melting time have considered that Stefan number is the only parameter influencing the melting. Although the Stefan number is a significant influencer, the other parameters also affect the melting time. In this connection, Rathod [6]  Eq.4 Eq.5 Eq.6 Cylindrical Bottom charged tube in tube system considered in this work resulted in the following expression. Paraffin and Water are used as PCM and HTF, respectively, in their work.
In extension to the above work, Rathod [7] included one more dimensionless parameter as a ratio of inlet temperature and initial temperature. Paraffin and Stearic acid are used as PCMs. However, this correlation (Eq.8) has resulted in a maximum error of 30.55%.
Kalapala [8] has conducted numerical simulations on a vertical top charged tube in tube system using lauric acid and water. Apart from St and Re, other parameters viz., Rayleigh Number (Ra), Length to Diameter (L/D) ratio, Ratio of Thermal Diffusivities of tube material and PCM(φ) and Tube Thickness to Diameter Ratio (σ) are also included in the Fo cm correlation(Eq.9).  Fan [10] has employed graphene-based nanocomposite PCMs in a vertical cylindrical system for various loadings, i.e., the weight percentage of nanoparticles. a, b, c values (Eq.11) vary slightly with the given weight percentage. For zero weight percentage, i.e., PCM without nanoparticles, the correlation is Eq.12. So far, various correlations available in the literature are summarized. Among these, Kalapala [8] did a detailed analysis on the influence of each dimensionless parameter on melting for the Top Charging case. In this present work, a similar attempt is made to understand the impact of dimensionless parameters in Bottom Charging and to establish the correlation for the Bottom charging case.

MF = a Fo
Kalapala [8] has considered correlation among Melt Fraction, St, Fo, Ra, Re, φ, σ, and L/D ratio. However, in the same work [8], it has been concluded that the parameters φ, and σ, had a negligible effect on total melting time. Considering this finding, in the present work, these two parameters are not considered. Eq.14

Numerical Model
As shown in Fig.1, the outer tube and inner tube are placed concentrically. The inner tube facilitates the HTF (water) flow. Space between tubes is filled with PCM.
For a cylinder in cylinder geometry, Longeon [12] has proved through experiments that the temperature does not significantly vary with the angle (θ) in any horizontal plane. Temperature variation is significant only in radial direction(r) and axial direction(z). So, in this work, the 2D axisymmetric domain is selected for performing the simulations, which reduces the computational load substantially compared to a 3D domain. HTF is made to flow in the opposite direction to gravity, as shown in Fig.1, to conform to bottom charging.
Thermophysical properties of the RT35 PCM (Table.1) used by Longeon [12] are adopted in this numerical study.

Solution Procedure
The governing mathematical equations are built based on the following assumptions.
• PCM is pure, homogeneous, isotropic, and is in thermodynamic equilibrium in solid and liquid form.
• PCM melting is two-dimensional and axisymmetric.
• The Boussinesq assumption approximates PCM density variation.
Melting and Solidification model in CFD software Ansys Fluent is used to simulate the numerical melting. The following governing equations are considered for solving the same by employing the Enthalpy and Porosity technique [13].

Continuity Equation:
1 ( ) Continuity Equation is adopted in the form of Eq.15 for this numerical model. Since density variation is modelled with Boussinesq Approximation, the continuity equation is considered neglecting the change in density and thus the corresponding transient term. As this is a 2D domain, theta (θ) terms are also ignored.

Momentum Equations:
r -momentum equation: Eq.16 z-momentum equation: In momentum equations, θ components are nullified owing to 2D simplification. As density ( ) variation is minimal w.r.t change in temperature, it can be taken out of the partial derivative. Similarly, since the PCM flow is inviscid, the viscosity term is also considered out of the derivative.
These are Navier stokes equations, except for the terms 'Au', 'Av' and 'S b .' These source terms are added to account for the phase change and buoyancy effect.
Among all these equations, 'A' is the crucial term accounting for the phase transition. 'A' will serve two purposes here. The first purpose is to represent the solid phase by A=0. To attain zero, the (1-ε) term is taken in the numerator as the porosity(ε) is 1 in the solid phase.
'A' (Eq.18) is essentially a Carman Kozeny equation [14] with a mushy constant C*. The second purpose is to mimic the momentum equations with the Carman-Kozeny equation of porous media flow, thus quantifying phase change.
Hameter [15], Kumar [16], in their studies, have concluded that the variation in the constant C* has a significant influence on the thermohydraulic behaviour of PCM. C* is a morphological constant (or mushy zone constant) that indicates the melting front's structure and pattern. Its value is chosen as 10 6 in this work, as the recommended range is between 10 4 to 10 7 [18]. 'b' with 0.001 value is added in the denominator to avoid the 'Undefined Number' error if ε equates to zero.

= gβ∆T
S b is a buoyancy source term incorporated into the z-momentum equation. This term is obtained based on the Boussinesq approximation, where ρ ref is the reference density, g is the gravitational acceleration, β is the coefficient of thermal expansion, and ΔT is the temperature difference.

Energy Equation
Eq.17 Eq.18 Eq.19 Eq.20 Energy equation (Eq.20) is written essentially in terms of sensible enthalpy. Density is considered constant owing to incompressible flows. Thermal conductivity is deemed to be constant as the PCM is assumed to be homogeneous and isotropic.
Latent enthalpy (Δh) is blended in the Source term S h (Eq.21). This is temperature-dependent property. If PCM temperature is less than Solidus, Δh value is zero. Whenever PCM temperature crosses liquidus, Δh is equal to PCM latent heat. In greater than solidus and less than liquidus range, Δh is a fraction of L, quantified as shown below (Eq.22).
The Finite Volume Method is employed in the CFD Analysis, while the SIMPLE algorithm does the pressure and velocity coupling. PRESTO and QUICK schemes are used for discretization of pressure and momentum, respectively. Under relaxation factors are chosen as 0.3, 1, 1, 0.7, 0.9, 1 for pressure, density, body forces, momentum, liquid fraction update, energy, respectively.

Model Independence
The numerical model built is checked for its Independence. Independence is established using three criteria, i.e., Mesh size, Time Step, and the Maximum number of iterations per Time Step.

Mesh Independence
Simulations are performed with three different mesh sizes with 10,025 nodes, 18,847 nodes, and 37,647 nodes. As shown in Fig.2, 10,025 nodes mesh varied by less than 7% from the other two for the mass fraction of 0.3 to 0.75. Mass Fraction is matching with each other in 18,847 nodes mesh and 37,647 nodes throughout the calculation. The total simulation run time is doubled in the case of 37,647 nodes compared to 18,847 nodes. So, to save time, the 18,847 nodes mesh is chosen in this work. Eq.21 Eq.22

Maximum Iterations per Time Step Independence
The converged solution is obtained only when the simulation is converged at every time step. In order to get the convergence at each time step, the Maximum number of Iterations per Time Step has to be chosen appropriately. So as to verify the influence of this, Five cases are tested viz 5,10,20,30 and 50, for 18,847 nodes mesh and 0.1-time step. As Fig.4 suggests, the Mass Fraction is identical in all cases. Keeping total simulation time and other parameter variations in perspective, 20 Maximum Iterations per Time Step is chosen.

Fig.5. Dimensions (in mm) of test domain and Temperature position in space
considered for experimental validation [12].
Longeon [12] conducted experiments for Bottom charging and plotted temperature evolution in RT35 PCM. As indicated in Figure 5, the temperature at point "a" is tracked. Fig.6. shows that the numerical temperature variation is consistent with the experiment. It is to be noted that the dimensions given in Fig.5 are for experimental validation purpose only.

Results and Discussion
To understand the thermo-hydraulic behaviour of PCM in concentric shell and tube heat exchanger, dimensionless parameters viz Rayleigh Number, Stefan Number, Reynolds Number, and L/D ratio are varied by varying respective physical properties. Table.2 shows how each dimensionless parameter is quantified. In this section, the effect of each parameter is analysed, and correlation among them is presented.   7 shows the Rayleigh Number Effect on melt fraction and dimensionless time. The plot suggests that the increase in Rayleigh Number decreases the dimensionless melting time. However, the rise in Rayleigh Number will increase the actual melting time, as the amount of PCM will increase with the increase in shell diameter. The Rayleigh Number trend points to the increased influence of buoyancy force over the viscous force, which is resulting in the reduction of non-dimensional melting time. It has been observed that the increase in Rayleigh Number from 3.04x10 5 to 65.75 x10 5 decreased the non-dimensional time by 45.79% (Fig.8). Stefan Number quantifies the sensible heat to latent heat proportion in a given PCM. In this work, Stefan Number is varied from 0.2 to 1 in steps of 0.1 by changing the HTF inlet temperature. Other parameters are kept constant at Re=1500 and L/D=9. Fig.9 shows the melt fraction evolution as the Stefan Number varies. But, due to the variation in HTF inlet temperature in each case, the Rayleigh Number also varies with Stefan Number. In order to counter-balance the Rayleigh variation, the expression Fo cm .Ra n is found to have values very close to each other at n=0.1215. Thus, Fo cm .Ra n Vs. MF plot (Fig.11) gives the exclusive influence of Stefan Variation on melting.

Variation of Stefan Number
From Fig.11, it is observed that a lesser Stefan Number will result in increased time for PCM melting. This is obvious because lesser the Stefan Number, lesser the HTF inlet temperature. Moreover, the lesser the HTF inlet temperature, greater the melting time as the PCM melts only due to the hightemperature source, i.e., HTF.  Reynolds number is solely related to HTF and physically represents its velocity. However, it impacts the forced convection heat transfer between HTF and the inner pipe. Thus, the melting of PCM can be influenced by Re. In this work, Re is varied from 600 to 3000 in steps of 300 by varying the mass flow rate of HTF. St=0.5, Ra=15.32x10 5 and L/D=9 are kept constant. Fig.12 shows that the increase in Reynolds Number decreased the dimensionless melting time. MF curves for Re=1500 to Re=3000 are very close to each other, indicating decresed Re influence on melting in these Reynolds Numbers. From Fig.13, it can be understood that the percentage reduction in non-dimensional time with the increase in Re is reducing with every step. Furthermore, the maximum percentage reduction observed is 16.18%. Thus, the influence of Re is not as significant as the influence of Ra and St on the melting of PCM.
Percentage reduction from Re=900 to Re=1200 increased from 5.14% to 8.29% with 3.1% increase. But from Re=2700 to Re=3000 increase is from 15.43% to 16.18% with 0.75% increase. This shows that with the increase in Re, the percentage reduction in melting time also decreases. From Fig.14, it can be observed that with an increase in the L/D ratio, the non-dimensional melting time has increased. This is obvious because, as L/D increases with D being constant, the volume of PCM increases. As PCM volume increases, the time for melting it also increases. L/D is different from the other three parameters. While L/D increase resulted in Fo cm increase, increase in Ra, St, Re resulted in Fo cm decrease. Fig.15 shows that a maximum increase of 92.98% is recorded when L/D raised from 2 to 15. In the same bar chart, from L/D = 3 to L/D = 9, the present case to previous case increase is greater than 6%. After L/D=9, the case after case increase is less than 6%. ( / ) Multi Linear Regression is performed to obtain correlations in the form of Eq23. As the plots (Fig.7, Fig.9, Fig.11, Fig.12, Fig.14) suggest, the variation of MF till 0.8 is following one particular trend. Whereas, after 0.8, till 1, the trend is completely changed. Based on this, for MF less than 0.8, one correlation is obtained, and for MF greater than 0.8, a similar correlation with the same variables and different values in the power of those variables is obtained (Table.3).

Eq.23
St* is modified Stefan number, which takes the initial temperature of PCM into consideration. To include this effect, correlations are also considered with St* replacing St. As Table. Fig.16 shows the comparison of a few numerical cases with correlation. It points to the good agreement of Correlation with the Numerical data. The R square values (Table.3) also confirm that the variation is minimal.  Eq.24 is the correlation for predicting the non-dimensional melting time, where G is the constant value and H, J, K, L are power values. Table.4 Eq.24 shows the values in the case of St and St*. Fig.17 indicates that the Fo cm computed values are close to the numerically obtained Fo cm values. R square values further validates the same.
It is to be noted that the Fo cm correlation is simpler compared to the MF correlation. This is expected because MF correlation has considered additional continuous variable MF. On the other hand, Fo cm correlation considered that MF is constant and equal to one. A similar tendency is evident with the correlations in the first section of this work.

Conclusions:
Parametric Analysis is conducted on PCM based concentric tube heat exchanger. A validated numerical model is used to predict PCM's melt fraction and dimensionless melting time.
Rayleigh number is varied from 3.04x10 5 to 65.75 x10 5 . A maximum decrease of 45.79% is observed in non-dimensional melting time (Fourier Number). Ra trend pointed to the increased influence of Buoyancy.
Stefan Number is varied from 0.2 to 1. "Fo cm .Ra 0.1215 " is used to neutralize the effect of simultaneous Rayleigh variation. While 78.22% of percentage reduction in Fourier Number is recorded, 73.54% is the exclusive effect of Stefan Number, and the remaining 4.68% is attributed to Rayleigh Number. Thus, Stefan Number is established as the primary factor influencing the melting time.
Reynolds Number is varied from 600 to 3000. This decreased the melting time only by 16.18%. Hence the influence of Re is not as great as the influence of Rayleigh Number and Stefan Number on the PCM melting.
L/D ratio is varied from 2 to 15. Among the parameters investigated, L/D is the only parameter that has witnessed an increase in the melting time with its rise.
The proposed correlations predictions are convincing for quantifying Melt Fraction and dimensionless melting time.