A numerical study of the motion of a neutrally
buoyant cylinder in two dimensional shear flow
Tsorng-Whay Pan111Corresponding author: e-mail: firstname.lastname@example.org, Shih-Lin Huang, Shih-Di Chen, Chin-Chou Chu,
Chien-Cheng Chang Department of Mathematics, University of Houston, Houston, Texas 77204, USA Institute of Applied Mechanics, National Taiwan University, Taipei 106, Taiwan, ROC Department of Mathematics and Taida Institute of Mathematical Sciences, National Taiwan University, Taipei 106, Taiwan, ROC
Abstract In this paper, we investigate the motion of a neutrally buoyant cylinder of circular or elliptic shape in two dimensional shear flow of a Newtonian fluid by direct numerical simulation. The numerical results are validated by comparisons with existing theoretical, experimental and numerical results, including a power law of the normalized angular speed versus the particle Reynolds number. The centerline between two walls is an expected equilibrium position of the cylinder mass center in shear flow. When placing the particle away from the centerline initially, it migrates toward another equilibrium position for higher Reynolds numbers due to the interplay between the slip velocity, the Magnus force, and the wall repulsion force. keywords: Shear flow; Neutrally buoyant particle; Equilibrium height; Fictitious domain/distributed Lagrange multiplier method; Finite element.
The problem of particle motions in shear flows is crucially important in many engineering fields such as the handling of a fluid-solid mixture in slurry, colloid, and fluidized bed. Segré and Silberberg [27, 28] experimentally studied the migration of dilute suspensions of neutrally buoyant spheres in a tube Poiseuille flow. The particles migrate away from the wall and the centerline to accumulate at about 0.6 of the tube radius from the centerline. The experiments of Segré and Silberberg [27, 28] have had a large influence on fluid mechanics studies of migration and lift of particles. Comprehensive reviews of experimental and theoretical works have been given by Brenner , Cox and Mason , Feuillebois  and Leal .
Among the theoretical studies of the neutrally buoyant particle migration in linear shear flow, Bretherton  found an expression for the lift force per unit length on a cylinder in an unbounded two-dimensional linear shear flow at small Reynolds number. Saffman’s lift force  on a sphere of radius in an unbounded linear shear flow with shear rate is where is the kinetic viscosity of the fluid, is the density of the fluid, is the particle Renolds number, and is the slip velocity of the sphere. In a bounded linear shear flow, Ho and Leal  examined the motion of a rigid sphere with inclusion of the inertia effects at small Reynolds numbers by a regular perturbation method. The sphere reaches a stable lateral equilibrium position which is the midway between the walls. Vasseur and Cox  also obtained the same stable lateral equilibrium position. Ho and Leal require that which is more restrictive than the one required by Vasseur and Cox where is the confined ratio, being the distance between two walls. Direct numerical simulations have been used for understanding particle motion in shear flows. Feng et al.  investigated the motion of neutrally buoyant and non-neutrally buoyant circular particle in plane shear and Poiseuille flows using a finite element method and obtained qualitative agreement with the results of perturbation theories and of experiments. The numerical results of a neutrally buoyant circular cylinder in a shear flow of have been discussed in details. The cylinder migrates back to the midway between two walls due to the wall repulsion at the small Reynolds number. They have suggested that that three factors, namely the wall repulsion due to a lubrication effect, the slip velocity, and the Magnus type of lift, are possible responsible for the lateral migration. Ding and Aidun  studied numerically the dynamics of a cylinder of circular or elliptic shape suspended in shear flow at various particle Reynolds number. They obtained the transient from being rotary to stationary as the particle Reynolds number is increased for an elliptic cylinder. For the cases of the circular cylinder, the effect of the two walls on the rotation speed has been studied. For the confined ratio , Ding and Aidun reported . Similar result,, has been observed experimentally by Zettner and Yoda .
In this paper, we first discuss the generalization of a distributed Lagrange multiplier/fictitious domain method (DLM/FD method) developed in  to non-spherical neutrally buoyant cylinders in two-dimensional shear flows and to the cases with zero angular velocity as a constraint. Similar DLM/FD methods for freely moving neutrally buoyant particle has been successfully applied, in [20, 23, 30], to simulate particulate flow in two and three dimensions with neutrally buoyant particles. But for the cases of a neutrally buoyant elliptic cylinder in two dimensional flows, the methodology has not been fully validated yet. We have validated the numerical method by comparing with the computational results in Ding and Aidun  for a cylinder of circular and elliptic shape and the experimental results in Zettner and Yoda  for a circular cylinder. On the wall effect on the angular velocity of the circular cylinder, we have obtained for the confined ratio which is in a good agreement with the results obtained by Ding and Aidun  and Zettner and Yoda . We have also studied the wall effect on the angular velocity of the elliptic cylinder which is more complicated due to the non-circular shape.
Concerning the equilibrium position of a neutrally buoyant circular cylinder in shear flow, recent studies by Cherukat, McLaughlin and Dandy  and Kurose and Komori  focus on lift and drag on a stationary sphere in unbounded linear shear flow. The equilibrium positions have not been studied in these works. Feng and Michaelides  have investigated the equilibrium positions of non-neutrally buoyant circular cylinders in two-dimensional shear flow. In their simulations, the density ratio between the solid and fluid is between 1.005 and 1.1. The equilibrium heights of their lightest circular cylinder (the density ratio of 1.005) are far below the centerline. We have obtained that the cylinder stays in the middle between two walls as expected when placing it there initially. But when the initial position of the mass center of a circular cylinder is away from the centerline, the equilibrium position depends on the particle Reynolds number and the confined ratio . For those placing away from the centerline initially, the circular cylinder migrates back to the centerline for where is the critical value. For , the equilibrium position is between the wall and the centerline. The critical particle Reynolds number is increasing when increasing the confined ratio. Concerning the Magnus lift effect on the equilibrium position, we have added a constraint, zero angular velocity, to the motion of a circular cylinder and obtained that the equilibrium position of the circular cylinder moving with zero angular velocity is lower than those of freely moving cylinder when both are away from the middle. These results show that the Magnus lift does play a role as expected. Also from the computed slip velocity of the circular cylinder, it shows that the circular cylinder lags the fluid, at least for , which means that there is a force pushing the cylinder toward the wall (see Fig. 1 for the setup of the boundary conditions). Hence the equilibrium position of the cylinder is up to the interplay between the slip velocity, the Magnus lift and the wall repulsion force. The content of this paper is as follows: In Section 2 we introduce a fictitious domain formulations of the model problem associated with the neutrally buoyant long particle cases; then in Section 3 we discuss the time and space discretization and in Section 4 we present and discuss the numerical results.
2. A fictitious domain formulation of the model problem
A fictitious domain formulation with distributed Lagrange multipliers for flow around freely moving particles and its associated computational methods have been developed and tested in, e.g., [12, 13, 24, 25]. For the cases of neutrally buoyant particles, similar methodologies have been developed in [20, 22, 23, 30]. But for the cases of a neutrally buoyant elliptic cylinder in two dimensional flows, the methodology has not been discussed and fully validated yet. In this paper, we first discuss the formulation for the case of a neutrally buoyant elliptic cylinder and then present the numerical results to validate the methodology.
Let be a rectangular region filled with a Newtonian viscous incompressible fluid (of density and dynamic viscosity ) and containing a freely moving neutrally buoyant rigid particle centered at of density . The flow is modeled by the Navier-Stokes equations and the motion of the particle is described by the Euler-Newton’s equations. We define
with , , and an inner product on which can be the standard inner product on . For simple shear flow, we have on the top wall and on the bottom wall. Then the fictitious domain formulation with distributed Lagrange multipliers for flow around a freely moving neutrally buoyant particle of the elliptic shape is as follows
where and denote velocity and pressure, respectively, is a Lagrange multiplier, is the translation velocity of the particle , is the angular velocity of , and is the angle between the horizontal direction and the long axis of the elliptic cylinder (see Fig.1). We suppose that the no-slip condition holds on . We also use, if necessary, the notation for the function .
The hydrodynamical forces and torque imposed on the rigid body by the fluid are built in (1)-(7) implicitly (see [12, 13] for details), thus we do not need to compute them explicitly in the simulation. Since in (1)-(7) the flow field is defined on the entire domain , it can be computed with a simple structured grid.
In (3), the rigid body motion in the region occupied by the particle is enforced via the Lagrange multiplier . To recover the translation velocity and the angular velocity , we solve the following equations
To investigate the effect of the Magnus type of lift on the lateral migration of the cylinder, we have considered the cases of the cylinder freely moving in shear flow with zero angular velocity. For this special consideration, we have the following modified formulation
with the modified multiplier space
The translation velocity is recovered via
3. Space approximation and time discretization Concerning the space approximation of problem (1)-(7) by a finite element method, we have used -- and finite elements for the velocity field and pressure, respectively (like in Bristeau et al. ). We approximate then , , and by the following finite dimensional spaces
A finite dimensional space approximating is defined as follows: let be a set of points covering ; the discrete multiplier space is defined by
where is the Dirac measure at . Then, we have the inner product defined by
and approximate by
Applying a first order operator splitting scheme à la Marchuk-Yanenko  to the equations (24)-(30) at each time step and the Euler backward method in time for some subproblems, we obtain (after dropping some of the subscripts ):
For , knowing , , , , and , compute and via the solution of
Then compute via the solution of
Next, compute via the solution of
Now predict the position and the orientation of the particle via:
Then set and . Now, compute , , , and via the solution of
and solve for and from
Finally, correct the position and the orientation of the particle via:
and set and . Finally, we set , , and . In above algorithm (31)-(43), we have , , is the region occupied by the particle centered at . Finally, and verify ; we have chosen and in the numerical simulations discussed later.
At each time step we have a sequence of subproblems (32), (33), (35) and (39). The degenerated quasi-Stokes problem (32) is solved by a preconditioned conjugate gradient method introduced in , in which discrete elliptic problems from the preconditioning are solved by a matrix-free fast solver from FISHPAK by Adams et al. in . The advection problem (33) for the velocity field is solved by a wave-like equation method as in [7, 21]. Problem (35) is a classical discrete elliptic problem which can be solved by the same matrix-free fast solver. To enforce the rigid body motion inside the region occupied by the particles, we have applied the conjugate gradient method discussed in, e.g., [20, 22, 23].
4. Results and discussions 4.1. The motion of an elliptic cylinder in linear shear flow
We have first considered the motion of a neutrally buoyant elliptical cylinder in linear shear flow studied by Ding and Aidun  for the validation purpose. Their results were obtained by the lattice Boltzmann equation. The domain of computation is , then the height is . The confined ratio is , and the aspect ratio is where is the length of the semi-major axis and is the length of the semi-minor axis. The density of the fluid is and the kinetic viscosity is determined by the specified value of the Reynolds number via . The shear rate is fixed at and (so the moving directions of the two walls are the same as those in Ding and Aidun ). The mesh size is and the time step is .
In Fig. 2, the computational results of the angle and angular velocity of the elliptic cylinder are in a good agreement with the Jeffery’s solution  for and the ones obtained by Ding and Aidun for and 0.25, respectively. Fig. 3 shows that at various particle Reynolds numbers from 0 to 7.5, our results of angular velocity match very well with the ones obtained by Ding and Aidun . The motion of the ellipse for is a periodic rotation with non-uniform angular velocity, while for , the ellipse does not rotate at all; instead it takes a stationary orientation in shear flow (see, e.g., Ding and Aidun  for further details). An example of the velocity field of around an elliptic cylinder with a stable orientation is shown in Fig. 6.
The minimum angular velocity decreases as the particle Reynolds number is increased with a nearly straight line relationship as shown in Fig. 6, where = 7.25 is the critical particle Reynolds number above which the rotational motion is stopped. Fig. 6 shows that the period of rotation increases rapidly as the particle Reynolds number is increased close to , and the results are in good agreements with those obtained by Ding and Aidun . The further studies of the motion of a cylinder of axisymmetric shape in simple shear flow will be reported in a forthcoming paper.
4.2. The rotation of a cylinder in a simple shear flow
In, e.g., Aidun and Ding , the angular velocity of a circular cylinder suspended in a simple shear flow has been studied via the direct numerical simulation. Experimental results have been also obtained in, e.g., Zettner and Yoda . In this section, we focus on the wall effect of the rotation speed of a cylinder freely suspended in a simple shear flow with various confined ratios. We have first considered the cases of a circular cylinder suspended in the middle between two walls initially. The domain of computation is with where is the circular cylinder radius. The density of the fluid is and the kinetic viscosity is . The confined ratios are 0.125, 0.25 and 0.5. We have varied the values of the shear rate to have different values of the particle Reynolds number . The initial position of the cylinder is at the midway between two walls. The mass center of the circular cylinder stays at the centerline between two walls in the simulations without giving any extra conditions for keeping it there. A typical velocity field (here, and =0.5) is shown in Fig. 6. At the zero Reynolds number, the rotation speed of a circular cylinder is from the Jeffery’s solution. At the small particle Reynolds numbers, the ratio is about to converge to 0.5 when decreasing the confined ratio from 0.5 to 0.125 as shown in Fig. 7. When 0.5, the normalized angular speed is found to be about 0.420 for . When 0.25, is about 0.482 for . Both are very close to those obtained in Ding and Aidun . When 0.125, is much closer to 0.5 for due to weaker effect from the walls.
When increasing the particle Reynolds numbers, the log-log plot of versus in Fig. 7 shows that for when 0.5. The one obtained in Ding and Aidun  is and the experimental results reported in Zettner and Yoda  is for the same confined ratio 0.5. For the smaller values of the confined ratios, 0.25 and 0.125, we have obtained and respectively. The results of 0.5 and 0.25 are in a good agreement with those obtained in Ding and Aidun .
For an elliptic cylinder suspended in linear shear flow, we have studied the effect of the aspect ratio . The mass center of an elliptic cylinder also stays at the middle between two walls if it is placed there initially. Since the elliptic cylinder has zero angular velocity when the particle Reynolds number is larger than the critical value , the behavior of the rotation of an elliptic cylinder is different from that of the circular cylinder. In Fig. 8, the minimal angular velocity decreases rapidly to about zero when the particle Reynolds number is closer but less than since the motion of the elliptic cylinder is about to transit into the one with a fixed orientation in linear shear flow. The minimal angular velocity decreases faster for the smaller aspect ratio shows that the slender shape is easier to reach a stable orientation in linear shear flow. But its maximal angular velocity shown in Fig. 8 have different behavior since the maximal angular velocity happens when the direction of the long axis is about perpendicular to the shear direction. For the small particle Reynolds numbers, the maximal and minimal values of the angular velocity are close to the Jeffery’s solution as in Fig. 8.
4.3. The equilibrium position
In Ho and Leal  and Vasseur and Cox , they concluded that the sphere reaches a stable lateral equilibrium position which is the midway between the walls for small particle Reynolds numbers. In Feng et al. , the cylinder migrates back to the midway between two walls at when placing it away from the middle between two walls. Feng et al. have suggested that that three factors, namely the wall repulsion due to a lubrication effect, the slip velocity, and the Magnus type of lift, are possible responsible for the lateral migration. In the previous section, we have obtained that the centerline is always (at least in the range we have studied in this paper) the equilibrium position for the cylinder of elliptic or circular shape when it is positioned there initially. When placing the mass center initially away from the middle between two walls, the mass center may not migrate back to the centerline.
For the cases of a circular cylinder of the confined ratio as in the previous section, the radius of the circular cylinder is , the density of the fluid is , and the kinetic viscosity is . To check the effect of the length of the channel, the initial height , the mesh size and time step on the equilibrium height of the mass center, we have considered different sets of parameter values as indicated in Fig. 9. The final equilibrium heights in Fig. 9 are almost the same for each value of considered here except those with zero angular velocity constraint. When 1 and 2, the mass center of the freely moving cylinder migrates back to the middle between two walls. But for higher particle Reynolds numbers, we have obtained an equilibrium height which is between the centerline and the wall. To find out the effect of the walls on the final equilibrium height, we have varied the distance between two walls. The initial height is always 0.1 unit below the centerline. The length of the computational domain is . The equilibrium height versus the particle Reynolds number is shown in Fig. 10. As expected, the wall repulsion force is stronger for the larger confined ratio. Hence the cylinder migrates back to the midway between two walls for for the range of considered here and the critical particle Reynolds number is increasing when increasing the confined ratio. We believe that the symmetric breaking shown in Figs. 9 and 10 is not a numerical artifact. To further validate this phenomenon, we have compared ours with the results obtained in Feng and Michaelides  where only the non-neutrally circular cylinders were considered. For the results presented in Fig. 11, the radius of circular cylinder is , the initial height is , the kinetic viscosity is , and the length and the height of the computation domain are and , respectively. The density of the non-neutrally circular cylinder is either 1.002 or 1.005. We have applied the DLM/FD method developed in, e.g., [12, 13] to obtain the computational results for the cases of the non-neutrally circular cylinder. Our results are in a good agreement with the results obtained in Feng and Michaelides . The transition of the equilibrium height from the one associated with a non-neutrally buoyant cylinder to that of a neutrally buoyant cylinder is consistent and support the existence of the symmetric breaking of the equilibrium height.
For the effect of the Magnus lift associated with the rotation, we have considered the cases of a circular cylinder of 0.25 moving with zero angular velocity. The equilibrium height of non rotating cylinder is lower than the one of the freely moving cylinder for as in Fig. 9. These results indicate that the Magnus lift does play a significant role as expected. The similar results concerning the Magnus lift have also been observed in Feng and Michaelides . When =1, the wall repulsion force is strong enough to push the cylinder back to the centerline even without the help from the Magnus lift.
We now focus on the effect of the slip velocity. Since the initial position of the cylinder is below the centerline, it is moving to the left in the lower region of the computational domain due to the given boundary condition (see Fig. 1). For getting the slip velocity, we first compute the fluid horizontal speed on the streamline through the mass center in front of the cylinder at the distance of the half of the computational domain width and then minus the horizontal speed of the disk to obtain the slip velocity. The negative slip velocity means that the cylinder speed to the left is slower than the fluid speed to the left since the both signs are negative. We then have that the cylinder lags the fluid. For =1, the slip velocity of the non rotating cylinder becomes negative after a short initial transition period as in Fig. 12. The one associated with free moving circular cylinder remains positive for, at least, the first 110 time units and then oscillates about zero. The cylinder with no rotation lags the fluid. But both cylinders migrate back to the middle between two walls due to that the wall repulsion force dominates the very weak slip velocity effect. We can see freely moving circular cylinder moves toward the centerline faster than the one with no rotation does in Fig. 12. When =5, the slip velocity becomes negative after the initial transition period as in Fig. 12. The cylinder for the both cases lags the fluid. The slip velocity of the one with no rotation is about two times larger than the other one. The rotating velocity of each case of free motion is about constant speed after a short initial transition period. Fig. 13 shows that the relative velocity field to the horizontal velocity of the cylinder mass center in which we can clearly see that the cylinder lags the fluid. Both cases have stronger slip velocity which creates force pushing the cylinder to the region with faster flow speed, which is the region next to the bottom wall. Without the extra help from the Magnus lift, the one without rotation is closer to the wall. Hence when the initial position is not at the middle of two walls, the balance between the effect of the slip velocity, the Magnus lift and the wall repulsion does play a role for determining the equilibrium position of the cylinder.
5. Conclusion We have investigated the motion of a neutrally buoyant cylinder of circular or elliptic shape in two dimensional shear flow of a Newtonian fluid by direct numerical simulation. The numerical results are validated by comparisons with existing theoretical, experimental, and numerical results, including a power law of the normalized angular speed versus the particle Reynolds number. The rapid slow down of the normalized minimal angular speed of an elliptic cylinder is totally different from the behavior of the circular cylinder since the motion of an elliptic cylinder can transit from rotating to a fixed orientation when the particle Reynolds number is increased beyond the critical value. The midway between two walls is an expected equilibrium position of the cylinder mass center in shear flow. But when placing the particle away from the centerline initially, it migrates toward another equilibrium position between the wall and the centerline for higher Reynolds numbers which is caused by the interplay between the slip velocity, the Magnus force, and the wall repulsion force. The further study of the motion of a cylinder of axisymmetric shape in simple shear flow will be reported in a forthcoming paper.
Acknowledgments. T.-W. Pan acknowledges the support by the US NSF under Grant No. DMS-0914788. S.-L. Huang, S.-D. Chen, C.-C. Chu, C.-C. Chang acknowledge the support by the National Science Council (Taiwan, ROC) under Contract Numbers, NSC97-2221-E-002-223-MY3, NSC99-2628-M-002-003 and NSC100-2221-E-002-152-MY3.
-  Adams J, Swarztrauber P, Sweet R. FISHPAK: A package of Fortran subprograms for the solution of separable elliptic partial differential equations (The National Center for Atmospheric Research, Boulder, CO, 1980).
-  Brenner H. Hydrodynamic resistance of particles at small Reynolds numbers. Adv. Chem. Engng. 1966;6:287-438.
-  Bretherton FP. The motion of rigid particles in a shear flow at low Reynolds number. J. Fluid Mech. 1962;14:284-304.
-  Bristeau MO, Glowinski R, Periaux J. Numerical methods for the Navier-Stokes equations. Applications to the simulation of compressible and incompressible viscous flow, Computer Physics Reports 1987;6:73-187.
-  Cherukata P, McLaughlin JB, Dandy DS. A computational study of the inertial lift on a sphere in a linear shear flow field. Int. J. Multiphase Flow 1999;25:15-33.
-  Cox RG, Mason SG. Suspended particles in fluid flow through tubes. Ann. Rev. Fluid Mech. 1971;3:291-316.
-  Dean EJ, Glowinski R. A wave equation approach to the numerical solution of the Navier-Stokes equations for incompressible viscous flow. C.R. Acad. Sci. Paris Série 1 1997;325:789-797.
-  Ding E, Aidun CK. The dynamics and scaling law for particles suspended in shear flow with inertia, J. Fluid Mech. 2000;423:317-344.
-  Feng, Z-G, Michaelides, EE. Equilibrium position for a particle in a horizontal shear flow. Int. J. Multiphase Flow 2003;29:943-957.
-  Feng J, Hu HH, Joseph DD. Direct simulation of initial value problems for the motion of solid bodies in a Newtonian fluid. Part 2: Couette and Poiseuille flows. J. Fluid Mech. 1994;277:271-301.
-  Feuillebois F. Some theoretical results for the motion of solid spherical particles in a viscous fluid, in Multiphase science and technology, edited by G.F. Hewitt, J.M. Delhaye and N. Zuber (Hemisphere Pub. Corp., New York, 1989), Vol. 4, 583-798.
-  Glowinski R, Pan T-W, Hesla T, Joseph DD. A distributed Lagrange multiplier/fictitious domain method for particulate flows. Int. J. Multiphase Flow 1999;25:755-794.
-  Glowinski R, Pan T-W, Hesla T, Joseph DD, Periaux J. A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: Application to particulate flow. J. Comput. Phys. 2001;169:363-427.
-  Glowinski R, Pan T-W, Periaux J. Distributed Lagrange multiplier methods for incompressible flow around moving rigid bodies. Comput. Methods Appl. Mech. Engrg. 1998;151:181-194.
-  Jeffery GB. The motion of ellipsoidal particles immersed in a viscous fluid, Proc. R. Soc. Lond. A 1922;102:161-179.
-  Ho BP, Leal LG. Inertial migration of rigid spheres in two-dimensional unidirectional flows. J. Fluid Mech. 1974;65:365-400.
-  Kurose R, Komori S. Drag and lift forces on a rotating sphere in a linear shear flow. J. Fluid Mech. 1999;384:183-206.
-  Leal LG. Particle motions in viscouS. Ann. Rev. Fluid Mech. 1980;12:435-476.
-  Marchuk GI. Splitting and Alternating Direction Methods, in Handbook of Numerical Analysis, edited by P.G. Ciarlet and J.L. Lions, (North-Holland, Amsterdam, 1990), Vol. I, 197-462.
-  Pan T-W, Chang C-C, Glowinski R. On the motion of a neutrally buoyant ellipsoid in a three-dimensional Poiseuille flow. Comput. Methods Appl. Mech. Engrg. 2008;197:2198-2209.
-  Pan T-W, Glowinski R. A projection/wave-like equation method for the numerical simulation of incompressible viscous fluid flow modeled by the Navier-Stokes equations. Computational Fluid Dynamics Journal 2000;9:28-42.
-  Pan T-W, Glowinski R. Direct simulation of the motion of neutrally buoyant circular cylinders in plane Poiseuille flow. J. Comput. Phys. 2002;181:260-279.
-  Pan T-W, Glowinski R. Direct simulation of the motion of neutrally buoyant balls in a three-dimensional Poiseuille flow. C. R. Mecanique, Acad. Sci. Paris 2005;333:884-895.
-  Pan T-W, Joseph DD, Glowinski R. Modeling Rayleigh-Taylor instability of a sedimenting suspension of several thousand circular particles in direct numerical simulation. J. Fluid Mech. 2001;434:23-37.
-  Pan T-W, Joseph DD, Bai D, Glowinski R, Sarin D. Fluidization of 1204 spheres: simulation and experiments. J. Fluid Mech. 2002;451:169-191.
-  Saffman G. The lift on a small sphere in a slow shear flow. J. Fluid Mech. 1965;22:385-400.
-  Segré G, Silberberg A. Radial particle displacements in Poiseuille flow of suspensions. Nature 1961;189:209-210.
-  Segré G, Silberberg A. Behavior of macroscopic rigid spheres in Poiseuille flow. Part I. J. Fluid Mech. 1962;14:115-135.
-  Vasseur P, Cox RG. The lateral migration of a spherical particle in two-dimensional shear flows. J . Fluid Mech. 1976;78:385-413.
-  Yang BH, Wang J, Joseph DD, Hu HH, Pan T-W, Glowinski R. Numerical study of particle migration in tube and plane Poiseuille flows. J. Fluid Mech. 2005;540:109-131.
-  Zettner CM, Yoda M. The circular cylinder in simple shear at moderate Reynolds numbers: An experimental study. Expts. Fluids 2001;30:346-353.