2D CFD Simulation of Dynamic Heat Transfer in an Open-Type Refrigerated Display Cabinet

Abstract

This study uses a simplified two-dimensional time-dependent computational fluid dynamic (CFD) model to study the influence of wall shelf configuration and outtake honeycomb inclination angle on the performance of the air curtain in an open-type semi-vertical refrigerated display cabinet. Numerical simulation results show that changing the angle, along with changing the shelf configuration, can reduce the temperature of the intake air grille and improve the warm/cold air ratio. The average temperature of the simulated food products decreased by 1.2 °C compared to the unmodified display cabinet and began to meet the requirements of the temperature category 3M0. The proposed CFD model was validated by the experimental results. The simulation results obtained are important for cabinet design to improve product storage quality and reduce energy consumption.

1. Introduction

Open-type refrigerated display cabinets are important parts of retail food stores used to ensure proper food storage and trading, as they help to attract buyers and increase sales [1]. Food products can be stored and displayed to consumers at different temperature levels in such cabinets [1]. The temperature of food products usually varies between −1 and +7 °C (3M2 temperature category) as described in BS EN ISO 23953-2 [2] for chilled food application. This temperature category is suitable for storing milk and dairy products, sandwiches, green salad, and cooked meat products, but is not appropriate for storing minced meat and raw meat. The storage of raw and minced meat requires a temperature between −1 and +4 °C (3M0 temperature category) as described in BS EN ISO 23953-2 [2].

From a commercial point of view, open-type refrigerated display cabinets have a great advantage because they do not have physical barriers (such as doors) between the buyer and the food, which allows consumers to have free access to the food [1]. However, due to the open design, food products stored inside the display cabinet are isolated from the environment only by one or more cold air curtains. Air penetration due to open design can contribute up to 67–77% of the total heat load of the open-type display cabinet [1,3]. Freezers with doors consume significantly less energy due to the door barrier [3,4].

Replacing open-type refrigerated display cabinets with doored ones becomes an alternative and provides many advantages including improved food quality, achievement of energy savings, and reduced temperature fluctuations. The authors of [5,6,7] experimentally and numerically investigated closed refrigerated display cabinets under various conditions: with/without doors and at various external air temperatures. However, although the doored display cabinets exhibited higher energy efficiency, open display cabinets are always preferred by the buyers, as open-type cabinets attract them and make it easier for them to reach the food products [8]. On the other hand, improvements can often be made to the design of display cabinets to improve air curtain performance. Numerical simulation is an economic way to better understand these systems and predict their thermal behavior more easily [9]. The authors of [10] presented results of 3D CFD simulation of air flows and temperatures in horizontal open-top display units. In this paper, the effect of the perforated plate was studied, and the porous-jump model is used. Rai et al. [11] numerically investigated the protective properties of the air curtain in a refrigerated truck. It has been found that an air curtain of optimal velocity can reduce energy consumption by about 50%. Sun et al. [1] compared the performance of open-type refrigerated display cabinets with and without air-guiding strips. The results of the CFD simulation showed that the air-guiding strips accelerated the air curtain vertically and the cooling capacity necessary to keep the food cooled was decreased by 34%. Tsamos et al. [8] proposed using cold shelf innovation with guiding strips to improve the performance of vertical multi-deck refrigerated display cabinets. The results showed energy savings of 16.7 kWh/24 h compared to the conventional cabinet. Nikitin [12] numerically studied refrigerated display cabinet designs and compared the performance of a single-jet air curtain with the performance of a double-jet one. The simulation results showed that the effective volume of the open-type refrigerated cabinet with a double-jet air curtain can be 13% larger compared to a single-jet air curtain. Many researchers studied the efficiency of open vertical display cabinets by addressing the parameters that have the greatest effect on keeping product temperatures within the correct range [9,13,14,15]. However, most of the CFD modeling studies presented in the literature were developed in 3D configuration and limited to steady-state.

This paper presents a simple, rapid, time-dependent 2D numerical simulation of air flows and temperature changes in an open-type semi-vertical refrigerated display cabinet carried out for different honeycomb inclination angles and shelf configurations. Heat transfer between ambient air, inside air, and food products, as well as flow turbulence, were simulated using commercial CFD software Comsol Multiphysics 5.5 (COMSOL AB, Stockholm, Sweden). The aim of this work is to compare unmodified and proposed refrigerated cabinets in terms of air curtain formation time, the temperature distribution inside the cabinet, the average temperature in the middle plane, intake and outtake air temperatures, and the temperature of the M-packets at the various places inside the cabinet.

2. Materials and Methods

2.1. Research Object

The open-type semi-vertical refrigerated display cabinet is schematically shown in Figure 1. This cabinet is used to store food at temperatures between −1 and +7 °C. The cabinet has 5 shelves for food products, including 4 wall shelves and the bottom panel. The overall dimensions of the cabinet are 1460 mm (height), 970 mm (depth), and 2500 mm (width). The spaces between the shelves are as follows: between the base and the first shelf, 280 mm; between the wall shelves 250 mm. The construction provides a 10 mm clearance between the rear edge of the shelves and the perforated back panel, allowing products to have better cold exchange with the air from the panel and reducing the potential risk of product freezing. Other dimensions and properties of the cabinet are presented in Figure 1 and

Table 1. Characteristics of the open-type display cabinet (Figure 1).

Parameter	Unmodified Version	Modified Version
Length without end walls, mm	2500	2500
Volume, m3	1.8	1.8
Exposition space, m2	6	5.375
Honeycomb angle, °	40	55
Shelf L1, mm	600	600
Shelf L2, mm	500	450
Shelf L3, mm	500	450
Shelf L4, mm	450	400
Shelf L5, mm	350	250
Dimensions of test product package, mm	100 × 100 × 50	100 × 100 × 50
Dimensions of DAG (width × height × depth), mm	2500 × 120 × 20	2500 × 120 × 20
Dimensions of RAG (length × width), mm	2500 × 70	2500 × 70
Number of propeller fans	4	4
Propeller fan rpm and flowrate, m3/h (back flow pressure 12 Pa)	1300; 400	1300; 400
Connection type	Hydroloop	Hydroloop

.

The required temperature inside the display cabinet is maintained by two cooled air streams. The first air stream moves down from a honeycomb (Figure 1) placed at the top of the display cabinet and creates an air curtain. The second cold air stream penetrates from the back through the perforated back panel and moves towards the front where these two air streams are mixed. The air is directed to a return air grille by 4 fans (Figure 1, right), then passes through the cooling coils to decrease the temperature, and returns to the display cabinet.

2.2. Mathematical Model

A two-dimensional (2D) CFD model was developed to study airflows and their influence on temperature distribution inside the unmodified and modified open-type display cabinet using time-dependent studies. The aim of this study was to examine the ability of the 2D k-ε model to reproduce the main phenomena and to predict the product temperatures and the speed of formation of the air curtain in the display cabinet. To simplify the numerical simulation and keep the basic characteristic of the process, the following assumptions are made in the present simulations:

• The thermal properties of the fluid in the air-on section (Figure 3) are in a steady-state.
• Mass transfer from load to air (weight loss) is not considered.
• The moisture transfer process is not considered.
• The evaporator and fans change the air-on stream.

The mathematical model of heat transfer in the fluid interface solves the following equation:

$$\rho C_p\left(\frac{\partial T}{\partial t}+u·\nabla T\right)+\nabla ·\left(q+q_r\right)={\alpha }_{p}T\left(\frac{\partial p}{\partial t}+u·\nabla p\right)+\tau :\nabla u+Q,$$

(1)

where ρ is the air density; C_p is the specific heat capacity at constant pressure; T is the absolute temperature; t is the time; u is the velocity; q is the heat flux by conduction; q_r is the heat flux by radiation; α_p is the thermal expansion coefficient; p is the pressure; τ is the viscous stress tensor; Q contains heat sources other than viscous dissipation.

For the standard k-ε model, the transport equations have the following forms:

Turbulence kinetic energy equation:

$$\frac{\partial }{\partial t}\left(\rho k\right)+\frac{\partial }{\partial x_i}\left(\rho u_{i}k\right)=\frac{\partial }{\partial x_j}\left(\left(\mu +\frac{{\mu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right)+G_k+G_b-\rho \epsilon -Y_M-S_k,$$

(2)

Equation of turbulence dissipation rate:

$$\frac{\partial }{\partial t}\left(\rho \epsilon \right)+\frac{\partial }{\partial x_i}\left(\rho u_i\epsilon \right)=\frac{\partial }{\partial x_j}\left(\left(\mu +\frac{{\mu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right)+C_{\epsilon 1}\frac{\epsilon }{k}\left(G_k+C_{\epsilon 3}{G}_b\right)-C_{\epsilon 2}\rho \frac{{\epsilon }^2}{k}+S_{\epsilon },$$

(3)

where k is the turbulent kinetic energy; ε is the rate of dissipation of turbulent kinetic energy; μ is the dynamic viscosity; μ_t is the turbulent viscosity; σ_k, σ_ε are the turbulent Prandtl numbers for k and ε, respectively; G_k is the generation of turbulence kinetic energy due to the mean velocity gradient; G_b is the generation of turbulence kinetic energy due to buoyancy; Y_m is the contribution of the fluctuating dilatation in compressible turbulence to the overall dissipation rate; S_k, S_ε are the additional source terms; C_ε1, C_ε2, C_ε3 are the constants.

For an ideal gas, the density is calculated using the ideal gas law, and C_p is a function of temperature:

$$p=\frac{w{p}_{abs}}{R_{0}T},$$

(4)

$$\partial h=C_p\partial T,$$

(5)

$$C_p=C_p\left(T\right),$$

(6)

where w is the molecular weight; p_abs is the absolute pressure; R₀ is the universal gas constant.

2.3. Numerical Procedure

Comsol Multiphysics software was used to simulate air flows and heat transfer between air and food products. Four different mesh sizes (M1–M4,

Table 2. Mesh details.

Mesh	Elements	Triangles	Quads	Edges	Vertex
M1	6605	6605	0	966	265
M2	10,726	7908	2818	1059	265
M3	13,973	10,363	3610	1253	265
M4	17,986	13,506	4480	1470	265

) were used in the initial study to ensure numerical accuracy. The total number of mesh elements for each 2D model is presented in Table 2. Figure 2 presents temperatures obtained at the same points with different meshes. Much finer meshes were particularly applied in the regions near the discharge air grille (Figure 1), return air grille, and solid/fluid interfaces. It was found that the M2 mesh with 10,726 elements was satisfactory since the numerical simulation results did not change as the mesh size decreased (Figure 2). This mesh (Figure 3, left) was used in all studies. The simulation was performed using a computer (64-bit operation system) with Intel(R) Core (TM) i7-7700K CPU running at 4.20 GHz and 32 Gb RAM. The computational time was about one hour using M2 mesh (Table 2).

2.4. Boundary Conditions

Figure 3 (right) shows the boundary conditions of the computational model. The values used for the boundary conditions are given in

Table 3. Boundary conditions.

	Temperature, K	Velocity, m/s	Pressure
Air-on	272.15	0.252
Air-off		0.95
Open boundary	298.15		0
Model	283.15

. Air-on and air-off streams are modelled in a closed-loop model. The velocity values are calculated by the fan flow rate and boundary length. The right boundary is defined as the open boundary with both inflow and outflow conditions supported. One of the most important tasks is to set an open boundary location in the horizontal distance to the display case. The distance significantly influences simulation results. Additional studies were conducted to set the location of the open boundary. Figure 4 shows the distance of the open boundary and the temperatures at the measuring points. The measuring points are marked in location 1 (Figure 4, temperature distribution in the middle plane after 10 s). It can be seen that the results are the same in the 2, 3, and 4 studies (Figure 4, locations 2, 3, and 4). Short distance produces worse results of the airflow temperature for the entire simulation time interval. Long distance produces worse results of the airflow temperature at the early air curtain formation time in that the initial model temperature is set to 10 °C. This value of the initial model temperature was chosen by measuring the temperature after the defrosting process of the refrigerated cabinet. The ceiling and floor are considered adiabatic walls. The air temperature is adopted at the one boundary of air-on. The air-on boundary is defined as the velocity inlet, and the air-off boundary is described as the velocity outlet. The products are defined as solid and non-slip boundary conditions. The thin-layer structure was used to model the honeycomb to provide a more accurate flow direction. The wall shelves, the bottom shelf, the back panel, and other parts of the cabinet, were made from steel sheets, through which heat transfer occurs. The wall shelves and the bottom shelf were modelled as thin steel walls with a thickness of 1.4 mm. Other inner parts were modelled as thin steel walls with a thickness of 0.7 mm. The boundary condition of the back panel is defined as a porous jump, with other surfaces being walls. Lastly, there is strong natural convection between the refrigerator and the ambient air affected by the density of air at different temperatures. To simulate this natural convection, gravity was activated in the models and the density of the air was calculated based on the ideal gas formulas.

2.5. Experimental Technique

Figure 5 presents a picture of the experimental refrigerated cabinet with thermocouples. The open-type display cabinet is loaded with test product packages. The M-packets consist of oxyethylmethylcellulose (dimension in width × height × length: 100 × 50 × 200 mm; according to ISO 15502, ISO 5155, ISO 7371, ISO 8187, ISO 8561 standards). The M-packets freezing point is −1 °C and thermal conditions are like lean beef. The refrigerated cabinet was sliced by three measuring planes: z = 0 m (middle plane); z = −1.0 m (left plane); and z = +1 m (right plane). Thermocouples were positioned at the geometric center of the test product packages. Carel Easy controllers with 3 probes were used in this study. The controllers were connected to the computer with Carel PlantVisor PRO monitoring software (CAREL S.p.A., Padova, Italy). The data collected was represented with Microsoft Excel 2016 (Microsoft Corporation, Redmond, WA, USA) software.

The temperature of the test products was measured in the middle plane of the open refrigerated cabinet, which was operated at a room temperature of 25 °C (humidity ratio of 0.003–0.004 kg/kg). Time-averaged temperatures were calculated over a 6 h cycle. The results of temperature measurements were compared with the numerical prediction of temperature distributions in the present study.

3. 2D Numerical Simulation Results and Discussion

3.1. Simulation of Air Curtain Formation Time and Its Shape

Figure 6 represents the complete formation of an air curtain within 60 s in the middle plane of the unmodified cabinet. The velocity of airflow out of the honeycomb is 0.6 m/s. The curtain air velocity is reduced to approximately 0.4 m/s near the return air grille (RAG, Figure 1) at the bottom of the cabinet. Due to the small area of the intake grille, the maximum intake air velocity value of 1.2 m/s is reached in this zone (Figure 6).

Figure 7 represents the complete formation of the air curtain in the middle plane of the modified cabinet with a changed honeycomb angle and length of wall shelves (Figure 1 (right), Table 1). The honeycomb was turned 15 degrees clockwise (Table 1, Figure 1). It can be seen from Figure 7 that unlike the case of an unmodified cabinet (Figure 6), the air curtain is formed in the outside contour of the cabinet, which is represented by a solid black line in Figure 6 and Figure 7. The air curtain is formed after 20 s only, that is, three times faster than in the unmodified cabinet case (Figure 6). The airflow velocity values obtained are similar to those obtained for the unmodified cabinet (Figure 6). A significant difference is observed only in the formation time and position of the air curtain. The ambient air infiltration ratio is less than in the case of an unmodified cabinet due to the shorter air curtain and the cold air travel time.

The warm air in the room has an impact on the air curtain of the refrigerated cabinet. In the case of the unmodified cabinet (Figure 6), the airflow radius is greater than in the modified cabinet (Figure 7). A larger radius means a weaker airflow and makes it easier for the warm air from the surroundings to infiltrate inside the cabinet. Under real conditions in retail stores, customers and hall workers who take out and put in food can also damage the airflow bubble. The air curtain may not even be restored during the opening hours of the retail store due to frequent interruptions. Therefore, the proposed design changes will create a stronger and shorter air curtain with a shorter formation time, which will also allow one to reduce energy consumption.

3.2. Simulation of the Temperature Distribution in the Display Cabinet

Figure 8 presents the air temperature distribution in the middle plane of the unmodified cabinet. A cold air bubble formed in the cabinet far from its front. This means that the refrigeration process is strongly influenced by environmental factors.

Figure 9 presents the air temperature distribution in the middle plane of the modified cabinet. A cold air bubble formed in the improved refrigerator, closer to the exterior contour line. This means that the refrigeration cycle is not influenced by environmental factors. In this configuration, cold air circulates inside the cabinet with less infiltration of ambient warm air.

The shape of the cold air temperature bubble has an impact on the operation of a refrigerated cabinet. In the case of the unmodified cabinet, about 1/3 of the cold air is outside the cabinet contour line (Figure 8). This means that the cooled air mixes with warm ambient air, which adversely affects the operation of the refrigerator. In the case of the modified refrigerator, the cold air bubble only slightly crosses the contour line (Figure 9). This air curtain configuration is more protected from warm ambient air infiltration, and the operating conditions of the cabinet are less dependent on environmental factors.

3.3. Comparison of the Intake and Outtake Air Temperature and Infiltration Ratio

Figure 10 presents the changes in honeycomb outtake air temperature versus time in the middle plane of the unmodified and proposed cabinet. After 20 min of operation, the outtake temperature was −0.9 °C in the case of the unmodified cabinet and −0.73 °C in the modified version. The temperatures are nearly the same, and the difference appears in the airflow resistance in the area.

Figure 11 presents the intake air temperature change graph obtained for the middle planes of the unmodified and modified cabinets in a time-dependent study. The simulation results obtained for the unmodified refrigerated cabinet show that the intake temperature after 20 min of operation is approximately 4.78 °C in the intake grille (RAG, Figure 1). However, in the modified one, the temperature is about 4.07 °C. The higher temperature of the intake air means that the evaporator works under worse conditions. This causes evaporator icing. On the other hand, the intake temperature has an impact on the front products of the bottom shelf, and it is necessary to keep the temperature as low as possible.

The infiltration ratio shows how effectively the air curtain works in the refrigerated cabinet. It can be calculated using the following formula:

$$IR=\frac{T_{in}-T_{out}}{T_{air}-T_{out}}100\%,$$

(7)

where T_in is the intake air temperature (return air grille); T_out is the outtake air temperature (discharge air grille or honeycomb); T_air is the ambient air temperature.

For an unmodified cabinet, the infiltration is 22% and in the case of the proposed design, the infiltration ratio decreases to 18.7%. The difference is 3.3%, therefore, the modified refrigerator works more efficiently using less energy. The main consumer of electricity in this case is the compressor. Lower energy consumption also means shorter operating times and longer service life.

3.4. Comparison of the M-Packets Temperatures in the Middle Plane of the Cabinet

Figure 12 presents a distribution of the M-packets (for a definition of the M-packets, see subchapter 2.5) temperatures in the middle plane of both refrigerated cabinets. Temperature values were obtained in the time-dependent study. The minimum temperature in the unmodified refrigerator reaches a value of 1.4 °C at location B1 (Figure 12). The highest value of 5.4 °C was obtained at location F6. The average food temperature is +3.38 °C. The air curtain cannot separate the volume of the refrigerator from the environment, so the most unfavorable temperatures are obtained at the front of the refrigerator. Additionally, the unmodified refrigerator meets only the requirements of the 3M2 class and can be used to store food products at temperatures ranging from −1 to +7 °C. In the case of the modified cabinet, the lowest temperature value of 0.5 °C was obtained at the same location B1 (Figure 12) and the highest value of 3.8 °C is observed at location F8. The average food temperature is +2.18 °C. The proposed cabinet satisfies the requirements of the 3M0 class (from −1 °C to +4 °C).

4. Experimental Results

Figure 13 presents the distribution of the M-packets temperatures measured in the middle plane of the modified refrigerated cabinet. The average cooling temperature for a 6 h cycle was −1 °C. In all M-packets, the temperature does not reach +4.5 °C (3M0 temperature category, tolerance +0.5 °C). Comparing the results of CFD and experimental results, the maximum difference in temperatures was 9.33% (

Table 4. Comparison of the temperatures of food products obtained in numerical and experimental studies.

Location (Figure 13)	Temperature Obtained in Numerical Study, °C	Average Measured Temperature, °C	Difference, %
F1	0.86	0.91	+5.81
F2	2.57	2.33	−9.33
F5	3.55	3.74	+5.35
F6	3.62	3.78	+4.42

).

5. Conclusions

A simplified time-dependent 2D heat transfer model was developed to analyze airflows in open-type refrigerated cabinets according to the BS EN ISO 23953-2 standard.

The results of the numerical simulation show that the new configuration of the shelves and the changed honeycomb angle decrease the average temperature of the food products from +3.38 to +2.18 °C and the intake air temperature from 4.78 to 4.07 °C. The infiltration ratio decreases from 22 to 18.7%, respectively.

Numerical simulation results show that the modified refrigerated cabinet meets the requirements of the 3M0 class while the unmodified cabinet satisfies only the requirements of the 3M2 class.

The difference between simulation and experimental results is less than 10%.