Transient Numerical Simulation of the Melting and Solidification Behavior of NaNO 3 Using a Wire Matrix for Enhancing the Heat Transfer

The paper presents the results of a transient numerical investigation of the melting and solidification process of sodium nitrate (NaNO3), which is used as phase change material. For enhancing the heat transfer to the sodium nitrate an aluminum wire matrix is used. The numerical simulation of the melting and solidification process was done with the enthalpy-porosity approach. The numerical analysis of the melting process has shown that apart from the first period of the charging process, where heat conduction is the main heat transfer mechanism, natural convection is the dominant heat transfer mechanism. The numerical investigation of the solidification process has shown that the dominant heat transfer mechanism is heat conduction. Based on the numerical results, the discharging process has been slower than the charging process. The performance of the charged and discharged power has shown that the wire matrix is an alternative method to enhance the heat transfer into the phase change material.


Introduction
The aim of the European Union is a reduction of the greenhouse gas emissions (referred to 1990) as well as an increase of the energy efficiency of 20% till 2020.To achieve this ambitious goal, it is necessary to increase the target share of renewable energies in the energy consumption as well as the energy efficiency.This relates to the situation that renewable energies, energy storage and energy distribution will be the most important and strongly growing markets for the next decades and hence also key topics for research and development.
The storage of thermal energy is a special form of energy storage.Beside the simple storage and release of heat, thermal energy storage systems can play an important role in electricity production in case of an integration into a thermodynamic cycle (Rankine, Brayton, etc.).To store thermal energy from, e.g., superheated steam, a storage system based on latent heat is necessary to store the heat of condensation in an exergetically optimal way.In such a system, the thermal energy transfer occurs when a so called phase change material (PCM) changes the aggregate state (phase) from, e.g., solid to solid, solid to liquid, liquid to gas, or solid to gas and vice versa.Initially, the temperature of PCMs increases like a conventional sensible storage material as they absorb heat until it reaches the melting point.Compared to sensible storage materials, a PCM absorbs and releases heat at a nearly constant temperature during the phase change and they can store 5-14 times more heat per unit volume than sensible storage materials such as water, sand, or rock [1].This results in a significant reduction of the storage volume using PCMs as storage material compared to sensible heat storage.That means Energies 2016, 9, 205 2 of 18 that the latent heat thermal energy storage system (LHTES) has a higher energy density compared to sensible storage systems.A further advantage of the latent heat thermal energy storage system is that the charging and discharging process takes place at small temperature differences.
The selection of the PCM material used in an application depends on the process parameters as well as the thermodynamic, physical, chemical, and economic aspects.By selecting a PCM material for a specific application the melting temperature must be selected close to the operation temperature.The volume based latent heat should also be as high as possible to reduce the size of the storage device.A further important specification on phase change materials is given to the thermal conductivity in both states, solid and liquid, which is desirable to improve the heat transfer into the phase change material.To maintain convenient container design for the phase change device the density variation during phase change from solid to liquid and vice versa as well as the vapor pressure at operation temperature should be low.
For the use in power cycles as well as industrial applications phase change materials in the medium and high temperature range are required.In these temperature ranges, a high number of materials are available like metals and metal alloys or inorganic salts and saline compounds.However, the most discussed storage materials for LHTES are salts and their eutectics.These materials are characterized by the low heat conductivity [2][3][4].A widely used method to overcome this heat transport limitation in LHTES is to increase the specific surface of the heat exchanger tubes.Various techniques, like fins [4][5][6][7][8][9][10][11], dispersed high-conductivity particles inside the PCM [12,13], metal and graphite-compound matrices [14][15][16], micro-encapsulation [17][18][19] as well as a combination of metal foam made of copper and sandwich structure fins [20] have been suggested and were analyzed.
Based on the low thermal conductivity of the PCMs the details of heat transfer and phase change in the presence of heat transfer enhancing techniques has become of interest.To analyze the melting and freezing process of a material under the influence of gravity, a higher number of analytical [21][22][23][24][25] as well as numerical methods [26][27][28] were developed in the last decades.The change with time of the position of the phase interface during the melting and/or solidification process is the main problem concerning the numerical simulation.As reported in [29], more than half of all investigations done on latent heat thermal storage systems were numerically and most of the used models are 2D (see e.g., [28,30]).The very slow laminar flow (natural convection) and the variation of the physical properties of the PCM during the phase change result in difficulties for the numerical simulation.But the main challenge of such problems is the presence of a moving liquid-solid interface involving a strong coupling of mass and heat transfer.To get accurate results, very small time steps as well as a fine structured computational grid is necessary.Therefore long computation time can be expected even for small 3D models [10].
In the present work the results of a transient numerical analysis on the melting and solidification behavior of sodium nitrate (NaNO 3 ), which is used as phase change material in a LHTES, will be presented.The heat transfer fluid flows inside a steel tube of the heat exchanger.The steel tube is surrounded by an aluminum wire matrix for enhancing the heat transfer into the PCM.The material aluminum was chosen based on the high thermal conductivity, which is approximately 10 times higher compared to steel, the low material density of approximately 2700 kg/m 3 and the applicable temperature range up to 330 ˝C.A further advantage of the selected material combination is that galvanic corrosion between aluminum, carbon steel and sodium nitrate is not critical, as mentioned in [31].Parts of the present study are presented in [32].

Mathematical Model of the Heat Exchanger Tube
In the next paragraphs, the physical and mathematical model as well as the thermo-physical properties for the materials will be presented.

Physical Model of the Wired Matrix
Figure 1 shows a sketch of the analyzed heat exchanger tube used in the present study for the heat transfer from the heat transfer fluid, which flows inside the tube, to the surrounding sodium nitrate Energies 2016, 9, 205 3 of 18 (NaNO 3 ), which is used as phase change material.The heat exchanger tube consists of a plain steel tube surrounded by an aluminum wire matrix.The analysis was done for a vertical arrangement of the heat exchanger tube.The main advantage of using a wire matrix for enhancing the heat transfer to the phase change material is the simple geometry as well as the simple assembling of the tube by inserting and pre-stressing of the wires, which are connected to the steel tube.Under operation temperature (up to 330 ˝C), the wire matrix is highly flexible and therefore the different linear expansion between steel tube and wire matrix can be easily compensated.
Energies 2016, 9, 205 3 of 18 arrangement of the heat exchanger tube.The main advantage of using a wire matrix for enhancing the heat transfer to the phase change material is the simple geometry as well as the simple assembling of the tube by inserting and pre-stressing of the wires, which are connected to the steel tube.Under operation temperature (up to 330 °C), the wire matrix is highly flexible and therefore the different linear expansion between steel tube and wire matrix can be easily compensated.

Mathematical Model
To obtain the numerical solution for the phase change problem the commercial computer code ANSYS FLUENT (ANSYS Inc., Canonsburg, PA, USA) was used.In ANSYS FLUENT the enthalpy-porosity approach [33,34] is used to calculate the melting and/or solidification process.
The enthalpy-porosity method depends on the liquid volume fraction γ which denotes the ratio of volume of liquid to the total volume in each cell.Therefore, the liquid fraction γ describes the fraction of a cell volume which is filled with liquid.For the calculation the computational domain in the PCM will be divided into three regions: solid, liquid, and mushy zone.With this technique the melting front is not calculated explicitly.Instead, the quantity liquid fraction γ is computed at each time step for every cell based on the enthalpy balance [33].Computational cells, where the phase of the material changes, are modeled as a pseudo porous medium with liquid fraction γ ranging between 0 and 1, for solid and liquid respectively, while 0 < γ < 1 denotes the mushy zone.Thus, the energy equation used here for the PCM system yields to: with the total enthalpy and the sensible enthalpy The enthalpy change during the phase change (latent heat) can be written as a function of the liquid fraction and the latent heat of fusion Hlat = γL.Based on a method proposed by [34], the liquid fraction γ can be determined by temperature.
Figure 1.Heat exchanger tube with wire.

Mathematical Model
To obtain the numerical solution for the phase change problem the commercial computer code ANSYS FLUENT (ANSYS Inc., Canonsburg, PA, USA) was used.In ANSYS FLUENT the enthalpy-porosity approach [33,34] is used to calculate the melting and/or solidification process.
The enthalpy-porosity method depends on the liquid volume fraction γ which denotes the ratio of volume of liquid to the total volume in each cell.Therefore, the liquid fraction γ describes the fraction of a cell volume which is filled with liquid.For the calculation the computational domain in the PCM will be divided into three regions: solid, liquid, and mushy zone.With this technique the melting front is not calculated explicitly.Instead, the quantity liquid fraction γ is computed at each time step for every cell based on the enthalpy balance [33].Computational cells, where the phase of the material changes, are modeled as a pseudo porous medium with liquid fraction γ ranging between 0 and 1, for solid and liquid respectively, while 0 < γ < 1 denotes the mushy zone.Thus, the energy equation used here for the PCM system yields to: with the total enthalpy and the sensible enthalpy The enthalpy change during the phase change (latent heat) can be written as a function of the liquid fraction and the latent heat of fusion H lat = γL.Based on a method proposed by [34], the liquid fraction γ can be determined by temperature.
Energies 2016, 9, 205 4 of 18 with T sol = 578.65K and T liq = 579.65K.The temperatures T sol and T liq determine the size of the mushy zone.As described in [35], an iteration between the energy balance Equation (1) and the liquid fraction Equation ( 4) is done for the calculation of the temperature.This method was suggested by [34] because a direct calculation of Equation ( 4) to update the liquid fractions results in poor convergence of the energy balance Equation (1).The continuity equation used can be written as: and the momentum equation for the enthalpy-porosity equation is given in Equation (6). with the porosity function A C (γ) to reduce gradually the velocities from a finite value in the liquid to zero in the solid.The value b in Equation ( 7) is a small computational constant (b = 0.001) used to avoid division by zero while the mushy zone constant C is supposed to reflect the structure of the melting front.
For the mushy zone constant the standard value of ANSYS FLUENT with C = 10 5 was used.The higher the value of the mushy zone constant, the steeper the damping term becomes, and the faster the velocity drops to zero as the material solidifies, see also [36].
The basic conservation equations of continuity, momentum, and energy were solved numerically, using the ANSYS FLUENT 14.5 software.For the pressure-velocity coupling, the COUPLED scheme, which is implemented in ANSYS FLUENT, was used for the present study.
The following assumptions are made for the numerical investigation: 1.Both the solid as well as the liquid phase is homogeneous and isotropic, and the melting process is symmetric within a segment.2. NaNO 3 in the liquid phase is considered to be an incompressible and Newtonian fluid.3. The volume change upon phase change is ignored.4. For the molten PCM laminar flow as Newtonian fluid is assumed.5.The solid is homogeneously distributed in the mushy region.6.It is assumed that the PCM has an ideal solidification behavior.Therefore, the subcooling effects are neglected and the solidification temperature is constant.7. A perfect contact between the surfaces of the single wires and also to the steel tube is assumed.
A normal downward-directed gravity field with the corresponding gravitational acceleration of 9.81 m/s 2 is considered.

Thermo-Physical Properties
The thermo-physical properties for the tube and wire matrix are taken from the ANSYS FLUENT properties data base.For the properties of the phase change material NaNO 3 a user defined function (UDF) was programmed.
Important values for the numerical simulation of phase change problems are the melting temperature, heat of fusion, specific heat capacity, heat conduction, dynamic viscosity and density.In the literature, different data for the melting temperature of the sodium nitrate are available.In the present work a constant melting temperature of T = 579.15K was used for NaNO 3 , which is based on the references [37,38].The heat of fusion with L = 176.256kJ/kg is taken from [38].It can be seen in the references [37,38], which have made a comparison with a higher number of published data, that the gradient of the specific heat capacity in the liquid region of NaNO 3 is approximately constant.In the solid region of the sodium nitrate, the specific heat capacity shows a slight upward trend starting from the environment temperature up to the melting temperature.At temperature levels of 523.15K < T < 573.15 K, which is close to the phase change temperature, the specific heat shows a local maximum (at approximately 543.15 K) as a result of a change of the crystal structure.In the region between approximately 573.15 K < T ď 579.15K the specific heat capacity is approximately constant and has a value close to the specific heat capacity of the liquid sodium nitrate.However, the local maximum of the specific heat capacity is not within the temperature region of the numerical simulation.Therefore, the specific heat capacity for the solid as well as liquid NaNO 3 was modeled within the present study with the constant value of c p = 1635 J/kgK (this is the value of the liquid NaNO 3 ) based on the references [37,38].
For the density of the sodium nitrate in the solid phase, the value of 2119.6 kg/m 3 was used according to [3,37,39], while for the liquid region the equation which is based on the data of [40], was developed.A comparison with available data in the literature shows that the polynomial used for the liquid density of NaNO 3 in this work is located between the data of [41,42].Following [43], the dynamic viscosity of the liquid PCM has been expressed as For the thermal conductivity in the liquid phase the correlation of [44] λ liq " 0.62 ´1.787 ¨10 ´4T (10) was used in the UDF.Based on the measurement data of [37], an equation for the thermal conductivity in the solid phase was derived.
If different values for a physical property at the solidus and the liquidus line are given, then a linear interpolation between these values is done to get the thermo-physical properties in the mushy zone.

Numerical Model of the Heat Exchanger Tube
The numerical simulation of the melting and solidification behavior of a material must be done in such a way that the grid size used in the computational domain as well as the used time step is very small.Otherwise, it is not possible to get realistic solutions for the problem under investigation.Therefore it is not possible to simulate the total heat exchanger tube of a LHTES with a three dimensional domain, because for large storage devices the single tube length in the heat exchanger can be up to several meters which will consume enormous computing resources.However, for a basic understanding of the physical mechanism small sized computational domains can be used.Therefore, for the three dimensional simulations made for the present study a height of 43.8 mm was used for the computational domain of the vertical arranged heat exchanger tube with wire matrix, which is presented in Figure 2. The computational domain consists of the steel tube (outer diameter is 33.7 mm) with an aluminum wire matrix and the PCM sodium nitrate.The cross section of the wire is a square with the dimension 3.65 mm ˆ3.65 mm.
In Figure 3, the boundary conditions and the direction of the gravity used for the numerical simulations are depicted.At the top and the bottom surfaces as well as at the outer shell of the NaNO 3 an adiabatic wall is defined.The adiabatic boundary conditions avoid a heat exchange with the surrounding through these surfaces.The consequence of the adiabatic boundary condition is that during the charging process no dissipation of the heat coming from the heat transfer fluid through the steel tube is given.In real thermal heat storage devices the container is completely insulated against the environment.The thickness of this insulation is chosen in such a way that the thermal loss is as low as economically feasible.Thus, nearly adiabatic conditions are given for such a storage unit.Therefore the assumption of an adiabatic wall is justified.The computational domain consists of the steel tube (outer diameter is 33.7 mm) with an aluminum wire matrix and the PCM sodium nitrate.The cross section of the wire is a square with the dimension 3.65 mm × 3.65 mm.
In Figure 3, the boundary conditions and the direction of the gravity used for the numerical simulations are depicted.At the top and the bottom surfaces as well as at the outer shell of the NaNO3 an adiabatic wall is defined.The adiabatic boundary conditions avoid a heat exchange with the surrounding through these surfaces.The consequence of the adiabatic boundary condition is that during the charging process no dissipation of the heat coming from the heat transfer fluid through the steel tube is given.In real thermal heat storage devices the container is completely insulated against the environment.The thickness of this insulation is chosen in such a way that the thermal loss is as low as economically feasible.Thus, nearly adiabatic conditions are given for such a storage unit.Therefore the assumption of an adiabatic wall is justified.At the side walls of the computational domain the symmetry boundary condition was used.For the numerical simulations, the initial temperature of the whole system (tube, wire and NaNO3) was 569.15K for the melting and 589.15K for the solidification process.The tube wall surface temperature Tw, which represents the boundary condition at the inner tube diameter surface, was 599.15K for the melting and 559.15K for the solidification process.The tube wall surface temperature was kept constant during the whole simulation.
The grid size used for the discretization of the computational grid was determined after a careful examination of the results of a grid refinement process.As a time step for the numerical At the side walls of the computational domain the symmetry boundary condition was used.For the numerical simulations, the initial temperature of the whole system (tube, wire and NaNO 3 ) was 569.15K for the melting and 589.15K for the solidification process.The tube wall surface temperature T w , which represents the boundary condition at the inner tube diameter surface, was 599.15K for the melting and 559.15K for the solidification process.The tube wall surface temperature was kept constant during the whole simulation.
The grid size used for the discretization of the computational grid was determined after a careful examination of the results of a grid refinement process.As a time step for the numerical simulation, a value of 0.05 s was used, which is also a result of preliminary tests [45].A further decrease of the time step size as well as of the size of control volumes of the computational grid did not show noticeable changes in the results.For the computational domain shown in Figures 2 and 3 a total number of 65,624 control volumes is used.The value for the volume ratio, which is defined as: Volume of the wire matrix Volume of the PCM (12) is 0.149 and the surface to volume ratio ε " Surface area of the wire matrix Volume of the wire matrix (13) is 0.965 m ´1.
The numerical simulations presented in the present article are done on an 8 core computer and the computation time for the melting process was approximately two weeks while for the solidification process the computation time was about three weeks.

Numerical Results and Discussion
In the following, the results of the numerical simulation of the melting process will be presented.

Results for the Charging Process
In Figure 4 the chronological sequence of the volume averaged temperature (red curve) and volume averaged liquid fraction (blue curve) for the PCM during the charging process using a wire matrix for improving the heat transfer is presented.The black dashed line in Figure 4 shows the melting temperature of the PCM.For comparison reasons a further time sequence of a volume averaged NaNO 3 temperature of a heat exchanger tube design with longitudinal fins (green colored) is included in Figure 4.
Energies 2016, 9, 205 7 of 18 simulation, a value of 0.05 s was used, which is also a result of preliminary tests [45].A further decrease of the time step size as well as of the size of control volumes of the computational grid did not show noticeable changes in the results.For the computational domain shown in Figures 2 and 3 a total number of 65,624 control volumes is used.The value for the volume ratio, which is defined as: Volume of the wire matrix Volume of the PCM  (12) is 0.149 and the surface to volume ratio Surface area of the wire matrix Volume of the wire matrix  (13) is 0.965 m −1 .
The numerical simulations presented in the present article are done on an 8 core computer and the computation time for the melting process was approximately two weeks while for the solidification process the computation time was about three weeks.

Numerical Results and Discussion
In the following, the results of the numerical simulation of the melting process will be presented.

Results for the Charging Process
In Figure 4 the chronological sequence of the volume averaged temperature (red curve) and volume averaged liquid fraction (blue curve) for the PCM during the charging process using a wire matrix for improving the heat transfer is presented.The black dashed line in Figure 4 shows the melting temperature of the PCM.For comparison reasons a further time sequence of a volume averaged NaNO3 temperature of a heat exchanger tube design with longitudinal fins (green colored) is included in Figure 4.The numerical simulation for the heat exchanger tube with longitudinal fins was done under the same conditions and with approximately same ratios φ and ε as for the wire matrix.The data for this curve is taken from [45].It can be seen in Figure 4 that the volume averaged temperature shows no distinct flattening in the curve around the melting temperature compared to the tube with longitudinal fins.This behavior is a result of the aluminum wire matrix.
The development of the liquid fraction of the PCM over the time is depicted in Figure 5.During the first phase of the charging process (the melted region is small) the heat conduction is the main heat transfer mechanism.Furthermore, it can be seen that the melting process starts at the top of the The numerical simulation for the heat exchanger tube with longitudinal fins was done under the same conditions and with approximately same ratios ϕ and ε as for the wire matrix.The data for this curve is taken from [45].It can be seen in Figure 4 that the volume averaged temperature shows no distinct flattening in the curve around the melting temperature compared to the tube with longitudinal fins.This behavior is a result of the aluminum wire matrix.
The development of the liquid fraction of the PCM over the time is depicted in Figure 5.During the first phase of the charging process (the melted region is small) the heat conduction is the main Energies 2016, 9, 205 8 of 18 heat transfer mechanism.Furthermore, it can be seen that the melting process starts at the top of the computational domain and grows synchronously towards the outer shell and the bottom.Therefore, an asymmetric melting behavior over the height is given.This asymmetric melting process is a result of the increasing influence of the natural convection (buoyancy-driven currents) on the charging process.The molten PCM rises up along the tube wall as it gets heated and flows downwards along the solid NaNO 3 where it cools down by transferring heat to the solid region.Thus, a vortex flow arises in the melted zone, which helps to enhance the heat transfer to the solid-liquid layer of the NaNO 3 .This effect is increased by the adiabatic boundary condition on the top of the computational domain.As mentioned above, no loss of heat to the environment during the charging process is given.Thus the natural convection flow in the melted region of the PCM forces heat transfer to the upper part of the computational domain.This intensifies the described asymmetric melting process.However, the results of experimental investigations, e.g., [46,47], show the same effect.The measured temperature at the top of the storage device with vertical heat exchanger tubes increases faster compared to the temperature at the bottom.This indicates that the melting front moves from the top to the bottom.This behavior was also observed experimentally by [46,47] as well as numerically by [9,10,48].
Energies 2016, 9, 205 8 of 18 computational domain and grows synchronously towards the outer shell and the bottom.Therefore, an asymmetric melting behavior over the height is given.This asymmetric melting process is a result of the increasing influence of the natural convection (buoyancy-driven currents) on the charging process.The molten PCM rises up along the tube wall as it gets heated and flows downwards along the solid NaNO3 where it cools down by transferring heat to the solid region.Thus, a vortex flow arises in the melted zone, which helps to enhance the heat transfer to the solid-liquid layer of the NaNO3.This effect is increased by the adiabatic boundary condition on the top of the computational domain.As mentioned above, no loss of heat to the environment during the charging process is given.Thus the natural convection flow in the melted region of the PCM forces heat transfer to the upper part of the computational domain.This intensifies the described asymmetric melting process.However, the results of experimental investigations, e.g., [46,47], show the same effect.The measured temperature at the top of the storage device with vertical heat exchanger tubes increases faster compared to the temperature at the bottom.This indicates that the melting front moves from the top to the bottom.This behavior was also observed experimentally by [46,47] as well as numerically by [9,10,48].A vector plot of the velocity including the liquid fraction of the PCM is represented in Figure 6.It can be seen that the velocity profile within the melted zone is complex based on the wire matrix.At the beginning of the charging process (see Figure 6a), where the melted region is small, we find a steep gradient of the phase front.This is, on the one hand, a result of the small distance to the steel tube and, on the other hand, a result of the heat conduction which is the main heat transfer mechanism during this period of the melting process.We can also see in Figure 6a that during this first phase of melting the radial directed wires restrict the start of natural convection.Therefore the melting process of the PCM on top of each wire is slower compared to the region beside the radial wire, where the buoyancy-driven flow is unrestricted and accelerated.This effect loses its importance with increasing liquid region of the PCM, see Figure 6b.A vector plot of the velocity including the liquid fraction of the PCM is represented in Figure 6.It can be seen that the velocity profile within the melted zone is complex based on the wire matrix.At the beginning of the charging process (see Figure 6a), where the melted region is small, we find a steep gradient of the phase front.This is, on the one hand, a result of the small distance to the steel tube and, on the other hand, a result of the heat conduction which is the main heat transfer mechanism during this period of the melting process.We can also see in Figure 6a that during this first phase of melting the radial directed wires restrict the start of natural convection.Therefore the melting process of the PCM on top of each wire is slower compared to the region beside the radial wire, where the buoyancy-driven flow is unrestricted and accelerated.This effect loses its importance with increasing liquid region of the PCM, see Figure 6b.A detailed three-dimensional view of the velocity field around the wire matrix is shown in Figure 7. Based on the flow restriction of the wire matrix the liquid PCM flows around the wire and increases locally its velocity.If the downward flow of the PCM along the solid-liquid layer is restricted by a tangential wire (a part of the tangential wire is surrounded by the solid PCM, see Figures 6 and 7) then the melting process is slowed down locally in the solid region behind the tangential wire.This is a result of a lower heat transfer from the wire to the solid region behind the tangential wire based on a small temperature difference between wire and solid material.
The time evolution of the averaged velocity of the molten PCM is presented in Figure 8.In Figure 8 it can be seen that during the first period of the charging process the buoyancy-driven velocity of the molten PCM is high.This is a result of the small liquid layer and the relatively high temperature difference between the solid PCM and the tube wall.With increasing thickness of the A detailed three-dimensional view of the velocity field around the wire matrix is shown in Figure 7. Based on the flow restriction of the wire matrix the liquid PCM flows around the wire and increases locally its velocity.A detailed three-dimensional view of the velocity field around the wire matrix is shown in Figure 7. Based on the flow restriction of the wire matrix the liquid PCM flows around the wire and increases locally its velocity.If the downward flow of the PCM along the solid-liquid layer is restricted by a tangential wire (a part of the tangential wire is surrounded by the solid PCM, see Figures 6 and 7) then the melting process is slowed down locally in the solid region behind the tangential wire.This is a result of a lower heat transfer from the wire to the solid region behind the tangential wire based on a small temperature difference between wire and solid material.
The time evolution of the averaged velocity of the molten PCM is presented in Figure 8.In Figure 8 it can be seen that during the first period of the charging process the buoyancy-driven velocity of the molten PCM is high.This is a result of the small liquid layer and the relatively high temperature difference between the solid PCM and the tube wall.With increasing thickness of the If the downward flow of the PCM along the solid-liquid layer is restricted by a tangential wire (a part of the tangential wire is surrounded by the solid PCM, see Figures 6 and 7) then the melting process is slowed down locally in the solid region behind the tangential wire.This is a result of a lower heat transfer from the wire to the solid region behind the tangential wire based on a small temperature difference between wire and solid material.
The time evolution of the averaged velocity of the molten PCM is presented in Figure 8.In Figure 8 it can be seen that during the first period of the charging process the buoyancy-driven velocity of the molten PCM is high.This is a result of the small liquid layer and the relatively high temperature difference between the solid PCM and the tube wall.With increasing thickness of the liquid layer and solid temperature the temperature difference and consequently the velocity decreases.In the time Energies 2016, 9, 205 10 of 18 period between approximately 200 s up to 800 s the graph for the wired matrix shows some velocity peaks.These peaks are a result of the transversal arranged wires which represent an obstacle during the charging process for the fluid flow of the melted PCM.With the development of a small melted PCM layer between the wire and the solid PCM, the downward moving flow of the liquid PCM is divided and flows through the small liquid layer and crosses the tangential arranged wire.This fluid flow through the layer shows a high velocity which results in the velocity peak.This behavior is intensified if two or more tangential arranged wires are in the same state.With increasing liquid layer behind the tangential wire the averaged velocity decreases.
Energies 2016, 9, 205 10 of 18 liquid layer and solid temperature the temperature difference and consequently the velocity decreases.In the time period between approximately 200 s up to 800 s the graph for the wired matrix shows some velocity peaks.These peaks are a result of the transversal arranged wires which represent an obstacle during the charging process for the fluid flow of the melted PCM.With the development of a small melted PCM layer between the wire and the solid PCM, the downward moving flow of the liquid PCM is divided and flows through the small liquid layer and crosses the tangential arranged wire.This fluid flow through the layer shows a high velocity which results in the velocity peak.This behavior is intensified if two or more tangential arranged wires are in the same state.With increasing liquid layer behind the tangential wire the averaged velocity decreases.Such velocity peaks cannot be observed after approximately 800 s after simulation start.The circumstance for this behavior is given by a general lowering of the averaged PCM velocity and the flattening of the melting front.When the melting front becomes more and more flat the number of transversal arranged wires which restrain the fluid flow decreases.
With increasing of liquid volume fraction the temperature difference within the molten region decreases.This results in a further decrease of the volume averaged velocity of the PCM.
The charged power per volume PCM is depicted in Figure 9.The curve shows the typical behavior of a LHTES during the charging process.During the first phase of the melting process the charged power is very high, based on the high temperature difference between the heat transfer fluid and the solid PCM.With the development of the first small liquid layer of PCM along the outer tube surface the temperature difference between the heat transfer fluid and the liquid PCM decreases and therefore, also the heat flux decreases.After this typical rapid decrease of power (see also measurement data in [46,47]), the charged power is approximately constant up to about 1000 s after simulation start (in this analyzed case).During this period of time the volume averaged temperature increases close to the melting temperature (see Figure 4).When the volume averaged temperature of the PCM achieves approximately the melting temperature, the charged power decreases.Such velocity peaks cannot be observed after approximately 800 s after simulation start.The circumstance for this behavior is given by a general lowering of the averaged PCM velocity and the flattening of the melting front.When the melting front becomes more and more flat the number of transversal arranged wires which restrain the fluid flow decreases.
With increasing of liquid volume fraction the temperature difference within the molten region decreases.This results in a further decrease of the volume averaged velocity of the PCM.
The charged power per volume PCM is depicted in Figure 9.The curve shows the typical behavior of a LHTES during the charging process.During the first phase of the melting process the charged power is very high, based on the high temperature difference between the heat transfer fluid and the solid PCM.With the development of the first small liquid layer of PCM along the outer tube surface the temperature difference between the heat transfer fluid and the liquid PCM decreases and therefore, also the heat flux decreases.After this typical rapid decrease of power (see also measurement data in [46,47]), the charged power is approximately constant up to about 1000 s after simulation start (in this analyzed case).During this period of time the volume averaged temperature increases close to the melting temperature (see Figure 4).When the volume averaged temperature of the PCM achieves approximately the melting temperature, the charged power decreases.Such velocity peaks cannot be observed after approximately 800 s after simulation start.The circumstance for this behavior is given by a general lowering of the averaged PCM velocity and the flattening of the melting front.When the melting front becomes more and more flat the number of transversal arranged wires which restrain the fluid flow decreases.
With increasing of liquid volume fraction the temperature difference within the molten region decreases.This results in a further decrease of the volume averaged velocity of the PCM.
The charged power per volume PCM is depicted in Figure 9.The curve shows the typical behavior of a LHTES during the charging process.During the first phase of the melting process the charged power is very high, based on the high temperature difference between the heat transfer fluid and the solid PCM.With the development of the first small liquid layer of PCM along the outer tube surface the temperature difference between the heat transfer fluid and the liquid PCM decreases and therefore, also the heat flux decreases.After this typical rapid decrease of power (see also measurement data in [46,47]), the charged power is approximately constant up to about 1000 s after simulation start (in this analyzed case).During this period of time the volume averaged temperature increases close to the melting temperature (see Figure 4).When the volume averaged temperature of the PCM achieves approximately the melting temperature, the charged power decreases.

Results for the Discharging Process
In the following the results of the numerical simulation of the solidification process will be presented.The boundary and initial conditions used for the numerical investigations are described in detail in Section 2.4.
At the beginning of the discharging process the complete sodium nitrate is in the state of liquid and the circulation of the NaNO 3 starts along the whole external shell of the computational domain with an averaged velocity close to zero (see Figure 10).This circulation is a result of the natural convection which is based on the temperature difference between the sodium nitrate (initial temperature is 589.15K) and the heat transfer fluid (represented by the boundary condition of the inner tube wall surface temperature T w = 559.15K).Based on the large temperature difference during this period of solidification the buoyancy-driven velocity increases rapidly.With decreasing temperature difference within the molten region and increasing thickness of the solid PCM the averaged velocity of the sodium nitrate decreases.

Results for the Discharging Process
In the following the results of the numerical simulation of the solidification process will be presented.The boundary and initial conditions used for the numerical investigations are described in detail in Section 2.4.
At the beginning of the discharging process the complete sodium nitrate is in the state of liquid and the circulation of the NaNO3 starts along the whole external shell of the computational domain with an averaged velocity close to zero (see Figure 10).This circulation is a result of the natural convection which is based on the temperature difference between the sodium nitrate (initial temperature is 589.15K) and the heat transfer fluid (represented by the boundary condition of the inner tube wall surface temperature Tw = 559.15K).Based on the large temperature difference during this period of solidification the buoyancy-driven velocity increases rapidly.With decreasing temperature difference within the molten region and increasing thickness of the solid PCM the averaged velocity of the sodium nitrate decreases.In Figure 11, the development of the liquid fraction of the sodium nitrate during the solidification process at different time steps is depicted.As we can see from Figure 11 the solidification process starts at the bottom of the computational domain and along the outer tube surface.This result is in good agreement with numerical results presented in [9,10,48,49] and experimental work in [46,47].In all references a single vertical arranged shell-and-tube heat exchanger is investigated.For improving the heat transfer longitudinal [9,10,46,47,49] and transversal [48,49] fins are used.The numerical investigations are done three-dimensionally, except [48], which was two-dimensionally.In all numerical simulations, the volume change of the PCM during phase change due to density changes is neglected, the outer shell of the computational region is considered as adiabatic, and the mushy zone constant used for the calculations is C = 10 5 .The numerical simulation presented in [9,10,49] are done for the same boundary conditions and with approximately the same ratios φ (between 0.164 and 0.098) and ε (between 1.19 and 0.79 m −1 ) as for the present study.In [48] the supply temperature and the outlet pressure for the heat transfer fluid are used as boundary condition for the heat transfer fluid.The analyzed heat exchanger tube in [48] shows with a volume ratio of φ = 0.016 and a surface to volume ratio of ε = 1.58 m −1 a larger deviation compared to the wire matrix.In Figure 11, the development of the liquid fraction of the sodium nitrate during the solidification process at different time steps is depicted.As we can see from Figure 11 the solidification process starts at the bottom of the computational domain and along the outer tube surface.This result is in good agreement with numerical results presented in [9,10,48,49] and experimental work in [46,47].In all references a single vertical arranged shell-and-tube heat exchanger is investigated.For improving the heat transfer longitudinal [9,10,46,47,49] and transversal [48,49] fins are used.The numerical investigations are done three-dimensionally, except [48], which was two-dimensionally.In all numerical simulations, the volume change of the PCM during phase change due to density changes is neglected, the outer shell of the computational region is considered as adiabatic, and the mushy zone constant used for the calculations is C = 10 5 .The numerical simulation presented in [9,10,49] are done for the same boundary conditions and with approximately the same ratios ϕ (between 0.164 and 0.098) and ε (between 1.19 and 0.79 m ´1) as for the present study.In [48] the supply temperature and the outlet pressure for the heat transfer fluid are used as boundary condition for the heat transfer fluid.The analyzed heat exchanger tube in [48] shows with a volume ratio of ϕ = 0.016 and a surface to volume ratio of ε = 1.58 m ´1 a larger deviation compared to the wire matrix.
The numerical simulations presented in [10] are done for the same heat exchanger tube design as used for the experimental work in [46,47].The experimental work was done at approximately constant heat transfer fluid temperature (difference between supply and return temperature of the heat transfer fluid is approximately 2 K).The length of the heat exchanger tube is 3 m.
With progress of the discharging process the solidified layer of the NaNO 3 increases and the region for the buoyancy driven flow will be restricted by the growing solid sodium nitrate layer.The decreasing buoyancy-driven fluid flow takes place between the outer shell and the solid NaNO 3 .
Based on the solidified NaNO 3 layer an additional heat resistance is given.Thus heat conduction takes place within the solidified NaNO 3 region.The heat conduction will be the significant heat transfer mechanism during the solidification process and the natural convection gets less relevant.The numerical simulations presented in [10] are done for the same heat exchanger tube design as used for the experimental work in [46,47].The experimental work was done at approximately constant heat transfer fluid temperature (difference between supply and return temperature of the heat transfer fluid is approximately 2 K).The length of the heat exchanger tube is 3 m.
With progress of the discharging process the solidified layer of the NaNO3 increases and the region for the buoyancy driven flow will be restricted by the growing solid sodium nitrate layer.The decreasing buoyancy-driven fluid flow takes place between the outer shell and the solid NaNO3.Based on the solidified NaNO3 layer an additional heat resistance is given.Thus heat conduction takes place within the solidified NaNO3 region.The heat conduction will be the significant heat transfer mechanism during the solidification process and the natural convection gets less relevant.
In Figure 12 the time evolution of the volume averaged liquid fraction as well as temperature of the sodium nitrate for the analyzed heat exchanger design is illustrated.The black dashed line in Figure 12 shows the melting temperature of the PCM.From Figure 12 it can be seen that during the first phase of solidification the volume averaged temperature drops rapidly down close to the melting temperature (first approximately 250 s after simulation start).During this time the volume averaged liquid fraction decreases only by approximately 8%.This overall behavior was also observed by experimental measuring, see e.g., [46,47].After that first solidification phase, the averaged NaNO3 temperature is approximately constant around the melting temperature for a longer period of time.This is also characteristic for the solidification process of vertical arranged heat exchanger tubes, see e.g., [46,47,49], and is in contrast to the melting behavior of the wire matrix.In Figure 12 the time evolution of the volume averaged liquid fraction as well as temperature of the sodium nitrate for the analyzed heat exchanger design is illustrated.The black dashed line in Figure 12 shows the melting temperature of the PCM.From Figure 12 it can be seen that during the first phase of solidification the volume averaged temperature drops rapidly down close to the melting temperature (first approximately 250 s after simulation start).During this time the volume averaged liquid fraction decreases only by approximately 8%.This overall behavior was also observed by experimental measuring, see e.g., [46,47].After that first solidification phase, the averaged NaNO 3 temperature is approximately constant around the melting temperature for a longer period of time.This is also characteristic for the solidification process of vertical arranged heat exchanger tubes, see e.g., [46,47,49], and is in contrast to the melting behavior of the wire matrix.With decreasing volume averaged liquid fraction the volume averaged temperature of the NaNO3 is also decreasing.
If we compare Figures 4 and 12, then it can be seen that the melting process is faster than the solidification process.This shows the great importance of the natural convection during the melting process.Based on this behavior, it can be concluded that the heat exchanger tube must be designed With decreasing volume averaged liquid fraction the volume averaged temperature of the NaNO 3 is also decreasing.
If we compare Figures 4 and 12 then it can be seen that the melting process is faster than the solidification process.This shows the great importance of the natural convection during the melting process.Based on this behavior, it can be concluded that the heat exchanger tube must be designed in such a way that no restriction should be given for the natural convection flow, see also [10].
The time evolution of the averaged discharged power per m 3 NaNO 3 during the solidification process is depicted in Figure 13.For reasons of comparison a further time sequence of a volume averaged NaNO 3 temperature of a heat exchanger tube design with longitudinal fins (green colored) is included in Figure 13.With decreasing volume averaged liquid fraction the volume averaged temperature of the NaNO3 is also decreasing.
If we compare Figures 4 and 12, then it can be seen that the melting process is faster than the solidification process.This shows the great importance of the natural convection during the melting process.Based on this behavior, it can be concluded that the heat exchanger tube must be designed in such a way that no restriction should be given for the natural convection flow, see also [10].
The time evolution of the averaged discharged power per m 3 NaNO3 during the solidification process is depicted in Figure 13.For reasons of comparison a further time sequence of a volume averaged NaNO3 temperature of a heat exchanger tube design with longitudinal fins (green colored) is included in Figure 13.Apart from the very first phase of the solidification process (first approximately 30 s), where the discharged power shows a maximum based on the high temperature difference between the tube wall and the sodium nitrate, the discharged power decreases approximately linear between 200 and 100 kW/m 3 .During this period the volume averaged liquid fraction reaches the value of approximately 40%.Compared to a heat exchanger tube with longitudinal fins, the discharged power per volume is higher during the first 750 s.
The study presented in this paper was done with a small computational height.The following question needs to be evaluated: how is the influence of the computational domain height on the numerical results?
The influence of the computational height on the simulation result depends upon whether a charging or discharging process is analyzed.The simulation results of the present study as well as Apart from the very first phase of the solidification process (first approximately 30 s), where the discharged power shows a maximum based on the high temperature difference between the tube wall and the sodium nitrate, the discharged power decreases approximately linear between 200 and 100 kW/m 3 .During this period the volume averaged liquid fraction reaches the value of approximately 40%.Compared to a heat exchanger tube with longitudinal fins, the discharged power per volume is higher during the first 750 s.
The study presented in this paper was done with a small computational height.The following question needs to be evaluated: how is the influence of the computational domain height on the numerical results?
The influence of the computational height on the simulation result depends upon whether a charging or discharging process is analyzed.The simulation results of the present study as well as other studies (experimental [46,47] and numerical [9,10,48,49]) have shown, that heat conduction is the main heat transfer mechanism during the discharging process.Therefore the height of the computational domain had no significant impact on the numerical results, especially with a constant temperature boundary condition at the tube wall.The main heat transfer takes place in this case in horizontal direction and only a smaller heat flux is given in the vertical direction.In contrast to the charging process the heat is not transferred to the upper part of the computational domain.
In case of the charging process the situation is different.The main heat transfer mechanism during the charging process is natural convection.As we have seen in the results of the numerical study the natural convection flow in the molten region of the PCM shifts heat to the upper part of the computational domain (this process is intensified by the adiabatic boundary condition as discussed above).The consequence of this behavior is that a larger tube (computational domain) height results in a decrease of the heat flux from the HTF to the PCM in the upper part of the tube based on the decreasing temperature difference between HTF and PCM in this region.Thus, non-uniform heat flux distribution along the vertical heat exchanger tube dimension in direction to the outer shell is given with passing time.With increasing charging time the heat flux to the phase change material decreases.The impact of this behavior is given in the overall melting time.With other words, an increase of the tube (computational domain) height by keeping constant all other (boundary and geometrical) conditions will result in a larger melting time.But, and that is the most important knowledge, the basic physical behavior will not be effected.
The performance of the charged and discharged energy has shown that under the analyzed conditions the wire matrix is an alternative method to enhance the heat transfer to the PCM.But for the numerical analysis an ideal contact between the surfaces of the single wire as well as to the steel tube was assumed.This cannot be guaranteed in a real storage application.Therefore a higher contact resistance between the single wires must be taken into account.Especially in case of changing the cross section of the wire from quadratic to circular will increase the contact resistance.This will result in a decrease of the charging/discharging power and in an increase of the charging/discharging time.
Considering the complexity of the problem, including the multi-dimensionality of the melting and solidification process, the non-uniform heating along the heat exchanger tube (based on the temperature decrease of the heat transfer fluid in case of a single phase flow), the transient character of the process, the irregular melting patterns, the subcooling of the PCM during solidification and the heat transfer to the surroundings, further investigations must be done in future for a better understanding of the complex overall processes.In a next step, the change in volume of the PCM will be included to the mathematical models for the simulation of the melting and solidification process and these results will be validated with measurements made with the test rig at IET.Additionally, the influence of the model height on the numerical results will be investigated for a better understanding of the buoyancy-driven flow.But all these measure will lead to a considerable increase of the used CPU time.

Conclusions
In the present work, a transient numerical analysis of the melting process of sodium nitrate for a heat exchanger tube design, which consists of a steel tube and an aluminum wire matrix, has been performed.For the numerical analysis, the enthalpy-porosity formulation was used to get quantitative information about the time evolution of the melting front within the phase change material.
The numerical analysis of the melting process has shown that during the first period of the charging process, where the melted region is very small, the heat conduction is the dominant heat transfer mechanism.With increasing melted region natural convection plays the key role for the melting process.This results in a faster melting progress at the top of a storage device compared to the bottom.The analysis of the heat exchanger tube design with wired matrix indicates that no flattening is given for the volume averaged temperature of the PCM around the melting temperature.
The numerical results of the solidification process have pointed out that the dominant heat transfer mechanism during solidification is heat conduction.The investigation has also shown that the discharging process has been slower than the charging process.
For the numerical analysis, an ideal contact between the surfaces of the single wire as well as to the steel tube was assumed, but this cannot be guaranteed in a real application.Therefore, a higher contact resistance between the single wires must be taken into account which results in a larger charging and discharging time.

Figure 2 .
Figure 2. Computational domain and geometry data of the analyzed heat exchanger tube with wired matrix.

Figure 2 .
Figure 2. Computational domain and geometry data of the analyzed heat exchanger tube with wired matrix.

Figure 3 .
Figure 3. Boundary conditions used for the analyzed heat exchanger configuration.

Figure 3 .
Figure 3. Boundary conditions used for the analyzed heat exchanger configuration.

Figure 4 .
Figure 4. Chronological sequence of the volume averaged temperature and liquid fraction for the phase change material during the melting process.

Figure 4 .
Figure 4. Chronological sequence of the volume averaged temperature and liquid fraction for the phase change material during the melting process.

Figure 5 .
Figure 5. Contours of the liquid fraction of the PCM during the melting process at different time steps.

Figure 5 .
Figure 5. Contours of the liquid fraction of the PCM during the melting process at different time steps.

Figure 6 .
Figure 6.Vector plot of the liquid velocity of the molten PCM at different time steps.

Figure 7 .
Figure 7. Vector plot of the liquid velocity of the PCM 225 s after simulation start.

Figure 6 .
Figure 6.Vector plot of the liquid velocity of the molten PCM at different time steps.

Figure 6 .
Figure 6.Vector plot of the liquid velocity of the molten PCM at different time steps.

Figure 7 .
Figure 7. Vector plot of the liquid velocity of the PCM 225 s after simulation start.

Figure 7 .
Figure 7. Vector plot of the liquid velocity of the PCM 225 s after simulation start.

Figure 8 .
Figure 8. Chronological sequence of the volume averaged velocity during the charging process.

Figure 9 .
Figure 9. Chronological sequence of the volume averaged charged power during the melting process.

Figure 8 .
Figure 8. Chronological sequence of the volume averaged velocity during the charging process.
solid temperature the temperature difference and consequently the velocity decreases.In the time period between approximately 200 s up to 800 s the graph for the wired matrix shows some velocity peaks.These peaks are a result of the transversal arranged wires which represent an obstacle during the charging process for the fluid flow of the melted PCM.With the development of a small melted PCM layer between the wire and the solid PCM, the downward moving flow of the liquid PCM is divided and flows through the small liquid layer and crosses the tangential arranged wire.This fluid flow through the layer shows a high velocity which results in the velocity peak.This behavior is intensified if two or more tangential arranged wires are in the same state.With increasing liquid layer behind the tangential wire the averaged velocity decreases.

Figure 8 .
Figure 8. Chronological sequence of the volume averaged velocity during the charging process.

Figure 9 .
Figure 9. Chronological sequence of the volume averaged charged power during the melting process.

Figure 9 .
Figure 9. Chronological sequence of the volume averaged charged power during the melting process.

Figure 10 .
Figure 10.Chronological sequence of the volume averaged velocity during the discharging process.

Figure 10 .
Figure 10.Chronological sequence of the volume averaged velocity during the discharging process.

Energies 2016, 9 , 205 12 of 18 Figure 11 .
Figure 11.Contours of the liquid fraction of the PCM during the discharging process at different time steps.

Figure 11 .
Figure 11.Contours of the liquid fraction of the PCM during the discharging process at different time steps.

Energies 2016, 9 , 205 13 of 18 Figure 12 .
Figure 12.Chronological sequence of the volume averaged temperature and liquid fraction for the PCM during the discharging process.

Figure 12 .
Figure 12.Chronological sequence of the volume averaged temperature and liquid fraction for the PCM during the discharging process.

Figure 12 .
Figure 12.Chronological sequence of the volume averaged temperature and liquid fraction for the PCM during the discharging process.

Figure 13 .
Figure 13.Chronological sequence of the volume averaged discharged power per m 3 NaNO3 during the solidification process.

Figure 13 .
Figure 13.Chronological sequence of the volume averaged discharged power per m 3 NaNO 3 during the solidification process.