The Development of a Mathematical Model of an Evaporative System with Film Flows

A spatial mathematical model of the evaporation system using film flows has been developed. A numerical study of several modes of film flow in an evaporation cell at a heater temperature below the boiling point has been carried out. The characteristics of the film flow and the features of the coupled heat and mass transfer processes in the gas-film-metal substrate system have been studied. Good agreement between the calculation and the experiment on the main parameters of the wave motion of the film and the heat transfer characteristics has been obtained.


Introduction
Cooling of microelectronic equipment is one of the most important problems of thermal physics today. Existing electronic microchips emit heat fluxes at the level of 100 W/cm 2 . It is expected, that heat fluxes in the equipment hot spots will approach 1000 W/cm 2 in future due to the development of microprocessor technology. Therefore, the problem of heat dissipation is one of the keys in modern electronics, and the further increase in microprocessor performance depends on its solution in many respects. Traditional heat pipes combined with a fan allow dissipating the heat fluxes of about 200 W/cm 2 , which is obviously insufficient. In this regard, microchannel heat exchangers with characteristic flow size less than 1 mm are widely used. The first works demonstrating the advantage of microchannel heat exchangers were published by Tuckerman [1] and Swift [2] in the early eighties and proved that the heat flux density significantly increases even in single-phase mode when switching to microchannels. The high heat transfer efficiency, in this case, is due to an increase in the ratio of the heat transfer area to the volume of the heat exchanger. As a result of further research, single-phase microthermal heat exchangers with a heat flux density of up to 300 W/cm 2 were developed [3][4].
In addition, recently there has been a significant interest in two-phase heat transfer in microchannels, where extremely thin liquid coolant films with huge values of the heat transfer coefficient could be realized. It is worth distinguishing separately the processes of evaporation and boiling of liquids in microchannels, at which the record-breaking values of the heat flux density up to 1500 W/cm 2 are reached [5]. Interest in these applications has led to the appearance of a huge number of works performed on this topic in recent years. An example of such a two-phase device is a new efficient cooling system (Fig. 1), in which heat is removed due to intense evaporation of a thin film of liquid moving in a flat microchannel under the action of a gas or vapour flow [6][7][8]. 2 A spatial mathematical model of this evaporative system has been developed in this work. The model fully reproduces the real spatial geometry of the measuring section and easily adapts to the parameters of the experimental setup.

The experimental setup
A detailed description of the experimental setup is given in [6][7][8]. Here we describe only its main features. Schematic of design of the test section of experimental setup is presented in Figure 1. The main part of the test section is a thin stainless steel plate with a flush-mounted copper rod. At the working surface the rod has a 1×1 cm square head emulating surface of a computer chip. The rod is heated by a heating spiral coiled around its bottom part. Such construction of the heater provides the condition at the heater surface T = const. Gas, pumped by a compressor, flows through the channel and passes to the atmosphere at the end of the test section. Liquid, supplied from a thermostat, gets into the channel through a liquid nozzle and flows under the friction of the gas along the stainless steel plate as a film.
The liquid is accumulated at the bottom part of the test section and is returned into the thermostat. The stainless steel plate is kept at constant temperature by a 3 mm diameter channel through which the working liquid is pumped. Distance from the gas inlet to the liquid nozzle is 57 mm while distance from the liquid nozzle to the heater is 32 mm. This provides steady flows of gas and liquid at the moment they reach the heater. Channel width is 30 mm. 6thermocouples; 7copper rod; 8nichrome spiral; 9stainless steel plate; 10outlet of fluid; 11outlet of gas; 12glass cover/

The mathematical model
A spatial mathematical model of the evaporation system using film flows has been developed. The model fully reproduces the real spatial geometry of the measuring section and easily adapts to the parameters of the experimental setup. The computational domain is the inner part of the measuring section, which includes sections for supplying liquid and gas, a flat channel for the flow of film and gas, a metal substrate, part of a copper heater and the output section. The geometry of the measuring section is symmetrical relative to the central longitudinal plane of the channel, therefore, to save the estimated time using that symmetry conditions, half of the calculation area is considered. A numerical technique based on the Volume of Fluid (VOF) method is used (Hirt C.W., J. Comput. Phys., 1982) to describe the conjugate heat transfer in the film flow regime in an evaporation cell. This model is based on our works [9][10]. To simulate surface tension in the framework of the VOF method, the CSF algorithm is used (Brackbill, et. Al. Comput. Phys., 1992); it involves introducing an additional force into the equations of motion. This force is determined at the interface through the gradients of the volume fraction of the liquid in the computational cell. In the general case, the model takes into account the thermocapillary convection of Marangoni. The surface tension coefficient linearly depends on the temperature on the film surface. In addition, the effect of free convection in liquid and gas is taken into account. In this case, to describe the dependence of gas density on the temperature the ideal gas model is used, and to describe the dependence of liquid density on the temperature the Boussinesq model is used. A simplified model proposed by Gupta et al (2010) is used to take into account the evaporation of the film. According to this model, an additional term is introduced into the right side of the liquid volume fraction transfer equation, which is proportional to the difference between the saturation temperature and the surface temperature on the film. The mass transfer coefficient is determined according to the model proposed by Schrage et al (1953). The corresponding source terms, which are proportional to the flow rate of the evaporated liquid, are introduced into the energy and momentum conservation equations. In the general case, the model takes into account the dynamic behaviour of the contact angle. To determine the dynamic contact angle, the model proposed by Kistler (Kistler, S.F., Hydrodynamics of wetting, Wettability, 1993.) and based on the use of the equilibrium value of the contact angle and capillary number, is selected. The gradient adaptation technology of the computational mesh is used to resolve the phase boundary during film motion. Using this technology, the calculation mesh is automatically thickened in the region of large solution gradients during the calculation process. In this case, the gradient of the volume fraction of the liquid is used as a control parameter. The mesh cells in the interface region can be fourfold to sixteenfold smaller than in the original mesh. This technology allows obtaining a solution, acceptable in terms of accuracy on the available calculation meshes. Both laminar and turbulent flow regimes are considered. To simulate turbulent flow regimes, the LES method is used. To determine the subgrid viscosity in the LES method, the WALE submodel is applied.
The difference analogue of convective-diffusion equations is found using the finite volume method for structured multi-block grids, the application of which automatically provides the conservative scheme. The connection between the velocity and pressure fields is realized using the SIMPLEC procedure on combined grids. An implicit second-order scheme is used to approximate the unsteady terms of the transport equations. A second-order central-difference scheme is used to approximate the convective terms of the hydrodynamics and energy conservation equations. A TVD scheme with an HRIC flow restrictor is used to solve the convective transport equation for the volume fraction of liquid in the cell.
A series of methodological calculations are carried out in order to determine the effect of the detailing of computational grids and other numerical parameters on the results of numerical modelling. As a result, the total number of cells of the computational grid for the complete model of the evaporation cell in the process of dynamic adaptation is about 28 million. The value of the time step is determined from the condition CFL < 2. Methodological calculations show that such detailing is sufficient to obtain a grid-independent solution.

The simulation results
The numerical simulation of several regimes of film flow in the evaporation cell at a heater temperature below the boiling point was carried out. The temperature of the heater was 40°C. The inlet temperatures of the liquid and gas were 25°C. The evaporation from the film surface was negligible under these conditions. An investigation of the flow regimes of the film was carried out. The Reynolds number of the fluid Rel (Rel=l/l, where l is the specific mass flow rate of the fluid, l is the dynamic viscosity of the fluid) in the calculations ranged from 5 to 200. The Reynolds number of the gas Reg (Reg=Qg/νg, where Qg is the specific volume gas flow rate, νg is the kinematic viscosity of the gas) ranged from 400 to 3000. Five different film flow regimes, which were previously discovered in experiments [6][7][8], corresponded to these Reynolds numbers. Thus, the following regimes have been considered: channel flooding, smooth film, 2D-waves, 3D-waves and film rupture regimes. As a result of numerical simulation, the flow regimes, which previously discovered in the experiment, have been fully confirmed.
The visualizations of two different modes of film and gas flow, observed in calculations and experiment, are presented in Figures 2-3 as an example. Such flow parameters correspond to the 2D-

wave (Figures a) and 3D-wave (Figures b) flow regimes. It is shown that the model in whole
reproduces well all the features of the behaviour of the film in these regimes. Good agreement has been obtained between the calculation and the experiment on the shape of the surface and the thickness of the film, as well as on the wavelengths and frequencies of the ridges. It should be noted that the symmetry conditions were used in the calculations and only half of the measuring section was considered. Analysis of the simulation results shows that such approximation worked well for low Reg numbers (less than 1000). For higher Reynolds numbers, the flow became substantially threedimensional, and deviations from longitudinal symmetry were observed (see Figure 3a), especially near the outlet section, therefore it is necessary to use the full formulation of the simulation. This will be done in the future. The temperature distribution on the film surface and in the central longitudinal flow of the substrate is shown in Figures 4. The calculation allows studying the local behaviour of the film temperature for this device for the first time. As you can see, there is an obvious correlation between the film thickness and the temperature distribution in it. The film temperature is strongly inhomogeneous in the twodimensional waves and transition regimes. Significant overheating of the liquid is observed in the troughs of the film. The local non-uniformity of the film temperature is smoothed out by turbulent pulsations upon transition to pronounced turbulent regimes of the film flow (see Figure 4b). In addition, using numerical simulation, heat leakages from the heater into the surrounding substrate and the environment are investigated. Figures 4 clearly show that the heating of the film begins earlier than it reaches the heater due to the thermal conductivity of the substrate. In addition, the influence of the substrate continues downstream.
In addition, the developed mathematical model allows studying the local distribution of the heat flux density both on the surface of the heater and its surrounding substrate for the first time. The distribution of heat flux density for the three-dimensional wave regime is shown in Figure 5. As one can see, a significant non-uniformity of the heat flow is observed. Local maxima are clearly visible, the values of which, in this case, are about 120 W/cm 2 against the background of an average value of about 30 W/cm 2 . An analysis of the film thickness shows that these heat flow maxima coincide with the troughs of the waves on the film surface. The heater is found to "feel" the relief of the film above it very well. At high values of overheating, these places will be centres of dry spots formation. In addition to the spatial non-uniformity of the heat flux, its significant non-stationarity is observed. The dependence of the average heat flux density over the surface of the heater on the time for two different regimes of film flow is shown in Figure 5b. As can be seen, there are significant pulsations of the heat flux. They are especially pronounced in the regime of two-dimensional waves. An analysis of the simulation results shows that the period of these pulsations correlates with the period of passage of large waves. The amplitude of the heat flux oscillations, in this case, can be 30-40% of the average value. This fact is important in the design of film-based cooling devices.

Conclusions
A spatial mathematical model of the considered evaporation system has been developed in this work. The model fully reproduces the real spatial geometry of the evaporation system and easily adapts to the parameters of the experimental setup. The model is based on the Volume of Fluid method and takes into account thermocapillary effects and the process of film evaporation, which, in its turn, takes into account the dynamic angle of wetting. A numerical study of several regimes of film flow in the 6 evaporation cell at a heater temperature below the boiling point has been carried out to test the model. As a result of numerical simulation, the flow regimes, previously discovered in the experiment, have been fully confirmed. The model is shown to well reproduce all the features of the film behaviour in these modes. Good agreement between the calculation and the experiment on the shape of the surface and the thickness of the film, as well as on the wavelengths and frequencies of the ridges has been obtained.