Numerical investigation of the in ﬂ uence of mushy zone parameter A mush on heat transfer characteristics in vertically and horizontally oriented thermal energy storage systems

(cid:129) Numerical simulations were conducted to examine the e ﬀ ect of A mush on predicted heat transfer and melting of a PCM. (cid:129) The smallest values of A mush resulted in unrealistic predictions of the melt front development. (cid:129) The higher values of A mush corresponded to delayed melting of the PCM. (cid:129) The impact of A mush is generally less prominent in regions which conductive heat transfer dominates. The e ﬀ ect of the value used for the mushy zone parameter (A mush ) on predicted heat transfer and melting characteristics of a phase change material (PCM) Lauric acid, in both vertical and horizontal enclosures was studied. There is a lack of clarity regarding which value of this parameter should be used for accurate simula- tions of phase change heat transfer, addressing this will aid in accurate simulation and design of systems for LHTES (Latent heat thermal energy storage). The numerical analysis undertaken used a commercial CFD code ANSYS FLUENT 18.2 and the enthalpy-porosity formulation. The range of mushy zone parameter used was from 10 5 to10 7 . The predicted locations of the melt front were compared to published experimental data available in the literature. The simulations provided quantitative information about the amount of energy stored and the melt fraction and providing improved understanding of the heat transfer process. Comparison between predictions using di ﬀ erent values of A mush , and experimental data showed that correct selection of the value of A mush to be used in the momentum equations is an important parameter for accurate modelling of LHTES and has a sig- ni ﬁ cant in ﬂ uence on the solid-liquid interface shape and progression. The study reveals that increasing the value of A mush leads to a decrease in ﬂ uid velocity, decreasing convection and the rate of heat transfer, therefore, proper selection of the mushy zone parameter is necessary to accurately simulate LHTES systems and provide a better understanding of the phase change behaviour and heat transfer characteristics.


Introduction
With the diminishing reserves of easily accessible fossil fuels, the increasing global energy demand, increasing levels of greenhouse gas emissions, more energy-efficient environmentally friendly devices are essential to reach the "3E" objectives (Economic, environmental and energy security), and deliver a clean environment, a sustainable energy policy, a solid economy, and social development [1]. Latent heat thermal energy storage systems (LHTES) using phase change materials [3].
To effectively design a latent heat storage system for a specific process, the time required for melting and solidification and thus rates of charge and discharge and heat transfer during the phase change process has to be known [4]. Any numerical models used to predict PCM system performance should realistically simulate the processes of melting and solidification and be computationally efficient allowing the accurate simulation of LHTES systems within affordable computational resources. Understanding the transient characteristics of the heat transfer processes and the limitations of different modelling approaches for the phase transition process is important for the effective design, evaluation and optimization of LHTES systems [5].
Numerical modelling of thermal energy storage (TES) systems has attracted considerable attention recently because of the increased attractiveness of TES for a wide range of different applications. Although appropriate numerical procedures can vary widely depending on the system, the impetus for their use can almost always be attributed to cost and time constraints [6]. For example, experimental investigations to assess performance of TES systems can be expensive and more difficult to undertake compared to numerical modelling, cases considering high temperature applications or using hazardous materials or operational conditions, may be simulated more easily than investigated experimentally.
Analytical models can be subject to extensive simplifying assumptions that render models incomplete but are necessary since the exclusion of such simplifying assumptions can make equations overly complicated and cumbersome to evaluate. This has resulted in an increasing requirement for numerical models to accurately describe TES system behaviour, allowing computational methods to help solve governing equations. As a result, much of the research regarding this subject reflects the advancement of appropriate numerical analysis.
To simulate the phase change process involves a moving boundary between phases and is generally known as a Stefan problem. Melting and solidification are the two processes generally involved in latent heat storage systems. When heat is transferred to a solid PCM and melting occurs latent heat is stored. When heat is transferred from the PCM and solidification occurs latent heat is released. The PCM phase transition proceeds from a solid to a mushy state and then to a liquid state during the melting process and vice versa during the solidification process. The heat transfer mechanisms associated with the melting and solidification processes can be either conduction or convection (natural convection in most situations) controlled or simultaneous conduction/ convection controlled [7].
To maximize the performance of thermal storage systems using PCMs, knowledge of the thermal behaviour of the PCM that is employed is required, several research studies haves concentrated on the detailed thermal behaviour of PCMs during the phase change process [8]. Different modelling approaches for solid-liquid phase-change have been developed in recent decades. The enthalpy-porosity approach can adequately describe the natural convection effect in the melt region [9]. This approach is capable of accurately predicting both the transient position and shape of the melt front at different times with relatively modest computational requirements. This approach uses a fixed-grid for the analysis of solidifying and melting of materials. The enthalpyporosity formulation currently serves as the solidification and melting model within ANSYS FLUENT and COMSOL Multiphysics solvers [10]. This technique assumes phase change to occur over a finite temperature range; thus, generating an artificial mushy region in which the melt fraction of a fluid element varies from zero (solid phase) to 1 (liquid phase). Naturally, the velocity of the fluid element within the mushy region should also vary from zero (solid phase) to the natural convection velocity (liquid phase). The enthalpy-porosity formulation deals with this velocity transition by modelling flow within the mushy region as flow through a porous medium. A sink term, in the form of the Carman-Koseny equation, is added to the Navier-Stokes equations to mimic the effect of damping within the mushy region.
The mushy zone is considered to be a semi-solid existing as an interface between the melted and un-melted region of a PCM undergoing melting or freezing. The mushy zone parameter A mush measures the amplitude of the damping; the higher this value, the steeper the transition of the velocity of the material to zero as it solidifies. Very large values of A mush may cause the predicted solution to oscillate [11]. This region significantly influences the heat transfer and flow characteristics during the melting and solidification process. Simulating fluid flow resulting from natural convection within the mushy zone has been the subject of numerous discussions in recent years. Although the Carman-Koseny equation is used in the majority of available models, the effect of the A mush parameter is still the subject of frequent inquiry [12]. Values for A mush ranging from 10 3 to 10 8 have been suggested in commercial software guidelines and by different researchers. The correct prediction of the heat transfer behaviour within the mushy zone has been the subject of numerous discussions in recent literature on phase change heat transfer. Table 1 summarises previous research examining the effect of varying the value of A mush on the overall simulation results for melting and solidification of PCM in different geometries. The values to use are still the subject of discussion since recent studies reported in the literature observe a significant difference between experimentally determined behaviour and that predicted by modelling depending on geometry and/or PCM material used.
The objective of the present study is to evaluate the effect different values of A mush have on the overall simulation of PCM melting within vertically and horizontally oriented storage geometries and its effect on the melt fraction, vortex strength and amount of heat stored. The model predictions with different values of the A mush parameter are compared with the experimental data published by [8].

CFD simulation
CFD simulation was carried out using the commercial software ANSYS FLUENT 18.2. The physical model and computational procedure are discussed below.

Physical model
A schematic diagram of the two-dimensional simulation geometry is shown in Fig. 1. The dimensions of the container containing Lauric acid used in the simulations are based on those in the experimental study carried out by Kamkari B et al [36]. In the experimental study PCM melting was quantified by visually tracking the shape of the melting interface and how it changed with time. The container had a rectangular cross section with inside dimensions of 0.05 m in width, 0.12 m in height and 0.12 m in depth, the wall of the enclosure indicated in Fig. 1 was held at constant temperature 70°C (343.15 K). The other three walls of the enclosure were made from 0.025 m thick plexiglass with a thermal conductivity k of 0.043 W/m K.
The PCM used in the experiments was Lauric acid of 99% purity. Table 2 presents its thermophysical properties. The entire system was initially at a constant temperature is 298.15 K.

Computational approach
The numerical approach adopted makes it possible to predict the natural convection in the liquid PCM during the melt process that  occurs inside the container. In the approach used to simulate melting, the flow was considered to be unsteady, laminar, incompressible and two-dimensional. It was assumed that both solid and liquid phases are homogeneous and isotropic. The enthalpy-porosity approach [9] was adopted for simulation of the PCM. The governing conservation equations used for the PCM system are: Continuity Energy where ρ is the density, k is the thermal conductivity, μ is the dynamic viscosity, S is the momentum source term, V is the fluid velocity, T is the temperature, P is the pressure, g is the gravitational acceleration and H is the specific enthalpy. The latter is defined as a sum of the sensible enthalpy (h) and the latent heat (ΔH). where: where h ref is the reference enthalpy at the reference temperature T ref , C p is the specific heat. The content of the latent heat can vary between zero (for a solid) and L (for a liquid): where the liquid fraction (β) can be written by the following relations: The source term Sin the momentum equation, Eq.
(2), is given by: where A(β) is the ''porosity function" defined by Brent et al. [9]. The source term was used to describe the flow in the porous medium in the momentum equation, and it has to be zero in the liquid phase to allow for free motion, but it has to be large in the solid phase to force the velocity values to near zero values [13], while different functions fulfil this requirement, most often the Carman-Kozeny equation, which is derived from the Darcy law for fluid flow in porous media, is used in a modified form: where ɛ is a small computational constant used to prevent zero in the denominator (in this work ɛ = 0.001), and A mush is the mushy zone constant which measures how fast the fluid velocity approaches zero as it solidifies. The intensity of natural convection during melting is quantified by the Rayleigh (Ra) number, defined based on the characteristic length H, which is calculated by; where T w is the wall temperature, T m is the PCM average melting temperature and H is the height of the domain. For the present simulations the calculated Rayleigh (Ra) number varies from 3.51 × 10 6 to 2.54 × 10 7 .

Computational procedure
The simulations have been performed using the ANSYS FLUENT 18.2 software, run as two-dimensional double precision (2ddp) code. The pressure based coupled algorithm was used to solve the momentum and continuity equations. The gravity vector was set to −9.8 m/s 2 in the y-direction to allow prediction of the natural convection in the PCM during the simulation. A second order upwind scheme for the advection term, central differencing for the diffusion term and a second order implicit discretization scheme for the transient term were used. The PRESTO pressure interpolation scheme for the transient calculations, was implemented. The under-relaxation factors for density, momentum, pressure correction, thermal energy and melt fraction used were 1, 0.7, 0.3, 1 and 0.9, respectively. The mesh was generated using the mesh generator ICEM CFD.
Different grid resolutions were tested by comparison of the predicted melt fraction with time to ensure independency of the solution from the adopted grid density. The mesh used comprised of mapped quadrilateral cells in the Insulation, Plexiglas and PCM domains. The effects of time step size and grid size on the solution were carefully examined in preliminary calculations The mesh independence test was performed using mapped grids of 25 × 80, 50 × 120 and 75 × 180 cells within the PCM domain as shown in Fig. 2, from Fig. 3 it can be seen that the higher mesh density resulted in no significant change in the prediction of the average melt fraction development with time, mesh 2, 50 × 120 elements was deemed satisfactory for assessing the effect of varying A mush on the overall melt front behaviour. The maximum number of iterations for each time step was fixed at 300 and was found to be sufficient to satisfy the convergence criteria of 10 −6 for the velocity components and continuity, and 10 −11 for the energy equations. The time required to achieve full melting is a good measure of time step dependence. By comparing the melt fraction obtained with simulations using time steps of 0.1, 0.2 and 0.4 s, the time step selected for integrating the time derivatives was set to 0.2 s. In the earlier stages of the simulation that were dominated by conduction, faster convergence occurred compared to latter stages with a significant melt fraction that were dominated by natural convection.
Simulations were conducted on the Loughborough University Research High Performance Computing (HPC) cluster which consists of 7 compute nodes, each having two six-core Intel Westmere Xeon X5650 CPUs and 24 GB of memory. A typical simulation to achieve complete melting required from 30 h to 50 h of computing.

Numerical analysis and comparison with experiments
The numerical model used in this study was an approximate representation of the storage system presented by Kamkari and Shokouhmand [2], [8] and [36]. From the 3D physical geometry, a twodimensional mid-plane was simulated. It was assumed that the boundary effect at the end walls in the third dimension (z-direction), have a negligible effet on the simulation domain in the mid-plane. In this study, the mushy zone constant was varied from 10 5 to 10 7 . Fig. 4 shows the predicted melt fractions for vertical and horizontal rectangular enclosures with the different values of A mush used. For A mush = 10 7 the melting rate is very slow, for A mush = 10 5 the melting rate increases significantly, being greater than that observed experimentally. For the vertical orientation, an A mush value of 5 × 10 5 was found to give the best agreement with the experimental results and the most physically realistic results. In comparison, for the horizontal orientation for a value of A mush of 2 × 10 5 simulation results agreed well with the experimental results. During melting, as previously stated, the heat transfer is strongly influenced by convection. An increase in the value of A mush leads to a decrease in the predicted level of convection with a consequent reduction in the melting rate [37].
Additional validation was performed by comparing the predictions from the model to those of the experimental benchmark data of [8] for T w = 55 and 60°C using an A mush value of 5 × 10 5 and is shown in Fig. 5. Fig. 6 shows the experimentally measured solid liquid interface boundaries during the melt process for the vertically oriented enclosure at elapsed times of 10,20,30,40,60, and 80 min and those predicted for four mushy zone values, 10 5 , 5 × 10 5 , 10 6 and 10 7 . From this Figure, it can be seen that, the effect of A mush on the melting rate of the PCM is now clearly more pronounced and noticeable especially in the later stages of melting.
The resulting shape of the solid PCM for A mush = 5 × 10 5 does match the one obtained in the experiments well, especially towards the end of the process. This Figure shows that melt interface positions best correlate with the experimental results at the early stages of the melt process when the heat transfer is dominated by conduction (at an elapsed time of about 20 min or less).
Predictions using an A mush value of 10 5 show accelerated melting in the upper part of the enclosure and that the melt interface progresses faster than that observed in the experimental results, while an A mush value of 10 7 results in predictions in which the liquid/solid interface lags behind the experimental results and leads to a sharper curvature along the melt interface.
In conclusion, the present predictions show that higher A mush values result in a shaper curvature along the melt interface, with an A mush value of 5 × 10 5 , the melt interface curvature more accurately resembles that of the experimental results, this further supports the statement that an A mush value of 5 × 10 5 is the optimal value out of the five values tested for simulating this PCM system. Fig. 7 shows the experimental and the predicted solid liquid interface for the horizontally oriented enclosure for 4 different values of A mush, similar trends are observed, the effect of different values of A mush are less obvious in the early stages of the simulation where heat transfer is dominated by conduction. Larger values of A mush decrease the rate of melting in all simulated cases. The locations of peaks and troughs along the melt front vary in a random fashion (expect adjacent to the right and left side walls where wall effects dominate); From Fig. 7, it can be seen that the average height of the melt interface over time is affected marginally by A mush , the most likely reason for this behaviour is the reduced magnitude of the velocity at or near the liquid/solid interface in the horizontally oriented enclosure. Lower magnitudes of velocity reduce the damping forces within the momentum conservation equations (which, at a given melt fraction, is proportional to the fluid velocity, with A mush the proportionality factor). At the early stage, 10 min, of the simulated melting from the bottom surface in the horizontal enclosure, it can be seen that the solid-liquid interface is nearly flat and heat conduction is the dominant mode of heat transfer, after this time, the interface shape becomes increasingly undulating indicating that convection currents in the liquid PCM are now determining the interface shape. The undulating shape of the interface is a result of the formation of Bernard cell convection in the liquid PCM.
The predictions of an undulating interface between the solid and liquid PCM provides excellent trend wise agreement with the experimental results. The peaks and troughs along the melt interface do not align perfectly when comparing experimental observations and numerical simulations; however, the height of the interface is accurately predicated for the first 40 min of the simulation.
The predicted temperature distributions in the vertically and horizontally oriented enclosures for different value of A mush , are presented in Figs. 8 and 9 corresponding to the time periods presented in Figs. 6 and 7. Fig. 8 presents temperature contour plots in the vertically oriented container, it can be seen that in the first 10 min during early stages of melting conduction plays a significant role, and the value of A mush does not have any significant effect on predicted temperature contours which are in all cases nearly parallel to the hot wall. After 10 min the temperature and volume of the liquid PCM adjacent to the hot wall increases sufficiently so that buoyancy forces overcome the viscous    forces and hot liquid PCM rises along the vertical hot wall. It is clear from this Figure that the effect of A mush becomes more significant as convection becomes the dominant heat transfer mechanism. This appears to be caused by predicted PCM flow within the mushy zone that is greater than that occurring in the experiments, resulting in accelerated rates of heat transfer and melting being predicted. Fig. 9 presents predicted temperature contours in the horizontally oriented enclosure. It can be seen that, at 10 min the predicted temperature profiles are similar for all values of A mush the undulating temperature isotherms at the solid/liquid interface results due to the vertical convective flow structure in the liquid PCM, this is responsible for efficient heat transfer from the hot wall to the solid-liquid interface and effectively mixes the liquid PCM. In the later stage of melting after 60 min, for different values of A mush , the predictions are different and there is a greater dependence on the value of A mush over the entire range of values used.
Because convective flow in the liquid PCM affects the shape of the melt front, velocity contours and velocity vector maps during melting in the liquid PCM are presented in Figs. 9 and 10 for both vertical and horizontal arrangements at a simulated time of 20 min.
In the vertically oriented container, higher velocities are located next to the wall, the liquid PCM has a higher temperature and lower density than the adjacent material which causes the hot liquid to move upwards. The already melted Lauric acid in contact with the still solid PCM is cooler than the surrounding liquid PCM, which has a greater density, causing this liquid to move downward. Thus, the differences in density lead to rising flow near the heated vertical boundaries, as expected, for the vertical container the liquid PCM velocity along the solid/liquid interface increases as the A mush value decreases due to enhanced buoyancy-driven convection as shown in Fig. 10.
The predicted flow for the horizontal container at 20 min is presented in Fig. 11, cellular flow is predicted at different positions along the entire length of the PCM domain. Most of the flow cells are in pairs and rotate in the opposite directions, this is because the heated liquid, which is less dense rises transferring heat to the solid PCM by free convection.
The extent of the fluid flow is restricted by the solid PCM above it, the liquid produced by melting the solid PCM at the solid liquid interface cannot rise any higher because of the solid PCM and sinks towards the lower heated surface.
The predictions of stored energy with time for the vertical and horizontal enclosures using different A mush values during the charging process are shown in Figs. 12 and 13. It can be seen that, differences in the predicted amount of the energy stored with time are significant for different values of A mush .
During the melting process, energy transfer to the PCM in the horizontal container is faster than to that in the vertical container. The figures also show that the A mush values have a significant impact on the rate of heat transfer and the energy stored at different times. This is due to the effect of A mush on the predicted natural convection process, at a charging time of 150 min, for the vertical enclosure, it can be seen that the difference between the predicted amounts of energy stored for A mush = 10 5 and 10 7 is up to 37.5% as shown in Fig. 12. For the horizontal enclosure, at a charging time of 100 min, the difference between the predicted amounts of energy stored for A mush = 10 5 and 10 7 is up to 28.5% as shown in Fig. 13.
The overall effect of A mush on the melting rate was as expected, the Carman-Koseny equation acts as a damping term in the momentum equations, large values of A mush increase the volume within the mushy region where the fluid PCM) is static; this reduces natural convection and effectively decreases the overall heat transfer within the mushy region. The CFD predictions indicate that proper selection of the mushy zone constant is essential for accurate modelling of phase change heat transfer.

Concluding remarks
In the present work, numerical simulations were conducted to examine the effect of A mush (mushy zone parameter) on predicted heat transfer (conduction and natural convection) and melting of a PCM in vertically and horizontally oriented thermal storage enclosures heated from one side. The simulations were performed using ANSYS FLUENT 18.2, for different values of the mushy zone parameter ranging from 10 5 to 10 7 . The numerical predictions were compared with experimental work carried out by Shokouhmand and Kankari [2].
Based on this study, it was observed that A mush is an important parameter for accurately modelling phase change phenomena. The smallest values of A mush (< 10 5 ) used resulted in unrealistic predictions of the melt front development, higher values of A mush (10 6 ) resulted in  delayed prediction of melting in the PCM. From this work, it is clear that the impact of A mush is generally less prominent in regions which conductive heat transfer dominates. The effect of A mush was more pronounced farther away from the heated surface when more PCM had melted and heat transfer was dominated by natural convection. This was because with increasing A mush , predicted convective strength was observed to decrease, which decreases the heat transfer rate. Proper selection of the mushy zone constant is essential for the accurate prediction of heat transfer characteristics within a PCM. A further similar study is required for different PCMs and experimental work to determine the extent to which correct selection of A mush can be correlated to material and store properties.