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

: This study uses a simpliﬁed two-dimensional time-dependent computational ﬂuid dynamic (CFD) model to study the inﬂuence of wall shelf conﬁguration 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 conﬁguration, 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 unmodiﬁed 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.


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.

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.
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.   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.   Figure 1).

Parameter Unmodified Version Modified Version
Length without end walls, mm 2500 2500 Volume, m 3 1

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: 1.
The thermal properties of the fluid in the air-on section (Figure 3) are in a steady-state.

2.
Mass transfer from load to air (weight loss) is not considered.

3.
The moisture transfer process is not considered.

4.
The evaporator and fans change the air-on stream.
The mathematical model of heat transfer in the fluid interface solves the following equation: 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: Equation of turbulence dissipation rate: 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: where w is the molecular weight; p abs is the absolute pressure; R 0 is the universal gas constant.

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) 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).      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   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.     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. 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). 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.   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 Figures 6 and 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.   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 Figures 6 and 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.

Simulation of Air Curtain Formation Time and Its Shape
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. 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. 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 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.

Simulation of the Temperature Distribution in the Display Cabinet
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. Appl. Sci. 2022, 12, x FOR PEER REVIEW 1 Figure 8. Temperature (°C) distribution in the middle plane of the unmodified cabinet. Figure 9 presents the air temperature distribution in the middle plane of the mo cabinet. A cold air bubble formed in the improved refrigerator, closer to the exterio tour line. This means that the refrigeration cycle is not influenced by environment tors. In this configuration, cold air circulates inside the cabinet with less infiltrat ambient warm air.    Figure 9 presents the air temperature distribution in the middle plane of the mod cabinet. A cold air bubble formed in the improved refrigerator, closer to the exterior tour line. This means that the refrigeration cycle is not influenced by environmenta tors. In this configuration, cold air circulates inside the cabinet with less infiltratio ambient warm air.   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. 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. 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 10. Average outtake air temperature at the middle plane as a function of operating time. 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.  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.

Comparison of the Intake and Outtake Air Temperature and Infiltration Ratio
Appl. Sci. 2022, 12, x FOR PEER REVIEW 13 of 16 Figure 11. Average intake air temperature at the middle plane as a function of operating time.
The infiltration ratio shows how effectively the air curtain works in the refrigerated cabinet. It can be calculated using the following formula: where Tin is the intake air temperature (return air grille); Tout is the outtake air temperature (discharge air grille or honeycomb); Tair 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 mod- The infiltration ratio shows how effectively the air curtain works in the refrigerated cabinet. It can be calculated using the following formula: 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.   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).  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).  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). 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). Figure 13. Distribution of the temperature of food products in the middle plane of the refrigerated cabinet (experimental data). Figure 13. Distribution of the temperature of food products in the middle plane of the refrigerated cabinet (experimental data).

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%.

Conflicts of Interest:
The authors declare no conflict of interest.