
Direct Simulation of the Motion of Particles in Flowing Liquids
Computational Methods
Two separate lines of code development have been pursued by the principal
investigators to date, a generalized ALE Galerkin finite element method and an
embedding method. Both approaches have been initiated by us, for quite
different kinds of applications. At present, one scheme does not fit all
applications. Perhaps ultimately, a ``best'' universal scheme for moving
particles may evolve, but it is not presently prudent to make a bet.
An important issue to be addressed in developing the computational schemes is
the modeling of particleparticle and particlewall collisions, with proper
account taken of the frictional forces between the particles (or between a
particle and the wall) at the moment of contact. This is an open problem.
In experiments, particles are observed to touch one another and the retaining
walls. In slow flows, the contact is usually gentle, so rebound and friction
are probably not important. In faster flows, however, collisions are much more
vigorous, so it is imperative to develop a satisfactory theoretical and
computational model for them. There is no solution to this fundamental problem
in the literature. The best current codes refine the mesh as the particles
approach, but this can only be carried so far. These codes do seem to reproduce
the local dynamics correctly, even though collisions are not well
represented.
A crucial computational issue to be addressed is the efficient solution of the
various algebraic systems which arise in the schemes. These systems can be
extremely large for 3D problems, and their solution can consume up to 95 of
the CPU time of the entire simulation. It is therefore imperative to use
efficient iterative solution methods, with matrixfree preconditioners, and to
implement them on parallel architectures.
Recent computational approaches to solidliquid flows, possibly inspired by
molecular dynamics, are cellular automata and the lattice Boltzmann method.
These models can handle huge numbers of particles. However, they replace the
equations of motion with computer rules and do not deal with stagnation and
separation points, wakes, turning couples, drafting, kissing and tumbling, etc.
The interesting results produced by these methods are not yet sufficiently
reliable to be used in engineering practice.
It is possible to simplify the flow description considerably by ignoring
viscous effects completely (inviscid, potential flow) or by ignoring inertia
completely (Stokes flow). Potential flow simulations do lead to crossstream
arrays, but the wakes needed for the drafting part of the fundamental
rearrangement mechanism in a fluidized suspension are absent.
Brady, et al. have created good techniques for dealing with the motion
of many particles in Stokes flow. These simulations are appropriate for
colloids at very small Reynolds number and they appear to successfully capture
the hydrodynamic interactions. On the other hand, the inertial mechanisms which
turn long bodies broadside on and control the lateral migration of particles do
not exist in Stokes flow.
Direct simulations of many particles at finite Reynolds numbers have only just
begun. Unverdi and Tryggvason introduced a front tracking/finite difference
method for computing the unsteady motion of drops and bubbles. In their work,
the drop surface is tracked by separate computational points that are moved on
a fixed grid by interpolating their velocity from the grid. These points form a
front that is used to keep the density and viscosity stratification sharp and
to calculate surface tension forces. Esmaelli (and Tryggvason) have computed
the rise of 16 bubbles in three dimensions and 144 and 324 twodimensional
bubbles in a doubly periodic domain at a Reynolds number near 2. These
simulations have not been adapted to solid particles and rigid walls, but they
do give rise to drafting, kissing and tumbling as in the case of sedimenting
solid spheres.
Johnson (and Tezduyar) did a threedimensional simulation using a spacetime
finite element method advocated by Hughes and Tezduyar, et al. of the
sedimentation of five solid spheres in a tube at a Reynolds number of 100.
Their simulation gives rise to drafting, kissing and tumbling into a plane
perpendicular to the fall.
The Minnesota group has done direct simulations of initial value problems for
the twodimensional motion of circular and elliptical particles in sedimenting,
Couette and Poiseuille flows of a Newtonian fluid at particle Reynolds numbers
in the hundreds. Recently, Hu has simulated the motion of 400 circular
particles in sedimenting and shear flows at particle Reynolds numbers in the
hundreds.
The first direct simulation of the initial value problem for the
twodimensional motion of circular and elliptical particles settling in a
viscoelastic fluid (Oldroyd~B) has recently been done by the Minnesota group.
A generalized Galerkin finite element scheme which
incorporates both the fluid and particle equations of motion into a
single coupled variational equation has been developed by Joseph's
group and Hu for Newtonian fluids in 2D and 3D domains. The formulation can
be easily extended to viscoelastic fluids. The hydrodynamic forces and torques
acting on the particles are eliminated in deriving the combined variational
equation, so need not be computed explicitly.
An arbitrary LagrangianEulerian (ALE) moving mesh technique has been adopted
to deal with the motion of the particles. In our implementation, the nodes on
the particle surface are assumed to move with the particle. The nodes in the
interior of the fluid are computed using Laplace's equation, to guarantee a
smoothly varying distribution of nodes. At each time step, the grid is updated
according to the motion of the particles and checked for element degeneration.
If unacceptable element distortion is detected, a new finite element grid is
generated and the flow fields are projected from the old grid to the new
grid.
Initially, the particles are positioned randomly in the fluid, with zero
velocity. The fluid is either at rest or flowing steadily around the particles.
The particles are then released and the motion of the combined fluid/particle
system is computed using the procedure described. In this scheme, the positions
of the particles and grid nodes are updated explicitly, while the velocities of
the fluid and the solid particles are determined implicitly.
This generalized ALE Galerkin finite element formulation gives rise to a set of
nonlinear algebraic equations which is solved via a quasiNewton scheme; the
resulting linear system is shown below.
For 3D domains with large numbers of particles, these systems will be
extremely large, requiring iterative solvers and matrixfree
preconditioning strategies.
Extending this approach to 3D domains will also require a robust 3D mesh
generator and an associated projection scheme for moving the mesh. We plan to
develop a parallel automatic grid generator based on the GHS3D generator
developed at INRIA, France. In work to date, it is assumed that the particles
do not collide, but that there is always a thin layer of fluid between
approaching particles, and between particles and retaining walls. The grid
generator must therefore possess a local refinement capability to handle grid
in these thin layers, which is essential to correctly model the ``particle
collision'' process.
Finally, in computing flows with large numbers of particles, it is often
desirable to use periodic boundary conditions in one or more directions because
it is a perfectly absorbing boundary condition, because particles near a
periodic boundary can interact with particles near the opposite periodic
boundary rather than with a solid wall, so that the same local physics applies
uniformly to all particles, and because no mathematical approximation is
involved.
At the periodic boundaries particles frequently leave the computational domain
and enter at the opposite periodic boundary. We will design the grid generator
to handle any number of pairs of periodic boundaries without introducing
artificial cuts.
The above approach will be further extended to viscoelastic fluids. In this
case, the constitutive equation relating the stress to the strain rate is more
complicated. Various constitutive models (with Oldroyd~B principal parts) will
be considered, including those in which the viscosity function and the
relaxation time may depend on the shear rate. While the general plan of the
computational scheme will be essentially identical to that for Newtonian
fluids, the computation of the fluid flow will need special attention. In
particular, we will need very effective linear and nonlinear system solvers,
because the problem will become dramatically larger (due to the larger number
of flow variables) and more illconditioned (due the strong coupling of the
velocity, pressure and stress fields), and because of the existence of thin
stress boundary layers.
Over the last ten years, several numerical schemes have been developed for
handling viscoelastic flows in simple twodimensional geometries (see for
example. To prevent numerical instability due to convection dominance in the
constitutive equations, we plan to develop an elasticviscous split stress
Galerkin leastsquares (EVSSGLS) finite element formulation. In the GLS
formulation, griddependent terms are added to the standard Galerkin equations.
These terms are functions of the residuals of the governing equations evaluated
elementwise, and are designed to enhance the stability of the Galerkin method
without degrading its consistency and accuracy. The GLS formulation has yet to
be applied to simulation of viscoelastic flows.
We also propose to investigate the simulation of
particle motions using a class of methods based on the principle of
embedded or fictitious domains. There are several ways to
design highperformance particle movers by generalizing the application of this
principle to the problem of fluid flow around fixed obstacles. Though particle
movers based on embedding methods have yet to be developed, they appear to be
natural and show great promise for development into very fast highperformance
particle movers. Ultimately, embedding methods may be even more powerful than
methods based on efficient remeshing of unstructured grids.
Embedding methods, or fictitious domain methods were introduced by Golub, et
al.\ in, and developed by Glowinski, et al. The idea is to embed an
irregular computational domain into a larger, simpler domain, and to specify
simple boundary conditions on its boundary. Thus, since the larger domain
admits a uniform grid, we can use fast elliptic and fast Stokes
solvers when applying these methods to the timedependent NavierStokes
problem, for example. Fictitious domain methods have been applied to linear and
nonlinear 2D and 3D elliptic and parabolic problems, to the 2D NavierStokes
problem, and to the 3D NavierStokes problem. It is also suitable for
flowrelated shape optimization problems.
To apply the fictitious domain approach to the problem of fluid flow around
fixed obstacles, we may take, for the larger simpler domain, the fluid domain
plus the interiors of the obstacles. That is, the fluid flow is
computed as if the obstacles were filled with fluid. The noslip
boundary condition on the obstacle boundaries is enforced as a constraint
equation using well chosen source terms which are either associated with
singularities distributed on the obstacle boundaries, or with body forces
distributed inside the obstacles. In the approach developed by Glowinski, the
unknown distributions of singularities are Lagrange multipliers
associated with the boundary conditions on the obstacle boundaries. It is
therefore natural to apply powerful variational principles to both the spatial
discretization by finite elements and to iterative solution methods.
Recent results show that fictitious domain methods can be successfully applied
to the NavierStokes flow around moving obstacles (i.e., particles)
when the motion of these obstacles is prescribed in advance
(oscillating airfoil in a channel).
Results by Buzbee indicate that the interface between the actual and simpler
domains was consistent with the same mesh used for both domains. In more recent
publications, it was shown that this consistency is unnecessary. In particular,
it was shown that if the actual boundary conditions on the interface are
enforced using Lagrange multipliers, it is important to have the discrete
multiplier space be associated with the intrinsic geometry of the original
(irregular) domain and not with the mesh of the larger simpler domain. This may
be somewhat complicated to implement on parallel architectures, but once done,
the scheme can easily be applied to problems with moving boundaries, such as
the surfaces of moving particles.
Fictitious domain methods can accommodate fixed gridsindeed,
uniform fixed gridsa definite advantage in dealing with moving
obstacles, despite the advances achieved with moving grids, since we can use
fast elliptic and fast Stokes solvers when applying such methods to the
timedependent NavierStokes problem. Based on boundary layer thickness
considerations and on the large number of particles assumed to fill the flow
domain, we can expect that a fine uniform grid will need to be used at large
Reynolds numbers.
In order to apply fictitious domain methods to the problem of simulating the
flow of a fluid with free particles (i.e., obstacles with
unprescribed motion), the hydrodynamic forces and torques exerted by
the fluid on the particles must be accurately evaluated. Once this has been
accomplished, solving the coupled equations of fluid and particle motion will
be fairly easy through approximate time stepping methods using, for example,
operator splitting techniques, such as the #tex2html_wrap_inline963#scheme.
The thetascheme is a particularly good time stepping method. In fact, recent
evaluations and comparisons done at the University of Heidelberg under the
supervision of R.~Rannacher (one of the leading world experts on the
approximation of the unsteady NavierStokes equations) show that the Glowinski
thetascheme offers the best compromise for accuracy, stability, robustness and
ease of implementation. The thetascheme is well suited to fictitious domain
methods, It is also well suited to the simulation of viscoelastic flow, which
is one of our main objectives. Research team members have been using the
#tex2html_wrap_inline971#scheme extensively in both Houston and
Minneapolis.
Coupling the #tex2html_wrap_inline973#scheme for the flow simulation to the
timediscretization scheme for the equations of particle motion will be easy
once the hydrodynamic forces and the contact forces between particles are
known.
The major computational task in fictitious domain methods, coupled with
operator splitting, is the solution of symmetric indefinite systems of the form
in each time step. Here, both A and B are sparse matrices with
A being symmetric positive definite and B of maximal column rank.
Iterative schemes for solving such indefinite systems are discussed below.
Both of the proposed computational schemes require
the solution of large systems of linear and nonlinear algebraic equations. For
3D simulations involving a large number of particles, these systems can be
extremely large. In simulations conducted by research team members, up
to 95% of the computational time of the entire simulation is consumed by the
linear system solver alone.
In both the generalized ALE Galerkin scheme and the embedding scheme, we must
solve large sparse linear systems for which direct methods based on Gaussian
elimination are too expensive in terms of both time and storage cost. For this
class of problems, iterative methods must be used, possibly in conjunction with
direct methods when we consider preconditioning strategies. Important building
blocks in the library of parallel system solvers to be developed for this
project are preconditioned iterative schemes for solving indefinite linear
systems described above. >br>
Such systems may be solved using a variety of iterative schemes. Of particular
interest is the family of the Uzawa schemes and their variations. Also, for the
case of nonsymmetric matrices A, one may resort to Krylov subspace
methods such as GMRES, CGS and variations. The main drawback of those
nonsymmetric iterative solvers is lack of robustness without effective
preconditioners. Developing effective preconditioners that are also suitable
for parallel computers will be one important focus of this research. Row
projection schemes, accelerated via the conjugate gradient algorithm, provide a
robust alternative to Krylov subspace methods. More work needs to be done to
increase the rate of convergence, however.
In the embedding scheme, outlined above, using the thetamethod together with a
uniform grid results in matrices A that are both structured and
symmetric positive definite. In such a case, an important part of the library
will be the rapid elliptic solvers developed mainly by Golub, as well Sameh, in
addition to the parallel conjugate gradient algorithm with a variety of
preconditioners. For the ALE Galerkin formulation, however, the matrix A
is nonsymmetric. We plan to develop tools for rapid prototyping of parallel
preconditioned solvers for both types of systems. Since the size of such linear
systems can be extremely large when handling 3D problems, we propose to study
preconditioners that do not require explicit storage of coefficient matrices.
In other words, we need to find preconditioners that do not depend on obtaining
approximate factorization of the systems involved. Moreover, we will develop
iterative solvers which share the property of direct solvers in the sense that
once they are used for solving the linear systems described with one righthand
side then they can be used at a much less cost for another righthand side,
assuming of course that the matrix of coefficients does not change, or changes
slowly, from one time step to another. Another motivation for considering
multiple righthand sides is efficiency in conducting parametric studies.
One approach to preconditioning without using approximate factorization of
the linear system under consideration, is the use of an inner iterative process
to affect the preconditioning step. For reasons of efficiency, this is only
done approximately. The question then arises as to how the accuracy of the
inner iteration affects the convergence of the outer iteration. For Chebyshev
outer iterations, this question has been studied in detail by Golub and his
collaborators. Using a conjugate gradienttype method as the outer iteration in
combination with an inner iteration method remains an open problem.
A related inner iteration scheme is one in which the outer iteration provides
information about the eigenvectors of the troublesome eigenvalues while the
inner iteration solves the projected linear system in which the projector
deflates those eigenvalues. These projectors become more effective as the outer
iterations proceed.
Finally, another approach for matrixfree preconditioning is via Newton's
iteration for obtaining the inverse of a matrix. Choosing an initial
approximation for the inverse that commutes with the matrix in question leads
to a recursive function for the successive approximation of the solution
vector. This approach is still in the experimental stage, and further
investigation is needed to explore its numerical stability and suitability for
parallelism.
The efficient handling of multiple righthand sides proves to be a challenge
to most iterative methods because one of their main advantages is that not much
information needs to be stored and therefore very little information is saved
from the solution of one system which can be reused for the solution of
subsequent systems.
Consider the situation where we have two righthand sides b and
c. Our goal is to solve the second system A x=c using
information gathered from the first system A x=b. Here we
consider the case where A is an n x n symmetric matrix but the
basic ideas apply to more general situations. Consider using the Lanczos method
for solving the linear systems since it is well known that the conjugate
gradient method is mathematically equivalent to the Lanczos method. Our idea is
easier to express in the Lanczos framework, but we expect to be able to derive
a conjugate gradient and GMRES version as well.
Project Home
 AEM Home 
Institute of Technology 
 Academics 
Research  People 
Information  Contact
AEM 
 Major Projects 
by Specialty 
Other Programs 
SolidLiquid Flows,
Computational Methods
Last updated October 16, 2000 