7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Modelling Mixed Flows in Civil and Environmental Engineering:
A 1D Threephase Approach
Kerger Frangois1'2, Archambeau Pierre1, Dewals Benjamin J.1,2, Erpicum Sebastien', Pirotton
Michell
SLaboratory of Hydrology, Applied Hydrodynamics and Hydraulic Constructions, Liege University
Chemin des chevreuils, 1, B52, B4000Liege, Belgium
Email: fkerger@ulg.ac.be
2 Belgian Fund for Scientific Research FNRSF.R.S
Keywords: Preissmann slot, Driftflux model, Stratified flows, Hydraulics, Numerical scheme
Abstract
Mixed flows, characterized by a simultaneous occurrence of free surface and pressurized flows, are frequently encountered in
civil and environmental engineering. Single phase models coupled with a transport equation fail to account for
subatmospheric pressurized flows and airwater interactions, which are both likely to appear in mixed flows. This research is
thus dedicated to the development of a computational code that succeeds in simulating transient airwater flows in a unified
framework. In particular, this paper presents an original mathematical model developed for this purpose and its
implementation into a computational code. The present research shows that the driftflux model constitutes a reliable basis to
derive multiphase models, especially to describe airwater interactions in hydraulic structures and environmental flows. An
original onedimensional bilayer driftflux model is then obtained by areaaveraging the threedimensional driftflux model
over a flow cross section presenting a freesurface. Next, an original pressure term extends the applicability of the model to
pressurized flows as well. This original mathematical model unifies and improves the description of transient mixed flows
with airwater interactions. It also serves as basis to develop a computational code called IMPack (Integrated MultiPhase
Pack). IMPack relies on a finite volume discretization and a flux vector splitting. Validity of the computational code IMPack
is assessed by comparison with experimental, analytical and numerical results. IMPack increases the security and reduces the
costs in the design of hydraulic structures. The code also improves the reliability of environmental impact study.
Introduction
The accurate evaluation of all the dynamic phenomena
appearing in a mixed flow is an industrial necessity.
Characterized by a simultaneous occurrence of free surface
and pressurized flow, mixed flows are indeed frequently
encountered in civil and environmental engineering.
Evocative examples are the rapid filling/emptying of water
pipe networks, storm water storage pipes (Vasconcelos et
al. 2005), flushing galleries, sewer systems (Kerger et al.
2009),... Such a flow presents three dynamic features
difficult to predict: an hybrid behavior, airwater
interactions and subatmospheric pressure. First, a mixed
flow displays an hybrid behaviour, as the motion is
governed mainly by the gravity in the freesurface part and
by the pressure in the pressurized part. Consequently, the
wave celerity ranges from about Im/s to about 1000m/s. In
particular, the celerity sharply varies at the transition point.
This feature strongly affects the dynamics of the flow.
Second, airwater interactions may arise, especially at the
freesurface/pressurized transition bore. For instance, the
water inflow may squeeze the air phase initially present
within the pipe (Leon 2007). In large water network and
pipelines (Figure 1), large air pockets may also accumulate
at high points of the pipe. Presence of both dispersed air in
the water and pure air above the freesurface strongly
affects the flow behavior. Finally, subatmospheric
pressure may appear in the pipe. Water hammer or very
high velocity in low points of the pipe usually explains
such a negative pressure. It results in the dilatation of the
water and the contraction of pipe without apparition of a
freesurface. Subatmospheric pressurized flows may only
appear if the aeration device does not supply enough air to
enable the apparition of a freesurface. In the cases where
the pressure drops below the liquid vapor pressure,
cavitation is likely to appear. In conclusion, a correct
description of onedimensional transient flows must
account for mixed patterns, airentrainment and
subatmospheric pressure. Inaccurate design has indeed
caused severe damage to existing structures (Guo et al.
1990).
Figure 1: Large air pockets may accumulate at high points
in pipelines and forms mixed flows.
A lot of research since the 60's has been aiming at
establishing a onedimensional model for freesurface,
pressurized and mixed flows. With variables success,
different approaches tackled the problem of overcoming
the discrepancy between pressure gradients appearing in
freesurface and pressurized flow equations. What is more,
various single and multiphase models have been
proposed to take into account airwater interactions. The
three following paragraphs review pros and cons of these
methods.
Describing mixed flows requires overcoming the
discrepancy between pressure gradients appearing in
freesurface and pressurized flow equations. In this respect,
3 families of methods are worth a mention. First, the
socalled shocktracking approach consists in solving
separately freesurface and pressurized flows through a
different set of equations (Cardle et al. 1988; Politano et al.
2007). Transitions between the freesurface and the
pressurized flow are tracked explicitly and regarded as
internal boundaries across which appropriate jump
conditions are imposed. The advantage of the approach is
that the transition is computed as a true discontinuity
(infinite resolution). However, the associated algorithms
are very complex and become impossible to implement for
complex wave interaction and multidimensional problems.
Second, the rigid water column approach (Li et al. 1999)
treats each phase (air/water) separately on the basis of a
specific set of equations. If the description of the gas phase
is highly accurate, many assumptions underlie the
simplification of the water phase mathematical model.
Applicability of the method is thus limited on specific
cases. Finally, the socalled shockcapturing approach (by
analogy with shockcapturing numerical schemes)
computes pressurized and freesurface flows by using a
single set of equations. The main advantage of the method
is that there is no need to track explicitly the interface. To
authors' knowledge, only four models fall into this
category: the Preissmann slot (Preissmann 1961), the
twocomponent pressure approach (Vasconcelos et al.
2006), the dual model (Bourdarias et al. 2007) and the
kinetic model (Bourdarias et al. 2008). Despite many
advantages, shockcapturing models are known to suffer
various shortcomings (Politano, Odgaard et al. 2007):
simulating subatmospheric pressurized flows is not
straightforward and not a single model takes into account
the presence of air.
Different mathematic models are usually used to simulate
transport of air in civil and environmental engineering but
none of them gives full satisfaction. To date, the use of
multiphase mathematical model for this purpose has
remained circumscribed to very few attempts. As
mentioned by Spasojevic and Holly about sedimentation
engineering (ASCE 2008), "the twophase flow approach
seems promising" but "the formulation of governing
equations in flowsediment problems are still in their
infancy". To author's knowledge, no reference treats of the
air entrainment in civil engineering. Instead, NavierStokes
equations (NVS) constitutes the basis of most classical
models in hydraulics. A transport equation usually
complements the various pure water models derived from
the NVS. This approach aims at describing the evolution of
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
the concentration of air. Nevertheless, such models fail in
many respects to give an accurate description of transient
airwater flows. It does not describe rigorously the impact
of the dispersed phase on the water flow. It does not
characterize correctly the mechanism of dispersion within
the water flow. Finally, the definition of the concentration
is ambiguous.
Local Instant Formulation
Eulerian Time Averaging
IV
Introduction of Macroscopic Variables
Mixture properties Phase properties
DriftFlux Model: Twofluid Model:
3 field equations 4 field equations
interfaciall conditions
_ Homogeneous
Equilibrium Model
Figure 2: Driftflux model is obtained by timeaveraging
the Local Instant Formulation
Finding thorough twophase theories requires investigating
literature in chemical and mechanical engineering. Various
multiphase models are in frequent use in these fields. Most
of them rely on the Local Instant Formulation (LIF). and
all have in common a first step of averaging (Ishii et al.
2006). The most widelyused group of averaging methods
in continuum mechanics is the Eulerian averaging, because
it is closely related to human observations and most
instrumentations (Verloop 1995). In the timeaveraged
Eulerian framework, models can be further classified
(Figure 2) into two separate groups of approach, namely:
the DriftFlux Model (DFM), and the TwoFluid Model
(TFM). The driftflux model assumes the multiphase flow
is a singlephase flow of mixture variables, which refer to
the motion of the centre of mass of the system. The motion
of the dispersed phase is treated in terms of diffusion
through the mixture. The model includes a continuity
equation for the mixture, a momentum equation for the
mixture and a continuity equation for the dispersed phase.
The two latter equations contain additional terms
representing the effect on the mixture flow of the relative
velocity between phases. Governing equations,
implementation, and application of the model is available
in numerous articles (Masella et al. 1998; Romate 1998).
On the other hand, the multiple fluids model (or twofluid
model) treats each phase or component as a separate fluid
with its own set of governing balance equations (Chahed et
al. 2003; Ishii&Hibiki 2006). Each phase has its own
velocity and its own pressure. In this case, the velocity
difference does not require introducing a constitutive
equation. Computation exploits two different momentum
equations, one for each phase. However, difficulties arise
in specifying jump conditions at the interface as well as in
deriving constitutive equations (Issa et al. 2006; Ullmann
et al. 2006). In conclusion, the model should be more
appropriately take into account the dynamic interactions
between phases. It is then useful for applications which
require a high level of accuracy in the treatment of the
twophase interactions (De Henau et al. 1995; Gudmunson
2005). As already mentioned, none of these models has
been use in the field of environmental and civil
engineering. What is more, onedimensional multiphase
models only consider pressurized flows. They fail to
describe the dynamics of freesurface flows.
This research is dedicated to the development of a
computational code that succeeds in simulating transient
airwater flows in a unified framework. In particular, this
paper presents an original mathematical model developed
for this purpose and its implementation into a
computational code. In the rest of the text, the driftflux
theory is proven to constitute a reliable basis to describe
the motion of a dispersed phase in a midscale water flows.
Driftflux model addresses indeed most of the
shortcomings identified in the singlephase approach.
Second, computational effort is decreased by deriving an
original onedimensional bilayer driftflux model. This
model unifies and improves the description of freesurface
flows with airentrainment. Third, an original pressure
term extends the applicability of the model to pressurized
flows as well. This original mathematical model is then
discretised by a finite volume scheme and implemented in
a computational code called IMPack (Integrated
MultiPhase Pack). Finally, validity and applicability of the
IMPack is assessed by comparison with experimental,
analytical and numerical results. In this process, particular
attention is given to benchmarks relevant in civil and
environmental engineering. This model increases the
security and reduces the costs in the design of hydraulic
structures. The code also improves the reliability of
environmental impact study.
Mathematical Model: The Driftflux Theory
In principle, a rigorous multiphase flow solver should
solve at the same time the local instant variables describing
the behavior of each phase by means of the local
NavierStokes equations (holding for each phases
separately) as well as jump conditions (expressing the laws
of conservation across interfaces) (Ishii&Hibiki 2006).
Obtaining a solution this way is however beyond the
present computational capability for many engineering
applications. Various simplified models are thus used in
practice. As pointed in this section, the driftflux model
constitutes a remarkable theory to describe water flows
with airentrainment. In particular, the driftflux theory
addresses most of the shortcomings identified in the
singlephase approach classically used in civil engineering
(Navier Stokes equations + transport equation).
The threedimensional driftflux model is obtained by
timeaveraging the local NavierStokes equations and the
jump conditions. In this process, it is assumed that the
multiphase flow may be described as a single phase flow
of mixture variables which refer to the motion of the centre
of mass of the system. The motion of the dispersed phase is
then treated in terms of diffusion through the mixture
(Ishii&Hibiki 2006). The momentum equation for the
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
dispersed phase is neglected in favor of a constitutive
equation for the relative velocity between the centre of
mass and each phase. We refer the interested reader to the
following books (Ishii&Hibiki 2006; Kerger et al. 2010a)
for the whole demonstration. At the end of the day, the
threedimensional driftflux model is given by three field
equations:
Pm+V(mvm)= 0
at
%Pd+ V (adPdVm)+V(dPdVt)) Fd
a (1)
89pm + V(p vmv )
at
= VPm +v(Tm +TT +D )+Pmg+M,
where ad is the void fraction of the dispersed phase, pm is
the mixture density, Vm is the velocity of the centre of mass
of the system, which is different from both the water
velocity and the dispersed phase velocity. The drift
velocity Vd, is defined as the relative velocity with respect
to the mass center of the mixture. The mass source term Fd
accounts for the exchange of mass between the water and
the active dispersed phase and is usually given by a
casespecific correlation. The mixture pressure pm is a
primitive unknown. The mixture momentum equation
includes three kinds of stresses: the classical Newtonian
viscous stresses Tm, turbulent stresses zT and diffusion
stresses zD due to the relative velocity between phases.
Finally, the mixture momentum source term Mm represents
the effect of the surface tension on the mixture momentum.
The driftflux model presents several advantages over a
classical transport model that we discuss in this paragraph.
The model keeps the simple form of three partial
differential equations for three primitive unknowns. These
equations still translate the momentum and mass
conservation. The mixture density pm is variable even if
both phases are assumed incompressible. Density
variations result from the unsteadiness of the void fraction
ad as follows:
Pm = (1 a)Pw +aPd (2)
where pw is the water density and Pd is the dispersed phase
density (air in this case). Similarly, the mixture velocity Vm
is neither the velocity of the water Vw nor the velocity of
the dispersed phase Vd. It is the velocity of the center of
mass of the mixture formed by both the water and the air.
It is thus defined as:
A (1 a)pwV +adPdVd (3
vm = (3)
Pm
Consequently, the mixture continuity does not ensure the
continuity of the water, but the continuity of both the water
and the air taken as a whole. Adding the diffusion equation
ensures the conservation of the water and the dispersed
phase taken separately. Diffusion equation accounts for the
conservation of the dispersed phase only and governs the
evolution of the void fraction ad. This parameter is
rigorously defined as the probability to find the dispersed
phase at a given point (Ishii&Hibiki 2006; Kerger et al.
2010b):
ak Mddt = [At (4)
At] [At]
where Md is the density function whose value is one in
presence of the dispersed phase and zero in any other case.
The time interval is broken down into [At]=[At]w+[At]d
where [At]w is the set of time intervals in which the
characteristics of the water dominate and [At]d is the set of
time intervals in which the characteristics of the dispersed
phase dominate. Void fraction may be considered as a
concentration in dispersed phase. Besides, most of the
experimental gauges accord with this definition of the
concentration. This is a big advantage of the driftflux
model: the concentration is rigorously defined. Next,
driftflux model accounts for all form of diffusion and
relative velocities thanks to the notion of diffusion
velocity:
V A 1v. (5)
Vd, Vd Vm (5)
Diffusion velocity may have any form and must be given
by a casespecific constitutive equation. For this purpose,
many formulations have been derived in the literature from
experimental data or theoretical consideration on the
dispersed phase momentum equation (Hibiki et al. 2003).
Finally, the mixture momentum equation integrates the
contribution of all kinds of stresses that affect the flow
mixture. They may be conflated in three contributions:
The timeaveraged viscous stresses which integrate
the contribution of both phases:
,. = adTd + (ad )r, (6)
The stresses originating from the difference in
velocity between both phases:
Tr = dPdVdmVdm +(1 d)PwVwmV (7)
The value of the diffusive velocity in the momentum
and diffusion equations are given by the same
constitutive equation.
The turbulent stress due to the macroscopic effect
of the local instant variations in both phases:
TT T T
Td=td +(1 d)W (8)
The mixture momentum source Mm.that represents
the effect of the surface tension on the mixture
momentum, which is totally neglected in the
transport model.
Naturally, constitutive equations must determine these four
contributions to headlosses in multiphase flows. Most of
these laws originates from empirical data and remains
casespecific. It is then of no interest to make here an
extensive review about the literature treating of this topic.
In conclusion, the threedimensional driftflux model
constitutes a reliable alternative to singlephase models in
order to describe transport phenomena in midscale
applications. This model overcomes indeed most of the
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
shortcomings identified in classical models: it describes
accurately the transport of the dispersed phase, it takes into
account the effect of the air on the water flow, the
concentration is rigorously defined, and the transport
equation appears naturally. For the present project, this
conclusion means the driftflux theory may be used to
derive a onedimensional system of equations that
improves and unify the description of freesurface,
pressurized and mixed flows.
Mathematical Model: Onedimensional Model
Like in singlephase models, computation of a
threedimensional multiphase model frequently requires a
prohibitive computational effort. Difficulties arise from the
complexity of treating threedimensional grid and of
tracking the freesurface. In many practical applications in
civil and environmental engineering (river, channel,
pipe,...), onedimensional models are sufficient because
both the flow depth and width are smaller than the flow
length. In this respect, the single areaintegrated driftflux
model that is proposed in the literature only describes
pressurized flows (Wallis 1969; Ishii&Hibiki 2006). It
does not meet our objective to describe freesurface and
mixed flows. What is more, its application has remained
circumscribed to predict the gas void fraction in small
scale gasliquid flows. Civil and environmental engineers
never made use of this model (except a theoretical
contribution by (Wu et al. 2000)). Consequently, the
threedimensional driftflux model is used herein to derive
an original onedimensional bilayer model. The
threedimensional driftflux model is integrated over a
general crosssection presenting a freesurface (Figure 3).
The integration takes into account two layers: a lower one
and an upper one. The resulting model improves the
description of transient freesurface flows with
airentrainment.
< Upper layer: air only
V Area Q
y Lower layer: airwater mixture
Area ,Q
Figure 3: The integration domain consists in a general
crosssection presenting a free surface and two layers.
The lower layer contains a mixture of both water and
dispersed air. It is contained between the rigid bed and the
mobile freesurface. In the integration process of the
threedimensional driftflux model over this layer, we
make two further assumptions. First, we neglect
momentum equations along the flow depth and width. It
means that motion of the flows transversally to the main
stream have a negligible impact on the flow. Second, the
pressure distribution in the flow is taken hydrostatic. The
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
integration gives the following set of three partial
differential equations:
The mixture continuity equation:
8(Pm^ {Pm umRQ
+ = 0 (9)
at ax
The diffusion equation for the dispersed phase:
8(a,)} ~ p
t ax ( (10)
=(x)o
The mixture momentum equation:
a(pm)UmSmQ (pm) )UmUm +(pm )gP
at 8x
+ 9( (ad} dw fj dp (11
+ l(ad) (11)
=(Pm) gQ(So S)+(pm)gP g P Sr
Cax
where Q is the area of the lower layer and the bracket
symbol (f) designates the areaaverage of the function f:
(f)A J fdA (12)
"Q
Equations (9) to (11) account for the mass and momentum
conservation of both phases. The mixture continuity
equation has a very similar form to the continuity equation
in the SaintVenant model (Cunge et al. 1964). However,
the driftflux model describes the evolution of mixture
parameters. Continuity equation is thus expressed in terms
of mean density of the mixture (p) :
(Pm) (ad)Pd + (1 ())Pw (13)
and mean mixture velocity urn:
u. ((14)
Obviously, the mean mixture velocity is chosen as the
mixture density weighted areaaverage of the mixture
velocity. Because of this choice of variables, the continuity
equation ensures the conservation of both phases as a
whole but not separately. The conservation of each phase
on its own is ensured thanks to the diffusion equation (10).
This equation presents several advantages over a simple
transport equation. The concentration is precisely defined
as the mean value of the void fraction over the layer:
CA (ad). A single term depending on the areaaveraged
driftvelocity Udj governs the diffusion of the air in the
water flow. This term accounts for the relative velocity
between both phases. It is defined as:
1 A ((ud))(j) with ((ud)) A
( ad)
where (j) is the areaaveraged volumetric flux defined as
follows:
(j)= () ((ud)) + (1 (ad)) ((u)) (16)
Driftvelocity is intimately linked to diffusion velocities
but constitutes a most convenient parameter to determine
(Brennen 2005). It is given by a casespecific constitutive
equation (Manninen et al. 1996; Hibiki&Ishii 2003). This
formulation of the diffusion is obviously rigorous and
consistent to describe the transport of air in the water flow.
Finally, mixture momentum equation (11) relies also on the
mixture parameters defined above. The equation ensures
the momentum conservation of both phases, but not of
each phase taken separately. So is the topographic slope, a
traditional term in freesurface hydraulic. What is more,
three pressure terms appear in equation (11):
* The hydrostatic pressure term P, :
h,
P ()= (h )(x, )d
h,
Its value depends only on the pipe geometry and on the
lower layer crosssection 2. It can be easily computed
under an analytical form or by a numerical integration.
It results in a discretised relation between the area 2
and the pressure term PQ which is easily implemented
in the numerical scheme.
* The sidewall reaction to section variations gP, :
P(, () (h ) (x d
hx
It accounts for the nonuniformity of the pipe.
Consequently, it cancels in uniform section such as
considered in this paper. We refer the interested reader
to specialized literature for its treatment (Goutal et al.
2002).
S The pressure acting on the freesurface pr that is a
primitive unknown of the upper layer equations.
4ii
h
.................................
hb
...................................................
Figure 4: significance of the notations in the pressure
terms (17) and (18)
In equations (17) and (18), 1 is the lowerlayer freesurface
width, h, is the freesurface elevation of the lowerlayer
and hb is the bottom elevation of the pipe. The total height
h= h,+ hb and is the variable of integration (Figure 4).
Finally, internal and external friction terms are conflated in
two terms. SF accounts for the frictional head loss in the
lower layer. It is computed by means of a twophase
friction correlation such as the Homogeneous Colebrook or
the MartinelliLockart formulations (Lockhart et al. 1949;
Wallis 1969). Sr accounts for the friction at the interface.
This last is usually depends on the difference of velocity
between both layers.
On the other hand, the lower layer contains only pure air.
The layer is enclosed between the freesurface and the pipe
crown. In case of pressurization of the flow, this layer
vanishes. In the integration process of the
threedimensional driftflux model over this layer, we also
make further assumptions. As in the lower layer, we
neglect momentum equations along the flow depth and
width. Second, the void fraction is set at ad=l because the
flow is pure air. Finally, the pressure distribution in the
flow is assumed constant and equal to pr. Since the
diffusion equation cancels, the integration gives the
following set of two partial differential equations:
The pure air continuity equation:
8pgpg apg ((u)g o) g
t =0 (19)
ct Cx
* The pure air momentum equation:
at x gP 9x
= pgggSg +S,
where Qg is the area of the upper layer, and pg is the air
density. The pure air velocity ((ug)) is a primitive
unknown, as well as the pressure pr of the layer. Sg is the
head loss due to the friction in the upper layer. It can be
computed by means of the DarcyWeisbach equation.
Finally, the upper layer equations are closed by using a
equation of state defining the speed of sound in air denoted
a in this paper:
a'2 d(p) (21)
d(pQ)
In some approximation, the air sound speed can be
assumed constant a=330m/s.
In conclusion, equations (9) to (21) define an original
onedimensional bilayer model. This model improves the
description of transient freesurface flows with
airentrainment: it describes accurately the transport of the
dispersed phase, it takes into account the effect of the air
on the water flow, the concentration is rigorously defined
and the diffusion equation appears naturally. However, this
model is still unable to describe pressurized and mixed
flows. For this purpose, we will modify the pressure term
P, in the next section.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Mathematical Model: Extension of the
Freesurface Model to Pressurized Flows
The hydrostatic pressure term P, as defined in equation
(17) does not account for the occurrence of pressurized
flows. As pointed in this section, shockcapturing models
propose to adapt its formulation in order to create a unified
mathematical model for freesurface, pressurized and
mixed flows. However, classical methods does not account
for airentrainment. Consequently, we establish herein an
original mathematical formulation for P, that extends its
applicability to both pressurized flows and aerated flows.
Classical shockcapturing approaches for mixed flows,
such as the Preissmann slot, extend the applicability of
freesurface equations by making the most of the similarity
between the SaintVenant equations (for pure water
freesurface flows) and the Allievi equations (for pure
water pressurized flows). In these last one, the pressure is
linked with pQ by the pressure wave celerity c:
c2 = 2 d(pgP, (Q)) (22)
d(pQ)
The hydrostatic pressure can be thus written as follows:
g (h )l(x,)d if FS
gP. (0)= hb 1 2 (23)
gI .(maxFS)+g f if PP
where 4,a is the value of 2 at the pipe crown, the State
FS indicates that the flow is at freesurface and the State
PP that the flow is pressurized. If 2 > 4x, the flow is
obviously pressurized PP. If 2 < ,, the flow is chosen
either pressurized PP or freesurface FS according to the
aeration rate. Lack of aeration device prevents indeed that
a freesurface appears when the pressure drops below the
atmospheric pressure. Instead, it is observed a
subatmospheric pressurized flow. In classical
shockcapturing methods, the celerity is assumed constant
along the computation, as well as the section 12, such that
integral (23) becomes trivial (Vasconcelos, Wright et al.
2006):
Q c2
gQ i = c gc(Q mx) (24)
where co is the celerity in presence of pure water and at
atmospheric pressure.
In this paper, we propose to use a nonconstant value for
the pressurewave celerity. According to (Guinot 2001),
the celerity for gasliquid mixture is given by the following
formulation:
/ 1 / 1+i/
cm = Co/ 1+caoPm0co P/ P (25)
Pmo = aoPd +(1 )Pw
where the subscript 0 designates the reference state
characterized by a pressure p=101325Pa (atmospheric
pressure). / is a coefficient equal to 1.0 for isothermal
processes and 1.4 for adiabatic conditions. The reference
void fraction ao is the volume fraction of dispersed air at
reference pressure and the p,ro is the mixture density at
reference pressure. By analogy with (Guinot 2001),
integration of (23) by using the airwater formulation (25)
for the pressure wave celerity gives:
Q 2
p pmg f (26)
where the pressure p is is determined by solving iteratively
(Newton Raphson method) the following equation:
a0
PPo+ ( PO oPm oa0p
In conclusion, formulation (23) of the pressure term
coupled with the iterative equation (27) extends the
applicability of the freesurface driftflux model to
pressurized flows as well. The new formulation accounts
naturally for subatmospheric pressure. Originality of the
approach relies in the incorporation of air effect in this
formulation. The combination of this result with the
onedimensional drift flux model for freesurface flows
constitutes an original mathematical model that unifies and
improves the description of transient airwater flows. Since
this model does not offers analytical solution, governing
equations must be solved in a discretised fashion on a
computer. For this purpose, continuous equations must be
discretised and implemented in a computational code. This
is the topic presented in the next section.
Computational Code: IMPack
The original mathematical model proposed in this paper
does not have any analytical solution. However,
discretizing equations and implementing them in a
computational code enable to solve the model for practical
case. It results in the computational code IMPack
(Integrated MultiPhase Pack). As pointed in this section,
IMPack relies on a finite volume discretization and a flux
vector splitting.
Ui1 Ui Ui+,
Figure 5: Finite volume method is popular for solving non
linear hyperbolic problems.
Finite volume method is popular for treating nonlinear
systems of partial differential equations. Its capacity to
ensure the conservation of mass and momentum explains
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
partly this popularity. Its ability to treat discontinuities that
are likely to arise in nonlinear hyperbolic problems is
another reason. In this paper, the onedimensional
driftflux model is solved over a uniform grid containing
cells of length Ax,=x,+/2 x,11/2, i=1,N (Figure 5). The
explicit updating conservative formula is given by:
U"1 U" At FF (U) (U
(28)
At P(J")+AtS:
Ax1
where i is the spatial index, n is the time index and At is
the time step of the temporal integration. The vector of
conservative unknowns Ui is given by:
U, = {Pmum
{Pg g }1
The conservative numerical flux F, +/2 is computed with a
Flux Vector Splitting (FVS) originally proposed for
shallow water equations in (Dewals et al. 2008; Erpicum et
al. 2009a; Erpicum et al. 2009b) and extended here for the
original model:
F
I+
{Pm mm}t +{ Ypl2
{P + +}, 2
tfpgug 1 gj
'+ Y2,)~
where L designate de upstream reconstructed value and R
the downstream reconstructed value (Figure 4). Upstream
and downstream sides are determined by analyzing the
sign of the mixture velocity Um for the first three equations
(describing the lower layer) and the sign of the gas phase
velocity ug for the last two equations (describing the upper
layer). On the other hand, the nonconservative numerical
flux P1+1/2 is given by:
0
0
2, )(tp I
p 0.5*({iR {}R (31)
0
0.5* ({ y (p r 1 Pr I)
Finally, all the other terms are conflated into the source
term S1. Consistence, order of accuracy and stability of the
method has been extensively studied. Constant
reconstruction gives a 1st order of accuracy in uniform grid
(0th order for nonuniform grid). Linear reconstruction
increase the accuracy of one order. The extensive analysis
is beyond the scope of this paper.
Discretised equations are implemented in a homemade
computational code named IMPack for Integrated
MultiPhase Package. Fortran routines constitute the core
of IMPack. IMPack is an academic code that results from
various researches aiming at describing midscale
environmental flows with multiphase models. It concerns
not only airwater flows, but also sediment transport and
pollutant dispersion (Kerger, Dewals et al. 2010a).
In conclusion, the discretization of the driftflux model and
its implementation in IMPack constitutes a useful tool for
modelling transient airwater flows. In particular, it may be
used to design hydraulic structures and supply impact
study in environmental engineering. Validating IMPack is
the last step of development, presented in the next section.
Validation
The original mathematical model presented in this paper
and its implementation in computational code IMPack
have been validated in several single and twophase cases
by comparison with experimental, numerical and analytical
data. Benchmarks include freesurface, pressurized and
mixed flows. Whereas additional ones can be found in
publications in preparation (Kerger 2010), only three
validation cases are presented here: a twophase water
hammer, a stratified flow, and a mixed flow. They all
underline the validity and accuracy of IMPack.
Twophase Water Hammer
This standard benchmark consists in the simulation of a
waterhammer in airwater mixture. This case assesses the
capacity of IMPack to correctly describe the propagation
of shock and rarefaction waves. The evaluation of a correct
celerity is a crucial point here. In the following, we
describe the case and compare numerical results with
reference data (Leon et al. 2008).
A reported in (Leon, Ghidaoui et al. 2008), the present test
considers a horizontal frictionless pipe connected to an
upstream reservoir and a downstream valve. The duct is
10000m long and has a diameter of 1m. The upstream head
(in the reservoir) is kept constant at 200m while the initial
steadystate discharge is 2.0m3/s. The fluid is an airwater
mixture. At the reference pressure (101325Pa), the
concentration in air is 0.2%. The reference pressure wave
celerity is 1000m/s. A sharp transience is created in the
pipe by closing instantaneously the downstream valve. The
discontinuity propagates from the downstream end to the
upstream end of the pipe. It reflects at the reservoir. Since
the celerity depends on the pressure head, the water
hammer takes the form of either a rarefaction or a shock
wave according of the direction of propagation.
Numerical computation is performed with IMPack on a
1000cells uniform grid of 10m. Time integration is
performed with a 3 step RungeKutta scheme (RK31) and
the Courant number is kept constant at 0.5. The
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
reconstruction is constant such that the scheme is 1st order
accurate.
Reference data comes from (Leon, Ghidaoui et al. 2008). It
consists in the nearexact solution of a twophase water
hammer model (2 equations) computed using a Method Of
Characteristic (MOC) scheme for a very fine grid.
The timeevolution of the pressure head at the middle of
the pipe computed with IMPack shows a good accordance
with the NearExact Solution (Figure 6). In particular, the
rarefaction wave clearly appears for the propagation of the
depression. On the opposite, the creation of a shock wave
is clear when the pressure arises. What is more, diffusion
appears and dissipates the waves along the time. Since no
physical friction is considered in the case, it means the
numerical scheme is affect by numerical diffusion.
500
Pressure Head (m)
0Time (s) 140
Figure 6: The timeevolution of the pressure head at the
middle of the pipe computed with IMPack shows a good
accordance with the NearExact Solution of (Leon,
Ghidaoui et al. 2008)
500
Pressure Head (m)
IMPack
Numerical
Solution
250
*LEON
NearExact
Solution
0
0 Position (m) 1000
Figure 7: Comparison of the profiles of the pressure head
after 140s computed with IMPack and the NearExact
Solution of (Leon, Ghidaoui et al. 2008) shows that
numerical solution is partly smeared out.
Comparison of the profiles of the pressure head after 140s
computed with IMPack and the NearExact Solution shows
that numerical solution is partly smeared out. This is still a
consequence of the numerical diffusion that affects
IMPack. However, the graph clearly underline the ability
of the code to propagate shock and rarefaction waves at a
correct velocity.
Stratified Flow
A movement of the freesurface in a closed pipe should
lead to reaction of the air in the upper layer in terms of
pressure and discharge. Since experimental data for such
flows are not available, we propose here to illustrate such a
situation with an academic case.
The circular pipe is 100m long and its diameter is 2m. The
pipe is assumed frictionless and without slope. The initial
pure water discharge in the lower layer is 0.01m3/s and the
initial water height is H=1.703m. The void fraction in the
lower layer is set to zero. In the upper layer, the initial
pressure is 140861Pa because the air density is assumed
equal to 1.29349 kg/m3 and the celerity is equal to 330m/s.
The air is initially still. At the left boundary, we impose an
impermeable boundary for the air (ug=0m/s) and water
discharge is linearly increased from 0.01m3/s to 0.5m3/s
during 25s. From then the discharge is maintained constant.
At the right boundary, the pressure is fixed at 140861Pa
and the water height is assumed constant at 1.703m.
Numerical computation is performed on a 100cells
uniform grid of im. Time integration is performed with a 3
step RungeKutta scheme (RK31) and the Courant number
is kept constant at 0.6.
Computed pressure and discharge of the water flow (lower
layer) at the middle of the pipe show the evolution of the
freesurface (Figure 10). We clearly observe, in a first time,
a surge in height and discharge in the lower layer, followed
by a decrease in both height and discharge after reflection
on the right boundary. Such variations are classical in
freesurface hydraulics and could be computed with the
SaintVenant equations.
1
Discharge (m3/s)
18
Pressure head (m)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
the upper layer to the temporal variation of the freesurface.
We observe an increase in pressure up to 1.4atm followed
by a decrease up to 0.1atm. This is a consequence of the air
entrapment. The possibility of negative or subatmospheric
pressure is herein clearly observed. Such a flow dynamics
may only be describe by a twolayer approach as
introduced in this paper.
Discharge (mn/s)
01
0
01
21
Pressure head (atm)
14
02 1 0
0 Time(s) 100
Figure 9: Computed pressure and discharge of the air flow
(upper layer) at the middle of the pipe exhibits a temporal
evolution due to the variations of the freesurface
Mixed Flow: Fast Closing Gate
The simulation of the experimental test J of Cardle (Cardle
et al. 1983) aims at demonstrating the ability of IMPack to
simulate freesurface/pressurized transitions. The
experimental apparatus of Cardle (University of
Minnesota) consists in a PVC pipe of 48.77m length. Its
inner diameter is 16.6cm. The pipe presents a 0.05% slope
and a Manning coefficient of 0.011. Initial condition is a
steady flow whose discharge is Q=0.005097m3/s. The
downstream boundary condition is a pressure head of
H=0.1372m. The fast closing of the downstream gate
creates a sharp transition bore that propagates upstream.
When this interface reaches the middle of the pipe (24.4m),
the downstream gate is reopened. It reverses the direction
of the bore propagation. What is more, a second transition
bore forms at the downstream end of the pipe and
propagates upstream. The pipe is assumed perfectly
aerated such that a subatmospheric pressurized flow
cannot appear. The fluid flow is assumed to be pure water.
Water Height (m)
0 1 1 62
0 Time(s) 100
Figure 8: Computed pressure and discharge of the water
flow (lower layer) at the middle of the pipe show the
evolution of the freesurface.
On the Figure 9, computed pressure and discharge of the
air flow at the middle of the pipe exhibits the reaction of
0 Time (s) 50
Figure 10: Experiment J of Cardle Numerical data
computed with IMPack accords reasonably well with the
experimental data at the pressure transducer P1.
As pointed in Figure 10, numerical data computed with
IMPack accords reasonably well with the experimental
data at the pressure transducer P1 (9.9m from the
downstream end). The time of propagation of each bore is
correctly captured. The pressure at the transition bore is
also correctly assessed. However, the end of the simulation
differs notably from the experimental data. This
discrepancy is reported by other authors. It should be
caused by a depression in the air layer.
Conclusions and Perspectives
Mixed flows, characterized by a simultaneous occurrence
of free surface and pressurized flow, are frequently
encountered in civil and environmental engineering. Single
phase models coupled with a transport equation fail to
account for subatmospheric pressurized flows and
airwater interactions. This paper proves that a
onedimensional driftflux model succeeds in simulating
transient airwater flows in a unified framework. The
demonstration of this point relies on five sections. First, the
driftflux theory is shown to constitute a reliable basis to
describe the motion of a dispersed phase in a midscale
water flows. Driftflux model addresses indeed most of the
shortcomings identified in the singlephase approach.
Second, computational effort is decreased by deriving an
original onedimensional bilayer driftflux model. This
model unifies and improves the description of freesurface
flows with airentrainment. Third, an original pressure
term extends the applicability of the model to pressurized
flows as well. This original mathematical model is then
discretised by a finite volume scheme and implemented in
a computational code called IMPack (Integrated
MultiPhase Pack). Finally, validity and applicability of the
IMPack is assessed by comparison with experimental,
analytical and numerical results. In this process, particular
attention is given to benchmarks relevant in civil and
environmental engineering.
From a more practical point of view, results of this paper
seem to furnish a useful tool to describe environmental
flows. In particular, IMPack increases the security and
reduces the costs in the design of hydraulic structures.
Used in conjunction with an experimental model, it may
indeed provide very interesting results, especially the in
optimization stage. The code also improves the reliability
of environmental impact study. Physically based methods
are often more and more frequently used for risk
management and environmental preservation.
To conclude, the driftflux theory seems a promising
paradigm for the description of transport phenomena. With
a single mathematical model, it seems possible to describe
the transport of sediment, pollutant and any active or
passive dispersed phase. The mathematical model should
improve the accuracy of the description as well as provide
a unified basis for a computational code. It requires
however further research, especially for the establishment
of generic constitutive laws adapted to civil and
environmental engineering.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
References
ASCE Sedimentation engineering : processes,
measurements, modeling, and practice. Reston, Va.,
American Society of Civil Engineers(2008).
Bourdarias, C.&Gerbi, S. A Finite Volume Scheme for a
Model Coupling Unsteady Flows in Open Channels and in
Pipelines. Journal of Computational and Applied
Mathematics 209(1): 109131 (2007).
Bourdarias, C., Gerbi, S.&Gisclon, M. A kinetic
formulation for a model coupling free surface and
pressurized flows in closed pipes. Journal of
Computational and Applied Mathematics 218(2): 522
(2008).
Brennen, C. E. Fundamentals of Multiphase Flows,
Cambridge University Press(2005).
Cardle, J.&Song, C. Mathematical Modeling of
Unsteady Flow in Storm Sewers. International Journal of
Engineering Fluid Mechanics 1(4): 495518 (1983).
Cardle, J.&Song, C. Mathematical Modeling of
Unsteady Flow in Storm Sewers. International Journal of
Engineering Fluid Mechanics 1(4): 495518 (1988).
Chahed, J., Roig, V.&Masberat, L. EulerianEulerian
twofluid model for turbulent gasliquid bubbly flows.
International Journal of Multiphase Flow 29(1): 2349
DOI: Pii S03019 221j ,il1 234(2003).
Cunge, J. A.&Wegner, M. Integration numhrique des
equations d'6coulement de Barrt de Saint Venant par un
schema implicite de differences finies. La Houille Blanche
3339 (1964).
De Henau, V.&Raithby, G. D. A transient twofluid
model for the simulation of slug flow in pipelinesI.
Theory. International Journal of Multiphase Flow 21(3):
335349 (1995).
Dewals, B. J., Kantoush, S. A., Erpicum, S., Pirotton,
M.&Schleiss, A. Experimental and numerical analysis of
flow instabilities in rectangular shallow basins.
Environmental Fluid Mechanics 8(1): 3154 (2008).
Erpicum, S., Dewals, B. J., Archambeau, P&Pirotton, M.
Dambreak flow computation based on an efficient
fluxvector splitting. J. Comput. Appl. lfh i published
online) DOI: 10.1016/j.cam.2009.08.110(2009a).
Erpicum, S., Meile, T., Dewals, B. J., Pirotton,
M.&Schleiss, A. J. 2D numerical flow modeling in a
macrorough channel. Int. J Numer Methods Fluids
61(11): 12271246 (2009b).
Goutal, N.&Maurel, F. A finite volume solver for 1D
shallowwater equations applied to an actual river.
International Journal for Numerical Methods in Fluids
38(1): 119 (2002).
Gudmunson, R. L. 2005. A Numerical Study of The
TwoFluid Models for Dispersed TwoPhase Flow.
Numerical Analysis and Computer Science. Stockholm,
Royal Institute of Technology, Stockholm. Phd: 40.
Guinot, V Numerical simulation of twophase flow in
pipes using Godunov method. International Journal for
Numerical Methods in Engineering 50(5): 11691189
(2001).
Guo, Q.&Song, C. Surging in Urban Storm Drainage
Systems. Journal of Hydraulic Engineering 116(12):
15231537 (1990).
Hibiki, T.&Ishii, M. Onedimensional driftflux model
and constitutive equations for relative motion between
phases in various twophase flow regimes. International
Journal of Heat and Mass Transfer 46(25): 49354948
(2003).
Ishii, M.&Hibiki, T Thermofluid dynamics of twophase
flow, Springer Science, USA(2006).
Issa, R. I., Bonizzi, M.&Barbeau, S. Improved closure
models for gas entrainment and interfacial shear for slug
flow modelling in horizontal pipes. International Journal
of Multiphase Flow 32(1011): 12871293 (2006).
Kerger, F. 2010. Modelling Transient Airwater Flows.
ArGEnCo. Liege, University of Liege. PhD.
Kerger, F, Archambeau, P., Erpicum, S., Dewals, B.
J.&Pirotton, M. Simulation numerique des 6coulements
mixtes hautement transitoire dans les conduites
d'6vacuation des eaux. Houille BlancheRev. Int. 2009(5):
159167 (2009).
Kerger, F., Dewals, B. J., Erpicum, S., Archambeau,
P.&Pirotton, M. Modelling Flows in Environmental and
Civil Engineering: Unification of the Mathematical Theory
and Application to Practical Cases. NewYork, Nova
Science Publishers (2010a).
Kerger, F., Dewals, B. J., Erpicum, S., Archambeau,
P.&Pirotton, M. 2010b. Modelling Flows in Environmental
and Civil Engineering: Unification of the Mathematical
Theory and Application to Practical Cases. NewYork,
Nova Science Publishers 155.
Leon, A. 2007. Improved Modeling of Unsteady Free
Surface, Pressurized and Mixed Flows in StormSewer
Systems. Urbana, University of Illinois at
UrbanaChampaign. Phd: 194.
Leon, A. S., Ghidaoui, M. S., Schmidt, A. R.&Garcia, M.
H. Efficient secondorder accurate shockcapturing
scheme for modeling one and twophase water hammer
flows. Journal of Hydraulic Engineering 134(7): 970983
DOI: Doi
10.1061/(Asce)0733',42'li 'Xsl l 34:7(970)(2008).
Li, J.&McCorquodale, A. Modeling Mixed Flow in
Storm Sewers. Journal of Hydraulic Engineering 125(11):
11701180 (1999).
Lockhart, R. W.&Martinelli, R. C. Proposed correlation
of data for isothermal twophase, twocomponent flow in
pipes. Chemical Engineering Progress 45: 3948 (1949).
Manninen, M., Taivassalo, V.&Kallio, S. 1996. On the
mixture model for multiphase flow. V. t. u. (VTT). Espoo,
Technical Research Center of Finland: 67.
Masella, J. M., Tran, Q. H., Ferre, D.&Pauchon, C.
Transient simulation of twophase flows in pipes.
International Journal of Multiphase Flow 24(5): 739755
(1998).
Politano, M., Odgaard, A. J.&Klecan, W. Numerical
Evaluation of Hydraulic Transients in a Combined Sewer
Overflow Tunnel System. Journal of Hydraulic Research
133(10): 11031110 (2007).
Preissmann, A. F,.',ljr,.', des intumescences dans les
canaux et rivieres. First Congress of the French
Association for Computation, Grenoble, France(1961).
Romate, J. E. An approximate Riemann solver for a
twophase flow model with numerically given slip relation.
Computers & Fluids 27(4): 455477 (1998).
Ullmann, A.&Brauner, N. Closure relations for twofluid
models for twophase stratified smooth and stratified wavy
flows. International Journal of Multiphase Flow 32(1):
82105 DOI: DOI
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
10.1016/j.ijmultiphaseflow.2005.08.005(2006).
Vasconcelos, J.&Wright, S. Experimental Investigation
of Surges in a Stormwater Storage Tunnel. Journal of
Hydraulic Engineering 131(10): 853861 (2005).
Vasconcelos, J., Wright, S.&Roe, P. L. Improved
Simulation of Flow Regime Transition in Sewers : The
TwoComponent Pressure Approach. Journal of Hydraulic
Engineering 132(6): 553562 (2006).
Verloop, W. C. The inertial coupling force. International
Journal ofMultiphase Flow 21(5): 929933 (1995).
Wallis, G. B. Onedimensional Twophase Flow(1969).
Wu, W.&Wang, S. S. Y Mathematical Models for
LiquidSolid TwoPhase flow. International Journal of
SedimentReserach 15(3): 288298 (2000).
