Conjugate Heat Transfer and Fluid Flow Modeling for Liquid Microjet Impingement Cooling with Alternating Feeding and Draining Channels

Abstract

Liquid microjet impingement cooling has shown the potential to be the solution for heat removal from electronic devices such as very-large-scale integration (VLSI) chips. The post-impingement dynamics of the jet, specifically the interaction between the liquid fronts on the surface engendered by the jets is a critical criterion improving the heat transfer characteristics. While some seminally important experimental studies have investigated this attribute, the amount of accurate data and analysis is limited by the shortcomings of real-life experiments. In this article, numerical investigations into the fluid dynamics and heat transfer in microjet cooling systems are carried out. Specifically, this paper addresses the question regarding the necessary fidelity of the simulations. Different Reynolds-averaged Navier–Stokes (RANS) models are compared to the Large Eddy Simulations (LES) simulation and the potential fidelity of different eddy-viscosity-based closures is clearly shown. Recommendations are made regarding the RANS closures that should give the best performance. It is demonstrated that the transition Shear Stress Transport (SST) model and $k$-$\omega$ SST model both show excellent ability to predict the local or average Nu, and also local level pressure coefficient $f$ with less than 5% difference in the range of 30 < Re_d < 4000, compared with the reference LES model. For the experimental measurements in the range of 130 < Re_d < 1400, the LES model, transition SST model and $k$-$\omega$ SST model all show less than 25% prediction error. Moreover, it is shown that the validity of the unit cell assumption for the temperature and flow distribution depends on the flow rate.

1. Introduction

Direct microjet impingement cooling shows great potential in high power devices or high-performance systems [1,2]. The post-impingement dynamics of the jet, specifically the interaction between the liquid fronts on the surface engendered by the jets is a critical criterion improving the heat transfer characteristics. In general, most of the experimental studies are based on the two configurations with common return [3,4,5,6] and locally distributed outlets [2,7,8,9]. While some seminally important experimental studies have investigated this attribute, the amount of accurate data and analysis is limited by the shortcomings of real-life experiments. Therefore, numerical studies into the liquid microjets design and analysis are needed, which can provide invaluable guidance and important reference for future investigations.

While we focus on the Reynolds number, the underlying rationale for the increment in heat transfer arises due to the velocity fluctuations. This is best characterized in terms of the turbulence intensity parameter, usually considered to be proportional to ${\mathrm{Re}}_d$, n < 1 [10,11]. The jet Reynolds number ${\mathrm{Re}}_d$, is defined as ${\mathrm{Re}}_d=\frac{d_i·V_{in}}{\mathsf{\nu }}$, where $V_{in}$ is the mean inlet nozzle velocity, $d_i$ is the inlet nozzle diameter, and ν is the kinematic viscosity. For typical regions in impingement jet flow with a common outlet, the flow field exhibits laminar flow properties at ${\mathrm{Re}}_d$ ≤ 1000 based on the hydraulic diameter as representative length scale. At ${\mathrm{Re}}_d$ ≥ 3000 the flow has fully developed turbulent features. A transition region occurs with 1000 ≤ ${\mathrm{Re}}_d$ ≤ 3000 [12,13]. In literature, there are numerous articles concerning Computational fluid dynamics (CFD) numerical modeling of impinging jet cooling [14]. Narumanchi et al. [15] showed that the standard $k$-$\omega$ turbulence model can achieve less than 20% difference with experimental data for circular submerged jet configurations (Womac et al. [16]). For the submerged confined jet configuration (Garimella and Rice [17]), the difference can be as low as 10% over a wide range of Reynolds numbers. Isman et al. [18] showed that the overall performance of renormalization group (RNG) and standard k-ε models are better in comparison with other models by considering stagnation region and wall jet region. Esch and Menter [19] showed that the Menter’s Shear Stress Transport (SST) model predicted heat transfer rates within 5% of those predicted by Durbin’s $\mathrm{v}$²-f (V2F) model. In [20], John Maddox compared the transition SST and the V2F turbulence models, and finally selected the transition SST model for the 3 × 3 jet array with common outlet flow based on the computational cost considerations. Subrahmanyam et al. [21] used Large Eddy Simulations (LES) to study the unsteady flow and heat transfer characteristics of a single impinging jet at Reynolds number of 20,000 at four normalized nozzle-to-impinging plate distances (0.5 ≤ z/d ≤ 2). Sung [22] used the standard two-equation k-ε turbulence model to study the effects of the jet pattern on single-phase cooling performance of hybrid micro-channel/micro-circular-jet-impingement.

Moreover, a few review literatures related to CFD modeling of impinging flow cooling have been reported in the literature. Polat et al. [23] reviewed the available numerical techniques to predict laminar and turbulent impingement heat transfer on a flat surface. Zuckerman and Lior [12] reviewed the suitability of different Reynolds-averaged Navier–Stokes (RANS) models in predicting average Nusselt number distribution and location and magnitude of the secondary peak in Nusselt number. The comparisons showed that direct numerical simulation (DNS)/LES time-variant models can accurately predict both Nu distribution and the secondary peaks. Moreover, the V2F and SST models also showed better predictions of fluid properties in impinging jet flows while the standard k-ε and $k$-$\omega$ models result in poor predictions. Behnia et al. [24] performed a critical review of important parameters in LES, DNS, and RANS-based techniques for computation of impinging flows. They concluded that the V2F model agrees very well with the experiments while k-ε model highly overpredicts the rate of heat transfer and yield physically unrealistic behavior. Among all the models, LES model shows encouraging results and clarified the understanding for the unsteady flow and heat transfer characteristics of multiple impinging jets even though the high computing cost [25,26,27,28,29,30,31]. In [29] the objective and key findings of different LES studies dealing with impinging flows in recent times are reviewed. Cziesla et al. [30] demonstrated the ability of LES to predict local Nu under a slot jet within 10% of experimental measurements. Draksler et al. [31] carried out LES simulation to provide a detailed insight into unsteady flow mechanisms and the associated heat transfer process of multiple impinging jets. In summary, the computation cost of LES can be considerably reduced as compared to the DNS if sacrificing the accuracy with small-scale turbulence [31]. Therefore, the LES model is used as a benchmark for the comparison with other turbulence models in Section 3.

Due to the complex internal structure of the microjet coolers for the heat removal from electronic devices such as VLSI chips, the available experimental studies on the locally distributed outlet configuration are very limited. In general, most of the studies for microscale Si [8] or ceramic-based coolers [9] are focused on the laminar flow regime, while for module level coolers the focus is on the turbulent flow regime [32]. Brunschwiler et al. [8] fabricated and analyzed a silicon-based microjet cooler with nozzle diameters ranging from 31 μm to 103 μm, resulting in a Reynolds number between 5 and 900 for deionized water. The deviation between the experimental data and the laminar model amounts to approximately 3%. However, the cost of the silicon microjet cooler is very high due to the introduction of lithography and silicon etching. As a more cost-efficient alternative, polymer-based impingement jet coolers show a great potential since it is already proved that very small nozzle diameter values are not necessary due to the saturation of the thermal improvement for further scaling of the nozzle diameter [33]. In [2,33,34], low cost polymer impingement jet coolers were demonstrated by using mechanical micromachining and 3D printing technology with nozzle diameters in the range of 300 to 600 µm. The flow rate used in electronic cooling applications is typically in the range between 300 mL/min to 2 L/min [33], resulting in Reynolds numbers ranging between 300 and 2000 based on the jet nozzle diameter. Therefore, the expected flow regimes for this type of microjet coolers are expected to be in the laminar, transitional and low Re turbulent flow regime.

To the best of the authors’ knowledge, there are no available data and research studies on the turbulence model analysis with locally distributed outlets configuration. In this work, we will investigate the thermo-fluidic behavior for microscale jet impingement coolers, especially in the transitional region for the larger nozzle diameter values of 600 µm. The structure of this paper is organized as follows: we will first introduce the full cooler level and unit cell modeling approaches. Next, different turbulence models based on CFD RANS models will be compared for the case of jet cooling with locally distributed outlets. Further, LES unsteady modeling is used as a benchmark for the comparison. Moreover, the unit cell model is compared with the full cooler level model under different flow rate conditions. In the last section, experimental results are used to compare with the unit cell model and full cooler model with different numerical modeling methods.

2. Conjugate Heat Transfer and Fluid Flow Model

Figure 1 shows a schematic representation of a jet impingement cooling solution with locally distributed outlets, identifying the four typical functional layers: inlet plenum, outlet plenum, jet nozzle plate and the impinging cavity. This configuration with inlet nozzles array and locally distributed returns can eliminate the accumulating cross-flow effects. For the optimization of the complicated cooler geometry and detailed analysis of the thermo-fluidic behavior inside the cooler, a local–global level CFD model approach is used to cover the challenging range of length scales that is required for the analysis. For the local level modeling, a unit cell model is used since the cooling nozzle array contains repeated cooling unit cells, with a central inlet and four neighboring outlets for single unit cell. Therefore, a unit cell model can be used as a fast approach to optimize the inlet/outlet nozzles parameters (nozzle diameter, nozzle pitch, nozzle plate thickness). However, the unit cell model cannot provide system level information, such as the total pressure drop of the cooler. The full cooler level model can predict the detailed thermo-fluidic behavior inside the cooler and the flow distribution and temperature uniformity. However, the computational cost of the full cooler level model is higher than that of a unit cell model. Therefore, the unit cell model and the full cooler modeling approach are complimentary. Where the unit cell model can be used to provide optimal configurational and geometrical jet parameters during the design phase, the full model is used to check the overall flow uniformity and the related distributor and collector properties as well as the overall chip temperature distribution.

2.1. Miroject Cooling Test Case

In order to investigate the conjugate heat transfer and fluid flow aspects related to the jet impingement cooling, a multi-jet impingement cooler [33] with a 4 × 4 inlet nozzle array is chosen as the test case for this study. Figure 1 shows the design of this cooler and the internal geometry with inlet chamber, nozzle jets and outlet chamber. The liquid jets will impinge directly on the surface of the silicon die, resulting in a high convective cooling. A systematical parameter sensitivity study of the cooler geometry was conducted in our previous research [33]. The 14 mm × 14 mm × 8.7 mm cooler is mounted on the 8 mm × 8 mm thermal test chip [34]. The diameter of the inlet and outlet tube is 2 mm. The thickness of the inlet chamber is designed as 2.5 mm for uniform flow distribution. The cavity height is 0.7 mm. In the local level, the inlet and outlet nozzle diameters are designed as 0.6 mm, while the size of the nozzle array is designed to match the dimensions of the Si chip, in this case 8 × 8 mm². The chip thickness is 0.2 mm. In order to study the flow impact on the chip temperature distribution, the conjugate heat transfer model within the CFD approach should include the fluidic part of the cooler, as well as the chip heat conduction part. In the next part, the details of the conjugate CFD model will be introduced and discussed.

2.2. CFD Model Introduction

In order to study the system level behavior of the cooler, including the total pressure drop, flow uniformity, and full chip temperature distribution, a full cooler level CFD model is needed. Since the cooler is fabricated with a polymer with very low thermal conductivity, the heat conduction through the cooler solid wall can be neglected [2], and only the fluidic parts of the cooler are included in the model, along with the solid domain for the Si chip. As shown in Figure 2, the internal fluidic domain is extracted from the CAD model of the cooler, including the inlet chamber, outlet chamber, nozzle plate and the impingement cavity. Moreover, the nozzle array with 4 × 4 inlet nozzles and 5 × 5 outlet nozzles distribution is shown in the enlarged view of the nozzle plate in Figure 2. The nozzle distribution is shown as a scalable system of repeated unit cells. For every single unit cell, there is one inlet located in the center and four surrounding neighboring outlets.

For the second part of the modeling analysis, the unit cell model is introduced in order to learn the local heat transfer and fluid flow behavior as function of the design parameters. As shown in Figure 3, a 1/8 model with symmetry boundary conditions is used for the unit cell models. For the dimensions of the unit cell, the inlet/outlet nozzle diameter and cavity height are all consistent with the full cooler model. The length of the unit cell L corresponds to the nozzle pitch and is 2 mm for the 4 × 4 nozzle array in this paper. The solid chip thickness is kept as 0.2 mm. The nozzle thickness is 2 mm.

2.3. Grid Sensitivity Study

For the meshing of the CFD models, hybrid meshing is chosen. The fluid domain mesh is chosen as tetrahedron mesh cells. The first layer thickness for the boundary layer along the wall is 1 µm to make sure the grid is fine enough to get y⁺ in the viscous sublayer. For the solid domain mesh prism cells are used with a 20 µm mesh size. The grid convergence index (GCI) is used for the meshing sensitivity analysis.

The grid convergence order is defined as follows:

$$p=\frac{\ln \left(\frac{f_3-f_2}{f_2-f_1}\right)}{\ln r}$$

(1)

where $p$ is the order of computational method. These solutions ($f_3$; $f_2$; $f_1$) are computed over three different grid levels (${\overline{h}}_3$; ${\overline{h}}_2; {\overline{h}}_1$), which are subsequently refined according to a constant grid refinement ratio r, shown as ${\overline{h}}_1=\frac{{\overline{h}}_2}{r}$ = $\frac{{\overline{h}}_3}{r^2}$.

Once the order of convergence p is known, the GCI can be calculated by using two subsequent results. In particular, if $f_3$ and $f_2$ are used and the final reported result is $f_3$, the one on the coarsest grid is defined as below:

$$\mathrm{GCI}=\frac{F_s{r}^p}{r^p-1}\left|\frac{f_3-f_2}{f_2}\right|$$

(2)

where the $F_s$ is a safety factor. It is also important to be sure that the selected grid levels are in the asymptotic range of convergence for the computed solution. The check for asymptotic range is evaluated using the equation as below:

$$\frac{{\mathrm{GCI}}_{23}}{r^p{\mathrm{GCI}}_{12}}\approx 1$$

(3)

where GCI₂₃ and GCI₁₂ are the values of GCI computed by considering, respectively, $f_2$; $f_3$ and $f_1$; $f_2$.

The GCI₁₂ and asymptotic range of convergence are listed in

Table 1. Grid convergence index (GCI) meshing sensitivity analysis of full model.

Temperature	GCI12	Asymptotic Range of Convergence
Stagnation Temp	0.002	0.99
Averaged chip Temp	0.004	1.01

for both the unit cell and full cooler level model. Based on the GCI analysis, the final meshing details are shown in

Table 2. Meshing comparison for the full model and unit cell model.

Model	Full Model	Unit Cell Model Reynolds-Averaged Navier–Stokes (RANS)	Unit Cell Model Large Eddy Simulations (LES)
Elements	8.5 M	0.4 M	3 M
Minimal Grid size	80 µm	20 µm	1 µm
Computation time	24 h	2 h	12 h

. The details of the mesh for the full model and unit cell model are both shown in Figure 4.

2.4. Modeling Methodology

For the model assumptions, the density, viscosity and other material properties of the fluid/solid are assumed to be constant during the simulation. All cavities are assumed to be completely filled with the liquid coolant, without any the presence of air (submerged jets). As for the flow boundary conditions, a Dirichlet boundary condition is used which means the velocity of the liquid at all fluid–solid boundaries is equal to zero (no slip condition). The boundary condition for the cooler inlet is set as a constant uniform inlet velocity while the static pressure for the outlet is set to 0 Pa, as a reference pressure. This means all pressure data obtained are specified relative to the outlet pressure. As for the thermal boundary conditions, the coolant inlet temperature is 10 °C. Moreover, a constant heat flux of 37.5 W/cm² is applied on the chip bottom to represent the power generation in the heating elements of the test chip. The fluid and solid interface is set as a flow-thermal coupled boundary condition. Since the cooler material is plastic with low thermal conductivity, the boundary walls of the fluid domain are set as adiabatic walls.

The modeling studies are performed in ANSYS Fluent. For both the full model and unit cell model, the Semi Implicit Method for Pressure Linked Equations (SIMPLE) [12,34,35] algorithm is used as the solution method. For the numerical scheme, the Quadratic Upstream Interpolation for Convective Kinematics (QUICK) scheme is chosen for the discretization. The convergence criteria were set at 10⁻⁵ for continuity, 10⁻⁶ for energy and 10⁻⁶ for k, ω and momentum (x, y and z components of velocity), respectively. For all the simulations, the net imbalance of overall mass, momentum and energy is kept below 0.02% [33,34]. In order to characterize the thermal and hydraulic performance of the cooler, the overall thermal performance and pressure drop are defined as dimensionless Nusselt and friction numbers, Nu and k, respectively. They are primarily a function of the Reynolds number Re:

$${\mathrm{Re}}_d=\frac{\rho ·d_i·V_{\mathrm{in}}}{\mu }$$

(4)

where the ${\mathrm{Re}}_d$ is defined as the Reynolds number based on the jet inlet nozzle diameter, and the $V_{\mathrm{in}}$ is defined as the averaged inlet nozzle velocity, $d_i$ is the inlet nozzle diameter, $\rho$ is the density of the fluid, and $\mu$ is the dynamic viscosity of the fluid. Moreover, the characteristic length used in the unit cell model and full cooler level model is based on the jet inlet nozzle diameter $d_i$, rather than the inlet diameter of the collector.

$$N{u}_{\mathrm{avg}}=\frac{h_{\mathrm{avg}}·d_i}{k_f}$$

(5)

$$N{u}_0=\frac{h_{\mathrm{stag}}·d_i}{k_f}$$

(6)

where is $N{u}_{\mathrm{avg}}$ is the averaged Nusselt number, $N{u}_0$ is the local Nusselt number, while $k_f$ is the thermal conductivity of the fluid.

$$h_{\mathrm{avg}}=\frac{Q}{A·\left(T_{\mathrm{avg}}-T_{\mathrm{in}}\right)}$$

(7)

$$h_{\mathrm{stag}}=\frac{Q}{A·\left(T_{\mathrm{stag}}-T_{\mathrm{in}}\right)}$$

(8)

where the $h_{\mathrm{avg}}$ represents the averaged heat transfer coefficient of the heat source surface, while $h_{\mathrm{stag}}$ represents the local heat transfer coefficient based on the stagnation chip temperature. The total heating area is $A$ and $T_{\mathrm{in}}$ is the inlet liquid temperature. The total chip power is defined as $Q$.

The flow inside the channels is assumed hydraulically developed flow. The friction of the flow is given by:

$$f=\frac{\Delta P}{\left(\frac{1}{2}\rho ·V_{\mathrm{in}}{}^2\right)\left(\frac{t}{d_i}\right)}$$

(9)

where f is the friction loss coefficients, t the thickness of the nozzle plate, and $d_i$ the hydraulic diameter. $\Delta P$ is defined as the pressure drop between the inlet and outlet nozzle in unit cell level.

3. Numerical Modeling Analysis

3.1. Unit Cell Modeling Analysis

In the numerical modeling analysis, we first check this physical problem with an unsteady solver, which is unsteady RANS (URANS) simulation to analyze the unsteady state behavior. The time step, 10⁻⁷ s in this case, is calculated based on the cell size and inlet velocity. Two velocity points with the stagnation point and recirculation point are monitored during the URANS modeling. As shown in Figure 5, after 600 time steps, the flow is fully developed, and from then on, there is no velocity fluctuation observed, which reveals the steady phenomenon. Therefore, this flow problem is steady at Re = 2048. So, for all the following simulations, we choose the RANS solver instead of URANS solver.

Three-dimensional RANS models with different numerical modeling techniques are simulated and compared, including $k$-$\omega$ SST model, laminar model, $k$-ε model, Transition SST and Spalart–Allmaras (SA) One-Equation Model. The Reynolds number ${\mathrm{Re}}_d$ is ranging from 30 to 4000, which covers the practical range for this electronic cooling application. The simulation results for ${\mathrm{Re}}_d$ = 1024 are shown in Figure 6. It can be seen that the temperature distributions for different turbulence models are different even though the flow streamlines show similar behavior. The temperature distributions of the SA model and the $k$-ε model show lower temperatures than the other models. In addition, the temperature patterns for $k$-$\omega$ SST, transition SST and the laminar model show similar temperature distributions at ${\mathrm{Re}}_d$ = 1024. The strong temperature differences using one turbulence model or another will be explained in details later based on the $N{u}_{\mathrm{avg}}$-${\mathrm{Re}}_d$ and $f$-${\mathrm{Re}}_d$ correlation curve.

As mentioned in the introduction session, LES is regarded as one of the most promising approaches for the simulation of turbulent flows with given length scale $\delta$, which is often connected to the mesh size. In [12], a well-resolved three-dimensional LES impinging jet model shows a very good heat transfer coefficient prediction. Therefore, we extend the impingement CFD RANS simulations to an unsteady LES simulation as a benchmark. The time step is set to 10⁻³ s. As shown in Figure 7, the flow is fully developed after around 10 iterations. In order to compare with the RANS model, the mean values of the variable are calculated by time-averaging of instantaneous results from 0.1 s to 0.5 s. The velocity and temperature simulation results with the LES model are shown in Figure 7. By using the LES model, the smaller scale flow behavior can also be captured. The simulated $N{u}_{\mathrm{avg}}$ and $N{u}_0$ for different RANS models are compared with the LES model at ${\mathrm{Re}}_d$ = 1024, as listed in

Table 3. Conjugate flow and thermal modeling under ${\mathrm{Re}}_d$ = 1024.

Turbulent Model	Re	Nuavg	Nuavg Difference	Nu0	Nu0 Difference
LES model	1024	35.14	0	39.42	0
Laminar model	1024	34.92	0.6%	39.42	0.1%
$k$-$\omega$ SST	1024	34.84	0.9%	39.37	0.1%
Transition SST	1024	35.05	0.3%	39.53	0.3%
k-ε model	1024	64.94	84.8%	81.95	107.9%
SA model	1024	69.27	97.1%	98.28	149.3%

. It can be seen that, on the one hand, the laminar model, $k$-$\omega$ SST and Transition SST model can produce better results than any of the high-Re models, matching $N{u}_{\mathrm{avg}}$ and $N{u}_0$ within 1% compared to the LES model, however by reducing the calculation time by a factor of 6. On the other hand, $k$-ε model and SA model show large $N{u}_{\mathrm{avg}}$ prediction errors up to 80%, and the $N{u}_0$ prediction errors are above 100%.

Detailed comparisons between the $N{u}_{\mathrm{avg}}$–${\mathrm{Re}}_d$ correlation, $N{u}_0$–${\mathrm{Re}}_d$ and $f$–${\mathrm{Re}}_d$ correlation are shown in Figure 8, Figure 9 and Figure 10. It can be seen that there is no transition observed for the Nu–${\mathrm{Re}}_d$ correlation while there is a clear transition point visible for the pressure coefficient correlation $f$–${\mathrm{Re}}_d$. For the Nu–Re correlations shown in Figure 8 and Figure 9, it can be seen that the laminar model, $k$-$\omega$ SST model and transition SST model show maximum prediction errors of 5% compared with the LES model. For the SA model or $k$-ε model, the difference increases from 5% to 20% as the Reynolds number increases. The reason why the $k$-$\omega$ SST and $k$-ε models are so different is because the $k$-ε model is based on the wall function model, which is good for high Reynolds number case. For the $k$-$\omega$ SST and transition SST models, both are based on the low Reynolds near wall model, the calculation starts from the near wall cells. As for the $f$–${\mathrm{Re}}_d$ correlations in Figure 10, the laminar flow model shows a large difference compared with the LES model around the transition point ${\mathrm{Re}}_d$ = 650. On the other hand, the $k$-$\omega$ SST model and the transition SST model match very well with the LES model after the transition point. The reason is that the transition SST model in ANSYS Fluent extends the traditional SST k-ω transport equations by tracking two additional variables for intermittency and transition onset using empirical correlations developed by Menter et al. [36]. Various authors have shown that the k-ω SST model shows unsatisfactory performance for jets, both free jets [37] and impinging jets [38]. This arises due to the eddy-viscosity hypothesis used in two-equation turbulence models, that over-predict the mixing rate in the CFD simulation [39]. However, for integral quantities of interest like the heat transfer, the interaction between the liquid fronts on the surface engendered by the jets is a critical criterion. This integral quantity of interest is still well predicted by the k−ω SST model.

In summary, as for the unit cell model, the transition SST model and $k$-$\omega$ SST model both can predict the average chip temperature, the stagnation temperature on the chip, and also the pressure drop with less than 5% difference, compared with the reference LES model, when the ${\mathrm{Re}}_d$ is in the range between 30 to 4000.

3.2. Full Cooler Level Modeling Analysis

The flow and temperature distributions are shown in Figure 11 shows the results of the full cooler level CFD model for the flow distribution in the fluid domain and the temperature distribution in the solid domain of the Si chip, for different flow rates ranging from 100 mL/min to 1000 mL/min. The velocity in Figure 11 is the inlet tube velocity, while the ${\mathrm{Re}}_d$ here is based on the inlet nozzle diameter and velocity. The corresponding ${\mathrm{Re}}_d$ number is between 130 and 1400. The chip temperature distribution map is linked to the velocity field inside the cooler. It can be seen that the temperature of the chip edge is higher for low ${\mathrm{Re}}_d$ numbers, due to the flow nonuniformity. Moreover, the slight asymmetry of the temperature distribution, visible at higher flow rates, is due to the asymmetric placement of the outlet connection, resulting in an asymmetric flow behavior.

3.3. Unit Cell versus Full Level Model

In order to evaluate the validity of the unit cell model, the temperature modeling data are compared with the full cooler level model. In this comparison, the transition SST model is used for both the unit cell model and full cooler model since this model showed a good agreement with the reference LES model in the previous section and offers a good compromise between the model accuracy and computation cost, especially for the large simulation domain of the full cooler model.

Figure 12 shows the temperature distribution at the location of the heat sources in the Si chip, calculated by the full cooler level model for flow rate values from 100 mL/min to 1000 mL/min. The maximum, minimum and average chip temperature are extracted as a function of the flow rate. In general, the accuracy of the unit cell model depends on the symmetry of the flow and temperature patterns. The comparison between the unit cell model and full cooler model results in Figure 12 shows a higher flow non-uniformity at low flow rate values with a higher local relative flow rate in the central nozzles, resulting in higher temperatures at the chip corners. For a moderate flow rate of 650 mL/min, the unit cell model shows a good agreement with the full model. The full profile comparison between the unit cell and full model are shown in Figure 13, which provides more information on the usability of the unit cell model. It also shows where the unit cell assumption is valid and that it can change with the flow rate.

Moreover, the pressure drop of the unit cell model and full cooler level model are compared in Figure 14. It can be seen that the pressure drop of the unit cell from the nozzle plate contributes roughly 1/3 of the total pressure drop of the full cooler at the flow rate of 1000 mL/min.

4. Experimental Validation

4.1. Test Case Demonstration

In order to validate the numerical simulations of the turbulence models, a microjet cooler has been fabricated by using 3D printing technology, shown in Figure 15. The nozzle array is 4 × 4 with printed nozzle diameters of 0.58 mm, compared to the nominal design value of 0.6 mm. In our previous work [34], the impact analysis showed that the 5% diameter variation for the fabrication only results in a reduction of the average temperature and the minimum temperature of 4.7% and 4.3% for the same flow rate respectively.

4.2. Experimental Set-Up

Figure 16 shows a schematic of the experimental setup for the measurement of the hydraulic and thermal performance of the microjet cooler. The closed flow loop includes flow meter, pressure transducer, heat exchanger, chiller, pump, filter and temperature sensor. The closed flow loop can be operated in a controlled mass flow rate mode or a controlled pressure mode. The working range of the magnetically coupled gear pump is within a maximum flow rate of 3 kg/min and a maximum pressure of 11.5 bar. The accuracy of the mini Cori-FLOW mass flow meter is ± 0.2% for a range of 0.1 to 3 kg/min. A differential pressure gauge (EL-PRESS) is used to measure the pressure drop across the cooler. The accuracy of the pressure transducer is ± 0.5% of the full scale in the range between 0.2 and 5 bar. The particle filter with a mesh size of 25 μm is used. In addition, thermocouples with an accuracy of 2.2 °C are used to measure the coolant temperature before and after the cooler. A liquid–liquid heat exchanger is used to cool the coolant back to the set-point. In this work, DI-water is used as the coolant during the tests, with specified temperature at 10 °C and ambient temperature is kept at 25 ± 1 °C.

The chip temperature is measured in the thermal test chip [34] with a 32 × 32 array of temperature sensors and configurations of 832 programmable heater cells. The size of the single square cell is 240 × 240 μm². All these cells contain a diode in the center of the cell as a temperature sensor, resulting in a detailed temperature map measurement with 32 × 32 ‘thermal pixels’ across the die surface. The voltage drop across the diode for a constant current is used as the temperature sensitive parameter of the sensor. The 95% confidence interval of the calibrated sensitivity is −1.55 ± 0.02 mV/°C for a current of 5 μA in the temperature range between 10 and 75 °C. The chip temperature sensors allow absolute temperature measurements with an accuracy of ± 2 °C.

4.3. Results and Discussions

In Figure 17, the measured average chip temperature values are compared with the full CFD model and unit cell model results. The measured flow rate is ranging from 100 mL/min to 1000 mL/min, resulting in a ${\mathrm{Re}}_d$ number from 130 to 1400. The heat flux applied on the thermal test chip is 80 W/cm². Similar with the unit cell model analysis, different turbulence models are used for the full cooler level model, including laminar model, $k$-ε model, $k$-$\omega$ model, Transition SST model and SA model. It can be seen that the full CFD model with SA model overestimates the Nusselt number by a factor of 4 comparing with the experimental result. Moreover, the full model with $k$-ε model also shows very high prediction errors compared with the experiments. As expected, the LES model shows good agreement with the measurements. In general, the comparison shows that the laminar model, the $k$-$\omega$ model and the transition SST model show good agreement with the measured chip temperature, for the ${\mathrm{Re}}_d$ number from 130 to 1400.

5. Conclusions

This paper presents the conjugate flow and heat transfer modeling for microjet cooling with locally distributed outlets. We analyzed the turbulence models for their applicability in unit cell level models and full cooler level models for these liquid microjet impingement cooling applications. The results of the different turbulence models based on steady state CFD RANS models for a microjet unit cell are benchmarked with LES results. It is concluded that the transition SST model and $k$-$\omega$ SST model both can accurately predict the average chip temperature, the stagnation temperature on the chip, and the pressure drop with less than 5% difference, compared with the reference LES model. Moreover, the unit cell model is validated with the full cooler level model for different flow rate conditions. However, the usability of the unit cell model changes with the flow rate. A test case with a microjet cooler has been demonstrated by using 3D printing technology in order to validate the numerical simulations of the turbulence models. The experimental results are compared with the unit cell model and full cooler model with different numerical modeling methods.

In summary, the transition SST model and $k$-$\omega$ SST model both show excellent ability to predict the local or average Nu, as well as the local level pressure coefficient f with less than 5% difference in the range of 30 < Re_d < 4000, compared with the reference LES model. For the comparison with experimental measurements, the LES model, transition SST model and $k$-$\omega$ SST model all show less than 25% prediction error as the ${\mathrm{Re}}_d$ number ranging from 130 to 1400.