Mathematical Modeling and Numerical Simulation of a Fischer

Article pubs.acs.org/EF

Mathematical Modeling and Numerical Simulation of a Fischer− Tropsch Packed Bed Reactor and Its Thermal Management for Liquid Hydrocarbon Fuel Production using Biomass Syngas Tae Seok Lee and J. N. Chung* Department of Mechanical and Aerospace Engineering, University of Florida, Gainesville, Florida 32611, United States ABSTRACT: A mathematical modeling and numerical simulation study has been carried out for a liquid hydrocarbon production system using syngas from biomass gasification. The system is based on a Fischer−Tropsch packed bed reactor with an external heat removal sink. To ensure the accurate prediction of exothermic heat release, the chemical reaction kinetics is modeled with a comprehensive product distribution scheme based on a novel carbon number dependent chain growth concept and stoichiometric relationship between the syngas and hydrocarbons produced. The Fischer−Tropsch synthesis involves a three-phase phenomenon: gaseous phasesyngas, water vapor and light hydrocarbons; liquid phaseheavy hydrocarbon; solid phasecatalyst. A porous medium model has been used for the two-phase flow through an isotropic packed bed of spherical catalyst pellets. An Eulerian multiphase continuum model has been applied to describe the gas−liquid flow through the porous medium. Heterogeneous catalytic chemical reactions convert syngas into hydrocarbons and water. Intraparticle mass transfer limitation has also been considered in this model. Two independent validation cases have been demonstrated to show a successful implementation of the computational schemes. In this study, major attention has been paid to reactor temperature profiles because thermal management is highly important for the current exothermic catalytic reaction. We believe that our model could provide optimum operating conditions for achieving the maximum yield for desired products, liquid hydrocarbon fuels, without the thermal deactivation of the catalyst.

1. INTRODUCTION As the world faces significant energy supply and security challenges stemming from our dependence on petroleum and oil, the need for sustainable alternatives has been receiving great attention. To achieve energy security and independence in the near-future, and in the long run to prepare for postoil energy needs, the U.S. NSF−DOE Workshop report1 concluded that liquid biofuels produced from lignocellulosic biomass can significantly reduce our dependence on oil, create new jobs, improve rural economics, reduce greenhouse emissions, and ensure energy security. Further the report emphasized that the key bottleneck for lignocellulosic biomass-derived fuels is the lack of technology for the efficient conversion of biomass into liquid fuels. As a result, new technologies are needed to replace fossil fuels with renewable and sustainable energy resources. Reliable estimates of renewable and sustainable lignocellulosic forest and agricultural biomass and municipal solid waste (mostly biomass) tonnage in the U.S. range from 1.5 to 2 billion dry tons per year,1 so that these biomass resources could contribute 10 times more to our primary energy supply (PES) than they currently do. Another forecast claims that all forms of biomass and municipal solid waste have the potential to supply up to 60% of the total U.S. energy needs.1 Lignocellulosic biomass is the fibrous, woody, and generally inedible portions of the plants that are composed of hemicellulose and lignin. So, lignocellulosic biomass is nonfood based and does not compete with the production of food crops. The U.S. Natural Resources Defense Council has projected that an aggressive plan to make lignocellulosic biofuels in the U.S. could produce 7.9 million barrels of oil per day by 2050, or © 2012 American Chemical Society

more than 50% of current total oil use in the transportation sector.2 The synthetic biofuels produced by the Fischer−Tropsch process from lignocellulosic materials contain no sulfur, no particulates, no aromatics, and no nitrous compounds, thus making them very clean burning and reducing the production of acid rain. Because it has exactly the same chemical properties as fossil based diesel, it can be blended with regular diesel, stored, and distributed using the same infrastructure. Although chemically identical to fossil diesel, it has a higher cetane number and on a gallon for gallon basis contains 22% more energy. In general, the synthetic lignocellulosic diesel has up to 80% less combustion emissions, compared to petroleum diesel emissions, which include carbon dioxide, carbon monoxide, particulate matter, sulfur oxides, and hydrocarbons.3 Furthermore, these green biofuels are renewable, carbon-neutral, and sustainable. Therefore, a very promising route to liquid fuels, in particular, to synthetic lignocellulosic diesel, is the woody biomass gasification to synthesis gas (syngas; CO + H2), followed by the Fischer−Tropsch process to convert the syngas to hydrocarbon products. Fischer−Tropsch technology can be briefly defined as the means used to convert synthesis gas containing hydrogen and carbon monoxide to liquid hydrocarbon products. The hydrocarbons produced include oxygenated hydrocarbons such as alcohols. This technology was developed almost one century ago when Franz Fischer and Hans Tropsch were working to Received: October 26, 2011 Revised: January 8, 2012 Published: January 8, 2012 1363

dx.doi.org/10.1021/ef201667a | Energy Fuels 2012, 26, 1363−1379

Energy & Fuels

Article

one-dimensional, heterogeneous model for packed bed reactors with cobalt catalyst. Moutsoglou and Sunkara22 proposed a more comprehensive numerical simulation model to evaluate the Fischer−Tropsch (F-T) synthesis in a tubular multitube reactor packed with an iron-based catalyst. The effects of process parameters on product distribution were investigated. Three distinct regions of polymerization reactions were identified. Jess and Kern23 further developed a two-dimensional, pseudohomogeneous model with a pore diffusion limitation for a fixed packed bed reactor for both iron and cobalt catalysts, utilizing boiling water as the coolant. On the chemical reaction kinetics research, most of the very few previous studies were developed with a lumped kinetics model for syngas reactants. The lumped kinetics model has some inherent drawbacks. It cannot predict the exact amount of heat released by the exothermic reaction or the stoichiometric molar consumption ratio of hydrogen to carbon monoxide. From the thermal management point of view, the accuracy of the amount of heat release is of foremost importance. For example, even though the same 10 mol of carbon monoxide was consumed for each case, producing 1 mol of n-decane or 10 mol of methane will result in the release of 1560 kJ or 2060 kJ of heat under the standard pressure and temperature, respectively To predict the amount of heat released during the exothermic reaction, the current model obtained the product distributions and/or individual product generation rates based on the carbon number dependent chain growth probability theory from the carbide mechanism.4 Claeys and van Steen24 have reviewed the literature on the development of the F-T chemical synthesis model and summarized a general model that can predict the correct amount of heat production using the product distributions and/or individual product generation rates. The F-T synthesis converts syngas (hydrogen and carbon monoxide gases) into hydrocarbons and water in both gaseous and liquid phases (sometimes solid is also present but it is definitely the unwanted phase). The most wanted products are synthesis gasoline and/or diesel (this is why F-T synthesis was invented). Both products are definitely in the liquid phase under the standard conditions and may be so under most operating conditions. Moreover, none of these previous studies mentioned above included the two-dimensional flow model for the two phases. The main objectives of the current study are as follows: first, develop a rigorous computational thermal-fluid model with a sufficient chemical reaction kinetics that can predict heat release with a reasonable accuracy for the Fischer−Tropsch packed−bed reactor and, second, use this computational model to provide the thermal management criteria for the Fischer−Tropsch reactor design.

produce hydrocarbon molecules, from which fuels and chemicals could be made, using coal-derived gases in the 1920s.4 The Fischer−Tropsch synthesis has a long history and has been an active energy research area in recent years. As a result, large quantities of publications are available in the open literature. Two articles,5,6 in particular, provide an excellent and extensive review and discussion on the development and progress of the Fischer−Tropsch synthesis process and technology. The major problem in developing the Fischer−Tropsch (F-T) reaction kinetics is the complexity of its reaction mechanisms and the large number of species involved. Over the years, many apparently different mechanisms have been proposed, but common to them all has been the concept that a stepwise chain growth process is involved.7−13 This assumption is strongly supported by the fact that the carbon number product distributions calculated solely based on the probabilities of the chain growth matched the experimentally observed results obtained in different reactor types and sizes over a widely varying process conditions and with different catalysts.14

2. BACKGROUND AND OBJECTIVES There are four types of Fischer−Tropsch reactors in commercial applications at the present: the circulating fluidized bed reactor, the standard fluidized bed reactor, the fixed packed bed reactor, and the slurry-phase reactor. However, fluidized bed reactors are not suitable for producing liquid phase products because liquid phase products may cause catalyst agglomeration and a loss of fluidization.15 Fluidized bed systems are categorized as HTFT (high temperature Fischer−Tropsch) reactors, and the notable distinguished feature between the HTFT and the LTFT (low temperature Fischer−Tropsch) reactors is the absence of liquid phase in the HTFT reactor. The fixed packed bed and slurry-phase systems categorized as the LTFT reactors are appropriate for producing liquid phase products. The most common packed bed reactors are in shelland-tubes design. A high pressure operation due to narrow and long tubes, fine catalyst pellets, and high operating gas velocities will increase gas compression costs and also could cause disintegration of catalyst pellets. Catalyst loading and sometimes unloading will be difficult for narrow and long tubes. F-T synthesis reactions are exothermic, so heat removal is also an important factor. Reactor design considering only the kinetics aspect may have hot-spots causing catalyst deactivation, thermal runaway, or sintering. Although, with these drawbacks, the packed bed reactor is widely used for F-T process studies such as catalyst development, kinetics measurement, catalyst deactivation studies, and so on. Despite these drawbacks, a packed bed reactor has some benefits, such as easiness of its operation, no need for additional separation device, and straight scale-up for large scale reactors. Abundant experimental and modeling research efforts concerning slurry-phase F-T reactors are reported elsewhere.16 In contrast to slurry-phase F-T reactors, the literature on the packed bed reactor modeling and design is very limited. Atwood and Benett17 developed a one-dimensional plug flow, heterogeneous model to investigate the parametric effects on commercial reactors. A two-dimensional plug flow, pseudohomogeneous model without intraparticle diffusion limitations was developed by Bub and Baerns.18 Jess et al.19 developed a two-dimensional, pseudohomogeneous model for a nitrogen-rich syngas. A one-dimensional, heterogeneous model to account for the intraparticle diffusion limitations had been developed by Wang et al.20 De Swart et al.21 developed a

3. MATHEMATICAL MODEL An axisymmetric two-dimensional multiphase heterogeneous model has been developed to simulate a tubular fixed packed bed Fischer−Tropsch reactor. The system schematic and coordinates system are given in Figure 1. As shown in Figure 1a, the synthesis gas (syngas), a mixture of carbon monoxide and hydrogen, enters the inlet of the packed bed reactor with a mass flux, F, temperature, TIn and hydrogen to carbon monoxide molar ration, θH2/CO. The F-T reactor is housed inside the cylindrical tube in the form of a fixed packed bed of catalyst pellets that is assumed to be a porous material. The fluid flow 1364


Energy & Fuels

Article

Figure 1. (a) Schematic of the axisymmetric cylindrical packed bed F-T reactor and (b) external coolant flow configuration.

LHHW model was applied to obtain the product distribution and individual hydrocarbon production rates by relating each production rate with syngas consumption rate based on stoichiometric and chain growth probability from the carbide mechanism. 3.1. Gas−Liquid Hydrodynamics. 3.1.1. Continuity Equation. In the F-T synthesis, syngas is converted mainly to hydrocarbons in the form of a two-phase flow with both gaseous and liquid phases. Among the several possible flow regimes for the gas−liquid flow, it was pointed out in ref 26 that the flow type in the F-T packed reactor is in the form of liquid hydrocarbons trickling through the packed bed. So, for the numerical modeling, the mixture model, one of the Euler− Euler approaches, which consider each phase as an interpenetrating continuum, has been applied for the trickling liquid flows. This mixture model is a good, simplified alternative to the full Eulerian multiphase model because it solves each transport equation for the mixture and the volume fractions of all the phases. The continuity equations for the gaseous and liquid phases are shown below:

inside the tube is modeled as a two-phase flow through a porous medium that represents the packed bed of solid catalyst particles. Figure 1b shows the external cooling arrangement, where the coolant, with a temperature of TC, flows over the reactor tube perpendicularly in a so-called cross-flow condition. The detailed kinetics for producing linear-paraffin in both gaseous and liquid phases is derived based on the carbide mechanism chain growth probability and stoichiometry for non-Anderson−Schulz−Flory (ASF) distribution. The Eulerian mixture model has been adopted along with a semiempirical Ergun equation for modeling the porous medium. Mass transfer though the catalyst pellet pores is also considered by means of a general Thiele modulus. Main assumptions of the current model for the Fischer− Tropsch synthesis are listed as follows: (1) The flow field inside the tube is an axisymmetric and two-dimensional steady flow where catalyst pellets are packed inside and coolant flows outside the tube. (2) A steady state operation has been assumed; that is, there will not be any change over the time including catalytic activity, selectivity, and stability. (3) The two-phase flow is composed of gaseous (syngas, water vapor, and light hydrocarbon products) and liquid (heavier hydrocarbon) components. (4) Solid hydrocarbons in the form of wax have been neglected. (5) No VLE (vapor−liquid equilibrium) is assumed, as the liquid phase contains only heavier hydrocarbons with small mass fractions. (6) The packed bed is assumed to be statistically uniform, no channeling with isotropic hydrodynamic properties. (7) The production of oxygenates (alcohols, etc.) is neglected because of their small amounts compared to those of hydrocarbons. (8) Concerning the chemical kinetics model, it is assumed that the syngas consumption rate is primarily based on the lumped kinetics from the semiempirical Langmuir−Hinshelwood−Hougen−Watson (LHHW) model given by Yates and Satterfield.25 However, the

∂ (ρ ) + ∇· (ρmvm⃗ ) = 0 (1) ∂t m where ρm is the mixture density defined as ρm = χGρG + χLρL and χG and χL are the volume fractions of the gas and liquid phases, respectively, and vm⃗ is the mass-averaged velocity vector defined as vm⃗ =

∑k = G,L χkρkvk⃗ ρm

(2)

3.1.2. Momentum Equations. Similar to the continuity equation, the momentum equations for the general mixture model are accomplished by using mixture properties (density and viscosity) and mass-averaged velocity. It is assumed that the primary phase and dispersed phase do flow at the same velocities, which could be the actual case in the typical multiphase phenomena. However, as the multiphase mixture flows through the packed bed, the momentum equations need to be modified 1365


Energy & Fuels

Article

3.1.5. Species Transport Equation. The species transport equation for the gas or liquid phase is given by

as follows:

∂ (ρ vm⃗ ) + ∇· (ρmvm⃗ vm⃗ ) = −∇p + ∇·τm + SM ∂t m

∂ (ρ χ Yk , i) + ∇· (ρkχkvk⃗ Yk , i) ∂t k k = −∇· (χkJk⃗ , i ) + ∑ 9 jγjiMW, ji

(3)

where SM denotes a momentum sink vector term due to the packed bed. In this study, we have adopted a porous medium model that considers the solid phase as a momentum sink for the packed bed. A general momentum sink term for the homogeneous porous medium model is given by

⎛μ ⎞ δ SM, i = − ⎜ m vm, i + ρm|vm|vm, i⎟ ⎝ κ ⎠ 2

j

where Yk,i is the mass fraction of ith chemical species in kth phase, Jk,i means the diffusion flux of ith species in the kth phase as a result of concentration and temperature gradients. A total of N − 1 species transport equations will be solved if N is the total number of species inside the kth phase. The Nth equation for the kth phase, usually the most abundant chemical species in the corresponding phase, will be the sum of total mass fractions in the kth phase and it is unity; that is, ∑iYk,i = 1. 3.2. Fischer−Tropsch Catalytic Chemical Reaction Model. 3.2.1. Intrinsic Kinetics and Intraparticle Mass Transfer Limitation. Describing the F-T reaction kinetics is quite difficult due to its complicated reaction mechanisms and a large number of chemical species involved. Besides those problems, chemical kinetics studies are affected by many factors, considering that the F-T catalyst activity depends on its preparation method, metal loading, and support.28,29 F-T kinetics studies can be categorized by three different approaches: (1) Mechanistic proposals consisting of a sequence of elementary reactions among surface adsorbents and/or intermediates. (2) Empirical expressions of general power-law kinetics. (3) Semiempirical kinetics expression based on F-T mechanisms. In this study, we have adopted the third approach and chosen the widely accepted, well-known, semiempirical Langmuir−Hinshelwood equation for the kinetic expression proposed by Yates and Satterfield:25

(4)

where SM,i is the sink term for the ith momentum equation, |vm| is the magnitude of the mixture velocity, μm is viscosity of the mixture, ρm is density of the mixture,κ is the permeability of the porous medium, and δ is the inertial resistance factor. The momentum sink term is composed of two parts: a viscous loss term (dominant in a laminar flow, based on Darcy’s law) and an inertial loss term (dominant in a highvelocity flow). The permeability and inertial resistance factor can be obtained from the well-known semiempirical Ergun equation:27

κ=

Dp2

ε3 2

150 (1 − ε)

δ=

7 (1 − ε) 2 Dpε3

(5)

where Dp is the catalyst particle diameter and ε is the packed bed porosity. 3.1.3. Energy Equation. The conservation of energy equation for the mixture model is given below:

∂ [ ∑ χ ρ Ek ] + ∇· [ ∑ χkvk⃗ (ρkEk + p)] ∂t k k k k = ∇·(λ m∇T ) + SE

(6)

− rFT =

where λm is the mixture thermal conductivity, SE is the volumetric heat source from the exothermic reaction, and Ek is defined as follows:

⎧ hL liquid phase ⎪ Ek = ⎨ v2 p hG − + G gaseous phase ⎪ ⎪ ρG 2 ⎩

(7)

L‐phase j

i

(1 + KCCO)2

(10)

⎤ ⎛ − 37400 ⎞ ⎡ m6 ⎥ ⎟ ⎢ k = 0.4 exp⎜ ⎝ RT ⎠ ⎢⎣ kg mol s ⎥⎦ cat

(11)

⎛ 68500 ⎞ ⎡ m3 ⎤ ⎟ ⎢ ⎥ K = 5 × 10−9 exp⎜ ⎝ RT ⎠ ⎣ mol ⎦

(12)

For a packed bed reactor with a large number of catalyst pellets, the size range of the pellets may be wide, which will have an effect on the mass transfer of reactants through the catalyst pores. Therefore, the effective catalytic chemical reaction rate, in general, would be different from the intrinsic rate due to the scale effect. An effectiveness factor η, as defined below, is provided to take this effect into consideration:

∂ (χ ρ ) + ∇· (χLρLvm⃗ ) ∂t L L

∑ 9 j( ∑

kC H2CCO

We have selected the reaction constant k and adsorption constant K from Jess and Kern,23 among others.30−32 These constants are given below:

where hk is the enthalpy for phase k. 3.1.4. Volume Fraction Equation for the Liquid Phase. The volume fraction equation for the liquid phase can be obtained from the mass conservation equation for the liquid phase as follows:

=

(9)

γjiMW, ji) (8)

− rFT,eff = η( − rFT)

where 9 j denotes the jth reaction, γji and MW,ji represent the stoichiometric coefficient and molecular weight of the ith species in jth reaction, respectively. The sign convention for the stoichiometric coefficient is positive for products and negative for reactants, respectively.

(13)

where η is effectiveness factor defined as the ratio of the actual overall reaction rate to the reaction rate with no internal mass transfer limitation. It can be written as solely a function of the general Thiele modulus, regardless of the pellet geometry and 1366


Energy & Fuels

Article

3.2.3. Product Distribution Model. Based on the same general idea of the chain growth probability concept,24 we have developed the following product distribution model. The total carbon monoxide consumption rate is definitely related to the hydrocarbon production rate, since the F-T synthesis could be considered to consist of a lot of parallel reactions for each hydrocarbon formation reaction. The stoichiometric information is required for relating the molar change of hydrocarbon with that of carbon monoxide. Total carbon monoxide consumption rate, is equal to the summation of each individual production rate considering the stoichiometry condition that is represented by eq 18 below:

reaction order as follows:

η=

tanh φ φ

(14)

where ϕ is the general Thiele modulus given by

φ = Lc

k pseudoρpHH2 Deff RT

(15)

where Lc is the characteristic length, kpseudo is the pseudo-first order rate constant, ρp is the pellet density, HH2 is the Henry coefficient (approximately 20 000 Pa m3/mol) and Deff is the effective diffusion coefficient. The FLUENT code,33 a commercial software package employed to perform the numerical computation, does not have a built-in capability to evaluate the effectiveness factor with the general Thiele modulus nor heterogeneous chemical reaction between phases. Therefore, a C/C++ code has been written to generate the effectiveness factor, Thiele modulus and heterogeneous chemical reaction, those were submitted to FLUENT through a UDF (user defined function).34 The details of the FLUENT package used are discussed later. 3.2.2. General Carbon Number Dependent Chain Growth Probability and Deviation from Anderson−Schulz−Flory (ASF) Distribution. The hydrocarbon product distribution is inseparable from the syngas consumption. Furthermore, hydrocarbon product distribution plays a critical role in the thermal management of the F-T reactor, as we mentioned in the very beginning that the amount of heat produced is different when different hydrocarbons are formed, even though the amount of carbon monoxide consumed is the same. In addition to these, rigorous product distribution information is needed for modeling or analyzing the F-T reactor, in which the concentration profiles for reactants are changing with time and along the axial distance from the entrance. The most wellknown and simplest product distribution model is the Anderson−Schulz−Flory (ASF) distribution given below:

mn = (1 − α̅)α̅ n − 1

− rCO = r1 + 2r2 + 3r3 + ··· =

i

= =

rn rn − 1

= =

rdesorption, n rdesorption, n − 1 (1 − αn)αn − 1rchain growth, n − 2 (1 − αn − 1)rchain growth, n − 2

(1 − αn) αn − 1 (1 − αn − 1) (19) Using successive substitutions, we can relate any hydrocarbon desorption rate to methane production rate as follows: r1 r2 = (1 − α 2)α1 (1 − α1) =

(16)

r3 = (1 − α3)α 2α1

r1 (1 − α1)

rn = (1 − αn)αn − 1αn − 2···α1

(20)

r1 (1 − α1)

where the common term, r1/(1 − α1), could be interpreted as a chain initiation process, which facilitates the formation of an adsorbed methyl from the building block, CH2, by adding another hydrogen atom. Thus, it is defined as the chain initiation rate, rini, as follows:

rchain growth, n rchain consumption, n rchain growth, n rdesorption, n + rchain growth, n rchain growth, n rchain growth, n − 1

(18)

where ri denotes production rate for the hydrocarbon with a carbon number i. To achieve our objective of building a comprehensive two-dimensional heterogeneous numerical simulation capability, we have derived each hydrocarbon production rate with the general carbon number dependent chain growth probability instead of the carbon number independent constant chain growth probability as a non-ASF distribution model. From the definition of the general chain growth probability, eq 17, the ratio of hydrocarbon production rates between carbon number n and carbon number (n − 1) is given as

where mn is the mole fraction of the hydrocarbon with a carbon number n and α̅ is the carbon number independent constant chain growth probability. In the current approach, we adopted the carbon number dependent definition of the chain growth probability, as given in eq 17, that depends on the particular hydrocarbon produced.

αn =

∑ iri

r1 = rini (1 − α1) (17)

(21)

Any individual hydrocarbon production rate, hence, can be written as follows, using the chain initiation rate:

where αn is the carbon number dependent chain growth probability for the hydrocarbon with a carbon number n. In summary, in this work, we are using the specific carbon number dependent chain growth probability rather than a simple and carbon number independent chain growth probability.

rn = (1 − αn)αn − 1αn − 2···α1rini n−1

= (1 − αn)[ Π αk]rini k=1

1367

(22)


Energy & Fuels

Article

(FVM). FLUENT can faithfully discretize the governing equations and constitutive equations with corresponding initial and boundary conditions described in section 3 and then solve the resulting simultaneous equations with the maximum possible accuracy. It can handle the two-phase flow through a porous media and the heterogeneous chemical reaction in the packed bed reactor. The GAMBIT 2.4.6 package was used for grid generation. Because the reactor is cylindrical and packed with isotropic spherical catalyst beads, the computational domain was assigned as a two-dimensional axisymmetric cylinder. Therefore, quadrilateral computational cells were generated by the GAMBIT. A pressure-based solver employing the SIMPLE algorithm was used for the pressure velocity coupling scheme.39 As we described in the previous section, the Eulerian mixture model was applied for the two-phase flow through a porous media. Laminar flow is assumed as a result of very small length scales in the pathways created among small spherical catalyst pellets. As we have assumed no vapor−liquid equilibrium, the gaseous phase consists of carbon monoxide, hydrogen, water vapor, and light hydrocarbon up to C6, while heavier hydrocarbons, from C7 to C15, make up the liquid phase. From the thermodynamic properties of pure species, boiling points at 20 bar for n-hexane and n-heptane are 206.62 and 245.08 °C, respectively. The maximum temperatures of all the thermally stable cases investigated in the current paper are always less than 245 °C. F-T synthesis reactions are treated as interactions between the two phases in the presence of the fixed bed catalysts. However, FLUENT does not have built-in reaction kinetics such as eqs 10 and 24. Therefore, the implementation of the F-T synthesis reactions with rates given by eqs 10 and 24 has been accomplished by using a UDF (user defined function) in FLUENT (ANSYS FLUENT UDF Manual, 2009).34 Our UDF also provides intraparticle mass transfer limitation defined previously by eqs 13−15, as well as the F-T synthesis kinetics discussed in section 3. To achieve higher accuracy results, a second-order upwind discretization scheme was applied, except for the volume fraction equation where the QUICK scheme was applied.40,41 All materials (gas species, liquid species, and solid catalyst) are assigned appropriate properties from the literature, as well as from the FLUENT database. The properties of the gas species (density, viscosity, thermal conductivity, and specific heat) are allowed to vary with their respective temperatures. The thermodynamics properties of the gas phase mixture are calculated from their pure substance properties and local compositions (ideal gas law for density, ideal gas mixing law for viscosity and thermal conductivity, and mixing law for specific heat). It is noted that the volumetric heat generation term, SE, in eq 6 represents the heat released due to the exothermic chemical reaction. Figure 2 provides the standard heats of reactions42 during the production of n-paraffins and α-olefins in both gaseous and liquid phases. For the liquid phase heats of reactions between n-paraffin and α-olefin, we have computed the percent differences as 11.08%, 9.88%, 8.6%, 7.88%, 7.2%, 6.53%, 6.1%, 5.7%, and 5.3% for C7, C8, C9, C10, C11, C12, C13, C14, and C15, respectively. On average, the difference is 7.6%. From the point of total heat released for the thermal management analysis, by not distinguishing olefins from paraffins in the current model, and lumping them together for the heat release estimation, the error would be reasonably minor and the thermal runaway criteria predicted would be on the conservative side, as the heats of reaction are slightly higher for the paraffins than those of the olefins.

However, the measurement of the chain initiation rate could be extremely challenging. Without forming carbon dioxide, which is a reasonable assumption for the cobalt catalyst, all the adsorbed carbon monoxide molecules would be dissociated and the dissociated carbon successively gains hydrogen atoms to form a monomer, CH2, in a carbide mechanism. Then, this monomer or a building block should face two possibilities in its reaction process toward forming a hydrocarbon: act as an initiator or participate in the chain growth reaction with a pre-existing adsorbed alkyl group. Therefore, the chain initiation rate and all chain growth reaction rates should be balanced with the carbon monoxide consumption rate. In terms of a balance equation, the relationship between carbon monoxide consumption rate and chain initiation rate is as follows: N

−rCO = rini +

∑ rchain growth, i i=1

= rini + (α1rini + α1α 2rini + α1α 2α3rini + ···) N

= rini[1 +

i

∑ ( Π αj)] i=1 j=1

(23)

In the above, N is the highest carbon number species assumed in the products. Combining eqs 22 and 23 yields a specific hydrocarbon production rate involving the carbon monoxide consumption rate using individually assigned chain growth probability values as follows: n−1

(1 − αn)[ Π αk]( − rCO) rn =

k=1

p

1 + ∑N p = 1 ( Π αq) q=1

(24)

4. NUMERICAL SOLUTION METHOD AND VALIDATION 4.1. Numerical Simulation by FLUENT. As outlined in section 3, the synthesis process of Fischer−Tropsch (F-T) chemical reactions is complicated because the sophisticated heterogeneous catalytic reactions are intensively coupled with a two-phase flow in a packed bed, together with simultaneous heat and mass transfer. Therefore, solving this process accurately by numerical computations is a very challenging task. A numerical simulation of this process including the entire system infrastructure can provide some guidance for the design, scaleup, and optimization of a F-T reactor. A numerical simulation of the Fischer−Tropsch reactor has been accomplished using a commercial thermal-hydraulic code FLUENT (ANSYS-FLUENT 12).33 For the past several years, FLUENT has been widely accepted as the main computational software package for the numerical simulation of thermal hydraulic transport phenomena. For example, FLUENT was used for the numerical study of ignition/combustion process of pulverized coal.35 Jin and Shaw36 performed a computational modeling of n-heptane droplet combustion in an air−diluents environment under reduced-gravity using the FLUENT package. Chein et al.37 predicted the hydrogen production in an ammonia decomposition chemical reactor using the FLUENT software. Ho et al.38 used the FLUENT software to simulate the two-phase flow in a falling film microreactor. The discretization method in FLUENT is based on a Finite Volume Method 1368


Energy & Fuels

Article

Figure 3. Product distribution comparison between current model prediction and experimental results by Elbashir and Roberts;42 nonASF distribution, logarithmic value of normalized hydrocarbon product weight fraction versus carbon number.

Figure 2. Heats of reaction for the production of paraffin and olefin.

obvious that we have to use a carbon number dependent chain growth probability in our model to compare with the nonASF distribution experimental result. Even though Elbashir and Roberts43 found that the deviation from the ASF distributions is limited up to C5 hydrocarbon, we assigned the carbon number dependent chain growth probability values up to C7. The carbon number dependent chain growth probability values used in our model for comparison are given as follows:

Because a large-scale numerical simulation requires huge computational resources, therefore, the grid independence study was performed considering computational resource effectiveness aspects. Grid refinement has been performed until smaller grids do not significantly improve the accuracy. 4.2. Model Validation. Two independent approaches have been conducted to validate our model. The first validation was made by comparing the product distributions predicted by our model with the experimental work reported by Elbashir and Roberts.43 In their study, hydrocarbon product distributions are provided for both supercritical fluid and conventional gas-phase Fischer−Tropsch synthesis over a 15% alumina-supported cobalt catalyst, Co/Al2O3, in a high-pressure fixed packed bed reactor system. The operating condition for the presented case is for a 50 sccm syngas flow rate, a reaction temperature of 250 °C, a syngas partial pressure of 20 bar, and a H2/CO molar feed ratio of 2. Their measured data of products by two online gas chromatographs were used for the selectivity and conversion calculations. The sample analyzed in determining the product distribution was collected after the activity, and the selectivity of the cobalt catalyst showed a steady performance. Their typical results for the conventional F-T synthesis product distributions are shown in Figure 3, where the abscissa is the carbon number and the ordinate is ln(Wn/n). The n is the carbon number and Wn is the weight fraction of the hydrocarbon product with a carbon number of n. The logarithmic value of Wn/n is linear with the carbon number in a semilog graph for the constant chain growth rate, which is the well-known ASF-distribution. Their typical experimental results did not show a complete ASF distribution. Their results showed that the range of deviation from the standard ASF distributions (linear behavior) is limited to the light hydrocarbon (below C5) products region, while the distribution for heavier hydrocarbons above C5 follows well with the standard ASF distribution with a chain growth probability of 0.8.43 With a carbon number independent chain growth probability, the methane selectivity is underestimated and those for C2−C5 hydrocarbons are overestimated. It is

⎧ 0.292 ⎪ αn = ⎨− 0.0317n + 1.0362 ⎪ ⎩ 0.8

(n = 1) (2 ≤ n ≤ 7) (n > 8)

(25)

Chain growth probability values for higher hydrocarbons (above C8) are assigned as carbon number independent. Our FLUENT computational result for the validation is also plotted on Figure 3. The comparison between the current model and the experimental work is based on products distribution that is shown in Figure 3. In the experiment, Elbashir and Roberts43 had provided products distributions up to the carbon number 30. In our simulations, including the comparison case, the carbon number has been limited by our computing capacity. Eighteen chemical species are involved in our simulation, these are CO, H2, H2O, and C1−C15 hydrocarbons. Therefore 15 production rates are evaluated based on our model, eq 25. In the evaluation of 15 hydrocarbon production rates using eq 24, the highest carbon number, N, was set at 30 instead of 15 to increase calculation accuracy for a better comparison with experimental results. We found that, for the products distributions, there is a good agreement between our model and the experimental work by Elbashir and Roberts,43 as shown in Figure 3. On the basis of the close comparison, we believe that our model and the methodology are realistic and correctly implemented by the FLUENT software. 4.3. Model Comparison. The second approach was conducted to compare the temperature profiles from the current 1369


Energy & Fuels

Article

Figure 4. Reactor bed temperature profile comparison between current results and those of Jess and Kern.23

comprehensive model are very similar but the simplified twodimensional model generally predicts slightly higher temperatures than our comprehensive model. The percent deviations for the peak temperatures between the two simulations are 1.05%, 1.43%, 2.53%, and 5.70% for coolant temperatures of 200 °C, 205 °C, 210 °C, and 214 °C, respectively. The differences in the maximum temperatures are increasing with the cooling temperature, but the temperature runaways happen at the same coolant temperature of 215 °C. It is shown that six reactor bed temperature profiles using more refined coolant temperatures were calculated by the current model and are given in Figure 5 where we found that the thermally viable condition and the temperature runaway condition are divided between a narrow coolant temperature range between 214.1

model with those by the simplified two-dimensional pseudohomogeneous model developed by Jess and Kern.23 They have simulated an industrial scale multitubular Fischer− Tropsch reactor based on both one-dimensional and twodimensional pseudohomogeneous models, taking into account the intrinsic kinetics for two commercial catalysts of iron and cobalt, the intraparticle mass transfer limitations, and also the radial heat transfer inside the fixed bed and to the cooling medium outside. The differences between the two models and approaches are tabulated in Table 1. In Figure 4, we have Table 1. Comparison between Current Model and That of Jess and Kern23 Jess and Kern23 computational domain reactor kinetics products distribution heat of reaction component phases physical properties intraparticle momentum eq pressure concentration

this study

axisymmetric 2D

axisymmetric 2D

intrinsic kinetics for syngas N/A

intrinsic kinetics for syngas

lumped heat of reaction pseudohomogeneous

based on general chain growth probability and stoichiometry all specified individual heat of reactions heterogeneous

representative

function of composition

mass transfer considered not considered N/A

mass transfer considered porous material correction ideal gas

compared our reactor bed center-line axial temperature profiles with those reported by Jess and Kern23 in their twodimensional simulations using the cobalt catalyst for five different coolant temperatures. From the influence of the cooling temperature on the axial temperature profiles in the tubular packed bed reactor we can deduce that the trends predicted by the simplified two-dimensional model and current

Figure 5. Detailed reactor bed temperature profiles for various coolant temperatures. 1370


Energy & Fuels

Article

and 214.2 °C. Syngas conversions, XCO = (Ṅ CO,In − Ṅ CO)/ Ṅ CO,In, are also compared in Figure 6. In the simplified two-

numerically in an F-T reactor23 and observed in the experiment of gas-phase polymerization of ethylene44 and ethylene oxidation catalytic chemical reaction in a packed bed reactor.45 5.1. Baseline Case. For the scale-up and optimization of the syngas to liquid hydrocarbon F-T system, modeling and simulation are the first step of the process. As a result, the basic performance characteristics of the packed bed reactor with cobalt catalyst are examined under various system parameters. To facilitate a parametric study, the baseline case, identical to the one used in the second verification study above, is adopted here as a benchmark. Physical properties and operating conditions for the baseline case are tabulated in Table 2. Table 2. Physical Properties and Operating Conditions for the Baseline Case reactor type catalyst shape catalyst i.d. of single tube length of tube catalyst mean diam. catalyst density packed bed porosity outlet pressure inlet and coolant temp. syngas inlet mass flux inlet H2/CO molar ratio overall heat transfer coefficient

Figure 6. Syngas conversion comparison between current model prediction and results by Jess and Kern.23

dimensional model by Jess and Kern,23 carbon monoxide and hydrogen conversions are the same because the H2/CO feeding molar ratio and consumption ratio are identical, while the H2/CO consumption ratio is a variable in our comprehensive model, so that CO and H2 conversions are different in our comprehensive model. Figure 6 shows that, as expected, the conversion calculated by the two-dimensional pseudohomogeneous model is somewhere between the carbon monoxide conversion and hydrogen conversion predicted by our comprehensive model. It is noted that both models predicted the same trends that the conversions increases with the coolant temperature. In summary, we have developed hydrocarbon production rates for 15 species based on the general chain growth probability, stoichiometry, and intrinsic consumption rate for syngas reactants. Based on the above two independent comparisons, we believe that our model and the methodology are realistic and correctly implemented.

tubes-and-shell, focused on a single tube spherical cobalt based 4.6 cm 12 m 3 mm 1063 kg/m3 0.66 20 bar 214 °C 3.3 kg/m2 s 2 364 W/m2 K

Using the baseline case, the effects of varying the inlet H2/CO molar ratio, inlet and coolant temperatures, and inlet mass flux on the packed bed reactor performance are computed. Detailed simulation results for the baseline case are illustrated in Figures 7−11. Figure 7 shows the reactor bed centerline static pressure and temperature profiles along the axial direction. As shown in

5. RESULTS AND DISCUSSION Because the main objective of the current paper is to develop a numerical model that can predict the system performance characteristics for the Fischer−Tropsch synthesis reactor and its thermal management, we decided to provide the following parametric study for the purpose of illustrating some physical aspects of syngas conversion to hydrocarbons and thermal management criteria. The primary guiding principle that the thermal management is based on for an exothermic chemical reactor is to avoid the socalled “temperature runaway”. In a catalytic chemical reactor, the temperature runaway or thermal runaway is a process by which an exothermic chemical reaction goes out of control, often resulting in a temperature spike, deactivation of the catalyst, or catalyst sintering that usually causes the reactor to fail. Thermal runaway occurs when the reaction rate suddenly increases due to an increase in the reactor temperature, causing a corresponding increase in the heat production and, therefore, a further increase in the reactor temperature and a further increase in the reaction rate. This cyclic chain reaction is the cause of the temperature runaway. The phenomenon of temperature runaway has been predicted

Figure 7. (a) Pressure and (b) temperature profiles for the baseline case; pure syngas with mass flux of 3.3 kg/m2 s, H2/CO = 2, and inlet and coolant temperature at 214 °C. 1371


Energy & Fuels

Article

Figure 8. Temperature contours at three downstream locations for the baseline case; pure syngas with mass flux of 3.3 kg/m2 s, H2/CO = 2, and syngas inlet and coolant temperature at 214 °C.

For the temperature profiles, we selected three relatively upstream locations (z = 0.5, 1.5, and 2.5 m) because most heat transfer interactions take place there. For the upstream region, the flow is heated up by the heat released from the exothermic catalytic chemical reaction where the reactants concentrations are the highest. Some portion of the released heat is removed by the coolant flowing outside of the packed bed tubes, while the rest of the released heat facilitates the temperature increase of the reactor bulk flow. This causes the reactants to become more reactive, which accelerates the catalytic reaction process. Basically, the reactor bed temperature increases from 225 to 244 °C between z = 0.5 and 2.5 m in the axial direction due to the exothermic reaction. Also, the temperature gradients in the radial direction represent the driving force for the heat loss to the outside coolant through the wall and the heat loss rate increases as we proceed downstream as a result of the increase in the temperature difference between the reactor bed and the external coolant. At a certain point downstream, the flow temperature reaches the peak point and then starts decreasing as a result of increased convection heat loss to the outside coolant and lowered heat release rate from the chemical reaction because of the exhaustion of reactants. In Figure 9, the axial mass fraction profiles at the centerline for the gaseous phase are depicted in two different vertical scales (linear scale and log scale). As seen in Figure 9, the hydrogen and carbon monoxide mass fractions decrease linearly in the axial direction. The profiles of mass fractions for water and methane are also linear and so are the other hydrocarbons (C2−C6). It is also noted that, for this baseline case, the outlet mass fraction of gaseous phase (all gaseous species combined) is 0.8961, and the rest is mass fraction for liquid phase, which is 0.1039. The most abundant chemical species in the liquid phase, n-heptane, possesses a mass fraction of 0.1567

Figure 9. Mass fraction profiles at the centerline in the gaseous phase for the baseline case, (a) CO, H2, H2O, and CH4, and (b) all hydrocarbons in log scale; pure syngas mass flux of 3.3 kg/m2 s, H2/ CO = 2, inlet and coolant temperature 214 °C.

Figure 7, the pressure is linearly decreasing, and temperature is rapidly increasing in the beginning and asymptotically decreasing after the peak point. The two-dimensional temperature contours are depicted in Figure 8 for the axisymmetric FT chemical reactor, where the abscissa represents the axial coordinate and the ordinate represents the radial coordinate. 1372


Energy & Fuels

Article

in the liquid phase; therefore, the mass fraction of n-heptane in the total flow at the outlet is only 0.01629 (= 0.1567 × 0.1039).

More specifically, the carbon monoxide and water molar fraction contour plots are provided in Figures 10 and 11 for the baseline case. Three downstream locations of z = 2, 3, and 6 m

Figure 10. Contour plots for CO molar fractions at three downstream locations: z = 2, z = 3, and z = 6.

Figure 11. Contour plots for H2O molar fractions at three downstream locations: z = 2, z = 3, and z = 6. 1373


Energy & Fuels

Article

mass fluxes (F/Fbase = 0.5, 0.75, 1.0, 1.25, and 1.5; where Fbase denotes baseline case mass flux of 3.3 kg/m2 s) are evaluated, and their effects on the axial temperature profiles are given. Comparing with Figure 7b of the baseline case, which is the case of F/Fbase = 1, if the syngas mass flux is increased over that of the baseline case (cases F/Fbase = 1.25 and 1.5), the F-T catalytic reaction becomes much more intense, so that all the reactants are consumed in the first one-third of the reactor and the temperature runaway takes place, which is not a thermally viable case for the production of synthesis hydrocarbons. If the syngas inlet mass flux is decreased below that of the baseline case (cases F/Fbase = 0.5 and 0.75), the maximum reactor temperature is reduced, the position where the maximum temperature occurs is moved to a more upstream location, and the conversions of hydrogen and carbon monoxide to hydrocarbons are increased as a result of prolonged residence time. The carbon monoxide and hydrogen final percent conversions are tabulated in Table 3. It can be seen that the

are chosen. As shown in the figures, carbon monoxide molar fraction decreases along the axial direction, while water molar fraction increases with the axial direction. Unlike the temperature profiles in Figure 8, the two-dimensional molar contour profiles do not change the gradients significantly in the radial direction. The reason for this is due to the fact that heat is also removed radially by the external coolant but the mass (or chemical species) cannot be transported through the wall. 5.2. F-T Chemical Reactor Thermal Characteristics. Because the temperature distribution of the two-phase flow plays a crucial role in the performance of an exothermic Fischer−Tropsch catalytic reactor from the aspects of reactivity, selectivity, and stability, we have conducted a parametric study from the point of view of thermal management to assess the effects of syngas mass flow rate, syngas inlet and coolant temperatures, and H2/CO feed molar ratio on the F-T reactor internal temperature distributions. As background information, it should be pointed out that the temperature profile in the reactor bed mainly depends on the balance between the heat generated by the exothermic reaction and the heat removed by both the internal convection and the external coolant. The exothermic heat release is basically a function of the syngas inlet H2/CO molar ratio. The convective loss is primarily a function of the two-phase mass flow rate, and the heat loss to the external coolant is a function of the bed temperature, the coolant temperature, and the heat transfer coefficient between the tube outer surface and the coolant temperature. In the current analysis, the heat transfer coefficient is a fixed value of 364 W/m2 K; therefore, the heat loss to the external coolant is only varying with the coolant temperature, TC, which is equal to the syngas inlet temperature, TIn. As a result, in the following thermal management study, only syngas inlet H2/CO molar ratio, syngas gas inlet mass flux, and syngas inlet temperature (identical to the coolant temperature) will be varied. In Figure 12, the focus is on the effects of different syngas inlet mass fluxes while keeping all other system parameters equal to those of the baseline case. Five different syngas inlet

Table 3. Calculated Conversion Values for Selected Operating Conditions operating conditions

results

TIn/C [°C]

F/Fbase

H2/CO

XCO [%]

XH2 [%]

205

0.5

1.5 2.0 2.5 3.0 2.0 2.5 3 1.5 2.0 2.5 1.5 1.7 2.0 2.4 1.5 1.5 2.0 2.2 1.5 2.0 1.5 1.7 2.0

22.65 33.25 47.16 67.18 18.74 26.76 38.90 30.81 46.54 69.70 18.26 21.87 28.39 42.15 55.03 37.88 58.58 70.61 29.09 47.16 23.97 29.35 42.18

35.81 39.49 44.83 53.23 22.26 25.43 30.82 48.76 55.28 66.25 28.92 30.56 33.72 41.74 85.79 59.95 69.58 76.26 46.06 56.03 37.96 41.01 50.11

1.0

210

0.5

1.0

214

0.25 0.5

0.75 1.0

percent conversions for both hydrogen and carbon monoxide increase with increasing H2/CO inlet molar ratio and inlet syngas temperature but decrease with increasing syngas mass flux. Figures 13 and 14 show how the syngas mass flux affects the temperature profile for different hydrogen to carbon monoxide inlet molar ratios. The system conditions used in these figures are the same as those in Figure 12, except for the hydrogen to carbon monoxide inlet ratio (H2/CO = 1.5 in Figure 13 and H2/CO = 2.2 in Figure 14). When the syngas mass flux is increased, for all the hydrogen to carbon monoxide ratios, the peak temperature ascends, and its location is moved further downstream. So, the flow exit temperature keeps a continuously increasing trend with an increasing mass flux. A higher mass flux case also corresponds to a higher mass flow rate and a

Figure 12. Reactor bed temperature profiles for inlet and coolant temperature of 214 °C, H2/CO = 2.0, and different mass fluxes, F/Fbase = 0.5, 0.75, 1, 1.25, and 1.5. 1374


Energy & Fuels

Article

Figure 15. Reactor bed temperature profiles for inlet and coolant temperature of 214 °C, H2/CO = 2.5, and different mass fluxes, F/ Fbase = 0.25, 0.5, 0.75, and 1.

Figure 13. Reactor bed temperature profiles for inlet and coolant temperature of 214 °C, H2/CO = 1.5, and different mass fluxes, F/ Fbase = 0.25, 0.5, 0.75, and 1.

monoxide feed ratio increases to 2.2 while maintaining all other conditions the same, the cases F/Fbase = 0.75 and 1.0 result in the runaway of the reactor bed temperature that is considered thermally unviable. The thermally viable case, F/Fbase = 0.5, which yields 70.61% CO and 76.26% H2 conversions, respectively, reaches the maximum bed temperature of 244 °C and converts 17.2% input syngas into liquid products. Figure 15 shows the case of hydrogen to carbon monoxide feed molar ratio further increased to 2.5. In this case, even with the lowest mass flux case, F/Fbase = 0.25, the reactor becomes thermally unviable. For the lower mass flux case, the flow carries lower momentum when passing through the porous bed; so, if the heat released by the F-T reaction is accumulated rather than removed due to the lower flow rate, the reactor bed temperature will be increased, and this accelerates the exothermic reactions further, which results in the temperature runaway. In the actual experiment, this temperature runaway behavior may not be observed; instead, the deactivation of the catalyst would take place because the catalyst activity will be affected by the temperature, which is not considered in most simulations. The loss of catalytic activity due to high temperature is known as the deactivation of catalyst. The catalytic deactivation caused by high temperatures is also called catalyst sintering (also known as thermal degradation). The temperature runaway behavior in the simulation is still useful in providing a thermal management limiting condition for design consideration. If the heat released in a chemical reaction can be adequately removed then the temperature runaway can be avoided, so that a higher hydrogen to carbon monoxide feed ratio can be operated safely for a higher conversion. However, too much heat removal makes the reactor stay at relatively lower temperatures, in which the F-T reaction also cannot be activated. This is why the thermal management is important on an exothermic catalytic reaction. Enhancing the heat removal can be obtained by either increasing the heat transfer coefficient or increasing the temperature gradient using lower coolant temperatures. In this study, we only considered different coolant temperatures, but used a constant heat transfer coefficient.

Figure 14. Reactor bed temperature profiles for inlet and coolant temperature of 214 °C, H2/CO = 2.2, and different mass fluxes, F/ Fbase = 0.5, 0.75, and 1.

higher bulk velocity as we use a constant flow cross sectional area. This higher mass flux not only pushes the peak temperature further downstream but also makes the residence time shorter. Due to the shortened residence time, a higher mass flux case always results in a lower conversion of syngas despite a higher bed temperature (conversions are tabulated in Table 3). Although the syngas conversion is lower, the rate of total amount of syngas converted into hydrocarbon is higher for the higher mass flux case. For example, with the coolant temperature at 214 °C and H2/CO molar ratio at 1.5, let us compare two different syngas mass fluxes of F/Fbase = 0.25 and 1. On the basis of Table 3, the conversion for carbon monoxide, XCO = 23.97% for F/Fbase = 1 and XCO = 55.03% for F/Fbase = 0.25 gives the F/Fbase = 1 case a rate of total amount of syngas conversion that is 1.74 times higher than that of the F/Fbase = 0.25 case. As shown in Figure 14, when the hydrogen to carbon 1375


Energy & Fuels

Article

Figures 16−18 show how syngas inlet and coolant temperatures affect the reactor temperature profile as well as the conversions. In Figure 16, reactor temperature profiles for various H2/CO

Figure 18. Reactor bed temperature profiles for inlet and coolant temperature of 205 °C, syngas mass flux F = Fbase, and different H2/ CO ratios of 2.0, 2.5, 3.0, and 3.5.

directly proportional to the hydrogen molar fraction, but the carbon monoxide dependency is more complex than the hydrogen. Carbon monoxide acts as an inhibitor when its concentration is relatively high. When the carbon monoxide adsorption term in the denominator is greater than unity, 1 ≪ KCCO, then the entire denominator term could be approximated as (KCCO)2. In this case, F-T synthesis intrinsic kinetics is inversely proportional to the carbon monoxide concentration, which is a typical characteristic of the Langmuir−Hinshelwood kinetics. As we increase the hydrogen to carbon monoxide molar ratio at the inlet, the syngas consumption rate would be accelerated as a result of the increased hydrogen concentration so that more heat will be released, the reactor temperature will be higher and finally the exit conversion will be raised. By reducing the coolant temperature to 210 °C, the catalytic packed bed becomes thermally viable up to H2/CO = 2.4, which is shown in Figure 17. However, as mentioned previously, the bed temperature remains lower resulting in low syngas conversions. We found that the performance for hydrogen to carbon monoxide feed ratio of 2.0 with a coolant temperature of 214 °C is similar to that of hydrogen to carbon monoxide feed ratio of 2.4 with 210 °C coolant temperature. For the coolant temperature of 210 °C and H2/CO of 2.4, the CO conversion is 42.15%, the H2 conversion is 41.74%, and the mass converted into liquid phase is 10.13%. If the coolant temperature drops further, even higher H2/CO ratios can be thermally viable. The results for coolant temperature of 205 °C are illustrated in Figure 18, where the case of a H2/CO ratio as high as 3 is still thermally viable. In the case of H2/CO = 3.3, its temperature suddenly rises at almost half of the reactor length. Next, results are obtained using the same coolant temperatures as those in Figures 16−18, except the syngas inlet mass flux value is reduced by half. Temperature profiles depicted in Figure 19 are obtained under the same conditions as those in Figure 16, except the syngas mass flux. The temperature profiles illustrated in Figure 19 are very similar to those in Figure 16, except that the maximum H2/CO ratio for a thermally viable application is 2.2 instead of 2, as in Figure 16. Also, the downstream location where the peak temperature

Figure 16. Reactor bed temperature profiles for inlet and coolant temperature of 214 °C, syngas mass flux F = Fbase, and different H2/ CO ratios of 1.5, 1.7, 2.0, 2.2, 2.5, and 3.0.

Figure 17. Reactor bed temperature profiles for inlet and coolant temperature of 210 °C, syngas mass flux F = Fbase, and different H2/ CO ratios of 1.5, 1.7, 2.0, 2.4, 2.5, and 3.0.

ratios with F/Fbase = 1.0 and a syngas inlet and coolant temperature of 214 °C are depicted. Among these, the baseline case (H2/CO = 2) has the best performance among the thermally viable cases, and the higher hydrogen to carbon monoxide feed ratios become thermally unviable. It is not shown here, but the case of H2/CO = 2.1 also yields the temperature runaway behavior. The higher the H2/CO ratio gives the faster temperature rise. The temperature dependency on the H2/CO ratio can be explained by its intrinsic kinetics behavior, with respect to hydrogen and carbon monoxide concentrations. The intrinsic kinetics of F-T synthesis, eq 10, is 1376


Energy & Fuels

Article

Figure 19. Reactor bed temperature profiles for inlet and coolant temperature of 214 °C, syngas mass flux F = 0.5Fbase, and different H2/ CO ratios of 1.5, 2.0, 2.2, 2.3, 2.5, and 3.0.

Figure 21. Reactor bed temperature profiles for inlet and coolant temperature of 205 °C, syngas mass flux F = 0.5Fbase, and different H2/ CO ratios of 1.5, 2.0, 2.5, and 3.0.

occurs is closer to the inlet for lower mass flux case. As discussed previously, syngas conversion is higher for lower mass flux as a result of relatively longer residence time. For example, CO conversions are 42.18% for H2/CO = 2.0 in Figure 16 and 70.61% for H2/CO = 2.2 in Figure 19 as provided in Table 3. The half syngas mass flux version of Figure 17 is depicted on Figure 20. Findings similar to those in the coolant temperature

reactor design. We have prepared a thermal viability map given in Figure 22 to summarize the thermal management strategy.

Figure 22. Thermal viability map for a F-T reactor.

The thermal viability map is developed using the three integral parameters: syngas inlet mass flux, F; syngas inlet and coolant temperature, TIn = TC; and hydrogen to carbon monoxide feed molar ratio, θH2/CO. For each data point with the specific set of F and TC, it represents the maximum θH2/CO value that the F-T reactor is thermally viable (no reactor bed temperature runaway). For example, for the point with F/ Fbase = 1 and TC = 214 °C, the maximum θH2/CO without a reactor bed temperature runaway is 2. In other words, any point in the area under a specific curve represents a thermally viable case for the coolant temperature corresponding to that curve. Therefore, each curve can be considered as the critical threshold boundary for thermal viability. In general, the critical

Figure 20. Reactor bed temperature profiles for inlet and coolant temperature of 210 °C, syngas mass flux F = 0.5Fbase, and different H2/ CO ratios of 1.5, 2.0, 2.5, and 3.0.

of 214 °C case can be seen. Temperature profiles for the 205 °C coolant temperature case with the half mass flux are plotted in Figure 21. No cases are found to be thermally unviable for H2/CO ratios in the range from 1.5 to 3.0. 5.3. Thermal Management Analysis. As discussed above, it is apparent that the temperature profiles in the reactor bed and its thermal management hold the key for the optimal F-T 1377


Energy & Fuels

■

θH2/CO value increases with decreasing coolant temperature and decreasing syngas mass flux, F. It is worth mentioning that, for the case of TIn = TC = 205 °C and the lowest syngas mass flux, F/Fbase = 0.25, the reactor could function without a thermal runaway for θH2/CO as high as 3.9. For this particular case with a relatively high critical θH2/CO value, the limiting chemical species, carbon monoxide, is totally consumed within one-third of the reactor length from the entrance and the bed temperature has reached its maximum point at this location. After the depletion of the limiting chemical species, the bed temperature will stop rising and stabilize as a result of no more heat release. However, this point may not be the ideal operating condition because the rest of the reactor is not utilized. This map just provides the thermal viability; however, for the ideal operating condition, the thermal viability should be considered together with reactant conversion and product selectivity.

Article

AUTHOR INFORMATION

Corresponding Author

*Telephone: 1 352-392-9607. Fax: 1 352-392-1071. E-mail: [email protected].

■

ACKNOWLEDGMENTS This research was supported by the Florida Department of Agriculture and Consumer Services (FDACS). The partial support by the Andrew H. Hines Jr./Progress Energy Endowment Fund is also acknowledged. T.S.L. was also partially supported by the Korea Science and Engineering Foundation Grant funded by the Korea government (Ministry of Science and Technology) (2005-215-D00037).

■

6. CONCLUSION An axisymmetric two-dimensional numerical model that is based on two-phase gas−liquid flow in a heterogeneous porous medium has been developed to simulate a fixed packed bed tubular Fischer−Tropsch reactor. The detailed chemical kinetics for producing linear-paraffin in both gaseous and liquid phases has been derived based on carbide mechanism chain growth probability and stoichiometry for a non-ASF distribution. The thermal-fluid transport is modeled as a twophase liquid trickling flow going through a packed bed of porous material consisting of solid catalyst particles. An Eulerian−Eulerian mixture model has been used for the twophase flow simulation, together with the porous material assumed as momentum sinks in the fixed bed reactor. Mass transfer through the catalyst pellet pores is also considered by means of the general Thiele modulus. The FLUENT computer code that accommodates the F-T reactor model described was employed to provide the numerical solutions. Two comparisons have been made to validate our model. First, the products distribution predicted for a non-ASF distribution gives a satisfactory agreement between the current model predictions and the experiment results. The second comparison with a simplified model on the packed bed temperature variations and thermal management not only validated the current model but also proved that a comprehensive model is more useful and important when assessing the thermal viability of the reactor. The thermal characteristics of a F-T chemical reactor have been investigated with a focus on the effects of syngas mass flux, syngas inlet and coolant temperatures, and H2/CO feed molar ratio on the reactor temperature profiles. The simulation results indicate the following findings: (1) An increased syngas mass flow rate results in a shorter residence time, which causes a lower conversion, a higher peak temperature, and the location of peak temperature to move more downstream. (2) A higher hydrogen to carbon monoxide feed molar ratio makes the higher syngas flow rates thermally unstable. (3) A lower syngas inlet/coolant temperature will quench the reactions so that a higher mass flow rate or a higher H2/CO molar ratio can still be thermally viable. Among the mass flux range from 0.825 to 4.95 kg/m2 s, inlet/coolant temperature range from 205 to 214 °C, and H2/CO feed ratio range from 1.5 to 3.0, we have found that the case of F = 1.65 kg/m2 s, θH2/CO = 2.2, and TC = 214 °C is the optimum operating condition that gives the highest syngas conversion with no reactor bed temperature runaway.

NOMENCLATURE C = concentration (mol/m3). D = diameter of the reactor (m). Deff = effective diffusion coefficient (m2/ s). Dp = catalyst particle (or pellet) average diameter (μm). F = inlet mass flux (kg/(m2 s)). H = Henry coefficient (Pa m3/mol). h = enthalpy (J/kg). J = molar flux (mol/m2 s). K = CO adsorption constant (m3/mol). k = reaction rate constant (m6/(kgcat mol s)). m = molar fraction. Lc = characterization length (m). MW = molecular weight (g/mol). N = highest carbon number species assumed in the product. Ṅ = molar flow rate (mol/s). p = pressure (Pa). R = ideal gas constant ( J/mol K). r = catalytic chemical reaction rate (mol/kgcat s). SM = momentum sink term in momentum conservation equations (kg/m2 s2). SE = volumetric heat source term in energy conservation equation (W/m3). T = absolute temperature (K). t = time (s). v⃗ = velocity vector (m/s). X = conversion (% or nondimensional). Y = mass fraction.

Greek

α̅ = carbon number independent constant chain growth probability. αi = chain growth probability for hydrocarbon with a carbon number i. χ = volume fraction. δ = the inertial resistance factor. ε = porosity of the packed bed. ϕ = general Thiele modulus. γ = stoichiometric coefficient (positive for products and negative for reactant). η = effectiveness factor. κ = permeability (m2). λ = thermal conductivity (W/mK). μ = viscosity (kg/m2 s). θ = molar ratio. ρ = density (kg/m3). τ = shear stress tensor (kg/m2 s2).

Other

9 = heterogeneous reaction rate. 1378


Energy & Fuels

Article

Subscript

(29) Ribeiro, F. H.; von Wittenau, A. E. S.; Bartholemew, C. H.; Somorjai, G. A. Catal. Rev.-Sci. Eng. 1997, 39, 49−76. (30) Maretto, C.; Krishna, R. Catal. Today 1999, 52, 279−289. (31) Hamelinck, C. N.; Faaij, A. P. C.; den Uil, H.; Boerrigter, H. Energy 2004, 29, 1743−1771. (32) Philippe, R.; Lacroix, M.; Dreibine, L.; Pham-Huu, C.; Edouard, D.; Savin, S.; Luck, F.; Schweich, D. Catal. Today 2009, 147, S305−S312. (33) ANSYS FLUENT, version 12 software package; ANSYS Inc.: Canonsburg, PA. Available online: http://www.ansys.com. (34) ANSYS FLUENT UDF Manual Version 12.0. ANSYS Inc.: Canonsburg, PA, 2009. (35) Jovanovic, R.; Milewska, A.; Swiatkowski, B.; Goanta, A.; Spliethoff, H. Int. J. Heat Mass Transfer 2011, 54, 921−931. (36) Jin, Y.; Shaw, B. D. Int. J. Heat Mass Transfer 2010, 53, 5782− 5791. (37) Chein, R. Y.; Chen, Y. C.; Chang, C. S.; Chung, J. N. Int. J. Hydrogen Energy 2010, 35, 589−597. (38) Ho, C. D.; Chang, H.; Chen, H. J.; Chang, C. L.; Li, H. H.; Chang, Y. Y. 2011, 54, 3740-3748. (39) Patankar, S. V. Numerical Heat Transfer and Fluid Flow, Hemisphere: Washington, DC, 1980. (40) Versteeg, H. K.; Malalasekera, W. An Introduction to Computational Fluid Dynamics, The Finite Volume; Method Longman: Essex, England, 1995. (41) Ferziger, J. H.; Perić, M. Computational Methods for Fluid Dynamics, 3rd ed., Springer-Verlag: Berlin Heidelberg, 2002. (42) NIST Chemistry WebBook. http://webbook.nist.gov/chemistry/ (43) Elbashir, N. O.; Roberts, C. B. Ind. Eng. Chem. Res. 2005, 44, 505−521. (44) Lynch, D. T.; Wanke, S. E. Can. J. Chem. Eng 1991, 69, 332− 339. (45) Berty, J. M. Experiments in Catalytic Reactor Engineering. In Studies in Surface Science and Catalysis; Elsevier: New York, 1999; Vol. 124, pp 157−159.

base = baseline case. C = coolant. G = gas phase. H2/CO = hydrogen to carbon monoxide. In = inlet. ini = chain growth initiation. L = liquid phase. m = mass average. n = carbon number.

■

REFERENCES

(1) Huber, G. W. Breaking the Chemical and Engineering Barriers to Lignocellulosic Biofuels: A Research Roadmap for Making Lignocellulosic Biofuels a Practical Reality. NSF, DOE, and American Chemical Society Workshop, Washington, DC, June 25−26, 2007. (2) Greene, N. Growing Energy: How Biofuels Can Help End America’s Oil Dependence; Natural Resources Defense Council: New York, 2004; available online: http://www.nrdc.org/air/energy/biofuels/biofuels. pdf. (3) Sheehan, J.; Camobreco, V.; Duffield, J.; Graboski, M.; Shapouri, H. Life Cycle Inventory of Biodiesel and Petroleum Diesel for Use in an Urban Bus, NREL/SR-580-24089 UC Category 1503; Department of Energy, National Renewable Energy Laboratory: Golden, CO1998. (4) Steynberg, A. P.; Dry, M. E. Fischer−Tropsch Technology. Studies in Surface Science and Catalysis; Elsevier: New York, 2004. (5) Sousa, E. F.; Noronhac, F. B.; Faro, A. Catal. Sci. Technol. 2011, 1, 698−713. (6) Damartzis, T.; Zabaniotou, A. Renewable Sustainable Energy Rev. 2011, 15, 366−378. (7) Huff, G. A.; Satterfield, C. N. Ind. Eng. Chem. Process Des. Dev. 1984, 23, 696−705. (8) Sarup, B.; Wojciechowski, B. W. Can. J. Chem. Eng. 1989, 67, 62−74. (9) Kellner, C. S.; Bell, A. T. J. Catal. 1981, 70, 418−432. (10) Krishna, K. R.; Bell, A. T. J. Catal. 1993, 139, 104−118. (11) Lox, E. S.; Froment, G. F. Ind. Eng. Chem. Res. 1993, 32, 71−82. (12) Adesina, A. A. Appl. Catal., A 1996, 138, 345−367. (13) Teng, B. T.; Chang, J.; Zhang, C. H.; Cao, D. B.; Yang, J.; Liu, Y.; Guo, X. H.; Xian, H. W.; Li, Y. W. Appl. Catal., A 2006, 301, 39−50. (14) Dry, M. E. Appl. Catal., A 1996, 138, 319−344. (15) Deckwer, D. W. Oil Gas J. 1980, 78, 198−213. (16) Troshko, A. A.; Zdravistch, F. Chem. Eng. Sci. 2009, 64, 892− 903. (17) Atwood, H. E.; Bennett, C. O. Ind. Eng. Chem. Process Des. Dev. 1979, 18, 163−170. (18) Bub, G.; Baerns, M.; Bussemeier, B.; Frohning, C. Chem. Eng. Sci. 1980, 35, 348−355. (19) Jess, A.; Popp, R.; Hedden, K. Appl. Catal., A 1999, 86, 321− 342. (20) Wang, Y. N.; Xu, Y. Y.; Li, Y. W.; Zhao, Y. L.; Zhang, B. J. Chem. Eng. Sci. 2003, 58, 867−875. (21) De Swart, J. W. A.; Krishna, R.; Sie, S. T.Studies in Surface and Catalysis; Elsevier: New York, 1997; Vol. 107, pp 213−218. (22) Moutsoglou, A.; Sunkara, P. P. Energy Fuels 2011, 25, 2242− 2257. (23) Jess, A.; Kern, C. Chem. Eng. Technol. 2009, 32, 1164−1175. (24) Claeys, M.; van Steen, E. Basic Studies. In Studies in Surface and Catalysis; Elsevier: New York, 2004; Vol. 152, Chapter 8. (25) Yates, I. C.; Satterfield, C. N. Energy Fuels 1991, 5, 168−173. (26) Steynberg, A. P.; Dry, M. E.; Davis, B. H.; Breman, B. B. Fischer−Tropsch Reactors. In Studies in Surface and Catalysis; Elsevier: New York, 2004; Vol. 152, Chapter 2. (27) Ergun, S. Chem. Eng. Prog. 1952, 48, 89−94. (28) Iglesia, E.; Soled, S. L.; Fiato, R. A. J. Catal. 1992, 137, 212−224. 1379


Mathematical Modeling and Numerical Simulation of a Fischer

Recommend Documents