Numerical simulation of heat transfer in falling wavy liquid films on unsteady heat release

The mathematical model which allows the calculation of the wave surface profile as well as velocity and temperature fields has been presented. The numerical simulation of heat transfer in falling wavy films of liquid nitrogen has been performed. The dependencies of boiling expectation time and total local evaporation time on heat flux density for different inlet Reynolds numbers have been calculated. The regime map which describes the different mechanisms of film decay was obtained by summing up the simulation results. The results of numerical simulation are in satisfactory agreement with the experimental data.


Introduction
The liquid flow over a vertical surface is inevitably accompanied by development of wave perturbations on the free film surface which bring considerable contribution to the heat transfer processes and development of crisis phenomena [1][2][3].
The area of highly efficient heat transfer at the film flow of liquid is limited by the heat flux and depends significantly on the heat release law. When the system reaches the critical heat flux, this leads to decay of the film flow, formation of large-scale dry zones, reduction of heat transfer intensity and as a consequence to a drastic increase in temperature of the heat-releasing surface. All these can cause destruction of the heat-releasing element. The ability to calculate the critical heat fluxes and maximal times of their impact for a given system is required for the design of stably working heat exchangers.
The time of crisis beginning and scenario of film flow decay are characterized by such parameters as time of boiling expectation and time of total local evaporation of the falling liquid film. These values determine maximal times of heat flux influence on the film with certain parameters.
The mathematical model which allows the calculation of the wave surface profile as well as velocity and temperature fields has been presented in the current work. This model allows calculations of the wave characteristics and parameters of liquid film decay at different regime parameters. The results of numerical simulation performed in the framework of the suggested model have been compared with the experimental data obtained for the liquid nitrogen in [2].

Mathematical model
Simulation has been carried out within a 2D calculation area. The abscissa axis is directed along a vertical wall, downstream the liquid film flow. The ordinate axis is orthogonal to the abscissa axis and directed from the wall to the free film surface. The detailed description of the calculation area and the mathematical model has been presented in [4].

Hydrodynamic model
To describe the processes of liquid flow and evolution of the wave surface profile, in the framework of the current study the 2D hydrodynamic model of Kapitza-Shkadov has been used: (4) The method of liquid flow description by means of flow rate characteristic Q allows excluding the dependence on the ordinate and reducing the problem to the pseudo-singledimensional one. Velocity components in the whole calculation area can be obtained from h and Q by means of expressions (3) and (4). The system of equations (1 4) supplemented with the initial and boundary conditions allows calculation of flow dynamics and evolution of surface profile within the 2D calculation area.
Perturbations at the inlet of the calculation area were set through the addition to the flow rate of the undisturbed film under the condition of keeping the thickness undisturbed: Values h 0 and Q 0 set by expressions (5) and (6) are also the initial conditions. In other words, at zero time moment the film is assumed to be waveless and its thickness and flow rate are constant and determined by the conditions of the problem. The additive member in (7) consists of two multipliers: A is a random real number from 1 to 1 and R is the relative perturbations amplitude.

Heat transfer model
In natural variables the boundaries of the calculation area are not constant because of the unsteady character of the free surface. The following coordinate transformation converting the calculation area into the rectangular one has been used to solve this problem: (8) After the coordinates transformation the unsteady heat conductivity equation in the liquid film looks like the following: At the inlet of the calculation area and at the free film surface, the temperature was assumed equal to the temperature of liquid saturation. The following boundary condition was used on the heating surface: (10) To take into account evaporation from the free surface, equation (2) has been modified as follows:

Results and conclusions
The boiling expectation time was determined as the minimum time of reaching the given overheat by the heating surface. In the framework of the present work three threshold values of overheating for the liquid nitrogen under conditions of atmospheric pressure were used: two theoretically calculated values, ΔT b.e = 26 K, corresponding to 95% and ΔT b.e = 32 K, corresponding to 100% of probability of boiling under the conditions of homogeneous nucleation [5], and one approximation for experimentally retrieved dependence of overheat on heat flux density obtained for boiling-up under the condition of heterogeneous nucleation: The calculated boiling expectation times (Fig. 1) are in satisfactory agreement with the experimental data [2].
The regime map which describes the different mechanisms of film decay was obtained by summing up the results of numerical simulation (Fig. 2). Calculated curves 1, 2 and 3 divide the area of the investigated Reynolds numbers and heat flux densities into the set of regions. Region I is characterized by the absence of crisis; region II describes the set of the problem parameters, where liquid film decay will be performed by total local evaporation without liquid boiling-up; region III describes the set of parameters, where the liquid film decay will be performed by the film boiling-up. Calculated critical heat flux densities of the total local evaporation 3 agree satisfactorily with the values calculated for steady heat release conditions by the formula, which has been presented in [6] -curve 4.
According to the presented regime map, it is possible to make a conclusion that boilingup becomes the prevailing decay mechanism independent of the heat flux density when the inlet Reynolds number Re in exceeds 130 -220 depending on liquid and heater parameters.  The study was performed at Kutateladze Institute of Thermophysics SB RAS and was financially supported by the grant of the Russian Science Foundation (project No. 14-49-00010).