Numerical Study of Solidification in a Plate Heat Exchange Device with a Zigzag Configuration Containing Multiple Phase-Change-Materials

Latent heat thermal energy storage (TES) plays an important role in the advocation of TES in contrast to sensible energy storage because of the large storage energy densities per unit mass/volume possible at a nearly constant thermal energy. In the current study, a heat exchange device with a zigzag configuration containing multiple phase-change-materials (m-PCMs) was considered, and an experimental system was built to validate the model for a single PCM. A two-dimensional numerical model was developed using the ANSYS Fluent 14.0 software program. The energy fractions method was put forward to calculate the average Ste number and the influence of Re and Ste numbers on the discharge process were studied. The influence of phase change temperature among m-PCMs on the solidification process has also been studied. A new boundary condition was defined to determine the combined effect of the Re and Ste numbers on the discharging process. The modelling results show that for a given input power, the Ste (or Re) number has a significant impact on the discharging process; however, the period value of inlet velocity has almost no impact on it. Besides, the zigzag plate with m-PCMs has a good impact on the temperature shock as “filter action” in the discharging process.


Introduction
In recent decades, environmental pollution and energy shortages have had a significant impact on our lives, and addressing these problems is a great challenge for scientific researchers.Thermal energy storage (TES) plays a key role in addressing environmental pollution and energy shortages.TES systems can be divided into sensible heat thermal storage (SHTS) and latent heat thermal storage (LHTS) systems [1].However, compared with SHTS, LHTS is a particularly attractive technique because it provides a high energy storage density and a constant or near-constant operating temperature that corresponds to the phase change temperature of phase-change-materials (PCMs).Thus, LHTS systems using PCMs are one of the most attractive methods to store thermal energy in different engineering fields, such as electronic products [2,3], waste heat recovery [4,5], room temperature regulation [6][7][8], solar energy utilization [9] and off-peak electricity tariff implementation [10,11].Nevertheless, the PCMs loaded in the unit possess a relatively low thermal conductivity that drastically suppresses the energy charging/discharging rates of LHTS systems.There are several common ways to enhance the heat transfer of the LHTS systems, including improving the thermal conductivity of PCMs (using metallic foam [12], carbon nanotube (CNT) [13]), using m-PCMs [14,15], changing LHTS unit configurations (adding fins [16,17], using heat pipes [18,19], using enhanced tubes [20], using eccentric tube [21], etc.) and synthesizing those three enhancements [22,23].
Gong and Mujumdar [14] numerically studied the cyclic melting-freezing of m-PCMs in a series arrangement in an LHTS slab unit on the fixed wall temperature boundary condition.The numerical results indicate that the instantaneous heat flux can be greatly enhanced using m-PCMs slabs.Ismail et al. [24] studied the influence of the fin number, fin length, fin thickness and annular spacing on the freezing process; they found that the fin length and annular spacing, but not the fin thickness, had significant influence on the complete solidification time on the constant wall temperature condition.Blen et al. [25] studied the effect of the fin number on the melting and solidification process of CaCl 2 ¨6H 2 O as a PCM in a vertical two-concentric shell and tube unit with the fins placed inside the PCM, and they also studied the mass flow rate and inlet temperature on the melting and solidification process.It has been reported that the effect of the fins placed inside the PCM to the melting and solidification time is much more than the effects obtained by the flow rate.
Ismail et al. [26] experimentally and numerically studied the effect of different operation parameters, such as the Dean number (De " Re{ a R{r), fluid flow rate and inlet temperature of the heat transfer fluid (HTF), on the solidification of PCM around a curved cold tube.Al-Abidi et al. [27,28] numerically studied the solidification of PCM in a triplex tube heat exchanger with internal and external fins on the boundary condition of a constant wall temperature using a 2D numerical model.They reported that the case with full-size fins achieved complete solidification in a short time.Al-Abidi et al. [27] also experimentally studied the effect of the flow rate on the solidification of PCM in a triplex tube heat exchanger with fins.Mosaffa et al. [29] studied the effect of the design parameters of thermal storage units, such as the slab length, thickness and fluid passage gap, on the storage performance in the constant inlet temperature condition.The optimum parameters were determined to be air channel thickness of 3.2 mm, slab length of 1.3 m and slab thickness of 10 mm.
As discussed above, there have been very few reports about the solidification process of m-PCMs within a zigzag plate unit.The charging or discharging process of LHTS often use constant wall temperature or wall thermal flux or by heat transfer between HTF and PCM.However, the work conditions of simulation or experiment were focused on the inlet velocity (mass flow rates) or temperature boundaries with one of them fixed while the other varied.In this article, the arrangement of m-PCMs (reversed and forward) within a zigzag plate unit was studied, and different work conditions, such as the inlet temperature (Ste) and inlet velocity (Re), were also studied and the energy fractions was defined to calculate Ste number for the m-PCMs in arbitrary ratio.In addition, a new work condition, inlet power, was defined to determine the synergistic effect of Ste and Re on the discharging process of a zigzag plate in LHTS systems and the effect of "filter action" on the the temperature shock was also studied.

Physical Model
The physical configuration of the zigzag plate is shown in Figure 1.With reference to this figure, the geometries of the zigzag plate are given as: H 1 = H 2 = 2 mm, H 3 = H 4 = 4 mm and H 5 = 2 mm.The thickness of the plate is 1 mm.The PCM considered in this work is a salt mixture of sodium chloride and magnesium chloride at a mass ratio of 4:6.The thermophysical properties of the salt mixture are tabulated in Table 1.The arrangement of m-PCMs in the forward and reversed directions is shown in Figure 2. Figure 2a shows the three PCMs arrangement in reversed direction which means the phase change temperature of the three PCMs increases gradually in the HTF flow direction; Figure 2b shows the three PCMs arrangement in forward direction which means the phase change temperature of the three PCMs decreases gradually in the HTF flow direction.

Governing Equation
The mathematical formulation of the phase change problem includes mass, momentum and energy balance equations [30].The energy equation takes the following form: where H, ρ, , k, T and Sh are the enthalpy, density, velocity vector, thermal conductivity, temperature and energy source term, respectively.The enthalpy, H, includes both sensible enthalpy (h) and latent enthalpy (ΔH): with the sensible enthalpy given by: d where href, Tref, and cp stand for the enthalpy at a reference temperature, the reference temperature and the specific heat at a constant pressure, respectively.The latent enthalpy in Equation ( 2) is related to the fraction of liquid phase (β) during a phase change and the specific latent heat of the PCM (L) by:  The arrangement of m-PCMs in the forward and reversed directions is shown in Figure 2.  The arrangement of m-PCMs in the forward and reversed directions is shown in Figure 2. Figure 2a shows the three PCMs arrangement in reversed direction which means the phase change temperature of the three PCMs increases gradually in the HTF flow direction; Figure 2b shows the three PCMs arrangement in forward direction which means the phase change temperature of the three PCMs decreases gradually in the HTF flow direction.

Governing Equation
The mathematical formulation of the phase change problem includes mass, momentum and energy balance equations [30].The energy equation takes the following form: where H, ρ, , k, T and Sh are the enthalpy, density, velocity vector, thermal conductivity, temperature and energy source term, respectively.The enthalpy, H, includes both sensible enthalpy (h) and latent enthalpy (ΔH): with the sensible enthalpy given by: d where href, Tref, and cp stand for the enthalpy at a reference temperature, the reference temperature and the specific heat at a constant pressure, respectively.The latent enthalpy in Equation ( 2) is related to the fraction of liquid phase (β) during a phase change and the specific latent heat of the PCM (L) by:

Governing Equation
The mathematical formulation of the phase change problem includes mass, momentum and energy balance equations [30].The energy equation takes the following form: where H, ρ, Ñ v , k, T and S h are the enthalpy, density, velocity vector, thermal conductivity, temperature and energy source term, respectively.The enthalpy, H, includes both sensible enthalpy (h) and latent enthalpy (∆H): with the sensible enthalpy given by: where h ref , T ref , and c p stand for the enthalpy at a reference temperature, the reference temperature and the specific heat at a constant pressure, respectively.The latent enthalpy in Equation ( 2) is related to the fraction of liquid phase (β) during a phase change and the specific latent heat of the PCM (L) by: with β given by the following equation: where T s and T l stand for the solidification point and the melting point, respectively.This method is known as the enthalpy-porosity model, which treats the phase-changing region (mushy region) as porous, with a porosity equal to the liquid fraction.The source term of the energy equation is given as: Because of the formation of liquid during a phase change, fluid flow occurs in the liquid phase and the mushy zones, which are governed by the following set of momentum equations: B pρuq Bt `∇¨´ρ where u and v are the xand y-direction velocity, respectively; P is the pressure; µ is the viscosity; and S x , S y and S b stand for the source terms of the phase change zone in the x-direction, the y-direction and that resulting from gravitational acceleration (y-direction), respectively, and are given in the following equations: S y " where g is the gravitational acceleration, ε is a small number (e.g., 0.001) to prevent division by zero, and A is the mushy zone constant with a value depending on the structure of the porous mushy zone but normally assumed to be 10 ´5.The phase change can induce turbulent flows, and as a result the following turbulent sink term needs to be considered: where Φ stands for the turbulence parameter obtained from the turbulent flow model.The large eddy simulation (LES) model is used because of the complexity of the flow field.

Boundary Conditions
At the initial time, the m-PCMs were superheated in a liquid state and the entire calculation domain was set to a zero velocity (that was higher than the highest melting point of m-PCMs).Periodic boundary conditions were assumed for both the upper and lower computational domain.Adiabatic boundary conditions were applied to the front and back walls.No-slip boundary conditions were applied to the interior zigzag tooth-like wall surfaces and the inlet and outlet boundary condition were set as velocity inlet and outflow respectively in the FLUENT (a commercial software of computational fluid dynamics) environment.The boundary conditions can be pictorial description as follows: T in = T (constant or changed with time depending on the work condition); v in = v (constant or changed with time depending on the work condition); T o = 803.15K.

Numerical Modelling
The models were solved in the FLUENT environment.The Semi-Implicit Method for Pressure-Linked Equations (SIMPLE) algorithm was used for velocity-pressure coupling.The PREssure Staggering Option (PRESTO) scheme was used for pressure spatial discretization, whereas the bounded central differencing (BCD) scheme was used for the discretization of the momentum terms.The under relaxation value factors for the pressure, velocity, energy, and liquid fraction were 0.3, 0.2, 1, and 0.9, respectively.The computational domain was represented by a two-dimensional base grid, which had 49,163 mesh cells, with 19,560 cells in the PCM domain and 13,770 cells in the HTF domain.A finer grid was created by increasing the mesh density near the plate walls.Grid independence was studied using grid sizes of 0.3 mm, 0.4 mm and 0.5 mm, which resulted in 84,932, 63,410 and 49,163 total numbers of cells, respectively.There was little effect of the grid size on the time evolution of the liquid fraction.It was therefore decided to use a grid size of 0.5 mm with 49163 cells in all of the subsequent simulations.The independent time steps and grid sizes for the melt fraction for a zigzag plate, as well as the different time steps of 0.005 s, 0.01 s (chosen), and 0.02 s, were examined prior to simulation.

Experimental System and Procedure
Experimental work was carried out to verify the model.Figure 3 shows a schematic diagram of the experimental system, which consisted mainly of an electric heating furnace, a high temperature blower, and two phase change heat exchangers.In the charge process, air was heated to a preset temperature by the electric heating furnace and driven to circulate in the loop by a high temperature blower.In the charge step, the heated air flowed through the phase changer heat exchangers where heat was transferred to and stored within the phase change material.When fully charged (judged from temperature measurements of the outlet air temperature), the electric furnace and the high temperature blower were switched off to end the charging step.In the discharge step, the ambient air blower was switched on to deliver a preset rate of air flow to the charged phase change heat exchanger.Air was used as the HTF because of the high melting points of the phase change materials considered in this work.The two heat exchangers could be charged at the same time.Through the use of a combination of valves, it was also possible to charge one of the heat exchangers while the other was discharged.In this work, only one phase change heat exchangers was used, and the charge-discharge processes occurred cyclically.
air blower was switched on to deliver a preset rate of air flow to the charged phase change heat exchanger.Air was used as the HTF because of the high melting points of the phase change materials considered in this work.The two heat exchangers could be charged at the same time.Through the use of a combination of valves, it was also possible to charge one of the heat exchangers while the other was discharged.In this work, only one phase change heat exchangers was used, and the charge-discharge processes occurred cyclically.

Verification and Validation of the Simulation Models
Figure 5 compares the measured and calculated time evolutions of the temperature at the outlet and the centre of PCM (as indicated in Figure 4a).It can be seen that reasonably good agreement has been achieved, thus establishing confidence in the model.The deviation of thermocouple 1 is likely associated with the existence of pores in the PCM because the material was loaded into the plate heat exchangers in compacted powder form.Other reasons may include the temperature dependence of the PCM properties, an uneven distribution of the HTF flow across all of the plates, which were not considered in the modelling, or the location of the thermocouple.4a).It can be seen that reasonably good agreement has been achieved, thus establishing confidence in the model.The deviation of thermocouple 1 is likely associated with the existence of pores in the PCM because the material was loaded into the plate heat exchangers in compacted powder form.Other reasons may include the temperature dependence of the PCM properties, an uneven distribution of the HTF flow across all of the plates, which were not considered in the modelling, or the location of the thermocouple.

Verification and Validation of the Simulation Models
outlet and the centre of PCM (as indicated in Figure 4a).It can be seen that reasonably good agreement has been achieved, thus establishing confidence in the model.The deviation of thermocouple 1 is likely associated with the existence of pores in the PCM because the material was loaded into the plate heat exchangers in compacted powder form.Other reasons may include the temperature dependence of the PCM properties, an uneven distribution of the HTF flow across all of the plates, which were not considered in the modelling, or the location of the thermocouple.

Modelling Results and Discussion
Numerical simulations were performed on three PCMs that were contained between zigzag plates.The ranking of the phase change temperatures followed PCM-1 (T PCM-1 ) > PCM-2 (T PCM-2 ) > PCM-3 (T PCM-3 ).It is assumed that T PCM-1 ´TPCM-2 = T PCM-2 ´TPCM-3 = ∆T/2, where T PCM-2 = (T PCM-1 + T PCM-3 )/2 = 713.15K and ∆T = T PCM-1 ´TPCM-3 .Four sets of m-PCMs with different phase change temperatures were selected for the modelling study with ∆T = 20, 40, 60, and 80 K.The results were compared with ∆T = 0 K (a single PCM only).Unless otherwise indicated, the inlet velocity of the HTF and the mass ratio of m-PCMs were assumed to be unity.Dimensionless time:

Multiple Phase-Change-Materials in Different Arrangements
where Ste and Fo are the Stefan and Fourier numbers, respectively; α is the PCM thermal diffusivity given by α = k/(ρc p ), with k representing the PCM thermal conductivity; T in and T o are the inlet temperature of HTF and initial temperature of PCMs respectively; and t is elapsed time.Regarding m-PCMs units, Ste is changed with Ste (ave) which is calculated as: where ω 1 , ω 2 , and ω 3 are the energy fractions of PCM-1, PCM-2 and PCM-3, respectively.They are defined as follows: Because all of the physical properties except the phase change temperatures are the same, the energy fractions are equal to the mass fraction in this study.However, except for the phase change temperature, if other physicals of PCMs are the same, the Ste (ave) number can be calculated in an easier way [31].
Figure 6 is the liquid fractions of all of the PCMs with equal mass ratios as a function of dimensionless time for five different combinations of melting temperatures in reversed and forward arrangement respectively.One can see that it takes less time for bigger ∆T value to complete the solidification process in the reversed arrangement (Figure 6a), but more time for the forward arrangement (Figure 6b).The main reason for this is the arrangement order of the m-PCMs.Thus, the reversed arrangement of m-PCMs can shorten solidification time.
Figure 6 is the liquid fractions of all of the PCMs with equal mass ratios as a function of dimensionless time for five different combinations of melting temperatures in reversed and forward arrangement respectively.One can see that it takes less time for bigger ΔT value to complete the solidification process in the reversed arrangement (Figure 6a), but more time for the forward arrangement (Figure 6b).The main reason for this is the arrangement order of the m-PCMs.Thus, the reversed arrangement of m-PCMs can shorten solidification time.Figure 7 plots the dimensionless outlet temperature (θ(out)) of the heat exchanger with an equal mass ratio as a function of dimensionless time (τ) in the reversed and forward arrangements of m-PCMs, respectively.The dimensionless outlet temperature is defined as: (16) where Tin is the inlet temperature of the HTF. Figure 7 plots the dimensionless outlet temperature (θ (out) ) of the heat exchanger with an equal mass ratio as a function of dimensionless time (τ) in the reversed and forward arrangements of m-PCMs, respectively.The dimensionless outlet temperature is defined as: where T in is the inlet temperature of the HTF.As shown in Figure 7a, the dimensionless outlet temperatures are close to each other and decrease quickly in the initial period (τ < 50), then separate at a constant value in the middle period (50 <  < 200), and finally coincide with each other in the discharging process.The constant value of θ(out) increases with ΔT because the PCM-1 is in the downward arrangement and the TPCM-1 increases with ΔT value.As shown in Figure 7b, when the m-PCMs are in the reversed arrangement, the values of θ(out) are close to each other, except in the middle period (125 < τ < 225).However, this separation in the forward arrangement is much smaller than that in the reversed arrangement.In addition, the value of θ(out) in the middle period decreases with the ΔT because the PCM-3 is in the downward arrangement and the TPCM-3 increases with ΔT.From Figure 7, one can find that the value of θ(out) in the reversed arrangement is higher than that in the forward arrangement.Comparing the value of θ(out) between the two imaginary lines from Figure 7, the θ(out) value in Figure 7a is bigger than in Figure 7b.
Thus, as we can see, it is more advantageous for the zigzag plate energy storage unit to arrange the m-PCMs in a reversed fashion than in a forward fashion, so the reversed arrangement is chosen in the following studies.
Figure 8 shows the liquid fractions of PCM as a function of dimensionless time at different As shown in Figure 7a, the dimensionless outlet temperatures are close to each other and decrease quickly in the initial period (τ < 50), then separate at a constant value in the middle period (50 < τ < 200), and finally coincide with each other in the discharging process.The constant value of θ (out) increases with ∆T because the PCM-1 is in the downward arrangement and the T PCM-1 increases with ∆T value.As shown in Figure 7b, when the m-PCMs are in the reversed arrangement, the values of θ (out) are close to each other, except in the middle period (125 < τ < 225).However, this separation in the forward arrangement is much smaller than that in the reversed arrangement.In addition, the value of θ (out) in the middle period decreases with the ∆T because the PCM-3 is in the downward arrangement and the T PCM-3 increases with ∆T.From Figure 7, one can find that the value of θ (out) in the reversed arrangement is higher than that in the forward arrangement.Comparing the value of θ (out) between the two imaginary lines from Figure 7, the θ (out) value in Figure 7a is bigger than in Figure 7b.
Energies 2016, 9, 394 9 of 17 Thus, as we can see, it is more advantageous for the zigzag plate energy storage unit to arrange the m-PCMs in a reversed fashion than in a forward fashion, so the reversed arrangement is chosen in the following studies.
Figure 8 shows the liquid fractions of PCM as a function of dimensionless time at different value of ∆T under the condition of v in = 6 m/s and T in = 303.15K with m-PCMs at a mass ratio of 1:1:1 in the reversed arrangement.As shown in Figure 8, although T PCM-3 is the lowest, it takes the shortest time to complete the freezing process among the three PCMs and the solidification rate of PCM-3 increases while the value of ∆T decreases.This can be explained as following: PCM-3 is located in the inlet region, which results in the largest solidification rate of the three PCMs, and the phase change temperature of PCM-3 decreases as the value of ∆T increases.
values of θ(out) are close to each other, except in the middle period (125 < τ < 225).However, this separation in the forward arrangement is much smaller than that in the reversed arrangement.In addition, the value of θ(out) in the middle period decreases with the ΔT because the PCM-3 is in the downward arrangement and the TPCM-3 increases with ΔT.From Figure 7, one can find that the value of θ(out) in the reversed arrangement is higher than that in the forward arrangement.Comparing the value of θ(out) between the two imaginary lines from Figure 7, the θ(out) value in Figure 7a is bigger than in Figure 7b.
Thus, as we can see, it is more advantageous for the zigzag plate energy storage unit to arrange the m-PCMs in a reversed fashion than in a forward fashion, so the reversed arrangement is chosen in the following studies.
Figure 8 shows the liquid fractions of PCM as a function of dimensionless time at different value of ΔT under the condition of vin = 6 m/s and Tin = 303.15K with m-PCMs at a mass ratio of 1:1:1 in the reversed arrangement.As shown in Figure 8, although TPCM-3 is the lowest, it takes the shortest time to complete the freezing process among the three PCMs and the solidification rate of PCM-3 increases while the value of ΔT decreases.This can be explained as following: PCM-3 is located in the inlet region, which results in the largest solidification rate of the three PCMs, and the phase change temperature of PCM-3 decreases as the value of ΔT increases.As can be seen from Figure 8, there is little difference in the liquid fraction curves of PCM-2 for different values of ∆T.Explanations for this observation could be as follows: the first is that the heat transfer in m-PCMs is controlled by the thermal conductivity in the solidification process; the second one is that there is no relationship between the phase change temperature of PCM-2 and the value of ∆T.
One can see from Figure 8 that although the phase change temperature of PCM-1 is the highest, it takes the longest time to complete the solidification process in all cases, which forms the rate-limiting factor of the whole solidification process in the reversed arrangement.The main reason for this result is because of the arrangement of the three PCMs and the location of PCM-1.When the m-PCMs is in the reversed arrangement, the HTF carries "cool" energy in and absorbs heat from PCM-3 first, then PCM-2 and finally PCM-1 during discharging process.The temperature of the HTF increases during this process.As a result, before the PCM-2 and PCM-3 are fully discharged, the temperature of the HTF at the inlet is always higher than that of PCM-2 and PCM-3.As seen from the liquid fraction curve of PCM-1, the solidification rate increases as the value of ∆T increases.An explanation to this observation could be as follows.In the initial stage of discharge process, heat stored in the PCMs releases in the form of sensible heat to the HTF.Due to the buffering effect of PCM-2, the temperature of HTF entering PCM-1 region does not change much with time.As a result, a bigger value of ∆T (meaning a higher T PCM-1 ) yields a quicker phase change process for PCM-1.The temperature of HTF decreases as the solidification progressing, making PCM-1 start to freeze at high ∆T values.
To quantitatively compare charging performance enhancements through the m-PCMs in reversed arrangement, the following ratio is introduced: where Q refers to the total amount of thermal energy stored and η is the ratio of the total heat stored in the m-PCMs (∆T ‰ 0 K) to that in a single PCM (∆T = 0 K).η can therefore be regarded as the m-PCMs' effectiveness.Figure 9 shows the change of η as a function of dimensionless time for different ∆T.For a given ∆T, η decreases first with time, reaches a minimum τ = τ min , then increases with time, peaks at τ = τ max , and finally decreases to 1.Both the τ min and τ max are functions of ∆T.The τ min lies in the range of 25-35, whereas the τ max is between 160 and 200.Dimensionless time durations with ∆T larger than 1 are far longer than that with ∆T smaller than 1.The m-PCMs' effectiveness at τ min and τ max are also a function of ∆T; a large ∆T gives a lower η at the τ min but a higher η at the τ max .For a given ∆T, η is higher than 1 for most of the discharging process, leading to an overall process intensification through the use of m-PCMs.As shown in the inset figure of Figure 9, the average η (η (ave) ) over all τ increases with the value of ∆T, which means a higher ∆T gives a larger extent of overall process enhancement in the case of the reversed arrangement.
process for PCM-1.The temperature of HTF decreases as the solidification progressing, making PCM-1 start to freeze at high ΔT values.
To quantitatively compare charging performance enhancements through the m-PCMs in reversed arrangement, the following ratio is introduced: (17) where Q refers to the total amount of thermal energy stored and η is the ratio of the total heat stored in the m-PCMs (Δ 0 K) to that in a single PCM (Δ 0 K).η can therefore be regarded as the m-PCMs' effectiveness.Figure 9 shows the change of η as a function of dimensionless time for different ΔT.For a given ΔT, η decreases first with time, reaches a minimum τ = τmin, then increases with time, peaks at τ = τmax, and finally decreases to 1.Both the τmin and τmax are functions of ΔT.The τmin lies in the range of 25-35, whereas the τmax is between 160 and 200.Dimensionless time durations with ΔT larger than 1 are far longer than that with ΔT smaller than 1.The m-PCMs' effectiveness at τmin and τmax are also a function of ΔT; a large ΔT gives a lower η at the τmin but a higher η at the τmax.For a given ΔT, η is higher than 1 for most of the discharging process, leading to an overall process intensification through the use of m-PCMs.As shown in the inset figure of Figure 9, the average η (η(ave)) over all τ increases with the value of ΔT, which means a higher ΔT gives a larger extent of overall process enhancement in the case of the reversed arrangement.

Effect of Different Ste or Re Numbers
In this part, the set of ∆T = 80 K was chosen to study the effect of various Ste or Re numbers on the discharging process when one of them is fixed while the other varies and the initial temperature of PCMs (T o ) is 803.15K.As shown in Figure 10, when the Ste number varies between ´1.376 and ´1.679 on the condition of v in = 6 m/s, the value of the Ste number almost has no impact on the solidification rate (Figure 10a) or θ (out) (Figure 10b).This is because the heat transfer between HTF and m-PCMs is controlled by conductivity, not by convection, and the temperature difference between HTF and m-PCMs is large.However, when the Re number varies from 91 to 364.1 with T in = 303.15K, the value of the Re number has a significant impact on the solidification rate and θ (out) and the Re can be calculated as Equation ( 18): where D H is the hydraulic diameter of the zigzag channels [30].
As shown in Figure 11a, the solidification rate increases with the Re number; however, the effect of Re on solidification weakens with an increasing Re number.As shown in Figure 11b, θ (out) decreases with the Re number, and the effect of Re on θ (out) weakens with an increasing Re number.An explanation for this observation could be as follows.With increasing Re value, the more "cool" energy is that transferred to the heat exchanger, the more heat that could be extracted from the m-PCMs, resulting in a higher solidification rate.The larger value of Re, the more "cool" energy that is transferred downstream, resulting in θ (out) decreasing with Re.
As shown in Figures 10b and 11b, the dimensionless time dependence of θ (out) shows a three-stage process.With θ (out) decreasing with τ rapidly at first, then no effect of τ on θ (out) , and finally θ (out) decreasing with τ again.This dependence is because the temperature difference between m-PCMs and HTF is large and the m-PCMs releases heat energy in sensible and latent way at first, but mostly in a sensible form.As the discharging process continues, the m-PCMs starts to freeze and release heat energy mostly in the latent form, which leads to the second stage of θ (out) .Finally, as the m-PCMs release most of their latent heat, the heat energy released to the HTF is controlled by a sensible form; thus, the θ (out) decreases.However, the temperature difference between the HTF and m-PCMs is smaller than the first two stages, which results in a much slower decrease of θ (out) than during the first stage.

Effect of Different Ste or Re Numbers
In this part, the set of ΔT = 80 K was chosen to study the effect of various Ste or Re numbers on the discharging process when one of them is fixed while the other varies and the initial temperature of PCMs (To) is 803.15K.As shown in Figure 10, when the Ste number varies between −1.376 and −1.679 on the condition of vin = 6 m/s, the value of the Ste number almost has no impact on the solidification rate (Figure 10a) or θ(out) (Figure 10b).This is because the heat transfer between HTF and m-PCMs is controlled by conductivity, not by convection, and the temperature difference between HTF and m-PCMs is large.However, when the Re number varies from 91 to 364.1 with Tin = 303.15K, the value of the Re number has a significant impact on the solidification rate and θ(out) and the Re can be calculated as Equation ( 18): where DH is the hydraulic diameter of the zigzag channels [30].As shown in Figure 11a, the solidification rate increases with the Re number; however, the effect of Re on solidification weakens with an increasing Re number.As shown in Figure 11b, θ(out) decreases with the Re number, and the effect of Re on θ(out) weakens with an increasing Re number.An explanation for this observation could be as follows.With increasing Re value, the more "cool" energy is that transferred to the heat exchanger, the more heat that could be extracted from the m-PCMs, resulting in a higher solidification rate.The larger value of Re, the more "cool" energy that is transferred downstream, resulting in θ(out) decreasing with Re. m-PCMs is controlled by conductivity, not by convection, and the temperature difference between HTF and m-PCMs is large.However, when the Re number varies from 91 to 364.1 with Tin = 303.15K, the value of the Re number has a significant impact on the solidification rate and θ(out) and the Re can be calculated as Equation ( 18): where DH is the hydraulic diameter of the zigzag channels [30].As shown in Figure 11a, the solidification rate increases with the Re number; however, the effect of Re on solidification weakens with an increasing Re number.As shown in Figure 11b, θ(out) decreases with the Re number, and the effect of Re on θ(out) weakens with an increasing Re number.An explanation for this observation could be as follows.With increasing Re value, the more "cool" energy is that transferred to the heat exchanger, the more heat that could be extracted from the m-PCMs, resulting in a higher solidification rate.The larger value of Re, the more "cool" energy that is transferred downstream, resulting in θ(out) decreasing with Re.

Effect of Different Ste/Re under a Given Inlet Power
The effects of inlet temperature and velocity are considered for a given input power and expressed as: where A a is cross-sectional area, v in is the inlet velocity of the HTF and T b is the base temperature defined as T b " 1{2 ´r T in `To ¯.The base case of this set of modelling work is under the conditions of r v in " 6 m s , r T in " 303.15 K, T o = 803.15K, and the set of ∆T = 80 K is chosen.Different conditions are realized by changing the inlet temperature and the inlet velocity of the HTF under the given inlet power.
Figure 12 shows the time evolution of the liquid fractions at different Re/Ste numbers under the given inlet power condition.As shown in Figure 12, the solidification rate of all of the m-PCMs increases as the Re number increases under the given inlet power condition.Such an increase in the solidification rate does not seem to be associated with the overall heat transfer driving force (the temperature difference).It is plausible that an increase in the Re number results in an enhancement of the overall heat transfer between the HTF and the m-PCMs and hence an overall increase in the solidification rate.This increase in the solidification rate can be explained as follows.For a given input power, a decrease in the inlet velocity of the HTF (Re number) indicates a decrease in the inlet temperature of the HTF.This implies that at the inlet region the temperature difference between the HTF and m-PCMs remains high even after the complete freezing of the m-PCMs, leading to more "cool" energy stored in the inlet region as sensible heat and less "cool" energy flow to the downstream region to freeze the m-PCMs there, resulting in the observed low solidification rate.
expressed as: where Aa is cross-sectional area, vin is the inlet velocity of the HTF and Tb is the base temperature defined as 1/2 .The base case of this set of modelling work is under the conditions of 6 , 303.15 K, To = 803.15K, and the set of ΔT = 80 K is chosen.Different conditions are realized by changing the inlet temperature and the inlet velocity of the HTF under the given inlet power.
Figure 12 shows the time evolution of the liquid fractions at different Re/Ste numbers under the given inlet power condition.As shown in Figure 12, the solidification rate of all of the m-PCMs increases as the Re number increases under the given inlet power condition.Such an increase in the solidification rate does not seem to be associated with the overall heat transfer driving force (the temperature difference).It is plausible that an increase in the Re number results in an enhancement of the overall heat transfer between the HTF and the m-PCMs and hence an overall increase in the solidification rate.This increase in the solidification rate can be explained as follows.For a given input power, a decrease in the inlet velocity of the HTF (Re number) indicates a decrease in the inlet temperature of the HTF.This implies that at the inlet region the temperature difference between the HTF and m-PCMs remains high even after the complete freezing of the m-PCMs, leading to more "cool" energy stored in the inlet region as sensible heat and less "cool" energy flow to the downstream region to freeze the m-PCMs there, resulting in the observed low solidification rate.Figure 13 shows the time evolution of the outlet temperature at different Re/Ste numbers under the given inlet power condition.One can see that the outlet temperature decreases quickly first, then reaches a certain value, and finally decreases relatively fast again.This is because the temperature of the m-PCMs is higher than the phase change temperature at the initial period and the heat released to the HTF is in a sensible form.As the discharge progresses, the m-PCMs begin to freeze and release Figure 13 shows the time evolution of the outlet temperature at different Re/Ste numbers under the given inlet power condition.One can see that the outlet temperature decreases quickly first, then reaches a certain value, and finally decreases relatively fast again.This is because the temperature of the m-PCMs is higher than the phase change temperature at the initial period and the heat released to the HTF is in a sensible form.As the discharge progresses, the m-PCMs begin to freeze and release heat mainly in latent form; over time, the heat released to the HTF is mainly in a sensible form in the final stage because most of the latent heat has been released to the HTF already.
Energies 2016, 9, 394 12 of 16 heat mainly in latent form; over time, the heat released to the HTF is mainly in a sensible form in the final stage because most of the latent heat has been released to the HTF already.It is interesting to examine the time evolution of the ratio of the latent heat (QL) to the total stored heat (Qall) under a given input power.The ratio of the latent heat to the total stored heat is calculated as in Equation ( 20): It is interesting to examine the time evolution of the ratio of the latent heat (Q L ) to the total stored heat (Q all ) under a given input power.The ratio of the latent heat to the total stored heat is calculated as in Equation ( 20): Figure 14 shows that Q L Q all rises abruptly after an initial stage and reaches a peak which increases with the Re number; then they decrease slowly and intersect with each other in the latter half of the discharging process.The difference among the values of Q L Q all is small (0.54-0.57) at the moment that the m-PCMs are fully frozen.As expressed in Equation ( 20), the value of lim It is interesting to examine the time evolution of the ratio of the latent heat (QL) to the total stored heat (Qall) under a given input power.The ratio of the latent heat to the total stored heat is calculated as in Equation ( 20): Figure 14 shows that rises abruptly after an initial stage and reaches a peak which increases with the Re number; then they decrease slowly and intersect with each other in the latter half of the discharging process.The difference among the values of is small (0.54-0.57) at the moment that the m-PCMs are fully frozen.As expressed in Equation ( 20), the value of lim → equals 0.271, 0.391, 0.405 and 0.415 when the Re number equals 273, 319, 364 and 410, respectively, which leads to an intersection in the latter half of the discharging process.Figure 15 shows the liquid fraction of m-PCMs as a function of time at different periods of inlet velocity under the given inlet power condition.For a given inlet power, the inlet velocity (vin) and inlet temperature (Tin) are calculated as in Equations ( 21) and (22), respectively.Figure 15 shows the liquid fraction of m-PCMs as a function of time at different periods of inlet velocity under the given inlet power condition.For a given inlet power, the inlet velocity (v in ) and inlet temperature (T in ) are calculated as in Equations ( 21) and ( 22), respectively.
T in " T b ´Pin c p ¨ρ¨A a ¨vin where T is the period of the inlet temperature.B and C are the amplitude and a constant value, which equal 1.5 and 7.5, respectively.A a is the cross-sectional area of inlet.To find the impact of a period on the solidification of m-PCMs, two small periods (T = 30 s/60 s) and two large periods (T = 5 min/10 min) are chosen in this work.As one can see from Figure 15, the liquid fraction curves of PCM-1, PCM-2 and PCM-3 almost coincide with each other, which means that the value of the period (30 s < T < 10 min) has little or no effect on the solidification of m-PCMs in the zigzag plate unit.
Figure 16 shows the time evolution of the outlet temperatures at different T under the given inlet power condition.The amplitude of the inlet temperature is 80 K, but the amplitude of the outlet temperature is much smaller (<10 K) at different values of the period, which indicates that the zigzag plate with m-PCMs has a good impact on the temperature shock as "filter action" in the communication field.which equal 1.5 and 7.5, respectively.Aa is the cross-sectional area of inlet.To find the impact of a period on the solidification of m-PCMs, two small periods (T = 30 s/60 s) and two large periods (T = 5 min/10 min) are chosen in this work.As one can see from Figure 15, the liquid fraction curves of PCM-1, PCM-2 and PCM-3 almost coincide with each other, which means that the value of the period (30 s < T < 10 min) has little or no effect on the solidification of m-PCMs in the zigzag plate unit.Figure 16 shows the time evolution of the outlet temperatures at different T under the given inlet power condition.The amplitude of the inlet temperature is 80 K, but the amplitude of the outlet temperature is much smaller (<10 K) at different values of the period, which indicates that the zigzag plate with m-PCMs has a good impact on the temperature shock as "filter action" in the communication field.

Conclusions
An experimental system was built to validate the two-dimensional mathematical model that was established to describe the discharging behaviour of the device.A heat exchange device with a zigzag configuration containing m-PCMs was investigated in this work.A new inlet power boundary condition was defined.The following conclusions were obtained in this work: (1) It is more advantageous for the plate energy storage unit to arrange the m-PCMs in a reversed direction than in a forward direction.(2) If only the Ste (inlet temperature) or Re (inlet velocity) numbers were changed, there was almost no impact of the Ste number on the discharge process of the m-PCMs unit, but there Figure 16 shows the time evolution of the outlet temperatures at different T under the given inlet power condition.The amplitude of the inlet temperature is 80 K, but the amplitude of the outlet temperature is much smaller (<10 K) at different values of the period, which indicates that the zigzag plate with m-PCMs has a good impact on the temperature shock as "filter action" in the communication field.

Conclusions
An experimental system was built to validate the two-dimensional mathematical model that was established to describe the discharging behaviour of the device.A heat exchange device with a zigzag configuration containing m-PCMs was investigated in this work.A new inlet power boundary condition was defined.The following conclusions were obtained in this work: (1) It is more advantageous for the plate energy storage unit to arrange the m-PCMs in a reversed direction than in a forward direction.(2) If only the Ste (inlet temperature) or Re (inlet velocity) numbers were changed, there was almost no impact of the Ste number on the discharge process of the m-PCMs unit, but there

Conclusions
An experimental system was built to validate the two-dimensional mathematical model that was established to describe the discharging behaviour of the device.A heat exchange device with a zigzag configuration containing m-PCMs was investigated in this work.A new inlet power boundary condition was defined.The following conclusions were obtained in this work: (1) It is more advantageous for the plate energy storage unit to arrange the m-PCMs in a reversed direction than in a forward direction.(2) If only the Ste (inlet temperature) or Re (inlet velocity) numbers were changed, there was almost no impact of the Ste number on the discharge process of the m-PCMs unit, but there was a significant impact of the Re number; for a given inlet power inlet condition, it was more advantageous to increase the Re number than the Ste number to shorten the phase change time in the discharge process.(3) For a given inlet power condition with a periodically varying velocity, the value of the period (30 s-10 min) had almost no impact on the discharge process, which indicates that the zigzag plate with m-PCMs had a positive impact on the temperature shock.
Author Contributions: All the authors contributed to the manuscript.Peilun Wang performed the main research and wrote the manuscript.Yun Huang provide guidance and revised the paper, Yulong Ding, Zhijiang Peng, and Yi Wang provided supervision.Dacheng Li and Xingang Zheng contributed in the experimental and numerical parts of the paper, respectively.

Conflicts of Interest:
The authors declare no conflict of interest.

Figure 1 .
Figure 1.Physical configuration of the zigzag configuration.H: height; M: tooth length; PCM: phase-change-material; and HTF: heat transfer fluid.

Figure 1 .
Figure 1.Physical configuration of the zigzag configuration.H: height; M: tooth length; PCM: phase-change-material; and HTF: heat transfer fluid.

Figure 4 .
Figure 4. Zigzag plate heat exchanger unit.(a) Arrangement of the thermocouples; (b) surface of the zigzag plates; and (c) internal structure.

Figure 4 .
Figure 4. Zigzag plate heat exchanger unit.(a) Arrangement of the thermocouples; (b) surface of the zigzag plates; and (c) internal structure.

Figure 5
Figure 5 compares the measured and calculated time evolutions of the temperature at the outlet and the centre of PCM (as indicated in Figure4a).It can be seen that reasonably good agreement has been achieved, thus establishing confidence in the model.The deviation of thermocouple 1 is likely associated with the existence of pores in the PCM because the material was loaded into the plate heat exchangers in compacted powder form.Other reasons may include the temperature dependence of the PCM properties, an uneven distribution of the HTF flow across all of the plates, which were not considered in the modelling, or the location of the thermocouple.

Figure 5 .
Figure 5.Comparison between the experimental and modelling results for the center and the outlet.

Figure 5 .
Figure 5.Comparison between the experimental and modelling results for the center and the outlet.

Figure 6
Figure 6 plots the liquid fraction of total PCMs as a function of time and the three PCMs are of equal mass ratio.The temperature and time expressed in dimensionless forms defined as follows:Dimensionless time:

Figure 6 .
Figure 6.Dimensionless time evolution of the liquid fraction of all of the m-PCMs at different ΔT.(a) Reversed arrangement; and (b) forward arrangement.

Figure 6 .
Figure 6.Dimensionless time evolution of the liquid fraction of all of the m-PCMs at different ∆T.(a) Reversed arrangement; and (b) forward arrangement.

Figure 7 .
Figure 7. Dimensionless time evolution of the dimensionless outlet temperature.(a) Reversed arrangement; and (b) forward arrangement.

Figure 7 .
Figure 7. Dimensionless time evolution of the dimensionless outlet temperature.(a) Reversed arrangement; and (b) forward arrangement.

Figure 8 .
Figure 8. Liquid fractions of the three m-PCMs as a function of dimensionless time at different value of ΔT (empty, half-filled and filled symbols stand for liquid fractions of PCM-1, PCM-2 and PCM-3, respectively).

Figure 8 .
Figure 8. Liquid fractions of the three m-PCMs as a function of dimensionless time at different value of ∆T (empty, half-filled and filled symbols stand for liquid fractions of PCM-1, PCM-2 and PCM-3, respectively).

Figure 9 .Figure 9 .
Figure 9. Dimensionless time dependence of m-PCMs effectiveness η (inset is the average η(ave) overall η by integrating over the distribution).

Figure 10 .
Figure 10.(a) Liquid fractions of the m-PCMs and (b) dimensionless outlet temperature as a function of dimensionless time at different Ste numbers.

Figure 11 .Figure 10 .
Figure 11.(a) Liquid fractions of the m-PCMs and (b) dimensionless outlet temperature as a function of dimensionless time at different Re numbers.

Figure 10 .
Figure 10.(a) Liquid fractions of the m-PCMs and (b) dimensionless outlet temperature as a function of dimensionless time at different Ste numbers.

Figure 11 .Figure 11 .
Figure 11.(a) Liquid fractions of the m-PCMs and (b) dimensionless outlet temperature as a function of dimensionless time at different Re numbers.

Figure 12 .
Figure 12.Time evolution of liquid fractions at different inlet velocity (Re)/inlet temperature (Ste) (Re/Ste) numbers under the given inlet power condition.

Figure 12 .
Figure 12.Time evolution of liquid fractions at different inlet velocity (Re)/inlet temperature (Ste) (Re/Ste) numbers under the given inlet power condition.

Figure 13 .
Figure 13.Time evolution of the outlet temperature at different inlet velocity (Re)/inlet temperature (Ste) (Re/Ste) numbers under the given inlet power.

Figure 13 .
Figure 13.Time evolution of the outlet temperature at different inlet velocity (Re)/inlet temperature (Ste) (Re/Ste) numbers under the given inlet power.

Figure 13 .
Figure 13.Time evolution of the outlet temperature at different inlet velocity (Re)/inlet temperature (Ste) (Re/Ste) numbers under the given inlet power.

Figure 14 .
Figure 14.Time evolution of the ratio of the latent heat to the total stored heat at different inlet velocity (Re)/inlet temperature (Ste) (Re/Ste) numbers.

Figure 14 .
Figure 14.Time evolution of the ratio of the latent heat to the total stored heat at different inlet velocity (Re)/inlet temperature (Ste) (Re/Ste) numbers.

Figure 15 .
Figure 15.Liquid fractions of the m-PCMs as a function of time at different periods (T) under the given inlet power condition (empty, half-filled and filled symbols stand for liquid fractions of PCM-1, PCM-2 and PCM-3, respectively).

Figure 16 .
Figure 16.Time evolution of the outlet temperatures at different periods (T) under the given inlet power condition.

Figure 15 .Figure 15 .
Figure 15.Liquid fractions of the m-PCMs as a function of time at different periods (T) under the given inlet power condition (empty, half-filled and filled symbols stand for liquid fractions of PCM-1, PCM-2 and PCM-3, respectively).

Figure 16 .
Figure 16.Time evolution of the outlet temperatures at different periods (T) under the given inlet power condition.

Figure 16 .
Figure 16.Time evolution of the outlet temperatures at different periods (T) under the given inlet power condition.

Table 1 .
Thermophysical properties of the NaCl-MgCl2 salt for the validation experiments.

Table 1 .
Thermophysical properties of the NaCl-MgCl 2 salt for the validation experiments.

Table 1 .
Thermophysical properties of the NaCl-MgCl2 salt for the validation experiments.