Microfluidic interlayer cooling has been demonstrated as a practical solution for the vertical stacking of high power microelectronics. Although a considerable amount of studies has been presented for single phase cooling with this approach, the flow boiling features in more complex arrangements have not been as thoroughly studied. The embedded cooling of microelectronics is feasible with the use of dielectric refrigerants, which are ideally used in two-phase conditions in order to exploit the latent heat of vaporization. In the present investigation, the two-phase cooling in silicon microgaps is assessed under variable power and heat transfer surface area densities. The dielectric refrigerant HFE-7200 is used as the working fluid under flow boiling conditions, analyzing useful characteristics such as the two-phase flow regime, heat transfer, and pressure drop. The present investigation uses a numerical model that is capable of predicting the relevant features of flow boiling phenomena through a mechanistic phase-change model. The results from this study demonstrate that multiple hotspots with variable pin densities can be effectively controlled, with relatively uniform temperatures, under flow boiling conditions with dielectric fluids.

## I. INTRODUCTION

Microfluidic interlayer cooling in vertically stacked microelectronics is known for its capability of removing relatively large heat fluxes in small volumes, constrained by the interconnection length in such systems. Although remarkable contributions have been made in the last years for better understanding the heat transfer characteristics in microfluidic cooling applications, there are still a number of challenges to be addressed from both applied and fundamental perspectives. Single-phase cooling has been extensively studied in the last decades^{1–8} due to its inherent advantages such as stability and predictability; however, an important challenge associated with single-phase cooling is that as the mean fluid temperature increases downstream due to its sensible heat, the wall surface temperature is offset in a similar trend for fully developed conditions, leading to relatively high temperature gradients from the inlet to outlet. In addition to the relatively high heat fluxes and constrained volumes, another relevant challenge in the practical implementation of vertically stacked microelectronics is the nonuniform heat flux distribution (power map) that arises from the specific circuit architecture and load.

Recently, the concept of implementing localized surface enhancements (variable pin fins) has been analyzed as a solution to maintain a more uniform device temperature by increasing the heat transfer surface area in the flow direction.^{9,10} Rubio-Jimenez *et al.*^{9} studied the effects of variable fin density arrays with different fin geometric parameters, with results indicating a more homogeneous temperature distribution with reduced pressure drops when compared to uniform fin arrays. The dissipated heat flux for this study was 230 W/cm^{2} using water as the coolant with a pressure drop of ∼30 kPa and a maximum temperature difference across the device of 10 °C. In the investigation by Lorenzini-Gutierrez and Kandlikar,^{10} a similar concept of variable density with strip fins was assessed for different configurations in silicon microchannels with a height of 100 *µ*m. Results indicated the removal of heat fluxes up to 200 W/cm^{2} using water as the coolant with a pressure drop of 34 kPa and a maximum temperatures gradients on the silicon surface below 30 °C. In the study by Green *et al.*,^{11} a hybrid heat sink was presented in which two separate loops were used for addressing heterogeneous heat fluxes using water for hotspots zones and air for background zones. Results indicated the removal of heat fluxes up to 365 W/cm^{2} in hotspot regions and background heat fluxes of 20 W/cm^{2}, with pressure drops below 100 kPa. More recently, Wan *et al.*^{12} investigated the codesign of a 16-core architecture using a compact thermal model with liquid cooling through microgaps with pin-fins, employing a realistic localized heat flux input to the model.

The heat transfer phenomena associated with flow boiling in microsystems is not as well understood as it is for the case of single-phase cooling; significant challenges are related to flow instabilities, critical heat flux (CHF), flow regime transition, among others. Several reviews regarding the challenges and opportunities on flow boiling in microchannels^{13–17} have given some guidelines for future design and analysis of such systems. In the study by Krishnamurthy and Peles,^{18} the features of flow boiling of water in a silicon device were assessed in a 250 *µ*m high, 1.8 mm width microchannel with cylindrical pin fins of 100 *µ*m diameter. The operating conditions were varied from 20 to 350 W/cm^{2} and mass fluxes in the range of 346–794 kg/m^{2}s, resulting in wall temperatures ranging between 100 and 140 °C. Han *et al.*^{19} studied the flow boiling of water in a silicon microgap 200 *µ*m high and 10 mm width over a 1 cm^{2} heated area, using cylindrical pin fins of 90 *µ*m of diameter in a staggered array. The operating conditions were varied with heat fluxes ranging from 198 to 444 W/cm^{2} and mass fluxes from 1351 to 1784 kg/m^{2}s. Subcooled water was supplied to the device at subatmospheric pressures, with results indicating a maximum wall superheat of 60 °C at a heat flux of 444 W/cm^{2}. Although water has been the most used coolant in several studies due to its superior thermophysical properties, there are reliability concerns such as the direct contact of the coolant with electrical connections, prompting the use of dielectric refrigerants instead of water. In a recent study by Asrar *et al.*,^{20} the flow boiling of the dielectric refrigerant R245fa was investigated in a 200 *µ*m high and 10 mm width silicon microgap, with cylindrical pin fins covering an area of 1 cm^{2}. Operating conditions were ranged for heat fluxes up to 498 W/cm^{2}, mass flux values up to 7896 kg/m^{2}s, and inlet temperatures of 13 °C and 18 °C; results reported the observed two-phase flow regimes. Alternative modeling approaches for multiphase flows relevant to phase-change problems are reported in Refs. 21–23.

In the present work, the numerical results of flow boiling in a silicon microgap with variable pin fin density are presented and discussed. The dielectric refrigerant HFE-7200 is used as the working fluid, and the microgap is subjected to different heating conditions with a central hotspot. These results are obtained with a numerical model capable of predicting the relevant features of the flow boiling phenomena, with the use of a mechanistic phase-change model. The effect of having multiple hotspots will be assessed through this numerical model, indicating interesting trends and useful results for the design of two-phase cooling microsystems.

## II. COMPUTATIONAL MODEL

The Coupled Level Set Volume of Fluid (CLSVOF) method is a hybrid technique^{24} that combines the accurate spatial gradient calculation of the LS model and the volume-conservative nature of the VOF model.^{25} Such an approach was employed for the interface tracking, and the details of its implementation for the specific case of flow boiling in microsystems can be referred from previous work.^{26}

### A. Phase-change model

The modeling approach used for calculating the volumetric mass transfer due to boiling is described in this section. From the kinetic theory of gases, the interfacial mass flux for a vaporization process of a given liquid can be represented by the Hertz-Knudsen equation as

where *P*_{v} is the partial vapor pressure and *σ*_{e} is an evaporation coefficient that represents the portion of vapor molecules absorbed at the interface; Cammenga *et al.*^{27} and Maa^{28} have suggested a value of *σ*_{e} = 1 based on experimental evidence. The flow variables such as pressure and temperature may be related at the saturation condition through the Clausius-Clapeyron equation to obtain the interfacial mass flux as

The volumetric mass transfer rate $m\u0307lv$ required for the volume fraction equation can be obtained as the product of the phase-change mass flux by the interfacial area density. As suggested in the model by Lee and Verizoglu,^{29} and adapted for the CFD studies in Refs. 30–32, the constant quantities in Eq. (2) may be grouped to cast a simplified expression as

where *λ*_{l} is a relaxation parameter that may be adjusted as discussed in Ref. 32. Energy is transported with mass during the boiling process, and such an effect is accounted through the energy source term *S*_{e}, computed by multiplying the latent enthalpy of vaporization by the mass transfer rate due to phase change as

It is also important to point out that the critical heat flux (CHF) condition is not within the scope or reach of the present formulation, but rather this is a useful approach for saturated flow boiling without arriving to such conditions. The applicability and suitability for different scenarios, working fluid, and operating conditions can be referred from previous work.^{26,31–34}

### B. Computational domain

Figure 1 depicts the geometrical features of the analyzed portion of a silicon microgap, inspired by the microgap with a single hotspot that has been studied both numerically and experimentally in previous work.^{33,34} For the present investigation, the number of hotspots has been increased to three, which have been located at the inlet, center, and outlet of the computational domain, as indicated in Fig. 1. The main reason for choosing such a configuration is the assessment of multiple heterogeneous heat sources and their effect on the heat transfer, pressure drop, and two-phase flow regimes at the microscale, giving a more realistic insight for microelectronic systems that are known to operate in such conditions. The modeled portion of the silicon microgap has a height of 200 *µ*m height with a substrate thickness of 50 *µ*m. The domain is populated with cylindrical pin fins that connect the bottom with the top of the microgap (200 *µ*m height) classified into two different zones; the background array has pin fins of 150 *µ*m diameter, having equal longitudinal and transversal pitches of 225 *µ*m, while the hotspot array has pin fins with both the diameter and pitches reduced to half, 75 *µ*m and 112.5 *µ*m, respectively. The active heated length of the domain is 10 mm in extension, with a width of 500 *µ*m.

### C. Boundary conditions

The computational domain is subjected to heterogeneous heat flux conditions in the different background and hotspot zones at the bottom of the silicon surface, as indicated in Fig. 2. The side walls of the computational domain are modeled as periodic boundaries, while all of the remaining exterior surfaces are prescribed with adiabatic boundary conditions. The background heating zones are defined with a heat flux value of *q*_{B}″ = 30 W/cm^{2} (taken from the experimental values in Ref. 26), while the multiple hotspots are subjected to a range of heat fluxes (*q*_{HS}″) from 60 to 120 W/cm^{2}, a ratio of 2–4 times the background heat flux value. Two different cases are analyzed, in which the magnitude of the hotspot heat flux is varied according to the values shown in Fig. 2. The flow inlet is prescribed with a mass flow rate of 1.099 × 10^{−7} kg/s, corresponding to the 60 ml/min of the dielectric refrigerant HFE-7200 used in the experimental tests in Ref. 26 (entire domain). The inlet fluid temperature is 81 °C, which corresponds to the saturation temperature of HFE-7200 observed in the experimental measurements at an inlet pressure of 365 kPa.^{26}

The solid domain (silicon) is modeled with a temperature-dependent thermal conductivity, using the values reported in Ref. 35 with the following polynomial form (R^{2} = 0.99):

where the units for temperature are defined in Kelvin (K), with a range of applicability between 280 and 373 K. The thermophysical properties of silicon and HFE-7200 used for the present modeling are listed in Table I.

Property . | Symbol . | Silicon . | HFE-7200 (liquid) . | HFE-7200 (vapor) . |
---|---|---|---|---|

Density (kg/m^{3}) | ρ_{si}, ρ_{l}, ρ_{v} | 2330 | 1 325 | 10.13 |

Thermal conductivity (W/m-K) | k_{si}, k_{l}, k_{v} | Equation (5) | 0.079 66 | 0.014 36 |

Specific heat (J/kg-K) | c, c_{pl}, c_{pv} | 712 | 1 333 | 972.2 |

Dynamic viscosity (Pa-s) | −, μ_{l}, μ_{v} | … | 3.79 × 10^{−4} | 1.035 × 10^{−5} |

Surface tension coefficient (N/m) | σ | … | 0.012 99 |

Property . | Symbol . | Silicon . | HFE-7200 (liquid) . | HFE-7200 (vapor) . |
---|---|---|---|---|

Density (kg/m^{3}) | ρ_{si}, ρ_{l}, ρ_{v} | 2330 | 1 325 | 10.13 |

Thermal conductivity (W/m-K) | k_{si}, k_{l}, k_{v} | Equation (5) | 0.079 66 | 0.014 36 |

Specific heat (J/kg-K) | c, c_{pl}, c_{pv} | 712 | 1 333 | 972.2 |

Dynamic viscosity (Pa-s) | −, μ_{l}, μ_{v} | … | 3.79 × 10^{−4} | 1.035 × 10^{−5} |

Surface tension coefficient (N/m) | σ | … | 0.012 99 |

### D. Numerical procedure and mesh

The computational fluid dynamics and heat transfer (CFD-HT) solver software ANSYS® FLUENT® R15.0 was used for the numerical solution of the governing equations through the finite volume method; the built-in multiphase CLSVOF model in this package was used and coupled with the described phase-change model through the implementation of a user-defined-function (UDF) for the mass and energy source terms discussed in Sec. II C. The simulation is set up as a transient, conjugate heat transfer problem with laminar fluid flow, using a first-order explicit temporal discretization and a second order upwind spatial discretization.

The Pressure Implicit with Splitting of Operators (PISO) algorithm was selected for coupling the velocity and pressure fields, and the convergence criteria for each time step was defined for residual values reaching lower than 10^{−4} for the mass, momentum, and level-set equations and lower than 10^{−7} for the energy equation. The time step is automatically updated in the solution with values ranging from 10^{−8} to 10^{−6} s in order to keep the global Courant number below 1.0 for numerical stability. The simulation of the analyzed cases was performed using a High Performance Computing (HPC) cluster consisting of 80 compute cores in Intel® Xeon® E5-2680 v2 processors at 2.8 GHz scaled in a parallel configuration with a 4× QRD InfiniBand® Interconnection at 40 GB/s; the average solution time was approximately 100 h per case.

The computational domain was meshed using the proposed nonconformal scheme in Ref. 33, where hexahedral elements are used for the fluid domain and tetrahedral elements for the solid (silicon) domain, providing a better convergence and lower time solutions with respect to all tetrahedral meshes. The mesh sensibility was inspected in the computational domain by comparing the simulation outputs at different grid densities; the field variables used for comparison of results are the average bottom surface temperature and the two-phase pressure drop at the quasi-steady state solution. The identified, time-efficient meshing model used for the present analysis has a mean element sizing of 5 *µ*m, using the same type of boundary layer (inflation method) recommended in Refs. 31–34 with five layers around fluid/solid interfaces at a growth rate of 1.2. Such meshing combination provides a good balance between accuracy and the solution time, the latter being known to considerably increase when using smaller cell sizing due to the explicit nature of the calculation and the way the interface is advected in the CLSVOF method.

## III. RESULTS

It is important to point out that the described computational model has been successfully validated in previous work,^{26} where such type of approach was directly compared to experimental results corresponding to the silicon microgap with a central hotspot using the same dimensions of pin fins for both background and hotspot zones. The present work adds to investigating the effect of having not one but multiple hotspots of different magnitude, which is relevant to more closely approach the case of operational microelectronics; such results can provide useful insights regarding design considerations for two-phase cooling heterogeneous microsystems.

### A. Transient temperature behavior

One of the important advantages of having a validated numerical model is that multiple cases can be analyzed in a great level of detail, gaining better insights regarding the temporal and spatial behavior of the flow variables. Figure 3 depicts the time evolution of the temperature field (measured at the device’s base) for the two analyzed cases in which the position and magnitudes of the hotspot heat fluxes are varied. Recalling that the simulation is started at the specified saturation conditions of HFE-7200, it can be observed how the wall superheat starts elevating in each of the hotspot zones. As expected, the hotspot zones with highest heat fluxes start to spread out faster to the adjacent background zones, but it is also important to notice how the local temperatures at each hotspot are very similar regardless of the location as the time evolves for the first 10 ms of two-phase flow. One of the known advantages of two-phase cooling is that system temperatures are considerably more uniform than single phase cooling approaches due to the latent heat of vaporization; for the simulated cases, it can be noted how the background zones are maintained at a relatively uniform temperature throughout the entire time evolution.

For the present transient simulations, the temperature fluctuations from one time step to another become asymptotic to quasi-steady state values, in which the flow field variables are stabilized. As indicated in Fig. 3, for the present cases with 20 ms of flow simulation, such values are reached, in which the results from 17 ms to 20 ms and onwards are virtually the same and indicate the maximum values at which such configuration would operate. Comparing case 1 with case 2 at the quasi-steady conditions, it can be noted in Fig. 3 how the hotspot location (inlet, central, or outlet) does have an effect on the local temperatures. In specific for the highest local heat flux of 120 W/cm^{2}, the local temperature is slightly affected by the location; for case 1 (120 W/cm^{2} at the inlet hotspot), the average wall superheat is ∼16 °C, whereas for case 2 (120 W/cm^{2} at the central hotspot), the average wall superheat is ∼18 °C. For the case of the hotspot magnitude of 90 W/cm^{2}, it is interesting to observe that the effect of position is significantly smaller between both cases, with an average wall superheat of ∼15 °C for case 1 (90 W/cm^{2} at the central location) and ∼15.5 °C for case 2 (90 W/cm^{2} at the outlet location). Such results are mainly attributed to the effect local heat spreading, with the adjacent hotspots to inlet and outlet zones being favored due to larger planar extensions for spreading.

In order to better quantify the temperature distributions at the quasi-steady state, Fig. 4 presents a plot of the spatial variation of the average wall superheat across the 10 mm of the active heated length; the dashed vertical lines denote the location of the inlet, central, and outlet hotspots (1, 2, and 3, respectively). The aforementioned trends can be detected in this plot, such as the increased peak temperatures caused by moving the heat fluxes with higher magnitudes toward the outlet. Although the peak wall superheats for all of the hotspots are comparable in value with offsets within 2 °C, the trend is consistent and therefore suggests that lower peak temperatures may be obtained by directing the flow inlet to the areas with higher heat fluxes, with case 1 having better thermal results than case 2 even though the effective power input and other operating conditions are the same for both cases. Table II lists the numerical averages obtained for both cases, in which the local values are reported indicating the specific heat flux at each zone, providing further evidence on the effects of location and heat spreading in the quasi-steady state temperature distributions across the device.

. | Wall superheat (°C) @ heat flux (W/cm^{2})
. | |||
---|---|---|---|---|

. | Background . | Inlet hotspot (HS1) . | Central hotspot (HS2) . | Outlet hotspot (HS3) . |

Case 1 | 12.93 °C (30 W/cm^{2}) | 16.07 °C (120 W/cm^{2}) | 15.28 °C (90 W/cm^{2}) | 11.71 °C (60 W/cm^{2}) |

Case 2 | 12.38 °C (30 W/cm^{2}) | 9.07 °C (60 W/cm^{2}) | 18.20 °C (120 W/cm^{2}) | 15.19 °C (90 W/cm^{2}) |

. | Wall superheat (°C) @ heat flux (W/cm^{2})
. | |||
---|---|---|---|---|

. | Background . | Inlet hotspot (HS1) . | Central hotspot (HS2) . | Outlet hotspot (HS3) . |

Case 1 | 12.93 °C (30 W/cm^{2}) | 16.07 °C (120 W/cm^{2}) | 15.28 °C (90 W/cm^{2}) | 11.71 °C (60 W/cm^{2}) |

Case 2 | 12.38 °C (30 W/cm^{2}) | 9.07 °C (60 W/cm^{2}) | 18.20 °C (120 W/cm^{2}) | 15.19 °C (90 W/cm^{2}) |

Despite the minor variations detected in peak temperatures across different zones in the computational domain, another feature to be pointed out is the capability of the device and pin fin arrangement to handle relatively high ratios of hotspot to background heat fluxes (up to 4:1) in multiple locations and maintain relatively low wall superheats. Unlike previous investigations in which the hotspot was limited to a central zone, for the present study the simulation results encourage the possibility of selectively tuning the heat transfer surface area depending on the localized power load and circuit architectures.

Figure 5 depicts the variation of the heat transfer coefficient (calculated as the local heat flux divided by the local wall superheat) along the microgap length for the different analyzed cases with heterogeneous heating conditions. The variable behavior is associated with the nonuniform power distribution and resulting wall temperatures, having the highest values at the hotspot areas.

### B. Transient two-phase flow regimes

Figure 6 depicts the transient fluctuations of the vapor-phase distribution inside the domain for 20 ms of flow. The depictions correspond to a top view of the domain, where the vapor-liquid interface is denoted by the light blue surfaces, with the flow occurring from the left to right and the hotspot heat fluxes being indicated for each analyzed case. Recalling that the simulations are started with the steady-state solution of the single-phase flow field at the saturation conditions, during the first millisecond (ms) of the simulation, the temperature starts to slightly increase in the hotspot zones (as previously discussed in Fig. 4), causing the nearby fluid zones to start generating vapor structures. For the snapshot corresponding to *t* = 3 ms, a more noticeable amount of vapor bubbles can be observed at the hotspot zones, with larger amounts spotted at the zones where the local heat flux is the highest (120 W/cm^{2}). For the subsequent time step of *t* = 5 ms, the vapor phase distribution is more evolved and the background zones start to show the presence of grown bubbles at the flow stagnation points, which have been experimentally proven to be the preferred zones for bubble nucleation and growth^{26} and are in consistency with the modeling predictions. For *t* = 7 ms, a number of detached vapor structures may be observed flowing downstream and merging, leading to the formation of a churn-type two-phase flow regime past the central hotspot zone. After 9 ms of two-phase flow, the vapor distributions for both cases look virtually the same and the two-phase flow regimes look significantly evolved; after this point in time, the subsequent time steps look very similar in a qualitative visual inspection, which is quantified in Sec. III C by looking at the two-phase pressure drop, the flow variable that is directly affected by the vapor phase distribution. The similar amounts of vapor and their general distribution in the domain for time steps beyond ∼10 ms indicate that a balance has been achieved in terms of evacuated and generated vapor due to boiling, which is also an indication of the so-called quasi-steady state.

### C. Two-phase pressure drop

The two-phase pressure drop is a relevant flow variable that is closely related to the two-phase flow regime and its evolution; for such cooling approach, the pressure gradient for transporting a liquid/vapor mixture is considerably larger than that required for a single phase. Although some experimental correlations are available in the literature to estimate the two-phase pressure drop, these are significantly limited to a small number of conditions and simplified geometries; for the present cases with a relatively complex distribution of heterogeneous geometries and distributions, the best option is the use of a numerical model with proven predictions capabilities for experimental conditions.^{26} Figure 7 depicts the transient two-phase pressure drop behavior for the two analyzed cases, where the initial pressure drop of ∼35 kPa is the same for both cases as the geometries are exactly the same and corresponds to the value for single-phase, steady state conditions at the specified flow rate. It can be observed how the pressure drop steeply increases for both cases during the first 6–7 ms of two-phase flow due to the gradual increase in vapor structures inside the microgap and their buildup. The plotted lines are consistent with the two-phase flow regime depictions in Fig. 6, observing a peak around 8 ms of flow that is associated with the point of maximum vapor buildup inside the domain right before the continuous removal and generation processes reach a balance. For subsequent times beyond *t* = 10 ms, the trend is very similar for both cases with fluctuations around a mean value of ∼62 kPa (horizontal dashed line).

## IV. CONCLUSIONS

The present investigation numerically assessed the two-phase cooling features in a low-profile silicon microgap with a variable density of heat transfer surface area, which was subjected to heterogeneous heating in an effort to more closely approach the operation of microelectronic systems. The flow boiling of the dielectric fluid HFE-7200 was studied under saturated conditions, and the analysis of the results leads to the following main remarks:

Under two-phase cooling conditions, the background zones with constant heat flux conditions can be kept at a relatively uniform temperature due to the fluid’s latent heat of vaporization, regardless of the location of the multiple hotspots and their magnitude for the analyzed cases.

There is a slight effect of the location of the hotspot and its heat flux magnitude on the local wall superheat; such results are attributed to the heat spreading effect as the inlet and outlet hotspots are more favored to a sink effect due to the adjacent solid zones with no power input, whereas for the case of the central hotspot, the adjacent zones are subjected to the background heat flux. For such a condition, it is therefore suggested to place the hotspots with higher heat fluxes either at the inlet or exit and not at the central location.

The effects of the hotspot location and magnitudes are negligible on the two-phase pressure drop characteristics as long as the power input and geometries are kept the same.

The present results may be useful for the thermal design of two-phase microcooling systems under heterogeneous operating conditions as these have proven that multiple locations and heat flux magnitudes can be effectively addressed under relatively low wall superheats using the concept of variable pin fin clustering.

## ACKNOWLEDGMENTS

The authors acknowledge financial support for this work from the DARPA ICECool Fundamentals program.

## NOMENCLATURE

*c*,*c*_{p}specific heat (J kg

^{−1}K^{−1})*E*energy (J kg

^{−1})- $F\u2192$
volumetric surface tension force (N m

^{−3})- $g\u2192$
gravitational acceleration (m s

^{−2})*G*mass flux (kg m

^{−2}s^{−1})*h*_{lv}latent enthalpy of vaporization (J kg

^{−1})*k*thermal conductivity (W m

^{−1}K^{−1})- $m\u0307lv$
mass source term due to phase change (kg m

^{−3}s^{−1})*P*pressure (Pa)

*q*″wall heat flux (W m

^{−2})*S*_{e}energy source term (J m

^{−3}s^{−1})*t*time (s)

*T*temperature (K)

- $v\u2192$
velocity vector (m s

^{−1})