Flow boiling heat transfer in silicon microgaps with multiple hotspots and variable pin fin clustering

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.⁹ 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² 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,¹⁰ 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² 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.,¹¹ 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² in hotspot regions and background heat fluxes of 20 W/cm², with pressure drops below 100 kPa. More recently, Wan et al.¹² 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,¹⁸ 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² and mass fluxes in the range of 346–794 kg/m²s, resulting in wall temperatures ranging between 100 and 140 °C. Han et al.¹⁹ studied the flow boiling of water in a silicon microgap 200 µm high and 10 mm width over a 1 cm² 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² and mass fluxes from 1351 to 1784 kg/m²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². 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.,²⁰ 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². Operating conditions were ranged for heat fluxes up to 498 W/cm², mass flux values up to 7896 kg/m²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.

The Coupled Level Set Volume of Fluid (CLSVOF) method is a hybrid technique²⁴ that combines the accurate spatial gradient calculation of the LS model and the volume-conservative nature of the VOF model.²⁵ 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.²⁶ 

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.²⁷ and Maa²⁸ 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 $ṁlv$ 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,²⁹ 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

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.

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² (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², 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⁻⁷ 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.²⁶ 

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² = 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.

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⁻⁴ for the mass, momentum, and level-set equations and lower than 10⁻⁷ for the energy equation. The time step is automatically updated in the solution with values ranging from 10⁻⁸ to 10⁻⁶ 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.

It is important to point out that the described computational model has been successfully validated in previous work,²⁶ 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.

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², the local temperature is slightly affected by the location; for case 1 (120 W/cm² at the inlet hotspot), the average wall superheat is ∼16 °C, whereas for case 2 (120 W/cm² at the central hotspot), the average wall superheat is ∼18 °C. For the case of the hotspot magnitude of 90 W/cm², 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² at the central location) and ∼15.5 °C for case 2 (90 W/cm² 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.

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.

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²). 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²⁶ 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.

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.²⁶ 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).

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:

1. 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.

2. 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.

3. 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.

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

c, c_p

specific heat (J kg⁻¹ K⁻¹)

E

energy (J kg⁻¹)

$F→$

volumetric surface tension force (N m⁻³)

$g→$

gravitational acceleration (m s⁻²)

G

mass flux (kg m⁻² s⁻¹)

h_lv

latent enthalpy of vaporization (J kg⁻¹)

k

thermal conductivity (W m⁻¹ K⁻¹)

$ṁlv$

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

P

pressure (Pa)

q″

wall heat flux (W m⁻²)

S_e

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

t

time (s)

T

temperature (K)

$v→$

velocity vector (m s⁻¹)

Greek letters

α

volume fraction

μ

dynamic viscosity (kg m⁻¹ s⁻¹)

ρ

density (kg m⁻³)

σ

surface tension coefficient (N m⁻¹)

λ

relaxation time for mass transfer between phases (s⁻¹)

Subscripts/superscripts

int

interface

l

liquid phase

si

silicon

sat

saturation

v

vapor phase