Critical Heat Flux Condition and Post-Critical Heat Flux Heat Transfer of Carbon Dioxide at High Reduced Pressures in a Microchannel

Flow boiling in microchannels has been a topic of much fundamental and practical interest over the better part of the 21st century. However, there are several major roadblocks hindering its wide use in cooling systems for high-power electronics. Chief among them are flow boiling instabilities, the critical heat flux (CHF) condition, high-pressure drop, temperature overshoot, high superheat temperature at the onset of nucleate boiling (ONB), and the prevalence of laminar flow, which is associated with a low-grade heat transfer coefficient. As a plethora of systems require low pumping power (also low-pressure drop) and low fluid inventory, it is desired to leverage the latent heat of vaporization to its fullest extent by enabling high mass quality flows (i.e., reduced flowrate, and thus, reduced pressure drop and low fluid inventory). However, the CHF conditions impose a stringent requirement on the ability of these systems to operate at desired high mass qualities. As a result, more recent research efforts have started to explore means to operate microscale flow boiling systems at mass qualities as high as 0.95 [1]. As part of this endeavor, several studies have explored the flow boiling of carbon dioxide (CO₂) in microchannels at a reduced pressure of 0.5 and higher, as it offers several distinct advantages. Unlike many potential coolants, the critical temperature of CO₂ of 31.1 °C is well aligned with the temperature domain permissible by electronic systems, and high reduced pressure has been shown to mitigate flow boiling instabilities [2–7]. Carbon dioxide has a low viscosity, allowing microscale cooling systems to operate with the enhancement offered by turbulent flow. It also has a low surface tension, which minimizes the ONB, eliminates temperature overshot, and helps mitigate the rapid bubble growth instability. In addition, at high reduced pressures, the thermal conductivity of the vapor can be comparable to that of liquid, albeit still lower, helping mitigate some of the harsh penalties imposed by the critical heat flux condition (i.e., the transition to post-CHF might not be so harmful since the vapor can still support heating loads at acceptable surface temperatures). (It should be noted that CO₂ has a global warming potential of 1, which is an order of magnitude smaller than hydrofluorocarbons, such as R134a with a global warming potential of 1430.)

Equipped with the realization of these potential benefits, Parahovnik et al. [8–10] and Parahovnik and Peles [11–13] embarked on a journey to reveal the nature of flow boiling heat transfer of CO₂ in microchannels at high reduced pressures. Initially, subcooled flow boiling was examined, and the ONB [8], subcooled heat transfer coefficient [8,14–17], saturated flow boiling patterns [12,18,19], bubble dynamics [11], the piston effect [9], and CHF conditions were examined [12]. Subsequently, saturated flow boiling heat transfer was studied and compared with previous models and correlations [12]. In parallel, Adeoye et al. [20,21] and Adeoye and Peles [22] studied flow boiling heat transfer with micro-impingement jets and elucidated the heat transfer coefficient and critical heat flux. Additionally, Asadzadeh et al. [23] investigated heat transfer characteristics in a silicon microchannel with a NACA pin fin configuration.

The current study extends these previous efforts and provides an in-depth analysis of the transition through the CHF condition. It delineates and quantifies the transition from pre-CHF condition to post-CHF condition and the associated heat transfer toll related to the transition. Specifically, it provides a measure of the criticality of the critical heat flux condition (i.e., the severity of the temperature increase, or the drop in the heat transfer coefficient, during the transition). The coefficients and exponents of the Katto CHF correlation [24] and the Bishop et al. heat transfer coefficient correlation [25] were examined and subsequently modified to better predict the experimental results.

The microfluidic device was made from two parts: a substrate with resistive temperature detectors (RTDs) and heaters, termed part 1, and a second custom-made substrate, part 2, fabricated by Mindrum Precision^® with a 296-μm hydraulic diameter microchannel and an embedded 81-μm hydraulic diameter orifice with a length of 1.8 mm (i.e., l/d ratio of about 22). The width and height of the orifice were 500 μm and 45 μm, respectively, and the width and height of the microchannel were 1800 μm and 160 μm (Fig. 1(a)). (Note that the orifice had slightly rounded edges due to tooling and computer-aided design (cad) software and a measurement of the area and the perimeter out of there resulted in a hydraulic diameter of 81 μm rather than 82.6 μm.)

Part 1 was made of a 0.5-mm thick fused silica plate, which included the electrical circuit (i.e., the RTDs, heaters, electrical vias, and electrical pads). These were made by the deposition of metallic layers (i.e., titanium, platinum, and aluminum) and multiple silicon oxide layers for electrical insulation. The RTDs were custom-made from vapor-deposition layers. The heater and RTDs were made using similar microfabrication processes but with much different geometrical shapes to render their functionality. The rectangular heater had dimensions of 5.82 mm by 0.89 mm (Fig. 1(b)) with a low electrical resistance of about 70 Ω, generating heat when the DC power was turned on. The serpentine RTDs had footprints of 50 μm by 60 μm, with a resistance of about 2.5k Ω each. The RTDs were placed on top of the heater at locations specified in Fig. 1(b). The linear dependency of the metal resistivity of the RTDs on the temperature was leveraged to measure the local, temporal temperature of the heated wall.

Figure 1(c) depicts a top view of part 1, where the microchannel, inlet, outlet, and orifice were marked with a dashed white line. The flow is shown inside the dashed area, and the restriction in the middle represents the orifice. The heater is shown as a darker area inside the dashed region, split into two independently controlled zones (upstream and downstream of the orifice). (Only the downstream heater was used in this study.) The heater and the RTDs were made from metallic layers of titanium (7 nm) and platinum (30 nm) deposited on a fused silica wafer and later sliced into rectangular dies. These components were placed inside the microchannel and had intimate contact with the flow. The electrical path included electrical vias, electric pads, spring pins, custom-made printed circuit boards (PCBs) mounted on the package, and standard nine-pin connectors to sense and control the in-channel elements (i.e., the heater and the RTDs). The vias had an additional aluminum layer (∼1 μm) to significantly reduce their electrical resistance and render it negligible compared to the RTDs and heaters. For insulation, the metallic layers were covered with a silicon oxide layer with a thickness of ∼1.1 μm that was later selectively etched to expose the electrical pads that interfaced with the push pins. The power was supplied from a DC power supply (Kysight^®, E3645A^®) through electrical wiring, push pins, contact pads, and the vias.

The microfluidic device was enclosed in a special housing that kept its structural integrity and enabled sealing between the two pieces while supplying the flow at desired experimental conditions (Fig. 1(d)). Additionally, the package enabled to read the temperatures for three RTDs located at length to hydraulic diameter ratios of 3.2, 7.4, and 11.6, and to apply the desired thermal heat load. The two-phase flow conditions were achieved using a throttling process through the embedded micro-orifice and an external preheater (HeatBrisk^®). The orifice allowed to study of a wide range of mass qualities and pressures in the microchannel, and its operational design and thermodynamic states were described in detail by Parahovnik et al. [10]. Because of the low mass flow rates through the microfluidic device, the carbon dioxide was allowed to flow through the open loop experimental setup and was released to the ambient. The experimental flow setup was intended to provide maximum control and flexibility to study a range of conditions, such as pressure, inlet temperature, and mass flux. The experimental setup had pressure transducers (Omega^®, Omega Engineering Inc., Norwalk, CT), a mass flowmeter (Alicat^®), and a high-precision T-type Thermocouple that enabled to measure of pressure, mass flowrate, and fluid temperature at the inlet of the microdevice, respectively. Additional details about the pressure and mass flow control are provided by Parahovnik et al. [8,11].

where ΔT_sat = T_RTD–T_sat and q″_loss = (T_RTD–T_amb)/R_t.

The mass of the microfluidic device and the microdeposited components were small such that the RTDs and the heater's thermal time response were practically instantaneous. Nevertheless, data was recorded only 30 s after a new heat load was set, allowing ample time for the temperature to stabilize. Subsequently, the data were sampled for an additional 20 s with a sampling frequency of 1 kHz and were averaged to obtain the measured point. The standard deviation was also calculated and added to the zero-shift accuracy to estimate the uncertainty of the reading.

The experimental data points, given in Table 1, cover pressures from 1.7 MPa to 6.9 MPa (i.e., reduced pressures of 0.23–0.935) and mass fluxes from 1405 kg/m²s to 2776 kg/m²s (based on the microchannel's cross sections of A = 0.284 × 10⁻⁶ m²), corresponding to a Reynolds number based on the homogeneous two-phase flow model ranging from 6622 to 32,248. Three length-to-diameter ratios were studied, l/d, of 3.2, 7.4, and 11.6. As the flow exited the orifice, it gradually expanded while passing through the first two RTDs (i.e., l/d = 3.2 and 7.4) before becoming fully developed at the third RTD at l/d = 11.6. (Since the flow was turbulent, it is believed that around l/d∼10 the flow was fully developed [27].)

Figure 2 depicts the boiling curves for reduced pressures of 0.67 and 0.91 at three different length-to-diameter ratios. As the flow was still developing, the region closest to the orifice exhibited the highest CHF value, a trend that is also consistent with the model of Katto [24] (although outside of the correlation boundaries), and the highest heat transfer coefficient. (Regardless of the developing flow versus developed flow, CHF typically starts downstream at the trailing edge of the heater and then propagates upward, but since the mass quality didn't change much along the heater, it is believed that the higher CHF near the inlet was primarily because of the effect of the developing flow.) At low reduced pressures, the transition across the CHF condition can be rapid and accompanied by a large increase in surface temperature (Fig. 2(a)). As a matter of proper engineering practice, the CHF condition should be typically avoided. However, it was found that as the reduced pressure increased, the transition to post-CHF conditions (Fig. 2(b)) was gradually accompanied by only a moderate increase in the surface temperature. At the same time, the absolute CHF value decreased with pressure. The vapor-to-liquid ratios of the thermal diffusivity, dynamic viscosity, and density approach unity at high reduced pressure smoothen the transition from pre-CHF to post-CHF conditions.

The temperature increase following the transition to post-CHF condition (i.e., ΔT_CHF) was inversely proportional to the reduced pressure (Fig. 3(a)) with a linearity Pearson coefficient [28] of –0.92, –0.9, and –0.74 for l/d of 3.2, 7.4, and 11.6, respectively. The coefficient indicates a strong linear dependency between the parameters. The CHF above a reduced pressure of 0.5 gradually declined (Fig. 3(b)). As the difference between the thermophysical properties of the liquid and the vapor diminished with pressure, ΔT_CHF approached 0 K near the critical condition [29]. Although the ΔT_CHF asymptotically approaches zero at the critical pressure, the CHF did not, and its lower limit was about 50 W/cm². This is somewhat surprising since the surface tension and the differences between the properties of liquid and vapor diminish with pressure, which should eliminate any difference between the pre-CHF to post-CHF conditions at the critical state. This lower limit is hypothesized to be related to the piston effect, a unique thermalization mode of near-critical fluids [9].

Since the velocity distribution varied significantly from RTD₁ to RTD₃, the heat transfer process was considerably different. As such, the heat transfer coefficient and the CHF varied between RTD₁ and RTD₃ (Fig. 4). Since the velocity near the wall at RTD₁ was much higher than downstream, the heat transfer coefficient was considerably higher.

As shown in Fig. 5, for reduced pressures above 0.6 (i.e., density ratio above 0.15), the correlation underpredicted the critical heat flux condition for all RTDs. The correlation overpredicted the experimental results for reduced pressures below 0.6 for RTD₃. The MAE values for RTD₁, RTD₂, and RTD₃ were 20.7%, 16.9%, and 12.4%, respectively.

Since the exponents had the strongest impact on the accuracy of the correlation, a sensitivity test of the exponents f and g was performed and is shown in Fig. 6. Each exponent was independently varied while the other was kept constant at its optimum value, and the average MAE was obtain based on results for all three RTDs. The MAEs were sensitive to the coefficients, and a shift of 10% from the optimal values doubled the MAEs. Figure 7 compares the modified correlation (Eq. (3)) with the experimental results for the three length-to-diameter ratios. The CHF peaked around P_r = 0.5 and gradually decreased with pressure. (For comparison, Boyd [32] reported a peak value at a reduced pressure of 0.75 for water.) However, as stated earlier, the lower CHF at high reduced pressures resulted in a moderate decline in the post-CHF heat transfer coefficient.

The new correlation also better captured the experimental trends and predicted the experimental data at reduced pressures above 0.6. The MAEs values for RTD₁, RTD₂, and RTD₃ were 5.86%, 4.84%, and 6.69%, respectively, compared to 20.7%, 16.9%, and 12.4% for the original Katto [24] correlation.

For pre-CHF conditions, the liquid was the dominant phase that directly carried the heat away from the wall through vaporization or convection [34]. For post-CHF conditions, convection to the vapor phase was an essential part of the heat transfer process, and the liquid phase served to a large extent as a thermal capacitor [1,34]. As the reduced pressure increased, the liquid-to-vapor thermodynamic property ratios, such as the thermal conductivity, thermal diffusivity, asymptotically approach unity, and the surface tension diminished [36,37]. Thus, the post-CHF heat transfer coefficient was expected to be appreciably higher than for low pressures, which could potentially support moderate heat fluxes at acceptable surface temperatures.

The experimental results were compared to five correlations detailed in Song et al. [35] to obtain quantitative measures of the post-CHF heat transfer coefficient under the high reduced pressures examined in this study, given in Table 2. Through the maximum root-mean-square (RMS) defined as √((h_exp.–h_calc.)/h_exp.)²). The Miropolskiy [38] and the Groeneveld [39] correlations had similar terms (i.e., vapor-to-liquid density ratio, mass quality, vapor-based thermal conductivity and Reynolds number, and wall-based Prandtl number) but different functional dependencies (i.e., coefficient and exponents). The governing terms in the Bishop [25] correlation were similar (i.e., density ratios, Prandtl and Reynolds numbers, and thermal conductivity), but the terms were mainly evaluated based on the wall and vapor conditions. The Quinn [40] correlation accounted for a viscosity ratio (i.e., μ_v/μ_l), which wasn't considered in the other correlations. The Song et al. [35] correlation was formulated through coefficients, which depend on the Reynolds number, Prandtl number, and mass quality that are used as multipliers of the Nusselt number calculated through the Dittus–Boelter correlation [27] (where the Reynolds number is based on the homogeneous two-phase flow model [34] and the Prandtl number is evaluated based on the wall temperature).

The correlations were compared to the heat transfer coefficients at RTD₃ (i.e., the location corresponding to a fully developed flow). To verify that the correlations were compared with stable post-CHF conditions, the heat transfer coefficients were calculated for a wall temperature of at least 40 K above the saturation temperature. The correlations were compared to 258 data points collected from 22 experiments. The calculated MAEs were 75.2%, 91.1%, 36.2%, 74.7%, and 95.2% for the Miropolskiy [38], Groeneveld [39], Bishop et al. [25], Quinn [40], and Song et al. [35], respectively. The Bishop et al. [25] correlation best predicted the experiments, which can be because it evaluates the thermophysical properties based on the wall temperature, suggesting that the superheated vapor near the wall was the dominant heat transfer medium. All the other correlations were mainly based on the saturated properties, which failed to predict the experimental results.

A similar procedure to the one described in Sec. 3.2 was performed to examine the validity of the exponents in the Bishop et al. correlation. The exponents, Table 4, were examined for their MAEs while each exponent was varied independently.

Exponents m, n, and s varied between 0.6–0.9, 0.7–1.5, and 0.3–0.8, respectively. The MAE was invariant to the exponent of the Pr_w (i.e., n), with an MAE variation lower than 1.25%. Variations of the m and s exponents had a significant influence on the MAEs. However, contrary to the CHF expression, the exponents, in this case, were consistent with the original Bishop et al. [25] correlation and, therefore, were not changed, Fig. 9.

Equation 5 predicted the experiments well for a reduced pressure of 0.49 and a reduced pressure of 0.78, Figs. 10(a) and 10(b), while it underpredicted the subcooled flow boiling experiment of reduced pressures of 0.85, 0.89, 0.91, and 0.94. Regardless, Eq. (5) captured well the experimental trend. For experiment #3, the thermophysical properties of the vapor were much inferior for post-CHF conditions resulting in a lower post-CHF heat transfer coefficient, while for experiment #22, with a reduced pressure of 0.94, the transition to post-CHF was relatively smooth and resulted in comparably high post-CHF heat transfer coefficient (i.e., 30 kW/m²K at ΔT_w = 40 K). A possible explanation for the underprediction of the heat transfer coefficient for experiment #22 can be related to the presence of the piston effect, which for a reduced pressure of 0.94, is predicted to account for 25% of the heat transfer rate [9]. By assuming an average calculated heat transfer coefficient of 20 kW/m²K for experiment #22. By adding the poison effect into consideration (i.e., 25% of the heat was dissipated into the bulk fluid through the piston effect), the actual heat transfer coefficient should be around 25 kW/m²K, which is approximately the value obtained from the experiments, Fig. 10(c). Since the Bishop et al. [25] correlation doesn't account for near-critical condition effects, it underpredicted the near-critical heat transfer coefficients.

The center velocity of a flow exiting an orifice is the highest and narrowest at the orifice exit. As the flow propagates downstream, the center velocity declines, but the nonzero velocity region becomes wider, as shown in the cone sketched in Fig. 1(c). The enhancement factor was derived by minimizing the MAEs for the RTDs at l/d of 3.2 and 7.4, Fig. 11. The enhancement factor for RTD₁ (i.e., t₁) and RTD₂ (i.e., t₂) was found to be 2.3 and 0.8 with MAEs of 38.5% and 27.5%, respectively.

The critical heat flux condition and post-CHF heat transfer coefficient of carbon dioxide at high reduced pressures in a microchannel were experimentally studied and compared to available correlations. Based on the experimental results, two new correlations were proposed for the post-CHF heat transfer coefficient and the CHF condition. It was found that as the reduced pressure increased, the CHF decreased, and the post-CHF heat transfer coefficient increased. As a result, the typical rapid increase in the surface temperature following the transition to film boiling was moderate, unlike lower reduced pressures where the burnout phenomenon is common. Moreover, a micro-orifice and its associated developing flow were demonstrated to enhance the heat transfer coefficient and the CHF condition significantly. Additionally, it was found that the post-CHF heat transfer coefficient prediction was best when the thermophysical properties were evaluated at the wall temperature. Consequently, a new correlation for CHF condition, based on Katto [24] formulation, and a new correlation for fully developed post-CHF heat transfer coefficient were proposed, resulting in reduced MAEs. An enhancement factor was derived to predict the increase in the heat transfer coefficient in the developing region.

This work was funded by the Office of Naval Research (ONR) under award number N00014-21-1-2653. The authors would like to acknowledge the Cornell Nano-Scale Science & Technology Facility (CNF) for staff support, help, and use of facilities to fabricate the microdevices. The CNF is a member of the National Nanotechnology Coordinated Infrastructure (NNCI), supported by the National Science Foundation.

• Office of Naval Research (Award No. N00014-21-1-2653l Funder ID: 10.13039/100000006).

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

a =

CHF correlation coefficient, Eq. (2)

A =

channel's cross section (0.284 × 10⁻⁶ m²)

b =

CHF correlation exponent, Eq. (2)

C =

correction coefficient for the developing flow, Eq. (4)

CHF =

critical heat flux (W/cm²)

D_h =

hydraulic diameter (m)

e =

CHF correlation coefficient, Eq. (2)

Exp # =

experiment #, Table 1 

f =

CHF correlation exponent, Eq. (2)

F₁, F₂, F₃ =

Song et al. [31] coefficients

g =

CHF correlation exponent, Eq. (2)

G =

mass flux (kg/m²s)

h =

heat transfer coefficient (kW/m²K)

H =

enthalpy (J)

h_fg =

latent heat of vaporization (J/kg)

k =

thermal conductivity (W/mK)

l/d =

length-to-diameter ratio

m =

Bishop [25] exponent, Table 4 

MAE =

mean averaged value, Eq. (1)

n =

Bishop [25] exponent, Table 4 

N =

number of samples, Eq. (1)

Nu =

Nusselt number

ONB =

onset of nucleate boiling

P =

pressure (MPa)

Pc =

critical pressure (i.e., =7.337 MPa [33])

Pr =

Prandtl number

q″ =

heat flux (W/cm²)

q″_loss =

heat losses to the environment from the package (W/cm²)

R_t =

thermal resistance (K·cm²/W)

R² =

coefficient of determination

Re =

Reynolds number

RMS =

root-mean-square, Table 3 

s =

Bishop [25] exponent, Table 4 

t =

enhancement factor, Eq. (6)

T =

temperature (K)

x =

mass quality

x_m =

Song et al. [31] mass quality defined in Table 3 

z =

dimensionless distance from the inlet (e.g., z = l/d)

@RTD =

property measured at the location of the RTD (K)

Greek Symbols

ΔT =

temperature difference (K)

μ =

viscosity (Pa·s)

ρ =

density (kg/m³)

σ =

surface tension (N·m)

Subscripts

amb =

ambient

b =

bulk fluid

B =

Bishop correlation [25], Fig. 8 

c =

calculated

exp =

experimental\measured

fd =

fully developed

i =

summation index, Eq. (1)

in =

at the inlet of the microchannel after the orifice

Katto =

based on the Katto correlation [24]

l =

liquid phase

max =

Maximum

r =

reduced (e.g., reduced pressure)

RTD =

property measured at RTD_'s location

sat =

saturation

TP =

two-phase

v =

vapor phase

w =

wall

1, 2, 3 =

numbers of the RTDs