A national laboratory of the U.S. Department of Energy

Office of Energy Efficiency & Renewable Energy

National Renewable Energy Laboratory 

Innovation for Our Energy Future 

Numerical Simulations of 

Boiling Jet Impingement 

Cooling in Power Electronics 

S. Narumanchi, V. Hassani, and D. Bharathan 

National Renewable Energy Laboratory 

A. Troshko 

Fluent, Inc.

Presented at ITHerm 2006 (Intersociety Conference on Thermal and 

Thermomechanical Phenomena in Electronic Systems) 

San Diego, California 

May 30–June 2, 2006 

Conference Paper 

NREL/CP-540-39401 

December 2006 

NREL is operated by Midwest Research Institute ● Battelle Contract No. DE-AC36-99-GO10337 
NOTICE

The submitted manuscript has been offered by an employee of the Midwest Research Institute (MRI), a 

contractor of the US Government under Contract No. DE-AC36-99GO10337. Accordingly, the US 

Government and MRI retain a nonexclusive royalty-free license to publish or reproduce the published form of 

this contribution, or allow others to do so, for US Government purposes. 

This report was prepared as an account of work sponsored by an agency of the United States government. 

Neither the United States government nor any agency thereof, nor any of their employees, makes any 

warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or 

usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not 

infringe privately owned rights. Reference herein to any specific commercial product, process, or service by 

trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, 

recommendation, or favoring by the United States government or any agency thereof. The views and 

opinions of authors expressed herein do not necessarily state or reflect those of the United States 

government or any agency thereof. 

Available electronically at http://www.osti.gov/bridge

Available for a processing fee to U.S. Department of Energy 

and its contractors, in paper, from: 

U.S. Department of Energy 

Office of Scientific and Technical Information 

P.O. Box 62 

Oak Ridge, TN 37831-0062 

phone: 865.576.8401 

fax: 865.576.5728 

email: mailto:reports@adonis.osti.gov

Available for sale to the public, in paper, from: 

U.S. Department of Commerce 

National Technical Information Service 

5285 Port Royal Road 

Springfield, VA 22161 

phone: 800.553.6847 

fax: 703.605.6900 

email: orders@ntis.fedworld.gov

online ordering: http://www.ntis.gov/ordering.htm

Printed on paper containing at least 50% wastepaper, including 20% postconsumer waste
NUMERICAL SIMULATIONS OF BOILING JET IMPINGEMENT COOLING IN POWER 

ELECTRONICS 

Sreekant Narumanchi1*

, Andrey Troshko,2 Vahab Hassani,1 and Desikan Bharathan1

1. National Renewable Energy Laboratory, 

1617 Cole Blvd, 

Golden, CO, USA 80401 

*

Phone: (303)275-4062 

Fax: (303)275-4415 *

Email: sreekant_narumanchi@nrel.gov

2. Fluent Inc., 

10 Cavendish Court, Centerra Park Resource, 

Lebanon, NH, USA 03766 

ABSTRACT 

Boiling jet impingement cooling is currently being explored to 

cool power electronics components. In hybrid vehicles, 

inverters are used for DC-AC conversion. These inverters 

involve a number of insulated gate bipolar transistors (IGBTs), 

which are used as on/off switches. The heat dissipated in these 

transistors can result in heat fluxes of up to 200 W/cm2

, which 

makes the thermal management problem quite important. 

In this paper, turbulent jet impingement involving nucleate 

boiling is explored numerically. The framework for these 

computations is the CFD code FLUENT. For nucleate boiling, 

the Eulerian multiphase model is used. A mechanistic model 

of nucleate boiling is implemented in a user-defined function 

(UDF) in FLUENT. The numerical results for boiling water 

jets (submerged) are validated against existing experimental 

data in the literature. Some representative IGBT package 

simulations that use R134a as the cooling fluid are also 

presented. 

KEY WORDS: boiling, power electronics, numerical 

simulations, CFD, IGBTs 

NOMENCLATURE 

21 ,CC εε Constants 

CP Specific heat, J/KgK 

d Jet diameter, m 

D Target equivalent diameter, m 

f Bubble departure frequency, Hz 

F

r

 Force vector, N/m 

g

r Gravitational acceleration vector, m/s2 

G Term in the turbulent kinetic energy equation, kg/ms3

h Heat transfer coefficient, W/m2K 

IGBT Insulated-gate bipolar transistor 

k Thermal conductivity, W/mK 

k Turbulence kinetic energy, m2

/s2 

L Latent heat, J/kg 

m& Mass transfer, kg/m3

s 

n Nucleation site density 

p Pressure, N/m2

Pr Prandtl number 

, qq ′′ Heat flux, W/m2

q

r

 Heat flux vector, W/m2

Q Interfacial energy exchange, W/m3

R

r

 Interfacial drag force, N/m3

Re Reynolds number 

S Source term in energy equation, W/m3

k S Source term in turbulent kinetic energy equation, 

kg/ms3

ε S Source term in dissipation rate equation, kg/ms4

u,U Liquid velocity, m/s 

v Phase velocity, m/s 

Jet velocity at the nozzle exit, m/s 

v

r

 Velocity vector, m/s 

Greek symbols 

α Phase volume fraction 

ΔT Temperature difference, K 

ε Dissipation rate, m2/s3 

κ Thermal conductivity, W/mK 

μ Dynamic viscosity, Ns/m2

ν Kinematic viscosity, m2/s 

ρ Density, Kg/m3 

σ Surface tension, N/m 

τ Shear stress, N/m2

Subscripts 

avg Average 

CHF Corresponding to critical heat flux 

d Jet diameter 

f Corresponding to fluid 

l Corresponding to liquid 

1
lv Corresponding to interaction between liquid and 

vapor 

pq Corresponding to interaction between phases p and q 

q Corresponding to phase q 

sat Corresponding to saturation 

sub Corresponding to subcooling 

v Corresponding to vapor 

INTRODUCTION 

Single-phase liquid jets have been studied very extensively in 

the literature [1-4]. These studies include experiments, 

theoretical analyses, and numerical simulations. Considerable 

attention has also been focused on boiling jets. Boiling liquid 

jets take advantage of the latent heat of vaporization of the 

fluid to provide fairly high heat transfer coefficients (>20,000 

W/m2

K), which makes them attractive for electronic cooling 

applications. The boiling curve for a saturated liquid is shown 

in Fig. 1. 

Typically, for electronic cooling applications, nucleate boiling 

is the preferred regime of operation because a small increase 

in wall superheat is accompanied by a large increase in the 

wall heat flux dissipated. Also, in electronics, it may not be 

possible to afford very large temperature differences between 

the solid surfaces and the liquid—a characteristic essential for 

regimes such as film boiling. 

In the context of boiling liquid jets, extensive work has 

already been reported in the literature [5-9]. Many studies 

have been carried out with circular [10-17] as well as planar 

[18-21] jets in both free-surface and submerged 

configurations. This includes single and multiple jets [14, 22-

26]. 

In the nucleate boiling literature, most of the correlations are 

cited in the following form: 

 (1) m

sat Δ= TCq sat

′′

where C and m are determined by curve fit to the experimental 

data, ΔTsat = Twall – Tsat is the wall superheat, Tsat is the 

saturation temperature of the fluid, Twall is the wall 

temperature, and is the wall flux. qsat ′′

Most of the heat transfer data are cited in the form given in 

Eq. 1, which can be rewritten as: 

 m

sub C

q T

q h /1

⎟

⎠

⎞ ⎜

⎝

⎛ ′′ +Δ

′′ = (2) 

where h is the heat transfer coefficient, ΔTsub = Tsat – Tf , 

where Tf is the fluid temperature, is the amount of subcooling in the fluid. 

Nucleate boiling is governed by intense bubble motion and 

mixing, so it is a strong process that does not depend on many 

jet parameters, unlike single-phase jets. Jet diameter, jet 

orientation, number of jets, jet configuration (free-surface or 

submerged), and even jet velocity do not have much effect on 

the heat transfer in nucleate boiling [5]. 

The target surface plays a critical role in the bubble nucleation 

process [27]. In fact, much of the difficulty in obtaining truly 

non-dimensional correlations for nucleate boiling arises from 

this. Surface conditions, surface aging, and even the condition 

of the surface during the course of an experiment [10] all have 

a considerable impact on the heat transfer results. 

The other aspect that has been given considerable attention is 

the critical heat flux (CHF) (Fig. 1) [5, 8]. When CHF occurs, 

the temperature of the wall shoots up because of dry-out 

conditions in which no liquid is in contact with the surface to 

sustain boiling. A schematic of this phenomenon is shown in 

Fig. 2 [5]. The liquid sub-layer (Fig. 2) drawn from the main 

liquid jet supply sustains the boiling process. When liquid 

cannot be supplied to this sub-layer, dry-out occurs, and CHF 

is reached. 

Considerable work has been done to develop non-dimensional 

correlations that show the dependence of the CHF on other 

parameters [13, 25, 28-32]. Typically in the literature, 

empirical correlations are presented only for certain simple 

geometries. For thermal design that involves more 

complicated geometries, it is important to have CFD modeling 

capability. Although multiphase models have been presented 

in the literature [33, 34], CFD modeling of boiling jets is still 

in its infancy [35, 36]. Numerical studies of phenomena such 

as nucleate pool boiling [37] and film boiling [38-40] are just 

beginning to appear in the literature. 

A significant amount of experimental work has also been 

reported in the literature on spray cooling (e.g. [41-43]). Both 

sprays and jets have their advantages and disadvantages. For 

some applications, sprays might be more suitable, while for 

others, jets may be a better choice. 

This paper presents CFD modeling of jets involving nucleate 

boiling. It is not possible to model boiling jets in the 

commercially available version of FLUENT, so a user-defined 

function (UDF) was used to perform these simulations. For 

nucleate boiling, the Eulerian multiphase model is used. A 

mechanistic model of nucleate boiling is implemented in a 

UDF in FLUENT. The numerical predictions are validated 

against experimental studies on submerged jets involving 

nucleate boiling. These experimental studies involve water 

and R-113 as the fluids. After this, IGBT package simulations 

are reported with a submerged boiling jet of water. 

To the best of our knowledge, these validations and IGBT 

package simulations with boiling jets are being reported for 

the first time. A comparison between single-phase and boiling 

jets from the heat transfer viewpoint and in the context of 

cooling the IGBT package is also presented. 

2
Fig. 1 General boiling curve for saturated liquids 

Fig. 2 Mechanism by which critical heat flux occurs

CFD MODELING OF JETS IN THE NUCLEATE 

BOILING REGIME: VALIDATION WITH 

EXPERIMENTS 

In this section , we present CFD modeling of jets in the 

nucleate boiling regime. 

Eulerian multiphase model description. The following are 

equations for the conservation of mass, momentum, and 

energy for the Euler multiphase model [33, 34] implemented 

in FLUENT 6.2.

Mass conservation equation for phase q: 

()( ) ∑=

=⋅∇+

∂

∂ n

p

qq qqq mv pq t 1

& r ραρα (3) 

where α is the phase fraction of phase q, ρ is the density, v

r

is 

the velocity vector and is the mass exchange between 

phases p and q. 

m&

Momentum conservation equation for phase q: 

( )( )

( )( ) qvmqliftqqq

n

p

pq pqpq

qqq q qqq qqqq q

vmR FFF

v vv p g t

, ,

1

rrr r & r

r rr r

+++ ++

⋅∇+ +⋅∇+∇−=

∂

∂

∑=

ρα

ρα ρα ρατα (4) 

where τ is the shear stress, R

r

 is the interfacial drag force, 

Fq

r is the turbulent diffusion force, is the lift force, p is 

the pressure, and 

F ,qlift

r

g

r is the gravitational acceleration vector. 

Energy conservation equation for phase q: 

( )( )

∑( ) =

+++∇

−∇+

∂

∂ ⋅∇+ −= ∂

∂

n

p

qq pq pqpq

q q qqq qqqq q

hmQSq

v

t

p h hv

t

1

:

& r

r r ρα ρα τα

 (5) 

where h is the enthalpy, q

r

 is the heat flux vector, S is the 

source term, and Q is the energy exchange term between the 

different phases. 

The following are restrictions on interfacial mass, momentum, 

and energy interfacial exchange terms: 

= − = 0 , pq mmm ppqp &&& (6) 

=−= 0 , pq RRR ppqp

r r r (7) 

=−= 0 , pq QQQ ppqp (8) 

Most sub-cooled boiling flows are turbulent, so the mixture 

phase k-epsilon model is used: 

() ( )

km

mk

k

mt

m mm

S

kvk Gk

t

+−

+⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ ∇⋅−∇=⋅∇+

∂

∂

ερ

μ

ρρ ,

,

Pr

r

 (9) 

where k is the turbulent kinetic energy, μ is the viscosity, ε is 

the dissipation rate, G and S are source terms in the turbulent 

kinetic energy equation. 

() ( )

( ) ε ε ε

ε

ερ

ε

ε

μ

ερερ

SCGC

k

v

t

mk m

mt

m mm

+ +−

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ ∇⋅−∇=⋅∇+

∂

∂

,1 2

,

Pr

r

 (10) 

where are constants, and S 21 ,CC εε ε is a source term in the 

dissipation rate equation. 

What follows is a closure for interfacial terms in the subcooled boiling model. 

Mass conservation equation. 

Rate of vapor formation per unit of volume in Eq. 1 becomes: 

( )

/( ) ( ) 0,max

/ 1

wE pl ls

isllvlv

n

p

qp

TTCLAq

LATThmm

+ ′′ + −

∑ −== =

&& (11) 

where is the interfacial heat transfer coefficient calculated 

from the Ranz-Marshall correlation; 

lv h

(( ) ( )) i sv v vsv −−= ααα // d116A (11a) 

is the interfacial area density, where ( 250 ) sv v α = α .,min

(Kurul and Podowski [33]); is the evaporating heat flux 

calculated from the RPI model [33, 34]; is the 

latent heat per unit of mass; 

qE

′′

00 −= hhL lsvs

( w wxxA

rr δ −= ) is interfacial area 

density of wall surface; and d is the diameter of the 

3
secondary phase (vapor bubble). Here, subscripts l , v , and s

mean liquid phase, vapor phase, and saturation state, 

respectively. In discretized form, becomes the ratio of the 

cell face area constituting the wall to the cell volume. The rate 

of condensation is calculated from Eq. 11. 

Aw

mvl &

Momentum conservation equation. The interfacial drag 

force per unit of volume is calculated as:

lv vrrvld R 75.0 C / dvv rr r αρ ⋅⋅⋅⋅⋅= (12) 

where the drag coefficient is calculated as: Cd

),( vis

d

dis

d = d CCMINC (13) 

where and are known drag correlations calculated 

for distorted and viscous regimes accounting for the high 

concentration effect: 

dis Cdvis Cd

( ) ( ( ) 75.0 2 Re1.01 5.0, Re

24 −

b ⋅+= l

b

vis Cd MAX α ) (14) 

( )( ( ) ( )) 2/1 1 / 5.0, 3

2 − = − ⋅ vl l

dis v

d g MAX d C σρρ α

6000,)2.012.0( 1

6000,0767.0

36000/ 073/

e

C ee e

l

φ

φ

φ

φ φ

 (15) 

The lift-force coefficient is calculated as (Moraga et al. [44]): 

⎜

⎜

⎜

⎝

⎛

+≥−

−− +<<

≤

= − +

059.1,002.0

9. e 05 (16) 

where φ = ReRe vb . This lift coefficient combines the 

opposing actions of two lift forces. “Classical” aerodynamics 

lift force results from interaction between bubble and liquid 

shear, and lateral force results from interaction between 

bubbles and vortices shed by bubble wake. Here, 

lrvb Re vd /ν r = is the bubble Re, and llvv Re vd /ν 2 r ×∇= is 

the bubble shear Re. 

Turbulent diffusion force is calculated as [33]: 

lv vlTD kCFF ∇−=−= αρ rr (17) 

where turbulent dispersion coefficient = 0.1 was used. CTD

Energy conservation equation. In a number of scenarios, it is 

reasonable to assume that the temperature of vapor inside 

bubbles is close to saturation. Ideally, the energy equation for 

the vapor should be solved. However, as a simplification, the 

vapor energy conservation equation is not solved, and the 

temperature of vapor mass is set to saturation. The source term 

due to wall heating goes entirely into the liquid phase:

qAS wwl ⋅= ′′ (18) 

According to the RPI model, the total heat flux from wall to 

liquid is partitioned into three components: 

EQlw′′ = ′′ + ′′ + qqqq ′′ (19) 

which are liquid convective heat flux, quenching heat flux, 

and evaporative heat flux. Under sub-cooled boiling 

conditions, the wall surface is subdivided into portion 

, covered by nucleating bubbles, and portion 

, covered by fluid. Therefore, convective heat flux is 

expressed as: 

Ω ( ≤Ω≤ 10

′′ (1TThq Ω−⋅−⋅= ) cell

lwlwl )( (20) 

where , the single-phase heat transfer coefficient, is 

derived from either log law if flow is logarithmic or Fourier 

law if flow is laminar. Liquid phase properties must be used 

while calculating for either turbulent or laminar flow. 

hlw

hlw

Quenching heat flux Qq′′ models additional energy transfer 

related to liquid filling the wall vicinity after the bubble 

detachment: 

( ) ( ) cell qQ Ω= plll − TTCf lw ′′ − 5.0 5.0 2 ρκπ (21) 

where f is the bubble departure frequency, κ is the thermal 

conductivity, C is the specific heat, and ρ is the density. 

Evaporative heat flux is given by: 

LfndqE ρ vvw

π 3

6 ′′ = (22) 

where d is the bubble departure diameter, n is the nucleation 

site density, and L is the latent heat. 

In Eqs. 20–22, closure must be provided for wall boiling 

parameters. Correct prediction of bubble departure diameter 

 is very important because evaporation heat rate depends 

strongly on this parameter, according to Eq. 22. It can be 

calculated from the following relation (Unal [45], Wei and 

Morel [46]): 

dvw

( ) 709.05 5.0 1042.2 − − dvw ⋅⋅⋅⋅= bap θ (23) 

where all variables are in SI units. In Eq. 23: 

( )

π

κρ

ρ

spss

v

sw TT C

a

2

− = (24) 

( ) ( ) ( )⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ ′′ − − =

lpll

w

ls

lv . ρ UC

q b TTMAX

00650 , /12

1

ρρ

 (25) 

( 0.1 ,61.0/ ) θ = UMAX l (26) 

where subscript s denotes the solid material of the wall where 

boiling occurs, U is the free stream velocity, and p is the 

pressure. In Eqs. 25 and 26, liquid temperature and velocity 

are local values at cells next to the wall. Because FLUENT 6.2 

is an unstructured solver, the calculation of these parameters is 

impossible in the general case. Again, these correlations were 

originally developed for use in structured or one-dimensional 

codes. So instead of using local cross-section averaged values, 

local values at the cell next to the wall are used. This means 

that bubble departure diameter will be somewhat 

overestimated, but bubble departure diameter is not very 

sensitive to local sub-cooling, as follows from Eqs. 23 and 25. 

According to the experimental data of Prodanovic et al. [47], 

Eqs. 23–26 give reasonable agreement with experiments for 

water boiling at atmospheric pressure. 

Nucleation site density is given by the following relation 

(Podowski [34]): 

( ( )) 80.1 200 −= TTn satw (27) 

Bubble departure frequency is calculated as: 

⎟

⎟

⎠

⎞

⎜

⎜

⎝

⎛ − =

lvw

vl

d

g f

ρ

ρρ

3

(4 (28) )

1 Ω−

The effective wall area occupied by boiling sites is given by: 

4
( 0.1,25.0min ) 2 ⋅=Ω vwnd ηπ (29) 

where η (−⋅= Ja 80/exp8.4 ), and Jacob number is given by 

(Kenning and Victor [48]) ( ) ( )

−1 ρ ⋅−⋅= ρ vlslpl LTTCJa . 

Bubble diameter in free stream is given by either a constant 

value or by Unal’s correlation as a function of local subcooling (Kurul and Podowski [33]): sub −= TTT lsat

⎜

⎜

⎜

⎜

⎝

⎛

<⋅

<<⋅−⋅

>⋅

=

−

− −

−

0,105.1

5.130,10105.1

5.13,105.1

3

3 4

4

sub

sub sub

sub

v

T

KTT

KT

d (30) 

Turbulence equations. The conventional mixture k-epsilon 

model contains two additional terms that describe additional 

bubble stirring and dissipation: one in turbulent kinetic energy 

and one in dissipation rate equation:

k CS dv vrvld 75.0 / r 2

αρ ⋅⋅⋅⋅= (31) 

k

v

rd S

d

vC

CS

r 3 = εε 3 (32) 

where 45.0 (Troshko and Hassan [49]). Cε 3 =

It is important to validate models against experimental data. 

We believe not many validations-of CFD models with 

experimental data on boiling jets-have been reported in the 

literature. Therefore, the code is validated with two 

experimental studies on submerged boiling jets. The effect of 

aspects, such as jet orientation and nozzle diameter, is also 

explored numerically. These aspects are covered in the 

following sections. 

Validation With Experimental Study of Katto and 

Kunihiro [11]. First, we examined the experimental study of 

Katto and Kunihiro [11]. The domain and the boundary 

conditions are shown in Fig. 3. A water jet with 3°C subcooling at atmospheric pressure (i.e., with Tinlet = 97°C) 

impinges on a 10-mm-diameter disk with an inlet velocity of 2 

m/s. The baseline nozzle diameter is 1.6 mm, and the distance 

between the nozzle exit and the heated plate is maintained at 3 

mm. A heat flux is imposed on the hotplate surface, as shown 

in Fig. 3. This is a submerged jet configuration. An 

axisymmetric domain is established. The RNG k-epsilon 

model with standard wall functions is used. With the use of 

the standard wall function, the y+ close to the walls should be 

maintained above 11. 

All the results presented here are expected to be meshindependent to within 5%. The properties of water at 1 

atmosphere pressure are listed in Table 1. Figure 4 shows the 

comparison of the wall superheats from the experiments and 

CFD modeling. 

There is ambiguity about the exact temperature measured in 

the experiments. All indications are that the stagnation point 

temperature is reported in the experiments. Figure 4 shows the 

stagnation as well as the average wall superheat. The 

experimental data are close to the stagnation superheat (within 

10%), which is encouraging. Given the nature of this problem, 

uncertainties on the order of even 30% are acceptable. There is 

temporal fluctuation in all quantities—such as temperature, 

fluid volume fraction, mass flow rates, and energy transfer 

rates—obtained from the CFD simulations. A temporal 

average is reported here. 

Figure 5 indicates the vapor volume fraction for two cases: 

heat flux of 50 W/cm2

 (left) and 100 W/cm2 (right). For the 

100 W/cm2

, the vapor volume fractions are much higher, as 

expected. Figure 6 lists the target surface superheat for two 

cases of 50 W/cm2

 (left) and 100 W/cm2 (right). As expected, 

the superheat for the 100-W/cm2

 case is higher than that for 

the 50-W/cm2

 case. Also, the stagnation point is associated 

with the highest temperature, whereas away from the 

stagnation point, the temperature is almost uniform. Figure 7 

shows the liquid temperature in the domain for the two cases 

of 50 W/cm2

 (left) and 100 W/cm2 (right). The liquid 

temperature rises as it flows along the hotplate and reaches a 

maximum temperature almost at the edge of the plate. The 

stagnation point is not associated with any vapor formation 

and accumulation because the liquid is not heated up, and 

vapor condenses in the liquid. However, away from the 

stagnation point, as the liquid is heated up, the vapor 

formation is sustained. This is evident in Fig. 5. 

Table 1. Properties of water and R-113 at 1 atmosphere 

pressure (1.013e+05 Pa)

Water R-113 

Liq. Vap. Liq. Vap. 

Saturation 

temperature (°C) 

100 47.6 

Surface 

tension 

(N/m) 

0.059 0.014 

Latent heat 

(J/kg) 

2,257,000 144,000 

Density, 

kg/m3

958 0.6 1,507 7.5 

Specific 

heat, J/kg-K 

4,219 2,010 980 724 

Dynamic 

viscosity, 

Ns/m2

2.83e

-04 

1.23e

-05 

5.23e04 1.08e

-05 

Thermal 

conductivity

, W/mK 

0.68 0.025 0.074 0.01 

5
Fig. 3 Domain used for the Katto and Kunihiro [11] validation 

study 

Fig. 5 Vapor volume fraction for two cases: 50 W/cm2

 (left) 

and 100 W/cm2

 (right) 

Fig. 6 Wall superheat for two cases: 50 W/cm2

 (left) and 100 

W/cm2

 (right) 

Fig. 4 Comparison of boiling curve from experiments (Katto 

and Kunihiro [11]) and CFD modeling 

6
Fig. 7 Liquid temperature in the domain for two cases: 50 

W/cm2

 (left) and 100 W/cm2 (right) 

Effect of jet orientation. Typically, experiments have 

revealed that jet orientation does not affect nucleate boiling [5, 

13]. To study the effect of jet orientation with respect to 

gravity, in Fig. 3 [11] the direction of gravity was reversed. In 

one case, the direction of vapor rise is against gravity (the 

baseline case); in the other case, the direction of vapor rise is 

in line with gravity. The results are summarized in Table 2. 

For both cases (with 100 W/cm2

), there is not much difference 

in the average wall superheat and the stagnation superheat, 

which is in line with experimental observations. However, for 

the case in which the vapor rise is in line with gravity, the 

vapor volume fraction increases because gravity causes the 

vapor to get locked up near the walls. The small difference in 

liquid temperature may also be related to this effect.

Effect of nozzle diameter. The Katto and Kunihiro study [11] 

indicates that the jet diameter does not have much effect on 

the boiling curve. Here, the code is tested with this 

experimental observation. For the case of 100 W/cm2

, three 

nozzle diameters are used (Table 3). The 1.6- and 1.16-mm 

cases are comparable in terms of average wall superheat, 

liquid temperature, and vapor content. With the 1.6-mm 

nozzle, the velocity is 2 m/s; with the 1.16-mm nozzle, the 

velocity is 2.5 m/s; and with the 0.71-mm nozzle, the velocity 

is 2.4 m/s. It is important to note that mass flow rate for three 

cases is not the same. The CFD results for the three cases 

(Table 3) match well with the experimental data (within 20%). 

Table 2. Effect of gravity on the thermal predictions for the 

case of 100 W/cm2

 (Katto and Kunihiro [11]) 

Table 3. Effect of nozzle diameter on the thermal predictions 

for the case of 100 W/cm2

 (Katto and Kunihiro [11]) 

Vapor rise 

direction 

against 

gravity 

Vapor rise 

direction in 

line with 

gravity 

Average wall 

superheat, K 

12.6 12.6 

Average liquid 

temperature, K 

370.9 370.8 

Average vapor 

content 

0.39 0.49 

Stagnation 

superheat, K 

19 19 

Nozzle diameter 1.6 

mm 

1.16 mm 0.71 

mm 

Average wall 

superheat, K 

12.6 12.6 13.86 

7
Impact of mesh refinement. The original mesh contained 

2312 cells. The mesh was refined so that the total number of 

cells increased to 9248. The average wall superheat was found 

to be 12.7 K. For the original mesh, which had 2312 cells, the 

wall superheat was 12.6 K. This suggests that the results from 

the original mesh are mesh-independent to well within 2%. 

IGBT PACKAGE SIMULATIONS WITH BOILING 

JETS 

The previous section established some confidence in the CFD 

predictions of FLUENT by comparing them with experimental 

data from the literature. This section explores IGBT package 

simulations with boiling jets. 

The axisymmetric domain used in the simulation is shown in 

Fig. 8. The IGBT package consists of several layers. The 

silicon die is mounted on a copper layer, followed by an 

aluminum nitride layer, another copper layer, and finally the 

aluminum baseplate/heat sink. Further details on the IGBT 

structure can be found in reference [50]. These boiling 

simulations are fairly intensive and take several hours to 

converge on 64-bit Linux machines with 1.4 GHz processor 

speed-even with a small spatial mesh (i.e., small number of 

cells). 

In the automotive industry, R134a is the working fluid used in 

the airconditioning units. Here we explore the possible use of 

R134a as the cooling fluid for the IGBTs. The simulations 

presented here are performed with an R134a (properties given 

in Table 4) jet inlet temperature of 47°C. For the boiling 

simulations, the pressure is maintained at 1.318e+06 Pa (i.e. ~ 

13 atmospheres), at which the saturation temperature of the 

R134a is 50°C. So there is a 3°C sub-cooling in the R134a 

temperature at the inlet. A volumetric heat generation term is 

included in the silicon layer to simulate heat dissipate. The 

heat dissipation in the silicon die is 25 W/cm2

. This 

corresponds to a heat dissipation of 20.3 W. The jet inlet 

velocity is 2 m/s (Fig. 9). Figures 9 and 10 show the velocity 

and temperature contours, respectively, in the domain. A 

maximum temperature of 56°C is attained. The maximum 

temperatures are in the silicon die, as expected. The domain 

used for this simulation contains 5506 cells. A simulation was 

performed with 22024 cells also. The results for the maximum 

temperature in the domain within this refined mesh were 

within 2% of the results from the original baseline mesh The 

vapor volume fraction in the domain is shown in Fig. 11. The 

volume-averaged vapor fraction for this case is 17% (Table 5). 

Figure 11 also gives a hint that in this reduced IGBT structure, 

the vapor does not have a very clear removal path. It tends to 

get trapped to some extent because of the confining copper 

and aluminum walls.

 

Copper

Aluminum

Fig. 8 Axisymmetric domain used for the IGBT package 

simulation 

 

Fig. 9 Velocity contours in the domain 

 

 

Fig. 10 Temperature contours in the domain – 25 W/cm2

8
 

 

Fig. 11 Contours of vapor volume fraction – 25 W/cm2

It would be illustrative to compare the performance of singlephase and boiling jets in the context of the IGBT package 

cooling. Figure 12 shows such a comparison for the case of 25 

W/cm2

. The plot on the left in Fig. 12 is the same as the plot in 

Fig. 10, whereas the plot on the right in Fig. 12 is the case in 

which there is no boiling, and everything else is identical to 

the case mentioned in this section. Table 5 also has a summary 

of all the results with and without boiling. For the cases not 

involving any boiling, the bubble nucleation site density is 

equal to zero. In all other aspects, they are similar to the case 

involving boiling. This facilitates a direct comparison between 

the results of the two cases. The maximum temperature in the 

domain (silicon die) is 69.8 °C. The case involving boiling 

therefore gives maximum temperatures almost 13.8°C lower 

than the case not involving any boiling (Table 5). This is a 

significant difference. These results hint at the advantage that 

heat of the fluid, boiling may be yielding lower die 

temperatures than comparable single-phase flows. 

Of course, these numerical results need to be validated 

experimentally. The validations presented in Section 2 give 

some degree of confidence in the predictions from the code. It 

would be interesting to see how experimental results from the 

IGBT package involving boiling would match up with the 

predictions presented here. Overall, there is evidence of the 

benefits of boiling from a heat transfer standpoint. It is 

important to note that this model is applicable only in the 

bubbly flow regime. This implies that the liquid should be the 

dominant phase and the vapor bubbles should be the 

secondary phase. If there is intense boiling with a large 

amount of vapor formation, the basic model may not be 

applicable. Such regimes remain an area of research. Although 

experiments are essential, modeling can yield very useful 

information when used within their range of applicability. 

Tools such as those presented in the paper can be explored for 

thermal design involving the use of R134a or any other fluid 

for thermal management in power electronics as well as other 

applications. 

Summary & Conclusions 

The CFD model and code are validated against experimental 

studies involving submerged jets. To the best of our 

knowledge, this is the first time such validations are being 

reported. A reasonable match is found between the 

experimental boiling curves and those obtained by CFD. IGBT 

package simulations suggest that, for the case examined here, 

boiling jets are providing significant benefits over non-boiling 

jets. A tool has been established that can be used for thermal 

design of two-phase cooling systems in the power electronics 

context. 

Fig. 12: Temperature contours in the domain with boiling (left) 

and without boiling (right) – 25 W/cm2

9
Table 4. Properties of R134a at 1.318e+06 Pa, saturation 

temperature = 323.15K 

Liquid Vapor 

Density, kg/m3 1103.4 66.16 

Specific heat, J/kgK 1575 1219 

Thermal 

conductivity, W/mK 

0.0735 0.018 

Dynamic viscosity, 

kg/ms 

0.000167 0.0000138 

Table 5. Summary of results obtained from IGBT package 

simulations 

Acknowledgments 

We acknowledge the support and funding provided by Susan 

Rogers, the program manager for the advanced power 

electronics and electric machines project under the U.S. 

Department of Energy FreedomCAR Program. 

References

[1] Webb, B.W. and Ma, C.-F., 1995, "Single-phase liquid jet 

impingement heat transfer", Advances in Heat Transfer, 26, 

pp. 105-217. 

[2] Lienhard, J.H., 1995, "Liquid jet impingement", in Annual 

Review of Heat Transfer, pp. 199-270. 

[3] Garimella, S.V., 2000, "Heat Transfer and flow fields in 

confined jet impingement", in Annual Review of Heat 

Transfer, pp. 413-494. 

[4] Narumanchi, S.V.J., Amon, C.H., and Murthy, J.Y., 2003, 

"Influence of pulsating submerged liquid jets on chip-level 

thermal phenomena", ASME Journal of Electronic Packaging, 

125, pp. 354-361. 

[5] Wolf, D.H., Incropera, F.P., and Viskanta, R., 1993, "Jet 

impingement boiling", Advances in Heat Transfer, 23, pp. 1-

132. 

[6] Zhou, D.W. and Ma, C.F., 2004, "Local jet impingement 

boiling heat transfer with R113", Heat and Mass Transfer, 40, 

pp. 539-549. 

[7] Mudawar, I., 2000, "Assessment of high-heat-flux thermal 

management schemes", in International conference on 

thermal, mechanics and thermomechanical phenomena in 

electronic systems, Las Vegas, NV, ASME, pp. 1-20. 

[8] Mitsutake, Y. and Monde, M., 2003, "Ultra high critical 

heat flux during forced flow boiling heat transfer with an 

impinging jet", ASME Journal of Heat Transfer, 125, pp. 

1038-1045. 

[9] Liu, Z.-H., Tong, T.-F., and Qiu, Y.-H., 2004, "Critical 

heat flux of steady boiling for subcooled water jet 

impingement on the flat stagnation zone", ASME Journal of 

Heat Transfer, 126, pp. 179-183. 

[10] Copeland, R.J., 1970, "Boiling heat transfer to a water jet 

impinging on a flat surface (-1g)",Ph.D. thesis, Department of 

Mechanical Engineering, Southern Methodist University, 

Dallas, TX. 

[11] Katto, Y. and Kunihiro, M., 1973, "Study of the 

mechanism of burn-out in boiling system of high burn-out heat 

flux", Bull. JSME, 16(99), pp. 1357-1366. 

[12] Katto, Y. and Monde, M., 1974, "Study of mechanism of 

burn-out in a high heat-flux boiling system with an impinging 

jet", in 5th International Heat Transfer Conference, B6.2, pp. 

245-249. 

[13] Monde, M. and Katto, Y., 1978, "Burnout in a high heatflux boiling system with an impinging jet", Int. J. Heat Mass 

Transfer, 21, pp. 295-305. 

[14] Monde, M., 1980, "Burnout heat flux in saturated forced 

convection boiling with an impinging jet", Heat Transfer - 

Japanese Research, 9(1), pp. 31-41. 

[15] Ma, C.-F. and Bergles, A.E., 1983, "Boiling jet 

impingement cooling of simulated microelectronic chips", in 

Heat Transfer in Electronic Equipment, ASME, New York, 

HTD-Vol. 28, pp. 5-12. 

[16] Ma, C.-F. and Bergles, A.E., 1986, "Jet impingement 

nucleate boiling", Int. J. Heat Mass Transfer, 29(8), pp. 1095-

1101. 

[17] Ruch, M.A. and Holman, J.P., 1975, "Boiling heat 

transfer to a Freon-113 jet impinging upward onto a flat, 

heated surface", Int. J. Heat Mass Transfer, 18, pp. 51-60. 

[18] Katto, Y. and Ishii, K., 1978, "Burnout in a high heat flux 

boiling system with a forced supply of liquid through a plane 

jet", in Proceedings of the 6th International Heat Transfer 

Conference, 1, FB-28, pp. 435-440. 

[19] Miyasaka, Y., Inada, S., and Owase, Y., 1980, "Critical 

heat flux and subcooled nucleate boiling in transient region 

between a two-dimensional water jet and a heated surface", J. 

Chem. Eng. Jpn., 13, pp. 29-35. 

[20] Mudawar, I. and Wadsworth, D.C., 1991, "Critical heat 

flux from a simulated chip to a confined rectangular impinging 

jet of dielectric liquid", Int. J. Heat Mass Transfer, 34(6), pp. 

1465-1479. 

[21] Wadsworth, D.C. and Mudawar, I., 1990, "Cooling of a 

multichip electronic module by means of confined twodimensional jets of dielectric liquid", ASME Journal of Heat 

Transfer, 112, pp. 891-898. 

[22] Nonn, T., Dagan, Z., and Jiji, L.M., 1988, "Boiling jet 

impingement cooling of simulated microelectronic heat 

sources", in Winter Annual Meeting, Chicago, Illinois, ASME, 

88-WA/EEP-3, pp. 1-9. 

[23] Nonn, T., Dagan, Z., and Jiji, L.M., 1989, "Jet 

impingement flow boiling of a mixture of FC-72 and FC-87 

liquids on a simulated electronic chip", in National Heat 

Transfer Conference, HTD-Vol. 111, pp. 121-128. 

25 W/cm2

With 

boiling 

Without 

boiling 

Tmax (°C) 56 69.8 

Volume averaged vapor 

fraction 

0.17 0 

Volume averaged liquid 

temperature (°C) 

47.3 47.3 

10
[24] Monde, M., Kusuda, H., and Uehara, H., 1982, "Burnout 

heat flux in saturated forced convection boiling with two or 

more impinging jets", Heat Transfer - Japanese Research, pp. 

18-31. 

[25] Monde, M. and Mitsutake, Y., 1996, "Critical heat flux in 

forced convective subcooled boiling with multiple impinging 

jets", ASME Journal of Heat Transfer, 117, pp. 241-243. 

[26] Pais, M.R., Chow, L.C., and Mahefkey, E.T., 1993, 

"Multiple jet impingement cooling", Journal of 

Thermophysics and Heat Transfer, 7(3), pp. 435-440. 

[27] Rohsenow, W.M., 1952, "A method of correlating heattransfer data for surface boiling of liquids", Trans. Am. Soc. 

Mech. Engrs., Paper No. 51-A110, 74, pp. 969-976. 

[28] Sharan, A. and Lienhard, J.H., 1985, "On predicting 

burnout in the jet-disk configuration", ASME Journal of Heat 

Transfer, 107, pp. 398-401. 

[29] Katto, Y. and Yokoya, S., 1988, "Critical heat flux on a 

disk heater cooled by a circular jet of saturated liquid 

impinging at the center", Int. J. Heat Mass Transfer, 31(2), pp. 

219-227. 

[30] Monde, M., 1985, "Critical heat flux in saturated forced 

convective boiling on a heated disk with an impinging jet", 

Warme-und Stoffubertragung, 19, pp. 205-209. 

[31] Monde, M., 1987, "Critical heat flux in saturated forced 

convection boiling on a heated disk with an impinging jet", 

ASME Journal of Heat Transfer, 109, pp. 991-996. 

[32] Monde, M. and Inoue, T., 1991, "Critical heat flux in 

saturated forced convective boiling on a heated disk with 

multiple impinging jets", ASME Journal of Heat Transfer, 

113, pp. 722-727. 

[33] Kurul, N. and Podowski, M.Z., 1990, "Multidimensional 

effects in forced convection subcooled boiling", in Ninth 

International Heat Transfer Conference, Jerusalem, Israel, 1-

BO-04, pp. 21-26. 

[34] Podowski, R.M., et al., 1997, "A mechanistic model of 

the ebullition cycle in forced convection subcooled boiling", 

in Eighth International Topical Meeting on Nuclear Reactor 

Thermal-Hydraulics, Kyoto, Japan, 3, pp. 1535-1542. 

[35] Wang, D., Yu, E., and Przekwas, A., 1999, "A 

computational study of two phase jet impingement cooling of 

an electronic chip", in Fifteenth IEEE Semi-therm 

Symposium, IEEE, pp. 10-15. 

[36] Seiler, N., et al., 2003, "Modeling and computation of 

heat exchanges in the configuration of an impinging jet on a 

hot plate", in International conference on supercomputing in 

nuclear applications, SNA 2003, Paris, France, pp. 1-13. 

[37] Mukherjee, A. and Dhir, V.K., 2004, "Study of lateral 

merger of vapor bubbles during nucleate pool boiling", ASME 

Journal of Heat Transfer, 126, pp. 1023-1039. 

[38] Agarwal, D.K., et al., 2004, "Planar simulation of bubble 

growth in film boiling in near-critical water using a variant of 

the VOF method", ASME Journal of Heat Transfer, 126, pp. 

329-338. 

[39] Esmaeeli, A. and Tryggvason, G., 2004, "Computations 

of film boiling. Part I: numerical method", International 

Journal of Heat and Mass Transfer, 47, pp. 5451-5461. 

[40] Esmaeeli, A. and Tryggvason, G., 2004, "Computations 

of film boiling. Part II: multi-mode film boiling", International 

Journal of Heat and Mass Transfer, 47, pp. 5463-5476. 

[41] Estes, K.A. and Mudawar, I., 1995, "Comparison of twophase electronic cooling using free jets and sprays", ASME 

Journal of Electronic Packaging, 117, pp. 323-332. 

[42] Lin, L. and Ponnappan, R., 2003, "Heat transfer 

characteristics of spray cooling in a closed loop", International 

Journal of Heat and Mass Transfer, 46, pp. 3737-3746. 

[43] Pautsch, A.G. and Shedd, T.A., 2005, "Spray 

impingement cooling with single- and multiple-nozzle arrays. 

Part I. Heat transfer data using FC-72", International Journal 

of Heat and Mass Transfer, 48, pp. 3167-3175. 

[44] Moraga, F.J., Bonetto, F.J., and Lahey, R.T., 1999, 

"Lateral forces on spheres in turbulent uniform shear flow", 

International Journal of Multiphase Flow, 25, pp. 1321-1372. 

[45] Unal, H.C., 1976, "Maximum bubble diameter, maximum 

bubble-growth time and bubble-growth rate during the 

subcooled nucleate flow boiling of water up to 17.7 MN/m2

", 

Int. J. Heat Mass Transfer, 19, pp. 643-649. 

[46] Wei, Y. and Morel, C., 2002, "Prediction of parameters 

distribution of upward boiling two-phase flow with two-fluid 

model", in Proceedings of ICONE 10, Arlington, VA, USA, 

10, pp. 

[47] Prodanovic, V.A., 2002, "Bubble behavior in subcooled 

flow boiling of water at low pressures and low flow rates", Int. 

J. Multiphase Flow, 28(1), pp. 1-19. 

[48] Kenning, D.B.R. and Victor, H.T., 1981, "Fullydeveloped nucleate boiling: overlap of areas of influence and 

interference between bubble sites", Int. J. Heat Mass Transfer, 

24(6), pp. 1025-1032. 

[49] Troshko, A.A. and Hassan, Y.A., 2001, "A two-equation 

turbulence model of turbulent bubbly flows", International 

Journal of Multiphase Flow, 27, pp. 1965-2000. 

[50] Narumanchi, S.V.J., Hassani, V., and Bharathan, D., 

2005, "Modeling single phase and boiling liquid jet 

impingement cooling in power electronics", National 

Renewable Energy Laboratory, NREL/TP-540-

38787,http://www.nrel.gov/docs/fy06osti/38787.pdf. 

11
REPORT DOCUMENTATION PAGE Form Approved 

OMB No. 0704-0188

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 

gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this 

collection of information, including suggestions for reducing the burden, to Department of Defense, Executive Services and Communications Directorate (0704-0188). Respondents 

should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a 

currently valid OMB control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ORGANIZATION.

1. REPORT DATE (DD-MM-YYYY)

December 2006 

2. REPORT TYPE

Conference paper 

3. DATES COVERED (From - To)

 

5a. CONTRACT NUMBER

DE-AC36-99-GO10337 

5b. GRANT NUMBER

4. TITLE AND SUBTITLE

Numerical Simulations of Boiling Jet Impingement Cooling in Power 

Electronics 

5c. PROGRAM ELEMENT NUMBER

5d. PROJECT NUMBER

NREL/CP-540-39401 

5e. TASK NUMBER

FC06.7010 

6. AUTHOR(S)

S. Narumanchi, A. Troshko, V. Hassani, and D. Bharathan 

5f. WORK UNIT NUMBER

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

National Renewable Energy Laboratory 

1617 Cole Blvd. 

Golden, CO 80401-3393 

8. PERFORMING ORGANIZATION 

REPORT NUMBER

NREL/CP-540-39401 

10. SPONSOR/MONITOR'S ACRONYM(S)

NREL 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

11. SPONSORING/MONITORING 

AGENCY REPORT NUMBER

12. DISTRIBUTION AVAILABILITY STATEMENT

National Technical Information Service 

U.S. Department of Commerce 

5285 Port Royal Road 

Springfield, VA 22161 

13. SUPPLEMENTARY NOTES

14. ABSTRACT (Maximum 200 Words)

This paper explores turbulent boiling jet impingement for cooling power electronic components in hybrid electric 

vehicles. 

15. SUBJECT TERMS

boiling jet impingement; power electronics; numerical simulations; CFD; IGBTs; power electronic cooling applications; 

hybrid electric vehicles 

16. SECURITY CLASSIFICATION OF: 19a. NAME OF RESPONSIBLE PERSON 

a. REPORT 

Unclassified 

b. ABSTRACT 

Unclassified

c. THIS PAGE 

Unclassified

17. LIMITATION 

OF ABSTRACT

UL 

18. NUMBER 

OF PAGES 

19b. TELEPONE NUMBER (Include area code)

Standard Form 298 (Rev. 8/98) 

Prescribed by ANSI Std. Z39.18