Ind. Eng. Chem. Res. 1991, 30, 909-918
909
Diffusion in Random Particle Models for Hydrodemetalation Catalysts Olivier Mace’ and James Wei* Department of Chemical Engineering, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139
Recent data showed that it is more realistic to model the structure of hydrotreating catalysts by random spheres instead of cylindrical pores and to model metal depositions by isolated growing crystallites instead of uniform layers. We develop here four progressive models to describe the structure, the average pore radius, the diffusivity, and the relation between percolation threshold and pore plugging. They are (1) the random sphere model, (2) the random needle model, (3) discrete percolation, and (4) continuous random walk. 1. Introduction Heavy crude oils usually contain significant amounts of heteroatoms, predominently sulfur, nitrogen, nickel, and vanadium, which have to be removed by catalytic hydrotreating. The deposition of nickel and vanadium sulfides on the catalyst leads to rapid deactivation. When the oil has a high metal content of several hundred parts per million (ppm), hydrodemetalation should constitute the first hydrotreating step in a guard bed, to avoid rapid deactivation of the valuable hydrodesulfurization/ hydrodenitrogenation (HDS/HDN) downstream catalyst. Most of the modeling efforts on hydrodemetalation catalysts have been directed toward developing a better understanding of this complex deactivation phenomenon. The model published by Newson (1975) features a quantitative treatment of catalyst deactivation by uniform deposition of metal sulfide layers into a cylindrical pore. At the end of the catalyst life, pore mouth plugging by accumulation of deposit layers blocks the access to the pore interior, leading to the final fast step on the deactivation curves (Tamm et al., 1981). Spry and Sawyer (1975) proposed a simple correlation for the effective diffusivity of the metal-bearing molecules in the pores:
De = DmL( T 1-
):
4
where D , is the bulk diffusivity of the metal-bearing molecule in the feed, e is the void fraction of the catalyst, T is its tortuosity factor, and r, and r,, are respectively the molecule and pore radii. Dautzenberg et al. (1978) extended Newson’s constant diffusivity, cjrlindrical pores, and uniform deposition model by expressing the “relative catalyst age” as the fraction of pore radius blocked by deposits (rini- rp)/rini,where rhi is the initial pore radius. Both hydrodemetalation (HDM) and hydrodesulfurization activities are correlated as functions of this parameter. Hughes and Mann (1978) compare the relative rates of “poisoning”, coverage of the active sites by an inactive monolayer, and “fouling”,successive adsorption of several layers until pore plugging occurs. Assuming wedge-shaped coke deposits, they distinguish between two regimes of deactivation dominated by poisoning or fouling. Rajagopalan and Luss (1979) adapted Newson’s model by describing the catalyst pellet as a sphere drilled by cylindrical nonintersecting pores of various lengths and by using eq 1 to calculate diffusion coefficients. In a model by Nitta et al. (1979), both HDM and coking were taken into account. The hypothesis of a reversible coking reaction results in the prediction of decreased levels of coke as metal ‘Current address: Mobil Refinery,Notre Dame de Gravenchon,
B.P.#2, 76330 Notre Dame de Gravenchon, France.
sulfide deposition increases. Tamm et al. (1981), using electron microprobe techniques, observed a metal deposition gradient in the catalyst pellet. They emphasize the importance of diffusional limitations by using a “distribution parameter”, the ratio of the amount of deposit averaged over the whole pellet to its maximum value near the edge. A low distribution parameter can be interpreted as an underutilization of the pellet center. More recently, investigators have tried to take into account the complexity of the porous geometry: Toulhoat et al. (1987)used a biconical pore shaped model to simulate micropores as well as macropores. D. M. Smith (1986) calculated diffusion coefficients in the case of periodically constricted pores. Investigating the effect of pore radius definition, he distinguished between the radius measured by mercury porosimetry and the radius calculated as twice the ratio of pore volume to surface area. Smith and Wei (1990) observed that metal sulfide deposition is not uniform but forms a discrete crystallite pattern and described both catalyst and deposit by using random distributions of overlapping spheres. Melkote and Jensen (1989) compared the effects of uniform and discrete deposition on both active sites coverage and pore plugging and used percolation theory concepts to model the decrease of coordination number of a pore to its neighbors. Recent electron microscope pictures (Toulhoat et al., 1987; Smith, 1988) showed that the microporous geometry is a complex conglomerate of grains of diverse shapes, particularly needles, plates, and spheres, and is not accurately described by the cylindrical model with uniform deposition. A representation of both catalyst and deposit by statistical distributions of spherical particles (Smith, 1988) shows better agreement with the EM pictures and offers a better quantitative understanding of diffusion and deposition. 2. Diffusion and Reaction with a Random Spheres Model The random spheres model was described by Weissberg (1963) as a statistical medium of randomly located, freely overlapping spheres of uniform size. The model is entirely defined by two parameters, the radius a of the spheres and the volumetric density n of their centers. The spherical symmetry allows computation of the void fraction = e-(4/3)*na9
(2)
and of the surface area
z = c(4rna2)
(3) This computation can be extended to a distribution of sphere radii and densities. Van Eekelen (1973) reviewed the theory of random overlapping spheres and proposed a calculation of the mean pore size in this medium. Smith
0 1991 American Chemical Society 0888-5885f 91/263Q-Q9Q9$Q2.5O/Q
910 Ind. Eng. Chem. Res., Vol. 30,No. 5, 1991 I
Table I. Input File for the Pellet Model, with Base Case Values ~~~
catalyst characteristics R, pellet radius SM,initial surface area Z,initial surface area tM, initial void fraction pe, catalyst density Pd, deposit density leading to al, random sphere radius nl,random spheres density coking parameters hf&, final amount of coke coking time p& coke density oil characteristics po, oil density CR,reactant concentration in the feed d,, molecule diameter characteristicrates parameters k, intrinsic rate constant r, activity ratio I),,,,,, initial porous diffusivity run length
1 mm 176 m2/g 262 m2/cms 0.64 1.49 g/cm3 5.82 g/cms (nickel case) 32.7 A 3.04 X 1@ m-s 0.2 wt % (hardly any coke) 1 day 1.5 g/cms
.
0.9 g/cms
500 PPm 30 A (asphaltene case) 5 X 10-lo kg/(m2 s) 0.1 5 x IC7 cm2/s 1 year
and Wei (1990)adapted this concept to the description of an aged HDM catalyst. The fresh CoMo/A1203catalyst is described (see Figure la) by a distribution of spheres of radius al and density nl (both values derived from the experimental values of the initial cini and Zini by using eqs 2 and 3). On these computer-generated drawings, all spheres with centers lying in a sample of thickness 80 A are represented. The parameter values are given in Table I. No account is taken of cross-cuttings, so that all spheres are represented as disks of the same radius, but an impression of depth is given by plotting lighter shades for the spheres on top of the sample and darker shades for the spheres at the bottom. Smith (1988) found experimentally that, for a given catalyst, there is a fixed density of nucleation sites n2, which is related to the acidity of the catalyst, and the crystallite radius u2 grows with time as the catalyst ages. Thus, the deposit is modeled as crystallites, by a second distribution of spheres, of radius a2and density n2,given as darker spheres in parts b and c of Figure 1. The void fraction and surface area of the aged catalyst become
= e-t4/3)r(wla+w~9 Z = t[4?r(nlul2+ n$22)]
(4) (5)
To confirm the validity of the model, we study its rigidity and connectivity when c increases. The result is that, for void fractions lower than 0.65, there is an infinitely connected solid cluster leaving only about 5% unconnected spheres. But for a void fraction of 0.7, the probability of having an infinite connection is only 0.5 and the model does not realistically describe a solid catalyst. This model of double random spheres distribution is incorporated in the geometry of a spherical catalyst pellet. We designate by R the radius of the pellet and by r the radial position from 0 (center) to R (edge), and the species conservation balance in the catalyst is obtained by identifying the incoming flux to the first-order reaction rate
with boundary conditions
I
k
A f
l&
Figure 1. Random spheres model; (a, top) fresh catalysG (b, middle) 25 wt 7 ' 0 metal sulfide; (c, bottom) 100 wt ?% metal sulfide.
Ind. Eng. Chem. Res., Vol. 30, No. 5, 1991 911 The local rate of deposition follows the equation aM,/at = skc (9) This system of partial differential equations can be solved by an explicit finite differences scheme (Limbach and Wei, 1988),with the assumption of slow evolution of the catalyst geometry compared to the rate of diffusion and reaction. For each time step, local void fractions are derived from M , and eq 4, and the local values of the surface area are then calculated by using the random spheres model (eq 5 ) . As we have three parameters (eini, Zfi, and e) and three equations (eqs 2, 3, and 4) for four unknowns (al, n,,a,, and n2),the choice of a degree of freedom is made by assuming that the crystallites have constant density n2and variable radius. The base case in Smith and Wei (1990) has a deposit density that is 10 times smaller than the alumina particle density, so our choice of adjustable parameter will be n, = nl/lO. Figure l illustrates the simulated progressive aging of a catalyst sample of thickness 80 A, following this hypothesis and the parameters values given in Table I. The local diffusivities are calculated by using eq 1,in which the pore radius rp is assumed to be the equivalent radius of an infinitely long cylinder of inner volume t and inner surface Z rp = 2 e / 2 (10) When this pore radius becomes equal to the molecular radius rm, it is assumed that the catalyst pellet becomes completely inactive. The tortuosity is assumed to be constant with the value of 4 (Satterfield, 1970) although it should increase with deposition. The model contains a treatment of the coking of the catalyst, where the coke level reaches an equilibrium in a few weeks. Coke is assumed to be inert and to deposit uniformly on the whole catalyst surface, resulting only in an increase of the catalyst microsphere radius al-and in the subsequent decrease of diffusivity in eq 1,void fraction in eq 2, and surface area in eq 3. In the model, coke does not decrease the intrinsic activity of the surface that it covers. On the contrary, the metal-deposit surface is assumed to be less active, although not completely inactive, than the fresh catalyst surface. Defining r to be the ratio of the deposit crystallite rate constant to the catalyst intrinsic rate constant k, we correct the surface area I: = + 2, (sum of the uncovered alumina surface and of the crystallites surface) to take into account the difference of activity. The surface area corrected for activity is 8' = I:~ rz2 (11) In this study, we use the letter I: for the surface area per unit volume and S for the surface area per unit weight. Pellet deactivation is simulated assuming a constant temperature run, rather than the industrial case of constant conversion. k is therefore kept constant. The required input is given in Table I; the indicated values are used as a base case, to simulate a hypothetical catalyst able to retain more than 100 w t % metal sulfide without complete deactivation. Figure 2 illustrates the evolution of surface areas averaged over the pellet volume. Sl decreases due to coverage by crystallites, while S2increases sharply at the beginning before slowing down as the crystallite radius increases. The deposition gradient of Figure 3 corresponds to a final Tamm distribution parameter of 0.79 and is the direct result of the concentration gradient of porphyrin in solution in Figure 4. Maximal deposition occurs at the edge of the catalyst, as we do not take into account the existence
---------Sa
//--
/'
/
/'
0
100
400
300
200
DAYS ON STREAM Figure 2. Aging of a catalyst pellet; surface area of the alumina (Sl), of the crystallites (SJ, total (S1+ Sz), and corrected for activity (S,
+ rs2).
n w
fi L,
-
100
-
-1
2
..
4
d
e
i
E
_ _ _ _-.---
t=e
_-__---___________-------
50
t-3
-
_---*----:
__----___----
________________------------
every three months
t=o
0
+
I
,
,
,
,
400 h
E
a
_----
a
v
_----
__--__---
I t=12
0 '
0.0
"
"
I
0.5
'
"
'
1.o
RELATIVE RADIAL POSITION ( R = l ) Figure 4. Aging of the catalyst pellet; evolution of the concentration
gradient.
912 Ind. Eng. Chem.
Res., Vol. 30,No. 5, 1991
Figure 5. Needle, of radius a and half-length b, featuring spherical caps at the ends.
of reaction intermediates described by Agrawal and Wei (1984)and Ware and Wei (1985). The Thiele modulus, in this case
varies a t the pellet edge from 1.7 a t the beginning to 3 a t the end of the run, which qualifies this base case as less diffusion-limited than typical industrial values. The model can be adapted to a bimodal catalyst by using the macrospheres/microspheresmethod (Agrawal, 1980). Characteristic length scales are a millimeter for the catalyst pellet (Agrawal's macrosphere), a micromillimeter for the microspheres separated by the macropores, and a nanometer for the random spheres describing the microporous region into each microsphere. Reaction in the macropores is neglected due to the small surface area involved.
3. The Random Needles Model On the basis of transmission electron microscopy pictures, the alumina catalyst support has been reported as a conglomerate of long needles (Nourbakhsh et al., 1990) or plates. On the basis of the random spheres model, we develop a model dealing with randomly located and oriented overlapping needles of uniform size. The fresh catalyst is now represented by a distribution of these needles, while the deposit crystallites remain as spheres. Defining n to be the volumetric density of needle centers, a the needle radius, and b its half-length (or alternatively a the ratio b / a ) , we have the needle void fraction and surface area ~-~4/s~Mslar+la~2~+~l/2~la~~a-r~/~~l~~+~l/2~~~af/l~arl~~P+1~J~l
(13)
d z
Z, = en(4ma2)[a2+ - a{=] (14) to be compared to eq 2 and 3. As given in the Appendix, the computations of e, and Z, require a slight change in the needle shape: I t can be described as the shape obtained by cutting a sphere of radius b by a cylinder of radius a (see Figure 5). With this convention, the sphere model itself is included as the case a = 1. Figure 6a has the same initial void fraction and the same initial surface area as Figure la. The needle radius is 21.9 A,the half-
Figure 6. Random needles model; (a, top) fresh catalyst; (b, middle) 25 w t % metal sulfide; (c, bottom) 100 w t % metal sulfide (c). The deposits are spheres.
Ind. Eng. Chem. Res., Vol. 30, No. 5, 1991 913 1.o
+ w
0.5
0.0 0
2
1
n.v Figure 7. Rate of overlapping ( 6 versus m) for spheres and needles.
length is 219 8, (so a = lo), and the density of needles n = 6.5 X loB m-3. The increasing deposition of spherical crystallites in the random needles model is shown in Figure 6. These drawings are generated with the same conventions as in Figure 1. The size of the crystallites is 43.6 A in Figures l b and 6b and 73.8 8, in Figures ICand 6c. The relationship between the volume of particles and the void fraction is given in Figure 7. The product nu represents the total volume of particles, overlappings counted several times, per unit volume. Figure 7 shows that the rate of overlapping of spheres (a = 1) and infi) the same and that the least nitely long needles (a = ~ 1 is overlapping occurs for an ellipsoid-likeshape obtained for a N 1.35. As the rates of overlapping for spheres and for needles of CY = 10 are very similar, the void fraction will be the same if we compare the two models at same density n and same particle volume u. However, the surface of a needle is significantly larger than that of a sphere of the same volume (72% larger if a = 10). Consequently, the main difference between the two models is the larger needle surface area at the same void fraction and particle volume or, alternatively, the need of less numerous but bigger needles than spheres to describe a given e and 8 . This difference in size can be seen in Figures l a and 6a. For the same initial void fraction 0.64 and surface area 2.62 X 108 m2/m3,an individual needle of a = 10 has a volume about 5 times larger than the sphere. It is worthwhile to point out that a random coins model can be established. The sphere of Figure 5 can be cut by two horizontal planes at equal distance a from its center, and there remains a median part shaped like a coin. The void fraction and surface area of the random coins model can also be computed by using the probabilistic concepts described in the Appendix. 4. Diffusion by Percolation in a Discrete Cubic Lattice Next we investigate the relative blocking effectiveness of a few large clusters of deposits vs many smaller clusters with the same total volume and the networking of clusters to achieve blocking. These concepts are investigated with the method of random walks in a defective lattice in the framework of percolation theory. First we will deal with the simpler discrete cubic lattice and then with a continuum random sphere model. In a simple cubic lattice, random sites are progressively covered and made unavailable for diffusion, either one by
Figure 8. Compact clusters of 1, 7, and 33 sites.
one or by clusters. The clusters, as shown in Figure 8, are of 1 (simple site percolation), 7, and 33 sites, and their compact shape is used to simulate spheres of radius 2 or 3 times the deposit atomic diameter d . For a given size of clusters and the proportion of remaining sites (or void fraction), a computer program generates several “random spheres” lattices and then simulates a large number of random walks in each lattice. The random walk in this case is defined by the most intuitive rules (Stauffer, 1985): At the beginning of each time step, the diffusing molecule sits on one of the still available sites-the void fractionand one of the six directions in space is randomly selected. If the neighboring site in this direction is still vacant, the molecules moves there. If the neighboring site is covered by solid, the molecule stays where it is. After a random walk of T time steps, the program collects the traveled distance R T and initiates a new random walk. The relative diffusion coefficient D, = D,/D,, for this size of clusters and void fraction, can then be obtained by averaging the squares of traveled distances over all different random walks and lattices (RT2)
= D,d2T
(15)
where d is the unit bond length (note that T , the number of time steps, is a dimensionless number). However, eq 15 is only valid above the percolation threshold, i.e., for void fractions-or proportion of available sites-greater than a limit ep’ Below this value ep, the void fraction is no longer infinitely connected but is reduced to isolated void chambers, and the proportion between ( R T 2 ) and T is no longer valid. Therefore, we interpret the final rapid step of the HDM deactivation curve as being closely related to a percolation phenomenon, rather than being the actual plugging of a pore mouth. Figure 9 illustrates the concept of percolation threshold for solid clusters of 1 site (simple site percolation): Equation 15 becomes log ( R T ’ ) = log T + log (d2D,) (16) so a straight line of slope 1 and intercept log (d2D,)at T = 1 is expected above the percolation threshold. Below the threshold (case 1 - c = 0.7), any random walk starting in a given isolated void chamber will end in the same chamber, explaining the almost horizontal line for 1- e = 0.7. The traveled distance in this case becomes a characteristic length of the average void chamber-it is actually
914 Ind. Eng. Chem. Res., Vol. 30, No. 5, 1991 1000000 I-
' ' ' """ ' ' ' ""'I 1--E = 0.0 (0) (solid fracuon) 0.4 (x)
'
'
10000
0
9
4
+
+
+
w
loo[
1 1000
100000
10000
'"'9
0
-
0
-
1 1000000
TIME STEPS Figure 9. Average square distance traveled in random walk versus number of time steps, above and below the percolation threshold. Clusters of 1 site.
Figure 11. Random walk in the continuous random spheres model.
-
the average distance between two random points in the chamber-which is a constant when T increases. The average sequence distance vs time steps can be described by a straight line of slope y, so that we can write
lation threshold, where relative diffusivity goes to zero, as the asymptote l/D, m. Concerning the three compact clusters, the result is that bigger solid clusters are more favorable for diffusivity than small ones, by allowing percolation to occur at a higher solid fraction. This result can be intuitively understood: that, for the same amount of solid in the lattice, diffusion is easier if there is a well-defined segregation between big solid clusters and big void chambers in comparison with a more intimate mixing of solid and void. Analogously, the random oriented rod case is less favorable to diffusion than the 33-site sphere cases, as extended rods introduce more blockage and steric hindrance than compact shapes. The similarity between random oriented rod of 7 and of 33 sites in a row can be understood, especially at high solid fractions, as a result of overlapping and touching, causing the fusion of rods several times longer than the unit size of 7. The partition of the rods among the three axial directions is necessary to their higher diffusional hindrance, in comparison with the case of vertical rods, which leads to a more open structure. A practical application of this set of results is that the qualitative trend-better diffusion in the case of bigger spheres-is certainly valid for crystallites growth: For a given amount of metal sulfide in the catalyst, the diffusivity will be better retained if the crystallites are less numerous and bigger. Most importantly, the total catalyst deactivation by percolation will occur at a bigger amount of deposit, so any technique of increasing the size and decreasing the number of the crystallites would improve the metal capacity of an hydrodemetalation catalyst.
(R& a TY (17) with y decreasing continuously from 1above the percolation threshold to 0 below the threshold, which is consistent with Stauffer (1985). The approximated value of cp between 0.3 and 0.35 is also consistent with more precise computations of the percolation threshold for simple site percolation in a cubic lattice (~0.3117). In Figure 10, diffusion coefficients are evaluated by the random walk method for different cluster sizes and different values of the solid fraction 1- e. We also indicate the results of some additional simulations where pseudospherical solid clusters are replaced by elongated clusters or "rods", randomly located and oriented along the three axes. The use of l/D, instead of D,allows for each cluster shape and size an immediate identification of the perco-
5. Diffusion by Random Walks in a Continuum Random Spheres Model The preceding percolation study demonstrates qualitative results, but its quantitative value is limited by the cubic lattice approximation. As a next step, we adapt the random walk concept to the continuous random spheres model. The rules of random walk in ths context are illustrated in Figure 11. As long as the molecule does not hit a solid sphere, it travels the distance X per time step before changing direction randomly, with equiprobability of all space directions (point 1). If the molecule encounters a random sphere before having traveled A, the rule is that it adsorbs on the sphere surface until the end of the time step. An alternative to taking adsorption into account would have been to model an elastic collision with im-
VOID FRACTION 1.o
I 30
-
(E)
0.5
"
"
"
"
'
0.0
1 0
A apherea of 1 sib x spheres of 7 sites 0
a
apherea of 33 sibs
A
OrOdrOf7dtOS 0 rods of 33 sftea
20
-
0 V0-d
X
V
rodr Of 7 sites
0
' a
b
t
--
lo
A
t
0.0
V
a
9( 0 v
0.5
:
i 1.0
SOLID FRACTION (1-e) Figure 10. Diffusivity and percolationthreshold for diverse cluster sizes and shapes.
Ind. Eng. Chem. Res., Vol. 30, No. 5, 1991 915
B c=0.5. h=3A
0
E-1, rm=6.85A
0 r=0.5, h=5A
o
€50.5, r,=S.BSA
0 ~ 0 . 5 h=50A .
-
z 4
(rm-5.S5A in all cases) m
A
0 e-0.2, r,=30A
(h=SOA
. ' A 0
.
in all cares)
100000
A
10000 W
3 W
1
1000
0.0 1
10
100
1000
10000
100000
TIME STEPS Figure 12. Relative diffusivity, from a local to a global value.
mediate reflection. In the present method, the adsorbed molecule, having no memory where it came from, desorbs in a direction depending only on the surface orientation. The angle t9 between the desorption direction and the normal to the sphere at point 2 follows a cosine probability distribution, so that the most probable direction is the normal to the sphere (cos (0) = 1) and the least probable are tangent to the sphere (cos (ir/2) = 0). As in a percolation lattice, values of the diffusivity D, = De/Dmcan be obtained by averaging traveled distances squared for a given number of time steps T and a given random spheres medium (a; n) (RT2)
= D,X2T
(18)
The coefficient A is analogous to the mean free path in the kinetic theory of gases, where molecules travel a significant distance before hitting each other. When X a, the molecule only collides with solid spheres as in Knudsen diffusion. Here A is a parameter necessary to the adaptation of percolation theory concepts to a continuum, where a diffusing molecule may collide with either another molecule or the wall. Kirkaldy and Young (1987) have reported a gaslike diffusion in the case of liquid metals, with a mean free path proportional to the linear thermal expansion coefficient multiplied by the temperature and by the molecule diameter. Values of this mean free path are reasonably small-about 1/10th of the molecule radius-as can be expected in a closely packed state. However, there is no direct similarity between liquid metals and heavy crude oil, and this random walk technique will only be taken as a qualitative approach of liquid-phase diffusion in a HDM catalyst. It can nevertheless prove very useful to evaluate the properties of the random spheres porous medium, generated now without any approximation. As a starting point, we consider values of X much smaller than the random sphere radius of 32.7 A, the order of magnitude of 1A seeming a priori appropriate. However, qualitative trends are more obvious and less costly to demonstrate with larger values of A, which explains our repeated use of X = 50 A. The number of time steps needed to get an accurate value of the relative diffusivity D, increases when X decreases, as shown in Figure 12: When the number of time steps is small, so that the average traveled distance is of the order of magnitude of the sphere radius, the random walk is confined to a local
-
1
10
100
1000
TIME STEPS Figure 13. Average square distance versus number of time steps, above and below the percolation threshold.
region-mainly, the void chamber around the starting point-and the molecule on1 "sees" a locally high void fraction. The case of X = 50 in Figure 12 confirms that the lower global value of the diffusion coefficient exists and can be reached in less than 30 steps. As the traveled distance is proportional to XT1/2, the number of time steps necessary to reach the global diffusivity varies as 1/X2, and T = 10OOO for A = 3 8,compared to the much smaller value of T = 30 for X = 50 A. In the rest of this study, all diffusion coefficients are evaluated by using numbers of time steps TA equivalent-by the 1 / X 2 rule-to T = loo00 for X = 3 A. This results probably in a slight (less than 10%) overevaluation of D,. Figure 13 is the equivalent for a continuous medium of Figure 9, and the straight lines of slope 1 and 0 can be identified as respectively above and below the percolation threshold. In the lattice case where the molecule size is basically 1 unit site, the percolation threshold is a purely geometrical property of the medium. In the continuum case, the percolation threshold is highly dependent of the molecule size. For the exact same random spheres medium of spheres radius 32.7 A and void fraction 0.2, the percolation threshold has already been reached for r, = 30 A (0) but not yet for r, = 5.85 A (A). The first case is a very big asphaltenic molecule, while the second is the smaller etioporphyrin. Asphaltenes and porphyrins are the two main types of metal-bearing compounds in crude oil (Quann et al., 1988). The concept of percolation threshold is usable in a continuous medium, but it is now a function also of the ratio of the size of the diffusing molecule to the random sphere radius. Therefore, the void fraction of a porous medium considered as percolated for asphaltenes does not have to be composed of isolated holes: The holes may still be connected and the connections, although too narrow for asphaltenes, may still be large enough for porphyrins, as is the case for r, = 5.85 A in Figure 13. Figure 14 illustrates the difference in diffusivity decrease between the two molecular sizes. The random spheres media are generated with constant sphere radius 32.7 A and variable density. The percolation threshold for the 30-A asphaltene case occurs at a void fraction of about 0.35, while there is no evidence of a percolation threshold until very low void fractions for the etioporphyrin case. These results can be compared to those predicted by the Spry and Sawyer formula, in which we assume that pore
1
1.09,
I
I
I
,
I
,
I
I
P h=3A
0 h=5A
E2
‘\.
VI
2cr,
0
0 h-50A
(rm=5.85A)
A
‘9
8
‘\\\
0.5
-
0 h=lOA
‘9,.
-
-
A
“\,
SOLID FRACTION ( 1-E) Figure 14. Relative diffusivity versus solid fraction for different molecular radii.
VOID FRACTION
(E)
0.0
0.5
1.0
P rm=5.85A x r,=30A
( h = U in both cases)
A A
0.5
5 w
A
x
E
X
X
0.0 0.0
I
I
I
0.5
I
I
1 .o
SOLID FRACTION (1-E) Figure 15. Relative diffusivity versus solid fraction for different values of the parameter A.
plugging-or percolation-occurs when the two radii rm and rp become equal. By using eq 10 for the pore radius, pore plugging according to this criterion happens for t = 0.025 in the etioporphyrin case and about 0.5 in the case of rm = 30 A. The first value gives an explanation to our failure in reaching a percolation threshold for etioporphyrins, as percolation probably occurs only for very small remaining void fractions. In the asphaltene case, however, the value predicted by eqs 1 and 10 is significantly lower than the result of random walk experiments. The value of the percolation threshold, being directly related to the catalyst metal capacity, is a crucial data in the evaluation of the catalyst quality and in the prediction of the final deactivation step, and a difference of 0.2 in the ultimate void fraction is more than significant. Our random walk method uses the exact description of the random spheres model, whereas eq 10 merely translates the model
Figure 16. Conditions of coverage of
M by the needle (C; v).
into its cylindrical equivalent. Finally, Figure 15 presents the effect of A on the relative diffusivity, for the etioporphyrin case. The absolute diffusivity De is highly dependent on A itself, as the bulk diffusivities D, between the cases A = 50 A and A = 3 A are in the ratio 502/32= 278. For A S 10 A, the curves follow more or less the diagonal D, = t, with the case of A = 3 A being certainly overevaluated, possibly due to the phenomenon described in Figure 12. The case of A = 50 A gives significantly lower values of D, and seems to follow a form D, = e/?, where the tortuosity factor 7 increases with the solid fraction. However, we have to consider that, at large values of A, our simulation switches from normal “gaslike”diffusion to Knudsen diffusion. If A becomes very large, the characteristic unit length of the random walk is not A anymore but is the average distance between two spheres separated by a void. 6. Conclusions Most catalyst supports are porous material formed by cementing and conglomerating small spheres, needles, and plates precipitated from solution and are much more realistically represented by the random spheres or the random needles models. These random models are much more sophisticated and mathematically difficult than the pore models, and require stochastic theories and extensive computation. The void spaces available for guest molecules to diffuse have continuously changing width and connectivity and do not possess the simple parameter of pore diameter. We have solved the problems of surface areas and void volumes, as well as diffusivity and pore plugging, in these random models. The four models of porous solid and deposit morphology give a progressive approach to a more realistic description of metal deposition in HDM. For the random spheres model, the evolution of surface area and deposition profile with time is computed by using a pellet model. Diffusion in a discrete lattice is studied in the framework of percolation theory, to investigate the influence of removing large clusters of sites in comparison with simple site percolation. For the same number of sites removed, diffusion is faster for large cluster percolation than for simple site percolation. Spherical clusters are also more open to diffusion than rodlike clusters. The percolation threshold is related to pore plugging and the final rapid step of the deactivation curve. Diffusion in the continuous random spheres model is investigated by random walk with a parameter A, which is interpreted as a mean free path in gaseous diffusion. When A is much smaller than a characteristic pore length p, the diffusing molecule does not encounter a solid surface for small number of time steps, so that the diffusivity is close to the bulk diffusivity. At large number of time steps, the molecule encounters solid spheres and the diffusivity decreases to its average value in the medium. Below a percolation threshold, the void space becomes a number of isolated chambers, so that the average traveled distance approaches an upper limit, which is constrained by the characteristic chamber size. This percolation threshold is a function of the ratio of the molecular radius to the sphere radius, and a large molecule would cease to diffuse at a void fraction that would not deter a smaller molecule.
Ind. Eng. Chem. Res., Vol. 30, No. 5, 1991 917
Figure 17. “Right directions” for v, as defined by the solid angle
a.
Acknowledgment We are indebted to Dr. Ernst Gail for his critical contributions. The random walk simulations in the continuous random spheres medium were done on the MITSF CRAY-2 supercomputer. Olivier Mac6 benefits from a grant from the Ecole Nationale Sup6rieure du P6trole et des Moteurs (France). Financial support is received from Mobil Research and Development and from Chevron Research Company. Nomenclature a = random particle radius (m) b = needle half-length (m) C = reactant concentration (kgmet/kgoil) CR = reactant concentration in the feed (kgm,,/kgoil) d = unit bond length in the cubic lattice (m) d , = molecular diameter (m) De = effective diffusivity (m2/s) D, = relative diffusivity D, = bulk diffusivity (m2/s) k = intrinsic rate constant (kgoil/(m28)) M,= local metal amount (kgmet/kg,,,) Mcok= equilibrium coke amount (kgcok/kgmt) n = random particle density (m-9 r = radial position in the pellet (m) rm = molecular radius (m) = pore radius (m) = pellet radius (m) RT = traveled distance after T time steps (m) s = particle surface (m2) S = surface area (m2/kg,,) t = time (s) T = number of time steps u = particle volume (m3)
2
Creek Letters CY = ratio b/a in the random needles model y = exponent of T in eq 17 r = activity factor t = void fraction tp = percolation threshold 0 = desorption angle in Figure 11 (rad) = random walk parameter (m) pc = catalyst density (kgmt/m3) pwk = coke density (kgcok/m3) Pd = deposit density (kgdep/m3) po = oil density (kgoil/m3) B = surface area (m2/m3) 7 = tortuosity factor 3 = Thiele modulus Subscripts 1 = fresh catalyst 2 = crystallites s = random spheres model n = random needles model ini = initial
Appendix This appendix gives the calculation of void fraction and surface area for the random needles model. The position of a needle of radius a and half-length b (see Figure 6) is completely determined by the three coordinates of its center C and the direction v of its axis. If c, is the void fraction for a density n of needla of unit surface s, the total surface area-overlapping not taken into account-is m and the remaining surface area after overlapping is I;, = c,(ns). The value of I;, can therefore be derived directly from E,. The void fraction E , is the probability that a random point M belongs to the void, i.e., is not covered by any needle. Figure 16 gives the two possibilities of coverage by a needle: If the distance between M and C is smaller than a, M is always covered. If it is between a and b, M is covered only if v is in the “right direction” in Figure 16. Finally, if the distance between M and C is greater than b, M cannot be covered by the needle of center C. Conversely, the point M is not covered by any needle if there is not needle of center closer than a (event of probability P1)and, for every center between a and b, the associated vector v is not in the “right direction”, event of probability p,. The void fraction of the random needles model is therefore c, = PIP2 (19) P1,the probability of not having a center closer than a for a volumetric density n of centers, is well-known from the random spheres theory (Weissberg, 1963)
p1 = e-(4/3)mw3 (20) To calculate P2,let C be a center at distance r of M, with a Ir I b. The vector v associated with C is not in the right direction if v is not in the solid angle @ in Figure 17. As @ = 2741 - cos a),the probabiilty that v is not in the right direction is 1- @/2a,27r being the solid angle of the complete half-sphere. This probability is exactly cos CY = (r2- a2)ll2/r. The needle shape of Figure 6 is a condition for the validity of this expression at r = b. Using a classical technique of differential calculus, we consider the crown [r; r + dr]. The probability that in this crown there is no needle center associated with a vector in the right direction is
P(,)= (cos cu,)Nr
(21)
where N , = n(4?rr2)dr is the number of centers in this crown. So p(,)= (cos a,)4~nr% = e4*nRn((r2-a2)1/*/r)dr (22) and integrating over r
p 2 = fip(,, = eSnb4rn*((a*)l/2/r)* ria
(23)
We finally obtain c, = PIP2 = e - 4 * n I ( a 5 / 3 ) + S n b ~ ( r / ( r ~ 2 ) ’ ~ 2 ) d r l (24) directly leading to eq 13. Literature Cited Agrawal, R. Kinetics and Diffusion in Hydrodemetallation of Nickel and Vanadium Porphyrins. Sc.D. Thesis, MIT,Cambridge, MA 1980.
Agrawal, R.; Wei, J. Hydrodemetalation of Nickel and Vanadium Porphyrins (I). Ind.Eng. Chem. Process Des. Dev.1984,23 (3), 505-5 14. Dautzenberg, F. M.; van Klingen, J.; Pronk, K. M. A.; Sie, S. T.; Wijffels, J. B. Catalyst Deactivation through Pore Mouth Plugging
Ind. E n g . Chem. Res. 1991,30,918-922
918
during Residue Desulfurization. ACS Symp. Ser. 1978, 65, 255-267. Hughes, C. C.; Mann, R. Interpretation of Catalytic Deactivation by Fouling from Interactions of Pore Structure and Foulant Deposit Geometries. ACS Symp. Ser. 1978,65,201-213. Kirkaldy, J. S.;Young, D. J. Diffusion in the Condewed State; The Institute of Metals: London, 1987. Limbach, K. W.; Wei, J. Effect of Nonuniform Activity on Hydrodemetallation Catalyst. AZChE J. 1988,34 (2),305-313. Melkote, R. R.; Jensen, K. F. Models for Catalytic Pore Plugging: Application to Hydrodemetallation. Chem. Eng. Sci. 1989,44(3), 649-663. Newson, E. Catalyst Deactivation Due to Pore-Pluggingby Reaction Products. Ind. Eng. Chem. Process Des. Dev. 1975,14(l),27-33. Nitta, H.; Takatsuka, T.; Kodama, S.; Yokoyama, T. Deactivation Model for Residual Hydrodesulfurization Catalysts. AIChE Annual Meeting, 86th Houston, 1979,paper 34E. Nourbakhsh, N.; Smith, B. J.; Webster, I. A.; Wei, J.; Tsotsis, T. T. Metal Deposition in Porous Anodic Alumina Films under Hydrotreating Conditions. J. Catal. 1990,in press. Quann, R. J.; Ware, R. A.; Hung, C. W.; Wei, J. Catalytic Hydrodemetallation of Petroleum. Adv. Chem. Eng. 1988, 14,95-259. Rajagopalan, K.;Luss, D. Influence of Catalyst pore size on Demetallation Rate. Ind. Eng. Chem. Process Des. Dev. 1979,lB(3), 459-465. Satterfield, C. N.Mass Transfer in Heterogeneous Catalysis; MIT Press: Cambridge, MA, 1970.
Smith, B. J. Deactivation in Catalytic Hydrodemetallation. Sc.D. Thesis, MIT, Cambridge, MA, 1988. Smith, D. M. Restricted Diffusion through Pores with Periodic Constrictions. AZChE J. 1986,32 (6),1039-1042. Smith, B. J.;Wei, J. Deactivation in Catalytic Hydrodemetallation (111). J. Catal. 1990,in press. Spry, J. C.; Sawyer, W. H. Configurational Diffusion Effects in Catalytic Demetallization of Petroleum Feedstocks. 68th AIChE Annual Meeting, Los Angeles, 1975,paper 30C. Stauffer, D. Introduction to Percolation Theory; Taylor & Francis: London, 1985. Tamm, P. W.; Harnsberger, H. F.; Bridge, A. G. Effect of Feed Metals on Catalyst Aging in Hydroprocessing Residuum. Znd. Eng. Chem. Process Des. Deu. 1981,20(2),262-273. Toulhoat, H.; Plumail, J. C.; Houpert, C.; Szymanski, R.; Bourseau, P.; Muratet, G. Modelling RDM Catalysts Deactivation by Metal Sulfide Deposits. Prepr.-Am. Chem. Soc., Diu. Pet. Chem. 1987, 32 (2),463. van Eekelen, H. A. M. The Random-Spheres Model for Porous Materials. J. Catal. 1973,29,75-82. Ware, R. A.; Wei, J. Catalytic Hydrodemetallation of Nickel Porphyrins (I). J. Catal 1985,93,100-121. Weissberg, H.L.J. Appl. Phys. 1963,34 (2),2636.
Received f o r review May 1, 1990 Revised manuscript received July 26, 1990 Accepted August 7, 1990
Free-Surface Effects in Torsional Parallel-Plate Rheometry Robert
W.G.Shipman,? Morton M. Denn,* and Roland Keuningst
Center for Advanced Materials, Lawrence Berkeley Laboratory, and Department of Chemical Engineering, University of California a t Berkeley, Berkeley, California 94720
Free-surface distortions caused by the combined effects of gravity and surface tension have little effect on the measured torque in torsional parallel-plate rheometry, but they can affect the measured normal stresses for dilute polymer solutions. A static free surface correction to the normal force is in good agreement with finite-element calculations and can be used to correct experimental measurements. Introduction The measurement of rheological properties in shear is normally carried out in either a parallel-plate or coneand-plate geometry. The conventional analysis of torsional shear rheometers (Walters, 1975) is based on a flow field that is unaffected by edge effects at the liquid/air interface (although explicit assumptions must be made about the shape of the interface in order to determine the level of the isotropic pressure, which is required for the measurement of viscoelastic normal stresses). The shape of the interface is known to have an effect on measured rheological properties; see Kalika and co-workers (1989), for example, for an illustration of the pronounced effect of edge conditions on the measured viscosity. We explore here, through the use of finite-element simulation and asymptotic analysis, the effect of free-surface distortions caused by gravity, surface tension, and inertia on measured rheological properties in parallel-plate torsional rheometry. The effect of inertia on the free surface is negligible under normal conditions for parallel-plate rheometry. The effects of surface tension and gravity are sufficiently large that they need to be taken into account
* To whom correspondence should be addressed at the Department of Chemical Engineering. +Present address: 3M, St. Paul, MN 55144. Present address: UniG de MBcanique AppliquBe, Universit6 Catholique de Louvain, B-1348 Louvain-la-Neuve, Belgium. f
for the measurement of normal stresses in dilute polymer solutions, but they will have little effect on concentrated solutions and melts. These effects of surface tension and gravity are adequately described by a static analysis of the interface that can be used to correct experimental normal stress measurements.
Parallel-Plate Rheometer The parallel-plate rheometer is shown schematically in Figure 1. In the usual mode of operation, one plate is moved at a fixed rate, either steady or oscillatory, and the shear and normal stresses are determined from torque and thrust measurements on the stationary plate. In the absence of inertial terms, and assuming infinite disks, the flow field for an arbitrary liquid is given by ue =
rQ{, u, = u, = 0
(1)
where uo, u,, and u, are the velocity components in a cylindrical coordinate system, Q is the rotation rate, and { = 1 - z / H , where H is the plate separation. The torque measured on the top plate in steady shear for a Newtonian fluid with viscosity 9 is
We assume for discussion that the polymer solutions are described by the Oldroyd-B fluid, which is a continuum
OS8S-5885/91/2630-0918$02.50/0 0 1991 American Chemical Society
