Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
2D Steady State Numerical Model for Critical TwoPhase Flow in a ConvergentDivergent
Nozzle with New Critical Criterions
Martel Sylvain*, Dostie Michelt and Mercadier Yves*
t Energy Technology Laboratory, HydroQuebec, Shawinigan, Canada
Department of Mechanical Engineering, Universite de Sherbrooke, Sherbrooke, Canada
2500 Boul. University, Sherbrooke, J1K 2R1, Canada
sylvain.martel2@cusherbrooke.ca
Keywords: Critical twophase flow, Critical criterion, Numerical model
Abstract
The use of twophase ejectors to improve refrigeration systems encounters today a great interest. However, modeling of such
devices with low void fraction at the entrance of the motive nozzle, presents significant challenges. The choking conditions
and the discontinuities appearing in a twophase flow in a nozzle are not well documented and some works are needed to
better anticipated flow behavior under these conditions. This paper presents a steady state twophase flow model including
new choking criterions for twodimensional conservative systems. The present twofluid model remains general since no
assumption on pressure, temperature or velocity equilibrium is made. For this reason, it allows the use of pressure correlation
between phases compared to those in literature which are dealing principally with singlepressure model. As a first step, this
model is used to study the flow in the motive nozzle of an ejector.
1. Introduction
The steady state flow of compressible fluid through
convergentdivergent nozzles covers various important
flow phenomena like the occurrence of critical flow
conditions, transition from subsonic to supersonic flow or
the occurrence of flow discontinuities. For singlephase
gas, relatively simple algebraic solutions exist as described
in many gasdynamic textbooks. For twophase flow
conditions, iterative algebraic solutions can be derived
only for simple onedimensional case with homogeneous
assumption and thermal equilibrium between phases. For
the more general case of nonhomogeneous and
nonequilibrium conditions, the numerical resolution
present a major challenge.
In fact, several mathematical models have been published
in the literature for twophase flow but only a few have
addressed the steady state behavior in a nozzle, especially
with low void fraction, heat and mass transfer between
phases, and phase interaction. In addition, twophase
critical flow has been the subject of many analytical and
experimental investigations, mainly because of it
importance in the safety analyses of pressurized water in
nuclear reactors and convergingdiverging nozzles.
One phenomena still of interest is the choking condition in
a critical flow such as in supersonic ejectors. Supersonic
ejectors are widely used in a range of applications such as
aerospace, propulsion, refrigeration and many thermal
systems. In refrigeration applications, ejectors could be
used as thermocompressor or to recover part of the work
that would be lost in the expansion valve. The flow at the
entrance of the ejector coming from the condenser unit of a
refrigeration system is either a subcooled or a saturated
liquid. Inside the motive nozzle, the flow undergoes
important pressure variation resulting in a very fast change
of phase (flashing process) giving turbulent twophase
flow with thermodynamical nonequilibrium. Even if
condensation and evaporation have been studied for
several years, many uncertainties about low void fraction
twophase flow and critical conditions remain.
In this paper, a twodimensional compressible steady state
twophase flow model is presented with new choking
criterions that are directly related to optimal flux
conditions developed by [1]. As a first step, this model is
used to study the flow in the motive nozzle of an ejector.
2. Theoretical model
Twofluid twophase flow models are formulated around
the macroscopic separate balance equations for each phase
based on space and time averaging of the local
instantaneous phasic flow equations. Consequently, the
model can provide information only on the average flow
behavior, which assumes that sufficiently accurate
empirical correlations can be used to describe heat, mass
and momentum transfer processes at the phasic interface
and at the boundary walls.
2.1 Governing equations
Although twopressure models show interesting aspects of
twophase flows none of them has yet reached a state of
Paper No
Nomenclature
C, specific heat J/kgK
E energy flux kgm3/s3
f phase fluxes scaling
F force N/m3
G critical mass flux kg/m2s
h enthalpy J/kg
hfg heat of vaporization J/kg
J impulse flux kgm2/s2
h mass flow rate kg/m2s
M mass flux kgm/s
P pressure Pa
R perfect gas constant J/kgK
s slip ratio
T temperature K
i vector of velocity m/s
U Norm of the velocity m/s
V Total volume m3
Zc compressibility factor
Greek letters
p density kg/m3
ak volumetric concentration of 
phase k
source term representing kg/m3s
interfacial transport process kg/m2s2
kg/ms3
Subscripts
exp experimental
k index for phases
th theoritical
Superscripts
mean mean value
D interfacial property
M mass
J momentum
E energy
Q heat
0 initial value
maturity for broader applications to scientific or technical
problems. This is the reason why in most of the actual
twophase flow models the assumption of pressure
equilibrium is introduced. This seems to be justified for
many technical applications as long as surface tension
effects can be neglected. Nevertheless, the present model
remains general and each phase can have a different
pressure. Neglecting the influence of the bulk heat
conduction and bulk viscous effects compared with the
governing effects resulting from the interfacial heat and
mass transfer processes, the following balance equations
are obtained:
Mass:
V *VIVa. I Vak (1)
Momentum:
VVaC p k2 +Pk)=Vack O(k) (2)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Energy:
V Vak puk khk + 2 =Vak (akE) (3)
2
where olM, oJ, and ok are the interfacial source terms for
mass, momentum and energy, respectively. V is the total
volume and ak the volume fraction of each phase.
The interfacial source terms for mass, momentum and
energy are defined as follow:
kM = idA
Vak AD
J = FD + UkMUkmean
ak k Mrk (4)
[ mean )2
E Q + FD D M k
ak =a +k k uk k hfg + 2
where aM correspond to phase change mass flow rate by
unit volume, Fk is the interfacial force per unit volume,
Ukmean is the average velocity of the phase just after the
phase change, oaQ is the volumetric heat source, hfg is the
heat of vaporization for phase changes, and ikD is the
interfacial velocity.
These interfacial source terms on the righthand side of
equations (1) to (3) distinguish between the internal
contributions describing the interfacial exchange between
phases. Their formulations are presented in [2]. These
equations are equivalent to those often found in the
literature. The only difference is related to the righthand
side of the equations where the properties transferred
between phases are explicitly introduced. The conservation
principles for mass, momentum, and energy require the
following balances of the corresponding source terms:
a =o0, FD = ZD =o (5)
k=g,l k=g,l k=g,l
As a result of the space and time averaging procedures,
information on local transport phenomena at the interface
and at the wall have been lost. This information has to be
substituted by appropriate modeling the source terms in the
balance equations. Generally three types of additional
terms might be distinguished: nondissipative interfacial
coupling terms, algebraic terms representing interfacial
exchange processes for mass, momentum and energy, and
physically dissipative terms such as bulk heat conduction
and turbulent viscosity.
The interfacial forces are usually split into viscous part and
nonviscous part. The viscous part represents the
interfacial frictional drag force resulting from different
local phasic velocities. The nonviscous part includes
additional forces due to virtual mass effects, and interfacial
pressure forces. These nondissipative terms, introduced to
compensate for the loss of information at the interface as a
result of the averaging processes, are entirely described by
space derivative of the primitive variables. These terms
might have considerable influence on the Eigenspectrum
of the basic flow equations of a transient model.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
2.2 Closure laws
To close de system of equations (four balance equations for
each phase), four additional equations are needed since
there are six unknown for each phase (P1, P2, hk, Ukl, Uk2, t,
pk). Two equations of state are used:
Pk = Zc k R/Tk (6)
where Tk is temperature of phase k defined as a function of
enthalpy assuming constant specific heat:
hk =h +Ck (TkTk) (7)
where hkO and TkO are constants.
In addition, the following constraint for the volumetric
phase concentration is used:
a + a2 = 1 (8)
Finally, a relation related to the topology of the flow whish
was lost in the averaging process is needed. This unknown
relation is often replaced by a correlation linking pressure
of each phase:
P2 = f() (9)
This correlation can take different forms where the
simplest are: equal local pressure or constant pressure
jump.
2.3 Critical conditions
The critical conditions for systems described by
Eqs.(1)(3) have been shown by Dostie and Mercadier [1]
to be related to local optimum values of the conservative
variables, or fluxes, with respects to the variation of the
primitive variables. These critical criterions provide a
simple tool to see if the flow is choked or not and they
express a mixture property. They have been shown to be
singular points of the trajectories of the dynamical
differential system associated with Eqs.(1)(3).
Considering a system of n phases, the first criterion is
based on a global scaling of the following mass, impulse,
and energy mixture fluxes:
M1 =Z VakokUkl
k
J = Voakpkk2 + VakP (10)
k k
Using relation (7), the energy flux becomes:
E1== hko+ ZCRk Pk +Uk2 akPkUkl (11)
kH ZckRk k 2 j
where Uk is the norm of the velocity of the phase k and uk1
is the velocity component of phase k along the coordinate
1.
By scaling the phase fluxes with respect fluxes of the
phase 1, the following global scaling factors are obtained:
fMk akpkl
fjk=
=akPU12 + akP
alplU12 + Clp1
0 Cpk Pk Uk2
hk ZckRk Pk 2
0 +i +
fE z 1R1 2
zR, p, 2
Subsequently, fluxes are:
M1 = ValplUl fMk
k
J,= ap U2 + Va fJk
E,= 0+
Cpl l+ V UZ f EkfMk
ZclR, p,1 2 k
With these fluxes, an explicit expression for momentum
flux J, as a function of U1, Mi and E, can be obtained. The
critical conditions can be obtained from:
a
SJ,l(Ul,M,,E,)=0
au1
pU12 1 (14)
CPl
where zli is the compressibility factor, R1 is the perfect gas
constant and Cp is the specific heat for the phase of
reference.
For a singlephase critical condition of perfect gas, this
equation resumes to the usual speed of sound. For a
singlephase liquid condition, in which case R is almost
zero, the right hand side is equal to on. This shows that a
critical flow criterion that applies to compressible as well
as to incompressible fluid is obtained. It is a local criterion
valid over the entire flow domain.
The next family of critical conditions is obtained by
providing scaling functions for the two terms in the
impulse mixture flux and the terms in the energy mixture
flux.
There is two terms to be scaled in the impulse flux:
2
afk p ukkl
f Jk 2
Q1'1 (15)
J k k Pk
f P,
aj~k
And three terms in the energy flux:
Paper No
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
akpkukl Pk Ck
f Pk ZckRk
jE k p ^
alplUl 1 P1
lEakkUklUk2
2p z
Critical conditions are then defined
equation:
P1u2 _
k
k
by the folio
1
ZfEA
Euk c
k ZclR1I
Z fEk cp1
k
This equation shows a new family of singular points which
defines critical conditions by scaling the evolutions of the
other phases to a phase of reference.
The third family is obtained by introducing the condition
that the flow must fill a prescribed volume V. From this,
the sum of the volume of all phases is equal to the
prescribed volume:
Vk
This prescription has a special nature as it does not link
two phases or does not apply to a particular phase but is a
global constraint on the flow mixture. The last critical
conditions are then defined by the following equation:
P12 1 (19)
Pl Y.f Y .fFE k
a k
k
k Zcl R
YfE k Cpl
As any phase can be taken as phase 1, Eqs.(14), (17) and
(19) provide a bundle of critical conditions that applied to
the steady state system being studied. These critical
conditions are mixture properties determined by the state
of phase 1 and how this state scales with the state of the
other phases. They provide simple algebraic relations that
define if the flow is locally subcritical, critical or
supercritical with respect to a given phase.
The information provided by these critical criterions is
taken into account in the pressure correction scheme and
the relaxation strategy to ensure convergence and stability
as the flow undergoes a transition between subcritical and
supercritical states.
3. Numerical Model
The model includes the turbulent friction as well as heat
and mass transfer phenomena affecting the behavior of the
twofluid twophase flow in a convergentdivergent nozzle
such as the primary nozzle of a refrigerant ejector. The
(16) model is twodimensional in direction perpendicular to the
flow. Velocity and pressure fields are solved with the
SIMPLE (SemiImplicit Method for PressureLinked
Equations) type pressurecorrector method developed for
steady state flow conditions. Correlations for the friction
phenomena and for the exchange area between phases are
proposed as a function of void fraction. The complete
wing model formulation allows the simulation of the complete
void fraction range from liquid to gas flow including
bubbly and droplet flows.
3.1 Algorithm and grid
A new steady state compressible multiphase solver in
OpenFOAM is under development. This solver will allow
physical compressibility, turbulence, heat and mass
transfer, phase interaction and different phase velocity. A
new choking criterion developed by [1] is used to build the
numerical scheme and ensure the stability and
convergence.
3.2 Flow physic
For the present study, the nozzle geometric design is
depicted in Fig. 1 and based on the ASTAR nozzle
geometry [2] to simplify validation of the model since
validated numerical results are already known.
Figure 1: Typical ASTAR nozzle geometry.
3.3 Boundary conditions
At the inlet of the nozzle, homogeneous equilibrium flow
is assumed from reservoir. For this reason, each phase have
the same temperature and pressure. The inlet volume
fraction is imposed. For the outlet boundary condition a
constant exit pressure is applied. The initial conditions in
the nozzle are identical with the upstream reservoir
conditions which imply that the transient calculation starts
with a strong discontinuity at the nozzle exit.
4. Comparison of critical conditions
These new critical criterions can be used to compare
resulted from previous works. In this section, comparisons
with previous numerical results obtained by [2] and with
previous experimental results showed in [4] and [5] are
presented.
4.1 Comparison of numerical critical conditions
Table 1 illustrates the interest of the new development
Paper No
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
presented in section 2.3. Numerical data for airwater
critical flow are presented by [2]. These data are
summarized in the upper part of Table 1 for some locations
in the nozzle for a flow with fMk = 1. The data were read
from Fig. 9.46 of [2] and the error associated with this
reading should be a few percent.
Table 1: Critical air velocity and mass flux as predicted
using Eq.(19). Numerical data taken from Fig. 9.46 of [2]
, il = 1.
Throat Throat Sonic
Numerical results oc Shock
inlet outlet point
Air velocity (m/s) 247 275 382 482
Water velocity (m/s) 124 160 224 271
Pressure (kPa) 683 600 358 192
Air temperature (K) 389 387 373 359
Water temperature (K) 391 390 379 367
Mass flux (kg/m2s) 2 891 2 891 2 076 1 563
Mach number 0.65 0.68 1.00 1.28
Prediction
Critical air velocity (m/s) 283 274 274 275
Dis number () 0.87 1.00 1.40 1.75
Critical mass flux (kg/m's) 3416 2 933 1740 940
The lower part presents the critical air velocity at the same
point as predicted by Eq.(19), and the associated mass flux.
The Dis number is the ratio of the numerical air velocity to
the critical air velocity predicted by Eq.(19). It can be
observed that the flow is critical at the end of the throat
section (Dis = 1) while the flow is still subsonic (Mach <
1). Between the end of the throat and the point where the
sonic velocity is reached (Mach = 1), the flow is
supercritical (Dis > 1) while being subsonic. The predicted
mass flux using Eq.(19) are equal to those obtained by the
numerical model reported by [2].
The critical conditions presented in section 2.3 suggests an
explanation for the pseudochoking phenomena reported by
[2] as Eq.(19) being the proper criteria to determine the
critical conditions. Results in Table 1 also indicate that
critical conditions may not be linked to propagation
phenomena.
4.2 Comparison of experimental critical conditions
This part provides a comparison between experimental
mass fluxes and theoretical predictions based on Eq.(19).
These results have been presented in [4] and [5]. A more
detailed discussion and analysis can be found in [4]. For
these data:
f/M =1
f 2 ap2U2 (20)
and the criteria given by Eq.(19) can be simplified to:
u2o, 1
U,, 1 Cl
A 2) (21)
Pl (1 +fM22) (1+ fM 2s22) Zl (21)
1+ ) ( l+fM2 Cp2
p1 C1 22
where s2 is the slip ratio defined by:
In this section, it is assumed that phases are isothermal.
The critical theoretical mass flux Gth is provided by:
Gth = (1+M2 )al cl (23)
Data selected have enough information available to
evaluate the flow variables at the critical section to avoid
computing the flow evolution up to the critical section. The
error is defined as:
Gth Gexp
e=100( (24)
Gth
where Gxp it the experimental mass flux.
The results of four data sets of critical flow in nozzles are
presented in Figure 2 and summarized in Table 2. The
results are in very good agreement with experimental data
over a wide range of quality and volumetric fractions. The
mass flux range covers about two orders of magnitude. The
conclusion is that the velocity difference (slip ratio) and
quality are correctly taken into account in Eq.(21).
100
Vogrin
o Smith
S Muir
A Thang
25
0
25
50
0. . 0 o O [
S 0~. e "
0,0 0,2 0,4 0,6 0,8 1,0
al
Figure 2: Critical mass flux prediction error as a function
of air volumetric fraction in nozzles.
Table 2: Description of the critical airwater nozzle flow
data sets.
Nb. Critical Exp. Mean Std.
Data Mass flow Critical
Of pressure error dev.
set data (kPa) ratio mass flux %) (%)
data (kPa) (\l (%) (%)
Vogrin 21 135392 10178 6.715.5 1.8 11.3
Smith 16 103214 0.5053 0.66.6 1.0 8.3
Muir 8 110 10807700 12.031.5 7.9 4.7
Thang 6 78204 4372800 9.616.4 0.9 6.8
All
data 49 78392 0.507700 0.631.5 0.7 9.6
set
5. Summary and future perspective
Multiphase systems present significant scientific
challenges and the choking phenomenon has to be taken
into account properly to achieve an accurate model of two
phase ejectors. This paper presents a steady
twodimensional model including new choking criterion
based on a new critical fluxes analysis in conservative
Paper No
Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
systems. They provide the information needed to include
the choking phenomenon into the numerical scheme. These
new criterions are compared with numerical data obtained
by the ASTAR project [2] and with experimental data for
validation. It is shown that these new criterions provide a
satisfactory explanation of the prechoking phenomena
reported in [2]. It also shows that steady state models are
achievable for critical multiphase flow. Such an approach
is in some cases simpler that using transient models which
involve much more complex propagation phenomenon.
Acknowledgements
The work described here was supported by NSERC,
FQRNT and NSERC Chair in Industrial Energy Efficiency
(chairholder: N. Galanis) at Universit6 de Sherbrooke with
the support of HydroQuebec (Energy Technology
Laboratory, LTE), Rio Tinto Alcan and the CANMET
Energy Technology Center (CETCVarennes, Natural
Resources Canada).
References
[1] Dostie M., Mercadier Y: Critical fluxes in conservative
systems: application to multiphase flow. To be published in
2009.
[2] Stadtke H.: Gasdynamic aspects of twophase flow:
hyperbolicity, wave propagation phenomena, and related
numerical methods, WileyVCH, Weinheim, 2006.
[3] Stadtke et al.: Advanced threedimensional twophase
flow simulation tool for application to reactor safety
(ASTAR), Nuclear and Engineering Design 235(2005),
379400.
[4] Dostie et al.: Choking criterion for onedimensional
critical multiphase flows in ducts and open channels. CNS
Conf. Procs. 13t Reactor Simulation Symposium, Chalk
river, April 2728, 1987, CRNL4139, p.5068.
[5] Dostie M: Contribution a 1'6tude des debits critiques
multiphases et a la modelisation des ejecteurs, PhD.
Thesis, University of Sherbrooke, 1988.
