7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
A Coupled Finite Volume Solver for the Simulation of Disperse Multiphase Flows
Marwan Darwish, Amer Abdel Aziz and Fadl Moukalled
American University of Beirut, Faculty of Engineering and Architecture, Mechanical Engineering
Riad El Solh Street, Beirut, 1107 2020, Lebanon
darwish@aub.edu.Ib
Keywords: Multiphase Flow, Coupled Algorithm, Disperse Flow, VelocityPressure Algorithm, EulerianEulerian, Finite
Volume Method
Abstract
In this paper a coupled multiphase finite volume collocated variable arrangements pressure based solver developed within a
EulerianEulerian framework is presented. The coupled solver differs from pressure based segregated solvers in that it
accounts implicitly for the PressureVelocity and the interphase drag couplings that are present in disperse multiphase flows,
thus yielding a system of coupled equations with velocity fields and pressure. Two test problems for dense and dilute
disperse flows are show to yield substantial decrease in the number of iterations to convergence as compared to a pressure
based segregated multiphase solver.
Introduction
At the heart of Computational Fluid Dynamics (CFD) is
the velocitypressure coupling algorithm that drives the
convergence of fluidflow simulations. Over the last
decades efforts to devise velocitypressure algorithms that
are more robust, and more efficient have resulted in a
better understanding of the numerical issues affecting the
performance of these algorithms, such as the choice of
primitive variables (densitybased see Chorin (1967) and
Rizzi et al. (1995) versus pressurebased see Moukalled
(2000), Lien et al. (1994)), the type of variable
arrangement (staggered versus collocated arrangement),
and the chosen solution approach (coupled versus
segregated), to cite a few of the important issues. These
issues and their effects are now much better understood by
the research community. Recently there has been a
renewed interest into the development of coupled solvers,
see Webster (1996), Darwish (2009) partly because of the
tremendous increase in available memory in computers in
addition to the convergence problems experienced by
segregated solvers (see van Doormal et al. (1985)) when
used with large computational domains. This last problem
has been addressed successfully through parallel
processing and domain decomposition, but nevertheless the
convergence issue has not by itself been resolved. In it
worthwhile in this respect to point out that density based
Euler solvers have been using quite successfully the
coupled approach. In the single fluid flow coupled
approach the NavierStokes equations and sometimes
additional equations are discretized and solved as one
system of equations as opposed to the segregated approach
where the equation for each variable is solved separately
using a fixed, best estimate values of the other dependent
variables. This approach has been shown by the authors to
yield drastically improved convergence rates and
robustness see Darwish (2007) and Darwish (2009).
When dealing with multiphase flows the complexity of
the velocitypressure coupling is compounded by the need
to resolve the coupling to pressure of multiple velocity
fields, in addition to dealing with the intervelocity
coupling that arise from the interfluid drag terms. This
dual layer of nonlinearity has resulted in large
computational time needed to solve even smallscale
multiphase flow problems.
In this paper the coupled singlefluid flow algorithm of
Darwish (2009) is extended to multiphase flow with the
fluidic momentum equations and the pressure equation
being assemble into one matrix with all the relevant
coupling between the various velocity and pressure fields
resolved in an implicit manner. The tight velocitypressure
coupling is developed within the context of a collocated
unstructured grid, and the system of equations involving
velocity and pressure is solved simultaneously.
In what follows the governing equations for multiphase
flow are introduced and presented in their semidiscretized
form. The coupled algorithm is then presented and its
performance compared to the segregated multiphase flow
finite volume solver (see Moukalled (2004a) adn
Moukalled (2i"li4b) in terms of iterations to convergence
and computational time, for two particles laden flow
problems.
Nomenclature
ap, anb coefficients of the algebraic equation
g gravitational constant (in i,
drag coefficient (
P pressure (Nm2)
v velocity
B body force
D RhieC lio\ Coefficient
I Interphase drag term
S surface vector
Greek letters
a volume fraction
p viscosity (Pas)
p density
Q volume of element
Subscripts
P main control volume or element
F neighboring control volume or element
f refers to face f of the element
n refers to neighboring element
Superscript
(k) refers to phase k
(d) refers to disperse phase
(c) refers to continuous phase
(slip) refers to slip velocity
Governing Equations
In Eulerian multiphase models the local conservation
equations for mass and momentum are can be obtained
through a volume averaging procedure to yield for steady
state incompressible flow (see Drew (1983))
V. c((ck) (k)V ))
o (k) (k)
V (k)' 'k)V(k)) (k) ( P+B(k))
(k) g(mik) (v(m)
mk
+V (a(k) (k))+ I(k)
 V (k)
(k) = )(Vv(k) + V() ) (4)
where for the kth phase, a(k) is the volume fraction, p(k) is
the phasic density, vk) the phasic velocity, P stands for the
pressure, which is regarded as being shared amongst the
fluids, B(k) is the body force per unit volume of phase (k)
and Ik) is the momentum transfer to phase (k) resulting
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
from interaction with other phases.
The volume fraction is the phasic equations are governed
by the a geometric constraint written as
Sa (k)=l (5)
k phases
For the computation of the pressure field a pressure
equation is derived from the global continuity equation.
The global continuity is obtained from the sum of the
phasic continuity equations to yield
SV. (a(k)p(k(k)) = a(k)(k)= 0 (6)
k phases k phases
This equation will be used in deriving the pressure
equation.
Discretization
When discretizing the above equations using a segregated
algorithm the pressure gradient in the momentum
equations and the velocity fields in the global continuity
equation are treated explicitly. In the coupled algorithm
this deficiency is alleviated by treating both the velocity
and pressure terms implicitly, coupling the momentum
equations and the pressure equations through a set of
coefficients that represent the mutual influence of the
continuity and momentum equations on the pressure and
velocity fields. To that end the pressure gradient term in
the momentum equations is integrated over the surfaces of
the control volume and evaluated implicitly by expressing
the pressure at the face in term of the nodal values
straddling the interface. Considering the flow of a disperse
phase (d) in a continuity phase (c) as in the transport of
particles in air we get after integrating the momentum
equation
X(a (c) P (CvC) (c) CVVC ) Sf
+a(c)PS yp + Q, ( v (v () v) (7)
f m k
= a (c)QBc)c
further discretization yields an equation set of the form
ap()u(,)+ av(c)v(,)+ aPu(d) (d) a+p(c) +
F=NB(P) F NB(P)
+ I a;p P b'
Pc)F = bFb
FNB(P)
avy(c)v()+ a(c)U(c)+ (d) (d) ap() +
(c) (c) vu(c) (e)
FNB(P) FNB(P)
vp(c) p=b
S + I aP)F F
F NB(P)
In a segregated algorithm the underlined terms are treated
explicitly i.e. using values from the previous iteration and
added to the RHS of the equation. In the coupled approach
these terms are written as a function of the nodal value and
their influence on the other variables accounted for. This
include the pressure and drag terms that generate a velocity
pressure influence and an interphase velocityvelocity
influence. When treating certain boundary conditions
coupling between the various components of velocity is
also accounted for. A similar set of equations is also
derived for the disperse phase.
Substituting the RhieC 1ho" interpolation of velocity into
the global continuity and then expressing the averaged
velocity in term of the straddling nodal velocities, we get a
pressure equation in which the pressure is tightly coupled
to the phasic velocity fields
Pf (DfVp,) S
a s f + s 
k~phases f=nb(P) If Vf f (9
f =nb(P) (9)
= k
k~ phases
{PPf(W f)* Sf}
Expanding on the discretization yields for each element a
system of equation of the form
aP ap ap"
a a 0 a
a" 0 a" a"'
ap ap ap
PP0 0 a
a aB (
o o a"
a,, af
U P k
a p
Sa v
P U
Y C akp Vp
Sa X u +
ap
PP
0 a v b
rP' PP P
"a x = b
a' a v b
a a P< b_
These equations are assembled into a global matrix of a
form
[] []
L 1[ 6 1 [ 6
where each [] is a 5x5 matrix of coefficients and { } is a
5x1 array of the variables or sources.
The RhieChow Interpolation
The computation of the RhieC ho"\ coefficient for single
fluid flows is relatively simple (see ), starting with the
equation for the xmomentum equation at element P
written as
a up = an Q u + nb (12)
N=NB(P)x
for element F a neighboring element to P we write
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
aFUF = al F p  + b (13)
N=NB(F) F. F
Interpolating the two equation to the common face f we get
afuf u f I; +bf (14)
n=nb(f) f
yielding the RhieC 1ho" interpolation
(P
Uf =Uf d (15)
I
where
df Qf
af)
(16)
Uf n  +bf
n=nb(f) ) if
and the average value is given by
f1 = s + sfI (17)
where f denote the face common to elements P and F, f is
the a geometric interpolation factor, Q( is the volume of
the element and ap is the coefficient of the momentum
equation.
For a two phase flow simulation with interface drag the x
momentum equations for the velocities ul and u2 of the
two fluids at element "P" can be written as
'l1(1) P 1'
U 2 (1)
a211 j(2
N a j {(1
\a(1) 1 ( k)L b pl)
a(2) x, b1 2)
Sp P
A similar equation can be written for element "F"
a" l a "(1)
u2U1 2 2 (2),
a u
2_ a) 0 ]U uIN (19)
N NB(F) 0 a"2M2 (2)
Lit
a ( P\ b (2
Interpolating the two equation to the common face f yields
u u ,
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
11 12 i I)
a ll a 12 i u2
u2u1 2a 2 1)
a a u
where now
.2 2
a
uu u2U2 ulu2 U2U1
a a a a
a
ulu1 U2U2 aulu2 aU2U1
a a a a
.1U2
a
a a a a
aul 1
ulu1 U2U2 aUlU2 aU2U1
a a a a
(1) aP ble
_[ a'] N" (r [ ^1
(2) ax ) b'(2)
thus writing the RhieC lio'\ interpolation as
r(1) f ] r(1) 1
u u 2
I u m if
d"" d""' ao> aP
d"2"1 d"22 a(2) x f
I f
of, "2 p,
(22)
In the limit of no interphase drag with
a"u"2 = a21 = 0
we recover the single phase formulation
1
0
1 o1
Df= a f (23)
1
0
a
f
In disperse flow problems where the momentum
coefficients are more important the simplified RhieC Iho"
relation was used.
Slip Walls
.
\ P=0
Figure : SchematicL of hori l tt
Figure 1: Schematic of horizontal test cases
Results and Discussion
The performance of the coupled algorithm is assessed by
solving four simple cases involving the transport of dense
and dilute solid particles in air for a horizontal and vertical
channel. In the horizontal cases, see figure 1, the inlet
velocity for air and particles is 5 m/s and 1 m/s
respectively, only the drag force between the two phases is
considered as the gravitational force is set to zero. In the
vertical test cases, shown in figure 2, the inlet velocities are
100 m/s and 10 m/s for the air and particles respectively, in
addition to drag, gravity effects are accounted for by setting
g to 9.81 m/s2. The drag coefficient between the two
phases is given, for both test problems, by:
(km) 3 CD a(d)p(c)(sp) (24)
4 d '(
CD is set to 0.44 in all cases. The subscripts c and d in the
above relation denote the continuous (i.e. air) and disperse
(i.e. particle) phases respectively. The parameter d(d) is the
particle diameter.
V(shp) is the slip velocity and is given by:
V(sUip)= V (d) v(c) (25)
No diffusion term appears in the disperse phase momentum
equation.
v`) ( .'
Sg=10 m/s2
SSlip Walls
V(d) ,(d)
Figure 2: Schematics of vertical test cases
For the horizontal cases, analytical solution of the velocity
of the disperse phase, u(d), is given as a function of x, the
distance from the inlet, and the two phases inlet velocities
and properties as
i
3 p(c) CD
4 P X
4 p(d) d(d)
/ )V(C)
+lnVc) + nmlet
Sinlet ) (c) V(d)
mlet inlet
Three meshes where used in all test cases in order to
evaluate the performance of the coupled and segregated
solvers with increasing mesh density. A Underrelaxation
factor of 0.7, 0.7 and 0.3 were used with the segregated
solver for the volume fractions equations, the momentum
equations and the pressure equation respectively. For the
coupled solver a false transient step varying between 0.006
was used for the cases with no gravity and 0.007 for the
cases with gravity.
Horizontal dilute solid particle transport in air
For this case the boundary and initial values are set as
follows: at inlet are V(C) = 5m/s, V(d) = Im/s, d) = 105 and
ac) = 1ad). The other properties are d(d) = 2mm, p(c) = 1
kg/m3 and pc) = 2000 kg/m3. The kinematic viscosity of air
is p=105 Pa.s. Since the solid concentration is dilute, the
equilibrium velocity is nearly equal to that of the
continuous phase. Figure 3 shows the two phases velocities
as function of x. The numerical solution by coupled and
segregated algorithms falls on top of the analytical solution
which is obtained from equation (26). The settling velocity
in this case is 4.99996m/s. Table 1 shows the iteration
number to convergence and the time to convergence for the
segregated and coupled solver for three meshes of size
10,000, 30,000 and 50,000 elements. A false transient step
was used for the coupled solver and under relaxation factor
for the segregated solver
Mesh Cost Coupled Segregated
(xl03) Reduction Iter# time Iter# time
10 45.8% 51 355 501 654
30 21.3% 138 3,061 959 3,890
50 27.1% 197 7,188 1466 9,856
Table 1 Dilute horizontal AirSolid Particle Transport
4SI
II I l l i l l m i l l .
AN Cpit"G
>cMCd
*CM4PII)
;10 4 0
X
s t0
Figure 3: Dense air/solid horizontal particle flow
InV (C) U (d) V nlet
inlet vc) l(d)
inlet
S6
S
I/ ApMP
2 / 8s anmPtU
/  nft
5
.10 .5
a tI
Figure 4: Dense air/solid horizontal particle flow
Table 2 shows the iteration # and computational time
obtained using the coupled and segregated solvers
Mesh Cost Coupled Segregated
(xl03) reduction Iter# time Iter# time
10 43.4% 54 373 500 659
30 22.7% 136 2,982 958 3,860
50 26.8% 195 7,159 1465 9,783
Table 2: Dense horizontal AirSolid Particle Transport
Vertical dilute and dense solid particle transport
in air
For this cases the flow parameters are set as follows: at
inlet are Vc) = 100 m/s, V(d) = 10 m/s, Xd) = 106 and ac) =
1ad). The other properties are d(d) = 2mm, p(C) = 1 kg/m3
and pC) = 1000 kg/m3. The kinematic viscosity of air is
p=105 Pa.s. The large inlet value is to ensure that the solid
does not exit through the inlet due to gravitation. For the
dense case (d) is set to 102. Most of the inertia in the dense
case is carried by the disperse phase.
Figure 5 and 6 show the numerical solution for the coupled
and segregated algorithms. For the dense particle, the air
velocity dips slightly as expected due to the momentum
transfer to the disperse phase, in the dilute case the dip is
imperceptible. The increase in the velocity of the disperse
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Horizontal Dense Solid Particle flow in air
In this case the boundary and initial values are similar to
the previous case, except that, a(d) is set to 102. The
analytical solution of this case can be computed from
equation (26). The parameters in this case are just as those
in the previous one except that a(d)=102. Although the
disperse volume faction is low, the flow problem is termed
dense since the disperse phase carried most of the inertia.
This is because the mass fraction of the solid disperse
phase is 20 time larger than that of the continuous phase:
PC,,, 1 x 0.99 1
= (27)
Pdadn 2000x0.01 20
The equilibrium velocity in this case is 4.96m/s. figure 4
shows the velocity profiles of the two phases.
r  .. 
momentum transfer from the
phase is due to drag
continuous phase.
100 
90 .
to o
S3
SIQto
Figure 5: Dilute air/solid vertical particle flow
Tables 3 and 4 show the number of iterations and
computational cost to convergence for the dilute and dense
cases, respectively. The cost reduction is around 40% for
the 30k and 50 k meshes, and slightly higher at 50% for the
coarser 10k mesh. This better than for the horizontal
transport cases, this can be explained by the fact that
adding effects of gravity increases the cost per iteration by
a nearly constant amount for discretization and because of
the added nonlinearity. This year more advantage to the
coupled algorithm that requires a relatively low number of
iterations to convergence.
Mesh Cost Coupled Segregated
(xl03) reduction Iter# Solving Iter# Solving
time time
10 51.9% 41 283 465 588
30 39.3% 99 2,182 906 3,595
50 41.8% 144 5,357 1397 9,198
Table 3: Dilute vertical AirSolid Particle Transport
lap
aqc.~sua
u ImqcmUC.
 Da4ILE)
Y UdISU)
0 4 0
X
S 10
Figure 6: Dense air/solid vertical particle flow
Se(cedrupen
a2d (CCbdl
Ab(SIPLI)
 o (iPd( LIt
$66d (SOPLEP
  _
0 I US l *
 0
x
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Mesh Cost Coupled Segregated
Mesh Cost P
(xl03) Reduction Iter# Solving Iter# Solving
time time
10 53.3% 41 286 468 613
30 41.3% 99 2,167 912 3,691
50 43.7% 145 5,330 1396 9,462
Table 4: Dense vertical AirSolid Particle Transport
Conclusions
A coupled velocitypressure algorithm was extended to
twophase EulerianEulerian flows. The coupled algorithm
shows promising performance improvement as compared
to the segregated algorithm. In general performance better
that segregated algorithms with computational cost
reduction ranging from 50% down to about 20%.
Further refinement to the coupled algorithm are currently
under review.
Acknowledgements
This work was supported by the the University Research
Board of the American University of Beirut.
References
Chorin AJ., A Numerical Method for Solving
Incompressible Viscous Flow Problems, J. Comput. Phys,
Vol. 2, 1226 (1967)
Rizzi A., Eriksson LE., Computation of Inviscid
Incompressible Flow with Rotation, J. Fluid Mech., Vol.
153, 275312 (1985)
Drew D.A. Mathematical modeling of twophase flow,
Ann. Rev. Fluid Mechanics, Vol.15, 261291 (1983)
Darwish M., Sraj I., Moukalled F. A Coupled
Incompressible Flow Solver for the solution of
incompressible flows on unstructured grids, Num. Heat
Trans. part B: Fundamentals, Vol. 52, 353371 (2007)
Darwish M., Sraj I., Moukalled F. A Coupled
Incompressible Flow Solver for the solution of
incompressible flows on unstructured grids, J. Comput.
Physics, Vol. 228, 180201 (2009)
Lien FS., Leschziner MA. A General NonOrthogonal
Collocated Finite Volume Algorithm for Turbulent Flow at
All Speeds Incorporating SecondMoment Turbulence
Transport Closure, Part 1: Computational Implementation,
Comput. Methods Appl. Mech. Engrg., Vol. 114, 123148
(1994)
Moukalled F. & Darwish M., A Unified Formulation of the
Segregated Class of Algorithms for Fluid Flow at All
Speeds, Numerical Heat Transfer, Part B: Fundamentals,
Vol. 37, 103139 (2000)
Moukalled F, Darwish M., PressureBased Algorithms for
MultiFluid Flow At All Speeds Part I: Mass Conservation
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Formulation Num. Heat Trans., part B: Fundamentals, Vol.
45, 495522 (',21"14,
Moukalled F., Darwish M., PressureBased Algorithms for
MultiFluid Flow At All Speeds Part II: Geometric
Conservation Formulation Num. Heat Trans., part B:
Fundamentals, Vol. 45, 523540 (2 1141b)
Van Doormal J.P., Raithby G.D., An evaluation of the
segregated approach for predicting Incompressible fluid
flows, ASME paper 85HT9, National Heat Transfer Conf,
Dever, Colorade, (1985)
Webster R., An Algebraic Multigrid Solver for Navier
Stokes Problems in the Discrete SecondOrder
Approximation", Int. J. Num.Methods Fluids, Vol. 22,
11031123 (1996)
