Evaluation of a Cooling Performance of a Typical Full Loaded Cool Storage Using Mono-scale CFD Simulation

A mono-scale three-dimensional validated CFD model was developed to predict airflow, heat and mass transfer in order to evaluate the cooling performance of a typical full loaded cool storage. The bottom of the room was subjected to fresher and cooler air with relatively higher velocity, therefore the product’s temperature reduced from top to bottom (at least from 288 to 286.21 ?k) as well as from left to right (at least from 289.91 to 283.86 ?k). The averaged absolute deviation of product’s temperature varied versus cooling time and reached to the maximum of 2.1 ?C between the 7 and 9th hours of cooling time and 0.7 ?C at the end of cooling process. Between 5.5 to 9.5 ?C difference was observed between the hottest and coolest product’s temperature during the cooling time. The model predicted the average weight loss of 1.85% for apple during the six months period in cool storage.


Introduction
Freshness and quality of fruits and vegetables as well as products saleable weight depend on heat and moisture transfer rate during air-cooling in the cool storages.The heterogeneity of product temperature during the storage may causes either the fast deterioration due to the increased respiration at higher temperature or chilling injury at lower temperature.Since numerous design and operational parameters such as flow rate and configuration of cool air blowers, the arrangement of product containers, and physical properties of products would affect the cooling performance of the system, improving it may necessitate numerous experiments.Numerical simulation is usually an economic alternative to physical experiments, although it is imperative that the reliability of the model be established through validation against experimental data (Smale, et al., 2006).Computational fluid dynamics (CFD) models have been successfully applied for analyzing the effects of design and operational parameters on cooling characteristics of fruit and vegetables and good agreements between experiments and predicted results have been reported (Tassou and Xiang, 1998;Xu and Burfoot, 1999;Foster, 2002;Hoang, et al., 2000Hoang, et al., , 2003;;Moureh, et al., 2002Moureh, et al., , 2009;;Nahor, et al., 2005;Chourasia andGoswami, 2006, 2007;Delele, et al., 2009).Due to small dimensions of slots on the fruit container's wall and large dimension of cool storage enclosure, the Simulation of airflow, heat and mass transfer inside vented fruits containers stacked in a large industrial cool storage would be a challenge.Nahor et al. (2005) developed a CFD model of the airflow, heat and mass transfer in cool room.The walls of fruit packages considered to be perforated and due to high computational costs, only empty and partially loaded (with eight packages) were modeled.In the study carried out by Moureh et al. (2009), filled boxes were taken into account as a porous media, and a separated pressure drop was used for slotted box walls.Ferrua and Singh (2009) investigated different design parameters effective on cooling rate and homogeneity of strawberry packages.Due to small slots on the packages together with large domain, only a column of packages was simulated using symmetry boundary conditions.Delele et al. (2009) applied multi-scale CFD simulation for optimization in cool storage.In multi-scale simulation the knowledge from the detailed simulation of the fruit package is used to create a model of cool storage.
Mono-scale simulation with considering directly the slotted walls and the whole large physical domain of full loaded cool storage with no simplification assumptions could cause a very fine grid which increases computational costs.Considering the availability of more powerful computers, the purpose of this study was to evaluate the cooling performance of a full loaded cool storage by developing a mono-scale three-dimensional CFD model.The cooling performance of fruits and vegetables during cooling is evaluated based on the cooling rate and homogeneity of temperature inside the products (Ferrua and Singh, 2009).Predicted air velocity, product cooling rate and weight loss were validated against experiments by the same authors (submitted for publication) and good agreements between the simulations and experiments were reported.

Cool Storage -The Simulation Domain
Cool storage modeled in this study was based on the experimental conditions of the pilot cool storage of postharvest laboratory at Dept of Food Engineering and Post harvest Technology (Karaj, Iran) with the dimensions of 4 m length, 2.6 m width and 1.9 m height.The cooling unit was located at the top center of the storage and consisted of two axial fans of 0.25 m diameter for air circulation and a finned tube heat exchanger (Figure 1).Apples (Malus domestica Borkh cv.Golden Delicious) were packed in the vented containers of 30 kg weight.The wall slots were spread over the five faces of each container and allowed air to go through about 10% of the surface (Figure 2).The container's were made of high density polyethylene (HDPE) with 0.004 m thickness and the walls modeled as conducting walls.The bulk of apple in the containers was modeled as a porous media.The enclosure was loaded with 96 containers with the dimensions of 0.54 m length, 0.4 m width and 0.3 m height.The arrangement of apple containers is shown in Figure 1.Four containers stacked with a small vertical gap of 0.06 m and there were four stacks along Z-direction with 0.06 m horizontal gap.Six columns of stacks C1 to C6 positioned alongside X-direction with 0.1 m horizontal gap.All the gaps between the container and room's sidewalls were 0.1 m.The distance between containers and roof, bottom, back and front walls of room considered as 0.46, 0.1, 0.85 and 0.25 m; respectively, considering Guidelines of settings in a refrigerated warehouse (IIR, 1966;Woolrich and Hallowell, 1970).

Continuity and Momentum Equations
The following are the governing equations in the turbulent flow regime inside the cool storage as well as fans, heat exchanger and containers of apples (Anonymous, 2003).Equations are expressed in Cartesian coordinate system.

(
) The standard k-ε turbulence model was applied in the model in which the specific Reynolds stress term ' ' u u i j was approximated as follows (Veorsteeg and Malalasekera, 1995): Bulk of apple inside the containers was considered as porous media by adding momentum source term to the standard fluid flow equations.The source term is composed of two terms: viscous and inertial loss terms.This momentum source contributed to the pressure gradient in porous cells, creating a pressure drop that is proportional to the fluid velocity in the cells and was calculated using the following equation (Anonymous, 2003): To take in to account the effect of confining walls on the viscose and inertial resistance of the porous zone, Eisfeld and Schnitzlein (2001) modified the correlation of Reichelt (1972) for confined beds, with the following coefficients for Eq.(4): Where the apple equivalent diameter ( ) was estimated using Gaskell (1992).A w and B w are analytical expressions for considering the effect of the confining walls can be found in Reichelt (1972).
Fan and heat exchanger acted as a source of momentum.The heat exchanger resistance which introduced by Zukauskas and Ulinskas (1983) applied as a momentum source: Where p  is the pressure drop due to the resistance of the cooler tubes; C is a constant factor and is equal to 0.3; z is the number of tube rows; m is a geometry parameter estimated using: Here, D h is the diameter of the tube; D 1 and D 2 are transversal and longitudinal pitch of the tubes; respectively.The fan model was formulated based on the fan characteristic provided by the supplier: , where C1, C2 and C3 are experimental coefficients, obtained from the characteristic curve and was equal to 95.97, -26.24 and 1.738; respectively.

Species transport equation
The variations of water vapor mass fraction in the enclosure during the cooling time were used to estimate the weight loss of the products.The conservation equations for chemical species of water vapor used in the model as follow (Anonymous, 2003): where W J is the diffusion flux of species water vapor, which arises due to concentration gradients and can be calculated for turbulent flow as:  is the turbulent viscosity and D t is the turbulent diffusivity.S W is the rate of the water vapor production per unit volume of products and can be estimated using: . S mAs W  . Transpiration rate per unit area of product surface ( .m ) was calculated using a lumped model, neglecting the moisture diffusion inside the products and assuming a local equilibrium at the product surface and could be expressed as (Chourasia and Goswami, 2007;Delele, et al., 2009) where vapor pressure at product surface (p vp ) and vapor pressure of ambient moist air (p va ) were estimated using: Mass transfer coefficient in bulk of product (h m ) consists of resistance of the apple skin mass transfer coefficient (h s ) and the air film mass transfer coefficient (h a ) and was estimated using : (Delele, et al., 2009).The skin mass transfer coefficient for apple was considered 0.167 × 10 -9 kg m -2 s -1 Pa -1 (Chau, et al., 1987) and the air film mass transfer coefficient was estimated using Sherwood-Reynolds-Schmidt correlation (Chau et al., 1987) ).In different superficial velocities between 0.001 to 0.1 m s -1 , the air film mass transfer coefficient was much more than skin mass transfer coefficients.Therefore, the overall mass transfer coefficient would be almost equal to skin mass transfer coefficient.

Energy Equations
Assuming a local thermal equilibrium between the air and the porous solid inside the containers, the governing energy equation for porous media was (Nakayama and Kuwahara, 2005): Where the effective thermal conductivity of the bulk (k eff ) estimated using: (1 ) where ma k and p k are moist air and products conductivity; respectively.For the containers, the source term S h took into account of metabolic rate of heat generation during respiration as well as evaporative cooling due to moisture loss during the transpiration and was calculated as the following: .
where the rate of respiratory heat generation per unit mass of product (q p ) calculated using USDA reported correlation (USDA, 1986): , f and g were: 1.69×10 -6 and 2.59; respectively (USDA, 1986).
The energy equation for the air inside the cool storage had the following form: A lumped model was used to describe the heat exchange inside the heat exchanger.Convection exchanges were assumed to follow Newton's law of heat transfer.At the cooler, the heat was removed from the air and was calculated by the following expression: where surface heat transfer coefficient at finned tube heat exchanger (h h ) was estimated using heat exchanger j ).

Boundary Conditions
Since the fluid flow and domain geometry have symmetry along X-direction, only one-half of the domain was simulated.A lumped model was considered for fan boundary condition to predict the amount of flow.However, it did not consider the turbulence due to fan blades.Apples inside containers considered to behave as a porous media and generate heat through respiration and lose moisture due to the vapor pressure deficit.The walls of containers supposed to be conducting wall.No-slip boundary conditions and zero heat flux were considered as the appropriate boundary conditions at the walls of cool storage.Standard wall functions approach was used to avoid very fine grid near the solid walls (Anonymous, 2003).

Properties and Model Parameters
All the parameters related to the apple, fan and heat exchanger were extracted from several measurements or different references, shown in Table 1.True density, bulk density, porosity, moisture content and surface area of 100 apples samples were measured in order to calculate other parameters.Bulk specific area was estimated from the correlation developed by Dullien (1979).Apple thermo-physical properties were extracted from correlations developed by USDA (1986) based on the apple moisture content.Moist air properties extracted using psychometric table developed by ASHRAE (1993).Fan parameters obtained from the manufacturer pressure-flow curve.Heat exchanger parameters were obtained from measurements and also using correlation developed by McQuiston (1981).Apple bulk viscose and inertial resistances were extracted from correlation developed by Reichelt (1972) and Gaskell (1992).

Solution Procedure and Numerical Method
The model geometry and grid were created using preprocessor Gambit v.2.3.16software.After comparing three different grid resolution results and considering the CPU time of each calculation, a uniform tetrahedral mesh size of 0.02 m for the entire domain was selected.The grid then refined near the walls to satisfy the standard wall function conditions to bound y between 30 and 300.The refinements were done gradually to avoid large errors, which could be occurred in the calculation when adjacent cell size changes more than 50%.The total numbers of cells was about 2150000 (Figure 3).A finite volume code was used for the numerical implementation of the models.Standard k-ε turbulence model was enabled based on the previous studies (Nahor, et al., 2008).The overall accuracy of the model was selected as second order upwind.The pressure-velocity coupling was ensured using SIMPLE algorithm.Prior to enabling heat and species transport and porous zone properties, the model was solved for velocity field, using small time step 0.1 s by unsteady solver applying four parallel processors.Then, time steps gradually increased up to 10 seconds.After overall flow time of 300 seconds, porous zones, heat and species transport and user defined functions for heat exchanger, heat generation and moisture loss of apples enabled in the solver.Initial temperature of 25℃ and initial zero mass fraction for vapor, patched to the entire of the domain.Convergence criteria set to absolute value of residuals 10 -3 for continuity, velocity, k, ε and vapor mass faction, and 10 -6 for energy equation.Then, the time steps increased up to 100 s and data files were auto saved after each time step.Calculations were done using a Windows XP work station with Intel Core TM 2 Quad 2.83 GHz CPU and 3.25GB of RAM.A total CPU time of about 230 hrs were needed to simulate 24 hrs of real cooling process.

Air Velocity Distributions
The general flow pattern is illustrated in Figure 4.The airflow leaving the fan accelerates along X-direction and then reaches the ceiling and detours to the front wall then flows downward and returns to the cooling unit through the containers with a little more velocity close the floor.The flow pattern through the containers on the surface, crossing the center of the fan plan, and topmost and bottommost containers are shown in Figures 5 and 6; respectively.Most of the flow was passing through the vertical gaps.Little flow passed through the horizontal gaps especially at the bottom.The air velocity close to the sidewall of cool storage is relatively more.More uniform direction velocity vectors are seen on the bottom in comparison with top containers.This is due to the closeness of the fan and the effects of air jet to the top containers.
Average velocity magnitude inside the six columns and four rows of containers is shown in Figure 7.It could be observed that the averaged velocity increased from 0.0077 to 0.0307 m s -1 (about 300% increases) from left to right in columns of containers.The similar result with less intensity was observed inside the rows where the averaged velocity magnitude increased from 0.0154 to 0.0195 m s -1 (about 25% increases) from top to bottom in rows.The maximum, minimum, average and averaged absolute deviations of predicted volume averaged velocity magnitude inside the containers were predicted 0.0408, 0.0047, 0.016 and 0.0065 m s -1 respectively.The average velocity magnitude predicted for the free space of cool storage enclosure was predicted 0.526 m s -1 .

Cooling Rate
The cooling curves for the hottest and coolest container's temperature are presented in Figure 8 and compared with the cooling curves of all containers (based on the average temperature of all containers) and averaged air temperature of cool storage.The temperature of coolest container is less than the air temperature because of the evaporation heat loss of products.
Figure 9 shows predicted product temperature after 8.7 hrs of cooling (half cooling time: time necessary to reduce the temperature difference between initial and set point temperature by a half) from an initial uniform temperature of 25℃.About 9.5 C difference was observed between the hottest and coldest product temperature at half cooling time.This difference changes during the cooling time (Figure 10) and reaches to about 5.5℃ at the end of the cooling period when the average product temperature reached to ideal maintaining temperature for apple (0℃).The product's temperature heterogeneity can be calculated as the averaged absolute deviations of temperature ( 1 ) inside the cool storage.The averaged absolute deviations of temperature versus time of cooling is shown in Figure 10.The maximum of product's temperature heterogeneity was predicted 2.1 C between the 7 and 9 hours of cooling and reduced to 0.7 C at the end of the cooling.
Averaged product temperature of the six columns and four rows of containers at half cooling time is shown in Figure 11.It could be observed that the average product temperature decreased from 289.91 to 283.86 k (about 6 k) from left to right in columns.The similar result with less intensity was observed in rows where the average product temperature decreased from 288 to 286.21 k (about 2 k) from top to bottom in rows.The bottom of the room was subjected to the fresher and cooler air with relatively higher velocity (Figure 5) that is why the product temperature reduced from top to bottom and from left to right.

Weight Loss
Average of predicted vapor mass fraction during the cooling time can be used to estimate the weight loss of the products.Due to a large computational costs needed for each simulation, the weight loss of the product after a long times can be estimated by extrapolation based on the weight loss rate at the end of the cooling process, when the product temperature attained the steady state condition.After the steady conditions, the rate of moisture transferred from the products to the air is more affected by skin mass transfer coefficient than air film mass transfer coefficient; therefore, the difference in air velocity inside the bulk of the apple in containers did not show large difference in moisture loss.The result showed the average weight loss of 0.102gr H 2 o per kg apple during 24hrs which is equal to 1.85% weight loss by extrapolation to 6 months.

Conclusion
A mono-scale three-dimensional validated CFD model was developed to simulate airflow, heat and mass transfer in a typical full loaded cool storage.Dynamic behavior of the fan, heat exchanger as well as respiration and transpiration of apple were considered in the model.A very fine grid was generated for slot ventilated walls of apple containers as well as the entire of the large domain of the cool storage.Airflow pattern and cooling curves of products was investigated.The numerical results were used to evaluate the cooling rate and temperature heterogeneity inside the cool storage.The results showed that the average velocity inside the containers increased from 0.0077 to 0.0307 m s -1 (about 300% increased ) along X-direction from left to right The similar results with less intensity were observed from top to bottom along Y-direction where velocity magnitude increased from 0.0154 to 0.0195m s -1 (about 25% increased).The bottom of the room was subjected to fresher and cooler air with relatively higher velocity, therefore the product's temperature reduced from top to bottom (at least from 288 to 286.21 k at half cooling time) as well as from left to right (at least from 289.91 to 283.86 k at half cooling time).The averaged absolute deviation of product's temperature (product's temperature heterogeneity) varied versus cooling time and reached to maximum about 2.1 C between the 7 and 9th hours of cooling time.At the end of the cooling, the product's temperature heterogeneity reduced to 0.7 C. Between 5.5 to 9.5 C difference was observed between the hottest and coolest product's temperature during the cooling time.
The model predicted the average weight loss of 1.85% for apple during the six months period in cool storage.Apple skin mass transfer coefficient 1.67×10 -10 kg m -2 s -1 Pa -1 Apple heat capacity 3639.275J kg -1 C-1 Apple heat conductivity 0.560395 W m -1 C-1 Moist air density 1.2893 kg m 3 Moist air relative humidity 90% Moist air heat capacity 1006.408J kg -1 C-1 Moist air heat conductivity 0.02397 J kg -1 C-1 Moist air viscosity 1.72×10 -5 kg m -1 s -1 Water vapor diffusivity in the air 2.1 ×10 -5 m 2 s -1 Latent heat of water (at 0 C) 2495460 J kg