Modeling of the maleic anhydride circulating fluidized bed reactor

Modeling of the maleic anhydride circulating fluidized bed reactor

Accepted Manuscript Title: Modeling of the Maleic Anhydride Circulating Fluidized Bed Reactor Authors: Pranava Chaudhari, Sanjeev Garg PII: DOI: Refer...

2MB Sizes 3 Downloads 170 Views

Accepted Manuscript Title: Modeling of the Maleic Anhydride Circulating Fluidized Bed Reactor Authors: Pranava Chaudhari, Sanjeev Garg PII: DOI: Reference:

S0098-1354(17)30068-6 http://dx.doi.org/doi:10.1016/j.compchemeng.2017.02.012 CACE 5715

To appear in:

Computers and Chemical Engineering

Received date: Revised date: Accepted date:

20-5-2016 1-11-2016 6-2-2017

Please cite this article as: Chaudhari, Pranava., & Garg, Sanjeev., Modeling of the Maleic Anhydride Circulating Fluidized Bed Reactor.Computers and Chemical Engineering http://dx.doi.org/10.1016/j.compchemeng.2017.02.012 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Modeling of the Maleic Anhydride Circulating Fluidized Bed Reactor Pranava Chaudhari and Sanjeev Garg* Department of Chemical Engineering, Indian Institute of Technology, Kanpur 208016, India [email protected], [email protected]

* to whom the correspondence should be addressed

Highlights:     

A simple model is proposed for the maleic anhydride circulating fluidized bed. A dispersion model with variable gas density is used in the hydrodynamics model. Configurational complexities of the reactors are considered in the proposed model. The tuned reactor models are obtained by the data fitting exercise. A multi-objective model tuning is proposed for the multi component reactor system.

Abstract: Maleic Anhydride (MA) production by selective oxidation of n-butane in a Multi-Tubular Fixed Bed reactor is constrained by the flammability limits. The use of Fluidized Bed circumvents this constraint but suffers from high back mixing. The Circulating Fluidized Bed (CFB) reduces the back mixing of the

1

catalyst and enhances the catalyst activity by reactivating the catalyst. A simple reactor model with suitable hydrodynamic and kinetic models at the operating conditions of both the pilot and commercial MA CFB reactors is proposed. A dispersion model with variable gas density is used for modeling the hydrodynamics of these reactors. The reactors' configurational complexities are also accounted for in the model. The data fitting exercise of the proposed reactor model is performed. The proposed method simultaneously minimizes the total bias based on the individual biases for each component of the multi component reactor system and the total error in a multi-objective framework.

KEYWORDS: Circulating Fluidized Bed; Maleic Anhydride; Reactor Modeling; Variable Density Reactor; Data Fitting Exercise; Multi-Objective Optimization.

Nomenclature

 

ACS , i

cross-sectional area of the reactor section i m2

CT

total concentration in the gaseous phase kmol m3

DT , i

diameter of the reactor section i  m 

dp

mean particle diameter  m 

Fj

molar flow of the jth-component  kmol s 

f conv

conversion factor  kg ms 2bar 

HT

total height of the reactor

HT ,i

height of the reactor section i  m 

NL

total oxygen capacity of the lattice site in the catalyst bulk  mmolO2 kg of cat 

PT

total pressure in the gaseous phase  bar 

pj

partial pressure of the jth-component





 bar 

2

Rj

reaction rate of the jth-component  kmol kg of cat  s 

Rg'

gas constant  kJ kmol  K 

Rg

gas constant  m3bar K  kmol 

rk '

rate of the k ' th catalytic reaction  kmol kg of cat  s 

T

temperature  K 

v pz , t

terminal settling velocity of the single particle  m s 

Yj

mole fraction of the gaseous jth-component   

Greek symbols Greek letters

g

viscosity of the gaseous mixture  kg ms 

g

density of the gaseous mixture kg m3

cat

density of the catalyst kg m3



void fraction of the reactor bed  m3 of void m3 of reactor 

 OL

fraction of occupied lattice oxygen sites in the catalyst bulk   









Superscript/subscript

i

section of the reactor (1: transport bed; 2: riser)

k'

specific catalytic reaction step

k"

adsorption desorption step (1: O2; 2: Bu; 3: MA)

1. Introduction The various valuable chemicals and intermediates are produced in the process industries by selective oxidation of the small alkanes in the presence of specific catalysts (Schlogl, 2009). The Maleic Anhydride (MA) production by the selective oxidation of n-butane (Bu) in the presence of Vanadium3

Phosphorus-Oxide (VPO) catalyst is one of the most successful example (Ballarini et al., 2006) (Global market of $2.8 billion with projected 7.2% increase for 2014-2019; Source: report by researchandmarket ID: 3292451, June 2015). The Multi-Tubular Fixed Bed Reactor (MTFBR) is predominantly used in the industrial processes (Contractor, 1999). The limited heat removal in a MTFBR generally leads to hot spots within the reactor in case of exothermic reactions. A previous study (Chaudhari and Gupta, 2012) has shown that the optimal production of the MA in a MTFBR is controlled by constraints of the flammability limit and the auto ignition temperature (AIT). The fluidized bed reactors have much better heat transfer characteristics and uniform temperature control that also reduces the chances of hot spots in the fluidized reactor. Moreover, in the fluidized bed reactors the catalyst particles also act as flame arrestors by quenching the free radicals of the homogeneous combustion reactions (Contractor et al., 1987). The fluidized bed reactors can, thus, be operated with much higher Bu (~ 4%) in the feed as compared to the MTFBR (~ 1.8%) (Contractor, 1999). The lower concentration of Bu in the MTFBR is controlled by its flammability limit. In the open literature, three major problems associated with the scale-up of a fluidized bed are reported (Contractor and Sleight, 1987). Firstly, the “gulf streaming” that creates a wide gas residence time distribution that may lead to reduced conversions; secondly, the “significant solid and gas back-mixing” that may lead to the loss in the selectivity of the desired product and thirdly the “limited oxygen availability” because of high Bu fraction in the feed that may lead to the over reduction of the catalyst. The use of a MA Circulating Fluidized Bed (CFB) reactor has been proposed (Contractor and Sleight, 1987) by DuPont that not only retains all the advantages of a fluidized bed reactor over a MTFBR but also overcomes the three problems associated with the fluidized bed reactors mentioned above. In a CFB reactor, the extra oxygen is fed separately through the spargers that results in less over reduction of the catalyst. The separate oxygen feed in the MA CFB reactor allows the reactor to operate with higher Bu concentration (25%) (Contractor, 1999) as compared to the fixed bed and MTFBRs. The catalyst particles in a CFB reactor are cycled between the fuel enriched and the oxygen enriched

4

environments. The lattice oxygen present on the VPO catalyst contributes to the reaction in the fuel enriched environment and is previously reported as the Lattice Oxygen Contribution (LOC) of the catalyst (Patience and Bockrath, 2010). The Langmuir-Hinshelwood reaction kinetics reported in the reactor model of a MTFBR (Chaudhari and Gupta, 2012) is, thus, not applicable for the operating conditions of a MA CFB reactor. Several other kinetic models are reported in the open literature for the CFB reactors (Gascon et al., 2006; Lorences et al., 2004, 2003; Patience and Lorences, 2006; Shekari and Patience, 2013; Wang and Barteau, 2003, 2002, 2001; Xiao-Feng et al., 2002). Ballarini et al., (2006), in their review, conclude that despite being one of the most studied reactions in the selective oxidation of the small alkanes the final picture for MA production over VPO catalyst is still not clear. Several studies (Avidan and Yerushalmi, 1982; Bi and Grace, 1995; Chan et al., 2010; Grace et al., 1999; Gupta, 1998; Yerushalmi and Cankurt, 1979; Zhang et al., 2015) report the different operating hydrodynamic regimes of the different gas-solid systems. The different hydrodynamic regimes present in the CFB reactors with solid recirculation are, namely, Dense Suspension Upflow (DSU), Fast Fluidized Bed (FFB), Core Annular Flow Bed (CAFB) and Pneumatic Conveying (PC). Recently, Zhu and Zhu (2008a) have suggested another hydrodynamic regime present in the CFB reactors as Circulating-Turbulent Fluidized Bed (C-TFB). The superficial gas velocity, U 0 , and the solid mass flux, G pz , are the two main parameters for defining these different hydrodynamic regimes. A schematic

representation of these regimes and the criteria are shown in Fig. 1. These criteria are the critical values of U 0 ( U ff , U CA , U C TFB and U PC ) and G pz ( G pzSCC ; SCC: saturation carrying capacity) across which the transition from one regime to another regime takes place. The correlations for criteria of regime transition along with the corresponding references are outlined in Table 1. The correlations given in Table 1 use U 0 measured at the outlet i.e., U 0OUT , as reported in the open literature (Patience and Bockrath, 2009).

Pugsley et al., (1992) have proposed a MA CFB reactor model without the data fitting exercise with the pilot or commercial MA CFB reactor data. The hydrodynamic model reported with this reactor

5

model requires experimental pressure drop profiles. Pugsley and Berruti, (1996) have suggested a predictive core annular hydrodynamic model that did not require the experimental pressure drop profiles. The reported core annular hydrodynamic model considers a downward solid movement in the annular region near the reactor wall. This downward solid movement disappears in the CFB reactors operating in the DSU regime with high G pz values (Grace et al., 1999). Thus, the reported core annular hydrodynamic model is suitable for the combustion reactors that operate at lower G pz values and is not suitable for the catalytic reactors that operate at much higher G pz values (Berruti et al., 1995). The absence of a downward solid movement in the annular region of the catalytic reactors limits the applicability of many other core annular hydrodynamic models reported in the open literature (Bai et al., 1995; Gupta, 1998; Puchyr et al., 1997) for modelling the hydrodynamics of a MA CFB reactor. A generic fluidized-bed reactor hydrodynamic model is also reported using a dispersion model with dense and lean sections in terms of the solid hold up (Abba et al., 2003). The model is generic as it is applicable to bubbling fluidized bed (BFB), turbulent fluidized bed (TFB) and fast fluidized bed (FF) regimes. The limitation of this generic model is that most of the catalytic reactors operate in the DSU regime and this generic model does not consider the DSU regime. Thus, a reactor model is desirable with suitable hydrodynamic and kinetic models valid at the operating conditions present in a MA CFB reactor. Moreover, the desired reactor model should be simple and accurate and the same is proposed in this work. The motivation of the proposed work is described in Section 1. A brief description of a Pilot Plant (PP) and a Commercial Plant (CP) MA CFB reactor is given in Section 2. The development of the proposed reactor model is discussed in Section 3. A detailed data fitting exercise using the PP and CP reactors is discussed in Section 4. The results obtained from the data fitting exercise are discussed in Section 5. The conclusions of this work are given in Section 6. 2. Brief Description of the MA CFB reactors

6

The MA CFB PP reactor (Patience and Bockrath, 2010) comprises of a transport bed, a riser, a riser stripper, a regenerator and a regenerator stripper. The MA CFB CP reactor (Patience and Bockrath, 2010) has no regenerator stripper but instead a cyclone separator is present. The details of these two reactors are not given here for brevity and may be found elsewhere (Hutchenson et al., 2010; Patience and Bockrath, 2010). Here, the schematic of a MA CFB reactor is shown in Fig. 2. Fresh Bu and nitrogen (N2) along with the recycle gas, which contains Bu, oxygen (O2), carbon dioxide (CO2) and carbon mono-oxide (CO) (Patience and Bockrath, 2010) are fed to the bottom of the transport bed. The velocity of the gas mixture is sufficient enough to fluidize and lift the catalyst particles coming from the regenerator in the upward direction. Some fresh O2 is fed separately through the spargers to avoid the flammability limit at the inlet gas mixture. It is noted that the fresh O2 is fed at three different locations at three different heights in the transport bed to avoid the over reduction of the catalyst. The oxidation reactions take place on the catalyst particles as the gas and solid move through the transport bed and the riser section of the reactor. The catalyst is reduced by the consumption of its lattice oxygen (Patience and Bockrath, 2010). The reduced catalyst is separated from the gas mixture at the outlet of the riser and send to the regenerator. In the regenerator, the reduced catalyst is fluidized with the fresh oxygen to replenish the consumed lattice oxygen (Patience and Bockrath, 2010). In the MTFBR and fluidized bed reactors, the oxidation reactions over the catalyst particles and the oxidation of the reduced catalyst particles by the gaseous oxygen occur simultaneously, while in the MA CFB reactor these two steps occur in two different zones. The desired product, MA, is separated from the outlet gas mixture in the form of maleic acid and a significant portion of the remaining gas is recycled back to the reactor (Patience and Bockrath, 2010). Thus, as discussed above, the MA CFB reactor is a multi-component reactor system (MCRS) involving 8 components. 7 ( NGC ) of these are present in the gas phase (Bu, MA, O2, H2O, CO2, CO and N2) and 1 in the solid phase (the lattice oxygen as the fraction of occupied lattice oxygen sites in the catalyst bulk,  OL ). The regenerator is assumed to be ideal and at the steady state with complete

7

regeneration of the lattice oxygen in each cycle. In the proposed work, the modeling is done only for the transport bed and the riser of the MA CFB reactor. The riser stripper, cyclone separator, regenerator and regenerator stripper are not modeled. Thus, the inlet of the transport bed is considered as the inlet and the outlet of the riser is considered as the outlet in the proposed MA CFB reactor model. The various reactor and catalyst parameter values used in the modeling are given in Table 2. 3. Reactor Modeling 3.1. Kinetic Model The heterogeneous catalytic reaction mechanisms may broadly be divided into two categories based on the role of the catalyst. In the first case, the catalytic surface acts only as an adsorption platform (e.g., Langmuir-Hinshelwood, Eley-Rideal, etc.) and in the second case, the catalyst plays an active role in the reaction through the redox reaction (Mars-Van Krevelen). In the kinetic model development for these cases quasi steady state assumption (QSSA) for the catalyst sites is generally assumed. Several reported kinetic models applying QSSA on all the catalyst sites at the steady state conditions, applicable for the MA MTFBR, are summarized in the open literature (Xiao-Feng et al., 2002). These kinetic models are not applicable for the MA CFB reactors because even for the steady state operation of these reactors the catalyst particles are circulated between a fuel enriched environment and an oxygen enriched environment. Thus, the catalyst is never truly in a steady state condition. Several reaction mechanisms and transient kinetic models are reported (Gascon et al., 2006; Lorences et al., 2004, 2003; Patience and Lorences, 2006; Shekari and Patience, 2013; Wang and Barteau, 2002, 2001; Xiao-Feng et al., 2002) by simulating the experimental conditions similar to that of the MA CFB reactor. This is achieved by operating the experimental reactors in cycles with the fuel enriched and the oxygen enriched feeds. The transient state of the catalyst particles under such operating conditions is modeled in terms of variable catalyst sites. These variable catalyst sites are reported as the fraction of adsorption sites occupied by O2  O2  and the fraction of occupied lattice oxygen sites in the

8

catalyst bulk OL  by Gascon et al., (2006); as O2 ,  OL , the fraction of lattice surface sites occupied by





S oxygen and MA OS and  MA  by Xiao-Feng et al. (2002); the Vanadium sites V 5 ,V 4 and V C 4 by

Lorences et al. (2004, 2003 and Patience and Lorences (2006) and as  O by Shekari and Patience (2013) L

and Wang and Barteau (2002, 2001). It is noted that other catalyst sites are also involved in the reported reaction mechanisms but the QSSA is made for these catalyst sites in their respective kinetic models. The LOC information for both the MA CFB reactors (Patience and Bockrath, 2010) is given as the variable catalyst site,  OL . Though the kinetic models proposed by Wang and Barteau (2002, 2001) are based on  OL but the models are not applicable for the MA CFB reactors. In their model, either too much oxygen (21%) in the experimental feed (Wang and Barteau, 2002) or no oxygen in the experimental feed (Wang and Barteau, 2001) was reported whereas the CFB reactor is known to have 3.6 -12.9 % O2 (Patience and Bockrath, 2010) in the feed. Shekari and Patience (2013) reported a kinetic model for the MA CFB reactors operating under similar conditions and in terms of the  OL . In their kinetic model the selectivity of the MA is a ratio of the kinetic parameters which is inconsistent with the observed experimental data for the MA CFB reactors (Patience and Bockrath, 2010). Among the other transient kinetic models, the kinetic model proposed by Gascon et al. (2006) involves the least number of variable catalyst sites ( O2 and  OL ). The reaction mechanism (Gascon et al., 2006) involves other catalyst sites (

Bu , MA and  OS ) and QSSA is made for these sites. A modification of the Gascon’s kinetic model (Gascon et al., 2006) is proposed for the MA CFB reactors in this study. The proposed modified kinetic model uses a single variable catalyst site. Along with the reported QSSA on Bu , MA and  OS sites, QSSA on O2 is assumed. A fundamental approach is used to develop the modified kinetic model equations using the above assumption on O2 site and the reaction mechanism proposed by Gascon et al. (2006). The proposed kinetic model as given in Table 3 has the following three differences with the Gascon’s kinetic model

9

I. Instead of considering O2 as an independent variable, it is calculated from Eq. (1).

O  2

Keq ,1 pO2

(1)

1  Keq ,1 pO2  K eq ,2 pBu  K eq ,3 pMA

II. The extra term for the MA adsorption, Keq ,3 pMA , is present in the denominator of the expressions for r4 to r8 as given in Table 3 and represents the adsorption-desorption equilibrium of MA. III. The rate expression for r8 in Table 3 is proportional to the adsorbed MA and  OS as reported in the reaction mechanism but not present in the Gascon’s kinetic model (Gascon et al., 2006).





The rate constants,  kk '  and the equilibrium constants, Keq , k " as defined in the reactions

 k '  1, 2,..., Nrxn ; k "  1: O2 ; 2 : Bu;3: MA

given in Table 3 are expressed in the form of Arrhenius

equation and Van't Hoff equation given by Eq. (2) and Eq. (3), respectively

 E A, k ' kk ' T   kk0' exp  '   673Rg

  673    1   ; T   

 k '  1, 2,..., N rxn 

 H eq , k "   673   Keq , k " T   K eq0 , k " exp  1   ;  k "  1: O2 ; 2 : Bu;3: MA '  T    673Rg  

(2)

(3)

The reported values of the parameters used in these equations are given in Table 4. The two extra parameters, K eq0 ,3 and H eq ,3 , corresponding to the MA adsorption are introduced in the proposed kinetic model. The reaction rates, R j , for each of the NGC gaseous components and  OL are written in terms of the intrinsic reactions rates with the stoichiometric coefficients,  k ' j in Eq. (4). The stoichiometric coefficients are given in Table 5. N Rxn

R j   k ' j rk ' ; j  1, 2,..., NGC and OL

(4)

k '1

10

3.2. Hydrodynamic modeling The identification of the operating regimes for the transport bed and the riser of the MA CFB PP and CP reactors is necessary for the development of a suitable hydrodynamic model. The experimental values of the variables U 0,OUT i , G pz , i and  i , (i = 1: transport bed; i = 2: riser) are required to determine OUT the operating regimes of the given section i. The experimental values of U 0,2 and G pz ,2 for the riser

are calculated from the reported data (Patience and Bockrath, 2010) for both the MA CFB PP and CP reactors. The analysis based on these values, along with the correlations given in Table 1, conclude that the risers of both the MA CFB PP and CP reactors operate in the DSU regime. The experimental value OUT of U 0,1 for the transport bed cannot be calculated directly from the reported data (Patience and



OUT OUT Bockrath, 2010). An approximate value for U 0,1 is calculated by multiplying U 0,2 with ACS ,2 ACS ,1



assuming negligible variation in U 0 in the riser. The analysis based on these values, along with the correlations given in Table 1, conclude that the transport beds of both the MA CFB PP and CP reactors operate in the C-TFB regime. Another hydrodynamic study is reported (Patience and Bockrath, 2009) on the PP reactor only (not on the CP reactor) with the complete information for U 0,OUT i , G pz , i and  i , (i = 1: transport bed; i = 2: riser). The analysis based on this experimental data also results in the same conclusion for the MA CFB PP reactor, i.e., the transport bed and the riser of PP reactor operates in the C-TFB and DSU regimes, respectively. The hydrodynamic modeling of the gas and solid flow behaviors in a CFB reactor can be done using any of the following three approaches - Approach 1: Assuming ideal flow reactors, either as PFR or CSTR; Approach 2: Using the first principle models with Computation Fluid Dynamics (CFD); and Approach 3: Empirical modeling with non-ideal flows. Approach 1 is the simplest of the three approaches but lacks good prediction capabilities. Approach 2 is good in visualizing the flow behavior. Nonetheless, the prediction capabilities of the resulting reactor model may not be good until most of the fluid flow phenomena are appropriately modeled. Incorporating this constraint makes the model

11

computationally intensive. Approach 3 is further classified in two sub-approaches - Approach 3a: Modeling the flow complexities by approximating them to the assumed flow behaviors (e.g., the coreannular model or the cluster flow model) based on the experimental observations and Approach 3b: PFR with dispersion or tank in series models (i.e., “One-parameter models”). The single parameter in the PFR with dispersion model is the dispersion coefficient while it is the number of tanks for the tanks in series model. The experimental studies on the residence time distributions of the pulse inputs in both the gas (Liu, 2001) and solid phases (Andreux et al., 2008) of the CFB systems have shown significant breakthrough time which represents a behavior close to the PFR than the CSTR. The tail behavior in the output response, E(t), also shows deviation from the ideal plug flow behavior. In the reviews reported in the open literature by Bi (2004) and Breault (2006), the gas and solid hydrodynamics for CFB reactors are modeled as a PFR with dispersion. The same approach is used in this study. The radial and axial dispersions are modeled with only the effective dispersion coefficients (Schugerl, 1967) for both the gas

 D  and the solid  D  phases. A uniform radial distribution is assumed for both the gas velocity and ge

pe

the solid hold up as suggested by Liu (2001). 3.2.1. Gas phase modeling: Choong et al. (1998) have suggested that the dispersive molar flow rate of a gaseous species, FjDisp  x  , in a binary mixture is proportional to its mole fraction gradient analogous to the Fick’s law of

diffusion. The justification for the suggestion was that the sum of dispersion flow rates for both the components must be zero. This justification was made on the fact that the density variation or the total concentration variation across the length of flow is counter balanced with the variation in the gas velocity by the continuity equation in an open system. The molar flow rate of a gaseous species, Fj  x  , at a given dimensionless height, x  x  z HT  , with the dispersion in the MCRS (similar to the binary case of Choong et al., 1998) is given by Eq. (5). The total molar flow rate is only convective and is given by Eq. (6).

12

Fj  x   FjConv  x   FjDisp  x   ACSU 0  x  CT  x  Y j  x  

ACS Dge  x  CT  x  dY j  x  HT

dx

FT  x   ACSU 0  x  CT  x 

(5) (6)

The continuity equation for a given species in the MCRS with dispersion, variable gas density and reactions is given by Eq. (7). This equation in the limiting case of constant gas density simplifies to the expression given in Eq. (8) and is similar to the expression reported by Danckwerts (1953).

dY j  x   d CT  x  Y j  x U 0  x   1 d   cat 1    x  HT R j  x   0  Dge  x  CT  x   HT dx  dx  dx



Dge d 2C j  x  dx 2

HT

U0

dC j  x  dx





 cat 1    x  HT R j  x   0



(7)

(8)

The overall continuity equation (Eq. (9)) used in the proposed model considers the effect of the volumetric changes due to both the density variation (as a result of pressure drop) and the reactions. Note that Eq. (9) in the model is an overall material balance and hence only NGC  1 gaseous species (Bu, MA, O2, H2O, CO2 and CO) balances are independent and considered in the individual species balances given by Eq. (7). The mass transfer resistance over the catalyst particles and the diffusional resistance within the catalyst particles are neglected. The effective dispersion coefficient because of the non-ideal flow is assumed to be the same for each species and constant throughout the transport bed

 D  and the riser  D  . ge ,1

CT  x 

ge ,2





dU 0  x  dC  x  cat 1    x  HT  U0  x T  RT  x   0 dx dx CT  x 

(9)

Abba et al. (2002) in their variable-gas-density fluidized bed reactor model suggested the superficial gas velocity be defined as given in Eq. (10). Further, the species continuity equation and the reactions were proposed to be written in terms of Fj  x  for the MCRS. The expression for Fj  x  reported in their model considered only the convective molar flow and neglected the dispersive molar flow. A correct expression for Fj  x  , as given in Eq. (5), considers both the convective and dispersive molar

13

flows. The use of this expression, though, makes it quite difficult to express the reaction rates in terms of Fj  x  . The use of Eq. (7) and Eq. (9) in terms of Y j , CT and U 0 is proposed as a much easier way of developing the reactor model for the variable gas density systems.

U 0  x   U 0IN

FT  x  T  x  PTIN FTIN TTIN PT  x 

(10)

The pressure drops in the two sections (i = 1: transport bed; i = 2: riser) of the CFB reactor are calculated using Eq. (11). The first two terms represent the contributions of the gas and solid weights, respectively. The third term represents the contribution of the acceleration effects and is present only in the transport bed. This term is absent in the riser as the solid flow is developed in the transport bed itself. The last two terms represent the contributions of the gas-wall and the solid-wall frictions, respectively. The friction factors for the gas, f g and the solid, f s can be calculated from the correlations given in Bai et al. (1997). The effect of the various terms in Eq. (11) are analyzed on the reported hydrodynamic data (Patience and Bockrath, 2009) for both the sections of the MA CFB PP reactor, which are operating in the C-TFB and the DSU regimes, respectively. This analysis concludes that solid weight is the significant contributor (~98%) for the pressure drop. Eq. (12) is, thus, used in the proposed model for calculating the pressure drop. The void fraction  1 ( z )  for the C-TFB regime is reported to be axially uniform (Zhu and Zhu, 2008b) and, thus, a constant value, 1 ( z )   1 , is considered. The DSU operating regime for the riser is characterized with a significant axial variation only at the bottom section because of the acceleration effects, followed by a near uniform void fraction in the upper section of the riser work (Issangya, 1998). In MA CFB reactors, the bottom section of the riser operates in the C-TFB regime and hence the void fraction of the riser is characterized only by the upper section of the DSU regime and is considered constant  2 ( z )   2 . This is also in agreement with a recently reported experimental work (Wang et al., 2014). It is noted that this unique set up of the C-TFB operating regime at the bottom and the DSU regime at the top ensures the assumption of two different but constant void fractions (constant solid hold up) in the transport bed and the riser. It is noted that if a reactor is 14

operating in the DSU regime only, then the axial variation in the void fraction will be required to be considered in the hydrodynamic model. The compressibility factor at the reported experimental conditions (Grace et al., 1999; Issangya et al., 1996) is almost one. Thus, the gas phase is considered an ideal gas mixture and the total concentration is given by Eq. (13). 2   G pz 2 f g  i  gU 0,2 i HT ,1  HT    g g  i H T  Cat g 1   i  x  H T  H T ,1 Cat 1   1 DT , i   dPT   f conv   2 dx G pz   2 f s HT ,i    DT , i Cat 1   i  

(11)

dPT  x    f conv cat 1   i gHT i dx

(12)













CT  x  





PT  x  RgT  x 

(13)

The MA CFB reactors have excellent heat removal characteristics. The temperature rise in the transport bed and the riser of a MA CFB PP reactor are reported as 7 K and 4 K, respectively (Patience and Bockrath, 2010) irrespective of the very high exothermic reactions involved in the selective oxidation of butane. The equations in form of the constant temperature gradients, as given in Table 6, for both sections of the MA CFB PP reactor are used in the proposed model with the reported values for the temperature rise in each section. In the MA CFB CP reactor the temperature rise values are considered to be zero since the CP reactor operation is reported at isothermal conditions (Patience and Bockrath, 2010). 3.2.2 Solid phase dispersion model The deactivation of the catalyst is modeled in terms of the available lattice oxygen. The concentration of the available lattice oxygen at the dimensionless height, x, is given by Eq. (14). The lattice oxygen capacity of the catalyst (NL) is reported (Wang and Barteau, 2001) as 310 mmol O2/kg

15

and used in the proposed work. The molar flow for the lattice oxygen with both convective and dispersive molar flow rates is given by Eq. (15) and is similar to Eq. (5).





C L  x   N L cat 1   OL  x  O

(14)





 ACS N L cat 1   Dpe d L  O F L  x   FConv  FDisp  ACS G pz N LOL  x     x   L L O O O HT dx    

(15)

The constant solid phase effective dispersion coefficients in the transport bed and the riser

D

pe ,1

and Dpe,2  are considered similar to the gas phase effective dispersion coefficients. The molar

balance for the lattice oxygen across a differential segment of the reactor between heights x and x  dx is given by Eq. (16). The balance for  OL is given by Eq. (17) which is obtained by combining Eq. (15) and Eq. (14).

dF L  x  O

dx





 ACS cat 1   HT ROL  x 

HT G pz N L

(16)

 d L d  d L O  x   cat 1   Dpe N L  O  x   cat 1   HT2 ROL  x   0 dx dx  dx 









(17)

3.2.3 Summary of the basic model for the MA CFB reactor The complete set of model equations in terms of 2 NGC  3 simultaneous 1st order ODEs are summarized in Table 6. The first three ODEs are for the temperature, the total concentration and the continuity equation, respectively. The NGC 2nd order ODEs represent the dispersion model for NGC  1 gaseous species and  OL as given by Eq. (7) and Eq. (17), respectively. These NGC 2nd order ODEs are converted into a set of 2 NGC 1st order ODEs. 2 NGC  3 variables ( T , CT , U 0 ,  OL , dOL dx , Y j , and dY j dx ; j = 1, 2 ..., NGC -1) involved in the 2 NGC  3 1st order ODEs are termed as the model

variables. Thus, a total of 2 NGC  3 specified values are required to solve these model equations. 3.2.4 Hydrodynamic parameters of the model

16





The hydrodynamic parameters involved in the reactor model are the two void fractions  1 ,  2 and





the four effective dispersion coefficients Dge,1 , Dge,2 , Dpe,1 and Dpe,2 . A correlation for  i (i = 1: transport bed; i = 2: riser) of the MA CFB PP reactor is reported in the open literature (Patience and Bockrath, 2009) as given by Eq. (18). The correlation involves 4 parameters ( ai , bi , ci and d i ) in both the sections. The reported d 2 value is used with a sign change from positive to negative. This is required to fit the trend in the reported graph between C0,2 and U 0 (Wang and Barteau, 2001). The same sign change is also reported in the literature (Patience and Bockrath, 2009). In the absence of the “fines” data for the MA CFB PP and CP reactors (Patience and Bockrath, 2010), b1 is used as 0 in the proposed model. It is noted that for the riser section (i = 2) b2 is already reported as 0 (Keraron, 2009) and the same is used in the proposed model. The values of the parameters ( ai , bi , ci and d i ) used for the MA CFB PP reactor are given in Table 7. The same correlation given by Eq. (18) for   i  is used for the MA CFB CP reactor but the parameters are retuned. It is noted that the correlation given by Eq. (18) use U 0 measured at the outlet of the transport bed and the riser i.e., U 0,OUT i , as reported in the literature (Patience and Bockrath, 2009).

 i   mf  1   mf

U 0, i  U mf  G pz C0, i  U 0, i  U mf   cat 1   mf 

  V   g ,i

(18a)

where,

G pz

di C0, i  1  cU i 0, i

cat 1   mf 

Vg , i   ai d p ,50   Cat   g  2

g 1 18 g  fines bi

(18b)

(18c)

The four effective dispersion coefficients which represent the measure of gas and solid mixing are the other set of hydrodynamic parameters. Breault (2006), based on the reported correlations in the

17

literature, concluded that the different correlations applicable for different operating regimes result in different order of values for the gas and solid dispersion coefficients. Bi (2004), based on the experimental data from the open literature, reported that the Peclet numbers based on the effective diffusivity decrease with an increase in G pz both for the gas  Pege  and solid phases  Pe pe  . Moreover, after the onset of the DSU regime Pege and Pe pe become constant (Bi, 2004). The constant values are attributed to the disappearance of the downward flow in the annular region in the DSU regime. Wei et al., (1993) reported that the dispersion coefficient represented in terms of the Peclet numbers is proportional to the square root of the diameter of the reactor. Hence, the constant values of the Peclet numbers ( Pege and Pe pe ) could be different for the PP and CP reactors. The range of diameters of the experimental reactors used in the study by Bi (2004) is 0.03 m – 0.94 m. The diameter of the PP reactor riser (0.15 m) lies well within this range, whereas the diameter of the CP reactor riser (1.8 m) lies outside this range. Thus, the constant value of the Peclet numbers ( Pege and Pe pe ) for the PP reactor model is used as 10 as suggested by Bi (2004) whereas for the CP reactor model this value is tuned. The transport bed of both the MA CFB PP and CP reactors operate in the C-TFB regime. In the absence of a proper analysis of the gas and solid mixing in the C-TFB regime, the values of the Peclet numbers for the transport bed section ( Pege,1 and Pe pe,1 ) of both the PP and CP reactors are used as tuning parameters. Both the solid and gas effective dispersion coefficients are related to the corresponding Peclet numbers with the superficial gas velocities as reported by Bi (2004). It is again

 

noted that U 0,i used in Eq. (18) – Eq. (20) is measured at the outlet i.e., U 0,f i , as reported in the literature (Patience and Bockrath, 2009).

Dge, i 

U 0, i H T

Dpe, i 

U 0, i H T

Pege, i

Pe pe, i

(19)

(20)

18

3.3 Conditions for the reactor configuration complexities: The reactor model in the form of 2 NGC  3 1st order ODEs, given in Table 6, remains same in both the transport bed as well as the riser. The parameters and the model variables values change from the transport bed to the riser. The sudden change in the variable values also happen at the mixing points of the oxygen spargers in the transport bed. The model variable values also change due to the sudden change in the cross-sectional area between the transport bed and the riser. The relationships between the model variables at the interface of each section using the fundamental balances and the basic assumptions are discussed next. 3.3.1

Conditions at the mixing point of oxygen spargers:

The location of the fresh oxygen mixing point through kth (k = 1, 2 and 3) sparger is represented by O2 . The temperature, pressure, molar balances, etc., are conserved at each of the mixing points x  xIN ,k O2 O2 (i.e., at xIN and xIN ) for the steady state operation of the reactor. The fresh feed of oxygen ,k   ,k 

is assumed to be at the same temperature as the reaction mixture before and after the mixing point. The pressure at a given height, for assumed constant temperature, represents the constant total concentration at that height. The side additions of the fresh oxygen result in change of the relative concentration values of the different gaseous species so as to keep the total pressure/concentration constant before and after the mixing point. The molar balances for the lattice oxygen, ( NGC  1 ) gaseous species along with the total molar balance in the gas phase are given by Eq. (21) – Eq. (23). O2 O2 F L xIN , k     F L  xIN , k    O

(21)

O

Fj  x

O2 IN , k

O2  Fj  xIN , k         O2 IN , k  Fj  xIN , k     FO2 

O2 O2 IN , k FT xIN , k     FT  xIN , k     FO2

 for j  1, 2,..., NGC  1; j  O2     for j  O2 

(22)

(23)

19

It is assumed that the fresh oxygen feed through the spargers only affect the gaseous phase and not the lattice oxygen of the catalyst as evident from Eq. (21) and Eq. (22). This is valid if and only if QSSA is not assumed on  OL . The assumption makes both the convective FConv and dispersive components L O

FDisp of F L to be continuous before and after the mixing point, which are represented by Eq. (21) and L O

O

Eq. (24). The convective molar flow rates of the gaseous species other than the oxygen are not affected by the fresh oxygen feed as given in Eq. (25). The fresh O2 feed is assumed to be added only in the molar convective flow of O2 as given in Eq. (26a). It is noted that the kinetic parameters of the reactor model are further tuned in the proposed work. The reaction extent and thus the molar flow rate values predicted by the model depend on the values of the kinetic parameters in the tuning exercise. The values of the tuning parameters are varied in the tuning exercise. The specific values of these tuning parameters may lead to specific situations. The conditions applied in the model for the reactor configuration complexities should be applicable for such specific situations. For example, one such specific situation at the oxygen mixing points could be that the gaseous O2 is completely consumed before any of the oxygen mixing points. The specific situation of complete O2 consumption is achieved in Eq. (26a) when





O2 IN , k xIN the model predicted value of FOConv as compared to being greater than , k    becomes equal to FO2 2

it. After this situation Eqn. (26a) is not applicable as the convective molar flow just before the mixing point cannot be negative and should be zero. Thus, for this case the condition given by Eq. (26b) is used in the model in place of Eq. (26a).

FConv L L xINO2, k    FConv xINO2, k    O

(24)

O

O2 Conv FjConv xIN xINO2, k   , k     Fj

 for j  1, 2,..., NGC 1; j  O2 

FOConv xINO2, k    FOConv xINO2,k     FOIN2 , k 2 2 FOConv xINO2, k    0 2

if F x Conv O2

O2 IN , k

(25)

    FOIN , k  2

(26a) (26b)

20

The temperature and pressure continuity along with the Eq. (21) – Eq. (26) are the basis for 2 NGC  3 relationships between 2 NGC  3 model variable values just before and after the mixing points, given in Table 8. 3.3.2 Conditions at the interface of the transport bed and the riser: The interface between the transport bed and the riser is represented by x  xINT . The temperature and pressure are constant immediately before and after the interface (i.e., at xINT    and xINT    ), resulting in a constant total concentration at the interface according to Eq. (13). The molar flow rates are also conserved at the interface for the steady state operation of the reactor. The conservation of species molar balances for the lattice oxygen, the NGC  1 gaseous species along with the total molar balance in the gaseous phase are given by Eq. (27) – Eq. (29). There is no addition from the outside of reactor at the interface, and hence the mole fraction Y j of each of the ( NGC  1 ) gaseous species and  O is L

assumed to be continuous as given by Eq. (30) and Eq. (32). The temperature and pressure continuity along with Eq. (27) – Eq. (31) are the basis for the 2 NGC  3 relationships, given in Table 8 between 2 NGC  3 model variable values at the outlet of the transport bed and at the inlet of the riser.

F L xINT     F L xINT    O

(27)

O

Fj xINT     Fj xINT   ;

 j  1, 2,..., NGC 1

FT xINT     FT xINT   

Yj xINT     Yj xINT   ;

(29)

 j  1, 2,..., NGC 1

OL xINT     OL xINT    3.4

(28)

(30) (31)

Boundary conditions

(Abba et al., 2002) used the Danckwerts boundary conditions, as given in Eq. (32), in the dispersion model considering the variable gas density.

21

dC j  x  dx dC j  x  dx



HT C j  x   C IN  j  U Dge,1 IN 0

at x  0

0

(32)

at x  1

The modified Danckwerts boundary conditions for a variable density binary gas mixture were proposed by Choong et al. (1998). The modified conditions consider the dispersive molar flow rates of different species to be proportional to the gradients of their mole factions Choong et al. (1998). It is noted that the present reactor system is a MCRS with significant pressure drop. Thus, the modified Danckwerts boundary conditions given by Eq. (33) are used in this study. It is noted that for the limiting case of constant density, Eq. (32) are obtained from Eq. (33) but the reverse is not possible. A similar set of modified Danckwerts boundary conditions given by Eq. (34) for the solid phase component,  OL , is used. dY j  x  dx dY j  x  dx

 FjIN  H T FTIN  Y j  x   IN  ACS ,1CTIN Dge,1  FT 

at x  0

0

at x  1

G pz ,1 H T dOL  x   dx Dpe,1 cat 1  



d

L O



  FINL L O O  x    ACS ,1G pz ,1 N L   

 x  0

 for j  1, 2,..., NGC  1

(33)

at x  0 (34)

at x  1

dx

The modified Danckwerts boundary conditions have been used only for the NGC  1 components in the gaseous phase (Eq. (33)), thus, CT and U 0 at the inlet and outlet are related by Eq. (35).

ACS ,1U 0  x  CT  x   FTIN ACS ,2U 0  x  CT  x   FTOUT

at x  0 at x  1

(35)

The species molar flow rate of the NGC  1 gaseous species and the molar flow rate of  OL are given by Eq. (36) and Eq. (37) at the outlet. It is noted that the dispersive flow rates according to the modified Danckwerts boundary conditions are zero at the outlet. The pressure at the outlet of the MA CFB PP and

22

CP reactors is reported (Patience and Bockrath, 2010). The temperature is calculated at the outlet of both the reactors from the reported information (Patience and Bockrath, 2010). The total concentration at the outlet is then calculated by Eq. (5). FjOUT  ACS ,1U 0  x  CT  x  Yj  x 

at x  1

FOUT  ACS ,2Gpz ,2 N LOL  x  L

at x  1

 for j  1, 2,..., NGC 1

(36) (37)

O

4. Data Fitting Exercise The data fitting exercise is done in two steps, namely, model tuning and model validation. The available experimental data of both the MA CFB PP and CP reactors are divided into the training and validation sets. The model tuning is done by error minimization between the experimental data and the model predicted values using the training set. The error minimization results in the tuned values of the kinetic and hydrodynamic model parameters. The model validation is done by analyzing the accuracy of the reactor model predictions with the tuned parameters using the validation set. 4.1 Experimental data The experimental data for both the MA CFB PP and CP reactors are required for the data fitting exercise of the respective reactor models. The experimental data for both the reactors are reported in terms of the conversion of Bu ( X Bu ), the selectivity of MA ( S MA ), the mole fractions of some components at the inlet and/or of other components at the outlet, the lattice oxygen contribution (LOC) and the catalyst mass flow rate (Patience and Bockrath, 2010). The experimental data for both the reactors are divided in two parts: the training and validation sets. The ratio of these sets is kept at 2:1. The training and validation sets are made to have similar range of values for each of the reported information. To simulate a reactor-model, the information for all the NGC  1 components are required at either the inlet of the transport bed and/or at the outlet of the riser. Thus, the reported experimental data of both the reactors (PP and CP) are converted into the individual molar flow rates of each component at the 23

inlet ( FjIN ) and at the outlet ( FjOUT ). This is achieved by applying the material balances on a process flow sheet as shown in Fig. 2. It is again stressed that the modeling is done only for the two sections of the reactors, namely, the transport bed and the riser. Out of the reported information one information was relaxed to get a unique solution of the mass balances for the designed flow sheet. The LOC value was relaxed based on the assumption that the LOC measurements depend on the residence time distribution of the solid catalyst particles which is not reported in the experimental data (Patience and Bockrath, 2010). The reported ratio of the oxygen feed rates through the three spargers is 4:2:3 for the MA CFB PP reactor based on the reported number of nozzles at the each sparger (Hutchenson et al., 2010). The ratio of the oxygen feed rates used for the MA CFB CP reactor is 1:1:0.325 based on the reported information that the feed rates at the lower and the middle spargers of the CP reactor were of the order of 1000 kg/h, whereas at the upper sparger it varied between 250 – 500 kg/h (Hutchenson et al., 2010). 4.2. Model simulation 4.2.1 Numerical methods for the model simulation A total of 2 NGC  3 conditions are required to simulate the MA CFB reactor model consisting of 2 NGC  3 simultaneous 1st order ODEs given in Table 6. Many relations as given by Eqs. (33 – 37) and

discussed in the Section 3.4 may act as the boundary conditions. These conditions require the IN , Exp

information of the molar flow rates both at the inlet ( Fj

OUT , Exp

) and at the outlet ( Fj

) of the reactor.

The pressure and temperature are also available at the outlet of the reactor. It is noted that only some part of the experimental data, discussed in Section 4.1, of each experimental case can used as the required conditions for the model simulation because the remaining information is required for the error minimization in the model tuning step. It is also observed that the complete set of required 2 NGC  3 conditions at one point can only be obtained at the outlet of the reactor based on the available information and the use of Eqs. (33 – 37). The same is not possible at the inlet of the reactor. Hence, the

24

reactor model can be simulated as an initial value problem (IVP) only in the backward direction from the outlet to the inlet of the reactor with the initial conditions given in the first column of Table 9. It is evident, from Table 9, that the model simulation as an IVP uses the information of the molar flow rates only at the outlet ( FjOUT , Exp ) of the reactor. The reactor model can also be simulated as a boundary value problem (BVP) with the boundary conditions given in the second column of Table 9. It is evident from the second column of Table 9 that the model simulation as a BVP uses the information of the molar flow rates only at the intel ( FjIN , Exp ) of the reactor. Thus, different information in the form of molar flow rates, FjOUT , Exp and FjIN , Exp are required for the model simulation as an IVP and a BVP, respectively. Thus, the remaining information FjIN , Exp and FjOUT , Exp used in defining the error formulation are also different for the model simulations as an IVP and a BVP, respectively. The reactor model includes the interface between the transport bed and the riser along with the three oxygen mixing points as discussed above. The conditions at these points make the ODEs unstable at these points and thus a stiff ode solver ode15s of MATLAB® is used to simulate the model as an IVP. The shooting method is inappropriate for simulating such unstable ODEs as BVP (Shampine et al., 2000). Nonetheless, the algorithm bvp4c of MATLAB®, a collocation based method, is used to simulate the model as a BVP. The elemental profiles for all the four elements involved in the system, C, H, O and N, are observed to be constant in the model simulation both as an IVP and as a BVP. The constant elemental profiles validate the model equations and the conditions applied at the different interfaces given in Table 8. The simulated BVP and IVP results in the form of profiles for FjModel  x  are observed to be consistent in the presence of the appropriate conditions. The appropriate conditions, here, are the required boundary conditions for the BVP simulations obtained as a result from the IVP simulation and vice versa. The same consistent profiles for FjModel  x  are not observed when two different sets of 2 NGC  3 conditions, given in the two columns of Table 9, are used for the model simulation as an IVP and as a BVP as shown in Fig. 3 for j = Bu. This is because different sets of

25

information ( FjOUT , Exp and FjIN , Exp ) are used for simulating the model as an IVP and as a BVP. It is also observed that the profiles obtained from both the simulations are under predicting the corresponding experimental values but to different extents. The model simulation both as a BVP and an IVP have their own limitations, which makes both the methods suitable for different steps of the data fitting exercise. In the model simulation as a BVP the collocation based bvp4c algorithm of MATLAB® converts the system of ODEs into a system of nonlinear algebraic equations at different mesh points between the two boundary points (Shampine et al., 2000). To solve these nonlinear algebraic equations bvp4c algorithm requires an initial guess. The solution of the reactor model depends on the various model parameter values. The values of these parameters change at each iteration of the model tuning step while searching for the optimal values and, thus, at each iteration a different initial guess for the solution is required. Moreover, it is observed that bvp4c is not able to solve the reactor model with randomly chosen initial guess values. Thus, the model simulation as a BVP is found to be inappropriate in the iterative model tuning step of the data fitting exercise. The IVP code does not need any initial guess of the model solution and is also found to be almost 50 times faster as compared to the BVP code in terms of the CPU time for a given tolerance. Thus, the model simulation as an IVP is used in the iterative model tuning step of the data fitting exercise. This approach also has its own limitations. Some approximations are required because of the backward simulation which are discussed in the next section. 4.2.2 Extrapolation required in the model simulation as an IVP in model tuning exercise The backward simulation of the model as an IVP leads to an opposite trend in the variation of the mole fractions of the different components compared to the real trends observed with the forward flow of the gas and solid phases. The mole fractions of the components acting as the reactants (Bu, O2 and

 OL ) increase with the backward simulation as an IVP. On the contrary the mole fractions of the components acting as the products (MA, H2O, CO2 and CO) decrease. The extent of the decrease depends on the values of the model parameters. In model tuning step these parameter values keep

26

changing. The reactor model with a given set of parameter values sometimes over predict and sometimes under predict the extent of different reactions. The over prediction of the reaction extents leads to a faster decrease in the mole fractions of the products in the backward simulation of the model as an IVP in comparison to the experimental data. This leads to zero values for the mole fractions of some products before reaching the inlet of the reactor model. It is noted that if the model is simulated in the forward direction and that results in complete consumption of any of the reactants, the corresponding reaction rates automatically becomes zero and the concentration of the given reactant remains zero thereafter. In the backward simulation of the model if the mole fraction of any of the product becomes zero, it will not remain zero automatically for the remaining part of the reactor. This is because the reactants are still going to be present in the simulation of the reactor model. The specific rate expressions based on these reactant concentration values forces the simulation of reactor model to further decrease the amount of such products. The mole fraction of the products can be set to zero in the reactor model once they reach to zero value, but the model will not be able to penalize the over predictability of the reaction rates properly as it will give a fix value of the error irrespective of the extent of the over predictability of the reaction rates. Thus, the tuned model obtained from such a model tuning exercise will have an inherent bias for the given product. Thus, to penalize the over predictability similar to penalizing the under predictability of the reaction rates, the mole fractions of the products are kept decreasing in the model, even after they reach to zero value. This is just a mathematical exercise used in the model tuning step to achieve a reactor model with unbiased error in the prediction of its products and does not represent the real concentration values inside the reactor. These negative values for the products, which are not involved in rate expressions (H2O, CO2 and CO) as given in Table 3, does not lead to any infeasible reaction rates. On the contrary, the negative mole fraction of the product MA which also takes part in the series reaction leads to infeasible negative reaction rates. Thus, the reactor simulation as an IVP is stopped as soon as the mole fraction of the MA reaches to zero value at x  xconstraint ( YMA  xconstraint   0 ). As discussed earlier both over and under predictions are desired and,

27

thus, calculated to guide the model tuning exercise towards the optimal solutions. To achieve this the molar flow rates are extrapolated at the inlet of the reactor by using the slope of all the molar flow rates

Fj  xconstraint  at the reactor height xconstraint as shown in Fig. 4 for FMA . It is observed from Fig. 4 that FMA achieves the zero value before xconstraint when YMA  xconstraint  becomes zero. This is because of the dispersive term of FMA being considered along with the convective term. This constraint in the model simulation as an IVP is even possible with the converged tuned model parameter values. This is because the error is randomly distributed with positive and negative values for different experimental cases for obtaining a statistically better prediction. This constraint is not present in the simulation of the model as a BVP as shown in Fig. 4. This is because the molar flow rates in the BVP simulation are fixed at the inlet using the boundary conditions. This forces Fj  x  0  for all the products (MA, H2O, CO2 and CO) to have either positive or zero values. This ensures that in the BVP simulation the Fj  x  for products cannot achieve negative values as the reactions lead to further increase in Fj  x  0; j :products  . Thus, even though the model tuning exercise is performed by the IVP simulations, as it does not need any initial guess of the model solution and also takes significantly less time than the BVP simulations, the final analysis of the tuned model is done using the BVP simulations. Here, the IVP simulation provides the initial guess of the model solution for the BVP simulation. The model validation exercise is done by simulating the model as a BVP only with the obtained optimal model parameter values. 4.3 Model parameters for the model tuning exercise The model parameters in the MA CFB reactor model involve several kinetic and hydrodynamic parameters. The kinetic parameters are the pre-exponential factors ( k k0' and K eq0 , k " ; k ' =1, 2,…, N Rxn and

k " = 1: Bu, 2: O2, 3: MA), the activation energies ( E A, k ' ; k ' = 1, 2,…, N Rxn ) and the reaction enthalpies ( H eq , k " ; k " = 1: Bu, 2: O2, 3: MA) as given in Eq. (2) and Eq. (3). The reported values of these

28

parameters (Gascon et al., 2006) are given in Table 4. It is reported (Gascon et al., 2006) that these values are obtained using commercial VPO catalyst produced by Amoco in a fluidized bed reactor. On the other hand the experimental data used for the data fitting exercise are generated using the VPO catalyst produced by DuPont (commercially named as VPP catalyst) in both the PP and CP MA CFB reactors (Patience and Bockrath, 2010). Thus, the number of catalyst sites (per weight of the catalyst) on which the specific reactions are taking place might be different that may lead to differences in the preexponential factors. Thus, k k0' ( k ' = 1, 2,…, N Rxn ) and K eq0 , k " ( k " = 1: O2, 2: Bu) are used as the tuning parameters while the values of E A, k ' ( k ' =1, 2,…, N Rxn ) and H eq , k " ( k " = 1: O2, 2: Bu) are used as reported (Gascon et al., 2006) in the proposed reactor model. Keq0 , k "3:MA and H eq , k "3:MA are also used as two more tuning parameters in the absence of any reported values. The hydrodynamic parameters include a total of five parameters each for both the transport bed and the riser. Three of these parameters, ai , ci and d i are used in the correlations for the void fractions as given in Eq. (18) and the rest two are the Peclet numbers ( Pge, i and Pe pe, i ) representing the effective dispersion coefficients in the gas and solid phases. The hydrodynamic parameters of the MA CFB PP reactor, except the two Peclet numbers ( Pge,1 and Pe pe,1 ) for the transport bed, are already reported and discussed in section 3.2.2. 4.4 Objective formulation for the model tuning exercise The kinetic and hydrodynamic parameters in the MA CFB reactor model control the predicted reaction extent,  k (k = 1, 2, …, N Rxn ), of each reaction. In a MCRS with multiple reactions,

k

can be

related to the change in the molar flow rates for each components across the reactor Fj ( FjIN  FjIN ; j = 1, 2, …, NGC and  OL ) and the stoichiometry of the reactions. If the number of components in a MCRS are more than the number of the reactions involved, then each Fj can be expressed in terms of independent

 k (which are less than the number of components). If the number of reactions becomes

29

more than the number of the components, as observed in the proposed kinetic model (Table 3), then each Fj can be independently expressed in terms of the overall conversion of the reactants, X j  j  Bu, O2 and OL  , and the overall selectivity of the products, S j  j  MA, H2O, CO2 and CO  , as

given by Eq. (38) and Eq. (39), respectively. It is noted that the different definitions of selectivity and yields have been used interchangeably in the open literature (Carberry, 2001). The definition of S j used in this work is as defined by Coker, (2001) with Bu as the limiting reactant. The expressions for the experimental and the model predicted (both as an IVP and a BVP) X j and S j are given in Table 10.

Xj 

Sj 

FjIN  FjOUT FjIN FjOUT  FjIN FBuIN  FBuOUT

 j  Bu, O

2

and OL 

 j  MA, H2O, CO2 and CO 

(38)

(39)

It is noted that the expressions for experimental values of X j and S j ( X Exp and S Exp ) are the same for j j IVP or BVP model simulations but the expressions for the model values of X j and S j ( X Model and j ) are different. This is depicted by the two different profiles for Fj  x  obtained from the IVP and S Model j BVP model simulations as shown in Fig. 3. The input of the model simulations as an IVP and a BVP are FjOUT , Exp and FjIN , Exp , whereas the output of the model simulations as an IVP and a BVP are FjIN , Model

and FjOUT , Model , respectively. Thus, the calculated model values of X j , u and S j , u ( X Model and S Model ) j j given in Table 10 are different in the model simulations as an IVP and a BVP. The percentage error, E j , u , between the experimental value and the model value of X j , u and S j , u for a given component, j,

and experiment, u, is calculated using Eq. (40). It is noted that E j , u is not defined for the inert component (N2) as it will always remain zero. This makes the total number of responses, N Resp , to be equal to NGC with NGC  1 responses in the gas phase and one response in the solid phase. The objective for the model tuning, I Error , is defined as the average of the absolute values of E j , u for all the reactive

30

components and for all the experimental cases of the training set, as given by Eq. (41). It is noted that in a given experiment, u, E j , u value for one component (say, j = J) is usually linked with E j , u values of the other components (j ≠ J). Nonetheless, for any component E j , u value of one experiment (say, u = U) is completely independent to that of the other experiment (u ≠ U). This leads to a possibility that the model tuning done with the objective, I Error , may result in E j , u values for a given component averaged over all the experiments predominantly positive or negative. This average value for a given component, j, is proposed and defined as Bias j , given by Eq. (42). The target Bias j value for any component is zero for a good model prediction. Hence, to achieve this target, a second objective, I Bias , is defined in terms of the Bias j as given in Eq. (43). Thus, the model tuning exercise is done with these two objectives using a multi-objective genetic algorithm (GA) implemented in MATLAB® (gamultiobj). Model  X Exp j ,u  X j ,u  100  X Exp j  E j ,u  %    Exp Model 100  S j , u  S j , u  S Exp j,u 

I Error  Min

Bias j 

I Bias

 and OL      j  MA, H 2O, CO2 and CO  

 j  Bu, O

2

(40)

 N Resp N Exp  1     j ,u  %  N Resp  N Exp  j 1 u 1 

1 N Exp

(41)

N Exp

 E % u 1

(42)

j ,u

N  1  Resp  Min   Bias j  N Resp  j 1 

The default values of the GA parameters are used in gamultiobj. The population size of 10

(43)



number

of tuning parameters is used. Several GA optimization runs are made with the first run of 300 generations and subsequent runs of 150 generations till the convergence is reached with less than 1 % improvement in the objective values. In every GA run of 150 generations 35 % of the population obtained from the previous GA run is kept and the rest 65 % is randomly generated. The lower and

31

upper bounds for each of the 12 kinetic parameters are kept as 0.1 times to 10 times of the reported values for the first 300 generations. The bounds are relaxed in every 150 generations by dividing and multiplying by 3 in the minimum and maximum values of each kinetic parameter values obtained from the previous GA run, respectively. The lower and upper bounds for the two hydrodynamic parameters, Pge and Pe pe , for the transport bed of the PP reactor operating in the C-TFB regime are kept as 5 and

50. This is based on the reported values (Bi, 2004) for the risers operating with the G pz values less than 200 kg/m2s. The lower bound for the Peclet numbers for the CP reactor is decreased to 1 from 5 based on the square root ratio of the diameters of the PP and CP reactors (transport bed as well as riser) as discussed in Section 3.2.4. This makes the lower and upper bounds for Pge and Pe pe for the transport bed and the riser of the CP reactor as 1 and 50. The lower and upper bounds for the other hydrodynamic parameters ( ai , ci and d i ; i =1: transport bed, i =2: riser) are appropriately relaxed with each GA run based on the minimum and maximum values obtained in the previous run. 4.5 Procedure for the data fitting exercise The kinetic parameters are considered same for the PP and CP reactors whereas the hydrodynamic parameters are assumed to be different. Thus, a step by step data fitting exercise is proposed. In the Step 1 data fitting exercise, the PP data is used for the data fitting exercise of the PP reactor model. In the Step 2 data fitting exercise, the CP data is used for the data fitting exercise of the CP reactor model. Most of the kinetic (Gascon et al., 2006) and hydrodynamic parameters (Bi, 2004; Patience and Bockrath, 2009) for the PP reactor are reported except the two kinetic parameters ( Keq0 , k "3:MA and H eq , k "3:MA ) and the two hydrodynamic parameters ( Pge,1 and Pe pe,1 ). Two data fitting exercises as the

Step 1a and the Step 1b are performed. The Step 1a is performed on the PP data by relaxing the 12 reported kinetic parameters (as mentioned in Section 4.3) along with the two unknown hydrodynamic parameters as the tuning parameters. The Step 1b is performed on the PP data by relaxing all the 10 reported hydrodynamic parameters (as mentioned in Section 4.3) along with the two unknown kinetic

32

parameters as the tuning parameters. The tuned values of the kinetic parameters obtained from the Step 1 data fitting exercise on the PP data are then used without any further tuning in the Step 2 data fitting exercise on the CP data. The hydrodynamic parameters used in the Step 2 data fitting exercise on the CP data are all the hydrodynamic parameters ( ai , ci , d i , Pge, i and Pe pe, i ; i =1: transport bed, i =2: riser). Thus, in summary, in the proposed step by step data fitting exercise, the kinetic parameters of the CP reactor model are tuned from the PP reactor data on the PP reactor model. It is observed that the reported conversion and selectivity data trends in the PP and CP reactors are different. Hence, it is envisaged that the step by step data fitting exercise though logical may not be appropriate to fit all the data. Thus, another fitting exercise is performed as simultaneous data fitting exercise. This is done by fitting both the PP and CP data simultaneously on the PP and CP reactor models. The kinetic parameters are kept the same for the PP and CP reactor models whereas the hydrodynamic parameters are kept different for the two reactor models. In the end, the prediction capability of CP reactor models obtained from both the step by step and the simultaneous data fitting exercise are compare using the validation set of experimental data. 5. Results and Discussion: The MA CFB PP reactor model is first simulated as a BVP that does not involve any approximations using the reported kinetic and hydrodynamic parameters values both on the training and validation sets to get reference values for E j , u , Bias j , I Error and I Bias . The values of the two unknown hydrodynamic parameters, Pge,1 and Pe pe,1 are taken as the reported values of Pge, 2 and Pe pe,2 . The values for the two unknown kinetic parameters, K eq673,k "3 and H eq , k "3 are taken as the reported values of K eq673,k "2 and H eq , k "2 . Thus, the E j , u values obtained without any data fitting exercise are plotted against the

experimental values of X Exp and S Exp in Fig. 5a. This figure represents a comparison between the j j model predicted values and the reported experimental values. It is noted that multiplication factors (0.2 for S H 2O and 0.66 for SCO2 and SCO ) are used in the experimental values plotted on the horizontal axis

33

to shift the points horizontally. This exercise is done for a better visual analysis. The reference values of the Bias j for all the reactive components (j = Bu, O2,  OL , MA, CO2, CO and H2O) , I Error and I Bias without any data fitting exercise are shown in Fig. 6a and Fig 7a, respectively. The first observation made from Fig. 5a, Fig. 6a and Fig. 7a is that the model is behaving similar on the training and validation sets. This validates the division of the experimental data in the two sets and, thus, the model validation exercise based on such a partitioning can be accepted. The second observation made from these figures is that the model values are quite different from the experimental values without any data fitting exercise. The E j , u values are quite high (Fig. 5a), which makes the I Error values of the order of 50. The third observation made from Fig. 5a is that the E j , u values for any given component are either dominantly positive or negative for all the experimental cases. This observation results in significant Bias j values as shown in Fig. 6a. Positive Bias j values for all the reactants j = Bu and O2 and  OL in Fig. 6a represent that the model values of X j for all the reactants are being under predicted in comparison to the experimental values. Similarly, it can be said that model values of S j for j = MA and CO are being under predicted while for j = CO2 are being over predicted in comparison to the experimental values. It is noted that E j , u values for j = H2O are quite low. This can be attributed to the fact that the amount of H2O produced in the conversion of Bu to MA, CO2 or CO is almost same. The maximum value of the I Bias is equal to I Error . This equality of objectives happens when E j , u for any component j remain positively or negatively biased for all experimental cases ( u  1,2,..., N Exp ) as observed in Fig. 5a. Here, the significant Bias j values makes I Bias very close to the value of I Error . These observations suggest the need for the data fitting exercise of the model. Since the model tuning exercise is performed by simulating the model as an IVP, the reference values of the two objectives ( I Error and I Bias ) are also calculated by simulating the model as an IVP using the base values of the model tuning parameters. The values of the two objectives for both the training and validation

34

sets are shown as the first bars of Fig. 8a and Fig. 8b, respectively. It is observed by comparing the values from Fig. 7a and 8a that the objective values for the BVP and IVP simulations are different. This is attributed to the previous observation of different profiles Fj  x  for the BVP and IVP simulations shown in Fig. 3 and as discussed in section 4.4. The evolution of the model tuning exercise with the number of GA runs is shown in Fig. 9 for the Step 1a data fitting exercise. The number of generations shown in Fig. 9 represent the cumulative number of the generations after each GA run. The results are shown in the form of different converged Pareto sets after each GA run. A Pareto set represents a set of non-dominated solutions. Each solution of the Pareto set is non-dominated to each other because the improvement in any one of the objective (say, I Error ) is only possible at the cost of worsening the other objective (say, I Bias ). Usually a unique solution

is required to obtain a final set of the tuned model parameter values which represent the minimum I Error as well as minimum I Bias . In the absence of any process design constraints on which to make this decision, the solution closest to the origin in the Pareto set is chosen as the final solution of the model tuning exercise in this work. The normalized objective values between 0 and 1 are used in the final selection. The decision variables corresponding to the chosen solution represent the tuned model parameters. The model obtained with these model parameters is called the tuned reactor model. The tuned reactor model after each GA run are simulated both as an IVP and a BVP to calculate the values of the I Error and I Bias both on the training and validation datasets. The first observation with the convergence of the model tuning exercise is made based on the different Pareto sets shown in Fig. 9 and the objective values given in Fig. 8a corresponding to the selected solutions from the Pareto sets. It is observed that both the objective values decrease till 900 generations (5th iteration of the GA runs) and after that they are almost converged. The tuned model parameters are also observed to be converged after 900 generations with less than 1 % variation after these generations. The second observation is made on the guidance capability of the model tuning exercise. It is clear from Fig. 8a and Fig. 8b that the improvement in terms of the two objective values

35

on the training set with the number of generations is properly reflected on the validation set too. This observation is further verified by comparing Fig. 8a, Fig. 8b and Fig. 8c, Fig. 8d for the simulation of the model as an IVP and a BVP where similar trends are also observed. The use of an IVP simulation in the tuning exercise is validated based on these observations and the fact that an IVP simulation takes significantly less computational time than the BVP simulation. Though, it is noted that the BVP simulation is more accurate as it does not require any extrapolation while calculating I Error and I Bias . Thus, the two objective values using BVP simulations are used to further analyze the final converged results. The final converged results of Step 1a and Step 1b of the data fitting exercise on the PP data both for the training and validation data sets are shown in Fig. 7a. The I Bias values of the order of I Error signifies that the model with tuned parameters from Step 1b data fitting exercise still gives biased predications for most of the components. The Step 1a of PP data fitting exercise with twelve kinetic parameters and two unknown hydrodynamic parameters, reduces both the objective values significantly. This justifies the use of kinetic parameters for the model tuning in the Step 1a data fitting exercise on PP data. The converged values of E j , u and Bias j are also given in Figs. 5b and 6b. The comparison of Fig. 5a and 5b shows that the values of E j , u are significantly decreased and are also evenly distributed with positive and negative values after the (Step 1a) data fitting exercise. The same can be observed by comparing Fig. 6a and 6b as the Bias j values of each component significantly decreased after the data fitting exercise. A similar exercise is done on the CP data for getting the reference values before performing Step 2 data fitting exercise, as done on the PP data before the Step 1 data fitting exercise. This is done by simulating the CP reactor model as a BVP on the training and validation sets with the tuned kinetic parameters obtained from the Step 1a data fitting exercise. Since the analysis of the tuned model is done by simulating it as a BVP, the reference values are obtained by simulating the model as BVP only. The hydrodynamic parameters are also kept the same (as obtained for PP reactor model). E j , u values for the

36

CP reactor without any data fitting exercise are plotted against the experimental values of X Exp and j in Fig. 10a for both the training and validation sets. A similar exercise of horizontally shifting the S Exp j points for better visual analysis is implemented in this figure as done in Fig. 5 by using different multiplication factors (0.2 for S H 2O and 0.66 for SCO2 and SCO ). The simulated values of the Bias j values and the two objectives ( I Bias and I Error ) without any data fitting exercise on the CP data are given in Fig. 11a and Fig. 7b, respectively, for both the training and validation sets. It is observed in Fig. 10a that the E j , u values are significantly high and unevenly distributed with either positive or negative values for different components. This is also clear in Fig. 11a as the Bias j values of all the three reactants (j = Bu, O2 and  OL ) are negative while for products it is negative for MA and positive for CO2 and CO. The opposite trends in the Bias j values are observed in Fig. 6a and Fig. 10a. The PP reactor model with the use of the reported kinetic parameters for the PP data was under predicting the reactant conversions which is addressed with the use of the tuned parameters. It is observed that with the tuned parameter values the CP reactor model now over predicts the reactant conversions as shown in Fig. 9. It is noted that the hydrodynamic parameters of the PP reactor are used for the CP reactor and could be the reason for this observation. This validates the Step 2 data fitting exercise on the CP data by tuning all the 10 hydrodynamic parameters as mentioned in Section 4.5. The results of the step by step data fitting exercise (Step 2) for the CP data in terms of E j , u , Bias j and the two objective values calculated by the BVP simulation of the converged CP reactor model are shown in Fig. 10b, 11b and 7b, respectively. It is observed from Fig. 10b that even though the E j , u values decreased after Step 2 data fitting exercise but the values are either dominantly positive or negative for any given component. This is shown in Fig. 11b with significant positive and negative Bias j values. This makes the I Bias values of the order of the I Error values as shown by the third and

fourth bars in Fig. 7b. Thus, the reactor model for the CP reactor obtained from the step by step data fitting exercise shows significant bias in the form of over prediction for MA and under prediction for 37

CO2 and CO. The CP reactor tuned model with such inherent biases is not appropriate for optimization studies. The kinetic parameters tuned from the PP data fitting exercise are unable to fit the CP reactor data as evident from the previous discussion. This validates the simultaneous data fitting exercise using the PP and CP reactor data simultaneously on the respective models as discussed in section 4.5. The E j , u , Bias j and the two objective ( I Error and I Bias ) values for the PP data calculated using the PP reactor

model obtained from the simultaneous data fitting exercise are shown in Fig. 5c, 6c and the last two bars of Fig. 7a, respectively. Similarly, the E j , u , Bias j and the two objective ( I Error and I Bias ) values for the CP data calculated using CP reactor model obtained from the simultaneous data fitting exercise are shown in Fig. 10c, 11c and the last two bars of Fig. 7b, respectively. It is observed that the PP reactor models obtained from the step by step and the simultaneous data fitting exercise show similar behavior in terms of the E j , u and Bias j values. The CP reactor models obtained from the two different data fitting exercises show quite different behavior. It is observed by comparing Fig. 10b and 10c that the CP reactor model obtained from the simultaneous data fitting exercise gives the E j , u values evenly distributed with both positive and negative values, unlike the CP reactor model obtained from the step by step data fitting exercise. This results in much lower Bias j values as observed in Fig. 10c in comparison to the Bias j values shown in Fig. 10b. A similar observation is made in terms of the two objective values ( I Error and I Bias ). The comparison between the simultaneous and the step by step data fitting exercise in Fig. 7 shows significant improvement in the objective value, I Bias , for the commercial plant (Fig. 7b) at the cost of not so significant worsening of the objective values for the pilot plant (Fig. 7a). It is also noted that the model obtained from the simultaneous data fitting exercise shows better I Bias values for the validation set of the PP reactor data. Thus, the simultaneous data fitting exercise results in a CP reactor model with better prediction capability in terms of any inherent bias for any component. In summary, the tuned PP and CP reactor models from the simultaneous data fitting

38

exercise with lower I Bias values represent better prediction capability models than the models obtained from step by step data fitting exercise. The tuned kinetic parameter values from the simultaneous data fitting exercise are given in Table 4 and the hydrodynamic parameters are given in Table 11. It is observed that the tuned kinetic parameters are different than the reported values. This is attributed to the different catalyst used in the study for the reported kinetic model development than that used in the PP and CP MA CFB reactors as discussed in section 4.3. Ballarini et al. (2006) in their review reported that the activity of the VPO catalyst is affected by the preparation method. It is observed from Table 4 that the rate constants for the reactions involving the formation and consumption of  OS sites are predicted to be higher than the reported values. On the other hand, the rate constants for the reactions happening on the O2 sites are predicted to be less than the reported values. Moreover, the tuned values of K eq0 , k " and H k " values for k” = 3: MA adsorption - desorption equilibrium are of the order of the tuned values for k” = 1: O2 and 2: Bu adsorption - desorption equilibriums. This justifies the use of K eq0 , k " and H k " terms in the modified kinetic model. The tuned hydrodynamic parameters are given in Table 11. It is again noted that only two hydrodynamic parameters are tuned for the PP reactor whereas ten hydrodynamic parameters are tuned for the CP reactor. It is observed that the tuned value of Pege,1 for the PP reactor is near to its upper bound on the other hand the same for the CP reactor is on its lower bound. The tuned value of Pe pe,1 for the PP reactor is also greater than the same obtained for the CP reactor, but both the values are well within the bounds (1 – 50: CP reactor, as discussed in Section 4.4). This difference is attributed to the fact that the data fitting exercise tries to fit the differences observed in the PP and CP reactor data with the same kinetic and different hydrodynamic parameters as discussed before. The values of the Peclet numbers represent high gas and solid mixing in the transport bed of the CP reactor than the PP reactor. The high gas mixing in the riser of the CP reactor is also observed from the tuned value of Pege,2 for the

39

CP reactor. The solid mixing in the riser of the CP reactor is observed to be similar to that of the PP reactor as the tuned value of Pe pe,2 for the CP reactor is near to the reported value (Bi, 2004) for the PP reactor. The tuned values of Pege,1 and Pege,2 for the CP reactor are at their lower bound. The tuned values of Pege,1 and Pege,2 being at their lower bounds represent a significant gas mixing in the CP reactor. A sensitivity analysis is done by further decreasing these values (one by one and together). It is observed that 50 % reduction of these values resulted in less than 1 % change in the two objective values. The lower tuned values of a1 and a2 for the CP reactor in comparison to the reported values for PP reactor signify lower drift velocity Vg , i (Eq. 18c) in the CP reactor than that in the PP reactor. The lower tuned values of c1 and c2 along with the tuned values of d1 and d 2 for the CP reactor in comparison to the reported values for the PP reactor signify lower distribution coefficient C0,i (Eq. 18) in the CP reactor than in the PP reactor. This means the radial distribution of the solid void fraction and the axial gas velocity is more uniform in the CP reactor than that in the PP reactor. Several assumptions are made during the model development to keep it simple. The simplification of the expression for pressure drop from Eq. (11) to Eq. (12) is one of the such assumptions. The final tuned model is simulated with the pressure drop including all the frictional and acceleration term as given in Eq. (11). It is observed that the change in the two objective values ( I Error and I Bias ) are less than 1 % in the training and validation sets. This validates the simplification for the pressure drop expression. Negligible external and internal mass transfer resistances in the gas film and within the microspores of catalyst particles, respectively, are the other assumptions made in the model. Mears criterion (Mears, 1971) and Weisz-Prater criterion (Weisz and Prater, 1954) are generally used for checking the validity of these two assumptions, respectively. These criteria need the experimental rate data. The validity of these two assumptions in absence of such rate data is made by simulating the final tuned model with and without these assumptions. The R j  x  for different components considering the two mass transfer resistances are calculated from Eq. (44) in the place of Eq. (4). 40

R j  x  3

 r 1

rcat  0

2 cat



R j , rcat  x  drcat

(44)

R j , rcat  x  in Eq. (44) represents the rate for component j within the catalyst (present at reactor height, x)

at the radius rcat . A similar approach is followed to calculate R j , rcat  x  implementing the external and internal mass transfer resistances, as reported (Chaudhari and Gupta, 2012) in the model of the MTFBR . The correlation proposed by Gunn (1978) is used for the mass transfer coefficient in the CFB reactor as suggested by Breault (2006). No significant change is observed in the value of R j  x  calculated from Eq. (4) or Eq. (44). On the other hand, the simulation of the model considering these mass transfer resistances leads to significant increase (more than 100 times) in the computational time. This validates the model assumption of negligible external or internal mass transfer resistances for the catalyst particles. This further validates the use of kinetic parameter tuning exercise of the PP and CP reactor data. One of the salient feature of the proposed model is the consideration of variable gas density in the form of axial variation of CT  x  and U 0  x  . The axial profiles of CT  x  and U 0  x  predicted by the tuned MA CFB PP and CP reactor models are shown in Fig. 12 for one of the given experimental conditions of the PP and CP each. It is noted that the experimental conditions of the PP and CP reactors are different, which is also observed by two different profiles of CT  x  and U 0  x  predicted by the tuned PP and CP reactor models. It is observed from Fig. 12a that the decrease in CT  x  is of the order of 10 % in the transport bed and the riser of the PP and CP reactors. Similarly, it is observed from Fig. 12b that the increase in U 0  x  is of the order of 30% in the transport bed and 15 % in the riser of the PP and CP reactors. The sudden increase in U 0  x  at xINT is because of the sudden decrease in the crosssectional area between the transport bed and the riser. The linear variation of CT  x  shows that the first term in the expression for dCT  x  dx in Table 6 (due to pressure drop) is the dominant term. The

41

constant bed void fraction assumption results in a constant pressure drop and a constant dCT  x  dx ignoring the non-significant second term in the expression for dCT  x  dx in Table 6 (due to temperature drop). Two terms lead to the axial variation in the U 0  x  , the pressure drop and the overall reaction rate RT  x  , as discussed in Section 3.2.1. It is observed from the simulations of the tuned model that the

contribution of the reaction rate in the axial variation of U 0  x  is of the order of 10 % and the remaining contribution is due to the pressure drop. This is the reason why the profile of U 0  x  is observed to be almost linear in Fig. 12b. Thus, consideration of axial variations in CT  x  and U 0  x  in the model is validated from the tuned model simulations. The assumption of the uniform radial distribution both for the gas velocity and the solid hold up can be relaxed by considering radial dispersion along with the axial dispersion as suggested in the literature (Liu, 2001). The radial and axial dispersions in the reactor model result in partial differential equations (PDE). The same is attempted by Abba et al. (2003). A PDE model may represent the real reactor better than the proposed model but with significant increase in the computational time. Thus, with the given prediction capability of I Error under 17 % and I Bias under 6 %, the use of proposed simple MA CFB reactor model is recommended. 6. Conclusions A simple MA CFB reactor model is proposed in this work. The transport bed and the riser of MA CFB reactor are modeled assuming ideal operation of the regenerator at the steady state operation with the complete regeneration of the lattice oxygen in each cycle. The non-ideal plug flow, both for the gas and solid phases, is modeled with the help of dispersion coefficients ( Dge and D pe ) for both the phases in the proposed model. The reported kinetics Gascon et al. (2006) appropriate for the operating conditions of the MA CFB reactor is suitably modified. The transport bed and the riser of both the MA CFB PP and CP reactors are found to be operating in the C-TFB and the DSU regimes, respectively. The appropriate hydrodynamic correlations for the identified operating regimes of the MA CFB reactor are used in the proposed model. The density variation because of the significant pressure drop in the

42

DSU and the C-TFB regime is implemented in the dispersion model with a very simple approach. The proper conditions are implemented in the proposed reactor model for the configurational complexities of the reactor in the form of variation in its cross – section area between the transport bed and the riser of the MA CFB reactor. The appropriate conditions for the addition oxygen fed at the different heights of the transport bed are also implemented in the proposed model. The modified Danckwerts boundary conditions for the variable gas density are used in the proposed model. A detailed data fitting exercise is performed on the proposed MA CFB PP and CP reactor models using the reported PP and CP reactor data. The model simulation as an IVP is found to be more suitable for the model tuning exercise than the BVP in absence of the proper initial guess for the model solution. The final analysis of the tuned model is done using the model simulation as a BVP. The use of the two objectives, I Error and I Bias , is proposed for the model tuning exercise of a MCRS. The kinetic parameters are assumed to be same for both the PP and CP reactors whereas the hydrodynamic parameters are assumed to be different for the two reactors. Two data fitting exercises are performed, namely, the step by step and the simultaneous data fitting exercises. It is observed from the model validation exercise that the tuned PP and CP reactor models obtained from the simultaneous data fitting exercise show better prediction capability in comparison to the models obtained from the step by step data fitting exercise. Similar results obtained in all the model tuning and validation exercises also validate the use of the proposed two-objectives ( I Error and I Bias ) minimization for modeling the MA CFB reactors.

Funding This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

Literature Cited

43

Abba, I.A., Grace, J.R., Bi, H.T., 2002. Variable-gas-density fluidized bed reactor model for catalytic processes. Chem. Eng. Sci. 57, 4797–4807. Abba, I.A., Grace, J.R., Bi, H.T., Thompson, M.L., 2003. Spanning the flow regimes: Generic fluidized‐ bed reactor model. AIChE J. 49, 1838–1848. Andreux, R., Petit, G., Hemati, M., Simonin, O., 2008. Hydrodynamic and solid residence time distribution in a circulating fluidized bed: Experimental and 3D computational study. Chem. Eng. Process. Process Intensif. 47, 463–473. Avidan, A.A., Yerushalmi, J., 1982. Bed expansion in high velocity fluidization. Powder Technol. 32, 223–232. Bai, D., Issangya, A.S., Zhu, J.X., Grace, J.R., 1997. Analysis of the overall pressure balance around a high-density circulating fluidized bed. Ind. Eng. Chem. Res. 36, 3898–3903. Bai, D., Kato, K., 1995. Saturation carrying capacity of gas and flow regimes in CFB. J. Chem. Eng. Japan 28, 179–185. Bai, D., Zhu, J.-X., Jin, Y., Yu, Z., 1995. Internal recirculation flow structure in vertical upflow gassolids suspensions Part I. A core-annulus model. Powder Technol. 85, 171–177. Ballarini, N., Cavani, F., Cortelli, C., Ligi, S., Pierelli, F., Trifirò, F., Fumagalli, C., Mazzoni, G., Monti, T., 2006. VPO catalyst for n-butane oxidation to maleic anhydride: A goal achieved, or a still open challenge? Top. Catal. 38, 147–156. Berruti, F., Chaouki, J., Godfroy, L., Pugsley, T.S., Patience, G.S., 1995. Hydrodynamics of circulating fluidized bed risers: a review. Can. J. Chem. Eng. 73, 579–602. Bi, H.T., Fan, J., 1991. Regime transitions in gas-solid circulating fluidized beds. In: AIChE Annual Meeting. Los Angeles, California, pp. 17–22. Bi, H.T., Grace, J.R., 1995. Flow regime diagrams for gas-solid fluidization and upward transport. Int. J. Multiph. Flow 21, 1229 1236. Bi, X.T., 2004. Gas and solid mixing in high-density CFB risers. Int. J. Chem. React. Eng. 2, 1542– 6580. Breault, R.W., 2006. A review of gas–solid dispersion and mass transfer coefficient correlations in circulating fluidized beds. Powder Technol. 163, 9–17. Carberry, J.J., 2001. Chemical and catalytic reaction engineering. Dover Publications, New York. Chan, C.W., Seville, J.P., Parker, D.J., Baeyens, J., 2010. Particle velocities and their residence time distribution in the riser of a CFB. Powder Technol. 203, 187–197. Chaudhari, P., Gupta, S.K., 2012. Multiobjective optimization of a fixed bed maleic anhydride reactor using an improved biomimetic adaptation of NSGA-II. Ind. Eng. Chem. Res. 51, 3279–3294. Choong, T.S.Y., Paterson, W.R., Scott, D.M., 1998. Axial dispersion in rich, binary gas mixtures: model form and boundary conditions. Chem. Eng. Sci. 53, 4147–4149. Coker, A.K., 2001. Modeling of Chemical Kinetics and Reactor Design, 1st ed. Gulf Professional Publishing, Houston, Texas. Contractor, R.M., 1999. Dupont’s CFB technology for maleic anhydride. Chem. Eng. Sci. 54, 5627– 44

5632. Contractor, R.M., Bergna, H.E., Horowitz, H.S., Blackstone, C.M., Malone, B., Torardi, C.C., Griffiths, B., Chowdhry, U., Sleight, A.W., 1987. Butane oxidation to maleic anhydride over vanadium phosphate catalysts. Catal. Today 1, 49–58. Contractor, R.M., Sleight, A.W., 1987. Maleic anhydride from C-4 feedstocks using fluidized bed reactors. Catal. Today 1, 587–607. Danckwerts, P.V., 1953. Continuous flow systems: distribution of residence times. Chem. Eng. Sci. 2, 1–13. Gascon, J., Valenciano, R., Tellez, C., Herguido, J., Menendez, M., 2006. A generalized kinetic model for the partial oxidation of n-butane to maleic anhydride under aerobic and anaerobic conditions. Chem. Eng. Sci. 61, 6385–6394. Grace, J.R., Issangya, A.S., Bai, D., Bi, H.T., Zhu, J., 1999. Situating the high‐ density circulating fluidized bed. AIChE J. 45, 2108–2116. Gunn, D.J., 1978. Transfer of heat or mass to particles in fixed and fluidised beds. Int. J. Heat Mass Transf. 21, 467–476. Gupta, S.K., 1998. Modelling the hydrodynamics of large scale circulating fluidized bed risers. PhD Thesis, The University of Calgary. Hutchenson, K.W., La Marca, C., Patience, G.S., Laviolette, J.P., Bockrath, R.E., 2010. Parametric study of n-butane oxidation in a circulating fluidized bed reactor. Appl. Catal. A Gen. 376, 91–103. Issangya, A.S., 1998. Flow Dynamics in High Density Circulaing Fluidized Bed. THE UNIVERSITY OF BRITISH COLUMBIA. Issangya, A.S., Bai, D., Bi, H.T., Lim, K.S., Zhu, J., Grace, J.R., 1996. Axial solids holdup profiles in a high-density circulating fluidized bed riser. In: CFB5, 5th International Conf on Circulating Fluidized Beds. Beijing, pp. 1–7. Keraron, M., 2009. Conversion du méthanol en oléfines sur SAPO-34: modélisation cinétique et design de réacteur. Doctoral dissertation, École Polytechnique de Montréal. Liu, J., 2001. Particle and gas dynamics of high density circulating fluidized beds. Doctoral dissertation, The University of British Columbia. Lorences, M.J., Patience, G.S., Díez, F. V., Coca, J., 2003. Butane oxidation to maleic anhydride: kinetic modeling and byproducts. Ind. Eng. Chem. Res. 42, 6730–6742. Lorences, M.J., Patience, G.S., Díez, F. V., Coca, J., 2004. Transient n-butane partial oxidation kinetics over VPO. Appl. Catal. A Gen. 263, 193–202. Mears, D.E., 1971. Tests for transport limitations in experimental catalytic reactors. Ind. Eng. Chem. Process Des. Dev. 10, 541–547. Patience, G.S., Bockrath, R.E., 2009. Drift flux modelling of entrained gas–solids suspensions. Powder Technol. 190, 415–425. Patience, G.S., Bockrath, R.E., 2010. Butane oxidation process development in a circulating fluidized bed. Appl. Catal. A Gen. 376, 4–12.

45

Patience, G.S., Lorences, M.J., 2006. VPO transient oxidation kinetics. Int. J. Chem. React. Eng. 4, 1542–6580. Puchyr, D.M.J., Mehrotra, A.K., Behie, L.A., Kalogerakis, N.E., 1997. Modelling a circulating fluidized bed riser reactor with gas—solids downflow at the wall. Can. J. Chem. Eng. 75, 317–326. Pugsley, T.S., Berruti, F., 1996. A predictive hydrodynamic model for circulating fluidized bed risers. Powder Technol. 89, 57–69. Pugsley, T.S., Patience, G.S., Berruti, F., Chaouki, J., 1992. Modeling the catalytic oxidation of nbutane to maleic anhydride in a circulating fluidized bed reactor. Ind. Eng. Chem. Res. 31, 2652– 2660. Qi, M., 2012. Hydrodynamics and micro flow structure of gas-solid circulating turbulent fluidized beds. Doctoral dissertation, University of Western Ontario. Schlogl, R., 2009. Concepts in selective oxidation of small alkane molecules. In: Mizuno, N. (Ed.), Modern Heterogeneous Oxidation Catalysis: Design, Reactions and Characterization. Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, Germany, pp. 1–42. Schugerl, K., 1967. Experimental comparison of mixing processes in two- and three-phase fluidized beds. In: Drinkenburg, A. (Ed.), Proceedings Intern Symp on Fluidization. Amsterdam: Netherlands University, pp. 782–794. Shampine, L.F., Kierzenka, J., Reichelt, M.W., 2000. Solving boundary value problems for ordinary differential equations in MATLAB with bvp4c, Tutorial Notes, pp. 437–438. Shekari, A., Patience, G.S., 2013. Transient kinetics of n‐ butane partial oxidation at elevated pressure. Can. J. Chem. Eng. 91, 291–301. Shekari, A., Patience, G.S., Bockrath, R.E., 2010. Effect of feed nozzle configuration on n-butane to maleic anhydride yield: From lab scale to commercial. Appl. Catal. A Gen. 376, 83–90. Wang, C., Zhu, J., Barghi, S., Li, C., 2014. Axial and radial development of solids holdup in a high flux/density gas–solids circulating fluidized bed. Chem. Eng. Sci. 108, 233–243. Wang, D., Barteau, M.A., 2001. Kinetics of butane oxidation by a vanadyl pyrophosphate catalyst. J. Catal. 197, 17–25. Wang, D., Barteau, M.A., 2002. Oxidation kinetics of partially reduced vanadyl pyrophosphate catalyst. Appl. Catal. A Gen. 223, 205–214. Wang, D., Barteau, M.A., 2003. Differentiation of active oxygen species for butane oxidation on vanadyl pyrophosphate. Catal. Letters 90, 7–11. Wei, F., Lin, S., Yang, G., 1993. Gas and solids mixing in a commercial FCC regenerator. Chem. Eng. Technol. 16, 109–113. Weisz, P.B., Prater, C.D., 1954. Interpretation of measurements in experimental catalysis. Adv. Catal. 6, 143–196. Xiao-Feng, H., Cheng-Yue, L., Chen, B.H., Silveston, P.L., 2002. Transient kinetics of n-butane oxidation to maleic anhydride over a VPO catalyst. AIChE J. 48, 846–855. Yerushalmi, J., Cankurt, N.T., 1979. Further studies of the regimes of fluidization. Powder Technol. 24, 187–205. 46

Zhang, H.L., Degrève, J., Dewil, R., Baeyens, J., 2015. Operation diagram of circulating fluidized beds (CFBs). Procedia Eng. 102, 1092–1103. Zhu, H., Zhu, J., 2008. Gas‐ solids flow structures in a novel circulating‐ turbulent fluidized bed. AIChE J. 54, 1213–1223.

47

No solid circulation 𝑮𝒑𝒛 = 0 Fixed Bed

Solid circulation 𝑮𝒑𝒛 > 0 CAFB 𝑼𝟎 ↑↑

𝑼𝟎 > 𝑼𝑪𝑨

𝑼𝟎 ↑ Bubbling Fluidized Bed

𝑼𝟎 ↑

Turbulent Fluidized Bed

𝑼𝟎 ↑

𝑼𝟎 ↑

FFB 𝑼𝟎 ↑ 𝑮𝒑𝒛 ↑

𝑼𝟎 ↓

𝑼𝟎 ↑

𝑼𝟎 > 𝑼𝒇𝒇 𝑮𝒑𝒛 < 𝑮𝑺𝑪𝑪 𝒑𝒛

PC 𝑼𝟎 > 𝑼𝑷𝑪

𝑮𝒑𝒛 ↑ 𝑼𝟎 ↓ C-TFB

DSU 𝑮𝒑𝒛 ↓

𝑼𝟎 < 𝑼𝑪−𝑻𝑭𝑩

𝑼𝟎 < 𝑼 𝒇𝒇 𝑮𝒑𝒛 > 𝑮𝑺𝑪𝑪 𝒑𝒛

Fig. 1. Regime map for gas-solid system with and without solid circulation

48

𝑭𝑶𝑼𝑻 ; 𝒋 = 𝑴𝑨, 𝑯𝟐 𝑶 𝒋

Cyclone Separator Riser Stripper

𝑹𝒆𝒄𝒚𝒄𝒍𝒆 𝑹𝒂𝒕𝒊𝒐 = 𝑹𝑹

𝑭𝑶𝑼𝑻 𝜽𝑳 𝑶

Absorber Regenerator (Fluidized Bed)

Recycler

𝟏 − 𝑹𝑹 𝑭𝑶𝑼𝑻 𝒋 𝒋 = 𝑩𝒖, 𝑶𝟐 , 𝑪𝑶𝟐 𝑪𝑶, 𝑵𝟐

𝑭𝑶𝑼𝑻 ; 𝒋 = 𝑩𝒖, 𝑴𝑨, 𝑶𝟐 , 𝒋 𝑯𝟐 𝑶, 𝑪𝑶𝟐 , 𝑪𝑶, 𝑵𝟐 Riser

Regenerator Stripper

Sparger  𝒌 𝑭𝑰𝑵, 𝑶𝟐 𝒌=𝟑

𝒌=𝟐

Transport

𝒌=𝟏 𝑭𝑰𝑵 𝜽𝑳

Bed

Recycle Stream

𝒙 = 𝒙𝑰𝑵𝑻 𝑶 𝒙 = 𝒙𝑰𝑵,𝟐  𝟑

𝑹𝑹 𝑭𝑶𝑼𝑻 𝒋 𝒋 = 𝑩𝒖, 𝑶𝟐 , 𝑪𝑶𝟐 𝑪𝑶, 𝑵𝟐

𝑶

𝒙 = 𝒙𝑰𝑵,𝟐  𝟐 𝑶

𝒙 = 𝒙𝑰𝑵,𝟐  𝟏 𝒙=𝟎

𝑶

𝑭, 𝑰𝑵

𝑭𝒋

𝑭𝑰𝑵 𝒋 ; 𝒋 = 𝑩𝒖, 𝑴𝑨, 𝑶𝟐 , 𝑯𝟐 𝑶, 𝑪𝑶𝟐 , 𝑪𝑶, 𝑵𝟐 ; 𝒋 = 𝑩𝒖, 𝑵𝟐

Fig. 2. A schematic of MA CFB reactor

49

Fig. 3. Profile of FBu  x  obtained from the model simulation as an IVP and a BVP.

50

Fig. 4. Profile of FMA  x  simulating the same model as an IVP and a BVP. The extrapolation needed in IVP simulation when the model is over predicting the formation of MA.

51

Fig. 5. Residual plots of the model values as compared to the experimental values for the commercial plant; □: X Bu ;  : X O2 ;  : X  L ; y : SMA ; : SCO2 ; : SCO ;  : S H2O (blank: training set; filled: O

validation set) (a) Model values of PP reactor without any data fitting exercise on PP reactor data (b) Model values obtained from step by step data fitting exercise on PP reactor data (c) Model values obtained from simultaneous data fitting exercise on PP and CP reactor data combined.

52

Fig. 6. Bias j for different component and for both training and validation sets of MA CFB PP reactor. (a) obtained from PP reactor model without any data fitting exercise (b) obtained from PP reactor model using step by step data fitting exercise (c) obtained from PP reactor model using simultaneous data fitting exercise

53

Fig. 7. Results of the data fitting exercise in terms of the objective values ( I Error and I Bias ) calculated on both the training set (T. S.) and the validation set (V. S.) and by simulating the model as a BVP. (a) Results for the PP data (b) Results for the CP data. First two bars in both the figures represent objective values without any data fitting exercise. Next four bars for the PP data in (a) and two bars for the CP data in (b) represent objective values for the PP and CP reactor models obtained from the step by step data fitting exercise. Last two bars in both figures represent objective values for the PP and CP reactor models obtained from the simultaneous data fitting exercise.

54

Fig. 8. Evolution of tuning exercise as I Error and I Bias values calculated between the experimental data (training and validation) and the model simulation values (IVP and BVP) with the tuned parameters corresponding to the selected point from different Pareto sets converged after each GA run. Here, the number of generations represent the cumulative sum after each GA run.

55

Fig. 9. Evolution of tuning exercise as the converging Pareto sets with each iteration of optimization runs. Here, the number of generations represent the cumulative sum after each GA run.

56

Fig. 10. Residual plots of the model values as compared to the experimental values for the commercial plant ; □: X Bu ;  : X O2 ;  : X  L ; y : SMA ; : SCO2 ; : SCO ;  : S H2O (blank: training set; filled: O

validation set) (a) E j , u  %  values without any data fitting exercise on CP reactor data (b) Converged

E j , u  %  values obtained from step by step data fitting exercise on CP reactor data (c) Converged E j , u  %  values obtained from simultaneous data fitting exercise on PP and CP reactor data combined.

57

Fig. 11. Bias j for different components and for both the training and validation sets of the MA CFB CP reactor. (a) obtained from the CP reactor model without any data fitting exercise (b) obtained from the CP reactor model using the step by step data fitting exercise (c) obtained from the CP reactor model using the simultaneous data fitting exercise

58

Fig. 12. Axial profile of (a) CT  x  and (b) U 0  x  simulated using the tuned PP and CP MA CFB reactor models

59

Table 1 Correlations for the terms involved in the criteria for regime transition as regime map given in Fig. 1 References

U ff

(Gupta, 1998) (Grace et al., 1999) (Kim et al., 2004)

G

SCC pz

(Bi and Fan, 1991)

(Bai and Kato, 1995)

Correlations for the criteria of the regime transition mentioned in Fig. 1

 gU ff d p  G pz  12.55  g  v pz , t cat

  

0.55

  g  cat   g  gd p3    g2 

1.064   g g  cat   g  U ff  0.0113G1.192 pz  g  

U ff  v pz , t

 G pz  22.8   cat 1   v pz , t 





 G pzSCC   21.6    gU 0  gd p  

0.542

U0

G pzSCC d p

g

 U 0  0.125   gd p 

   

0.59

1.85

  g  cat   g  gd p3    g2 

U CA

(Bi and Fan, 1991)

U C TFB

(Qi, 2012) (Bi and Fan, 1991)

U PC  10.1 gd p 

U PC

 G pz   g

  

   

0.31

 dp     DT 

   

0.105

0.139

   

0.2

0.105

  g  cat   g  gd p3    g2 

  g  cat   g  gd p3 U CA  21.6 gd p    g2  UC TFB  0.0041Gpz  1.78 0.347

0.36

0.064

  g  cat   g  gd p3    g2 

   

   

   

0.63

 G pzSCC      gU 0 

 cat   g   g

  

0.542

  g  cat   g  gd p3    g2 

   

0.021

60

Table 2 Reactor and catalyst parameter values Reactor properties

PP reactor

CP reactor

Reference

HT ,1  m  , DT ,1  m 

11.15, 4.2

6, 0.3

(Patience and Bockrath, 2010)

HT ,2  m  , DT ,2  m 

28.5, 1.8

24, 0.15

(Patience and Bockrath, 2010)

zOIN2 ,1  m  , zOIN2 ,2  m  , zOIN2 ,3  m 

0.45, 1.9, 5.5 0.9, 2.1, 3.7 (Shekari et al., 2010)

 kg  3  m 

Cat 

1573

1573

(Patience and Bockrath, 2009)

d p  m

70

70

(Shekari et al., 2010)

N L  mmolO2 kg Cat 

310

310

(Patience and Bockrath, 2010)

61

Table 3 Modified Gascon’s kinetic model Oxidation of surface lattice site from adsorbed oxygen r1  k1O2 1  OS 

k O  2VS   2OS  V 1

2

Diffusion of lattice oxygen from subsurface lattice sites to surface lattice sites: r2  k2 OL  OS 

k OL  OS 2

Reactions happening at the adsorbed oxygen site: k3 6.5O2  C4 H10  g    4CO2  g   5H 2O  6.5V

r3  k3 pBu O2

k4 3.5O2  C4 H10  g   C4 H 2O3  g   4H 2O  3.5V

r4  k4 pBu O22

Reactions happening at the lattice oxygen site with adsorbed Bu and MA: 7  C4 H10  C4 H2O3  4H 2O  7 S O

k5

S V

k6 9OS  C4 H10   4CO  5H 2O  9VS  V

k7 13OS  C4 H10   4CO2  5H 2O  13VS  V

k8 6OS  C4 H2O3   4CO2  H 2O  6VS  V

r5 

k5 pBu OS 

2

1  Keq ,1 pO2  K eq , 2 pBu   K eq ,3 pMA  k6 pBuOS r6  1  Keq ,1 pO2  Keq , 2 pBu   K eq ,3 pMA 

r7 

k7 pBuOS

1  Keq ,1 pO2  Keq , 2 pBu  Keq ,3 pMA  k8 pBuOS r8  1  Keq ,1 pO2  Keq , 2 pBu  Keq ,3 pMA 

62

Table 4 Reported and tuned values (simultaneous data fitting exercise) of the kinetic parameters

Kinetic parameter†

Reported values (Gascon et al., 2006)

Final Tuned Values (Simultaneous data exercise fitting)

Kinetic parameter††

EA,1 104  kJ kmol 

Reported values (Gascon et al., 2006)

k10 106  kmol kg  s 

5.52

8.81

k20 106  kmol kg  s 

3.65

13.44

EA, 2 104  kJ kmol 

13.40

k30 106  kmol kg  bar  s 

2.33

0.57

EA,3 104  kJ kmol 

13.81

k40 106  kmol kg  bar  s 

14.3

1.41

EA, 4 104  kJ kmol 

5.03

k50 106  kmol kg  bar  s 

5.82

19.93

EA,5 104  kJ kmol 

3.24

k60 106  kmol kg  bar  s 

8.15

6.51

EA,6 104  kJ kmol 

10.61

k70 106  kmol kg  bar  s 

2.77

2.55

EA,7 104  kJ kmol 

8.20

k80 106  kmol kg  bar  s 

1.45

6.69

EA,8 104  kJ kmol 

13.79

Keq0 ,1 1 bar 

55.85

34.52

H1 104  kJ kmol 

-3.06

Keq0 ,2 1 bar 

94.213

5.50

H 2 104  kJ kmol 

-11.43

Keq0 ,3 1 bar 

-

7.97

H 3 104  kJ kmol 

-

-2.42

3.09

† Kinetic parameters which are used as tuning parameters in data fitting exercise †† Kinetic parameters which are used as reported and not used as tuning parameters in data fitting exercise

63

Table 5 Rate expressions for different components N Rxn

Rj

 k '1

r

k' j k'

RC4 H10

r3  r4  r5  r6  r7

RC4 H 2O3

r4  r5  r8

R O2

r1  6.5r3  3.5r4

R H 2O

5r3  4r4  4r5  5r6  5r7  r8

RCO2

4r3  4r7  4r8

RCO

4r6

R N2

0

ROL

r2

64

Table 6 Model equations as simultaneous 1st order ODEs dT  x  dx



dCT  x  dx

dU 0  x  dx

HT OUT Ti  Ti IN  Hi









f conv cat 1   i gHT Rg T  x 

U 0  x  dCT CT  x  dx

 x 





PT  x  dT  x 

RgT 2  x 



cat 1   i HT CT  x 

dx

RT

 d  dY j 1 dY j  x  dCT  x  U 0  x  H T   x     dx  dx CT  x  dx dx Dge, i 

G pz , i HT  d  dOL  x    dx  dx  cat 1   i Dpe, i











 dY  x  HT cat 1   i  j  Y j RT  x   R j  x   dx CT  x U 0  x  







  ;  j  1, 2,...NGC  1  

 d L   1   i HT  O  x   cat ROL  x    dx  G pz , i N L  

dY j d Yj  x    x  ;  j  1,2,...NGC  dx dx d L d L O  x   O  x  dx dx

65

Table 7 Parameter values used in Eq. (18)

ai

bi

ci

di

i = 1: transport bed 2.3 0

2.7 1

i = 2: riser

4.5 -1.7

1

0

66

Table 8 Conditions at the oxygen mixing points and at the interface of the transport bed and the riser Conditions at x =

Conditions at x = xINT

O2 xIN ,k

T    T   

T    T   

PT     PT   

PT     PT   

  FOIN2 , k O2  U 0  xIN  U  1   ,k  0   ACS ,1U 0    CT     

U 0      ACS , 2 ACS ,1 U 0   

OL     OL   

OL     OL    dOL   

dOL d L    O    dx dx

Yj   

dx

FT   

FT     FOIN2 , k

YO2    

dY j dx

  

FT     F

IN , k O2

dY j dx

dx

ACS , 2

2

pe , 2

dx

if FOConv     FOIN2 , k 2

Yj     Yj    ;  j  1,2,..., NGC 

if FOConv     FOIN2 , k 2

   ;  j  1, 2,..., NGC ; j  O2 

dYO2 dYO2

  1    D

ACS ,1 1   1 Dpe,1 dOL   

Y j    ;  j  1, 2,..., NGC ; j  O2 

FT    YO2     FOIN2 , k 0



dx



    dY    O 2

dx

if FOConv     FOIN2 , k 2

FOIN2 , k  U 0    HT  Conv IN , k  YO2      if FO2     FO2 Dgz ,1  FT    

dY j    dx



ACS ,1 Dge,1 dY j    ACS , 2 Dge, 2

dx

;  j  1, 2,..., NGC 

67

Table 9 Boundary and Initial conditions for simulating the model as a BVP and an IVP, respectively Initial Conditions for IVP

Boundary Conditions for BVP

(In terms of output molar flow rates)

(In terms of input molar flow rates)

PTOUT CT 1   Rg T OUT

PTOUT CT 1   Rg T OUT

T 1   T OUT

T 1   T OUT



U 0 1  

Y j 1  

ACS , 2CT 1 

dx dOL 1  dx

U 0  0  

dOL  0 

ACS , 2G pz , 2 N L

dx

O

0

0

 j  1,2,...,6 

FTIN

ACS ,1CT  0 

dY j  0  U 0  0  HT  dx Dge,1

FOUT L



dY j 1 

j  1,2,...,6

FTOUT

 1   L O

FTOUT

FjOUT





dY j 1  dx dOL 1  dx



 FjIN   Y j  0   IN  j  1,2,...,6 FT  

  FINL L  O O  0    ACS ,1G pz ,1 N L  1   1  

G pz ,1 HT Dpe,1 cat

0





 j  1,2,...,6 

0

68

Table 10 Expressions for experimental and model (for simulation as IVP and BVP) predicted X j and S j Experimental

X

Exp j



S Exp  j

FjIN , Exp  FjOUT , Exp FjIN , Exp FjOUT , Exp  FjIN , Exp FBuIN , Exp  FBuOUT , Exp

Model (IVP)

X

Model j



S Model  j

Model (BVP)

FjIN , Model  FjOUT , Exp FjIN , Model FjOUT , Exp  FjIN , Model FBuIN , Model  FBuOUT , Exp

X

Model j



S Model  j

FjIN , Exp  FjOUT , Model FjIN , Exp

FjOUT , Model  FjIN , Exp FBuIN , Exp  FBuOUT , Model

69

Table 11 Tuned values (simultaneous data fitting) of the hydrodynamic parameters of the PP and CP MA CFB reactor model Hydrodynamic Parameter MA CFB reactor

PP

CP

PP

Final tuned values

47.4

1.06

20.7

Pege,1

Pe pe,1

Pege, 2

Pe pe, 2

a1

c1

d1

a2

c2

d2

CP

CP

CP

CP

CP

CP

CP

CP

CP

6.09

1.52

10.73

1.86

0.77

1.16

0.69

0.22

-1.79

70