Numerical analysis of a cooling of a hot plate by an array of microjets

. In this paper heat transfer during impingement of an array of 8 ¥ 8 microjets on a hot surface was investigated. The influence of a ratio of a distance between a nozzle and hot plate ( H / d ) and microjet diameter-based Reynolds number (Re d ) on the temperature and heat transfer coefficient (HTC) on the hot plate were numerically studied. The numerical model which was based on the steady-state compressible Navier-Stokes equations and SST k - w turbulence model was developed and applied for the analysis. During simulations the ratio of the distance between the nozzle and hot plate to the microjet diameter was H / d = 3.125, 25 and 50, while the microjet diameter-based Reynolds number was equal to Re d = 690, 1100 and 1510. The ratio of the microjet pitch to the microjet diameter was s / d = 31.25. It was found that both the H / d ratio and Re d significantly influenced flow patterns in the gap between the nozzle and hot plate as well as the temperature and HTC on the surface of the hot plate. With increase of the H / d ratio a more uniform distributions of the plate temperature and HTC were observed, while the rise in the Re d intensified heat transfer on the hot plate.


Introduction
The microjet impingement is gaining attention as it has large potential for applications in cooling technologies for miniaturized devises [1], e.g., in electronics. It is very efficient transport phenomenon which offers very high heat and mass transfer rates on an impinged surface. However, to cover large surfaces, many microjets are joined in an array [2].
Many investigations were conducted for conventional jet impingement to gain knowledge about these phenomena and to assess their applications and effectiveness in many fields of industry such as electronics cooling, drying of paper, glass, lumber and textiles, cooling of turbine blades, food processing and annealing of metals [1,2]. Still, due to many factors which influence on heat transfer and fluid flow during this phenomenon, e.g., a type of fluid, thermophysical properties of the fluid, flow rate and geometrical parameters of the cooling system, further investigations are required.
Many experimental and numerical investigations have been carried out for conventional jets impinging upon the surface to find their fluid and heat flow characteristics. Various types of geometries of the jet which impinged the surface, e.g., single round jet [3][4][5][6], twin slot jets [7], single row of jets [8] as well as elliptic, circular and rectangular jet arrays [9] were considered. Influences of various geometrical parameters, e.g., the ratio of the jet height to the jet diameter [3,4,6,7,9], the inclination angle of the jet to the impingement surface [7], the ratio of the jet to the jet separation distance [7,8] and the ratio of the cross-flow velocity to the jet velocity [8] were investigated. The Reynolds number based on the jet diameter in these analyses varied approx. from 1000 up to 100000. Different ways of turbulence modelling were tested, e.g., standard [6] and low Reynolds number k-e models [6,9], RNG k-e model [6,7], realizable k-e model [6], standard k-w model [6], SST k-w model [6][7][8], RANS k-w model of Wilcox [3], hybrid RANS/LES models with different subgrid scale eddy viscosity in the LES mode [3,4] and direct numerical simulation [5].
In the case of impingement of miniaturized jets or microjets several experimental and numerical works were conducted. These investigations concentrated on both characterization of heat transfer and fluid flow during microjet impingement [10] as well as investigation of performance of devices based on the microjet impingement technology, e.g., cooling system of gas turbine blade [11], silicon-based heat sink for electronic applications [12], modular heat exchanger [13] and cooling system integrated within an electronic devices [14].
In this paper numerical parametric analysis of the cooling system based on 64 microjets, arranged in the rectangular array of 8¥8, and impinged on the hot surface was carried out. At first considered microjet based cooling system as well as developed mathematical and numerical models were described. Then numerical test carried out were presented and the results obtained discussed. Finally, conclusions were formulated.  2 Thermal model of an array of microjets impingement

Geometrical model
The considered cooling system which based on microjets impingement is presented in Fig. 1. It consisted of an airgap between the nozzle and hot plate where air microjets expanded, impinged on the heater surface and then flow away, through lateral openings, cooling hot surface. Microjets were generated in orifices, i.e., outlets from circular micro-channels in the nozzle. The nozzle consisted of 64 micro-channels of diameter d = 0.16 mm which were arranged in the rectangular array of x and y pitches of s = 5 mm. For this configuration the ratio of the microjet pitch to the microjet diameter was s/d = 31.25. Because considered geometry has at least two symmetry planes only its quarter denoted by red colour in Fig. 1A was considered during preparation of the geometrical model. The computational geometry and assumed boundary conditions (i.e., microjets outlets, lateral air outlets and symmetry planes) are presented in Fig. 1B. This geometry corresponded to the one-fourth of the heater's surface. In this way, the whole space above the heater was modelled. The computational domain had the shape of cuboid of dimensions: 20¥20¥H mm, where three heights of the airgap between the nozzle and hot plate were assumed, i.e., H = 0.5, 4 and 8 mm (see Fig. 1C). These heights corresponded to the following values of the H/d ratio: H/d = 3.125, 25 and 50, respectively. Distances from outer rows of microjet outlets to symmetry planes and lateral air outlets were equal to 2.5 mm.

Computational grids
Meshes were generated applying the sweep meshing technique in which the face mesh (i.e., created on the nozzle surface) was extruded in the direction perpendicular to the nozzle and heater surfaces. Meshes were also refinement in flow sensitive regions, i.e., close to microjets outlets and impingement surface. Moreover, in the airgap the fine mesh was generated close and in zones defined by cylinders which corresponded to the microjet outlets. These cylinders embraced microjets. In this way meshes were refined in areas of the highest shear stresses where microjets were generated and expanded. Despite relatively simple geometry the obtained number of elements were very large and equal to approx. 37.22, 36.75 and 43.23 millions for H/d = 3.125, 25 and 50, respectively, due to mesh refinement. For generated grids mesh quality parameters were the following: maximum values of the aspect ratio and skewness were equal to 42.86 and 0.737, 180.0 and 0.749 as well as 184.6 and 0.699 for H/d = 3.125, 25 and 50, respectively. Moreover, y + values obtained on the hot plate where microjets impinged were below y + < 2. Tests carried out showed that these values of mesh quality factors were within acceptable ranges and ensured grid independent solutions.

Governing equations
Air flow in micro-channels in the nozzle was generated by significant pressure difference between its inlet and outlet. Resulting velocities of microjets at outlets from micro-channels and at inlets to the airgap between the nozzle and hot plate were very high (i.e., approx. 87, 130 and 169 m/s for Red = 690, 1100 and 1510, respectively). Therefore, steady-state compressible air flow was assumed in the airgap and was described by the following continuity, momentum, energy and SST k-w turbulence model equations [15][16][17]: where: i, j, k = 1, 2, 3 are indices of Cartesian axes, Dwcross-diffusion term, E -total energy, Gk -generation of k due to mean velocity gradients, Gw -generation of w, k -turbulence kinetic energy, p -pressure, Ttemperature, u -velocity component, Yk and Ywdissipation of k and w, respectively, due to turbulence, Gk and Gw -effective diffusivity of k and w, respectively, d -Dirac delta function, l and lt -molecular and turbulent thermal conductivity, respectively, µ -dynamic viscosity, r -density, w -specific dissipation rate. The air parameters were calculated from ideal gas model for which the total energy was given by the following formula [15][16][17]: where cp is the specific heat. The viscous (Reynolds) and deviatoric stress tensors were found by applying the Boussinesq hypothesis from following equations [15][16][17]: where: u¢ denotes velocity fluctuation and µt is the turbulent (eddy) viscosity which was calculated using solution of the SST k-w turbulence model.

Boundary conditions and data assumed in simulations
Boundary conditions assumed in simulations are shown in Fig. 1B and C. Velocity, turbulence and temperature profiles set at each microjet outlet were generated in separated simulation of air flow through the single micro-channel in the nozzle. In all numerical simulations the same mass flow rate of air in each micro-channel was assumed. Total volumetric flow rates of air at temperature 300 K and supplied to the considered cooling system were assumed equal to 5, 8 and 11 dm 3 /min which corresponded to constant mass flow rates of 1.512◊10 -6 , 2.420◊10 -6 and 3.327◊10 -6 kg/s, respectively, through each of 64 micro-channels. These mass flow rates referred to Red = 690, 1100 and 1510, respectively. On the surface of the heater the fixed value of heat flux density equal to 3750 W/m 2 was assumed which corresponded to 6 W of heat dissipated in the whole cooling device. Through lateral outlets air flowed away from the cooling system to the surroundings at 300 K and 1 atm. Thermophysical properties of air used in simulations were following: the specific heat at constant pressure was equal to cp = 1006.43 J/kg/K, the thermal conductivity of l = 0.0242 W/m/K, the dynamic viscosity of µ = 1.7894◊10 -5 kg/m/s and the molecular weight of M = 28.966 kg/kmol.

Numerical solution
The set of governing equations, eq. (1)-(5), were solved in the Finite Volume based ANSYS Fluent software. Settings of the density-based solver were the following: second order upwind schemes for approximation of convective terms in the Navier-Stokes, energy, turbulent kinetic energy and specific dissipation rate equations, eq.
(1)-(5), as well as the implicit solution's formulation. Moreover, the least square cell based gradient reconstruction method was used. To avoid excessive generation of the turbulence energy in the vicinity of stagnation points production limiter was applied to limit the production term in the turbulence kinetic energy equations, eq. (4). Simulations were performed utilizing a multi-processor workstation. Due to very large number of mesh elements times of simulations applying parallel computing (i.e., 12 nodes per case, Intel® Xeon® CPU E5-4620 v4 @ 2.10 GHz) approached one month for each case for convergence criteria on the level of 10 -5 for all equations. The RAM usage exceeded 200 GB.

Results
Numerical simulations were performed for three values of the ratio of the distance between the nozzle and hot plate to the microjet diameter, i.e., H/d = 3.125, 25 and 50 as well as for three values of the microjet diameterbased Reynolds number, i.e., Red = 690, 1100 and 1510. Results of calculations, i.e., distributions of temperature and HTC on the hot plate are presented on Fig. 2 and 3, respectively, for three values of the H/d ratio and three values of the Red. Bottom left corners of maps of temperature and HTC refer to the centre of the heater plate, whereas their top and right edges to lateral air outlets. In these plots the areas of the microjet impingement and flow stagnation might be easily identified -see blue and red colour regions in Fig. 2 and red and blue colour ones in Fig. 3. The temperature of the hot plate varied from about 300 K up to 375 K -see Fig. 2 The lowest average temperature of the heater surface was found for H/d = 25 and was equal to 322.8, 315.5 and 314.6 K for Red = 690, 1100 and 1510, respectively -see Table 1 with average temperatures and HTCs on the hot plate. This means that for H/d = 25 the best compromise between the cooling effectiveness and the size of affected area by microjets was achieved. Further increase of the H/d ratio resulted in the rise of the average heater temperature and the decrease of the cooling efficiency. The degradation of the velocity profile with the increase of the microjet height was responsible for this effect. Moreover, rising in the H/d ratio resulted in the influence of the microjet on adjacent microjets. For H/d = 3.125 and 50 average temperatures of the hot plate were approx. 13-17 K and 1-3 K higher, respectively than for H/d = 25. Values of the lowest plate temperature were also depended on the Red. If the Red grew, both the lowest and average temperature of the heater surface fell regardless the height of the microjet.
Distributions of the HTC presented in Fig. 3 were calculated assuming that the reference temperature was equal to 300 K. Local values of the HTC varied from several up to 3910 W/m 2 /K. Similarly as for temperature, in Fig. 3 the microjet impingement areas are denoted by red and green colour which corresponds to the highest local HTCs on the heater plate, while flow stagnation zones are denoted by blue colour which corresponds to the lowest local HTCs on the plate. Values of the lowest and highest local HTCs were depended on both the H/d ratio and Red. If the H/d ratio increased, the maximum value of the HTC fell. Again, the dependency between the influence area by the microjet and its height was observed. In terms of heat transfer, the best ratio of the distance between the nozzle and hot plate to the microjet diameter was found to be H/d = 25. For this value the highest average HTCs on the heater surface were observed -see Table 1 The effect of the lateral air flow in the airgap between the nozzle and hot plate on the temperature and HTC on the hot plate is also clearly visible and was the most significant for higher values of the H/d ratio and Red. With increase in the distance between the nozzle and hot plate the effect of lateral flow of air increased as well -see Fig. 2 and 3. What is more, the outer microjets cooled the hot surface more efficiently than inner ones. Due to the lateral flow in the airgap the relative velocities for outer microjets were slightly higher than for inner ones (i.e., the lateral flow was perpendicular to microjet axes) which resulted in the more intensive heat transfer rates close to lateral air outlets. This effect might be also related to the suction of air from inner microjets by the outer ones which resulted in their higher impinging velocities. The influence of the lateral air flow in the airgap between the nozzle and hot plate on the shape of impingement area is more visible for variation of the Red than for change of the H/d ratio. Rising of the Red resulted in more intensive flow between microjets and therefore led to increase in heat transfer intensity and decrease of hot plate temperature in these regions.

Conclusions
In this work heat and fluid flow phenomena during impingement of an array of 8¥8 microjets on a hot surface were numerically investigated. The influence of two parameters, i.e., the distance between the nozzle and hot plate (H/d) as well as the microjet diameter-based Reynolds number (Red) on both the local and average temperature and HTC on the surface of hot plate were studied. During development of the computational model the steady-state compressible air flow was assumed. To account for turbulence the SST k-w turbulence model was applied. Numerical simulations were performed for three values of the H/d ratio, i.e., H/d = 3.125, 25 and 50 and for three values of the Red, i.e., Red = 690, 1100 and 1510. It was found that both parameters have significant affect on the heat and fluid flow phenomena during impingement of the array of microjets on the hot surface.
In the case of variation of the H/d ratio, its optimal value for which the highest heat transfer intensity occurred was found. From analysed cases the ones for H/d = 25 were characterized with the lowest average temperatures and the highest average values of the HTC on the heater surface. In the case of variation of the Red, its increased resulted in the decrease of the average temperature as well as the increase of the average HTC on the heater surface regardless the height of the microjet. Moreover, with the rise in the H/d ratio, more uniform distribution (i.e., lower differences between maximal and minimal values) of temperature and HTC on the hot plate were observed. For the small H/d ratio, microjets with air at high speed impinged on the hot surface and the impingement spots were small. Rising in the H/d ratio resulted in the expansion of microjets and degradation of their velocity profiles. Velocity of the air in microjets decreased and impingement area greatly enlarged covering almost the whole hot plate and affecting on adjacent microjets. The H/d ratio and Red had also impact on the lateral flow in the airgap between the nozzle and hot plate. For small values of the H/d ratio and Red impingement spots were not influenced by flow of air which left the cooling system by lateral openings. Increasing of the H/d ratio and Red resulted in interactions of the lateral flow of air with microjets. Moreover, outer microjets cooled the hot surface more efficiently than inner ones due to their higher relative velocities.
This work has received funding from the the Faculty of Power and Aeronautical Engineering of Warsaw University of Technology in the framework of statutory activity.