7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Efficient Numerical Solution of OneDimensional Governing Equations for
Evaporating Flow in a Tube
C.T. Veje, S. Madsen, and M. Willatzen
Mads Clausen Institute, University of Southemn Denmark
Alsion 2, DK6400 Sldnderborg, Denmark
veje~mci.sdu.dk, sma~mci.sdu.dk, willatzen~mci.sdu.dk
Keywords: Twophase flow, evaporation, heattransfer, numerical methods
Abstract
We present results of the computational stability and efficiency in simulating twophase flow and heat transfer in
a horizontal tube. Two structurally different models of the process are presented and compared for two dynamic
computational scenarios; dynamic response for a constant and changing number of flowzones, respectively. The first
model is a simple generic lumped model of the moving boundary type and the second is a fully distributed model
in one dimension. The moving boundary model is solved as a set of conventional ordinary differential equations
whereas the distributed model is solved using a highlyefficient numerical scheme following Kurganov & Tadmor
(2000). The governing equations are formulated in terms of mass flow, enthalpy, pressure, and tube wall temperature
based on the continuity, momentum and energy equations for the fluid. A wall heatexchange equation accounts for
convective heat exchange between the fluid and the ambient air. The comparison demonstrates a \ignilk~.lill difference
in the results obtained using the two models, especially in the case where the number of fluid zones is changing. This
emphasizes the importance of using either more advanced moving boundary models or fully distributed models in the
simulation of systems comprising twophase flows. In contrast to conventional belief we find that the finite difference
model is indeed both stable to discontinuities in boundary conditions; changing number of fluid zones; high frequency
perturbations and also at the same time numerically fast enough for realtime simulation and control.
Introduction
Modeling and simulation of twophase phenomena is
important in the development and optimization of a wide
range of industrial applications and products. Among
others these include essential components of, e.g., the
HVAC industry, power plants and the growing efforts
within energy storage technologies. As requirements for
model accuracy of such systems increase there is a grow
ing if not compelling need for detail in the modeling
of the basic phenomena. Such increased detail usually
comes at the expense of computation time. The central
question is then if a reasonable compromise between ac
curacy (detail) and computational speed can be met. The
present work will address this.
A main component in the aforementioned systems
is the evaporative heat exchanger subject to twophase
flow and heat transfer. Traditionally, the lumped com
ponent models for evaporators are solved using NTUe
 methods. The need to investigate for example evap
orator control strategies demands realistic modeling of
parameters such as the superheat which has given rise
to a widespread use of movingboundary models, see
e.g., He et. al. (1998); Willatzen et. al. (1998); Zhang
& Zhang (2006); McKinley & Alleyne (2008). These
models show an advantage in demonstrating dynami
cal behavior of the component in relation to capacities,
temperature, and pressure levels although formulated as
lumped models. This has enabled evaluation of system
performance and analysis with respect to varying condi
tions and different control strategies. The drawback is
loss of detail in the modeling and the handling of fluid
zone switching which may result in numerical instability
and increased computation time.
In a general effort to investigate the dynamic behav
ior of evaporative flow systems several works have been
carried out on fully distributed models solving the gov
emning equations, see e.g., Jia et. al. (1995, 1999); Val
ladares et. al. (2004); Madsen et. al. (2009). In this pa
per, we present models and comparative simulation re
sults for two implementations of a horizontal tube with
an evaporating fluid; a conventional moving boundary
I > z
Mixed phase Super heat
gas phase
I II I
Figure 1: Schematic representation of moving boundary model.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
mif t) h, (t) I I
and a distributed model. Two different cases of dynamic
response are investigated.
In the following we will give a brief review of the
model equations. First, a set of ordinary differential
equations is developed based on a conventional moving
boundary formalism following the works of He et. al.
(1998) and secondly a full set of governing equations
is developed following Madsen et. al. (2009). In both
cases we make the following assumptions:
 constant tube cross section,
no gravity effects,
onedimensional flow,
 twophase inlet condition,
no axial heat conduction in the fluid,
no viscous heating or viscous pressure loss
 fixed temperature heat source in the airside,
 constant airside heat transfer coefficient and area.
The working medium is isobutane (R600a) and thermo
dynamic quantities such as T(h, P) and p(h, P), satura
tion values and fluid properties are evaluated using Engi
neering Equation Solver, EES (Klein (2009)) or Refrig
erant Equations v. 6.1 (Skovrup (2001)) by generating a
lookup table for fast access to these relations.
1 Moving Boundary Model
The system under consideration is shown in figure 1. It
constitutes a straight horizontal tube of length L and in
ternal diameter Di with fluid flowing from left to right
while it exchange heat with the tube wall.
A schematic view of the simple moving boundary
model is also indicated in figure 1. We treat the tube
wall as a single lumped block. On the fluid side the
model constitutes a twophase part L, P, and a single
phase gas part; the superheat region L~SH. Clay 081 =
L, P + L~SH and the position of the dryout point L, P is
determined by the balance between the inlet mass flow,
outlet volume flow, the pressure and heat transfer rates
through the tube wall. The air to wall heat flow is;
lw=fini7Do1a(TA TW)7,
where to, Do, ffin, TA, TW are the heattransfer coef
ficient between evaporator wall and ambient, outer tube
diameter, fin factor, ambient temperature, and wall tem
perature, respectively. The outer tube area is reduced by
a factor y LTP/L. (as a correction) due to reduced
heat transfer over the superheat region given by the liq
uid volume V, = D/ L,2 ~P(1 a), total volume V and
void fraction a,
V( a (2)
The wall energy equation is given by:
dTw
T~ww d =(QW Qf), (3)
where M~wCw is the thermal capacity of the wall and
Qf is the heat flow into the tube fluid
Qf = Gevap + sh,
where
Qevap = i7DiLa4(TW Te)7. (5)
Te is the evaporation temperature and as is the internal
heat transfer coefficient which is assumed constant and
estimated based on the correlation used in the finite dif
ference model described later. Qsa is evaluated as
Qsh = ro2(ho hg),
where h, is the dewpoint enthalpy and ho, ro2 is the tube
outlet enthalpy and mass flow respectively. The outlet
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
These equations are recast in terms of new variables
(pressure, mass flow, and enthalpy) P, m = pAw, and h
8~p BP 8~p 8& 1 84m
+h +p (13)
B3P 8t 8&h 8 A 8z
A 8 8 8z (14)
88A h BP dz BP
8 pA 8z\ 81 pA 8zm P
4 zr (Tw T), (15)
Di
temperature is evaluated using a functional form for the
superheat similar to He et. al. (1998). Finally the change
in liquid mass is given by
dMI
= (1 x)As Reap,(7)
where x4 and ris is the inlet quality and mass flow, re
spectively, and mevap evap/A/* iS the evaporation
flow with al being the evaporation enthalpy at the rel
evant pressure. The outlet mass flow is simply assumed
to follow the inlet mass flow and the change in liquid
content and thus:
keeping in mind that p p(P, h), T T(P, h).
() The derivative of p with respect to P equals the inverse
square root of the speed of sound and the latter quantity,
tion
in the case of twophase flow, is derived using the model
of Nguyen et. al. (1981). ish evaluated by numerical
heat 8
differentiation. To close the system the energy equation
for the wall is given by
mo = rji
(9)
An important variable in the model is the determine
of the pressure level which results as a combinatioI
the in and outlet flow conditions and the overall 1
transfer. The pressure is calculated as
dP BP 8 p
di 8p 08 t '
8 Tw
(CwpwAw)8
aiixoi(T Tw)
+aox7Doffin (TA TW) + XA~W 2
i8z2
where dpis evaluated through the continuity equation
for the total gas volume of the tube V~(LSH TiP) L
Thus, the pressure dynamics is assumed to be deter
mined only by the balance of inlet and outlet mass flow
and not by variations in temperature since temperatures
vary slowly compared to mass flows.
The moving boundary model is solved with appro
priate boundary conditions using a conventional ODE
solver. For an evaporator tube of length L we have an
inlet mass flow ris at z 0 and an outlet volume flow
Vi mo/po at z L. Pipe wall ends (z 0 L) are
insulated and the fluid enthalpy is specified at the inlet
z = 0
2 Distributed Model
For the distributed model we assume that variations in
physical properties basically take place along one coor
dinate direction (the length direction of the evaporator).
Then, following Madsen et. al. (2009) we have
where pw, Aw, and A are the wall density, wall cross
sectional area, and wall heat conductivity, respectively.
As in the case of the moving boundary model we use a
constant coefficient of heat transfer for the airside. For
the fluid side we use the correlation due to Kandlikar
(1991) for the twophase flow up to a vapor quality of 0.7
and the DittusBoelter correlation for the single phase
flow regime. A simple linear interpolation is used be
tween the two correlations. A correlation for the viscous
pressure loss could have been included but is assumed
to be zero for the sake of comparison with the moving
boundary model.
The computational domain as indicated in figure 1 is
discretized in N fixed finite elements as opposed to the
two variable elements for the moving boundary case.
The above formulation is essentially a homogeneous
flow model since we have formulated the equationset
in terms of overall density and enthalpy neglecting the
separation in gas and liquid phases as is more conven
tional in modeling twophase flows Wallis (1969).
The boundary and initial conditions are similar to
those for the moving boundary model and the relations
are given by:
8p 8 (pw)
+ 8z= 0,
8 (pW ) dP
8z 8z '
agi(Tw T),
Di
8 (pw)
+ w
8~z
8&
(
p 8
m(0, t) = rjs ,
h(0, t) = he ,
rj(L,t) = p(L,t)V,
aTw
d0
where w is the onedimensional velocity along the evap
orator length coordinate.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
A conventional central finite difference method ap
plied to Eqs. (13)(16) typically shows spurious oscil
lations which quickly destroy the numerical solution. In
the fortunate case of numerical convergence, which is
highly unlikely for a dynamic simulation, the compu
tation time is unreasonably long and dynamic control
of a real time simulation is not feasible. As demon
strated in Madsen et. al. (2009), the KurganovTadmor
(KT) highresolution scheme can avoid these difficulties
and allows possible discontinuities in the solutions while
maintaining computational speed and stability. We have
applied the KT scheme in the second order semidiscrete
form to Eqs. (13)(15) and implemented it in c++. The
heat equation (16) is discretized using straightforward
second order finite differences. Further details can be
found in Kurganov & Tadmor (2000) and Madsen et. al.
(2009).
Variable/Parameter
Value
6 103
8 103
10
(D
385
8.96 10"
4((D D")
386
2100
100
300
2.5 10"
1.75 103
1.001 10"
263
273
h;
m
m
[]
m
J/K
kg/m3
m
W/(mK)
W/(m2 K)
W/(m2 K)
K
J/kg
m3/s
Pa
K
K
J/kg
and initial data
A
as (MB model)
TA
h;
1"
P(t =
T(t =
Tw (t
h(t
= t
0, z)
0, z)
S z)
0, z)
Table 1: Parameters, boundary values,
used in the simulations.
Figure 2: Steady state temperature T(2) at the tube out
let for different discretization levels N. Inset
shows the result for the entire length of the
tube.
Figure 2 shows the steady state solution for the tem
perature at the tube outlet for different values of dis
cretization level N. Already for N 50 not much is
gained in increasing to N 100 and the relative differ
ence in temperature T,, between N 100 and N 400
is as low as 5 104 K. The inset shows the result for
the entire length of the tube. Temperature gradients are
converged to 104 K/s after 10 seconds and steady state
is obtained by solving the system for 100 s. For the re
maining part of this work we use N = 100. Compu
tation time is O(N2), but even for N 100 we only
use 1.17 s per simulated second on a 2 GHz desktop.
In other words, realtime simulation is definitely achiev
able using a distributed model!
3 Simulation Results
In the following, we present results for the moving
boundary (MB) and the finite difference (FD) models
MB model
FD model
290
~285
S280
2 46 8
Figure 3: Steadystate solution for the fluid temperature
as function of position z in the tube. Re
sults are recorded after 100 seconds simula
tion. Data for the MB model is generated us
ing the exponential form of He et. al. (1998).
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
using the basic parameters, boundary values, and ini
tial data shown in Table 1 and otherwise stated in the
text. In figure 3, we show the fluid temperature T as a
function of position z for both methods. The results are
steady state values after 100 seconds simulation time.
By roughly tuning as in the MB model we find approx
imate agreement between the two models. Notice the
acceleration pressuredrop for the FD method. The tem
perature drops slightly until dryout at z ~ 7 m where
after the temperature rises quickly in the single phase
gas zone.
To test the dynamic response and numerical stability
of the model we next investigate the model response to
steps in inlet mass flow. After simulating the system for
100 s to a steady state we increase the inlet mass flow
ris 7.5 103 kg/s by 10% at t 100 s and decreases
it again to the initial value at t 150 s. Figure 4 shows
simulation results of rjs, rn0 (top), To (middle) and Ly P
(bottom) for both models as function of time. The mod
els behave as expected and are numerically stable. Dif
ferences in initial transients for t < 10 s reflect differ
ences in initial conditions between the models. For the
response of the outlet mass flows reasonable agreement
between the models is demonstrated, however, the tem
perature dynamics is clearly different for the two mod
els. This is to a large extent due to the lumped character
of the wall energy equation in the MB model. Since
outlet temperatures are typical signals in the control al
gorithms for such systems this points to the necessity of
increasing the detail of the modeling.
The computation times are of the order of one second
for the MB method depending on the choice of integra
tor. For the FD model it is of the order 100 s. For the FD
model, we use N 100 and the timestep is determined
by the maximum velocity criterion as given by the KT
method. This roughly corresponds to a Courant num
ber based on the maximum flow and sound speed and
is therefore orders of magnitude smaller than the typical
timestep of the MB method.
When the system is perturbed by an instantaneous
change in, e.g., inlet mass flow as in figure 5, a pressure
wave is formed shuttling back and forth in the evapora
tor tube. Figure 5 shows a detailed spatiotemporal view
of rj(z, t) near t 100 s when the inlet mass flow iS
increased as in the case of figure 4. The wave is formed
at the inlet z = 0 m and moves with the flow towards
the outlet. At t ~ 100.1 s, the wave reaches the outlet at
z=10 m and reflects back into the system. The speed
of the wave is not constant over the evaporator tube due
to the changing speed of sound in passaging mainly the
twophase region. Obviously, this phenomenon cannot
be captured by the MB method. However, although
these effects may not be important on timescales rele
vant in the heat transfer process the resolution of these
_m
mo (MB)
mo (FD)
L
so loo
Time [s]
298
296
S294
S292
E 290
S288
5 286
O
284
282
38n 
MB
150 200
50 100
Time [s]
i ('
50 100
Time [s]
150 200
Figure 4: ri, T, and L, P as a function of time. Top:
outlet mass flow mo response to a step in in
let mass flow ris at time t 100 s. Middle:
Outlet fluid temperature. Bottom: Position
of dryout point L, P. Results are shown for
both MB and FD models. Differences in ini
tial transients for t < 10 s reflects differences
in initial conditions.
I I I I I lU
0 24 6 810
z [m]
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
100.4
100.3
100.2
x 103
10
m
mo (MB)
mo (FD)
115 120
100.1
100
99.9
105 110
Time [s]
 MB
FD
294
a 292
7.4 a 290
E
~ 288
S286
Figure 5: Spacetime dependence of th2(2, t) as the in
let mass flow thi is increased by a 10% at
t 100 s. A pressure wave is seen to rapidly
shuttle through the system.
20900 105 110 115 12
Time [s]
in the FD method points to the stability of the model.
For the integration of ODE's in the MB model the
main part of the computation time is spent at events of
discontinuity such as when the inlet mass flow changes
abruptly or when the number of flow zones changes. A
typical situation in superheat control of evaporative flOW
systems involves the response of the inlet mass flow tO
external perturbations. Situations may occur where there
is a repeated change of fluid zones if the inlet mass flow
is modulated. Figure 6 (top) shows this situation. Start
ing from a similar situation as in figure 4 we now in
crease the inlet mass flow and modulate it with a higher
frequency to an extent that the number of fluid zones
changes with each period.
For this highly perturbed situation the simulation re
sults reveal large differences between the models. Fig
ure 6 (top) shows the mass flow responses and the miss
ing dynamic terms in the MB formulation is now clearly
reflected in the dynamics of the model. The correspond
ing results for the temperature and dryout point in fig
ure 6 (middle and bottom plots) reflect the same short
MB"
1 5 120
105 110
Time [s]
Figure 6: rih, T, and Lr P response to a modulation of
inlet mass flow thi. Top: resulting outlet mass
flow tho. Middle: Outlet fluid temperature.
Bottom: Position of dryout point L~TP is now
pushed to the exit of the tube. Results are
shown for both MB and FD models.
103rn [kg/s]
In 8.4
0
98
96
94
~92
9
"88
8 6
84
82
10
9 9
9 8
97
9 6
9 5
9 4
9 3
9 2
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Conclusions
100.4
100.3
100.2
Through implementation of two numerical models for
the solution of the twophase flow and heat transfer in
a horizontal tube we have found \ignilk~.llll differences
in the results of the two models. Especially in the case
9.8
where the number of fluid zones changes there are sub
stantial deviations. This stresses the importance of using
either more advanced moving boundary models or fully
distributed models in the simulation of systems compris
9.6 ing twophase flows. The FD model implemented here
is an explicit formulation of the complete set of govern
ing equations including dynamic effects such as propa
gation of sound waves and acceleration pressure losses.
9.4 Despite this level of detail the FD model is found to be
numerically stable and sufficiently efficient that it facil
itates realtime simulation and feedback control of such
systems.
9.2
References
A. Kurganov, E. Tadmor, New highresolution semi
9 discrete central schemes for HantiltonJacobi equations,
Journal of Computational Physics 160, 720742 (2000).
100.1
100
99.9
z [m]
Figure 7: Spacetime dependence of t2(2, t) as the in
let mass flow the is modulated with a high
frequency and the number of fluid zones is
changing. Pressure waves are seen shuttling
back and forth through the system.
comings of the MB model. The discrete levels for the
FD model in the bottom figure is due to the discretiza
tion of the tube but could be easily avoided using inter
polation. A dryout position of 10 m corresponds to a
situation where the tube contains only two phase flow
and liquid appears at the tube exit.
An interesting test of numerical stability and effi
ciency of the FD model is shown in figure 7 where
the spatiotemporal plot of the mass flow in the tube is
shown for the case where the modulation frequency of
figure 6 is increased to a period of 0.1 s. The dynamic
pressure waves are reflected at the tube ends and crosses
while the number of fluid zones changes periodically in
a similar way as in figure 6. The computation time of
the results in figure 7 is 0.58 s which clearly facilitates
realtime simulation and control of such a system even
in the case where detailed distributed models are used.
S .Madsen, C.T. Veje and M. Willatzen, knplementa
tion Of The KurganovTadmor High Resolution Semi
Discrete Central Scheme For Numerical Solution Of The
Evaporation Process h2 Dry Expansion Evaporators, in
Proceedings of the 50'th SIMS Conference, (2009).
X.D. He, S. Liu, H. Asada, and H. Itoh, Multivariable
control of vapor compression systems, HVAC&R Re
search 4, 205230 (1998).
M. Willatzen, N.B.O.L. Pettit, and L. PlougSedrensen,
A general simulation model for evaporators and con
densers in refrigeration. part I: Moving boundary for
mulation of twophase flows with heat exchange, Int. J.
Refrig. 21, 398403 (1998).
W.J. Zhang, C.L. Zhang, A generalized moving
boundary model for transient simulation of dry
expansion evaporators under larger disturbances, Int. J.
Refrig. 29, 11191127, (2006).
T. L. McKinley, A. G. Alleyne, An advanced nonlinear
switched heat exchanger model for vapor compression
cycles using the mo vingboundary method, Int. J. Refrig.
31, 12531264, (2008).
X. Jia, C.P. Tso, P.K. Chia and P. Jolly, A distributed
model for prediction of the transient response of an
evaporator, Int. J. Refrig. 18, 336342, (1995).
103rn [kg/s]
n 1 0
I
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
X. Jia, C.P. Tso, P. Jolly and Y.W. Wong, Distributed
steady and dynamic modeling of dryexpansion evapo
rators, Int. J. Refrig. 22, 126136, (1999).
O. GarciaValladares, C.D. PerezSegarra and J. Rigola,
Numerical simulation of doublepipe condensers and
evaporators, Int. J. Refrig. 27, 173182, (2004).
S.A. Klein, Engineering equation solver (EES), Aca
demic Commercial V8.4 McGraw Hill (2009).
M. Skovrup, Thermodynamic and thermophysical prop
erties of refrigerants, 3.10, Department of Energy Engi
neering, Technical University of Denmark (2001).
D.L. Nguyen, E.R.F. Winter and M. Greiner, Sonic Ve
locity in ETwoPhase Systems, Int. J. Multiphase Flow 7,
311320, (1981).
S.G. Kandlikar, A Model for Correlating Flow Boiling
Heat Transfer in Augmented Libes and Compact Evap
orators, J. of Heat Transfer 113, 966972, (1991).
G. Wallis, Onedimensional LT'ophase Flow, McGraw
Hill (1969).
