7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Fourway coupling of dense particle beds of black powder in turbulent pipe flows
C. Narayanan, M. Labois, D. Lakehal
ASCOMP GmbH
Technoparkstr. 1, Zurich, CH8005 Switzerland
Email: lakehal@ascomp.ch
Keywords: Twophase flow, EulerianLagrangian method, black powder, fourway coupling, VLES, LES
Abst ract
The accumulation of black powder in pipelines may be very harmful for hydrocarbon transportation installation. Recent
engineering studies are heavily investing in computerbased predictive strategies for this phenomenon to anticipate
hypothetical blackpowder slug formation and develop fast and efficient removal techniques to operate in time. But the
modelling of this phenomenon is one of the most challenging one in multiphase flow, because the physics is multifaceted and
rather complex: turbulence of the carrier phase; particleturbulence interaction, particlewall interaction, particleparticle
interaction, two and four way coupling, particle agglomeration, deposition and resuspension. We will discuss these issues and
present new routes for particle collision stress modelling. Practical examples including sand and oil slugs' formation up slope
pipelines will be presented and discussed. The model employed is based on denseparticle formulation accounting for
particleturbulence interaction, particlewall interaction, particleparticle interaction via a collision stress. The model solves the
governing equations of the fluid phase using a continuum model and those of the particle phase using a Lagrangian model.
Interparticle interactions for dense particle flows with high volume fractions (> 5%) have been accounted for by mapping
particle properties to an Eulerian grid and then mapping back computed stress tensors to particle positions. Turbulence within
the continuum gas field was simulated using the VLES (Very LargeEddy Simulation) and full LES, which provides sufficient
flow unsteadiness needed to lift up the particles and move the deposited bed. The main issues and limitations will be discussed
in the paper.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
This work aims at studying the mechanisms of turbulent
dispersion, deposition and transport of solid particle in
dense packed beds formed in 3D pipes. For the purpose,
use has been made of the CMFD solver TransAT (2010), a
multiphaseflow dedicated computational fluid dynamics
code. The physical model is introduced below. The
mechanism of particle dispersion requires the turbulence
in the carrier fluid phase to be modelled with a more
sophisticated approach than RANS. In TransAT, we
promote the use the VLES and LES approaches,
depending on the Reynolds number.
The Physical model in TransAT
2.1. EulerianLagrangian for dense particle systems
The EulerianLagrangian formulation for dense particle
flows featuring mildtohigh volume fractions (a > 5%) in
incompressible flow conditions is implemented in the
CMFD code TransAT as follows (Eulerian mass and
momentum conservation equations for the fluid phase and
Lagrangian particle equation of motion):
a, (ap) + 8, ap, )= (1)
(2)
t, [x,,(r)] )+ g + F,,it
d,(v, ) = f 9ft (st,,
2p,,dy,
fd=1+0.15Re;
where a is the incell volume of fluid (a+a,= 1), u is the
velocity of the carrier phase, u, is the velocity of the
carrier phase at the particle location, v, is the particle
velocity, II is the sum of viscous stress o and pressure p, I
is the turbulent stress tensor (depending whether RANS ,
VLES or LES is employed). Sources terms in (2) denote
body forces, Fb, and the rate of momentum exchange per
volume between the fluid and particle phases, F,. In (3),
F,,on denotes the interparticle stress force. The momentum
equation (2) presented here does not neglect viscous and
turbulent diffusion mechanisms in the fluid phase. The
interphase drag model in (3) is set according to Gidaspow
(1988).
2.2. Fourway coupling modeling
The particle volume fraction is defined from the particle
distribution function (4) as
ag,= V,, ~dV,, dp,i du, (4)
The interphase momentum transfer function per volume
in the fluid momentum equation (2) is
FP JTVP [A]d'Vp d/p, du : (5)
with A standing for the particle acceleration due to
Introduction
Accurate prediction of particle dispersion in the
atmosphere and in various engineering applications is
very important in the assessment of human safety and in
the analysis of risk. One of the engineering segments
which is very much interested in acquiring accurate
prediction tools for particle dispersion is the Oil and Gas
industry, for instance in relation to Black Powder
deposition and transport in pressurized gas pipelines.
Black powder (Smart, 2007) is toxic, with a complex
constitution, including radioactive elements. It consists of
iron compounds such as magnetite and iron sulfide and
includes sand and clay, salt, weld spatter and metallic iron.
It is generated during gas production or in wet gas
pipelines when hydrogen sulfide, carbon dioxide or
oxygen is present in the gas, by bacterial corrosion of the
steel, or from construction when lines are not cleaned
adequately. It influences the flow performance of gas
pipelines, impairs the function of valves and metering
systems, and leads to severe accidents during transport.
In denseparticle beds systems, the flow behaves in a very
subtle way, with very complex physical mechanisms near
the wall, where the powder accumulates. A number of
simplified analytical solutions to determine the conditions
of particlebed removal in pipes and channels have indeed
been proposed, but with limited success due the
simplifications implied in these models. Today intensive
research is devoted to understand the conditions for dense
particlebed formation and removal, in hydrocarbon and
in many other related areas, like chemical engineering,
but the difficulties encountered in measurements and flow
visualization have hindered this progress. There are
various incentives to explore the use of advanced
prediction methods for this class of flows, featuring
Lagrangian particle tracking including fourway coupling,
instead of average EulerEuler formulations, LargeEddy
Simulation (LES) instead of RANS, and transient rather
than steadystate simulations.
In particleladen pipelines, the particles tend to be
transported through the pipeline by gas flow under
specific conditions. The velocity required to move the
particles could in some cases be estimated, made based on
pipeline diameter, gas pressure, and particle size and
density (Tsochatzidis & Maroulis, 2007; Smart, 2007;
Davies, 1987). When black powder moves, it shatters and
becomes very small in size, in the range of one micron or
less, making it difficult to filter and possibly easier to
move. Deposition of black powder will occur if there are
solids in the pipeline fluid and the velocity is not high
enough to drag the particles along by viscous flow forces.
Sediment deposits can lead to blockage of the line,
especially during piggino, while flowing powder can
damage compressors, plug filters and damage user
equipment. In some extreme cases, the piping could be
half full of black powder, causing shutdown of the
compressor and up to 60 tons of black powder could
subsequently be removed from the piping.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Immersed Surfaces Technology (IST) developed by
ASCOMP GmbH and implemented in TransAT (Labois et
al., 2010). The idea is inspired from ITM for twophase
flows, in that the solid is described as the second 'phase'
with its own thermomechanical properties. The technique
has the major advantage to solve conjugate heat transfer
problems. The solid is first immersed into a cubical grid
covered by a Cartesian mesh. The solid is defined by its
external boundaries using the solid level set function. Like
in fluidfluid flows, this function represents a distance to
the wall surface; is zero at the surface, negative in the
fluid and positive in the solid. The treatment of viscous
shear at the solid surfaces is handled very much the same
way as in all CFD codes. To better resolve boundary
layers, IST is complemented by the BMR (block mesh
refinement) technique, where additional refined
subblocks are automatically generated around solid
surfaces; with dimensions made dependent on the
Reynolds number (based on the boundary layer thickness)
and desired y+ for wall treatment (low Re model or wall
functions). This combined method can save up to 75%
grid cells in 3D, since it prevents clustering grids where
unnecessary.
VLES versus LES
The idea of VLES (i.e. VeryLarge Eddy Simulation; see
Speziale, 1999) is to combine RANS and LES for a
specific flow, where the size of the most important scales
can be identified. Here the flow is decomposed into
resolved and subscale part, the latter being independent
from the grid (in contrast to the subgrid scale modelling
in LES), but is dependent on the flow, and thus the flow
characteristic length scale. Larger scales than this specific
scale are directly resolved, meaning actually no model is
included, but the subscale part is modelled, though with a
more refined statistical turbulence model than the
Smagorinsky one, because the subscales of turbulence
are not necessarily isotropic or independent of the
boundary conditions, as speculated in LES. Typically,
twoequation or Reynolds Stress models can be applied.
The approach assumes that the Kolmogorov equilibrium
spectrum applies to the subfilter flow portion, which
justifies the use of RANS models. The VLES used in this
study uses the k Emodel as a subfilter model. The
filter width is no longer related to the grid size, but is
made related to a characteristics lengthscale of the flow.
Increasing the filter width beyond the largest length scales
will lead to predictions similar to standard RANS models,
whereas in the limit of a small filterwidth (approaching
the grid size) the model predictions should tend towards
those of LES. If the filter width is smaller than the length
scale of turbulence provided by the RANS model, then
larger turbulent flow structures will be able to develop
during the simulation, depending on the simulation
parameters (e.g. grid, time stepping and order and
accuracy of the schemes employed). The VLES theory as
currently used was proposed by Johansen et al. (2004).
aerodynamic drag (1st term in the RHS of Eq. 3), i.e.
excluding body forces and interparticle stress forces (2nd
and 3rd terms, respectively). The pressure gradient
induced force perceived by the solids is not accounted for.
The fluidindependent force Fooll is made dependent on the
gradient of the socalled interparticle stress, n, using
Feat = VE/ #;,a, (6)
Collisions between particles are estimated by the isotropic
part of the interparticle stress (its offdiagonal elements
are neglected.) In most of the models available in the
literature n is modelled as a continuum stress (Harris &
Crighton, 1999), viz.
r;rb'=' (7)
max [acy, a,; E(1 a,)]
In TransAT, the particle field is predicted in a Lagrangian
way first, which enables defining the particle volume
fraction (4) and interphase momentum transfer function
(5), then highorder accurate interpolations are resorted to
map the Eulerian field (to estimate n), then back again to
the Lagrangian system to determine Fooll). The constant Ps
has units of pressure, and O', is the particle volume
fraction at close packing, and constant P is set according
to Auzerais et al. (1988). The original expression by
Harris & Crighton (1994) was modified to remove the
singularity at close pack by adding the expression in the
denominator (Snider, 2001). The E is a small number on
the order of 10 7. The particle stress is unaffected by the
modification except when the volume fraction approaches
or exceeds close pack limit, which is somewhat arbitrary
and depends on the size, shape, and ordering of the
particles. This limit can actually be physically reached or
even slightly exceeded.
The Numerical Approach in TransAT
The CMFD code TransAT@ developed at ASCOMP is a
multiphysics, finitevolume code based on solving
multifluid NavierStokes equations. The code uses
structured multiblock meshes. OpenMP and MPI parallel
based algorithm can be used. The grid arrangement is
collocated and can thus handle more easily curvilinear
skewed grids. The solver is pressure based (Projection
Type), corrected using the KarkiPatankar technique for
lowMach number compressible flows. Highorder time
marching and convection schemes can be employed; up to
third order Monotone schemes in space and 5th Order
RungeKutta time marching schemes. Multiphase flows
can be tackled using interface tracking for both laminar
and turbulent flows (Level Set, VOF with interface
reconstruction, and Phase Field), the phase averaged
homogeneous mixture model (with Algebraic Slip), and
the Lagrangian particle tracking (one, two, and fourway
coupling, including with heat transfer).
To mesh complex geometries, use is made of the
VLES of the flow and particles transport in pipe
5.1. Problem setup and modeling
Particle transport from bed
TransAT (ASCOMP. Zurkh)
Figure 1. Initialization of the particle bed; the particles are
colored by their size
The objective of this section is to estimate/predict the
critical velocity of dust transport in gas pipelines based on
calculations. From the modelling point of view, the goal is
to come up with the correct modelling requirements
which allow the estimation of the critical velocity. The
critical velocity is defined as the gas velocity needed to
transport 10% of the initially injected dust mass from the
starting point to the filter separator. Gas composition, gas
temperature and dust characteristics being kept constant.
Simulations were performed for one experimental
condition (Narayanan & Janssen, 2009) of gas flow at
system pressure of 10 bar and Gas temperature of 200C.
The simulations involved changing the gas superficial
velocity. The other flow parameters considered are:
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
transported further in the pipe with increasing gas velocity,
still without reentrainment. (II.1InI'_II' the VLES model
parameterization has not helped liftup the particles.
Length of dust layer: Im.
Dust layer thickness: 12mm.
Dust density: 2650kg/m .
Dust size: 200 400pm.
Initial dust mass: 20gm.
The gas flow field was solved using a VLES approach
with unsteady inflow conditions, as introduced previously.
The dust mass was modelling using Lagrangian particle
tracking (Eqs. 13). The full (4way) mass and momentum
coupling between the dust and gas phases was accounted
for based on the model of Snider (2001); Eqs. (47). In
this model the interparticle collision is not directly
simulated, instead a collisionpressure stress is introduced
that indirectly takes into account the closepacking of
particles. The dust bed was initialized as a smooth bump
as specified (see Fig. 1). Particles of 4 different diameters
were used to represent the particle size distribution: 225,
275, 325, 375pm. Inflow gas velocities of 2, 3, 5, 10 m/s
were used. Each simulation was run for 10'000 iterations,
with different time steps to adhere to the CFL criterion,
for real times of 25, 16, 9, and 4.2 seconds, respectively.
5.2. Discussion of the results
Figure 2 shows that particles are indeed moved by the
flow along the pipe. However, all the particle transport
takes place along the wall of the pipe. Very little particle
entrainment into the core flow is observed. Particles are
Figure 2. Snapshots of the simulations at 4 different times for
4 different gas velocities.
The critical velocity of transport was estimated by
extrapolating the rate of change of particles in the domain
to 3min, from the existing simulation data. This is
presented in Fig. 3. It can be inferred from this graph that
the critical velocity (defined as 10% particle transport in
3min) is around 3m/s or slightly lower. In the graph, the
dashed part of the lines is linearly extrapolated results
from the simulation. It was assumed that the simulations
were run long enough to reach a quasisteady state of
particle removal. Even with the particles not being
liftedup in the core flow, the predicted value of the
critical velocity matches exactly the measured value.
The size distribution of particles at the exit of pipe for
critical transport conditions is shown in Figs. 4 & 5. All
the values are taken at the end of the simulation period,
which is different for different inlet velocities.
50000
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
LES of the flow and particles transport in channel
Figure 5 in particular shows that at low velocities, a larger
fraction of smaller particles are transported, whereas at
higher velocities an equal fraction of all particle sizes is
transported. A larger difference in relative transport would
be evident if a larger range in particle sizes is considered.
In the simulated cases, the 4 diameter values considered
are not that different (225375pm). The particles size
distribution was assumed to be Gaussian with a mean
diameter of 300pm and a standard deviation of 50pm. The
particles were randomly initialised in the initial particle
bed volume as shown in Fig. 1. In summary, it seems that
the rate of particlestripping from the bed is actually well
simulated together with the critical transport velocity,
even if the particle cloud has not lifted up.
Figure 6. Particle transport in channel (4way coupling)
We present here first attempts to resolve the flow using
LES of the gas flow in the pipe, instead of VLES.
PUTSuing the flow analysis in this direction is motivated
by the perception that smallscale turbulence (modelled
by VLES) may be responsible for liftingup the particle
ClOud. The Lagrangian particle module is kept the same as
explained previously; 4way coupling. The subgrid scale
(SGS) model for fluid flow turbulence was treated with
the Nicoud & Ducros (1999) model. No SGS model has
been used for particle motion equation. A channel flow
was simulated, with initial and fluid flow conditions
exactly as in Narayanan et al. (2003). Periodic boundary
conditions were set. The initialization of the particle bed
was done as in the VLES case.
Snapshots of the flow obtained with the model are
presented in Fig. 6 below, depicting particle concentration
in the bed (coloured by their size) and the flow developed
through the interaction with the carrier phase. Turbulence
structures generated subsequent to the interaction are
better highlighted with the vertical velocity contours.
Again, it seems that even with full LES, the particle
reentrainment in the core flow is not predicted.
Conclusions
The paper presents a simulation campaign of a turbulent
gas pipeflow laden with solid particles of different size,
under different system pressure levels. The rate of
particlestripping from the bed is well simulated, and the
prediction of the critical transport velocity is in line with
the experiment, even if the speculated resuspension of
solid particles is predicted, neither with VLES nor with
LES. The reasons for the inability of capturing particle
entrainment could be the following:
The particle bed is too thin (12mm implying
45 dust particle layers for the given maximum
particle size of 4004), requiring very high grid
resolution.
The particles are relatively large and gravity has
a strong impact on their motion bringing them to
the pipe floor quickly. Also, in the case of
particles being ejected from the bed, since they
2m/s~
3 m/s
/
, 0ms
....~~~~~~~~~~~~~~~~ .I...............................
\
I
60000
58000
E56000
54ooo
52000
Figure 3.
30 60 90 120 150
t (sec)
Number of particles remaining in the domain
Diameter (pLm)
Figure 4. Number of particles still remaining in the pipe at the
end of the simulation
.Oo~l c~10 m/s
0.0900 250 300 350 400
Diameter (pm)
Figure 5. Percentage of particles transported out based on the
initial number of particles
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
S. Lain, M. Sommerfeld & J. Kissin, Experimental studies
and modelling of fourway coupling in particleladen
horizontal channel flow, Int. Journal qf Heat and Fluid
Flow, Vol. 23, 647656, (2002).
D. Lakehal, LEIS for the Prediction of Turbulent
Multifluid Flows Applied to Thermal Hydraulics
Applications. Nucl. Eng. Design, (available online:
(doi: 10.1016/j.nucengdes.2009. 11.030), 2010.
C. Narayanan, D., Lakehal, L. Botto & A. Soldati,
Mechanisms of particle deposition in a fullydeveloped
turbulent open channel flow, Phys. of Fluids, vol. 15(3),
763, (2003).
C. Narayanan & J. Janssen, Simulation of Particle
Transport/Entrainment in Gas Pipelines, Internal Report
GAS0109, ASCOMP GmbH, (2009).
F. Nicoud & F. Ducros, Subgridscale stress modelling
based on the square of the velocity gradient tensor. Flow,
Turbulence and Combustion, Vol. 62, 183200, (1999).
J. Smart, Determining the Velocity Required to Keep
Solids Moving in Pipelines, The Journal of Pipeline
Engineering, Vol. 6 No. 1 (1st Qtr, 2007).
D. M. Snider, P. J. O'Rourke, & M. J. Andrews, Sediment
flow in inclined vessels calculated using multiphase
particleincell model for dense particle flow, hIt. J.
Multiphase Flow 24, 1359 (1998).
D. M. Snider, An Incompressible ThreeDimensional
Multiphase ParticleinCell Model for Dense Particle
Flows, JCP, Vol. 170, 523549, (2001).
C.G. Speziale, Turbulence modelling for timedependent
RANS and VLES: A review. AIAA Journal, Vol. 36(2),
173184, (1998).
TransAT User Manual, (2010). www.ascomp.ch/transat.
N.A. Tsochatzidis & K.E. Maroulis, Methods help remove
black powder from gas pipelines, Oil and Gas Journal,
March (2007).
are already very close to the wall (%/ mm), they
do not develop any vertical velocity due to
gravitational acceleration. This means that wall
collisions might not contribute to resuspension.
For particles that are much heavier that the fluid
the Saffman lift force is known to be negligible.
This was also tested as part of the model
development effort.
Accounting for wall roughness or dust
nonsphericity might be important to obtain
sustained particle suspension (Lain et al., 2002).
Particle rotation is not accounted for in the
simulations. If the particle rotation is high, then it
can affect the wall rebound characteristics and
also produce the additional Magnus lift force.
Acknowledgements
This work is partly sponsored by GASUNIE (now
KEMA), The Netherlands.
References
F. M. Auzerais, R. Jackson, & W. B. Russel, The
resolution of shocks and the effects of compressible
sediments in transient \c'llline~. J. Fluid Mech. Vol. 195,
437 (1988).
J.T. Davies, Calculation of Critical velocity to maintain
solids in suspension, Chent. Eng. Science, Vol 42(7), 1987
D. Gidaspow, Hydrodynamics of fluidization and heat
transfer supercomputer modeling, Appl. Mech. Rev. Vol.
39, 1, (1986).
S. E. Harris & D. G. Crighton, Solutions, solitary waves
and voidage disturbances in gasfluidized beds, J. Fluid
Mech. Vol. 266, 243 (1994).
S.T. Johansen, J. Wu & W. Shyy, 2004. Filteredbased
unsteady RANS computations. hIt. Journal of Heat and
Fluid Flow, Vol., 25, 1021, (2004).
M. Labois, C., Narayanan & D. Lakehal, A New CFD
approach for urban flow and pollution dispersion
simulation. In Proc. Comp. wind Eng. CWE10, Chapel
Hill, NC, USA, May 2327 (2010).
