3D numerical investigation on flow boiling characteristics in large length-to-diameter ratio microchannel

To meet the demand of thermal management with large power and high heat flux in the future engineering applications, more and more attentions have been draw to develop the two-phase microchannel cooling technology. In this work, a three dimensional numerical investigation is performed to study the influence of mass flux and heat flux on the two phase flow boiling characteristics in a single microchannel with large length-to-diameter ratio. The related calculation results show that the Volume of Fluid (VOF) model coupled with Lee model can predict the two phase pressure drop and heat transfer coefficient of the microchannel well. And the flow patterns in the microchannel are classified as bubbly flow, slug flow and annular flow, which can be captured by this numerical model precisely. Besides, the transient pressure drop augments with the increase of mass flux. But the effect of mass flux on the amplitude of quasi-steady pressure drop and the peak value of time-averaged heat transfer coefficient along the microchannel is small. At last, the amplitude of quasi-steady pressure drop increases with the augmentation of heat flux. The maximum and fully developed time-averaged heat transfer coefficients along the microchannel increase with the heat flux.


Introduction
Microchannel cooling technology is a promising scheme to solve the bottleneck of thermal management with large heat power and high heat flux in various industrial applications like photovoltaic cell, electronic equipment, solar panels and spacecraft [1].Compared to that with single-phase Mechanical Pumped Fluid Loop (MPFL) at significant heat loads, microchannel integrated with Mechanical Pumped Two-phase Loop (MPTL) thermal management systems can reduce the pump power and meanwhile improve the heat transfer performance [2].However, the designer of MPTL need to master the flow boiling characteristics of the microchannel heat sink so as to prevent pre-mature Critical Heat Flux (CHF) and flow instabilities [3].To relieve the problems of two phase flow and heat transfer in the microchannels mentioned above, some recent studies have been centered on developing Computational Fluid Dynamics (CFD) approach for prediction of two phase flow pattern, void fraction, pressure drop, heat transfer coefficient and CHF [4][5][6][7][8][9][10][11][12][13][14][15][16].
Jeyaraj and Pankaj [4] conducted a 2D numerical study to evaluate the prediction performance of the Mixture Model (MM), VOF model and Eulerian Model (EM) on flow boiling in a microchannel.
The results showed that the MM was the best model to predict the void fraction compared to the other two models.Yang et al. [5] used the VOF model and the modified Lee model [6] to study the two phase flow and heat transfer characteristics for NH 3 /H 2 O mixture in a horizontal two-dimensional single microchannel.The results indicated that at some cases, it was benefit to used NH 3 /H 2 O mixture as an alternative coolant to prevent local drying of electronic components and maintain the temperature within their limits.Alugoju et al. [7] used the VOF model to numerically simulate the two-dimensional expanding microchannels with uneven heat flux.The performance of straight channel and expansion channel under similar heat flux boundary condition were compared.They reported that in the expanding channels, more small-diameter bubbles could be captured.However, fewer large-diameter bubbles in the straight channel could be observed.Besides that, the nucleation time in straight tube was shorter.
In order to predict the flow boiling characteristic in microchannel precisely, more and more researcher have been conducted their studies on two phase flow in microchannel by 3D numerical simulation.Ferrari et al. [8] used the VOF method to conduct a numerical study on the slug flow and heat transfer characteristics in a microchannel with square cross section.It was illustrated that the channel geometry had a great effect on the bubble dynamics and heat transfer characteristics.Magnini and Matar [9] systematically analyzed the influence of channel structure on bubble dynamics and heat transfer at flow boiling state.The dynamics interface of liquid phase and vapour phase were tracked by VOF model.The results emphasized that the thermal performance of various channel geometry is closely connected with the circumferential liquid film distribution around the ultra-long bubble.Luo et al. [10] used the VOF model to conduct a basic numerical study on the hydrodynamics and thermal performance of two phase flow in a large width-to-depth ratio microchannel with rectangular cross section.The results showed that the thickness of liquid film between the heated wall and the interface decreases with the increase of inlet quality or heat flux.Ariyo and Bello-Ochende [11] used CFD method of non-equilibrium boiling model to simulate the subcooled flow boiling.It was reported that the microchannel could operate exceed the optimized heat flux.Yan et al. [12] used VOF method and Lee model to study the influence of aspect ratio and non-uniform heat flux boundary condition on flow boiling characteristics.The results showed that the flow patterns in the microchannel could be classified as dispersed bubbly flow, bubbly flow and restricted bubbly flow.The distribution of bubble shape had some relationship with the ratio of width to height.Lorenzini and Joshi [13] presented a CFD approach using a three-dimensional VOF model.This model is used to evaluate the two-phase cooling performance of silicon microchannels with uneven heat flux.The results showed that this model could be used for the trend prediction of flow boiling in microchannels, with remarkable capabilities for two-phase flow visualization, evolution and tracking.Lin et al. [14] studied the growth of slug bubbles in expanding and uniform cross section microchannels by using the VOF model coupled with the phase change model provided by Hardt and Wondra.The results showed that the diverging of the cross section could accelerate the movement of bubbles to the outlet.As a result, the flow instability and channel chocking were reduced.Alugoju et al. [15] used the VOF model to conduct a three-dimensional numerical investigation of expanding microchannel heat sink.The results illustrated that compared with the straight microchannel, the expanding microchannel had a better performance in suppressing the flow instability.However, all the numerical studies mentioned above are focused on small length-to-diameter ratio microchannel compatible with electronics devices.None of them are about heat sinks consisted of several parallel large length-to-diameter microchannels, which used for future engineering applications that could dissipate large power waste heat from concentrated heat sources.
In fact, to the author's knowledge, only Lee et al. [16] studied the flow boiling characteristics in a large length-to-diameter ratio microchannel by CFD simulation, while the model is two-dimensional.
The primary motivations for this study are indicated as below: 1) Develop 3D CFD model for simulating flow boiling heat transfer process in a large length-to-diameter ratio microchannel.
2) Explore the two phase flow boiling characteristics in a large length-to-diameter ratio microchannel by 3D numerical simulation.

Physical model
The 3D flow boiling heat transfer process in a single rectangular microchannel is simulated.The structure and geometric parameters of the studied microchannel used for numerical simulation are illustrated by Fig. 1.The microchannel is made of AlSi10Mg with a length (L) of 330mm.The inlet width (W f ), height (H f ) and wall thickness (δ s ) are 0.76mm, 0.9mm and 0.2mm, respectively, so the hydraulic diameter (D h ) of the microchannel is about 0.82mm.As a result, the length-to-diameter ratio (L/D h ) is about 402, which is much larger than the conventional microchannels.The ammonia, widely used in aerospace heat pipe, is adopted as the Heat Transfer Fluid (HTF).The material property of AlSi10Mg and ammonia are shown in Table 1.
. Table 1.Material property of AlSi10Mg and ammonia

Numerical model
The three dimensional continuity, momentum and energy governing equations are indicated as below: Continuity equation: where ρ, t and u represent the density, time and instantaneous velocity respectively.Momentum equation: where p, μ, F and g are the static pressure, dynamic viscosity, surface tension force and gravity, respectively.
Energy equation: where E, T and k eff are the energy per mass, temperature and effective thermal conductivity, respectively.S is the source term, which equals the product of mass transfer rate of unit volume and vaporization heat.
The VOF model coupled with Lee model is adoped to study the two-phase flow boiling in the microchannels.The tracking of the interphase interface is realized by solving the volume fraction equation for the secondary phase of vapor ammonia: Volume fraction equation: where α is the volume fraction.The subscript v and l represents the vapor phase and liquid phase, respectively.
The volume fraction of liquid ammonia is calculated according to the following formula: The Lee model is expressed as follows: where r e , r c , T v , T l and T sat are the evaporation empirical coefficient, condensation empirical coefficient, vapor temperature at the interface, liquid temperature at the interface and the saturation temperature, respectively.
Based on the assumption that the interface temperature is equal to the saturation temperature, the empirical coefficients r e and r c are set as 100s −1 so as to keep the interface temperature within T sat ± 1K during the simulation.
Taking effect of surface tension in this work, Continuum Surface Force (CSF) model established by Brackbill et al. [17] is adopted in our simulation.
The boundary condition for the inlet is the mass flux in the subcooled state.The boundary condition for outlet is the pressure boundary condition at saturation state.The top wall is adiabatic, while the double sides of the walls are symmetry conditions.And no slip boundary conditions are used for the inner wall.
The Finite Volume Method (FVM) is used for discretizing the governing equations.The PISO algorithm is implemented for coupling velocity and pressure field.The second upwind scheme and QUICK scheme are adopted for discretizing the convective terms in momentum and energy equations, respectively.The pressure is discretized by PRESTO scheme.In volume fraction calculation, the geometric reconstruction scheme is adopted.The first-order implicit scheme is used to solve the unsteady term.The convergence criterion for the mass, momentum and energy equations is that the residuals are less than 10 -3 , 10 -5 and 10 -8 , respectively.For numerical stability, the global Courant number is maintained less than 1.0.As a result, the time steps ranged from 10 -7 s to 10 -6 s.

Data reduction
The microchannel module's pressure drop, △p, is defined as follows: where p in and p out are area-averaged pressure at inlet and outlet of the microchannel, respectively.The heat transfer coefficient of the microchannel, h, is defined as follows: where q w , T f and T w are the heat flux based on the inner wall of the microchannel, the mass-averaged temperature of the fluid and the area-averaged temperature of the inner wall, respectively.

Grid independence checking
As shown in Figure 2, the whole microchannel is meshed by orthogonal hexahedral grids.In the streamwise direction, the mesh is uniform.However, the mesh near the wall is refined gradually so as to capture the liquid film.Grid sensitivity analysis is done by increasing the number of cells by successive refinements from one case to another until a grid-independent criterion of one or more related variables is obtained.The selected operating conditions of mass flux (G), heat flux (q) and saturation temperature are 300kg/(m 2 •s), 20000W/m 2 and 303.15K, respectively.It is important to note that the contact angle of AlSi10Mg with liquid ammonia is 10.8º , which was measured experimentally.When the flow time t=1000ms, the quasi-steady state condition was achieved.
The criteria of grid independence checking is |(N j -N j+1 )/N j |<10 -2 , where N j and N j+1 represent the concerned field variable for the current grid and next grid, respectively.The pressure drop along the microchannel and the average heat transfer coefficient of the inner wall were selected for comparison.
Table 2 shows the results of grid independence checking.So as to balance the solution accuracy and the calculation time, the third grid (No.3) was adopted for our numerical study.

Model validation
The reliable of VOF model with Lee phase change model in predicting flow boiling characteristics in microchannel has been validated by a lot of studies.In this paper, the validate calculation of our two phase flow boiling model are performed for a stainless steel circular tube with D h =1.224mm, L=245mm and uniform heat flux by comparing the numerical results with empirical correlation [18][19][20] at G=100kg/(m 2 •s), q=10341~30750W/m 2 , T sat =316.15K and θ=19.3º[21].As Table 3 shows, the relative errors of ammonia flow boiling pressure drop between present results and empirical correlation provided by Maqbool et al. [18] are within 15.9%.And the relative errors of ammonia flow boiling heat transfer coefficient between the numerical results and empirical correlation proposed by Cooper [19] and recommended by Maqbool et al. [20] are within 18.98%.Considering that the uncertainties of Maqbool correlation and Cooper correlation are 30% and 20%, respectively, the presented numerical model can be used to simulate the flow boiling performance in microchannel.

Table 3.
Comparison of the numerical results with previous studies.

Results and discussion
In this section, the influence of mass flux and heat flux on ammonia two phase flow boiling characteristics in a large length-to-diameter ratio microchannel will be discussed.

The effect of mass flux
Fig. 3 illustrates the transient inner wall temperature field and two phase flow regime at various mass flux with a subcooled degree (∆T sub ) of 3K.It can be seen that as the inner wall temperature is superheated, which reveals the local fluid temperature achieves the saturation temperature, the flow transforms from subcooled to saturated flow boiling (Fig. 3(a-i)).As time goes on, the bubbles emerge, then rapid growth, expansion and coalescence in the axial direction, which sequentially lead to the bubbly, slug and annular flow (Fig. 3(g-i)).However, as the mass flux increases, the annular flow region gradually decreases until it disappears (Fig. 3(c, f, i)).And the distance between the inlet and the location of bubbly flow emergence increases (Fig. 3(a, d, g)).The reason is that the heat absorbed by the fluid entering from the inlet is mainly related to the flow time, and the local flow pattern is closely dependent on the amount of heat absorbed by the fluid.The larger mass flow rate, the longer distance traveled at the same time interval, which will delay the evolution of the flow pattern in the axial direction.Fig. 4 presents the transient pressure drop at different mass flux.It is obvious that as the fluid at the inlet starts moving, the pressure drop rapidly increases to its maximum value, then decreases quickly and followed remain almost unchanged.Up to now, it is single-phase laminar flow.While at t=370ms,     It illustrates that the peak and quasi-steady two-phase pressure drop both increase with the increase of heat flux.And it should be noted that as heat flux is 10000W/m 2 , the two-phase pressure drop does not decrease after reaching its peak value, but reaching a quasi-steady state.This indicates that the gradual decrease in two-phase pressure drop from the peak point for other cases is likely related to the transition from bubbly flow to slug flow and slug flow to annular flow.Besides, we can also see that the amplitude of quasi-steady pressure drop increases with the increase of heat flux.G=200kg/(m 2 •s), T sat =303.15K,∆T sub =3K.

Conclusions
In this work, a three dimensional numerical simulation is performed to investigate the influence of mass flux and heat flux on the flow boiling characteristics in a single large length-to-diameter ratio microchannel.The related calculation results show that the VOF model coupled with Lee model can predict the two-phase pressure drop and heat transfer coefficient well.The flow patterns in the microchannel are classified as bubbly flow, slug flow and annular flow, which can be captured precisely by our numerical model.Besides, the transient pressure drop increases with the augmentation of mass flux.But the influence of mass flux on the amplitude of quasi-steady pressure drop and the peak value of time-averaged heat transfer coefficient along the microchannel is small.At last, the amplitude of quasi-steady pressure drop augments with the increase of heat flux.The maximum and fully developed time-averaged heat transfer coefficients along the microchannel augments with the heat flux.We will combine more experimental data to improve the numerical model and carry out a further systematic investigation on the two phase flow boiling mechanism in the microchannel with large length-to-diameter ratio.

Figure 1 .
Figure 1.The 3D structure and geometric parameters for the studied microchannel.

Figure 2 .
Figure 2. Grids in the cross section.
.1088/1742-6596/2599/1/012024 7 the pressure drop gradually increases, indicating that the fluid begins to vaporize.Until t=530ms, the pressure drop gradually decreases and finally reaches the quasi-steady state.In the entire process, the pressure drop increases with the increase of mass flux.The pressure drop oscillation mainly occurs after t=530ms, and the effect of mass flux on amplitude is small.Fig. 5 presents the axial distribution of time-averaged heat transfer coefficient with different mass flux.It reveals that the peak value of time-averaged heat transfer coefficient is less affected by the mass flux.And as the mass flux increases from 100kg/(m 2 •s) to 200kg/(m 2 •s), the time-averaged heat transfer coefficient in axial fully developed region increases slightly.However, as the mass flux increases from 200kg/(m 2 •s) to 300kg/(m 2 •s), there is no significant change for the time-averaged heat transfer coefficient of axial fully developed region.The axial distribution of time-averaged heat transfer coefficient mainly due to the entrance effect and the evolution of flow pattern.

Figure 3 .
Figure 3.The transient inner wall temperature field and two phase flow regime at different mass flux: T sat =303.15K,∆T sub =3K, q=20000W/m 2 .

Fig. 6
shows the transient inner wall temperature field and two phase flow regime at t=1574ms with various heat flux.It indicated that as the heat flux is 10000W/m 2 , the flow regime in the microchannel is bubbly flow.As the heat flux augments to 20000W/m 2 , the flow patterns along the axial direction in the microchannel are bubbly flow and slug flow, and the location where bubbly flow emergence moves towards the inlet.When the heat flux augments to 30000W/m 2 , the starting point of flow boiling further moves towards the inlet, and the annular flow appears after the slug flow.Fig.7 presents the transient pressure drop at various heat flux.

Fig. 8
Fig. 8 illustrates the axial variation of time-averaged heat transfer coefficient with different heat flux.It shows that the maximum and fully developed time-averaged heat transfer coefficient along the microchannel increase with the increase of heat flux.However, the axial variation trend of time-averaged heat transfer coefficient under different heat flux is similar, although the flow patterns near outlet varies under different heat flux.