Compurers
Pergamon
them.
Engng
Vol. 19. Suppl., Copyright
009%.1354(95)00098-4
Printed
in &eat
0
pp. S351-S356,
1995 Elsevier Britain.
All
009%1354/95
1995
Science
rights
Ltd
reserved
$9.50
+ 0.00
SIMULATION OF TUBULAR REACTORS PACKED WITH LARGE-PORE CATALYSTS WITH SPHERICAL GEOMETRY Rosa M. Quinta-Ferreira and Paul0 M. Simbes Department of Chemical Engineering, University of Coimbra, 3000 Coimbra, Portugal Alirio E. Rodrigues Laboratory of Separation and Reaction Engineering, University of Porto, Porto, Portugal
ABSTRACT The effect of intraparticle convection on the behavior of fixed-bed reactors with spherical catalyst particles is analyzed. In the model equations for the particles, both radial and angular coordinates are accounted for. The study concerning the number of discretization points to be used in the numerical solution of the intraparticle profiles is presented. For the reactor, a comparison between unidimensional heterogeneous models with and without intraparticle convection and with spherical and slab geometry for the solid particles is also performed.
KJxwoRDs
Intraparticle convection; spherical particles; catalytic reactors
INTRODUCTION Nir and Pismen (1977) and Nir (1977) presented pioneering results showing the enhancement of the catalyst efliciency when the additional intraparticle convective transport is allowed inside isothermal slab, spherical and cylindrical catalyst particles. Also at the particle level, Lopes et al. (1993) studied the efficiency of non-isothermal catalysts with spherical geometry. In the sequence of previous studies concerning the behavior of packed-bed reactors, where the catalyst particles promoting intraparticle convection were assumed with slab geometry (Rodrigues and Quinta Ferreira,l990; Quinta Ferreira et al. 1992a,b), the present work performs the analysis of such systems considering now the more realistic spherical shape for the solid particles. In this case a two-dimensional problem arises since internal gradients in both radial and angular particle coordinates are developed; special attention was given to the boundary conditions of the mathematical model. The selective oxidation of o-xylene to phthalic anhydride was used as an example of an exothermic reaction. In this paper we have assumed an irreversible first order kinetics and a complex reaccional scheme will be considered in a future work. Since in the first case an analytic solution can be achieved for the intraparticle concentration profiles of isothermal particles (Nir and Pismen, 1977), it was possible to test the computational program by comparing both numerical and analytical results, having in mind its further use on situations which can only be solved numerically
S352
European Symposium on Computer Aided Process Engineering-5
MODEL EQUATIONS The dimensionless model equations for the unidimensional heterogeneous intraparticle di&sion/convection model with spherical particles, HTk,e are presented in Table l., being the mass and thermal balances for the fluid phase expressed by eqs. (l)-(3). In the modelization of the solid particles, the internal thermal gradients were neglected since the Prater thermicity factor was low. Therefore, the temperature of the solid in each axial position of the reactor is obtained through eq. (4), which represents the thermal balance to the fluid/particle interface. Table 1 - Dimensionless Equations for the heterogenous intraParticlediflkionkonvection aalyst Particles, HTA, p
model with spherical
Fluid Phase MassBalance
(11
Thermal Balance
(2) z*=o:
Boundary Conditions
; @,,=I
fb’l
(3:
Fluid - Particle Interface Nfi [@,-@b]=Nf B [f,,-fs]
(4)
CatalystParticle MassBalallc42
VZfp = x,
( 1-c12 1!!lL+,: r
V2fP =-y(
with , Boundary Conditions r=l,
p~+km--
fp
(51
ap
r’~)++$F”)~)
Ip)=l,Vr:
&L: Bim[fb-fs]=IJ$ifs
afP =-&=O
(6) (7)
(8)
;
l=o)
p=o:
$0
(91
Due to the existence of the intraparticle convective flux, the concentration profiles are noi symmetric in relation to the center of the sphere as in the case of the conventional catalysts where diision is the only mechanism taken into account. So, two spatial coordinates must be considered, the radial position, r, and the polar angle 0; in respect to the angle cp there are no concentration gradients. The mass balance is then described by an elliptical partial ditferential equation of second order, eqs. (5) and (6), being the angular position represented by u, with u = cos 0. The two boundary conditions for the angular coordinate were based in the symmetry of the concentration profiles when 8=0 and x, or Iul=l, eq. (7) (Fig. 1). One of the boundary conditions for the radial coordinate was defined through the concentration at the particle surf&e, which was assumed equal for all the values of 8, being given by the mass balance to the tilm surrounding the solid particle, eq. (8). The other boundary condition was obtained by looking at the concentration gradients at the particle center. The concentration profiles and their derivatives are continuous functions along the radial coordinate for any angular position, 8 and 8+x, Fig. 1. Since r varies between zero at the particle center, and one at the particle surface , it will come eq. (lOa) for r-0. Moreover, one can also state for the particle center that, due to symmetry, the derivative in respect to r will be equal for both x+0 and x-0, eq. (lob). So, accounting for eq.
European
Symposium
on Computer
Aided Process
Engineering-S
s353
(IOa), the derivatives for 8 and x-8 have opposite signals, eq. (10~). As a consequence, it exists an angular position where the derivative is null which will happen for 0=x/2. Therefore, one finally gets the second boundary condition for r=O and B=x/2 (u=O), by using eq. (9) of Table 1. r = 0:
=- af, ar +0 t =i
-%P a
1 n-0
(lo@ 1n-0
Fig. 1 - Diagram of the spherical Particle
In each reactor axial position, the intraparticle concentration profiles along the radial and angular coordinates were obtained by solving eqs. (5)-(g) with the package SLDGL, which implements the finite differences method with variable order and non-uniform grid (Glotz et al., 1984). Since this code is used for integration of elliptic equations defined in rectangular regions, the boundary conditions for the sphere referred above were adapted to a rectangular domain. The knowledge of the intraparticle profiles allowed the calculation of the effectiveness factor for the particle, 1, defined as the ratio between the observed reaction rate inside the catalyst and the intrinsic rate which would occur at the catalyst surface conditions, eq. (1 la). Since in this case, there are only gradients in r and 8, q will be given by eq. (1 lb). This efficiency was calculated by using the Simpson rule for the integrals.
JJIv,‘wpm dvp ‘= JJJv$(C,,T,)dVp
J:gfr rdrd8 (‘la)
’
‘=
fn
(1 lb)
For the heterogeneous model which neglects the intraparticle convective fluxes, model HTde, there are only concentration gradients along the radial coordinates; therefore, there is no need to account for the angular coordinate. However, the model equations presented above can be extended for this situation by putting hm=O, which eliminates the intraparticle convective transport. The concentration and temperature profiles in the fluid phase where finahy Catculated by solving the corresponding system of ODE’s of first order (eqs. (l)-(3)) through the GEAR package for initial value problems which makes use of the Adams method of variable order and variable step (Hindmarsh, 1974).
COMPUTERRESULTS In the first part of this work we have analyzed the behavior of the spherical catalyst particles separately, in order to optimize the numerical parameters for the computational solution of the intraparticle concentration profiles. We have tested different numbers of discretization points in both radial and angular coordinates, by comparing the values of the effectiveness factor calculated numerically with those obtained through the analytical solution. For low and high Peclet numbers, Fig. 2(a) (&=3) and Fig. 2(b) (Xrn=30) show the relative error of the numerical solution and the CPU time when a VAX cluster 8530 is used, as a function of the number of discretization points in both coordinates: Nr, for r and Na for u; these results were obtained for the solid temperature of 609 K. When Na increases from 5 to 11 a decrease in the relative error
s354
European
Symposium
on Computer
Aided Process
Engineering-5
of about 8% can be observed for km=3 and 20% for km=30. In both cases, a further increase of Na until 21 doesn’t improve the final solution, being , however, the computing time about four times higher than when 11 points are used. On the other hand, the relative error is approximately maintained for Nr, higher than 11. One can still observe that for low Peclet numbers a decrease of N, leads to lower errors, since the profiles are almost symmetric being then well represented with a low number of discretization points. Neverthless, one can not be very confident in the numerical solution for the cases with high h, values when low Nr are used. In fact, when the intraparticle convective transport increases, the assymetry of the concentration profiles is also more pronounced, requiring a high number of discretization points for their representation. Figures 3 compare the numerical profiles (curves 1) for h=3 and 30 when Nr=5, in (a), and Nil 1, in (b), with those calculated analytically (curves 2), when the pellet temperature is 625K. 4.0 2
A
4.7 .
6 r
iP
D 4.3
B w
4.5
F
4.1 0
5
10
15
20
25
40
% 20
1
ri 10 3
30 0
k
0
5
(a)
IO
Nr
a
IS
20
25
0
10 @I
20 Nr
Fig. 2 - Relative errors and CPU time as a function of Nr for ditIerent N, values: A-5, B-7, C-l 1, D-2 1 (a) Xm=3 ; (b) &=30.
30
-I
a.5 @)
0
0.5
I
r
Fig. 3 - Intraparticle concentration protiles for different &,:A-3, B-30 [l-nun., 2-anal.]. (a) Nr=5 ; (b) Nr=ll
Allowing the entrance of more reactant by convection, by using higher Peclet numbers, the intraparticle mass gradients are reduced and the catalyst efficiency will increase, as can be observed in Fig. 4, where the effectiveness factors calculated for diierent hm values and different solid temperatures are shown. Moreover, for a same Xm, higher temperatures will increase the reaction rates which will lead to higher internal gradients, decreasing then the catalyst efficiency. The results presented before were obtained with a first order method in both coordinates, which revealed to be the more adequate compromise between the final solution and the computing time. In fact, by using an order equal to 3 and 5 one could only get an improvement of ahout 2% on the catalyst efficiency with CPU times twice higher. At the reactor level, quite good results were obtained with Ni7, Na=l 1 and a 1st order method.
The analysis of the reactor behavior when large-pore catalysts are used, was performed by comparing the predictions of the complete heterogeneous intraparticle diflbsionkonvection model with spherical particles, HTd ce (krn=lO), with those obtained for the heterogeneous model considering only intraparticle dkirsion HTde. In both cases, the results were compared with the situations where slab geometry was used for the catalyst particles, models HTdc,p and HTdp (Quinta Ferreira et al., 1992).
European
Symposium
on Computer
Aided Process
s355
Engineering-5
.. 02
00 (a)
Fig. 4 - Efficiency of the spherical particles for different temperatures and h,: A-O; B-10; C-25; D-50.
0:4 z*
06
0
02 (b)
0.4
0.6
0.8
1
2”
Fig.5 - (a) Axial temperature profiles and (b) o-xylene concentration profiles in the bulk phase of the reactor for different models A-PIT: B-H&+ C-m&,,’ D-m&; E-m&
A comparison with the simple pseudo-homogeneous model, PH, which neglects all inter and intraparticle resistances was also obtained. Figures 5 show in (a) the temperature profiles in the bulk phase and in (b) the o-xylene concentration, for all those models, when the inlet and wall temperatures were equal to 660 K and the o-xylene partial pressure at the reactor entrance was 0.014 atm. The additional intraparticle convective transport will allow more reactant going inside the catalyst, which will increase the reaction rate. As a consequence, when the same operating conditions are used, higher temperatures and lower o-xylene concentrations, i.e., higher final conversions will be achieved than when the internal convection is neglected either for spherical or slab geometry for the particles. However, due to the enhancement of the catalyst efficiency when intraparticle conversion is taken into account, less severe operating conditions can be used in order to get the same final conversions as the ones obtained with the conventional catalysts. Consequently, lower hot-spots are developed inside the catalytic bed being then possible to increase the catalyst life. The profiles predicted in the cases where the catalyst particles are assumed with slab geometry are more abrupt than when spherical geometry is considered which can be observed for both models by comparing curves B and C with curves D and E in Figs. 5. Moreover, the predictions of the PH model are very different from the ones obtained with the heterogeneous models showing that the heterogeneity of the system has to be taken into account in the analysis of such strongly exothermic systems.
CONCLUSIONS An analysis of fixed-bed reactors packed with spherical particles was presented by modelling the system through unidimensional heterogeneous models accounting for intraparticle diision and convection. The mass balances to the catalyst particles led to elliptical PDE’s where the radial and angular polar coordinates were included. The numerical profiles for the particle were compared with the analytical solution since a pseudo-irreversible first order reaction was assumed. The results obtained at the reactor level were compared with the ones obtained with those heterogeneous models where only intraparticle diffusion was considered and also with the simple pseudo-homogeneous model. For the same operating conditions, when the intraparticle phenomenon is taken into account, the final reactant conversion is higher than when this transport is neglected . By comparing also these results with the ones obtained with the reactor model considering the catalyst particle with slab geometry, we observed that this model predicts higher hot spots and higher conversions. The use of large-pore catalysts promoting the additional transport by convection inside the solid seems to be an interesting way of increasing the efficiency of such industrial heterogeneous processes which justifies the need for the identification of their behavior .
S356
European Symposium on Computer Aided Process Engineering-S
ACKNOLEWDGMENT The authors are gratefblly acknowledged to Professor Maria da Graqa Rasteiro for the advise on the computational program and for the discussion about the numerical solution of the spherical particles.
NOTATION
temperaturerise
Da - Damkoler number (Da=WoI JP .JP n De - effective & sw~ty
f - dimensionless reactant concentration -&I - heat of reaction h - film heat transfer coefficient k - reaction kinetic constant k - film mass tranfer coefftcient t - reactor length No - number of discretization points on the angular coordinate N,. - number of discretization points on the radial coordinate Nf- nyber of film mass transfer umts O$-=“-“kppJrp~$ NJ~ - number of film heat transfer
11 - effectiveness factor referred the bulk conditions
units (N,h=2Ut/aRrpFp$ PO- o-xylene partial pressure at the reactor entrance R - ideal gas constant R, - reactor radius % - reaction rate rP -,radil coordinate in the particle r- dunensionless radial coordinate in the particle (r=r-4$J T - absolute temperature U - global wall heat transfer coefficient * - interstitial velocity v- intraparticle convective velocity 2 - reactor axial coordinate I* - dimensionless reactor axial coordinate (~‘=a&) Greek Svmbols y - Anhenius ntmrber ()=E&TJ E- bulk porosity r,~- effectiveness factor referred to the catalyst surface
(5 =Itkfl&W&d 6 - angular coordinate 0 - dimensionless tcrnperatwe (S=TfTJ Am- mass intraparticle Peclet mm&r (& =vR/D& p- angularcoordinate @=cos6) p - bulk density pr- fluid density r - space time (t=Uu) dr - Thiele modulus (&&&i?=
)
SUbSClitiS b -bulk conditions d - diffusion dc - ditfusionlconvection p - catalyst particle s - particle surface w-wall
REFERENCES Glotz, G., K. Raith and W. Schoenauer, Programpackage-SUXX, Univ. Karlsruhe, Germany, 1984. Hindmarsh, A. (1974). GEAR-ordinary differential equations systems solver.Lawrence Livermore Laboratory, Report UCID 30001, Livermore. Lopes, J.C., M.M. Dins, V.G. Mata and A.E. Rodrigues (1994). Flow fieldand non-isothermal effectson difibsion,convectionand reactionin catalystparticles. Ind. Eng. Chem. Res.. accepted. Nir, A. (1977). Simnkaneous intraparticle forced convection, diffusion and reaction in a porous catalyst-II. Chem. Eng. Sci., 32,925930 Nir, A. and L. Pismen (1977). Simultaneous intraparticle forced convection, diffusion and reaction in a porous catalyst-I. Chem. Eng. Sci., 32, 35-41. Quinta Ferreira. R.M., M.M. Marques, M.F. Babo and A.E. Rodrigues (1992a). Modelling of the methane steam reforming reactor with large-pore catalysts. Chem. Eng. Sci., 47,2909-2914. Quima Ferreira, RM., A.C. Costa and A.E. Rodrigues (1992b). Dynamic behavior of fixed-bed reactors with large-pore catalysts: a bidimensional heterogeneous difhtsionkonvection model. Camp. Chem. Eng., l&, 721-751. Rodrigues, A.E. and R.M. Quinta Ferreira (1990). Effect of imraparticle convection on the steady-state behavior of fixed-bed catalytic reactors. Chem. Eng. Sci.. 45,2653-2660.