Abstract
A comprehensive and detailed computational fluid dynamics (CFDs) modelling of air flow and heat transfer in an open refrigerated display cabinet (ORDC) is performed in this study. The physical-mathematical model considers the flow through the internal ducts, across fans and evaporator, and includes the thermal response of food products. The air humidity effect and thermal radiation heat transfer between surfaces are taken into account. Experimental tests were performed to characterize the phenomena near physical extremities and to validate the numerical predictions of air temperature, relative humidity, and velocity. Numerical and experimental results comparison reveals the predictive capabilities of the computational model for the optimized conception and development of this type of equipments. Numerical predictions are used to propose geometrical and functional parametric studies that improve thermal performance of the ORDC and consequently food safety.
1. Introduction
Nowadays, the energy spent during the commercialization of food products stored at refrigerated temperatures is about 50% of the total consumption of energy of a typical supermarket [1]. The growing of energy consumption in the commercial sector is due to the increasing demand of large quantities of perishable food products in urban areas and to an effective rigorous regulation of the sector as well as quality and food safety requirements of consumers. In past decades, the investments and costs associated with food refrigeration increased with the installation of refrigeration equipment to maintain perishable products in perfect conditions for consumption [2]. This trend continues as the worldwide demand for commercial refrigeration equipment that is projected to rise 5.2% annually through 2014. Beverage equipment demand will post the fastest gains among products. Reach-in and walk-in coolers and freezers are expected to post solid gains due to their widespread use in all of the major markets. Display cabinets will benefit from a growing middle class, which will spur gains in the food and beverage retail segment [3]. A large part of the refrigeration equipments installed in supermarkets, and retail stores are vertical open refrigerated display cabinets (ORDCs). This merchandising solution contributes to the largest part of electrical energy consumption related with refrigeration on these sites [4]. This type of appliance is equipped with a recirculating air curtain that establishes an aerothermodynamics barrier between the conservation space and the environment, without physical restrictions for the consumer. Thus, the product to be acquired can be seen and handled without inconvenience.
The air curtain reduces the infiltration of exterior air at higher dry bulb temperature and specific humidity. The effectiveness of this aerothermodynamics barrier changes with thermal and mass-diffusive effects that affect thermal entrainment. These effects, among others, depend on flow instabilities and boundary effects. These conditions lead to a decrease in conservation quality of food products and greater energy consumption and costs. In what concerns this topic, it must be highlighted that the ambient air infiltration load is around 67% to 77% of ORDC cooling load as exposed by Gaspar et al. [5]. The thermal entrainment of hot ambient air in empirical recirculating air curtains is due to inbalance of air distribution between discharge air grille (DAG) and perforated backpanel (PBP), and consequently to air curtain temperature, velocity, and thickness values as described by D’Agaro et al. [6] and Gray et al. [7]. Also, refrigerated air leakage from the bottom part of the frontal opening (cold leg effect) increases the energy loss to surroundings. Those facts, with the remaining heat gain components [1], lead to an increase of thermal load and consequently higher energy consumption. ASHRAE [8] points out that thermal load reduction is the first step for a better energy efficiency of refrigeration equipments that may be accomplished by air curtain optimization, reducing the thermal entrainment with ambient air, as well as the temperature of air returned to the evaporator. Therefore, most of research in this field, both experimental and numerical, to evaluate the thermal performance of ORDC focus on air curtain. Due to similarity, the research on nonrecirculating air curtains is helpful. It is mainly related with heat and mass transfer studies, considering the impinging jet of several angles, different initial velocities, temperatures and thicknesses, and generation/suppression of turbulence inside air stream [9–14]. Although, due to specific characteristics of air curtains installed in ORDC, researchers carried out two- (2D) and three- (3D) dimensional computational fluid dynamics (CFD) parametric studies. Cortella et al. [15] and Navaz et al. [16] evaluated the influence of DAG velocity in thermal performance, quantifying the air infiltration through the frontal opening. Axell and Fahlén [17] developed a CFD-parametric study to evaluate the influence of air curtain height/width ratio and inlet velocity on the thermal performance. Navaz et al. [18] evaluated the optimum operating condition based on entrained air amount, taking into account the jet width and velocity and inlet turbulence intensity. Foster et al. [19] developed 3D CFD models to analyse the effect of changing size and position of the evaporator coil, width and angle of DAG, and inserting baffle plates into the upper duct. D’Agaro et al. [6] carried out 2D and 3D CFD-parametric studies to evaluate the influence of: longitudinal ambient air movement; display cabinet length, and air curtain temperature on the extremity effects and how it reflects on ORDC performance. Chen [20] developed CFD-parametric studies to evaluate the thermal barrier performance of air curtains, adjusting the length/width ratio and its discharge angle, the height/depth ratio of the cavity, and the dimension and position of shelves. Ge and Tassou [21] developed correlations for heat transfer across air curtain with reasonable agreement with experimental data in steady-state conditions. Ge et al. [22] developed an ORDC model integrating CFD and cooling coil models. The airside inputs of the cooling coil model are the outputs of the CFD model, and, inversely, the airside outputs from the cooling coil model are used to update the boundary conditions of the CFD model. The use of the validated integration model of the ORDC allows several analyses such as the optimal designs of the geometrical structures of ORDC, curtains and coils, and further of alternative control strategies and operating states. However, this type of simulation is computationally expensive.
Other research works were just experimental, such as the study developed by Chen and Yuan [23] to evaluate the influence of ambient air temperature and relative humidity, indoor airflow, DAG velocity, PBP airflow, and night-covers application, on the performance of an ORDC. Gray et al. [7] also conducted an experimental study to evaluate the effect of perforation pattern of PBP on the airflow distribution. Among the experimental techniques used by researchers are thermocouple thermometry, hot-wire/film anemometry, laser Doppler anemometry, digital particle image velocimetry, hygrometry, tracer gases, and infrared thermography. Other method that can be used to evaluate the thermal performance of an ORDC, based on the thermal barrier provided by the air curtain, consists in the thermal entrainment factor (TEF) calculation [6, 18, 23, 24]. This formulation was adopted by Gaspar et al. [5] to evaluate air curtain TEF for different environmental conditions (air temperature, relative humidity and velocity, magnitude, and direction) and how it influences the thermal performance of the ORDC.
The experimental techniques described above are reliable and provide results with a high degree of confidence, but its use is labour intensive and time consuming, involving high costs, being the results dependent on ORDC geometry, DAG parameters, and ambient air conditions. By other hand, the development of a complex simulation model requires large computational resources, being the computational results untrustworthy, until reliability and accuracy are checked by comparison with experimental data. Although the significant progress of CFD codes in geometry creation, data transfer from CAD/CAE packages, improvement of meshing tools, numerical robustness, and well-validated physical models for specific classes of physical phenomena, made them desirable in the design phase as an expedite analysis method. However, none of CFD research works cited above has modelled the several devices that make part of the refrigeration system. Hence, this study provides a detailed CFD model that combines the characteristics of aforementioned works. It considers a 2D CFD simulation of an ORDC taking into account the air flow through internal ducts, across fans, evaporator and grilles, and the thermal response of food products. The objective is to provide a complete CFD design tool for this type of equipment, allowing fast and efficient development of parametric studies devoted to predict how operational and geometrical modifications can improve the ORDC performance. This paper presents some results derived from experimental testing to define boundary conditions and to validate numerical predictions, the physical-mathematical and numerical model formulations, the CFD-modelling parameters and methodology for complete simulation of the air flow and heat transfer of an ORDC. The nonisothermal turbulent flow and temperature field characteristics are predicted, and several parametric studies to improve the ORDC performance are designed.
2. Experimental Testing
2.1. Experimental Apparatus
Figure 1(a) shows the midplane cross-section of the open refrigerated display cabinet (ORDC) experimentally tested and numerically modelled on heat and mass transfers. The air is drawn by fans located in front of the evaporator. The air passing through the evaporator is cooled below the conservation temperature of the perishable products exposed in equipment’s shelves. This air is convected to the rear duct, where part of it is discharged inside the conservation space at low velocity through the perforated backpanel (PBP). The other part of this air mass flow rate will supply the air curtain, which develops vertically between discharge (DAG) and return (RAG) air grilles.
(a) Display cabinet configuration and dimensions
(b) Conservation space measuring locations (CSML)
(c) Air curtain measuring locations (ACML)
The authors have already performed extensive experimentation on this ORDC model [5]. It is a self-contained system, in which the condensing unit and controls are built into the cabinet structure (beneath the cabinet, taking up the entire lower part). Its dimensions are 1900 × 796 × 1911 mm (L × W × H). It has four shelves and a well tray, being its frontal opening height, , 1209 mm (see Figure 1(a)). The evaporator is housed under the well tray and integrated with refrigerant feed and return lines. Temperature regulation and defrosting are digitally conducted by an electronic thermostat, being the cooling system fan-assisted (four forward fans that blow air through the coil).
The experimental test was performed in a climate chamber Aralab Fitoclima 650000 EDTU. The test probes described in Table 1 and placed inside the cabinet as shown in Figure 1(b) (conservation space measuring locations (CSMLs)) and Figure 1(c) (air curtain measuring locations (ACMLs)) were connected to a data acquisition system Intab PC-Logger 3100. A probe-positioning system is used to evaluate the 3D effects of thermal entrainment on air curtain and properties variations along length and height of the conservation space [25]. The probe-positioning system was settled on each shelf of the equipment to measure air temperature, relative humidity, and velocity for three positions across the air curtain width and eight vertical cross-sections along the equipment’s length. The positioning system moved the test probes in 240 mm increments for the 1800 mm length of shelves, taking 1 min to move between positions to reduce flow perturbation. Also, the properties values were acquired 1 min after reaching each position to ensure flow stabilization. The experimental results obtained with this point measuring technique along the air curtain and conservation space coordinates show a similar behaviour to the experimental results obtained by Gray et al. [7], Chen and Yuan [23], and Evans et al. [26]. It was found that the average variation of air temperature was 0.4 K, and of air-relative humidity was 4.5% which are not particularly significant. Based in these results, the probes were distributed in the midplane of the equipment’s length. In order to consolidate and extend the measurements data, a thermoanemometer AM 4003 was used to measure air temperature and velocity near DAG and RAG, at the refrigeration unit (compressor + condenser) and at several heights of PBP. Pressure loss was measured using a micromanometer Air Instruments Resources MP3KDS near inlets and outlets (DAG, RAG, and PBP). The temperature of internal surfaces was obtained with K-type thermocouple contact probe using a digital thermometer FLUKE 51 (location n.°9). Additionally, several sets of experimental tests were performed to evaluate the influence of surrounding air in air temperature, velocity, and humidity inside the conservation space [5].
2.2. Testing Procedure and Results
The experimental testing follows the procedure defined by EN-ISO Standard 23953 [27] for test room climate class n.°3, which considers air temperature ( = 298 K; air relative humidity, = 60%; air velocity magnitude, = 0.2 m , with direction parallel to the frontal opening plane of the ORDC, that is, = 0°) and for M-package temperature class M1 (272.15 K ≤ ≤ 278.15 K). The experimental test followed the procedure described in Gaspar et al. [5]. Table 2 contains the average values of the parameters measured during this period and their statistical parameters: maximum and minimum values, standard deviation (), standard error of mean (), and variance (). Notice that overall values of air temperature and relative humidity inside conservation space (probes location n.°0 to n.°4) consider the average values along equipment’s length and height.
3. Mathematical Formulation
The turbulent air flow and nonisothermal heat transfer process are modelled by a 2D (length midplane neglecting extremity effects) steady-state mathematical model. The basic equations governing transport phenomena in an ORDC are continuity, momentum, and energy [28]. For application on the computational domain, the airflow modelling is coupled to products modelling, where within the latter region the convection term is neglected.
The air is considered as an ideal gas. Considering that flow can be driven by buoyancy forces in specific zones of the domain, the Boussinesq approximation is applied. Thus, density has a constant value in all solved equations, except for the buoyancy term in the momentum equation where it is determined as function of temperature differences. So, reference density is set to = 1.17 kg m−3 for the operating conditions defined by test room climate class n.°3 [27].
The energy equation is developed as function of temperature in steady state with constant specific heat. Further simplifications are accomplished neglecting viscous dissipation due to flow characteristics.
The turbulence is modelled by the RNG model [29] since in previous works [30–33] its ability to model turbulence was assessed in comparison with other models (Spalart-Allmaras, , standard ). The comparison of numerical results shows the applicability of models type for simulation of a wide range of flows with minimum coefficients adjustment and also for its relatively simple formulation. The numerical predictions obtained with standard and RNG models were similar. The latter model was used due to its constants revaluation and additional term in dissipation rate equation, to improve the precision of flow simulation.
The set of model equations is suitable for fully turbulent flow. To account for viscous effects and high gradients in proximity of walls, the turbulence model equations are used in conjunction with empirical wall functions. The complete description and implementation details of wall functions in turbulence models can be found in Rodi [34] and Launder and Spalding [35].
The influence of ambient air-relative humidity is considered by making use of a species transport model. The fluid is considered as a mixture of dry bulb air and water vapour. The heat gain of ORDC by thermal radiation is one of most important cooling load components [5]. A surface-to-surface radiation model (based in surfaces view factors calculation) is used to take this heat gain component into account.
4. Numerical Model
The mathematical model is a set of coupled nonlinear partial differential equations, describing mass, momentum, and energy conservation which can be simultaneously and interactively solved. The equations set was solved using the pressure-based, nonstructured grid, finite volume method CFD code FLUENT.
4.1. Geometry and Computational Mesh
The 2D geometry closely follows a mid-cross-section of ORDC. An automatic orthogonal unstructured mesh generator (Gambit) was used to develop the computational grid. The mesh was refined in the internal ducts, near fans, evaporator, DAG, and RAG vicinities and across air curtain, where velocity and temperature gradients are expected to be higher. The grid size was refined to predict accurately heat transfer by conduction inside the conservation products. Grid dependence tests were carried out for models with different grid size, that is, increasing number of control volumes (cells): (a) 17 939 cells, (b) 72 350 cells, (c) 110 029 cells, and (d) 350 533 cells. From the comparison with experimental results, the computational mesh of model (c) provides numerical results independent of grid size both in solid and fluid regions. The computational domain, grid, and boundary conditions location are shown in Figure 2. Notice that connected but independent mesh zones were defined near DAG, RAG, fans, evaporator, and internal ducts in order to speed up the process of developing future parametric studies of geometrical and/or functional modifications and the analysis of their influence on the ORDC overall performance.
4.2. Discretization of the Partial Differential Equations
The computational procedure is based on a numerical iterative process using the pressure-implicit with splitting of operators (PISOs) algorithm [36] for pressure-velocity coupling. This algorithm was derived from semi-implicit method for pressure-linked equations (SIMPLEs) algorithm [37], but it has higher performance, and it is more efficient as described by Jang et al. [38].
The equations were discretized in the control volume form using MUSCL differencing scheme. The monotone upstream-centered schemes for conservation laws (MUSCLs) scheme proposed by Van Leer [39], derives from central differencing scheme (CDS) and second-order upwind (SOU) differencing scheme) as described by Patankar [37]. It is a differencing scheme with low numerical diffusion; that is, it shows higher spatial precision for all types of computational grids and for complex flows. In this study, the comparison with experimental results shows that numerical predictions computed with MUSCL scheme are more realistic and precise. The models run on a server Intel Xeon DualCore 2.33 GHz (4 MBytes internal cache) with 16 GBytes RAM.
4.3. Boundary Conditions
Boundary conditions (BCs) of common practice in numerical simulations, defined for climatic class n.°3 of EN-ISO Standard 23953 [27], are imposed in the computational domain.
4.3.1. Ambient Boundary Condition
The ambient boundary is simulated by an “opening” type BC, that is, a constant pressure boundary which allows both inflow and outflow. The pressure value is considered to be the total pressure based on the normal component of the velocity when flow enters the domain and static pressure when it leaves the domain. The ambient air temperature is supposed to be = 298.15 K, and the water vapour mass fraction for = 60% is = 11.80 gv kgm−1. The radiative black body temperature is assumed to be = 298.15 K for the algebraic calculation of radiative view factors as described by Modest [40]. A black body emissivity value is assumed in this BC.
Free stream values for turbulent kinetic energy and its dissipation rate are assumed on the free boundary at this fixed-pressure condition set in terms of turbulence intensity, = 10%, (the worst situation considers the effects of consumers passage in front of equipment, of air-conditioning system operation and the influence of pressure perturbations) and hydraulic diameter of equipment’s frontal opening to ambient, = 1.2 m, as being the characteristic turbulence length scale. The values of parameters specified at fixed pressure BC are shown in Table 3.
The influence of the distance at which this free boundary is defined is not consensual as it affects the predicted solutions. The numerical studies of ORDC developed by Cortella et al. [15], and George and Buttsworth [41] consider this BC at and , respectively, of the equipment width. Bhattacharjee and Loth [9] stated that this BC should be imposed far away from zone of interest in the computational domain, likewise at a distance of 40, being the air curtain width. Gaspar et al. [42] tested different distances of this free BC to ascertain the effect on the numerical predictions of temperature and velocity within an ORDC. The following cases were tested: (1) W; (2) W; (3) W; and (4) W, being the width of the ORDC. Based on comparison between numerical predictions and experimental values, the precision of numerical results increases with the distance at which is defined the BC. Thus, the fixed pressure BC was defined at a distance equal to W from the frontal opening of the equipment since it allows modelling the flow with heat and mass transfers more precisely.
4.3.2. Wall Boundary Condition
(1) Shear Boundary Condition
Wall boundary conditions are used to bound fluid and solid regions. At the walls a nonslip BC (zero velocity) is considered.
(2) Thermal Boundary Conditions
Heat Flux Boundary Conditions
An adiabatic BC is defined for walls not considered in heat transfer calculation. However, a heat flux BC is used to simulate heat generated by the illumination (85% for fluorescent lamp OSRAM L58W/20) and heat flux through conduction across material layers that compose the equipment walls. The heat flux across them is determined by Fourier Law using a global heat transfer coefficient determined by the conductive thermal resistances of each wall material. The experimental values of surface temperature on the internal and external sides of the equipment are used. Table 4 shows the heat flux values fixed as BC.Species Boundary Conditions
A zero-gradient condition for water vapour mass fraction is assumed at walls.Radiation Boundary Conditions
It is necessary to specify the emissivities of different surfaces to use a surface-to-surface radiation model. Constant emissivities are fixed for the internal surfaces of the equipment = 0.9 and for external ground, = 0.7. It is considered the black body emissivity, = 1, for the external enclosure surfaces.
4.3.3. Product Load (Solid Region)
Following EN-ISO Standard 23953 [27] and ISO 15502 [43], the products simulators are made of tylose, in which thermal characteristics are similar to meat. Considering values given by ASHRAE [44], the equivalent solid thermal characteristics to simulate food products are shown in Table 5.
4.3.4. Air Pressure Drop through the PBP: Porous Medium
The detailed flow simulation across perforation of the backpanel requires a very high grid refinement and consequently a huge computational effort. Thus, PBP was modelled as a porous medium to simplify the numerical model. The thickness of PBP is 1.6 mm. The porous medium is composed by 520 grid cells (4 grid cells along thickness and 130 grid cells along height): = 0.4 mm and = 20 mm. The porous medium is modelled by the addition of a momentum source term corresponding to Forchheimer law to the standard fluid flow equations [45]. The source term is composed of two components: a viscous loss term (Darcy’s law, the first term on the right-hand side of (1)), and an inertial loss term (the second term on the right-hand side of (1)) due to high flow velocities [45] and low thickness to hole diameter ratio [46]. We have
In (1), is the average value of fluid velocity through the surface normal to flow, is the medium permeability, and is the inertial resistance factor. The permeability of perforated surfaces is defined by Bear [45] and Tang et al. [47], in which the thickness of the porous medium, δ, has an analogical correspondence to the perforation diameter, as described in (2). We have where is the porosity (void fraction), that is, the volume fraction of fluid within the porous region (i.e., the open volume fraction of the medium) given by the ratio between the open volume and total volume of porous medium (see (3)). One finds that
The constant value parameter provides a correction for inertial losses in porous medium. This constant can be viewed as a loss coefficient per unit length along flow direction, thereby allowing the pressure drop to be specified as a function of dynamic head (see (4)). As this factor is specified for fully open porous cells, the loss coefficient, K, must be converted into dynamic head loss, considering the same flow rate per unit length of porous region [48, 49]. So we have where and are respectively the fluid velocities of the fully open porous medium and of flow passage through the porous medium. The loss coefficient is calculated as proposed by Idel’Cik [50]. The values imposed at PBP are shown in Table 6.
4.3.5. Air Pressure Drop through DAG and RAG: Porous Medium Model Simplification
The air pressure drop is experimentally measured at DAG and RAG. Equation (5) represents a one-dimensional simplification of the porous media model defined at the faces that simulate these grilles. The thin porous medium has a finite thickness, δ, over which the pressure change is defined as a combination of Darcy’s Law and Forchheimer inertial loss term. We obtain that
The fixed parameters on these BC based in geometrical characteristics of grilles and experimental measures of air temperature and relative humidity are shown in Table 7.
4.3.6. Heat and Mass Transfer and Pressure Drop Modelling across the Evaporator
The ORDC contains a wavy fin and tube heat exchanger (evaporator). At the evaporator, it is specified a heat exchanger BC type in which are defined both pressure drop and heat transfer coefficient as functions of velocity in direction normal to the heat exchanger. The BC is specified at an infinitely thin face, and the pressure drop through the heat exchanger is assumed to be proportional to the dynamic head of the fluid. The friction factor, = 0.0078, is calculated through the empirical correlation (6) proposed by Wang et al. [51]. We have
The pressure drop, = 8.33 Pa, is calculated by (7) proposed by Kays and London [52], being contraction, , and expansion, , loss coefficients given by McQuiston and Parker [53]. We have
The heat transfer modelling through a heat exchanger BC type in the CFD code requires the specification of overall heat transfer coefficient, U, and refrigerant (R404A) temperature, , at the main entrance of evaporator. The variation of refrigerant temperature across evaporator (with phase change process) is reduced [54], so its value can be considered constant. It is assumed no superheating of refrigerant, and all heat transfer is on the two-phase zone at constant temperature. The value determined for the overall heat transfer coefficient, U = 124.16 W m−2 K−1, is referred to the air temperature downstream of heat exchanger, . Based in the equation of energy conservation and considering the experimental values of air properties upstream and downstream the heat exchanger, the value of refrigerant temperature at heat exchanger entrance, = 271.05 K, is determined by a trial and error-iterative procedure. These parameter values are based in the heat balance of the evaporator, determined with the average experimental values.
The heat flux from heat exchanger to surrounding fluid, = 124.31 W m−2, is determined by (8). Taking into account the total surface area of the evaporator, , this value is consistent with the ORDC’s cooling load determined by Gaspar et al. [5]. We have
However, since the evaporator is modelled as an infinitely thin BC, it is necessary to obtain the equivalent heat flux based on surface areas ratio given by (9). We obtain that
The determined equivalent overall heat transfer coefficient, = 910.42 W m−2 K−1, is assumed to be similar to a convective heat transfer coefficient. This simplification was based on the evaluation of empirical correlations for Colburn factor and friction factor available for heat exchangers with corrugated fins.
The mass transfer modelling through a heat exchanger BC type consists in the specification of a Dirichlet BC for water vapour mass fraction, considering the air dehumidification process. By psychometric analysis, the water vapour mass fraction specified at evaporator exit is = 3.3722 gv .
4.3.7. Fan Boundary Condition Modelling the Pressure Rise across the Fan
A discontinuous pressure rise across fans (infinitely thin face) is specified as a function of air velocity. The empirical characteristic curve which governs the relationship between head (pressure rise) and flow rate (velocity) across a fan element (obtained at the manufacturer: EBM Papst series 4500) is converted into a 4th-order polynomial relationship (10) and specified as BC. We obtain that
4.4. Solution Monitoring and Control Techniques
The linear relaxation method is used to reduce high variation of dependent variables during the iterative process of calculation. Linear relaxation values ranged from 0.3 for pressure to 0.8 for momentum. The convergence monitoring was done by the sum analysis of absolute residuals of mean field variables. The iterative procedure run until a prescribed convergence criterion for absolute residuals is met.
5. Results and Discussion
The CFD modelling of air flow and thermal patterns inside the conservation space of a vertical ORDC has an operation temperature varying from 273.15 K to 278.15 K. This type of equipment design still requires studies due to the major influence of the frontal ambient opening on the heat and mass transfers.
The numerical simulations allow the evaluation of air temperature, relative humidity, and velocity distributions within the equipment in order to propose possible paths for technical evolution. Only most significant results concerning flow properties and thermal behaviour are discussed here. The numerical predictions solution required approximately 4 h 30 m performing 10 000 iteration sweeps.
5.1. Comparison with Experimental Data
The validation of numerical predictions of the 2D ORDC model is accomplished by its comparison with experimental measurements data.
Figures 3, 4, and 5 show the comparative profiles, for different planes, of the experimental average values of air temperature, relative humidity, and velocity. The predicted steady-state air flow and heat transfer inside the ORDC both present a reasonable quantitative agreement. The experimental average points are represented by ○, and numerical predictions points are represented by ■. For the validation points set (90 validation points), the minimum absolute deviations from experimental measures for air temperature, relative humidity, and velocity were respectively: (y/W = 0.34, z/H = 0.39) = 0.08 K (conservation zone at 3rd shelf height), (y/W = 0.52, z/H = 0.88) = 0.26% (conservation at 4th shelf height), and (y/W = 0.52, z/H = 0.88) = 0.002 m s−1 (near RAG). The highest quantitative discrepancies are found outside the air curtain limits on the external ambient side (y/W 1), thus without relevant significance for the study.
Figure 6 shows the comparison of numerical predictions of product core average temperatures with experimental data obtained by Gill et al. [55] (ten retail stores were tested), Foster et al. [19] (difference between minimum and maximum package temperature is around 8 K), Evans et al. [26] (single average temperature value for all M-packages through its standard deviation), Gray et al. [7] (ORDC with inclined shelves), and Lu et al. [56] (evaporator in the back of the cabinet and possessing two air curtains). The predictions of core products temperature are compared with results from other studies because a uniform initial value for all products core temperature was not ensured. The products core temperatures are in the range of the results obtained by the abovementioned studies.
Considering the range of air temperature, relative humidity and velocity variations measured during experimental testing, the deviation between experimental data and numerical predictions is acceptable for this type of engineering application. The nonoverlapping between experimental and numerical results is due to experimental errors (measurements precision, physical phenomena perturbation, etc.), physical-mathematical model assumptions (2D, steady state, turbulence model, PBP and evaporator formulations, open boundary BC definition, etc.), numerical model (the purpose of the MUSCL scheme is to reduce the higher-order terms to first-order eliminating instabilities. However, its use for incompressible flows can damp the physical solution in some extend, providing nonphysical predictions especially in the vicinity of mixing layers and where instabilities occur), i.e., in the air curtain (a phenomenon also observed by D’Agaro et al. [6] and Hammond et al. [57]). Nevertheless, the combined analysis of experimental data and numerical predictions shows that the computational model follows the physical phenomena occurring in the real equipment.
5.2. Numerical Results Analysis and Discussion
The numerical predictions of air flow and heat transfer for the ORDC model subjected to climatic class n.°3 of EN-ISO Standard 23953 [27] are analysed in this section. This analysis will allow the development of numerical parametric studies to improve the global performance of this type of equipments.
5.2.1. Velocity Field Predictions
The numerical predictions of air velocity field are shown in Figure 7. It can be observed the main characteristics of air flow, low air velocity between shelves in contrast with high velocities in ducts and at DAG exit. The thermal mixture between refrigerated air curtain and ambient is detected by the increase of air curtain thickness as it moves toward RAG.
A substantial spillage of refrigerated air to ambient can be observed in RAG vicinity due to increased mass flow rate as consequence of entrainment from air curtain. It is also related with the air curtain momentum weakening, being the air velocities quite smaller when compared to those near DAG. These effects convert into energy losses.
Flow details are analysed locally. Figure 8 shows the velocity vectors predictions near DAG and RAG. Several eddies that promote thermal mixture are predicted in DAG vicinity as shown in Figure 8(a). A parametric study where fan velocity and DAG angle influence on thermal performance are studied can provide valuable information to an improved design of the ORDC.
| (a) DAG region, (m s−1) |
| (b) RAG region, (m s−1) |
Some recirculations around food products are predicted, as well as the “plug flow” from PBP carrying heat towards air curtain. Other parametric study that can provide an insight about this condition is related to the PBP holes diameter and their distribution, homogeneous or not, in the panel.
The velocity vectors field near RAG is shown in Figure 8(b). Spillage to exterior ambient and eddies formation are predicted at this location. These conditions can be further analysed in a parametric study where different values of fan velocity and RAG angle are tested. Higher air velocities inside rear duct, smaller ones through PBP, and the lowest air velocities inside the conservation space are predicted. Also, vortices are predicted on duct curves, which consequently decrease the air velocity at DAG. Thus, a parametric study can be developed to analyse the inclusion of air guides or deflectors inside ducts to improve air flow by reducing friction. A recirculation region is predicted due to interaction between low velocity air in the conservation zone and high velocity of air curtain. This recirculation reduces air curtain momentum, decreasing its performance as aerothermodynamics barrier to ambient air.
5.2.2. Temperature Field Predictions
The homogeneity degree of temperature distribution in the refrigerated display space is shown in Figure 9. A temperature gradient rising along ORDC’s width and height is predicted, suggesting that thermal entrainment occurs. The radiation heat transfer affects mostly the surfaces “viewed” from the external environment. The maximum value of air temperature is predicted in the well tray. The internal temperature of products is not constant, as products located at shelves front are more exposed to thermal entrainment through air curtain and by thermal radiation. These results are in concordance with numerical predictions obtained by Cortella [58] and with experimental results obtained by Gill et al. [55], Foster et al. [19], and Evans et al. [26].
(a) computational domain
(b) products conservation zone (reduced scale)
5.2.3. Relative Humidity Field Predictions
The numerical results of air relative humidity are coherent with psychometric analysis, that is, the air-relative humidity decreases with air temperature increase. At the interface region between air curtain and ambient, the relative humidity increases due to thermal interaction promoted by eddies development.
5.2.4. Mass Flow and Heat Transfer Rates across the Air Curtain Predictions
A region that virtually encloses air curtain boundaries is considered as shown in Figure 10. These boundaries enclose five control volumes (VC). The numerical results of mass flow (11) and heat transfer (12) rates per unit length are determined for each one. One has
Figures 11 and 12 show the values and directions of mass flow and heat transfer rates across air curtain per unit length. These results indicate that air curtain gains thermal energy from ambient and conservation region. The latter gains are related to thermal radiation, illumination, and heat conduction through equipment walls.
(a) West faces (w)
(b) East faces (e)
(a) West faces (w)
(b) East faces (e)
The mass flow and heat transfer rates predictions show that recirculated air curtain gains energy from external and internal borders, W m−1. It gains energy from ambient at higher heights (East faces of VC n.°1, 2, 3, and 4) and losses energy near RAG (East face of VC n.°5) at the external border. Simultaneously, it loses energy to the conservation space (West faces of VC n.°1, 2, and 3) and gains it at lower heights (West faces of VC n.°4 and 5). These predictions show that the air curtain is a dynamic system, gaining and losing energy, depending on air velocity, turbulence, air mass ratio between discharge grille and perforated backpanel, among other factors. An air curtain that provides full protection is an unfeasible condition in ORDC. However, the analysis of CFD predictions, such as those provided by the CFD model described along this paper, can help to improve the aerothermodynamics blockage efficacy.
6. Conclusions
A CFD model of an ORDC, including ducts, grilles, perforated backpanel, evaporator, and fans, has been developed to simulate air fluid flow and heat transfer. The aim of the CFD model is to cover the simulation of both air curtain and air flow inside the ducts, which is not common way of using CFD in this type of application. The main objective consists in the development of a detailed CFD model which allows an expedite simulation of design improvements aimed to increase the thermal performance and to reduce the energy consumption of ORDC.
The characterization of air flow and heat transfer allows identifying parameters that can be adjusted and hence reducing the impact of thermal entrainment and improving global performance of ORDC. The main thermal load in this type of equipments is the infiltration of ambient air. The determination of parameters that influence the overall air curtain efficacy, likewise mass flow and heat transfer rates across it, can provide valuable information to performance improvement. The progressive downward thermal entrainment into the air curtain is very dependent on the development of eddies which trigger the mixing. The momentum reduction decreases air curtain stability. This condition promotes a nonuniform air temperature distribution inside the ORDC and influences temperature differences inside the products. Also, the predicted spillage of cold air to surroundings near RAG and the temperature value at this location affect the energy performance of the equipment.
The agreement between numerical predictions and experimental results is quite good and adequate to these type of engineering problems although the authors recognize that this model can be improved to perform transient regime simulations, and the possibility of considering three-dimensional effects will provide additional achievements. Nevertheless, this CFD model is suitable to investigate the effect of a number of modifications that affect the equipment’s overall performance, that is, aiming the reduction of thermal entrainment rate and maintaining stable air curtain momentum until it reaches RAG.
Based on the numerical predictions analysis, several parametric studies can be developed to evaluate the influence of some parameters in thermal entrainment and temperature distribution homogeneity inside products conservation region, as these factors are directly related with food safety and energy consumption. As determined by the numerical predictions analysis, these parameters involve fans velocity; hole diameters and their distribution on the backpanel perforation; DAG and RAG angles; guides and deflectors inside ducts, among others. In the future, parametric studies of the abovementioned parameters based on this detailed CFD model will be developed with the aim to analyse design modifications that improve the overall performance of the open display cabinet.
Nomenclature
| : | Area, (m2) |
| : | Total surface area, (m2) |
| : | External tube surface area, (m2) |
| : | Air curtain width, (m) |
| : | Inertial resistance coefficient, (m−1) |
| : | Specific heat, (J kg−1 K−1) |
| : | Fin collar outside diameter, (m) |
| : | Hydraulic diameter, (m) |
| : | Energy, (J) |
| : | Force, (N) |
| : | Friction coefficient |
| : | Fin pitch, (m) |
| : | Fin spacing, (m) |
| : | Gravitational acceleration, (m s−2) |
| : | Mass flux of the air based on minimum flow area, (kg s−1 m−1) |
| : | Height, (m) |
| : | Turbulence intensity, (%) |
| : | Loss coefficient |
| : | Turbulent kinetic energy, (m2 s−2); permeability, (m2); thermal conductivity (W m−1 K−1) |
| : | Viscous resistance coefficient, (m−2) |
| : | Length, (m) |
| : | Mass flow rate, (kg s−1) |
| : | Number of longitudinal tube row |
| : | Pressure, (Pa) |
| : | Waffle height, (m) |
| : | Longitudinal tube pitch, (m) |
| : | Transverse tube pitch, (m) |
| : | Heat flux, (W m−2) |
| : | Heat transfer rate, (W) |
| : | Radius, (m) |
| Re: | Reynolds number |
| : | General source term |
| : | Standard error of mean |
| : | Temperature, (K) |
| : | Global heat transfer coefficient, (W m−2 K−1) |
| : | Average velocity, (m s−1) |
| : | Width, (m) |
| : | Spatial coordinates (along length, width, and height), (m) |
| : | Projected fin pattern length for one-half wavy length, (m) |
| : | Water vapour mass fraction, (kgv kgm−1). |
| 0: | Total surface |
| : | Air |
| amb: | Ambient |
| bb: | Black body |
| : | Air curtain; contraction; cross-sectional; fin collar |
| comb: | Honey comb |
| cons: | Conservation |
| : | Expansion |
| eq: | Equivalent |
| evap: | Evaporator |
| e, w, n, s: | Control volume faces identification (East, West, North, and South) |
| : | Fin |
| : | Hydraulic |
| : | Component of cartesian directions according to , , and |
| in: | Input; upstream |
| m: | Mixture |
| out: | Output; downstream |
| ref: | Reference |
| refrig: | Refrigerant fluid |
| sup: | Surface |
| : | Total; turbulent |
| : | Water vapour; void. |
| : | Density, (kg m−3) |
| : | Dissipation rate of turbulent kinetic energy, (m2 s−3); emissivity; porosity |
| : | Relative humidity, (%); general dependent variable |
| : | Dynamic viscosity, (kg m−1 s−1) |
| : | Standard deviation; ratio of the minimum flow area to frontal area |
| : | Variance |
| Absolute humidity, (kgv kga−1) | |
| Thickness, (m) | |
| Convergence criterion. |
| ACML: | Air curtain measuring location |
| BC: | Boundary condition |
| CAD: | Computer-aided design |
| CAE: | Computer-aided engineering |
| CSML: | Conservation space measuring location |
| DAG: | Discharge air grille |
| ORDC: | Open refrigerated display cabinet |
| PBP: | Perforated backpanel |
| RAG: | Return air grille |
| VC: | Volume of Control |
| TEF: | Thermal entrainment factor. |
Acknowledgments
The authors wish to acknowledge the support of University of Beira Interior-Engineering Faculty-Electromechanical Engineering Department, Covilhã, Portugal (http://www.ubi.pt/) and of JORDÃO Cooling Systems, Guimarães, Portugal (http://www.jordao.com/). The authors do not have any conflict of interest is with the content of the manuscript.