7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Atomization of Stratified Turbulent Flow
E.W. Smith*, R.M. Robertst and M. J. McCready*
Department of Chemical and Biomolecular Engineering, University of Notre Dame, Notre Dame, IN 46556, USA
t Chevron Energy Technology Company, Houston, TX 77002, USA
mjm@nd.edu and roberrml@chevron.com
Keywords: stability, atomization, turbulence, stratified, wavy, OrrSommerfeld, kepsilon, oilgas pipeline
Abstract
A fully computational approach to interfacial stability in stratified turbulent flows is presented. The approach is valid for a
wide range of conditions from laminar to highly turbulent (Re = 107). The basestate is modelled using a kepsilon turbulence
model that has been tailored to stratified flow including waveroughening effects. The stability analysis is performed by
solving a modified OrrSommerfeld equation with a ChebyshevTau method. The determination of atomization and other
phase transitions comes from calculating the growth rates of long and short waves as they provide the mechanisms for
splashing, droplet shearing, and growth of nonlinear waves and slugs. Calculations are in good agreement with data from 1"
channels at atmospheric conditions and 2" pipes at 2 atm with the working fluids being water and air. But more than that, the
calculations determine the growth rate of detectable waves in experiments with finite residence times.
Introduction
Multiphase flows occur commonly in the petroleum
production and transportation infrastructure. Despite many
experiment studies and complete knowledge of the
fundamental governing equations, it is still not possible to
predict with confidence pressure drops, slugging and
atomization conditions and solids carrying capacity for
many important operating ranges. This paper discusses
our new computational efforts to predict the onset of
atomization for gasliquid flows. In addition to altering
the pressure drop flow rate relation, atomized drops can
wet the top of a pipe and hence reduce topline corrosion.
Our approach is to calculate the basestate (turbulent gas,
laminar or turbulent liquid) for both a smooth surface and
waveroughened conditions. Stability analysis is done for
both basestates so that secondary transitions, such as the
onset of roll waves, can be determined. Atomization has
been associated with roll wave onset in some previous
studies and atomization can also occur on thin films if
solitary waves are present. Our results consider horizontal,
inclined and declined flows. This work goes beyond
previous studies in some important respects. First it is
recognized that when comparing calculated and measured
basestates for gasliquid flows, the interfacial roughness
must be considered carefully. Second, the onset of long
waves, since these generally occur at conditions when short
waves are unstable, need to be calculated from this already
rough basestate.
The first question to be addressed is the accurate prediction
of smooth and wavy basestates. A finite difference
method is used to solve a kepsilon turbulence model for the
gas and liquid flows. This approach will be shown to match
available data for smooth and rough interfaces very well.
The physical phenomenon that occurs in the presence of
short waves is a vertical shift of the maximum velocity
away from the interface. This occurs because the wavy
interface dampens momentum, which can be modeled as an
increase in the eddy viscosity at the interface. These
considerations have not been made in the literature
pertaining to stability and yet they are the key initial step in
the investigation of long wave formation.
A stability analysis was performed using an
OrrSommerfeld equation that was modified to include the
basestate velocity and the eddy viscosity. It was found that
using the ChebyshevTau method on the modified
OrrSommerfeld equation at higher Reynolds numbers and
lower wavenumbers became increasingly difficult to resolve
numerically. To check the validity of the solutions being
obtained, a longwave analysis was implemented. The
longwave analysis gives an efficient means of calculating
the neutral stability lines while the ChebyshevTau method
provides information about the growth rates of short waves,
which can be correlated to the magnitude of interfacial
roughness giving a completely theoretical prediction of the
onset of atomization. The stability boundary for long
wave predictions has been compared to data and
consistently underpredicts the roll wave transition as shown
in the middle figure. However, this should be the case
because the data are from channels of lengths no more than
10m whereas the stability analysis produces the neutral
stability curve, which corresponds to waves that would
require an infinite length to develop. A contour plot of
growth rates reveals that the roll wave transition matches
roughly with a predicted rate of about 0.1/s which is
nominally the rate necessary for roll waves to form in the
10m channels where the measurements were taken. Also
in this case, atomization is seen to occur with the presence
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
of roll waves.
Nomenclature
cx Turbulence model constants
c Wave speed (m.s'1)
d Channel height or pipe diameter
f Smoothing turbulence model function
F Froude number
g Gravitational constant (m.s 2)
h Perturbation height (m) (in stability)
h Holdup (in basestate)
m Viscosity ratio
p Perturbation pressure (Pa)
P Pressure (Pa)
r Density ratio
R Reynolds number at interface (Uo.dL.n1)
Re Reynolds number using avg phase velocity
U Basestate velocity (m.s1)
v* Friction velocity (m. s')
T Weber number
x Position in the flow direction
Greek letters
a Wave number (m'1)
y Perturbation variable
F Total viscosity (Pa.s)
E Turbulent dissipation (m2.s3)
k Turbulent kinetic energy (m2/s2)
1T Eddy viscosity (Pa.s)
p Dynamic viscosity (Pa.s)
p Density (kg.m3)
p Stream function (m2.s')
T Shear stress (N.m 2)
numerical considerations.
Because we are interested in the stability of the turbulent,
potentially roughened basestate, we need an approach for
which stability can be implemented. Methods such as
Direct Numerical Simulation or Large Eddy Simulation are
not only extremely time consuming, but it is doubtful that
information from these calculations could readily transfer to
a stability analysis. Alternatively, highly empirical
methods such as zero equation models or mixing length
models are not flexible enough or lack the empirical values
to use with high pressure oilgas flows. Two equation
models were chosen for their application to a wide range of
conditions and also in the flexibility afforded by them to
adjust for the effects of roughness and turbulent effects at
low Reynolds numbers.
The kepsilon, twoequation models are a way of
determining the eddy viscosity, which is part of the
Boussinesq approximation. By modelling the Reynolds
stress as a second viscosity, the momentum equation
simplifies to
dy(a )d ) = Pdx
ay ~ d d jP Slr)
where p, is p, = pjc,, k2/E .
The modeling
parameters c, and f, depend on version of the turbulence
model and Z is the isotropic part of the turbulent
dissipation. The variables, k and E are determined by
the continuum equations
d Ty dk y Pj (E+ D)
kdy k dy
Subscripts
i Interface
j Phase index (L or G)
k Turbulent kinetic energy
L Liquid phase
G Gas phase
e Turbulent dissipation
p Eddy viscosity equation
Superscript
+ Dimensionless friction variable
Basestate
Theory
Here we will describe the approach used to model wavy
flows and give an overview of the numerical methods. A
stability analysis is only a tool that provides information
regarding a particular basestate. Thus without an accurate
model of the base flow there is simply no relevance of
hydrodynamic stability results to a desired physical situation.
To model a stratified turbulent flow, a single phase
turbulence model must be adjusted for the roughening effect
of interfacial waves and also increased turbulence in the
liquid at Reynolds numbers well below the empirical
transition (Re=2300). In addition, there are several
d + LT\ dE\ Lt2 dU
=. r C ) cE2 f2 2 C f r PT L )
dy \ E dy )k k (dy)
PjE,
where again cx, fx, E, and D vary depending on the
model. The various models of Chien (1982), Dutoya and
Michard (1981), Hoffman (1975), Lam and Bremhorst
(1981), Launder and Sharma (1974), and Reynolds (1976)
were tested. The version chosen for this work was based on
the Launder and Sharma model, but with f, taken from
Chien.
Results
The maximum velocity of a fluid will move away from a
wavy interface or a rough wall because there is more
resistance to motion in this region. This is modeled simply
as a nonhomogeneous boundary condition for the turbulent
kinetic energy. The value of the condition is determined
from matching data, which includes profile data when
available or by matching flow rates and pressure drops. It
was found that in general k+ = k/v*2 = 4 gave a good fit,
where v* is the friction velocity. An example of the
success of this approach is shown in figure 1.
/dU \2
PT 
(dy)
0.6
0.4
0
0.2
0 2 4 6 8
Velocity [mrs]
Figure 1: Velocity profile from Cohen and Hanratty (1968)
versus the model of wavy flow
Further evidence of the viability of the twoequation model
is given in table 1 where the pressure drop data calculations
are compared to data from Bruno (1988). While the
condition, k = 4 is useful, it is not expected to be a
completely general result. It could probably be improved
by using the growth rate of the smooth basestate to estimate
a wave height.
ReG ReL Pressure Drop Pressure Drop error
Data [Pa/m] Calculation
9000 720 27.4 28.2 3%
7670 720 20.5 19.3 6%
6100 720 14.7 14.9 1%
11000 745 26.2 22.8 13%
10050 745 23.6 20.5 13%
9000 745 19.0 17.7 7%
7670 745 13.9 12.9 7%
6100 745 9.7 9.4 3%
4750 745 6.4 6.4 0%
11000 255 18.7 15.6 17%
10050 255 16.3 14.4 12%
9000 255 12.6 11.5 9%
7670 255 9.1 8.8 3%
6100 255 6.5 6.2 4%
4750 255 4.3 4.1 4%
Table 1: Basestate comparison with data from 1" airwater
channel data
Before discussing stability, it is worth reiterating the
importance of the basestate. One goal of this work is to
correctly predict flow pattern transitions. This involves
correctly calculating the stabilizing and destabilizing
mechanisms that include gravity, shear stress, pressure, and
surface tension. If the basestate velocity or eddy viscosity
profiles are not correct then the relative magnitudes of these
mechanisms will give erroneous transition lines. The
transitions of interest are that of stratified to wavy, the onset
of entrainment, the initiation of slug flow, and maximum
entrainment. A turbulent, linear, stability analysis provides
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
crucial information for all of these transitions. Obviously,
nonlinear techniques are required for a study of slug flow,
but again such an undertaking is not possible without a
correct linear analysis. Figure 2 illustrates this by showing
how much the wave speeds of long waves depend on the
details of the linear analysis and the underlying basestate.
This has not been available and is why this work poses great
opportunities in producing far more accurate and
fundamentally sound nonlinear analyses. Ultimately, it is
the necessary first step in reassessing flow pattern
boundaries, a subject that has seen little improvement since
Taitel and Dukler over 30 years ago.
0.4
0.3
I0.2
Turbulent
0.1 Laminar
Inviscid
0.0
0,01 0.1 1 10 100 1000
Wavenumber Il/m]
Figure 2: Comparing wave speeds of simple laminar and
inviscid flows versus a complete turbulence model and
stability analysis
Linear Stability Analysis
Theory
The OrrSommerfeld equation is an expansion of the
laminar terms of the NavierStokes equations. Because the
Boussinesq approximation is employed, the expansion can
be repeated with the eddy viscosity term. First, the
velocity is recast by the Reynolds decomposition,
u(x, y, t) = U(y) + yu(x, y, t) and v(x,y,t) = V(y) +
yv(x, y, t), where y is the expansion parameter and is
much smaller than one. The perturbation velocity is
assumed to have the form of a traveling wave, un(x, y, t) =
u (y) exp(ia (x c t)). Use of a stream function means
the perturbation velocities can be written as
M(x,y, t) = 0'(y) exp(ia(x ct))
v(x, t) = ia 0(y) exp(ia(x ct)).
The Boussinesq approximation remains in the momentum
equations, but the eddy viscosity is not decomposed into
basestate and perturbation viscosities. The reasons for this
are twofold, the first is that this term is presumably small
even compared to the other perturbation values (the correct
phase shift of the shear stress confirms this). Also, there is
little reason to believe that the perturbation of a modeling
convenience such as the eddy viscosity accurate as it may
be for the basestate would hold any physical relevance in
the stability analysis. The pressure is also represented as a
sum of a basestate and normal, travelingwave, perturbation,
p = P + yp(y) exp(ia(x c t)) Upon substituting
these expansions the order y' terms for the horizontal and
vertical momentum equations are
dy
dy (r(y)i "(y)) 2r' (y)
rk
= aiR R (p(y) c Y'(y) P(y)U(y)
mk
+ )'(y)U(y))
and
a(a22(y)P(y) '(y)'(y) f(y)v"(y))
= i r (a2C 2 (y) p(y)
mk
a2 U(y)P(y)),
respectively. Differentiating the horizontal momentum
equation with respect to y, multiplying the vertical
momentum equation by a, and adding the two equations,
the pressure perturbation can be eliminated giving
d2 d
d ((y)"(y)) 22 d~ (r(y)('(y)) + a4r(y)(y)
dy2 dy
iRk ((U(y) c)(a2 (y)
mk
"(y)) + U"(y)(y)),
which is a modified OrrSommerfeld equation. Unlike the
standard form, it includes a basestate eddy viscosity as well
as the basestate velocity profile.
In addition to an expansion of the differential equations, the
boundary conditions must also be reevaluated to include the
eddy viscosity. The no slip and zero flux conditions at the
walls remain homogeneous:
4=0
(' =0.
Disturbances in the flow produce perturbations at the
interface. Therefore, the boundary conditions do not occur
precisely at y = 1, but on an undulating interface of the
form h(x, t) = yh exp(ia(x c t), where h is a complex
scalar. The displaced velocities are written as a Taylor
series expansion with respect to h around y = 1 giving
du(1)
u(1 + h) = u(1) + h d
dy
dv(l)
v(1 + h) = +h
dy
Substituting the normal wave expansion into the velocity
matching condition gives
PL = (G
P'L(y) + UL(y)h = b (y) + U'(y)h.
Matching the shear stress requires both horizontal and
vertical velocity components of the Cauchy stress tensor
r = F (du/dy + dv/dx) Expanding and equating this
term across the interface gives
FrL(L' h U (' + aPL) = m2F (PG + h UL' + a2PL).
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Finally, the normal stress is required to match at the
interface. In addition to the velocity gradients, a
perturbation pressure, the capillary pressure produced from
surface tension, and hydrostatic pressure created by an
uneven interface must be included. The normal stress
condition is
(FL dvL
2 Ly
dy
m y R (PL PG)
dy
dh2 (, t)
STRd2
dx2
R F (1 r2)h,
where T is the Weber number, and F is the Froude
number. The pressure for phase j is
P= i (a2p + 1 (3)
ryR
+ a(U + c ().
Substituting the normal wave expansion, the pressure, and
collecting y' terms gives the normal stress condition
3ia2
T a3h aF h (1 r2) + ( ; m2FG')
+ a(#PU' Uj(P rT2ct 2 r222
1 d d
Rd (L L') m2 y (FG )
= ca(P' r2 9).
The system is still underdetermined. There are 8 boundary
conditions, 2 fourth order equations, and an unspecified
interfacial perturbation variable h The kinematic
condition provides the necessary condition for solvability
h = c h.
In practice, one of the boundary conditions is dropped by
solving for h using one of the above conditions. While
much of the literature uses the kinematic condition to
determine h, this leads to an eigenvalue in the denominator
and thus the matrix formulation does not work. Therefore,
the value of h determined by the continuity of tangential
velocity is substituted into all the boundary conditions. It
is possible to keep h as a variable in the vector of
unknowns, and use all 9 boundary conditions. In this case
the vector would include the value of the series coefficients
followed by the perturbation height.
Results
The solution of the eigenvalue problem is to factor the
modified OrrSommerfeld equation into two second order
differential equations and apply the ChebyshevTau Spectral
method. The application of this method is not
straightforward. The difficulty is using enough ( Ii.bI sheIc\
polynomials to resolve the basestate profiles and their
derivatives, but few enough that numerical error remains
insignificant. For moderate flows this condition is easily
met. Figure 3 shows clearly defined growth rate and wave
speed curves for the interfacial modes of a moderate flow of
water and air in a 5" channel.
10 102
wavenumber [1/rnm]
100 102
wavenumber [1/m]
Figure 3: Growth rate and wave speeds of interfacial modes
at ReG=31,000 and ReL=3900. The interfacial modes
correspond to the downstream (solid) and upstream (dashed)
traveling waves
Increasing the inertia necessitates a finer grid with more
points in the basestate, which in turn requires more modes in
the stability problem. The consequence is a narrow or
nonexistent range of modenode combinations that produce
realistic profiles. To broaden the space between poor mode
expansions and machine precision, several numerical tricks
were implemented as described in the thesis of Smith (2010).
The success of these numerical techniques is shown in
figure 4.
When data are unavailable, there are other implicit means of
checking the validity of results. For instance, it is not
sensible to say that upstream traveling long waves will grow
while downstream traveling waves are stable, which is what
the upper plot of figure 4 indicates. For some conditions
well above stability a linear analysis may reveal that
upstream waves are unstable, but this is only valid if the
downstream modes are unstable first. Intuition can serve
as a check in other ways. An example is that the speed
curve should reach a minimum where the growth rate is
highest, this is a wellknown energy tradeoff between
growing versus propagating the wave. Additionally, the
interfacial modes should be mirror images of one another, a
mathematical and physical truism of wave behavior. The
speeds should also increase where the restoring force of
surface tension (large wavenumbers) and gravity (small
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
wavenumbers) are applicable. The bottom curves of figure
4 pass all of these tests and are easily defended as valid
growth rates and wave speeds.
U)
E0.5
,0.5
1
2
E0.5
_0
.5
1
100 102
wavenumber [1/m]
100 10/
wavenumber [1/rm]
Figure 4: For highly turbulent flows (ReG=3.2x105,
ReL=4000) the solution becomes exceedingly difficult. By
judiciously implementing a number of numerical techniques
the mode interactions in (a) can be eliminated (b).
Having developed a stability method that is at the very
minimum defensible (the accuracy will be defended in the
next section), this method can be used to probe parameter
space and determine where instabilities might arise due not
only to flow rates but flow parameters such as viscosity
ratio, density ratio, holdup, pipe diameter, incline, etc. An
example of this type of analysis is given in figure 5. In the
first plot the increase of pressure, while maintaining the
'N"
/
same Gas Reynolds numbers, decreases stability and
ultimately leads to the formation of long waves. Similarly,
changing the flow from a decline to an incline produces a
long wave transition.
P=1atm
 P=2atm
P=1.5atm
01
0.212
102
S0.3 r 9=+0.03
0.2 
0.1
D o0 /
0.1
0.2
10 10 102
a
Figure 5: Long wave transitions produced by increasing a) the
system pressure and b) the inclination
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
was calculated and the resulting growth rate contours are
compared with Bruno's long wave transition in figure 6.
10000
5000
500 1000
ReL
1500 2000
Figure 6: Bruno's experimental long wave transition
(dashed line) as determined by a spectral analysis of wave
tracing is compared to the growth rate contours of long
waves as determined by a stability analysis
We can use growth rate contours to find neutral stability for
highly inertial flows. A program has been developed that
handles such conditions and produces contour plots. Some
data are available that determine the onset of atomization by
visually confirming the presence of droplets in air. Since
these experiments were done in relatively short pipes it
should be expected that growth rates corresponding to onset
would be higher than that of Bruno's spectral analysis. A
plot of growth rate contours and the onset of atomization for
a waterair system in a 2" by 20' pipe is given in figure 7.
The kinks in the lines result from the tradeoff between
precision and speed between automation and a laborious
adjustment of numerical parameters. If need be, trend
lines can be fitted to smooth contours which is shown in
figure 8.
Phase Transition and Atomization
0.1
A motivation of this work is to study the onset of
atomization. A mechanism for this is the shearing of liquid
from the crest of roll waves. Therefore a prediction of
long wave instability should lead to waves that will correlate
with the presence of atomization. Here we compare the
stability calculations of long waves with experimental
values for the transition. The difficulty in measuring the
existence of roll waves is that laboratory experiments are
performed in channels or pipes no longer than about 10 m,
and thus detectable waves are limited by growth rates and
finite residence times. Furthermore, outside of some
heuristics, there is no quantitative method of determining
the growth rate of these realizable waves. A stability
analysis provides a quantitative context for experimental
observations.
At conditions or flow distances when roll waves are not
visible, they can still be detected by a spectral analysis of
the wave tracings. Bruno and McCready (1988) performed
this analysis for waterair experiments in a 1" inch channel.
These data are ideal for testing the stability analysis because
the flow rates are moderate and the analysis provides the
best experimental estimate of neutral stability. The
stability of long waves at various points in ReLReG space
0.01
0.001
0.0001
Use
Figure 7: Growth rate contours of a highly inertial flow and
the experimental onset of Mantilla (2008)
~6)
~ bC
.5 3
.2 \
0.1
.I . . . I
01
001
0.001
0.0001
 \
\ 3.0
0.5 1
).1 1 10 100
USG
Figure 8: Fitted growth rate contours
The stability analysis also includes the growth rates of short
waves. One reason to doubt the applicability of the long
wave mechanism to all systems is the evidence in data that
surface tension has an effect on the transition line. Since
surface tension only appears in one boundary condition and
at order a3 where a is the wavenumber, it is likely that
long wave instability is not the only criterion for
atomization. Therefore, the short wave growth rate
contours can be compared to data with the purpose of
finding a growth rate threshold that corresponds to the
observed onset of atomization. An example of this
analysis is given in figure 9. The short wave analysis has
only recently been considered, but it is presented here as an
example of information that is only available with a
comprehensive approach to turbulent stability analysis.
I
0.1
0.01
0.001
0.0001
Figure 9: Short wave growth rates
experimental onset line (o)
compared to the
Another possible application of this code is determining the
maximum entrainment. The onset of entrainment is the
point where the first liquid drop becomes airborne. At
maximum entrainment there is no net increase in the liquid
content of the gas phase. At this point the interface is once
again stable with respect to long waves. The progression
from onset to maximum entrainment can be posed as a
stability problem applied to a mass balance. In other
words, given two flow rates the percent entrained
determines the gas density, the liquid height, and the
stability. Simply moving from zero to 100% entrained will
produce a line such as the one found in figure 10, which
determines the location of the onset and maximum
entrainment for those flow rates.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
x 10
.o0
W 2 %00
0
ct' 0
1 0 Unstable
0 Stable Below Onset
O Stable Max Entrainment
0 50 100 150 200 250 300 350
Re,
Figure 10: Entrainment line for a liquid mass flow rate of
0.53kg/min and a gas mass flow rate of 0.10kg/min. The
50 points span entrainment values from 0 to 99%
Summary
This paper briefly describes a new stability procedure for
turbulent gasliquid flows. A kepsilon model is shown to
provide good predictions of wavy stratified flows. The
onset of atomization is predicted by the instability of long
waves that occur on roughened twolayer basestate. The
criterion is not neutral stability but a finite growth rate
sufficient to allow such waves to grow in finite length
systems. The conditions of maximum entrainment are
predicted when the longwaves again become stable because
the liquid layer is getting thinner as the entrainment gets
larger. This boundary is calculated by reducing the liquid
flow rate and adjusting the gas flow density and viscosity to
account for the entrained drops.
Acknowledgements
The authors would like to thank the Chevron Energy
Technology Company for their funding.
References
Bruno, A Study of Interfacial Waves in GasLiquid Flows,
o University of Notre Dame (1988)
Bruno K., McCready M.J., "Origin of roll waves in
horizontal gasliquid flows" AIChE J. 34, 9, 14311440
(1988)
Chien K.Y, Predictions of Channel and BoundaryLayer
Flows with a LowReynoldsNumber Turbulence Model,
AIAA Journal, 20, 3338, (1982)
Cohen L.S., Hanratty T.J., "Effect of Waves at a GasLiquid
Interface on a Turbulent Air Flow" J. Fluid Mech. 31, 3,
467479 (1968)
Dutoya D., Michard P., "A Program for Calculating
Boundary Layers along Compressor and Turbine Blades",
Numerical methods in Heat Transfer, John Wiley & Sons,
New York (1981)
Hoffman G.H., "Improved Form of the LowReynolds
Number kepsilon Turbulence Model", Physics ofFluids, 18,
309312 (1975)
S20 40
0 '
0
5 "
. . . . . . . .
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Lam C.K.G., Bremhorst K.A., "Modified Form of the
kepsilonModel for Predicting Wall Turbulence", Journal
of Fluids Engineering, 103, 456460 (1981)
Launder B.E., Sharma B.I., "Application of the
EnergyDissipation Model of Turbulence to the Calculation
of Flow Near a Spinning Disc", Letters in Heat and Mass
Transfer, 1, 131138 (1974)
Mantilla I., "Mechanistic modeling of liquid entrainment in
gas in horizontal pipes". University of Tulsa (2008)
Reynolds W.C., "Computation of Turbulent Flows", Annual
Review ofFluid Mechanics, 8, 183208 (1976)
Smith E.W., Thesis, University of Notre Dame (2010)
Unpublished
Taitel Y., Dukler A.E., "A Model for Predicting Flow
Regime Transitions in Horizontal and Near Horizontal
GasLiquid Flow", AIChE J, 22, 1, 4755 (1976)
