Experimental and Computational Investigation of Flow Boiling in a 52 μm Hydraulic Diameter Microchannel Evaporator With Inlet Restrictions and Heat Spreading

As advances in electronics production have resulted in increasingly miniaturized and higher performance devices, thermal management strategies must improve commensurately to facilitate safe, reliable, and predictable operation. Power semiconductor devices, such as laser diode arrays, logic and switching transistors, spacecraft and avionic electronics, and other Micro-electro-mechanical systems (MEMS) electronics, are known to experience large, often nonuniform heat fluxes in localized regions within the devices which can exceed 1 kW/cm² [1,2]. Effective thermal management for these devices must address the large heat fluxes while maintaining a small form factor. Microchannel flow boiling is an effective high heat flux thermal management strategy due to its high heat transfer coefficients [3], large surface-area-to-volume ratios, improved device temperature uniformity [4], and lower pumping power requirements due to decreased mass flow rates relative to single phase cooling. Heat fluxes of greater than 1.1 kW/cm² have been successfully rejected using two-phase cooling [5]. However, microchannel flow boiling is complex due to the strong couplings between the fluid flow dynamics, capillary and surface tension forces, and different heat transfer mechanisms [6].

Many studies have been devoted to elucidating the fundamental mechanisms in microchannel flow boiling for predicting or improving its performance under a given set of operating conditions. These studies have attempted to classify the prevailing heat transfer mechanism as nucleate boiling dominant, convective boiling dominant, or some combination of both, with comprehensive reviews given in Refs. [7] and [8]. Correlations and databases have been compiled from these studies for predicting the heat transfer coefficients under a range of fluid conditions, heat fluxes, mass fluxes, and outlet vapor qualities [8–13]. Unfortunately, because the experimental design, setup, fluid, and geometric parameters tested vary greatly in the literature, it is difficult to draw specific conclusions as to what may be a universal understanding of the dominant heat transfer mechanism in microchannel flow boiling. Many of the experimental procedures and assumptions can introduce large uncertainties into the correlations which may explain the differential results and highlight shortcomings in empirical studies. In addition to understanding heat transfer, others have explored the impact of myriad flow instabilities which diminish performance, with an extensive review given in [14]. Several studies have successfully mitigated these flow instabilities and improved flow uniformity by including pressure drop elements at the inlet of the flow channels which greatly restrict the flow area, with performance generally improving as area reduction ratio increases [15–17]. However, including these inlet restrictions also impacts heat transfer and fluid dynamics due to the sudden expansion into a larger channel. These features may result in local flow separation, recirculation zone formation, and turbulation of otherwise laminar flow through the microchannel [18].

Computational modeling of microchannel flow boiling regimes can be utilized to gain greater insights into the underlying mechanisms and behaviors, particularly for complicated devices which do not easily lend themselves to simplified analytical treatments. However, computational modeling of two-phase flows can require vast amounts of computing power and is heavily dependent on the proper boundary conditions, scope of the domain, and experimental validation for its results to be considered valuable. Nevertheless, combining an experimental study with a computational model, particularly in a complex device, may allow for a better understanding of the heat transfer and fluid dynamics at work and facilitate improved design of thermal management strategies aimed at cooling high heat flux devices.

Szczukiewicz et al. [19] studied microchannel flow boiling with and without inlet restrictions with R245fa. They found that while orifices stabilized the flow, larger orifice contraction ratios more negatively affected heat transfer coefficients at outlet vapor qualities below 0.2, while the opposite heat transfer coefficient trend was observed at greater than 0.2 outlet vapor quality. This was thought to be due to the impact of the increased fluid velocity out of the inlet orifice, differentially affecting the liquid–vapor film thickness and entrainment of liquid into the vapor at different vapor qualities. Halon et al. [20] also studied microchannel flow boiling with R245fa with and without inlet restrictions of ∼2:1 area restriction and found that heat transfer coefficient was lower with restrictions at all tested vapor qualities. This was attributed to the fact that the orifice pressure drop reduces the saturation temperature into the channels and increases the wall superheat, thus lowering heat transfer coefficient. Oudah et al. [21] studied microchannel flow boiling in water with up to five inlet orifices per channel, with area restriction ratios from 4.6% to 22.95%, and found that heat transfer coefficients with orifices initially lagged the corresponding cases without orifices until heat fluxes exceeded 80 W/cm² and mass fluxes reached 144 kg/m²s. This was attributed to two factors: increased flow inertia improved the turbulent mixing in recirculation zones which were formed just downstream of the orifices, and reduced flow instabilities at the higher heat fluxes mitigated local dryout. The results of these studies suggest that inlet orifices reduce flow instabilities but impact local heat transfer and fluid flow dynamics.

While experimental measurements can often face challenges in microscale devices, numerical simulation has its strength in resolving fine details that are impossible to observe in experiments. Most computational studies rely upon some form of computational fluid dynamics (CFD) modeling coupled with models for interface tracking and boiling/condensation mass transfer. However, finite element analysis (FEA) models have also been constructed utilizing reduced order models for pressure drop and heat transfer coefficient from the literature. Burk et al. [22] validated a comsol Multi-Physics FEA model along with literature correlations for two-phase heat transfer coefficient on a small hydraulic diameter (<100 μm) silicon microchannel evaporator subjected to large heat fluxes exceeding 1 kW/cm². It was found that conjugated heat transfer effects and heat spreading were important, and that improved two-phase heat transfer correlations are needed.

Computational fluid dynamics simulations can help develop more accurate heat transfer correlations when validated with experimental results. Some of the early efforts in CFD of microchannel boiling can be found in several review articles [23–25]. One popular method for the numerical simulation of multiphase flow is the volume of fluid (VOF) method. Gorle et al. [26] performed a validation study for VOF simulations and conjugate heat transfer of subcooled R245fa boiling in a 100 μm square silicon microchannel with an inlet orifice and a uniform applied heat flux down the channel length (test data obtained from Szczukiewicz et al. [19]); however, the orifice was neglected in the model. The multiphase mass transfer was modeled using the Lee model with a reduced-order relaxation coefficient which was tuned to match experimental data. It was found that a spatially-variable reduced order constant most appropriately fit the data, but it was emphasized that the specific values used are case-specific. Chen et al. [27] studied flow boiling of low-pressure water in a variety high aspect ratio (spanning 1:1 to 10:1 height to width, with hydraulic diameters down to 167 μm) copper microchannels using the VOF model coupled with conjugate heat transfer and the Lee mass transfer model. A uniform base heat flux of 150 W/cm² was applied down the channel length. The Lee relaxation parameter used was assumed to be constant at 0.1 but was not parametrized. The model was validated by comparing it to literature correlations. It was found that a sharp interface model resulted in more nucleate boiling-dominant behavior due to the formation of bubbly and slug flow while a dispersed model smoothed out the bubbles and resulted in more convection-dominant boiling.

While the previous studies investigated uniform applied base heat fluxes down the channel length, Lorenzini and Joshi [28] presented a numerical investigation of flow boiling with water in a 200 μm wide by 100 μm deep silicon microchannel with a nonuniform heat flux applied based on a heat flux map taken from a CPU processor. They also utilized the VOF method combined with the Lee model for phase change tracking, noting that the relaxation parameter in the source term must be tuned to experimental data but used a value of 0.1. Reasonable results compared with standard literature correlations for two-phase heat transfer coefficient and pressure drop were obtained. They found that for cases where the nonuniform heating was high near the channel outlet, large temperature gradients due to a large accumulation of vapor under bubble confinement in the channel were found, and that reversing the direction of the flow (changing the channel inlet side) greatly reduced the intensity of temperature gradients by improving the nucleate boiling heat transfer coefficients. In addition, liquid subcooling reduced temperature gradients by reducing the accumulation of vapor in the channel and subsequently reduced the channel pressure drop. These results thus indicate that nonuniform heat fluxes can greatly alter local fluid behaviors and result in exacerbated flow instabilities and variable heat transfer coefficients compared with the more commonly studied uniform channel heating.

These prior studies utilized the VOF model and compared results to either experimental data or literature predictions. Other studies have investigated the performance of different phase tracking models. Lorenzini and Joshi [29] performed CFD modeling of flow boiling in silicon microgaps and compared two different phase-tracking methods: VOF and coupled level set-VOF techniques. It was concluded that the CLSVOF method offered a sharper interface reconstruction than the standard VOF method by predicting bubble nucleation and departure mechanisms more closely to experimental observations. Soleimani et al. [30] performed a validation study on subcooled two-phase flow boiling with HFE7100 in a copper microchannel heat sink. Five different numerical models including VOF, Eulerian–Eulerian and Mixture model were compared with experimental data. VOF was recommended as the best model to simulate flow boiling in microchannels at the working conditions, as it most accurately captured the nucleate boiling behavior seen in the channels while attaining the most accurate average temperatures. From these studies, a VOF-type method offers the best resolution of the liquid–vapor interface and captures the behaviors seen in microchannel boiling most closely. In addition, the CLSVOF model seems to improve on the performance and may offer superior interface tracking capabilities.

While a wealth of knowledge on microchannel flow boiling has been assembled from experimental studies, the understanding of the heat transfer and fluid dynamics is limited by nonlocalized sensor placement, 1D assumptions, the presence of flow stabilizing features such as inlet restrictions, and nonuniform heating that complicate dynamics. Computational studies can help advance the understanding of microchannel flow boiling, especially in more realistic devices which often have unique features not often captured by fundamental research studies. Guidance from these simulations can influence future experimental design and assist in understanding how more realistic devices may operate with microchannel flow boiling. Furthermore, models can be tuned against real data to further improve their fidelity and thus improve the predictability of techniques used in designing two-phase microchannel heat sinks.

Previous CFD studies have not explored the influence of an inlet orifice on the subcooled flow boiling heat transfer in small (<100 μm) hydraulic diameter microchannel evaporators. Most studies utilized water as the coolant and used the Lee mass transfer model but kept the relaxation parameter as a constant at 0.1, with only one study using R245fa showing that refrigerants clearly require different values for this parameter. In addition, while nonuniform heat fluxes have been investigated, they have been focused on hotspot targets from CPUs in small areas, and results were not linked to experimental data. Laser diode bars and power electronics operate under relatively uniform high heat loads and high heat fluxes across a larger device area, and the effects of heat spread from a heat sink with a larger footprint than the device may have slightly different impacts than in hotspot regions on CPUs. The aim of the current investigation is to develop a 3D computational conjugate heat transfer model coupled with a CLSVOF-based multiphase model and the Lee mass transfer model. The model will be validated with experimental hardware data generated on a multimicrochannel evaporator with channel inlet restrictors and a broad area heater with nonuniform heat spreading at the channel level. The results will be used to investigate the localized fluid and heat transfer behavior in such a device that is representative of a high heat flux laser diode bar cooled with flow boiling in microchannels. Flow boiling behavior will be investigated, and heat transfer coefficients will be computed and assessed across flow regime and other typical parameters such as heat flux, mass flux, and vapor quality.

In the present study, a computational model will be developed and compared with existing experimental data. The microchannel flow boiling study was conducted over a range of test conditions using test pieces and a testing facility built and tested in previous investigations [5,31].

The test section investigated in this study was fabricated from a 300 μm wafer of silicon using standard Micro-Electro-Mechanical Systems (MEMS) microfabrication techniques. The details of the test section design and fabrication are given in Ref. [32]. The test section is shown in Fig. 1. Fluid enters and exits through vias which contain supporting silicon features that distribute and collect fluid from the channels. There are 142 parallel rectangular channels which are nominally 30 μm wide, 200 μm tall (D_h = 52 μm), and 2 mm long. Each channel has an inlet restrictor to promote flow stability and uniformity. Since the etching process etched the narrower (17.5 μm average width) inlet orifices at a slower rate than the channels, the orifices are shallower than the channels by ∼68 μm. The fluid is hermetically sealed within the device with a 0.5 mm thick borosilicate glass cap anodically bonded to the etched side of the wafer. On the back side of the device a thin film platinum heater was deposited over the center 1 mm of the channels to provide Joule heating, though the heater dimensions compared to the channel length create a nonuniform heating effect due to heat spreading in the silicon substrate.

The experimental facility used was constructed for prior experiments using similar evaporators, with full details given in Ref. [32]. Some features of the facility were modified for an investigation using transient heat loads [31]; an updated schematic is shown in Fig. 2, and a brief description follows. A positive displacement pump circulates R134a refrigerant while mass flow rate is monitored by a Coriolis meter. The degree of inlet subcooling is set using a plate heat exchanger upstream of the test section and is measured using a K-type thermocouple. The evaporator inlet pressure is set using a bladder accumulator by adjusting the compressed air side pressure and is measured by an upstream pressure transducer. The heat load is applied by varying the voltage to the test section heater with a programable power supply, and test section voltage and current are measured using a National Instruments module and shunt resistor. The test section average heater temperature is measured using an in situ calibrated infrared (IR) camera. A thin layer (∼25–50 μm) of high emissivity (∼0.9 total hemispherical emissivity) black paint is applied to the heater and facilitates more accurate IR measurements, but the effect of emissivity was calibrated out, as discussed in Sec. 4.3.3. The test section pressure drop is measured using additional pressure transmitters. Another plate heat exchanger downstream of the evaporator condenses the refrigerant. A National Instruments cDAQ-7184 running a custom labview program logs measurements at a sampling rate of 0.33 Hz.

The experimental procedure and resulting experimental data are presented and reduced here for inclusion and comparison with the computational model.

Prior to testing, a system steady-state is first established; it is attained when fluctuations are less than 2 mg/s for mass flowrate, 0.1 °C for evaporator fluid inlet temperature, and 0.6 kPa for inlet pressure. Then, the heat load for the given test is initiated as a step function and held while the test section average heater temperature stabilizes (<0.2 °C variation over a period of 30–60 s). Once the heater reaches a steady temperature, the temperatures, mass flowrate, and pressure data are then collected and averaged with the DAQ.

The inlet saturation temperature (and inlet subcooling) deviated by <3 °C, and the flow rates deviated by up to 35% from the nominal due to system perturbations at the onset of boiling in the test section.

A 30% uncertainty was assessed for the pressure drop calculated from the correlations and is propagated through accordingly in the following section.

Microchannel flow boiling is a multiphase flow problem with liquid-to-vapor phase change at the heated channel wall. Based on the earlier literature search, the VOF model was found to be more numerically stable and less computationally expensive than the Eulerian–Eulerian method. However, the major drawback of the VOF model is that it can only predict the location but not the surface normals and curvature of the interface. ansysfluent [40] provides several reconstruction schemes for calculating the interface normal vector. While the geometric reconstruction scheme using a piecewise-linear approach is considered the most accurate among them, it is still difficult to compute accurate local curvature from volume fraction because the volume fraction is not continuous near the interface. To alleviate the drawbacks of the VOF model, the level set method [41] was developed to simulate two-phase flow. The coupled level-set and VOF (CLSVOF) method was shown to be superior to either VOF or the level-set method alone as mentioned earlier. Therefore, the CLSVOF model was chosen for the current study.

The VOF model only solves a single set of momentum and energy equations regardless of number of phases. The fields for all variables and properties are shared by the phases and represent volume-averaged values. The interface between the phases is resolved by solving the continuity equation for the volume fraction of each phase. The governing equations are as follows [40]:

• Continuity Equation for volume fraction of the secondary phase (vapor)

  $∂αvρv∂t+∇·(αvρvu)=slv−svl$

  (8)
• Momentum Equation for Mixture

  $∂(ρu)∂t+∇·(ρuu→)=−∇p+∇·(μ(∇u+∇uT))+ρg+Fσ$

  (9)
• Energy Equation for Mixture

  $∂(ρE)∂t+∇·(ρu(ρE+p))=∇·(k∇T)+Se$

  (10)

The fluid saturation temperature is important in microchannel boiling heat transfer. Up to 75 kPa pressure drop was recorded during the experiments. Based on the REFPROP database [46], the saturation temperature for the worst case will decrease about 4 °C from 22.8 °C at the inlet to 18.6 °C at the outlet when the pressure drops from 623 kPa to 548 kPa. This drop in saturation temperature will impact the local heat transfer behavior as both fluid temperature and properties change. Constant saturation temperature was used as a first pass estimate based on the average pressure at the inlet and outlet. The model was run, and an updated axial pressure distribution was subsequently computed. To improve the accuracy of the simulation, a pressure-dependent saturation temperature was implemented for the high pressure-drop cases. The saturation temperature-pressure curve for R134a was determined from REFPROP, and the pressure distribution was used to calculate the local variation in saturation temperature using a UDF. Due to the computational demand of the simulations, only several iterations of this saturation temperature update were performed.

As the physical device consists of 142 parallel microchannels, only a half-channel/half-fin symmetry unit cell was modeled to save on computational demand. The symmetry boundary condition was justified from previous studies on similar evaporators which showed a small temperature variation perpendicular to the flow wise direction of the device, indicating flow and temperature uniformity due to the inlet restrictions [5,22,45]. In addition, a representative thermal image is shown in Fig. 3. The temperature variation across the channels is <2 °C away from the edges where there is some heat spread due to the thicker edge supports. However, the variation is still <4 °C and thus reflects reasonable flow uniformity across the device due to the inlet orifices. Therefore, appropriate symmetry boundary conditions were applied to the planes that divide the geometry into the half-unit cell. The heater was modeled as a surface heat flux since the thickness of the heater was ∼200 nm. The fins were considered insulated at the tips, which is reasonable because of the low thermal conductivity of the glass. Flow rate and pressure boundary conditions were set at the channel entrance and exit, as well as temperature at the inlet. The remaining boundaries were considered insulated based on the negligible heat loss mentioned in Sec. 4. The unit cell computational domain is shown in Fig. 4. A mesh sensitivity study was conducted on a simplified version of case 6 (the highest outlet vapor quality case) using 55,296, 108,528, and 145,656 high orthogonal quality hexahedral cells, with mesh refinement such that y+ < 5 near the wall. The resulting heater temperature and pressure drop were used as the comparison metrics. The results of this study are shown in Fig. 5 and demonstrate mesh convergence since the average heater temperature differs by 0.5 °C and the pressure drop differs by 0.1 kPa between the intermediate and fine meshes. Ultimately, the 145,656-cell mesh was chosen from a balance of high computational cost of running the models as well as the best agreement with the experimental data.

The ten different experimental cases outlined in Sec. 4.3.2 were simulated using ansysfluent 2021. The k–ε realizable turbulence model was added to the conservation equations to capture the turbulence, creating two additional transport equations to solve. The standard PISO algorithm was used for pressure–velocity coupling. The transient flow solver was used with adaptive time-stepping, with a global Courant number of one chosen as the time-step limiter. The resulting time steps ranged from 1.0 × 10⁻⁷ to 2.1 × 10⁻⁷ s.

The pressure drop comparison is also included for validation. It was mentioned earlier that correlations do not exist for all the conditions and complex geometries in the microchannel device studied here. In addition, many correlations have larger error bounds of 30% or more. Furthermore, at these scales, small discrepancies in channel dimensions have a large impact on pressure drop. When considering the geometric uncertainties as well as the outflow pressure drop uncertainty, which have a larger impact with higher flow rates, the computed pressure drop agrees reasonably well and is <30% difference in all cases.

The fluid dynamics and heat transfer behavior from the models will be shown using selected contours and distributions extracted from the simulation results to understand the boiling behaviors.

To understand the flow patterns, the volume fraction of vapor was examined. Figure 6 displays the contour plots of vapor volume fraction for all ten cases. Based on the parameter ranges, the cases can be grouped into three categories, which are detailed in Table 4. Four cases in group 1 have low vapor quality of less than 0.2. Boiling occurs near the leading edge of the heater as the combination of higher mass flux and lower heat flux in these cases increases the single-phase liquid region by inhibiting bubble nucleation until the fluid absorbs more heat underneath the heater footprint. The boiling initiates at the channel floor with bubbly flow which develops into slug or churn flow toward the end of the channel in cases 1, 5, and 7. Group 2 includes three cases with the lowest mass flux in the current study, and group 3 has the highest heat fluxes. Both group 2 and 3 indicate boiling occurs right after the exit of the orifice near the bottom of the channel, indicative of significant enough heat spread to initiate boiling; however, the flow patterns are different for these two groups.

Group 2 is characterized with bubbly flow at the onset of boiling which coalesces into slug flow, churn flow, or annular flow. Cases 2 and 3 display similar flow regimes of slug or churn flow since they have similar heat flux, mass flux, and vapor quality. Case 6 differs as vapor occupies a large fraction of the channel due to the lowest mass flux in the study, which may cause some localized dryout since the thin film appears to have nearly vanished from the bottom of the channel. In contrast with group 2, group 3 exhibits a separated flow with the vapor at the bottom and subcooled liquid on the top which develops bubbles and eventually coalesces with the vapor at the bottom to form churn flow. In group 3, the heat flux is increasing from case 8 to 10 and the mass flux is decreasing by case, so the vapor quality is increasing in each case. Due to this, bubbles nucleate further upstream in the separated liquid region moving from case 8 to 10. Based on these groupings and flow regimes, it appears that, in general, nucleate boiling is the dominant flow regime due to the bubbly and slug flows which extend across much of the channel. However, in case 6 convective boiling appears to be dominant as a churn or possibly annular flow comprises a large fraction of the channel length.

The liquid phase velocity streamlines can be used to understand trends in fluid dynamics within the channel. Figure 7 details these streamlines at the symmetry center plane by case. Fluid jets out of the orifice and creates a recirculation zone on the channel floor immediately downstream of the orifice, characteristic of a backward-facing step behavior due to boundary layer detachment and flow separation. This recirculation zone results in high localized velocity streamline values due to the mixing. The blank areas preclude the existence of vapor-only regions. In the higher outlet quality cases, this area is filled almost entirely with vapor and matches well with the vapor fraction contours in this region. In the higher heat flux cases, it is probable that this flow separation and eddy formation results in separated boiling flow, where the lower velocity floor at the outlet of the orifice allows boiling to commence more readily, but the increased flow inertia of the subcooled bulk liquid exiting the inlet orifice toward the channel top inhibits nucleation for longer.

The cross-sectional averaged pressure distribution down the channel length was obtained for each case. Figure 8 compares these pressure distributions. It is evident that the pressure gradient through the orifice is large, indicative of the much smaller cross section through which fluid must flow. The pressure drop through the orifice is quite steep for the higher flowrate cases and is a substantial portion of the overall pressure drop particularly for the lower vapor quality cases. In the pressure profiles a flow recirculation zone is also apparent in several of the cases due to the sudden expansion of the fluid into the channel, which acts as a pressure recovery element. The pressure recovery effect is most pronounced in the higher flowrate cases 5, 7, and 8–10, where up to roughly a third of the pressure drop in the orifice is recovered. However, for the cases where the recirculation zone is not apparent, it is possible that the vapor has pressure recovery since the streamlines shown are for the liquid only. In addition, the eddy at the bottom of the channel is due to vertical expansion out of the orifice. The fluid also expands horizontally out of the orifice, which will also result in pressure recovery that may not be apparent from this center-plane view. Cases 5 and 7 have additional pressure recovery which occurs further downstream of the orifice. These cases also have the highest flowrate, so it is likely that the flow expansion may occur for a longer duration down the channel due to increased flow inertia.

The slope of the pressure drop curves in the channel grows more negative as the vapor quality throughout the channel increases, as inertial and frictional pressure losses grow, and as the flow moves downstream of the circulation region. This transitional region occurs between 0.5 and 1.25 mm downstream of the orifice inlet and is dependent on local vapor quality and flowrate. Cases 1–4 have similar total pressure drop and general pressure distributions (except that case 4 is shifted higher up due to the higher inlet pressure) due to lower flow rates and vapor qualities. Interestingly, case 6 has low pressure drop despite having the largest outlet vapor quality; however, this is likely because case 6 also has the lowest flowrate. Figure 7 also shows relatively low liquid velocities toward the outlet of the channel for case 6 compared with cases 8–10 which have vapor qualities above 0.2 but also have much higher flow rates.

Because of the nonuniform heating distribution where the 1 mm heater is centered over the 2 mm channels, 1D assumptions typically used to calculate experimental heat transfer coefficients are not applicable. This statement is confirmed in the following figure. The axial perimeter-averaged wall heat flux contours down the channel (not including the orifice) are shown in Fig. 9. The heater lies between 0.65 and 1.65 mm, and nonuniform heating occurs outside this region due to both heat spreading and the difference between single phase and two-phase heat transfer performance. In some cases, the fluid is still subcooled partially down the length of the heater, reducing the heat fluxes in this region. The two-phase region starts with an abrupt increase in dissipated heat in these cases. Regardless, the heat fluxes are still highest under the heater footprint.

The axial heat transfer coefficients (HTCs) are plotted underneath the heater footprint in Fig. 10. The results will be cross-referenced against Fig. 6 to understand the impact of boiling regime on heat transfer coefficient.

In cases 1, 4, and 5, the fluid is still mostly or fully subcooled at the leading edge of the heater, and the heat transfer coefficients increase rapidly at the onset of boiling. Boiling initiates on the channel floor as bubbly flow, and heat transfer coefficients then slightly decrease before a period of stagnation. As the bubbly regime fills the channel, heat transfer coefficients generally tend to increase and then stagnate again. Cases 2 and 3 experience nucleation on the channel floor directly downstream of the orifice in either a slug or plug flow. The heat transfer coefficient trends for these two cases are nearly identical except that case 3's is shifted higher due to its higher heat flux. The heat transfer coefficients initially are stagnant but then increase gradually around 0.8 mm as the liquid jet from the orifice disrupts the flow and results in an abrupt transition to churn flow. The case 7 HTC resembles cases 2 and 3 but looks more like cases 1, 4, and 5 with the flow morphologies. This is likely because the boiling initiates slightly upstream of the heater on the channel floor, and the slight decline after the first onset of boiling occurs just past the leading edge of the heater. Then, the discrete bubbly regime persists before transitioning to churn flow, after which the HTC slightly increases again near the downstream edge of the heater.

The remaining cases are the highest vapor quality cases. In case 6 (which also has the lowest flowrate), because the bubbly/plug flow initiated immediately downstream of the orifice, the HTC peaks just past the leading edge of the heater. Then, an abrupt transition to high vapor fraction churn flow with some possible dryout results in a gradual decline of the HTC down the length of the heater. Cases 8–10 look similar across both HTC and flow patterns. The peak HTC shifts further upstream going from case 8 to 10 as the bubbly regime occurs earlier due to the commensurately increasing heat fluxes. In case 8, larger plugs exist further down the heater, which results in the HTC increasing at a slightly increased rate on the trailing edge of the heater in comparison with case 9 and 10.

and are plotted against heat flux, mass flux, and outlet vapor quality in Fig. 11. In general, the HTCs increase with heat flux, with the exception being from case 3 to case 4. Comparing cases 4 and 5, although the mass flux of case 5 is two times higher than case 4, the HTC increases only 15%. Similarly, the mass flux of case 7 is about three times higher than case 6 at nearly identical heat flux, but the HTC increases by ∼39%, indicating the HTC is only weakly dependent on the mass flux. Comparing cases 5 and 6, which have the maximum and minimum mass flux and vapor quality, respectively, the HTCs are nearly identical, confirming the observation that vapor quality or mass flux has a weak influence on the HTC.

Case 4 has the lowest heat transfer coefficient despite its average parameters in all categories. Case 6 has the second lowest HTC despite an average heat flux, but this is likely due to its highest vapor quality and lowest mass flux. Cases 3 and 4 have a similar vapor quality (0.21 and 0.19, respectively). Although case 4 has slightly higher heat flux and mass flux than case 3, its HTC is about 50% lower than case 3. This suggests that there may be a flow pattern change from case 3 to case 4. Recalling the vapor fraction contours in Fig. 6, it is evident that boiling initiates much later in the channel and thus reduces the heat transfer coefficient due to the larger amount of subcooling in the channel. The larger saturation temperature and subcooling in case 4 relative to the other similar cases may also play a role. As mentioned earlier, the relationship between case 4 and 5 suggests a weak dependence of HTC on mass flux. However, case 4 and 5 show slightly different boiling regimes, with bubbly and slug flow as the more dominant regimes in case 5, which, along with a higher heat flux, increases the HTC slightly relative to case 4. It is possible that the rapid change in flow regime in case 6 may lead to some localized dryout which results in the lowest average heat transfer coefficient in this investigation.

Cases 1–3 have the lowest heat fluxes in the study, while cases 8–10 have the highest; however, the HTCs of cases 1–3 are still large and comparable. Cases 2 and 3 share similar overall boiling flow behavior, subcooling, and parametric trends, while case 1 has lower vapor quality and slightly higher mass flux despite nearly identical subcooling and saturation temperature. However, case 1 has a lower heat transfer coefficient than both cases 2 and 3. Case 1 displays a different flow regime, with boiling initiating much further down the channel and with initially bubbly-dominant, then slug-dominant flow. The lower vapor quality of case 1 along with higher heat fluxes of cases 2 and 3 likely results in their higher HTC behavior. Overall, it appears that HTC increases to some extent with heat flux but is more weakly dependent on both vapor quality and mass flux. Case 5 likely has a vapor quality too small for efficient two-phase heat transfer, while cases 4 and 6 may suffer from local dryout. A combination of high heat flux, mass flux, and vapor quality result in high heat transfer coefficients as does a combination of lower heat flux, mass flux, and vapor quality. It appears that flow regime impacts results as well, with bubbly, slug, or low vapor quality churn flow resulting in a higher HTC, while higher vapor quality churn flow results in lower HTC. In addition, the flow behavior due to the orifice likely also impacts results since a large pressure drop occurs in the orifice and differentially affects the local saturation properties at different flow rates. The fluid expansion and recirculation zones just downstream of the orifice most likely play a role in HTC and flow regime as well.

In summary, this investigation created and validated a 3D conjugate heat transfer CFD model using experiments on a small hydraulic diameter microdevice under high heat fluxes. The model was validated with the experimental data and facilitated more in-depth analysis into the thermal-fluidic behavior of microchannel boiling in this evaporator. The following conclusions can be drawn:

• The inlet orifice was found to impact the flow and heat transfer dynamics downstream due to the formation of eddies and flow separation. Pressure drop and velocity streamline profiles confirmed a large pressure gradient within the orifice followed by a pressure recovery zone immediately downstream of the restrictor. Separated two-phase flow was observed in several cases.

• Two-phase heat transfer coefficients extracted from the simulations generally tended to be most sensitive to heat flux and to a lesser extent mass flux and outlet vapor quality. Bubbly, slug, and churn flows usually resulted in the highest heat transfer coefficients while dryout and high vapor volume fraction churn flow resulted in the lowest. However, saturation properties and dynamics from the orifice also likely impacted heat transfer coefficients due to the high velocity liquid jetting which disrupted downstream flow.

• Qualitatively, nucleate boiling appeared to be dominant in most of the simulation cases due to the large length of bubbly and slug flows in the channels.

Future efforts should focus on obtaining more controlled experimental results to investigate the impact of other parameters on the boiling behavior more carefully. Specifically, these results are limited to the conditions tested on this specific evaporator geometry. Higher outlet vapor qualities should be investigated to explore further trends in heat transfer coefficient and flow regimes. This computational model can then ultimately be used to develop improved predictive heat transfer correlations for small hydraulic diameter evaporators with nonuniform heat fluxes and inlet restrictions, which can guide the design of microchannel evaporators requiring aggressive peak heat fluxes with low peak temperatures. Furthermore, assessments should be made comparing the current heat transfer coefficient calculations with those from correlations in the literature. Finally, orifice design should be carefully explored to optimize flow regimes which maximize heat transfer performance and retain flow stabilizing features with the lowest penalties to pressure drop.

The authors would like to acknowledge Lawrence Livermore National Lab for the test section fabrication. This paper is dedicated in memoriam to Nicholas P. Niedbalski of the Air Force Research Lab, Dayton, Ohio, for his role in forming this collaborative effort.

• Air Force Research Laboratory (Award No. FA8650-19-F-2907; Funder ID: 10.13039/100006602).

• U.S. Office of Naval Research (Award No. N00014-18-1-2198; Funder ID: 10.13039/100000006).

D_h =

hydraulic diameter (m)

d =

distance (m)

E =

energy (J)

$F$ =

force (N)

F =

evaporation/condensation flux (kg·m⁻²·s⁻¹)

G =

channel mass flux (kg·m⁻²·s⁻¹)

$g$ =

gravitational acceleration (m·s⁻²)

H =

Height (m)

h =

heat transfer coefficient (W·m⁻²·°C⁻¹)

k =

thermal conductivity (W·m⁻¹·°C⁻¹)

L =

length (m)

L_vap =

latent heat of vaporization (J·kg⁻¹)

$m˙$ =

mass flow rate (kg·s⁻¹)

N =

number of

$n$ =

normal vector

$n̂$ =

unit normal vector

P =

pressure (Pa)

Q =

heat load (W)

$q″ =$ =

heat flux (W·m⁻²)

r =

Lee model proportionality coefficient (s⁻¹)

R_sv =

surface-volume ratio (m⁻¹)

S_e =

energy source term (J·m⁻³·s⁻¹)

s =

volumetric mass transfer rate (kg·m⁻³·s⁻¹)

T =

temperature (°C)

t =

thickness (m)

U =

uncertainty

$u$ =

velocity (m·s⁻¹)

W =

width (m)

x =

coordinate (m)

$Y¯$ =

generic dependent variable mean

Y =

generic parameter/independent variable

y+ =

non-dimensional wall distance for first boundary layer element

z =

axial variable (m)

Subscripts

a =

applied

avg =

average

base =

substrate floor between channels and heater

ch =

channel

f =

frictional

fin =

fin/wall between channels

h =

heater

i =

index in sum

in =

inlet

l =

liquid

le =

leading edge of heater

meas =

measured value

minor =

minor loss (pressure drop)

ori =

inlet orifice

out =

outlet

outflow =

downstream region after channels to pressure transducer

pred =

predicted

sat =

saturated

Si =

silicon

sim =

simulation

sub =

subcooled

tp =

two-phase

v =

vapor

w =

wall

Superscripts

* =

interface

i =

index in sum

T =

transpose

Greek Symbols

α =

volume fraction

α₀ =

void fraction

β =

accommodation coefficient

$Δ$ =

difference or change

δ =

Dirac delta function

ε =

half-thickness of the interface (m)

$κ$ =

interface curvature (m⁻¹)

μ =

dynamic viscosity (Pa·s)

ρ =

density (kg·m⁻³)

$σ$ =

surface tension (N·m⁻¹)

$θ$ =

angle (deg)

$υ$ =

specific volume (m³·kg⁻¹)

$χ$ =

thermodynamic vapor quality