7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Paper No
Evaluation of a
Lagrangian Discrete Phase Modeling Approach for Resolving
Cluster Formation in CFB Risers
S.Cloetet, S. Johansent, M. Braun*, B. Popoff*,S. Amini
"f Flow Technology research group, SINTEF Materials and Chemistry, Alfred Getz vei 2, NO7465 Trondheim,
Norway
*Ansys Germany GmbH, Birkemveg 14a, D64295 Darmstadt, Germany
Shahriar.aminii~isintef.no
Keywords: Lagrangian particle tracking, Circulating fluidized bed, Riser, Particle size distribution
Abstract
A computational fluid dynamics (CFD) model was developed to simulate the hydrodynamics of gassolid flows in a
circulating fluidized bed (CFB) riser using a Lagrangian discrete phase approach. This new development extends
the standard discrete phase model (DPM) to account for particleparticle interaction effects and to include the
dispersed phase volume fraction in the flow field, thereby making it applicable to denser particulate flows. The
model was evaluated by comparing its predictions with experimental results reported for a CFB riser (10 m high
and 76 mm ID) operating at solids circulation rate of 100 kg/m2 S and a gas superficial velocity 3.5 m/s. The model
was capable of predicting the main features of the complex gassolids flow, including the cluster formation of the
solid phase along the walls. Cluster formation could be resolved on substantially coarser grids than would be
possible with the traditional twofluid model. Results compared well to experimental data even when using grid
cells as wide as 60 particle diameters. When compared
computational time of almost one order of magnitude.
Introduction
Numerical modelling of riser flows has been a topic of
great interest in recent years. Traditionally, an
EulerianGranular approach, based on the kinetic theory of
granular flows (KTGF) has been implemented (Lun
Savage et al. 1984: Gidaspow, Bezburuah et al. 1992:
Syamlal, Rogers et al. 1993). In this approach, the random
particle motion is likened to the thermal motion of
molecules in a gas. A granular 'temperature' is introduced
as a central concept to quantify the energy contained
within these random particle motions and is subsequently
used to estimate fluid properties such as pressure and
viscosity for the granular phase. Thus, no actual particles
are modelled and both the gas and granular phases are
assumed to be a continuum. This approach is conunonly
referred to as the two fluid model (TFM).
On the other side of the spectrum, discrete element
methods (DEM) have also recently been used to model
riser flows (Zhang, Chu et al. 2008: Zhao, Cheng et al.
2010). Within this methodology, the particles are tracked
individually within a Lagrangian framework according to
Newton's laws of motion, while the surrounding gas phase
is modelled as a continuum. Particleparticle interactions
are also directly simulated, making this approach much
more fundamental than the TFM. It is, however, very much
limited by the amount of particles that can be tracked with
present computational capacities.
A third method, lying somewhere between the TFM and
DEM approaches, is also available. This method, known as
to the two fluid model, this can lead to a decrease in
the dense discrete phase model (DDPM), also tracks
particles through the domain in a Lagrangian fashion, but
does not simulate particleparticle interactions directly.
Particleparticle interactions are estimated by the KTGF
approach. This method allows the simulation of a much
larger quantity of particles than the DEM and has some
distinct fundamental advantages over the TFM.
The wide particle size distribution often present in
riser flows can be implicitly accounted for.
The phenomenon of 'delta shocks' (Passalacqua
and Fox 2009) observed in the simulation of
crossing dilute jets with the TFM is
fundamentally avoided.
Wall interactions, another topic of substantial
uncertainty in the TFM, can be more
fundamentally accounted for.
The DDPM approach as applied to granular flows is still in
an early stage of development, but holds a large amount of
promise due to these advantages.
Nomenclature
C Coefficient
d Diameter (m)
epp Particleparticle restitution coefficient
jForce rector per unit volume (N/nr')
g Gravity vector (ms1i'
g, Radial distribution function
K Interphase exchange coefficient
Paper No
ke, Diffusion constant
#2 Mass transfer rate per volume (kg/s/m )
P Pressure (Pa)
Re Reynolds number
t Time (s)
Greek letters
a Volume fraction
Rate of energy exchange (W/m )
Te> Dissipation rate of fluctuating energy (W/m )
PU Viscosity (Pa.s)
8 Granular temperature (nr1 5)
P Density (kg/m )
z Stressstrain tensor
tJ Velocity vector (m/s)
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
The particle equation of motion is solved for each particle
in the form:
duG 1VI+,(~)
dt p,
(1.3)
g9 /))
+ F + n;~eeract on
The right hand side terms represent the pressure force, drag
force, gravitational force, any additional force and the
particleparticle interaction force.
The drag force is calculated as follows:
18p CD Rep
F, = '
Dp~d 24
(1.4)
Subscripts
DPMl Discrete phase model
P Phase p or Particle/Solids
4 Phase q
D Drag
Abbreviations
DDPM Dense discrete phase model
DEM Discrete element method
DPM Discrete phase model
KTGF Kinetic theory of granular flows
pd Particle diameters
TFM Two fluid model
with the drag coefficient modelled according to (Wen and
Yu 1966):
(1.5)
The interaction force is estimated from the solids pressure:
1
EInteraction V ~p
(1.6)
where (Lun, Savage et al. 1984):
P,=aYp + 2Pp (1+ e ) a go
(1.7)
The granular temperature used in the KTGF is calculated
aS follows:
dt nplmse.
St <=1
s
(1.1)
3 6a~,,~~~~
2p St V' + V.k) y, + ,
(1.8)
ap~ +
V [ayp (V6, + v ~+ a~pY f +F +
C K, v ~v ) + mg2 6 mY2 ) +
Here, the right hand side terms represent the generation of
fluctuating energy by the solid stress tensor, the diffusion
of fluctuating energy, the collisional dissipation of
fluctuating energy (Lun, Savage et al. 1984) and the energy
exchange between the fluctuating particles and any
additional phases (Gidaspow, Bezburuah et al. 1992). This
equation was solved algebraically by neglecting the
convection and diffusion terms to improve numerical
stability.
(1.2)
K DPTI (DPTI 6 V) + DPTI,eplcit
The conservation equations are not solved for the
particulate phase, but the appropriate volume fraction or
velocity values are taken directly from the particle field.
C = 1+ 0.15a R
G, R
Simulations
Model Equations
In its standard form, the discrete phase model (DPM) does
not account for the volume fraction of the discrete phase
particles. The DDPM formulation (Popoff and Braun 2007)
overcomes this limitation by solving a set of conservation
equations for multiple phases (generalized form written
below for phase p).
Paper No
The solids stress tensor in equation (1.8) is written as
follows:
Here, the solids pressure and the bulk viscosity iS
calculated according to (Lun, Savage et al. 1984) and the
shear viscosity according to (Syamlal, Rogers et al. 1993).
Within these formulations, the radial distribution function
is calculated according to (Ogawa 1980).
Computational Domain
A 2D, rectangular section of the riser was modelled under
the assumption of fully developed flow. The section waS
0.076 m in width and 0.4 m in height and was meshed with
completely square structured grid cells according to the
specific simulation being processed.
Boundary and Initial ConditionS
The top and bottom boundaries of the domain were
designated as translationally periodic. A pressure gradient
was specified across the domain to drive the flow against
gravity. This pressure gradient was constantly adjusted by
means of a negative feedback controller to keep the gas
superficial velocity moving through the periodic domain
within + 1% of the given set point.
Side boundaries of the domain were designated as walls
with a noslip boundary condition for the gas phase and a
standard Johnson and Jackson (Johnson and Jackson 1987)
boundary condition for the particles. The specularity
coefficient used in the Johnson and Jackson boundary
condition was set to 0.5 and the discrete phase particle
reflection coefficients were set to 0.9 both in the tangential
and normal directions.
The solution was initialized with a constant volume
fraction and particle field as well as a constant vertical
velocity for both phases. These values were calculated
according to the specific simulation run in question.
Flow Solver
The commercial CFD solver package, FLUENT 12.1 was
used in the calculations. The phasecoupled SIMPLE
algorithm was selected for pressurevelocity coupling,
while the second order scheme was employed for
discretization of all remaining equations, except for the
volume fraction which was discretized according to the
QUICK scheme. Second order implicit timestepping was
used.
Experimental Data
Experimental data for comparison was found from the
literature (Yan and Zhu 2004; Yan and Zhu 2005). Time
averaged void and velocity profiles were collected in a tall
(10 m) riser in which the flow could generally be seen as
fully developed above a riser height of 5m. The average
particle diameter used was 67 pLm and the particle density
was 1500 kg/m3. Experimental radial volume fraction
profiles were collected at five measuring heights above 5
m. These were subsequently averaged to form the best
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
estimation of the void fraction and solids velocity profiles
in the fully developed region. This procedure results in the
void and velocity profiles shown in Figure 1 for the case
considered in the present study (flux 100 kg/m s and
superficial velocity 3.5 m/s).
0,18 6
0,16
0,14
.p 0,12 4 4
S0,1
0,0 solids volu me fraction a
0,06  Solidsvelocity 2,
S0,04 *
0,02
0 0,2 0,4 0,6 0,8 1
Normalized radial distance
Figure 1: Experimental time averaged data (Yan and
Zhu 2004; Yan and Zhu 2005).
Data was collected in a 3D cylindrical geometry, but
modelled on a 2D planar geometry. In order to adjust for
this difference, the radial profiles for void fraction and
velocity were averaged in 2D space to provide an adjusted
value of the average void fraction and the superficial
velocity. The superficial velocity was calculated from the
solids velocity curve under the assumption of a constant
slip velocity across the radius of the riser. When looking at
Figure 1, this adjustment would increase the superficial
velocity and decrease the average void fraction since the
central areas of the riser will now weigh the same as the
wall areas. For example, the void profile in Figure 1 would
result in an average voidage of 5.65% when averaged in
3D space and 3.72% when averaged in 2D space. The
target flux can also be adjusted in this fashion to a value of
99 kg/m s. Also note that, for ease of discussion, the terms
'radial profile' and 'radius' will be used throughout the text
even though they are not technically applicable to a 2D
planar geometry.
Results and Discussion
Two sets of numerical experiments were performed. Since
this is a new approach, the first set was a basic
independence study establishing solution independence in
terms of grid size, timestep size, iterations per timestep and
the number of particles. The second set of experiments
focused solely on determining the maximum grid size on
which cluster formation could be adequately resolved.
Results obtained from these experiments could then be
compared to simulations using the TFM.
Run Granular
Flux tmeaue
(kg/nrs) teprtr
(I I
Centre pint 91.68 0.00449
Grid: 14.14 pd87.20 0.00547
Grid: 28.28 pd85.59 0.00363
Timestep:2e5 93.49 0.00456
Timestep: 8e5 unstable unstable
Iter pr TS: 24 90.65 0.00467
Iter pr TS: 6 125.13 0.00416
Part prcell: 64 91.46 0.00384
Part prcell: 16 98.59 0.00562
When looking at the flux results, outputs from the
simulation compare well with the target flux of 99 kg/m s.
Flux results from all simulations also reside within +10%
from the centre point, except for the experiment with the
low number of iterations per timestep. The solids flux is a
very important measure to be extracted from the model and
is the most important dependent variable for which to
establish numerically independent results. Thus, when
looking at the flux, the solution can be seen to be
numerically independent even at the coarsest levels of all
the independent variables, except for the number of
iterations per timestep. This can be confirmed by looking
at the time averaged radial profiles of solids volume
fraction and velocity shown in Figure 2.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
0,18
SExperimental
0,16 Centrepoint
*Fine Grld
0,14 Small Timestep
0,1 oo Many Iterations per TS
Many particles per cell
o,as
O ,06
O,04
O,02
0 0,2 0,4 0,6 0,8
Normalized radius
10
SExpenimental
.. Centre point
.. *Fine Grld
Small Timestep
*~Many Iterations per TS
_~ Many particles per cell
P4
ll 0,2 0,4 0,6 O,8
Normalized radius
Figure 2: Comparison of radial solids volume fraction
and velocity between the centre point and the finer
numerical settings. Experimental data from (Yan and
Zhu 2004; Yan and Zhu 2005).
The averaged granular temperature, on the other hand,
shows grid and particle per cell dependent behaviour even
at the finest levels investigated. These dependence can be
explained by noting that both a larger amount of particles
and fewer cells will smooth out the particle distribution,
leading to somewhat smaller gradients in the particle
volume fraction distribution. Granular temperature is
primarily generated in the denser regions, more of which
will be present in a more segregated volume fraction field.
Thus it would be expected that smaller cells and a lower
amount of particles per cell will lead to a larger production
of granular temperature. It seems, however, that the failure
to reach numerical independence in terms of granular
temperature is not reflected directly on the more important
parameter of solids flux. This is encouraging and allows
for investigations into the performance of the DDPM on
even coarser grids.
Paper No
Independence study
The four independent variables under consideration were
evaluated around a centre point specified as follows:
Grid size: 20 times the average particle diameter
Time step size: 4e5 S
Iterations per timestep: 12
*Number of particles: 32 particles per cell
A high and low point of each of these variables was
evaluated by doubling or halving the value at the centre
point. For the grid size the value was adjusted such that the
number of cells was doubled or halved.
Results were quantified in terms of the average flux
through the domain and the average granular temperature
in the domain. An average value of these quantities over
the domain was collected once every ten timesteps and
subsequently averaged over time as well. Results are
displayed in Table 1.
Table 1: Results
experiments.
from the independence numerical
Grid size (k s t rtaulrre
14.14 pd87.20 0.00547
20 pd91.68 0.00449
28.28 pd85.59 0.00363
40 pd100.90 0.00204
56.57 pd95.67 0.00121
80 pd115.48 0.00081
It is shown that the overall solids flux is predicted with
satisfactory accuracy up to a grid size of 80 particle
diameters. The granular temperature, however, shows
strongly grid dependent behaviour once again, but still
makes no impact on the solids flux. The adequate capturing
of the average solids flux is quite encouraging since the
periodic nature of the simulations will expose any
discrepancy much more clearly than would be observed in
a full riser simulation. If the momentum interaction
between solids and gas is underpredicted in a real riser the
volume fraction would only increase slightly, creating a
larger surface area for the drag force to work on and
therefore increasing the momentum interaction again. The
average volume fraction is fixed in the periodic domain,
meaning that this countering effect is not possible.
Radial void and velocity profiles for the cases in Table 2
are shown in Figure 3. Here, it is seen that the correct time
averaged system behaviour is maintained up to a grid size
of 56.57 particle diameters. The solids velocity profile at a
grid width of 80 particle diameters is also adequate, but the
solids void fraction profile starts to become too flattened
out. It can therefore be concluded that the overall time
averaged behaviour of the system is predicted correctly up
to a grid width around 60 particle diameters.
o,oo
0 0,2 0,4 0,6 0,8
Normalized radius
10 a Expenmental
1414 pd
a *20 pd
...28 28 pd
O,*40 pd
6 FPI\1 56 57 pd
P4
S0,2 0,4 0,6 0,8
Normalized radius
Figure 3: Comparison of radial solids volume fraction
and velocity for the different grids in Table 2.
Experimental data from (Yan and Zhu 2004; Yan and
Zhu 2005).
Paper No 7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Coarse grid simulations
0,18
SExpenimental
Following the grid independence of the overall solids fluX 14 14 pd
as well as the radial solids void and velocity profiles, some ,6*20 pd
simulations on coarser grids were carried out as well. The 0,14 2828pd
grid sizes used in these simulations were increased from = *40 pd
oP 0,12 5657 pd
the coarsest grid investigated in the independence study so *0p
that each successive simulation would have half the 0,1
number of grid cells of the previous. Simulation results are
presented in Table 2. p o,as
( O,06
Table 2: Results from the coarse grid numerical u,
experiments. o,o4 .. *
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
0.58
0.56
0.54
0.52
S 060.50
0.49
0.47
0.45
0.43
0.41
0.40
0.38
0.36
0.34
0.32
0.31
0.29
0.27
0.25
0.23
0.22
0.20
0.1 8
0.1 6
0. 14
0. 13
0.1 1
0.09
0.07
0.05
0.04
0.02
0.00
Figure 5: Cluster resolution with the DDPM on
progressively coarser grids.
When directly compared to a TFM simulation performed
on the same grid, the advantage of using the DDPM is
clear. Figure 6 shows a comparison of the cluster
resolution attained on the 40 pd grid. The difference in the
results from the two modeling approaches is evident in that
the DDPM resolves the clusters much more clearly than
the TFM on this relatively coarse grid.
The main reason for the reduced grid dependence of the
DDPM approach is the Lagrangian tracking of the particles.
It provides a mechanism to use the accuracy of particle
tracks in the context of a multifluid formulation for the
conservation equations. Grid dependence will increase
with a decrease in particle diameter because the grid based
gas velocity is more dominant for small particles. On the
other hand, for particles with a large diameter and higher
inertia, even larger grid independence advantages can be
expected.
Paper No
Cluster resolution
Granular simulations start to fail on coarser grids because
clusters are not being sufficiently resolved. When solved
on a sufficiently fine grid, particles would cluster together
to form dense masses that fall downwards much more
readily, thereby reducing the flux through a periodic
domain. This effect was investigated recently in the same
flow system by using the TFM and quantified as illustrated
in Figure 4.
Figure 4: The response of solids flux to changes in grid
width (GW) and cell aspect ratio (AR) (Cloete, Amini et
al. 2010). The axes are given in coded variables. For
GW, 1.414 5 pd, O 10 pd and 1.414 20 pd. For
AR, 1.414 1 O 2 and 1.414 4.
Figure 4 shows that reasonable grid independence in terms
of solids flux is attained at a grid width of 10 particle
diameters (GW = 0) and an aspect ratio of 2 (AR = 0).
When the aspect ratio of the cells are set to 1 (AR =
1.414), a grid width of 20 particle diameters (GW =
1.414) can be used to give approximately the same flux.
This situation (grid width of 20 pd at an aspect ratio of 1)
was used as the centre point of the independence study
reported earlier. It is also shown in Figure 4, however, that
some strong grid dependent behavior is shown at grid
widths and cell aspect ratios above the aforementioned
points. Since the dense masses of particles are not resolved
on coarser grids, the now dispersed particle field is easily
swept upwards by the rising gas stream and the flux
increased accordingly.
The DDPM approach seems to do much better in this
regard by creating regions of dense particle packing even
on very coarse grids. Naturally, the finer details of
individual clusters are lost on coarser grids, but the actual
effect that cluster resolution has on the overall flow
situation is maintained. This is demonstrated in Figure 5.
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
When looking at the timeaveraged void and velocity
profiles in Figure 7, both plots for the TFM behave as
would be expected if clusters were not resolved. The void
profile plot is flattened out since the dense clusters are not
resolved as clearly towards the walls of the riser. More
solids are therefore simulated in the central, high velocity
regions of the riser and the overall flux through the system
is overpredicated. Despite this flattened out void profile,
the velocity profile is actually more peaked. This is
because the resolution of clusters will have a strong
diffusive effect on the overall momentum transfer,
resulting in flatter, more turbulentlike time averaged
velocity profile such as the one seen for the DDPM.
Computational time
Computational time is the largest problem in solving riser
flows. The fine spatial and temporal resolutions required to
adequately resolve the physics makes even the simulation
of relatively small scale 2D systems rather expensive.
The DDPM is presently significantly slower than the TFM
when solved on the same grid resolution. In addition to the
continuous phase calculations, the discrete phase particles
also have to be tracked, adding some load to the
computation. In addition, the DDPM is presently less
stable than the TFM, requiring substantially lower
underrelaxation factors and therefore more iterations per
timestep to converge. All in all, the DDPM is probably
about 3 times slower than the TFM on a fixed grid.
As shown in the present work, however, the DDPM can
provide adequate results even at a grid size 3 times greater
than the limit for the TFM. In 2D, this translates to 9 times
fewer cells in the domain. Additionally, the timestep size
C8H 81So be increased by a factor of 3. Therefore, when all
is considered, the DDPM can be seen to be able to provide
adequate results in almost 10 times less time than the
Although these simulations were done for mono sized
particles, most CFB risers are operated with particle size
distributions. When a particle size distribution is included,
the gain in computational efficiency will be even greater.
In the TFM, every particle diameter class needs to be
simulated as a separate phase, thereby increasing the
number of transport equations to be solved. Since the effort
of the DDPM approach scales only with the number of
tracked parcels, additional diameter classes can easily be
considered at a low cost.
Paper No
0.56
0 5
0.43
0.41
0.40
0.41
0 4
0.22
0.20
0.18
0 13
0 2
0 1
Fiue6 Cmaio o lserrsltonwt h
DDMan h TMo a4 d rd
0,18
0,16
0,14
P0,12
0,1
0 06 8
o,04
0.02
0 0,2 0,4 0,6
Normalized radius
0,8 1
Normalized radius
Figure 7: Time averaged particle void and velocity
profiles for the DDPM and TFM on a 40 pd grid.
Experimental data from (Yan and Zhu 2004; Yan and
Zhu 2005).
7th International Conference on Multiphase Flow
ICMF 2010, Tampa, FL USA, May 30June 4, 2010
Yan, A. and J. Zhu (2l1***4. "ScaleUp Effect of Riser
Reactors (1) Axial and Radial Solids
Concentration Distribution and Flow
Development." Industrial & Engineering
Chemistry Research 43(18): 58105819.
Yan, A. and J. Zhu (2005). "ScaleUp of Riser Reactors:
Particle Velocity and Flow Development." AIChE
Journal 51(11): 29562964.
Zhang, M. H., K. W. Chu, et al. (2008). "A CFDDEM
study of the cluster behavior in riser and downer
reactors." Powder Technoloev 184(2): 151165.
Zhao, Y., Y. Cheng, et al. (2010). "EulerianLagrangian
simulation of distinct clustering phenomena and
RTDs in riser and downer." Particuoloev 8(1):
4450.
Paper No
ConclusionS
The dense discrete phase modeling approach to riser flows
was evaluated in this work. In this approach, the particulate
phase is tracked in a Lagrangian framework and the
particleparticle interactions are modeled according to the
kinetic theory of granular flows. The total flux as well as
time averaged radial solids velocity and volume fraction
profiles calculated with this method matched well with
experimental results.
A general independence study was done to establish that
the solution was independent of the grid size, timestep size,
the number of injected particle parcels and the number of
iterations per timestep. It was found that the solution
retained grid independency at resolutions substantially
coarser than that achieved by the two fluid model.
Further investigations showed that the dense regions of
clusters were resolved even up to grid sizes in excess of 60
particle diameters. The resolution of clusters is the most
important factor determining the momentum interaction in
risers. This finding suggests that accurate results can be
attained in the DDPM framework using a computational
load of roughly ten times lower than that required by the
two fluid model.
References
Cloete, S., S. Amini, et al. (2010). "On the Effect of
Cluster Resolution in Riser Flows on Momentum
and Reaction Kinetic Interaction." To be
submitted.
Gidaspow, D., R. Bezburuah, et al. (1992). Hydrodynamics
of Circulating Fluidized Beds, Kinetic Theory
Approach. 7th Engineering Foundation
Conference on Fluidization 7582.
Johnson, P. C. and R. Jackson (1987).
"FrictionalCollisional Contitutive Relations for
Granular Materials, with Application to Plane
Shearing." Joumnal of Fluid Mechanics 176:
6793.
Lun, C. K. K., S. B. Savage, et al. (1984). "Kinetic
Theories for Granular Flow: Inelastic Particles in
Couette Flow and Slightly Inerlastic Particles in a
General Flow Field." Joumnal of Fluid Mechanics
140: 223256.
Ogawa, S. U., A.; Oshima, N. (1980). "On the Equation of
Fully Fluidized Granular Materials." Joumnal of
Applied Mathematics and Physics 31: 483.
Passalacqua, A. and R. O. Fox (2009). Multiphase CFD for
GasParticle Flows: Beyond the TwoFluid Model.
Seventh Intemnational Conference on CFD in the
Minerals and Process Industries. CSIRO,
Melboumne, Australia.
Popoff, B. and M. Braun (2007). A Lagrangian Approach
to Dense Particulate Flows. 6th Intemnational
Conference on Multiphase Flow. Leipzig,
Germany.
Syamlal, M., W. Rogers, et al. (1993). MFIX
Documentation: Volume 1. Theory Guide.
Springfield, National Technical Information
Service.
Wen, C. Y. and Y. H. Yu (1966). "Mechanics of
Fluidization." Chemical Engineering Progress
Symposium Series 62: 100111.
