Airflow simulation in a horizontal closed refrigerated display cabinet

Vítor Hugo Moura Guedes Airflow simulation in a horizontal closed refrigerated display cabinet Master’s Dissertation in Mechanical Engineer Work developed under supervision of Professor Doutor Pedro Alexandre Moreira Lobarinhas Professora Doutora Senhorinha de Fátima Capela Fortunas Teixeira July 2019 ii DIREITOS DE AUTOR E CONDIÇÕES DE UTILIZAÇÃO DO TRABALHO POR TERCEIROS Este é um trabalho académico que pode ser utilizado por terceiros desde que respeitadas as regras e boas práticas internacionalmente aceites, no que concerne aos direitos de autor e direitos conexos. Assim, o presente trabalho pode ser utilizado nos termos previstos na licença abaixo indicada. Caso o utilizador necessite de permissão para poder fazer um uso do trabalho em condições não previstas no licenciamento indicado, deverá contactar o autor, através do RepositóriUM da Universidade do Minho. Licença concedida aos utilizadores deste trabalho Atribuição CC BY https://creativecommons.org/licenses/by/4.0/ iii ACKNOWLEDGMENTS To my tutors, Professor Pedro Lobarinhas, for all the knowledge and support provided, to Professor Senhorinha Teixeira, for the help and make me think a step ahead. To Master João Vasconcelos, for following this project since its beginning and all the help especially with the software ANSYS – FLUENT. To my family: my mother, my father and my brother, for making this project possible and all the support during it. To my friends and colleagues that help me in many ways. iv DECLARAÇÃO DE INTEGRIDADE Declaro ter atuado com integridade na elaboração do presente trabalho académico e confirmo que não recorri à prática de plágio nem a qualquer forma de utilização indevida ou falsificação de informações ou resultados em nenhuma das etapas conducente à sua elaboração. Mais declaro que conheço e que respeitei o Código de Conduta Ética da Universidade do Minho. v RESUMO Os equipamentos de refrigeração tornaram-se essenciais em várias vertentes, especialmente na conservação de alimentos. Nos últimos anos os equipamentos de refrigeração fechados têm ganho destaque devido ao seu consumo ser muito inferior comparativamente com os abertos. O estudo centra-se na análise do escoamento de ar frio que percorre o equipamento, tendo em conta as condições exteriores. Para este efeito recorreu-se ao uso do programa CFD ANSYS Fluent. Este usa a técnica de volumes finitos na solução das equações de conservação da massa, momentum e energia. Na primeira parte do trabalho fez-se um conjunto de simulações 2D em que foram estudados vários parâmetros. Primeiramente, foi analisada a dimensão da malha mais adequada (com apenas elementos quadriláteros); de seguida, foi feita uma comparação entre diferentes tipos de saídas e modelos de turbulência. Para terminar o estudo 2D foram colocados produtos de diferentes dimensões, atuando como obstáculos ao escoamento. Quanto à influência dos parâmetros, pode concluir-se que o mais predominante é a temperatura na zona inferior junto ao evaporador. O fluxo de calor junto à parede e a influência das lâmpadas na geometria com produtos afetam pontualmente, não se refletindo noutras zonas geometria de uma forma evidente. O aumento da área perfurada nas costas na zona inferior permite uma maior ventilação nessa zona, contudo a temperatura aumenta ligeiramente. Quanto à geometria com produtos, pode concluir-se através dos resultados que a geometria 2D utilizada é bastante limitadora. Palavras-Chave: Vitrine fechada, simulação CFD, escoamento de ar, ANSYS Fluent. vi ABSTRACT Refrigeration equipment has become essential in several areas, especially in food preservation. In the last few years, closed refrigeration equipment becomes more popular because its consumption is much lower compared to the open ones. The work of this study focuses on the analysis of the cold air flow that traverses the equipment, knowing the external conditions. For this purpose, CFD ANSYS Fluent program was used. This uses the finite volumes technique to solve the equations of conservation of mass, momentum and energy. In the first part of the work, it was made a set of 2D simulations in which several parameters were studied. Firstly, the most appropriate mesh dimension was analysed (using only quadrilateral elements), then a comparison was made between different types of exits and turbulence models. To finish the 2D study products were introduced in the geometry with different dimensions, acting as obstacles to the fluid flow. Concerning the influence of the parameters, it can be concluded that the most predominant is the temperature in the lower zone near the evaporator. The heat flux near the wall and the influence of the lamps on the geometry with products affect punctually, not being reflected in other geometry zones significantly. The increased perforated area in the lower back allows greater ventilation in the lower back, however the temperature increases slightly. In the geometry with products, it can be concluded that the 2D geometry used is very limiting. Keywords: Closed display cabinet, CFD simulation, airflow, ANSYS Fluent vii CONTENTS Acknowledgments ................................................................................................... iii Declaração de integridade ....................................................................................... iv Resumo .................................................................................................................... v Abstract .................................................................................................................. vi Contents ................................................................................................................ vii List of figures ........................................................................................................... ix List of tables ........................................................................................................... xii Nomenclature ........................................................................................................ xiii List of symbols ........................................................................................................ xiv 1. Introduction ...................................................................................................... 1 1.1. Closed showcases .............................................................................................. 2 1.2. Objectives .......................................................................................................... 3 1.3. Dissertation structure ........................................................................................ 4 2. State of art ........................................................................................................ 6 2.1. Problem description ........................................................................................... 7 2.2. Cold equipment classification ............................................................................ 8 2.3. Literature review ................................................................................................ 9 3. Computational model ...................................................................................... 15 3.1. Mathematical model ........................................................................................ 17 3.2. Turbulence models .......................................................................................... 24 3.3. Finite volumes .................................................................................................. 27 3.4. SIMPLE.............................................................................................................. 29 viii 3.5. Residual values and convergence .................................................................... 33 4. Case study ....................................................................................................... 36 4.1. Strategies adopted ........................................................................................... 36 4.2. Mesh definition ................................................................................................ 38 4.3. Mesh optimization ........................................................................................... 40 5. Results and discussion ..................................................................................... 46 5.1. Study of some parameters influence ............................................................... 46 5.1.1. Boundary condition on the bottom boundary: 0°C ..................................... 46 5.1.2. Heat flux 50 and 150W/m2K......................................................................... 49 5.1.3. Back entrances double size bottom ............................................................. 50 5.2. Geometry with products inclusion and lamps influence ................................. 52 6. Conclusion and future work ............................................................................. 56 Conclusions ...................................................................................................... 56 Future work ...................................................................................................... 57 7. References ...................................................................................................... 58 ix LIST OF FIGURES Figure 1.1 - example of food preservation equipment .................................................... 1 Figure 1.2 - cold showcase description ............................................................................ 3 Figure 2.1 - showcase's full scheme, including the cold production zone and the storage zone .................................................................................................................................. 6 Figure 2.2 – Real steam compression system .................................................................. 7 Figure 2.3 - Temperature field using multiple air curtains ............................................. 10 Figure 2.4 - mean temperature and standard deviations of the display cabinet with doors (a) and without doors (b)................................................................................................ 11 Figure 2.5 - Iced drink cabinet description, .................................................................... 12 Figure 2.6 - results at the end of each stage: (a) flow field, (b) temperature field ........ 13 Figure 2.7 - Pathlines for plate evaporator, 1 - with a finned surface, 2 – without a finned surface ............................................................................................................................ 14 Figure 3.1 - Continuity on control volume, .................................................................... 18 Figure 3.2 - Conservation of momentum, ..................................................................... 19 Figure 3.3 - Conservation of momentum in x ................................................................ 20 Figure 3.4 - Energy conservation, ................................................................................... 23 Figure 3.5 – SST k-𝝎 model definition ............................................................................ 26 Figure 3.6 -Fluxes in a 2D element according to finite volumes method ...................... 28 Figure 3.7 - Nodes and vector designation in a 2D element according to finite volumes method ........................................................................................................................... 29 Figure 3.8 - Control volume ............................................................................................ 31 Figure 3.9 - Residual values definition on ANSYS - FLUENT ........................................... 35 Figure 4.1 - Showcase's geometry with no products in the shelves .............................. 36 2 also a greater effort on the part of the manufacturers, who must guarantee the correct functioning of the equipment with the minimum power dissipated. For that reason, closed display cabinets become more and more popular over the years[1]. The most popular way to produce cold in these devices is using refrigerants, that will transport the heat from the cold source to the hot source. Depending on which source is the one is interested (to refrigerate is the cold source and to heating is the hot source), there are two different classifications: heat pumps (heating) and refrigeration machines (cooling). That the operation is the same for both, and there are many types of equipment that can do both tasks, such as air conditioners. 1.1. Closed showcases Commercial refrigeration equipment is indispensable for the preservation of food in such establishments, namely butchers, bakeries, etc. A schematic drawing of the cabinet is presented in figure 1.2. There is an upwards airflow on the back that will be split in two: a fraction of this flow is entering on the equipment by the back perforated plated and the rest of it is entering on the top, near the glass on the front, forming an air curtain. This type of equipment is divided into two zones: cold production and storage (figure 1.2). The cold production zone is usually located at the bottom of the equipment. This zone contains insulating material to prevent energy loss to the outside and there is only one connection to the storage zone to enable the airflow. 3 Figure 1.2 - cold showcase description Air enters the storage area through the back of the showcase (through a perforated plate) and through an existing channel on the top. Afterwards, the air is exhausted from the bottom of the storage zone, entering the cold production zone again, starting a new cycle. The existence of a perforated plate is crucial to ensure uniformity of flow in the storage zone. This area is bounded by a glass on the face that is customer-facing because the glass is a transparent material to facilitate customer choice. If there is a freezing zone, it is not visible, and it is attached to the cold production zone. It is also common that in the storage zone there are different temperatures intended for different foods, even in equipment without freezing. 1.2. Objectives As referred before, the main goal is to optimize the showcase in order to reduce its consumptions without compromising its purpose. There are 2 methods to evaluate this kind of situations: computational and experimental. The computational method is an extension of the mathematical, but its faster and it is easier to visualise the problem. The main reason for using CFD as a tool to optimize the equipment is the reduced cost Top air entrance Perforated plate glass storage zone outlet Inlet 3 1 – evaporator 2 – conpressor 3 – condenser 4 – expansion valve 4 2 1 4 and time when compared to experimental studies. So, the project should follow this schedule: • Modelling a closed display cabinet • Study of the flow in a closed display cabinet in order to determine the following parameters: - The temperature range inside it, so that this range of values is the smallest possible and that these temperatures allow the preservation of food. - The flow velocity so that it does not exceed a certain value (so that it does not produce excessive noise and does not dehydrate the food) - Study the influence of different parameters and the relative importance of each one. - Influence of the heat flux from all sources, namely heat transfer from the exterior through the boundaries and dissipated power from electric devices (for example the lamps). - Influence of the products in the refrigerating process. 1.3. Dissertation structure Considering the objectives mentioned before, in the next chapter the problem will be contextualized, explaining the origin of the cold air entering the storage zone. Some classifications of this type of equipment will be presented according to different entities, and finally a brief literature review including CFD and experimental work. Chapter 3 will introduce numerical simulation and the software ANSYS-Fluent used for this project, considering some calculation strategies included on it and its functionalities. In chapter 4, the case study will be presented, including the geometry, mesh definition and mesh optimization. It will be explained how to optimize the mesh to the study. In chapter 5, the influence of several parameters will be studied. The last case studied will include products, for which a new geometry and mesh will be needed. For each model, there will be a small explanation of why this model was made and the next step of the study. 5 Finally, the last chapter will draw some conclusions and propose some future work that can continue this research. 6 2. STATE OF ART In this chapter, the problem will be described, and the equipment will be presented and its working process. Then it will be presented the legislation to this kind of equipment and its classification considering the properties and operation parameters. Lastly, some scientific papers and other works similar to this one will be reviewed. Figure 2.1 presents the full showcase, with the cold production zone and the storage zone. Figure 2.1 - showcase's full scheme, including the cold production zone and the storage zone The next subchapter explains how the cold is produced in this kind of equipment. However, in the present study, this process will be ignored and it will be considered only the storage zone. 4 3 2 1 1 – evaporator 2 – compressor 3 – condenser 4 – expansion valve storage zone 7 2.1. Problem description The most common system used in the showcases to produce cold is by vapour compression (figure 2.2) Figure 2.2 – Real steam compression system[2] At point 1, the refrigerant is in a state of superheated vapour (in an ideal cycle, it would be saturated steam, but in the real case it is not to guarantee that there is no liquid in the compressor). Then there is a compression in order to raise the temperature and pressure of the refrigerant before the condenser, which will remove heat from the fluid. At this stage, the energy from the refrigerant can/should be used to warm the environment or another process in order to monetize all the energy (most important step of the heat pump). The refrigerant at the end of the condenser is liquid (similarly the liquid is subcooled to guarantee there is no gas on the expansion). Then there is an expansion valve before the refrigerant enters in the evaporator when it will be heated by the refrigerated space, restarting another cycle. This refrigerated space will be the study object of this project. 8 There are many other systems used in certain situations that will not be analysed because they are not relevant in this case[3]. Therefore, the main goal of this project is the improvement of the efficiency in a showcase, that is reducing the consumptions without compromising its purpose. The consumption reduction is achieved by lowering the mass flow rate and/or increasing the air inlet temperature (reducing the power of the refrigeration system). In order to optimize the equipment, the temperature field will be evaluated inside the showcase, maintaining the inlet conditions. 2.2. Cold equipment classification Thus, several classifications were assigned according to their application and operating parameters. According to the European Community PRODCOM (from the French 'PRODuction COMmunautaire'), the cold store can be fitted in the code 'CPA 29.23.13: Refrigerating and freezing equipment (table 2.1) and heat pumps, except household type equipment', which can be divided into several subgroups [4]: Table 2.1 - Classification of refrigerating equipment according to PRODCOM Code Description Characteristics 29.23.13.33 Refrigerated show-cases and counters incorporating a refrigerating unit or evaporator for frozen food storage Refrigerated and frozen food 29.23.13.35 Refrigerated show-cases and counters incorporating a refrigerating unit or evaporator (excluding for frozen food storage) Refrigerated good 29.23.13.40 Deep-freezing refrigerating furniture (excluding chest freezers of a capacity ≤ 800 litres, upright freezers of a capacity ≤ 900 litres) Frozen food, small equipment 9 29.23.13.50 Refrigerating furniture (excluding for deep-freezing, show-cases and counters incorporating a refrigerating unit or evaporator) Frozen food, large equipment 29.23.13.90 Other refrigerating and freezing equipment The equipment under analysis is a commercial showcase that does not contain frozen food storage, so, considering this classification, this furniture fits in group 29.23.13.35. There are many other entities using other classifications, such as EUROVENT, US Department of Energy, Energy Star Program requirements, etc. that are more precise, considering the following specifications: • Volume • Operating temperature: - Refrigeration (temperature higher than 0°C) - Freezing (temperature bellow 0°C) - Combined refrigeration and freezing • Vertical or horizontal • Self-contained or remote compression • With or without doors, material (solid doors or glass) • Opening type: sliding or pivoting. The normative EN ISO 23953 classifies equipment according to their orientation (vertical, semi-vertical or horizontal), shape (line-up, islands, etc), operating temperature, service type (self-service or assisted service). 2.3. Literature review Ribeiro et al (2016) [5] developed a numerical study (2D) using 3 different models of an open vertical display cabinet: the first one using a simple air curtain, the second one using deflection blades and the third one using multiples air curtains. These curtains 10 are refrigerated by the evaporator on the bottom. The air velocity was defined at the exit of the evaporator. This study surprisingly obtained better results on the first model (figure 2.3); however, it refers that the second one has more potential if the deflection blades were optimized. According to the author, the third model was not successful because of the increment of complexity. These results have natural limitations due to being performed in 2D simulations. Figure 2.3 - Temperature field using multiple air curtains[5] Chaomuang et al. (2019) developed an experimental study of a closed display cabinet, where the main goals were analysing the factors that influence the temperature inside the equipment and its variation, comparing these results with the same cabinet but with the doors opened (figure 2.4) [6]. The authors observed that 70% of the heat transfer was due to thermal bridges and gaps (for example the air infiltration between the doors). Comparing the temperatures with and without doors (when the cabinet was empty), in the storage zone there was an increment from 0.4 to 3.1°C and from 8.2 to 9.5°C in the shelf edges (near the air curtain). The influence of ambient temperature was 11 also studied: the authors verified that there is an almost linear relationship between the air temperature inside the equipment and the external temperature. Figure 2.4 - mean temperature and standard deviations of the display cabinet with doors (a) and without doors (b), [6] Wang et al. (2015) developed a numerical 3D study of an iced drink refrigeration cabinet, where the temperatures inside fluctuate between -7°C and -5°C, using an automated switch on and automated switch off [7]. The cabinet has an enclosing perforated plate that separates the evaporator and the storage room, which is simulated by the porous jump model. In the coolant air chamber, there is a fan in the bottom to assist the air circulation (figure 2.5). 18 Figure 3.1 - Continuity on control volume, [9] Adding the transient term and dividing all terms by δxδyδz: 𝜕𝜌 𝜕𝑡+𝜕𝜌𝑢 𝜕𝑥 +𝜕𝜌𝑣 𝜕𝑦 +𝜕𝜌𝑤 𝜕𝑧 =0 3.3 The last equation represents the mass conservation of three-dimensional transient flow in a compressible fluid. The first term represents the term transient while the others represent the term convective Momentum conservation law (momentum equation) According to Newton's 2nd law F = ma, the momentum in the control volume particle is equal to the sum of the forces in the particle. The variation of momentum in the x, y, and z axes per unit volume in the control volume is given by: 𝜌𝐷𝑢 𝐷𝑡 𝑖𝑛 𝑥 ; 𝜌𝐷𝑣 𝐷𝑡 𝑖𝑛 𝑦; 𝜌𝐷𝑤 𝐷𝑡 𝑖𝑛 𝑧 3.4 19 There are two types of forces in the control volume: • Surface forces - Forces due to pressure - Forces due to viscosity • Body forces - Forces due to gravity - Centrifugal forces - Coriolis forces - Electromagnetic forces Surface forces will be considered in terms of momentum and the body forces contained in the source term. In a three-dimensional case, the stress state in the fluid element is defined in terms of pressure and nine components of viscosity stresses as shown in Figure 3.2. The pressure is identified by the letter p and the viscous stress by the letter τ. Figure 3.2 - Conservation of momentum, [9] Considering that the forces on the components in x resulting from the pressure p and the surface stresses τxx, τxy and τxz, (Figure 3.3). The resulting force in x is the sum of all the components acting in that direction of the element. 20 Figure 3.3 - Conservation of momentum in x, [9] The resulting forces in x on the faces (E, W) are: [(𝑝−𝜕𝑝 𝜕𝑥1 2𝛿𝑥)−(𝜏𝑥𝑥 −𝜕𝜏𝑥𝑥 𝜕𝑥 1 2𝛿𝑥)]𝛿𝑦𝛿𝑧 + [−(𝑝−𝜕𝑝 𝜕𝑥1 2𝛿𝑥)−(𝜏𝑥𝑥 −𝜕𝜏𝑥𝑥 𝜕𝑥 1 2𝛿𝑥)]𝛿𝑦𝛿𝑧 =(−𝜕𝑝 𝜕𝑥+𝜕𝜏𝑥𝑥 𝜕𝑥 )𝛿𝑥𝛿𝑦𝛿𝑧 3.5 In faces (N, S) in the x-axis results: −(𝜏𝑦𝑥 −𝜕𝜏𝑦𝑥 𝜕𝑦 1 2𝛿𝑦)𝛿𝑥𝛿𝑧+(𝜏𝑦𝑥 −𝜕𝜏𝑦𝑥 𝜕𝑦 1 2𝛿𝑦)𝛿𝑥𝛿𝑧=𝜕𝜏𝑦𝑥 𝜕𝑦 𝛿𝑥𝛿𝑦𝛿𝑧 3.6 Lastly, the forces on direction x on faces (T, B) are given by: −(𝜏𝑧𝑥 −𝜕𝜏𝑧𝑥 𝜕𝑧 1 2𝛿𝑧)𝛿𝑥𝛿𝑦+(𝜏𝑧𝑥−𝜕𝜏𝑦𝑥 𝜕𝑧 1 2𝛿𝑧)𝛿𝑥𝛿𝑦=𝜕𝜏𝑧𝑥 𝜕𝑧 𝛿𝑥𝛿𝑦𝛿𝑧 3.7 21 The total force by volume unit of the fluid resulting from the surface tensions is obtained by the sum of equations mentioned above, which divided by δxδyδz results: 𝜕(−𝑝+𝜏𝑥𝑥) 𝜕𝑥 +𝜕𝜏𝑦𝑥 𝜕𝑦 +𝜕𝜏𝑧𝑥 𝜕𝑧 3.8 The equation of momentum in component x is obtained by the equality of the momentum variation with the resultant of the force due to the surface tensions of the equation (previous equation) plus the term source in x. 𝜌𝐷𝑢 𝐷𝑡 =𝜕(−𝑝+𝜏𝑥𝑥) 𝜕𝑥 +𝜕𝜏𝑦𝑥 𝜕𝑦 +𝜕𝜏𝑧𝑥 𝜕𝑧 +𝑆𝑀𝑥 3.9 Similarly, on the y component: 𝜌𝐷𝑣 𝐷𝑡=𝜕(−𝑝+𝜏𝑦𝑦) 𝜕𝑦 +𝜕𝜏𝑦𝑥 𝜕𝑥 +𝜕𝜏𝑧𝑦 𝜕𝑧 +𝑆𝑀𝑦 3.10 And z component: 𝜌𝐷𝑤 𝐷𝑡 =𝜕(−𝑝+𝜏𝑧𝑧) 𝜕𝑧 +𝜕𝜏𝑦𝑧 𝜕𝑦 +𝜕𝜏𝑧𝑥 𝜕𝑥 +𝑆𝑀𝑧 3.11 Energy conservation law (Energy equation) This equation is derived from the first law of thermodynamics which states that the internal rate of change of a fluid particle is the sum of the variation of heat and the variation of the work. Saying that the variation rate 𝜌 is given by: 𝜌𝐷𝐸 𝐷𝑡 3.12 22 Work is equal to the product of force by the component of velocity in the direction of the force. Making the product of equation 3.9 with the velocity in the component in x: 𝜕𝑢(−𝑝+𝜏𝑥𝑥) 𝜕𝑥 +𝜕(𝑢𝜏𝑦𝑥) 𝜕𝑦 +𝜕(𝑢𝜏𝑧𝑥) 𝜕𝑧 3.13 The same happens with y and z components, obtaining respectively: 𝜕𝑣(−𝑝+𝜏𝑦𝑦) 𝜕𝑦 +𝜕(𝑣𝜏𝑥𝑦) 𝜕𝑥 +𝜕(𝑣𝜏𝑧𝑦) 𝜕𝑧 3.14 𝜕𝑤(−𝑝+𝜏𝑧𝑧) 𝜕𝑧 +𝜕(𝑤𝜏𝑦𝑧) 𝜕𝑦 +𝜕(𝑤𝜏𝑥𝑧) 𝜕𝑥 3.15 Isolating the terms that contain pressure p: −𝜕(𝑢𝑝) 𝜕𝑥 −𝜕(𝑣𝑝) 𝜕𝑦 +𝜕(𝑤𝑝) 𝜕𝑧 =−𝑑𝑖𝑣(𝑝𝑢) 3.16 The total work variation resulting from the superficial forces is given by: [−𝑑𝑖𝑣(𝑝𝑢)]+ [ 𝜕(𝑢𝜏𝑥𝑥) 𝜕𝑥 +𝜕(𝑢𝜏𝑦𝑥) 𝜕𝑦 +𝜕(𝑢𝜏𝑧𝑥) 𝜕𝑧 ] + [ 𝜕(𝑣𝜏𝑥𝑦) 𝜕𝑥 +𝜕(𝑣𝜏𝑦𝑦) 𝜕𝑦 +𝜕(𝑣𝜏𝑧𝑦) 𝜕𝑧 ] + [ 𝜕(𝑤𝜏𝑥𝑧) 𝜕𝑥 +𝜕(𝑤𝜏𝑦𝑧) 𝜕𝑦 +𝜕(𝑤𝜏𝑧𝑧) 𝜕𝑧 ] 3.17 Equation 3.17 represents the balance of all x, y, and z components of work on the fluid particle resulting from forces acting on surfaces. 23 The balance of the heat flux through the faces of the control volume is represented by the vector of heat flow q in the three components qx, qy and qz, according to Figure 3.4. Figure 3.4 - Energy conservation,[9] [(𝑞𝑥−𝜕𝑞𝑥 𝜕𝑥 1 2𝛿𝑥)−(𝑞𝑥+𝜕𝑞𝑥 𝜕𝑥 1 2𝛿𝑥)]𝛿𝑦𝛿𝑧=−𝜕𝑞𝑥 𝜕𝑥 𝛿𝑥𝛿𝑦𝛿𝑧 3.18 The same happens for the components y and z. −𝜕𝑞𝑦 𝜕𝑦 𝛿𝑥𝛿𝑦𝛿𝑧 3.19 −𝜕𝑞𝑧 𝜕𝑧 𝛿𝑥𝛿𝑦𝛿𝑧 3.20 Dividing the equations above by δxδyδz: −𝜕𝑞𝑥 𝜕𝑥 −𝜕𝑞𝑦 𝜕𝑦 −𝜕𝑞𝑧 𝜕𝑧 =−𝑑𝑖𝑣 𝑞 3.21 24 According to the Fourier law for the heat flow by conduction: 𝑞𝑥−𝐾𝜕𝑇 𝜕𝑥; 𝑞𝑦−𝐾𝜕𝑇 𝜕𝑦;𝑞𝑧−𝐾𝜕𝑇 𝜕𝑧 3.22 That can be expressed by: 𝑞=−𝑘𝑔𝑟𝑎𝑑 𝑇 3.23 Combining the two previous equations, the final energy balance equation is obtained resulting from the conduction heat transfer process. −𝑑𝑖𝑣 𝑞=𝑑𝑖𝑣 (𝑘 𝑔𝑟𝑎𝑑 𝑇) 3.24 3.2. Turbulence models Most flows studies are turbulent studies. A flow is defined as laminar, intermediate or turbulent through the Reynolds number. This number is given by the ratio of inertial forces to viscous forces (UL / ν). Since U is the velocity of the fluid, L is the characteristic length of the flow, and ν is the kinematic viscosity. For low values of the Reynolds number the flow is laminar, above the critical value (transient Reynolds value) the flow becomes turbulent. The ANSYS fluent software uses an approach to the numerical simulation names RANS (Reynolds Averaged Navier-Strokes Simulation) [10]. This technique solves timeaveraged Navier-Strokes equations and offers several methods of turbulence for all turbulent length. This is the most common approach in the industry due to his low cost compared to other ones (DNS, LES). There are many turbulent models that can be used in this approach, which must be chosen considering the Reynolds number, fluid properties and its application (Table 3.1). 25 Model Behaviour and usage Spalart-Allmaras Economical for large meshes. Performs poorly for 3D flows, free shear flows, flows with strong separation. Suitable for mindly complex (quasi-2D) external/internal flows and boundary layer flows under pressure gradient (e.g. airfoils, wing, airplane fuselages, missiles, ship hulls Standard k𝜺 Robust. Widely used despite the known limitations of the model. Performs poorly for complex flows involving severe pressure gradient, separation, strong streamline curvature. Suitable for initial screening or alternative designs, and parametric studies Realizable k𝜺 Suitable for complex shear flows involving rapid strain, moderate swirl, vortices, and locally transitional flows (e.g. boundary layer separation, massive separation and vortex shedding behind bluff bodies, stall in wide-angle diffusers, room ventilation). RNG k𝜺 Offers largely the same benefits and has similar applications as Realizable. Possibly harder to converge than realizable. Standard k𝝎 Superior performance for wall-bounded boundary layer, free shear, and low Reynolds number flows. Suitable for complex boundary layer flows under adverse pressure gradient and separation (external aerodynamics and turbomachinery). Can be used for transitional flows (thought tends to predict early transition). Separation is typically predicted to be excessive early SST k𝝎 Offers similar benefits as standard k𝜔. Dependency on wall distance makes this less suitable for free shear flows. Table 3.1 - Turbulence models description 26 There will be studied two different models, Realizable k-𝜀 and SST k-𝜔, which can be the most suitable for this study and its variations. Turbulence model Realizable k-𝜺 (RKE) This model, like the Standard k-𝜀 (SKE) model, is suitable for most situations but has some benefits compared to the previous one, such as: • Accurately predicts the spreading rate of both planar and round jets. • Also, likely to provide superior performance for flows involving rotation, boundary layers under strong adverse pressure gradients, separation, and recirculation. It is expected to have some rotation in this case study, that is the main reason the SKE was excluded for this situation. Turbulence model SST k-𝝎 k-𝜔 models have better performance for boundary layer flows and low Reynolds number flows. SST k-𝜔 specifically, is a mix of a k-𝜔 model near the wall and k-𝜀 in the freestream (Figure 3.6). Figure 3.5 – SST k-𝝎 model definition[10] RSM Physically the most sound RANS model. Avoids isotropic eddy viscosity assumption. More CPU time and memory required. Tougher to converge due to close coupling of equations. Suitable for complex 3D flows with strong streamline curvature, strong swirl/rotation (e.g. curved duct, rotating flow passages, swirl combustors with very large inlet swirl, cyclones). 27 This model is a good compromise between k-𝜔 and k-𝜀 models. Near wall treatment Four different models can be chosen in order to define the flow behaviour near the wall: Standard wall functions, Non-equilibrium wall functions, Enhanced wall treatment, User-Defined Wall functions. As it was done in the turbulence models, standard wall functions and enhanced wall treatment models will be considered. Standard wall functions model is designed for high Re attached flow, the nearwall region is not resolved, the near-wall mesh is relatively coarse. The enhanced wall treatment model is used for low Re flows or flows with nearwall complex phenomena. Usually requires a very-fine near-wall mesh capable of resolving this region but can also handle coarse near-wall mesh. 3.3. Finite volumes The finite volume method uses as a starting point the integral form of the conservation equation. The domain of the solution is divided into a number of control volumes, the above equation is applied to each of them (Figure 3.7). At the central point of each volume control is located a computational node, in which the values of the desired variables are calculated, as it happens in the boundaries of the volume control. The values of the variables at the boundary are obtained through interpolations as a function of nodal values (in the centre). Volume and surface integrals are approximated using appropriate quadrature formulas. As a consequence, an algebraic equation is obtained for each volume control, in which the value of the variable to be studied appears in this node and in the neighbours. This method can be used in any type of mesh, so it can be used in very complex geometries. The mesh itself defines all the boundaries of the volume control; however, it does not need to be related to the coordinate system. This method is 34 𝑎𝑝∅𝑝=∑ 𝑎𝑛𝑏∅𝑛𝑏 +𝑏 𝑛𝑏 3.43 Where 𝑎𝑝 is the coefficient of the centre of the control volume, 𝑎𝑛𝑏 the contribution of the neighbouring cell coefficients and b the contribution of the constant to the source term. 𝑎𝑝=∑ 𝑎𝑛𝑏 −𝑆𝑝 𝑛𝑏 3.44 The calculated value of the residual by software according to the Pressure-based model is given by the unbalance of equation (3.42) summed in all cells, this is a nonstaggered residue and is written as follows: 𝑅∅=∑ |∑ 𝑎𝑛𝑏∅𝑛𝑏+ 𝑏− 𝑎𝑝∅𝑝 𝑛𝑏 | 𝑐𝑒𝑙𝑙𝑠 𝑃 3.45 The software uses two models of staggered residues, representative of the flow ∅ through the domain. These are global scaling and local scaling. They are respectively defined by the following equations: 𝑅∅=∑|∑𝑎𝑛𝑏∅𝑛𝑏 + 𝑏− 𝑎𝑝∅𝑝𝑛𝑏 | 𝑐𝑒𝑙𝑙𝑠 𝑃 ∑| 𝑎𝑝∅𝑝| 𝑐𝑒𝑙𝑙𝑠 𝑃 3.46 𝑅∅=√∑(1 𝑛 𝑛 𝑐𝑒𝑙𝑙𝑠 )(∑𝑎𝑛𝑏∅𝑛𝑏 + 𝑏− 𝑎𝑝∅𝑝𝑛𝑏 𝑎𝑝)2 (∅𝑚𝑎𝑥 −∅𝑚𝑖𝑛)𝑑𝑜𝑚𝑎𝑖𝑛 3.47 Residual values are a good indicator of the convergence of problems. By default, ANSYS Fluent activates the global scaling option. There are some models to evaluate the 35 convergence of results. The definition of residual values in some cases is sufficient to ensure convergence but may be misleading in other cases. It is, therefore, good practice to assess convergence not only by the value of residues but also by monitoring some specifications, such as the heat transfer coefficient. Figure 3.10 shows the configuration panel of the residuals, and in the image, the values of the residuals were assigned by default by the program. Figure 3.9 - Residual values definition on ANSYS - FLUENT 36 4. CASE STUDY In this chapter, the main strategies adopted in the study are presented. The physical parameters, geometry and mesh are presented as well as the optimization of the mesh 4.1. Strategies adopted Analysing this type of equipment, it can be realized that the initial geometry is very simple. The first models are tests that do not include some conditions of the last models. As the simulations go on, these conditions will be considered one by one until the final model is accomplished. This strategy allows us to detect errors and their origins and correct them before the final model is accomplished. The last simulations include products which will transpose in more complex geometry. The first decision to make is choosing between a 2D or a 3D geometry. Showcases usually have a transversal section that is uniform, so the 3D geometry would probably not accomplish better results compared with the 2D, but the complexity and the computational resources would be much higher. For these reasons, this study would be focused on a 2D geometry (figure 4.1) Figure 4.1 - Showcase's geometry with no products in the shelves Zoom A 37 The first concern is to approximate the geometry to reality and simplify it whenever is possible without compromise the reliability of the results. With that said, a velocity inlet will be considered, with air at 0°C and ignore the cold production zone, focusing only on the storage zone. However, some care is needed when transposing the real geometry to 2D because there are some details that cannot be perfectly represented in a 2D geometry. In this showcase is not possible to design the air entrances as they are, because the lateral section is not constant in all the equipment. In the showcase the air intakes are holes but, in this geometry, they will be represented by a line, which means a rectangular rip at the full length. In order to not compromise the results and simplify the geometry, this line represents a rip that has the same area as the sum of the holes. The perforation area ratio (perforated area/total area) in the back is around 7.5%. The first analysis after the geometry was defined was the mesh refinement. The simulation was started with a coarse mesh that would necessarily result in a bad solution. In the next simulations, the mesh was gradually thinner. At some point the solution does not depend on the mesh that was used, that means it is not necessary to refine more the mesh. In the next studies, a mesh that has roughly the same element size will be used. Boundary conditions delimitate different zones of the domain. This limitation can be physical and easy to define, but sometimes the physical phenomenon is hard to translate in a boundary condition. Despite that, it is critical to have a boundary definition that has the best approach to reality in order to minimize the differences between the model and the reality. Most of the boundary conditions were defined according to other experimental and numerical studies mentioned previously on the state of art in chapter 1.4. In this case study different boundary conditions were used and also different initial values, which were optimized during the study: 38 • Inlet velocity = 2 m/s, according to other studies, the velocity is usually bellow 3m/s in order to reduce the noise. This value varies according to the refrigerated volume and inlet area[5], [11], [12]. • Pressure inlet = 0 Pa, the pressure on the equipment should be equal to the atmospheric pressure, this value is used by default on other studies similar to this one. • Outlet: outflow, other studies may use pressure outlet but based on previous studies made for this showcase, outflow obtains better results[13]. • Inlet temperature = 0°C, the objective is to maintain the temperature at 5°C or lower in the product zone, but also higher than 0°C to avoid freezing. The inlet temperature is the lowest on the showcase because there are many heat sources during until the air reaches the outlet (environment temperature, products temperature, lamps, etc) and there are no temperatures below 0°C. • Glass heat transfer: Q(glass) = 8.29 W/m2, this value was used on previous studies for this showcase, but other values will be evaluated on further simulations[13]. • Initial temperature = 15oC, this temperature should be equal to the temperature of the room where the showcase is, usually equal or lower than 25°C [12], [14]. • Turbulence model: RKE, according to the turbulence models available on the software (some of them analysed on chapter 3.2), this should be one of the most suitable models for room ventilation and standard situations like this one[10]. 4.2. Mesh definition Mesh consists in dividing the domain into small areas/volumes (in 2D or 3D simulations). In the 2D model were used only quadrangular elements. In order to accomplish a good mesh to further, the mesh definition was made in two steps. The first step was to choose the type of elements used (the choice was to use quadrilateral elements whenever is possible) and the size of each one related to the others. Smaller elements should be used in areas where stronger variations were expected. There are also some factors given by the program that were considered, like element quality, 39 minimum orthogonal quality, skewness, maximum aspect ratio and y+. These factors consider the element quality based on its shape (aspect ratio, orthogonal quality and skewness). The size variation between neighbour elements should also be verified. The element size variation should be smooth because of calculation stability. The simulation starts only after the mesh definition. The second step is to get the intended precision, this can be accomplished by using a coarse mesh (that may present low-quality results due to the low number of elements) and to increase the number of elements gradually. This raise of the number of the elements must be proportional in order to maintain a similar mesh between smaller elements, that way generates better results. At some point, as the number of elements were increased, the meshes start to converge and should present the same results regardless of the number of elements. The best mesh in this group is the one that has the least number of elements and presents the same results of all these with a higher number of elements. The showcase’s geometry is presented in figure 4.2. As it is shown in figure 4.2 the elements that are closer to the wall are smaller because it is expected to have a higher variation of velocity in this area (boundary layer). This is because the air velocity close to the wall is zero (due to the friction) and quickly increases as moving away from the wall. Figure 4.2 - Mesh definition near the glass and wall, element size = 1.65mm (mesh 4, approximately 230,000 elements) glass wall Zoom A 40 4.3. Mesh optimization As mentioned in the last chapter, the first step is to optimize the mesh quality based on the factors given by the program. Some of these factors are presented in table 4.1. Table 4.1 - Mesh quality report Researching for a better mesh, some refinements were done based on the first mesh (mesh with 61,000 elements). The second and third meshes (with 93,000 and 117,000 elements, respectively) were generated increasing the number of elements in all the domain, because the velocity results reveal some unexpected vortices (figure 4.4, meshes 1, 2) and was expected that the temperature variations were smoother (figure 4.3, meshes 1 and 2). The fourth mesh (with 230,000) was created increasing the number of elements in all zones of the showcase because the temperature field still changed, despite the bigger changes were below the bottom shelf (figure 4.3, mesh 3). The velocity field remains similar, except bellow the top shelf the vortices generated are different (figure 4.4, mesh 3). The fifth mesh (with 415,000 elements) was created increasing the number of elements in all the domain to confirm the results obtained for the fourth mesh (figures 4.3 and 4.4, meshes 4 and 5). Mesh number Mesh report Quality Number of elements Average orthogonal quality Average skewness Average aspect ratio 1 0.92 0.18 1.85 61225 2 0.92 0.18 1.98 93198 3 0.92 0.18 2.21 117247 4 0.93 0.17 2.18 232141 5 0.93 0.17 2.42 415295 41 To the mesh study, there were used many different refinements, resulting in consecutively more elements: roughly 61,000 (mesh 1), 93,000 (mesh 2), 117,000 (mesh 3), 230,000 elements (mesh 4), 415,000 elements (mesh 5). The lines in black represent the main streams of air. The results in figure 4.3 and 4.4 represent 60 seconds in real-time for the temperature and velocity fields, respectively. 42 1 2 3 4 5 Figure 4.3 - Temperature field after 60s in the meshes 1 to 5, T0 = 15oC, v0 = 2 m/s, T(inlet) = 0oC, Q(glass) = 8.29 W/m2K Mesh 1 – 61.000 elements Mesh 2 – 93.000 elements Mesh 3 – 117.000 elements Mesh 4 – 230.000 elements Mesh 5 – 415.000 elements Temperature 43 1 2 3 4 5 Figure 4.4 - Velocity field after 60s in the meshes 1 to 5, T0 = 15oC, v0 = 2m/s, T(inlet) = 0oC, Q(glass) = 8.29W/m2K Mesh 1 – 61.000 elements Mesh 2 – 93.000 elements Mesh 3 – 117.000 elements Mesh 4 – 230.000 elements Mesh 5 – 415.000 elements Velocity 50 Figure 5.6 - temperature field after 120s, T0 = 15oC, v0 = 2 m/s, T(inlet)=0oC, Q(glass)=8.29 W/m2K, T(bottom) = 0oC As it can be seen on figures 5.5 and 5.6, the impact of this heat flux is only local for both 50 and 150 W/m2K. The maximum temperature in the highest heat flux is not admissible. Although, as it can be seen in figure 5.6, those values only appear near the glass, which will not affect eventual products. Although the air curtain is at a higher temperature than the rest of the air inside the showcase, the hot air will be exhausted in a short period of time, preserving the conditions on the zone products. On the next simulations, the heat flux will be 50 W/m2K. 5.1.3. Back entrances double size bottom In previous simulations was evident that there was an issue below the bottom shelf because the velocity was too low, and that could cause some problems to the quality of the products, depending on which products are refrigerated. To avoid this, it was decided to increase the size of the holes in this zone in order to increase the mass flow (figures 5.7 and 5.8). Figure 5.5 - temperature field after 120s, T0=15oC, v0 = 2 m/s, T(inlet) = 0oC, Q(glass) = 8.29 W/m2K, T(bottom) = 0oC 51 On this study, it was decided to use only mesh 4 because the previous simulations have the same results for mesh 4 and 5, so there is no point simulate with both meshes. As it can be seen in figure 3.14, the velocity was increased on the pretended zone, although there is still a small area where the velocity remains low. However, this also caused a reduction of the mass flow on the rest of the showcase (because the inlet Figure 5.8 - temperature field after 120s, T0 = 15oC, v0 = 2 m/s, T(inlet) = 0oC, Q(glass) = 8.29 W/m2K, T(bottom) = 0oC, with double sized entrances on the bottom of the perforated plate Figure 5.7 - velocity field after 120s, T0 = 15oC, v0 = 2 m/s, T(inlet) = 0oC, Q(glass) = 8.29 W/m2K, T(bottom) = 0oC, with double sized entrances on the bottom of the perforated plate 52 velocity/mass flow was maintained), which caused a slight increase of the temperatures overall (Figure 5.8). This is not critical in this particular case, but this model is not yet completed because the products and some other parameters that would affect negatively the performance were not considered yet. Then, it can be concluded that this variation is not the best in a general case, but it can be a better solution if the products stored need ventilation. 5.2. Geometry with products inclusion and lamps influence In this final stage, the products in the geometry will be introduced. These products have different dimensions and the bigger products are at the bottom because there is more space to place them (figure 5.9). Figure 5.9 - geometry with products On the bottom, the dimensions are 200 x 150 and 250 x 150 mm, respectively. On the shelves, the products are 70 to 80 mm wide, and 100 to 210 mm height, except the smaller one on the top shelf that is 50 mm height. The products have similar dimensions that some refrigerated food and drinks (packed products, water bottles, 53 soda cans, etc) and the gaps between have 20 mm or less. The point is to evaluate different situations where the space between the products are refrigerated or not. The heat flux in the products will not be studied so there is no mesh in them, they act only as obstacles to the flow. 5.2.1. Mesh definition The mesh definition is different from the previous ones because the introduction of products makes it impossible to produce the same mesh. The first change was the type of mesh: before it was a quadrilateral mesh and on this one is hybrid (with quadrilateral and triangular elements). In this geometry, the same technique shown in Figure 5.2 could not be used, so there are triangular elements near the glass (Figure 5.10). On the other hand, the elements on the rest of the geometry are quadrilateral and the principals followed on the mesh construction are the same, such as the element size being smaller near the wall because of the boundary layer. Figure 5.10 - mesh definition near the glass In order to have the same precision, the element size is roughly the same, and not the total elements. This mesh has 70% fewer elements because there is no mesh on the products, which represents 30% of the total area. 5.2.2. Results Figures 5.11 and 5.12 represents the simulation without and with the lamps’ influence, respectively. The lamp is collocated near the air entrance on the top. Figure 54 5.13 represents the velocity field of both situations (the velocity results are the same, as expected). Figure 5.12 - geometry with products, velocity field after 120s, T0 = 15oC, v0 = 2m/s, T(inlet) = 0oC, Q(glass) = 8.29 W/m2K, T(bottom) = 0oC, with double sized entrances on the bottom of the perforated plate As it can be seen in Figures 5.11 and 5.12 the influence of the lamp is local, as it happened on the heat flux through the glass. The only propagation of the heat generation by the lamp is to the air curtain. As can also be seen on these figures, between the products there is almost no air flow, except on the top shelf between the first and third product. That means that in this zone the refrigeration will be reduced. Figure 5.12 - lamps influence, geometry with products, temperature field after 120s, T0 = 15oC, v0= 2 m/s, T(inlet) = 0oC, Q(glass)= 8.29 W/m2K, T(bottom) = 0oC, with double sized entrances on the bottom of the perforated plate Figure 5.111 - geometry with products, temperature field after 120s, T0 = 15oC, v0 = 2 m/s, T(inlet) = 0oC, Q(glass) = 8.29 W/m2K, T(bottom) = 0oC, with double sized entrances on the bottom of the perforated plate 55 So, it can be concluded that 20 mm between the products is not enough to ventilate this space, a larger space was needed as there is on the top shelf, which is confirmed by figure 5.13 (the velocity is nearly 0 in the 20 mm gap). However, there are a lot of limitations in the 2D simulations. Firstly, the gap between the products all the way across cannot be represented. That means that the air could be flowing the products, contrary to what is represented in oud 2D results. Secondly the inlet and outlets, so as the gaps in the back of the showcase, are represented as lines in 2D, which would be transposed to a rip all the way across. In the real geometry, the inlets, outlets and back perforation are often circular holes or oblongs, which can be only designed in a 3D geometry. These limitations in the 2D geometry can leave us to some conclusions that are not true. 56 6. CONCLUSION AND FUTURE WORK In this chapter some conclusion will be presented and the future work that would best complete this one. Conclusions This works contains some simulations with different conditions, where the first one is the simplest because it contains only the geometry itself, and the next simulations increase one more condition that were taken into account. The point of analysing the results this way is to visualize which condition has more impact on the temperature and/or velocity field and it becomes easier to see the errors on definition. If only the last case was made, it was not possible to verify the influence of each parameter. The first case helped us to define the mesh definition for the following ones. At the beginning were used more than 1 mesh because it was not certain that the complexity of the next simulations would have different results, but it was concluded later that was not necessary. The heat flux through the wall and the influence of the lamps’ studies are similar in the way that it only affects locally for the values used. On the other hand, defining the temperature as 0°C on the bottom wall (near the evaporator) affected a large area compared with the heat flux definition. The last step on the empty showcase was to double the size on the bottom of the back perforated plate in order to increase the velocity on this area. It was concluded that by increasing the mass flow on the bottom, it will be reduced on the other entrances, causing to reduce the velocity of the airflow near the glass (air curtain). This might be interesting if depending on the application: if the products are packed, the first option is better because the temperature is lower overall, if there are products that need to be ventilated to maintain is properties, such as fruits and vegetables, it is better to use the second option. 57 Lastly, the geometry of the product is very limited because the products are represented by boxes, which have the same width as the showcase, that is not the real case, where the products are placed side-by-side, having enough space to the air circulate all the way across. This causes that the velocity is zero between the products, which is not true neither in width nor in length. For those reasons, where the products are included is suitable to develop a 3D simulation in order to get reliable results. Future work The future work to complement this one would be a 3D simulation that can represent better the reality and it would include more geometry changes, such as air entrances and exits, product dimensions and glass configuration, in order to optimize the configuration of the final equipment. The velocity inlet could also be analysed in order to minimize to reduce costs, always guaranteeing the temperature field remain acceptable. Another very important parameter that should be analysed is considering the products as a thermal load and not just an obstruction to the flow. As a final task, it could be done a prototype to validate all the results. 58 7. REFERENCES [1] N. Chaomuang, D. Flick, and O. Laguerre, “Experimental and numerical investigation of the performance of retail refrigerated display cabinets,” Trends Food Sci. Technol., vol. 70, no. May 2017, pp. 95–104, 2017. [2] P. Lobarinhas, “Ciclos Frigoríficos,” in Ciclos frigoríficos, 2018. [3] Leitão, B.M.R., “Equipamentos de congelação industrial de produtos alimentares perecíveis: análise comparada de apoio à decisão,” Instituto Superior de Engenharia de Lisboa, 2015. [4] PRODCOM, “NACE 13.10 : Mining of iron ores PRODCOM LIST 2006.” p. 322, 2006. [5] S. C. Ribeiro, “Estudo numérico de escoamentos em equipamento de refrigeração vertical aberto,” Universidade do Minho, 2016. [6] N. Chaomuang, D. Flick, A. Denis, and O. Laguerre, “Experimental analysis of heat transfer and airflow in a closed refrigerated display cabinet,” J. Food Eng., vol. 244, no. April 2018, pp. 101–114, 2019. [7] L. Wang, L. Zhang, and G. Lian, “A CFD simulation of 3D air flow and temperature variation in refrigeration cabinet,” Procedia Eng., vol. 102, pp. 1599–1611, 2015. [8] J.M. Belman-Flores, “Analysis of the flow and temperature distribution inside the compartment of a small refrigerator,” Appl. Therm. Eng., vol. 106, pp. 743–752, 2016. [9] H. K. Versteeg and W. Malalasekera, versteeg_computational_fluid_dynamics (2).pdf. 1996. [10] ANSYS, “ANSYS FLUENT Theory Guide.” 2013. [11] P. D. Gaspar, A. Miranda, and R. A. Pitarma, “Estudo numérico do escoamento de ar e do desempenho térrmico de equipamentos e expositores refrigerados,” Métodos Numéricos Ing. V, p. 19, 2002. [12] D. Parpas, C. Amaris, J. Sun, K. M. Tsamos, and S. A. Tassou, “Numerical study of air temperature distribution and refrigeration systems coupling for chilled food processing facilities,” Energy Procedia, vol. 123, pp. 156–163, 2017. [13] V. Guedes and D. Brás, “Simulação do escoamento numa vitrine de frio,” 2018. [14] J. Sun, K. M. Tsamos, and S. A. Tassou, “CFD comparisons of open-type refrigerated display cabinets with/without air guiding strips,” Energy Procedia, vol. 123, pp. 54–61, 2017. [15] F. Casa, “No Title,” O vidro e o isolamento térmico, 2019. [Online]. Available: https://forumdacasa.com/discussions/?PostBackAction=Download&Attachment ID=1069: . [Accessed: 21-Oct-2019].