3154
Ind. Eng. Chem. Res. 1997, 36, 3154-3162
Studies of Transport of Asphaltenes through Membranes Using Hindered Diffusion Theories for Spheres and Spheroids V. S. Ravi-Kumar, Theodore T. Tsotsis,* and Muhammad Sahimi Department of Chemical Engineering, University of Southern California, Los Angeles, California 90089-1211
The results of ongoing efforts by this group to model the transport of asphaltene molecules through model membranes are presented. A model is described which aims to capture the effect of the polydisperse nature of asphaltene molecules on their transport properties. The asphaltene structure is generated stochastically using Monte Carlo techniques. Individual asphaltene molecules are approximated as spheroids for the purpose of calculating their hindered diffusivities. Continuum hydrodynamic theories and boundary element methods are used to calculate the diffusion coefficients. A number of analytical expressions, scaling relationships, and approximations utilized in the literature are evaluated. 1. Introduction It is well-known that the diffusivity of large macromolecules in narrow pores can deviate significantly from bulk solution diffusivity. Diffusion of such molecules in restricted porous systems is described as hindered and is of significance during the transport of liquid macromolecules into catalyst particles and membranes and in various chromatographic separations. Key to the estimation of hindered diffusivities is the calculation of the Stokes resistance reflecting an increased drag in the pore either from “exact” (i.e., numerical) (Haberman and Sayre, 1958; Chen and Skalak, 1970; Paine and Scherr, 1975; Lewellen, 1982; Davidson and Deen, 1988; Kim and Karilla, 1991; Tullock et al., 1992; Phan-Thien and Graham, 1992; Ravi-Kumar et al., 1994; Yang and Kim, 1995) or approximate asymptotic (Wakiya, 1957; Anderson and Quinn, 1974; Bungay and Brenner, 1973; Brenner and Gaydos, 1977; Mavrovouniotis and Brenner, 1988; Nitsche and Balgi, 1994) solutions of the Stokes equations. Of the above, only Wakiya (1957), Chen and Skalak (1970), Davidson and Deen (1988), and Ravi-Kumar et al. (1994) address nonspherical particles. Chen and Skalak (1970), for example, used singularity techniques to solve the Stokes flow problem in a cylindrical tube for spheroidal particles moving along the centerline. Equation 1 below is the expression they
D ) 6πµb0[KU(Λ, 0, Θ h )U - KV(Λ, 0, Θ h )V]
(1)
used for calculating the hydrodynamic viscous drag D of spheroidal particles of spheroidal radius to pore radius ratio b0 and spheroidal dimension Λ translating axisymmetrically with velocity U in a cylindrical pore filled with a fluid of viscosity µ and moving with velocity V. KU and KV are the normalized drag coefficients for diffusion and convection, respectively, and Θ h is the axisymmetrical orientation. Due to the complexity of the resulting expressions, Chen and Skalak tabulated the values of KU(Λ, 0, Θ h ) and KV(Λ, 0, Θ h ) for only a few pairs of (Λ, b0) for b0 values in the range 0 < b0 < 0.8. The present group (Ravi-Kumar et al., 1994) has extended the calculations of Chen and Skalak over a much broader range of b0 values. Averaging over the pore cross section is needed to approximate the observed “average” hydrodynamic hindrance factor (or osmotic reflection coefficient) (Anderson, 1981). In the past, the centerline value has been used as an approximation in lieu of the averaging over the cross section. The validity of this approximation for S0888-5885(96)00598-2 CCC: $14.00
spherical particles has been discussed by Deen (1987). The averaging can be carried out either numerically (Anderson and Quinn, 1974; Lewellen, 1982) or by using asymptotic techniques (Brenner and Gaydos, 1977; Mavrovouniotis and Brenner, 1988; Nitsche and Balgi, 1994; Nitsche, 1995). Of the above studies, only Anderson (1981) and Nitsce (1991, 1995) address nonspherical particles. The overall effective diffusivity includes a partition coefficient since the effective diffusivity of a given molecule through the pore is reduced not only as a consequence of the steric hindrance and the hydrodynamic viscous drag resistance due to the pore walls but also as a result of the partitioning of the particles between the bulk phase and the porous medium. The coefficient describing the partitioning is called the partition coefficient Φ, and the coefficient reflecting the steric hindrance and hydrodynamic resistance is called the hydrodynamic viscous drag coefficient Kd (Deen, 1987). The diffusivity through the pore Dh is given as
Dh ) D∞ΦKd
(2)
where D∞ is the corresponding bulk diffusivity value. Analytical and numerical (Monte Carlo) calculations of partition coefficients of nonspherical molecules have already appeared in the literature (Giddings et al., 1968; Limbach et al., 1989). Hindered diffusion is of particular significance in the hydrocracking and upgrading of crude oil and coalderived liquids. These materials contain macromolecular species whose dimensions often approach the size of typical catalyst pores; as a result, their transport is severely hindered. Understanding the phenomena involved during the hindered transport of such macromolecules is the motivation behind the experimental and modeling transport investigations of this group (Sane et al., 1992; Ravi-Kumar et al., 1994, 1995) and a number of prior investigators (Baltus and Anderson, 1983, 1984). A model system often utilized in studies of transport in lieu of the crude oil itself is asphaltenes. Asphaltenes are a solubility class of compounds found in the original crude characterized by a high heteroatom content and a macromolecular colloidal nature. The choice of asphaltenes as a model system has attracted considerable debate (Bunger and Cogswell, 1981; McKay et al., 1978; Koots and Speight, 1975; Wiehe, 1992; Sane et al., 1994), related to the issue of whether asphaltenes bear simi© 1997 American Chemical Society
Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997 3155
larity to compounds found in the original crude or are simply artifacts of the isolation technique. The view commonly held today is that asphaltenes do indeed exist in the original crude, but to uniquely differentiate them from other crude fractions, one needs more than one property. At least two properties are needed to characterize crude oil fractions (molecular weight and polarity are common choices) and three for coal liquids. Interest in asphaltenes dates back to the 1940s. The literature on their chemistry and structure is voluminous, with a wealth of information on the subject (Strausz et al., 1994; Yen, 1994). The technical details are many and often obscure, but overall, a consistent (but not necessarily universally accepted) picture is starting to emerge. It is apparent that asphaltenes are not simple molecules but rather complex mixtures of components (micelles) of various sizes and shapes consisting of assemblies of smaller particles, which in turn result from the clustering of lower molecular weight components, all in a state of dynamic exchange strongly affected by the presence of solvents, temperature, concentration, and fluid mechanical conditions. The type and exact nature of the forces resulting in the clustering of the lower molecular weight (MW) components are still being studied; the detailed chemical composition of such compounds is still not completely delineated. Sahimi and co-workers (Rassamdana et al., 1996; Rassamdana and Sahimi, 1996), for example, have recently used concepts from nonequilibrium growth phenomena to propose a simple aggregation model of asphaltene particle formation. An important contribution in this area is by Klein and co-workers (Savage and Klein, 1989; Trauth et al., 1994). They have taken the vast knowledge of the chemical structure of asphaltenes and have succeeded in instilling in it a degree of mathematical formalism. 2. Model Development Statistical Model of Asphaltene Transport. A model of asphaltene transport, which is consistent with the literature on the detailed chemistry and structure of asphaltenes, has been previously developed by this group (Ravi-Kumar et al., 1994). This model couples the mathematical formalism of Klein and co-workers, which is capable of encompassing the vast knowledge on the chemical nature and structure of asphaltenes, with a detailed description of the transport of asphaltene clusters through cylindrical pores. The model in its present form does not explicitly account for the cluster agglomeration/delamination phenomena, which characterize the dynamic nature of asphaltene molecules and result in the formation of large MW and diameter (100-300-Å) species defined as micelles in the asphaltene literature, and the omission is deliberate. Though a number of groups (Rassamdana et al., 1996; Rassamdana and Sahimi, 1996; Dabir et al., 96; Sheu et al., 1990, 1992) have recently begun to look at agglomeration/delamination behavior, still little concrete knowledge exists in this area. The mathematical framework developed by Klein and co-workers (Savage and Klein, 1989; Trauth et al., 1994) can, in principle, utilize available chemical and structural data (Speight, 1972; Speight and Moschopedis, 1981; Yen et al., 1984a,b) to describe the structure of asphaltenes and other resid fractions. The technique requires, however, the knowledge of the probability distribution functions for the number of sheets in the asphaltene particles, the number of aromatic and satu-
rated rings, and the number/length of paraffinic chains. Though direct analytical techniques for identifying all these structural elements are available, the approach requires considerable experimental resources and is time consuming; data of this nature are not currently available in the open literature. The potential of the technique of Klein and co-workers lies instead in its proposed use as a tool in an iterative stochastic modeling approach, whereby one tries to determine the probability distribution functions of the structural elements by matching measured macroscopic properties like average molecular weight, density, viscosity, diffusivity, etc. For the calculations presented in this paper, the probability distributions of the various structural elements like the number of unit sheets, the aromatic and saturated rings, and the length of alkyl chains were assumed to be two-parameter log-normal distribution functions, widely used to describe polymer molecular weight distributions:
p)
exp[-(ln x - ln xj)2/2σ2] xx2πσ
(3)
where p is the probability, x is the number/size of the structural element, xj is the mean, and σ is the standard deviation. The choice of the log-normal distribution is one of convenience. Other choices are also available. For example, Trauth et al. (1994) utilized γ distribution functions to develop a representative structure for petroleum resids, while Rassamdana and Sahimi (1996) and Dabir et al. (1996) proposed more complex distributions based on the similarity between asphaltene formation and aggregation phenomena. By adjusting the parameters of the probability distribution functions and using an iterative stochastic structure determination procedure, one can generate asphaltene molecules whose MW distribution and structure match those experimentally measured and reported in the literature. Further details are given elsewhere (Ravi-Kumar et al., 1994, 1995). Calculation of Transport Coefficients of Asphaltenes. There are three approaches available, when one wishes to model the hindered resistance of asphaltene particles through membrane pores and the partitioning between the bulk and the membrane structure, each with its intrinsic advantages and disadvantages. Molecular dynamic calculations of the motion of the particle would, of course, be the preferred technique, but its implementation for asphaltene particles is presently numerically intractable. Monte Carlo simulations of transport and partitioning are another technique, and some preliminary calculations have been carried out. They are the preferred option for calculating partition coefficients but are not without shortcomings in the calculation of transport coefficients, since one is forced to introduce ambiguous parameters for the probability of motion of the molecules and hindered drag resistance. Moreover, the relation between the Monte Carlo time and the real time is not clear. In this paper, we use the continuum hydrodynamic theory to simulate the transport coefficients of asphaltene clusters, which to a first approximation are assumed to be rigid, neutrally buoyant, and uncharged particles moving through a solvent continuum in a porous membrane assumed to consist of nonintersecting, straight, chemically inert, cylindrical pores. The effective diffusivity of solutes in such a system as described
3156 Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997
constructing an “equivalent” spheroidal asphaltene molecule are given by Ravi-Kumar et al. (1994). The partition coefficient between the bulk phase and the pores, defined as the equilibrium ratio of the solute concentration in the pores to that in the bulk solution, is calculated in the model using the corrected (Limbach et al., 1989) statistical mechanics formulae of Giddings et al. (1968) discussed in the Introduction. Only steric interactions are considered for the calculation of the partition coefficient. For an axisymmetrically shaped particle, Φ can be calculated as the fraction of the orientationally averaged pore volume accessible to the center of the molecule. The partition coefficient is given by (Limbach et al., 1989) Figure 1. Schematic diagram of the ellipsoidal particle. ζ is the ratio of a radial coordinate to the ellipsoidal radius; θ is the angle between the pore and the molecule axis; Ra and Rb are the ratios of the ellipsoidal radius and half-axial length to pore radius.
4 π
Φ)
(
1-
in the Introduction is reduced as a result of the partitioning of the solute molecules between the bulk phase and the porous medium and as a consequence of the steric hindrance and the hydrodynamic viscous drag resistance due to the pore walls. The diffusivity Dh is described by eq 2. For the calculation of Φ and Kd, the asphaltene particles have been approximated as spheroids (RaviKumar et al., 1994), following the original suggestion of Baltus and Anderson (1983). The molar volume of the entire asphaltene particle was not used in the calculation of the hindered diffusivities and partition coefficients, since it was found (Tsai et al., 1991) that the hindered diffusivity values so calculated are usually overestimated. For hindered diffusivity and partition coefficient calculations, the critical dimensions of the condensed asphaltene portion have been found to be more suitable. As previously discussed (Ravi-Kumar et al, 1994), the asphaltene particle is thought to consist of a “stack of condensed polyaromatic sheets”, which is referred to as the “condensed portion” and “paraffinic chains” extending away from the condensed portion. The condensed portion is characterized by two dimensions. The first is the “height” of the stack, which can be directly calculated from the average inter-sheet distance, reported in the literature (Yen et al., 1961), while the second dimension is the critical diameter, which is the largest linear dimension of any of the sheets found in the stack. The volume of the condensed portion is taken to be that of the spheroid with these two dimensions. To this volume, the volume of the paraffinic chains is added, which is calculated from the van der Waals radii and the inter-atom distances between carbon-carbon atoms and carbon-hydrogen atoms. Since the paraffinic chains are reported to lie in the plane of the polyaromatic sheet they are attached to (Yen et al., 1961), it is assumed in this model that the “height” of the entire asphaltene particle is the height of the condensed portion. The diameter of the entire asphaltene particle is then calculated based on the total volume, which is the sum of the condensed portion and the molar volume of the paraffinic chains. The diameter of the asphaltene particle is the spheroidal diameter, and it is the major axis for an oblate spheroid and the minor axis for a prolate spheroid. The height of the asphaltene particle is the height of the spheroid, and it is the minor axis for an oblate spheroid and the major axis for a prolate spheroid. Further details about the structure of asphaltene molecules and the technique for
(
∫01dζ∫0π/2dθ sin θ(1 - ζ2)1/2 × (Ra2 - Rb2)cos2 θ + Rb2
{Ra2ζ2 + [(Ra2 - Rb2)cos2 θ + Rb2](1 - ζ2)}1/2
1-
)
×
Ra2[(Ra2 - Rb2)cos2 θ + Rb2] {Ra2ζ2 + [(Ra2 - Rb2)cos2 θ + Rb2](1 - ζ2)}3/2
)
(4) where ζ is the ratio of a radial coordinate to the ellipsoidal radius, θ is the angle between the pore and the molecule axis, and Ra and Rb are the ratios of the ellipsoidal radius and half-axial length to pore radius, respectively (see Figure 1). For the partition coefficient of a sphere in a cylindrical pore, the above equation gives the following well-known expression:
Φ ) (1 - λ)2
(5)
where λ is the ratio of the sphere diameter to pore diameter. The average hydrodynamic viscous drag coefficient Kd is dependent upon the steric and hydrodynamic viscous drag effects due to the presence of the pore wall. As discussed in the Introduction, expressions of Kd are readily available for spherical particles. The most comprehensive expressions covering the entire range of particle-to-pore size ratios have been provided by Bungay and Brenner (1973). They are ad hoc expressions obtained by combining asymptotic centerline and perturbation expansions. According to Bungay and Brenner (1973), the average hydrodynamic viscous drag coefficient of spherical solutes in cylindrical conduits is given by the following equation, commonly used in the literature:
Kd )
6π Kt
(6)
where Kt is given by the equation 2 4 9 an(1 - λ)n + an+3λn Kt ) π2x2(1 - λ)-5/2 1 + 4 n)1 n)0 (7)
[
∑
]
∑
where the coefficients for an are given in Table 1. As previously discussed, other authors have also derived expressions for Kd by radially averaging over the entire pore (Brenner and Gaydos, 1977; Mavro-
Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997 3157 Table 1. Expansion Coefficients for the Hydrodynamic Functions Kt n
an
1 2 3 4
-73/60 77 293/50 400 -22.5083 -5.6117
n
an
5 6 7
-0.3363 -1.216 1.647
vouniotis and Brenner, 1988; Nitsche and Balgi, 1994). For example, Nitsche and Balgi (1994) derived the following analytical expression for the hindered diffusivity, which is accurate to order λ2:
Kd ) (1 - λ)-2[1 + (9/8)λ ln λ 1.539λ + 1.2λ2 + o(λ2)] (8) However, the above expression can only be used for λ < 0.1. Mavrovouniotis and Brenner (1988) derived the following expression for Kd for λ g 0.9:
Kd ) 0.984(1 - λ)9/2
(9)
For arbitrarily shaped rigid molecules (with scaled vector dimensions Λ) in cylindrical pores, Anderson and Quinn (1974) developed an expression for Kd which involves an integral of the normalized hydrodynamic viscous drag coefficient KU(Λ, β, Θ) over all radial positions β and all orientations Θ:
Kd(Λ) )
∫Θ∫0βh(Λ,Θ)βKU-1(Λ, β, Θ) dβ dΘ ∫Θ∫0βh(Λ,Θ)β dβ dΘ
(10)
In the above equation, βh (Λ, Θ) is the radial position, where the molecule first touches the pore wall for any given orientation. To calculate Kd from the above equation, the normalized hydrodynamic drag coefficient KU(Λ, β, Θ) for given values of Θ and β must be computed by solving the Navier-Stokes equations. For creeping Newtonian flow of incompressible viscous fluids, the continuity and Stokes equations of motion are
-∇P + µ∇2V ) 0
(11)
∇‚V ) 0
(12)
Once the Stokes equations are solved, the drag of the molecule can be calculated and then KU can be derived. Due to the random nature of the particles’ motion for a given Θ and β, the boundary conditions are not clearly defined. However, assuming that the random fluctuations in the motion of the molecule can be averaged over a long time and noting that the balance between the driving force and the hydrodynamic drag leads to a steady particle motion, the problem can be reformulated by assuming that the particle is moving axially with a steady arbitrary velocity U and then writing the corresponding boundary conditions. Once the Stokes equations have been solved, the drag on the entire molecule is known, and hence, KU can be calculated. For nonspherical particles, the rigorous calculation of KU(Λ, β, Θ) for all β and Θ is quite formidable, however. As discussed in the Introduction, only a handful of authors have attempted such calculations. Chen and Skalak (1970) proposed a singular perturbation technique to calculate KU for spheroidal particles which translate
axisymmetrically along the centerline. Their technique was utilized by this group (Ravi-Kumar et al., 1994) to calculate Kd over a wide region of particle sizes. The validity of the centerline approximation for such molecules has never been tested before. In this paper, KU is calculated by direct numerical integration of the original set of Stokes equations. Boundary element methods (BEM) will be used to calculate the hindered diffusivity of spheres and spheroids in cylindrical pores. The BEM formulations of Tullock et al. (1992), Kim and Karrila (1991), and Yang and Kim (1995) have been extended to study the diffusion of spheroidal particles. Calculations have been performed for various radial positions within the pore but only for aligned orientations, i.e., for the axisymmetric axis of the spheroid being parallel to the pore axis. How close the Kd calculated based on this orientation comes to the true Kd for the particle motion within the pore is still under investigation. One expects the calculations to be more accurate for prolate spheroids in narrow pores. The numerical solution of the Stokes equations for configurations other than axisymmetric are currently not available. In addition, our model assumes dilute solutions in which various asphalteneasphaltene interactions are absent. The purpose of these calculations is 2-fold. The first goal is to extend the region of conditions for which values of Kd have been calculated. This is important if one is to incorporate such values in the statistical model of asphaltene transport, since asphaltenes are polydisperse entities with a broad MW distribution. The second goal is to test the validity of the centerline approximation for spheroidal particles. The latter has implications which go beyond the narrow confines of the model of asphaltene transport. The newly calculated Kd values will be utilized, furthermore, to test an assumption commonly utilized in reactor design involving petroleum and coal macromolecules. The hindered diffusivity of such compounds is typically expressed in terms of their Stokes-Einstein radius, often following a simple exponential dependence. The model presented here does not support the validity of such an assumption over a broad range of conditions. Boundary Element Method. The Stokes equations can be solved, at least in principle, by using either BEM or finite element methods (FEM). The advantages of the BEM over FEM for this class of problems have been discussed in a number of classical papers (Youngren and Acrivos, 1975; Zick and Homsy, 1982). In BEM, the set of Stokes equations are first recast into a boundary integral equation (BIE) for the velocity field (Tullock et al., 1992; Kim and Karrila, 1991; Yang and Kim, 1995). Subsequently, one divides the boundary into a series of M elements over which the geometry, velocity, and traction are approximated by piecewise polynomials. The BIE written in a discretized form is applied to a series of collocation points over the boundary, leading to a set of linear algebraic equations. Given appropriate traction and velocity boundary conditions, these equations are solved (by Gaussian elimination) for the remaining boundary values. In this paper, the completed double-layer boundary integral equation method, presented in Kim and Karrila (1991), is used. In their approach, the integral equations are derived by constructing the velocity field as a superposition of the relevant Green function responses to point forces in the cylinder (Yang and Kim, 1995). Each spheroid was descretized using a 320-element
3158 Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997
Figure 2. 1/KU vs the sphere radial position β for different values of the scaled diameter of the sphere λ.
Figure 3. 1/KU vs the oblate spheroid radial position for different values of the scaled diameter of the oblate spheroid λ.
mesh. The calculations were performed on a parallel computer (IBM SP2) with parallel environment (IBM AIX). The IBM AIX Parallel Environment (PE) program (IBM, 1994), which uses a distributed memorypassing system, is used to schedule the tasks. A singleprogram multiple-data (SPMD) model was used to assign the calculation of each spheroid to a different processor. The SPMD program executes a number of individual, but related, parallel tasks on a number of system processor nodes. 3. Results and Discussion Calculations were peformed for both oblate and prolate spheroids for various radial positions in a cylinder. To verify the accuracy of the numerical technique, similar calculations were also performed for spheres. The spheroidal particles are characterized by two dimensionssthe spheroidal thickness (a) and the spheroidal diameter (b). For such particles, the scaled diameter λ is defined as the ratio of the spheroidal diameter to the pore diameter. For spheroidal particles, the calculations have been performed for the case in which the axisymmetric axis of the spheroid stays parallel to the pore axis. To calculate the average hydrodynamic viscous drag coefficient Kd, one must calculate the inverse of the normalized hydrodynamic drag coefficient 1/KU for different radial positions β and radially average over the entire pore (see eq 10). How different the value of Kd so calculated is from the value based on the centerline approximation (where one assumes that the particle only translates along the centerline) will, of course, depend on how strongly 1/KU depends on β. Such data are shown in Figures 2 and 3. In Figure 2, the inverse normalized hydrodynamic drag coefficient (1/KU) for spheres is plotted vs the radial position (β) for different values of the scaled diameter of the sphere. The value of 1/KU deviates significantly from the centerline value only for values of β close to the wall. In Figure 3, 1/KU for an oblate spheroid is plotted vs its radial position in the tube for different values of its scaled diameter (λ). As in the case of the sphere (Figure 2), 1/KU deviates from its centerline value significantly only for positions close to the tube wall. Similar behavior is also observed with prolate spheroids. The value of Kd was found by integrating eq 10 using a closed-type Newton-Cotes 11-point formula (Abramowitz and Stegun, 1964). Since data for spheres are readily
Figure 4. Plot of average hydrodynamic drag coefficient Kd vs the ratio of particle diameter to pore diameter λ for spherical particles. Comparison of BEM results with literature data: (O) BEM, centerline; (b) BEM; (+) centerline, Paine and Scherr (1975); (-‚-) Nitsche and Balgi (1994); (s) centerline, Bungay and Brenner (1973).
available, the values of the average hydrodynamic viscous drag coefficients of spheres in cylindrical pores are calculated to test the accuracy of the BEM technique and are shown in Figure 4 as a function of the ratio of sphere diameter to pore diameter. The solid circles are calculated using BEM and averaging over the pore radius. The BEM data for spheres moving along the centerline are shown as the open circles in Figure 4. For comparison, the literature data are also plotted in the figure. The solid line is the analytical expression derived by Bungay and Brenner (1973) using asymptotic analysis (see eqs 6 and 7). The + symbols are the numerical data of Paine and Scherr (1975) for the centerline approximation. The agreement between the BEM calculations and the literature data for the centerline approximation is excellent. The dotted line is the analytical expression of Nitsche and Balgi (1994) for small values of λ. Again, it very closely follows the BEM calculations averaging over the pore radius. Figure 5 shows the calculated Kd for both the centerline case and for the radially averaged case for spheroidal particles. In this figure, the results are shown for two values of the aspect ratio a/b of 0.5 and 2. The centerline case overestimates the value of Kd. The deviation is not significant, however, over the range of λ values of relevance to the asphaltene-transport cal-
Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997 3159
Figure 5. Plot of average hydrodynamic drag coefficient Kd vs λ for spheroids with different aspect ratios.
culations. The degree of deviation, furthermore, between the centerline and radially averaged values does not appear to be affected by changes in the spheroidal aspect ratio a/b. Figure 6 compares Kd based on the centerline approximation using BEM with the results of Chen and Skalak (1970) for spheroids. In Figure 6, Kd is plotted vs a/b for different values of λ. For comparison, the published data of Paine and Scherr (1975) for spheres (aspect ratio of 1) are also included in the figure. The BEM data are in good agreement with literature data, indicating the accuracy of the present BEM calculations and, of course, of the prior literature asymptotic calculations. The transport behavior of spheroids and other particles of arbitrary shape has often been approximated by the behavior of spheres with an equivalent StokesEinstein radius. Baltus and Anderson (1984) used the mean projected radius to relate the Stokes-Einstein radius to the partitioning behavior of ellipsoidal molecules. In the literature, the Stokes-Einstein radius as of an equivalent sphere to a given spheroid is calculated using the following expressions (Happel and Brenner, 1991):
as )
kT 6πηD∞
Figure 6. Plot of average hydrodynamic drag coefficient Kd vs aspect ratio; a, spheroid axial length; b, axisymmetric diameter; λ, scaled axisymmetric diameter; comparison of BEM results with literature data.
Figure 7. Partition coefficient Φ vs λs, the ratio of the equivalent Stokes-Einstein radius to the pore radius, for various spheroid sizes.
(13)
1 2 D∞ ) D∞| + D∞⊥ 3 3
(14)
where D∞| and D∞⊥ are the diffusivities of the molecule in the directions parallel and perpendicular to the orientation of the axis of rotation, respectively. Equations 13 and 14 are used to calculate as for spheroids. For a prolate spheroid, the equation is
{
[
]
as 8 22 - 1 + (2 - 1)1/2 ) ln + a (2 - 1)3/2 - (2 - 1)1/2 2(22 - 3) (2 - 1)3/2
}
-1
ln[ + (2 - 1)1/2]
(15)
where ) b/a, and for an oblate spheroid,
as ) (1 - 2)1/2[sin-1(1 - )1/2]-1 b
(16)
where b is the major axis of the spheroid and, in this case, ) a/b.
Figure 8. Average hydrodynamic drag coefficient Kd vs λs, the ratio of the Stokes-Einstein radius to the pore radius, for various spheroid sizes.
The calculated values of Φ are shown in Figure 7 vs the ratio λs of the Stokes-Einstein radius to the pore radius, for different values of the spheroidal aspect ratio a/b. Figure 8 shows the calculated values of Kd vs λs over the spheroidal size ranges of interest in this paper.
3160 Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997 Table 2. Parameters of the log-Normal Distribution Functions for the Various Structural Elements of an Asphaltene structural element
mean
SD
no. of sheets length of chains no. of aromatic rings
3.58 6.77 10.50
0.26 0.23 0.30
The lines corresponding to a/b equal to 1 are the result calculated for a sphere. As expected, all the curves converge to unity for small spheroidal sizes. Note that in addition to its dependence on λ, the drag coefficient is also a function of the spheroidal aspect ratio. This points out the difficulty in trying to approximate the transport behavior of spheroids with that of spheres of the same equivalent Stokes-Einstein radius. The same is true for the partition coefficients, as is obvious from Figure 7. Asphaltene Diffusivity. The hindered diffusivity data generated for spheroidal particles are now utilized to calculate the diffusivity of asphaltene molecules. To generate the structure of the asphaltenes, model results were averaged over five simulation runs, each involving 200 000 particles. The MW distribution of the asphaltene molecules is calculated by counting the number of particles that are within a particular range of molecular weights. The structural element probability distributions of the number of unit sheets of aromatic and saturated rings and of the length of alkyl chains were all assumed to be two-parameter log-normal distribution functions (see eq 3). The parameters of the distribution functions are shown in Table 2 and correspond to the asphaltenes studied by Savage and Klein (1989). The generation of the asphaltene molecules has been described in detail elsewhere (Ravi-Kumar et al., 1994). After one defines the various probability distribution functions, the number of unit sheets in every asphaltene particle and the number of aromatic rings and saturated rings in every sheet are determined. The aromatic and saturated rings are stochastically arranged with a condensation index of 0.8. The condensation index, as defined by Savage and Klein (1989), depends on the number of peripheral carbons. The aromatic and saturated rings in the sheet are packed as closely as possible. The packing is done by starting with one ring at the center and then adding rings adjacent to it in a clockwise manner. The molecule is now peri-condensed. Then the peripheral rings are arranged randomly so that the condensation index, calculated from the number of peripheral rings, is close to 0.8. The number of peripheral carbons and the number of peripheral carbons that can be substituted with aliphatic chains are then calculated. The peripheral carbons are randomly substituted with alkyl chains, whose lengths are determined from the cumulative distributions. The resulting diffusion coefficients are plotted vs the MW by averaging over all the molecules within a range of MWs. Figure 9 shows how the diffusion coefficient changes with the MW and the pore size. The diffusion coefficients decrease with increasing MW and decrease for more constricted pores. The calculated asphaltene diffusivities are of the same order of magnitude as the experimental results in the literature (Baltus and Anderson, 1983; Sane et al., 1992). Note that ln(Dh/ D∞) scales linearly with MW1/3; i.e.,
Dh ∼ exp(-KMW1/3) D∞
(17)
Figure 9. ln(Dh/D∞) vs the molecular weight MW of a simulated asphaltene through membranes of varying pore diameter. The straight lines are best fit lines in the figure.
Figure 10. Plot of ln(Dh/D∞) as a function of λs, the ratio of Stokes-Einstein radius to the pore radius, for varying pore sizes. Plotted on the figure are also the best-fit exponential line (-‚-) and the Renkin (s, eq 18) and power-law (- - -, eq 19) expressions.
The scaling coefficient K, however, is a function of the cylindrical pore size, as can be seen in Figure 9. In Figure 10, ln(Dh/D∞) is plotted vs the ratio of the average Stokes-Einstein radius to the pore radius. The calculated data for a broad range of pore sizes and λs all follow a smooth curve. However, this is not a simple exponential function as one often encounters in the engineering literature. The best exponential fit, exp(-5.13λs), is also shown in Figure 10. It deviates significantly from the calculated data for λs > 0.3. As can be seen in Figure 10, the other two relationships commonly utilized in the engineering literature, namely, the Renkin relationship (Pappenheimer et al., 1951) given by
Dh ) (1 - λs)2(1 - 2.1044λs + 2.089λs3 - 0.948λs5) D∞ (18) and the power-law relationship (Deen, 1987)
Dh ) (1 - λs)4 D∞
(19)
are much more accurate approximations of the numerical data. The Renkin relationship, in particular, provides a very satisfactory fit of the results.
Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997 3161
4. Conclusions A statistical model of asphaltene transport, previously developed by our group, has been utilized to study unique features of the asphaltene transport through cylindrical pores. The model utilizes the statistical molecular structure models developed by Klein and coworkers to generate the asphaltene molecules. Hydrodynamic continuum theories have been used for calculating the partition coefficients and hindered resistance for asphaltene molecules by approximating them as spheroids. Boundary element methods have been used to calculate the hindered diffusivity of spheres and spheroids in cylindrical pores for particles translating along the centerline and away from the centerline. The diffusivity of a generic asphaltene molecule has also been calculated. The diffusivities of asphaltene components are plotted vs the average MW and the ratio of the Stokes-Einstein radius to the pore radius, calculated by averaging over narrow ranges of MWs. Though an exponential scaling relationship exists between Dh/D∞ and MW1/3, i.e., Dh/D∞ ∼ exp(-KMW1/3), the exponent K is not universal but shows a strong dependence on the pore size. When hindered diffusivities are plotted vs λs, the average ratio of the StokesEinstein radius to the pore radius, the data for a broad range of pore radii and λs seem to form a smooth curve. However, the relationship is not a simple exponential one as it is commonly utilized in the engineering literature. The model described here is at an early stage of development. For a more realistic description of asphaltene transport, one must incorporate into the model the physics of the agglomeration/delamination process (Rassamdana et al., 1996; Rassamdana and Sahimi, 1996; Dabir et al., 1996). The idea of approximating an asphaltene particle as a rigid spheroid, though useful for making the calculations numerically tractable, is, of course, still an oversimplification. The validity of any model (and its further development) must eventually be judged by agreement of its predictions with experimental data. To critically evaluate the asphaltene-transport model presented here, one must have structural information and diffusion data over a wide range of conditions. Though asphaltene diffusivity data through membranes are available (Baltus and Anderson, 1983; Sane et al., 1992), they are not comprehensive enough at this time to allow meaningful comparisons. Work in this area is continuing in our group. Acknowledgment This paper is written to commemorate the 65th birthday of Professor Gilbert Froment, a valued friend of many years. We acknowledge the support of the National Science Foundation and the U.S. Department of Energy. We also thank Professor Sangtae Kim of the University of Wisconsin-Madison for making available to us the BEM code and for helpful discussions. Nomenclature as ) Stokes-Einstein radius b0 ) ratio of spheroidal radius to pore radius D ) hydrodynamic viscous drag Dh ) hindered diffusivity D∞ ) bulk diffusivity D∞| ) bulk diffusivity of molecules in the direction parallel to the axis of rotation
D∞⊥ ) bulk diffusivity of molecules in the direction perpendicular to the axis of rotation Kd ) hydrodynamic viscous drag coefficient KU ) normalized drag coefficient for diffusion KV ) normalized drag coefficient for convection MW ) molecular weight p ) probability P ) pressure T ) temperature V ) velocity Greek Symbols Ra ) ratio of ellipsoidal radius to pore radius Rb ) ratio of half axial length to pore radius β ) radial position ζ ) ratio of radial coordinate to the ellipsoidal radius η ) kinematic viscosity Θ ) angle between pore and molecule axes; see Figure 1 λ ) ratio of molecular diameter to pore diameter λs ) ratio of Stokes-Einstein radius of solute to pore radius Λ ) spheroidal dimension µ ) viscosity σ ) standard deviation Φ ) partition coefficient
Literature Cited Abramowitz, M.; Stegun, I. A. Handbook of Mathematical Functions; Dover: New York, 1964. Anderson, J. L. Configurational Effect on the Reflection Coefficient for Rigid Solutes in Capillary Pores. J. Theor. Biol. 1981, 90, 405. Anderson, J. L.; Quinn, J. A. Restricted Transport in Small Pores: A Model for Steric Exclusion and Hindered Particle Motion. Biophys. J. 1974, 14, 130. Baltus, R. E.; Anderson, J. L. Hindered Diffusion of Asphaltenes through Microporous Membranes. Chem. Eng. Sci. 1983, 38, 1959. Baltus, R. E.; Anderson, J. L. Comparison of GPC Elution Characteristics and Diffusion Coefficients of Asphaltenes. Fuel 1984, 63, 530. Brenner, H.; Gaydos, L. J. The Constrained Brownian Movement of Spherical Particles in Cylindrical Pores of Comparable Radius. J. Colloid Interface Sci. 1977, 58, 312. Bungay, P. M.; Brenner, H. The Motion of a Closely-Fitting Sphere in a Fluid-Filled Tube. Int. J. Multiphase Flow 1973, 1, 25. Bunger, J. W.; Cogswell, D. E. Characteristics of Tar Sand Bitumen Asphaltenes as Studied by Conversion of Bitumen by Hydropyrolysis. In Chemistry of Asphaltenes; Bunger, J. W., Li, C. N., Eds.; Advances in Chemistry Series; American Chemical Society: Washington, DC, 1981. Chen, T. C.; Skalak, R. Stokes Flow in a Cylindrical Tube Containing a Line of Spheroidal Particles. Appl. Sci. Res. 1970, 22, 403. Dabir, B.; Nematy, M.; Mehrabi, A. R.; Rassamdana, H.; Sahimi, M. Asphalt Flocculation and Deposition. III. The Molecular Weight Distribution. Fuel 1996, 75, 1633. Davidson, M. G.; Deen, W. M. Hydrodynamic Theory for the Hindered Transport of Flexible Macromolecules in Porous Membranes. J. Membrane Sci. 1988, 35, 167. Deen, W. M. Hindered Transport of Large Molecules in Liquidfilled Pores. AIChE 1987, 33, 1409. Giddings, J. C.; Kucera, E.; Russell, C. P.; Myers, M. N. Statistical Theory for the Equilibrium Distribution of Rigid Molecules. Exclusion Chromatography. J. Phys. Chem. 1968, 72, 4397. Haberman, W. L.; Sayre, R. M. Motion of Rigid and Fluid Spheres in Stationary and Moving Liquids inside Cylindrical Tubes. Technical Report, David Taylor Model Basin. Report No. 1143. U.S. Navy, Washington, DC, 1958. Happel, J.; Brenner, H. Low Reynolds Number Hydrodynamics: with Special Applications to Particulate Media; Prentice-Hall: Englewood Cliffs: NJ, 1991. IBM. IBM AIX Parallel EivironmentsOperation and Use. Technical Report, Release 2.1, 2nd ed., 1994. Kim, S.; Karrila, S. Microhydrodynamics. Principles and Selected Applications; Butterworth-Heinemann: Boston, 1991.
3162 Ind. Eng. Chem. Res., Vol. 36, No. 8, 1997 Koots, J. A.; Speight, J. G. Relation of Petroleum Resins to Asphaltenes. Fuel 1975, 54, 179. Lewellen, P. C. Hydrodynamic Analysis of Microporous Mass Transport. Ph.D. Thesis, The University of WisconsinMadison, 1982. Limbach, K. W.; Nitsche, J. M.; Wei, J. Partitioning of Nonspherical Molecules Between Bulk Solution and Porous Solids. AIChE J. 1989, 35, 42. Mavrovouniotis, G. M.; Brenner, H. Hindered Sedimentation and Dispersion Coefficients for Rigid, Closely Fitting Brownian Spheres in Circular Cylindrical Pores Containing Quiescent Fluids. J. Colloid Interface Sci. 1988, 124, 269. McKay, J. F.; Amend, P. J.; Cogswell, T. E.; Harnsberger, P. M.; Erikson, R. B.; Latham, D. R. Petroleum Asphaltenes. Chemistry and Composition. Analytical Chemistry of Liquid Fuel Sources: Tar Sands, Oil Shale, Coal, and Petroleum; Uden, Peter C., Jensen, Howard B., Eds.; Advances in Chemistry Series; American Chemical Society: Washington, DC, 1978; Vol. 170, p 128. Nitsche, J. M. Hydrodynamic Coupling and Non-equilibrium Distribution in Pore Diffusion of Non-spherical Particles. Part. Sci. Technol. 1991, 9, 135. Nitsche, J. M. Pore Diffusion of Nonspherical Brownian Particles. Ind. Eng. Chem. Res. 1995, 34, 3606. Nitsche, J. M.; Balgi, G. Hindered Brownian Diffusion of Spherical Solutes within Circular Pores. Ind. Eng. Chem. Res. 1994, 33, 2242. Paine, P. L.; Scherr, P. Drag Coefficients for the Movement of Rigid Spheres Through Liquid-filled Cylindrical Pores. Biophys. J. 1975, 15, 1087. Pappenheimer, J. R.; Renkin, E. M.; Borrero, L. M. Filtration, Diffusion and Molecular Sieving through Peripheral Capillary Membranes. Am. J. Physiol. 1951, 167, 13. Phan-Thien, N.; Graham, A. Translation and Rotation of Spheres Settling in Square and Circular Conduits: Experiments and Numerical Predictions. Int. J. Multiphase Flow 1992, 18, 1061. Rassamdana, H.; Sahimi, M. Asphalt Flocculation and Deposition. II. Formation and Growth of Fractal Aggregates, AIChE J. 1996, 42, 3318. Rassamdana, H.; Dabir, B.; Nematy, M.; Farhani, M.; Sahimi, M. Asphalt Flocculation and Deposition I. The Onset of Precipitation. AIChE J. 1996, 42, 10. Ravi-Kumar, V. S.; Tsotsis, T. T.; Sahimi, M.; Webster, I. A. Studies of Transport of Asphaltenes Through Porous Membranes: Statistical Structural Models and Continuum Hydrodynamic Theories. Chem. Eng. Sci. 1994, 49, 5789. Ravi-Kumar, V. S.; Tsotsis, T. T.; Sahimi, M.; Webster, I. A. Inorganic Membranes for the Study of Transport and Reaction of Petroleum and Coal Liquids. Proceedings of the 3rd International Conference on Inorganic Membranes, Worcestershire, MA, 1995. Sane, R. C.; Tsotsis, T. T.; Webster, I. A.; Ravi-Kumar, V. S. Studies of Asphaltene Diffusion and Structure and Their Implications for Resid Upgrading. Chem. Eng. Sci. 1992, 47, 2683. Sane, R. C.; Tsotsis, T. T.; Webster, I. A.; Ravi-Kumar, V. S. Studies of Asphaltene Diffusion. Implications For Asphaltene Structure And Optimal Upgrading Reactor Design. In Asphaltenes and Asphalts; Developments in Petroleum Science; Yen, T. F., Chilingar, G. V., Eds.; Elsevier-Science: Amsterdam, 1994; Vol. 1.
Savage, P. E.; Klein, M. T. Asphaltene Reaction Pathways. V. Chemical and Mathematical Modeling. Chem. Eng. Sci. 1989, 44, 393. Sheu, E. Y.; Liang, K. S.; Sinha, S. K.; Overfield, R. E. Polydispersity of Asphaltenes in Toluene. Prepr. Pap.sAm. Chem. Soc., Div. Pet. Chem. 1990, 35, 813. Sheu, E. Y.; De Tar, M. M.; Storm, D. A.; DeCanio, S. J. Aggregation and Kinetics of Asphaltenes in Organic-solvents. Fuel 1992, 71, 299. Speight, J. G. The Application of Spectroscopic Techniques to the Structural Analysis of Coal and Petroleum. Appl. Spectrosc. Rev. 1972, 5, 211. Speight, J. G.; Moschopedis, S. E. On the Molecular Nature of Petroleum Asphaltenes. Chemistry of Asphaltenes; Bunger, J. W., Li, C. N., Eds.; Advances in Chemistry; American Chemical Society: Washington, DC, 1981; Vol. 195. Strausz, O. P.; Lown, E. M.; Mojelsky, T. W. The Molecular Structure of Asphaltene. Prepr. Pap.sAm. Chem. Soc., Div. Fuel. Chem. 1994, 39, 199. Trauth, D. M.; Stark, S. M.; Petti, T. F.; Neurock, M.; Klein, M. T. Representation of the Molecular Structure of Petroleum Resid through Characterization and Monte Carlo Modeling. Energy Fuels 1994, 8, 576. Tsai, C. H.; Massoth, F. E.; Lee, S. Y.; Seader, J. D. Effects of Solvent and Solute Configuration on Restrictive Diffusion in Hydrotreating Catalysts. Ind. Eng. Chem. Res. 1991, 30, 22. Tullock D. L.; Phan-Thien, N.; Graham, A. L. Boundary Element Simulations of Spheres Settling in Circular, Square and Triangular Conduits. Rheol. Acta 1992, 31, 139. Wakiya, S. J. J. Phys. Soc. Jpn. 1957, 12, 1130 (as cited by Happel and Brenner (1991)). Wiehe, I. A. A Solvent-Resid Phase Diagram for Tracking Resid Conversion. Ind. Eng. Chem. Res. 1992, 31, 530. Yang, H.; Kim, S. Boundary Element Analysis of Particle Mobilities in a Cylindrical Channel: Network-based Parallel Computing with Condor. Comput. Chem. Eng. 1995, 19, 683. Yen, T. F. The Overall Properties of Asphaltenes as Deduced from Structure. Prepr. Pap.sAm. Chem. Soc., Div. Pet. Chem. 1994, 39, 198. Yen, T. F.; Erdman, J. L.; Pollack, S. P. Investigation of the Structure of Petroleum Asphaltenes by X-ray Diffraction. Anal. Chem. 1961, 33, 1587. Yen, T. F.; Wu, W. H.; Chilingar, G. V. A Study of the Structure of Petroleum Asphaltenes and Related Substances by Infrared Spectroscopy. Energy Sources 1984a, 7, 203. Yen, T. F.; Wu, W. H.; Chilingar, G. V. A Study of the Structure of Petroleum Asphaltenes and Related Substances by Proton Nuclear Magnetic Resonance. Energy Sources 1984b, 7, 275. Youngren, G. K.; Acrivos, A. Stokes Flow Past a Particle of Arbitrary Shape: A Numerical Method of Solution. J. Fluid Mech. 1975, 69, 377. Zick, A. A.; Homsy, G. M. Stokes Flow Through Periodic Arrays of Spheres. J. Fluid Mech. 1982, 115, 13.
Received for review September 30, 1996 Revised manuscript received March 11, 1997 Accepted March 17, 1997X IE960598C X Abstract published in Advance ACS Abstracts, June 15, 1997.
