Effect of Particle Shape on Catalyst Deactivation during 2-Butene and Isobutane Alkylation of Liquid Phase in Fixed-Bed Reactor Using Particle-Resolved CFD Simulation

Zhang Sizhen; Zhu Zhenxing; Xin Feng; Chu Menghan

(1. School of Chemical Engineering and Technology, Tianjin University, Tianjin 300350, China; 2. SINOPEC Research Institute of Petroleum Processing Co., Ltd., Beijing 100083, China)

Abstract: How catalyst shape affects its deactivation is a crucial issue for quickly decaying catalysts such as zeolite in 2-butene and isobutane alkylation. In this work, steady simulations are used to determine the temperature and species distribution in fixed beds filled with particles of four shapes. Subsequently, unsteady simulations are used to study the deactivation behavior of the catalysts based on the steady simulation results. We describe the deactivation rate and type of catalyst deactivation by defining a local internal diffusivity, which is affected by catalytic activity. The results reveal that the internal diffusion distance of the catalyst determines the deactivation rate, whereas the local internal diffusivity determines its deactivation type.

Key words: alkylation; catalyst deactivation; particle shape; fixed bed; particle-resolved CFD

1 Introduction

The alkylation of isobutane and 2-butene is one of the most important petrochemical processes for forming longer chain alkylates, including branched alkanes such as trimethylpentanes (TMPs), which have high octane numbers and low Reid vapor pressures and contents of compounds containing sulfur and nitrogen, alkenes, and aromatics. This process, therefore, has been viewed as an environmentally friendly production of gasoline blending components[1-3].

Traditional alkylation processes in refineries are catalyzed by liquid acids such as concentrated sulfuric acid (H2SO4) or hydrofluoric acid (HF), which can induce severe environmental pollution and safety risks owing to high acid consumption, corrosiveness, and toxicity. Therefore, researchers have been focusing on developing alternative alkylation catalysts in recent decades. Solid acid catalysts are replacing liquid acids in commercial production due to their lower cost and operational safety. However, rapid catalyst deactivation remains a big challenge for solid acids, such as in the zeolite-catalyzed alkylation of isobutane and 2-butene[4], and on-line regeneration is essential. Among the heterogeneous solid acid catalysts, zeolites have been studied in great detail because of their suitable structure and pore size, as well as their relatively strong acidity, which enhances the desired product selectivity. However, zeolites deactivate fast due to the formation of coke.

Two strategies are normally used for slowing the deactivation of zeolite catalysts: one is to ensure the high concentration of strong Bronsted acid sites governed by the aluminum content in the zeolite and exchange procedure; another is to enlarge the molar ratio of isobutane to 2-butene and use continuous stirred tank reactors (CSTRs) to lower 2-butene concentration, whose oligomerization is the main cause of catalyst deactivation. However, the separation of solids and liquids in the CSTR slurry complicates the reactor structure and operation. Commercial plants with large production capacity usually adopt a fixed-bed reactor due to its simple design and operation[5]. Zuazo et al.[6]argued that the key issue of isobutane/butene alkylation in the fixed-bed tank is the rapid deactivation of the catalyst, whose frequent regeneration has to be ensured in industry. This situation has hindered the steady-state operation of commercial fixed-bed tanks and increased their running cost. In view of its economy, the catalyst operation period must be prolonged to reduce regeneration.

The fixed-bed reactor is filled with larger catalyst particles than those in the CSTR. Most parts of the above two strategies are suitable for the fixed-bed reactor, where the catalyst must tolerate a higher concentration of 2-butene near the reactor inlet and has a longer lifetime. The catalyst shape and size determine the resistances of fluid flowing through the catalyst bed and diffusing inside the catalyst particle. Herein, we attempt to determine the relationship between catalytic activity, lifetime, and catalyst shape for optimizing its geometry and activity. Catalyst shapes for the fixed-bed tank include sphere, cylinder, ring, and complex shapes such as multilobe or multichannel cylinder. These particles are used to construct the fixed beds, whose conventional simulations are based on pseudo-homogeneous and heterogeneous models. However, all of these models are too macroscopic to describe local transport phenomena like local mixing and the near-wall effect[7]. The chemical and physical processes in both the bed and the catalyst scales must be understood more fully for better designing and operating the reactor.

In order to precisely characterize the micro-scale reaction and transport occurring in the fixed-bed reactor, the particle-resolved CFD method has emerged for simulating the complex flow around the particles and providing detailed local information inside the bed[8]. Most of the particle-resolved CFD simulations reported in the literature are for spherical or cylindrical particles, and few research the effect of complex particle shape on reactor performance. Karthik et al.[8]investigated the effect of particle shape on catalyst deactivation using five different complex shapes of particles and showed that particles of lower surface area favored propane dehydrogenation due to slower catalyst deactivation. Shen et al.[9]concluded that a certain amount of holes in cylindrical catalysts potentiate heat transfer to obtain high CH4yield in the Sabatier reaction.

For the alkylation reaction in the fixed reactors, much attention has been paid to the catalyst type and reaction mechanism, whereas the selection of catalyst shape was ignored. In our previous work[10], we studied the effect of particle shape on the fixed-bed reactor performance for the alkylation reaction and optimized the catalyst shape by considering the fluid flow, heat transfer, and internal diffusion. However, the study was limited to only the alkylation reaction without catalyst deactivation. Apart from the strong influence of internal diffusion, the study of catalyst deactivation type and rate is also an extremely important issue in the solid-acid-catalyzed alkylation. Our work aims to establish an efficient method to simulate the rates of the reaction and catalyst deactivation during C4alkylation in the fixed-bed reactor and optimize catalyst particle shapes. Herein, we numerically simulate catalyst deactivation of the alkylation reaction in fixed beds filled with particles of four different shapes by using the particle-resolved CFD method. The mutual influence between the reaction rate and internal diffusivity variation is analyzed in detail, and we compare the simulated results with experimental data. Moreover, the type of catalytic deactivation and the effect of catalyst shape on deactivation rate are studied to better understand the solid-acid-catalyzed alkylation.

2 Reaction Kinetics

Understanding the relevant reaction kinetics is essential for designing the reactor and optimizing its operation. A kinetics model proposed by Taylor et al.[11]considered the alkylation reaction as isobutane and 2-butene converting into trimethylpentane.

Under a high isobutane to 2-butene ratio, an apparent and pseudo-first-order reaction rate law in an ultra-stable Y zeolite catalyst of 1.6 diameter spherical pellet was modified by Chu et al.[10]to extrapolate an intrinsic kinetic model in Eq. (1), and their parameters are listed in Table 1:

Table 1 Intrinsic kinetic parameters

wherea(t) andkiare the deactivation function and intrinsic rate constant, andki0andEintare the intrinsic frequency factor and intrinsic activation energy, respectively. Moreover, the deactivation rate constantkdis obtained from Reference [11] and expressed as follows:

It is noteworthy that the effect of deactivation on internal diffusion is not considered when deriving the intrinsic rate law. In Section 4.2.3, we compensate for the mutual interaction of the deactivation and internal diffusion with an effective diffusivity.

3 Methodology

Our numerical simulation consists of three steps: geometry building, meshing, and numerical calculation. First, a random stacked fixed bed made of cylindrical catalysts is generated by the rigid body dynamic (RBD) method, and the position information of every particle is extracted out. Similarly, the rest of the three fixed beds filled with different particle shapes is generated by maintaining the same packing arrangement for rational comparison. Then, the domain is meshed in both fluid and solid regions. Finally, the governing equations are solved in the domain with proper boundary conditions to quantify velocity, concentration, and temperature fields by the commercial STAR-CCM+ software package.

3.1 Geometry

Three different particle shapes—ring, trilobe, and quadrilobed—are constructed from the cylindrical particle. The three particle shapes are created with equal particle voidage (0.4) based on the reference particle volume of the cylindrical particle. The random packing algorithm based on RBD embedded in open-source software Blender was applied by Partopour et al.[12]to generate packed beds for different particle shapes. The flow behavior of the generated packed bed was consistent with the experiments, which proved the feasibility of the method. In our work, we use the same method and preset the parameters as shown in Table 2.

Table 2 Parameters of RBD

In order to save computational time, a quarter volume of the fixed bed was calculated under an assumption of symmetry[13-14]. The final packing arrangement is shown in Figure 1. The particles are packed in a cylindrical tube with a total height of 35 mm and a diameter of 16 mm, whose tube-to-particle diameter ratio is 8.

Figure 1 The particle packing arrangement and fixed bed geometry

3.2 Meshing

The gaps method is used to void the bad cell quality in the contact regions[9,15]. The particle volume is shrunk by 0.5% to deal with the issue of contact region in the present work. Polyhedral meshes are applied in both fluid and porous solid regions. Moreover, to capture the temperature and species gradients at the particle surface, boundary layer prism cells are set at the fluid–solid interface and the tube wall. A mesh independence study is performed with three different mesh resolutions for the fixed bed filled with cylindrical particles. A suitable mesh size setting is determined, and the parameters for meshing are listed in Table 3.

Table 3 Parameters used for mesh generation

3.3 Governing equations

For an incompressible laminar flow in the liquid phase, the governing equations (including the continuity, momentum, species transport, and energy conservation equations) can be found in the literature[9,16]. In the porous solid phase, we only consider the species transport and energy conservation equation due to the absence of convection. The species transport and energy conservation equations are given by Eqs. (3) and (4), whereSiandSTare the species and heat source terms due to chemical reactions in the solid phase:

The effective diffusivity (Deff,i) depends on the porosity (ε), and tortuosity (τ) is calculated by Eq. (5), where the molecular diffusivity (Di) is calculated by the Wilke–Chang equation[17]. The effective heat capacity (Cp,eff) and thermal conductivity (λeff) in the porous medium region are expressed in Eqs. (6) and (7):

whereCp,fandCp,sare heat capacity for fluid and solid,λfandλsare the thermal conductivity for fluid and solid.Cp,effandλeffare the weighted values of mass fractions of all components, whose individual values are obtained from the NIST (https://webbook.nist.gov) and fitted by polynomials versus temperature, as shown in Table 4.

Table 4 Physical properties of each component

3.4 Boundary conditions and simulation parameters

We set boundary conditions for the velocity at inlet and pressure at outlet to simulate laminar flow over the catalyst particles in a fixed bed; the wall temperature of the reactor is constant, and all solid surfaces are regarded as non-slip. Both the fluxes of mass and heat are conserved on the fluid–solid interface, and the boundary conditions are as follows:

The initial parameters of the simulation and boundary conditions are shown in Table 5. To compare the deactivation of catalysts with different particle shapes in fixed-bed reactors, the mass space velocity of 2-butene (OSV) is set to be a constant of 0.8 kg2-butene/kgcat/h according to industrial operating conditions. The OSV of the reaction is low in value to avoid a fast deactivation of the solid catalyst.

Table 5 Boundary conditions and simulation parameters

3.5 Numerical method

The governing equations of liquid and porous solid regions are solved separately and the data are exchanged with each other to achieve flux conservation, as shown in Figure 2. First, the leading simulation in the fluid phase updates its physical property after a certain number of iterations considering the fluid flowing. Second, the information of mass fraction and heat flux on the boundary is read out and mapped to the boundary of the lagging simulation. Third, the lagging simulation in the porous solid phase iterates certain steps and updates the physical properties considering the reaction-diffusion process. Finally, the mass flux of each species and temperature on the boundary are solved for the continuous simulation until achieving convergence.

Figure 2 Schematic for solving multi-region governing equation

The simulations implemented in our research are carried out in the commercial CFD software STAR-CCM+12.02 using the finite volume method to solve the governing equations. Convergence is achieved by tracking the residual curves of the overall mass and energy. We only consider an unsteady simulation in the solid domain, where the catalyst deactivation occurs. A variable time step is set in our unsteady simulation to reduce computational cost and maintain calculation accuracy. Based on the complexity of the calculation, the simulations are performed in parallel using multiple CPU cores (60–120).

4 Results & Discussion

4.1 Steady simulations

Steady simulations are used to investigate the mass and heat transfer on the fresh catalyst. Additionally, in this section we calculate the internal diffusion effectiveness factor and 2-butene conversion to screen the catalyst shapes.

The mass fraction profiles on the two planes and three planes of both fluid and solid regions in Figure 3 show strong mass fraction gradients inside particles due to internal diffusion limitations, and lower mass fraction gradients are found in the beds with trilobes, quadrilobes, and rings compared with the bed with cylinders due to their better mass transferring characteristics.

Figure 3 Effect of particle shape on mass fraction profile: (a) two planes, (b) three planes, (i) cylinders, (ii) trilobes, (iii) quadrilobes, and (iv) rings. 323 K, w2B=0.00116, OSV=0.8 kg2-butene/kgcat/h

The temperature profiles on the two planes and three planes of both fluid and solid regions in Figure 4 show that the bed temperature gradually increases along the flow direction due to the exothermic reaction, and we find a small increase in temperature in all four fixed beds due to the low reactant concentration and small reaction heat. Moreover, lower temperature gradients are observed in the beds with trilobes, quadrilobes, and rings compared with the bed with cylinders.

Figure 4 Effect of particle shape on temperature profile (a) two planes, (b) three planes, (i) cylinders, (ii) trilobes, (iii) quadrilobes, and (iv) rings. 323 K, w2B=0.00116, OSV=0.8 kg2-butene/kgcat/h

The internal diffusion effectiveness factorη, defined as Eq. (10), is the ratio of the volume-averaged reaction rate to the surface-averaged reaction rate of a catalyst particle.

whereVpandApare the pellet volume and pellet surface, respectively. Theηis calculated to quantify the extent of internal diffusion limitations for different particle shapes. The ring shape is found to have the highest value ofη, as shown in Figure 5, because the ring has the shortest diffusion distance from its external surface to its center, leading to easier access to the interior of the catalyst for the reactants. In addition, the 2-butene conversion in Figure 5 is positively associated with the value ofηunder the same OSV. The order of 2-butene conversion rate for the catalyst bed arrangements is ring > quadrilobe > trilobe > cylinder.

Figure 5 Effect of particle shape on internal diffusion effectiveness factor (η) and conversion of 2-Butene. 323 K, w2B=0.00116, OSV=0.8 kg2-butene/kgcat/h

4.2 Unsteady simulations

After simulating the reaction on fresh catalyst, we input the outcomes of temperature and mass profiles as initial conditions for unsteady simulation. The catalyst deactivation process is simulated to observe the timevariations of both catalytic activity and reaction rate profiles inside the particles, and the impact of different particle shapes on the deactivation is also quantized. In addition, the simulation results are compared with the experimental results.

4.2.1 Reaction on deactivated catalyst

We take the cylindrical catalyst as an example to describe the deactivation process. The variation of catalytic activity on the three planes of the solid region is shown in Figure 6(a), where the overall catalytic activity decreases with time because of the continuous deactivation of catalyst. At the initial stage, the descending degree of the activity near the fluid–solid interface is more significant than that in the inner zone, which can be attributed to the higher reactant mass fraction in this area. The catalytic activity in the inner zone begins to decrease when the activity in the shell zone is at a low level, which means that the reactants diffuse into the inner zone of the particles, and catalyst deactivation occurs gradually from the outside to the inside.

The variation of alkylation reaction rate on the three planes of the solid region is shown in Figure 6(b), where the reaction rate tends to decrease with time due to the catalyst deactivation. In addition, the distribution of the reaction rate is obviously non-uniform due to the uneven distribution of reaction mass fraction and catalytic activity. At the initial stage, a higher reaction rate is observed in the near-wall region due to the high reactant mass fraction and catalyst activity, whereas a lower reaction rate is observed within particles due to the internal diffusion limitation. As time progresses, the catalytic activity in the near-wall region decreases to a low level, which almost stops the reaction and make more reactants diffuse to the inner zone. In this stage, the reaction rate inside the particles is higher than that in the near-wall region. Finally, the reaction rate decreases to zero in the whole particle zone, as the catalyst is completely deactivated.

To quantify the effect of different particle shapes on deactivation, the time-variations of the volume-averaged activity [Eq. (11)] and area-averaged activity [Eq. (12)] with different particle shapes are shown in Figure 7. A faster decline of area-averaged activity than volumeaveraged activity is observed due to faster deactivation on the outer surface of the catalyst. Furthermore, the volume-averaged activity of the ring is found to decrease faster than other particle shapes, which means a faster deactivation inside the particles. This is because the ring has the shortest diffusion length and a better access for reactants to the particle inner zone, leading to an acceleration of deactivation processes related to the reactant concentration.

Figure 6 Time variations of cylindrical catalyst for (a) catalytic activity and (b) reaction rate. 323 K, w2B=0.00116, OSV=0.8 kg2-butene/kgcat/h

The variableX, defined as Eq. (13), is used to quantify the difference of catalytic activity between the inside and outside of the particles. The values ofXfor different particle shapes are shown in Figure 8. The value of outer surface area to volume ratio (STV), which is normalized by being divided by the corresponding value of the cylinder particles, is also shown in Figure 8. A larger value of STV corresponds to better transferring characteristics and access for reactants to the catalyst interior. Therefore, the difference in the degree of deactivation between the outer surface and inner zone of the catalyst becomes smaller, leading to a smaller value ofX.

Figure 7 Time-variations of the volume-averaged activity and area-averaged activity for (a) cylinders, (b) trilobes, (c) quadrilobes, and (d) rings. 323 K, w2B=0.00116, OSV=0.8 kg2-butene/kgcat/h

Figure 8 Effect of particle shape on catalyst activity deviation and surface area-to-volume ratio (STV). 323 K, w2B=0.00116, OSV=0.8 kg2-butene/kgcat/h

4.2.2 Comparison with experimental result

To verify the efficiency of the simulation result, a comparison is made between the experimental and simulated results for the spherical catalyst. The time evolution of the experimental and simulated 2-butene conversion is given in Figure 9. It should be noted that the reaction conditions remain the same between experiment and simulation. In addition, we capture a significant difference on the decline of the 2-butene conversion between Taylor’s work and simulated result. Under the same reaction condition, the simulated result shows a slower decrease for the reactant conversion. It reveals the necessity to revise our model. We find the main reason for the deviation is that the deactivation rate expression was obtained on the basis of a 1.6-mm-diameter spherical pellet without eliminating the effect of deactivation on internal diffusion. Therefore, the model is improved in the following section.

Figure 9 Time-evolution of 2-Butene conversion over catalyst pellets. (♦) Taylor et al., dp=1.6 mm; (●) Simulation, dp=1.6 mm. 360 K, OSV=0.2 kg2-butene/kgcat/h

Figure 10 Time variations of cylindrical catalyst for (a) catalytic activity and (b) reaction rate n=1, 323 K, w2B=0.00116, OSV=0.8 kg2-butene/kgcat/h

4.2.3 Modified diffusivity

Catalyst deactivation is closely related to the internal diffusivity inside particles. It is generally believed that irreversible adsorption of heavier hydrocarbons on the acid site is the major cause of catalyst deactivation for the alkylation reaction catalyzed by a solid acid. Moreover, the carbonaceous compounds adsorbed on the poisoned sites will lead to pore plugging. Here, the effective diffusivity (Deff,d) of the deactivated region is modified as a function of local activity of the catalyst pellet and can be expressed as follows:

wherenis an index representing the effect of the activity on the local effective diffusivity and ranges from 0 to 2.

Takingn=1, for example, the variations of the activity and reaction rate distribution on different planes inside the cylindrical particles are shown in Figure 10. At the initial stage (0.1 h), a higher reaction rate is observed in the near-wall region due to the high reactant mass fraction and catalyst activity, as shown in Figure 10(b). At the same time, a rapid decline of activity occurs in the shell zone [Figure 10(a)], which reduces the diffusivity and limits the access for reactants to the inner zone of the particles. In spite of a high activity inside the particles, the reaction rate is almost zero throughout the whole region at last (20 h).

The time-variations of the volume-averaged activity and area-averaged activity for cylindrical particles with different values ofnare shown in Figure 11. We can recognize that the descending degree of the area-averaged activity is similar for different values ofn, which is characterized by a fast decrease to zero. However, the rate of decrease for the volume-averaged activity significantly slows down with an increased value ofn. The volumeaveraged activity is observed to remain at a high level whennis bigger than 1. This is because the deactivation of the particle surface reduces the effective diffusivity heavily, thereby leading to the underutilization of catalyst inside the particles. In other words, the type of catalyst deactivation is pore plugging.

The evolution ofXfor different values ofnand different particle shapes is shown in Figure 12. TheXvalue increases withnbecause of the stronger effect of deactivation on internal diffusion. For the samen, the value ofXfor the cylinder is the largest, corresponding to the most serious pore plugging. In addition, when the value ofnis 2,Xis large for all catalyst shapes, indicating that serious pore plugging is occurring.

To evaluate the impact ofnon the reactor performance, the time-variations of the reactant conversion with different particle shapes at differentnvalues are shown in Figure 13. Irrespective of the particle shapes, the descending degree of conversion is more significant at highernvalues because the deactivation exerts an unfavorable impact on the effective diffusivity, leading to a lower 2-butene conversion. At this point, the fastest decrease of conversion is observed for the ring at the samenvalue. The reason for this phenomenon is that a shorter diffusion path makes it easier for reactants to diffuse to the interior of the catalyst, leading to a larger deactivation rate constant.

Simulations with differentnvalue are performed to compare with the experimental results. From Figure 14, we find that the simulation results agree well with the experiment results whenn=1, which indicates that the modified effective diffusivity makes the simulation result reasonable. Under the condition ofn=1, the timevariations of the reactant conversion with different particle shapes are shown in Figure 15. The fixed beds filled with rings, trilobes, and quadrilobes are found to have a higher 2-butene conversion rate compared to the cylinders during the whole deactivation stage, which also proves that catalyst shape optimization is of great importance in dealing with the deactivation.

Figure 11 Time-variations of the volume-averaged activity and area-averaged activity for cylindrical catalysts with different values of n. 323 K, w2B=0.00116, OSV=0.8 kg2-butene/kgcat/h

Figure 12. Evolution of the activity deviation (X) with the value of n for different particle shapes. 323 K, w2B=0.00116, OSV=0.8 kg2-butene/kgcat/h

Figure 13 Time-variations of reactant conversion with time for different particle shapes at different n values. 323 K, w2B=0.00116, OSV=0.8 kg2-butene/kgcat/h

Figure 14 Time-evolution of 2-butene conversion over catalyst pellet of dp=1.6 mm. 360 K, OSV=0.2 kg2-butene/kgcat/h

Figure 15 Time-variations of the reactant conversion with time for different particle shapes at n=1. 323 K, w2B=0.00116, OSV=0.8 kg2-butene/kgcat/h

5 Conclusion

We numerically simulate the catalyst deactivation of an alkylation reaction in fixed-bed reactors filled with different particle shapes by using particle-resolved CFD and conclude the following:

The steady simulations show that the fixed bed filled with rings exhibits lower gradients of species due to its shorter diffusion distance compared with beds filled with cylinders, trilobes, or quadrilobes. In addition, 2-butene conversion increases with the increased internal diffusion effectiveness factor under the same mass space velocity of 2-butene.Based on the steady simulation results, unsteady simulations are carried out to study the catalytic activity considering the catalyst deactivation. The effective internal diffusivity is modified by comparing the simulation results with referenced experimental results. Pore plugging behavior is revealed by observing the catalytic activity difference between the outer surface and the inner zone of the catalyst. Moreover, the effect of particle shape on the deactivation of catalyst is investigated. The results show that the fixed bed of rings exhibits a higher 2-butene conversion rate during the whole deactivation stage compared with the beds of the other three shapes due to the shorter diffusion distance of the ring.

Acknowledgements:We acknowledge financial support from National Engineering Research Center for Petroleum Refining Technology and Catalyst (RIPP, SINOPEC, Grant No. 33600000-20-ZC0607-0009).