On modeling of the initial stage of nonstationary nucleate boiling for the high heat fluxes

During recent years, there have been made significant achievements in the numerical description of the bubble boiling, particularly, in the calculation of the bubble growth dynamics, the nucleation density, and the bubble boiling threshold [1, 2]. However, the numerical prediction is mostly based on the empirical correlations, the accuracy of which does not exceed 20%. Boiling is a complicated process, where each parameter affects not only the general outcome but other parameters, too. The non-stationary heat release is most difficult for modeling, because many of the existing researches are based on analytical expressions for the fluid temperature distribution. The basic stages of the explosive boiling are: (i) heating of the wall to the nucleation temperature; (ii) nucleation of the isolated bubbles; (iii) merging of bubbles and coverage of the entire surface by the steam phase; (iv) heating of the surface to the Leidenfrost temperature and transition to the boiling crisis. This paper presents results of the experimental study for the initial stages of the explosive boiling (i) and (ii), as well as an attempt to simulate them in order to clarify whether the existing approaches can be extended to the case of the nonstationary heat release.


Introduction
The present research considers the initial stage of the bubble boiling without bubble agglomerations. In the case of high heat fluxes the subcooled liquid is located at the small distance from the heater surface. Figure 1 shows the comparison of the experimental bubble diameter with the theoretical value of the overheated layer thickness, √ , where α is the thermal diffusivity and t is the time passed since the wall achieved the saturation temperature. Analysis of this figure shows that the bubble diameter can be estimated on the basis of the overheated layer thickness, taking into account the liquid bulk temperature, wall overheat, and the heat flux. However, it was revealed in [3] that no single approximation of this kind can cover all possible regimes of non-stationary heat release. We propose the following approach: the ratio between the overheated layer thickness and the bubble maximum diameter is assumed to be constant [2], and the nonstationary overheated layer thickness is obtained numerically by taking into account not only the conductive heat transfer but also the evaporative heat transfer. The coupling between the temperature field and the averaged bubble diameter in homogeneous flow is modeled with the help of COMSOL Multiphysics (ver. 5.3) software.
The calculation results are compared with the experimental data obtained in a vertical channel with a cylindrical heater (120 mm long, 12 mm in diameter and 1 mm thick). The stepwise heat flux was used with the maximum up to 1.0 MW/m 2 .

The heat balance based model for nucleate boiling
The expression for the maximum diameter, given in [5,6], is based on determining the heat balance for a bubble, part of which is submerged in a superheated liquid layer and another part is surrounded by the liquid subcooled to the saturation temperature: , k is the thermal conductivity. The factor m is the fraction of the bubble surface area covered by the subcooled liquid at the temperature T0, as proposed in [5]. Following [4], φ (characterizes the effect of liquid velocity) and C (a function of system pressure) for the present condition are equal to 1 and 65, correspondingly. The cooling area fraction m of 0.3 was indicated as the best one for approximating the experiment data in [4] and 0.5 for the data in [4]. By definition, m is related to the location of Ts isotherm as follows: where xs is the thickness of the superheated liquid layer. For the boiling stages of short duration, we can assume that the ratio between the thickness of the superheated layer and the maximum bubble diameter is a value of small variability. It means that during nonstationary boiling simulation we can calculate the temperature field and easily obtain the bubble diameter via Eq. (7). In [3] we employed the homogenous representation of the liquid wall layer without taking into consideration the influence made by the nucleate boiling on the temperature distribution. It was shown that in this case the value of m stayed within the range 0.5-0.9. In the present work we enhance the existing axisymmetric model of the nonstationary heat release into the boundary liquid layer with the dependencies for the main heat fluxes arising from bubble boiling.

Numerical simulation
To calculate the heat balance, it is required to know the description of the temperature field. According to the measured wall temperatures and bubble diameters, the temperature gradient near the heater surface exceeds 0.5K/µm. Such detail is possible only if we use numerical modeling. The water temperature profile, the heater surface temperature, and the radial heat flux were calculated with the help of the COMSOL software.
We used a 2D axisymmetric model (Fig.2), which comprised a 55 mm long steel heater with inner diameter of 10 mm and outer diameter of 12 mm, the surrounding water layer was 3 mm thick. To obtain a finer resolution for the superheated layer, we applied a mesh of about 180000 elements with 2 boundary layers near heaterwater interface: the first layer is 0.05 mm thick, and has the step of 0.001 mm, the second layer is 0.27 mm thick, and has the step of 0.005 mm (Fig. 3). The remaining part of the model is meshed with the step of 0.1 mm. In order to obtain an overheated layer with a sufficient resolution, we used a mesh with a total number of elements of about 180,000 and two boundary layers near the interface between the heater and water regions: the first 0.05 mm thick layer has a mesh step of 0.001 mm and the second 0.27 mm thick layer has a mesh step of 0.005 mm. The rest of the model is covered with a 0.1 mm mesh.
The water flows into the channel with uniform velocity of uin=0.2 m/s and flows out at atmospheric pressure. No-slip boundary conditions were set at the outer wall and the heater surface. To obtain the temperature field relevant to experimental bubbles data, the resistive heat release inside the heater was calculated by solving the problem of electrical conductivity with the boundary conditions based on the voltage history of each individual experiment. The bulk temperature was specified at the water inlet and at the outer wall. The zero normal heat flux condition was used at the water outlet.
The laminar fluid flow problem coupled with the heat transfer in fluid problem was coupled with the direct electric current problem The changes in the temperature profile that took place beyond the first 50 mm of the heater length and the first 3 mm of the heater thickness, as well as the gravitation effect, were checked numerically and found negligible. The heater surface temperature measured by thermocouples showed good agreement (<5%) with numerical modeling until the moment of the fully developed nucleate boiling. In addition to the heat fluxes described by the heat transfer equations (10-11), this paper also introduces evaporative heat transfers associated with the life cycle of bubbles (Fig. 4). In this case, the bubbles, as well as the vapor cavities of a specific localization, are not considered, and the sequence of events of the bubble life cycle is not modeled. Heat fluxes are averaged over the bubble lifetime , 1 ≈ 0.2 , and over the area of its influence of size ± in the plane of the heater surface. Since the axisymmetric model cannot be used for representation of an individual bubble, we model a solid ring of bubbles at the height of 45 mm from the beginning of the channel. Modeling the ring of bubbles, rather than a continuous coating of the heater surface, makes it possible to estimate the effect of convection and thermal conductivity along the channel on the temperature field, and to determine the dimensions of the long-term area of influence over the entire period of the nucleate boiling, which lasted about 2.5 ms. The limitations of the axisymmetric model do not allow for the distribution of heat in the angular direction. Due to this fact, in this work we obtained a high estimate of the influence of evaporative heat transfer on the change in the thickness of the superheated layer and the coefficient m.
In addition to estimating the local heat fluxes in the area of bubble influence, it would be interesting to measure the contribution of all bubbles to the total heat transfer from the heater surface. To do this, we need to know the nucleation centers density. When processing video data, for each frame we calculated the corresponding time from the start of the heating and wrote out the coordinates and sizes of all the bubbles in the frame. Based on the scale of the frames (0.01 mm/pixel along the x axis and 0.02 mm/pixel in the y axis) and the angle of inclination between the heater surface and the filming direction, we obtained the area of the frame (8.5 mm 2 ) and the nucleation centers density, n, in each frame. The obtained values are approximated by the equation = 0.243 * exp (( − 6.14)/0.77) where the time from the start of heating is given in milliseconds. The nucleation centers density at the end of the nucleate boiling period reaches about 5 mm 2 . Therefore, the distance between the centers of nucleation is about 0.44 mm. With a maximum observed bubble diameter of this period up to 0.2 mm, the distance between the centers of the bubbles is approximately 2.2 diameters. When studying the sequence of frames, it can be seen that  Steel q-the bubbles appear in the same placesnucleation centers, the number of which becomes bigger as the temperature of the heater increases. Throughout their life, bubbles are spherical and touch the surface of the heater. They do not slide along the surface, but collapse, or decrease, remaining at the nucleation point. For each nucleation center, a sequence of frames corresponding to the local maxima of the bubble diameters was chosen. The volume-averaged maximum diameters, for all the nucleation centers in the field of view, were approximated by the formula = 0.45 − 0.17 * + 0.017 * 2 (16) where the time from the start of heating is given in milliseconds, whereas the maximum bubble diameter is measured in millimeters. For each measurement of the local maximum of the bubble diameter, its lifetime was determined as the period of time elapsed since the previous measurement of the local maximum of the diameter in the same nucleation center. We also measured the inverse of the lifetime valuethe nucleation frequency. The average nucleation frequency was = 4.67 kHz. In the model, the dependency of the nucleation frequency on the start of the heating was represented as Here the frequency is expressed in kHz, and the time from the start of heatingin milliseconds. The approximations presented here for the nucleation centers, the nucleation frequency, and the bubble diameter, as well as the results of calculations of the temperature field with allowance for evaporative heat transfer, were used in the same experiment with an initial temperature of = 90 ° and an average growth rate of the heater surface temperature of 5.6 K/ms. The size of the influence area of an individual bubble can be estimated based on its lifetime. Due to forced convection with a maximum flow rate of 0.2 m/s, the temperature distortions associated with the formation and collapse of bubbles move 40 mkm. However, it should be taken into account that at the distance of one bubble diameter from the wall, the water flow velocity does not exceed 0.05 m / s. In this case, the estimation of the influence distance due to forced convection is reduced to 10 μm. The distance of influence due to thermal conductivity, calculated by the formula √ , amounts to about 5 μm for water at the saturation temperature of 102 ° C and to δ-=30 μm for steel. Therefore, the maximum diameter of the influence area of a bubble during its lifetime is stipulated by the thermal conductivity of steel, and is equal to D-= D+=0.23 mm, which corresponds to the influence area − = + = ± 2 4 = 4.15 • 10 −8 2 .
Judging by the size of the temperature perturbations, which are visible on the frames following the bubble collapse, the thickness of the influence area is about one and a half maximum diameters, + ( ) = 3 2 ( ).
The heat flux per unit surface of the heater is considered to be composed of four components: the bubble and the microlayer evaporation heat, both averaged over the lifetime and the influence area; the thermal conductivity through the microlayer; the thermal conductivity through the superheated layer.
This is the local heat flux in the influence area of a bubble. When estimating the contribution of bubbles to the overall heat transfer, it is necessary to use data on the density of nucleation centers: The heat transfer through the dry spot and the rewetting effect are neglected, since it is known that in this series of experiments bubbles were spherical and the diameter of a dry spot was much smaller than that of a bubble. The heat flux of the initial evaporation, , , is determined as the evaporation heat released during the bubble growth, averaged over the influence area and the bubble lifetime The absorption of heat by a bubble due to evaporation of the microlayer is determined by the volume of the microlayer, The heat flux due to the thermal conductivity of the microlayer is assumed to be constant within the radius of the microlayer 2 from the center of a bubble, and is averaged over the whole bubble influence area.
In accordance with Fig. 4, the radius of the microlayer is defined as follows: As shown in some papers, including [6][7][8], the microlayer thickness reaches its maximum value at the stage of bubble growth. Further, the area and thickness of the microlayer decrease due to its evaporation. Therefore, the calculations for a constant-sized microlayer that were carried out in this work are the estimates aimed at evaluation of the relative contribution of each component of the heat flux. The microlayer thickness, similar to the formula obtained in [8], is determined by the expression Here, the bubble growth time is estimated as 1/3 of the nucleation period, and the coefficient 0.25 is chosen so that the thickness of the microlayer does not exceed 1.5 μm.
The heat flux through the superheated layer is considered proportional to the wall overheating ( − ), but the layer thickness varies from to . This heat flux is integrated over the heater area under the corresponding part of a bubble and is averaged over the entire influence area.
The obtained value ′′ is the heat flux, additionally transferred from the heater into the liquid. Heat transfer is realized in the form of a volumetric heat source, ′′′ . With respect to z, this source is determined in the range 45 − ± 2 < < 45 + ± 2 . In the thickness of the heater, the heat is absorbed to the depth − . In a liquid, heat is released linearly at the distance up to δ_ + = 1.5 Dm and reaches its maximum at the center of a bubble. The heat release integral ′′′ is equal to zero over all r .

Conclusions
Taking into account the contribution of interphase transitions to heat exchange, the change in the thickness of the superheated layer is proportional to the change in the maximum bubble diameters. The fraction of the bubble protruding into the region of the subcooled liquid, as a result of the introduction of heat fluxes associated with the bubble into the model, decreased by no more than 0.1, remaining substantially higher than the values indicated by other authors (Fig. 6).