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

Abstract

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 .

Keywords:

• CFD simulation
• fixed-bed reactor
• GTL technology
• hot spot formation
• Fischer–Tropsch synthesis
• optimal conditions

Public Interest Statement

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.

1. 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 (Li et al., 2009). 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₁ to C₇⁺) 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; Wan et al., 2007). 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₇⁺ 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.

Figure 1. A schematic diagram of fixed-bed type for FTS process.

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 (Wang et al., 2003). 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 pseudo-homogeneous 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₂/CO increase.

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.

2. 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):(1) $$\text{Paraffin formation}:n\text{CO}+(2n+1){\text{H}}_2\iff {\text{C}}_n{\text{H}}_{2n+2}+n{\text{H}}_2\text{O}$$(1) (2) $$\text{Olefin formation}:n\text{CO}+2n{\text{H}}_2\iff {\text{C}}_n{\text{H}}_{2n}+n{\text{H}}_2\text{O}$$(2) (3) $$\text{Water-gas shift}(\text{WGS}):\text{CO}+{\text{H}}_2\text{O}\iff {\text{CO}}_2+{\text{H}}_2$$(3)

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:(4) $$2{\text{s}}_1+{\text{H}}_2\iff 2{\text{Hs}}_1$$(4) (5) $${\text{Hs}}_1+\text{CO}\iff {\text{Hs}}_1\text{CO}$$(5)

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).(6) $${\text{Hs}}_1\text{CO}+{\text{H}}_2\iff {\text{Hs}}_1\text{C}+{\text{H}}_2\text{O}$$(6) (7) $${\text{C}}_n{\text{H}}_{2n+1}{\text{s}}_1+\text{CO}\Rightarrow {\text{C}}_n{\text{H}}_{2n+1}{\text{s}}_1\text{CO}$$(7) (8) $${\text{C}}_n{\text{H}}_{2n+1}{\text{s}}_1\text{CO}+{\text{H}}_2\iff {\text{C}}_n{\text{H}}_{2n+1}{\text{s}}_1\text{C}+{\text{H}}_2\text{O}$$(8) (9) $${\text{C}}_n{\text{H}}_{2n+1}{\text{s}}_1\text{C}+{\text{H}}_2\iff {\text{C}}_n{\text{H}}_{2n+1}{\text{s}}_1{\text{CH}}_2$$(9)

Third step is chain propagation (Yang et al., 2003):(10) $${\text{Hs}}_1\text{C}+{\text{H}}_2\iff {\text{Hs}}_1{\text{CH}}_2$$(10) (11) $${\text{CH}}_3{\text{s}}_1+{\text{Hs}}_1\Rightarrow {\text{CH}}_4+2{\text{s}}_1$$(11)

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):(12) $${\text{C}}_n{\text{H}}_{2n+1}{\text{s}}_1+{\text{H}}_2\Rightarrow {\text{C}}_n{\text{H}}_{2n+2}+{\text{Hs}}_1$$(12) (13) $${\text{C}}_n{\text{H}}_{2n+1}{\text{s}}_1\iff {\text{C}}_n{\text{H}}_{2n}+{\text{Hs}}_1$$(13)

where s₁ 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 listed in Table 1. According to water-gas shift reaction, the rate equation for CO₂ formation is as follows:(14) $$R_{{\text{CO}}_2}=\frac{k_v\left(\frac{P_{\text{CO}}{P}_{{\text{H}}_2\text{O}}}{{P_{{\text{H}}_2}}^{0.5}}-\frac{P_{{\text{CO}}_2}{P_{{\text{H}}_2}}^{0.5}}{K_p}\right)}{1+K_v\frac{P_{\text{CO}}{P}_{{\text{H}}_2\text{O}}}{{P_{{\text{H}}_2}}^{0.5}}}$$(14)

Table 1. Kinetic parameters used in the present study

Parameter	Value	Unit
k5,0	7.88 × 10^3	mol g^−1 s^−1 bar^−1
E5	75.52	kJ mol^−1
k7M,0	2.01 × 10^6	mol g^−1 s^−1 bar^−1
E7M	97.39	kJ mol^−1
k7,0	1.10 × 10^6	mol g^−1 s^−1 bar^−1
E7	111.48	kJ mol^−1
k8,0	8.79 × 10^3	mol g^−1 s^−1
E8	97.37	kJ mol^−1
kv,0	3.42	mol g^−1 s^−1 bar^−1.5
Ev	58.43	kJ mol^−1
k-8	2.77 × 10^−5	mol g^−1 s^−1 bar^−1
kV	2.76 × 10^−2	bar^−0.5
K1	2.59	bar^−1
K2	1.67 × 10^−3	bar^−1
K3	8.34 × 10^−2	[–]
K4	1.21	bar^−1
K6	0.10	[–]

where(15) $$K_p=\frac{5078.0045}{T}-5.8972089+13.958689\times {10}^{-4}T-27.592844\times {10}^{-8}{T}^2$$(15)

The reaction rate for methane, paraffins, and olefins formation are given:(16) $$R_{C{H}_4}=\frac{k_{7M}{K}_4{K}_6{K}_3^{\prime }\frac{P_{{\text{H}}_2}^3{P}_{\text{CO}}}{P_{{\text{H}}_2\text{O}}}}{{\left[\begin{array}{c}1+\sqrt{K_4{P}_{{\text{H}}_2}}+K_1{P}_{\text{CO}}+K_3^{\prime }\frac{P_{{\text{H}}_2}^2{P}_{\text{CO}}}{P_{{\text{H}}_2\text{O}}}+K_1{K}_2{P}_{\text{CO}}{P}_{{\text{H}}_2}\hfill \\ +K_6{K}_4^{0.5}{K}_3^{\prime }\frac{P_{{\text{H}}_2}^{2.5}{P}_{\text{CO}}}{P_{{\text{H}}_2\text{O}}}+K_3^{\prime }\frac{P_{{\text{H}}_2}^2{P}_{\text{CO}}}{P_{{\text{H}}_2\text{O}}}\left(1+K_6\sqrt{K_4{P}_{{\text{H}}_2}}\right){\sum }_{i=2}^n{\prod }_{j=2}^i({\mathit{\alpha }}_{ij})\hfill \end{array}\right]}^2}$$(16) (17) $$R_{{\text{C}}_n{\text{H}}_{2n+2}}=\frac{k_7{K}_4{K}_6{K}_3^{\prime }\frac{P_{{\text{H}}_2}^3{P}_{\text{CO}}}{P_{{\text{H}}_2\text{O}}}{\prod }_{j=2}^n{\mathit{\alpha }}_j}{{\left[\begin{array}{c}1+\sqrt{K_4{P}_{{\text{H}}_2}}+K_1{P}_{\text{CO}}+K_3^{\prime }\frac{P_{{\text{H}}_2}^2{P}_{\text{CO}}}{P_{{\text{H}}_2\text{O}}}+K_1{K}_2{P}_{\text{CO}}{P}_{{\text{H}}_2}\hfill \\ +K_6{K}_4^{0.5}{K}_3^{\prime }\frac{P_{{\text{H}}_2}^{2.5}{P}_{\text{CO}}}{P_{{\text{H}}_2\text{O}}}+K_3^{\prime }\frac{P_{{\text{H}}_2}^2{P}_{\text{CO}}}{P_{{\text{H}}_2\text{O}}}\left(1+K_6\sqrt{K_4{P}_{{\text{H}}_2}}\right){\sum }_{i=2}^n{\prod }_{j=2}^i({\mathit{\alpha }}_{ij})\hfill \end{array}\right]}^2}(n\ge 2)$$(17) (18) $$R_{{\text{C}}_n{\text{H}}_{2n}}=\frac{\left[\frac{K_8^{+}(1-{\mathit{\beta }}_n)K_3^{\prime }{P}_{{\text{H}}_2}^2{P}_{\text{CO}}}{P_{{\text{H}}_2\text{O}}{\prod }_{j=2}^n{\mathit{\alpha }}_j}\right]}{{\left[\begin{array}{c}1+\sqrt{K_4{P}_{{\text{H}}_2}}+K_1{P}_{\text{CO}}+K_3^{\prime }\frac{P_{{\text{H}}_2}^2{P}_{\text{CO}}}{P_{{\text{H}}_2\text{O}}}+K_1{K}_2{P}_{\text{CO}}{P}_{{\text{H}}_2}\hfill \\ +K_6{K}_4^{0.5}{K}_3^{\prime }\frac{P_{{\text{H}}_2}^{2.5}{P}_{\text{CO}}}{P_{{\text{H}}_2\text{O}}}+K_3^{\prime }\frac{P_{{\text{H}}_2}^2{P}_{\text{CO}}}{P_{{\text{H}}_2\text{O}}}\left(1+K_6\sqrt{K_4{P}_{{\text{H}}_2}}\right){\sum }_{i=2}^n{\prod }_{j=2}^i({\mathit{\alpha }}_{ij})\hfill \end{array}\right]}^2}(n\ge 2)$$(18)

where(19) $${\mathit{\alpha }}_n=\frac{K_5{K}_3^{\prime }\frac{P_{{\text{H}}_2}^2{P}_{\text{CO}}}{P_{{\text{H}}_2\text{O}}}}{K_5{K}_3^{\prime }\frac{P_{{\text{H}}_2}^2{P}_{\text{CO}}}{P_{{\text{H}}_2\text{O}}}+k_7{K}_4{K}_6{P}_{{\text{H}}_2}+K_8^{+}(1-{\mathit{\beta }}_n)}(n\ge 2)$$(19) (20) $${\mathit{\beta }}_n=\left(\frac{k_8^{-}}{k_8^{+}}\right)\left\{\frac{P_{{\text{C}}_n{\text{H}}_{2n}}}{{\mathit{\alpha }}_A^{n-1}{K}_3^{\prime }{P}_{\text{CO}}\frac{P_{{\text{H}}_2}^2}{P_{{\text{H}}_2\text{O}}}+\frac{k_8^{-}}{\frac{K_5{K}_3^{\prime }{P}_{{\text{H}}_2}^2{P}_{\text{CO}}}{P_{{\text{H}}_2\text{O}}}+k_7{K}_4{K}_6{P}_{{\text{H}}_2}+k_8^{+}}\times {\sum }_{i=2}^n\left({\mathit{\alpha }}_A^{i-2}{P}_{C_{(n-i+2)}{H}_{2(n-i+2)}}\right)}\right\}(n\ge 2)$$(20) (21) $${\mathit{\alpha }}_A=\frac{\frac{K_5{K}_3^{\prime }{P}_{{\text{H}}_2}^2{P}_{\text{CO}}}{P_{{\text{H}}_2\text{O}}}}{\frac{K_5{K}_3^{\prime }{P}_{{\text{H}}_2}^2{P}_{\text{CO}}}{P_{{\text{H}}_2\text{O}}}+k_7{K}_4{K}_6{P}_{{\text{H}}_2}+k_8^{+}}$$(21)

The rate constant in the above equations can be defined according to Arrhenius relation, which is consistent with what CFD solver needs:(22) $$k_i(T)=k_{i,0}{T}^{n}exp\left(-\frac{E}{RT}\right)$$(22)

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.

3. 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.

3.1. Momentum conservation equation

The conservation of momentum is expressed by:(23) $$\mathrm{\nabla }·\left(\mathit{\rho }\overrightarrow{v}\overrightarrow{v}\right)=-\mathrm{\nabla }·P+\mathrm{\nabla }·\left(\overline{\overline{\mathit{\tau }}}\right)+\left(\mathit{\rho }\overrightarrow{g}\right)+\overrightarrow{F}$$(23)

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.(24) $$\left(\frac{\mathit{\eta }}{k}+Q\right)u=\mathrm{\nabla }·\left[-pI+\frac{1}{\mathit{\epsilon }}\left\{\mathit{\eta }\left(\mathrm{\nabla }·u+{(\mathrm{\nabla }·u)}^T\right)-\left(\frac{2}{3}\mathit{\eta }-k\right)(\mathrm{\nabla }·u)I\right\}\right]+F$$(24)

3.2. Mass conservation equation

General equation for conservation of mass, continuity equation, is defined as follows:(25) $$\frac{\mathit{\partial }}{\mathit{\partial }t}(\mathit{\epsilon }\mathit{\rho })+\mathrm{\nabla }·(\mathit{\rho }\overrightarrow{v})=Q$$(25)

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.(26) $$\frac{1}{r}\frac{\mathit{\delta }}{\mathit{\delta }r}(\mathit{\rho }r{v}_r)+\frac{\mathit{\delta }}{\mathit{\delta }Z}(\mathit{\rho }{v}_Z)=0$$(26)

3.3. Energy conservation equation

Energy conservation equation valid for porous media including conduction flux terms is given by:(27) $$\mathrm{\nabla }·\left(\overrightarrow{v}\left({\mathit{\rho }}_f{E}_f+p\right)\right)=\mathrm{\nabla }·\left[k_{\text{eff}}\mathrm{\nabla }T-\left(\sum _j{h}_j{j}_j\right)+\left(\overline{\overline{\mathit{\tau }}}·\overrightarrow{v}\right)\right]+S_f$$(27)

The term on the left-hand side of Equation 27 is the heat transfer due to convection. On the right-hand 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):(28) $$k_{\text{eff}}=\mathit{\epsilon }{k}_f+(1-\mathit{\epsilon })k_s$$(28)

3.4. Species transport equation

The equation governing species transport through the reactor is of course well known:(29) $$\left(\frac{\mathit{\partial }{C}_{\mathit{\alpha }}}{\mathit{\partial }t}+v_r\frac{\mathit{\partial }{C}_{\mathit{\alpha }}}{\mathit{\partial }r}+\frac{v_{\mathit{\theta }}}{r}\frac{\mathit{\partial }{C}_{\mathit{\alpha }}}{\mathit{\partial }\mathit{\theta }}+v_z\frac{\mathit{\partial }{C}_{\mathit{\alpha }}}{\mathit{\partial }z}\right)=-\left[\frac{1}{r}\frac{\mathit{\partial }}{\mathit{\partial }r}(r{J}_{\mathit{\alpha }}^{*})+\frac{1}{r}\frac{\mathit{\partial }{J}_{\mathit{\alpha }\mathit{\theta }}^{*}}{\mathit{\partial }\mathit{\theta }}+\frac{\mathit{\partial }{J}_{\mathit{\alpha }z}^{*}}{\mathit{\partial }z}\right]+r_{\mathit{\alpha }}$$(29)

Assuming steady-state and irrotational flow (θ direction), Equation 29 diminishes such:(30) $$\left(v_r\frac{\mathit{\partial }{C}_{\mathit{\alpha }}}{\mathit{\partial }r}+v_z\frac{\mathit{\partial }{C}_{\mathit{\alpha }}}{\mathit{\partial }z}\right)=-\left[\frac{1}{r}\frac{\mathit{\partial }}{\mathit{\partial }r}\left(r{J}_{\mathit{\alpha }}^{*}\right)+\frac{\mathit{\partial }{J}_{\mathit{\alpha }z}^{*}}{\mathit{\partial }z}\right]+r_{\mathit{\alpha }}$$(30)

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.

4. 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.

Figure 2. Mass fraction of CO along the reactor length obtained for two different grids.

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.

Figure 3. Geometry and grid with 84430 structured quadrilateral elements (the lower boundary is axis of symmetry).

5. 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.

6. 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.

Table 2. Characteristics of fixed-bed reactor and catalyst

System-type parameter	Value
Reactor dimension (m)	7 × 0.05
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

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:(31) $$F_i=-\left(\sum _{j=1}^3{D}_{ij}\mathit{\mu }{v}_j\right)$$(31)

The first 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:(32) $$F_i=-\left(\frac{\mathit{\mu }}{\mathit{\alpha }}{v}_i\right)$$(32)

where α is permeability and C₂ 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).(33) $$\frac{\mathrm{\Delta }P}{L}=\frac{150\mathit{\mu }}{D_p^2{\mathrm{\varnothing }}_p^2}\frac{{(1-\mathit{\phi })}^2}{{\mathit{\phi }}^3}v$$(33)

Or in simpler form:(34) $$\frac{\mathrm{\Delta }P}{L}=\frac{\mathit{\mu }}{\mathit{\alpha }}v$$(34)

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:(35) $$\mathit{\alpha }=\frac{D_p^2}{150}\frac{{\mathit{\phi }}^3}{{(1-\mathit{\phi })}^2}$$(35)

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² and 12064 1/m, respectively, for viscous and inertial resistances. It must be said that viscous resistance is inversely associated with permeability.

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).¹ 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.

Table 3. Species mass fraction at the inlet

Component	Mass fraction
H2O	0.0014
H2	0.0224
CO	0.3245
CO2	0.001
C5H12	0.0004
CH4	0.0006
N2 and other components	0.6497

7. 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%.

Table 4. Comparison between CFD simulation results with experimental data

Run number	T (K)	P (bar)	H2/CO	CO conv. CFD (%)	CO conv. exp. (Yang et al., 2003) (%)	Relative error (%)
1	570.17	20.2	1.03	35.71	34	5.027
2	570.17	15	2.07	35.11	37	−5.113
3	556.27	15	1.03	21.8	22	−0.917

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.

Figure 4. Temperature profile along reactor length under different pressures.

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.

Figure 5. Temperature profile along reactor length for various feed temperatures.

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 fraction for CO and H₂ are illustrated in Figure 6, showing that mole fraction of CO and H₂ 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⁻¹, 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.

Figure 6. Contours of H₂ and CO mole fraction.

Predicted concentration of products through the axis of fixed-bed reactor is shown in Figures 7^{Figure 8Figure 9Figure 10Figure 11}–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.

Figure 7. Contours of C₂H₄ mole fraction.

Figure 8. Contours of C₂H₆ mole fraction.

Figure 9. Contours of C₃H₆ mole fraction.

Figure 10. Contours of C₃H₈ mole fraction.

Figure 11. Contours of C₄H₈ mole fraction.

Figure 12. Contours of C₄H₁₀ mole fraction.

Figures 7^{Figure 8Figure 9Figure 10Figure 11}–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₇⁺ 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.

8. Conclusion

Two-dimensional CFD technique was employed to predict the conversion of synthetic gas (CO and H₂) and distribution of products through fixed-bed reactor in the FTS process. Results show that conversion of CO and H₂ and the selectivity of C₇⁺ as it seen in 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.

Figure 13. Contours of C₇⁺ mole fraction.

Nomenclature
g	=	acceleration due to gravity (m s−2)
F	=	surface tension force components and also body force (N)
Q	=	mass source term (kg m−3 s−1)
C	=	molar concentration (kg mol m−3)
cp	=	heat capacity (J mol−1 K−1)
h	=	fluid enthalpy (kJ kg−1)
hj	=	enthalpy of species j (kJ kg−1)
Jj	=	diffusion flux of species j (kg m−2 s−1)
k	=	conductivity of fluid phase (W m−1 K−1)
α	=	permeability (m2)
keff	=	effective thermal conductivity (W m−1 K−1)
kf	=	fluid-phase thermal conductivity (W m−1 K−1)
ks	=	solid medium thermal conductivity (W m−1 K−1)
P	=	total pressure (bar)
ra	=	the et rate of production of species (mol g−1 s−1)
Sf	=	heat source (kJ m−3 s−1)
T	=	temperature of the reaction system (K)
t	=	time (s)
u	=	superficial velocity of liquid phase (m s−1)
ρf	=	density of fluid phase (kg m−3)
$\overrightarrow{v}$	=	velocity vector (m s−1)
τ	=	shear stress (N m−2)
ε	=	void fraction (–)
ε	=	porosity of the medium (–)
$\mathit{\eta }$	=	dynamic viscosity (N s m−2)
$\mathit{\gamma }$	=	specific weight (kg m−2 s−2)
$\mathrm{\varnothing }$	=	porosity of porous layer (–)
Ef	=	total fluid energy (kJ)

Additional information

Funding

Funding. The authors received no direct funding for this research.

Notes on contributors

Hamed Aligolzadeh

Hamed Aligolzadeh received his MSc in Chemical Engineering from Iran University of Science and Technology in Tehran, I.R. Iran. His research interests are GTL technology and CFD simulation.

Ali Jebreili Jolodar

Ali Jebreili Jolodar received his MSc in Chemical Engineering from Mazandaran University of Science and Technology, I.R. Iran. His research interests are GTL and numerical simulation. He manages above 10 national projects in the Research Institute of Petroleum Industry, I.R. Iran.

Rasool Mohammadikhah

Rasool Mohammadikhah is assistant professor of Chemical Engineering and a scientific member of the Research Institute of Petroleum Industry, I.R. Iran. He received his doctorate in Chemical Engineering from University of Tehran in Tehran (2014), I.R. Iran. His academic publications include up to 30 technical papers. His research interests are molecular thermodynamics, molecular dynamic simulation, and CFD simulation.

Notes

1. “Incompressible ideal gas” in the FLUENT 6.3 ed, FLUENT Inc., USA.