Application of CFD Simulation in the Development of a New Generation Heating

This paper deals with the application of Computational Fluid Dynamics simulation in the development of a new generation cooking appliance in Gorenje concern. As the oven is multifunctional, radiation, conduction, natural and forced convection mechanisms of heat transfer are used. The Discrete Ordinate (DO) model is used for radiation. The density of air is described by incompressible ideal gas equation in a natural convection model. The intention was to create the best possible baking conditions for different heating systems. Several discrete models were created. The influence of geometry change and boundary conditions variations to the velocity and temperature field distribution in the oven cavity was analyzed. The results of numerical simulations are validated with measurements taken from an oven prototype. The agreement was good. After successfully passing the standard tests, the oven came into serial production and was launched on the market.


INTRODUCTION
Consumer society, fashion trends, new materials and demands for efficient energy use require an increasingly rapid development of new products.This also applies to ovens, which have long been an indispensable piece of kitchen equipment.
In the past, heating ovens used solid fuel or gas, but today the main source of energy is electricity.Electric ovens are generally combined, which means that there is a choice of three different baking methods: classic, fan and grill.In the traditional mode, the electric heaters are located on the floor and on the ceiling of the oven.Heat is transferred [1] and [2] by conduction, natural convection and radiation.In the fan mode, the air is heated by electric heaters and it circulates in the oven as a result of a fan.The main mechanism of heat transfer is forced convection.In the grill mode the electric heater on the oven ceiling is heated to such high temperatures that it glows.Heat is transferred by radiation.Each of these heating methods is specific.In developing the new oven we tried to provide optimal conditions for all the three methods of heating.
Computational fluid dynamics [3] to [5] is very useful for predicting temperature and velocity fields in the oven.With the rapid progress of parallel computing and development of commercial codes, CFD has become an indispensable tool in the development of new products.Together with the experiment it allows for a significant reduction of time from the first design through the prototype to the finished product.Therefore, CFD is increasingly used in the food industry and in developing appliances used in food preparation.This is proved by a growing number of publications in recent years.
In their study, Amanlou et al. [6] present the use of CFD for design cabinets for fruit drying.They investigated the effect of different geometries on the efficiency of the cabinet dryer.They also tried to find the best construction with regard to the uniform distribution of temperature and air flow.Comparing the experimental and CFD, data have revealed a very good correlation coefficient for drying air temperature and air velocity in the drying chamber.
Chhanwal et al. [7] use a CFD model of electric oven for baking bread, where the radiation is the dominant mechanism of heat transfer.They employed three models of radiation: discrete transfer radiation model (DTRM), surface to surface (S2S) and discrete ordinates (DO).The results were compared with the experiment.They reported a good agreement for all three models.Boulet et al. [8] made CFD model of the bakery pilot oven.All three heat transfer mechanisms are considered and coupled with turbulent flow, for which they used k-ε realizable model whereas the S2S model simulates the radiation.The model predictions show a good qualitative agreement with the experimental measurements.
In their article Verboven et al. [9] discuss the application of CFD modeling and validation of the isothermal air flow in a forced convection oven.The governing fluid flow equations were expanded with a fan model.They used a standard and renormalization group (RNG) version of the k-ε turbulence model.The calculated speeds were validated by measurements with a hot-film velocity sensor.The calculation error was 22% of the actual velocity.The reason for the relatively high error was within the limits of turbulence modeling and numerical grid density.
Therdthai et al. [10] performed a threedimensional CFD simulation of temperature profiles and flow patterns during a continuous industrial baking process.They used a moving grid approach.Their model could predict the dynamic responses.The model was also used to investigate the oven operating conditions, which could produce optimum baking conditions.The numerical results were consistent with the actual measurements.
Williamson et al. [11] presented the development, testing and optimization of the design of a novel gasfired radiant burner for industrial tunnel ovens.CFD simulation was used to model the burner and baking chamber environment, and in particular to predict radiation heat fluxes incident on the top surface of the food.CFD model has been validated by measurements with thermocouples.The agreement was within 10%.The simulations have indicated that the new burner is capable of delivering irradiation to a traveling conveyor more uniformly than existing radiant burner designs.
Wong et al. [12] used the 2D CFD model of the oven chamber of an industrial continuous breadbaking oven in order to better understand the process of baking.They used the sliding mesh technique to simulate the continuous movement of dough/bread during the whole traveling period.The CFD modeling was proven to be a useful approach in studying the unsteady heat transfer in the oven as well as the heating history and temperature distribution within dough/bread.
Ghani et al. [13] numerically simulated the natural convection heating of canned liquid food for axisymmetric case.They predicted transient flow patterns and temperature profiles.The density and viscosity of the fluid were dependent on temperature.It was found that the slowest heating zone moves towards the bottom of the can.The shape and size of this area depends on the used liquid.
The optimization of new generation ovens NGKA3 was performing a numerical simulation of temperature and velocity fields.Numerical simulation of such a complex case was made on mini computational cluster Supermicro (5 nodes, 14 processors, 38 cores, 104 GB RAM) with the ANSYS Fluent software.
The numerical simulation was carried out in several steps.The solutions of transport equations enabled the study of pressure, velocity and temperature fields under steady conditions.The second step was the validation [14] of the model.The agreement between the numerical simulation results and the measured data from the prototype was checked.After the numerical model was validated, the optimal temperature distribution in the oven cavity by changing the geometry (heaters shape, shape of the fan cover) and the boundary conditions (heater temperature, ventilation system flow rate) at the third step were looked for.

Governing Equations
Temperature distribution and air circulation in the oven cavity are governed by conservation equations.The system of partial differential equations is: • momentum conservation: • energy conservation: • radiation intensity conservation: where  v is the velocity vector, τ is the stress tensor, ρ  g and  F are the gravitational body force and external body forces respectively, e is the energy, h is the sensible enthalpy, ρ is the density, p is the pressure, T is the temperature, λ is the thermal conductivity, c p is the specific heat and S h is the volumetric heat source (0 in our case), I is the radiation intensity,  r is the position vector,  s is the direction vector,  s ' is the scattering direction vector, s is the path length, a is the absorption coefficient, n is the refractive index, σ s is the scattering coefficient, σ is the Stefan-Boltzman constant.
As we were only looking for a steady state, the governing equations are simplified since all time derivatives are zero.
The Reynolds number varies between 230 (inlet) and 880 (outlet) in the case of forced convection, while the Rayleigh number is 6•10 8 for natural convection.These numbers indicate [2] a transition zone between laminar and turbulent flow.We have decided to treat flow as a laminar since Reynolds averaged (RANS) turbulence models available in Fluent do not satisfy the y+ wall conditions for such a case.
The pressure variations are small enough, so the incompressible ideal gas law is used to express the relationship between air density and temperature: where p op is the operating pressure (101.325Pa), R is the universal gas constant and M w is the molecular weight of the air.This means that full buoyancy model was used.
The following material properties presented in Table 1 were used in the numerical model.

Choosing a Radiation Model
At high temperatures, the fourth-order dependence of the radiative heat flux on temperature implies that the radiation will dominate and the radiative heat transfer should be included in the simulation.ANSYS Fluent provides five radiation models: Discrete Transfer (DTRM), P-1, Rosseland, Surface-to-Surface (S2S) and Discrete Ordinates (DO).To choose the most accurate radiation model a test simulation was done.
Maximum temperature of the biscuit surface was compared with the measured temperature, Table 2. Quite large differences indicate that not all models are appropriate for our simulation.The optical thickness aL is a good indicator of which model to use, where L is the length scale of the domain.In our case (L = 0.4 m) the optical thickness is <<1.This comparison proves that the P-1 and Rosseland models are not suitable for optically thin problems.Although DTRM and S2S models can be used in these kinds of problems, they do not support semi-transparent walls (glass doors).This is manifested in over-prediction of biscuit surface temperature.The only appropriate radiation model for our analysis is DO model which works across the full range of optical thicknesses, and it also allows the solution of radiation at semitransparent walls.

Influence of Biscuit Surface Emissivity
The oven inner walls are enameled, while heaters are made from steel.The emissivities of enamel and steel are well known and were supplied by manufacturers.On the contrary, the emissivity of biscuit is not well known.Computations with ε in the range from 0.4 to 0.8 were conducted to study the influence of biscuit wall emissivity to maximum biscuit surface temperature, Table 3.As it can be seen, the biscuit surface emissivity has no significant effect on the surface temperature.

Discrete Models
Four different geometrical models were created in GAMBIT [15] using the combination of various geometrical parameters, see    The manual repair by merging of problematic cells was needed.The maximal skewness was 0.76 for all grids after mesh repair.This means that the quality of the surface mesh was very good.The initial volume grids were generated on the basis of surface meshes.They needed additional refinement and smoothing to achieve the grid quality within the recommended limits.
For a more effective calculation (smaller number of cells, faster and more stable convergence) the tetrahedral cells were converted to polyhedral cells in Fluent [17], see Fig. 5. Table 5 shows the size of computational grids for the four models after conversion of tetrahedral cells.The quality of polyhedral grid was checked [18] with Fluent before starting the numerical simulation.In addition to the above parameter skewness, the criterion of aspect ratio (a measure of the stretching of the cell) and squish index (using the vector from the cell centroid to each of its faces and the corresponding face area vector) are used.Table 6 shows the parameters of cell quality for all four models.Maximum values are much lower from the recommended upper values, which means that the spatial discretization is made very well.
The mesh consistency test was not performed due to the short deadline.The choice of computational mesh size was based on experience with similar cases and available computer resources (memory and CPU time).

Boundary Conditions
The correctness of the numerical calculation results is largely dependent on boundary conditions.If these do not correspond to the actual situation, the solution of equations does not give the expected values.Temperature boundary conditions for oven outer walls, infrared heaters and fan cover were obtained by the measurements on the oven prototype in the Electrothermics Laboratory of Cooking Appliances department in Gorenje to approach the real situation, see Fig. 6.Thermostat was set to 200 ˚C.
The measured average temperatures at the oven outer walls in the case of conventional heating are shown on Fig. 7.The influence of the heater switching with time period of 10 minutes on the top and bottom wall can be observed.After 30 minutes oscillations of temperature are minimal on the walls away from heaters.A similar scenario was observed in an oven with convection heating.Steady state conditions are established after 60 minutes of heating and they are used for temperature boundary condition at oven walls, see Table 7.
The maximum temperatures of the last three periods are used for the top and bottom wall boundary conditions.The oven has a fan to provide hot air circulation in the case of convection heating system.The air is sucked from the oven cavity and after it is heated it travels through the slots back to the oven cavity.Due to the complexity, the fan is not modeled in details.
At the perforated area (air exit from oven cavity) the pressure drop is prescribed.It is proportional to the fluid dynamic pressure: where the resistance coefficient [19]    Fan cover slots, where hot air enters into the oven cavity, are divided into 16 segments, see Fig. 8. Eleven different combinations of boundary conditions have been prescribed for segments, where some were closed and others open.The measured air velocity and air temperature were 0.4 m/s and 250 °C, respectively.

Numerical Simulation
Four numerical simulations have been done.Only heat transfer with radiation of hot walls (infrared heater, oven cavity walls, glass door and tray with biscuit) was considered in the first case (Model 0).Forced convection was included in other simulations (Models 1, 2 and 3).
The commercial CFD code ANSYS Fluent was used for the numerical simulation.The code uses a pressure-based solver with SIMPLE [20] method for velocity-pressure coupling.The relaxation factors were 0.3 for pressure, 0.7 for velocity, and 1 for energy and radiation.The velocity and temperature fields were discretized with a second order upwind scheme, whereas the pressure field was discretized with a PRESTO!Scheme.The convergence criteria for residuals of continuity and momentum equations was 10 -4 and 10 -6 for energy and radiation equations.

Validation
Validation of the numerical model was made for a conventional heating system.Fig. 9 shows the positions of five measuring probes of temperature on a tray with biscuit.The calculated temperature field on the biscuit surface is shown in Fig. 10.

Fig. 9. Positions of measuring probes
A comparison of measured and calculated temperatures at selected points is shown in Table 8.Peak temperature in the middle of the tray (measuring point 3) matched very well with the measurement within 1%.A slightly higher difference in temperature, up to 15%, was measured at the front region (measuring points 1 and 2), which is probably a result of inexact boundary condition for the glass doors.The deviation of temperatures at measuring points 4 and 5 was up to 10%.

Optimum Conditions in the Oven Cavity
Several factors have an influence on the conditions inside the oven, such as: the shape and power of heaters, fan rotor rotational speed, thickness and quality of insulation, the design of oven doors, etc. Optimum conditions can be approached through the variation of these factors.The change in the oven cavity geometry and the variation of the boundary conditions was used in this case.We varied the inclination of the back wall and the width of the air inlet slots.Boundary conditions were changed by opening or closing the segments of the air inlet slots.The numerical simulation results for the convection heating system variants are presented below.

Geometry Change
The change in the form of oven cavity has a strong influence on the velocity field, see Fig. 11.The inlet velocity of 0.7 m/s is high because of a 5 mm slot width in the initial geometry of Model 1.The straight back wall directs hot air perpendicular to the side walls, resulting in hot air recirculation in the back corners and its rapid suction from the oven cavity due to the fan proximity.The zone of cooler and slower air appears at the front region.Such a situation means worse baking conditions.By changing the inclination of the back wall and increasing the width of the inlet slot, the conditions are significantly improved.Due to a better flow guidance the movement of hot air is assured throughout the oven cavity.

Change of Boundary Conditions
Figs. 12 and 13 show the velocity vectors and temperature distribution in the vertical plane in the middle of the oven for the initial Model 2 (variant 00, all segments are open) and improved conditions (variant 10, closed segments 1 to 4 and 5 to 9).
Variant 00 shows the above vortex, which causes the flow of hot air to be directed down in the middle of tray, from where the fan sucked it from the oven cavity.Therefore, the temperature field around the biscuit is non-uniform, which means poor baking at the front where there is already high heat transfer due to the glass doors (radiation, convection).Temperature is about 20 °C higher than the desired 200 °C in a great portion of the oven cavity.
Due to the top and bottom air inlet slots closure in variant 10, the vortex does not appear.Hot air flows from the glass doors along the tray towards the fan.Therefore, the air temperature in the vicinity of biscuit is constant, which results in even baking.The temperature throughout the oven cavity is only a few degrees higher than 200-°C.
Figs. 14 and 15 show the velocity vectors and the temperature distribution in the horizontal plane at the middle of the oven cavity for variants 00 and 10 of the Model 2. The air circulation at the rear of the oven can be observed again in the case of initial model.This means higher temperature in that region and consequently, uneven baking.
At the variant 10, the top and the bottom air inlet slots are closed, therefore the air flow is increased at the side slots.Hot air travels along oven side walls in the direction of doors, where it turns back towards the fan on the back wall.As can be seen from the picture of isotherms, the temperature field is much more uniform as in the case of variant 00.
Fig. 16 show streamlines, where color represents temperature, and biscuit surface temperature for the initial and the optimized variant of the fan cover slots boundary conditions.

CONCLUSION
The velocity, pressure and temperature field in the oven cavity have been dealt with the use of CFD simulation.Natural convection, forced convection, conduction and radiation mechanisms of heat transfer were used.The comparison of five radiation models showed that only appropriate model is DO.Biscuit surface emissivity didn't show significant effect on the surface temperature.Four different oven geometries and eleven variations of velocity boundary conditions were used.
Control of the circulating hot air velocity in the oven cavity has emerged as the decisive impact factor on the baking quality.Predictions from numerical simulations are confirmed by the positive results from the validation and the functional testing of the oven prototypes and they allowed fast optimization of the temperature field in the oven cavity.

Fig. 1 .
Fig. 1.Conventional oven1.5 Computational MeshDiscretization of the geometrical models was performed in the program Tgrid [16] by using tetrahedral control volumes.The grid quality parameter skewness (the difference between the shape of the cell and the shape of an equilateral cell of the equivalent area) of the automatically generated surface mesh exceeded a value of 0.85 in some cells.The manual repair by merging of problematic cells was needed.The maximal skewness was 0.76 for all grids after mesh repair.This means that the quality of the surface mesh was very good.The initial volume grids were generated on the basis of surface meshes.They needed additional refinement and smoothing to achieve the grid quality within the recommended limits.For a more effective calculation (smaller number of cells, faster and more stable convergence) the tetrahedral cells were converted to polyhedral cells in Fluent[17], see Fig.5.Table5shows the size

Fig. 7 .
Fig. 7. Measured average temperatures at the oven outer walls

Fig. 8 .
Fig. 8. Segments of the fan cover slots

Fig. 6 .
Fig. 6.Temperature measurement locations at oven outer wall and infrared heater

Table 1 .
Physical properties of the materials in use

Table 3 .
Maximum biscuit surface temperature at different emissivities

Table 2 .
Maximum biscuit surface temperature for different radiation models

Table 5 .
Number of tetrahedral and polyhedral cells

Table 7 .
Temperature boundary conditions at oven walls

Table 6 .
Parameters of cell quality for polyhedral mesh

Table 8 .
Comparison of measured and calculated temperature