University of Massachusetts Amherst

ScholarWorks@UMass Amherst
Masters Theses 1911 - February 2014

2013

Optimization of Mixing in a Simulated Biomass

Bed Reactor with a Center Feeding Tube
Michael T. Blatnik
University of Massachusetts Amherst

Follow this and additional works at: [Link]

Part of the Catalysis and Reaction Engineering Commons, and the Energy Systems Commons

Blatnik, Michael T., "Optimization of Mixing in a Simulated Biomass Bed Reactor with a Center Feeding Tube" (2013). Masters Theses
1911 - February 2014. 1108.
Retrieved from [Link]

This thesis is brought to you for free and open access by ScholarWorks@UMass Amherst. It has been accepted for inclusion in Masters Theses 1911 -
February 2014 by an authorized administrator of ScholarWorks@UMass Amherst. For more information, please contact
scholarworks@[Link].
OPTIMIZATION OF MIXING IN A SIMULATED BIOMASS
BED REACTOR WITH A CENTER FEEDING TUBE

A Thesis Presented
by

MICHAEL T. BLATNIK

Submitted to the Graduate School of the

University of Massachusetts Amherst in partial fulfillment
of the requirements for the degree of

MASTER OF SCIENCE IN MECHANICAL ENGINEERING

September 2013

Mechanical and Industrial Engineering

 OPTIMIZATION OF MIXING IN A SIMULATED BIOMASS
BED REACTOR WITH A CENTER FEEDING TUBE

A Thesis Presented

by
MICHAEL T. BLATNIK

Approved as to style and content by:

Stephen M. de Bruyn Kops, Chair

David P. Schmidt, Member

T. J. Mountziaris, Member

Donald L. Fisher, Department Chair

Mechanical and Industrial Engineering
 ACKNOWLEDGMENTS

This work was funded by the National Science Foundation Office of Emerging Frontiers

in Research and Innnovation (EFRI) grant number 0937895 and greatly benefited from
discussions with Prof. George Huber.
High performance computing resources were provided by the U.S. Department of Defense

High Performance Computing Modernization Program at the Army Research Laboratory

(ARL) and at the Army Engineer Research and Development Center (ERDC).

iii
 ABSTRACT

OPTIMIZATION OF MIXING IN A SIMULATED BIOMASS

BED REACTOR WITH A CENTER FEEDING TUBE

SEPTEMBER 2013

MICHAEL T. BLATNIK
[Link]., LYNCHBURG COLLEGE
M.S.M.E., UNIVERSITY OF MASSACHUSETTS AMHERST

Directed by: Professor Stephen M. de Bruyn Kops

Producing gasoline-type fuels from lignocellulosic biomass has two advantages over pro-

ducing alcohol-type fuels from plant sugars: gasoline has superior fuel characteristics and
plant lignin/cellulose does not compete with human food supplies. A promising technology
for converting lignocellulose to fuel is catalytic fast pyrolysis (CFP). The process involves in-

jecting finely ground biomass into a fluidized bed reactor (FBR) at high temperatures, which
reduce the biomass to gases that react inside the catalyst particles. This entails complex
hydrodynamics to efficiently mix a stream of biomass into a catalyst bed that is fluidized

by a separate stream of inert gas. Understanding the hydrodynamics is complicated by the

fact that the entire process occurs inside a heavily insulated, opaque, reactor vessel. Nu-
merical simulations offer a promising approach to understanding, predicting, and optimizing

hydrodynamic mixing in a CFP biomass reactor.

iv
 The purpose of this research is to understand the simulation techniques and statistical

measures appropriate for quantifying mixing in a CFP biomass reactor. The methodology
is validated against the canonical configuration of a non-reacting, single-inlet fluidized bed.
A new finding is that the minimum bubbling velocity may be predicted by a significant

increase in temporal variance of the pressure drop. The methodology is then applied to a
non-canonical FBR in which biomass is injected into the catalyst bed via a vertical center
tube. Since no hydrodynamic mixing data exist from laboratory experiments, mixing is

inferred from the aromatics yield from the laboratory reactor. Flow configurations with
which simulations demonstrate the best mixing have the highest aromatic yields in the
experiments. The simulations indicate that when the bed is in the bubbling regime, the

gasified biomass from the center tube is efficiently mixed radially throughout the catalyst
bed. If the flow rate of inert gas is insufficient to bubble the bed, then the gasified biomass
exits the center tube, reverses direction, and flows upward along the tube’s outside wall.

Provided the bed is bubbling due to the inert gas stream, the upper limit on the flow
through the center tube, and thus the aromatic yield potential, has yet to be determined.

v
 TABLE OF CONTENTS

Page

ACKNOWLEDGMENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii

ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iv

LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix

LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . x

CHAPTER

1. BACKGROUND . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1

1.1 Fluidized Beds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.2 Flow Regimes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 Geldart Classification of Particles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.4 Gas Introduction Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.5 Deviations from Canonical Beds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

1.5.1 Horizontal Tubes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

1.5.2 Vertical Tubes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2. THEORY AND MODEL . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.1 Governing Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2.2 Closure Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.2.1 Granular Temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2.2.2 Solids Pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2.3 Viscous Stress . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2.4 Solids Thermal Conductivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2.5 Frictional Stress . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2.6 Radial Distribution Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

vi
 2.2.7 Drag Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.3 Semi-empirical correlations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.3.1 Minimum fluidization velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.3.2 Minimum Bubbling Velocity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.3.3 Bed expansion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3. SIMULATION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

3.1 Modeling methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

3.1.1 Lagrangian Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

3.1.2 Eulerian Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.2 Simulation of Geldart Type A Particles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.3 Physical Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
3.4 Geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.5 Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

3.5.1 Walls . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.5.2 Axis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.5.3 Inlets . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.5.4 Outlet . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.5.5 Internal field . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

3.6 Numerical Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

3.7 Quantifying Mixing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

3.7.1 Tracer Gases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

3.7.2 Mixture Fraction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
3.7.3 Volume-weighted Statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.7.4 Statistical Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.7.5 Volume-weighted Mean . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.7.6 Volume-weighted Variance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
3.7.7 Bed Mean and Variance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

4. RESULTS AND DISCUSSION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

4.1 Standard Configuration Fluidized Bed Simulations . . . . . . . . . . . . . . . . . . . . . . . . . 39

4.1.1 Determination of the Minimum Fluidization Velocity . . . . . . . . . . . . . . . . . 39

4.1.2 Determination of the Minimum Bubbling Velocity . . . . . . . . . . . . . . . . . . . 42

vii
 4.2 Effect of Gas Distribution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

4.2.1 Visualization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.2.2 Statistical Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.2.3 Statistical Analysis of the Catalyst Volume Fraction . . . . . . . . . . . . . . . . . 47
4.2.4 Statistical Analysis of the Mixture Fraction . . . . . . . . . . . . . . . . . . . . . . . . . 52
4.2.5 Bed Statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

4.3 Effect of Total Gas Flow Rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

4.3.1 Visualization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.3.2 Statistical Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
4.3.3 Statistical Analysis of the Catalyst Volume Fraction . . . . . . . . . . . . . . . . . 56
4.3.4 Statistical Analysis of the Mixture Fraction . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.3.5 Bed Statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

5. CONCLUSIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

APPENDIX: GOVERNING EQUATIONS AND CLOSURE MODELS . . . . . 66

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

viii
 LIST OF TABLES

Table Page

3.1 Helium, argon, and nitrogen density and dynamic viscosity vs.
temperature . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

3.2 Summary of 2D axisymmetric meshes of the reactor with and without the
tube . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24

3.3 Flow rate configurations for 3 sets of CFD simulations: a set without a tube
in which the fluidizer velocity is varied, a set of runs with the tube and a
near-constant Qtotal as a function of Qf t , and a set of tube simulations
with a constant Qf t /Qtotal as a function of Qtotal . . . . . . . . . . . . . . . . . . . . . . . . 26

3.4 Under-relaxation factors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

3.5 Ideal N2 tracer gas mass fraction and mixture fraction for each simulation
with the feeding tube . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

4.1 Flow rate configurations and abbreviations for simulations with a

near-constant Qtotal as a function of Qf t . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

4.2 Flow rate configurations and abbreviations for simulations with a constant
Qf t /Qtotal as a function of Qtotal . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54

A.1 Governing equations and closure models: I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

A.2 Governing equations and closure models: II . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

A.3 Governing equations and closure models: III . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

ix
 LIST OF FIGURES

Figure Page

3.1 CFD meshes at ∼ 30 cells/cm mesh resolution for meshes without a tube
(left) and with a tube (right) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

4.1 Minimum fluidization velocity determination via pressure drop vs. fluidizing
velocity for Fluent simulations (bars indicate 2 standard deviations above
and below the mean), the experimental reactor, and literature models . . . . . 41

4.2 Catalyst volume fraction profiles for (top, from left to right) vf l = 0.2, 0.3,
0.4, 0.6, and 0.8 cm/s and (bottom, from left to right) vf l = 0.9, 1.0, 1.5,
2.0, and 3.0 cm/s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

4.3 Four measures of the transition to the bubbling regime: (left) ε and Ω,
(right) χ and ψ for tubeless Fluent simulations, where
h−hmf Vαs ≤0.15 hα02s i
ε = αhα si
s,max
, Ω = hmf
, χ = Vbed
, and ψ = α02
. . . . . . . . . . . . . . . . . . . . . . 43
s,max

4.4 Time snapshots of contour plots of (a) the catalyst volume fraction αs , (b)
the feed tube tracer gas mass fraction YAr , and (c) the catalyst axial
velocity magnitude vs for Fluent simulations with a feed tube with a
constant total flow rate Qtotal ≈ 55 cm3 /s and different ratios of
Qf t /Qtotal (left to right: 8%, 17%, 33%, 51%, and 68%) . . . . . . . . . . . . . . . . . . 46

4.5 Statistical convergence of bed mixing parameters (ε, ψ, θ, and Φ) vs.

sampling period (TS ) for simulations run with Qf t /Qtotal ≈ 8%, 17%,
3
33%,D51%, E and 68% withD aEconstant QtotalD ≈ 55E cm /s, where
02
ε = αs /αs,max , ψ = αs02 /αs,max , θ = ξ(r) /ξideal , and
D E D E2
Φ = ξ 02 / ξ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

4.6 The effect of the gas distribution between the fluidizer and the feed tube
(Qf t /Qtotal (8%, 17%, 33%, 51%, and 68%; Qtotal ≈ 55 cm3 /s) on the
mixing of the catalyst . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

x
4.7 The effect of the gas distribution between the fluidizer and the feed tube
(Qf t /Qtotal (8%, 17%, 33%, 51%, and 68%; Qtotal ≈ 55 cm3 /s) on the
mixing of the feed tube tracer gas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

4.8 Bed catalyst mean and variance, bed mixture fraction mean and variance
(left axis, solid lines), and aromatic product yield (left axis, dashed line)
vs. the ratio of gas distribution in the feed
D Etube (Qf t /Qtotal
D ) with
E a
02
constant Qtotal ≈ 55 cm3 /s, where ε = αs /αs,max , ψ = αs02 /αs,max ,
D E D E D E2
θ = ξ /ξideal , and Φ = ξ 02 / ξ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

4.9 Time snapshots of contour plots of (a) the catalyst volume fraction αs , (b)
the feed tube tracer gas mass fraction YAr , and (c) the catalyst axial
velocity magnitude vs for Fluent simulations with a feed tube with a
constant Qf t /Qtotal ≈ 33% and different ratios of Qtotal /Qmf :
(left-to-right: 3, 4, 6, and 8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

4.10 Statistical convergence of bed mixing parameters (ε, ψ, θ, and Φ) vs.

sampling period (TS ) for simulations run with Qtotal /Qmf ≈ 3, 4, 6, and
8 with a constant Qf t /Qtotal ≈ 33% . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

4.11 The effect of increasing the total gas flow rate (Qtotal /Qmf ≈ 3, 4, 6, and 8;
Qf t /Qtotal ≈ 33%) on the mixing of the catalyst . . . . . . . . . . . . . . . . . . . . . . . . 59

4.12 The effect of increasing the total gas flow rate (Qtotal /Qmf ≈ 3, 4, 6, and 8;
Qf t /Qtotal ≈ 33%) on the mixing of the feed tube tracer gas . . . . . . . . . . . . . . 60

4.13 Bed mean and variance (solid lines) with aromatic yield (dotted lines) as a
function of Qtotal /Qmf . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

xi
 LIST OF SYMBOLS

β Interphase drag coefficient (kg/(m3 ·s))

g Gravity constant (m/s2 )
qk Heat flux of phase k (W/m2 )
vk Velocity of phase k (cm/s)
τ̄¯k Viscous stress tensor of phase k (N/m2 )
δ Mesh size, uniform in both x and y directions (m)
γs Dissipation of granular energy due to inelastic particle collisions (kg/(m3 ·s))
κdil Dilute phase granular conductivity (kg/(m·s))
κs Solids thermal conductivity (kg/(m·s))
λk Bulk viscosity of phase k (Pa·s)
µg Gas dynamic viscosity (Pa·s)
µs Solids shear viscosity (Pa·s)
µf Frictional viscosity (Pa·s)
Ω Bed expansion, percentage above hmf (dimensionless)
D E D E2
Φ Statistical parameter related to variance of ξ: ξ 02 / ξ (dimensionless)
φ Angle of internal friction/azimuthal angle coordinate position (rad)
D E
02
ψ Statistical parameter related to variance of αs : αs02 /αs,max (dimensionless)
ρk Density of phase/species k (kg/m3 )
Θ Granular temperature (m/s2 )
D E
θ Statistical parameter related to mean of ξ: ξ(r) /ξideal (dimensionless)
D E
ε Statistical parameter related to mean of αs : αs /αs,max (dimensionless)
αk Volume fraction of phase k (dimensionless)
αs,max Maximum solids volume fraction packing limit
αs,min Minimum solids volume fraction for frictional model

xii
ξ Mixture fraction (dimensionless)
ξideal Ideal mixture fraction in a well-mixed case (dimensionless)
A Frictional model empirical constant
CD Coefficient of drag
D Bed depth (m)
ds Particle diameter (m)
e Coefficient of restitution
Fv Viscous drag force (N)
fl Abbreviation for the bottom fluidizer inlet
Fr Frictional model empirical constant (Pa)
ft Abbreviation for the feed tube inlet
g0 Radial distribution function
H Reactor height (m)
h Height of the catalyst bed (cm)
hk Specific enthalpy of phase k (J/kg)
HT Height of feed tube above bottom fluidizer (cm)
hmf Bed height at minimum fluidization velocity (dimensionless)
Js Dissipation of granular energy due to velocity fluctuations (kg/(m3 ·s))
K Empirical constant (N)
m Total mass of the catalyst (g)
n Frictional model empirical constant
O Arbitrary, cell-centered quantity (dimensionless)
P Pressure (Pa)
p Frictional model empirical constant
Ps Solids pressure (Pa)
Pf Normal frictional stress (Pa)
Qf l Fluidizer flow rate (cm3 /s)
Qf t Feed tube flow rate (cm3 /s)
Qtotal Total flow rate Qf l + Qf t (cm3 /s)
R Radius of the reactor (cm)
r Radial coordinate (m)
RT Radius of the feed tube (cm)

xiii
T Temperature (K)
t Time (s)
TS Sampling period for statistical analysis (s)
uk x-component of vk for phase k (cm/s)
V (r, t) Volume occupied by an axial slice of the catalyst bed a radial location r and time t
(cm3 )
V (z) Volume occupied by a radial slice of the reactor at a height z (cm3 )
Vbed (t) Volume occupied by the catalyst bed at time t (cm3 )
vf l Feed tube velocity (cm/s)
vf l Fluidizer velocity (cm/s)
vk y-component of vk for phase k (cm/s)
vmf Minimum fluidization velocity (cm/s)
x Cartesian position coordinate (m)
Xi Mole fraction of species i (dimensionless)
y Cartesian position coordinate (m)
Yi∞ Free-stream mass fraction of species i (dimensionless)
Yiideal Ideal, well-mixed mass fraction of species i (dimensionless)
Yi Mass fraction of species i (dimensionless)
z Axial coordinate (m)
f# Abbreviation for a simulation run with Qf t /Qtotal = #%
Re Relative Reynolds number
t# Abbreviation for a simulation run with Qtotal /Qmf = #

xiv
 CHAPTER 1

BACKGROUND

The goal of this research is to numerically simulate and statistically quantify the hydro-

dynamics of a non-standard fluidized bed using catalytic fast pyrolysis (CFP) for the produc-
tion of gasoline-grade fuel from cellulose. The reactor is under development by the Emerging
Frontiers in Research and Innovation (EFRI) group at both the University of Massachusetts

Amherst and the University of Wisconsin−Madison. From a downward-facing center feed-

ing tube, cellulose or other wood feedstocks are injected into a fluidized bed partially filled
with zeolite (ZSM5) catalyst particles. The process is maintained at high temperatures (773

K), which causes the biomass to pyrolize into oxygenated vapors. As the vapors enter the
catalyst pores, they are converted to aromatics, the desired product, as well as olefins, CO,
CO2 , water, and other by-products. The aromatics are desirable because they can be con-

verted to a gasoline-range hydrocarbon fuel. Advantages of this emerging technology include

short mixing residence times (2−10 s), low cost of the catalyst, an abundance of feedstock,
a ready-to-use fuel product, and high oil prices.

Numerical simulations are performed using computational fluid dynamics (CFD) codes in
c Fluent, version 14.0. We neglect chemical reactions since
the commercial software ANSYS
the hydrodynamic mixing is the study of interest. To assess the accuracy of the simulations,

we compare our results with data from the laboratory reactor and literature models.

1
1.1 Fluidized Beds
A typical fluidized bed is a cylinder, containing a bed of particles supported by a bottom
plate through which a fluidizing gas is introduced [30]. At sufficient velocities, the gas will
suspend or lift the particles, causing them to move throughout the bed. The gas velocity

required to support the weight of the bed is called the minimum fluidization velocity vmf .
The value is empirically determined by increasing the gas velocity until the pressure change
(∆P ) across the height of the bed equals the weight of the bed over the inlet surface area.

Momentum losses due to collisions between the particles and the walls or bottom plate cause
∆P > bed weight/area. The fluctuations are primarily due to the center of mass of the bed
accelerating up and down.

After minimum fluidization velocity, increases in gas velocity result in “bubbles,” which
are voidage regions with low particle density that rise to the top of the bed. The behavior
of such fluidization depends upon the properties of the solids and of the gas, the fluidizing

velocity, and the geometry of the bed.

1.2 Flow Regimes

Fluidization flow regimes are divided into 4 categories: expanding, bubbling, slugging,
and turbulent [30]. The expanding regime has no bubbles. Increases in velocity yield a higher

bed height and thus a lower average solids volume fraction within the bed. Bed expansion
does not always occur before bubbling but is typical for smaller particles [29].
The bubbling flow regime occurs when the fluidizing velocity rises above vmb , the mini-

mum bubbling velocity. At this point, bubbles begin to form that grow with bed height but
typically do not reach full bed width, except in narrow, tall beds. Higher fluidizing velocities
generally lead to larger bubbles. A bubbling bed is an essential component of mixing within

the bed. In industry, fluidized beds are used to facilitate chemical processes such as coal

2
gasification, combustion, and biomass conversion. Ideal mixing maximizes the stoichiometric

amount of desired products [18].

The slugging flow regime occurs at high fluidizing velocities when bubbles (“slugs”) in
deep, narrow beds grow to the diameter of the bed. Slug flow can further be divided into

axisymmetric slugs, wall slug, plugs, and continuous slug flow [20]. At even higher velocities,
the bed can transition from either a bubbling or slugging flow to a turbulent regime. In this
scenario, the bubbles or slugs break up into smaller bubbles and the particles form clusters

or streamers. At this point, the solids and the gas phase become more indistinguishable and
the flow more dilute [20]. Every fluidized bed does not necessarily go through each of these
four stages. The uniqueness of fluidization behavior depends upon the particles’ properties.

1.3 Geldart Classification of Particles

Geldart classified particles into 4 different categories (Geldart Type A, B, C, and D),
depending on the particle’s size, density, fine content, and cohesiveness [29]. Particles of the
same type tend to progress through the same fluidization regimes.

Type A particles, such as fine catalysts, have small diameters and low densities (ρs < 1400
kg/m3 ). Type A particle beds undergo expansion before bubbling, and thus vmb /vmf > 1.
In addition, bubble growth is limited to a fixed maximum bubble size [20], which is relevant

whether or not slugs can form.

Geldart Type B particles exist within a range of a mean particle diameters (40 < ds < 500
µm) and a range of densities (4000 kg/m3 < ρs < 1400 kg/m3 ). For particles of this type,

bubbling begins at minimum fluidization velocity, i.e., vmb /vmf = 1.

Particles belonging to Type C consist of fine, cohesive powders that are difficult to fluidize.
In small-diameter reactors, the particles lift as plugs. With larger beds, channeling occurs

in which the gas bypasses a majority of the particles [30]. Plugging and channeling make
Type C particles non-ideal for mixing situations. Lastly, Type D contains large or very dense

3
particles (ds > 1000 µm), which are prone to spouting. This causes an abrupt transition
from the bubbling to the turbulent flow regime.

The catalyst particles used in the physical reactor have a density of 1750 kg/m3 and the
particle diameter has a bimodal distribution with peaks between 40−80 and 120−160 µm.

Type A particles are preferable for the experimental reactor because they lack cohesion and
facilitate bed expansion. For our numerical simulations, we use a diameter of 100 µm, which
means that the particles can be classified as either Geldart Type A or Geldart Type B. We

show in Chapter 4 that the simulations capture the particle behavior that occurs within the
physical reactor with regard to vmf .

1.4 Gas Introduction Methods

There are several commonly used bed configurations for introducing gas or other feed-
stocks into the bed besides a bottom plate distributor. One method is to introduce an

upward-facing central jet of higher velocity in combination with the uniform velocity plate
distributor maintained at vmf . Feeding tubes are often used for inputting gases into the
bed. Experimental [54, 7, 31, 68] and computational [46, 45, 58] studies of horizontal side

feeding tubes are common in the literature. A central, downward-facing vertical feeding tube
is utilized in our reactor to force gas and cellulose particles deep into the bed for catalytic
fast pyrolysis. The use of central feeding tubes is less common in the literature, and thus

the tubes’ effects on the hydrodynamics of the bed is less predictable.

1.5 Deviations from Canonical Beds

Even if no gas is injected through the feeding tube, the presence of the tube will alter
bubble formation and bed fluidization due to the obstruction caused by the tube. Additional
complications of the flow can result from high flow rate through the tube.

4
1.5.1 Horizontal Tubes

Immersed horizontal tubes within fluidized beds have been shown empirically and numer-
ically to alter fluidization and bubbling behavior [56, 35, 36, 5, 44]. Hull et al. [35] performed
numerical simulations to determine the effects of horizontal tubes on the average bubble size

and velocity. The authors compared their work to experimental studies [81] for verification.
The tubes broke the rising bubbles into smaller bubbles, which then coalesced back into
larger bubbles after rising above the tubes. The effect of the tubes is the retardation of

growth and speed of the bubbles, resulting ultimately in smaller, slower bubbles.
Further studies by Hull et al. [36] considered the effect on mixing times with and without
submerged horizontal tubes by tracking the mixing residence times of an added tracer gas

jet. The location of the tracer gas inlet was a significant determining factor of mixing time,
with a minimum mixing time occurring when the tracer gas was located in the center of
the bed. For slower fluidizing velocities, tracer gas inlets located on the fringes of the bed

led to mixing times on the order of two times longer than when the inlet was in the center
of the bed. The authors concluded that local minima in mixing residence times occurred
when both the tracer gas was injected centrally and when tubes were located either near the

distributor plate or near the top of the solids.

The horizontal immersed tube studies confirm that bubble behavior, and therefore mixing
times, are sensitive to not only the presence but also the location of obstructions relative to

gas inlets. Mixing is quickest near the center of the bed, indicating that the walls have a
strong influence in deterring mixing.

1.5.2 Vertical Tubes

Studies involving immersed vertical tubes and baffles have also shown that the presence of
such obstructions within fluidized beds hinders bubble growth, which helps prevent slugging

[63, 16, 57, 43, 66].

5
 Ozawa et al. [57] employed a neutron radiography method on a fluidized bed to visualize

and quantify void fractions for a bed with and without a vertical tube bank. The authors
considered Geldart Type A and B particles with ds = 62 µm and 218 µm, respectively. They
found that the presence of the tubes reduced the bubble diameters for both cases, and that

the tubes served as a dividing “guider” for the bubbles. Bubbles in a freely bubbling bed
“wandered” throughout the bed in their trajectory; with the tubes, however, bubbles were
“guided” vertically. Additionally, bubble diameter reduction was more prevalent in beds

with a lower height-to-depth ratio and the change in bubble reduction was attenuated at
higher ratios (2 ≤ H/D ≤ 3).
Our reactor includes one vertical tube through which gas is injected into the bed, and,

because the reactor is enclosed and insulated, its effect on fluidization, bubble size, and
mixing are unknown and require CFD simulations.

6
 CHAPTER 2

THEORY AND MODEL

Although fluidized beds have been used in industry for decades, only recently, with

advancements in computing power, has computational fluid dynamics (CFD) come to fuller
fruition in predicting the behavior of fluidized beds. In the past, fluidized beds could only
be modeled semi-empirically, theoretically, and with coarse-mesh CFD simulations. The

behavior of such beds is dependent upon the fully closed Navier-Stokes equations. However,
these exact equations cannot be solved without closure models to account for particle stress
and strain, gas drag, and friction between particles. Developments and improvements in

computing power in recent years have led to an extensive number of verified numerical
studies of fluidized beds. These simulations and the methods employed are now accepted
and commonplace for predicting bed behavior for a fluidized bed in standard configuration.

2.1 Governing Equations

c Fluent are derived
The governing two-phase conservation equations used in ANSYS

by Ishii [37] and applied to gas-solids flows by Enwald et al. [24]. We utilize the versions of
these equations as reproduced by van Wachem et al. [75]. A summary of all the governing
equations and closure models employed in this study is presented in the Appendix in Tables

A.1, A.2, and A.3. All symbol and variable definitions can also be found in the List of
Symbols. Further information about all models available in Fluent can be found in the
Fluent Theory Guide [3].

7
 The continuity equation for the gas phase g with velocity vg is

∂αg
+ ∇ · (αg vg ) = 0; (2.1)
∂t

for the solid phase s with velocity vs , the continuity equation is

∂αs
+ ∇ · (αs vs ) = 0. (2.2)
∂t

Here, αg represents the gas volume fraction and αs represents the solids (catalyst) volume
fraction. By conservation of mass, αs + αg = 1.
The two-phase gas/solids conservation of momentum equations are given by Enwald et

al. [24], where

 
∂vg
ρg αg + vg · ∇vg = −αg ∇P + ∇ · αg τ̄¯g + αg ρg g − β(vg − vs ) (2.3)
∂t

is for the gas phase, and

 
∂vs
ρs αs + vs · ∇vs = −αs ∇P + ∇ · τ̄¯s − ∇Ps + αs ρs s + β(vg − vs ) (2.4)
∂t

is for the solids phase. Here, τ̄¯k is the viscous stress tensor for phase k, β is the interphase
momentum transfer coefficient, P is the pressure shared by both phases, and Ps is the solids

pressure due to particle-particle interactions described by kinetic theory.

The conservation of energy equation for phase k involves the specific enthalpy hk and is
given in the Fluent Theory Guide [3] as

∂ ∂Pk
(αk ρk hk ) + ∇ · (αk ρk vk hk ) = αk + τ̄¯k : ∇vk − ∇ · qk , (2.5)
∂t ∂t

8
where qk is the heat flux of phase k. For our simulations, there are no chemical reactions or
heat source terms and thus qk = 0.

2.2 Closure Models

Closure of the momentum equations requires a description of the stresses on the particle.
Chapman and Cowling [13] developed the kinetic theory of granular flow (KTGF) and this
theory was furthered by Jenkins and Savage [39] and Lun et al. [52]. For our purposes, we

list only the closure models relevant to our simulation of a densely packed bed. Additional
information on other models can be found in the comparative studies of van Wachem et al.
[75] and the Appendix.

2.2.1 Granular Temperature

Particle velocity variations are measured by the granular temperature Θ, where

1 D 02 E
Θ= v . (2.6)
3 s

The granular temperature for the solid phase is analogous to the thermodynamic temperature

for the gas phase.

The energy balance for granular energy is

 
3 ∂  
(αs ρs Θ) + ∇ · (αs ρs Θvs ) = −Ps I¯ + τ̄¯s : ∇vs + ∇ · (κs ∇Θ) − γs − Js . (2.7)
2 ∂t

Here, κs is the solids thermal conductivity, γs is the dissipation of granular energy from

inelastic particle collisions, and Js is the dissipation of granular energy due to velocity fluc-
tuations.

9
 The dissipation term is given by Lun et al. [52] as

3
αs2 ρs g0 Θ 2
γs = 12(1 − e2 ) √ . (2.8)
ds π

The coefficient of restitution e is usually close to 1 (e = 0.9 in our case); ds and g0 are the
particle diameter and the radial distribution function, respectively.
Dissipation from velocity fluctuations is given as Js = β(vs0 · vg0 − vg0 · vs0 ). Closure of
these terms is given by Louge et al. [50], where

βds (vg − vs )2
 
Js = β 3Θ − √ . (2.9)
4αs ρs πΘ

Following the work of Syamlal et al. [71] and Boemer et al. [9], the energy balance
for dense flows can be simplified by assuming that the granular energy is in steady state
∂
( ∂t (αs ρsΘ) = 0), the granular energy is dissipated locally (Js = 0), and the convection
∇·(αs ρs Θvs ) and diffusion terms ∇·(κs ∇Θ) can be neglected. The energy balance simplifies
to
 
0 = −Ps I¯ + τ̄¯s : ∇vs − γs . (2.10)

Here, the first term represents the generation of fluctuating energy due to solids shear. For
dense flows, the remaining generation and dissipation terms dominate, so Equation (2.10) is
a reasonable estimate [75]. Solving Equation (2.10) for the granular temperature gives,

¯ )
−(K1 αs + ρs )tr(D̄ s
Θ=
2αs K4
q (2.11)
¯ ) + 4K α [2K tr(D̄
(K1 αs + ρs )2 tr2 (D̄ ¯ )]
¯ 2 ) + K tr2 (D̄
s 4 s 3 s 2 s
+ ,
2αs K4

10
with the abbreviations

K1 = 2(1 + e)ρs g0
4 2
K2 = √ ds ρs (1 + e)αs g0 − K3
3 π 3
 √    (2.12)
ds ρs π 2 8αs
K3 = 1 + (1 + e)(3e − 1)αs g0 + √ g0 (1 + e)
2 3(3 − e) 5 5 π
2
12(1 − e )ρs g0
K4 = √ .
ds π

When the flow is dilute, Equations (2.10) and (2.11) do not suffice and Equation (2.7)
must be solved. The average solids volume fraction within our bed is between 60−100% the
maximum catalyst packing value. Since the bed is primarily dense rather than dilute, we

employ the algebraic granular temperature formulation (Equation (2.11)).

2.2.2 Solids Pressure

The solids pressure Ps , given by Lun et al. [52], incorporates the granular temperature
into the momentum equations with

Ps = αs ρs Θ + 2g0 αs2 ρs Θ(1 + e). (2.13)

The two terms in the solids pressure represent the respective kinetic and collisional contri-

butions due to particle-particle interactions.

2.2.3 Viscous Stress

The viscous stress tensor for phase k is

 
¯ 2 ¯ )I,
¯
τ̄¯k = 2µk D̄k + λk − µk tr(D̄k (2.14)
3

11
where
¯ = 1 [∇v + (∇v )> ].
D̄ (2.15)
k k k
2
The solids bulk λs and shear µs viscosities must be modeled. Lun et al. [52] gives the solids

phase bulk viscosity as r

4 Θ
λs = αs ρs ds g0 (1 + e) . (2.16)
3 π
The solids shear viscosity is divided into two parts, the collisional contribution µs,col
and the kinetic contribution µs,kin , which together yield µs = µs,col + µs,kin . According to
Gidaspow [30], the collisional part is given as

r
4 Θ
µs,col = αs2 ρs ds g0 (1 + e) , (2.17)
5 π

and the kinetic part,

1√ 1√ 10 √ ρs ds
µs,kin = πΘρs ds g0 (1 + e)αs2 + πΘρs ds αs + πΘ . (2.18)
15 6 96 (1 + e)g0

2.2.4 Solids Thermal Conductivity

The solids thermal conductivity κs is considered when solving the full PDE for Θ. The
Gidaspow model [30] is utilized when relevant:

 2 r
2 6 2 Θ
κs = 1 + (1 + e)g0 αs κdil + 2αs ρs ds g0 (1 + e) , (2.19)
(1 + e)g0 5 π

where the dilute phase granular conductivity is given by

75 √
κdil = ρs ds πΘ. (2.20)
384

12
2.2.5 Frictional Stress

Frictional stresses occur at high solids volume fractions due to prolonged contact between
¯f is written as
particles. In Newtonian form, the frictional stress tensor σ̄

¯f = Pf I¯ + µf ∇v + (∇v)> .
 
σ̄ (2.21)

Frictional stresses are only considered when the solids volume fraction packs above the

critical point αs,min . The stresses are added to the solids pressure and solids bulk viscosity,
as described by kinetic theory in Equations (2.13), (2.17), and (2.18), such that

Ps = Pkinetic + Pf (2.22)

and

µs = µkinetic + µf . (2.23)

The two most commonly used frictional stress models are those of Johnson and Jackson
[40] and Schaeffer et al. [67] (furthered developed and implemented by Syamlal et al. [71]).
We use the first of these models, in which the normal frictional stress Pf is modeled semi-

empirically as
(αs − αs,min )n
Pf = F r , (2.24)
(αs,max − αs )p
where F r, n, and p are empirically determined material constants and αs,max is the maximum
packing solids volume fraction (typically between 0.63−0.65). Using the linear law proposed
by Coulomb [17], the frictional shear viscosity is related to Pf by

µf = Pf sin φ, (2.25)

where φ is the particle’s angle of internal friction. The empirical constants within this
frictional stress model have not been determined for a wide range of particles. In addition,

13
although Johnson and Jackson [40] only give values for glass and polystyrene beads (ds =
1800 µm and 1000 µm, respectively), the constants were not robustly determined, and as a

result, van Wachem et al. [75] recommend that frictional stress models be used with caution.
We use the Johnson and Jackson [40] model, with typical values F r = 0.05, n = 2, p = 5,

φ = 25◦ , and αs,min = 0.6 [75, 4]. Instabilities tend to occur when αs increased well above
αs,min and so we set αs,min closer to the packing limit of 0.63. Friction played the largest
role when vf l < vmf , i.e., when the bed is very densely packed. For this reason, we use the

frictional model to better match the experimental behavior of the bed below vmf .

2.2.6 Radial Distribution Model

The radial distribution function g0 increases the probability of particle collisions when

αs approaches the packing limit. The radial distribution function tends towards infinity at
maximum packing and unity when the solids volume fraction approaches zero [75]. The Lun
and Savage [51] radial distribution model is given as

 −2.5αs,max
αs
g0 = 1− . (2.26)
αs,max

2.2.7 Drag Model

The Gidaspow drag model [30] for the interphase momentum transfer coefficient β is

used for densely packed beds. It combines the Ergun equation [25] and the Wen and Yu
drag model [79], such that


2
150 αs µ2g + 7 ρg αs |vg −vs |


αg d 4 ds
for αs ≥ 0.2
s
β= (2.27)
 3 Cd ρg αs αg |vg −vs | αg−2.65


4 ds
for αs < 0.2.

Rowe [65] relates the drag coefficient CD to the relative Reynolds number Re with

14
 
 24 [1 + 0.15(Re)0.687 ] for Re < 1000


Re
CD = (2.28)

0.44
 for Re ≥ 1000,

where
αg ρg |vg − vs |ds
Re = . (2.29)
µg

2.3 Semi-empirical correlations

Before the prevalence of CFD models, correlations based on theory and empirical studies
were used to predict fluidization in fluidized beds. These include correlations for the mini-
mum fluidization velocity, the minimum bubbling velocity, and bed expansion. The empirical

correlations discussed below provide a basis of comparison for CFD simulations, especially
where experimental data is not available. It should be noted that all of the empirical cor-
relations are derived for a standard cylindrical fluidized bed configuration with a uniform

fluidizing gas flow, without the central feeding tube employed in our reactor. For this reason,
we first run simulations to validate our numerical model with a standard-geometry, tubeless
fluidized bed.

2.3.1 Minimum fluidization velocity

The minimum fluidization velocity equation originates from Stokes [70] for viscous resis-
tance of a spherical particle of diameter ds in a fluid with viscosity µg and velocity v. The

viscous drag is
Fv = 3πµg vds . (2.30)

πd3s
Setting the force equal to the effective gravitational force, 6
(ρs − ρg )g, we solve for the
terminal falling velocity of a single spherical particle,

15
 d2s (ρs − ρg )g
v0 = . (2.31)
18µg

Following the same logic for a single particle’s terminal velocity, Robinson [64] suggested
that a collection of uniform particles in a fluidized bed could be suspended with a minimum

fluidization velocity of
d2 (ρs − ρg )
vmf = K g, (2.32)
µg
where K is an empirically determined constant. Typically, subsequent correlations for vmf

follow the form of Equation (2.32).

Davies and Richardson [21] averaged experimental data for various catalysts (55 ≤ ds ≤
142 µm and 945 ≤ ρs ≤ 1230 kg/m3 ) to obtain a value for K. The Davies and Richardson

[21] correlation for minimum fluidization velocity is

d2 (ρs − ρg )
vmf = 0.0078 g. (2.33)
µg

Another correlation that is often used for fine Geldart Type A particles was derived by
Baeyens [6]. For 48 gas/solid systems with 20 ≤ ds ≤ 70 µm and 1117 ≤ ρs ≤ 3920 kg/m3 ,

Abrahamsen and Geldart [1] confirmed the correlation as the best-fit equation for vmf ,

0.0009(ρs − ρg )0.934 d1.8 g 0.934

vmf = . (2.34)
µ0.87
g ρg
0.066

Many other authors [55, 65, 79, 27, 62, 12] have similar correlations, based on different
experimental data. A discussion and summary of various vmf correlations can be found in

[80] in which Wu et al. study the effect of temperature on the minimum fluidization velocity.
The gas viscosity and density are both functions of temperature, and thus vmf can vary
dramatically depending on the operating temperature.

16
 Our research operates with particles most similar to those of Davies and Richardson [21],

Baeyens [6], and Abrahamsen and Geldart [1] – catalysts with ds = 100 µm and ρs = 1750
kg/m3 . To fully understand the fluidization process, we compare these correlation models,
our experimental data, and the computational data generated.

2.3.2 Minimum Bubbling Velocity

The minimum bubbling velocity is defined as the velocity at which bubbles begin to form
within a fluidized bed [29]. The traditional equation for vmb applied only for air at ambient

conditions and was given by Geldart as

vmb = 100ds . (2.35)

Abrahamsen and Geldart [1] updated this equation to account for gases with different den-
sities and viscosities. Based on their data, they found the correlation

ds ρg0.06
vmb = 2.07 exp (0.716F45 ) , (2.36)
µ0.347

where F45 is the fraction of particles less than 45 µm.

2.3.3 Bed expansion

Bed expansion with respect to vmf and vmb is empirically derived by Abrahamsen and
Geldart [1] as
 0.22
Hmb vmb
= . (2.37)
Hmf vmf
The minimum fluidization velocity, the minimum bubbling velocity, and the bed expan-

sion correlations provide comparison points for our simulations of a standard configuration
fluidized bed.

17
 CHAPTER 3

SIMULATION

3.1 Modeling methods

The computational fluid dynamics (CFD) of a fluidized bed are typically simulated using
either the Lagrangian or the Eulerian method. The Lagrangian method applies Newtonian
physics to each particle, tracking its position, velocity, momentum, and how particles collide

and interact with each other and the walls. The Eulerian method, on the other hand, is a
control-volume approach, which treats each phase as interpenetrating continua [30].

3.1.1 Lagrangian Method

Examples of Lagrangian simulations include Hoomans et al.’s simulations of 40,000 par-

ticles treated as hard spheres [32], Bokkers et al.’s numerical and experimental study with
30,000 particles [10], Chu and Yu’s particle-fluid flows with 3000−20,000 particles [14], and
more recently, Wang et al. simulations with 42,000 particles [78]. For further information on

Lagrangian method simulations, see reviews written by Deen et al. [22] and Zhu et al. [82].
The discrete particle method (DPM) is the state-of-the-art, Lagrangian particle-tracking
method to discretely account for all but interparticle cohesive forces each particle encounters.

The method is widely accepted as an accurate way to simulate fluidized beds [22, 83, 82], but
it is severely limited due to the computational expenses involved in tracking a large number
of particles.

Most discrete tracking simulations only consider beds with up to 50,000 particles. Our
lab-scale reactor uses 90 grams of ZSM5 catalyst particles, which is on the order of 107

18
particles, a prohibitively expensive number of particles to track. A more practical approach

is to instead use the faster Eulerian method.

3.1.2 Eulerian Method

The Eulerian method discretizes the problem into a mesh consisting of cells within which

conservation equations (mass, momentum, and energy) are applied. The method requires
less computing power and is therefore more frequently applied for solving dense, multiphase
flows. Because both the gas and solids phase are treated as liquids, this method is often

referred to as an Euler-Euler (or Euler-Granular, due to the granular modeling of the solids
phase) approach, as opposed to an Euler-Lagrange approach in which a control volume is
used for the gas but not the particles.

In theory, the finer the mesh, the more accurately the Eulerian approach can describe
multiphase flow. Often, mesh refinement studies are included within the literature [53, 48, 33]
to ensure that for a given mesh size, a finer mesh still tends towards the same solution.

Inaccurate bed expansion values for unrefined grids is a primary problem for the simulation
of fine particles. In the following section we discuss studies on the limitations of the Eulerian
approach when applied to fine particles.

3.2 Simulation of Geldart Type A Particles

Good agreement with experimental data has been found for Geldart Type B and D par-

ticles for numerical simulations using the Eulerian approach [60, 34, 72]. However, numerical
simulations of Type A particles continue to challenge scientists and engineers. In the past,
simulations tended toward unrealistic solutions with overestimation of bed expansion on the

order of two times the experimental bed height [26, 2, 53, 84]. Based on the work of Parmen-
tier et al. [59] and Di Renzo and De Maio [23], Wang et al. [76] noted that the cause of this
overestimation may be due to insufficient resolution of subgrid structures. These structures,

19
such as small bubbles, could be resolved by decreasing the mesh size such that cell lengths

δ are 2−4 times the particle diameter [76].

Another theory on expansion overprediction is that structures begin to cluster due to
interparticle forces (IPFs) [28, 53, 84, 74]. Mckeen and Pugsley [53] conducted numerical

studies with fluid catalytic cracking (FCC) particles (ds = 75 µm) and found that by scaling
down the drag coefficient by a factor of about 0.2−0.3, their simulations were able to match
experimental data. They argued that the particle clusters at a subgrid level necessitate this

scaling factor. The agglomerate cluster size they find to be ∼135−170 µm, about two times
the actual particle size.
Following the work of Mckeen and Pugsley [53], other authors validated Geldart Type

A simulations with experimental data through the use of drag scaling factors [19, 33, 48,
47]. Lindborg et al. [48] created a varying scaling factor specifically designed to match
the pressure vs. fluidizing velocity measured experimentally in their laboratory. Similarly,

Hosseini et al. [33] found an ideal drag scaling factor of 0.1 in order to find agreement
with experimental bed expansion, but they note that the scaling factor varies with fluidizing
velocity and a general scaling factor model has not been fully established for Type A particles.

It has been shown that for Type A particles of size ds = 50−105 µm, IPFs are negligible
since the fluidization characteristics can be fully described with frictional stresses [49, 69,
73]. To validate this, Wang et al. [78] performed simulations using the DPM to explore

IPFs at different particle diameters (50 ≤ ds ≤ 150 µm and ρs = 1500 kg/m3 ). They
found that for Type A particles larger than 55 µm, IPFs are essentially negligible, and
a refined mesh resolves the subgrid bubbling structures needed to properly predict bed

expansion. The authors stressed, however, that their studies only involved cold flow under
atmospheric pressure, and that the role of inter-particle cohesive forces could also be a
function of temperature.

20
 Wang et. al [77] also performed an Eulerian-Eulerian CFD study to determine the effect

of mesh size on determining vmb and compared their results to the semi-empirical vmb model
derived by Abrahamsen & Geldart [1] (Equation (2.36)). The authors found that with a
course mesh (δ = 10 × ds ), vmb was overpredicted by over 77% times the literature value due

to lack of resolving any bubbles. When the mesh size was changed to δ = 2.7 × ds , vmb was
determined to be within 20% times the empirically determined vmb . Wang et. al considered
there to be no bubble phase in a simulation if within every cell in the catalyst bed, αs > 0.2.

Conversely, the authors declared a simulation to have a bubble if within any cell in the bed
αs < 0.15. We use a similar criteria for defining the bubble phase and finding vmb .
In recent years, the advancement of computing capability makes meshes of minute size

realizable. Here, we use a mesh resolution of ∼ 30 cells/cm. This equates to about 3.4 times
the particle diameter for ds = 100 µm, which fits well within the 2−4 range for sufficient
resolution.

3.3 Physical Parameters

In our simulations, we approximate the ZSM5 catalysts to have a mean particle diameter
of ds = 100 µm, a density of ρs = 1750 kg/m3 , and the total catalyst bed to have a mass of
m = 90 g. Since the actual distribution of particle sizes in the reactor is bimodal, we can

expect some discrepancy between the simulations and the experimental reactor.
The fluidizing gas properties, in particular ρg and µg , vary according to the operating
temperature. Due to the high temperature necessary for CFP, the experimental reactor is

operated at T = 773 K. Table 3.1 shows the effect of temperature on the density (based on
the ideal gas law) and dynamic viscosity of helium and the argon and nitrogen tracer gases
[61, 42, 11].

The minimum fluidization velocity is also a function of temperature as well (see Equations
(2.32), (2.33), and (2.34)). Since the difference in densities (ρs − ρg ) is dominated by the

21
 T ρHe ρAr ρN2 µHe (10−5 µAr (10−5 µN2 (10−5
(K) (kg/m3 ) (kg/m3 ) (kg/m3 ) kg/(m·s)) kg/(m·s)) kg/(m·s))
300 0.163 1.623 1.138 1.991 2.261 1.788
400 0.122 1.217 0.854 2.436 2.883 2.210
500 0.098 0.974 0.683 2.847 3.426 2.579
600 0.081 0.811 0.569 3.235 3.914 2.911
773 0.063 0.630 0.442 3.863 4.676 3.418
800 0.061 0.609 0.427 3.956 4.779 3.494

Table 3.1: Helium, argon, and nitrogen density and dynamic viscosity vs. temperature.

solids density, which is not a function of temperature, the primary parameter in evaluating
the dependence of vmf on temperature is the gas viscosity. Between room temperature (300
K) and the operating conditions (773 K), the viscosity of helium gas increases by a factor

of 1.94. Thus we can expect the minimum fluidization velocity to be inversely proportional
to this ratio, i.e., 0.51 times less. Experimentally, the gas is initially at room temperature
when it begins to proceed down the tube. However, the gas is rapidly heated by the inside

wall of the feeding tube. We assume that by the time the gas proceeds down the tube, it
has reached the full 773 K, and as such isothermal conditions are met.

3.4 Geometry
The geometry of our simulated model is based on the experimental reactor, which is
60.96 cm (24 in) tall and has a 5.08 cm (2 in) outer diameter, a 4.925 cm (1.939 in) internal

diameter, and a 1.27 cm (0.5 in) diameter feeding tube down the center axis of the reactor
extending to a depth of 1.524 cm (0.6 in) above the distributor plate. Due to computational
expenses, we simulate a height of H = 12.7 cm (5 in) with a radius of R = 2.54 cm (1 in) on

a 2D axisymmetric mesh. Because the flow becomes constant above the catalyst bed and the
catalyst bed height is only about 4−6 cm, the height restriction is a reasonable assumption.
A mesh resolution of 29.528 cells/cm (∼ 30 cells/cm = 75 cells/R) is used, which introduces

22
some slight rounding on the tube locations. The feed tube is 19 cells wide with a radius of

RT ≈ 0.643 cm and extends 46 cells above the bottom fluidizer inlet: HT ≈ 1.558 cm. The
meshes with and without the feeding tube can be seen in Figure 3.1. In the figure, the mesh
is mirrored over the axis of symmetry. The meshes are also described in Table 3.2.

Figure 3.1: CFD meshes at ∼ 30 cells/cm mesh resolution for meshes without a tube (left)
and with a tube (right). Gas inlets, the feeding tube, and the axis of symmetry are indicated.

23
 Cells/Radius Grid Size (µm) δ/ds Total Cells Tube
75 338.67 3.39 28125 No
75 338.67 3.39 21938 Yes

Table 3.2: Summary of 2D axisymmetric meshes of the reactor with and without the tube.

3.5 Boundary Conditions

c Fluent 14.0 is utilized for the simulations. The
The commercial software ANSYS

following subsections detail the boundary conditions for each wall, inlet, and outlet used for
these simulations. Information regarding the boundary conditions and options available in
Fluent is available in the Fluent User’s Guide [4].

3.5.1 Walls

We consider one (two) wall(s) for the mesh without (with) a feed tube. Each wall is set
to the boundary type wall. For both walls, a free-slip boundary condition is prescribed for

the two phases due to the uncertainty of boundary layers in such a two-phase flow. The free-
slip condition is applied by setting the Specified Shear for both the x- and y-components
to a constant 0 Pa. The walls are set to be stationary with temperature = 773 K, wall

thickness = 0 m, and heat generation = 0 W/m3 . One can define the material of the walls
by creating a steel material from the Fluent database of materials. However, since heat-
transfer is insignificant in an isothermal reactor, the material of the walls does not matter

in these simulations. Lastly, for the two tracer gases, the wall boundary condition is set to
Zero Diffusive Flux.

3.5.2 Axis

The internal axis of symmetry is defined by selecting the type axis. No further specifi-
cation is needed beyond selecting both 2D in the Fluent start-up menu and Axisymmetric

24
in the general options. When a 2D axisymmetric setting is used, derivatives with respect

to φ are set equal to 0, where in a cylindrical coordinate system (r, z, φ), φ is the angular
coordinate, r is the radial coordinate, and z is the axial coordinate.

3.5.3 Inlets

Each of the inlets is set to type velocity-inlet for both phases. This ensures Dirichlet
boundary conditions [15], a constant inlet velocity (non-zero value for gas, zero value for
solids), which we prescribe in the axial direction. The inlet solids volume fraction is set to 0

as well, and the inlet granular temperature is set to 0.001 m2 /s2 to account for minor velocity
fluctuations. For the mixture of the two phases, the Supersonic/Initial Gauge Pressure
is set to 0 Pa.

For simulations without a feed tube, there is only one inlet to consider. The fluidizer
inlet gas velocities are set to one of the following, depending on the simulation: vf l = 0.1,
0.2, 0.3, 0.4, 0.6, 0.8, 0.9, 1.0, 1.5, 2.0, or 3.0 cm/s.

Two sets of simulations are performed with the feeding tube based on flow rate configura-
tions used with the experimental reactor, the results of which are in development for a paper
[41]. In the first set, five simulations were run to study the effect of gas distribution between

the feed tube and the bottom fluidizer. This was achieved by changing the percentage of
the flow from the feed tube Qf t (Qf t /Qtotal ≈ 8%, 17%, 33%, 51%, or 68%) while keeping
the total gas flow rate constant (Qtotal ≈ 55 cm3 /s), where Qtotal = Qf t + Qf l . For this set

of simulations, Qtotal ≈ 8 × Qmf , where Qmf is the fluidizer flow rate at the experimentally
determined minimum fluidization velocity (Qmf = 6.917 cm3 /s, vmf = 0.363 cm/s). Ex-
perimentally, the best aromatic yield was found when Qf t /Qtotal = 34%. The second set of

simulations are used to study the effect of changing the total flow rate (Qtotal /Qmf ≈ 3, 4,
6, or 8) while holding constant the ratio Qf t /Qtotal ≈ 34%.

25
      
cm

cm3 cm3 cm3 Qtotal cm

Qf t
vf l s
Qf t s
Qf l s
Qtotal s Qmf
vf t s Qtotal
(%)
0.1 - 1.91 1.91 0.28 - -
0.2 - 3.81 3.81 0.55 - -
0.3 - 5.72 5.72 0.83 - -
0.4 - 7.62 7.62 1.10 - -
0.6 - 11.43 11.43 1.65 - -
0.8 - 15.24 15.24 2.20 - -
0.9 - 17.15 17.15 2.48 - -
1 - 19.05 19.05 2.75 - -
1.5 - 28.58 28.58 4.13 - -
2 - 38.10 38.10 5.51 - -
3 - 57.15 57.15 8.26 - -
     
Qf t cm3 cm3 cm3 Qtotal cm

cm


Qtotal
(%) Qf t s
Qf l s
Qtotal s Qmf
vf t s
vf l s
7.74 4.33 51.63 55.96 8.09 3.33 2.55
17.48 9.75 46.01 55.76 8.06 7.49 2.27
33.21 18.41 37.04 55.45 8.02 14.16 1.83
51.11 28.16 26.94 55.10 7.97 21.65 1.33
68.24 37.37 17.40 54.77 7.92 28.73 0.86
     
Qtotal cm3 cm3 cm3 Qf t cm

cm


Qmf
Qf t s
Qf l s
Qtotal s Qtotal
(%) vf t s
vf l s
3.15 7.24 14.57 21.82 33.21 5.57 0.72
4.21 9.66 19.43 29.09 33.21 7.43 0.96
6.31 14.49 29.14 43.63 33.21 11.14 1.44
8.41 19.32 38.86 58.18 33.21 14.85 1.92

Table 3.3: Flow rate configurations for 3 sets of CFD simulations: a set without a tube in
which the fluidizer velocity is varied, a set of runs with the tube and a near-constant Qtotal
as a function of Qf t , and a set of tube simulations with a constant Qf t /Qtotal as a function
of Qtotal .

26
 The experimental gas velocities were used for the simulations, but as these were calculated

using the experimental internal diameter (4.925 cm) rather than the 5.08 cm diameter used in
simulation, upon back-calculation of the flow rates, we find that the flow rate ratios deviate
slightly from the experimental runs. This introduces a relatively small error, which bears no

difference in the interpretation of the results. As noted in the summary of flow conditions
in Table 3.3, Qtotal differs at most by 2.17% in the first set of runs.

3.5.4 Outlet

The outlet is set to pressure-inlet, in which a constant pressure (0 Pa above atmo-

spheric) is prescribed and the gas and solids velocities compensate appropriately (known

as a Neumann boundary condition [15]). Due to the constant flow rates at the top of the
reactor, we rarely observed backflow occurring in the simulations. The backflow conditions
are as follows: αs = 0, T = 773 K, YAr = 0, YN2 = 0, and Θ = 0.0001 m2 /s2 , where YAr and

YN2 are the respective argon and nitrogen mass fractions.

3.5.5 Internal field

The initial uniform solids volume fraction in the interior of the bed is determined by the
height of the bed and the initial mass of zeolite catalysts. For a bed mass m, the following

equation must be satisfied within a standard (i.e. no tube) fluidized bed configuration:

m
αs = , (3.1)
πR2 hρs

where h is the catalyst bed height. A maximum-packed bed occurs when αs = αs,max = 0.63
and corresponds to h = 4.242 cm in our case. Above this bed height, αs = 0 and thus, αg = 1.
A loosely packed bed (αs = 0.44) is commonly used as an initial condition corresponding to

h = 5.8 cm within a tubeless bed. When the tube is considered, the bed height for a given
mass will be higher due to the volume occupied by the tube.

27
 The initial solution for the tube simulations comes from a simulation by Saba Almalkie,

a previous researcher for the EFRI group. Due to differences in these simulations, slight
changes had to be made such as the tube size (from RT = 0.457 cm to 0.643 cm), the slip
condition on the walls (from no-slip to free-slip), the frictional model (from none to the

model described in Section 2.2.5), the particle diameter (from ds = 50 µm to 100 µm), the
mesh size (from ∼ 39 cells/cm to ∼ 30 cells/cm), and the operating temperature (from 900
to 773 K). These changes were primarily made by mapping the solution from the older mesh

to the newer. The internal temperature was adjusted by patching a new temperature into
each cell, and subsequently, the properties of the gas were changed according to Table 3.1.
Lastly, because we mapped the previous solution onto a mesh with a larger feed tube, we

had to add the appropriate amount of mass back to the system to conserve the 90 grams of
catalyst. This mass was added well above the bed and was allowed to fall into the bed. After
the new catalyst mixed with the rest of the bed, gas velocities were changed according to

the tube configurations in Table 3.3. This method of solution mapping was utilized due to
difficulties in initializing the solution caused by instabilities from the frictional stress model.
The internal initial conditions for the cases without the tube were set to an operating

temperature of 773 K. All these cases were run using the vf l = 1.0 cm/s simulation as the
initial solution.

3.6 Numerical Method

c Fluent gives a wide range of selection for solution methods, numerical schemes,
ANSYS

and solvers. These were selected based on simulation consistency, accuracy, and stability.
The multiphase model in Fluent uses the Semi-Implicit Method for Pressure-Linked Equa-
tions (SIMPLE) method. Further information for this solving method can be found in the

Fluent Theory and User’s Guides [3, 4] and in the CFD book by Chung [15]. We typically
set the number of SIMPLE iterations to 60−100 iterations per time-step, with time-steps

28
 p ρ Body Forces Momentum αs Θ Energy Tracer Gases
Tubeless runs 0.2 0.5 0.6 0.4 0.2 0.2 0.5 -
Tube runs 0.2 0.4 0.8 0.4 0.2 0.1 0.6 0.8

Table 3.4: Under-relaxation factors.

varying between 10−5 and 10−4 s. We used modified, smaller under-relaxation factors for
better stability. These are listed in Table 3.4.

To initialize the flow field for tubeless cases, the numerical solvers were first order implicit
in time and first order upwind for solving the momentum equation. After about 100 time-
steps, we switched to a second order implicit solver in time and a higher order QUICK

scheme. The QUICK scheme is especially designed for quadrilateral meshes and uses a
solution-dependent weighted average between the second-order upwind scheme and central
interpolation [3]. QUICK was also used for tube simulations and for solving all variables

(momentum, volume fraction, granular temperature, energy, and tracer gases). Additionally,
a least-squares, cell-based method was used for solving gradients. Generally, the time-step
was increased to 10−4 s after a quasi-steady state was reached unless simulation stability

required us to reduce the time-step.

The simulations were run on 4−32 processors each, and took, on average, 200 h to
simulate 10 s of simulation time. We noted that the benefits of parallel processing tapered

off after about 8 processors due to the relatively low number of total cells (∼ 28, 000 cells).

3.7 Quantifying Mixing

Once a fluidized bed has reached a quasi-steady state, we study its long-term behavior
and quantify the fluctuating quantities by using time-averaged statistics such as the mean
and variance. A quasi-steady state is reached after about 2−4 s of simulation time. This

29
is determined by performing statistical analyses on different simulation durations until the

results are statistically stable. These are discussed furthermore in Chapter 4.

3.7.1 Tracer Gases

Using Fluent, tracer gases are introduced at each inlet to track mixing. The gas phase

is primarily helium (He) but at the fluidizer inlet, a nitrogen (N2 ) gas is introduced at 2%
mole fraction and at the feed tube inlet, an argon (Ar) gas is introduced, also at 2% mole
fraction. The mole fraction Xi of species i can be converted to a mass fraction Yi via the

following relationship:
M
Yi = Xi × , (3.2)
Mi
where M and Mi are the average molecular weights of the mixture and the species i, respec-
tively. It follows that within each cell, YHe + YAr + YN2 = 1.

The free-stream mass fraction is pertinent to the inlet from which the tracer gas is
introduced and represents the maximum possible mass fraction of that tracer gas. We
∞
calculate the free-stream mass fractions of Ar and N2 , defined as YAr and YN∞2 , respectively,

∞ M (4.0026 + 39.948) kg/kgmol

YAr = XAr × = 0.02 × = 0.169217
MAr 39.948 kg/kgmol
(3.3)
M (4.0026 + 28.0134) kg/kgmol
YN∞2 = XN2 × = 0.02 × = 0.124981.
MN2 28.0134 kg/kgmol

In addition, the free-stream mass fractions of each tracer gas are the same for every case,

regardless of the velocity at the inlets.

The well-mixed case occurs when the tracer gases and the fluidizing helium gas mix
proportional to the volume or mass flow rate of each component. The outlet mass fractions

of each tracer gas represents such a well-mixed case, and we denote these as the ideal mass

30
 ideal
fractions YAr , YNideal
2
ideal
, and YHe . We use YN2 to demonstrate how YNideal
2
can be found. An

alternate form of the mass fraction can be defined in terms of mass flow rates:

ṁoutlet
N2
YNideal = YNoutlet =
2 2
ṁoutlet
mix
(3.4)
ṁfNl2
= ,
ṁfmix
t
+ ṁfmix
l

where, f t and f l are the respective feed tube and bottom fluidizer inlets. In order to calculate
the mass flow rate, we first find the density of the mixture at each inlet:

!−1 −1
YN∞2 Y fl YN∞2 1 − YN∞2

ρfmix
l
= + He = + (3.5)
ρN2 ρHe ρN2 ρHe
!−1 −1
∞ ft  ∞ ∞
YAr YHe YAr 1 − YAr
ρfmix
t
= + = + , (3.6)
ρAr ρHe ρAr ρHe

and then multiply each these by the respective volume flow rate to get the mass flow rates:

ṁfmix
l
= Qf l × ρfmix
l
(3.7)

ṁfmix
t
= Qf t × ρfmix
t
(3.8)

ṁfNl2 = YN∞2 × ṁfmix

l
. (3.9)

These mass flow rates depend upon the volume flow rates at the inlets (see Table 3.3)

and the gas densities (see Table 3.1). We tabulate these ideal mass fractions in the following
subsection (see Table 3.5) and discuss their significance in determining the mixture fraction
and the ideal mixture fraction.

31
3.7.2 Mixture Fraction

The mixture fraction ξ is our main measure for quantifying the mixing of the feed tube
tracer gas with the remainder of the bed. The mixture fraction is a conserved scalar typically
used in combustion problems in turbulent mixing cases [8, 3]. Here, there are no chemical

reactions occurring, simply diffusive mixing between the two tracer gases (Ar and N2 ) and
the dominant He fluidizing gas. We define the mixture fraction such that ξ = 1 when the
maximum possible concentration of the Ar feed tube tracer gas component is attained and

ξ = 0 when the component concentration is 0:

YN∞2 − YN2
ξ= , (3.10)
YN∞2

Thus, at the feeding tube inlet, where the argon concentration is at its maximum and the
nitrogen concentration is negligible, ξ = 1; at the fluidizer inlet, where YAr = 0, ξ = 0.

Note that in the well-mixed case, the mixture fraction will not be equal to 1. Combining
Equation (3.10) with Equations (3.4), (3.9), (3.8), and (3.7) yields the following equation for
the ideal mixture fraction, denoted with ξideal :

YN∞ ×ṁfmix
t

YN∞2 − 2
Qf t ×ρfmix
t
+Qf l ×ρfmix
l
ξideal =
YN∞2
Qf t × ρfmixt
= 1−
Qf t × ρfmix
t
+ Qf l × ρfmix
l
(3.11)
Qf l × ρfmix l
=
Qf t × ρfmix
t
+ Qf l × ρfmix l
!−1
Qf t × ρfmixt
= +1 .
Qf l × ρfmixl

ξideal and YN2 are tabulated in Table 3.5 for each simulation with a feed tube.

32
 Qf t
Qtotal
(%) YNideal
2
ξideal
7.74 0.115 0.081
17.48 0.102 0.182
33.21 0.082 0.344
51.11 0.059 0.524
68.24 0.038 0.693

Qtotal
Qmf
YNideal
2
ξideal
3.15 0.082 0.344
4.21 0.082 0.344
6.31 0.082 0.344
8.41 0.082 0.344

Table 3.5: Ideal N2 tracer gas mass fraction and mixture fraction for each simulation with
the feeding tube.

3.7.3 Volume-weighted Statistics

The following subsections describe the method of calculating the volume-weighted mean

and variance for an arbitrary, cell-centered quantity O(r, z, t). Spatial averages are denoted
D E
with an overbar O while time averages are denoted as O . Since we are simulating a 3D
fluidized bed using a 2D axisymmetric slice, simulation outputs are dependent only on r, z,

and t, the radial distance from the axis of symmetry to the cell center, the axial distance
from the fluidizer to the cell center, and the simulation time, respectively. Statistics on the
solids volume fraction and the mixture fraction will be performed as a function of bed height,

bed radius, and over the catalyst bed as a whole.

In order to inspect the mixing quality within the bed, we employ primarily 4 statistical
quantities applied over the entire space of the catalyst bed and over a statistically steady

time period (4 s). These quantities are scaled to show relevant information. The first is
D E
ε = αs /αs,max , the time-ensemble, volume-weighted mean within the catalyst bed of the
solids volume fraction, normalized by the maximum packing limit. Secondly, we analyze

33
 D E
02
ψ = αs02 /αs,max , the normalized variance of the catalyst. The maximum possible variance
02
αs,max is determined by considering the variance between two samples, one at the lowest
possible value (αs = 0) and the other at the highest possible value (αs = 0.63). It follows
02
that αs,max = (0.63 − 0)2 /2 = 0.19845. The parameter ψ shows how much the catalyst
changes inside the bed, and is a good measure of the amount of bubbling within the bed. A
slugging bed, for instance, will exhibit a high value of ψ as large voidages rise over the width

of the bed followed by the catalyst, filling the space in the wake.
Although ε and ψ are useful in probing the mixing of the catalyst, we are essentially
interested in how well the feed tube tracer gas (representing the gasified cellulose) mixes

throughout the bed. The time-ensemble, volume-weighted mean of the mixture fraction
D E
scaled by the ideal mixture fraction θ = ξ(r) /ξideal provides a good measure of how well
this tracer gas is mixed. A well-mixed bed should thus have θ ≈ 1.

The most important statistic we use in quantifying the mixing potential of the feed
D E D E2
tube gas is Φ = ξ 02 / ξ , the time-ensemble, volume-weighted variance of the mixture
fraction, which we scale by the mean mixture fraction squared. We use this scaling in order

to acquire a dimensionless quantity that signifies not just a low variance, well-mixed case
but also one with a significant amount of tracer gas present. The mixing parameter Φ is
equivalent to the squared coefficient of variation, where the coefficient of variation is the

standard deviation over the mean. In our case, the lower Φ is, the better we consider the
case for the mixing of the feed tube gas. It is this value that we compare to the aromatic yield
since stoichiometrically, better mixed cellulose gas and catalyst produces more aromatics.

This set of quantities {ε, ψ, θ, Φ} is also calculated vs. bed height (ε(z), ψ(z), θ(z), and
D E D E2
Φ(z)) and vs. bed radius (ε(r), ψ(r), θ(r), and Φ(r)). For instance, Φ(r) = ξ(r)02 / ξ(r) ,
D E D E2
Φ(z) = ξ(z)02 / ξ(z) , and so forth. Two other statistical quantities that are relevant
h−hmf
especially to the cases run without the tube are the bed expansion Ω = hmf
and the volume
Vαs ≤0.15
of the bed occupied by the bubble phase χ = Vbed
(see the bubble phase definition from

34
Section 3.2). Thus the set of quantities we analyze for simulations without a tube is {ε, ψ,

Ω, χ}

3.7.4 Statistical Convergence

Before calculating any time-averaged statistics, we must first gauge a statistically steady

time period (TS ) over which to perform our analyses. We look specifically at the convergence
of the four mixing parameters ε, ψ, θ, and Φ with respect to a range of sampling time lengths
(TS = 0.1, 0.2, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0 s) from the latest time for each

simulation. If, for a specific case the mixing parameters all have converged to a statistically
significant value, then we use the converged sample time for all statistical analyses. Because
some cases mix better than others, the sampling time can be shorter for some cases (TS ∼ 2

s) than for others (TS ∼ 4 s).

3.7.5 Volume-weighted Mean

We weigh each cell by its respective volume before averaging, since the 2D axisymmetric
assumption means the outer cells within the cylinder occupy a greater volume than the inner
cells. Horizontal slice averages convert the three-dimensional (two dimensions in space, one
D E
dimension in time) value O(r, z, t) into a one-dimensional vector as a function of z, O(z) .
For each z location, we integrate over the radius and time, dividing by the volume of the
slice and the time over which the average is taken, ti to tf :

Ztf ZR
D E 1 2π∆z
O(z) = O(r, z, t) r dr dt, (3.12)
tf − ti V (z)
ti 0

where the volume of the horizontal slice V (z) is dependent on whether or not z ≥ HT :

V (z) = π(R2 − RT2 (z))∆z, (3.13)

35
and 

0
 for z < HT (below feeding tube)
RT (z) = (3.14)

RT
 for z ≥ HT (above feeding tube).

Since we utilize discrete time-steps, the integral can be written as multiple summations.
The discrete time-ensemble, volume-weighted radial mean of a quantity O is therefore

N R
D E 1 2∆r X X
O(z) = O(r, z, t) r. (3.15)
N R2 − RT2 (z) i=1
r=RT (z)

Vertical slice averages do not require a strict volume-weighting since cells in the axial
direction have the same volume. We follow the formulation above to derive an equation
for the time-ensemble, volume-weighted axial mean. The averages must be taken from the

bottom of the bed to the height of the catalyst at each radial point, and so we must have
some knowledge of what the bed height h is at each radial point r and at each time t, denoted
with h(r, t). This bed height, is found by scanning from the top of the reactor down until

αs (r, z, t) > 0.2, at which point h(r, t) = z.

The time-ensemble, volume-weighted axial mean of O(r, z, t) is, then,

t h(r,t)
D E 2πr∆r Z f Z O(r, z, t)
O(r) = dz dt, (3.16)
tf − ti V (r, t)
ti 0

where
V (r, t) = 2πr∆rh(r, t). (3.17)

For discrete data with N samples in time, Equation (3.16) simplifies to

E 1 X N h(r,t)
D X O(r, z, t)
O(r) = . (3.18)
N i=1 z=0 h(r, t)

36
3.7.6 Volume-weighted Variance

The volume-weighted variance can be derived in similar fashion to the volume-weighted

mean, except using the square of the difference between the current value O(r, z, t) and the
D E D E
respective time-ensemble mean, O(z) or O(r) . The time-ensemble, volume-weighted
D E
radial variance of O(r, z, t) is denoted as O02 (z) and is given by

Ztf ZRh
D E 1 2π∆z D Ei2
O02 (z) = O(r, z, t) − O(z) r dr dt, (3.19)
tf − ti V (z)
ti 0

which simplifies in discrete form to

N R
D E 1 2∆r X X h D Ei2
O02 (z) = O(r, z, t) − O(z) r. (3.20)
N R2 − RT2 (z) i=1
r=RT (z)

For vertical slices, the time-ensemble, volume-weighted axial variance of O(r, z, t), de-
D E
noted as O02 (r) , is

t h(r,t)
h D Ei2
D E 2πr∆r Z f Z O(r, z, t) − O(r)
O02 (r) = dz dt, (3.21)
tf − ti V (r, t)
ti 0

which simplifies to

h D Ei2
D E 1 N h(r,t) O(r, z, t) − O(r)
X X
O02 (r) = . (3.22)
N i=1 z=0
h(r, t)

3.7.7 Bed Mean and Variance

In order to quantitatively compare cases with different inlet velocities, we average statis-

tics over the entirety of the bed. To calculate the bed average, we perform a time-ensemble,

37
volume-weighted average over all space within the bed and over all time such that

Ztf ZRT h(r,t)

2π O(r, z, t) r
D E Z
O = dz dr dt (3.23)
tf − ti Vbed (t)
ti 0 0

where the volume of the bed at each time t is

ZRT h(r,t)
Z
Vbed (t) = 2π r dz dr, (3.24)
0 0

or, as a summation,
RT h(r,t)
X X
Vbed (t) = 2π∆z∆r r. (3.25)
r=0 z=0

Equation (3.23) is converted discretely to give the time-ensemble, volume-weighted bed

mean of O(r, z, t):
N PRT Ph(r,t)
z=0 O(r, z, t) r
D E 1 X
r=0
O = PRT Ph(r,t) (3.26)
N i=1 r=0 z=0 r

The time-ensemble, volume-weighted bed variance of O(r, z, t), is similarly given as

Ztf ZRT h(r,t)

Z 
2
D E 2π O(r, z, t) − O r
O02 = dz dr dt, (3.27)
tf − ti Vbed (t)
ti 0 0

whose discrete form is

N PR T Ph(r,t) 
2
D E 1 X r=0 z=0 O(r, z, t) − O r
O02 = PRT Ph(r,t) . (3.28)
N i=1 r=0 z=0 r

The time-ensemble, volume-weighted variance within the catalyst bed, henceforward referred

to as simply the bed variance, is the main metric used for determining mixing within the
fluidized bed. The next chapter includes the simulation results and data analyses following
the time-ensemble, volume-weighted statistical methods described above.

38
 CHAPTER 4

RESULTS AND DISCUSSION

In this chapter, we first present no-feeding-tube simulation results to determine the min-
imum fluidization velocity, fluidization behavior, and the minimum bubbling velocity. Next,

we show the results of simulations run with the tube (see Section 3.5.3 for setup details).
The goal of these simulations is to relate the set of statistical quantities {ε, ψ, θ, Φ} with
the experimental product yields. The simulations provide visualization of the internal re-

actor flows, which helps illustrate the effect of different flow rate configurations on product
yields. Finally, based on what we observe from the simulations, we offer suggested flow rate
configurations for optimal mixing in a fluidized bed reactor.

4.1 Standard Configuration Fluidized Bed Simulations

Ten simulations were run in a standard, i.e. tubeless fluidized bed configuration with

a mesh resolution of ∼ 30 cells/cm at fluidizing velocities between 0.1−3.0 cm/s. These

simulations were run to calculate the minimum fluidization velocity, the minimum bubbling
velocity, and to analyze the internal flow behavior without the feed tube. The results are

compared to both experiment and literature models for {ε, ψ, Ω, χ}, the set of parameters
discussed in Section 2.3.

4.1.1 Determination of the Minimum Fluidization Velocity

Experimentally, the minimum fluidization velocity was determined by inserting a pressure

probe at the top and another at the bottom of the catalyst bed and recording the average

39
pressure drop across the bed given a constant fluidizing velocity without any gas flowing

from the feeding tube [41]. The fluidizing gas velocity was both increased and decreased
to find the time-averaged pressure drop at that velocity. Different pressure drop curves
were observed for the increasing and the decreasing velocity cases due to hysteresis of the

fluidization process (see Figure 4.1). The minimum fluidization velocity vmf was taken to
be the fluidizing velocity at which the slope of the pressure curve changed significantly. The
pressure drop found is measured between single points in space at the bottom and top of

the bed, and averaged in time.

c Fluent were initialized with a fluid velocity well above
The simulations run in ANSYS
vmf (vf l = 1.0 cm/s) and then adjusted to the higher or lower velocities. After the pressure

at the inlet readjusted to changes in velocity, we time-averaged the inlet pressure across
the entire inlet. Additionally, we calculated the standard deviation in time to investigate
the fluctuations in pressure about the mean. Figure 4.1 shows the experimental results for

increasing and decreasing the fluidizing velocity, the inlet pressure from Fluent, and vmf as
determined by the Abrahamsen & Geldart [1] and the Davies & Richardson [21] literature
models.

From Figure 4.1, we note that when increasing the fluidizing velocity in laboratory ex-
periments, the pressure drop behavior is marked by a sharp peak followed by a short trough,
which is indicative of hysteresis in the process of fluidization since the particles are ini-

tially packed closely and held together by friction and gravity. Thus, the decreasing velocity
experimental runs are generally used as the more accurate measurement. We observe a
significant change in the slope of the pressure vs. velocity in the decreasing-velocity case

between vf l = 0.3−0.4 cm/s, consistent with the experimentally determined minimum flu-
idization velocity vmf ≈ 0.363 cm/s [41]. The value of the minimum fluidization velocity for
the CFD results is also in good agreement with the literature value of vmf ≈ 0.3−0.4 cm/s.

In addition, we note that the average inlet pressure decreases while the standard deviation

40
 500

400

Pressure Drop (Pa)

300 Experiment (decreasing velocity)
Experiment (increasing velocity)
Abrahamsen & Geldart 1980
Davies & Richardson 1966
200 Fluent Simulations

100

0
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5
Fluidizing Velocity (cm/s)

Figure 4.1: Minimum fluidization velocity determination via pressure drop vs. fluidizing
velocity for Fluent simulations (bars indicate 2 standard deviations above and below the
mean), the experimental reactor, and literature models.

increases dramatically with higher velocity (vf l > 1.0 cm/s). There is thus a marked differ-
ence between simulations run below and above vf l = 1.0 cm/s. The standard deviation bars
demonstrate that the inlet pressure always fluctuates close to the experimentally determined

pressure (∼ 435 Pa), and the averaged inlet pressure at vf l = 3.0 cm/s is only 3.3% lower
than the experimental pressure at this velocity.
The Fluent simulations include a frictional model that accounts for frictional forces be-

tween particles at close packing (αs,min = 0.6). The maximum packing value of the catalyst
volume fraction is set to αs,max = 0.63. Without the use of this friction model, the pressure
drop curves were constant around 435 Pa for all Fluent simulations, including those below

41
the minimum fluidization velocity. This illustrates that the frictional model played a large

role in producing pressure curves that agree with the experimental pressure curves.

4.1.2 Determination of the Minimum Bubbling Velocity

To visualize fluidization and analyze the qualitative differences between flows with vf l >

1.0 cm/s and vf l < 1.0 cm/s, we look at contour plots of the catalyst volume fraction for
different velocities (see Figure 4.2). First, we note that at vf l < 0.4 cm/s, there are no
changes in bed height or the packing of the catalyst, since vf l < vmf . With vf l between 0.4

and 0.8 cm/s, fairly uniform bed expansion occurs with no bubbling. When vf l ≥ 0.9 cm/s,

Figure 4.2: Catalyst volume fraction profiles for (top, from left to right) vf l = 0.2, 0.3, 0.4,
0.6, and 0.8 cm/s and (bottom, from left to right) vf l = 0.9, 1.0, 1.5, 2.0, and 3.0 cm/s.

42
bubbles begin to rise throughout the bed. It is commonly known that for Geldart Type A

particles, such as this catalyst, vmb /vmf > 1 [29]. This indicates that we should expect bed
expansion before bubbling, as the simulations do indeed show.
In order to determine vmb from the simulations, we look at the parameter χ, the per-

centage of the bed occupied by the bubble phase. Additional parameters we would like to
consider are the scaled average catalyst volume fraction ε within the bed, the bed expansion
Ω, and the scaled bed variance of the catalyst ψ. These are plotted in Figure 4.3.

The bed expansion Ω and ε follow very similar trends when one is viewed on an inverted
scale from the other. This is expected, since, in general, h ∝ αs−1 (see Equation 3.1). There
is no significant bed expansion until vf l > 0.3 cm/s, followed by a roughly linear increase in

Ω until vf l = 1.0 cm/s, and then a less steep slope of Ω vs. vf l for vf l > 1.0 cm/s. These

0.0 0.10 0.25

ε − Catalyst mean χ − Bubble phase %
Ω − Bed expansion 1.0 ψ − Catalyst variance
vmb Abrahamsen & Geldart 1980
0.2 0.08 vmb Geldart 1973 0.20
0.8

0.4 0.06 0.15

0.6
Ω

ψ
χ
ε

0.6 0.4 0.04 0.10

0.2
0.8 0.02 0.05

0.0

1.0 0.00 0.00

0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5
vf l (cm/s) vf l (cm/s)

Figure 4.3: Four measures of the transition to the bubbling regime: (left) ε and Ω, (right)
h−h V ≤0.15
χ and ψ for tubeless Fluent simulations, where ε = αhα si
s,max
, Ω = hmfmf , χ = αVsbed , and
hα i
02
ψ = α02 s .
s,max

43
observations agree with those noted from the catalyst contours in Figure 4.2 and, in addition,

serve to quantify the bed expansion vs. fluidization velocity.

The dimensionless catalyst variance parameter ψ and the fraction of the bed occupied by
the bubble phase χ exhibit very similar trends. In particular, both parameters dramatically

increase with velocity for vf l > 1.0 cm/s. We find that vf l = 0.9 cm/s is the lowest velocity
at which bubbling occurs, indicating that vmb ≈ 0.9 cm/s. The Abrahamsen & Geldart
[1] model, which accounts for changes in gas properties, predicts vmb significantly lower

(vmb ≈ 0.6 cm/s) than either the CFD determined vmb or the simpler literature model
described by Geldart [29]. The simulation and the simpler literature model, however, are in
very close agreement with vmb ≈ 1.0 cm/s.

Although small bubbles do form at vf l = 0.9 and 1.0 cm/s, it is not until vf l = 1.5 cm/s
that a significant amount of the bed (2%) is occupied by the bubble phase in addition to a
significant increase in catalyst variance. The change in internal bed bubbling at and above

this fluidizing velocity provides a strong explanation for the observed changes in expansion
rates and in the variance of the catalyst. The catalyst variance parameter serves to reinforce
the fact that more bubbling means a higher variance in the catalyst.

Lastly, we compare the literature model for hmb /hmf , given by Equation (2.37), with the
CFD results. For this exercise, we declare vmb = 0.9 cm/s and vmf = 0.4 cm/s. This gives
the following relationship:

hmb 5.207 cm
= = 1.187
hmf 4.382 cm
 0.22  0.22 (4.1)
vmb 0.9 cm/s
≈ = = 1.195.
vmf 0.4 cm/s

We find very good agreement with the empirical correlation of Abrahamsen & Geldart. Our

simulations yield a hmb /hmf ratio that is only 0.7% lower than the empirical model predicts.

44
This relationship held true even though the value we found for vmb was about 50% higher
than the corresponding literature model predicted.

4.2 Effect of Gas Distribution

To test the effect of the gas distribution between the feed tube and fluidizer inlets, five
simulations were run with varying ratios of the feed tube gas flow to the total gas flow rate
(Qf t /Qtotal ≈ 8%, 17%, 33%, 51%, and 68%) and a constant total flow rate (Qtotal ≈ 55

cm3 /s) (see Tables 3.3 and 4.1).

4.2.1 Visualization

The catalyst volume fraction αs , the Ar feeding tube tracer gas mass fraction YAr , and

the axial velocity of the catalyst vs after the tracer gas is well-mixed are shown in Figure
4.4. From Figure 4.4a, we note that the Qf t /Qtotal ≈ 8%, 17%, and 33% cases (hereafter
f8, f17, and f33, respectively) all have similar αs profiles, bed heights, and bubble sizes.

The Qf t /Qtotal ≈ 51% and 68% cases (hereafter f51 and f68, respectively), however, show a
marked difference with a lower bed height, a denser bed, and spouts of particles along the
feeding tube wall. The feeding tube tracer gas profiles appear to show higher variance in

these latter two cases (see Figure 4.4b), particularly with a high Ar concentration along
 3
Qf t Qtotal
Case abbr. Qtotal (%) Qtotal cms cm cm



Qmf
vf t s
vf l s
f8 7.74 55.96 8.09 3.33 2.55
f17 17.48 55.76 8.06 7.49 2.27
f33 33.21 55.45 8.02 14.16 1.83
f51 51.11 55.10 7.97 21.65 1.33
f68 68.24 54.77 7.92 28.73 0.86

Table 4.1: Flow rate configurations and abbreviations for simulations with a near-constant
Qtotal as a function of Qf t .

45
 (a) αs

(b) YAr

(c) vs

Figure 4.4: Time snapshots of contour plots of (a) the catalyst volume fraction αs , (b) the
feed tube tracer gas mass fraction YAr , and (c) the catalyst axial velocity magnitude vs for
Fluent simulations with a feed tube with a constant total flow rate Qtotal ≈ 55 cm3 /s and
different ratios of Qf t /Qtotal (left to right: 8%, 17%, 33%, 51%, and 68%).

46
the feed tube. In Figure 4.4c, we note different flow behavior with the catalyst having a

downward flow along the feeding tube wall for f8, f17, and f33, and with the catalyst having
a high upward velocity along the feeding tube wall for f51 and f68.

4.2.2 Statistical Convergence

As discussed in Section 3.7.4, we must find TS for each case, the time period over which the
four mixing parameters {ε, ψ, θ, Φ} are in a statistically steady state. The four parameters
are plotted vs. sample length in Figure 4.5. The ε vs. TS plot shows that for all five cases,

full convergence occurs after about 1.5 s. The normalized variance of the catalyst volume
fraction converges after 3 s. The mean mixture fraction over the ideal mixture fraction θ
converges after about 2 s for f8, f17, f33, and f51 and after about 4 s for f68. Lastly, Φ shows

a similar set of convergence times with all five cases in full convergence for a sampling time
period of ∼ 4 s. For consistency, we use the same time period of TS = 4 s for all 5 of these
cases and for all the following statistical analyses.

4.2.3 Statistical Analysis of the Catalyst Volume Fraction

The normalized axial mean of αs , ε(r), shown in Figure 4.6, reveals the difference between
the low feed tube flow rate cases (f8, f17, and f33) and the high feed tube flow rate cases

(f51 and f68) with a significantly lower ε(r) along the feed tube wall (RT = 0.643 cm) in the
latter two cases due to more gas passing along the tube wall. We note a reversal of slope in

ε(r) vs. r between the two sets of cases, which results in a less dense bed for the lower vf t
cases and a more packed bed for the high vf t cases. When plotted against the bed height,
the packing of the bed is more evident for f51 and f68 with the latter packing to about 80%

the maximum possible packing vs. ∼ 70% for the former. Additionally, we note that the bed
height is defined by a sharp drop in ε(z) vs. z for f8, f17, and f33 while the bed height is less
sharply defined for the other two cases, with longer, gradual drops in ε(z) with increasing

47
 1.0 1.0

0.8 0.8

0.6 0.6
ε

θ
0.4 0.4
8% 8%
17% 17%
0.2 33% 0.2 33%
51% 51%
68% 68%
0.0 0.0
0 1 2 3 4 5 0 1 2 3 4 5
Sample Length TS (s) Sample Length TS (s)
0.5 1.0
8% 8%
17% 17%
0.4 33% 0.8 33%
51% 51%
68% 68%
0.3 0.6
Φ
ψ

0.2 0.4

0.1 0.2

0.0 0.0
0 1 2 3 4 5 0 1 2 3 4 5
Sample Length TS (s) Sample Length TS (s)

Figure 4.5: Statistical convergence of bed mixing parameters (ε, ψ, θ, and Φ) vs. sampling
period (TS ) for simulations run with QDf t /Q
Etotal ≈ 8%, 17%,
D E 33%, 51%, andD 68%E with a
02
constant Qtotal ≈ 55 cm3 /s, where ε = αs /αs,max , ψ = αs02 /αs,max , θ = ξ(r) /ξideal ,
D E D E2
and Φ = ξ 02 / ξ .

48
z. This gradualness of slope vs. z indicates a less even bed height for f51 and f68, which can

also be noted in the contour plots in Figure 4.4.

Despite the fact that there are more bubbles in the low feed tube flow rate cases, the

catalyst variance parameter ψ(r) peaks to a significantly higher value for the high feed tube
flow rate cases near the feed tube (see Figure 4.6). The catalyst variance ψ(z) also peaks
higher around the top of the catalyst bed for the high feed tube velocity cases. For the three

lower vf t cases, we note a nearly linear increase in variance of the catalyst with bed height
above the feeding tube (HT = 1.558 cm), which is indicative of bubble growth vs. bed height.
Such a linear growth is not seen for the high feed tube flow rate cases due to a general lack

of bubbling within the bed. Similar to the sharp drop in ε(z) for the low vf t cases, we note a
sharp drop in ψ(z) at the top of the catalyst bed and, in f51 and f68, wide trailing edges of
the peaks in variance. We again attribute these responses to the spouting of particles along

the feed tube wall above the catalyst bed, which cause variations in the bed height across
the radius of the bed.
The occurrence of the trough in ε(RT ) and the peak in ψ(RT ) for the high vf t cases can

be linked to the direction of flow along the feed tube wall, which reverses from negative
to positive when Qf t /Qtotal > 33%. The fluidizing velocity, and particularly how close
it is to vmb , is one of the main differences between cases run above and cases run below

Qf t /Qtotal = 33%. The three cases run at Qf t /Qtotal ≈ 8%, 17%, and 33% have respective
fluidization velocities of vf l = 2.55, 2.27, and 1.83 cm/s, which are well above vmb ≈ 0.9
cm/s. Meanwhile, the other two cases, where Qf t /Qtotal ≈ 51% and 69%, respectively have

vf l = 1.33 and 0.86 cm/s, much closer to vmb . Although vf l = 0.9 cm/s is the first velocity
associated with bubble formation, vf l = 1.5 cm/s was the lowest fluidizing velocity at which
both the bed catalyst variance ψ and the percentage of the bed occupied by the bubble phase

49
 1.0 1.0
8%
17%
0.8 0.8 33%
51%
68%
0.6 0.6

ε(z)
ε(r)

0.4 0.4
8%
17%
0.2 33% 0.2
51%
68%
0.0 0.0
0.0 0.5 1.0 1.5 2.0 2.5 0 2 4 6 8 10 12
r cm z cm

8% 8%
17% 17%
0.5 0.5
33% 33%
51% 51%
0.4 68% 0.4 68%
ψ(z)
ψ(r)

0.3 0.3

0.2 0.2

0.1 0.1

0.0 0.0
0.0 0.5 1.0 1.5 2.0 2.5 0 2 4 6 8 10 12
r cm z cm

Figure 4.6: The effect of the gas distribution between the fluidizer and the feed tube
(Qf t /Qtotal (8%, 17%, 33%, 51%, and 68%; Qtotal ≈ 55 cm3 /s) on the mixing of the cat-
alyst. Axial (left) and radial (right) time-ensemble means (top) and variances
D (bottom)
E of
the solids volume fraction from Fluent simulations with a feed tube: ε(r) = αs (r) /αs,max ,
D E D E D E
02 02
ε(z) = αs (z) /αs,max , ψ(r) = αs02 (r) /αs,max , and ψ(z) = αs02 (z) /αs,max .

50
 2.0 2.0
8% 8%
17% 17%
33% 33%
1.5 51% 1.5 51%
68% 68%

θ(z)
θ(r)

1.0 1.0

0.5 0.5

0.0 0.0
0.0 0.5 1.0 1.5 2.0 2.5 0 2 4 6 8 10 12
r cm z cm
0.5
8% 8%
17% 1.2 17%
0.4 33% 33%
51% 51%
1.0
68% 68%
0.3
0.8
Φ(z)
Φ(r)

0.6
0.2

0.4
0.1
0.2

0.0 0.0
0.0 0.5 1.0 1.5 2.0 2.5 0 2 4 6 8 10 12
r cm z cm

Figure 4.7: The effect of the gas distribution between the fluidizer and the feed tube
(Qf t /Qtotal (8%, 17%, 33%, 51%, and 68%; Qtotal ≈ 55 cm3 /s) on the mixing of the feed tube
tracer gas. Axial (left) and radial (right) time-ensemble means (top) and variances
D (bottom)
E
of the mixture fraction from Fluent simulations with a feed tube: θ(r) = ξ(r) /ξideal ,
D E D E D E2 D E D E2
θ(z) = ξ(z) /ξideal , Φ(r) = ξ 02 (r) / ξ(r) , and Φ(z) = ξ 02 (z) / ξ(z) .

51
χ were significantly high (see Figure 4.3). Thus, the two high vf t cases have a vf l less than
the significant bubbling velocity (1.5 cm/s).

4.2.4 Statistical Analysis of the Mixture Fraction

The mixing of the catalyst does not directly show how the feeding tube gas mixes within

the bed. Rather, we look at the mixture fraction statistics. Within the majority of the bed,
we note from Figure 4.7 that θ(r) and θ(z) are close to unity (nearly ideally mixed) for f8,
f17, and f33. The cases f51 and f68 do not show as good of mixing except near the feeding

tube. These high vf t cases show a linear increase in mixing (θ(z)) with bed height rather
than a well-mixed case throughout the whole bed.
The variance of the mixture fraction over the squared mean of the mixture fraction vs.

radius (Φ(r)) and vs. bed height (Φ(z)) show that the f33 case exhibits the lowest values of
Φ(r) and Φ(z). Thus, this case has the best mixing of the feed tube gas with the catalyst
bed. This is in strong agreement with experimental highest aromatic yield case (when 34% of

the total gas flow is injected through the feeding tube). The other four cases show increased
variance vs. r and z with f51 and f68 having the worst mixing of the feed tube tracer gas.

4.2.5 Bed Statistics

The mixing parameters (ε, ψ, θ, and Φ) described in Section 3.7.3 are plotted in Figure
4.8. First, we note that ε is nearly constant for the three low vf t cases (f8, f18, and f33)
and increases slightly for the other two cases (f51 and f68). The bed variance of the catalyst

additionally is similar for all five cases at ψ ≈ 0.18. Again, we analyze the mixture fraction
to determine the mixture quality and compare to the product yields. The bed mean of the

mixture fraction over the ideal mixture fraction is close to unity for the lower three vf t cases.
For the two higher vf t cases, θ is far from ideal and closer to 0.7 and 0.6.

52
 0.0

30

0.2
25

Carbon Yield %
ε − Catalyst mean
0.4 ψ − Catalyst variance 20
θ − Mixture fraction mean
Φ − Mixture fraction variance
Aromatics 15
0.6

10

0.8
5

1.0 0
0 10 20 30 40 50 60 70
Qf t /Qtotal (%)

Figure 4.8: Bed catalyst mean and variance, bed mixture fraction mean and variance (left
axis, solid lines), and aromatic product yield (left axis, dashed line) vs. the ratio of gas
distribution in the feed tube (Qf t /Qtotal ) with a constant Qtotal ≈ 55 cm3 /s, where ε =
D E D E D E D E D E2
02
αs /αs,max , ψ = αs02 /αs,max , θ = ξ /ξideal , and Φ = ξ 02 / ξ .

The mixture fraction variance over the mean mixture fraction squared Φ within the
catalyst bed is minimized at f33, which has the same ratio of Qf t /Qtotal that a maximum of
the aromatic yield was found experimentally (Qf t /Qtotal = 34%). Similarly, cases above and

below the Qf t /Qtotal = 34% point agree with the experimental data: the further the case
is from this optimum operating condition, the greater the drop in mixing quality (drop in
aromatic yield). The value of Φ is close to 2 times higher for the worst mixed case (Φ ≈ 0.3

for the 68% case) than for the ideally mixed case (Φ ≈ 0.15 for the 33% case).

53
4.3 Effect of Total Gas Flow Rate
The second set of simulations tests the effect of increasing the total flow rate (Qtotal /Qmf ≈
3, 4, 6, and 8; Qmf = 6.917 cm3 /s) while keeping a constant ratio of the feed tube gas flow
rate to the total inflow gas flow (Qf t /Qtotal = 33.21%). In this section, mixing parameter

definitions are excluded from figure captions as these can be found in prior figures, sections,
and the List of Symbols. A summary of simulations run can be found in Tables 3.3 and 4.2.

4.3.1 Visualization

As in Section 4.2.1, time snapshots of αs , YAr , and vs are shown after the tracer gas mass
fraction at the outlet is roughly constant (see Figure 4.9). From Figure 4.9a, we note that at

low total flow rate Qtotal /Qmf ≈ 3 and 4 (hereafter t3 and t4, respectively), there is a lack
of bubbling within the majority of the bed while bubbling is evident at the higher two total
flow rates, Qtotal /Qmf ≈ 6 and 8 (hereafter t6 and t8, respectively). For t3 and t4, vf l = 0.72

and 0.96 cm/s, respectively, and for t6 and t8, vf l = 1.44 and 1.92 cm/s, respectively. Thus,
the lack of bubbling within the bed is consistent with the lack of significant bubbling below
vf l = 1.5 cm/s. The tracer gas from the feeding tube shows a lack of mixing within the

majority of the bed for t3 and t4 as the gas rises from the feeding tube and passes quickly
through the system along the feeding tube outer wall. The t6 and t8 cases demonstrate that
the Ar tracer gas disperses radially before proceeding out of the bed.
 3
Qtotal Qf t
Case abbr. Qmf Qtotal cms cm cm



Qtotal
(%) vf t s
vf l s
t3 3.15 21.82 33.21 5.57 0.72
t4 4.21 29.09 33.21 7.43 0.96
t6 6.31 43.63 33.21 11.14 1.44
t8 8.41 58.18 33.21 14.85 1.92

Table 4.2: Flow rate configurations and abbreviations for simulations with a constant
Qf t /Qtotal as a function of Qtotal .

54
 (a) αs

(b) YAr

(c) vs

Figure 4.9: Time snapshots of contour plots of (a) the catalyst volume fraction αs , (b) the
feed tube tracer gas mass fraction YAr , and (c) the catalyst axial velocity magnitude vs for
Fluent simulations with a feed tube with a constant Qf t /Qtotal ≈ 33% and different ratios
of Qtotal /Qmf : (left-to-right: 3, 4, 6, and 8).

55
 The axial velocity component of the catalyst in Figure 4.9c reveals that for t3 and t4,

the catalyst rises quickly along the feeding tube. There is a marked shift in the higher total
flow rate cases with downflow of the catalyst along the feeding tube and a higher upward
flow in between the feeding tube and the reactor wall. This is very similar to the behavior

noted in the first set of simulations, and it can be noted that in both sets of simulations, the
change in flow occurs when vf l transitions from below ∼ 1.5 cm/s to above ∼ 1.5 cm/s, the
so-called significant bubbling velocity.

4.3.2 Statistical Convergence

As with the first set of simulations, we ensure that statistical convergence of the set of
mixing parameters {ε, ψ, θ, Φ} has been achieved before performing further statistics. The

sampling time period of statistical convergence TS is determined from Figure 4.10. All four
statistical quantities converge after about 4 s. The variance of the mixture fraction statistic
Φ fluctuates around the converged value for t3 and t4, but these fluctuations have relatively

small oscillations. Regardless of any small-scale variations, after 4 s, all the statistics per-
taining to the four different cases maintain the same qualitative behavior with varying flow
rates. We thus again establish the same TS = 4 s as was used for the first set of simulations.

4.3.3 Statistical Analysis of the Catalyst Volume Fraction

Plotted in Figure 4.11 are the catalyst volume fraction statistics as functions of radius
(ξ(r)) and of bed height (ξ(z)). The first of the plots reveals the difference between the two
lower total flow rate cases (t3 and t4) and the two higher total flow rate cases (t6 and t8).

Higher ε(r) is observed along the tube wall for t3 and t4 due to the gas bypassing along the
tube wall. The inverse slope of the catalyst volume fraction vs. radius for t6 and t8 shows a
marked difference in catalyst distribution. Ultimately, t6 and t8 have a higher void fraction

56
 1.0 1.0

0.8 0.8

0.6 0.6
ε

θ
0.4 0.4

3 3
0.2 4 0.2 4
6 6
8 8
0.0 0.0
0 1 2 3 4 5 0 1 2 3 4 5
Sample Length TS (s) Sample Length TS (s)
0.5 1.0
3 3
4 4
0.4 6 0.8 6
8 8

0.3 0.6
Φ
ψ

0.2 0.4

0.1 0.2

0.0 0.0
0 1 2 3 4 5 0 1 2 3 4 5
Sample Length TS (s) Sample Length TS (s)

Figure 4.10: Statistical convergence of bed mixing parameters (ε, ψ, θ, and Φ) vs. sampling
period (TS ) for simulations run with Qtotal /Qmf ≈ 3, 4, 6, and 8 with a constant Qf t /Qtotal ≈
33%.

57
within the majority of the bed, as is also to be expected with a higher fluidizing velocity.

The high void fraction with high Qtotal is also evident in the second plot, wherein the bed
height increases with Qtotal .
The catalyst variance parameter ψ(r) is highest along the feed tube wall due to large

bubbles passing along the tube for t3 and t4. The lack of internal bubbling causes these cases
to have a lower catalyst variance across the majority of the bed than for t6 and t8. The bed
height is characterized by a high catalyst variance, as demonstrated by ψ(z) in Figure 4.11.

From the width of the trailing edge of this peak, we can again ascertain whether or not the
bed height is flat in time and whether particles are blown off the top leading to uneven bed
heights. This trailing edge is elongated for t3 and t4 while it is a steep drop for t6 and t8.

The contour plots of αs in Figure 4.11 confirm the same uneven bed height behavior with
the particles blown high above the height of the bed along the feeding tube wall for t3 and
t4 and a more even dispersion spread across the radius of the bed in the other two cases.

4.3.4 Statistical Analysis of the Mixture Fraction

Figure 4.12 shows that for t6 and t8, θ(r) is close to unity across the radius of the bed
after the feed tube. For t3 and t4, however, the feeding tube tracer gas does not disperse

well radially into the bed with θ(r) ≈ 0.5−0.7. For t3 and t4, θ(z) increases with bed height,
only reaching ideal mixing above the catalyst bed (∼ 5−6 cm). The t6 and t8 case have
θ(z > HT ) ≈ 1. For z < HT , the feeding tube tracer gas is poorly mixed for all cases since it

must reach deep into the bed and against the flow of the fluidizing gas in order to mix with
this portion of the bed.
The mixture fraction variance parameter Φ(r) is high when r ≤ RT for t6 and t8. We

focus mostly on the variance across the radius in the majority of the bed (r > RT ), wherein
Φ(r) decreases steadily for these two cases (see Figure 4.12). This is indicative of better

mixing away from the feeding tube wall. We also note a higher Φ(r) for t3 and t4 with a

58
 1.0 1.0
3
4
0.8 0.8 6
8

0.6 0.6

ε(z)
ε(r)

0.4 0.4

3
0.2 4 0.2
6
8
0.0 0.0
0.0 0.5 1.0 1.5 2.0 2.5 0 2 4 6 8 10 12
r cm z cm

3 3
4 4
0.5 6 0.5 6
8 8
0.4 0.4
ψ(z)
ψ(r)

0.3 0.3

0.2 0.2

0.1 0.1

0.0 0.0
0.0 0.5 1.0 1.5 2.0 2.5 0 2 4 6 8 10 12
r cm z cm

Figure 4.11: The effect of increasing the total gas flow rate (Qtotal /Qmf ≈ 3, 4, 6, and 8;
Qf t /Qtotal ≈ 33%) on the mixing of the catalyst. Axial (left) and radial (right) time-ensemble
mean (top) and variance (bottom) of the solids volume fraction from Fluent simulations with
a feed tube.

59
 2.0 2.0
3 3
4 4
6 6
1.5 8 1.5 8

θ(z)
θ(r)

1.0 1.0

0.5 0.5

0.0 0.0
0.0 0.5 1.0 1.5 2.0 2.5 0 2 4 6 8 10 12
r cm z cm
0.5
3 3
4 1.2 4
0.4 6 6
8 8
1.0

0.3
0.8
Φ(z)
Φ(r)

0.6
0.2

0.4
0.1
0.2

0.0 0.0
0.0 0.5 1.0 1.5 2.0 2.5 0 2 4 6 8 10 12
r cm z cm

Figure 4.12: The effect of increasing the total gas flow rate (Qtotal /Qmf ≈ 3, 4, 6, and
8; Qf t /Qtotal ≈ 33%) on the mixing of the feed tube tracer gas. Axial (left) and radial
(right) time-ensemble mean (top) and variance (bottom) of the mixture fraction from Fluent
simulations with a feed tube.

60
peak of variance in between the feed tube wall and the outer wall. For all four cases, above

the feeding tube inlet, Φ(z) decreases with height, and it is significantly lower for t6 and t8
than for the other two cases. The high peaks of Φ(z) for t3 and t4 can be explained by the
fact that the Ar tracer gas is not mixing radially across the bed.

4.3.5 Bed Statistics

Figure 4.13 compares the statistical quantities ε, ψ, θ, and Φ to the aromatics yield
percentages. We observe that with increasing total flow rate, ε steadily decreases while ψ

increases. This finding simply serves to confirm that by increasing Qtotal , the bed becomes
less dense and the bubbling more vigorous.

0.0
40

35
0.2
30

Carbon Yield %
ε − Catalyst mean
0.4 ψ − Catalyst variance 25
θ − Mixture fraction mean
Φ − Mixture fraction variance
20
Aromatics
0.6
15

10
0.8

1.0 0
2 3 4 5 6 7 8 9
Qtotal /Qmf

Figure 4.13: Bed mean and variance (solid lines) with aromatic yield (dotted lines) as a
function of Qtotal /Qmf

61
 The mixture fraction is significantly closer to the ideally mixed case for t6 and t8 than

for t3 and t4: θ ≈ 0.9 for t6 and t8 and θ ≈ 0.55 for t3 and t4. By analyzing Φ, we note
that the feed tube tracer gas is increasingly better mixed as the flow rate increases, up to
Qtotal /Qmf ≈ 6, where a slight minimum of Φ is found. The high variance in the low total

flow rate cases can be traced back to the feeding tube gas rising along the feeding tube outer
wall, which causes a large volume of the gas to bypass the majority of the bed. The shallow
minimum of the mixture fraction variance over the mean-squared at t6 is not a complete

match with the significant peak in aromatic yield, noted at Qtotal /Qmf ≈ 4. Although the
CFD simulations do not exhibit the same exact point of best mixing as the experimental
reactor, the simulations do still show a subtle, best-mixed case vs. Qtotal at vf l ≈ 1.5 cm/s,

the point of significant internal bubbling. This shows that the bubbling point is a critical
point for mixing within the fluidized bed.

62
 CHAPTER 5

CONCLUSIONS

c Fluent were used

Computational fluid dynamics (CFD) simulations run on ANSYS
to simulate a three-phase fluidized bed reactor (FBR). As in a traditional FBR, the bed
is fluidized by a gas stream entering through the bottom of the reactor. An additional gas

stream enters downward through a vertical tube in the center of the reactor. The overarching
objective is to optimize the mixing of this second gas stream with the particles in the bed.
The simulations are validated using experimental data and theory for the traditional two-

phase configuration. Good agreement was found with the minimum fluidization velocity and
it was noted that in the simulations, bubbling within the bed did not occur until a fluidizing
velocity of 0.9 cm/s, the minimum bubbling velocity. Significant bubbling did not occur until

vf l = 1.5 cm/s. Ultimately, the significance of these findings manifested in the simulations
run with a feeding tube as it was noted that heavy feeding tube gas bypassing occurred along
the feeding tube wall when the fluidizing velocity was below 1.5 cm/s. The feed tube tracer

gas dispersed radially into the majority of the bed for cases run with vf l & 1.5 cm/s. The
flow of the gas along the feed tube wall is thus important in determining whether a flow rate
configuration will be well-mixed or not. The feed tube disperses radially into the bed only

if the bed consists of a more dilute, bubbling bed. This could be the primary explanation
for why there is a peak in aromatic yield when adjusting both the fluidizing and the feeding
tube velocities simultaneously.

In the first set of simulations in which the feeding tube flow rate is varied with a constant
total flow rate, we find very good agreement with the experimental reactor. Here, we find

63
 D E D E2
a minimum of Φ = ξ 02 / ξ , and thus a maximum of mixing of the feeding tube tracer
gas with the catalyst bed, at the same Qf t /Qtotal ≈ 34% that produced the highest aromatic

yield. In the second set of simulations, where Qf t /Qtotal was held constant and Qtotal was
varied, there is also a minimum in Φ at Qtotal /Qmf ≈ 6. This minimum is close to but not
at the same Qtotal /Qmf = 4 that a maximum aromatic yield was found experimentally.

It should be noted that the CFD simulations were run with multiple assumptions, in-
cluding a 2D-axisymmetric assumption, no cellulose or chemical kinetics, and a constant
particle diameter (instead of the actual bimodal distribution). Thus, it is not expected that

the CFD simulations will fully match the behavior of the experimental fluidized bed. As the
fluidized bed with the tube is enclosed and insulated, it is difficult to confirm the fluidization
behavior within the bed such as the flow direction along the walls or if internal bubbling

occurs. One experimental test could provide validity to the CFD simulations if vmb is noted
to be around 1.0−1.5 cm/s. If this is indeed the case, then, as the simulations show, it is
likely that the radial dispersion of the feeding tube gas (and thus the radial dispersion of

the gasified cellulose) is best when the internal bed is bubbling.

From the simulations, it was shown that the bed is the most well-mixed in configurations
right above the bubbling border. Thus, we conclude, that the highest aromatic yields should

occur when both vf l ≥ vmb and there is a significant amount of feeding tube flow rate
(Qf t /Qtotal ≈ 34 % in this case). A further test of this theory would be to determine vmb
experimentally, and then to fix vf l to vmb and adjust the feeding tube flow rate, noting

the aromatic yields. Based on what we observe in our simulations, vmb can be determined
experimentally in the enclosed reactor by analyzing the standard deviation in time of the
inlet pressure vs. fluidizing velocity, as was done in Figure 4.1, and noting when the standard

deviation increases significantly. Running cases with these flow rate configurations (vf l ≈

64
vmb ) and adjusting vf t is our suggestion for both further experimental and further simulation
work.

Additional future work to extend from this thesis is running simulations to test how the
mixing is affected by the scale-up of the reactor to a full industrial size, running a 3D model
simulation in order to verify the 2D-axisymmetric model by comparing the volume-weighted

statistics for each, and simulating chemical kinetics. Scale-up, 3D simulations, and chemical
reactions provide computational challenges, however, particularly because simulation run

times were not significantly improved when run on 32 processors instead of 8 processors.
Further tests of the effect of the mesh size and time step size could help alleviate the problem
of long run times by using a more coarse mesh or larger time steps.

The computational fluid dynamics simulations of this non-standard fluidized bed provide
a challenge in both numerical modeling and computational effort. The simulations operate
as both confirmation for the experimental reactor and a stand-alone testing of hydrodynamic

mixing within an unconventional biomass reactor feeding mechanism through the use of a
central, downward-facing feeding tube. The statistical analyses and methods described and
implemented in this thesis can further serve to quantify mixing within fluidized beds of any

configuration.

65
 APPENDIX

GOVERNING EQUATIONS AND CLOSURE MODELS

Continuity
∂αg
+ ∇ · (αg vg ) = 0
∂t [38]
∂αs
+ ∇ · (αs vs ) = 0
∂t
Momentum
 
∂vg
ρg αg + vg · ∇vg = −αg ∇P + ∇ · αg τ̄¯g + αg ρg g − β(vg − vs )
∂t
  [38]
∂vs
ρs αs + vs · ∇vs = −αs ∇P + ∇ · τ̄¯s − ∇Ps + αs ρs s + β(vg − vs )
∂t
Granular Temperature Θ
1 D 02 E
Θ= v [52]
3 s
Granular Energy
 
3 ∂
(αs ρs Θ) + ∇ · (αs ρs Θvs ) = − Ps I¯ + τ̄¯s : ∇vs


2 ∂t [52]
+ ∇ · (κs ∇Θ) − γs − Js
Granular Temperature (Algebraic) Θ
¯ )
−(K1 αs + ρs )tr(D̄s
Θ=
2αs K4
q [71]
(K1 αs + ρs )2 tr2 (D̄ ¯ 2 ) + K tr2 (D̄
¯ ) + 4K α [2K tr(D̄ ¯ )]
s 4 s 3 s 2 s
+
2αs K4

Table A.1: Governing equations and closure models: I

66
Granular Temperature Algebraic Abbreviations K1 , K2 , K3 , & K4

K1 = 2(1 + e)ρs g0
4 2
K2 = √ ds ρs (1 + e)αs g0 − K3
3 π 3
 √
[71]
  
ds ρs π 2 8αs
K3 = 1 + (1 + e)(3e − 1)αs g0 + √ g0 (1 + e)
2 3(3 − e) 5 5 π
2
12(1 − e )ρs g0
K4 = √
ds π

Collision Disspipation γs

αs2 ρs g0 3
γs = 12(1 − e2 ) √ Θ2 [52]
ds π
Fluctuation Dissipation Js

βds (vg − vs )2
 
Js = β 3Θ − √ [52]
4αs ρs πΘ
Solids Pressure Ps
Ps = αs ρs Θ + 2g0 αs2 ρs Θ(1 + e) [52]
Solids Bulk Viscosity λs
r
4 Θ [52]
λs = αs ρs ds g0 (1 + e)
3 π
Viscous stress tensor of phase k, τ̄¯k
 
¯ 2 ¯ )I¯
τ̄¯k = 2µk D̄k + λk − µk tr(D̄ k
3 [75]
¯ = 1 [∇v + (∇v )> ]
D̄ k k k
2

Table A.2: Governing equations and closure models: II

67
Solids Thermal Conductivity κs
 2 r
2 6 Θ
κs = 1 + (1 + e)g0 αs κdil + 2αs2 ρs ds g0 (1 + e)
(1 + e)g0 5 π [30]
75 √
κdil = ρs ds πΘ
384
Solids Shear Viscosity µs
√  2
10ρs ds Θπ 4
µs = 1 + g0 αs (1 + e) αs [30]
96αs (1 + e)g0 5

Frictional Viscosity µf

(αs − αs,min )n
µf = F r sin φ [40]
(αs,max − αs )p

Radial Distribution Model g0

 −2.5αs,max
αs
g0 = 1− [51]
αs,max

Drag Model β

αs2 µg 7 ρg αs |vg − vs |
β = 150 + (αs ≥ 0.2);
αg d2s 4 ds
3 ρg αs αg |vg − vs | −2.65
β = Cd αg (αs < 0.2)
4 ds [30]
 
24 0.687
CD = 1 + 0.15(Res ) (Res < 1000);
Res
CD = 0.44 (Res ≥ 1000)

Table A.3: Governing equations and closure models: III

68
 BIBLIOGRAPHY

[1] Abrahamsen, R., and Geldart, D. Behavior of gas-fluidized beds of fine powders. 1.
homogeneous expansion. Powder Technology 26 (1980), 35–46.

[2] Agrawal, K., Loezos, P. N., Syamlal, M., and Sundaresan, S. The role of meso-scale
structures in rapid gas-solid flows. Journal of Fluid Mechanics 445 (2001), 151–185.

[3] ANSYS, Inc. ANSYS FLUENT 12.0 Theory Guide, 2009.

[4] ANSYS, Inc. ANSYS FLUENT 12.0 User’s Guide, 2009.

[5] Asegehegn, T. W., Schreiber, M., and Krautz, H. J. Numerical simulation and ex-
perimental validation of bubble behavior in 2D gas-solid fluidized beds with immersed
horizontal tubes. Chemical Engineering Science 66, 21 (2011).

[6] Baeyens, J. Heat transfer in gas fluidized beds. PhD thesis, University of Bradford,
1973.

[7] Bauer, W., and Werther, J. The role of gas-distribution in fluidized-bed chemical reactor
design. Chemical Engineering Communications (1982).

[8] Bilger, R. W. Turbulent flows with nonpremixed reactants. In Turbulent reacting flows.
Springer, 1980, pp. 65–113.

[9] Boemer, A., Qi, H., Renz, U., Vasquez, S., and Boysan, F. Eulerian computation of
fluidized bed hydrodynamics—a comparison of physical models. In Proceedings of the
20th International Conference on Fluidized Bed Combustion (1995).

[10] Bokkers, G. A., van Sint Annaland, M., and Kuipers, J. A. M. Mixing and segregation
in a bidisperse gas-solid fluidised bed: a numerical and experimental study. Powder
Technology 140 (2004), 176–186.

[11] Bonilla, C. F., Brooks, R. D., and Walker, Jr., P. L. The viscosity of steam and nitrogen
at atmospheric pressure and high temperatures. In Proceedings of General Discussion
on Heat Transfer (1951).

[12] Broadhurst, T. E., and Becker, H. A. Onset of fluidization and slugging in beds of
uniform particles. American Institute of Chemical Engineers Journal 21 (1975), 238–
247.

69
[13] Chapman, S., and Cowling, T. S. The mathematical theory of non-uniform gases. Cam-
bridge University Press, 1970.

[14] Chu, K. W., and Yu, A. B. Numerical simulation of complex particle-fluid flows. Powder
Technology 179 (2008), 104–114.

[15] Chung, T. J. Computational Fluid Dynamics. Cambridge University Press, 2002.

[16] Coronella, C. J., Yee, S. Y., and Seader, J. D. Minimum slugging velocity in fluidized
beds containing vertical rods. Fuel 73, 9 (1993), 1537–1543.

[17] Coulomb, C. A. Essai sur une application des régles de maximis et minimis á quelques
problémes de statique, relatifs á l’architecture. Mémoires de Mathématique & de
Physique présentés á l’Académie Royale des Sciences par divers Savans, & lûs dans
ses Assemblées 7 (1776), 343–382.

[18] Crowe, C., Sommerfeld, M., and Tsuji, Y. Multiphase Flows with Droplets and Particles.
CRC Press, 1998.

[19] Das Sharma, S., Pugsley, T., and Delatour, R. Three-dimensional CFD model of the
deaeration rate of FCC particles. American Institute of Chemical Engineers Journal 52
(2006), 2391–2400.

[20] Davidson, J. F., Clift, R., and Harrison, D., Eds. Fluidization, 2nd ed. Academic Press,
1985.

[21] Davies, L., and Richardson, J. F. Gas interchange between bubbles and the continuous
phase in a fluidised bed. Transactions of the Institution of Chemical Engineers 44
(1966), T293–T305.

[22] Deen, N. G., van Sint Annaland, M., van der Hoef, M. A., and Kuipers, J. A. M. Review
of discrete particle modeling of fluidized beds. Chemical Engineering Science 62 (2007),
28–44.

[23] Di Renzo, A., and Di Maio, F. P. Homogeneous and bubbling fluidization regimes
in DEM-CFD simulations: hydrodynamic stability of gas and liquid fluidized beds.
Chemical Engineering Science 62 (2007), 116–130.

[24] Enwald, H., Peirano, E., and Almstedt, A. E. Eulerian two-phase flow theory applied
to fluidization. International Journal of Multiphase Flow 22 (1996), 21.

[25] Ergun, S. Fluid flow through packed columns. Chemical Engineering Progress 48, 2
(1952), 89–94.

[26] Ferschneider, G., and Mege, P. Eulerian simulation of dense phase fluidized beds. Rev.
Inst. Fr. Pet. 51, 2 (1996), 301–307.

70
[27] Frantz, J. S. Minimum fluidization velocities and pressure drop in fluidized beds. Chem-
ical Engineering Progress Symposium Series 62 (1966), 21–31.

[28] Gao, J., Chang, J., Xu, C., Lan, X., and Yang, Y. CFD simulations of gas solid flow in
FCC strippers. Chemical Engineering Science 63 (2008), 1827–1841.

[29] Geldart, D. Types of gas fluidization. Powder Technology 7 (1973), 285–292.

[30] Gidaspow, D. Multiphase flow and fluidization. Academic Press, San Diego, 1994.

[31] Hong, R. Y., Li, H. Z., Li, H. B., and Wang, Y. Studies on the inclined jet penetration
length in a gas-solid fluidized bed. Powder Technology 92, 3 (1997), 205–212.

[32] Hoomans, B. P. B., Kuipers, J. A. M., Briels, W. J., and Van Swaaij, W. P. M. Discrete
particle simulation of bubble and slug formation in a two-dimensional gas-fluidized bed:
a hard sphere approach. Chemical Engineering Science 51, 1 (1996), 99–118.

[33] Hosseini, S. H., Rahimi, R., Zivdar, M., , and Samimi, A. CFD simulation of gas-solid
bubbling fluidized bed containing FCC particles. Korean Journal Chemical Engineering
26, 5 (2009), 1405–1413.

[34] Hosseini, S. H., Zhong, W., Esfahany, M. N., Pourjafar, L., and Azizi, S. CFD simulation
of the bubbling and slugging gas-solid fluidized beds. Journal of Fluids Engineering 132,
4 (2010).

[35] Hull, A. S., Chen, Z., and Agarwal, P. K. Influence of horizontal tube banks on the
behavior of bubbling fluidized beds: 1. bubble hydrodynamics. Powder Technology 103
(1999), 230–242.

[36] Hull, A. S., Chen, Z., and Agarwal, P. K. Influence of horizontal tube banks on the
behavior of bubbling fluidized beds: 2. mixing of solids. Powder Technology 111 (2000),
192–199.

[37] Ishii, M. Thermo-Fluid Dynamic Theory of Two-Phase Flow. Direction des Etudes et
Recherches d’Electricitò de France. Eyrolles, 1975.

[38] Jackson, R. Locally averaged equations of motion for a mixture of identical spherical
particles and a newtonian fluid. Chemical Engineering Science 52 (1997), 2457.

[39] Jenkins, J. T., and Savage, S. B. A theory for the rapid flow of identical smooth, nearly
elastic spherical particles. Journal of Fluid Mechanics 130 (1983), 187–202.

[40] Johnson, P. C., and Jackson, R. Frictional collisional constitutive relations for
antigranulocytes-materials, with application to plane shearing. Journal of Fluid Me-
chanics 176 (1987), 67–93.

71
[41] Karanjkar, P., Coolman, R. J., Blatnik, M. T., Almalkie, S., de Bruyn Kops, S. M.,
Mountziaris, T. J., Conner, W. C., and Huber, G. W. Hydrodynamics of a fluidized bed
reactor for catalytic fast pyrolysis of cellulose to make green aromatics. In preparation,
2013.

[42] Kestin, J., Ro, S. T., and Wakeham, W. A. Viscosity of noble gases in the temperature
range 25-700o C. The Journal of Chemical Physics 56, 8 (1972), 4119–4124.

[43] Law, C. L., Tasirin, S. M., Daud, W. R. W., and Geldart, D. Effect of vertical baffles
on particle mixing and drying in fluidized beds of group d particles. China Particuology
1, 3 (2003), 115–118.

[44] Li, T., Dietiker, J., Zhang, Y., and Shahnam, M. Cartesian grid simulations of bubbling
fluidized beds with a horizontal tube bundle. Chemical Engineering Science 66 (2011),
6220–6231.

[45] Li, T., Pougatch, K., Sacludean, M., and Grecov, D. Mixing of a secondary gas injection
in a bubbling fluidized bed. Chemical Engineering Research and Design 87 (2009),
1451–1465.

[46] Li, T., Pougatch, K., Salcudean, M., and Grecov, D. Numerical simulation of horizontal
jet penetration in a three-dimensional fluidized bed. Powder Technology 184, 1 (2008),
89–99.

[47] Lindborg, H., and Jakobsen, H. A. Sorption enhanced steam methane reforming process
performance and bubbling fluidized bed reactor design by use of a two-fluid model.
Industrial and Engineering Chemistry Research 48 (2009), 1332–1342.

[48] Lindborg, H., Lysberg, M., and Jakobsen, H. A. Practical validation of the two-fluid
model applied to dense gas-solid flows in fluidized beds. Chemical Engineering Science
62 (2007), 5854–5869.

[49] Loezos, P. N., Costamagna, P., and Sundaresan, S. The role of contact stresses and wall
friction on fluidization. Chemical Engineering Science 57 (2002), 5123–5141.

[50] Louge, M. Y., Mastorakos, E., and Jenkins, J. T. The role of particle collisions in
pneumatic transport. Journal of Fluid Mechanics 231, 8 (1991), 345.

[51] Lun, C. K. K., and Savage, S. B. The effects of an impact velocity dependent coefficient
of restitution on stresses developed by sheared granular matrials. Acta Mechanica 63
(1986), 15–44.

[52] Lun, C. K. K., Savage, S. B., Jeffrey, D. J., and Chepurniy, N. Kinetic theories for
granular flow - inelastic particles in couette-flow and slightly inelastic particles in a
general flowfield. Journal of Fluid Mechanics 140 (1984), 223–256.

72
[53] Mckeen, T., and Pugsley, T. Simulation and experimental validation of a freely bubbling
bed of FCC catalyst. Powder Technology 129 (2003), 139–152.

[54] Merry, J. M. D. Penetration of a horizontal gas jet into a fluidised bed. Transactions
of the Institution of Chemical Engineers and the Chemical Engineer 49, 4 (1971), 189.

[55] Miller, C. O., and Logwinuk, A. K. Fluidization studies of solid particles. Engineering
and Process Development 43, 5 (1961), 1220–1226.

[56] Olsson, S. E., Wiman, J., and Almstedt, A. E. Hydrodynamics of a pressurized fluidized
bed with horizontal tubes: Influence of pressure, fluidization velocity and tube-bank
geometry. Chemical Engineering Science 50, 4 (1995), 581–592.

[57] Ozawa, M., Umekawa, H., Furui, S., Hayashi, K., and Takenaka, N. Bubble behavior
and void fraction fluctuation in vertical tube banks immersed in a gas-solid fluidized-bed
model. Experimental Thermal and Fluid Science 26 (2002), 643–652.

[58] Papadikis, K., Gu, S., and Bridgwater, A. V. Geometrical optimzation of a fast pyrolosis
bubbling fluidized bed reactor using computational fluid dynamics. Energy Fuels 24
(2010), 5634–5651.

[59] Parmentier, J. F., Simonin, O., and Delsart, O. A numerical study of fluidization
behavior of Geldart B, A/B and A particles using eulerian multifluid modeling approach.
In Ninth International Conference on Circulating Fluidized Beds (2008), J. Werther,
K. E. Wirth, and W. Nowak, Eds.

[60] Patil, D. J., Annaland, A. V., and Kuipers, J. A. M. Critical comparison of hydrody-
namic models for gas-solid fluidized beds - Part II: freely bubbling gas-solid fluidized
beds. Chemical Engineering Science 60, 1 (2005), 73–84.

[61] Petersen, H. The properties of helium: density, specific heats, viscosity, and thermal
conductivity at pressures from 1 to 100 bar and from room temperature to about 1800
K. Tech. rep., Danish Atomic Energy Commission: Research Establishment Risö, 1970.

[62] Pillai, B. C., and Rao, R. M. Pressure drop & minimum fluidization velocities in air-
fluidized beds. Indian Journal of Technology 9 (1971), 77–86.

[63] Ramamoorthy, S., and Subramanian, N. Axial solids mixing and bubble characteristics
in gas-fluidized beds with vertical internals. The Chemical Engineering Journal 22
(1981), 237–242.

[64] Robinson, C. S. Some factors influencing sedimentation. Industrial and Engineering

Chemistry 18 (1926), 869–871.

[65] Rowe, P. N. Drag forces in a hydraulic model of a fluidized bed: II. Transactions of the
Institution of Chemical Engineers 39 (1961), 175.

73
[66] Rüdisüli, M., Schildhauer, T. J., Biollaz, S. M. A., and van Ommen, J. R. Bubble char-
acterization in a fluidized bed with vertical tubes. Industrial & Engineering Chemistry
Research 51 (2012), 4748–4758.

[67] Schaeffer, G. Instability in the evolution-equations describing incrompressible antigran-

ulocytes flow. Journal of Differential Equations 66, 1 (1987), 19–50.

[68] Song, X. Q., Grace, J. R., Bi, H., Lim, C. J., Chan, E., Knapper, B., and McKnight,
C. A. Gas mixing in the reactor section of fluid cokers. Industrial & Engineering
Chemistry Research 44, 16 (2005), 6067–6074.

[69] Srivastava, A., and Sundaresan, S. Role of wall friction in fluidization and standpipe
flow. Powder Technology 124 (2002), 45–54.

[70] Stokes, G. G. Mathematical and physical papers. Cambridge: University Press, 1880.

[71] Syamlal, M., Rogers, W., and J., O’Brien T. Mfix Documentation Theory Guide. U.S.
Department of Energy, Office of Fossil Energy Technical Note, 1993.

[72] Taghipour, F., Ellis, N., and Wong, C. Experimental and computational study of gas-
solid fluidized bed hydrodynamics. Chemical Engineering Science 60 (2005), 6857–6567.

[73] Tsinontides, S. C., and Jackson, R. The mechanics of gas-fluidized beds with an interval
of stable fluidization. Journal of Fluid Mechanics 255 (1993), 237–274.

[74] van Wachem, B., and Sasic, S. Derivation, simulation and validation of a cohesive
particle flow CFD model. American Institute of Chemical Engineers Journal 54 (2008),
9–19.

[75] van Wachem, B. G. M., Schouten, J. C., van den Bleek, C. M., Krishna, R., and
Sinclair, J. L. Comparative analysis of CFD models of dense gas-solid systems. American
Institute of Chemical Engineers Journal 47, 5 (2001), 1035–1051.

[76] Wang, J., van der Hoef, M. A., and Kuipers, J. A. M. Why the two-fluid model fails to
predict the bed expansion characteristics of Geldart A particles in gas-fluidized beds:
A tentative answer. Chemical Engineering Science 64, 3 (2009), 622–625.

[77] Wang, J., van der Hoef, M. A., and Kuipers, J. A. M. CFD study of the minimum
bubbling velocity of Geldart A particles in gas-fluidized beds. Chemical Engineering
Science 65, 12 (2010), 3772–3785.

[78] Wang, J., Van der Hoef, M. A., and Kuipers, J. A. M. The role of scale resolution versus
inter-particle cohesive forces in two-fluid modeling of bubbling fluidization of Geldart
A particles. Chemical Engineering Science 66 (2011), 4229–4240.

[79] Wen, C. Y., and Yu, Y. H. Mechanics of fluidization. Chemical Engineering Progress
Symposium Series 62 (1966), 100.

74
[80] Wu, Y. Effect of operating temperature on minimum fluidization velocity. Powder
Technology 67, 2 (1991), 217–220.

[81] Yates, J. G., Ruiz-Martinez, R. S., and Cheesman, D. J. Prediction of bubble size in a
fluidized bed containing horizontal tubes. Chemical Engineering Science 45, 4 (1990),
1105–1111.

[82] Zhu, H. P., Y., Zhou Z., Yang, R. Y., and Yu, A. B. Discrete particle simulation of
particulate systems: A review of major applications and findings. Chemical Engineering
Science 63 (2008), 5728–5770.

[83] Zhu, H. P., and Zhou Z. [Link] Yu, A. B. Discrete particle simulation of particulate
systems: Theoretical developments. Chemical Engineering Science 62 (2007), 3378–
3396.

[84] Zimmerman, S., and Taghipour, F. CFD modeling of the hydrodynamics and reaction
kinetics of FCC fluidized-bed reactors. Industrial & Engineering Chemistry Research
44 (2005), 9818–9827.

75