Printer Friendly

A comparative study of the effect of initial turbulence on the performance of an open vertical refrigerated multi-deck--a numerical study and its experimental validation.


Open refrigerated display cases (ORDC) are common used in supermarkets to maintain the food products at setting temperatures. In typical design of ORDC, having a correct temperature setting and an appropriate controlling entrainment or infiltration rate are the most crucial deign parameters as far as energy saving is concerned. The suitable design parameters can be obtained either numerically or experimentally. As far as cost is concerned, CFD is regarded as the best design tool because it is generally fast and reliable. Hence, there were many researchers adopt CFD tools to optimize the design parameters and to improve the performance of open display cases [1, 3-6, 8-14].

Amid the CFD examinations of relevant design parameters under vertical design situation, most of them were performed in two-dimensional conditions [1, 3-5, 9-14], and only very few literatures (D'Agaro et al. 2006; Foster et al. 2005) were using three-dimensional simulation. The simulations [1, 5-6, 8-14] were generally compared with experimental tests but its agreement with the experimental data was usually qualitative rather than quantitative. For the horizontal refrigerated display cases, the simulation by Cui and Wang (2004) is the only one that is in line with measured temperature profile at the outlet of the air curtain. For the vertical refrigerated display cases, unfortunately, none of the existing literatures can provide a quantitative agreement with the experimental data. Notice that in vertical refrigerated display case the momentum force of air jetting from DAG and buoyancy force counteracts with each other, resulting in imbalance during simulations. Yet the importance of buoyancy force can be made clear from the numerical simulation carried out by Bhattacharjee and Loth (2004) and was confirmed experimentally by Field and Loth (2006). Without considering the influence of buoyancy force, it was not surprised that the simulation results [9-14] were not consistent with the measurements.

There were some studies associated with the influence of buoyancy force. Cortella et al. (2001, 2002) included the buoyancy force in the simulation but the simulation was unable to extend to 3D situation for the model was based on stream-vorticity formulation. D'Agaro et al. (2006) considered both buoyancy force and 3D effects, and showed that the cabinet performance was highly dependent on 3D flow structures. However, they mentioned that it was quite difficult to achieve full agreement with experimental data due to the uncertainty in the experimental boundary conditions, especially in the velocity distribution at the curtain outlets, small imperfection in the actual cabinet geometry, or modeling issues in radiative heat transfer.

The forgoing survey suggests that a 3D simulation model taking into account the effect of buoyancy force accompanied with and the exact boundary conditions would resolve the present inaccurate agreements between simulations and experiments. Hence it is the purpose of this study to include these effects for simulation. The simulation is made in a typical ORDC having six decks (Fig. 1). In addition, the experimental verification is also carried out to compare with the simulation result.



The present 3D formulations of the mathematical modeling consist of equations of mass, momentum, and energy conservation with turbulence being expressed by the k-[epsilon]-E model. The Boussinesq approximation for the buoyancy force is adopted in momentum equation. Schematic of the front view and side view of open refrigeration display cabinet is shown in Fig. 1, and the computational domain is shown in dashed purple line.The simulation is conducted via a commercially available CFD package (CFD-ACE+2007 from ESI group). Details of the formulation are described as follows.


The governing equations of turbulence flow with heat transfer are solved by the mass conservation, averaged momentum equation, and averaged energy equation with the k-[epsilon]-E two-equations turbulence model which is based on the Boussinesq eddy viscosity assumption.

* the mass conservation equation,

[[partial derivative]/[[partial derivative][x.sub.j]]]([rho][u.sub.j]) = 0 (1)

* the averaged momentum equation,

[[partial derivative]/[[partial derivative][x.sub.j]]]([rho][u.sub.i][u.sub.j]) = -[[[partial derivative]p]/[[partial derivative][x.sub.i]]] + [[partial derivative]/[[partial derivative][x.sub.j]]][[mu]([[partial derivative][u.sub.i]]/[[partial derivative][x.sub.j]] + [[partial derivative][u.sub.j]]/[[partial derivative][x.sub.i]])] + [[partial derivative]/[[partial derivative][x.sub.j]]](-[rho][bar.[u'.sub.i][u'.sub.j]]) + [F.sub.buoyancy] (2)

[F.sub.buoyancy] is the buoyancy force from the difference in density subject to temperature change.

* the averaged energy equation

[[partial derivative]/[[partial derivative][x.sub.j]]]([u.sub.j]T) = [[partial derivative]/[[partial derivative][x.sub.j]]]([alpha][[partial derivative]T]/[[partial derivative][x.sub.j]]) + [[partial derivative]/[[partial derivative][x.sub.j]]](-[[bar.u'].sub.j][theta]) (3)

* k-[epsilon] model

Two-equation k-[epsilon] model are used. Turbulence kinetic energy k and dissipation rate [epsilon] are computed through transport-diffusion equation:

[[partial derivative]/[[partial derivative][x.sub.j]]]([rho][u.sub.j]k) = [[partial derivative]/[[partial derivative][x.sub.j]]]([[mu].sub.t]/[[sigma].sub.k][[partial derivative]k]/[[partial derivative][x.sub.j]]) + [[mu].sub.t]([[partial derivative][u.sub.i]]/[[partial derivative][x.sub.j]] + [[partial derivative][u.sub.j]]/[[partial derivative][x.sub.i]])[[[partial derivative][u.sub.i]]/[[partial derivative][x.sub.j]]] - [rho][epsilon] (4)

[[partial derivative]/[[partial derivative][x.sub.j]]]([rho][u.sub.j][epsilon]) = [[partial derivative]/[[partial derivative][x.sub.j]]]([[mu].sub.t]/[[sigma].sub.[epsilon]][[partial derivative][epsilon]]/[[partial derivative][x.sub.j]]) + [C.sub.1][[mu].sub.t][[epsilon]/k]([[partial derivative][u.sub.i]]/[[partial derivative][x.sub.j]] + [[partial derivative][u.sub.j]]/[[partial derivative][x.sub.i]])[[[partial derivative][u.sub.i]]/[[partial derivative][x.sub.j]]] - [C.sub.2][rho][[[epsilon].sup.2]/k] (5)

The generalized Boussinesq eddy viscosity is adopted to represent the Reynolds stress equation, -[rho][bar.[u'.sub.i][u'.sub.j]], and the Reynolds heat flux equation, -[bar.[u'.sub.j][theta]] as shown in Eq. (2) and (3). Thus, the Reynolds stress (-[rho][bar.[u'.sub.i][u'.sub.j]]) and the Reynolds heat flux (-[bar.[u'.sub.j][theta]]) are modeled on the basis of the eddy viscosity and eddy diffusivity models, respectively. In this model, the Reynolds stress and the Reynolds heat flux can be written as:

-[rho][bar.[u'.sub.i][u'.sub.j]] = [[mu].sub.t]([[partial derivative][u.sub.i]]/[[partial derivative][x.sub.j]] + [[partial derivative][u.sub.j]]/[[partial derivative][x.sub.i]]) - [2/3][[delta].sub.ij][rho]k (6)

-[bar.[u'.sub.i][theta]] = [[alpha].sub.t][[[partial derivative]T]/[[partial derivative][x.sub.i]]] (7)

where [[mu].sub.t] is the turbulent eddy viscosity and [[alpha].sub.t] is turbulent eddy thermal diffusivity, and can be expressed as

[[mu].sub.t] = [rho][C.sub.[mu]][[k.sup.2]/[epsilon]] (8)

[[alpha].sub.t] = [[C.sub.[mu]]/[Pr.sub.t]][[k.sup.2]/[epsilon]] (9)

Where [Pr.sub.t], ranging from 0.8 to 1.3 in this study, is the turbulent Prandtl number. The five empirical constants of Eqs. (4-9) are

[C.sub.[mu]] = 0.09, [C.sub.1] = 1.44, [C.sub.2] = 1.92, [[sigma].sub.k] = 1.0, [[sigma].sub.[epsilon]] = 1.3 (10)


For problem having small or medium density gradient, the governing equations can be simplified by invoking the Boussinesq approximation (Chen and Jaw, 1998). The Boussinesq approximation for buoyancy force [F.sub.buoyancy] in Eq. (2) takes the form

[F.sub.buoyancy] = [beta][rho][[vector].g](T - [T.sub.amb])[[delta].sub.ij] (11)

In Eq. (11), [beta] is the volume expansion coefficient and can be expressed as

[beta] [approximately equal to] -[1/[rho]][[[DELTA][rho]]/[[DELTA]T]] = -[1/[rho]][[[[rho].sub.amb] - [rho]]/[[T.sub.amb] - T]] (at constant p) (12)

As is well known, the Reynolds number (Re) is the most important dimensionless parameter for forced convection whereas the Grashof number characterizes the natural convection condition. The Grashof number represents the ratio of the buoyancy force to the viscous force and is termed as

Gr = [[g[beta](T - [T.sub.amb])[L.sub.C.sup.3]]/[v.sup.2]] (13)

The Reynolds number is defined as

Re = [[[rho]v[L.sub.c]]/[mu]] (14)

In mixed convection both modes may be comparable to each other. The relative importance of each mode of heat transfer is determined by the value of Gr/Re. The width of the DAG (w = 70 mm [2.8 in.], in Fig. 1a) is chosen as the characteristic length ([L.sub.c]).


The velocity of DAG is prescribed using the same measured velocity profile, yet the turbulent kinetic energy (k) of DAG can be estimated from:

k = [3/2][(TI x u).sup.2] (15)

Where TI is the turbulence intensity, and it is defined as the ratio of the root-mean-square of the fluctuation velocity, u', to the mean flow velocity, u

TI = u'/u (16)

The dissipation energy [epsilon] of DAG is initially guessed as 0.01k.

Except the door region in Fig. 1(b), the adiabatic (ambient temperature T = 298K [76.7[degrees]F]) and no-slip boundary conditions (u = v = w = 0) are given on the wall surfaces that are shown in purple dashed line in Fig. 1 Fixed pressure (p = [p.sub.a]) and the ambient temperature T = 298K (76.7[degrees]F) are given at the outlet boundary of the door region. In this study, the air channel from the evaporator to the discharge air grill (DAG) is neglected. To amend this minor difference, a negative pressure is given at the outlet of the RAG (Fig. 1b). The magnitude of the pressure is then calculated when the flow rate is balance between the DAG and RAG.


For the open refrigerated display case, the cold air is provided through an inlet jet called the discharge air grille (DAG) located at the top front of the unit as appeared in Fig. 1b. The cold air jet, termed as air curtain, acts as an invisible barrier between cabinet and outside warm air. The continuous flow of outside warm air into the air curtain and its subsequent mixing with cold air is called entrainment, and some of the warm air becomes infiltrated. When outside warm air has infiltrated through the return air grill (RAG), it will increase the air temperature and impose a cooling load on the refrigeration system. Obviously, to keep a minimum infiltrated rate from outside warm air is a critical issue in energy conservation. The infiltration rate, [[alpha].sub.e], was usually used to indicate the performance of the open refrigerated display case, and is defined as (Bhattacharjee and Lothe, 2004; Navaz et ak, 2005; D'Agaro et al., 2006)

[[alpha].sub.e] = ([T.sub.capture] - [T.sub.DAG])/([T.sub.amb] - [T.sub.DAG]) (17)

where [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. In the ideal air curtain, the inflitration rate would be zero ([a.sub.e] = 0).


The vertical open fronted chilled multi-deck cabinet is used for experimental verification. The size is 0.91 m (2.99 ft) L x 0.6 m (1.97 ft) W x 1.9 m (6.23 ft) H. Tests were conducted according to the EN441 Standard requirements. An environment chamber (5 x 3.6 x 2.8 m [16.4 x 11.81 x 9.19 ft] in size) maintains the ambient conditions at 25[degrees]C (77[degrees]F) and at a relative humidity of 60%. The corresponding air curtain velocities are measured by a Sierra-600 hot wire anemometer having an accuracy of 1% of the measured value in the range of 0-2 m/s (0-6.56 ft/s). The air curtain temperature is measured by T-type thermocouple (accuracy of [+ or -]0.1[degrees]C [[+ or -]0.18[degrees]F]) with calibrated range from -20[degrees]C to 50[degrees]C (-4[degrees]F to 122[degrees]F).

Figure 1 is the schematic of the side view and the front view of the open cabinet. A thermocouple and an anemometer are installed in the middle of discharge air grill (DAG) to record the time history diagram of temperature and velocity. Six thermocouples are also installed in the middle of each shelf, and the positions are shown designated as S1 to S6 in Fig. 1.


The computational domain and grid system in x-y plane, x-z plane, and y-z plane are plotted in Fig. 2(a)-(c). The total cells number is 174864. The effect of evaporator area and the air channel are neglected in simulation, so that the cold air is provided from the discharge air grille (DAG) located at the top front of the unit (Fig. 1b), and the inlet velocity profile of DAG is given from experimental measurement. The measured velocity profile along the x-direction is shown in Fig. 3 and this is adopted as the entrance velocity profile for calculation input with z-direction being assumed uniformly. The average velocity is 1.3 m/s (4.27 ft/s) and the Reynolds number Re is equal to 6314. The temperature of the measured cold air of DAG is 0[degrees]C (32[degrees]F) and the Grashof number is near 3.5x[10.sup.6].




Figure 4a is the temperature distribution on the media vertical section (x-y plane at z = 0.42 m [1.38 ft]) with turbulence intensity being 5%, 10%, and 25%, respectively. By checking the temperature distribution from the different shelves, it is found that a significant rise of temperature accompanied with higher turbulence intensity (Fig. 4a). The width of the cold air curtain becomes thicker as the turbulence intensity is increased from 5% to 25% (Fig. 4b). Figure 4c illustrated the velocity magnitude distribution of the x-z plane at y = 0.62 m (2.03 ft) from the bottom, which is located near the sixth shelf. It is obvious that the air is infiltrated or entrainment from the outside of the cabinet, even though the turbulence intensity is as low as 5%. The infiltrated or entrainment is increased with the turbulence intensity further increased, so that the infiltration rate is increased from 0.4 to 0.5 as turbulence intensity is increased from 5% to 25% (Fig. 5).




As mentioned in the introduction, the three dimensional effect may play a decisive role. Hence, it is interesting to examine the qualitative difference amid the two-dimension (2D) and three-dimension (3D) simulation. The relevant comparisons for the temperature and velocity along the open front of DAG to RAG (position L1 to L6) associated with different turbulence intensities are plotted in Fig. 6. The y' (0 m [0 ft]) to 1.2 m [3.94 ft]) is the distance from RAD (y'=0 m [0 ft]) to DAG (y'=1.2 m [3.94 ft]) along the L line appeared in Fig 1(b). The temperature rise for TI = 5% from DAG (y' = 1.2 m [3.94 ft]) to RAG (y' = 0 m [0 ft]) is 5[degrees]C (9[degrees]F) via applying the 2D simulation (Fig. 6a), yet it is increased to around 6[degrees]C (10.08[degrees]F) by the 3D simulation. Analogous results are seen for the case of TI = 25%, the temperature rise of 3D simulation is also significantly higher than that of 2D simulation. The difference is attributed to the swirled flow structure caused by the 3D flow field in x-z plane as shown in Fig. 4c, through which one can expect a better mixing. The 2D model assumes a symmetry boundary (or infinitely) in z-direction that eventually eliminate the formation of the swirled flow structure in x-z plane.


For the velocity distribution at the open front of DAG to RAG (Fig. 6b), the flow inertia dwindled as it flows from DAG to RAG. It is not surprised that the rise of turbulence intensity gives way to larger energy dissipation, thereby leading to a sharp decline of velocity from DAG to RAG. Yet the phenomenon becomes more pronounced at a higher turbulence intensity (TI = 25%). This phenomenon is observed in both 2D and 3D simulation, but the effect of 3D simulation considerably exceeds that of 2D, indicating the effect of outside warm air entrainment or infiltration from the 3D simulation is reinforced by the swirled flow structure. This can be further made clear from Fig. 5 where the difference in infiltration rate is generally above 30%.

The temperature distribution along the center of the shelves (the line S in Fig. 1) is plotted in Fig. 7. The y' (0 to 1.2 m [0 to 3.94 ft]) is the distance from the bottom shelf (S6, y'=0 m[0 ft]) to top shelf (S1, y'=1.2 m [3.94 ft]) (Fig. 1a). Although the temperature profile along the center of the shelves shows a step shape (kept at the same in the same shelf), but the value is rather similar to it along the open front of DAG to RAG. It is notable that the correct food storage temperature is an important parameter for an open vertical refrigerated display case, and the temperature inside the shelf is significantly related to air curtain. Thus, the design of efficient air curtain is very important for ORDC. Based on the present simulation, it can be found that reducing the turbulence intensity of DAG is an effective method to minimize the temperature rise at the bottom shelf. According to the comparison of 2D and 3D simulation (Fig. 7), the temperature rise calculated by 3D is higher than 2D. The 3D effect can not be neglected in the ORDC simulation.



Figure 8 shows the comparison of the averaged temperature at position S1 to S6 between the 3D simulation model and measured data. The positions S1 to S6 are in the middle positions inside each shelves as shown in Fig. 1. Since the outside warm air infiltration or entrainment from air curtain, the measured temperature inside the shelves are increased from the first shelf to the sixth shelf. The temperature difference between the first shelf and the sixth shelf is near 8[degrees]C (14.4[degrees]F). The computed temperature agreed well with the experimental data at a turbulence intensity of 15%. Invoking the effects of 3D and buoyancy force into the present k-[epsilon]-E turbulence model can well evaluate the performance of ORDC.



This study adopts a three-dimensional k-[epsilon]-E turbulence model in association with the effect of buoyancy force to examine the performance of an open refrigerated display case having vertical installation. The simulation is compared with the experimental measurements. It is found that the effect of turbulence intensity or the influence of 3D flow structure is rather important. The infiltration rate is increased as the turbulence intensity is increased from 5% to 25%, yet the infiltration rate for 3D simulation exceeds that of 2D results by approximately 20~30%. This is due to the influence of swirled flow structure. This can be made clear from the temperature difference between shelves. In summary of the comparison of the temperature, velocity distribution, and the infiltration rate, the effect of 3D can not be ignored for it provides better flow mixing and gives rise to higher temperature. To reduce the temperature difference between the first and sixth shelf, the present calculation suggests that reducing the turbulence intensity would be a very effective method. When the turbulence intensity is reduced from 25% to 5%, the temperature difference between the first and sixth shelf can be reduced from 11[degrees]C to 5[degrees]C (19.8[degrees]F to 9[degrees]F) and the entrainment rate can be reduced by 25%. In addition to the numerical simulation, the calculated temperature is validated by experimental measurements. The calculations with a turbulence intensity of 15% are in line with the measurements.


The authors would like to express gratitude for the Energy R&D foundation funding from the Bureau of Energy of the Ministry of Economic of Taiwan.


[C.sub.p] = heat capacity, kJ/kg*K (Btu/lb*[degrees]F)

[F.sub.buoyancy] = buoyancy force

g = gravitational acceleration, m/[s.sup.2] (ft/[s.sup.2])

Gr = Grashof number

k = turbulence kinetic energy, k

[k.sub.t] = thermal conductivity, W/m*K (Btu*ft/h*[ft.sup.2]*[degrees]F)

[L.sub.c] = characteristic length equals the width of the DAG, w

p = mean pressure, N/[m.sup.2] ([lb.sub.f]/[ft.sup.2])

[Pr.sub.t] = turbulent Prandtl number

Re = Reynolds number

T = mean temperature, K ([degrees]F)

TI = turbulence intensity

u' = root-mean-square of the fluctuation velocity,

u' = [square root of [[[i = 3].summation over (i = 1)][u'.sub.i.sup.2]]]

u = mean flow velocity, m/s (ft/s)

[u.sub.i] = mean velocity, m/s (ft/s)

[u.sub.i]' = velocity fluctuating value, m/s (ft/s)

w = the width of the DAG, m (ft)

[mu] = viscosity, kg/ms (lb/ft*s)

v = kinetic viscosity, [m.sup.2]/s ([ft.sup.2]/s)

[[mu].sub.t] = turbulent eddy viscosity

[alpha] = thermal diffusivity,

[[alpha].sub.e] = entrainment rate

[[alpha].sub.t] = turbulent eddy thermal diffusivity

[beta] = volume expansion coefficient, 1/K (1/[degrees]F)

[epsilon] = dissipation rate

[theta] = temperature fluctuating value, K ([degrees]F)

[rho] = density, kg/[m.sup.3] (lb/[ft.sup.3])

[[sigma].sub.k] = effective Prandtl numbers of kinetic energy

[[sigma].sub.[epsilon]] = effective Prandtl numbers of dissipation rate


amb = ambient

DAG = discharge air grill

RAG = return air grill


Bhattacharjee, P., and E. Loth. 2004. Entrainment by a Refrigerated Air Curtain Down a wall. ASME J. of Fluids Engineering 126: 871-79.

Chen C. J. and S. Y. Jaw. 1998. Fundamentals of turbulence modeling. Chapter 6.3 Boussinesq approximation, London: Taylor & Francis book, Ltd.

Cortella G. 2002. CFD-aided retail cabinets design. Computers and electronics in agriculture 34: 43-66.

Cortella G., M. Manzan, G. Comini. 2001. CFD simulation of refrigerated display cabinets. Internal Journal of Refrigeration 24: 250-260.

Cui J., S. Wang. 2004. Application of CFD in evaluation and energy-efficient design of air curtains for horizontal refrigerated display cases. International Journal of Thermal Sciences 43:993-1002.

D'Agaro P., Cortella G., and G. Croce. 2006. Two- and three-dimensional CFD applied to vertical display cabinets simulation, International Journal of Refrigeration 29:178-190.

Field B. S. and E. Loth. 2006. Entrainment of refrigerated air curtains down a wall, Experimental Thermal and Fluid Science 30:175-184.

Foster A. M., M. Madge, and J.A. Evans. 2005. The use of CFD to improve the performance of a chilled multi-deck retail display cabinet Internal Journal of Refrigeration 28:698-705.

Navaz, H.K., R. Faramarzi, D. Dabiri, M. Gharib, and D. Modarress. 2002. The application of advanced methods in analyzing the performance of the air curtain in a refrigerated display case. ASME J. of Fluids Engineering 124: 756-64.

Navaz, H.K., B.S. Henderson, R. Faramarzi, A. Pourmovahed, and F. Taugwalder. 2005. Jet entrainment rate in air curtain of open refrigerated display cases. International journal of Refrigeration 28:267-275.

Navaz, H. K., M. Amin, D. Dabiri, and R. Faramarzi. 2005. Past, present, and future research toward air curtain performance optimization. ASHARE Transactions: Symposia: 1083-88.

Navaz H. K., M. Amin, S.C. Rasipuram, and R. Faramarzi. 2006. Jet entrainment minimization in an air curtain of open refrigerated display case. International Journal of numerical methods for heat & fluid flow 16(4):417-30.

Stribling, D., S.A. Tassou, and D. Marriott. 1997. Optimisation of the design of refrigerated display cases using computational fluid dynamics, AIRAH Journal, May 34-44.

Stribling, D., S.A. Tassou, and D. Marriott. 1997. Two-Dimensional CFD Model of a Refrigerated Display Case, ASHARE Transactions:Research, 89-94.

Y.F. Chen, PhD

H.W. Lin, PhD

W.D. Hsieh, PhD

J.Y. Lin, PhD


C.C. Wang


Y.F. Chen and H.W. Lin are researchers, W.D. Hsieh is an engineer, J.Y. Lin is a senior researcher, and C.C. Wang is a senior lead researcher in the Department of Energy and Environment Research Laboratories, Industrial Technology Research Institute, Hsinchu, Taiwan.
COPYRIGHT 2009 American Society of Heating, Refrigerating, and Air-Conditioning Engineers, Inc.
No portion of this article can be reproduced without the express written permission from the copyright holder.
Copyright 2009 Gale, Cengage Learning. All rights reserved.

Article Details
Printer friendly Cite/link Email Feedback
Author:Chen, Y.F.; Lin, H.W.; Hsieh, W.D.; Lin, J.Y.; Wang, C.C.
Publication:ASHRAE Transactions
Article Type:Report
Geographic Code:1USA
Date:Jul 1, 2009
Previous Article:A closer look at [CO.sub.2] as a refrigerant.
Next Article:An experimental evaluation of HVAC-Grade carbon dioxide sensors--part I: test and evaluation procedure.

Terms of use | Privacy policy | Copyright © 2019 Farlex, Inc. | Feedback | For webmasters