Simulation of vapor diffusion in anisotropic particulate deposits

Chemical Ewinm Printed in Great

Science, Britain.

Vol. 47, No. 2. pp. 421-443,

1992.

ow-2So9p2 S5.00 + 0.00 Q 1991 Fkrgamon Press plc

OF VAPOR DIFFUSION IN ANISOTROPIC PARTICULATE DEPOSITS+

SIMULATION

MENELAOS TASSOPOULOS Yale University, Chemical Engineering Department, Laboratory,

(First

received

and DANIEL

E. ROSNER*

High Temperature Chemical New Haven, CT 06520-2159, U.S.A.

2 Jnnunry

1991; accepted in revisedfirm

Reaction

Engineering

18 April 1991)

Abstract-In this study we determine the effective diffusivity tensor of computer-generated, particulate, anisotropic deposits. Deposits of polydispersed, spherical particles were generated by several algorithmic models, which simulate different deposition conditions. Depending on the nature of the arriving particle trajectories (diffusive/deterministic), the angle of incidence and the subsequ&t motion after they establish contact with a.n already deposited particle, both the deposit anisotropy and 6nal solid fraction can be varied. Monte Carlo simulations of vapor molecules difIiing through cubic samples extracted from the generated deposits provided accurate estimates of their effective diffusivity tensor and associated tortuosities. The effective cliffusivities were directly determined over the entire Knudsen number range, assuming self-diffusion in a constant pressure pure-gas that consists of hard-sphere elastic molec& &ith exl& nentially distributed free-paths. In the continuum limit the effective diffusivities of the Mrous solids conside& here could be closely correlated with porosity alone. Archie’s power law [eq. {17)] with an exponent m = 1.45 fitted our computational results accurately. In the Knudsen limit on the other hand the details of the local microstructure became important and so an accurate prediction/correlation of the effective ditiivity required higher order information. An interpolation based on the experimentally measured effective diffusivities in the two Knudsen limits (Scott and Dullien, 1962) was used to correlate the diffusivities in the transition regime. For the deposits analyzed here the diffusivity properties became anisotroGc onlv in the hi& Knudsen reaime. that is. when diffusinn molecule/solid interactions dominate.

In the I&& limit we-have also fo6d s&nificak “rectification;’effects. &ecifically, the transmission probability as determined by the test-molecule method could be ut, to 25% hiher when the molecules diffused away from the tar&t than when they diffused towards the target. Such%ctification” effects are attributed to the tree-like microstructures of the deposits.

1. INTRODUCRON Mass transport in porous media is a well-studied problem since it occurs in many engineering appli-

cations and in materials science. Gas and vapor transport in porous solids can be very important when considering catalytic and noncatalytic heterogeneous reactions (Sat&-field and Sherwood, 1963; Satterfield, 1970), certain composite material fabrication techniques like chemical vapor infiltration (CVI) (Melkote and Jensen, 1989; Sotirchos and Tomadakis, 1989), as well as in particle filtration. Related examples include the fabrication of optical wave guides and various coatings, as well as in combustion turbines, and heat exchangers where particle deposition often occurs in the presence of condensible vapors. In the latter case both the temporal evolution of the deposit, as determined For example by the fraction of impinging particles that stick (Rosner and Nagarajan, 1987), and the deposit final properties (transport or optical properties, chemical activity and uniformity, resistance to wear, etc.) can depend on the rate of vapor uptake by the growing deposit, clearly a function of the local effective diffusivity tensor, &. One approach for predicting the effective transport coefficients through porous media is based on theoretical models. An important example is the dusty-gas

+HTCRE No. 180. *Author to whom correspondence

should be addressed. 421

model (Mason and Malinauskas, 1983) where the porous medium is treated as an extra component of the gas phase, consisting of very heavy, immobile molecules. The effect of the solid on the gas is then determined from the kinetic theory of gases applied to this disparate molecular weight gas mixture. A-nother route is provided by using network models (see, Sahimi et al., 1990 for a recent review; also Jackson, 1977 and Sotirchos, 1989) where a simplified model of the actual microstructure is proposed, the desired quantities are obtained at the microscopic level (say for a single capillary tube) directly from the governing equations, and the macroscopic effective quantities are then expressed as appropriate averages (Whitaker, 1967; Slattery, 1972; Sahimi et al., 1990). Still in the same spirit, the actual porous solid is sometimes approximated by a periodic structure, thus, reducing the problem to the determination of the property of interest for a single unit cell (see Nitsche and Brenner, 1990 for a review). Finally, the effective property of interest can also be estimated from rigorous bounds that incorporate explicit information about the local microstructure (see Torquato, 1987 for a review). A completely different and rather “brute-force” approach involves actual Monte Carlo (MC) simulations of the diffusive motion of the gas molecules through many realizations of the porous solid. Such random walk techniques, although well-known for many years (Hammersley and Handscomb, 1964), have become practical only recently, thanks to the

422

ME-LAOS

TASSOP~IJLOS and DANIEL E. RUNNER

advancement of computers and the development of efficient algorithms. The main advantage of such MC methods when compared to the more traditional approaches outlmed earlier is that the actual rnicrostructure of the porous sample can be represented in the computer to any desired degree of accuracy, thus, enabling an accurate estimation of the effective property of interest, here the diffusivity tensor, provided, of course, a large enough number of molecules is followed for large enough times (“large” will be quantified in Section 3). Evidently Evans and coworkers (Evans et al., 1980, Abbasi et al., 1983) were the first to report Ddf results of model porous solids obtained by MC simulations. Using a method based ou the fraction of molecules transmitted through a porous slab (see Section 3), they determined Den for assemblages of penetrable spheres randomly arranged in space in both limits of Knudsen and ordinary (bulk) diffusion. Later, Smith and his collaborators (Huizenga and Smith, 1986; Smith, 1986; Olague et al., 1988) measured the effective diffusivity in the free molecule limit (Kn = A/d, %- 1 where 2 is the prevailing gas meanfree-path and d, is the mean particle diameter) of various gases diffusing through ordered and disordered packings of small spherical particles and compared their experimentally obtained tortuosity factors (defined below) to predictions from simulations they performed using Evans’ technique. They reported very good agreement between the two. Schwartz and Banavar (1989), and Schwartz et al. (1989) have also used random walk simulations to determine the electrical properties of insulating consolidated porous solids saturated with a conducting pore fluid. Due to the mathematical equivalence of the electrical and the diffusion problem in the continuum limit (Ku -+ 0), their formation factor, F = ko/keff (where ketf is the effective conductivity and k, is the reference conductivity), can be readily related to an appropriately nondimensionalized effective diffusivity. Kim and Torquato (1990, 1991) have recently developed an efficient random walk algorithm based on first-passage-time probability theory that enabled them to accurately determine the effective thermal conductivity of two-phase media for arbitrary ratio of the conductivities of the two phases. Finally, in a series of papers Burganos and Sotirchos (1988,1989a, b) considered the problem of Knudsen diffusion in networks of capillary structures; see, also Melkote and Jensen (1989), who treated the related problem of molecules diffusing through a substrate of random overlapping fibers. In this study we report on the effective diffusivity tensor, O,,, , of anisotropic media formed by particle deposition from the gas phase, and its dependence on the local Knudsen number. Such porous solids, apart from a practical interest due to their relevance in many engineering applications, exhibit quite different microstructures depending on the prevailing deposition conditions (Section 2). Thus, they vary from the very open “tree-like” structures associated with diffusion-limited aggregation (DLA), to the

columnar microstructures obtained under “ballistic” deposition conditions (solid fraction, Q N O-15), to dense configurations approaching the so-called random loose packing of spherical particles (4 CL0.60) if post-collision rolling events occur. D,, was determined, assuming self-diffusion in a constant pressure pure gas that consists of hard-sphere elastic molecules moving at constant speed. The diffusing vapor molecules depart from vapor molecule collisions isotropicdly and are assumed to be reflected from the particle surface mrding to the laws of diffuse reflection. To our knowledge, this is the first time that “exact” results on the combined effects of local microstructure and Knudsen number are reported, apart from the work of Melkote and Jensen (1990) who, as we recently found out, have been independently and concurrently extending their work to treat transition diffusion through fibrous structures. Finally, we should mention the elegant work of Strieder (1971) and Ho and Strieder (1979,1980), who, assuming the same model as here for gaseous diffusion in the transition regime, derived variational upper-bounds on D crf that incorporated limited information on the local microstructure. The paper is organized as follows. In Section 2, we briefIy describe the models we use to generate particulate deposits, and report on some of the relevant microstructural quantities (solid fraction, radial distribution function, distribution of normal contacts, etc.) that characterize them. In Section 3, we present the methods available for obtaining the diffusivity tensor as a function of Knudsen number and deposit anisotropy while in Section 4 we give some details on certain computational aspects and verify our code in the limits of Kn --c 0 and Kn -D a0 . In Section S we report representative results on the dependence of &r and the related formation- and tortuosity-factors on deposit microstructure and Knudsen number, with conclusions and generalizations presented in Section 6.

2 CONSTRUCXiON

OF “MODEL”

DEFQSITS

Predicting the microstructure of a deposit as a function of the relevant deposition parameters can be a difficult proposition, because of the large number of interactions that must be accounted for. In general, we expect that the microstructure of a growing’deposit will depend on: (I) The

particle

size distribution

(PSD)

of the

actually depositing aerosol. Note, that in most cases because of the size-dependent deposition efficiency (Friedlander, 1977) the PSDs of the aerosol in the gas main stream and at the gas/target interface are quite different (Rosner and Tassopoulos, 1989, 1990). Particle shape, i.e. dense primary particles vs aggregates (Vicsek, 1989) or spherical vs elongated particles (Tassopoulos et al., 1989a) can further affect the morphology of the resulting deposit.

Simulationof vapor difrusionin anisotropicparticulate&posits (2) The nature of the incoming particle trajectories, as determined by the interaction of the particles with the various fields (host fluid velocity, and temperature field, as well as radiation and electrical fielda) existing in the fluid boundary layer (3L) adjacent to the target/deposit. As we further discuss below, especially for small enough particles, the final deposit morphology is highly dependent on the relative importance of deterministic (convective/phoretie) and ditrusive motion (Tassopoulos et al., 1989a). (3) The kinetic energy of the particles prior contact with the target or the bed of the deposited particles. It is well-known that larger particles, depositing, say, due to inertia or interception, arrive at the target with a finite velocity, and h&e the final position where they come at rest, (provided, of course, they do eventually stick), is determined by the local impact dynamics as well as the strength of the adhesive and dissipative forces. The former are most often electro-

static in nature, while dissipative

forces can he

423

phoresis. In passing we note that it is possible to simulate deposits grown by particles of arbitrary Peclet number, and thus generate microstructures “between” the DLA and ballistic, by simply letting them perform a biased random walk (Tassopoulos et al., 1989x Houi and Lenormand, 1984) All such models further assume that there is no background fluid affecting the particle motion, the particles arrive at the target one at a time, and upon contact with the target surface *or an already deposited particle, they stick irreversibly. Moreover, the deposit grows in the “frozen” limit, according to which particles already deposited are not allowed to move at subsequent times. We generated three-dimensional off-lattice diffusion-limited and ballistic deposits on targets 45 x 45 and more typically 32 x 32 particle diameters large. Throughout the simulations we employed periodic boundary conditions in the x- and y-directions in order to minimize finite target size effects. Each deposit was grown until it reached an average height of

about 100 particle diameters (all particles assumed spherical), at which point we “extracted” cubic

due to nonelastic deformations or, in the presence of condensible vapors,.a case of interest here, are dominated by viscous dissipation (see, e.g. Barnocky and Davis, 1988). (4) Changes in the microstructure of the existing deposit, because of displacement of the already deposited particles by the incoming particles. (5) Longer term changes due to sintering, condensation of vapors, chemical surface reactions that can lead to pore shrinkage, etc.

samples with an edge equal to the target size to be used in the diffusion simulations. Since it is wellknown that both the part of the deposit close to the target surface and the part close to the evolving top have different structural properties than the “bulk” (see, e.g. Family and Vicsek, 1985). we never sampled from regions within 25 particle diameters from the deposit ‘edges” (in the z-direction). In order to demonstrate ‘more clearly the very different microstruc-

Although a deposition simulation that would account simultaneously for all the above mentioned effbets is still impractical, in recent years many models of varying complexity have been proposed that do capture several features of the microstructures associated with deposits grown in different environments. It is beyond the scope of this communication to review the various approaches, but in a general sense it is possible to distinguish between two classes of deposits, namely deposits formed by essentially Brownian particles “walking” towards the target, and by “ballistic” particles (higher Stokes number) that have appreciable inertia. Deposition of small particles, that arrive at the target with small velocity and so have unity sticking probability, is readily modelled by the &@&on limited aggregation (DLA) (Witten and Sander, 1981) and ballistic aggregation models (Void, 1963; Meakin et al., 1986). The DLA model is valid when the deposit grows under isotropic diffusion-limited conditions, i.e. the depositing particles undergo equilibrium Brownian motion (the dimensionless Peclet number, Pe, (see e.g. Rosner, 1988), that can be viewed as the ratio of a characteristic time for di&sion to a characteristic time for convection, goes to zero), while according to the ballistic model the, particles follow straight line deterministic trajectories (Pe + CO), that locally at least can simulate thmnophoresis or electro-

dimensional lattices. The first one, generated under pure DLA conditions, is characterized by very open tree-like fractal structures; the other two generated under pure ballistic conditions, with deposit lb formed from monodisperse particles arriving normal to the target and deposit lc from particles arriving at a 60” angle as measured from the target normal. For larger particles, that arrive at the target with nonvanishing velocity, the sticking coefficient is almost never unity, but rather the outcome of the impact event as well as the subsequent motion of the impinging and target particles, as determined by the actual collision dynamics. The postcollision particle motion gives rise to denser deposits, and can be still modelled within the frozen limit approximation, provided the impact velocities are not high and the particle interarrival times are longer than the typical time required for the formation of strong “sinter” type bonds (Tassopoulos et al., 1989b; Konstandopoulos and Tassopoulos, 1990; Diemer, 1989). In the opposite extreme where the particle kinetic energy’at impact is very high and/or the adhesive forces between the particles are weak, the deposits approach a possibly anisotropic (Thornton and Barnes, 1986) random loose packing configuration (Berryman, 1983), irrespective of the details of the particle trajectories and particle-particle interactions. Finally, in intermediate cases where the energies associated with the adhesive

CBS

47:2-S

tures associated with DLA and ballistic deposits we show in Fig. 1 three such deposits generated on two-

424

MENELAOS T~ssor=ou~os and DANI@L E. Rosm~

forces acting between the particles and the kinetic energy of the arriving particles at impact are of the same order, it will be necessary to adopt the powerful Distiltct Element method (Cundall and Strack, 1979; Satake and Jenkins, 1988);originally used in problems related to aeolian saltation and granular flows (for applications in the present context see, Rosner et al., 1991). This method can readily handle the simultaneous motion of many particles, but because it is computationally intensive we are not yet able to grow the large deposits required for determining structural and transport properties.

In this work, we generate denser deposits using a simple algorithmic model based on some recent impaction simulations (Tassopoulos et al., 1989b). The aforementioned three dimensional (3D), ballistic type, off-lattice impaction simulations were performed in the frozen deposit limit. Incident particles were modelled as hard, nearly elastic, frictional spheres that arrived at the target with prespecified linear and angular velocities. The outcome of the binary collisions between the incident and target particles were based on the stereomechanical theory of impact (Goldsmith, 1960), while sticking was modelled with

Fig. 1. (a)-(b).

Simulation of vapor diffusion in anisotropk particulate deposits

425

(4 Fig. 1. Typical 2D deposits without particle restructuring. Deposit (a) was generated under pure diffusion limited conditions (DLA), deposit (b) under ballistic conditions with particles arriving normal to the target and deposit (c) with particles arriving at a 60” angle with respect to the target normal.

the help of a local collision

energy balance (Dahneke,

1971; Wall and John, 1988) that allowed for adhesive barriers in the normal and tangential directions. The

essential features of the current algorithmic model, which gives microstructures similar to the ones generated by the impaction simulations (as determined. by similar final solid fractions, distributions of coordination numbers, radial distribution functions, etc.) are the following: l

l l l

We perform three-dimensional off-lattice simulations, with periodic boundary conditions. The incident particles are started at random points above the existing &posit, and are assumed to travel on straight lines always parallel to the xzplane and inclined at an angle 8 from the normal to the target &y-plane). The particles are hard spheres, sampled from any appropriate particle size distribution. Already deposited particles cannot be displaced by subsequent particle arrivals (frozen limit). A particle that hits the target surface sticks immediately. If it hits another particle, then it rolls down the plane defined by the line along the two particle centers and the velocity vector (coinciding here with the original direction of incidence) of the moving particle. The rolling motion is continued until it encounters a new (third) particle, in which case it roils down the new plane defined by the velocity vector and the new contact normal. Each new contact between the rolling particle and another fixed particle, is considered as a new “rolling event”. The rolling is continued until the number of rolling events specified at the beginning of the simulation is completed or the particle reaches a position of local minimum potential energy, or it reaches the original target surface.

In Fig. 2 we show three two-dimensional, off-lattice deposits, generated by particles that arrive normal to the target and undergo one rolling event (deposit 24, two (deposit 2b), or three rolling events. The effect of the rolling on the deposit microstructure is apparent, especially when compared to the pure ballistic deposit of Fig. lb. Note further, that the difference in the final solid fraction, #J, between deposits 2b and 2c is very small, an. effect we also observed in the 3D simulations. On the other hand the high degree of crystallinity apparent in these 2D deposits is completely absent in three dimensions, where cubic or hexagonal configurations are rare. This is demonstrated by the radial distribution functions shown in Fig. 3. The radial distribution function, g(r), gives the probability of finding a particle center in a small volume element about r given that there is a particle at the origin, and it is usually normalized with respect to the average particle number density so that g(r) + 1 in the limit of larger r. Figure 3(a) shows the radial distribution function in the z-direction (norma to the target surface), gZ. for a ballistic deposit (solid line) and a denser deposit (dashed line) formed by allowing each particle to roll to the position of minimum local potential energy. Both deposits were generated by particles arriving normal to the target. The similarity of the two curves suggests that the relative short-range positioning of the particles in the z-direction is independent of the local porosity. Note further the absence of peaks normally associated with simple and body-centered cubic solids. In Fig. 3(b) we plot the corresponding distribution functions in directions parallel to the target (q-plane), gxg. Here we have averaged over the azimuthal angle since deposits formed by particles arrivingsormal to the target possess azimuthal symmetry. For the denser deposit g, and gXr are similar, indicating that these deposits are nearly isotropic. On the contrary, the ballistic deposit gXr is flat, i.e. there

MEUELA~~ TASSOFQULOS and DANIEL E. ROSNER

(4

(4 Fig. 2. Typical two-dimensional deposits with particle restructuring. All deposits were generated by particles following straight line trajectories, normal to the target. In deposit (a) particles undergo one “rolling event” after contact with an already deposited particle, in deposit (b) two, and in deposit (c) three rolling events. is absence of short range order in planes parallel to the target surface. We should also point out that similar restructuring mechanisms have been proposed in the past, either in the context of the statistics of random packings of spheres and their association to the properties of hard-sphere liquids (Bernal, 1964, Finney, 1970; Bennett, 1972; Visscher and Bolsterli, 1972; Reyes and Iglcsia, 1991), or in relation to deposits formed by the chemical vapor deposition process (Henderson et al., 1974,.Dirks and Leamy, 1977) and by the sedimentation of spheres in a viscous fluid (Meakin and Jullien, 1987; Jullien and Meakin, 1987). The main difference between our model and previous ones, evident only in three dimensions, is in the motion of the rolling particle after it establishes contact with a third particle, The previous models, trying in most cases to maximize the final solid fraction assumed that the moving particle after the completion of what we term a “rolling event” would roll in the direction of steepest descent in contact with borh new particles. This may

be a plausible scenario if one is modelling surface diffusion or packings of cohesionless particles, but is not supported by the “full” impaction simulations mentioned earlier. In Fig. 4 we plot deposit solid fraction vs number of rolling events for different angles of incidence. The empty symbols correspond to solids generated with the current restructuring mechanism, while the tilled symbols are obtained using the restructuring model of Visscher and Bolsterli (1972), see also Jullien and Meakin (1987). As expected, our model results in less dense deposits, especially in the limit of many rolling events. On the other hand, the dependence of the final deposit solid fraction on the angle of incidence is similar between the two models, in the sense that the higher the angle of incidence, as measured from the target normal, the more open the deposits are. For comparison, we further indicate on Fig. 4, the local solid fraction of a deposit grown under pure diffusion limited (DLA) conditions. This solid fraction corresponds to a distance of about 40 particle diameters away from the target.

Simulation of vapor diffusion in anisotropic particulate deposits

427

3.0

. 0.55

sc

2.0

-----_

L c _

1

1

.

2.5

0.85 0.55

gt.5 3 1 .o 0 7

0

0.25

0.5

0.0 0.0

1.0

2.0

.

3.0

Angle of incidence CCOCQ 0 dgr OmM 30 dgr

4.0

3.0

Fig. 4. Effect of number of rolling events on deposit solid fraction. Variation with angle of incidence and restructuring model used.

7

2.5

rank tensor proportional 2.0

----_-

L * _

to (Satake,

1985):

F -

0.85 0.55

“0

(2)

where II denotes the contact normal vector (the angular brackets in eq. (2) represent an ensemble average).

51.5 3 1a

3. ESTIMATION 0.5

0.0 @I

OF DEPOSIT

EFFECTIVE

DIFFUSI-OMPUTATIONAL

at.0

-

3.0

1

:‘“/

D,

a 4.0

5.0

Fig. 3. Radial distribution functions along the direction normal to the target, g,(r), and parallel to the target, Q_, of a three-dimensional ballistic deposit (# = 0.15) and a denser deposit (4 = 0.44) obtained by allowing the particles to roll to the position of local minimum potential energy.

The angle of incidence of the arriving particles affects the final deposit microstructure in two ways, namely by increasing the deposit porosity, i.e. Fig. 4, and by determining the principal directions of anisotropy. For example, CVD experiments (see, e.g., Dirks and Leamy, 1977) indicated that deposits formed by particles that arrive at a relatively large angle with respect to the target normal, say 8,, have an appropriate defined mean direction of growth, which was initially approximated by the “tangent rule”, tan 0, = 3 tan Bi

(1)

where t?# is the angle formed by the direction of growth and the target normal. Although a complete description of an arbitrary anisotropic microstructure may require a representation in terms of an infinite series of increasing-order tensors (see e.g. Onat and Leckie, 1988), the principal directions of anisotropy of deposits formed by discrete unconsolidated particles are usually defined to be the eigenvectors of thefnbric tensor. The fabric tensor, F, is a symmetric second

METHODS

3.1. Test-Molecule method The test-molecule method was proposed by Evans et al. (1980) and we refer the interested reader to the original paper for a more compIete description. Suppose that we want to compute, D,,, that is the effective diffusivity in the z-direction. Consider a slab of thickness L, in the z-direction, and “infinite” (by using periodic boundary conditions) in the x- and y-directions. The method is based on generating many molecules on one side of the slab and estimating the fraction transmitted, fr, through the porous slab to the other side. The molecule trajectories within the slab are determined by the local Knudsen number, and it is further assumed that molecules undergo diffuse reflections when impinging on the solid; the law of diffuse reflection is consistent with microscopially rough particles. D,, is then obtained from the equation (Evans et al., 1980): D,, =

lim ‘A r.r-oo 4

E

(3)

where E = (8RT/xM)‘12 is the mean thermal speed of the gas molecules, T is the temperature, R is the universal gas constant and M is the molecular weight of the gas. According to eq. (3) the effective difhtsivity in any one direction can be obtained only in the limit of large penetration L, or, equivalently, by using a sufficiently wide “porous slab”. In practice we satisfy this requirement by using appropriate boundary conditions at the edges (boundaries) of the sample. In this study we have consistently used periodic boundary conditions

428

MENELAOS TASSOPOULOS and DANIEL E. ROSNER

in the x- and y-directions, since the deposits themselves were generated with periodic boundary conditions. At the sample edges in the z-direction we have tried both periodic and specular reflection boundary conditions (Burganos and Sotirchos, 1988, 1989a). A specular reflection boundary condition is equivalent to replacing the infinite porous medium by identical cells arranged in such a way that any two are mirror images of each other. When considering porous samples that are not inherently periodic, e.g. our samples in the z-direction, both approaches have advantages and disadvantages. Periodic boundary conditions preserve the direction of anisotropy of the void phase, but they also lead to an apparent decrease of the void area at the cell edges; thus while for any random-plane cut through the sample the area fraction free of solid is on average equal to E, the two planes that define the cell edges have a void coverage equal to 6’. Specular reflection boundary conditions, although not having the latter problem, they do not preserve the direction of anisotropy of the sample, unless the direction of anisotropy coincides with the direction in which the specular reflection boundary condition is applied. We have tried both types of boundary conditions on samples obtained from deposits formed by particles arriving normal to the target (note that in this case the specular reflection boundary condition is applied in the z-direction which is also the direction of anisotropy). We did not see any statistically significant difference in D,, obtained by the two methods, probably due to the fact that we deal with relatively hi& porosity beds of spherical particles. Since one of the goals of this study was to determine anisotropic effects we have thus adopted periodic boundary conditions for the z-direction as well. In the original implementation of the test-molecule method, in order to obtain the effective diffusivity in all three directions, say JC,y and z, it was necessary to perform three independent simulations, one for each direction. Furthermore, one can argue (see, e.g. Burganos and Sotirchos, 1989a) that because most of the molecules penetrate only a short distance into the slab before they are reflected (exit through the “entrance” face), they do not adequately sample the whole deposit. It turns out, that both these Shortcomings can be easily circumvented if we apply the testmolecule method as follows: Consider a large cubic sample, obtained from the porous solid of interest, that is rendered “infinite” by application of periodic boundary conditions (PBCs) in all three directions. Generate test-molecules, one at a time, at random positions r,, in the void-space, and assign to them initial velocities vO, chosen again randomly. If L, denotes the maximum distance we allow the molecules to penetrate in the ith direction, then clearly a molecule will have successfully transmitted if at some time t, its ith coordinate satisfies

while it will have been reflected if 5

c%(t) -

xm1 < 0.

By following each molecule until it is transmitted/reflected in all three directions, a single run yields the transmission probabilities associated with all three directions, with significant savings in computation time. 3.2. Mean-square displacement method The mean-square displacement method is based on the fact that the effective orientation averaged diffusivity of a porous medium is given by:

(see, e.g. Chandrasekhar, 1943; for an application in the present context see, also, Burganos and Sotirchos, 1989a, b) while the di&sivity in the ith direction by D,

=

lim

(st>

where denotes the mean-square displacement of a large number of molecules, the ith component of the mean-square displacement and t the travel time. Here we should point out that eqs (4) and (5) as given are valid for regular Brownian diffusion, i.e. when the diffusing molecules (random walkers) take steps of constant length (in our context this is equivalent to assuming that the step size at each discrete time is equal to the mean-free-path). On the other hand, for gas kinetic theory diffusion, that is, when the freepaths are exponentially distributed, cf. eq. (12), the right-hand side of eqs (4) and (5) must be multiplied by a factor of 2. In a fashion similar to the extended testmolecule method, molecules are generated at random positions in the void space, they are assigned random initial velocities, and from their trajectories and , i = 1, 3 are computed. The effective ditTusivities are then determined from the slope of mean-square displacement vs travel time plots, in the limit of large time. 3.3. Anisotropic difisivity and “rectification” effects The porous solids considered here can be quite anisotropic depending on many factors, including the angle of incidence of the arriving particles and subsequent particle restructuring. In Section 2, we presented one way of determining the principal directions of anisotropy in the solid phase, from the distribution of the contact normals, . On the other hand, vapor diffusion occurs through the voids, and so we expect that the eigenvalues and eigenvectors of , which provide one descriptor of the anisotropy of the solid phase, will not necessarily coincide with those of the diffusivity tensor, D,, . Moreover, it is obvious that anisotropic effects will be more pronounced in the high Kn regime (see, also, Satterfield, 19701, since there molecules move in straight lines

429

Simulation of vapor diffusion in anisotropic particulate deposits between successive collisions with the solid phase and so “blocking” and “streaming” effects should he maximixed. This is further suggested by effective-medium type analyses, which in the continuum limit predict that solid-solid interactions do not enter to first order with 4, while in the Knudsen liiit they must be always accounted for, cf. eq. (10). In order to determine in a more quantitative way the effect of the anisotropy of the solid phase on Den as a function of local Knudsen number, we propose here a way to determine the complete diffusivity tensor. If the principal directions of the diffusivity tensor are a priori known, we can choose a coordinate system coincident with these directions, so ‘that the D,,‘s determined from either the test-molecule or the meansquare displacement method, yield directly the full diisivity tensor in diagonal form. In fact, the principal directions are known a priori in the case of deposits formed by particles arriving at normal incidence to the target; they coincide with the axes of the Cartesian frame of reference, with properties along x and y being the same. When particles arrive at an angle with respect to the target normal there are no symmetries that fix the principal directions, and so they must be determined by diagonalizing the complete ditIusivity tensor. The off-diagonal components of 0, D,, with i, j = 1, 3; i # j can be readily determined using the mean-square displacement method, if we keep track of the coordinates of the diffusing molecule with respect to three frames of reference. Of course, since D is symmetric, D, = D,,. Suppose, for example, that we want to determine Dlz_ Let (x, y, z) be the original frame of reference, and construct a new frame (x’, y’. z’) by rotating the original one about the z-axis by an angle B in the counter-clockwise direction. During the simulations we record the molecule position with respect to both frames without sign&ant computational overhead, since given the coordinates in the unprimed frame, the coordinates in the primed one are immediately obtained as a function of cos 8 and sin 8. Thus, at the end of the simulation we determine Dli, Dz2, and the corresponding primed quantities, 0; 1, Dz2. We further know, that the diffusivity tensors with respect to the two frames of reference satisfy: U = Q=DQ where Q is the rotation matrix associated with the rotation described above. From this relation, solving for D,,, we find that D 12 = =

D;,

- D,, ~0~~8 -D,, 2cosBsin0

sin2 e

-~D;,+D,,~os20+D2rsin28 2cos0sin8

(6)

By considering similar rotations about the x and y axes we compute D,, and D,, thus constructing the full diffusivity tensor. Another point of interest, again illustrated by the particular porous solids considered here, is whether

the direction of diffusion, namely towards the target or away from the target, has any effect on the resulting transmission probability, fr. This was suggested to us by the often “tree-like” shape of the deposit microstructures, cf. Fig. 1, indicating that under appropriate conditions the short-time effective diffusion coefficient in directions away from the target might be higher than Dcff in the opposite direction. We have explored such possible “rectification” effects using the testmolecule method, which can inherently distinguish between opposite directions. 4. CODEDRSCRIPTION

AND

VERIFICATION

4.1. Knudsen limit In the Knudsen regime, the molecule mean-freopath (MFP), i.e. the average distance a gas molecule travels between successive collisions with other molecules, is much greater than the characteristic length scales of the deposit, such as the mean particle diameter or an appropriately defined mean pore size. It is therefore accurate to assume that the diffusing molecules move in straight-line trajectories between successive collisions with the solid. From a code implementation point of view, this type of motion can be efficiently followed in a computer because, given the current position and velocity vector of a molecule, the time, t,, required to reach the ith particle located at some known position can be conveniently obtained from the solution of a quadratic equation. Of course, if there is a list of potential collision partners, the molecule will actually collide with the particle it will reach first, that is the particle with the smallest positive ti. In order to further accelerate the simulations, the computational domain is subdivided in cubic cells of size equal to the diameter of the largest particle present, so that at any point during an individual run, we search for the next collision partner out of a very small number of particles. Depending on the solid fraction of the deposit considered and particle polydispersity, at any point during a molecule trajectory we may have to consider 10 to 50 potential collision partners, while the total number of particles in the deposit can be up to 4O,fK!O. Throughout the simulations reported in this communication, the length scale is set by the particle radius, R,, and the time scale is R,,G, where F is the gas mean-thermal-speed. In Fig. 5, we plot meansquare displacement vs time of travel (s being the total distance travelled) for molecules diffusing in the Knudsen limit. The three curves correspond to the same solid, but were obtained by varying the total number of molecules considered in each run (see the figure legend), and the total diffusion time (travel time) allowed to each molecule. Note that in all cases we obtain the expected linear behavior, and furthermore the diffusivities computed from the slopes of the three curves agree to within 4%. We found that in most cases it was adequate to consider 2,000 molecules per run, with each molecule diffusing for 5000 to 8000 dimensionless time units. The statistical uncertainty in the results reported in this communication ranges

T~sso~ou~osand

MENELAOS

DANIELEROSNER

results on the Knudsen effective diffirsivity for porous media are often presented in terms of the dimensionless torruosity factor. The tortuosity factor, 7, is defined by the equation:

&L#S

D CFf

10000

C

D

Fig. 5. Mean square. displacement vs time curves obtained for Knudsen diffusion.

0.50

.

I

.

.

.

.

.

I

.

I

.

L

0.40

.

=

,

.

I

.

I

I

0.1s

-

=

-

7

where E is the porosity, D,, is some reference Knudsen diffusivity and D,, is the actual el?ective diffusivity. The tortuosity factor is essentially a correction to the reference Knudsen diffusivity, DKm, that accounts for the on-average longer travel path of a molecule when compared to the actual distance between two points along.the flow, constrictions of the pores, dead ends, etc. Although there is no agreement in the open literature on the most appropriate definition of the mean pore diameter and hence the reference Knudsen diBusivity, .especially when dealing with solids that have broad pore size distributions (see, e.g., Melkote and Jensen, 1989 for a discussion), for beds consisting of unconsolidated spherical particles the mean pore diameter, db, is usually defmed as four times the ratio of pore volume to pore area (Dullien, 1979, see, also, Section 5.1). Note, that if the bed further consists of equal size particles of radius R,, then

Given &,,, the reference Knudsen to be

D

Fig. 6. Fraction of molecules transmitted vs inverse of penetration distance. Effect of deposit porosity and direction of diffusion. between 3 and lo%, depending on the sample solid fraction and local Knudsen number. The anisotropic results have the highest uncertainty; the values for D, given represent averages obtained by performing the diffusion simulations on 3 to 5 different deposit realizations. In Fig. 6, we plot typical curves of fraction of molecules transmitted vs reciprocal penetration distance, corresponding to samples of two different porosities. Note, that for the higher porosity sample linearity is attained only after high penetrations, the final results being quite sensitive on the portion of the curve used to determine the slope. Because of this uncertainty associated with the testmolecule method, we have primarily used the meansquare displacement method to determine D,, , while the former was used to assess possible “rectification” effects. In order to verily the code in the Knudsen regime, we compared our ditIusivities for regular packings of spherical particles with experimental data obtained by other investigators. In the engineering literature

=”

-

ditfusivity is taken

fEFQ

(9)

the Knudsen diflusivity in a straight cylindrical pore of dieter equal to the mean pore diameter. We have performed Knudsen diffusion simulations in simple cubic (SC) solids and body-centered (BC) cubic solids, and have determined the associated tortuosity factors. In Table 1, we present the tortuosity factors we computed for these regular arrays, along with related experimental data obtained by Olague et al. 1988. Note the relatively good agreement between the various tortuosity values, especially when realizing that the real solids considered in Table 1, although of similar porosity (the lowest order microstructural descriptor), are combinations of simple cubic and bodycentered packings. A perhaps more stringent test of the Knudsen diffusion code can be devised by using it to directly compute a theoretically known bound. Specifically, Strieder and Prager (1968) have applied variational principles to formulate a rigorous upper bound on the effective Knudsen diffusivity of a porous medium. They have shown that for Knudsen diffusion in an isotropic material the following inequality is always satisfied: Deff (

1

6 ‘T

b*>

-

q

z <(,,r _ ,,).@,f _ ,,)>

=

Deft.

rmr

(10)

Simulation of vapor diffusion in anisotropic particulate deposits

431

Table 1. Comparison of tortuosity factors obtained in this work for cubic arrays of particles and some experimental values reported for regular arrays of spheres

Solid Simple cubic Body centered Regular array Regular array Regular array

&

‘F

Source

0.476 0.320 0.339 0.356 0.348

1.25 1.56 1.68 1.95 1.75

This work This work Olague et al. (1988) Olague et 01.(1988) Olaguc et al. (1988)

Table 2. Evaluation of an upper bound on the effective Knudsen diffusion coefficient for packings of mndom interpenetrating spheres, and comparison to the theoretical value E

0.90 0.70 0.5

Code [eq. (lo)]

mry

11.15 2.82 1.20

Ceq- (lU1 11.08 2.87 1.23

where m = (r’ - r) is the displacement of a molecule between successive particle collisions, II and n’are the local normals pointing towards the void space at r and r’ respectively, and S, is the total pore-surface area (including the area associated with disconnected pores) per unit volume. Next, for a pore system generated by randomly overlapping spheres they were able to analytically evaluate the ensemble averages appearing in eq. (lo), to obtain the dimensionless upper limit of the effective diffusion coefficient:

Deff,mar

48 =

-

&= --

13 s,

(11)

If we further note, that in order to generate the trajectories of the diffusing molecules, we compute both m and n, II’(the local normal is used to determine a new velocity after a diffuse reflection), it becomes apparent that the ensemble averages of eq. (10) and thus the upper bound itself can be evaluated directly from the simulation. We have used the Knudsen code to numerically evaluate the upper bound, D.ff mlx, for porous samples generated by random addition of monodisperse overlapping spheres, and in Table 2 we compare our results with the predictions of the theory, cf. eq. (11). Note that the agreement between the two is excellent, further verifying the Knudsen code. 4.2. Continuum limit In the continuum limit the molecule MFP becomes very small, and so the molecule trajectories can be represented by stochastic random walks, of ideally vanishing (A -, 0) step size. By letting many molecules “wander” through the porous medium and in a manner completely analogous to the case of Kn + 00, the slope of the mean-square displacement vs time curve yields Dem. As we have already indicated in the introduction this approach has been used quite exten-

sively to accurately determine various transport prop erties of porous materials (see references in Section 1). Implementation of the random walk simulations has been approached by two different ways: (1) each molecule performs a complete random walk, sometimes also termed discrete random walk, by generating a random direction, taking a step of constant length in that direction, generating a new direction isotropically distributed, taking another step, etc. When dealing with multiphase materials where every phase has its own conductivity, this discrete approach can be computationally intensive since before taking a step one must check whether this new step will lead the random walker into a different phase. This type of searching is potentially the most CPU-intensive part during off-lattice simulations. In the present context, where the solid phase is assumed nonconducting, if it is found that the next step will lead to solid penetration it is customary to use the so-called “blind-ant” boundary condition, according to which the molecule is not moved from its current position but the clock counter, which keeps track of the total number of steps taken so far, is incremented by one. Both the “blind-ant” boundary condition and the law of diffuse reflection satisfy the constraint that the normal flux at an interface between conducting and insulating regions must be zero. (2) use a continuum, sometimes also called Pearson or “floating” random walk (Brown, 1956), to generate the diffusive motion of the molecule. According to this method, the next location of the diffusing molecule is determined as follows: (a) inscribe the imaginary sphere that is centered at the current location of the molecule, and tangent to the closest solid particle, (b) select a random point on the surface of this imaginary sphere and move the molecule to this new position, and (c) update the time counter either by sampling the actual first passage time distribution, i.e. the distribution of times taken by a molecule starting at the origin of the sphere to reach its surface as determined from the solution of the diffusion equation (Siegel and Langer, 1986; Zheng and Chiew, 1989) or by assuming that the molecule always takes the mean time (R2/Do for a step size equal to unity) to reach the sphere surface (Torquato and Kim, 1989), where D, the intrinsic diffusion coefficient in the conducting medium and R the radius of the inscribed sphere. It turns out that the computed D,, is insensitive to whether the time counter is incremented at each step, by always using the

432

MENELAOS TASSOPOLILOSand DANIEL E. ROSNER

mean first passage time, or by sampling the first passage time distribution (see e.g., Siegel and Langer, 1986). A valid concern when using the solution of the diffusion equation to determine the time required to reach the surface of the inscribed sphere (see e.g., Reyes et al., 1989), is that the diffusion equation is valid only in the limit of R/A + 00, where R is the sphere radius, and 1 the fixed step size (equivalent to the mean-free-path). It turns out that for R/A E O(1) the theoretical mean first passage time N R’/l.) underpredicts the actual one, (_TTThWV as indicated from Fig. 7, where we plot the T Simulrtiall~ ratio of the two vs dimensionless sphere radius. Fslmu,a,l,,awas obtained by averaging the time 100,000 diffusing molecules took to reach the surface of spheres of radius R/A. In the continuum limit (Kn -D 0) simulations we have used a hybrid technique, where the trajectories of the molecules are determined from discrete random walks when closer than four mean free paths to the solid particles, and from the inscribed sphere technique using the mean first passage time otherwise (Reyes et al., 1989). We have also performed discrete random walk simulations, that is, without using the concept of the first passage time (see Section 4.3 for further discussion and a description of an efficient algorithm to perform such simulations) with a step size of the random walk motion in the range 1O-2-1O-3 particle diameters (the corresponding Knudsen numbers are lo-’ and lo-‘, respectively) with no statistically significant effect on D,, (see, also, Fig. 8 below). A similar conclusion was reached by Schwartz and Banavar (1989), and so we have used a the “continuum” value of Kn = 0.01 throughout simulations. In Fig. 8, we compare effective diffusivities obtained from discrete random walks to Doa’s obtained from the Pearson random walk. Note, the exceptionally good agreement between the discrete and Pearson random walk simulations. Our main interest being vapor diffusion through the gas phase in the pores, we have also considered

P

IF \

s

‘2

S i

Il-

‘Z

Fig. 7. Ratio of computed to theoretical mean first passage time vs dimensionless sphere radius.

Full

Random

Walk

((D)

x

10’)

Fig. 8. Comparison

of effective diffusion coefficient obtained using a discrete random walk approach abscissa and using the Pearson (continuum) random walk.

diffusion with a variable free-path. In a uniform gas that consists of rigid elastic spherical molecules that move at constant velocity, the probability that a molecule will traverse a distance greater than 1 is eerja, where L is the mean-free-path. More elaborate calculations that account for a Maxwellian velocity distribution provide nearly the same result (see e.g., Chapman and Cowling, 1960). If we further assume that the exponential free path distribution is not altered by the presence of the rough particles (there exist both experimental indications and theoretical arguments that suggest that this might be the case; see Fowler, 1955 for a discussion) then the free-path between successive molecul~molecule collisions is given by: I = I In (l/U)

(12)

where U is a random variable uniformly distributed between zero and unity, U = U(0, 11. Using eq. (12) to sample the free-path distribution, we found that in the continuum limit, with Kn Q 0.01, assuming a tied or variable free-path did not have any statistically significant effect on the estimated diffision coefficients. A similar conclusion was reached by Abbasi et al. (1983). In Fig. 9 we plot mean-square displacement vs time obtained from continuum simulations performed on the same solid and involving different total number of molecules. Note that in all three cases the data points lie on a straight line. Most of our data where obtained from runs of 2000 molecules, with each molecule diffusing for a total dimensionless travel time equal to 8000. For Kn = 0.01, this travel time is equivalent to taking 400,000 discrete steps per molecule. In the case of continuum diffusion we are able to verify the computer code by comparing our simulation results to theoretical results, available for regular arrays of cubic solids. In Fig. IO we plot formation factor vs porosity for simple cubic solids. The formation factor, sometimes also termed electrical resistivity is defined as the ratio of the factor, F( = k,/k,,), conductivity of the pore filling fluid to the effective

Simulation of vapor diffusion in anisotropic particulate deposits

40-

a*--G-rooo. .mmm- N”,.p2000. *.a.. #q&-JOOD.

D -3.2zE-03 D -3.15E-03 D -3.25E-05

3D-

0

&

Fig. 9. Mean-square displacement vs time curves obtained for continuum diffusion (Kn = 0.01).

Fig. Il. Variation of formation factor with deposit porosity for body-centered cubic solids. Comparison between theory and simulation results.

cubic lattice had a size equal to 2 x 10m3, in R, units, corresponding to Kn = 10m3). Note, the very good agreement between the predictions from the simulations and the theory. In Fig. 11, we plot formation factor vs porosity for body-centered cubic solids. The solid line corresponds to the theory (Sangani and Acrivos, 1983) while the data points are obtained from simulations involving the hybrid technique. Again the agreement between theory and simulation is very good.

0.5 1 0.25

I 0.50

0.75

1 1.oo

&

Fig. 10. Variation of formation factor with deposit porosity for simple cubic solids. Comparison between theory and simulation .results.

conductivity of the two phase material. Note, that keff is related to the effective diffusion coefficient through Einstein’s relation, namely k - DB.s. For diffusion in free space, the dimensionless bulk diffusion co&icient, D,, directly proportional to k, since the porosity is unity, is equal to Kn/3 for a random walk of fixed step length, and equal to 2 Knj3 for a “kinetic theory” random walk. The solid line in Fig. 10, corresponds to the theoretical results obtained by McPhedran and McKenzie (1978) who extended a technique originally developed by Lord Rayleigh; according to their method they compute the conductivity associated with a unit cell by taking into account tbe effect of multipoles of arbitrarily high order. The symbols in the figure correspond to data obtained from: (a) discrete random walk simulations using a variable freepath, (b) off-lattice simulations using the hybrid technique (discrete and Pearson random walk) and (c) some on-lattice simulations we performed on a discretized unit-cell of the cubic array (each cell of the

4.3. Transition regime If we again-invoke a description consistent with the simple kinetic theory of gases, the free-path distribution in the transition regime remains exponential, cf. eq. (S), with the mean-free-path, J.,related to the local Knudsen number through Kn = L/d,, d, being the reference particle diameter; d, = 2 in our simulations. It is known @at in the transition regime, both the effective mean-free-path of the diffusing molecules and the resulting effective dXusivities depend on moleculemolecule and molecule-particle interactions flokunaga, 1985; Scott and Dullien, 1962), tith the molecule-particle interactions becoming more important with increasing Kn. Thus, during the simulations in the transition range it is necessary to allow for a variable free-path, in order to properly account for the interactions between the molecule and the solid phase. This can be readily seen if we consider diffusion of a gas, molecule through a tube of fixed diameter (see Pollard and Present, 1948). In the limit of Kn + 0, where in the present context Kn is defined with respect to the tube diameter, assuming fixed or exponentially distributed free-paths has no effect on D cff since in both cases the diffusing molecules “see” the tube walls only rarely. On the contrary at intermediate Kn values molecules diffusing with an exponentially distributed free-path will see the tube walls more often than molecules taking fixed-length steps and so they will approach the Knudsen limit diffusion faster. This asserted behavior has been also confirmed

MENELAOS TASBXVUL~~ and DA-L

434

by our simulations, which predict a smooth transition between the continuum and Knudsen diffusion limits only if the free-path is allowed to vary. In the transition regime we use the following reference diffusion coefficient, D,, to define the formation and tortuosity factors (see also Section 4.1): 1 -=

Do

‘+& DK,

(13) Br

where D,, the reference Knudsen dflusivity, cf. eq. (9), and DB the dimensionless reference self-diffusion coefficient, which in the context of a “gas kinetic theory” random walk in free space is equal to, D, = 2Kn/3. Scott and Dullien (1962) derived eq. (13) based on momentum transfer considerations, and have shown that it is strictly valid for self-diffusion or equimolal counterdiffusion of a gas at constant pressurti. This “harmonic average” form of D,,, can also be derived from probabilistic considerations on the effect of the porous medium on the free-path distribution (Tokunaga, 1985). In Section 4.2 we have already indicated that offlattice slmulations of discrete random walks are computationally very expensive, and so an &bent algorithm is required. It turns out that such an algorithm can be obtained as follows: Consider a molecule, which at some initial time is located somewhere in the void space. Construct the largbt iriscribed sphere that does not overlap any particles;in a manner completely analogous to the Pearson algorithm. Next generate a random direction, sample the free-path distribution to obtain the length of the next step and move the molecule accordingly. If the new molecule position remains within the inscribed sphere, it is autoniatically guaranteed that no pcnetration has occurred, and the molecule can take a new step. If at some point a step takes the molecule outside the inscribed sphere, then before moving the molecule, we check whether a ptirticle blocks its trajectory. If this is the case, then it moves to the particle surface undergoes a diffuse reflection and takes a ntw random step; otherwise it moves to the new position, we construct the new largest inscribed sphere and the whole process is repeated. In Table 3, we compare the computing (CPU) time 1 to follow 100 required by a Sun Sparcstation Table 3. Typical code execution Algorithm

E. ROSNES

molecules, each one diffusing over a total of 5000 dimensionless time units, when using the computer codes developed for this study. We present typical execution times of the code that performs diffusion simulations in the Knudsen limit (Section 4.1), of the “hybrid” algorithm based on first passage time concepts (Section 4.2), of the efficient discrete random walk algorithm described here, and of an “inefficient” discrete random walk algorithm where we check for possible collisions with the neighboring particles at each step. In the continuum simulations it was assumed that Kn = 0.01 and data are reported for two sample porosities. Note that all codes become slower with decreasing deposit porosity, with the “Knudsen code” performance being the most sensitive. Note further that the new efficient discrete random walk algorithm, which can accurately treat the transition regime (within the inherent model assumptions), is always slower than the mean first passage time code, but still up to 5 times faster than the inefficient discrete random walk algorithm that checks for possible collisions with the neighboring particles at each step. 5. RESULIX AND DISCUSSION 5.1.

Orientation avernged results

We have analyzed porous samples obtained from deposits that were formed by spherical particles arriving with angles of incidence equal to 0”, i.e. normal to the target, as well as 30,45 and 60”. The resulting solid fractions, 4 ( = 1 - E), along with subsequent particle restmcturing effects (rolling events) have been shown in Fig. 4. The effective diffusivity values reported here represent averages obtained from 3 to 5 different samples with a resulting statistical uncertainty between 3 and 5%. In each independent run we typically followed 2000 mokcules diffusing for 6000 to 8000 dimensionless time units. In this section we will present the orientation averaged results, while anisotropit and rectification effects will be explored In the following subsection. Given the components of the diffusivity tensor along the three principal directions of anisotropy, D,,, the orientation averaged diffusivity can be evaluated from:

= h-&f

times’ for the various algorithms study

’_X

E

Knudsen wde, Kn + cn Knudsen code. Kn + 00 First passage &me, Kn = 0.01 First passage time, kn = 0.01 Efiicient discrete random walk, Kn = 0.01 E%Ecient discrete random walk, Kn = 0.01 Inefficient discrete random walk, Kn = 0.01 Inefficient discrete random walk, Kn = 0.01 tCPU used by a Sun Sparcstation dimensionless time units.

X:K 0.844 0.558 0.844 0.558 0.844 0.558

2 per 100 molecules

= %h,

used in this CPU (min) 6.3 37.3 5.5 14.4 48.0 66.5 159.9 321.6

Musing

SW0

(14)

Simulation of vapor diffusion in anisotropic particulatedeposits where we use Einsteins summation convention, and I is the unit tensor. We may further define an orientation averaged tortuosity factor (for an anisotropic medium the tortuosity factor itself is a second rank tensor; see, e.g., Dullien, 1979), (T),

(15)

2.0

*

435

. Knudsen

1.8

Diffusion

_

-

T

where, as before, D, is the reference diffusivity in the absence of the solid, calculated at the appropriate Knudsen number. In the Knudsen limit we define the reference diffusivity with respect to a circular tube of diameter equal J 1.0 1 to the hydraulic pore diameter. The hydraulic dia0.7 0.E 0.9 0.4 0.5 0.6 0.3 & meter, &,( = 4V/S), is defined as four times the inverse of pore area per unit pore volume. From stereolFig. 13. Variation of tortuosity factor with deposit porosity ogy it is further known (see, e.g., Underwood, 1970) in the Knudsen limit. Effect of the local deposit microthat 4 V/S is also the “mean intercept length”, defined structure,as determined by the angle of incidence of the as the average unobstructed pore surface to pore arriving particles. surface distance. Evans et al. (f980) have pointed out that an estimate of the mean intercept length can be they are related to the macroscopic anisotropy, espeobtained from the molecule travel distance between cially apparent in the less dense deposits. Further successive collisions with the pore walls in the support to this claim is provided by the fact that we Knudsen limit. Indeed Fig. 12 shows the good agreeobserved a similar discrepawy for simple cubic solids ment between the estimated me_an pore diameter even at maximum solid fraction (4 = 0.524). Based on based on the hydraulic diameter, dp, n_oti _,, and the the molecule trajectories we measured a mean-pore one computed from the molecule trajectories in the diameter equal to 1.21, while 4V/S gives 1.31. Again, Knudsen regime, dpsMdemurrd. In Fig. 12, we further we believe that the larger mean-free-path of the moleindicate the sample porosities corresponding to the cules is due to the completely open channels associmean-pore diameters, cf. eq. (10). In the limit of very ated with the simple cubic packing, while 4V/S is high porosities we have consistently observed deviabased on a random isotropic two-phase material. tions between the predicted and measured pore diaNote, that for body-centered solids, where such open meters, the measured ones being slightly higher than channels are absent the agreement between the two the hydraulic diameter. We have considered samples was very good. up to 60 x 60 x 60 particle diameters, and have let In Fig. 13, we plot the orientation-averaged molecules diffuse upto 20,000 dimensionless time untortuosity factor [empty symbols), obtained in the its without any difference in the obtained results. It limit of Kn 4 co, for all the different deposits conthus seems to us that the deviations are not due to finite size or short diffusion time effects, but rather’ sidered here vs sample porosity, E. The respective angles of incidence are indicated in the figure legend. Let us first concentrate on the data points at the very right of the figure (highest porosity, labeled 1) corres, 20.00 0.60 0.75 0.86 0.66 0.90 ponding to ballistic deposits without subsequent particle restructuring. Note the significant differences between the tortuosities associated with these open deposits, almost independently of sample porosity; the porosity of these deposits varies from ca 15% for normal incidence to about 11.5% for particles arriving at an angle of 60” with respect to the target normal. This almost 30% variation in the effective diffusivity of essentially similar, and, for that matter, high porosity deposits, clearly demonstrates the importance of the detailed deposit microstructure on the resulting effective transport properties. It is obvious that, at least in the free molecule limit, sample porosity alone is not an adequate descriptor of the deposit 0 2 46 a 10 12 &m.omtico~ microstructure. The next set of data points (labeled 2) corresponds to deposits generated by letting 50% of Fig. 12. Comparison of mean-pore diameter obtained from the arriving particles to undergo one rolling event, the molecule mean travel distance between 8uccessive colpoints labeled 3 correspond to all particles undergolisions with the solid particles and the mean hydraulic ing one rolling event and data points 4 to all particles diameter.

436

MENEMOS

TASXPOULO~

undergoing two rolling events. Note that for higher deposit solid fractions, becomes a linear function of porosity. In fact, in this limit and independently of angle of incidence we find that (7) can be approximated by (7)

N 3.4 - 3E,

0.55 < & < 0.75.

(16)

In this case the angle of incidence just determines the range of porosities for which eq. (16) remains valid. A linear dependence ,of Knudsen tortuosity ,on sample porosity has been also reported by Evans et al. (1980). In Fig. 13 we have also plotted typical Knudsen tortuosities (filled symbols), obtained experimentally by Smith and his collaborators. The filled triangles correspond to random packings of spherical particles (C u 0.363 f 0,030; Huizenga and Smith, 1986) while the Nled circles correspond to spherical packings that exhibited a high degree of ordering (Olague et al., 1988). In both cases the packings where unconsolidated and the reference Knudsen diffusivity was defined with respect to the hydraulic radius. Note that it is again the actual channel geometry which determines the resulting resistance to Knudsen diffusion as opposed to the absolute porosity. In Fig. 14 we plot orientation-averaged tortuosity vs sample porosity in the continuum limit. The empty data points correspond to our simulation results while the filled ones to experimental data. The first observation to be made is that in the Kn + 0 limit and for not highly consolidated porous media the effect of deposit microstructure on resistance to bulk diffusion can be adequately described as a function of porosity alone. If the porosity is very low, say, due to consolidation, certain void regions may become disconnected, and then one must use the accessible porosity. In any event the solids considered here are in a regime where such percolation problems are completely absent. If we take into account the statistical uncertainty of our data (about 7%), we find that they are quite well

2.0 -

Continuum 1.9lC \

..

*

Diffusion

Hanhin-Shtrikman Loww Bound _x Archian’ Law. m-1.45 Brqgamon. 1935 -v cutie. I SW Hoaplchcl9.n. ,915 .99..wang Turner. 1979

_

et al.. tm4

1.6 -

1.4 -

1.2 -

Fig. 14. Variation of tortuosity factor kith deposit porosity in the continuum limit. Effect of the local deposit microstructure, as determirsed by the angle of incidence of the arriving particles.

and DANIEL E. ROSNER described by an empirical relationship that links sample conductivity and porosity, known as Archie’s law: k cfF

=

ak,e”’

(17)

where k, the conductivity of the pore filling liquid (the solid matrix is assumed insulating) and E the porosity. a and m are empirical parameters, that, in general, may vary between different porous solids. For not too dense packings of consolidated particles (.s in the range between 0.2 and 0.42) it has been found that values of Q = 1 and m = 3/2 predict very well expcrimentally obtained formation factors (Wong et nl., 1984; Guyon et al., 1987). As can be seen from Fig. 14, Archie’s law with m = 1.45 describes our orientation averaged data quite well, indicating that it remains valid in the high porosity regime as well. Although the microstructures considered here, and as indicated, for example, by the distribution of the contact normals, are anisotropic the continuum limit diffusion simulations yielded effective diffusivities that were nearly isotropic independent of angle of incidence and deposit porosity. We have verified this finding independently of the diffusion simulations by evaluating a two-point microstructural tensor which appears in a perturbation expansion of the conductivity tensor of a two-phase anisotropic material (Torquato and Sen, 1990, and references therein). Note, that a two-point tensor is the lowest-order microstructural descriptor that contains information about the local anisotropy. For the special cast of an anisotropic medium that possesses azimuthal symmetry, e.g. our “normal-incidence” deposits, Torquato and Sen were able to express this tensor [eq. (50) of their paper] in terms of a single scalar quantity, Q [eq. (Sl)], which involves a volume integral of the twopoint anisotropic correlation function, &(r) [for a homogeneous statistically anisotropic medium S,(r) gives the probability of finding two points with relative displacement r in the ith phase]. Q is defined in such a way that for a statistically isotropic medium it is exactly equal to l/3, while deviations from this value indicate anisotropy (a similar integral but in a different context arises also in the analysis of Onat and Leckie, 1987). We have evaluated Q and the associated anisotropic upper bounds on the effective conductivity for representative deposits. For a ballistic deposit (4 N 0.15) we estimated Q = 0.35 and thus computed the effective conductivity upper hounds (k,,,Yy/ko)m_ = 0.89 and (k,,/ko)_. = 0.88. Similarly, for a denser deposit (4 = 0.39) we found Q = 0.37 and CLJkOLx = 0.67 and (k,,/kO),., = 0.68. The fact that the measured Qs are almost equal to l/3 is consistent with the diffusion simulation result of a nearly isotropic D,, in the continuum limit. With these considerations in mind, we show on Fig. 14 the lowest permitted tortuosity for any given porosity, as obtained from the Hashin and Shtrikman (1962) upper bound on k,,, for a two phase isotropic material. For insulating inclusions the Hashin and

Simulation of vapor diffusion in anisotropic particulate deposits

Shtrikman

bound becomes

Note, that this is the best, i.e. most tight, bound that can be obtained given information on the relative conductivity between the two phases and the respective volume fractions. The fact that our computed tortuosities remain always above the HS bound is a further consistency check for our continuum simulations. On Fig. 14, we have also plotted the predicted T values based on Bruggeman’s (1935) effective medium theory, who found that in the limit of e -+ 1 k -=eff k,

3E -

1

-’ 2

It turns out that Bruggeman’s formula signiRcantly overestimates the tortuosities. We have also determined D,, for porous samples “extracted” from DLA deposits, but explicit results are not given here since these very high porosity solids (a B 0.95) followed the low solid fraction expressions, cf. eqs (17) and (18), closely. The filled squares correspond to experimentally obtained bulk diffusion tortuosities for beds of smooth uniform spheres (Hoogschagen, 1955, Currie, 1960), while the filled circles represent electrical resistivity measurements of slightly consolidated samples formed from three particle sixes. Note that for any given porosity, the consolidated samples yielded higher tortuosities, most probably because of the larger constrictions (smaller “pore throats”) formed between contacting particles. Finally the filled triangles are taken again from resistivity measurements of fluid&d beds of nonconducting spherical particles (Turner, 1976). Resistivity measurements of fluidized beds have been the only experimental data we where able to find on the resistivity/tortuosity of high porosity packings of spheres. Apart from the scatter in the limit of unit porosity, attributed to experimental difficulties in obtaining a constant measurement and/or anisotropic effects as E + 1, they follow quite closely Maxwell’s prediction (1873) which for insulating inclusions coincides with the bound of Hashin and Shtrikman. Typical variations of the orientation-averaged tortuosity with local Knudsen number are shown in Fig. 15 for various deposits. Note, that for the specific class of porous solids considered here and when using eq. (13) to detine the reference diffusivity in the transition regime, the estimated tortuosity varies with Kn, especially as the deposit solid fraction increases. In agreement with previous investigators we also find that the transition regime approximately extends over Kn between lo- ’ and 50. In Fig. 16 we show a scatter plot of the ratio of in the Knudsen and continuum limit vs porosity. The straight line shown is a linear fit to the data, which gives <‘LX.> N - 0.348 + 1.41,

0.55 6 & d 0.90.

(20)

Fig. 15. Typical variation of tortuosity factor with Knudsen number, for deposits of different porosities.

b.,iilil 0.6

0.7

0.6

0.9

&

Fig. 16. Scatter plot of the ratio of the tortuosities in the Knudsen and continuum liits vs deposit porosity for all the microstructures considered in this study.

Of course, we must stress that “correlations” like eq. (20), based on sample porosity alone are, at best, oversimpliReations, particularly for the high porosity regime, where Q” becomes highly dependent on details of the microstructure, 6. Fig. 16. In the engineering literature it is customary to treat the tortuosity factor as a constant independent. of Knudsen number. Indeed, there have been cases when z was measured to be approximately the same in the two Kn limits, provided that the right characteristic pore size had been used when defining the reference Knudsen difhusivity. For example, Huizenga and Smith (1986) measured effective diffusivities of various gases ditising through pellets of small spherical particles (e = 0.363 f 0.030) and using a reference Knudsen diffusivity based on the hydraulic diameter found that their computed tortuosity, ranging between 1.45 and 1.51, agreed well with bulk diffusion tortuosities of similar porous solids, reported earlier

438

MENELAOS TASSOPOULOSand DANIEL E. ROSNER

by Hoogschagen (1955, z = 1.42) and Currie (1960, z = 1.48). On the other hand, in another communication (Smith, 1986) they considered Knudsen and bulk diffusion in constructed pores and found that the tortuosity factors obtained in the two limits agreed only if they defined D,, with respect to the minimum pore diameter, which one can presumably obtain from mercury porosimetry or adsorption isotherm analysis. These considerations along with our new results clearly demonstrate that it is not possible to know (z priori which is the appropriate pore length scale that will give tortuosities that are independent of the local Knudsen number. Of course, a tortuosity factor that is independent of Kn is quite desirable from a practical point of view, because then any single diffusivity measurement at a convenient Knudsen number in conjunction with eq. (13) provide a complete description of gas diffusion through this particular solid. Scott and Duilien (1962), in order to test the formula they derived for gas diffusion at arbitrary Kn (es. 9, in their paper), performed certain experiments where they measured effective diffusivities associated with four types of industrial catalyst supports ranging in porosity between 0.614 to 0.347. Although they had an estimate of the mean-pore diameter obtained from mercury porosimetry measurements, which could be used to define a Knudsen reference diflusivity, they instead took the following approach: From the efiective diffusivity they measured in the continuum limit, D arscFI and the self-diffusion coefficient of the gas at the same pressure and temperature, D&, they estimated the bulk diffusion tortuosity (in their paper they actually present the analysis in terms of formation factors)

Then they assumed that the same tortuosity factor, cf. eq. (21), was applicable in the Knudsen limit, and so they obtained their reference Knudsen diffusivity from

cff the effective Knudsen diffusivity was measured experimentally, and z is the tortuosity factor obtained from eq. (21). Using this scheme, they were able to describe their experimental results very successfully. This approach provides a good fit to the experimental data since it essentially guarantees the correct behavior in the two Knudsen limits. We tested this method using our data, and it was again very successful. This is shown in Fig. 17, where we plot normalized effective diffusivities vs Kn for two deposits of different porosities. The data points correspond to the actual simulation results, the solid line corresponds to the estimate we would obtain if we knew the complete dependence of z on Kn, with the reference Knudsen diffusivity based on the hydraulic pore diameter, and the dashed line gives the prediction obtained when we define the reference Knudsen

where b,

Fig. 17. Normal&d effectivediffusivityvs Knudsen number. Comparison of the simulation data and the prediction based on the method of Scott and Dullien (1962).

diffusivity as suggested by Scott and Dullien (1962), which just requires the knowledge of Deff in the two limits of pure Knudsen and pure bulk diffusion. From these observations it seems to us that whenever D,, is known in the two limits, this last approach probably yields the most accurate estimate of the effective diffusion coefficient in the transition regime. Let us also point out in passing, that the characteristic pore diameter that we backcalculate given DRn, namely = 3D,/E, falls approximately in the middle be4. alf tween the value estimated from the hydraulic radius and the value we would have obtained from say mercury porosimetry, estimated here to be about 75% of 4V/S. A similar
Simulation of vapor diffusionin anisotropicparticulate deposits

D 33_

D 33

D 11

D Z1

N 1.25 k 0.04,

KR --+ m

439

2.0

deposits generated by particles arriving normal to the target (23a)

1.5

while in the opposite limit D AZ=-D11

0”

D33 N 0.98 f 0.05, D 22

Kn + 0

(23b) I.0

where Dir, i = 1, 3, the diffusivity along the ith principal direction. Note, that for deposits generated by particles arriving normal to the target, D,, would be parallel to the z-axis. The effect of the Knudsen number on the anisotropic properties of the diEusivity tensor is further shown in Fig. 18 where we plot 2D,,/(D_ + Dg,,),a measure of the anisotropy of D, vs for two deposits formed by particles arriving Kn. normal to the target and thus possessing azimuthal symmetry. For both deposits, b is isotropic in the continuum limit, i.e. the diffusivities along the three principal directions are equal. Anisotropy appears only in the transition regime which extends approximately between Kn = 10-l and Kn = 50, and becomes most apparent in the pure Knudsen diffusion regime (Kn > 1000). Note also, that for the denser deposit (triangles in Fig. 18) the difference in the diflusivity along the z-direction and along the x- and y-direo tions is smaller, since particle restructuring tends to block (till) the elongated pores formed in ballistic type deposits (Fig. 1). When particles arrive at an angle with respect to the target normal, the principal directions of D do not necessarily coincide with the cartesian frame of reference (x, y, z) or with the principal directions of anisotropy of the solid, and so they must be directly determined from the simulations (Section 3.3). In the following discussion we limit ourselves to Knudsen diffusion, and like before let D,, i = 1, 3, be the directional diffusivities along the principal directions. Without loss of generality, we take D,,

1.5.

-

I.3

, . .‘.‘.,,

b

,

D,,

B

, ,,.,.,,

_

Normal Incidence

-

-L AAbAA.3 -

(24)

DI,.

,

, ,,,,,,,

,

0.66 0.56

1

C?

+ v2 \

. 1.1 41

.

*.

A .

:

a

*

0.9 -

cl” c\I

, ,,fiTT

0.7 -

0.5

IO -

..".I ' 10 -1

"'..'I

. . ..'..I

10

cLu* 101

IL

Fig 18. Effect of Knudsennumberon the anisotropic prop erties of the difflMivity tcnPlor. CEB

47:2-K

Knudsen Diffu~on Empty Symbols: d - 0 !=DIad Symbdm: d - 60 “? l.S

0.6

0.7

(dar)

0.6

,

&

Fig. 19. Directional tortuosities vs deposit porosity for deposits formed by particlesarrivingat a normal incidence and at 60” with respect to the targetnormal.

We may also define the directional tortuosity, zii, by

(25) where DRn the rcferetice Knudsen diffusivity obtained from eq. (9). In Fig. 19 we plot tii. i = 1, 3, vs deposit porosity, E, for particles arriving normal to the target (empty symbols) and at a 60” angle with respect to the target normal (filled symbols). The former deposit possesses azimuthal symmetry as indicated by the equal tortuosities along the 11 and 22 directions. On the other hand, the deposit formed by particles arriving at a 60” angle is anisotropic with respect to all three directions. Note also, that the computed tortuosity along the 33 direction is less than unity, indicating that diffusion along this direction is “faster” than the reference D,, based on an isotropic material of equal porosity. We have also investigated the relationship between the principal directions of the Knudsen diffusivity and the fabric tensors (cf. eq. 2). Representative results are provided in Table 4, where we report the angle (degrees) formed between each principal direction of F and of D with the Cartesian axes x, y and z. These angles are denoted as a, /3and y respectively; the ( - ) before an entry indicates that the angle is formed with the corresponding negative axis. Angle of incidence of the arriving particles, or, and deposit porosity, e, are also provided. Note further that these results are not averages, but were obtained from particular deposit realizations. The first two sets of data correspond to a high and a lower porosity deposit, both formed by particles arriving normal to the target, the next two correspond to I!+ = 60 and the final two to 0, = 45 degrees. As already stipulated the eigenvectors associated with F,, are aligned with the z-axis (respective y = 1.4 and 1.1); the eigenvectors associated with D,, give y = 10.7 and 16.4. In a rather qualitative sense [we have not been able to obtain quantitative statistics on the variation between different realizations,

440

and DANIEL E.ROSNER

MENELAOSTASSOFOULOS

Table 4. Principal directions of the fabric and di&sivity

tensors for typical deposits D

F a

Y

u

B

Y

89.3 1.4

78.2 ( - )12.1 63.4

( - )15.3 ( - )78.8 79.7

80.4 85.4 10.7

10.2 79.8 ( - )89.8

89.9 89.0 1.1

34.7 ( - )56.3 82.6

55.4 38.4 ( - )75.4

87.7 73.7 16.4

53.9 82.3 ( - )37.2

77.9 12.1 ( - )89.5

38.7 80.7 52.9

37.6 ( - )88.1 ( - )52.5

( - )809 17.6 ( - )75.1

53.9 72.5 41.4

0.696

89.6 46.1 ( - )43.9

(-)1.4 88.6 89.3

88.6 43.9 46.1

25.3 80.1 ( - )67.0

( - )80.0 10.2 71.7

67.1 87.8 23.1

4s

0.869

62.5 52.9 ( - )49.4

( - )38.3 51.7 89.0

65.6 59.9 40.6

38.7 78.8 ( - )53.6

( - )75.0 15.0 ( - )88.7

55.3 80.1 36.4

45

0.611

71.2 38.2 ( - )58.1

( - )22.5 67.5 89.6

78.0 60.9 31.9

24.3 77.6 ( - )69.5

74.2 ( - )17.4 83.0

72.1 78.1 21.8

ii

I3

E

11 22 33

0

0.858

11 22 33

0

0.558

( - )79.8 80.2 ( - )89.9

11 22 33

60

0.079

11 22 33

60

11 22 33 11 22 33

B

59.5 ( - )30.6 89.4

88.9

All angles reported in degrees. The ( - ) before an entry indicates that the reported angle is formed with

the corresponding negative axis.

due to the relatively small number of configurations (3 to 5) considered] we can say that for the deposits analyzed in this work we have always found that the principal directions of F and D although not coinciding were still strongly correlated typically forming angles between themselves of 5 to 20”. The dynamic nature of the solid-molecule interactions in the Knudsen limit is further demonstrated by the “rectification” effects we observed in this regime. We have explored such effects using the test molecule transmission technique, which can readily account for diffusion in the positive and negative direction. Specifically, we determined the transmission probability, JT for molecules that started with II, > 0, that is for molecules diffusing away from the target, and similarly for molecules diffusing towards the target, i.e. with u, -C 0 (see also Section 3.1). Plots of fraction of molecules transmitted, fx , vs inverse of Penetration distance, l/L,, are shown in Fig. 6, for a ballistic particles

and a denser deposit, both generated by arriving normal to the target. The ( + )

symbols correspond to molecules diffusing towards the target. The estimated fTLI values, are also shown in the figure legend. For the open deposits we found that rectification effects where typically very small (less than about 2%). On the contrary, for denser deposits the transmission probability in the positive z-direction was up to 25% higher than in the negative z-direction, as demonstrated in Fig. 6. The observed rectification effects, reminiscent of the Maxwell demon [see, e-g_ the discussion in Feynman et al. (1963)], would

mostly affect the short-time diffusion coefficient of vapor molecules through such deposits and will be the subject of future work from this laboratory. 6. CONCLUSIONS

AND

GENERALIZATIONS

In this study we determine tensor of computer-generated,

the effective diffusivity particulate,

anisotropic

Deposits of polydispersed, spherical particles were generated by several algorithmic models, which simulate different deposition conditions. Depending on the nature of the arriving particle trajectories (diffusive/deterministic), the angle of incidence and the subsequent motion after they establish contact with an already deposited particle both the deposit anisotropy and final solid fraction can be varied. Monte Carlo simulations of vapor molecules diffusing through cubic samples extracted from the generated deposits provided accurate estimates of the effective diffusivity tensor and associated tortuosities. The effective diffusivities were determined over the whole Knudsen number range, assuming self-diffusion in a constant pressure pure gas that consists of hardThe diffusing vapor sphere elastic molecules. molecules depart from vapor molecule collisions isotropically and are assumed to be reflected on the particle surface according to the laws of diffuse reflection. The free-paths between successive molecular collisions were sampled from an exponential distribution. with the probability of the free-path being greater than I given by exp ( - l/n), where 1 is the prevailing mean-free-path. Although this is admittedly deposits.

one

of

the simplest

molecular

models

of gaseous

Simulation

of vapor diffusion

in anisotropic

diffusion, it has been used quite extensively especially in analytical treatments of the transition regime (see e.g. Pollard and Present, 1948; Present, 1958; Strieder, 1971). For the near-continuum and transition regimes we have developed an efficient algorithm that has allowed us to determine Defr from discrete random walk simulations, thus avoiding the use of first passage time concepts which become questionable whenever the local Knudsen number is not much less than unity. The algorithm developed was efficient enough that we are able to perform the simulations reported 1. In the continuum limit here on a Sun Sparcstation we found that the effective diffusivities of the porous solids considered in this study could be closely correlated with porosity alone. Indeed, Archies’ law [eq. (17)j with the exponent m = 1.45 fitted all of our data quite accurately (Fig. 14). In the Knudsen limit on the other hand the details of the local microstructure became important and so an accurate prediction/correlation of the effective diffusivity requires higher order information of the type suggested by eq. (10). The measured orientation-averaged tortuosity factors, r, based on a Knudsen reference diffusivity defined with respect to a tube of diameter equal to the hydraulic diameter, varied with the local Knudsen number up to 25% indicating that z cannot be a priori assumed constant. Still, an interpolation based on the measured effective diffusivities in the two Knudsen limits, i.e. Scott and Dullien (1962), prediied the diffusivities in the transition regime very successfully. For the deposits analyzed here the diffusivity tensor becomes anisotropic only in the high Knudsen regime, that is when difIusing molecule/solid interactions dominate. In this case the principal axes of anisotropy of the diffusivity tensor and the principal axes of anisotropy of the solid as quantified here by the fabric tensor [cf. eq. (2)7 form angles between approximately 5 and 15”. We have not been able to obtain more quantitative results due to the relatively small number of realizations considered in this study. Finally, in the Knudsen limit we have found that rectification effects can be quite significant. Specifically, the transmission probability as determined by the test-molecule method could be up to 25% higher when the molecules diffused away from the target than when they diffused towards the target. This effect was attributed to the tree-like microstructures of the deposits. In conclusion, we should point out that the “continuum” tortuosities presented here can be readily used to obtain accurate estimates of the fluid permeability of these deposits using the procedure outlined by Schwartz et al. (1989). Similarly, the Knudsen tortuosities can be also thought of as estimates of the radiative transport properties of these deposits, at least in the limit of nonconducting particles that are very large compared to the incident radiation wavelength. In the near future we will use the technique described in this communication to determine the effective diffusivity tensor of more realistic deposits generated by a dynamic deposition model

particulate

deposits

441

(Konstandopoulos and Tassopoulos. 1990). According to this model the particles are assumed to be frictional, hard, elastoplastic, adhesive spheres and outcome of the binary colliiions between the incident and target particles is based on conservation of linear and angular momentum. Acknowledgements-This research was supported,in part, by the U.S. Department of Energy-Pittsburgh Energy Technology Center via Grant DE-FG22-90 PC 90099 and the I!?90 Industrial AlEliates of the Yale High Temperature Chemical Reaction Engineering (HTCRE) Laboratory: Shell Corp., SCM-Chemicals, DuPont, Union Carbide and General Electric. The authors would also like to thauk A.G. Konstandopoulos and Professors P. Garcia-Ybarra, J. Castillo, S. Torquato and J. A. O’Brien for their comments and helpful discussions throughout this work

NOTATION

diffusivity tensor scalar diffusivity mean pore (hydraulic) diameter reference particle diameter fabric tensor formation factor fraction of molecules transmitted radial distribution function unit tensor thermal conductivity Knudsen number thickness of porous slab gas free path gas molecular weight displacement vector exponent in Archie’s law (eq. 17) contact normal (eq. 2) and local Pe

Q R R* : S

V V <>

normal

(es. 10) Peclet number rotation matrix universal gas constant reference particle radius position vector total pore surface area displacement distance of travel two-point correlation function total pore surface area per unit volume gas temperature time random variable uniformly distributed between 0 and 1 total void volume velocity vector ensemble average

Greek letters porosity ; angle of incidence A gas mean-free-path tortuosity factor solid fraction

442

MENELAOS

TAS.WKKJLOS and DANIEL

Subscripts efkctive e ff B i ii

direction direction

Kn Br

Knudsen Brownian

XY

parallel

along

along

of growth of incidence

the ii-direction,

i -

1, 2, 3

to the xy-plane the z-direction

intrinsic. reference

:

REFERENCE

Abbasi M. H., Evans, J. W. and Abramson, I. S.. 1983, Diffusion of gases & porous solids: Monte-Carlo simulations in the Knudsen and ordinary diffusion regimes.

A.I.Ch.E. J. 29,617-624. Bamocky. G. and Davis, R. H., 1988, ‘Elastohydrodynamic collisi& and rebound of spheres: experim&tal ~erification. Phys. Fluids 31, 1324-1329. Bennett, C. H., 1972, ‘Serially deposited amorphous aggregates of hard spheres. J. appl. Phys. 43,2727-2734. Bemal, J. D., 1964, The st-&tu& of liquids. Proc. R. Sot. A Zso, 299-322. Berryman, I. G., 1983, Random close packing uf hard spheres and disks. Phys. Rev. A27, 1053-1061. Brown, G. W., 1956. Monte Carlo methods, in Modern Mathematics for the Engineer, (Edited by E. F. Beckenbach), McGraw-Hill, New York. Bruggeman, D. A. G., 1935, Bcrcchnung Versciedener Physikalisher Konstanten von Heterogenen Substanzen I. Dielektrizitltskonstanten und I..eitT”ahiglreiten der Mischkiirper aus Isotropen Substanzen. Ann. Physik. 24,

636-679. Burganos, V. N. and Sotirchos, S. V., 1988, Simulation of Knudsen diffusion in random networks of parallel pores. Chem. Engng Sci. 43,1685-1694. Burganos, V. N. and Sotirchos. S. V., 1989a, ‘Knudsen diffusion in parallel multidimensional or randomly oriented capillary structures.’ Chem. Engng Sci. 44, 2451-2462. Burganos, V. N. and Sotirchos, S. V., 1989b, Effective diffusivities in cylindrical capillary-s~herical~vity pore structures. Che&. Engng Sci 44, i659-2637. _Chandrasckhar, S.. 1943, Stochastic problems in _ physics and _ astronomy. Rev. mod.Phys. 15, l-89. Chapman, S. and Cowling, T. G., 1960, The Mathematical Theory of Non-Un(fo& Gases. Cambridge University Press, Cambridge. Cundalj, P. A. and Strack, 0. D. L., 1979, A discrete numerical model for granular assemblies. Gevtechnique 29,

47-65. Currie, J. A., 1960, Gaseous diffusion in porous media II: dry ucanular materials. &it. J. amd. Phvs. 11, 192-200. D&neke, B., 197i, The cap& of-aerosoi particles by surfaces. J. Colloid Interjkce Sci. 37, 342-353. Diemer, B. R., 1989. A simple model for sticking coefficient prediction, Paoer 132d. AIChE 1989 Annual Meeting.__ San Fran&co,-CA, November 5-10. Dirks, A. G. and Leamy, H. J., 1977, Columnar microstructure in vapor-deposited thin films. Thin Solid Fibns 47, 219-233. Dullien, F. A. L., 1979, Porous Media Fluid YlYun.sportand Pore Structure. Academic Press, New York. Evans, J. W.. Abbasi. M. H. and Satin, A., 1980, A Monte Carlo simulation of the diffusion of gases in wrous solids. J. them. Phys. 72, 2967-2973. a Family, F. and Vicsek, T., 1985, Scaling of the active zone in the eden process on percolation networks and the ballistic deposition model. J. phys. A: Math. Gen. 18, L75-L81. Feynman, R. P., Leighton, R. B. and Sands, M., 1963, The

E. ROSNER

Feynman Lectures on Physics. Vol. 1, pp. 46.1-46.9, Addison-Wesley, Reading, MA. Finney, J. L., 1970, Random packings and the structure of simple liquids I. The geometry of random close packing. Proc. R. Sot. Land. A319,479493. Fowler, R. H., 1955, Statistical Meckanlcs, 2nd edition, Cambridge University Press, Cambridge. Friedlander, S. K., 1977, Smoke, I)lrst and Haze-Fundamentals ofAerosol Behavior. Wiley, New York. Goldsmith, W., 1960, Impact, 1st Edition, Arnold, London. Guyon, E., Oger, L. and Plona, T. J., 1987, Transport properties in sintered porous media composed of two particle sizes. J. Phys. D: Appl. Phys. 20, 1637-1644. Hammersiey, J. M. and Handscomb, D. C., 1964. Monte Carlo M&hods. Methuen, London. Hashin, 2. and Shtrikman, S.. 1962, A variational approach to the theory of the effective magnetic permeability of multiphase materials. J. appl. Phys. 33, 1514-1517. Henderson, D., Brodsky, M.-H. &d Chaudhari, P., 1974, Simulation of structural anisotropy and void formation in amorphous thin films. Appl. Phyi: Lett. 25,641~643. Ho, F. G. and Strieder, W., 1979, Asymptotic expansion of the porous medium, effective diffusion coe&ient in the Knudsen number. J. chetn,Phys. 70, 5635-5639. Ho, F. G. and Strieder, W., 198b, Numerical evaluation of the porous medium effective diffmivity between the Knudsen and continuum limits. J. &em. Phys. 73, 62966300. Hoogscagen, J., 1955, Diffusion in porous catalysts and adsorbents. Ind. Engng Chem. 47,906-909. Houi, D. and Lenormand, R., 1984, in Kinetics of Aggregation and Gel&ion (Edited by F. Family and D. P. Landau), Elsevier, Amsterdam. Huizenga, D. G. and Smith, D. M., 1986, Knudsen diffusion in random assemblages of uniform spheres A.I.ChE;. J. 32, l-6. Jackson, R., 1977, Transport in Porous Catalysts&t Edition, Elsevier, Amsterdam. Jullien, R. and Meakin, P., 1987, Simple three-dimensional models for ballistic deposition with restructuring. Europhys. Lett. 4, 1385-1390. Kim, I. C. and Torquato, S., 1990, Determination of the conductivitv of heteroneneous media bv effective Brownian motion similation. J. appl:Phys. 68,3892-390% Kim, I. C. and Torouato, S.. 1991, Effective conductivity of s&pensions of h&d s&e&s by Brownian motion s&ulation. J. a&. Phys. 69, 2280-2289. Konstandop&los, A. G. and Tassopoulos, M.. 1990, Simulation of particle impaction: sticking- probability and . microstru&re evolu&on in the ‘frozen’ deposit limit. Paper 6E.6, AAAR 1990 Annual Meeting, Philadelphia, Pi, June 18-22; (Manuscript prepared for submiss& to J. Colfoid Interface Sci.). Mason, E. A. and Malinauskas,

A. P., 1983, Gas Transport in Porous Media: Tke Dusty-Gas Model, 1st Edition. Elsevier,

Amsterdam. Maxwell, J. C., 1873, Treatise on Electricity and Magnetism, 1st Edition, Clarendon Press, Oxford. McPhedran, R. C. and McKenzie, D. R., 1978, The conductivitv of lattices of spheres I. The simDle cubic lattice. Proc. R. &c. torui. A35$45-63. Meakin, P. 1983, ‘Diffusion controlled deposition on fibers and surfaces,’ Phys. Rev. A 27,2616-2623 Meakin, P. and Jullien, R., 1987, Restructuring effects in the rain model for random deposition. J. Physique 48, 1651-1662. Meakin, P., Ramanlal, P., Sander, L. M. and Ball, R. C.. 1986. Ballistic dew&ion on surfaces. Phvs. Rev. A34.

5091-5103.

-

Melkote, R. R. and Jensen, K. F., 1989, Gas diffusion in random-fiber substrates. A.Z.Ck.E. J. 35, 1942-1952. Melkote, R. R. and -Jensen, K. F., 1990, Simulation of transition and ordinary regime diffusion in fiber media,

Simulation

of vapor diffusion in anisotropic

paper 155g AIChE 1990 Annual Meeting, Chicago, IL. November 11-16. Nits&e, L. C. and Brenner, H., 1989, Eulerian kinematics of flow through spatially periodic models of porous media

Archs ration. Mech. Analysis 107,225-292. Olague, N. E., Smith, D. M. and Ciftcioglu, M., 1988, Knudsen diffusion in ordered sphere packings. A.I.Ch.E.

J. 34, 1907-1909. Onat, E. T. and Leckie,

F. A., 1988, Representation of mechanical behavior in the presence of changing internal structure. J. OPDI.AJech. ASME. 55, l-10. Pollard, W. G. &d Present, R. D., 1948, On gaseous selfdiffusion in long capillary tubes. Phys. Rev. 73, 762-774. Present, R. D., 1958, Kinetic Theory ofGases. McGraw-Hill, New York. Reyes, S. C. and Iglesia, E., 1991, Monte Carlo simulations of structural properties of packed beds. Chem. Engnu - _ Sci. 46, 1089-1099. Reyes, S. C., Iglesia, E. and Chiew, Y. C., 1989, Effective diffusion coe%zients of randomly-assembled sintered porous structures, paper 54j AIChE 1989 Annual Meeting, San Francisco, CA, November 5-10. Rosner, D. E., Konstandopoulos, A. G., Tassopoulos, M. and Mackowski, D. W., 1991, Deposition dynamics of combustion-generated particles: summary of recent studies of particle transport mechanisms, capture rates and resulting deposit microstructure/properties. Paper presented at 1991, Engineering Foundation Conference: Inorganic

Transformations

nnd Ash Deposition during Combustion,

Palm Coast FL (March 10-15) (in press). Rosner, D. E., 1988, Transport Processes in Chemically Reacting Flow Systems. Butterworth Publishers, Stoneham, MA. Third printing 1990. Rosner, D. E. and Nagarajan, R., 1987, Toward a mechanistic theory of net deposit growth from ash-laden flowing combustion gases self-regulated sticking of impacting particles and d
particulate

deposits

443

homogeneous continuum systems. Phys. Rev. B40, 9155-9161. Scott, D. S. and Dullien, F. A. L., 1962, Diffusion of ideal gases in capillaries and porous solids. A.1.Ch.E. J. 8,

113-117. Siegel, R. A. and Langer, R., 1986, A new Monte Carlo approach to diffusion in constricted porous geometries. J. Colloid Interjkce Sci., 109, 426440. Slatterv. J. C., 1972. Momentum. Enerav and Mass Transfer in . Con&&a, ist Edition. McGraw-&i, New York. Smith, D. M., 1986, Knudsen diffusion in constricted pores: Monte Carlo simulations. A.I.Ch.E. J. 32, 329-331 Sotirchos. S. V.. 1969. ‘Multicomwnent diffusion and convec&on in. capillary structu~s,’ A.1.Ch.E. J. 35, 19X%1961. Sotirchoq S. V. and Tomadakis, M. M., 1989, Modelling transport, reaction, and pore structure during densification of cellular or fibrous structures. MRS 1989 Fall Meeting, Boston, MA. Strider. W. C., i971. Gaseous self-diffusion through a porous medium, J. them. Phys. S4,4050-4053. Strieder. W. C. and Prager, S., 1968, Knudsen flow through a porous medium. Phys. Fluids 11,2544-2548. T&sopodos, M.. O’B&en, J. A. &d Rosner, D. E., 1989a, Simulation of microstructure/mechanism relationships in particle deposition. A.Z.Ch.E. J. 35, 967-980. Tassopoulos, M.. O’Brien, J. A. and Rosner, D. E., 1989b, Simulation of d-sits generated bv the canture of oartitles that follow dete&inistic tr;jectoriei paper. 6Oj AIChE 1989 Annual Meetinn. San Francisco. , CA. I November 5-10. Thornton, C. and Barnes, D. J., 1986, Computer simulated deformation of compact granular assemblies. Acta Mech. I.

64.45-61. Tok&aga, T. K., 1985, Porous media gas diffusivity from a free path distribution model. J. them. Phvs. 82.5298-5299. Torquato. S., 1987, Thermal conduct&y or disordered heterogeneous media from the microstructure. Rev. Chem. Engng 4, 151-204. Torquato, S. and Kim, I. C.. 1989, Efficient simulation technique to compute effective properties of heterogeneous media. Appl. Phys. L&t. 55, 1847-1849. Torquato, S. and Sen, A. K., 1990, Conductivity tensor of anisotropic comwsite media from the microstructure. J. appl. Phjs. 67, 1‘145-1155. Turner, J. C. R.. 1976, Two phase conductivitv. The electrical conductance~ of l&id-fi-uidized beds of -spheres, Chem. Engng Sci. 31,487-492 Underwood, E. E., 1970, Quantitative Stereology. AddisonWesley, Reading, MA. Fractd Growth Pherwmerta. World Vicsek, T., 1989. Scientific, Teaneck, NJ. Visscher, W. M. and Bolster& 1972, Random packing of equal and unequal spheres in two and three dimensions. Nature 239,504-507. Void, M. J., 1963, Computer simulation of floe formation in a colloidal suspension. J. Coil. Sci. lS, 684-695. Wall, S. and John, W.. 1988, Impact adhesion theory applied to measurements of particle kinetic energy loss. J. Aerosol

Sci. 19, 789-792. Whitaker, S., 1967, Diffusion and dispxsion in porous media. A.Z.Ch.E. J. 13, 420-427. Witten, T. A. and Sander, L. M., 1981, Diffusion limited aggregation, a kinetic critical phenomenon. Phys. Rev. L&t. 47, 1400-1403. Wang, P., Koplik, J. and Tomanic, J. P.. 1984, Conductivity and permeability of rocks. Phys. Rev. B30, 6606-6614. Zheng, L. H. and Chiew, Y. C., 1989, Computer simulation of diffusion-controlled reactions in dispersions of spherical sinks. J. them. Phys. 90, 322-327.