CFD analysis of hot spot formation through a fixed bed reactor of Fischer-Tropsch synthesis

One of the interesting methods for conversion of synthesis gas to heavy hydrocarbons is Fischer–Tropsch process. The process has some bottlenecks, such as hot spot formation and low degree of conversion. In this work, computational fluid dynamics technique was used to simulate conversion of synthetic gas and product distribution. Also, hot spot formation in the catalytic fixed-bed reactor was investigated in several runs. Simulation results indicated that hot spot formation occurred more likely in the early and middle part of reactor due to high reaction rates. Based on the simulation results, the temperature of hot spots increased with increase in the inlet temperature as well as pressure. Among the many CFD runs conducted, it is found that the optimal temperature and pressure for Fischer– Tropsch synthesis are 565 K and 20 bar, respectively. As it seems that the reactor shall work very well under optimal conditions, the reaction rates and catalyst duration would simultaneously be maximum . Subjects: Chemical Processing & Design; Process Control Chemical Engineering; Reaction Engineering


Introduction
Nowadays, the energy crisis is a real challenge to human society, since oil sources in the world are declining. A promising alternative approach in this regard is to directly convert natural gas to liquid fuel by means of gas-to-liquid (GTL) technology . Fischer-Tropsch synthesis (FTS), with a complex network of reactions in parallel and/or in series, converts synthesis gas into valuable liquid hydrocarbons, clean fuels, and chemical feedstock (Li, Zhang, Zhang, & Cao, 2007;Ma, Sun, Ying, & Fang, 2009;Pour, Housaindokht, Tayyari, & Zarkesh, 2010a;Pour, Zare, & Zamani, 2010b;Takassi, Salooki, & Esfandyari, 2011). In this process, a wide range of paraffins and olefins (from C 1 to C 7 + ) can be synthesized, depending on the operating conditions, catalyst characteristics, and type of the process (Feyzi & Jafari, 2012). Because of escalation in the price of crude oil, air pollution, global warming, and resulting climate change, the process is of great interest as a vital alternative for producing clean fuels (Dai & Yu, 2008;Forghani, Elekaei, & Rahimpour, 2009;Hao, Bai, Xiang, & Li, 2009;James, Mesubi, Ako, & Maity, 2010;. The fixed-bed Fischer-Tropsch process plays a unique role in FTS industrial units, so that a variety of hydrocarbons ranging from methane to wax are produced by this process (Wang, Xu, Li, Zhao, & Zhang, 2003). Figure 1 shows a typical FTS plant equipped with fixed-bed reactor. The kinetic scheme used in this paper has some novelties; firstly, the scheme is not a kinetic totally rendered earlier. The kinetic relations are selective from different models. Secondly, we select this kinetic scheme according to our pilot plant and experimental data which is published earlier (Aligolzadeh, 2011). In this paper, we focus on the production of C 7 + hydrocarbons from reactants, which are plentifully available and compatible with national resources in Iran (natural gas). Hence, the process and operation conditions are set, so that the waxes formation is low.
A heterogeneous fixed-bed reactor was modeled/simulated, so as to evaluate the performance of the reactor, in which study, there was good agreement between simulation and experimental data . Later on, an experimental study on the milli-structured fixed-bed reactor for FTS process was performed, along with mathematical modeling (Knochen, Güttel, Knobloch, & Turek, 2010). Rafiq, Jakobsen, Schmid, and Hustad (2011) carried out experimental analysis as well as mathematical modeling of a fixed-bed reactor for FTS process fed with biosynthesis gas. A two-dimensional pseudohomogeneous model was developed by Rafiq et al. to evaluate the performance of the FTS reactor, the results of which showed that the model predicts experimental data successfully. Computational fluid dynamics (CFD) technique as a powerful device was applied for three-dimensional simulation of hydrogenation of carbon monoxide in a fixed-bed FTS reactor (Miroliaei, Shahraki, Atashi, & Karimzadeh, 2012). They discovered that operational conditions extremely affect the product distribution; and, the more temperature increases, the more conversion increases. Moreover, it was found that the selectivity parameter and mass fraction of the products increase in value, as a result of temperature or the ratio of H 2 /CO increase.

Figure 1. A schematic diagram of fixed-bed type for FTS process.
In this work, two-dimensional simulation of GTL fixed-bed reactor is carried out using CFD technique. The commercial CFD solver, FLUENT.6.3, is used for numerically solving the governing equations. Effects of operating pressure and feed temperature on the hot spot formation are investigated. A detailed kinetic scheme is used to estimate the rate of materials generated/consumed during the process once started. Mass fraction of present components, evaluation of heat of reactions, and analysis of hot spot formation and its location, all ascertain the necessity of using a detailed kinetic scheme. In fact, the novelty of present work is limited for incorporating a detailed kinetic scheme into the CFD solver and looking for hot spot formation.

Kinetic scheme
In this work, we use the kinetic which best fitted the experimental data at the outlet of GC. Kinetic scheme of Yang et al. is used to predict the rate of reactions and added to the CFD solver (Yang et al., 2003): We can consider four steps for these reactions, the following reactions are considered as leading FTS reactions (Yang et al., 2003).

First step is adsorption of reactants:
Second step is chain initiation. We assumed that chain growth is initiated by the adsorption of CO on the active sites (Yang et al., 2003).
Third step is chain propagation (Yang et al., 2003): Fourth step is chain termination and desorption of products, termination of the chain growth can occur via several routes for FTS (Yang et al., 2003): where s 1 is a vacant catalytic site, on which hydrocarbon can be formed. The reaction rate equations for carbon dioxide, methane, paraffins, and olefins are as follows, and the kinetic parameters are (1)

Paraffin formation:
nCO Olefin formation: Water-gas shift (WGS): Hs 1 + CO ⇔ Hs 1 CO Hs 1 CO + H 2 ⇔ Hs 1 C + H 2 O C n H 2n+1 s 1 + CO ⇒ C n H 2n+1 s 1 CO C n H 2n+1 s 1 CO + H 2 ⇔ C n H 2n+1 s 1 C + H 2 O C n H 2n+1 s 1 C + H 2 ⇔ C n H 2n+1 s 1 CH 2 (10) CH 3 s 1 + Hs 1 ⇒ CH 4 + 2s 1 C n H 2n+1 s 1 + H 2 ⇒ C n H 2n+2 + Hs 1 C n H 2n+1 s 1 ⇔ C n H 2n + Hs 1 listed in Table 1. According to water-gas shift reaction, the rate equation for CO 2 formation is as follows: where The reaction rate for methane, paraffins, and olefins formation are given: where The rate constant in the above equations can be defined according to Arrhenius relation, which is consistent with what CFD solver needs: where n is temperature exponent (dimensionless). The aforementioned kinetics has been developed, when using industrial manganese-promoted iron catalyst, including 0.5 wt% Mn and 1.5 wt% Fe. All reactions are defined volumetric type with generalized laminar finite-rate formulation.

Governing equations
There are a variety of commercial CFD solvers coded based on the finite element, finite difference, and finite volume methods. In the present study, finite volume method (FLUENT 6.3) was utilized for solving the governing equations.

Momentum conservation equation
The conservation of momentum is expressed by: The equation is transformed to so-called Brinkman equation, which is valid for porous media and cylindrical coordinates, similar to the conditions of FTS catalytic reactor.

Mass conservation equation
General equation for conservation of mass, continuity equation, is defined as follows: This equation may be simplified if we make the assumption that there is no variation of fluid density with time (steady-state condition), so the instability term (partial derivative of density with respect to time) in the left-hand side of Equation 25 becomes zero. The right-hand side term is likewise zero.

Energy conservation equation
Energy conservation equation valid for porous media including conduction flux terms is given by: The term on the left-hand side of Equation 27 is the heat transfer due to convection. On the righthand side, the terms are conduction heat transfer, species enthalpy due to diffusion flux of, for example, species j, heat dissipation, and heat source, respectively. Through the porous medium, conduction flux is calculated using an effective conductivity, as Equation 28 shows. The effective thermal conductivity in the porous medium, k eff , is calculated as the volume-weighted average of the fluid and solid conductivity (Shahhosseini, Alinia, & Irani, 2009):

Species transport equation
The equation governing species transport through the reactor is of course well known: Assuming steady-state and irrotational flow (θ direction), Equation 29 diminishes such: The intrinsic kinetics scheme is independent of boundary layer effects. But in the species transport equation, we have the terms, including the local velocity, which affects on the mass transfer and conversion, and the component concentration. When the species concentration changes with velocity profile, thus the rate of all reactions alter according to local concentration, as the reaction rate equations dictate. This is why we should consider boundary layer effects on the calculations.

Mesh generation and grid dependency
A two-dimensional axisymmetric system (0.046 × 7 m) is considered as the domain of fluid, through which the governing equations should numerically be solved. Two different grids for the mentioned geometry are generated using Gambit software 2.3.2 edition; they have a total number of 84330 and 112440 structured quadrilateral elements, respectively. Primary simulation is carried out by commercial CFD solver, FLUENT (Fluent Inc., 2003) 6.3 edition, for checking grid dependency. Mass fraction of CO versus reactor length for both grids is shown in Figure 2.
It is clear that the results obtained for both grids are coincident; therefore, the results of coarser grid are as good as finer grid. Consequently, choosing the coarser grid, rather than finer one, would enable solving the problem more quickly. Even the coarser grid is fine enough to capture near the wall treatment (wall effects). The height of first cells adjacent to the wall boundary is fine enough and the height of next cells grows up with a density ratio of 1.2 to capture near the wall treatment. The mesh is designed to be denser near the inlet and outlet boundaries, as one can see in Figure 3.

Boundary conditions
The geometry discussed above is subjected to the following boundary conditions: (1) The flow velocity was specified as inlet boundary condition (0.55 m/s).
(2) Outflow boundary was considered as outlet boundary condition.
(3) Wall boundary condition was set for other boundaries with zero diffusive flux for all species.
The solver is selected to be pressure-based because of solving incompressible flow. Among several pressure-velocity-coupling algorithms, the well-known SIMPLE algorithm is utilized. Moreover, the second-order upwind scheme is employed to discrete the governing equations due to its high accuracy. The turbulence modeling is passed over throughout the domain because of laminar flow.

Operating condition
The temperature at the inlet of the reactor is 553 K and the system works under different pressures (see Table 4). The characteristics of the reactor and catalyst are listed in Table 2.

Figure 3. Geometry and grid with 84430 structured quadrilateral elements (the lower boundary is axis of symmetry).
Bed characteristics are estimated through rigorous porous media relations having basic information about size and shape of catalyst. Particle solid is considered spherically 3 mm in diameter and the porosity of the bed including particles is assumed as 0.3894. The term F in Equation 24 symbolizes the momentum source term induced by body forces or here by porous media which is identical to: Thefirst term on the right-hand side is engaged in F as viscous loss and the second as inertial loss, where F i refers to source of ith (x, r or z) momentum equation. The notation |V| shows velocity magnitude everywhere, and D and C are prescribed matrices being able to switch to scalar in case of isotropic environment. For homogeneous porous media, above equation is changed to: where α is permeability and C 2 is the inertial resistance factor and shows itself as diagonal elements of D ij matrix, where the rest of elements are all zero. In fact, porous media theory is nothing more than a momentum sink on the right-hand side of momentum conservation equations. Other important subject in the porous media formulation is the pressure drop modeling; somehow, there are a few theoretical, empirical, or semi-empirical models for evaluating the pressure drop across porous media. The most used model is the Ergun's equation, a semi-empirical equation, which correlates the pressure drop with velocity magnitude, as demonstrated beneath (Ergun, 1952).

Or in simpler form:
Permeability and inertial loss coefficient may theoretically be determined by physical properties of the catalyst and bed porosity. The relevant relations for the aforesaid quantities are: Effect of sphericity on the porous pressure drop is not included in Equation 35 holding it to elsewhere (Equation 33); hence, both permeability and inertial loss coefficient are shape independent. Viscous and inertial resistances for the bed are calculated with the help of equations. We get the values of 1.052385e + 08 1/m 2 and 12064 1/m, respectively, for viscous and inertial resistances. It must be said that viscous resistance is inversely associated with permeability. (31) (1 − ) 2

Number of tubes 1
Reactor internal diameter (m) 0.046 Bed density (kg/m 3 ) 1100 Diameter of catalyst particles (mm) 3 Bed porosity 0.3894 In the present work, all physical properties of feed are assumed to be constant and temperature independent, except for density and viscosity that they are computed by (incompressible flow of an ideal gas). 1 Mass diffusivity is also calculated by kinetic theory (FLUENT6.3., XXXX).
The species mass fraction at the inlet is considered, as one can see in Table 3. The system is insulated and the heat flux from the wall to environment is zero.

Results and discussion
In this section, the results of CFD simulation are compared with experimental data reported by Yang et al., as one can see in Table 4, showing that the simulation results are in good agreement with experimental data (Yang et al., 2003). The last column shows that the maximum relative error obtained is below 6%.
Temperature profile as a function of reactor length under different total pressures is displayed in Figure 4. According to this figure, the temperatures of hot spots and products at the outlet rise as pressure increases. It is observed that the hot spot location occurs almost near the reactor inlet due to high reactants concentration there. The locus of hot spot formation shifts gently to the inside, as a result of decreasing pressure. This conclusion is completely according to what Chapman-Enskog equation represents; penetration rate is directly proportional to the temperature, penetration, and molar flux increase caused by rising temperature, leading to shift in hot spot location to the inside of the reactor (Miroliaei et al., 2012). Therefore, considering the key role of the temperature on the product selectivity, it can be expected that rising reactor bulk temperature caused by bulk pressure increase will lead to the production of undesirable heavy linear hydrocarbons.
In the following, four different temperatures were considered to study the effect of inlet temperature on the hot spot formation. The characteristics of flow inlet was reported in literature (Yang et al., 2003). Figure 5 shows the effect of inlet temperature on the hot spot formation.
As obvious, the temperature related to hot spot formation increases as inlet temperature rises. It is ascribed to increase in reactants enthalpy in addition to raising reaction rate. Contours of mole  Figure 6, showing that mole fraction of CO and H 2 decreases along the length of the reactor. Increasing pressure and temperature affects to increasing heavy hydrocarbons. We can rationalize this idea by the activation energy of Olefin formation, which is 97.37 kJ mol −1 , much smaller than that of paraffin's (Yang et al., 2003), which can interpret the general fact that rising reactor bulk temperature caused by bulk pressure will increase the production of heavy linear hydrocarbons on the Fe-Mn catalysts than the other iron-based catalysts.
Predicted concentration of products through the axis of fixed-bed reactor is shown in Figures 7-12. All contours are extracted from FLUENT6.3 and then depicted by MATLAB 6.5Ed software. In all figures, mole fractions change was depicted on vertical axis and reactor length on horizontal axis. Figures 7-12 indicate that the concentration of products increases by reaction progress along the length of the reactor. It can be attributed to the inhomogeneity in the flow around a catalytic particle, which strongly affects the concentration distribution in the reaction zone (Miroliaei et al., 2012).
Although temperature rise has a positive effect on the reactor performance, further increase will cause catalyst to damage because of hot spot formation. Using CFD simulation, optimal operational temperature and pressure for fixed-bed reactor of FTS were obtained as 565 K and 20 bar, respectively. The optimal conditions were found through several CFD runs; however, the reactor works under optimal conditions as long as the pressure and temperature are 565 K and 20 bar. Simulation results demonstrated that C 7 + productivity and the rate of production of linear hydrocarbons increase as temperature and pressure rise. In addition, more paraffins with equal carbon chains are produced in comparison with olefins. It is ascribed to higher reaction rate of paraffins compared to olefins.

Conclusion
Two-dimensional CFD technique was employed to predict the conversion of synthetic gas (CO and H 2 ) and distribution of products through fixed-bed reactor in the FTS process. Results show that conversion of CO and H 2 and the selectivity of C 7 + as it seen on Figure 13 increase along the length of reactor. As a consequence of high reaction rates in the early part of the reactor that decrease along the length of the reactor, hot spot formation occurred more likely in the early part and less in the middle part of the reactor. Results show that any increase in the temperature leads to higher reaction rate, but it causes deactivation of catalyst due to hot spot formation. For this reason, there must be some optimal conditions for FTS process, under which the rate of desirable reactions and catalyst duration shall simultaneously be maximum. Among the many CFD runs done, optimal temperature and pressure for fixed-bed reactor of FTS process were acquired as 565 K and 20 bar, respectively.