Comparison of the Volume of Fluid and CLSVOF Methods for the Assessment of Flow Boiling in Silicon Microgaps

The increasing power density and miniaturization of electronic systems of recent years has motivated the research of novel and improved thermal management technologies at the microscale. A clear example is the development of functional 3D ICs, which offer several advantages over planar architectures, such as significant increases in speed of signal transmission through the use of reduced-length vertical interconnections, commonly referred to as through-silicon-vias (TSVs) [1]. One of the main challenges in the integration of such technology is the small volume that is available for heat removal, constrained by the length of the TSVs, which is typically lower than 300 μm. Since the use of conventional heat-sink hardware is not suitable for this type of devices, new methods and technologies have been sought in recent years. The use of microfluidic cooling in interlayer arrangements has been proposed as an effective alternative to achieve relatively large heat transfer coefficients through either single-phase [2–5] or two-phase cooling [6–8].

Krishnamurthy and Peles [6] investigated the features of flow boiling of water in a silicon device with a flow passage 250 μm high and 1.8 mm width, populated with cylindrical pin fins of 100 μm diameter. The device was heated at the bottom with heat fluxes ranging from 20 to 350 W/cm² and mass fluxes in the range of 346–794 kg/m²s. The reported wall temperatures ranged between 100 °C and 140 °C, due to the saturation properties of water above atmospheric pressures. More recently, Han et al. [7] experimentally investigated the flow boiling of water in a silicon microgap 200 μm high and 10 mm width, with cylindrical pin fins of 90 μm diameter in a staggered array. The device was tested for heat fluxes ranging from 198 to 444 W/cm², and mass fluxes from 1351 to 1784 kg/m²s. In an effort to reduce the wall temperatures, subcooled water was supplied to the device at subatmospheric pressures. Results indicated a maximum wall superheat of 60 °C at the heat flux of 444 W/cm². Water has been the preferred coolant in several studies due to its superior thermophysical properties, such as high thermal conductivity and latent heat of vaporization, when compared with other fluids. However, reliability concerns, such as the direct contact of the coolant with the device in interlayer configurations, have prompted the use of dielectric refrigerants instead of water. Asrar et al. [8] experimentally studied the flow boiling of the dielectric refrigerant R245fa in a silicon microgap 200 μm high and 10 mm width, with pin fins of hydrofoil cross-sectional shape. The cooling performance was assessed for heat fluxes up to 326 W/cm² and maximum outlet quality of ∼13%, while also discussing on the observed two-phase flow regimes.

Despite the encouraging experimental evidence that suggests two-phase cooling as a potential candidate for the 3D-stacking of ICs, the detailed thermal modeling of such has lagged behind due to the complexity of the involved phenomena. Although this type of modeling is computationally intensive, the development and refinement of new and current models can eventually lead to more cost-effective tools for optimizing and analyzing such systems; considering the expensive, time-consuming, and complex processes related to the microfabrication and testing of new prototypes. Few investigations in this area have been reported in recent years. Wu et al. [9] and Yang et al. [10] used computational fluid dynamics (CFD) tools to model the flow boiling of the refrigerant R141B in an evaporator coiled tube by coupling the VOF technique with a mechanistic phase-change model, originally proposed by Lee [11]. Their results indicated a good capability to predict experimentally observed two-phase flow regimes and their evolution, encouraging the use of this type of modeling. Fang et al. [12] presented a CFD study of the flow boiling of water in a microchannel formed by three solid walls and a vapor-venting membrane on top, with the main objective of reducing the accumulation of vapor and related flow instabilities; the employed model used a similar phase-change model to that in Refs. [9] and [10]. More recently, Magnini and Thome [13] presented the numerical simulation of saturated flow boiling in a circular microchannel by using the VOF method with extra functions to improve the surface tension calculation and to include the phase change, based on the kinetic theory of gases. Recent modeling efforts by Lorenzini and Joshi [14,15] have focused in the transient analysis of flow boiling for cooling layers in 3D ICs applications, where the VOF method has been coupled with a mechanistic phase-change model based on the same principles as those used in the listed literature [9–13], while also solving the heat conduction problem in the silicon domain.

The VOF model has been preferred over other alternative techniques (for this type of analysis) due to its conservative formulation in which the interface is advected without violating the mass conservation. However, the main disadvantage of the VOF model is the resolution of the interface when compared with other available models. Other alternatives for a sharper interface reconstruction are the level-set (LS) and phase-field (PF) methods; however, they present the disadvantage of being not-conservative in mass. A more detailed discussion regarding the features of these methods may be consulted in Refs. [16] and [17]. In an effort to reduce the associated issues with the available models for the analysis of multiphase flows, Sussman and Puckett [18] proposed and developed an algorithm more accurate than either the LS or the VOF techniques; the so-called CLSVOF method. It was shown that the CLSVOF could calculate sharp interfaces as the LS method, with the important advantage of conserving mass within a fraction of a percent, while the standard LS method could lose up to 20% of the mass.

In the present investigation, the VOF and CLSVOF are compared for the analysis of flow boiling in a silicon microgap with cylindrical pin fins for 3D ICs applications. These interface calculation techniques will be coupled with a mechanistic phase-change model for the mass and energy transfer in order to account for the saturated boiling of the dielectric refrigerant R245fa, and the implications of implementing both approaches will be assessed and discussed. The transient flow regime evolution and quasi-steady state solutions for the specified applications will be presented for both methods, and will be compared with reported experimental data regarding bubble nucleation and departure from this type of microstructures.

where the term $ρ/((ρl+ρv)/2)$ is used to refine the capability of the CSF model for high-density ratio interfaces (e.g., liquid refrigerant and its vapor); this method is referred as density scaling of the CSF, and helps to keep a constant interface thickness for tracking [20]. The interface is reconstructed by using the piecewise-linear-interface-calculation (PLIC) algorithm developed by Youngs [21], where the interface between phases is represented as a linear slope within each element for calculating the fluid advection through the cell faces.

where $vint$ is the characteristic velocity of the interface, calculated as $vint=v+mlv″/ρ$⁠. In the case of no phase change, the velocity of the interface therefore corresponds to the velocity field. The term $mlv″$ represents the phase-change mass flux between phases due to boiling, and its calculation is discussed in Sec. 2.3.

In order to reduce the mass loss problem of the standard LS method, both the LS and VOF function values are used for the interface reconstruction, where the solution of the volume fraction equation is used to get the size of the cut in the cell where the interface is advected, and the gradient of the LS function provides the direction of the interface [22]. The PLIC algorithm [21] is used for the interface reconstruction, and after this step is completed, all of the possible distances from a given point to the front-cut segments are minimized to reinitialize the LS function [22].

where $λl$ is referred to as the relaxation time, which is a parameter that may be adjusted in order to match experimental data. A reasonable agreement with experimental correlations for saturated flow boiling of water in microchannels was reported in Ref. [14], and the effect of varying this parameter on the flow boiling regimes of R245fa in a silicon microgap was discussed in Ref. [15].

The variation of the saturation temperature as a function of the local pressure is enabled in the present model through the available thermodynamic data for the refrigerant R245fa from the EES^® software [25].

Recalling that the ultimate application is the thermal control for the 3D integration of ICs, the computational domain is based on the dimensions for the interlayer cooling of such devices. The standard dimensions of planar dies for this application are 10 mm by 10 mm, which set the available footprint area in which surface enhancements are usually incorporated. Since the available volume for heat removal is constrained by the length of the vertical interconnections (TSVs), the use of closely spaced pin fins is an effective alternative since it provides a significant increase of heat transfer area, while also a medium through which TSVs cross. Considering the complexity and transient nature for modeling flow boiling in such application, the simulation of the entire domain and its number of surface features results computationally restrictive, even using high-performance computing, and therefore, for the present investigation just a portion of the silicon microgap is considered.

Figure 1 depicts the analyzed domain consisting of a microgap populated with cylindrical pin fins of 150 μm diameter in a staggered array with equal longitudinal and transversal pitches of 225 μm; these pin fins are distributed on a length of 10 mm, and a width of 2 mm. Adiabatic entrance and exit sections of 0.6 mm each are also defined in the computational domain. The pin fin and microgap height are equal, with a value of 175 μm and the substrate thickness is defined as 50 μm, which are relevant dimensions to this application. Although just the silicon domain is depicted in Fig. 1, the fluid domain is known to fill out the microgap.

Regarding boundary conditions for the present modeling, heating in the form of a constant heat flux of q″ = 100 W/cm² is defined at the bottom surface of the computational domain in the area that is populated by the cylindrical pin fins only, while the bottom surfaces below the inlet and exit sections and top wall of the microgap are prescribed as adiabatic surfaces. The lateral surfaces of the computational domain are defined as symmetry boundary conditions, considering that ten rows of pin fins may provide useful insight for the analysis of the two-phase flow regime evolution, and the objective of the present study is comparing the interfacial resolution and implications of using the VOF and CLSVOF methods for this application. The flow inlet is specified with a constant mass flux condition using the dielectric refrigerant R245fa with a value of G = 1000 kg/m² s at a saturation temperature of T_sat = 30 °C. The fin structures located in the exit section serve the purpose of splitting the flow in different outlets in order to mitigate flow reversal issues in the simulation and hence, help its convergence; these outlet faces are specified as pressure outlets. The internal walls where the fluid and solid domains interact are prescribed with no-slip boundary conditions for the liquid phase, while also specifying a temperature continuity condition at the solid/fluid interface.

with the units for temperature are defined in Kelvin (K), and the range for applicability is 280–373 K. The thermophysical properties of silicon and R245fa used for the present modeling are listed in Table 1.

The CFD solver ANSYS^®FLUENT^® 15.0 was used for the numerical solution of the governing equations through the finite volume method. The built-in multiphase models for the VOF and also CLSVOF techniques in this package were used and coupled with the described phase-change model through the development of a user-defined-function (UDF) for mass and energy source terms. The problem is treated as transient laminar with a first-order explicit temporal discretization, and a second order upwind spatial discretization.

The coupling between pressure and velocity fields is done through the use of the pressure implicit with splitting of operators (PISO) algorithm, and the solution for each time step was considered converged when the residual values reached to a value of less than 10⁻⁴ for the mass and momentum equations, and less than 10⁻⁷ for the energy equation. The time step was automatically adjusted with values in the range from 10⁻⁸ to 10⁻⁶ s in order to maintain the global Courant number below a value of 1.0 for stability.

The transient simulation was run through high performance computing (HPC), in which a total of 80 cores of Intel^® Xeon^® E5-2680 v2 processors at 2.8 GHz are used through a parallel configuration with a 4× QRD InfiniBand^® Interconnection at 40 GB/s.

The microgap cooling layer was meshed through the use of a hybrid method where hexahedral elements are used for the fluid domain and tetrahedral elements for the solid (silicon) domain. The number of solid elements is reduced through a nonconformal meshing scheme, which is acceptable due to relatively low temperature gradients in this type of zone when compared with those arising in the fluid domain. Boundary layer meshing (five inflation layers with a growth rate factor of GR = 1.2) is defined at every fluid/solid interface, where the steep gradients in both velocity and temperature require more resolution. The parallel meshing of the computational domain indicated in Fig. 1, under this specified settings, was achieved in ∼2 h for a model with ∼2.5 × 10⁶ elements using the ANSYS^® meshing software with four cores of an Intel^® Xeon^® E5-1507 processor, consuming ∼4 GB of RAM memory. Figure 2 shows a zoomed-in top view picture of the mesh details for the solid (pin fins) and fluid domains.

A mesh sensibility analysis was conducted to ensure grid-independency; the field variables used for comparison were the average bottom surface temperature and the two-phase pressure drop for the quasi-steady state solution. Three mesh models with different element densities were created: the first with 912,234 cells, the second with 2,446,686 cells, and a third mesh with 3,895,640 cells. While the first mesh presented stability issues due to its coarse resolution, the results for the second and third meshes were in a reasonable agreement with an error between solutions of 1.5% for absolute temperature, and 7.5% for the mean two-phase pressure drop. The second meshing scheme was therefore selected in order to balance the solution accuracy with the computational time, which is directly increased by reducing the element size due to the Courant number definition in such explicit formulations. The computational time for the solution to reach the quasi-steady state with the second mesh was ∼80 h while the third mesh took ∼200 h for the same problem using the HPC resources aforementioned.

It is important to note that currently there is a lack of a comprehensive experimental database or correlations for the flow boiling of R245fa in cooling layers for 3D IC applications, specifically in the described configuration of micro pin fins. Hence, at this point, an objective validation for the presented work may not be carried out, and therefore, it is disclosed here. It is also pointed out that the main objective from the present investigation is the comparison of the two-phase flow regimes with the VOF and CLSVOF methods, by also making some qualitative comparison of the predicted physics for bubble nucleation with behaviors reported in the literature. As a reference, a similar form of the presented flow boiling model has been used in our previous work [14] for the analysis of flow boiling of water in a silicon microchannel, and the results were compared with an experimental correlation obtaining a reasonable agreement for a range of heat fluxes in a flow boiling curve.

It is also relevant mentioning that in order to start the transient solution for the flow boiling case, the problem requires to be initialized with the single-phase flow field at the specified mass flux conditions (1000 kg/m² s for this investigation); the initial flow field also requires being isothermal at the specified saturation condition. Failure to do this may result in divergence of the solution.

For the present investigation, the behavior of two-phase flow regimes and their tracking through the specified methods are compared by fixing the operating conditions. It was discussed in the phase-change model section, the presence of the tunable parameter $λl$⁠, referred to as the relaxation time, that is directly related to the amount of generated vapor phase within the computational domain. In our preliminary work [15], the effect of varying this parameter on the evolution of the two-phase flow regime, temperature, and pressure fields was assessed. Since this parameter has not been calibrated yet for the specific application of interest due to experimental work in via of development, it is fixed here to a value of $λl=50$⁠, which corresponds to a midvalue of the range that was studied in Ref. [15]. Hence, the present study is of qualitative nature but useful in terms of providing insight about the capability predictions of the flow boiling physics and the interfacial resolution and implications of using the specified techniques.

Once the saturated, single-phase flow field is loaded, the heat flux condition at the bottom of the microgap is enabled and vapor formation will start. According to the conditions defined in the phase-change model and the implemented UDF, bubble nucleation will occur at fluid cells in the fluid–solid interface that exceed the saturation condition. Although the real process involves bubble nucleation at wall crevices with entrapped vapor, these effects cannot be captured through a CFD based-model considering the large number of features in the computational domain. However, it has been observed that such approach is reasonable and comparable to the experimental bubble nucleation processes reported in the literature for similar configurations, which are discussed in Sec. 3.2.

The evolution of the two-phase flow regime for the saturated flow boiling conditions in the silicon microgap is depicted in Fig. 3 from a top view (+z) perspective, where the liquid-vapor interface is represented by a blue isosurface. At first glance, it may be supposed that both methods provide similar results in terms of the flow regime evolution inside the microgap. For the first millisecond (ms) of flow, a few bubbles start to nucleate at the rear part of the pin fins, and these can be barely noticed in the corresponding subset of Fig. 3. As the flow simulation time increases, at t = 3 ms, it can be noted a significantly higher amount of vapor inside the microgap as the silicon temperature is increasing. For latter times, the two-phase flow regime keeps evolving as a result of the tracking capabilities of both the VOF and CLSVOF that allow vapor merging. The trend of vapor accumulation within the microgap continues until quasi-steady state conditions are reached, which are defined as the point in which the two-phase flow regime and variable fields such as temperature and pressure are not varying significantly with time. For the present cases, such point was observed approximately after 10 ms of two-phase flow conditions. It can be observed in Fig. 3 how the vapor structures between t = 7 ms and t = 9 ms are very similar, resulting from the balance between the applied heat flux and mass fluxes, boiling conditions through the phase-change model, and the continuous vapor generation and removal from the computational domain.

Although the overall trend in vapor generation and evolution between both methods is similar, there are important differences regarding the interfacial reconstruction when the details of each solution are analyzed. Figure 4 shows the comparison between both methods for the quasi-steady state solution at t = 10 ms, where the zoom-in windows provide further details of the vapor structures and their behavior near the inlet, center, and outlet of the analyzed silicon microgap. Starting from the inlet, and recalling that flow is in the positive x-direction, it may be observed that bubble nucleation starts from the first column of pins and that the reconstructed interfaces actually differ between both methods; for the standard VOF method, these bubbles have an oval shape and some of them are somewhat distorted, while the CLSVOF clearly offers a sharper interfacial reconstruction with rounded and smooth bubbles that better resemble the physics of the problem (see Sec. 3.2 for further discussion). When comparing the midsection details between both solutions (where there is a considerable amount of vapor due to buildup), it is observed that the VOF solution is formed by a set of vapor blobs varying in shape and size, and that the interface resolution is not optimal. In contrast, the CLSVOF solution is capable of providing a significantly better interpretation for the same problem through smoother vapor films that continuously break up at the front pins and remerge as they flow downstream. Such features are possible through the use of the LS function, which is continuous across the interface and helps for improved interfacial reconstruction and therefore, a more accurate surface tension calculation than the standard VOF method. Regarding the outlet section, it may be noted that with the VOF solution the section is nearly occupied in its entirety by vapor, while the CLSVOF technique handles the breaking and merging of vapor structures in a more physical manner.

In the present section, the bubble nucleation and detachment is discussed for a region of pin fins in the three first columns of the microgap, since vapor buildup is not significant here and allows a better interpretation of this mechanism and its representation by the implemented UDF. It is also important to remark that the present type of modeling has a number of advantages for the detailed analysis of two-phase flow, such as the three-dimensional representation of the liquid-vapor interface with access to any point within the computational domain. The CFD results for two-phase flow were postprocessed in this study with the software ANSYS^® CFD-Post^® 15.0, which allowed the representation and tracking of the interface in a variety of useful ways. Figure 5 shows a set of detailed views for a selected array of pin fins near the inlet section for analyzing the different predictions for bubble nucleation between both studied methods. In order to facilitate the interpretation of results, the isometric view and top view of the liquid-vapor interface for an array of eight upstream pin fins is shown; a top view of the cross section plane that depicts the volume fractions contours is also included, where the red color denotes the vapor phase, the blue color corresponds to the liquid phase, and the in-between range of colors indicates the interface. It is noted that there are important differences between the VOF and CLSVOF solutions; the isometric view for the VOF case shows that the bubbles are nucleated and detached from thin vapor films that originate approximately 90 deg from the stagnation point of the pin fins. There is also some bubble breakup into slightly distorted shapes, which may be attributed to the discontinuous VOF function across the interface and its impact on the surface tension calculation. In contrast, the CLSVOF solution provides a different representation in which bubbles form at the rear part of the pin fins and eventually detach as they grow and are removed by the surrounding flow; the interface curvature is considerably better resolved recalling that both solutions were obtained with the same mesh.

Overall, the predicted behavior of bubble growth and detachment is found to be in agreement with the experimental observations reported by Isaacs [26], who investigated the two-phase flow features of R245fa in a silicon microgap with a similar arrangement. The points with the lowest pressure in a cylindrical pin fin array correspond to those where boundary layer separation occurs, which for laminar flows is approximately 90 deg from the stagnation point; these points also have the lowest saturation temperature condition and hence promote bubble nucleation at wall crevices in the experiment. A similar process is approached by the implemented CFD model, where these flow field conditions are also predicted and hence, promote the bubble nucleation at computational cells in the solid–fluid interface meeting the phase-change criteria. Isaacs [26] reported that vapor bubbles originated at the mentioned location, and then, these were moved due to inertial forces to the recirculation zone (rear part of the pin fin, at 180 deg from the stagnation point), growing to sizes up to 50–75 μm before detachment. The last column of Fig. 5 depicts the bubble nucleation mechanism observed by Ref. [26], which is comparable to the predictions by the CLSVOF model. In support to this discussion, Fig. 6 shows a comparison of the transient bubble growth process at one of the pin fins adjacent to the flow inlet (upstream). It can be observed that under the same mesh and operating conditions, the CLSVOF model predicts smoother bubbles than the standard VOF model, where the interface resolution and surface tension calculation lead to early bubble departure and distortion for the latter.

It was mentioned in Sec. 3.1 that for the implemented phase-change model, bubble nucleation occurs in the fluid–solid interfaces that exceed the saturation condition. However, for the case of flow boiling in a silicon microgap, vapor generation actually comes from pre-existing vapor nuclei at the crevices of the solid surface (heterogeneous nucleation), which requires some degree of liquid superheat for the onset of bubble nucleation. A plausible approximation for this may be obtained with the employed phase-change criteria. The required liquid superheat for the onset of nucleation in this CFD approach may be controlled by tuning the so-called relaxation parameter $λl$⁠, which provides control on the vaporization rate (i.e., a low value of $λl$ would require large values of liquid superheat for bubble nucleation, and vice versa). As discussed in Sec. 3.1, this value was fixed in order to compare the two-phase flow features between the VOF and CLSVOF methods due to the scope and reach of the present investigation, but further details about its variation may be consulted in Ref. [15]. Figure 7 depicts a zoom-in comparison of the temperature field (cross section contour) around an array of pin fins, in an upstream location of the microgap, at a flow time of 6 ms for both analyzed methods. It can be observed that the liquid is superheated (T > T_sat) in the fluid periphery and wake zone of the pin fins, with the bubble nucleation process that has been discussed.

Noting that thermal control is one of the ultimate objectives sought through this type of modeling, the present section is devoted to discuss the implications on the temperature field solution by using the two described methods in this investigation. Figure 8 depicts the predicted temperature contours by the VOF and CLSVOF techniques for the quasi-steady state solution at t = 10 ms, where the temperature contours (displayed at the bottom surface of the microgap) look similar in terms of the qualitative behavior. However, recalling that the implemented phase-change model includes a source term for the energy equation that is directly proportional to the vaporization rate, it is expected a difference between both approaches due to the described two-phase flow features in Secs. 3.1 and 3.2.

In order to provide a quantitative comparison of the predicted temperature field by both methods, Fig. 9 shows a plot of the average temperature variation for each case across the heated area of the computational domain; the average temperature is plotted by equally dividing the heated area into 20 zones in the x-direction. It is observed that both solutions are in agreement until the first millimeter of heated length, and after that the CLSVOF method predicts higher temperature values than its counterpart; this discrepancy is attributed to the interfacial reconstruction and tracking capabilities of each method. The fact that temperature curves depart in the downstream regions is consistent with the higher amounts of vapor-phase predicted by the VOF method as discussed in Sec. 3.1 for Figs. 3 and 4, which therefore provide a larger interfacial area and energy sink through the defined source terms. This behavior is related to the intrinsic capabilities for curvature resolution of these techniques, which dictate the surface tension calculation and how the vapor structures form, break up, and interact in the computational domain. The largest temperature difference between both methods occurs near the outlet section, where the vapor amount difference is the highest; the steep temperature decrease is due to the heat spreading to the nonheated zones after the cylindrical pin fin section.

The variation of the two-phase pressure drop as a function of time was also assessed for both techniques and it is plotted in Fig. 10. The two-phase pressure drop is calculated as the difference in the mass-flow-weighted average pressure values between the inlet and outlet sections. As expected, the pressure drop increases as vapor structures are formed within the microgap, and follows this trend in a nonlinear fashion as the two-phase flow regime evolves. The plotted data is consistent with the discussed results of Fig. 3, where the vapor continues to build up until the point that it is balanced by the surrounding flow conditions leading to vapor evacuation. After t = 5 ms, the pressure drop begins to stabilize for both methods, and hence, it is assumed that the two-phase flow regimes are evolved. Noting that the pressure drop values oscillate due to the continuous generation/evacuation of vapor structures, a mean value is achieved at the so-called quasi-steady state. The mean pressure drop for the evolved two-phase flow regimes is 161.88 kPa and 161.58 kPa for the VOF and CLSVOF methods, respectively. Although both values are very similar, a difference in behavior should be noted between both techniques, where the CLSVOF solution presents larger pressure drop oscillations than the VOF results. This is again associated to how both methods reconstruct the vapor structures and the corresponding calculation of the surface tension forces.

The VOF and CLSVOF methods have been coupled to a mechanistic phase-change model for the analysis of the flow boiling behavior in a silicon microgap for 3D stacking of ICs. The main features and capabilities from each technique have been discussed, and the results for a case with fixed operating conditions have led to conclude that the CLSVOF technique offers a significantly better interface reconstruction that better approaches the physical behavior of the problem. The main advantage of using the LS function, continuous across the interface, is the capability of predicting smoother interfaces and therefore a more accurate surface tension calculation, which is of critical importance in microscale two-phase flows. The bubble nucleation process for both techniques was discussed and compared to experimental observations from the literature for a similar case, observing the early breakup of bubbles and its distortion for the standard VOF model due to surface tension calculation issues. In contrast, the CLSVOF offered a better interpretation of this process and a reasonable agreement in the qualitative behavior of the experimental reports.

Although the validation with experimental data for the specific configurations is pending for this work, it provides evidence of the capabilities to model flow boiling in a relatively complex domain with a number of pin fins for relevant microscale applications, such as the 3D-stacking of ICs. This type of modeling is lacking in the literature, and its further development and refinement may lead to a better understanding of this phenomena and its optimization for enhanced thermal control in emerging technologies.

The authors acknowledge financial support for this work from the DARPA ICECool Fundamentals program. Daniel Lorenzini was partially supported by CONACYT Fellowship No. 382817 for doctoral studies at Georgia Tech. The authors also acknowledge the valuable comments and suggestions from the reviewers that have improved the overall quality of the paper.

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⁻¹)

$hlv$ =

latent enthalpy of vaporization (J kg⁻¹)

$k$ =

thermal conductivity (W m⁻¹ K⁻¹)

$m˙lv$ =

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

$P$ =

pressure (Pa)

q″ =

wall heat flux (W m⁻²)

$Se$ =

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

$t$ =

time (s)

$T$ =

temperature (K)

$v$ =

velocity vector (m s⁻¹)

Greek Symbols

$α$ =

volume fraction

$λ$ =

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

μ =

dynamic viscosity (kg m⁻¹ s⁻¹)

ρ =

density (kg m⁻³)

$σ$ =

surface tension coefficient (N·m⁻¹)

Subscripts/Superscripts

int =

interface

$l$ =

liquid phase

si =

silicon

sat =

saturation

$v$ =

vapor phase