Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010




Transient solution for flow of evaporating fluid in parallel pipes using the drift flux
model



Yehuda Taitel and Dvora Barnea

Tel-Aviv University, School of Mechanical Engineering
Tel-Aviv 69978, ISRAEL
taitel@eng.tau.ac.il and dbarnea@eng.tau.ac.il


Keywords: Two-Phase flow, Evaporation, Parallel pipes



Abstract

Maldistribution may occur in evaporating liquid flowing in parallel pipes with a common inlet and outlet manifolds. This is
an undesired phenomenon that takes place even for the case of equal heating of all pipes. This mal flow rate distribution is
caused by the awkward behavior of the pressure drop vs. the flow rate in evaporating lines. In single phase flow the pressure
drop is a monotonic increasing function of the flow rate. Once liquid is evaporating there is a range where the pressure drop
decreases with increasing flow rate.
Simplified models for the prediction of the flow rate distribution among the parallel heated pipes at steady state conditions
were presented by Minzer et al. (2" 1'4, 2006).
In this work we present a new transient model for the prediction of the system behavior at unsteady state conditions. The
model is based on the drift flux formulation and it is solved numerically. The pipes are subdividing into numerical cells and
heat and momentum balances are applied to each cell.
The transient model is used to analyze the stability of the steady state solutions to finite disturbances and to check the system
response to transient changes of various operational variables such as changes in the input flow rate and changes in the
heating power.
The result of this study may be useful to analyze the operation of parabolic trough solar power plants with direct steam
generation (DSG). At the present such current technology uses oil as the heated fluid that produces steam via secondary heat
exchangers.


Introduction

Flow in parallel pipes with evaporation is a common
occurrence in heat exchangers. Recently it is also
associated with solar power plants based on the trough
technology where solar radiation is focused on a pipe
heating the fluid within the pipe. At the present synthetic
oil is used as the heated fluid and steam is generated in a
secondary heat exchanger. Since the purpose is to generate
steam it would be more efficient to evaporate water
directly in the heated pipes. In this case, however
maldistribution of the flow rate within the parallel pipes
may occur.
Flow in real large-scale experiments of a single pipe with
Direct Steam Generation (DSG) operation is described by
Zarza et al. (2002), Eck and Steinmann (2002), Price et al.
(2002), and Eck et al. (2003). These papers bring to our
attention the up-to-date status of the Direct Solar Steam
(DISS) test facility at the Plataforma Solar de Almeria
(PSA). Apparently, the flow in parallel pipes is still not
well analyzed. The operation of parallel rows under steady
state and transient conditions, remain one of the main open
questions concerning Direct Steam Generation (Zarza et al.
2002, Eck et al., 2003).


Flow in microchannels is also an important area primarily
related to cooling of electronics devices (Thome et al.,
2004, Chen et al., 2002). State of the art review can be
found in Thome (2006) and Yarin et al. (2009).
Steady state solutions for the case of flow in 2 parallel
pipes was performed by Natan et al (2003) for which the
pressure drop was calculated using flow pattern analysis as
well as drift flux analysis. Minzer et al. (21 4) extended
this work and added stability analysis to demarcate among
stable and unstable steady states. Later Minzer et al. (2006)
also added a simple model to account for the transient
analysis. In their work the pipe is subdivided into 3 parts,
the subcooled liquid section near the inlet, evaporation
section and superheated section near the exit. Average
pressure drop and temperature are calculated for each
section.
In this work each pipe is subdivided into n differential
sections and the drift flux analysis is applied to each
section. The analysis yields solutions of the concentration,
C, the enthalpy, H, the pressure, P, and the flow rate G as a
function of time and position along each of the parallel
pipes. Some simplified assumptions are applied in order to
make the solution quick, robust and simple, yet without
loosing the main physical phenomenon.





7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010


a
A
b
D
c
C
Co
f
g
G
H
L
n
P
Q
Qmax
t
T
U
VGm
VGJ
W
x
X


Momentum equation of the mixture:


a 8a G2 f, C
(AG)+- A- +- A Vp
at 8x p 8x i-C


4A G\G
D 2p


coefficient (eq. 13)
cross section
coefficient (eq. 14)
pipe diameter
coefficient (eq. 15)
mass concentration
distribution parameter
friction factor
gravitational constant
mass flow rate per unit cross sectional area
enthalpy
pipe length
number of parallel pipes
pressure
heat absorbed (per unit length)
heat generation (per unit length)
time
temperature
proportionality factor for heat transfer
diffusion velocity
drift velocity
mass flow rate
axial coordinate
quality


pgA sin A
t8x


H,)= (4)
pA


Where


PGPL
p=LC PG -C)


C is the mass concentration of the gas.
FG is the gas phase generated per unit cross sectional area
and pipe length owing to evaporation:


Greek letters
a void fraction
p angle of inclination (positive for upflow)
FG gas generation per unit area and length
p density

Subsripts
G gas
in at the manifold inlet
L liquid
out at outlet
S at the position of the power generation
oo at the surrounding


Analysis

The transient behavior of the system is analyzed using the
drift flux model which is composed of the 4 balance
equations.

Continuity of the two phase mixture:


a a
(Ap) + -(AG)
at a


SA(H QH
F A(H HL)


The diffusion velocity, Vcm, is given in terms of the drift
velocity, VGJ:

VG= PG + C(pL-PG) V (C -1)G (7)
PG + C(p p) CO) (

It is assumed that the gas and fluid properties PG, PL, HG
and HL are constants in the evaporation range. Inserting eq.
(5) and rearranging yields the following 4 equations for C,
G P and H as a function of x and t.


ac [cpL+p -c)] 2 OG
t PPG(PL-P) =0
dt pLpG (pL Gp)d


aC G C avGm CVGm 9p 7C
ct p a VmC Gm pC r (9)

C 9VG 9G FG
OG )a p


Continuity of the gas phase

-C G C 1 8 (C V,,:
S-+--+-- CPAVm)
at p ax pA tx


Paper No


Nomenclature


and the mixture energy equation::

8H G H 1 8
-+--+--[AVGPC(HG
dt p Ox pA Ox





7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010


V2 1 2VGmpC OVGm
(I -C)2 1 C ac a

+ VcmC p2 C
Gm 2


2G 2VGmpC avGm G 8P
-+ -pC- -G- -+- (10)
p C aG ax 9x
4 G\G\
= -f pg sin P
D 2p


w
C(H H,)- VG + (HG H


p 8C
OVm OG OH GH OH
+C(HG H, ) +-+--
8G 8x 8t p 8x pA

Inserting eq. (8) into eq. (9) yields,

S+ C m VGm + CVGm n p + PC
p 8C p 8C 8x

S[Cp, +pG(l- C)2 VG. cG FG
PLPG(PL -P) aG x p


2b V -1 2VGmpC aVGm
b = yV p +
I2 C 1C

Gm V C_ ) ac
c, IC P2 OC


c=C(HG -H) VG +VGm(HG

+Gm (HG H )C p
P aC


HL)


Equations (13, 14 and 15) are the balance equations for
any single pipe in the set of parallel pipes.
In order to handle multiple parallel pipes we integrate the
(11) momentum equation (14) along the pipe..


lb BC d+G L +
i -&+L+P^


]Pn)


I D -2p pgsin dx
0 1 D 2p


For n parallel pipes with common inlet and outlet
manifolds the pressure drop in all n pipes is the same.
(12) Summing eq (16) for n parallel pipes yields

P,, -Po't


We further assume that the mass flow rate, G, is not a
function of x (quasi equilibrium assumption), thus
aG / ax = 0. In this case the equations are:

Eq. (12) yields,


ac
a-
a


Eq. (10) yields,


bC a G
b + +
ax at


D4 G2G
D 2p


-pgsin p (14)


And eq. (11) yields


ac 8H G H Q
Cx -+-+-
Dx at p ax pA


Where


1L cc 9C ?L 4 GG\ in P 1
Jfb +[f D 2p +pg sin 0x
Jo & L 4 GGl (17)i



+ ldG
n dt

Where G1m is the total mass flow rate per unit area entering
the inlet manifold. Note that for Gm=constant dGm/dt=0.
Equation (17) is used to calculate the pressure at the inlet
of the parallel pipes. Once Pm is known aG / at in each
pipe is calculated using eq. (16).
Equations (13, 14, 15) for each single pipe together with
equation (17) for n parallel pipes are used to predict C, P
and H as a function of x and t and G as a function of t.


Solution

The finite difference form of eqs. (13) and (15) is:


G 8VGm CVG 8a
a =-+C + VGm -
p aC p aC


C -C
lk k+1 1,k+l
Ax


Paper No


aG
+-+
at





Paper No


C,k+l Cl ,kl H,+l Hk
k Ax At

GK Hk+1-Hi-1,k+1 Q
pl. Ax pA ,,


Where i is the spatial index and k is the time increment.
Using (18) we calculate Cl,k+1 in each differential section of
every pipe, then using (19) we calculate H,k+1. Next we
calculate the inlet pressure, Pin, using (17) in which the
integration is performed numerically. The flow rate in each
pipe is then calculated at the next time step using eq. (16)
yielding the next Gi,k+l

L 4 G\G\

G,k+l = G,k + (20)

t -P ) J aC dx


Note that the flow rate along each pipe is constant at any
time, thus Gi,k+l =Gk+l.
The fluid inlet flow rate and properties and the outlet
pressure are specified as boundary conditions.
In our experimental system the heat flux Qmax is generated
in the electrical resistance located between two layers of
insulating materials (see figure 1). Only part of the
generated heat is absorbed by the fluid due to heat losses.
The heat absorbed depends on the temperature of the fluid
in each section:

Qmx =U(T T)+ U (T T,) (21)

where the first term on the r.h.s. is the heat absorbed by the
fluid, Q, and the second term is the heat loss to the
surroundings. T and Ts are the fluid temperature and the
resistance temperature and U and Uo are the equivalent
heat transfer coefficients to the fluid and to the
surroundings respectively.
Eq. (21) is used to calculate the resistance temperature Ts
and the heat absorbed by the fluid. Note that as the fluid
temperature increases more of the generated heat is lost to
the surroundings and less is absorbed by the fluid.
Eventually the fluid reaches a maximum temperature.

Distribution parameter and drift velocity.
Since during evaporation the flow pattern is mostly in the
form of annular flow with high void fraction the
correlations for annular flow are used for the distribution
parameter and the diffusion velocity (Ishii and Hibiki,
2006).

1-c
Co =1+ (22)
a+4 -PG PL


G -D(p a)
VGJ = (23)
0.015p,


7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

Since the void fraction for annular flow is close to unity Co
can be taken as 1 and the drift velocity is close to zero.
The above equations are applicable for the range of the
evaporation. Obviously the flow also contains sections
were the flow is subcooled liquid at the entrance and it
may be superheated vapor near the exit. The equations for
the single phase, though used in the calculations, are
omitted in this paper.

Results:
The solutions are performed for the operational conditions
of our experimental system presented in Minzer et al.
(2006). The evaporating fluid is water and the conditions
are as follows:


Pipe length
Pipe diameter
Outlet pressure
Inlet temperature
Surrounding pressure
Power
Heat coefficient
Heat coefficient


L
D
P out
Tm
To
Qmax
U
U


6m
5 mm
0.1 MPa
300C
250C
1200 W/m
8 W/mC
2 W/mC


Steady state solution, single pipe
Figure 2, 3 and 4 show the variation of C, X, a and T along
a single pipe at steady state. C is the local mass
concentration of the gas, X is the ratio of the local mass
gas flow rate (WG/W), a is the local void fraction and T is
the temperature. Figure 2 show the solution for the
homogeneous model (VGJ=0.0). Figure 3 shows the case of
a very low flow rate and figure 4 the case where the drift
velocity is 0.025 m/s.
As can be seen the concentration, C, and the quality, X,
increases linearly with distance at the evaporation zone. C
and X are the same in figures 2 and 3 which corresponds to
the homogeneous model. The temperature increases
linearly with distance in the subcooled liquid, it is constant
at the evaporation region and it increases dramatically at
the superheated gas region. Note that the rate of increase of
the temperature in the gas zone decreases as the
temperature increases. This is due to the losses to the
surroundings that increases as the water temperature
increases. In the case of a longer pipe, or a low flow rate
the subcooled and the two phase regions are short. The
superheated section is long and the gas temperature may
reach a maximum value. At this temperature no heat is
absorbed by the fluid and the generated heat is lost to the
surroundings. In this case T=Ts (see eq. 21). Figure 3
represents the case of a very low flow rate and that a
constant maximum temperature is reached along the pipe
at the exit region.
The void fraction, c, increases very quickly to a value
close to 1. Thus most of the two phase flow region is in
annular flow. Note that the results for the quality, void
fraction and temperature are not affected by the drift
velocity. The concentration, however is affected by the
drift and it is lower than the value of the quality.
In figure 5 the steady state inlet pressure (the black curve)
as a function of the flow rate in a single pipe is plotted.
Typical to the variation of Pin vs. flow rate W for
evaporating system there is a range of flow rate where the
pressure decreases with the increase of the flow rate.






Paper No


Steady state solution, two parallel pipes
Since the multiple pipes are connected by common inlet
and outlet manifolds the solution requires the same
pressure drop in all parallel pipes. Scanning now all the
possible pressure drops along the parallel pipes yields all
possible flow rate distribution among the pipes for a given
total inlet flow rate and given heat fluxes that satisfy the
requirement of equal pressure drop in all pipes. Figures 6
and 7 show the results for 2 pipes in parallel with identical
heating. Figure 6 present the flow rate ratio in each pipe as
a function of the inlet flow rate. As can be seen multiple
solutions are possible. For high or low flow rates a single
even flow rate distribution is the only solution. Yet for
intermediate inlet flow rate 3 and 5 solutions are possible.
Linear stability analysis (see Minzer et al., 2006) allows to
identify the unstable steady states solutions shown in the
figure by the red lines while the stable solutions are the
blue lines. It is interesting to see that there is a wide range
of inlet flow rate that equal distribution of the flow is
unstable while non even distribution is stable.
Figure 7 shows the steady state solution on a Pin vs. Wm
chart. Again the blue lines are the stable solution and the
red ones are the unstable, meaning they will not exist in
practice.

Transient simulation
Transient simulations are used to investigate the system
response to finite disturbances in the state variables such as
the flow distribution ratio as well as to simulate transient
responses to time dependent input in the operational
conditions such as the total flow rate or the heating power.
Three examples of transient processes are calculated
(1) Starting at an hypothetical unstable solution (point Ai
in figures 6 and 7) and follow the response to a finite
disturbance that ends at points B and C in figures 6,7.
(2) Starting at a high inlet flow rate (point A2), decrease
the inlet flow rate and end up again at B and C.
(3) Starting at a low flow rate A3 and increase the inlet
flow rate again to end up at B and C.
Figure 8 shows the transient response to a finite
disturbance starting from a hypothetical even, but unstable,
flow distribution (point A1). It is obvious that the solution
will depart from this unstable steady state towards a stable
steady state for the same inlet flow rate, Win. yielding an
uneven flow rate distribution (points B & C). During this
process the pressure decreases The trajectories of this
process is also shown on figures 5.
Figures 9 to 14 shows the response of the system to a
change in the operational conditions.
Figure 9 starts at a stable steady state condition where
Wm=0.08 kg/s and the flow is equally distributed between
the two pipes as subcooled liquid (point A2 in figure 6).
The total inlet flow rate is deceased gradually to Win=0.24
kg/s (the black line in figure 9).
Initially the path is identical for the two pipes (the green
and the red curves). Then the flow rate in each pipe reaches
the uneven steady and stable flow rate distribution (points
B & C). The path of this process is shown also in figure 10
(A2 to B and C). The transient behavior of the inlet
pressure is shown by the blue line on figure 9. It decreases
to a minimum and then overshoots before reaching the new
steady state. Note (figure 10) that the inlet pressure does


7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

not follow the steady state pressure line and it is lower than
the steady state pressure. This is due to momentum losses
during the decrease of the flow rate and lesser inlet
pressure is needed to drive the flow rate in each pipe
during the transient change.
Figure 11 and 12 show the same change as in figures 9 and
10 but the change of the inlet flow rate is slower. In this
case the transient trajectory is closer to the steady state
pressure line.
Figures 13 and 14 shows the transient response of the
system for increasing the inlet flow rate from Win=0.006
kg/s to Wm=0.024 (the black curve on figure 13). As a
result the solution changes with time from the stable steady
state (point A3) where the flow rate is equally distributed to
the points B and C (as before). Note that during this
transient the pressure increases substantially before
deceasing towards the final steady state.

Summary and Conclusions
The drift flux model is used to calculate the transient
variation of the flow rate distribution and the inlet pressure
in evaporating fluid flowing in parallel pipes.
The results are demonstrated for conditions similar to our
laboratory system that use electrical power to heat water in
parallel pipes, 6 m long and 5 mm diameter.
The result show how unstable solution of even flow rate
spontaneously changes to a stable unequal flow
distribution between 2 pipes.
Slow Change of the inlet flow rate with time will follow
the steady state curve. A faster change in the inlet flow rate
result in a path different than the steady state path and may
yield higher and lower peaks in the inlet pressure.


Acknowledgements

The financial support from the Israeli Science Formation is
greatly appreciated.


References

Chen, W.L., Twu, M.C., Pan, C., Gas-liquid two-phase
flow in micro-channels. Int. J. Multiphase Flow Vol. 28,
1235-1247 (2002)

Eck, M. Steinmann, W.D., Direct steam generation in
parabolic troughs: first results ofDISS project. Transactions
of the ASME, J. of Solar Energy Eng., Vol. 124, 134-139.
(2002)

Eck, M., Zarza, E., Eickhoff, M., Rheinlknder, J.,
Valenzuela, L., Applied research concerning the direct
steam generation in parabolic troughs. Solar Energy Vol.
74, 341-351 (2003)

Ishii, M. and Hibiki, T, Thermo fluid dynamics of
two-phase flow. Springer (2006)

Minzer, U., Barnea, D., Taitel, Y., Evaporation in parallel
pipes splitting characteristics. Int. J. Multiphase Flow,
Vol. 30, 763-777 (2~ri14).






Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010

Minzer, U., Barnea, D., Taitel, Y., Flow Rate Distribution in
Evaporating Parallel Pipes Modeling and Experimental.
Chem. Eng. Science, Vol. 61, 7249-7259 (2006). 1.2 350
W=0.0016 kg/s /C 3n
Natan, S., Barea, D., Taitel, Y., Direct steam generation in 1 ,- W0.0016 kg s 3
parallel pipes. Int. J. Multiphase Flow, Vol. 29, 1669-1683 / 250
(2003). 0.8- a

Price, H., Liipfert, E., Kearney, D., Zarza, E., Cohen, G., X 0.6 -15
Gee, R., Mahoney, R., Advances in parabolic trough solar 04 -1
power technology. Transactions of the ASME, J. of Solar 1.4
Energy Eng. Vol. 124, 109-125 (2002). 0.2 50

Thome J.R., Dupont V., Jacobi A.M., Heat transfer model o -
for evaporation in microchannels, Part 1. Presentation of the o 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
model. Int. J. Heat Mass Transfer, Vol. 47, 3375-3385 x/L
("'14). Figure 2: Variation of C, X, a and T along the pipe,
negligible drift velocity
Thome J.R., State of the art overview of boiling and two
phase flow in microchannels. Heat Transfer Eng. Vol. 27
(9) 4-19 (2006). 12 700
T
Yarin, L.P, Mosyak, A, Hetsroni, G, Heat transfer and 1 600
boiling in micro channels. Springer (2009). Ca x5
500
08
Zarza, E., Valenzuela, L., Leon, J., Weyers, H. D., Eickhoff, a W=0.0001 kg/s 400
M., Eck, M.., Hennecke, K.,. The DISS project: direct steam x 06
generation in parabolic trough systems. Operation and u 300
maintenance experience and update on project status. 04 200
Transactions of the ASME, J. of Solar Energy Eng. Vol.
124, 126-133 (2002). 02 10oo

0o II-------------0o
0 01 02 03 04 05 06 07 08 09 1
x/L
Figure 3: Variation of C, X, a and T along the pipe for low
flow rate

1.2 350

max 1 W-0.0016 kg/s 300
VGJ=0.025 m/s 250
0.8

SX C 200
X 0.6
T To X 150s
0.4
U Ts U.,
0.2 50


0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x/L
Figure 1: Geometry of the heating system Figure 4: Variation of C, X, a and T along the pipe, with
drift velocity







Paper No



030
028
026
024
m 022
0 20
01 018


7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010


000 001 002 003 004 005
W (kg/s)

Figure 5: Steady state solution for the inlet pressure and
transient path of the flow rate in a system of 2 parallel pipes







09 B
STABLE STEADY
0 8 STATE
07
= 06
A3 A, A2
05
04
03
0 2 UNSTABLE STEADY
STATE
01 C
0
0 001 002 003 004 005 006 007 008 009 01

Win (kg/s)

Figure 6: Flow rate ratio for steady state solution and linear
stability for 2 parallel pipes.


03
028
026
024
022


91 018 A2
016
014


012
0 001 002 003 004 005 006 007 008 009 01

Wi, (kg/s)

Figure 7: Inlet pressure for steady state solution and linear
stability for 2 parallel pipes.


03
028
026
024
" 022
-
02
.E 018
016
014
012
01


0025


0 02


0015 '


001


0 005


0


0 01 02 03

t (sec)

Figure 8: Transient spontaneous response of inlet pressure
and flow rate in 2 pipes.





03 009
0 28 0 08
026 007
0 24 Win0.08 0.024 kg/s 006
S022
02 A2 05
004
0 18
Olo \002
016 B 0 03
014 W1 002
012 W2 C 001
01 0
0 05 1 15 2

t (sec)

Figure 9: Transient solution of the flow rate in a system of 2
parallel pipes during decrease of inlet flow rate in 1 sec


000 001 002 003
W (kg/s)


004 005


Figure 10: Steady state solution for the inlet pressure and
transient path. Decrease of inlet flow rate in 1 sec.






Paper No



03 -
028
026
024
S022
02
0 018
016
014
012
01


Figure 11: Transient solution of the flow rate in a system of
2 parallel pipes during decrease of inlet flow rate in 2 sec


000 001 002 003
W (kg/s)


004 005


Figure 12: Steady state solution for the inlet pressure and
transient path. Decrease of inlet flow rate in 2 sec.


0 09
008

0 07
Win=0.08 0.024 kg/s 006

005
004

W1 B 003
0 02

001
W2 C

05 1 15 2 25 3
t (sec)


000 001 002 003 004 005
W (kg/s)

Figure 14: Steady state solution for the inlet pressure and
transient path. Increase of inlet flow rate in 1 sec


7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30-June 4, 2010


03 003
028 Wm=0.006 0.024 kg/s
026 -B 0 025
0 26
024 W 002



018 /
016 001
014XA3 0005



0 05 1 15 2
t (sec)

Figure 13: Transient solution of the flow rate in a system of 2
parallel pipes during increase of inlet flow rate in 1 sec