Boiling heat transfer simulation in rectangular mili-channels

Due to the high heat transfer coefficient and compactness of a system, mili-channel-based cooling and heating techniques are greatly expected to be distributing high heat flux from the electronic devices. In terms of cooling performance, the two-phase evaporating flow of boiling flow in mini and mili-channels is more effective than the single-phase flow due to the inclusion of latent energy in the process. In this study, a numerical model was proposed to simulate the boiling heat transfer of multiphase flow in a channel using different boundary conditions in the channel surfaces. The fluid volume approach regulating the hydrodynamics of the two-phase flow was used. Source terms of the energy and mass transfer that were taken into account at the interface of liquid and vapor were included in the management equations for the conservation of energy and vapor quality. A 3D Ansys-Fluent © simulation model was developed and numerical simulations were conducted for four different boundary conditions. A mili-channel with a length of 140 mm was used. The liquid and gas phases that were used in the model were liquid water and vapor; the total mass flux at the inlet was varied at 118–126 kg/m 2 s. In order to realize thin film annular flow over the boiler surface, employed specific boundary conditions in the 3D simulation model were obtained by means of one dimensional Matlab © simulation code. By means of utilizing the evaluated numerical results, distribution of heat transfer coefficient, vapor quality and dimensionless temperature over the heat transfer surfaces were reported and compared to experimental results. Numerically evaluated results are in agreement with experimentally measured results. For the studies cases an average value of 23600 W/m 2 .K was obtained for the heat transfer coefficient. Cite this article as: Koca A, Khalaji MN, Sepahyar S, Boiling heat transfer simulation in rect-angular mili-channels. J Ther Eng 2021;7( 6 ):1432–1447.


INTRODUCTION
The miniaturization of electronic components and the simultaneous performance increases have led to increased volumetric heat dissipation requirements and therefore to the need for more compact and efficient thermal management systems.Investigation of heat transfer performance of the mili-channels are available in the literature [1].There are three different methods have been used to investigate the heat transfer in the mili-channels: Microchannel flow boiling, spray cooling and jet jamming have received great attention in recent years.Many previous researchers have investigated heat transfer of the mili-channels by means of experimental measurements, CFD simulations and predicted heat transfer correlations.Tuckerman and Pease [1] first demonstrated single-phase microchannel cooling using water in 1981 and dispersed heat flows up to 790 W cm -2 .Unfortunately, single-phase cooling requires very high flow rates and therefore produces very high pressure drops that require high pumping power.This also requires large liquid stocks that do not lend themselves to compact thermal management systems.In another study of electronic cooling, the permissible boiling is a limitation on surface temperatures.For example, typically 75-85°C is the maximum allowable temperature for boiling surfaces connected to cold plates used in electronic cooling [2].
Surface heat transfer mechanisms and fluid flow are regulated by bubbles due to their ability to occupy the entire cross-section of the channel.High heat removal by using low refrigerant flow rate and temperature homogeneity in the ducts are two positive features of the two-phase cooling technique.Based on the International Semiconductor Technology Roadmap estimates in 2011, assuming today's average chip die sizes, we can estimate the average heat flow of up to 450 W cm -2 for microprocessor chips by 2026 and the highest heat flows up to 4.5 kWcm -2 [3].
Ganapathy et al. [4] reported a numerical simulation of micro-channels for condensation heat transfer.The fluid interface used monitoring volume and reported that the pressure drop and the Nusselt number deviated on average by 8.1% and 16.6% by empirical correlation.Bevis and Bandhauer [5] showed the distribution of heat fluxes up to 1.1 kW cm -2 by flow boiling in a series of parallel micro-channels.Spray cooling and jet shock demonstrated the ability to use higher heat flows (12 kW cm -2 [6] and 18 kW cm -2 [1], respectively), but the complexity of creating robust and durable systems with large heat transfer areas is the most are great obstacles.
Kandlikar [7] experimentally studied multi-channel boilers and observed the flow of water.Experimental setup poor six parallel channels with a cross-section of 1 * 1 mm 2 .In some channels, localized flow change was observed due to large fluctuations in pressure drop.Kandlikar [8] studied the boiling flow in micro chemicals and proposed two dimensionless groups of K1 and K2; wherein K1 is the ratio of the evaporation moment force to the inertial force and K2 is the ratio of the evaporation moment force to the surface tension force.Nucleate reported that due to the periodic flow of vapor and liquid deposits formed in the boiling process, the boiling source is reduced due to the conventional process.
Mili flow boiling heat exchangers, on the other hand, can easily be stacked in robust compact systems.The additional advantages of relatively low pumping power and small fluid inventory requirements enable the mili-channel flow to boil one of the most promising thermal management strategies for the future.Many flow boiling studies have been done and many things have been learned, but it is necessary to start applying the acquired knowledge to real industrial situations.Most of the boiling experiments were carried out with homogeneously applied heat flux along the length of the duct.This allowed the heat sink material to reasonably assume the conduction of heat during the analysis of the results, thus allowing many researchers to obtain a duct wall heat flux equal to that of the applied flow and calculate duct wall temperatures [9][10][11][12][13][14][15][16][17].Hoffman and Stephan [18] used thermodynamic crystals to measure temperatures below an evaporating meniscus and observed a decrease in temperature due to the high evaporating heat flux around the micro domain.Bogojevic et al. [19] conducted experiments on smooth cross-section micro and mili-channels, and flow boiling instabilities were investigated.Mukherjee and Kandlikar [20] presented a numerical analysis of micro-channels with an inlet contraction to suppress instability in flow boiling.They observed that the bubble caused by the formation of vapor tends to move towards the limitless end.In many studies, the researchers had obtained accurate results for condensation of steam in millimeter-scale ducts and tubes with various geometries and boundary conditions.In these studies, there is a focus on heat transfer in terms of theoretical, experimental and numerical modeling.In their studies, they investigated the phase change heat transfer mechanism in order to increase performance of the thermal devices such as heat exchangers, flat plate solar air collectors and other electronic devices, which play an important role in the industry today [21][22][23][24][25][26][27][28][29][30][31][32].In this regard, they carried out comprehensive studies on laminar or turbulent flow-forced convection in mini and micro channels of various shapes, such as; vertical, horizontal or inclined.
From the above review, it has been found that numerous experimental studies of multi-phase flows in the mini and mili-channel focus on boiling and condensing heat transfer.In this study, a numerical model for evaporation heat transfer in the national channel was proposed and the results were compared from the experimental studies.
The motivation of this research can be explain as follows: The conventional two phase cooling technologies used for critical applications such as: aviation and aerospace applications include pumped two-phase loops, heat-pipes and loop heat-pipes are all inactive but reliable systems.Though these systems will not able to meet future challenging cooling needs, since they have inherent limitation of capillary pumping mechanism in regards to transport distance and heat transfer.
In spite of that, functionality problems with boiling and condensing flows arise because, at the millimeter and micrometer scale hydraulic diameter and other operating conditions of interest, shear/pressure forces dominate over gravitational forces and give rise to [33]: • Thermally ineffective and hydro-dynamically problematic liquid-vapor configurations -such as long lengths of non-annular (plug/slugs, etc.) regimes • Inability to reduce the significantly higher sensitive coupling between vapor and liquid motions • High consumption of pumping power for small hydraulic diameters.• Inability for removing large amount of heat for the flows where the gravitational forces dominate other forces or the other forces dominate gravitational forces.
For this reason, innovative system designs for removing large amounts of heat from small spaces (as in electronic cooling, etc.) and resolving the issues, aforementioned above, for critical applications are dependent on breakthroughs in the development for phase-change devices such as new millimeter-scale flow boilers and condensers.
A key innovation of the proposed innovative boiler is the operations can successfully use re-circulating vapor flows to ensure that thermally and hydrodynamically efficient annular flows are realized over most of the devices' heat-exchange surfaces for the boiling flows.These surfaces are therefore continuously irrigated by thin liquid films [13].Besides, the ability to properly control the recirculating vapor flow rates associated with the proposed innovative boiler allows these devices to operate completely in the annular flow regimes and design flexibility.To achieve this, the inlet liquid (properly guided and at near saturation temperature) flow rate is chosen to be consistent with the heat load, and inlet vapor flow rate is chosen to lie within a range that ensures annular boiling [13].
In this regard, for the design of millimeter scale boiler, one dimensional Matlab © simulation code were developed by using the existing theories and the correlations in the literature.This simulation code allowed to design boilers for different areas of usage and optimizing their operational conditions for different boundary conditions.In order to realize thin film annular flow over the boiler surface, employed specific boundary conditions in the 3D simulation model were obtained by means of the Matlab © simulation code (Theories behind the Matlab © code can be find in [33]).After then, three dimensional simulations were conducted for four cases and the results were compared with the experimental results.

THEORY
The model consists of a rectangular mili-channel with actual correlations.Differences in the vapor quality of vapor at different locations along the centerline of the mili-channel for different wall boundary conditions are presented, and FLUENT version 16.0 (ANSYS, Inc.) was used to perform the simulations.Continuity, momentum, and energy conservation equations were solved by numerical model for multiphase flow are given in Equations (1), (2) and (3), respectively.
Continuity equation: Momentum Equation ( ) ( ) Energy conservation equation: ρ and µ, fluid density and viscosity can be written in (ρ = ρ l α l + ρ v α v ) and (µ = µ l α l + µ v α v ), respectively.Herein and prey are the liquid and vapor fraction, respectively.Sh is the energy source term associated with the liquid-vapor phase change determined as follows: And effective thermal conductivity (k eff ) is given by: The finite volume method is the solution of partial differential equations with algebraic equations was used, which is very similar to the finite difference method.Finite volume method, volume integrals containing deviation term in partial differential equations are converted to surface integrals by divergence theorem.These terms are considered as flows on the surfaces of each finite volume.This method is especially used in computational fluid mechanics problems.Most computational studies used computational fluid dynamics (CFD) modeling using level set (LS) or Volume of Fluid (VOF) techniques.Detailed descriptions of these techniques can be found elsewhere [34].In this study the VOF method was used to model the two-phase flow (Hirt [35]).Continuous surface force (CSF) technique was used to model the surface tension forces (Brackbill [36]).It is not trivial to form developing fluid/vapor interface regions that flow simultaneously along the length of a channel.Because of this complexity, most CFD studies have studied single bubbles and single nucleation sites to learn about bubble formation, biphasic flow dynamics, and heat transfer.The surface tension is due to the difference in molecular positions between the two contact phases which Where, (h st ) is the total local enthalpy at the point of measurement, (h g ) and (h 1 ) are the enthalpies of vapor and liquid respectively.

NUMERICAL ANALYSIS
The three-dimensional rectangular mili-channel is shown in Figure 1.The characteristic diameter of the mili-channel is 5 mm high and 10 mm wide and 140 mm long.The number of mesh used in the simulation had approximately 6*10 6 cells and the time step was 1*10 -2 .Depending on the internal pressure of the system at the inlet of the duct, the inlet condition of 106-113 kPa and liquid temperature 371K and the vapor temperature of 373K are given.Also, a total mass flow inlet condition of 118-126 kg/ m 2 s is specified for the liquid and vapor phase and is given as an outlet pressure at the outlet of the duct.In each simulation, different boundary conditions were determined on the channel walls that kept the other parameters the same.The effects of wall conduction were negligible and therefore neglected.The number of meshes to capture liquid films during annular streams was made high near the channel walls.At the inlet of the channel, a sprinter was used to separate the liquid phase from the vapor phase.Water (liquid and vapor) was used as the working fluid which obtained thermophysical properties at saturation temperature corresponding to 106-113 kPa working pressure from the FLUENT © material database.The properties of the working fluids used in numerical simulations at the saturation temperature corresponding to the system pressure of 106-113 kPa are shown in Table 1.
In addition, constant temperatures as the values of in between 383-397 K were given from the bottom wall surface.In the simulation, the mass flow of the incoming liquid was kept constant at a flow rate of 118-126 kg/m 2 s having an inlet temperature of 371K for the liquid and inlet temperature of 373 K for the vapor and an atmospheric outlet pressure.Mesh validation study was performed in numerical simulation as shown in Figures 2 and 3.The SIMPLE algorithm was used for pressure-velocity matching and the QUICK scheme for convective flows and the first-order implicit scheme for time separation.Calculations were performed using the CFD package FLUENT 16.0.User-defined functions have been developed to include terms of terms from the equations in (3) and (7).Inlet liquid water and water vapor, system pressures and constant temperatures, and the entering steam quality ratios are given in Table 2 for each analyses cases (Cases 1-4).
Table 2 shows the input and output conditions leading to numerical CFD and experimental analysis were performed under these conditions.These boundary conditions were taken from the experimental studies of Sepahyar [39], where the details of the experimental set-up, methodology, calculation of the parameters and uncertainty analysis can be found.cause a pressure jump across the interface when the interface is curved.The volume fraction of phase (i) in each cell is represented by (α i ), the sum of the volume segments of all phases adding up to 1, as given in Equation (6): In a liquid-vapor two-phase flow, when a cell is filled with liquid, the liquid volume ratio is (α l = 1) and the corresponding vapor volume ratio is (α v = 0), then (α l = 0) and (α v = 1), a cell is filled with vapor.If the cell is partially filled with liquid and with the remaining vapor, 0,1 and (α v ) are between 0 and 1 indicating that the interface passes through the cell, and the volume of the Fluid model can be written in an equation for the volume fraction given by α: ( ) The very high latent heat of evaporation due to many fluid phases changes greatly reduces the problems associated with the distribution and pressure drop in the flow source using a relatively low flow rate compared to single-phase systems.To account for the mass and heat transfer due to phase change, the mass and heat sources in the equations were defined as [37,38].The liquid vapor phase change at the interface was accompanied by mass transport between the two phases and the release (condensation) or absorption (evaporation) of the latent heat: For liquid mass source: For vapor mass source: Here (f o ) is an adjustable parameter used to reduce the temperature difference between T and T w to negligibly small values.It should be noted that S 1 and S v should be summed to zero for mass conservation and the heat source is then obtained by equation ( 4).Using the mass and area of the velocity and the enthalpy depending on the equation ( 9) can also be written the local vapor quality of the refrigerant system: and To solve a specific problem, it is necessary to create geometry, assign material, apply boundary conditions to exterior surfaces, create a mesh defining element shape and size, and then the software must solve it repeatedly until convergence.The final step before the solution is the construction of a suitable network defining the element size, shape and number.The final mesh shown in Figure 2 is produced by a quadratic equation at one end of the pattern, which is applied under the same conditions along the channel.A thinner and more precise mesh can be applied to the input/output and lower boundary wall regions without loss of precision.A mesh validation study was performed to ensure the accuracy of the solution because the mesh     quality increases as the fluid moves along the length of the mili-channel as time passes, which may provide sufficient time to evaporate a greater amount of water.Constant temperatures in Table 2, the vapor quality change over the time (at 120, 420 and 720 seconds) because of the inherent matter of the transient analysis.When the simulations reached the steady state conditions (Fig. 5c) final vapor quality distributions were obtained.
The average temperature distribution of each plate is shown in Figures 10-14 for different cases at different intervals of time (120 seconds, 420 seconds and 720 seconds).As can be seen from the figures, the temperature increases along with the vapor quality throughout the channel length.The temperature, and vapor quality are the highest at 120 seconds and this decreases at 420 and 720 seconds, and the temperature and vapor quality are almost indistinguishable in the last seconds and the system becomes stable.The inlet water temperature was maintained near the saturation temperature (2-3°C lower).
Figure 15 represents the experimentally measured [39] and numerically calculated vapor quality distributions.As can be seen in the figure, numerically calculated vapor quality values are compatible with the experimental results.Where, in Case 4 steam quality is lower than the other cases which is due to the higher liquid mass vapor flow rate, vice versa Case 1 has more steam quality values than the other cases.
The change in local dimensionless temperature along the mili-channel in the thermally developing region is given in Figure 16.Both numerical and experimental dimensionless temperatures decrease along the mili-channel.In the experimental study, it is important to define the level of heating and the "method of heating".For non-uniform temperature-controlled heating, a specific "method of heating, " whereas for uniform temperature heating that specific function is expected to be θ w (x) = 1 over 0 ≤ x ≤ L in ideal cases [33].According to experimental studies, the dimensionless temperature θ w (x), which is within ± 5%, is shown in Figure 16 based on all data obtained and reported in Table 2.
design has a major impact on both calculation time and accuracy.After each iteration during a solution, it may take from a few seconds to a few hours or days, depending on the number of elements.A mesh with a very large element can solve in a few seconds, but it can give very wrong results.As the optimum mesh count increases and the element size decreases, the calculation time and accuracy increase.
The quality criteria that are taken into consideration for the element quality while creating the solution network structure are; skewness and orthogonal qualities and non-dimensional wall distance parameter (y + ).For the boundary layer (inflation), the parameter (y + ) is also examined, specifying the zone of Wall Law where the governing equations are solved since the fluid shear-stress and the heat transfer greatly affected by the near wall refinement.The most important parameter affecting the element quality is the skewness value.This value shows how close it is to the ideal geometry (equilateral triangle or square).However, if this value is as low as possible, it will improve the element quality.In this regard, for the chosen mesh configuration, maximum skewness value of 0.28, minimum orthogonal quality of 0.8 and y + value of 4.7 were provided.
The most basic and accurate method for evaluating mesh quality is to achieve a critical result, ie the mesh we use is good until the solution is stable or the results do not change significantly with each thinning.For example, the outlet temperature and vapor quality of a 3D spindle channel model are shown in Figure 3 for Case 3. In this case, the number of meshes affect the results, but after a certain number of meshes, results are not affected by the mesh structure.

RESULTS AND DISCUSSIONS
Simulations were performed according to the predefined inlet and outlet conditions of mass flow inlet, system pressure operating pressure, constant temperature defined from the bottom wall (383-397 K) and the isothermal condition of the surrounding walls.Figure 4 shows the placements of the planes along with the channel, which were used to evaluate the cross-sectional average vapor quality at different locations of the channel.
Figure 5 shows the contour of the vapor quality distribution at different locations of the channel in different regions.At certain locations, the vapor quality near the bottom wall is high and then decreases towards the middle of the mili-channel.Figures 6, 7, 8 and 9 show the variation of the vapor fraction along the length and in the mili-channel centerline (averages of the values of each planes) for different cases indicated in Table 2.As the mixture comes into contact with the higher temperature walls, gradual evaporation occurs within the mili-channel, and as a result, as the liquid advances along the spindle channel, the formation of vapor increases and consequently the vapor quality increases with respect to the length.The change in vapor       Local enthalpy value (h st ), is the important parameter in the calculation of the steam quality (Eq.10). Figure 17 shows the experimentally, numerically obtained and calculated local enthalpy values for each case.
Important themo-hydrolic parameters are given in Table 3.In the table, experimentally [39] and numerically calculated parameters are compared.Whereby, the most important parameters in the table are; calculated outlet vapor qualities (X out ) and the average wall heat fluxes (q′ w ).In the experimental study, the integration of local 1-D energy equation as discussed in Ch 2 [39], yields reliable estimates of quality X(x) variations -and this is shown in Figure 15 and Table 3 for the representative cases in Table 2. Table 3 shows that a representative outlet quality (X out ) not       2) and outlet quality X out .
Figure 18 shows the velocity contour and the vector for the Case 1 model based on the different time steps.As   [CrossRef]

Figure 1 .
Figure 1.(a) General view of the experimental set up and (b) Mili-channel schematics and dimensions [39].

Figure 2 .
Figure 2. Mesh construction of the mili-channel.

Figure 4 .
Location of the predefined planes.

Figure 6 .
Figure 6.Vapor quality along the length of case-1 mili-channel.

Figure 7 .
Figure 7. Vapor quality along the length of Case-2 mili-channel.

Figure 8 .
Figure 8. Vapor quality along the length of Case-3 mili-channel.

Figure 9 .
Figure 9. Vapor quality along the length of Case-4 mili-channel.

Figure 15 .
Figure 15.Experimental and numerical Vapor Quality (X) values along with the channel length for cases 1 to 4 marked in Table1.

Figure 16 .
Figure 16.Experimental and numerical Non dimensional Temperature θ Exp-distance (m) curves for cases 1 to 4 marked in Table 1.

Figure 17 .
Figure 17.Experimental, numerical, and calculated local enthalpy at the outlet for cases 1 to 4 marked in Table1.

Figure 18 .
Figure 18.Velocity vectors and contours along the mili-channel length according to Case-1 conditions at different times.

Table 1 .
Properties of fluids used in numerical simulations

Table 2 .
Input and output conditions used in numerical simulations[39]

Table 3 .
Comparison results of the experimentally and numerically evaluated parameters