Flow Boiling in Flexible Polymer Microgaps for Embedded Cooling in High-Power Applications

Remarkable advancements have been made in the last years to develop structurally flexible, multifunctional microsystems, which have also been enabled by breakthroughs in the materials science and microfabrication techniques [1,2]. The trend of increased functionality in smaller form factors (thus higher heat fluxes) in such microsystems has posed a challenge for thermal solutions, where two-phase cooling has been explored as a feasible option due to the latent heat of vaporization in stable boiling conditions. Single-phase cooling in microchannels and microgaps using DI water as the coolant has been shown to provide significantly high heat transfer coefficients [3–5]; however, the use of electrically conductive liquids is a concern for micro-electronic systems. Dielectric refrigerants are a more suitable option for applications in which the coolant is in direct contact with electrical connections, with the tradeoff of having significantly inferior thermal properties than water. The flow boiling of dielectric refrigerants has been explored in recent years as a promising alternative for embedded cooling in microsystems [6–8]; although significant advances have been made in this field, the majority of investigations have been performed in rigid microchannels/microgaps made of silicon (or metals) and therefore are not suitable for flexible microsystems. Polydimethylsiloxane (PDMS) has gained attention in the last years for microsystems fabrication due to its elasticity, transparency, biocompatibility, and its ability to create conformal contact with many types of surfaces; such combination of properties make this material an attractive option for developing flexible microcooling layers in embedded systems. Korniliou et al. [9] investigated the boiling heat transfer in a PDMS microchannel, indicating its potential as an effective thermal control method in flexible microsystems, while also reporting on challenges such as the heat transfer coefficient variation due to flow instabilities. In recent years, the enhancement of boiling heat transfer has been explored through different methods such as selectively tuning the surface properties [10–12], inducing flow mixing/disturbances [13,14], and enhancing nucleation sites and mechanisms [15,16]. The addition of wicking microstructures to microgaps operating in boiling conditions has also been investigated and proposed as a promising alternative for the thermal control of such microsystems [17]. Dai et al. [18,19] have reported that the implementation of wicking microstructures leads to a significantly larger thin-film extension of the liquid meniscus when compared to other approaches, thus resulting in a higher heat transfer coefficient.

Despite its many desirable properties, the main challenges in the adoption of PDMS-made cooling layers are its relatively low thermal conductivity (0.15 W/m-K) and its hydrophobic surface properties (contact angle of ∼ 117 deg with water) [2]. Peter et al. [20] and Chen et al. [21] have proposed the surface wettability tuning through methods such as plasma treatments and chemical modification, which have resulted in a significantly increased performance in two-phase cooling conditions.

Considering the aforementioned literature and challenges, there are many opportunities for the development of this emerging technology. In the present investigation, a flexible PDMS microgap is proposed and assessed in an effort to prove its viability for thermal management in the aforementioned applications. The proposed PDMS microgap is constituted by a wick of micropillars in an effort to enhance the local rewetting mechanism with capillary-driven flow; such approach results in the formation of thin liquid films that enhance the boiling mechanism. Noting that inefficient vapor removal (such as the stagnation of vapor slugs) is negative to the liquid supply, the proposed device integrates a number of dedicated vapor removal lanes; such pathways are proven to assist in the efficient removal of vapor from the microsystem and mitigation of early critical heat flux (CHF) issues. The dielectric refrigerant HFE-7100 is used as the working fluid under flow boiling conditions, reporting on the two-phase flow regime, heat transfer, and pressure drop. In addition to experimental results, the numerical modeling of the relevant features of flow boiling is explored with the use of a mechanistic phase-change model that is proven to accurately predict the flow variables, while also being compatible with general-purpose CFD-HT codes and constitutes a valuable tool in the analysis and design of such microsystems. The results from this study demonstrate that this approach is feasible for the removal of relatively high heat fluxes which are comparable to metallic-based or silicon microchannels, with the added advantage of structural flexibility while also providing a stable two-phase cooling mechanism.

In an effort to address the aforementioned challenges, a PDMS microgap cooling layer was designed and microfabricated for its experimental assessment under boiling conditions with the dielectric coolant HFE-7100. The PDMS material was chosen due to its elasticity, ease of microfabrication, and ability to create conformal contact with different surfaces, which are key features for developing flexible microcooling layers in embedded systems. Figure 1 depicts a sample photograph of the microfabricated PDMS cooling layer sheet with multiple zoom-in boxes for a better understanding of the micropillar array and its features/dimensions. The PDMS sheet is populated with an inline array of cylindrical micropillars with a diameter of 50 μm and 70 μm height, which is also the height of the cooling layer; the micropillar array has the same longitudinal and transversal pitches of 75 μm. The micropillar array and dimensions used for the proposed cooling layer design are similar to those used in the investigation by Cho et al. [22], resulting in a porosity of 0.57. Figure 1(a) depicts an isometric CAD view of section of the PDMS microgap with the minimum transversal extension used for its modeling with symmetry conditions. Figure 1(b) shows a scaling comparison of the PDMS sheet with a fingertip, indicating 10 mm of extension, which has 18 vapor removal lanes, each with a width of 100 μm (as indicated in the top view of Fig. 1(d)). Figure 1(c) shows a scanning electron microscope (SEM) photograph of the PDMS sheet, corroborating the correct microfabrication of proposed design and its vapor removal lane (or vapor pathways) features. Recalling that vapor accumulation is one of the main issues leading to early CHF and poor thermal performance in two-phase cooling systems, the proposed design with multiple vapor removal lanes is an effort to address such issues observed in conventional microchannel layouts and enhance the heat transfer performance.

The PDMS microgap cooling layer was microfabricated from Polydimethylsiloxane (PDMS, purchased from Dow Inc., Midland, MI), where the lithography process was used in a cleanroom facility was used to create the microstructure layout described in Fig. 1. The first step of the microfabrication process was the creation of a silicon mold through the deep reactive ion etching (DRIE) process and subsequently modified with the chemical 1H, 1H, 2H, 2H-Perfluorooctyltriethoxysilane (purchased from Sigma–Aldrich, St. Louis, MO). The PDMS was prepared with a base-elastomer to curing agent ratio of 10:1, and subsequently cast to the silicon mold with a curating time of 2 h at a temperature of 100 °C.

The microfluidic device shown in Fig. 2 was formed by bonding a portion of the flexible PDMS microfabricated sheet to a Pyrex glass substrate (purchased from University Wafer Inc., Boston, MA) treated with an oxygen plasma process. The reason for using a transparent glass substrate is to form the microgap while also allowing the visualization of the two-phase flow behavior when operating in flow boiling conditions. Heating is provided to the device through a serpentine microheater made of aluminum film and deposited on the bottom side of the glass substrate before the described bonding process. The aluminum heaters are also employed as resistance temperature detectors (RTDs), in which the electrical resistance is correlated to the average temperature; such thin film heaters were calibrated in an isothermal oven, obtaining the linear correlation between these two parameters. Figure 2(a) depicts a cross section view to better understand how the microfluidic device works, in which the inlet and outlet ports were etched from the PDMS cooling layer of 200 μm thickness and the microgap is formed by the bonding with the glass substrate of 500 μm thickness. Figure 2(b) shows a top-view photograph of the microfluidic device, indicating its dimensions and features, such as the heater geometry (serpentine) and heated area of 10 mm by 20 mm.

After the microfabrication and assembly procedures were completed, each of the contact pads on the microfluidic device was wired for their connection to the power supply and data acquisition unit. The experimental setup consists of three subsystems as shown in the schematic of Fig. 3, noting an optical imaging system, an open coolant loop, and a data acquisition system.

A pressurized nitrogen tank was connected to a stainless steel reservoir in order to pump the dielectric refrigerant HFE-7100 through the open loop, which was degassed before the experimental measurements. Pressure and temperature values at the inlet and outlet were measured using type-K thermocouple probes (TMQSS-062G, OMEGA) and solid-state pressure transducers (PX219-300A5V, OMEGA), respectively. The system flow rate was monitored and recorded using a Krhone Optimass 3300c flowmeter. The measured values for these variables were processed by a data acquisition unit (Agilent 34972 A, Keysight Technologies Inc., Santa Rosa, CA), while the electrical power to the microheater was supplied by a digital programable power supply (XLN10014, BK-PRECISION). The optical imaging system was customized by coupling a high-speed camera (Phantom V7.3, VISION Research, Inc., Charlottetown, Canada) with a digital microscope (BX-51, Olympus). The heat loss was assessed by performing single-phase energy balances, comparing the electrical power input to the heat absorbed by the liquid, resulting in less than 4% of the input power; such value was subtracted from the power supply input in order to correct the effective heat flux. More details on the heat loss calculation may be consulted in authors' previous experimental work [23]. The Kline-McClintock methodology [24] was used to estimate the experimental uncertainties; the calculated measurement uncertainties of flow rate, pressure, voltage, current, temperature, and microfabrication resolution are ±1%, ±1.5%, ±0.5%, ±0.5%, ±1 °C, and 3 μm, respectively. The experimental data were acquired after quasi-steady-state conditions were reached, which were assumed once the fluctuations in temperature in the microheater and pressure transducers were within ±0.1 °C and 0.5 kPa, respectively, during a time period of 10 min.

In addition to the experimental measurements that are key to demonstrate the technology feasibility for practical applications, the complementation with a numerical model is advantageous for a better insight into the occurring physics. Computational fluid dynamics and heat transfer (CFD-HT) modeling is valuable tool that allows the simulation of multiple variations in both geometry and operating parameters; although multiphase simulations are computationally intensive, an adequately-validated model with experimental data is one of the most cost-effective tools for the detailed analysis and optimization of two-phase microsystems.

where $vint$ is the characteristic velocity of the interface, calculated as $vint=v+mlv¨/ρ$⁠, and the term $mlv¨$ denotes the phase-change mass flux due to vaporization (see Sec. 4.2).

where $h$ is the grid spacing and $α=1.5h$⁠. With this approach, a more accurate surface tension calculation than the VOF method can be obtained, which is desirable for microscale applications where such effects are important.

Despite the surface representation advantages of the LS method over the standard VOF technique, the former has the disadvantage of not being conservative in its mass formulation, leading to calculation errors in this variable up to 20% [25]. In contrast, with the CLSVOF method developed by Sussman and Puckett [25], mass is conserved within a fraction of a percent. With this hybrid approach, both LS and VOF functions are employed for the interface reconstruction, using the volume fraction equation to get the cut size in the cell where the interface moves, while the LS function is used to get the direction of such interface [28]. The Piecewise-Linear-Interface-Calculation (PLIC) algorithm [30] is used for the interface reconstruction, using linear slopes within each cell for computing the fluid advection through the faces, and then, the distances from a given point to the front-cut segments are minimized and the LS function is reinitialized [28].

The energy for each phase (⁠$El$ and $Ev$⁠) is obtained with the corresponding specific heat and temperature, and the mixture thermal conductivity $k$ is computed similarly as Eqs. (5) and (6).

It is important to acknowledge that one of the main limitations of CFD modeling of boiling processes is the prediction and representation of active nucleation sites (which are known to randomly vary in space and time, and occur at wall crevices with entrapped vapor). Although other approaches such as heterogeneous boiling representation with adjusted nucleation temperatures are of significant interest, implementation of such in large computational domains and time-explicit calculations is a challenge beyond the scope of the present investigation. The use of the proposed model is justified since the results obtained not just in this investigation, but also in multiple complex cases in authors' previous works [5,7,8,34–37] have shown that the model results are in reasonable agreement with in-house experimental observations.

It is important to notice that the microfluidic device consists of an array of numerous, closely-spaced micropillars of small dimensions, and therefore the task of simulating the entire device of Fig. 2(b) becomes nearly impossible and not cost-effective. However, given the symmetry conditions that the array of micropillars exhibits (see Fig. 1(a)), it is possible to accurately predict the fluid and thermal behavior by just analyzing a portion of the aforementioned cooling microgap (as long as the symmetry and operating conditions are correctly implemented). Figure 4 shows the 3D computational domain consisting of a portion of the microgap structure with the required symmetry conditions, which contains an exact representation of the micropillars in their total flow length (20 mm).

Regarding boundary conditions, heating is defined as a constant heat flux at the indicated surface of the glass component in Fig. 4. The numerical formulation entails a conjugate heat transfer setup with continuity and nonslip conditions at the solid/fluid interfaces, the effect of phase change and latent heat of vaporization is included with the appropriate source term given in Eq. (18). The transversal surface cuts are modeled as periodic boundary conditions, while the remnant exterior surfaces are defined as adiabatic surfaces.

The heat flux value (q_B″) was varied in the range of 1–15 W/cm² in order to compare with the generated experimental data. The operating conditions for the numerical simulations correspond to those from experimental measurements with a device mass flow rate of 0.25 kg/h (2.315 × 10⁻⁴kg/s for the computational domain) of HFE-7100 at an inlet temperature of 18 °C. Table 1 lists the detail of thermophysical properties used to model the solid (PDMS) and liquid (HFE-7100) domains. Although most thermophysical properties have temperature dependence, it has been shown in authors' previous works [5,7,8,37] that such properties may be modeled as constant within temperature ranges relevant to micro-electronics (i.e., no significant variation) and still get a reasonably accurate prediction with a stable numerical solution.

The commercial CFD-HT software, ansys^®fluent^® (R16.2) was used as the solver in the present study the numerical solution of the presented governing equations using the finite volume method. In regards to the multiphase formulation, the built-in CLSVOF method within FLUENT was configured and the phase-change modeling was included through a user-defined-function (UDF) that accounts for the mass and energy source terms that were introduced in Section 4.2. The numerical simulations were setup as a transient, laminar flow cases (Re = 1991 for the analyzed mass flow rate) using a first-order explicit temporal discretization, whereas the spatial discretization uses a second-order upwind scheme.

The Pressure Implicit with Splitting of Operators (PISO) algorithm was used for the velocity and pressure coupling as it provides a more stable calculation than other available methods. The convergence criteria are considered when the residual values reach values of less than 10⁻⁴ for the mass, momentum, and Level-Set equations, and less than 10⁻⁷ for the energy equation. Time step adjustment is a key parameter to consider in such types of calculations, it is recommended to define an automated time-step adjustment in the range from 10⁻⁸ to 10⁻⁶ s with a maximum change factor of 2, with the objective of maintaining the Global Courant number below values of 2, being 1 the most common choice since limiting the time-step will, of course, reduce the risk of divergence.

The numerical simulations using the described method were run on a HPC cluster due to the considerable amount of geometric features resulting in a mesh of ∼30 million elements, such cluster consists of six blades of dell^®poweredge^® servers resulting in a total of 384 computing cores of amd^®epyc^® 7601 processors with a turbofrequency of 3.2 GHz, interconnected in a parallelized configuration with an infiniband^® Interconnection of 40 GB/s; the average solution time per case was in average 250 h due to the large amount of resulting cells in the computational domain coming from the many micropillar features.

The computational domain meshed with a nonconformal scheme, in which hexahedral elements are used for the fluid domain and tetrahedral elements for the solid (silicon) domain. Figure 4 shows the details of such meshing scheme with the computational domain, where boundary layers (or inflation layers) are defined around each pin fin structure. The mean element sizing in the fluid domain is 5 μm and the boundary layer (growth) ratio has a value of 1.2.

The mesh independence analysis was carried out by inspecting the effect of a range of element sizing and comparing the average bottom surface temperature for the quasi-steady-state solution. Three mesh models with different element densities were analyzed: the first with 8,020,989 cells, the second with 18,031,084 cells, and a third mesh with 30,122,791 cells. The first mesh resulted in convergence problems for the mass conservation equation because of its coarse resolution, while the results for the second and third meshes were in a reasonable agreement within 2.5% for absolute temperature. The second meshing model was selected for the present investigation in order to obtain a balance between the solution accuracy and computational time, which is known to significantly increase when reducing the element size due to the Courant number definition in explicit calculations. The selected grid density (mean fluid element sizing of 5 μm) provides a cost-effective balance in terms of solution accuracy and the total HPC time. With the simulation framework setup, the time-marching simulation of two-phase cooling requires to be initialized by loading a data solution file which contains a steady-state point with an isothermal and adiabatic flow field; more details on such implementation can be found in Ref. [37].

The mesh independence analysis was carried out by inspecting the effect of a range of element sizing and comparing the average bottom surface temperature for the quasi-steady-state solution. Three mesh models with different element densities were analyzed: the first with 8,020,989 cells, the second with 18,031,084 cells, and a third mesh with 30,122,791 cells. The first mesh resulted in convergence problems for the mass conservation equation because of its coarse resolution, while the results for the second and third meshes were in a reasonable agreement within 2.5% for absolute temperature. The second meshing model was selected for the present investigation in order to obtain a balance between the solution accuracy and computational time, which is known to significantly increase when reducing the element size due to the Courant number definition in explicit calculations.

After describing the experimental and numerical modeling methods used in the present investigation this section presents the relevant results from these approaches in order to demonstrate the feasibility of the PMDS microgap cooling layer and the CFD-HT modeling technique. Figure 5 reports key results of the present investigation, since it shows the model validation with the generated experimental data for the microfluidic device. The experimental average temperature is indicated by the red squares, while the surface-based average from the CFD-HT model is represented in blue dots. It can be noted that the results reported in Fig. 5 correspond to both single-phase and two-phase cooling conditions, this is because the fluid inlet temperature is below its saturation point and the required amount of power input to transition to boiling conditions; under the selected operating conditions, such transition occurred at ∼ 10 W/cm². As demonstrated in Fig. 5, the CFD-HT model results are in a very good agreement with the experimental data generated for the selected range of heat flux input, therefore indicating the useful capability of this model to predict the thermal performance in the relatively complex microfluidic device. It is important to point out that for the studied device, there is only one continuous heater deposited on the surface, therefore getting only one average value for resistance/temperature in the data acquisition system which is compared to the average value from the CFD simulation. Nonetheless, this numerical tool has been used in a great variety of two-phase scenarios in authors' previous works with multiple heaters operating in heterogeneous conditions [34–37] and was found to accurately predict the flow variables, therefore constituting a valuable tool in the analysis and design of two-phase cooling solutions in microsystems, including the proposed device.

One of the key hypothesis of the present investigation is that the addition of auxiliary vapor removal lanes may enhance the cooling performance and device stability by mitigating the vapor accumulation that leads to early CHF. Figure 6 shows the comparison of the predicted flow boiling mechanism between the CFD-HT model (Fig. 6(a)), and the experimental visualization (Fig. 6(b)) for a device mass flow rate of HFE-7100 of 0.25 kg/h, and a heat flux input of 15 W/cm². Consistent with authors' previous observations of flow boiling with micropillar structures [34–37], vapor bubbles begin to nucleate at the rear section of the micropillars due to the lower pressure generated by the boundary layer separation; such conditions also lead to a localized lower saturation temperature, therefore promoting the phase-change. The vapor bubbles continue to grow at the localized points until reaching conditions for their detachment from the micropillar surface, being dragged by the shear forces induced by the surrounding flow field. The proposed micropillar layout induces a suction effect that alleviates the vapor buildup by redirecting the generated bubbles to the vapor removal lanes, as it can be observed in both numerical predictions and experimental measurements, which are in excellent agreement. Because of the high density of vapor bubbles migrating to the vapor removal lanes, such structures start to coalesce and generate elongated vapor bubbles (or slugs), that are evacuated from the domain in an efficient and sustained mode. Figure 6(c) shows a velocity magnitude contour at a midplane (xy) cross section that provides additional information to better understand the vapor removal mechanism, it can be observed how the proposed layout results in the formation of several lanes that differ in magnitude, ranging from an average velocity of ∼ 6 m/s between adjacent rows of micropillars to an average of ∼ 15 m/s in the wider vapor removal lanes. The incurred velocity gradient between these sections contributes to the aforementioned vapor suction effect, while also removing such vapor structures faster than they advance in the narrower lanes. Recalling that the CLSVOF model is chosen due to its advantages in reconstructing a more physical resemblance with actual vapor interfaces than the standard VOF model, a comparison between both models in such conditions is of interest. A thorough comparison between both methods in microscale boiling conditions can be consulted in authors' previous work [35].

Figure 7 shows a photograph from the customized optical system for visualization (as described in Sec. 3), that allows a remarkable observation of the flow boiling and vapor removal mechanism, while also providing additional evidence that the proposed design and approach is relevant for microfluidic cooling applications. The darker strip zones observed in Fig. 7 are a shade projection of the aluminum microheater (deposited on the top glass surface) on the PDMS microcooling layer; such method of visualization was devised because visualization with such clarity is not possible through the PDMS because of its optical properties (partially opaque).

Having justified the selection of transparent glass as the mating surface for the PDMS microfluidic device for merely visualization purposes, it is clear that such material constitutes a highly restrictive thermal resistance to the heat transfer system due to both its relatively large thickness (500 μm) and low thermal conductivity (1.005 W/m-K), therefore leading to a significantly low thermal performance and limitation to experimentally increase the heat flux before reaching limit temperatures that would damage the microheater (>140 °C). Figure 8 shows the temperature stratification across the device with a heat flux input of 8 W/cm², where the significantly high-temperature gradient (>40 °C) across the glass thickness can be easily noted to constrain the heat transfer system.

It is important to clarify that the proposed use of this technology is one in which the PDMS microgap cooling layer is in direct contact with the heat source, therefore eliminating all un-necessary interfacial heat transfer resistances to the system, while also providing the structural flexibility features exhibited by the PDMS. Among the many potential applications for the investigated cooling technology, a notable example is the embedded cooling concept of electrical motors proposed by Joshi and Li [38], in which the slot liners of the copper windings are replaced by the PDMS microfluidic cooling layer. In this approach, the copper windings are directly cooled, thus eliminating major thermal resistances found in conventional electric motors and allowing higher motor performance by increasing torque and reducing weight. The main advantage of this cooling technology in electric motors is the structural flexibility of the cooling layer, which wraps around the windings, acts as an electrical insulator, and efficiently transfers the generated heat to the dielectric coolant. In such application, the copper windings constitute the mating surface that forms the microgap with the PDMS cooling layer, therefore minimizing the overall thermal resistance.

Figure 9 shows the temperature comparison between glass and copper as the PDMS mating surface, where it can be noted the pronounced effect on the resulting temperature due to the thermal conductivity of the compared materials. For the range of analyzed heat fluxes, the comparatively low thermal resistance of the copper/PDMS system results in device temperatures near ambient, indicating the potential of such approach for the removal of much higher high heat fluxes. Calculating the Biot numbers for both configurations, the glass/PDMS system resulted in Bi > 10 for the analyzed range of heat fluxes, while for the copper/PDMS system results yield Bi < 0.1 for the explored ranges; such results indicate the effect of thermal conductivity in the resistance to transmission by conduction in both systems.

In order to explore the copper/PDMS configuration thermal performance at more challenging conditions, the heat flux was ramped up to values of 50 and 100 W/cm² at the same mass flow rate of dielectric coolant (HFE-7100) of 0.25 kg/h. Figure 10 shows the resulting temperature distribution for such operating conditions and configuration, in which the thermal gradient across the copper surface is, as expected, negligible in comparison with the previously analyzed glass configuration, where the dominating thermal resistance is now that due to convection. The temperature scale bar in Fig. 10 is separated for both subsets in order to observe the small temperature gradients across the solid, in accordance with the previous discussion regarding Biot numbers and dominant resistances. Results indicate an outstanding capability of the proposed configuration to remove significantly high heat fluxes, which are comparable to metallic-based or silicon microchannels [5]. Even for the case in which the heat flux input is 100 W/cm², the surface temperature (101 °C) is still within acceptable ranges that also indicate a considerable potential for use in micro-electronic applications, with comparable values to those reported in Ref. [7], where a silicon microgap in two-phase cooling operation with HFE-7100 was maintained at ∼91 °C. Although such results have yet to be experimentally demonstrated in the application of electric motors, where the PDMS microgap cooling layer is indeed taking the shape of the copper windings, numerical predictions presented so far and confidence in such approach indicates the sizable potential of the proposed technology in many applications requiring a structurally flexible cooling solution.

In the present investigation, a flexible PDMS microgap is proposed and assessed in an effort to prove its viability for thermal management in applications that require microfluidic cooling and structural flexibility. The flow boiling operation with the dielectric fluid HFE-7100 was assessed in a range of heat fluxes, with the following main findings:

1. The proposed vapor removal lane layout results in a high density of vapor bubbles migrating to the vapor removal lanes, such structures start to coalesce and generate elongated vapor bubbles (or slugs), that are evacuated from the domain in an efficient and sustained mode, therefore mitigating vapor buildup.

2. The incurred velocity gradient between flow lane sections contributes to the vapor suction effect, while also removing such vapor structures faster than they advance in the narrower lanes.

3. The proposed approach can control relatively high heat fluxes which are comparable to metallic-based or silicon microchannels, with the added advantage of structural flexibility while also providing a stable two-phase cooling mechanism.

4. Applications that may significantly benefit from the use of the proposed embedded cooling technology are flexible microelectronics and high-performance electrical motors.

Besides demonstrating these practical features through the presented results, the employed CFD-HT model for the two-phase cooling analysis was compared with the thermal behavior observed in the experiments, observing a remarkable capability for predicting such with a high level of detail. The model was capable to predict with excellent agreement both thermal results and two-phase flow regimes. The experimental demonstration and validated model from the present work constitute as well a contribution for the thermal design and optimization of microgap cooling layers operating under two-phase flow conditions at the microscale.

The information, data, and work presented herein was funded by the Advanced Research Projects Agency-Energy (ARPA-E) of the U.S. Department of Energy, under the Award Number DE-AR0001023. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof. The microfluidic devices were fabricated in the cleanroom of the Institute of Electronics and Nanotechnology (IEN) of the Georgia Institute of Technology; the SEM images were taken at the IEN Microscopy Center.

• Advanced Research Projects Agency (Funder ID: 10.13039/100009224).

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

c, c_p =

specific heat (J kg⁻¹ K⁻¹)

E =

specific energy (J kg⁻¹)

$F$ =

volumetric surface tension force (N m⁻³)

$g$ =

gravitational acceleration (m s⁻²)

h_lv =

latent enthalpy of vaporization (J kg⁻¹)

k =

thermal conductivity (W m⁻¹ K⁻¹)

L =

length (m)

M =

molar mass (kg kmol⁻¹)

$n̂w$ =

wall surface normal vector

P =

pressure (Pa)

q" =

heat flux (W m⁻²)

$Sl,v$ =

mass source term due to phase change (kg m⁻³ s⁻¹)

$Se$ =

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

t =

time (s)

T =

temperature (K)

u =

velocity component in the x-direction (m s⁻¹)

v =

velocity component in the y-direction (m s⁻¹)

w =

velocity component in the z-direction (m s⁻¹)

x,y,z =

cartesian coordinate components

Greek Symbols

$α$ =

phase volume fraction

$ϕ$ =

level-set function

μ =

dynamic viscosity (kg m⁻¹ s⁻¹)

$σ$ =

surface tension coefficient (N m⁻¹)

ρ =

density (kg m⁻³)

Subscripts/Superscripts

B =

background

f =

fluid

HP =

hotspot

int =

interfacial

l =

liquid phase

Si =

silicon

SiO₂ =

silicon dioxide

v =

vapor phase

wz =

wall zone