CFD–DEM modeling of gas–solid flow and catalytic MTO reaction in a fluidized bed reactor

Computers and Chemical Engineering 60 (2014) 1–16

Contents lists available at ScienceDirect

Computers and Chemical Engineering journal homepage: www.elsevier.com/locate/compchemeng

CFD–DEM modeling of gas–solid flow and catalytic MTO reaction in a fluidized bed reactor Ya-Qing Zhuang a , Xiao-Min Chen a , Zheng-Hong Luo b,∗ , Jie Xiao c a b c

Department of Chemical and Biochemical Engineering, College of Chemistry and Chemical Engineering, Xiamen University, Xiamen 361005, PR China Department of Chemical Engineering, School of Chemistry and Chemical Engineering, Shanghai Jiao Tong University, Shanghai 200240, PR China College of Chemistry, Chemical Engineering and Materials Science, Soochow University, 215123, PR China

a r t i c l e

i n f o

Article history: Received 15 February 2013 Received in revised form 12 August 2013 Accepted 14 August 2013 Available online xxx Keywords: Fluidization Granulation Multiphase reactors Mathematical modeling DEM MTO process

a b s t r a c t The methanol-to-olefins (MTO) process is currently being implemented successfully in fluidized bed reactors (FBRs) in China. Characterizing the gas–solid flow is crucial in operating MTO FBRs effectively. In this work, a combined discrete element method (DEM) and computational fluid dynamics (CFD) model is developed to describe the gas–solid flow behavior in an MTO FBR. In this model, the particles are modeled using DEM, and the gas is modeled using Navier–Stokes equations. The combined model incorporates the lumped kinetics in the gas phase to achieve the MTO process. Moreover, the combined model can characterize the heat transfer between particles as well as that between the gas and the particles. The distinct advantage of the combined model is that real-time particle activity can be calculated by tracking the motion history of the catalyst particle with respect to heat transfer. The simulation results effectively capture the major features of the MTO process in FBR. Moreover, the simulation results are in good agreement with the classical calculation and experimental data. The particle motion pattern and distributions of a number of key flow-field parameters in the reactor are analyzed based on the validated model. The effects of operating conditions on FBR performance are also investigated. The simulation results show that the particle motion exhibits a typical annulus–core structure, which promotes excellent transfer efficiency. The results also demonstrated that the feed temperature, inlet gas velocity, and feed ratio of water to methanol significantly affect reaction efficiency. © 2013 Elsevier Ltd. All rights reserved.

1. Introduction A computational fluid dynamics (CFD) model was recently proposed to describe the gas–solid flow in fixed-bed methanol-toolefins (MTO) reactors (Zhuang, Gao, Zhu, & Luo, 2012). Interesting flow behavior, such as coke deposition, as well as component and temperature distributions was observed in fixed-bed reactors. The coke deposition and the component distributions were also simulated in a fixed-bed reactor as a function of feed temperature, composition, and space velocity (Zhuang et al., 2012). The laboratory-scale reactor filled with catalyst particles was assumed to be a continuous porous medium. Given that the continuous Eulerian scheme is incapable of characterizing discrete solid particles (Wu, Cheng, Ding, & Jin, 2010a), some types of important solid flow behavior, such as particle motion pattern, temperature, and velocity, were not studied. Consequently, several important flow features in reactors, particularly in fluidized bed reactors (FBRs), such as excellent solid mixing, as well as heat and mass transfer

∗ Corresponding author. Tel.: +86 21 54745602; fax: +86 21 54745602. E-mail address: [email protected] (Z.-H. Luo). 0098-1354/$ – see front matter © 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compchemeng.2013.08.007

properties, were not analyzed (Zhu, Zhou, Yang, & Yu, 2008). Meanwhile, considering that the MTO reaction is highly exothermic, FBR is a more appropriate reactor for the MTO process compared with a fixed-bed reactor because of the excellent heat transfer capability of FBRs (Levenspiel, 2002). The FBR-based MTO process is currently being industrialized in China (Xing, Yue, Zhu, Lin, & Li, 2010; Zhu et al., 2010). To the best of our knowledge, the application of CFD to the MTO process in FBRs has not been reported because of the complexity of particle motion in FBRs. Particle motion, particularly its pattern, is reported crucial for CFD modeling (Chu, Wang, Xu, Chen, & Yu, 2011; Chu & Yu, 2008; Wu, Cheng, et al. 2010; Wu, Yan, Jin, & Cheng, 2010b; Zhao, Ding, Wu, & Cheng, 2010a; Zhao, Cheng, Wu, Ding, & Jin, 2010b). Thus, the effect of particle motion pattern on the flow-field parameters of the reactor requires further investigation to achieve the proper scaling up and design of reactors (Deen, Annaland, van der Hoef, & Kuipers, 2007; van der Hoef, Annaland, Deen, & Kuipers, 2008; Zhu et al., 2008). Understanding FBRs in the MTO process requires a fundamental knowledge of the gas–solid flow behavior while considering the motion of discrete particles in the gas phase. CFD approaches that can simulate gas–solid flow in FBRs are generally classified into two categories: Eulerian–Lagrangian and

2

Y.-Q. Zhuang et al. / Computers and Chemical Engineering 60 (2014) 1–16

Notation Cj,r Cp Cp,c Cp,g Cp,p CDS CAMF dp Di,m E fs F FD,i F /F n,ij

g hgs hi Hf,298 Hc Hj I kc keff kf km ks /kg kn /kt ki Kc kf,r Kfc mi Mi Mm Mw,i  ij n Nu NuD p Pr q Qgs,i Qi,k Qsg R S Sct Si rp Rep Ri Rˆ i,r t tos tij

t,ij

molar concentration of species j in reaction r, kmol/m3 heat capacity at constant pressure, J/kmol K heat capacity of molten salt, kJ/(kmol K) heat capacity of gas, kJ/(kmol K) heat capacity of particles, kJ/(kmol K) standard drag coefficient feed rate of methanol, gMeOH/gCat particle diameter, m diffusion coefficient for species i in gas mixture, m3 /s activation energy, kJ/kmol static friction coefficient interaction term between gas and particles, N/m3 drag force from gas to ith particle, N normal (tangential) contact force between the ith and jth particles, N gravitational acceleration, m/s2 solid-to-gas convective heat transfer coefficient, W/(m2 K) species enthalpy of formation, kJ/kmol enthalpy at 298 K, kJ/kmol contact conductance, W/K heat of reaction j, J/kg inertia of moment, kg m2 molten salt thermal conductivity, W/(m K) effective thermal conductivity of the medium, W/(m K) fluid phase thermal conductivity, W/(m K) mixture thermal conductivity, W/(m K) solid (gas) thermal conductivity, W/(m K) normal (tangential) spring constant, N/m reaction rate constant initial rate constant for coke formation forward rate constant for reaction r total heat transfer coefficient, W/(m2 K) mass of the ith catalyst particle, kg molecular weight, kg/mol molecular weight of gas mixture, kg/kmol molecular weight of species i, kg/kmol normal unit vector between two contacted particles Nusselt number in the tube Nusselt number outside the tube pressure, kPa Prandtl number heat flux, W/m2 heat transferred from gas to ith particle, W heat transferred from kth particle to ith particle, W heat source per unit volume in gas phase, W/m3 universal gas constant, kJ/kmol K mean strain rate, 1/s turbulent Schmidt number, dimensionless mass source of species j, kg/m3 s particle radius, m Reynolds number reaction rate, kg/(m3 s) arrhenius molar rate of creation of species i in reaction r, kmol/(m3 s) time, s time on stream, s tangential unit vector between two contacted particles

T Tg Tp,i Tij g u  p,i u V n,ij /V t,ij Vcell

vT

WHSV XA Xi Xw Yi ˛f ˛c ˇ ı n /ı t ε g c t ˜  i,r 

i,r n /t  j,r 

j,r b g s c ␭ ¯¯ ϕi ϕc

temperature, K gas mixture temperature, K temperature of ith particle, K torque due to tangential contact force, N m gas velocity, m/s ith particle velocity, m/s normal (tangential) relative particle velocity between particle i and j, m/s volume of computational cell, m3 transpose of velocity vector, m/s weight hourly space velocity, gMeOH/gCat h methanol mole fraction on a water-free basis outlet mole fraction of component i mole fraction of water in the feed mass fraction of species i fluid side convective heat transfer coefficient at wall, W/(m2 K) salt bath side convective heat transfer coefficient at wall, W/(m2 K) inter-phase momentum transfer coefficient, kg/(m3 s) normal (tangential) displacement vector, m turbulent kinetic energy dissipation, m2 /s3 viscosity of gas phase, Pa s viscosity of molten salt, Pa s turbulent viscosity, Pa s turbulent kinematic viscosity, m2 /s stoichiometric coefficient for reactant i in the rth reaction stoichiometric coefficient for product i in the rth reaction damping coefficient, N s/m rate exponent for reactant species j in the rth reaction rate exponent for product species j in the rth reaction bulk density of bed, kg/m3 density of gas mixture, kg/m3 density of solid, kg/m3 density of molten salt, kg/m3 thermal conductivity of steel, W/(m K) shear stress of gas phase, Pa medium porosity deactivation rate constant carbon deposition rate constant

Eulerian–Eulerian methods (Ranade, 2000; Lim, Zhang, & Wang, 2006; Utikar & Ranade, 2007; Vaishali, Roy, & Mills, 2008). Continuity and momentum equations are used in the Eulerian–Eulerian method, because this method considers the full interpenetration of continual subjects, suggesting that this method requires less numerical effort compared with the Eulerian–Lagrangian method, in which the motion of particles is calculated individually (Zhang, Lim, & Wang, 2007; Mahecha-Botero, Grace, Elnashaie, & Lim, 2006, 2009; Yan, Luo, Lu, & Chen, 2012). Given the high computational cost required to calculate the motion for a large number of particles, the application of Eulerian–Lagrangian method is restricted to low particle density. Most of the models proposed in the literature to describe gas–solid flows were based on the Eulerian–Eulerian method (Chen, Luo, Yan, Lu, & Ng, 2011; Gao, Chang, Lu, & Xu, 2008; Gao et al., 2009; Mahecha-Botero, Grace, Elnashaie, & Lim, 2006; Shi, Luo, & Guo, 2010). Although capable of

Y.-Q. Zhuang et al. / Computers and Chemical Engineering 60 (2014) 1–16

providing quantitative agreement with the limited experimental findings, these CFD models cannot physically describe the particle motion in the solid phase. Thus, the gas–particle, particle–particle, and particle–reactor wall interactions must be fully considered at the micro-scale to describe particle motion effectively (Xu, Ge, & Li, 2007; Zhou, Pinson, Zou, & Yu, 2011; Zhu, Zhou, Yang, & Yu, 2007; Zhu et al., 2008). The Eulerian–Lagrangian method is preferable for this endeavor because in this method, the particles are tracked individually by solving Newtonian equations of the motion, and the gas phase is treated as a continuum that is also coupled to the motion of particles through an interphase interaction term. The particle motion in FBRs has recently received considerable attention. Several models based on the Eulerian–Lagrangian method, which were pioneered by Tsuji, Kawaguchi, and Tanaka (1993), have been proposed to describe the gas–solid flows in FBRs. In addition, numerous groups (including Cheng’s group (Xu et al., 2007; Wu, Cheng, 2010; Wu, Yan, et al. 2010b; Zhao, Ding, et al., 2010; Zhao, Cheng, et al., 2010b), Li’s group (Ouyang & Li, 1999a, Ouyang & Li, 1999b), Yu’s group (Feng, Xu, Zhang, Yu, & Zulli, 2004; Zhang, Chu, Wei, & Yu, 2008)) have also made significant contributions (please refer to the reviews by Yu’s group (Zhu et al., 2007, 2008)). For instance, Zhang et al. (2008) used CFD and DEM to investigate the particle clustering behavior in a riser/downer reactor. Results showed that the existence of particle clusters, which were unique to the solid flow behavior in such reactor, can be predicted by the first-principles approach. Tsuji, Yabumoto, and Tanaka (2008) used CFD–DEM coupling simulation to investigate the flow structures induced by bubbles formed in 3D shallow rectangular gas-fluidized beds. A numerical code was parallelized, and over 4.5 million particles were tracked by 16 CPUs. Zhao et al. (2010a) simulated the gas–solid flows in a 2D downer (height: 10 m; width: 0.10 m) using CFD–DEM method, in which the motion of the particles was modeled using DEM, and the gas flow was modeled using by Navier–Stokes equations. The simulation results showed a rich variety of flow structures in the downer under different operating conditions. Geng and Che (2011) proposed a transient three-phase numerical model for simulating multiphase flow, heat and mass transfer and combustion in a bubbling fluidized bed of inert sand. The gas phase was treated as a continuum and solved using CFD approach, whereas the solid particles were treated as two discrete phases with different reactivity characteristics and solved using extended DEM based on individual particle scale. A new char combustion sub-model with inhibitory effects was also developed to describe the combustion behavior of char particles in FBRs. Karimi, Mansourpour, Mostoufi, and Sotudeh-Gharebagh (2011) simulated a gas-phase polyethylene reactor using CFD–DEM method. A comprehensive kinetic mechanism was used to evaluate the reaction rate of ethylene and 1-butene copolymerization. Li and Guenther (2012) applied an open-source MFIX-DEM code to simulate the effect of gas volume change resulting from chemical reactions on the flow hydrodynamics in bubbling FBRs. In their work, 2D simulations of ozone decomposition and the reverse reaction were conducted for systems with different particles. Previous efforts prove that CFD–DEM method can be used to describe particle motion and to capture particle activity. Although pressure drop evolution, particle mixing, gas and particle velocity vectors, and temperature with respect to fluidization time were characterized in previous studies by employing CFD–DEM method, the bubble behavior in FBRs was not revealed. To date, no attempt has been made to study the gas–solid flow in FBRs for the MTO process using CFD–DEM method. Reports on the MTO process are largely available (Aguayo, Mier, Gayubo, Gamero, & Bilbao, 2010; Alwahabi & Froment, 2004; Chen, Rebo, Gronvold, Moljord, & Holmen, 2000; Chen, Gronvold, Moljord, & Holmen, 2007; Gayubo, Benito, Aguayo, Castilla, & Bilbao, 1996; Gayubo, Aguayo, del Campo, Tarrio, & Bilbao, 2000;

3

Lwahabi & Forment, 2004; Park & Froment, 2001a, Park & Froment, 2001b). However, these reports focused on the catalyst and the kinetics in the MTO process. For instance, Schoenfelder, Hinderer, Werther, and Keil (1995) developed a lumped kinetic model which was then incorporated into a reactor model of circulating fluidized bed (CFB) for the MTO process. Soundararajan, Dalai, and Berruti (2001) also simulated the MTO process in a CFB reactor. In the simulation, the kinetic model with SAPO-34 as catalyst was combined with a core-annulus type hydrodynamic model. Fatourehchi, Sohrabi, Royaee, and Mirarefin (2011) investigated the MTO reaction in a differential fixed-bed reactor using acidic SAPO-34 molecular-sieve as catalyst. In addition, a mechanism for this process based on Langmuir–Hinshelwood formulation was also proposed, and the kinetic parameters were evaluated as functions of temperature. Mier, Gayubo, Aguayo, Olazar, and Bilbao (2011) studied the joint transformation of methanol and n-butane fed into a fixed-bed reactor on an HZSM-5 zeolite catalyst. In their work, the proposed kinetic scheme of lumps integrated the reaction steps corresponding to the individual reactions (cracking of n-butane and MTO process at high temperature) and considered the synergies between the steps of both reactions. These reactors consist of both CFB and fixed-bed reactors. However, no reports were found on the solid-solid flow in the MTO FBRs based on the CFD–DEM method. Nevertheless, related studies can still provide some insight into this work. In this study, a combined CFD–DEM model, which incorporates the heat transfer between the particles and between the gas and the particles, as well the lumped kinetics in the gas phase for the MTO process, is applied to study the flow behavior in an MTO FBR. The distinctive advantage of the proposed model is that the realtime particle activity can be calculated by tracking the history of the particle motion with respect to heat transfer. This presumed advantage is validated by classical calculation and experimental data. The combined model is also used to identify the gas–solid flow field, particularly the particle flow structure and the MTO reaction efficiency in MTO FBR. Simulations are then performed to understand the effects of the operating conditions on MTO FBR performance. 2. Mathematical model and numerical simulations Under the Eulerian-Lagrangian scheme, the governing equations of the gas–solid flows include the Navier–Stokes equations for gas flow and the DEM for particle motion, thereby comprising the CFD–DEM model. Given the wide applications of this approach in multiphase flow, the CFD–DEM model is well documented in the literature (Hoomans, Kuipers, Briels, & van Swaaij, 1996; Xu & Yu, 1997; Ouyang and Li, 1999a,b; Zhou et al., 2011). For brevity, only a number of core governing equations, such as the equations for gasphase flow, particle motion, and the interaction between gas and particles, are discussed and listed in Table 1 For example, the gas phase is modeled by the conservation equations of mass, momentum, energy, species, and the turbulence model (i.e., Realizable k − ε turbulent model (Shih, Liou, Shabbir, Yang, & Zhu, 1995)); the particle movement is described by the translational and rotational motions of a particle at any time, which follows the Newton’s laws of motion. The equation shown in Table 1 on particle movement is presented as follows mi

 p,i du dt

= FD,i + mi g +

ni 

(Fn,ij + Ft,ij )

(1)

j=1

In practice, the instantaneous particle velocity is calculated using Eq. (1), which is determined by the drag force (FD,i ), the gravity (mi g ), and the contact force (Fn,ij , Ft,ij ). In addition, the gravity and the drag force act on the mass center of the ith particle,

4

Y.-Q. Zhuang et al. / Computers and Chemical Engineering 60 (2014) 1–16

Table 1 Main governing equations of the DEM–CFD model. Governing equations for gas phase flow g ) = 0 Gas phase continuity equation: ∂∂t ( g ) + ∇ · ( g u  g ) + ∇ · ( g u g u  g ) = −∇ p − F + ∇ · ( g ) + g g Gas phase momentum equation: ∂∂t ( g u Gas phase equation of state: p = g RTg /Mm  g T ) = ∇ · ( kg ∇ T ) + Qsg Gas phase energy equation: ∂∂t ( g Cp,g T ) + ∇ · ( g Cp,g u   g Yi ) = −∇ · ( Ji ) + Ri , i ∈ [1, Nr ], where, Ji = − g Di,m + Species conservation equation: ∂∂t ( g Yi ) + ∇ · ( g u

⎡ ⎞⎤ ⎛  + j,r ( ) Nr Nr j,r    

⎢  −  ⎜k ⎟⎥ Mass rate of reaction: Ri = Mw,i Cj,r ⎣ i,r i,r ⎝ f,r ⎠⎦ r=1

d up,i

= FD,i + mi g +

dt

ni 

(Fn,ij + Ft,ij ) and I i

i dw dt

=

ni 

j=1

n

−

i=1

FD,i

Vcell

and FD,i = FD0,i −(ˇ+1) , FD0 =

   g − f dp u up,i 

Re =

 Tij where, Fn,ij = −kn ı n − n V n,ij , Ft,ij =

1  C 2 f DS

2

dp

−kt ı t − t V t,ij

  −fs Ft,ij  tij

4





ij

ij

t,ij

ij

from the mass center of the particle to the contact point. Thus, the equation governing the rotational motion of the ith particle is

 i dw Tij , = dt



Re ≤ 1

24/Re (0.63 + 4.8Re−1/2 )

2

, where,

Re > 1

the gas phase to the ith particle (Ranz, 1952; Kaneko, Shiojima, & Horio, 1999): Qgs,i = dp2 hgs (Tk − Ts,i ),

(2)

j=1 →

where wi is the angular velocity, and Ii is the moment of inertia of the ith particle that can be obtained using Ii = (2/3)mi rp2 . Thus, the magnitude and the vector of the velocity of the ith particle are determined by Eqs. (1) and (2). The equations for calculating these forces are also listed in Table 1 (Di Felice, 1994). For the gas–solid MTO reacting flow, the hydrodynamic model incorporates the sub-models describing the heat transfer and the chemical reactions, which are presented in the following sections. 2.1. Heat transfer sub-model Borkink and Westerterp (1992) proposed a seven-step mechanism to describe intraparticle heat transfer. Considering the physics in a fluid catalytic cracking (FCC) process, Wu et al. (2010a) used two of the seven mechanisms (i.e., the thermal conduction through the contact area between two particles and the heat transfer by the fluid convection) to construct the heat transfer sub-model, which was then incorporated into the CFD–DEM model to simulate the gas–solid reacting flows in the FCC processes. Given that the intraparticle heat transfer characteristic of the MTO process is similar to that of the FCC process, the heat transfer sub-model suggested by Wu et al. (2010a) is also adopted in this work. The main equations for the heat transfer sub-model are listed as follows The heat flow between the ith and kth particles can be calculated as follows (Vargas & McCarthy, 2001): Qik = HC (Tk − Ti ),

(3)

where the contact conductance can be expressed as follows (Vargas & McCarthy, 2001):

3Fn ·r ∗ 1/3 4E ∗

= 2˛.

(5)

where hgs is the heat transfer coefficient than can be obtained using Eq. (6):

ni

=2

 ij × Ft,ij Tij = rp n

g

ij

HC ks

         Ft,ij ≤ fs Ft,ij  Ft,ij  > fs Ft,ij 

 p,i  (u g − u  p,i ), ˇ = 3.7 − 0.65exp[− (1.5 − lgRe)2 /2], CDS = g − u 2 u

whereas the contact force acts at the contact point between the ith and jth particles, thereby generating a torque (Tij ) and causing the ith particle to rotate. For a spherical particle with the radius →  × F , where n is a vector running of rp , T is given by T = rp n

Ii

∇ Yi

j=1

Interaction between gas and particle, F៝ F៝ =



j=1

Governing equations for particle motion mi

t Sct

hgs · dp Nu = = kg



0.03Re1.3 p

Rep ≤ 100

1/2 2.0 + 0.6Rep Pr 1/3

Rep > 100





,

(6)

 p,i  /g , and Pr = Cp,g g /kg . g − u where Rep = g dp u In addition, the volumetric heat source in a computational cell of the gas phase can be expressed as follows: −

Qsg =

n i=1

Qgs,i

Vcell

.

(7)

Meanwhile, the particle–particle and particle–gas heat transfers determine the temperature of each particle. The temperature of each particle is also determined by the heat released from the reaction because the particles are used as catalysts. mi Cp,p

dTp,i dt

=



Qik + Qgs,i .

(8)

k

The temperature within each particle is assumed to be uniform. This assumption is valid based on the Biot number (i.e., Bi), which represents the ratio of the intraparticle heat transfer resistance to the external heat transfer resistance of the particles (Deliang, Watson, & McCarthy, 2008). The rate of convective heat transfer at the inlet conditions was estimated as follows: First, Rep was obtained using Eq. (9): Rep =

g dp uinlet 0.54 × 1.0 × 10−3 × 2.8 = = 8.79 < 100, g 1.72 × 10−4

(9)

Therefore, the convective heat transfer coefficient can be expressed as follows: Nu = 0.03Re1.3 p = 0.51,

(10)

(4)

For the particle–gas convective heat transfer, the following equations are used to calculate the convective heat transfer from

hgs =

Nu · kg 0.51 × 0.0254 = = 12.95 W/(m2 K), dp 1 × 10−3

(11)

Y.-Q. Zhuang et al. / Computers and Chemical Engineering 60 (2014) 1–16

5

Fig. 1. Lumped-species reaction scheme: ((a) detailed scheme and (b) simplified scheme)).

Therefore, Bi =

hgs · dp 12.95 × 1 × 10 = 0.156 ks

−3

= 8.3 × 10−2 << 1,

(12)

In addition, the value of Bi based on the solid conduction can be determined using the following equation: Bi =

HC HC 2˛ << 1. = =

r ks r ks A/r

(13)

Accordingly, the sum of Bi that incorporates convective and conductive heat transfers is significantly lower than 1, suggesting that the intraparticle heat transfer resistance is much smaller than the external heat transfer resistance of the particles. This finding proves that the temperature within one particle is uniform. According to the above equations, the heat released from the reaction causes the higher temperature of the particles, which then heat up the surrounding gases. Thus, the temperature of the gases and the particles is increased by the reaction. 2.2. MTO reactions sub-model Chen, Gronvold, Moljord, Fuglerud, and Holmen (1999), Chen, Rebo, Moljord, and Holmen (1999), Chen et al. (2000, 2007) reported an elaborate coking rate equation, which showed good agreement with the actual data. Therefore, Chen et al.’s kinetic model is applied in this work. Similar to Chen et al.’s work, SAPO34 is used as the MTO reaction catalyst in this study. The adopted lumped-reaction scheme for the MTO process is shown in Fig. 1a. C5 is assumed to represent C5 and C6 together because of the low fraction of both C5 and C6 as macromolecular compound and byproduct in the lumped-species scheme. In addition, methane is neglected because of its low concentration. Ethane and propane are selected to represent paraffin. Therefore, the reaction network shown in Fig. 1a can be simplified as Fig. 1b (Zhuang et al., 2012). The reaction model is given by → dXi = Ri . d(W/FM0 )

(14)

Moreover, Chen et al.’s kinetic equations for the reaction network, which are presented in Fig. 1b, are also adopted in this study

and presented in Table 2, where the subscripts 1 to 5 represent ethene, propene, butene, C5 , and (ethane + propane), respectively. The coke formation equation is listed as follows: C=

1 [1 − exp(−ϕc kc CAMF(1 − XA ))], ϕc

Kc = 1.27 exp

 −6707  8.314T

,

(15)

(16)

ϕc = 1.07T 2 − 0.00181T + 0.7954,

(17)

CAMF = WHSV × tos .

(18)

2.3. Numerical simulation The CFD–DEM coupled model was successfully solved using the commercial CFD code FLUENT 6.3.26 (Ansys Inc., US). The equations of gas flow were solved using the SIMPLER algorithm and the DEM model was solved using the finite difference algorithm. In addition, the DEM model and the other sub-models were coded with the C++ language as their user defined functions (UDFs), and the UDFs were implemented in FLUENT. Meanwhile, the velocity of the particle was computed by the DEM model, and the gas velocity at the position of the corresponding particle was computed by the FLUENT solver. The velocities of particle and gas were compared. If the velocity of the gas blown from the bottom of the reactor was greater than that of the particle, then the particle moved toward the top of the reactor, and vice versa. Then, the positions, the velocities and the other data related to the particles were expressed using the Tecplot software. Thus, not only the particles trajectories, but also the temperature and the components of the gas surrounding the particles, could be investigated using the CFD–DEM coupled model. All simulations were performed in 2D domains for a laboratory scale FBR (width: 0.16 m; height: 0.8 m). The computational conditions and the additional parameters are given in Table 3. In addition, model parameters related to the gas and catalyst properties are listed in Tables 4 and 5 (Reid, Prausnitz, & Poling, 1987; Shackelford, Alexande, & Pork, 2002). These parameters, expect otherwise noted,

Table 2 The kinetic equations and their parameters (Chen, Gronvold, Moljord, Fuglerud, & Holmen, 1999; Chen, Rebo, Moljord, & Holmen, 1999; Chen et al., 2000, 2007). Kinetic expression (ri )

Reaction rate constant (ki /(mol g−1 kPa−1 s−1 ))

Deactivation rate constant (ϕi )

r1 = k1 (1 − ϕ1 C)pA r2 = k2 (1 − ϕ2 C)pA r3 = k3 (1 − ϕ3 C)pA r4 = k4 (1 − ϕ4 C)pA r5 = k5 (1 − ϕ5 C)pA

k1 = 209 · exp(− 38, 400/RT) k2 = 40 · exp(− 27, 000/RT) k3 = 15 · exp(− 26, 900/RT) k4 = 17 · exp(− 49, 800/RT) k5 = 181 · exp(− 59, 600/RT)

ϕ1 = −3 × 10−7 · T2 + 0.0006 · T − 0.2365 ϕ2 = −6 × 10−7 · T2 + 0.0011 · T − 0.4143 ϕ3 = −9 × 10−7 · T2 + 0.0014 · T − 0.5148 ϕ4 = −8 × 10−7 · T2 + 0.0013 · T − 0.44532 ϕ5 = 2 × 10−6 · T2 − 0.0032 · T + 1.247

6

Y.-Q. Zhuang et al. / Computers and Chemical Engineering 60 (2014) 1–16

Table 3 Computation conditions and additional parameters. Parameters

Table 5 Model parameters (Shackelford et al., 2002; Zhuang et al., 2012). Value

Particles Shape Diameter (m) Density (kg/m3 ) Initial temperature of particle (K) Number of particles DEM–CFD model parameters Particle normal stiffness coefficient (N/m) Particle tangential stiffness coefficient (N/m) Particle normal damping coefficient Particle tangential damping coefficient Time step for particle phase (s) Time step for gas phase (s) Computation mesh for gas phase (m × m)

Sphere 1.0 × 10−3 640 300 8000 ∼ 40,000 20,000 10,000 0.5 0.3 5.0 × 10−7 2.5 × 10−6 0.004 × 0.004

Descriptions

Value

Thermodynamic and physical parameters Gas mixture Cp,g (kJ/(kmol K)) g (kg/m3 ) g (Pa s) kg (W/(m K)) Di,m (m2 /s) Solid phase CP,s (kJ/(kmol K)) ks (W/(m K))

“Mixing law” “Idea gas” 1.72 × 10−4 0.0254 2.393 × 10−5 774 0.156

were used in the following simulations. Three boundary conditions: (1) the uniform mass flow inlet, (2) the pressure outlet, and (3) the non-slip and zero heat flux for the wall, were applied in the simulations. The under-relaxation factors are set to 0.9 for the equations of all species to achieve better convergence. A two-stage calculation was implemented. First, the flow field was simulated without the MTO reaction until a fully fluidized flow field was obtained under a cold-flow condition. Thereafter, the reaction process was triggered by incorporating the DEM part and other sub-models.

3. Results and discussion In this section, the grid sensitivity study for the gas phase model is first investigated. The coupled model is then verified using the data calculated from the classical methods (see Eqs. (14)–(16)) and the experimental data. A number of basic particle flow characteristics are predicted using the coupled model without the reaction sub-models (i.e., under cold-flow conditions). In addition, the main reacting flow behavior is also obtained using the complete model.

3.1. Grid sensitivity study Simple grid sensitivity was conducted to determine the minimum number of cells for the gas phase to reach the convergence of the model. Five cases were designed for the grid sensitivity analysis. The simulation domain for the FBR was divided into 12,800, 19,200, 25,600, 32,000 and 38400 cells, respectively. The effect of cell number on the radial velocity of the gas at a height of 0.1 m was tested in the simulation. As shown in Fig. 2, the radial velocity converged when the number of cells exceeded 12,800. Accordingly, 19200 cells were selected for the gas phase in the following simulations. The instantaneous radial gas velocity data were plotted against the width, as shown in Fig. 2, and the data show asymmetry.

Fig. 2. Radial velocity distribution using different grids.

3.2. Particle flow characteristics 3.2.1. Pressure drop and model verification under cold-flow conditions Pressure drop, which can reflect particle flow behavior, is one of the most important parameters to achieve proper scaling up and reactor design. The pressure drop in a fluidized bed can always be described by the buoyant weight of the suspension, which can be expressed as follows (Chai & Zhang, 2004; Lettieri, Felice, Pacciani, & Owoyemi, 2006; Shi et al., 2010): Ps = (s − g )(1 − )gL.

(19)

In this work, the gas phase density is 1.4 kg/m3 . The effect of gas phase weight on the pressure drop is described by Pg = g gL.

(20)

The total pressure drop of FBR (i.e., the pressure drop from the FBR inlet to its outlet) can be calculated using Eq. (20). Pt = Ps + P g .

(21)

Table 4 The thermodynamic data of the components in the reaction system (Reid et al., 1987). Components

Hf,298 × 10−5 (J mol−1 )

Gf,298 × 10−5 (J mol−1 )

CH3 OH (g) C2 H4 (g) C3 H6 (g) C4 H8 (g) H2 O (g) C2 H6 (g) + C3 H8 (g) C5 (g)

−2.013 0.523 0.204 0.00126 −2.420 −0.847 −1.262

−1.618 0.682 0.628 0.716 −2.286 −0.328 −0.17

Cp = A + B · T + C · T2 + D · T3 (j mol−1 K−1 ) A

B × 102

C × 105

D × 108

21.150 3.806 3.710 −2.994 32.240 19.250 90.487

7.092 15.660 23.450 35.320 0.192 5.409 33.13

2.587 −8.348 −11.600 −19.900 1.055 17.800 −11.08

−2.852 1.755 2.205 4.463 −0.360 −6.938 −0.282

Y.-Q. Zhuang et al. / Computers and Chemical Engineering 60 (2014) 1–16

Fig. 3. The bed pressure drop as fluidization proceeds (space velocity = 2.8 m/s).

The corresponding pressure drop data calculated by these equations are shown in Fig. 3, which shows the bed pressure drop profile as a function of flow time in FBR. As shown in Fig. 3, the calculated values of the bed pressure drop are slightly higher than the data calculated from classical equations. This discrepancy may be caused by the disregard given to the pressure drop resulting from friction and the particle collision in the classical calculation. Nevertheless, Fig. 3 proves that the simulated pressure drop data are in agreement with the classical data. Moreover, as shown Fig. 3, three typical regions in the MTO FBR at cold-flow conditions can be identified before the formation of a stable fluidization state; these regions consist of the bed expansion and slug formation stage (t ≤ 0.2 s), the gas–solid fluidization establishment stage (0.2 s < t < 0.5 s), and the stable fluidization stage (t ≥ 0.5 s). The maximum bed pressure drop in the first stage is higher than those in the other stages because the interparticle locking is overcome. Evident bed expansion and slugging can be observed at this stage. After the bed expansion and slugging, the interparticle locking is broken. Hence, the fluctuation of the bed pressure drop shifts from an irregular to a regular pattern (Fig. 3). At the third stage (t ≥ 0.5 s), the bed pressure drop fluctuates regularly and the stable gas–solid fluidization state is then established. The gas–solid fluidization in FBRs is known to be significantly affected by the inlet gas velocity, which causes the existence of different pressure drop profiles. According to a well-known fluidization theory (Jin, Zhu, Wang, & Yu, 2001), there exist four typical fluidization states that correspond to different gas velocities in a typical FBR. The effect and the four states can be captured by the current numerical model. As shown in Fig. 4, the four typical states in FBR can be identified as follows: fixed-bed state (from A

Fig. 4. The effect of the space velocity on the bed pressure drop.

7

to B, A → B), complete fluidization state (from C to D, C → D), initial fluidization state (from C to E, C → E), and fluidization vanishing state (from E to F, E → F). At the fixed-bed state, the pressure drop increases with the increase in inlet gas velocity. In practice, particle interactions result in a particle bridge and particle scarfing. Accordingly, a greater driving force is needed to loosen the bed, that is, to form a high bed pressure drop. At the maximum pressure drop (i.e., position B in Fig. 4), the bed begins to loosen and its voidage starts to increase. At the complete fluidization state (C → D), continuous bubbles come into being, and the bed pressure drop is kept unchanged. With the decrease in the inlet gas velocity (C → E), the fluidization degree in FBR and the pressure drop decrease, and the gas velocity reaches the minimum/initial fluidization velocity, which is approximately 1.8 m/s. When the gas velocity is lower than 1.8 m/s (E → F), the particles in the bed cannot be fluidized, that is, the bed is at the fluidization vanishing state. The comparisons are only qualitative, because the present simulations were not conducted under the same conditions as in the theory (Jin et al., 2001).

Fig. 5. Particle configurations at different fluidization time positions (space velocity = 2.8 m/s).

8

Y.-Q. Zhuang et al. / Computers and Chemical Engineering 60 (2014) 1–16

Fig. 6. Velocity vector profiles at t = 0.312 s ((a) particle-phase velocity vector profile and (b) gas-phase velocity vector profile).

The information shown in Figs. 3 and 4 suggest that the developed CFD–DEM model can capture the major particle flow characteristics in an MTO FBR at cold-flow conditions. The pressure drop analysis, together with detailed information provided in the following sections, can result in a better understanding of the flow structure and the mechanisms in MTO FBR. 3.2.2. Gas–solid flow pattern As previously described, the gas–solid flow pattern in the gas phase, particularly the solid flow pattern, cannot be obtained using the Eulerian–Eulerian method (Mahecha-Botero et al., 2009; Yan et al., 2012; Zhang et al., 2007). However, the current CFD–DEM model can successfully capture the gas–solid flow pattern, as shown in Fig. 5. The different particles are labeled using different

colors at 0 s to distinguish these particles and to investigate their mixing (Fig. 5). During the bed-expansion and slug-formation stage (t ≤ 0.2 s, see Fig. 5), these particles move upward, and the spouted jet zone forms near the FBR inlet because of the drag acting on these particles at the bottom of FBR. Meanwhile, a number of top particles in FBR move downward to the bottom to fill the voids in FBR (Fig. 5a–c). In addition, at this stage, the bubbles come into being and develop over time. Accordingly, the bed expands and its height increases over time. After the first stage, the process proceeds to the second stage (i.e., the gas–solid fluidization establishment stage, 0.2 s < t < 0.5 s, Fig. 3), these bubbles aggregate to form two big bubbles, which move toward the two walls of the FBR (Fig. 5d). These particles near the bottom of FBR are blown upward and a convex forms with respect to time (Fig. 5d–f). At the stable

Fig. 7. Velocity vector profiles at t = 0.52 s ((a) particle-phase velocity vector profile and (b) gas-phase velocity vector profile)).

Y.-Q. Zhuang et al. / Computers and Chemical Engineering 60 (2014) 1–16

9

Fig. 8. Velocity vector profiles at t = 1.04 s ((a) particle-phase velocity vector profile and (b) gas-phase velocity vector profile).

fluidization stage (t ≥ 0.5 s, Fig. 3), the films between two bubbles become thinner over time (Fig. 5g–i) and ultimately rupture. These gases in the bubble move circularly from left to right, causing the particles adjacent to the bubble to move. Furthermore, at t = 1.3 s (Fig. 5m–r), the particles in FBR move circularly around the bubble under the influence of the inlet gas. In summary, Fig. 5 proves that the particles in the middle of FBR move upward, and the particles near the two walls move downward, generating a typical annulus-core structure in FBR. Strong particle back mixing and internal circulation exist for this flow structure in FBR. This type of particle flow pattern promotes the thorough mixing of particles, resulting in high heat transfer efficiency. 3.2.3. Gas–solid velocity vector The gas and the solid phase velocity profiles at three typical time instants are obtained using the proposed model to verify the flow pattern further (Figs. 6–8). Fig. 6 shows that the particles in the middle of FBR move up, and a convex forms because of the blowing up of the inlet gas along the FBR axis at 0.312 s (Fig. 6a). This phenomenon is identical to that observed in Fig. 5e. Meanwhile, the gas velocity in the axial middle of FBR decreases with the increase in the bed height (Fig. 6b). In addition, the particles at the top of FBR are initially thrown upward and subsequently fall down toward the two walls. Given that the gas velocity near the two walls is lower than that in the middle, the fallen particles along the two walls must move down to the dense phase zone of particles. These fallen particles, as well as the other particles in the dense phase zone, then move slowly downward until these particles are rolled into the spouted jet zone and eventually sprayed out again. Therefore, a typical annulus-core structure for the solid particles in FBR as described in Section 3.2.2 is formed, as shown in Fig. 6a. Fig. 6b describes the velocity vector of the gas phase in FBR at 0.312 s. As shown in Fig. 6b, some gases near the bottom of the FBR diffuse toward the two walls and form two small bubbles, with the gas flowing into FBR (Fig. 5e). This phenomenon can be clearly observed on the magnified images. As shown in the bottom-right image in Fig. 6b, one portion of the gases move circularly around the bubble, resulting in the circular movement of

some particles in the FBR surrounding the bubbles. Meanwhile, another portion of the gases, along with the entering particles, move upward from the gap between the two bubbles, as illustrated in the local magnified image in Fig. 6b. Fig. 7 describes the gas and solid velocity vector profiles at 0.52 s. The annulus-core structure for the solid particles in FBR can still be observed. However, the right bubble becomes larger and would ultimately break (Figs. 5g and 7b), suggesting that the particle movement is primarily determined by the entrainment of the gas. The gas and solid velocity vector profiles at 2.6 s are listed in Fig. 8, which corresponds to the gas–solid flow pattern shown in Fig. 5r. As shown in Fig. 8a, a large bubble can be formed, and the particles move circularly around the bubble. The same result is also obtained from the gas velocity vector profile, as shown in Fig. 8b. In summary, the information shown in Figs. 6–8 is consistent with that presented in Fig. 5, further proving that a typical annuluscore structure in the MTO FBR can be formed. 3.2.4. Particle trajectory Compared with the Eulerian method, DEM can track the history of the particle motion in real time. In this study, the particle trajectory in the fluidization process is simulated using the coupled model, as shown in Fig. 9. The selected particle is one of the particles entering FBR at 0 s. The particle temperature is described in Fig. 9, because the heat transfer sub-model is incorporated into the coupled model. In addition, the particle temperature profile with respect to time is shown in Fig. 10. As shown in Fig. 9, the particle is first blown up by the inlet gas. The particle then moves circularly from left to right in FBR as the fluidization proceeds. Meanwhile, the particle is heated up continuously, resulting in the increase in particle temperature over time. However, based on the local magnified image on the right side of Fig. 9, the particle temperature is always low when the particle moves across the right bottom of FBR because of the low gas temperature at this position. Fig. 10 shows that the particle temperature increases quickly and reaches a plateau value, which is the gas temperature. In practice, the fluctuation of particle temperature is due to the circular movement of the particle, as shown in Fig. 9. In summary, the information shown in Figs. 9 and 10 suggest that the

10

Y.-Q. Zhuang et al. / Computers and Chemical Engineering 60 (2014) 1–16

Fig. 9. Single-particle trajectory profile (space velocity = 2.8 m/s).

developed CFD–DEM can also capture the particle trajectory when the heat transfers between the particles and between the gas and the particles are considered. Furthermore, Figs. 9 and 10 demonstrate that excellent heat transfer efficiency can be achieved in the FBR. 3.3. Reacting flow behavior 3.3.1. Reacting flow field distribution and model verification During the preliminary work (Zhuang, 2012b), methanol conversion was found to reach approximately 100% during the initial fluidization stage (the simulated profiles are not listed because of the limitation in space). Given that the SAPO-34 catalyst is easily deactivated, the gas and particle temperatures, the coke content, and the product concentrations in the reactor would change as the reaction and the fluidization proceed. Therefore, the distributions of these parameters in the reactor should be investigated. In this study, four typical time instants during the process (0.052, 0.312, 0.520, and 2.600 s) are selected as representatives, which correspond to the simulations mentioned in Sections 3.2.1.

Fig. 10. Single-particle temperature profile (space velocity = 2.8 m/s).

Fig. 11. Main reaction parameter distribution profiles in the MTO FBR at t = 0.052 s ((a) gas-phase temperature, (b) particle temperature, (c) coke content, (d) ethene mole concentration, (e) propene mole concentration, and (f) butene mole concentration; space velocity = 2.8 m/s, inlet feed temperature = 723 K, feed ratio of water to methanol = 0).

Fig. 11 describes these parameter distributions in FBR at 0.052 s. Meanwhile, according to Fig. 3, the reacting flow is in the bed expansion and the slug formation stage (t ≤ 0.2 s). Given that the inlet raw gases have the highest concentration, when the raw gases come in contact with the particles, the particles are heated up by the exothermic heat released from the reaction. The gases surrounding the particles are then heated up by the heat transfer from the particles. Therefore, the temperature of the gases and the particles are highest at the inlet. The gas with the highest temperature diffuses upward to both sides, as shown in Fig. 11a. Accordingly, as shown in Fig. 11b and c, both the particle temperature and the coke content profiles exhibit a distribution similar to that of the gas temperature in FBR because of the close positive relations among the three parameters, as shown in Fig. 11b and c. Meanwhile, product concentrations (Fig. 11d–f).exhibit distributions similar to those of temperature profiles, thus proving that the change in gas/particle temperature is determined by the exothermic heat of the MTO reaction. The parameter distributions at 0.208 s described in Fig. 12 suggest that the reacting flow is in the gas–solid fluidization

Y.-Q. Zhuang et al. / Computers and Chemical Engineering 60 (2014) 1–16

Fig. 12. Main reaction parameter distribution profiles in the MTO FBR at t = 0.208 s ((a) gas-phase temperature, (b) particle temperature, (c) coke content, (d) ethene mole concentration, (e) propene mole concentration, and (f) butene mole concentration; the space velocity = 2.8 m/s, inlet feed temperature = 723 K, feed ratio of water to methanol = 0).

establishment stage (0.2 s < t < 0.5 s, Fig. 3). As shown in Fig. 12a, the gas with the highest temperature continuously diffuses toward the outlet of FBR. Given that the total temperature in FBR is primarily controlled by the reaction heat, the zone with the highest gas temperature should be the place stacked with the catalyst particles (Fig. 5d). Meanwhile, the gas temperature in the zone with the bubbles (i.e., the “bubble” zone) is lower, as shown in Figs. 5d and 12a. For the particle temperature, the highest temperature is still at the gas inlet of the bed, as shown in Fig. 12b. However, a relatively uniform temperature distribution exists in the other zones in FBR except in the “bubble” zones. Meanwhile, the coke content distribution is similar to the particle temperature distribution described in Fig. 12c, that is, the coke contents in the “bubble” zones of the bed are low. However, the total coke content increases over time based on comparison between Figs. 11c and 12c. The products also continuously diffuse along the axial outlet direction in the bed as the fluidization and the reaction proceed (Figs. 12d–f). Meanwhile, the product concentrations in the “bubble” zones of the bed remain low.

11

Fig. 13. Main reaction parameter distribution profiles in the MTO FBR at t = 0.52 s ((a) gas-phase temperature, (b) particle temperature, (c) coke content, (d) ethene mole concentration, (e) propene mole concentration, and (f) butene mole concentration; space velocity = 2.8 m/s, inlet feed temperature = 723 K, feed ratio of water to methanol = 0).

At the stable fluidization stage (t ≥ 0.5 s), the bubbles in FBR undergo breakage and aggregation to form a big bubble. In addition, the particles are at the entire fluidization state (Fig. 3). In this study, two series of parameter distribution profiles at two time instants (0.52 s and 2.6 s) during this period are presented in Figs. 13 and 14. At 0.52 s, the left bubble breaks, and the fluidization in the bed is just at the entire fluidization state. The highest gas temperature in the “bubble” zone is approximately 750 K, which is lower than that at 0.208 s (Fig. 12) because of the increase in coke content during the reaction and the fluidization. The products diffuse into the entire FBR, and all their concentrations are uniform because of the entire fluidization. During the fluidization at 2.6 s, both the gas and solid temperatures in the bed are lower than those at 0.52 s because of due to the increase in coke concentration. When the fluidization reaches a stable state, the MTO reactions are at the steady state. In this case, the product concentrations in the outlet of FBR can be used to verify the coupled model of the reacting flow. In addition, the selectivity of ethylene and

12

Y.-Q. Zhuang et al. / Computers and Chemical Engineering 60 (2014) 1–16

2010), indicating that the coupled model of the reacting flow is capable of simulating the MTO process in FBR.

Fig. 14. Main reaction parameter distribution profiles in the MTO FBR at t = 2.6 s ((a) gas-phase temperature, (b) particle temperature, (c) coke content, (d) ethene mole concentration, (e) propene mole concentration, and (f) butene mole concentration; space velocity = 2.8 m/s, inlet feed temperature = 723 K, feed ratio of water to methanol = 0).

propylene, as well as the methanol conversion, is almost identical in any type of reactor (including microbalance reactor, fixed-bed reactor, and FBR) during the MTO process, suggesting that the bulk kinetic behavior in a microbalance reactor is almost identical to that in an FBR under the reaction conditions of high methanol concentrations (Hu, Cao, Ying, Sun, & Fang, 2010; Zhuang et al., 2012). Although deduced at the laboratory scale, the kinetic model can be applied for the macro-reactor because the influence of reactor type on the bulk reaction mechanisms is negligible. The experimental data collected from a practical FBR can be used to validate the simulated data. Table 6 shows that the simulation data are in good agreement with the experimental data (Qi et al.,

3.3.2. Effects of key operation parameters The SAPO-34 catalyst becomes deactivated because of coke formation. The coke covers the active sites of the catalyst, thereby reducing the free space in the cavities of the catalyst and ultimately influencing the MTO reaction. The catalyst was completely deactivated when the MTO reaction proceeds to 3600 s. Simulating the complete MTO reaction is time consuming. However, the catalyst was found to be deactivated rapidly at the beginning of the reaction, and the change in reaction parameters is evident at 26 s. Therefore, investigating the effects of the operating conditions on the MTO reaction from the beginning up to ∼26 s is sufficient. The coupled model can be used to evaluate the effects of feed temperature, ratio of water to methanol in the feed, and feed velocity on the reactor performance. The results are shown in Figs. 15–17. Fig. 15a and b display the predicted average bed and particle temperature profiles at different feed temperatures, respectively. Both the average bed and particle temperatures initially increase and then decrease at certain feed temperatures. Moreover, both parameters increase with the increase in feed temperature, because the MTO process is highly exothermic, that is, the exothermic heat increases when methanol is totally converted and decreases because of the deactivation of catalysts. The increase in feed temperature also results in the increase in bed temperature. In practice, the bed temperature can be reflected by the particle temperature. The average outlet mole ratio of ethane to propene and the coke content profiles at different feed temperatures are presented in Fig. 15c and d, respectively. As shown in Fig. 15c, the average ratio of ethane to propene increases with the increase in feed temperature. This finding is consistent with that reported by Wu, Michael, and Rayford (2004). In addition, the average outlet mole ratio of ethane to propene decreases with catalyst deactivation at certain feed temperatures. These results can be explained using the MTO reaction sub-model. When the catalyst is deactivated, the exothermic heat decreases, and then the bed temperature decreases as well, resulting in the decrease in ratio, as shown in Fig. 15c. All of these findings indicate that the MTO reaction rate indirectly affects all of these changes, which are related to catalyst deactivation, as shown in Fig. 15d. Moreover, the catalyst decays faster at higher temperatures, suggesting that a higher feed temperature results in a higher coke content. Fig. 16 shows the simulated effect of the feed ratio of water to methanol on the average bed parameters, which are in contrast to those described in Fig. 15. With the increase in feed ratio of water to methanol, the simulated average bed and particle temperatures, as well as the average outlet mole ratio of ethane to propene, decrease, whereas, coke content increases. For the bed and particle temperatures, the addition of water initially slows down the average bed temperature and then reduces the partial pressure of methanol and the heat released from the reaction. Therefore, the lower bed and particle temperatures resulting from the increase in the feed ratio of water to methanol must result in a decrease in the ratio of ethane to propene and an increase in coke content (Fig. 15c and d). Fig. 17 shows the simulated effect of space velocity on the average bed parameters, and the obtained results are similar to those described in Fig. 16. Fig. 17a to c prove that all of the

Table 6 Comparison between the practical data and the simulation data (the space velocity = 2.8 m/s, the inlet feed temperature = 723 K and the feed ration of water to methanol = 0.) (Qi et al., 2010).

Practical data Simulation data

Methanol conversion (%)

Ethene (%)

Propene (%)

Butene (%)

C5 (%)

99.13 100

42.58 40.2

38.63 38.3

10.96 13.2

3.47 2.15

Y.-Q. Zhuang et al. / Computers and Chemical Engineering 60 (2014) 1–16

13

Fig. 15. Effects of feed temperature on the reacting flow field parameters ((a) the average bed temperature, (b) the average particle temperature, (c) the ratio of ethene to propene, and (d) the coke content; space velocity = 2.8 m/s, feed ratio of water to methanol = 0).

Fig. 16. Effects of the feed ratio of water to methanol on the reacting flow field parameters ((a) the average bed temperature, (b) the average particle temperature, (c) the ratio of ethene to propene, and (d) the coke content; inlet feed temperature = 673 K, space velocity = 2.8 m/s).

14

Y.-Q. Zhuang et al. / Computers and Chemical Engineering 60 (2014) 1–16

Fig. 17. Effects of the space velocity on the reacting flow field parameters ((a) average bed temperature, (b) average particle temperature, (c) ratio of ethene to propene, (d) coke content; inlet feed temperature = 673 K, feed ratio of water to methanol = 0).

simulated average bed and particle temperatures, as well as the ratio of ethane to propene, decrease with the increase in space velocity. By contrast, coke content increases with the increase in space velocity, as described in Fig. 17d. These changes are also related to the reaction rate, that is, a lower space velocity (or a higher conversion with longer contact time) results in a faster reaction rate and the release of more heat, resulting in higher bed and particle temperatures. 4. Conclusions A coupled CFD–DEM model, which integrates the sub-models that describe the heat transfer behavior between particles and that between gas and particles as well as the lumped kinetics in the gas phase for the MTO process, was developed to characterize the gas–solid flow behavior in an MTO FBR. This model successfully captured the important flow features in an MTO FBR; such features include the bed pressure drop profiles at cold-flow conditions and the outlet product concentrations at reacting-flow conditions. In addition, the particle motion pattern and key flow-field parameter distributions in the reactor were analyzed. The effects of the key operation parameters on the reacting flow field were also analyzed, and the findings are summarized as follows (1) The evolution of three typical regions in the FBR and four typical fluidization states with different gas velocities was investigated using CFD–DEM model. In addition, the gas–solid flow pattern, particularly the solid flow pattern, which cannot be obtained using the Eulerian–Eulerian method, was also successfully captured by the CFD–DEM model. A typical annulus-core structure for particle motions was observed in the simulation results under cold-flow conditions. This type of particle flow pattern promotes excellent particle mixing and heat transfer efficiency, which is essential to the MTO reactor because the MTO process is exothermic.

(2) The simulation by the CFD–DEM model under reaction-flow conditions were conducted, and the distributions of temperature, coke content and product mass fractions in FBR at different regions were obtained. The results showed that for the MTO process, the parameter distributions at different regions have significant differences. The temperature distribution in FBR is determined by the exothermic heat of the MTO reaction, and the product fraction distributions in FBR are uniform because of the excellent mass transfer capability of FBR. (3) The feed temperature, the feed ratio of water to methanol, and the space velocity have significant effects on the reactor parameter distributions. All of the average bed and particle temperatures, as well as the average ratio of ethane to propene, increase with the increase in feed temperature. However, a higher feed temperature results in higher coke content. Moreover, the feed ratio of water to methanol and the space velocity have opposing effects on the reactor performance. Acknowledgments The authors thank the National Natural Science Foundation of China (No. 201076171), the National Ministry of Science and Technology of China (No. 2012CB21500402) and the State-Key Laboratory of Chemical Engineering of Tsinghua University (No. SKL-ChE-13A05) for supporting this work. References Aguayo, A. T., Mier, D., Gayubo, A. G., Gamero, M., & Bilbao, J. (2010). Kinetics of methanol transformation into hydrocarbons on a HZSM-5 zeolite catalyst at high temperature (400–550 ◦ C). Industrial and Engineering Chemistry Research, 49, 12371–12378. Alwahabi, S. M., & Froment, G. F. (2004). Conceptual reactor design for the methanolto-olefins process on SAPO-34. Industrial and Engineering Chemistry Research, 43, 5112–5122. Borkink, J. G. H., & Westerterp, K. R. (1992). Influence of tube and particle diameter on heat-transport in packed-beds. AIChE Journal, 38, 703–712.

Y.-Q. Zhuang et al. / Computers and Chemical Engineering 60 (2014) 1–16 Chai, C. J., & Zhang, G. L. (2004). Chemical engineering fluid flows and heat transport. Beijing: Chemical Engineering Press (Chinese). Chen, X. Z., Luo, Z. H., Yan, W. C., Lu, Y. H., & Ng, I. S. (2011). Three-dimensional CFD–PBM coupled model of the temperature fields in fluidized bed polymerization reactors. AIChE Journal, 57, 3351–3366. Chen, D., Gronvold, A., Moljord, K., Fuglerud, T., & Holmen, A. (1999). The effect of crystal size of SAPO-34 on the selectivity and deactivation of the MTO reaction. Microporous and Mesoporous Materials, 29, 191–203. Chen, D., Rebo, H. P., Moljord, K., & Holmen, A. (1999). Methanol conversion to light olefins over SAPO-34. Sorption, diffusion, and catalytic reactions. Industrial and Engineering Chemistry Research, 38, 4241–4249. Chen, D., Rebo, H. P., Gronvold, A., Moljord, K., & Holmen, A. (2000). Methanol conversion to light olefins over SAPO-34: Kinetic modeling of coke formation. Microporous and Mesoporous Materials, 35–36, 121–135. Chen, D., Gronvold, A., Moljord, K., & Holmen, A. (2007). Methanol conversion to light olefins over SAPO-34: Reaction network and deactivation kinetics. Industrial and Engineering Chemistry Research, 46, 4116–4123. Chu, K. W., & Yu, A. B. (2008). Numerical simulation of the gas–solid flow in threedimensional pneumatic conveying beeds. Industrial and Engineering Chemistry Research, 47, 7058–7071. Chu, K. W., Wang, B., Xu, D. L., Chen, Y. X., & Yu, A. B. (2011). CFD–DEM simulation of the gas–solid flow in a cyclone separator. Chemical Engineering Science, 66, 834–847. Deen, N. G., Annaland, M. V., van der Hoef, M. A., & Kuipers, J. A. M. (2007). Review of discrete particle modeling of fluidized beds. Chemical Engineering Science, 62, 28–41. Deliang, S., Watson, L. V., & McCarthy, J. J. (2008). Heat transfer in rotary kilns with interstitial gases. Chemical Engineering Science, 63, 4506–4516. Di Felice, R. (1994). The voidage function for function for fluid–particle interaction systems. International Journal of Multiphase Flow, 20, 153–162. Fatourehchi, N., Sohrabi, M., Royaee, S. J., & Mirarefin, S. M. (2011). Preparation of SAPO-34 catalyst and presentation of a kinetic model for methanol to olefin process (MTO). Chemical Engineering Research and Design, 89, 811–816. Feng, Y. Q., Xu, B. H., Zhang, S. J., Yu, A. B., & Zulli, P. (2004). Discrete particle simulation of gas fluidization of particle mixtures. AIChE Journal, 50, 1713–1728. Gao, J. S., Chang, J., Lu, C. X., & Xu, C. M. (2008). Experimental and computational studies on flow behavior of gas–solid fluidized bed with disparately sized binary particles. Particuology, 6, 59–71. Gao, J. S., Lan, X. Y., Fan, Y. P., Chang, J., Wang, G., Lu, C. X., & Xu, C. M. (2009). CFD modeling and validation of the turbulent fluidized bed of FCC particles. AIChE Journal, 55, 1680–1694. Gayubo, A. G., Aguayo, A. T., del Campo, A. E., Tarrio, A. M., & Bilbao, J. (2000). Kinetic modeling of methanol transformation into olefins on a SAPO-34 catalyst. Industrial and Engineering Chemistry Research, 39, 292–300. Gayubo, A. G., Benito, P. J., Aguayo, A. T., Castilla, M., & Bilbao, J. (1996). Kinetic model of the MTG process taking into account the catalyst deactivation: Reactor simulation. Chemical Engineering Science, 51, 3001–3006. Geng, Y. M., & Che, D. F. (2011). An extended DEM–CFD model for char combustion in a bubbling fluidized bed combustor of inert sand. Chemical Engineering Science, 66, 207–219. Hoomans, B. P. B., Kuipers, J. A. M., Briels, W. J., & van Swaaij, W. P. M. (1996). Discrete particle simulation of bubble and slug formation in a two-dimensional gas-fluidized bed: A hard sphere approach. Chemical Engineering Science, 51, 99–118. Hu, H., Cao, F. H., Ying, W. Y., Sun, Q. W., & Fang, D. Y. (2010). Study of coke behaviour of catalyst during methanol-to-olefins process based on a special TGA reactor. Chemical Engineering Journal, 160, 770–778. Jin, Y., Zhu, J. X., Wang, Z. W., & Yu, Z. Q. (2001). Fluidization engineering principles. Beijing: Tsinghua University Press. Kaneko, Y., Shiojima, T., & Horio, M. (1999). DEM simulation of fluidized beds for gas-phase olefin polymerization. Chemical Engineering Science, 54, 5809–5821. Karimi, S., Mansourpour, Z., Mostoufi, N., & Sotudeh-Gharebagh, R. (2011). CFD–DEM study of temperature and concentration distribution in a polyethylene fluidized bed reactor. Particulate Science and Technology, 29, 163–178. Lettieri, P., Felice, R. D., Pacciani, R., & Owoyemi, O. (2006). CFD modelling of liquid fluidized beds in slugging mode. Powder Technology, 167, 94–103. Levenspiel, O. (2002). Chemical reaction engineering. New York, U.S.A: John Wiley & Sons, Inc. Li, T. W., & Guenther, C. (2012). MFIX-DEM simulations of change of volumetric flow in fluidized beds due to chemical reactors. Powder Technology, 220, 70–78. Lim, E. W. C., Zhang, Y., & Wang, C. H. (2006). Effects of an electrostatic field in pneumatic conveying of granular materials through inclined and vertical pipes. Chemical Engineering Science, 61, 7889–7908. Lwahabi, A. S. M., & Forment, G. F. (2004). Single event kinetic modeling of the methanol-to-olefins process on SAPO-34. Industrial and Engineering Chemistry Research, 43, 5098–5111. Mahecha-Botero, A., Grace, J. R., Elnashaie, S. S. E. H., & Lim, C. J. (2006). Comprehensive modelling of gas fluidized-bed reactors allowing for transients, multiple flow regimes and selective removal of species. International Journal of Chemical Reactor Engineering, 4, A11. Mahecha-Botero, A., Grace, J. R., Elnashaie, S. S. E. H., & Lim, C. J. (2009). Advances in modeling of fluidized-bed catalytic reactors: A comprehensive review. Chemical Engineering Communications, 196, 1375. Mier, D., Gayubo, A. G., Aguayo, A. T., Olazar, M., & Bilbao, J. (2011). Olefin production by cofeeding methanol and n-Butane: Kinetic modeling considering the deactivation of HZSM-5 zeolite. AIChE Journal, 57, 2841–2853.

15

Ouyang, J., & Li, J. (1999a). Particle-motion-resolved discrete model for simulating gas–solid fluidization. Chemical Engineering Science, 54, 2077–2083. Ouyang, J., & Li, J. (1999b). Discrete simulations of heterogeneous structure and dynamic behavior in gas–solid fluidization. Chemical Engineering Science, 54, 5427–5440. Park, T. Y., & Froment, G. F. (2001a). Kinetic modeling of the methanol to olefins process. 1. Model formulation. Industrial and Engineering Chemistry Research, 40, 4172–4186. Park, T. Y., & Froment, G. F. (2001b). Kinetic modeling of the methanol to olefins process. 2. Experimental results, model discrimination and parameter estimation. Industrial and Engineering Chemistry Research, 40, 4187–4198. Qi, Y., Liu, Z. M., Lv, Z. H., Wang, H., He, C. Q., & Zhang J. L., Wang, X. G. (2010). Method for producing light olefins form methanol or/and dimethyl ether. US. 20100063336AI. Ranade, V. V. (2000). Computational flow modeling for chemical reactor engineering. London: Academic Press. Ranz, W. E. (1952). Friction and transfer coefficients for single particles and packed beds. Chemical Engineering and Processing, 48, 247–256. Reid, R. C., Prausnitz, J. M., & Poling, B. E. (1987). The properties of gases and liquids (4th ed., pp. 125–137). New York, US: McGraw-Hill Press Inc. Schoenfelder, H., Hinderer, J., Werther, J., & Keil, F. J. (1995). Methanol to olefins—prediction of the performance of a circulating fluidized-bed reactor on the basis of kinetic experiments in a fixed-bed reactor. Chemical Engineering Science, 49, 5377–5390. Shi, D. P., Luo, Z. H., & Guo, A. Y. (2010). Numerical simulation of the gas–solid flow in fluidized-bed polymerization reactors. Industrial and Engineering Chemistry Research, 49, 4070–4079. Soundararajan, S., Dalai, A. K., & Berruti, F. (2001). Modeling of methanol to olefins (MTO) process in a circulating fluidized bed reactor. Fuel, 80, 1187–1197. Shackelford, J., Alexande, W., & Pork, J. S. (2002). CRC materials science and engineering handbook. New York, US: CRC Press Inc. Shih, T. H., Liou, W. W., Shabbir, A., Yang, Z., & Zhu, J. (1995). A new k–ε eddy viscosity model for high Reynolds number turbulent flows model development and validation. Computers & Fluids, 24, 227–235. Tsuji, T., Yabumoto, K., & Tanaka, T. (2008). Spontaneous structures in threedimensional bubbling gas-fluidized bed by parallel DEM–CFD coupling simulation. Powder Technology, 184, 132–140. Tsuji, Y., Kawaguchi, T., & Tanaka, T. (1993). Discrete particles simulation of twodimensional fluidized bed. Powder Technology, 77, 79–87. Utikar, R. P., & Ranade, V. V. (2007). Single jet fluidized beds: Experiments and CFD simulations with glass and polypropylene particles. Chemical Engineering Science, 62, 167–183. van der Hoef, M. A., Annaland, M. V., Deen, N. G., & Kuipers, J. A. M. (2008). Numerical simulation of dense gas–solid fluidized beds: A multiscale modeling strategy. Annual Review of Fluid Mechanics, 40, 47–59. Vaishali, S., Roy, S., & Mills, P. L. (2008). Hydrodynamic simulation of gas–solids downflow reactors. Chemical Engineering Science, 63, 5107–5119. Vargas, W., & McCarthy, J. (2001). Heat conduction in granular materials. AIChE Journal, 47, 1052–1063. Wu, C. N., Cheng, Y., Ding, Y. L., & Jin, Y. (2010). CFD–DEM simulation of gas–solid reacting flows in fluid catalytic cracking (FCC) process. Chemical Engineering Science, 65, 542–549. Wu, C. N., Yan, B. H., Jin, Y., & Cheng, Y. (2010). Modeling and simulation of chemically reacting flows in gas–solid catalytic and non-cataytic processes. Particuology, 8, 525–530. Wu, X. C., Michael, G. A., & Rayford, G. A. (2004). Methanol conversion on SAPO34: Reaction condition for fixed-bed reactor. Applied Catalysis A: General, 260, 63–69. Xing, A. H., Yue, G., Zhu, W. P., Lin, Q., & Li, Y. (2010). Recent advances in typical technology for methanol conversion to olefins II. Development of chemical process. Modern Chemical Industry, 30, 18–25 (in Chinese). Xu, B. H., & Yu, A. B. (1997). Numerical simulation of the gas–solid flow in a fluidized bed by combining discrete particle method with computational fluid dynamics. Chemical Engineering Science, 62, 2785–2809. Xu, M., Ge, W., & Li, J. H. (2007). A discrete particle model for particle-fluid flow with considerations of sub-grid structures. Chemical Engineering Science, 62, 2302–2308. Yan, W. C., Luo, Z. H., Lu, Y. H., & Chen, X. D. (2012). A CFD–PBM–PMLM integrated model for the gas–solid flow fields in fluidized bed polymerization reactors. AIChE Journal, 58, 1717–1732. Zhao, Y. Z., Ding, Y. L., Wu, C. N., & Cheng, Y. (2010). Numerical simulation of hydrodynamics in downers using a CFD–DEM coupled approach. Powder Technology, 199, 2–12. Zhao, Y. Z., Cheng, Y., Wu, C. N., Ding, Y. L., & Jin, Y. (2010). Eulerian–Lagrangian simulation of distinct clustering phenomena and RTDs in riser and downer. Particuology, 8, 44–50. Zhang, M. H., Chu, K. W., Wei, F., & Yu, A. B. (2008). A CFD–DEM study of the cluster behavior in riser and downer reactors. Powder Technology, 184, 151–165. Zhang, Y., Lim, E. W. C., & Wang, C. H. (2007). Pneumatic transport of granular materials in an inclined conveying pipe: Comparison of computational fluid dynamics-discrete element method (CFD–DEM), electrical capacitance tomography (ECT), and particle image velocimetry (PIV) results. Industrial and Engineering Chemistry Research, 46, 6066–6083. Zhuang, Y. Q., Gao, X., Zhu, Y. P., & Luo, Z. H. (2012). CFD modeling of methanol to olefins process in a fixed-bed reactor. Powder Technology, 221, 419–430.

16

Y.-Q. Zhuang et al. / Computers and Chemical Engineering 60 (2014) 1–16

Zhuang, Y. Q. (2012). CFD modeling of methanol to olefins process in the fixed-bed and fluidized-bed reactor. Xiamen, PR China: Xiamen University. Master Dissertation (In Chinese). Zhu, H. P., Zhou, Z. Y., Yang, R. Y., & Yu, A. B. (2008). Discrete particle simulation of particulate systems: A review of major applications and findings. Chemical Engineering Science, 63, 5728–5770. Zhu, J., Cui, Y., Chen, Y. J., Zhou, H. Q., Wang, Y., & Wei, F. (2010). Recent researches on process from methanol to olefins. CIESC Journal, 61, 1674–1684 (in Chinese).

Zhou, Z. Y., Pinson, D., Zou, R. P., & Yu, A. B. (2011). Discrete particle simulation of gas fluidization of ellipsoidal particle. Chemical Engineering Science, 66, 6128–6145. Zhu, H. P., Zhou, Z. Y., Yang, R. Y., & Yu, A. B. (2007). Discrete particle simulation of particulate systems: Theoretical developments. Chemical Engineering Science, 62, 3378–3396.