Dynamical Simulation of the Flow of Suspensions: Wall-Bounded and

6312

Ind. Eng. Chem. Res. 2002, 41, 6312-6322

Dynamical Simulation of the Flow of Suspensions: Wall-Bounded and Pressure-Driven Channel Flow C. Pozrikidis* Department of Mechanical and Aerospace Engineering, University of California, San Diego, La Jolla, California 92093-0411

A boundary-element method is implemented to simulate the motion of two-dimensional particles with arbitrary shapes suspended in a viscous fluid under conditions of Stokes flow. Numerical discretization of an integral equation of the first kind for the particle surface traction results in a system of linear algebraic equations for the components of the traction over boundary elements distributed along the particle contours, as well as for the velocity of translation and the angular velocity of rotation of the particles about designated centers. The linear system is solved by the method of successive substitutions based on a physically motivated iterative method that involves the decomposition of the influence matrix into diagonal blocks consisting of physical particle clusters and then carrying out updates by matrix-vector multiplication using the inverses of the diagonal blocks. The iterations are found to converge as long as the minimum gap between the particles is larger than a threshold that depends on the particle shape and level of numerical discretization. To improve the numerical efficiency, the stiffness of the governing equations due to lubrication forces developing between intercepting particles is removed by preventing the particles from approaching one another to within less than a specified distance. Simulations of singly periodic suspensions of circular or elliptical particles in semi-infinite shear flow bounded by a plane wall and in pressure-driven flow in a channel bounded by two parallel plane walls are carried out for an extended period of time. The results confirm the occurrence of particle migration due to an effective hydrodynamic diffusivity and illustrate the dependence of the suspension effective viscosity and dynamics of the microstructure on the solid-phase areal fraction and particle aspect ratio. 1. Introduction Considerable progress has been made in the past 15 years in describing the rheological and transport properties and the dynamics of the microstructure of suspensions of solid particles by numerical simulation at vanishing and nonzero Reynolds numbers. The most popular simulation algorithm for Stokes flow, coined the method of Stokesian dynamics for particulate Stokes flow, is built on the idea of composite expansions. In this approach, hydrodynamic interactions between wellseparated particles are computed in terms of multipole expansions implemented by Faxen’s laws, and lubrication forces developing between neighboring particles are taken into account by local solutions for squeezing and lubrication flow.1-4 Simulation algorithms based on standard or novel implementations of the finite-element method,5-9 the lattice Boltzmann method,10,11 and the finite-difference method12 have been developed by other authors, but the Stokesian dynamics formulation has been favored in the majority of numerical investigations.13-15 In practice, suspended particles appear in a variety of shapes, ranging from platelets encountered in photographic emulsions to slender particles and fibers encountered in liquid crystals and composite materials. Although a generalization of the Stokesian dynamics method to nonspherical shapes is possible,16,17 analytical and computational complexities, combined with high computational costs, have prevented dynamical simulations. * E-mail: [email protected]. Internet: http://stokes.ucsd.edu/c_pozrikidis. Tel.: (858) 534-6530. Fax: (858) 534-7078.

An alternative problem formulation hinges on integral equations of the first or second kind for Stokes flow originating from the standard or generalized boundaryintegral representation.18-22 One significant advantage of the integral-equation method is its ability to account for the presence of boundaries readily by the use of appropriate Green’s functions. Integral equations of the first kind for the distribution of the traction over particle surfaces arise from the single-layer representation descending from the standard boundary-integral formulation. Dynamical simulations of multiparticle flows based on this approach have been discouraged by two concerns: nonuniqueness of the solution of the integral equations resulting from the presence of eigenfunctions defined over the individual particle surfaces and prohibitive computational costs required to compile and solve the system of linear equations arising from the boundary-element discretization. Integral equations of the second kind for the density of a double-layer potential descending from generalized double-layer representations are generally amenable to iterative solutions by the method of successive substitutions, but the extra cost required for evaluating and integrating higher-dimensional and strongly singular kernels is a practical disadvantage. Recently, it was shown that the two main difficulties of the single-layer formulation and associated boundaryelement methods stated in the previous paragraph can be circumvented in a way that preserves the simplicity and efficiency of the boundary-element implementation.21,22 The multiplicity of solutions can be eliminated by deflating the spectrum of the single-layer operator

10.1021/ie010878e CCC: $22.00 © 2002 American Chemical Society Published on Web 05/30/2002

Ind. Eng. Chem. Res., Vol. 41, No. 25, 2002 6313

or by projecting the linear system arising from the boundary-element discretization onto the orthogonal complement of the eigenvectors of the single-layer potential, and the solution can be computed efficiently using a physically motivated iterative method based on successive substitutions. Other features of the numerical method include algorithms for dynamical particle clustering that determine the splitting of the influence matrix upon which the iterations rely and algorithms for preventing the particles from approaching one another by less than a specified minimum distance with the objective of eliminating the stiffness of the equations governing the particle motion and thus allowing for affordable simulation times. Simulations of doubly periodic suspensions of particles with circular and elliptical shapes in simple shear flow have demonstrated that, properly implemented, the boundary-element method allows the study of two-dimensional systems at low and moderate areal fractions.22 In this paper, we focus our attention on semi-infinite shear flow and pressure-driven, Hagen-Poiseuille flow in a channel confined between two parallel plane walls. Previous authors simulated the flow of suspensions of spherical particles arranged in a monolayer using the method of Stokesian dynamics and found that the particles slowly migrate away from the walls and toward the channel centerline.23 Lateral migration results in a particle velocity profile that is blunted with respect to the parabolic profile of the unperturbed Poiseuille flow, in agreement with some, but not all, laboratory observations dating back to Karnis et al.,24 as reviewed in ref 25. Simulations of concentrated suspensions at particle volume fractions higher than 0.30 have shown that steady state is achieved after an equilibration time on the order of (a/〈u〉)(H/a)3, in agreement with theoretical predictions based on the concept of hydrodynamic diffusivity (e.g., ref 26); in this expression, a is the particle radius, H is the channel half-width, and 〈u〉 is the average suspension velocity. Subsequent simulations have illustrated the effect of gravitational settling in horizontal channel flow.27 More recently, it was demonstrated that particles in low-frequency, quasisteady, oscillatory Stokes flow can migrate toward the walls instead of the centerline for a certain range of mean displacements over a half-cycle of the oscillation, in agreement with laboratory observations.28,29 Our main objective in this paper is to illustrate the effect of the suspended-phase areal fraction and particle aspect ratio on the dynamics of the microstructure and effective viscosity of the suspension in wall-bounded flow. In section 2, we present the problem statement and mathematical formulation in the general framework of flow in a domain with arbitrary geometry, assuming that the corresponding Green’s function is available in analytical or readily computable form. In section 3, we discuss the method of solution and outline the numerical method. In section 4, we report simulations on the expansion of a particle bed in simple shear flow with the objective of demonstrating the occurrence of particle migration due to an effective hydrodynamic diffusion mediated by particle interceptions. In section 5, we consider pressure-driven channel flow in a channel bounded by two parallel walls and discuss the dynamics of the microstructure and the suspension effective viscosity. In section 6, we recapitulate and comment on theoretical modeling as an alternative to numerical simulation.

Figure 1. Schematic illustration of (a) an isolated swarm of particles in semi-infinite shear flow above a plane wall, (b) a periodic collection of particles in semi-infinite shear flow above a plane wall, (c) a periodic collection of particles in a channel bounded by two parallel plane walls, and (d) an isolated or doubly periodic collection of particles in unbounded simple shear flow.

2. Problem Statement and Mathematical Formulation Consider the flow of a suspension of two-dimensional rigid particles in one of the configurations illustrated in Figure 1. The Reynolds number of the flow around each particle is assumed to be so small that the motion of the fluid is governed by the linear equations of Stokes flow, including the Stokes equation and the continuity equation

∇p ) µ∇2u + Fg, ∇·u ) 0

(2.1)

where p is the pressure, µ is the fluid viscosity, u is the fluid velocity, F is the fluid density, and g is the acceleration of gravity (e.g., ref 20). Neglecting the inertia and buoyancy of the particles, we find that, in the absence of an externally induced force or torque other than gravity, the velocity of translation and the angular velocity of rotation of each particle about a designated center are determined by the requirement that the hydrodynamic force and torque exerted on the particle vanish. We begin setting up the mathematical formulation by decomposing the fluid velocity into an incident component u∞ and a disturbance component uD due to the presence of the particles. For shear flow past a nonperiodic or periodic collection of particles in a semi-infinite domain bounded by a plane wall located at y ) w, as illustrated in Figure 1a,b, the incident flow is given by

u∞x ) k(y - w) +

δ (y - w)2, u∞y ) 0 2µ

(2.2)

where k is the shear rate at the wall and δ is the pressure gradient; the associated pressure field is p∞ ) -δx. For shear- or pressure-driven flow of a nonperiodic or periodic collection of particles in a channel bounded by two parallel plane walls located at y ) (H, as illustrated in Figure 1c, the incident flow is given by

δ y+H + (H2 - y2), u∞y ) 0 2H 2µ (2.3)

u∞x ) U1 + (U2 - U1)

where U1 and U2 are the velocities of the lower and upper walls, respectively; H is the channel semi-width;

6314

Ind. Eng. Chem. Res., Vol. 41, No. 25, 2002

and δ is the pressure gradient. Again, the associated pressure field is p∞ ) -δx. For simple shear flow past a doubly periodic suspension, as illustrated in Figure 1d, the incident flow is given by

u∞x ) ky, u∞y ) 0

∫C Gij(x,x0) ni(x) dl(x) ) 0 q

(2.4)

where k is the constant shear rate; the associated pressure field is uniform. In all cases, the incident flow causes the kth particle to translate with linear velocity U(k) and to rotate about a designated center x(k) c around the z axis with angular velocity Ω(k). Using the standard boundary-integral formulation of Stokes flow (e.g., ref 20), and requiring no-slip and nopenetration boundary conditions at the particles surfaces, we find that the distribution of the surface traction, f ≡ σ‚n, where σ is the Newtonian stress tensor and n is the unit normal vector pointing into the fluid, satisfies the integral equation of the first kind

1

Now, mass conservation for the flow due to a point force situated at the point x0 requires the integral identity

for q ) 1, ..., Np. This identity suggests that any solution of the integral eq 2.5 over the contour of each particle can be enhanced with an arbitrary multiple of the unit normal vector. Thus, the integral equation has an infinite number of solutions parametrized by Np arbitrary constants. A symmetry property allows us to switch the arguments of the Green’s function, provided that we also switch the indices (e.g., ref 20), which is another way of saying that the single-layer operator is self-adjoint. Using this property, we recast the identity in eq 2.9 into the form

∫C Gij(x,x0) nj(x0) dl(x0) ) 0 q

Np

∑ ∫C Gij(x,x0) fi(x) dl(x) ) Rj(x0)

4πµ q)1

q

(2.5)

∫C Rj(x0) nj(x0) dl(x0) ) 0 q

Np is the number of particles, Cq is the contour of the qth particle, l is the arc length along the particle contour, the subscript m designates the mth component of the vector difference enclosed by the preceding parentheses, and the point x0 lies on the contour of the kth particle. The kernel of the single-layer potential, Gij, is the Green’s function of two-dimensional Stokes flow corresponding to the particular geometry of the domain of flow under consideration. In the case of a doubly periodic suspension undergoing simple shear flow, the Green’s function conforms to the instantaneous periodicity of the flow. The x and y components of the force and the z component of the torque exerted on the kth particle are given by

F(k) i )

∫

∫

f dl, T(k) z ) zjm C i k

f (x - x(k) c )m dl (2.7) C j k

for i ) x, y. In the case of the inertialess, neutraly buoyant, force-free, and torque-free particles presently considered, the distribution of the traction must be such that the left-hand sides of eqs 2.7 vanish for all k. Equations 2.5 and 2.7 provide us with a system of integral/algebraic equations for determining the particle surface traction f and the linear and angular particle velocities U(k) and Ω(k). Once the velocities are available, the positions of the designated particle centers x(k) c and the particle inclination angles θ(k) can be updated using a standard method for integrating the ordinary differential equations

dx(k) dθ(k) c ) U(k), ) Ω(k) dt dt

(2.8)

subject to a specified initial condition corresponding to an assigned initial particle distribution.

(2.10)

which suggests the solvability condition

where (k) (k) Rj(x0) ≡ u∞j (x0) - U(k) j - jpmΩp (x0 - xc )m (2.6)

(2.9)

(2.11)

where Rj is defined in eq 2.6. The condition in eq 2.11 is satisfied for any incompressible incident flow. To render the solution of eq 2.5 unique, we can remove the multiple eigenfunctions of the single-layer integral operator by adding to the left-hand side a deflating term, obtaining

1

Np

∑ ∫C Gij(x,x0) fi(x) dl(x) +

4πµ q)1

q

∫C ni(x) fi(x) dl(x) - bm] ) Rj(x0)

nj(x0) [

m

(2.12)

where the point x0 lies on the surface of the mth particle, and bm is an arbitrary constant. Multiplying eq 2.12 by nj(x0), integrating the product over the surface of the mth particle, and using the identity in eq 2.10, we find

∫C ni(x) fi(x) dl(x) ) bm m

(2.13)

which can be substituted into eq 2.12 to produce eq 2.5. A more general method of regularization applicable to the discretized form of eq 2.5 will be described in section 3. 3. Numerical Method In the numerical method, each particle contour is discretized into boundary elements with straight, circular, or elliptical shapes, and the components of the traction are approximated with constant functions over the elements. Requiring that the integral eq 2.5 be satisfied at the midpoint of each element and enforcing the integral constraints in eqs 2.7, we obtain a system of linear equations of the form A‚w ) b, where A is the “master influence matrix”; the right-hand side b contains the x and y components of Rj defined in eq 2.6, followed by three zero entries corresponding to the lefthand sides of eqs 2.7. The majority of the entries of the master influence matrix are defined in terms of integrals of the Green’s function over the boundary ele-


ments. A smaller number of entries are defined in terms of the simple kernels of the integral constraints in eqs 2.7. The integrals over the elements were evaluated by the six-point Gauss-Legendre quadrature. Over the singular elements, the logarithmic singularity of the Green’s function was subtracted and then integrated analytically to improve the accuracy. The vector of unknowns w contains the x and y components of the traction at the boundary elements around each particle contour, followed by the x and y components of the particle velocity of translation and the scalar angular velocity of rotation about the z axis. 3.1. Uniqueness of Solution. The discrete counterpart of the integral identity in eq 2.9 is A(q)‚e(q) ) 0, where e(q) is the discrete right eigenvector for the qth particle, containing the x and y components of the normal vector at the midpoints of the boundary elements around the particle contour. The matrix A(q) is the diagonal block of the master influence matrix corresponding to the qth particle, reduced by the three last columns and three last rows associated with the particle linear and angular velocities and with the integral constraints of eqs 2.7. All entries of the matrix A(q) are defined in terms of integrals of the single-layer potential over the boundary elements of the qth particle. When the particle contour is polygonal, in which case the normal vector is constant over each element, the discrete identity is satisfied up to numerical error due to numerical integration. More generally, the discrete identity is satisfied up to the boundary-element discretization error. In either case, the master linear system is precisely or nearly singular, and normalization or preconditioning is required for a robust numerical method. 3.2. Singular Preconditioning and Regularization. The discrete version of the solvability condition in eq 2.11 is v(q)‚R(q) ) 0, where v(q) is the discrete left eigenvector of the influence matrix for the qth particle, containing the x and y components of the normal vector at the midpoints of the boundary elements around the particle contour multiplied by the element arc length. In the subsequent discussion, we shall tacitly assume that v(q) has been normalized so that its Euclidean norm is equal to unity. The elements of the vector R(q) contain the x and y components of the right-hand side of eq 2.6 evaluated at the collocation points. If the particle contour is polygonal, in which case the normal vector is constant over each element, and if the vector R(q) is also constant over each element, the left eigenvector v(q) satisfies the discrete solvability condition to machine precision. The discrete counterpart of the integral identity in eq 2.10 integrated with respect to x over the boundary elements is v(q)‚A(q) ) 0. If the particle contour is polygonal, this identity is satisfied up to the numerical error incurred by numerical integration. More generally, this identity is satisfied up to the error incurred by the particle contour discretization. One way to regularize the master linear system is by applying “singular preconditioning”, implemented by premultiplying both sides with a preconditioning matrix P defined as T

T

T

P ≡ (I - V(1)V(1) )·(I - V(2)V(2) )· ... ·(I - V(Np)V(Np) ) (3.1) where V(q) is the extended discrete left eigenvector for the qth particle and the superscript T denotes the

transpose. All elements of this vector are equal to zero, except for the elements residing in the block corresponding to the qth particle, which are equal to the corresponding elements of the normalized discrete left eigenvector v(q). The equations of the master linear system corresponding to the integral constraints in eqs 2.7 remain unchanged. After singular preconditioning, the master linear system is transformed into A ˆ ‚w ) b ˆ, where A ˆ ≡ P‚A and b ˆ ≡ P‚b. Note that V(q) is an exact left eigenvector of the singular matrix A ˆ corresponding to the vanishing eigenvalue and the solvability condition ˆ ) 0 is satisfied to machine accuracy. Thus, V(q)‚b singular preconditioning guarantees the exact satisfaction of both the discretized form of the integral identity in eq 2.9 and the solvability condition in eq 2.11, while introducing the mildest possible perturbation to the discretized governing equations. The master linear system can now be regularized by discarding one linear equation corresponding to an arbitrary collocation point over each particle and replacing it with the discrete form of the integral constraint in eq 2.13. Regularization removes the Np degrees of freedom in the solution of the integral equation and allows us to solve a perfectly well-posed problem. 3.3. Iterative Solution of the Master Linear System. The regularized master linear system was solved using a physically motivated iterative method based on the equation

D·w(k+1) ) N·w(k) + b

(3.2)

where the matrix D consists of diagonal blocks corresponding to particle clusters, as will be discussed later in this section, and N ≡ A ˆ - D. The iterations are carried out by computing the inverses of the diagonal blocks comprising the matrix D using Gauss elimination and then explicitly solving for w(k+1) by carrying out matrix-vector multiplications. When the particles are well separated and each diagonal block contains one particle, only a few iterations are necessary to obtain the solution with accuracy that is comparable to the round-off error. When, however, particle clusters develop, achieving rapid convergence, or even convergence at all, requires that the diagonal blocks contain particle clusters. Solving the master linear system by an iterative method is crucial for the practical feasibility of the simulations. For example, in simulations with 36 particles and 16 boundary elements along each particle contour, corresponding to 35 scalar unknowns per particle, the size of the master linear system is 1260 × 1260. Solving this system by Gauss elimination requires on the order of 109 operations. By comparison, with the iterative method, computing the inverse of 36 diagonal blocks requires 104 operations, and carrying out the iterations requires an additional 12602 × m ) 1.6 × 106 × m operations, where m is the number of necessary iterations, which is on the order of 10. Thus, issues of numerical stability aside, the iterative method allows computational savings by a factor of 102. 3.4. Particle Clustering. As two particles approach one another under the influence of the incident flow, the convergence of the iterations based on eq 3.2 degrades. When the particles nearly touch, the iterations do not converge, and the simulation fails. To ensure rapid convergence, the particles are arranged into dynamically defined clusters, and the iterations are carried out with respect to diagonal blocks of the master

6316


linear system corresponding to the clusters. In the numerical implementation, the particles are tabulated into a temporary cluster list immediately before the execution of each time step, so that successive ordered groups in the list consisting of a recorded number of particles form a cluster. The original particle labels are registered into an evolving label list, and the particle center coordinates and velocities are unscrambled to be recorded after the completion of each time step. 3.5. Cluster Definition. Cluster definition is done using an algorithm that is similar to the bubble-sort algorithm for arranging a group of numbers in ascending order, as follows. First, a symmetric pairwise association matrix L is defined such that Lij ) 1 if particles labeled i and j, and their periodic images when appropriate, are sufficiently close to form a cluster, and Lij ) 0 otherwise; the diagonal elements of L are set equal to zero. The computation of Lij will be discussed in section 3.6. Second, the cluster list is initialized to the permanent particle labels. Third, the first particle is compared with all other particles by inspecting L1,j, j ) 2, ..., Np, with the objective of identifying whether the first and the jth particle belong to the same cluster. If the kth particle is found to form a cluster with the first particle, then the labels of the second and kth particle are switched in the cluster list, and the corresponding rows and columns of the association matrix L are transposed. Fourth, the first and second new particles are compared with all subsequent particles by inspecting Li,j, where i ) 1, 2 and j ) 2, ..., Np. If the kth particle is found to form a cluster with the first or second particle, then the labels of the third and kth particle are switched in the label list, and the corresponding rows and columns of the association matrix L are transposed. The process then continues until no additional particles can be found that belong to the first cluster. When this has been achieved, the algorithm is repeated for the remaining particles that do not belong to the first cluster. The cluster-definition algorithm is similar, but not identical, to cluster-labeling algorithms used by previous authors.30-32 The latter is implemented by moving along a tree that emanates from a chosen particle, marking along the way all particles that belong to a cluster. Marked particles are overlooked in the next pass, and the interrogation continues until all particles have been assigned a label. Test computations have confirmed that the cluster-definition algorithm produces results that are identical to those generated by the cluster-labeling algorithm. The implementation of the cluster-definition algorithm, however, is simpler in the absence of selfcalling subroutines, a nonstandard feature of Fortran. 3.6. Computation of the Pairwise Association Matrix. To compute the pairwise association matrix, the contours of a pair of particles in the xy plane are described by the equations f1(x,y) ) 0 and f2(x,y) ) 0, where i ) 1, 2, such that fi(x,y) < 0 or fi(x,y) > 0, respectively, in the interior or exterior of the ith particle. To identify the point of minimum separation on the second particle, we compute the minimum of the function f1(x,y) with respect to the natural parameter describing the contour of the second particle, where the point (x, y) is located on the second particle. First, a rough minimum is computed by direct search, and then the minimum is refined using Newton’s method. If the value of the function f1 at the location of the minimum, denoted by xmin 2 , is negative, then the two particles

overlap. A similar method is used to compute the point of minimum separation on the first particle, denoted by min - xmin defines the vector of xmin 1 . The distance x1 2 minimum separation. To determine whether two particles belong to a cluster, the particle contours are expanded by the factor cl while the designated particle centers are held fixed, and testing for overlap is done on the expanded contours and their periodic images when appropriate. If overlap is detected, the corresponding entry of the association matrix Lij is set equal to unity. 3.7. Removal of Stiffness. In the course of the motion, two particles might approach one another and remain close for a substantial length of time on the order of several local inverse shear rates, forming long-lived doublets, triplets, chains, and clusters. Physically, the spontaneous development of a resilient microstructure is due to lubrication forces developing between the particle surfaces in the absence of surface mobility. The system of equations governing the motion of the particles becomes stiff during these interceptions, and an unaffordably small time step is required to prevent artificial particle overlap. Such difficulties are not encountered in the case of suspensions of liquid drops at moderate viscosity ratios and in the absence of surfactants where the interfaces are able to slide over each another propelled by the ambient shear flow. In previous Stokesian dynamics simulations of suspensions, the particles were allowed to overlap, but this did not seem to affect the performance of the numerical method seriously and was tolerated as an inevitable but potentially serious side effect of the numerical method.14 Glowinski et al.8 implemented a finite-element method in which circular particles were prevented from colliding by an artificial repelling force activated at close range. Algorithms for detecting contact and overlap of ellipsoidal particles and for tracking the minimum distance between two-dimensional particles with arbitrary shapes have been developed by other authors.33,34 In the present simulations, particle overlap causes the iterative method to diverge, thereby halting the simulation. To prevent particle overlap without reducing the size of the time step during the randomly occurring periods of cluster formation, we implemented two algorithms: one based on elastic collisions of expanded particle contours and the second based on geometrical exclusion. In the first algorithm, artificial elastic forces and accompanying torques are exerted on each member of a particle pair when the minimum distance between the particles is less than a preset threshold. To compute the elastic forces, the contours of two test particles are expanded by the factor els, and testing for overlap is done on the expanded contours as discussed in section 3.6. If the expanded contours or their periodic images are found to overlap, then a neutral force doublet is applied at the midpoint of the vector of minimum separation defined previously in this section. One force and the corresponding torque with respect to the particle center are applied to the first particle, and an equal and opposite force and accompanying torque are applied to the second particle. Both forces are parallel to the vector of minimum separation, and their magnitude is equal to the minimum separation multiplied by a modulus of elasticity. The force exerted on each particle is applied at the midpoint of the minimum separation vector and is directed toward the particle surface. Superficially, each force could be translated to


the respective point of minimum separation to give a separated point force doublet. If this is done, however, the net torque exerted on the particle pair will not necessarily vanish, and the particles will no longer be force-free and torque-free, as required. It is thus imperative that the force doublet be neutral, that is, that the constituent forces be applied at the same location. In the second algorithm, the particles are prevented from approaching one another to within less than a specified distance by geometrical exclusion. In the numerical implementation, the particle contours are expanded by the factor sep, and testing for overlap is done on the expanded surfaces. If the expanded particles or their periodic images are found to overlap, then the particle centers are displaced either symmetrically or unsymmetrically along the vector of minimum separation, so that the expanded contours become nearly tangential. Symmetric displacement presents the risk of infinite bouncing of particles that belong to one cluster or to neighboring clusters. To prevent this occurrence, unsymmetrical displacement, in which the position of the first particle is held fixed and the neighboring particles are displaced, was used in the majority of the simulations. Computational experimentation showed that the method of geometrical exclusion performs better than the method of elastic collisions, especially for noncircular shapes. The rate of convergence of the iterations is a strong function of the value of the expansion factor sep determining the minimum particle separation; the larger the value of sep, the higher the rate of convergence. When the time step is sufficiently small, geometrical exclusion produces visually smooth particle trajectories, while allowing the iterations to converge under a wider range of conditions. Accordingly, all simulations discussed in sections 4 and 5 were carried out using this method. 4. Expansion of a Particle Bed Next to a Wall Before proceeding to a discussion of pressure-driven channel flow, it is instructive to illustrate the concept of hydrodynamic diffusivity and its dependence on the particle aspect ratio by considering the expansion of a particle bed next to a plane wall under the action of an overpassing semi-infinite shear flow, as illustrated in Figure 1b. In the numerical simulations, 36 particles were distributed randomly within each periodic cell of a particle bed whose thickness was equal to the period L. The initial particle positions and orientations were produced by a random-number generator, with care taken to avoid interparticle and particle-wall overlap. The contour of each particle was discretized into 24 circular or elliptical elements, and eqs 2.8 describing the particle motion were integrated in time using the second-order Runge-Kutta method with a fixed time step of k∆t ) 0.05, where k is the shear rate of the overpassing shear flow. The Stokes flow Green’s function for semi-infinite flow is available in closed form and was computed by direct evaluation.20 A typical value of the numerical parameter sep determining the particle contour expansion for geometrical exclusion discussed in section 3.7 was 1.05 or 1.075. Particle cluster identification was done on particle contours expanded by a factor of 1.20 or 1.25. Each time step required a few minutes of CPU time on a 600-MHz Intel processor running Linux. The total simulation time, spanning

several thousand time steps for each case, was on the order of 1 week. Figure 2 shows stages in the expansion of a bed of circular particles of radius a ) 0.042 05D. The label above each frame indicates the corresponding dimensionless time kt since the beginning of the simulation. The initial areal fraction of the particulate phase inside the bed is φ0 ≡ 36πa2/(DL) ) 0.20, where D is the initial bed thickness. The results show that the initially homogeneous particle distribution develops a structure that is dominated by particle pairs, triplets, and clusters. The thickness of the bed slowly increases in time as a result of net particle displacements after random interceptions in a processes that can be loosely described as hydrodynamic diffusion. Figure 3 shows corresponding stages in the evolution of a bed of elliptical particles with an axis ratio of 3:1 for the same initial areal fraction φ0 ) 0.20. A comparison with Figure 2 clearly shows that the bed of elliptical particles expands more slowly than the bed of circular particles as a result of smaller net displacements after interception. Migration away from the walls, however, appears to be somewhat more pronounced in the case of the elliptical particles. Particles with very high aspect ratios rotate to align themselves with the streamlines of the simple shear flow, and are then convected with the local velocity field without introducing a perturbation flow that is responsible for bed expansion and particle migration away from the wall. 5. Channel Flow Consider now the pressure-driven flow of a periodic suspension, as depicted in Figure 1c. Each fundamental flow cell contains Np identical particles, each of area Ap, repeated in the x direction with period L. The areal fraction of the solid phase is φ ) NpAp/(2HL), where H is the channel half-width. We assume that the pressure drop across each period of the channel is held constant, independent of time, and that the flow rate changes in response to the unsteady motion of the suspension. The consideration of flow subject to a constant flow rate and varying pressure gradient, studied by previous authors,23,27 requires a simple substitution for the Green’s function. 5.1. Computation of the Green’s Function. Significant improvements in efficiency arise when the Green’s function for channel flow is computed by trilinear interpolation with respect to y0, y, and x - x0, where (x0, y0) are the coordinates of one point force in the periodic array and (x, y) are the coordinates of the evaluation point. Interpolation is done for a desingularized part defined by

G′ij ≡ G2D-1P-WW (x,y,x0,y0) - G2D-1P-W (x,y,x0,y0;+h) ij ij 2D-1P-W FS (x,y,x0,y0;-h) + Gij (x,y,x0,y0) (5.1) Gij where G2D-1P-W (x,y,x0,y0;w) is the singly periodic ij Green’s function of Stokes flow in a semi-infinite domain bounded by a plane wall located at y ) w and GFS ij (x,y,x0,y0) is the free-space Green’s function (Stokeslet). The smooth behavior of the components of the desingularized Green’s function G′ij near the point force and in the vicinity of the walls considerably improves the accuracy of the interpolation. The comple-

6318


Figure 2. Stages in the expansion of a bed of circular particles next to a plane wall under the action of an overpassing simple shear flow.

mentary components expressed by the last three terms on the right-hand side of eq 5.1 are available in closed form and were computed by direct evaluation. Test computations revealed that interpolation through a table of size 32 × 32 yields the Green’s function accurate to the fifth decimal place. 5.2. Effective Viscosity of the Suspension. In the absence of particles, the axial flow rate of the rectilinear Hagen-Poiseuille flow, denoted by Q0, is related to the negative of the pressure gradient δ ≡ -dP/dx by the

Hagen-Poiseuille law

Q0 )

2δH3 3µ

(5.2)

In the case of a singly periodic suspension, the axial flow rate is given by

Q ) Q0 -

H2

Np

∑∫C

2Lµ i)1

( )

i

1-

y2

H2

fx dl

(5.3)


Figure 3. Stages in the expansion of a bed of elliptical particles with an aspect ratio of 3:1 next to a plane wall under the action of an overpassing simple shear flow.

where δ ≡ -∆P/L is the pressure drop over each period and fx is the streamwise component of the traction.35-38 The relative effective viscosity of the suspension can now be defined in terms of the Hagen-Poiseuille and actual flow rate as

η≡

µeff µ

≡

Q0 Q

)1-

3

Np

∑∫C

4δLH i)1

( )

i

1-

y2

H2

fx dl (5.4)

To establish a point of reference, we note that the effective viscosity of an infinitely dilute suspension of circular particles in homogeneous simple shear flow is given by39

η ) 1 + 2φ

(5.5)

To make physical sense, infinite dilution should be interpreted as the absence of interparticle hydrodynamic interactions. 5.3. Numerical Simulations. Four main simulations were conducted for L/H ) 2 and Np ) 36 circular particles or elliptical particles with an aspect ratio 3 suspended within each period at areal fractions of φ ) 0.10 and 0.20. For the suspension of circular particles, φ ) 0.10 and 0.20 correspond, respectively, to a/H ) 0.084 10 and 0.059 47, where a is the particle radius. For the suspension of elliptical particles, φ ) 0.10 and 0.20 correspond, respectively, to a/H ) 0.048 55 and 0.034 34, where a is the major particle axis. In the main

simulations, the contour of each particle was discretized into 16 circular or elliptical elements, and eqs 2.8 describing the particle motion were integrated in time using the second-order Runge-Kutta method with a fixed time step of ∆t′ ) 0.05, where t′ ≡ tδH/µ is a dimensionless time. The other parameters of the simulation are the same as those discussed in section 4 for semi-infinite flow. Figures 4 shows the distribution of circular and elliptical particles initially and at an advanced stage of the simulation for φ ) 0.20. Viewing the results in animation reveals a multitude of elementary motions, including net migration away from the walls, pairwise and multiparticle interceptions, formation of transient clusters, and tumbling of the elliptical particles near the walls. A comparison of the top and bottom panels in Figure 4 suggests that the rate of migration of the elliptical particles is significantly more pronounced than that of the circular particles, in agreement with the results of section 4 for semi-infinite flow. This is also evident in Figure 5a, which shows the evolution of the RMS value of the y position of the particle centers. The RMS value of the centers of the circular particles does not show a significant change over the simulation period. If a significant change occurs, then it must take place over a period of time that is much longer than that required for a particle traveling with the mean fluid velocity to translate downstream by a distance equal to 100 times the clearance of the channel. In contrast,

6320


Figure 4. Snapshots of the distribution of circular (top) and elliptical (bottom) particles at the beginning of the simulation and at an advanced stage of the motion for an areal fraction of φ ) 0.20.

Figure 5. (a) Evolution of the RMS value of the y position of the centers of circular (solid line) and elliptical (dashed line) particles for φ ) 0.20. (b) Profile of the particle number density distribution (circles) and solid-phase velocity (squares); the solid lines are for the circular particles, and the dashed line are for the elliptical particles.

the RMS value of the centers of the elliptical particles shows a monotonic decline and appears to tend to an asymptotic limit at long times corresponding to statisti-

cal equilibration. A similar decline was observed in previous simulations of suspensions of spherical particles arranged in a monolayer at the higher areal fraction φ ) 0.40.23 The time scale for equilibration can be used to define an effective particle hydrodynamic diffusivity. The results presented in Figure 5a suggest that, for the same areal fraction, the effective diffusivity of the elliptical particles is significantly higher than that of the spherical particles. This conclusion, however, contradicts the results presented in section 4 on the expansion of a particle bed. The difference serves to underline the subtlety of particle motions in a nondilute suspension and warns against elementary phenomenological modeling. The circles in Figure 5b correspond to profiles of the particle number density distribution across the channel, normalized by H-1, defined such that the integral with respect to y across the clearance of the channel is equal to unity. The squares correspond to profiles of the solidphase velocity normalized by -∆pH2/(µL), where ∆p is the pressure drop over each period. These results were computed by averaging instantaneous values over the final simulation period of duration ∆t′ ) 30. The solid lines correspond to the circular particles, and the dashed lines correspond to the elliptical particles. The number density distribution for the elliptical particles is notably low in the vicinity of the walls as a result of the particle migration illustrated in Figure 4. The profiles of the solid-phase velocity are blunt with respect to the parabolic profile of single-phase Hagen-Poiseuille flow. It is especially interesting to note that the velocity of the elliptical particles is higher than the velocity of the circular particles; this suggests that a suspension of elliptical particles is transported more efficiently under the action of a constant pressure gradient. Figure 6a,b shows the evolution of the suspension relative effective viscosity defined in eq 5.4 for circular and elliptical particles. In both cases, the solid and dotted lines correspond, respectively, to areal fractions of φ ) 0.10 and 0.20. The solid line of limited duration in Figure 6a, hovering over the dotted line, represents the results of a simulation with twice the number of


in the bottom two frames of Figure 4. As the elliptical particles are convected by the flow, they spend a longer time oriented parallel, rather than normal to the walls, and this reduces the disturbance flow due to the solid phase. In the theoretical limit of infinite aspect ratio, the particles reduce to flat plates that tend to align with, and be convected parallel to, the streamlines of the unperturbed Hagen-Poiseuille flow without introducing a disturbance to the flow and thus affecting the effective viscosity of the suspension. 6. Discussion

Figure 6. Evolution of the suspension relative effective viscosity for a suspension of (a) circular, and (b) elliptical particles with an aspect ratio of 3:1. The solid and dotted lines correspond, respectively, to areal fractions of φ ) 0.10 and 0.20.

boundary elements, NE ) 32. This simulation took 8 days of CPU time on the aforementioned computer. A comparison with the dashed line suggests that the coarse discretization is sufficient for describing the main features of the motion. In all four cases illustrated in Figure 6, the relative effective viscosity undergoes an initial transition and then settles on a fluctuating path with a well-defined mean value. The mean value for φ ) 0.10 is in good agreement with the theoretical estimate (eq 5.5) for simple shear flow, which predicts η ) 1.2. For the higher areal fraction φ ) 0.20, the relative effective viscosity of the suspension of circular particles after equilibration, 〈η〉 ≈ 1.8, is substantially higher than the theoretical estimate obtained using eq 5.5. Previous simulations of monolayer suspensions of spherical particles predicted 〈η〉 ≈ 1.3 at φ ) 0.15 and 〈η〉 ≈ 1.6-1.9 at φ ) 0.30.23 These values are in good overall agreement with the present results. In contrast to the case of circular particles, the relative effective viscosity of the suspension of elliptical particles shows an overall increase from the initial value corresponding to the random distribution and exhibits a monotonic decline at longer times. The reason for this behavior is evident from the particle profiles depicted

We have discussed the implementation of a boundaryelement method for simulating the motion of suspensions of two-dimensional particles with arbitrary shapes and studied the expansion of a particle bed next to a wall under the influence of an overpassing shear flow, as well as the pressure-driven flow of a suspension in a channel bounded by two parallel walls. The results revealed that a bed of elliptical particles expands significantly more slowly than a suspension of circular particles, but that the circular particles appear to migrate away from the wall more quickly. The effective viscosity of a suspension in pressure-driven flow is reduced substantially as the particle aspect ratio becomes higher while the areal fraction is held constant. Previous authors have developed theoretical models to describe the dynamics of the mean velocity profile and particle concentration field in suspensions of spherical particles, as reviewed in ref 40. Theoretical models for suspensions of nonspherical particles have not yet been developed. In the “local model”, the flux of spherical particles is expressed as the sum of one contribution due to the spatial variation of the rate of particle interception and an opposing contribution due to the spatial variation of the effective viscosity associated with the generally nonuniform particle distribution.41,42 In the nonlocal model, the dynamics of the motion is described in terms of the intensity of the particle velocity fluctuations, termed the suspension temperature.23 In the case of tube flow, the local model predicts that the particle concentration profile has a cusp at the centerline, whereas the nonlocal model predicts a smooth profile with finite curvature. The sensitivity of the predicted flow behavior to the particular model underscores the importance of numerical simulation as a venue for probing the nature of the motion of inhomogeneous and wall-bounded suspension flows. Current efforts aim at extending the simulation algorithms described in this paper to three-dimensional flow. Acknowledgment This research was made possible by a grant provided by the National Science Foundation. Literature Cited (1) Brady, J. F.; Bossis, G. The rheology of suspensions in shear flow. J. Fluid Mech. 1985, 155, 105. (2) Brady, J. F.; Bossis, G. Stokesian Dynamics. Annu. Rev. Fluid Mech. 1988, 20, 111. (3) Sierou, A.; Brady, J. F. Accelerated Stokesian dynamics simulations. J. Fluid Mech. 2001, 448, 115. (4) Ichiki, K. Improvement of the Stokesian dynamics method for systems with a finite number of particles. J. Fluid Mech. 2002, 452, 231. (5) Hu, H. H.; Joseph, D. D.; Crochet, M. J. Direct simulation of fluid particle motions. Theor. Comput. Fluid Dyn. 1992, 3, 285.

6322


(6) Hu, H. H. Direct numerical simulation of flows of solidliquid mixtures. Int. J. Multiphase Flow 1996, 22, 335. (7) Maury, B.; Glowinski, R. Fluid-particle flow: A symmetric formulation. C. R. Acad. Sci. Paris Ser. I 1997, 324, 1079. (8) Glowinski, R.; Pan, T.-W.; Hesla, T. I.; Joseph, D. D. A distributed Lagrange multiplier/fictitious domain method for particulate flows. Int. J. Multiphase Flow 1999, 25, 755. (9) Patankar, N. A.; Singh, P.; Joseph, D. D.; Glowinski, R.; Pan, T. W. A new formulation of the distributed Lagrange multiplier/fictitious domain method for particulate flows. Int. J. Multiphase Flow 2000, 26, 1509. (10) Ladd, A. J. C. Sedimentation of homogeneous suspensions of non-Brownian spheres. Phys. Fluids 1996, 9, 491. (11) Aidun, C. K.; Lu, Y.; Ding, E.-J. Direct analysis of particulate suspensions with inertia using the discrete Boltzmann equation. J. Fluid Mech. 1998, 373, 287. (12) Ho¨fler K.; Schwarzer, S. Navier-Stokes simulation with constraint forces: Finite-difference method for particle-laden flows and complex geometries. Phys. Rev. E 2000, 61, 7148. (13) Revay, J. M.; Higdon, J. J. L. Numerical simulation of polydisperse sedimentation: Equal-sized spheres. J. Fluid Mech. 1992, 243, 15. (14) Dratler D. I.; Schowalter W. R. Dynamic simulation of suspensions of non-Brownian hard spheres. J. Fluid Mech. 1996, 325, 53-77. (15) Dratler D. I.; Schowalter W. R.; Hoffman, R. L. Dynamic simulation of shear thickening in concentrated colloidal suspensions. J. Fluid Mech. 1997, 353, 1. (16) Clayes, I. L.; Brady, J. F. Suspensions of prolate spheroids in Stokes flow. Part 1. Dynamics of a finite number of particles in an unbounded fluid. J. Fluid Mech. 1993, 251, 411. (17) Clayes, I. L.; Brady, J. F. Suspensions of prolate spheroids in Stokes flow. Part 2. Statistically homogeneous dispersions. J. Fluid Mech. 1993, 251, 443. (18) Power, H.; Miranda, G. Second-kind integral formulation of Stokes flow past a particle of arbitrary shape, SIAM J. Appl. Math. 1987, 47, 689. (19) Kim, S.; Karrila, S. Microhydrodynamics: Principles and Selected Applications; Butterworth-Heineman: Boston, 1991. (20) Pozrikidis, C. Boundary Integral and Singularity Methods for Linearized Viscous Flow; Cambridge University Press: Cambridge, U.K., 1992. (21) Pozrikidis, C. A spectral-element method for particulate Stokes flow. J. Comput. Phys. 1999, 156, 360. (22) Pozrikidis, C. Dynamical simulation of the flow of suspensions of particles with arbitrary shapes. Eng. Anal. Bound. Elem. 2001, 25, 19. (23) Nott, P. R.; Brady, J. F. Pressure-driven flow of suspensions: Simulation and theory. J. Fluid Mech. 1994, 275, 157. (24) Karnis, A.; Goldsmith, H. L.; Mason, S. G. The kinetics of flowing dispersions. I. Concentrated suspensions of rigid particles. J. Coll. Interface Science 1966, 22, 532. (25) . Phan-Thien, N.; Fang, Z. Entrance length and pulsating flows of a model concentrated suspension. J. Rheol. 1996, 40, 521.

(26) Davis, R. H. Hydrodynamic diffusion of suspended particles: A symposium. J. Fluid Mech. 1996, 310, 325. (27) Morris, J. F.; Brady, J. F. Pressure-driven flow of a suspension: Buoyancy effects. Int. J. Multiphase Flow 1998, 24, 105. (28) Morris, J. F. Anomalous migration in simulated oscillatory pressure-driven flow of a concentrated suspension. Phys. Fluids 2001, 13, 2457. (29) Butler, J. E.; Majors, P. D.; Bonnecaze, R. T. Observations of shear-induced particle migration for oscillatory flow of a suspension within a tube. Phys. Fluids 1999, 11, 2865. (30) Hoshen, J.; Kopelman, R. Percolation and cluster distribution. I. Cluster multiple labeling technique and critical concentration algorithm. Phys. Rev. B 1976, 14, 3438. (31) Gawlinski, E. T.; Stanley, H. E. Continuum percolation in two dimensions: Monte Carlo tests of scaling and universality for noninteracting disks. J. Phys. A: Math. Gen. 1981, 14, L291. (32) Chang, C.; Powell, R. L. Self-diffusion of bimodal suspensions of hydrodynamically interacting spherical particles in shearing flow. J. Fluid Mech. 1994, 281, 51. (33) Perram, J. W.; Rasmussen, J.; Praestgaard, E.; Lebowitz, J. L. Ellipsoid contact potential: Theory and relation to overlap potentials. Phys. Rev. E 1996, 54, 6565. (34) Maury, B. A many-body lubrication model. C. R. Acad. Sci. Paris, Ser. I 1997, 325, 1053. (35) Zhou, H.; Pozrikidis, C. The flow of suspensions in channels: Single files of drops. Phys. Fluids A 1993, 5, 311. (36) Zhou, H.; Pozrikidis, C. The flow of ordered and random suspensions of liquid drops in a channel. J. Fluid Mech. 1993, 255, 103. (37) Zhou, H.; Pozrikidis, C. Pressure-driven flow of suspensions of liquid drops. Phys. Fluids 1994, 6, 80. (38) Li, X.; Pozrikidis, C. Wall-bounded and channel flow of suspensions of liquid drops. Int. J. Multiphase Flow 2000, 26, 1247. (39) Belzons, M.; Blane, R.; Bouillor, J.-L.; Camoin, C. Viscosite’ d’une suspension dilue’e et bidimensionnelle de sph’eres. C. R. Acad. Sc. Paris Ser. II 1981, 292, 939. (40) Subia, S. R.; Ingber, M. S.; Mondy, L. A.; Altobelli, S. A.; Graham, A. L. Modelling of concentrated suspensions using a continuum constitutive equation. J. Fluid Mech 1998, 373, 193. (41) Leighton, D.; Acrivos, A. Measurement of shear-induced self-diffusion in concentrated suspensions of spheres. J. Fluid Mech. 1987, 177, 109. (42) Phillips, R. J.; Armstrong, R. C.; Brown, R. A.; Graham, A. L.; Abbott, J. R. A constitutive equation for concentrated suspensions that accounts for shear-induced particle migration. Phys. Fluids 1992, 4, 30.

Received for review October 23, 2001 Revised manuscript received April 5, 2002 Accepted April 8, 2002 IE010878E

Dynamical Simulation of the Flow of Suspensions: Wall-Bounded and

Recommend Documents