Research Article pubs.acs.org/journal/ascecg
Predictive Model for Particle Residence Time Distributions in Riser Reactors. Part 1: Model Development and Validation Thomas D. Foust,*,†,§ Jack L. Ziegler,†,§ Sreekanth Pannala,‡ Peter Ciesielski,† Mark R. Nimlos,† and David J. Robichaud† †
National Bioenergy Center, National Renewable Energy Laboratory, 15013 Denver West Parkway, Golden, Colorado 80401, United States ‡ Corporate Research and Development, SABIC Americas, 14100 Southwest Freeway, Sugar Land, Texas 77478, United States S Supporting Information *
ABSTRACT: In this computational study, we model the mixing of biomass pyrolysis vapor with a solid catalyst in circulating riser reactors with a focus on the determination of solid catalyst residence time distributions (RTDs). A comprehensive set of 2D and 3D simulations were conducted for a pilot-scale riser using the Eulerian−Eulerian two-fluid modeling framework with and without subgrid-scale models for the gas−solid interaction. A validation test case was also simulated and compared to experiments, showing agreement in the pressure gradient and RTD mean and spread. For simulation cases, it was found that for accurate RTD prediction, the Johnson and Jackson partial slip solid boundary condition was required for all models, and a subgrid model is useful so that ultra high resolutions grids, which are very computationally intensive, are not required. We discovered a 2/3 scaling relation for the RTD mean and spread when comparing resolved 2D simulations to validated unresolved 3D subgridscale model simulations. KEYWORDS: Biomass pyrolysis, Catalytic upgrading, Riser reactor, Multiphase flow simulation, Catalyst residence time distribution
■
INTRODUCTION Many important chemical conversion processes are carried out in fluidized catalytic cracking (FCC) reactors, which allow for short and controlled contact times between vapors and particulate catalyst. FCC also provides a convenient way to enable continuous regeneration of the catalyst, which can become deactivated (e.g., due to coking) over time. One specific example of an FCC process of intense current interest is catalytic upgrading1,2 of biomass pyrolysis vapors, or catalytic fast pyrolysis (CFP) as it is commonly referred to, is a promising strategy for improving the compatibility of bio-oil with the conventional petroleum processing infrastructure. Upgrading pyrolysis vapors with acidic zeolites such as HZSM5 can convert the complex set of oxygenated compounds found in pyrolysis vapors into a much simpler mixture of compounds.3 The undesirable properties of the raw bio-oil, including instability, viscosity, corrosivity, and immiscibility, are greatly reduced by the CFP process. FCC reactors allow for short, and controlled, contact time with the catalyst while enabling continuous regeneration, which will be critical for performing CFP at the commercial scale4,5 because HZSM-5 and other catalysts that appear effective for CFP tend to deactivate quickly, presumably through coking.2,6−9 Development of CFP to the commercial scale will require methods for accurately predicting and controlling residence time distribution of catalyst particles such as HZSM-5 © 2017 American Chemical Society
since this is a key factor affecting performance because yields, product suites, and coking and charring rates are direct functions of catalyst contact times. Using bulk flow calculations as described in the Supporting Information is not accurate for predicting particle RTDs due to clustering, wall effects, and other flow factors that have significant effects on particle RTDs. Particle/catalyst residence time distributions (RTDs) in FCC risers have been studied experimentally for more than 30 years.10−17 Some of the first papers to discuss RTDs from both experimental and modeling perspectives are Ambler et al.12 and Rhodes et al.17 A more recent paper that studies a nonreactive air case using simulation and experiments with crystallized NaCl tracers is Andreux et al.18 Harris et al.13 also determined RTDs experimentally using phosphorescent particles with photomultiplier detection. Several groups have worked on methods to computationally predict particle RTDs. Hua et al.19 used Fluent incorporating the energy minimization multiscale (EMMS) drag model with a coarse 3D grid to see if they could match the Andreux experimental measurements. The RTDs were calculated using the time elapsed diffusion of tracers. Using a solids flowReceived: November 1, 2016 Revised: January 11, 2017 Published: January 16, 2017 2847
DOI: 10.1021/acssuschemeng.6b02384 ACS Sustainable Chem. Eng. 2017, 5, 2847−2856
Research Article
ACS Sustainable Chemistry & Engineering dependent diffusion parameter for the tracer yielded good agreement with the experimental RTD. Liu et al.16 conducted a modeling study of varying the riser height to control RTDs. Using the commercial Barracuda software,20 they found the RTD to be proportional to height and the spread to increase with increasing height. Other notable references11,14 also simulated RTDs in risers for the reactive and nonreactive petroleum cases at the industrial scale. Affum et al.11 modeled the reactive case, and Lan et al.14 used the multiphase particle in cell (MP-PIC) method to compare nonreactive RTDs for the fast-fluidization and dilute-phase transport regimes. Most of these modeling efforts were able to qualitatively confirm the experimental results indicating an RTD distribution with a single peak and a long tail or in the case of Ambler et al.12 a possible smaller secondary peak in the tail. However, a significant limitation of these works was that coarse grids were utilized, and hence, the flow was not resolved. Because clustering, back mixing, random solids fluctuations, boundary conditions, and other aspects of the flow field could have a dramatic impact on RTDs, it is highly desirable to resolve the flow field for accurately determining RTDs, and hence, this will be the focus of this work. Given the magnitude and complexity of this undertaking, we divided this effort in a two-part staged project. In the first stage (this paper), we focus on the simulations formulation and methodology. We show results to validate our modeling approach compared to experimental RTD measurements from Andreux et al.18 In the second stage, part two of this study, we assess the functional utility of using this validated model to assist in the development of CFP processes in FCC reactors to a commercially viable state.
ing scale models. If the fluid−particle interactions are approximated with drag models, while still resolving the motion and collisions of single particles, more particles can be included, and this is the approach used by the discrete element method (DEM)25,26 and the NGA-multiphase flow software.27 The next level of approximation is the particle in cell and cluster modeling methods such as MP-PIC, which is implemented in the commercial software Barracuda20 and MFIX.23 However, this approach is still too computationally intensive for pilot-scale reactors so the TFM approach was chosen for this study given its compromise to accurately model RTDs while providing the computational efficiency to allow an engineering-scale model of the CFP process. The PDEs for the TFM model are shown below. There are continuity, momentum, and mass conservation equations for the gas mixture and the single solid phase. All transport equations in MFIX are solved in the nonconservative form. A reaction term is included in this study for the gas species but not the solid phase because the CFP kinetic formulation that will be detailed in part two of this study does not require reaction terms for the catalyst phase. In these equations, there are no additional body forces aside from gravity. Radiative heat transfer is neglected. The flow is modeled as incompressible; hence, the pressure terms in the gas-phase energy equation are also neglected. The fluid-phase continuity equation is supplemented by an equation of state, where we use the ideal gas law with perfect species.
COMPUTATIONAL APPROACH Gas−particle flows in riser reactors occur over a wide range of time and length scales. To have a simulation that can cover this spectrum, an Eulerian−Eulerian (E-E) continuum and filtered equation approach is warranted. The two-fluid model (TFM) equations for gas−particle flows, which have been developed and analyzed extensively over the past five decades, are able to model these flows in a robust manner.21 These continuum partial differential equations (PDEs) model the multicomponent gases and solids through continuity, momentum, and energy conservation for each of the mixture-averaged gas and solid species. Interaction terms appear in the momentum and energy PDEs in addition to an equation of state for the gas mixture and a granular energy PDE for the solids. After an extensive evaluation of the computational tools available for this work, we selected the Multiphase Flow Interactions with Exchanges (MFIX) open source research software. MFIX solves the kinetic theory model equations with a finite-volume method and was extensively developed and documented over decades of extensive effort primarily by the U.S. Department of Energy, National Energy Technology Laboratory (NETL).22−24 Continuum computational fluid dynamics (CFD) models differ in accuracy and scale, starting with the more accurate direct numerical simulation (DNS) methods and ending with the more efficient, larger grid-scale methods like TFM.24 DNS methods resolve all fluid and solids scales, finding an exact solution of the coupled Navier−Stokes and particle Newtonian and collision dynamics equations. This is too computationally intensive for typically more than ∼1000 particles for engineer-
(2)
∂ ∂ (ϵg ρg ) + (ϵg ρg Ugi) = ∂t ∂xi
■
Ng
∑ R gn n=1
∂ ∂ (ϵmρm ) + (ϵmρm Umi) = 0 ∂t ∂xi
(1)
∂ ∂ (ϵg ρg Ugi) + (ϵg ρg UgjUgi) ∂t ∂xj = −ϵg
∂Pg ∂xi
+
∂τgij Θm ∂xj
− Igmi + ϵg ρg g
(3)
∂ ∂ (ϵmρm Umi) + (ϵmρm UmjUmi) ∂t ∂xj = −ϵm
ρg =
∂τmij ∂Pm − Igmi + ϵmρm g + ∂xj ∂xi
(4)
PgM w RTg
(5)
⎡ ∂Tg ∂Tg ⎤ ϵg ρg Cpg ⎢ + Ugi ⎥ ∂xi ⎦ ⎣ ∂t =−
⎡ ∂pg ∂pg ⎤ + γgm(Tm − Tg) − ϵg ⎢ + Ugi ⎥ − ΔHg ∂xi ∂xi ⎥⎦ ⎢⎣ ∂t
∂qg
(6)
⎡ ∂T ∂q ∂T ⎤ ϵmρg Cpm⎢ m + Umi m ⎥ = − m − γgm(Tm − Tg) ∂xi ∂xi ⎦ ⎣ ∂t 2848
(7)
DOI: 10.1021/acssuschemeng.6b02384 ACS Sustainable Chem. Eng. 2017, 5, 2847−2856
Research Article
ACS Sustainable Chemistry & Engineering where the subscripts “g” and “m” stand for gas and solid phases, respectively; ε stands for volume fraction; U stands for velocities; ρ stands for density; Rgn stands for rate of production of the nth chemical species in the fluid phase (this is not utilized for part 1 of this paper but is included here since it is used in part 2 of this work); P stands for pressure; t stands for time; T stands for thermodynamic temperature; R stands for universal gas constant; Mw stands for average molecular weight of gas; γgm stands for fluid−solids heat transfer coefficient corrected for interphase mass transfer; and ΔHg stands for heat of reaction in the fluid phase. In the gas-phase momentum balance equation, τgij is the gas-phase stress tensor; Igmi is the momentum transfer between the gas and solid phases, and ϵgρgg is the body force due to gravity. In the solid-phase momentum balance equation, τmij is the solid-phase m stress tensor (note for this simulation the catalyst phase is the only solid phase). The internal energy balances are presented here in terms of the temperatures, where
∂qg ∂xi
and referred to as the transport form of the granular energy (TGE): Pm = ϵsρm Θm[1 + 4ηϵmg0]
(8)
⎡ ∂Θ Umj∂Θm ⎤ 3 ϵmρm ⎢ m + ⎥ ∂xi ⎦ 2 ⎣ ∂t =
∂U ∂ ⎛ ∂Θm ⎞ ⎜Κ m ⎟ + τmij mi + Π m − ϵmρm Jm ∂xi ⎝ ∂xi ⎠ ∂xj
(9)
In the expression above, Jm is the collisional dissipation term given by M
Jm =
3/2
Θ 48 η(1 − η) ∑ (εng0,mn) m π dp n=1
1+e 2 and Πm is the exchange term given by η=
in the first term on the
right-hand side is the fluid-phase conductive heat flux, and the second term is the fluid−solids interphase heat transfer. Here, qm in the first term on the right-hand side of the solids energy equation is the solid-phase conductive heat flux, and the second term is the fluid−solids and interphase heat transfer. In the formulation of internal energy equations, the irreversible rate of increase of internal energy due to viscous dissipation and interphase momentum transfer has also been neglected. Einstein summation convention implied only on subscripts i and j. An additional PDE is required to solve the solids conservation of momentum as the solids stress tensor, τmij, and solids pressure, Pm, depend on the granular temperature, Θm (pseudothermal energy), of the solids. The form of the solids pressure and granular temperature PDE is shown below
Π m = −3β Θs +
(10) (11)
81εmμg2 |u g − u m|2 g0.mmd p3ρm π Θm
(12)
where β stands for the coefficient for the interphase force between the fluid phase and the solid phase, dp is the diameter of the particles, μg is the molecular viscosity of the vapor phase, and g is acceleration due to gravity. The implementation of the detailed granular energy equation in MFIX is still under development for simulations requiring complex boundary conditions with the Cartesian cut cell method. Therefore, all 3D simulations with nonrectangular domains use an algebraic approximation of the granular temperature and energy (AGE), which assumes that the granular energy is dissipated locally
2 2 2 2 1/2 ⎫2 ⎧ ⎪ − K1m ϵm Dmii + {K1m(Dmii ) ϵm + 4K 4m ϵm [K 2m(Dmii ) + 2K3m(Dmij Dmij )]} ⎪ ⎬ θm = ⎨ ⎪ ⎪ 2ϵmK4m ⎩ ⎭
where Dm stands for the rate of strain tensor for the solid phase, and the values for k1−4 are given in the MFIX documentation.28 In the simulations in this study, unless indicated, the Wen and Yu29 drag model was used, as the subgrid filters have been calibrated with this particular drag model, and therefore, the subgrid-scale (SGS) and non-SGS simulations’ accuracies were studied while using the equivalently designed drag models. The Gidaspow drag model30 is considered more accurate for flows containing high solids volume fractions and, therefore, was also used for comparison. A supercomputer31 that is available at the National Renewable Energy Laboratory (NREL) with ∼15,000 available cores and nodes with two 12-core Intel Xeon 2.4 GHz CPUs was utilized. The use of a hybrid model using both distributed (domain decomposition) and shared memory (threading of loops) parallelism was crucial for efficient, multidimensional, multicomponent simulations. More details on the use of the hybrid parallelization can be found in Syamlal et al.32 and Pannala et al.33 The multidimensional pseudocolor plots shown have been constructed with the open source VisIt tool.34 These numerical tools, by design, average the flow features smaller than the computational cell and do not represent single
(13)
particles. In the simulations, clusters correspond to the regions of highly dense solids. Johnson−Jackson Boundary Conditions. The Johnson−Jackson (J-J) boundary partial-slip condition is used in this study.35,36 This BC equates the tangential force per unit area exerted on the boundary by the flowing particles to the corresponding stress within the particles flow approaching the boundary. The force exerted on the boundary by particles is the sum of the collisional and frictional contributions during collisions. This is a function of the specularity coefficient, ϕp, which ranges from 0 to 1 (smooth to rough) and the wall coefficient of restitution. Here, ϕp = 0 is similar (but not equivalent) to a free-slip wall (FSW) for the solids. The J-J BCs appear as −n × qPT = D + usl × Scb
(14)
where usl is the slip velocity of the particles at the wall (given by usl = u − uw), n is the unit normal vector of the wall, qPT is the fluctuation energy flux, and D the dissipation rate of fluctuation energy given by 2849
DOI: 10.1021/acssuschemeng.6b02384 ACS Sustainable Chem. Eng. 2017, 5, 2847−2856
Research Article
ACS Sustainable Chemistry & Engineering D=
transient multidimensional wall-bounded flows that resolve the small inertial scales and wall effects are needed. To resolve the catalyst particle clusters and bubble-like voids at all length scales, extremely fine spatial resolution is necessary. This approach is very computationally intensive, for example, the 2D simulations of this study were requiring 1000 cores for several weeks of computer time. For a 3D simulation including the CFP kinetics for part two of this study, the extremely fine grid resolution approach would be computationally intractable. An alternative less computationally intensive approach would be to augment coarse multidimensional simulations with validated SGS models that capture the unresolved physics of gas−solids−wall interactions. Clusters are spatiotemporally evolving structures that would at any instant be asymmetrical with significant variations in the axial direction and thus would equivalently impact conversion and selectivity. These dynamics would in turn alter the stationary (time-averaged behavior) in terms of the core− annulus structure, developed flow region, RTD of gas and catalyst particles, contact times, and resulting conversion/ selectivity. Subgrid-Scale Models. We choose to use the Igci et al.49 SGS model in the context of the Johnson−Jackson solids BC.35 This approach allows modeling of the flow details with empirical correlations without requiring extremely fine resolution grids. The derivation of the filtered multiphase flow equations involves the separation of the resolvable filtered scales from the actual scales leaving the nonresolvable subgrid contributions, which is approximated with empirical correlations. Note that the best way to do this is still an open research problem. In Igci et al.49−51 and Milioli et al.52 for 2D, the subgrid residuals are found for nonreactive flows in the absence of walls by conducting empirical fits to fully resolved solutions with periodic domains yielding the subgrid core model. Particle/Catalyst RTD Computations. Calculation of the RTDs involves tracing particle trajectories. Because we are using a continuum and nonsteady modeling approach, this requires a special treatment. With an alternative Lagrangian particle method, for instance, the DEM, each particle is traced as it enters and exits the reactor. With a continuum Eulerian model, particle seed points can be designated (weighted by the local solids concentration) in space and time and tracked according to the solids velocity field as the simulation evolves. A consideration when calculating the RTDs in this manner is confirming that the simulation has reached a statistical stationary state (SSS) before the RTD is calculated. The approach to determining SSS was measured by calculating the variation of the mass inventory (mass accumulation of solids) in time. In all cases, this was found to be oscillatory, with an obvious mean value once an SSS was reached. In our analysis, we use the probability distribution (the standard RTD frequency plot) and the cumulative frequency distribution. We obtain the RTD, f(t), from the following relation
1 3Θ πρ d p3Θ(1 − e w2 ) 4 p d p[(αp,max /αp)1/3 − 1] ×
1 d p2(αp,max /αp)2/3
(15)
where ew is the normal restitution coefficient of particle−wall collisions, αp is the solids volume fraction, αp,max is the solids volume fraction at a closely random packing state, Θ is the granular temperature, and Sbc is stress exerted on the wall due to particle collisions given by Scb =
ϕ 3Θ πρp αp|usl| 6αp,max[1 − (αp/αp,max )1/3 ]
(16)
Combined with eq 14, the J-J BCs was implemented as ϕ 3Θ πρp αp|usl| usl × σc × n =0 + |usl| 6αp,max[1 − (αp/αp,max )1/3 ]
(17)
The solids boundary condition is thusly controlled by the parameters ϕp and ew and also depends on the local flow variables of solids volume fraction, density, granular temperature, and velocity. Clustering. In risers with Geldart type A particles,24,37 typically with a diameter of 20−100 μm, density of
