7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Comparison between 2D CFD and 1D code for Wave Growth Simulations
S. Kalogerakos*, M. Gourma* and C. P. Thompson*
Department of Applied Mathematics and Computing, Cranfield University, Cranfield, Beds. MK43 OAL, UK
s.kalogerakos@cranfield.ac.uk
Keywords: 2D CFD, wave growth, multiphase flow, VOF
Abstract
Limitations in the predictive capability of onedimensional simulations of pipeline multiphase flow have led to
a recent increase in the use of threedimensional commercial CFD codes in the oil industry. On the other hand,
issues arise in 3D simulations when setting initial conditions and with simulation time lengths, especially when long
pipelines are involved. An alternative is to use a twodimensional CFD code with some modifications. In this paper
various methods are discussed that allow swift simulations of pipeline twophase flow with 2D code. Moreover
results from 2D commercial code simulations are compared with results obtained with inhouse developed 1D code
EMAPS (OmgbaEssama 2004). Comparisons include analysis of the wave growth problem, and measurements of
the wave rate of growth (Valluri et al 2008). Recommendations are given in cases where 2D CFD can be used with
no particular disadvantage compared with 3D CFD, and with comparable running times to 1D code. Specific cases
involving simulations where slug flow is observed (Ujang et al 2008) are also discussed.
Nomenclature
Roman symbols
g gravitational constant (ms 1)
j complex part
k wavelength (m)
p pressure (Nm 2)
u velocity (m/s)
Greek symbols
aL Liquid volume fraction
ac Gas volume fraction
S viscosity (kg/(sm))
p density (kg/(ms2))
w frequency (1/s))
Introduction
Commercial software Ansys FLUENT is a wellknown
CFD software solution that is being used by various
companies and in this paper use of FLUENT was aimed
at simulations of twophase flow in horizontal circular
pipes of constant diameter. On the other hand, the one
dimensional software called EMAPS (Eulerian Multi
phase Adaptive Pipeline Solver) was written in Cranfield
University and is based on Fortran F90.
Validation of FLUENT results was initially carried
out by investigating the wavegrowth problem, which
consists in introducing a perturbation at the inlet and
analysing its propagation and wave growth rate. In this
case, a full 8m long sine wave was introduced in a 30m
long horizontal pipe with constant diameter of 78mm,
at 15m from the inlet and constant superficial velocities
for liquid 0.4 m/s and gas 2.0m/s (water and air) were
imposed; these conditions were implemented by writing
custom user code (UDFs) in C++ for FLUENT. Initial
conditions are shown in Fig. 1.
Figure 1: Initial sine wave introduced in the 30m pipe.
I Xt '
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Simulations were carried out using both viscous
model kepsilon and Reynolds stress model. There was
an issue with the fact that a further perturbation appears
to stem from the inlet (in the FLUENT simulation) and
then propagates along the channel, quickly neutralising
the sine wave. In the EMAPS onedimensional code on
the other hand, the sinewave is propagating but keeps
its general sine shape well beyond 10 seconds.
It was thus decided to attempt to counteract the forma
tion of the perturbation at the inlet by first carrying out
a steady state simulation in FLUENT, and then applying
the sinewave to the steady state solution and starting
the transient simulation from that point. Gas was set to
be incompressible, in line with EMAPS settings. The
steadystate simulation converged, and the gas and liq
uid velocity profiles are given in Fig. 2.
0.99
0.94
S0.9
0 .84
079
074
Figure 3: Wavegrowth as
different cells.
10 15 20
Time(s)
calculated by
Figure 2: Velocity profile after steady state simulation
of wavegrowth.
A transient simulation was then started from the
steady state simulation result, and graphs of the wave
growth at different times are shown in Fig. 6 and Fig. 7
(contours of volume fraction). Using Matlab, it was pos
sible to estimate the rate of growth. The rate of growth in
EMAPS simulation was estimated to be 0.31 with resid
ual 0.013, while in FLUENT 2D simulation it was esti
mated to be 0.34 with residual 0.089. These values were
estimated from the slope of log(max liquid holdup) vs
log(time). As shown in Fig. 4, the time after which the
wave starts to grow is around Is and is different com
pared to EMAPS (around 11s, as shown in Fig. 3), but
the wavegrowth rate is similar as it only concerns the
rate of growth rather than its position.
It was observed that mesh refinement is essential:
with coarse grid, wavegrowth rate was 0.137 with 0.064
residual. In table 1 are the details of the meshes used.
A recent research carried out by Valluri et al (2008)
using research code TRIOMPH in Imperial College,
London concentrated on the wave growth problem us
Wavegrowth of sine wave, incompressible flow FLUENT 2D
0076
0074
0 072
0 07
0 068 +
0 066
0 064
0 062 +
006
0 2 3 4 5 6 7
Time (s)
Figure 4: Wavegrowth as simulated using incompress
ible FLUENT 2D
Table 1: Comparison of meshes used for wavegrowth in
FLUEND 2D
Max face area Cells Faces Nodes
Coarse Mesh 3.608e02 186,760 377,626 190,867
Refined Mesh 2.070e02 336,076 679,504 343,429
ing compressible flow, and thus it was decided to repeat
the previous simulation also using compressible flows
in FLUENT. The graph of maximum liquid heights vs.
time is shown in Fig. 5 and it starts occurring at around
3s, thus marginally closer to EMAPS results as com
pared with the incompressible flow results. Using Mat
lab, wavegrowth rate was found to be 0.32 with residual
0.062, thus again a better result. Therefore it appears
that using compressible flow in FLUENT 2D for wave
): 3
) '
EMAPS with
Wavegrowth of sine wave, compressible flow FLUENT 2D
0078
0076
0074
0072
007
0068
0066
0064
0062
0 5 10 15 20
Time (s)
Figure 5: Wavegrowth as simulated using compressible
FLUENT 2D
Figure 6: Contour of liquid volume fraction using Flu
ent 2D simulation of wavegrowth after 8.9s
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
growth analysis does indeed give results that are very
close to the ones obtained in 1D code, even though the
model is different.
Mathematical validation of rate of growth
In order to carry out a validation on the number obtained
in the previous section for the rate of growth, a pertur
bation analysis of the VOF model was completed. A
perturbation was introduced in the continuity and mo
mentum equations, and the dispersion frequency w was
eventually calculated, after solving the dispersion equa
tion D(k, w) = 0. Initially any second order terms were
neglected.
Continuity eqn. : dOtL + dOxLU 0
Momentum eqn. : Otpu + V(pd u + p)
a0 (0au) + y (, i) + pg+ FP
It was assumed that all timedependent functions can
be expressed as = , I '. Moreover (Fluent
2006) the source term F can be expressed as:
F= i (1)
(F (Pi + Pj)
where i and j are the two phases, Tij is the surface ten
sion coefficient, and K, is the curvature at the surface
where the surface tension is calculated. After various
simplifications, the real part and imaginary part of the
dispersion relation are shown below. It is assumed that
the pressure can be expressed as p = ( 1)cvTp = Fp,
where 7 is the ratio of specific heats and cv is the spe
cific heat coefficient at constant volume. Also Ap =
Pi Pj.
Real Part
D,(k, w) = k2ou(w 
w aoAp,9 + kuo
kpo
Imaginary Part
Dj(k,w)
uokaouAp + upo(w 
Spow 2 + w(uokaoAp
+ 2(2, ..,, , .,,, .., 
Figure 7: Contour of liquid volume fraction using Flu
ent 2D simulation of wavegrowth after 9.8s
kuo) kaouApg = 0
(2)
kaouAPkF
2 ..) kuo
 3pokuo)+
aoApF) 0
(3)
This is just a secondorder equation of the form ax2 +
bx + c, where the solution is given by bb2 ac
the previous equation, the terms are:
b2 +,qI12Ap2+,,.,9 O .12a O
"..,, a5Ap + ,,  : ",.,,A 0 an
4ac 8k2 2 2 ?11? ,,',A,) 4k2 opoApF
S\a() 2oa 1/2
b2 4ac kuopo (1 +c (Ap 2 2ao AP +4 4
A po P 0o 0 Apo }
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
When using the real part solution for the dispersion
relation as per eqn. 2, and substituting the variables used
for FLUENT simulation of wavegrowth, the number ob
tained is 0.31, which is remarkably close to the value
0.34/0.32 obtained with the simulation. This would in
dicate that as a magnitude check the procedure is consis
tent.
Slug formation and growth
After having shown that the wavegrowth problem can be
simulated in a satisfactory manner using 2D Fluent, thus
resulting in a simulation time that is approximately 10
times smaller compared with 3D Fluent, it was important
to investigate and compare results of simulations that
deal with twophase flow. In particular transient flows
with initial conditions that are known to produce slugs
of certain frequencies have been chosen for the current
investigation and comparison. Manolis cases (Manolis
1995) have been used in order to compare results with
known experimental data. The values of the initial con
ditions can be seen in table 2. It is important that the fre
quency obtained satisfies mesh independence, therefore
for each simulation three different meshes were used,
each with a higher number of cells as shown in table 3.
The mesh refers to a horizontal channel of length 30m
and diameter 78mm. The term cell configuration refers
to number of cells in the y direction and number of cells
in the x direction, so for example 203000 refers to 20
cells in the diameter direction and 3000 cells in the pipe
length direction. An example of a graph of liquid holdup
vs. time for case 36 is shown in Fig. 8.
Table 2: Manolis Cases used for hydrodynamic slug
flow simulation
Case VG (m/s) VL (m/s) RL
22 4.016 0.519 0.670
36 1.548 0.519 0.808
38 2.058 0.498 0.766
where VG is the gas superficial velocity, VL is the
liquid superficial velocity, and RL is the initial liquid
holdup.
As shown in Fig. 9 all three simulations of Mano
lis cases seem to converge to frequency results that
are within 4% of the experimental results. This agree
ment gives more support to the possibility of using two
dimensional Fluent in order to carry out pipe simulations
for twophases. Extrapolations to more phases and/or
more complicated shapes are not straightforward, and
more tests will have to be carried out in order to confirm
that. In more complicated shapes with less symmetry
compared with a pipe, the approximation of a pipe using
Manolis Case 36 Mesh609000
Liquid Holdup vs Time
Time (s)
Figure 8: Liquid holdup vs. time for Manolis case 36
Table 3: Slug frequencies in FLUENT 2D simulations
of Manolis cases 22, 36, 38
Case Cells 2D Fluent Experimental Manolis
Frequency (Hz) Frequency (Hz)
22 203000 0.0255 0.1333
22 406000 0.1286 0.1333
22 609000 0.1300 0.1333
36 203000 0.197 0.2444
36 406000 0.236 0.2444
36 609000 0.238 0.2444
38 203000 0.1467 0.2167
38 406000 0.2000 0.2167
38 609000 0.2250 0.2167
Mesh Convergence of
Slug Frequencies for Manolis cases
0.2500
0.2000 i"a' '
0.1500 Case
0.1000
ShCase 38
0.0500
0.0000
203000 406000 609000
Mesh used
Figure 9: Mesh convergence for slug frequencies calcu
lated for Manolis cases
a twodimensional channel will probably not be appro
priate.
7th International Conference on Multiphase Flow,
ICMF 2010, Tampa, FL, May 30 June 4, 2010
Conclusions
The use of 2D code gave satisfactory results that are
in agreement with EMAPS 1D code and with Manolis
data, and also mesh independence was reached. Vali
dation tests carried out include analysis of wave growth
and slug development. One of the advantages of this use
is knowing that the friction factors that are widely used
in many 1D codes are not satisfactory in slug regimes,
especially in the front region of a slug. In these regimes
it is difficult experimentally to estimate stresses at the
interface between the liquid film and gas (Hewitt 2009).
In fact in these regions slug superficial velocity is a
quantity that does not enter in any dynamics. As field
discontinuities are important, friction factors cannot be
approximated by continuous functions (Taitel & Duk
ler 1977). Simulation of transient flows using FLUENT
2D after an initial steadystate simulation, offers also the
possibility of calculating these stresses and these values
can be used to calibrate friction factors in 1D code, and
moreover simulation times are approximately 10 times
smaller compared with FLUENT 3D.
References
FLUENT 6.3 User Guide, Fluent Inc., 2006
Hewitt G. F., Private communication, TMF4 Meeting,
Imperial College of Science, Technology and Medicine,
UK, 2009
Manolis I.G., High Pressure GasLiquid Slug Flow,
Ph.D. Thesis, Department of Chemical Engineering
and Chemical Technology, Imperial College of Science,
Technology and Medicine, UK, 1995
OmgbaEssama C., Ph.D. Thesis, Cranfield University,
UK, 2004
Taitel Y. and Dukler A. E., A model for predicting slug
frequency during gasliquid flow in horizontal and near
horizontal pipes, Int. J. Multiphase Flow, 1977
Ujang P.M., Pan P., Manfield PD., Lawrence C.J., He
witt G.F., Prediction of the translational velocity of liq
uid slugs in gasliquid slug flow using computational
fluid dynamics. Multiphase Science and Technology,
Vol. 20. pp. 2579, 2008
Valluri P, Spelt PD.M., Lawrence C.J., et al. Numer
ical simulation of the onset of slug initiation in laminar
horizontal channel flow, INT J MULTIPHASE FLOW,
2008, Vol: 34, Pages: 206 225, ISSN: 03019322
