#### Abstract

Physical and mathematical model has been developed to predict the two-phase flow and heat transfer in a microchannel with evaporative heat transfer. Sample solutions to the model were obtained for both analytical analysis and numerical analysis. It is assumed that the capillary pressure is neglected (Morris, 2003). Results are provided for liquid film thickness, total heat flux, and evaporating heat flux distribution. In addition to the sample calculations that were used to illustrate the transport characteristics, computations based on the current model were performed to generate results for comparisons with the analytical results of Wang et al. (2008) and Wayner Jr. et al. (1976). The calculated results from the current model match closely with those of analytical results of Wang et al. (2008) and Wayner Jr. et al. (1976). This work will lead to a better understanding of heat transfer and fluid flow occurring in the evaporating film region and develop an analytical equation for evaporating liquid film thickness.

#### 1. Introduction

Over the last decade, micromachining technology has been increasingly used to develop highly efficient heat sink cooling devices due to advantages such as lower coolant demands and smaller machinable dimensions. One of the most important micromachining technologies is the ability to fabricate microchannels. Hence, the studies of fluid flow and heat transfer in microchannels which are two essential parts of such devices have attracted attention due to their broad potential for solving both engineering and medical problems. Heat sinks are classified as either single-phase or two-phase according to whether liquid boiling occurs inside the microchannels. The primary parameters that determine the single-phase and two-phase operating regimes are the heat flux through the channel wall and the coolant flow rate. For a fixed heat flux (heat load), the coolant may maintain its liquid state throughout the microchannels. For a lower flow rate, the liquid coolant flowing inside the microchannel may reach its boiling point, causing flow boiling to occur, resulting in a two-phase heat sink [1–3]. Since the initial conception of microheat pipes in 1984 by Cotter [4], a number of analytical and experimental investigations have been conducted. Most of these investigations have concentrated on the capillary heat transport capability. With rapid advancements in microelectronic devices, their required total power and power density are significantly increasing. Thin-film evaporation plays an important role in these modern highly efficient heat transfer devices. When thin-film evaporation occurs in the thin-film region, most of the heat transfers through a narrow area between the nonevaporation region and intrinsic meniscus region as shown in Figure 1. The flow resistance of the vapor phase during thin-film evaporation is very small compared with the vapor flow in the liquid phase in a typical nucleate boiling heat transfer configuration. In addition, the superheat needed for the phase change in the thin-film region is much smaller than that for a bubble growth in a typical nucleate boiling, in particular, at the initial stage of the bubble growth. The heat transfer efficiency of thin-film evaporation is much higher than the nucleate boiling heat transfer. Because thin-film evaporation occurs in a small region, increasing thin-film regions and maintaining its stability are very important for heat transfer enhancement. It is known that the thermodynamic properties of the liquid thin-film region are very different from those of the macroregion. The effect of the intermolecular forces between the liquid thin film and wall can be characterized by the disjoining pressure, which controls the wettability and stability of liquid thin film formed on the wall. A better understanding of the evaporation mechanisms governed by the disjoining pressure in the thin-film region, especially for high heat flux, is very important in the development of highly efficient heat transfer devices. As early as 1972, Potash Jr. and Wayner Jr. [5] expanded the Derjaguin-Landau-Verwey-Overbeek (DLVO) theory [6] to describe evaporation and fluid flow from an extended meniscus. Following this work, extensive investigations have been conducted to further understand mechanisms of fluid flow coupled with evaporating heat transfer in thin-film region. Stephan and Busse [7] developed a mathematical model based on the theoretical analysis presented by Wayner Jr. et al. [5, 8] to investigate the heat transfer coefficient occurring in small triangular grooves and found that the interface temperature variation plays an important role in the thin-film evaporation. Schonberg and Wayner Jr. [9] developed an analytical model by ignoring capillary pressure and found an analytical solution for the maximum heat evaporation from the thin film. Ma and Peterson [10] studied the thin-film profile, heat transfer coefficient, and temperature variation along the axial direction of a triangular groove. Hanlon and Ma [11] found that when particles become smaller, thin-film region can be significantly increased. More recently, Shao and Zhang [12] considered the effect of thin-film evaporation on the heat transfer performance in an oscillating heat pipe. Wang et al. [13, 14] established a simplified model based on the Young-Laplace equation and obtained an analytical solution for the total heat transfer in the thin-film region. In the current investigation, a mathematical model is established and its analytical solutions are obtained to evaluate the heat flux, total heat transport per unit length along the thin-film profile, thin-film thickness, and location for the maximum heat evaporation in the thin-film region and maximum heat transfer rate per unit length by thin-film evaporation.

#### 2. Theoretical Analysis

Figure 1 illustrates a schematic of an evaporating thin film formed on a wall. For the current investigation, it is assumed that fluid flow in the thin-film region is two-dimensional and pressure in the liquid film is a function of the -coordinate only. The wall temperature, , is greater than the vapor temperature, . The momentum equation governing the fluid flow in the thin film can be found by where Solving (1) with boundary condition (2), we get where is the slip coefficient. If is equal to zero, then no slip boundary condition is obtained.

So we can find where where is the limiting slip length and represents the critical value of the shear rate. The slip coefficient under the condition of this study turns out to be approximately m. From (4), the mass flow rate at a given location can be found as Taking a derivative of (6), But so the net evaporative mass transfer can be obtained as The heat transfer rate by evaporation occurring at the liquid-vapor interface in the thin-film region can be determined by But at the same time in (10) is equal to the heat transfer rate through the liquid thin film; that is, From expanding the Clausius-Clapeyron equation, we can find The pressure difference between vapor and liquid, , at the liquid-vapor interface is due to both the capillary pressure and disjoining pressure and is expressed using the augmented Young-Laplace equation: The disjoining pressure for a nonpolar liquid is expressed as where is the dispersion constant and is the film thickness. The capillary pressure is the product of interfacial curvature and surface tension coefficient : Substituting (13) into (11) the heat flux can be rewritten as Substituting (14) into (16) the heat flux can be rewritten as Substituting (19) into (10) gives We can find by differentiating (14) with respect to , and assume uniform vapor pressure, , along the meniscus: Substituting (21) into (20) yields For the evaporating thin-film region, the disjoining pressure is one dominant parameter, which governs the fluid flow in the evaporating thin-film region. And in the evaporating thin film region, the absolute disjoining pressure is much larger than the capillary pressure especially when the curvature variation along the meniscus is very small. In order to find the primary factor affecting the thin-film evaporation in the evaporating thin-film region, it is assumed that the capillary pressure, , is neglected.

Then, (22) becomes From (15) we find and ; Substituting (15) and (24) into (23) yields We can rewrite the heat flux as Simplifying (25) gives We need to solve (27) to find .

By multiplying 2 sides by and integrating two sides with respect to we get , are integration constants.

So Let To find constant we apply the boundary condition So we get From (26) and (27), So Then Note that The optimum thickness of the evaporating thin film, , where the heat flux reaches its maximum, can be found by taking a derivative of for heat flux equation (34): So Then For the nonevaporating film region, the heat flux is zero. Clearly, the interface temperature is equal to the wall temperature. The equilibrium thickness can be readily found by So From (31), (33), and (34) we get The total heat transfer rate per unit width along the meniscus, , can be calculated by From (34) So Dimensionless thickness of the evaporating region can be defined as where is the thickness of the nonevaporating region.

A dimensionless position also can be defined as For dimensionless heat flux, where is the heat flux at the interface temperature equal to the vapor temperature.

Substituting (36), (34), and (50) in (49) we found or By considering , we can determine the maximum dimensionless heat flux . The local heat flux through the evaporating thin film reaches its maximum when . Letting in (52), the maximum dimensionless heat flux, , can be found as Equation (53) indicates that the maximum heat flux occurring in the evaporating thin-film region is not greater than 0.473 times of the characteristic flux heat.

Consider the heat transfer coefficient of To solve (31) we have two methods: (1)exact solution or analytical solution,(2)numerical solution.

We are going to solve (31) by using both methods

##### 2.1. Analytical Solution

We can rewrite (31) in the form where Since so we can rewrite (55) as We introduce and to get Since , we can make the following approximation: To get we have Let Because , we neglect Therefore, the solution is accurate for as The general solution of (63) by using the initial condition is And from (65) we can evaluate an analytical solution for as

##### 2.2. Numerical Solution

The solution can be readily obtained using the fourth order Runge-Kutta method for the evaporating thin film profile. The governing equation (31) is solved with the use of a Runge-Kutta (4) method; the solution procedure is iterative. As a first guess , the values from the previous step are used, and the calculated values are returned from the Runge-Kutta solver and compared to the guess values. A comparison of the guess and calculated values is performed and looped until reaching convergence criteria for both film thicknesses. The numerical solver is coded in MATLAB. With these initial conditions, we have

#### 3. Results and Discussion

##### 3.1. Comparison of Analytical Solution and Full Model

As presented above, a mathematical model for predicting evaporation and fluid flow in thin-film region is developed. Utilizing dimensionless analysis, analytical and numerical solutions are obtained for the heat flux distribution, total heat transfer rate per unit length, location of the maximum heat flux, and ratio of the conduction to convection thermal resistance in the evaporating film region. In order to verify the analytical solution derived herein, results predicted by Wang et al. [13] and numerical solution by Wayner Jr. et al. [8] are used. Figure 2 shows the comparison of analytical and numerical results of the total heat transfer rate through thin-film region with results presented by Wang et al. [13] and Wayner Jr. et al. [8]. Total heat flux is presented in function of the superheat temperature. Our numerical and analytical solutions are compared to the ones of Wang and Schonberg and Wayner. It shows good agreement for moderate superheat temperatures; however, our analytical solution tends to underestimate total heat flux at large superheat temperatures. In addition, the model presented herein can be used to predict analytically and numerically all of the equilibrium film thickness, heat flux distribution, film thickness variation of evaporating film region, maximum total heat transfer rate through the evaporating film region, and ratio of the conduction to convection thermal resistance. The following calculations and predictions are based on the thermal properties and operating conditions shown in Table 1.

##### 3.2. Comparison of Analytic Equation for with Previous Studies

It was seen from Table 2 that, in many studies on a wide range of time, we did not find any researcher who found the analytical equation for even only approximately, but we found that they have numerical studies. So according to this, we can say that (65) is the first analytical equation for or at least approximately.

##### 3.3. Distribution of Evaporative Film Thickness

Figures 3, 4, 5, 6, 7, and 8 compare our numerical and analytical predictions for the dimensionless film thickness as function of the dimensionless position for 0.1 K, 0.5 K, 1 K, 2 K, 5 K, and 10 K superheat temperatures. They clearly show that the error of the analytical approximation is decreasing with increasing superheat temperature. We can see that for large positions the solution is dominated by exponential growth.

Figures 9, 10, 11, 12, 13, and 14 compare our numerical and analytical predictions for the dimensionless heat flux as function of the dimensionless position for 0.3 K, 0.7 K, 1 K, 2 K, 5 K, and 10 K superheat temperatures. Similar property can be seen; that is, the error of the analytical approximation is decreasing with increasing superheat temperature. The maximum of the heat flux is predicted correctly for all values of the superheat temperature analytically; however, the location of the maximum is predicted with significant error at low superheat temperatures.

Figures 15, 16, 17, and 18 show that nondimensional heat flux is presented as function of the nondimensional position for different values of : 0, 0.5, 1 and . As gets larger, the error in the analytical heat flux approximation is getting larger.

Figures 19, 20, 21, and 22 compare our numerical and analytical predictions for the dimensionless film thickness as function of the dimensionless position for different values of the slip coefficient, : 0, 0.5, 1 and . Comparing the figures indicates similar conclusions as previously mentioned: smaller slip coefficient yields smaller error in the analytical approximation.

#### 4. Conclusions

This paper presents a mathematical model for predicting evaporation and fluid flow in thin-film region. The thin-film region of the extended meniscus is delineated. Utilizing analytical solutions were obtained for heat flux distribution, total heat transfer, and liquid film thickness in the evaporating film region. The mathematical model also develops an analytic equation for . So according to this, we can say that (65) is the first analytical equation for or at least approximately. Also we can conclude that there is small effect of slip condition on the thin-film evaporation for two-phase flow in microchannel. The dimensionless heat flux through thin-film region is a function of dimensionless thickness. Also the results showed the assumption that neglecting the capillary pressure is acceptable.

#### Highlights

(i)Analytical two-phase flow for microchannel heat sink is studied.(ii)Numerical two-phase flow for microchannel heat sink is studied.(iii)New evaporating film thickness equation is developed.

#### Nomenclature

: | Dispersion constant (J) |

: | Heat of vaporization (J/kg) |

: | Conductivity (W/m·K) |

: | Interface net evaporative mass transfer (kg/(s)) |

: | Mass flow rate (kg/ms) |

: | Capillary pressure (N/) |

: | Liquid pressure (N/) |

: | Vapor pressure (N/) |

: | (N/m) |

: | Heat flux (w/) |

: | Characteristic heat flux (w/) |

: | Dimensionless heat flux |

: | Maximum dimensionless heat flux |

: | Total heat transfer rate per unit width (W/m) |

: | Temperature (K) |

: | Velocity along -axis (m/s) |

: | -coordinate (m) |

: | Dimensionless -coordinate |

: | -coordinate (m). |

*Greek Symbols*

: | Film thickness (m) |

: | Equilibrium film thickness or characteristic thickness (m) |

: | Dimensionless film thickness |

: | Optimum thickness corresponding to the maximum heat flux (m) |

: | Dynamic viscosity (N s/) |

: | Kinematic viscosity (/s) |

: | Density (kg/) |

: | Surface tension (N/m) |

: | The slip coefficient (m). |

*Subscripts*

: | Time rate of change |

: | Liquid |

: | Maximum quantity |

: | Total |

: | Vapor |

: | Wall. |

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgments

This work is supported by Universiti Teknikal Malaysia Melaka (UTeM) and Thi Qar University.