Numerical Simulation of Solidification Around Staggered Tube Arrangement With Convection – Dominated

This work describes and analyses thermal storage system as a phase change problem which involves a fluid flowing inside cooled tubes in a staggered arrangement installed in a rectangular duct surrounded by a phase - change material (water). The temperature of the fluid inside the tubes is below the freezing temperature of the PCM which causes ice formation around each tubes. The problem is modeled as, two-dimensional, time dependent and convection–dominated phenomena .A finite volume numerical approach is developed and used to simulate the physical details of the problem .This approach is based on the enthalpy method which is traditionally used to track the motion of the liquid – solid front and obtain the temperature and velocity profiles in the liquid –phase. The study gives an instruction on the presentation of ice – on – coil storage tank. Results of solidification experiments are used to assess and evaluate the performance


Introduction
Rectangular power consumption peak in summer causes a serious problem as the air-conditioning load increases rapidly.It is an important issue to save energy and distribute a heat-driven airconditioning system to reduce power consumption peak.An energy saving technology , i.e. ice storage system that utilizes sensible and latent coming from phase change has been receiving increased interest recently.According to the "ASHRAE Applications ( 2000)",the use of the thermal storage systems may be an economically attractive approach for heating or cooling loads when they are of short time duration and occur cyclically ;when there is an incompatibility with the available energy source ,or its supply is limited ;when the energy costs are timedependent ; and when economical incentives are provided for the use of load-shifting equipment.
In the refrigeration field, ice banks have been employed as an energy storage solution to reduce the electrical power required during the hours of peak in the demand , shifting the consume curve to a more favorable situation.It must be noted that in many countries, the cost of industrial electricity is time dependent during the day, having higher values in the periods of high demand.It can be also mentioned the contractual demand, resulting into a fine to the user in case the energy consumption is exceeded.Aiming attempted to develop systems that employ the storage of latent heat, despite its disadvantages.Among those, it can be mentioned the inherent irreversibilities due to heat transfer, the decrease in the efficiency of the thermodynamic cycle due to the decrease in the evaporation temperature in the air-conditioning installations, the increase in the physical space that is necessary to accommodate the ice bank, and uncertainties in the evaluation of the heat transfer between the PCM and the system working fluid.the air-conditioning installations, the increase in the physical space that is necessary to accommodate the ice bank, and uncertainties in the evaluation of the heat transfer between the PCM and the system working fluid.Water in the solid phase has been employed as the storage substance due to its properties and characteristics: high melting latent heat (335 kJ/kg), melting temperature of 0˚C at the atmospheric pressure, accessibility, low cost, no pollutant effects, and chemical stability.The ice ice bank consists of thermally insulated tanks that are inserted into a refrigeration loop.In such systems, a limit must be impose to the ice thickness formed on the piping, since if it Examples of works that investigated the phase change occurring in the outside of the tubes are "Ho and Chen (1986)", Zhang and Faghri (1996)","Abugderah and Ismail exceeds a given value there will be a decrease in the heat transfer during the discharge process, caused by the ice thickness and the flow obstruction.

(2000)","Stampa et al.(2002)".
The current work is a latent "cold storage" modulus, which consists of a bank of tubes containing fluid inside, subjected to a crossed external flow To simulate the phase change process, this paper proposes a simple numerical scheme based on the enthalpy method.The numerical results are verified with the experimental data.The spatial and temporal discriminations are achieved in the context of the finite volume method and the fully scheme, respectively.

Description of The Problem
The Thermal Energy Storage TES unit being investigated consists of a copper tubes arranged in staggered manner installed in the vertical position in a rectangular duct as surrounded by phase change material (water).The tubes were arranged in three rows with three tubes in each row, as shown in figure ( 1 ) .
In freezing process, the temperature of the fluid (glycol) running through the tubes is below the freezing temperature of the PCM, causing solid layers to form and propagate from the tubes outside surfaces, until steady-state condition is reached.Where the steady-state condition is defined as stoppage of ice formation process resulting from heat flow balance at the ice-water interface.The freezing process starts, when at time = 0, the wall temperature of the cooled tube drops to T w = - Eng.&Tech.Vol.26,No.2,2008Numerical Simulation of Solidification Around Staggered Tube Arrangement With Convection -Dominated .191 10ºC, while the temperature range of the flowing water are T ∞ = 2 -3ºC, and the freestream velocity range of water upstream of tubes are U ∞ = 0.006 -0.04 m/s.The problem is modeled as conjugate phase-change, convection-dominated using enthalpy formulation for two-dimensional transient, with solidification / melting of a PCM having a single fusion / melting temperature, T m .In the numerical runs, the solution obtained for both transient process and a steady-state condition.

CFD Model
The following model is adopted in the numerical analyses : 1.The fluid is a Newtonian fluid, and Boussinesq approximation in which density is considered as a function of temperature and density in the other terms including the convective is assumed to be constant.2.The flow is two dimensional , laminar , and incompressible.3.The irreversible dissipation from kinetic energy to thermal energy is neglected.4.The momentum field is subjected to no-slip boundary conditions at the wall.The thermo physical properties (C P , k , µ )of the material are constant, but may differ for the liquid and solid phases.

Governing Equations
When the assumptions in the CFD model are applied, the governing equations for the transient convective heat transfer are continuity equation, momentum equation, and energy equation as follows.The momentum equation is divided into x and y directional components : In this method, the absorption and evolution of the latent heat during phase change leads to the modification of the energy equation because the interface is not tracked, and thus the interface conditions are not imposed explicitly.The method relies on the enthalpy formulations, which employs the enthalpy as a dependent variable in the energy equation rather than the temperature.The enthalpy formulation defines the liquid mass fraction ( f ) as enthalpy will be : [ 8 ] cT fL h + =

………. ( 5 )
The ratio of the liquid mass to the total mass in a given computational cell.If ( h s and T m ) are set to the reference enthalpy and Temperature , respectively , the specific The heat capacity ( c ) may vary with the phase.The liquid mass fraction can be obtained from the enthalpy :

Numerical Methodology
For the integration of the governing equations, the finite volume method is applied , as describe in " Patanker (1980) " .
The governing equations is discretize by a Cartesian uniform and fixed grid with staggered arrangement.In the staggered grid, the velocity components are computed for the points lying on the control-volume faces.All other variables including pressure are calculated at the grid points "Vasilios etal ( 2003)".A computational mesh 51× 31 was employed after grid independence tests were performed.In the equation , the diffusion coefficient depends whether the considered region is the liquid , interface or solid .The systems of equations resulting from the integration of the differential equations is solved by the TDMA method.

Numerical Method
The discretized energy equation for two dimension in the finite volume formulation (Patankar,1980) can be expressed as: The term relating to the liquid fraction separate the non-linear behavior associated with the phase change into a source term.The superscript denotes the value at the previous time step.
If the discretized equation ( 10) is solved at the n-th iteration step the enthalpy obtained with physical properties assumed at the n-th iteration step satisfies the energy conservation.Then, the enthalpy and the liquid fraction can be obtainedfrom equations (5&6), respectively.This procedure enables the energy contained in the cell to redistribute so that the excessive (or deficient) energy can be stored into (or retrived from) latent heat rather than spurious sensible heat.Also, the temperature can be re-estimated according to the new mass fraction.
Which gives a better estimation of the temperature for the next iteration.
The whole set of the governing equations do not need to be solved during this procedure.It should be noted that during the inner iteration, the temperatures at neighboring cells are not changed but the updated temperature is a good estimation for the next iteration.The proposed algorithm seems to be very effective due to its simplicity and low computational cost.Further more, this can be readily adaptable to any numerical scheme designed for computational efficiency.This algorithm always ensures the energy conservation at the phase changing cell.Inner iteration are stopped when the residual is reached by a factor set to 10 -10 for the u ,v , h, and pressure correction equation The convergence criterion for outer iterations is based on a total number of outer iteration rather than a prescribed value for the residual.The choice of total number of outer iterations is based on the monitoring of the outer iteration errors (ratio of the norms of the difference between two consecutive iterates | Ф n+1 -Ф n | and the difference between the first two iterates.Where Ф stands for any of u, v ,h or P). .The outer iterations are stopped when the residual is reached a factor set to 5×10 -3 .

Numerical Parameters
The three dimensionless parameter describing the problem are the Reynolds number (Re), Prandtl number (Pr), and Stefan number (Ste) defined as: Where ∆T = T m -T W is the difference between the tube wall temperature T W and the melting temperature.In the above definition u, d, v, α, c , L denote the velocity of water following through the duct (TES), diameter of tube, kinematics viscosity, thermal diffusivity, specific heat of constant pressure, latent heat of fusion.A modified Boussinesq approximation has been used, that is a non-linear density variation has been considered.The density was approximated using a fourth order polynomial obtained by fitting it to the data range of [0ºC-20ºC] given by : ρ = 999.840281+ 0.0673268.T -0.00894484.T 2 + 8.7846287×10 -5 .T 3 -6.6213979×10- 7 .T 4 …………(15) Where the temperature T is given in degrees Celsius and ρ in kg/m 3 .
The remaining physical properties of the liquid are assumed to be functions of temperature as given by "Karlekar & Desmond (1982)".Eng.&Tech.Vol.26,No.2,2008Numerical Simulation of Solidification Around Staggered Tube Arrangement With Convection -Dominated .

193
For the temperature range investigated (i.e. from 0ºC to 20ºC) this resulted in the following variations:( see table 1)

Results and Discussion
The problem chosen for demonstrating the present technique is that of forced convection of water flowing around tubes bank in a staggered arrangement with tubes wall maintained at a temperature below the freezing point of water.
Radiation and free convection effects were neglected.
Two sample results are chosen to be presented in this study for small and large Reynolds number( Re d = 500,θ c = 4.3& Re d = 1400,θ c = 5.0).respectively.These results in form of velocity vectors and solid-liquid temperature contours.The results are presented for different times intervals ,i.e.t = 40min, 100min, 190min, and steady-state condition; its required 360 -720min to obtain a steady-state condition in this range of the present study.The predicted behavior during the transient process is summarized below.

Solid-Liquid Temperature Plots
( 1 ) At t = 40min, the profile of the temperature contours of ice was concentric around each individual tube as shown in figures (2a &3a).
A close observation revealed that the temperature profile and the ice thickness were effected by a stagnant flow resulting from the existence of neighboring stagnation region and the local ice thickness at the forward stagnation region become thick.This phenomenon was due to a decrease of heat transfer rate at the stagnation region.
( 2 ) at t=100 min, as seen in figures (2b,3b) are the result obtained under conditions that produce thicker ice layers than figure (2.a), because the water temperature decreases as it flowed across the cooled tubes for longer time.
( 3 ) at t =190 min, figure (2.c) and figure (3.c) the tubes at the same row down stream were linked by the ice layer, and the amount of linked ice between the tubes were comparatively small and a flow separation occurred on the linked ice surface between the tubes, as will be shown later.( 4 ) at steady-state condition, the figures (2.d & 3.d) show the completely linked ice layer between the tubes at the same row.Under the conditions that produce thicker ice layer, a flow separation occurred on the linked ice surface between the tubes.The ice-water interface of figures (2. d & 3.d) shows the sinuous shape of a wavelength corresponding with the tube pitches.The ice layer become thicker downstream, which implies that the ice ice pattern, the symmetric ice layer was layer was also effected by a temperature boundary layer in the water stream.For the linkedalways observed.The freezing process is required a longer period with increasing Re d , and the ice layer thickness is decreasing with increasing Re d .It suggests that the ice growth rate in the freezing process decreases due to larger heat transfer by the separated and reattached flows.

Velocity Vectors
( 1 ) 40 min ≤ t ≤ 100 min , during this time , the profile of the solid is approximately cylindrical shape.As can be seen from figures 4 (a&b), 5(a&b), the mathematical model for the ice layer forming around each tubes gives velocity equal to zero.The velocity vectors adjusted by the initial flow area become approximately parallel and their profile develops after the first tube of the symmetry axis.There is retardation and solidification occurring in the proximity of the tube surface.The resulting fully developed flow profile is not of the parabolic laminar shape but it is flatter along the core.This may be the result of the thickening of the solidification band near down stream tube; i.e. the phase-change effects near the outlet are strong even in the core, where a mushy region exists.
( 2 ) 160 min ≤ t ≤ steady-state.See figures (4c&d , 5c&d) , during this time, the high axial flow velocity dominate leading to a more parabolic velocity profile.The core flow has accelerated to roughly four times its initial axial velocity value.In this case, an obvious two-dimensional profile is present with considerable negative radial components(i.e. in y-direction).These facts would suggest that the rate at which the ice thickness increases with time at the tubes downstream, would decrease on nearing blockage due to high convection there.Looking along constant radial values, vectors near the solid at the entry are very small and are hence easily absorbed into the ice.This process takes longer further downstream, above last tube the outlet vectors in the solidifying band are very high in comparison, indicating a thick, steady -state mushy region.

Validationofthe Proposed Methodology
The ultimate test for any mathematical model is to compare its predictions with results from accurate physical experiments.In this case the model was run until steady state was reached for the laminar, flow of water, in order to compare predictions with data found in Hirata and Matsui [ 4 ] .In these experiments ( for which a schematic view of the apparatus is given in figure ( 6), an anastigmatic lens were used to photographed the ice layer to measure its thickness and the coordinates of the ice-water interface were measured too.The predicted results were compared with the experimental results, for the conditions as stated in table ( 3 ).

Conclusion
The convection-dominated solidification phenomena around the cooled tubes in staggered arrangement installed in rectangular duct with water flowing across them is investigated numerically on the basis of the fixed-grid formulation.The phase change process and the velocity suppression is modeled by the enthalpy method.A general mathematical model has been developed for predicting the freezing process applicable to any geometrical or heat flow conditions.It appears that the enthalpy method is general and conforms well with general controlvolume techniques for fluid-flow and heat transfer simulations.The present model is verified with the experimental data of Hirata and Matsui .There exists an a fairly good agreement between the two set of results.

References
Eng.&Tech.Vol.26,No.2,2008Numerical Simulation of Solidification Around Staggered Tube Arrangement With Convection -Dominated .194 Fig. (6) Schematic illustration of experimental apparatus[4 ]. Figure ( 7 ) shows the experimental and numerical results of mean ice thickness ( δ m ).It may be seen that for both low and high Reynolds numbers the present method shows very fair agreement for any practical application.Figure ( 8 ) shows the comparison for temperature contour for (Re d =500 and θ c =4.3) at steady -state condition for both computational and experimental, the symmetric ice layer was observed.

Table ( 1
) Physical properties of water.