Abstract
A direct numerical simulation of the boiling phenomena is one of the promising approaches in order to clarify their heat transfer characteristics and discuss the mechanism. During these decades, many DNS procedures have been developed according to the recent high performance computers and computational technologies. In this paper, the state of the art of direct numerical simulation of the pool boiling phenomena during mostly two decades is briefly summarized at first, and then the nonempirical boiling and condensation model proposed by the authors is introduced into the MARS (MultiInterface Advection and Reconstruction Solver developed by the authors). On the other hand, in order to clarify the boiling bubble behaviors under the subcooled conditions, the subcooled pool boiling experiments are also performed by using a high speed and high spatial resolution camera with a highly magnified telescope. Resulting from the numerical simulations of the subcooled pool boiling phenomena, the numerical results obtained by the MARS are validated by being compared to the experimental ones and the existing analytical solutions. The numerical results regarding the time evolution of the boiling bubble departure process under the subcooled conditions show a very good agreement with the experimental results. In conclusion, it can be said that the proposed nonempirical boiling and condensation model combined with the MARS has been validated.
1. Introduction
Despite extensive research efforts regarding the boiling heat transfer, the mechanism of the nucleate boiling phenomena is still unclear, and therefore its mechanistic model without any empirical correlation has not been developed. Since the tempospatial scales of boiling phenomena vary from the molecular scale motion in the heat conduction to the convective motion governed by the macroscopic scale in the bulk fluid, the scales of the phenomena of interest should be clear when a mechanistic boiling model is developed. Especially, the molecular motion is very fast in the nucleation process, thus, the nonthermal equilibrium state must play a very important role in this time scale.
On the other hand, a direct numerical simulation (DNS) of the boiling phenomena is one of the promising approaches in order to clarify the heat transfer characteristics of boiling phenomena and discuss their mechanism. With the advances in recent years of high performance computational technology, several DNS procedures to solve the conservation equations of mass, momentum, and energy for liquid and vapor phases simultaneously when an interface is continuously evolving at and near a heated surface have been developed: to see the recent review [1] by one of the authors.
This study focuses mostly on the thermal equilibrium processes of boiling phenomena, which means that the nucleation theory can be applied to determine the size of vapor embryo, and then its growth, departure, and condensation processes will be considered later. In this paper, a brief review of the recent DNS on the boiling phenomena is firstly summarized, and then in order to validate the DNS data, the results [2–6] obtained by the DNS based on the MARS which is the interface volume tracking method [7] are compared to the visualization data obtained by using a very high speed video camera and a Cassegrain telescope [8].
2. Direct Numerical Simulation of Nucleate Bubble Growth, Departure, and Condensation
With the advances in recent years of high performance computational technology, several DNS procedures to solve the conservation equations of mass, momentum, and energy for liquid and vapor phases simultaneously when an interface is continuously evolving at and near a heated surface have been developed [1]. The solutions provide not only the detailed physics associated with thermal and hydrodynamic processes, but also the evolution of the interface shape. Most of the numerical investigations focus on either the nucleate boiling at moderate heat flux (isolated bubbles) or the stable film boiling. One of the key problems is the numerical description of two-phase flow. Typically, there are several methods which are used for the simulation of boiling heat transfer as follows.
Marker and Cell (MAC) Method [9]. The interface is marked by many massless particles that are convectively transported with the Eulerian velocity field and can be used to reconstruct the interface position on a fixed mesh. As far as the author knows, the first boiling simulation based on MAC method was done by Madhaven and Mesler in 1970 [10].
Front Tracking (FT) Method [11, 12]. This method seems to be an advanced version of MAC method and is thus based on massless particles that follow the motion of the interface driven by the Eulerian velocity fields. The interfaces are tracked by a Lagrangian frame, and the interface position can be tracked precisely, but it is necessary to introduce some artificial rules to treat the break-up interface and also coalescence of interfaces during computation. Esmaeeli and Tryggvason used the model mainly for the simulation of film boiling [13].
Arbitrary Lagrangian-Eulerian (ALE) Method [14–17]. This is based on a dynamic mesh that follows the motion of the interface. Thus, the interface coincides with a boundary of the computational domain at all times. According to this feature, the center of the control volume may not coincide often with the center of mass, and eventually the introduction of large artificial viscosity will be inevitable. An important advantage of the ALE method and boundary-fitted meshes in general is the possibility of treating the liquid-vapor interface as a boundary of the computational domain. This facilitates the calculation of the heat flux at the interface and therefore of the evaporation rate.
Level-Set (LS) Method [18, 19]. This method uses a field that contains information about the distance of a numerical cell to the interface and which is convectively transported with the Eulerian velocity fields. The interface is represented by the zero-contour line of the level-set field. The criticism regarding LS method is a transport feature of the distance function because the distance is not a measure of physical property. Recent works were reviewed in [1].
Volume of Fluid (VOF) Method [4, 20–29]. The most popular volume capturing procedure is based on the transport equation of VOF (volume of fluid). This uses a field that contains information about the volume fraction of one of the phases in a numerical cell and which is convectively transported with the Eulerian velocity fields. The volume fraction field has a piece-wise value at the position of the interface. There are some problems with VOF method, such as the interface reconstruction between neighborhood cells and the precise calculation of interface curvature. Some recent works on boiling simulation based on VOF methods were listed in [1]. However, none of the aforementioned models based on VOF method include any submodel for evaporation at the 3-phase contact line.
2.1. Governing Equations Based on MARS
The governing equations of the MARS consisted of the continuity equation for multiphase flows, the momentum equation based on a one-fluid model, and the energy equation with an external work done by a phase change phenomenon as follows: where is VOF fraction, the suffix denotes the th fluid or phase, is time, is pressure, is temperature, is gravity, is viscous shear stress, is body force due to a surface tension based on the continuum surface force (CSF) model [30], is density, is specific heat at constant volume, is thermal conductivity, is heat source, and a bracket <> denotes an average of thermal properties. In order to satisfy the conservation of , the third term of the continuity equation (1) must be included. The second term of the right hand side of the energy equation (3), that is, the Clausius-Clapeyron relation, was considered as the external work done by the phase change, for example, a bubble oscillation caused by the expansion and contraction with the bubble growth and condensation processes. The interface volume tracking technique [7] was applied to the continuity equation in the MARS. The projection method [31] was applied to solve the momentum equation and the pressure Poisson equation was solved by the Bi-CGSTAB [32].
2.2. Nonempirical Boiling and Condensation Model
The boiling and condensation model in the MARS for the subcooled nucleate boiling phenomena consisted of both a nucleation model and a bubble growth-condensation model based on the temperature-recovery method [2, 3]. This model was applied only to the interfacial cells which have the VOF fraction between 0 and 1. A density-change between water and vapour was considered as a volume-change by a phase change ratio expressed as where is specific heat at constant pressure, is degree of superheat or subcooling (), is latent heat, and the suffixes of and denote gas and liquid phases, respectively. Equation (4) represents the ratio of the sensible heat to the latent heat at the interfacial cell. To conserve the total volume, the following constraint condition must be satisfied:
The improved enthalpy method consists of (4)-(5) based on the assumptions of a zero-thickness interface and means a “rapid” phase change from “State 1 (Water)” to “State 2 (Vapour)” or vice versa due to the local quasi-thermal equilibrium hypothesis. This hypothesis requests another contradictory assumption that is a “very slow” phase change from “State 1” to “State 2.” In order to satisfy this requirement, a finite thickness of interface must be considered and both from/to “very slow” to/from “rapid” changes will occur simultaneously in the phase change process. This process can be described by a relaxation or waiting time for consuming the latent heat at the interface region, that is, the unsteady heat conduction in the finite interface region as the “very slow” change process: the relaxation time can be defined as a time when the phase change front is passing through a fictitious interface thickness , so that can be expressed by using the thermal diffusivity of the water as follows:
On the other hand, a well-known thermal penetration length for the unsteady heat conduction in a semiinfinite slab with a constant boundary temperature was approximated by the following expression:
If the thermal penetration depth can be assumed the same as the fictitious interface thickness, can substitute into (7): . As a result, an invariant relation between the thermal penetration length and the fictitious interface thickness can be obtained as follows:
According to this relationship, the fictitious interface thickness in around 30% of the thermal penetration depth corresponds to the “very slow” phase change process. In other words, the rapid phase changed volume during is 70% of the thermal penetration depth, not 100%. In this study, the relaxation time can be considered to control the VOF function as a phase change limiter. For example, the relaxation time for both phase-fronts is assumed to be 15% at either edge of the interface such as evaporation or condensation front as follows:
This proposed model has been named the “nonempirical boiling and condensation model” by the authors [5].
3. Visualization Experiments
3.1. Experimental Setup
In order to verify the boiling and condensation model, the visualization experiments of subcooled pool boiling by using ultrahigh speed video camera with a highly magnified telescope system have been conducted [8]. Figure 1 shows the schematic of experimental apparatus. Test section was a rectangular section of which wall was made of polycarbonate. The size of the test chamber was 150 mm (Length) mm (Width) mm (Height). Test section was filled with the purified water. Water in a pool was open to the air, so the system pressure was an atmospheric pressure. The distillated water was used as a working fluid and the degassing was done by using the subheater. The water temperature, that is, the degree of subcooling, was controlled by the subheater located at the bottom of the pool. Test heater was a platinum wire of 0.1 mm in diameter located on the flat polycarbonate plate to prevent the natural convection. Direct electric current was supplied to the platinum, and boiling bubbles were generated somewhere on the wire. The degree of subcooling was set from 0 to , and heat flux was also set to 0.25 W/mm2.
The bubble behavior was visualized by using an ultrahigh speed video camera (Phantom 7.1, Vision Research Co.) with the Cassegrain optical system (Seika Co.) which works as the highly magnified telescope. The flame rate of recording was 10,000–60,000 frames per ms.
3.2. Experimental Results
Figure 2 shows the time variation of bubble volume in subcooled pool boiling at various degrees of subcooling by using simple image analyses [6, 8]. The bubble volume was estimated based on an assumption that the bubble was a spheroid. The thickness of the bubble edgeline was , that is, 6–10 per pixel, so the error of the bubble volume based on the above spatial resolution was ±2–13%. Every generated bubble showed almost the same behavior: the bubble size, the deformation history, the interval of bubble generation, and so forth were almost the same and were periodically observed. The open symbols depict the experimental data before the bubble departure, and the solid symbols depict the data after the bubble departure. In the case of , that is, almost the saturated boiling, the bubble volume is monotonically increasing, and then the bubble departs from the heating surface. With increase of the degree of subcooling, especially, , the bubble growth rate is decreased and the departure time of bubble is also reduced. It seems that the bubble departure is accelerated.
In order to estimate the time variation of bubble shape, Figures 3(a), 3(b), and 3(c) show the bubble height H, width W, and the aspect ratio H/W at various degrees of subcooling, respectively. The open and solid symbols were defined the same as in Figure 2. In Figure 3(a), the height of the bubble is monotonically increasing and the growth rate of the bubble height is almost independent on the degree of subcooling until the bubble departure from the heating surface. In contrast, the width of the bubble decreases with the increase of the degree of subcooling as shown in Figure 3(b). In Figure 3(c), the aspect ratio of the bubble increases with the increase of the degree of subcooling; that is, the bubble shape becomes more vertically elongated just before the bubble departure from the heating surface. A part of this phenomenon might be related to the condensation because the subcooled water entrains to the root of growing bubble due to the inertia force caused by the surface tension.
(a)
(b)
(c)
From the results of the visualization experiments, it is found that the processes of bubble departure from the heating surface in subcooled pool consisted of the following processes:(1)the bubble is generated and is growing with the horizontally elongated shape in subcooled pool, then,(2)the bubble becomes of vertically elongated shape because of the recoil force due to the surface tension, and finally,(3)the bubble is departed from the heating surface due to both the inertia force and the condensation due to the entrainment of subcooled water.
4. Comparison of Numerical Results to Visualization Results
4.1. Computational Domain
In order to validate the numerical simulation, the bubble departing behaviors obtained in Section 3 were compared to the 3D numerical simulations based on the MARS coupled with the nonempirical boiling and condensation model described in Section 2.2. Figure 4 shows the computational domain. The computational domain size was . The computational cell size was 20–100 in and directions, respectively, and 20 in direction. The number of computational cells was , and the time increment in the computation was set to 1. The periodic boundary conditions were imposed at the and directions, respectively. The nonslip velocity boundary conditions were applied to all the walls, and the upper boundary condition in direction was set to a constant hydraulic pressure condition. The computational conditions were basically the same as the experiments. The initial pressure was set to an atmospheric pressure and the degrees of subcooling in the water pool were set to 5.1, 10.3, 15.4, and , the same as the experimental conditions. The heating platinum wire of 0.1 mm in diameter used in the experiment was located at the bottom of computational domain, and the volumetric heat source at the center of the wire was set to 0.2 W which was corresponding to the heat flux on the heating surface in the experiment: 0.25 MW/m2. The initial temperature of heating wire was set to which was estimated by using a waiting time obtained from both the experimental results and the analytical solution of unsteady heat conduction equation for the thermal boundary layer on the heating wall. The solid heat conduction from the center of the heating wire to its surface was numerically considered. In order to evaluate the effect of the wettability on the heated wire surface, the static equilibrium contact angle between the wire surface and the liquid was set to 20, 45, 60, and .
Since this study focuses on the bubble departing behavior from the heated wire surface, the initial bubble shape was assumed at the maximum bubble size in a horizontal direction obtained from the experimental data. The initial bubble diameter was thus set to the values as shown in Table 1, and the initial bubble as a hemispherical shape was put at the center of the heating surface as shown in Figure 1. The initial bubble temperature was also set to the saturated temperature under the atmospheric pressure. Although the velocity and temperature fields around the growing bubble were existed in the experiments, these initial fields in this numerical simulation were assumed to be both stationary and homogenous temperature fields with the degree of subcooling, respectively, because the experimental data corresponding to these fields could not be obtained.
The grid dependency of the computation is very important to calculate the curvature of the bubble which strongly depends on the number of grids. In general, a sphere can be approximated by 24 subdivisions, so the calculation of bubble curvature must be needed more than 24 grids. Since the minimum initial bubble diameter as shown in Table 1 was () = (0.88 mm, 0.70 mm) in case of , the ratio of the bubble diameter to the grid size corresponded to () = (44, 35). All the ratios of the other cases were more than 24. This means that the number of grids for the bubble in this study could be sufficient for the bubble growth process, because the bubble size increased with the increase of time. However, the condensation process was different, because the bubble size decreased with the increase of time. The grid dependency in the condensation process will be discussed in Section 4.4.
4.2. Bubble Departure Behavior
Figures 5(a) and 5(b) show the time variation of bubble departing behavior observed in the experiments (upper) and the numerical simulation (lower) at and , respectively.
The results of numerical simulation show the bubble shapes as an isosurface of VOF = 0.5, the temperature contours, and the velocity vectors. At (Figure 5(a)), the numerical results retrieved the experimental ones that the bubble changed from the flattened shape in the superheated layer to the vertically elongated one in the saturated or subcooled liquid layer before the bubble departure from the heated wire surface. At (Figure 5(b)), the bubble shape became more vertically elongated with the increase of the degree of subcooling. From the results of numerical simulation at , it is found that a large upward velocity like a jet from the bottom of bubble to the top appears when the bubble becomes the vertically elongated shape.
In order to validate the time variation of bubble shape changes obtained from the numerical results, Figure 6 shows a quantitative comparison of the time variation of the bubble aspect ratio (H/W= bubble height/bubble width) between the experimental data and the numerical simulation results at various degrees of subcooling. The open symbols depict the experimental data before the bubble departure, and the solid symbols depict the data after the bubble departure. The solid lines depict the numerical results before the bubble departure. As the results show, the time variations of the bubble aspect ratio obtained from the numerical results were in good agreement with the experimental data: with the increase of the degree of subcooling, the bubble became more vertically elongated shape, and the time interval from the bubble nucleation to its departure from the heated surface is decreasing.
However, it can be seen that the numerical results over predicted the time interval from the bubble nucleation to its departure from the heated surface, compared with the experimental data. This reason might be considered so that the initial condition of numerical simulations could not be completely consistent with the experimental conditions. Especially, it seems that the boiling bubble behavior is sensitive to its initial temperature field. Therefore, the effect of initial temperature field for the bubble departure behavior was investigated. In this study, the temperature field after one bubble departure was used as initial condition in the simulation; that is, the initial temperature field with the bubble departure history, was considered.
Figures 7 and 8 show the comparison of the bubble shapes just before the bubble departure from the heated surface and the bubble aspect ratio between experiments and numerical results using the initial temperature field with or without bubble departure history at various degrees of subcooling, respectively. As the results show, both the bubble shape and bubble aspect ratio predicted by numerical simulations with bubble departure history rather than without bubble departure history are shown to be in good agreement with the experiments. Consequently, we can say that the present numerical simulation based on the MARS with the nonempirical boiling and condensation model can predict the bubble departure from the heated surface in subcooled pool boiling behaviors as experimentally observed.
4.3. Effects of Wettability on Departure Behavior
According to the previous study [33], the critical heat flux is increasing with the improvement of surface wettability. It means that the departure of bubble from the heating surface could be accelerated by the better surface wettability. Figure 9 shows the bubble departing behavior with the various static equilibrium contact angles of 20, 45, 60, and at [5]. As a result, it is found that the bubble departure from the heating surface becomes difficult with the deterioration of surface wettability. Therefore, it can be said that the present numerical results agree with the previous study [33].
Figure 10 shows a quantitative comparison of the time variation of the bubble aspect ratio (H/W= bubble height/bubble width) between the experimental data [8] at , 15.4, and C, respectively [5]. The open symbols depict the experimental data before the bubble departure, and the solid symbols depict the data after the bubble departure. The solid lines depict the numerical results before the bubble departure. According to the results, it is found that the bubble aspect ratio of did not show the vertically elongated shape regardless of the degrees of subcooling and that of was shown to be in good agreement with the experimental data except in the case of . According to these results, it can be seen that the time variation of aspect ratio at was overpredicted, compared to the experimental data. This reason might be considered so that the initial condition of numerical simulations was not completely consistent with the experimental conditions: temperature and velocity fields. Especially, it seems that the boiling bubble behaviour of high subcooling condition is sensitive to its initial temperature field. Therefore, the temperature field experiencing one bubble departing history was used as the initial condition for . The dotted line in Figure 10 denotes the bubble aspect ratio of using the initial temperature field with one bubble departing history at . As a result, the time variations of bubble aspect ratio with one bubble departing history rather than without bubble departing history were shown to be in good agreement with the experiments.
4.4. Bubble Condensation Behaviors
Figure 11 shows the computational domain and the boundary conditions for the condensation process. The domain size was , and the uniform computational grids of 50 in , , and directions were used. The radius of bubble departure at various degrees of subcooling obtained from the experiment as shown in Table 2 was used as an initial vapor bubble, and then it was put as a sphere bubble in the grid nearest the wall. The ratios of the bubble diameter to the grid size were 22 for , 14 for , 12 for , and 8 for . The initial temperature field calculated until the temperature at the heated surface became larger than . The time increment in the computation was set to 1. Other conditions were the same as the bubble growth process simulation.
Figure 12 shows the comparison between experimental and computational results at . Despite the usage of rather small number of grids compared to the growth process, the computational results showed good agreement with the experimental data. On the other hand, in case of , the computational results by using 50 grids showed a slower condensation compared to the experimental results as shown in Figure 13. The 50 grid did not show enough spatial resolution. Thus, the smaller grid sizes as the grid refinement test by using 25 um and 10 um applied to the same condition were examined. As a result, the results by using 10 grid showed fairly good agreement in the early stage of the condensation process with high subcooling . The 10 grid size corresponds to the ratio of 40, so the grid ratio could be sufficient.
5. Conclusions
The numerical simulations based on the MARS with the nonempirical boiling and condensation model considering the relaxation time based on the local quasi-thermal equilibrium hypothesis were conducted for the bubble departing behavior from the heated wire surface in the subcooled pool boiling. The results of numerical simulations were compared with the experimental results, especially for the bubble shapes and the aspect ratios. As the results show, the numerical results for the bubble departing behavior were shown to be in very good agreement with the experimental results as follows.(1)The time interval from the bubble nucleation to its departing from the heated surface decreased with the increase of the degree of subcooling.(2)The aspect ratio of the bubble shape increased with the increase of the degree of subcooling; that is, the bubble shape became more vertically elongated before the bubble departure from the heated surface.(3)The initial temperature field with the bubble departure history has a profound effect on bubble departure behavior.(4)The departure of bubble from the heating surface is accelerated by the better surface wettability. In the case of , that is, the poor surface wettability, the bubble could not depart from the heating surface. On the other hand, the time variations of bubble aspect ratio in case of predicted by numerical simulation were shown to be in good agreement with experimental data.(5)The initial temperature field experiencing one bubble departing history has a profound effect on the bubble departing behaviour.
Resulting from the numerical simulations on the subcooled pool boiling phenomena, the numerical results obtained by the MARS are validated by comparing to the experimental ones and the analytical solutions. The numerical results regarding the time evolution of the boiling bubble growth/departure processes under the subcooled conditions were shown to be in a very good agreement with the experimental results. In conclusion, it can be said that the proposed nonempirical boiling and condensation model combined with the MARS has been almost validated. The remaining issue needs to consider a microlayer model [34] underneath the nucleate bubble close to the saturated condition. However, according to the present DNS results, the influence of the microlayer on the heat transfer might be limited in the highly subcooled conditions.
Nomenclature
| : | Specific heat at pressure (J/(kg·K)) |
| : | Specific heat at constant volume (J/(kg·K)) |
| : | Gravity (m/s2) |
| : | Volume of fluid (VOF) fraction (–) |
| : | Body force due to an interface tension (N/m3) |
| : | Bubble height (m) |
| : | Latent heat (J/kg) |
| : | Pressure (Pa) |
| : | Heat source (W/m3) |
| : | Wall heat flux (W/m2) |
| : | Bubble semimajor axis (m) |
| : | Bubble semiminor axis (m) |
| : | Time (s) |
| : | Temperature () |
| : | Heating surface temperature () |
| : | Velocity (m/s) |
| : | Bubble width (m) |
| : | , , and direction. |
| : | Thermal diffusivity (m2/s) |
| : | Thermal penetration length (m) |
| : | Fictitious interface thickness (m) |
| : | Phase change ratio (–) |
| : | Degree of subcooling () |
| : | Relaxation time (s) |
| : | Distance to the heating surface (m) |
| : | Thermal conductivity (W/(m·K)) |
| : | Static equilibrium contact angle (°) |
| : | Density (kg/m3) |
| : | Interface tension coefficient (N/m) |
| : | Viscous share stress (N/m2). |
| : | Gas phase |
| : | Liquid phase. |
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.