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 particle-particle and particle-wall 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 3-D 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 matrix-free preconditioners, and to implement them on parallel architectures.

Existing Literature

Recent computational approaches to solid-liquid 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 cross-stream 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 two-dimensional 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 three-dimensional simulation using a space-time 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 two-dimensional 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 two-dimensional motion of circular and elliptical particles settling in a viscoelastic fluid (Oldroyd~B) has recently been done by the Minnesota group.

Generalized ALE Galerkin Finite Element Method

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 2-D and 3-D 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 Lagrangian-Eulerian (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 quasi-Newton scheme; the resulting linear system is shown below.


For 3-D domains with large numbers of particles, these systems will be extremely large, requiring iterative solvers and matrix-free preconditioning strategies.

Extending this approach to 3-D domains will also require a robust 3-D mesh generator and an associated projection scheme for moving the mesh. We plan to develop a parallel automatic grid generator based on the GHS-3D 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 ill-conditioned (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 two-dimensional geometries (see for example. To prevent numerical instability due to convection dominance in the constitutive equations, we plan to develop an elastic-viscous split stress Galerkin least-squares (EVSS-GLS) finite element formulation. In the GLS formulation, grid-dependent terms are added to the standard Galerkin equations. These terms are functions of the residuals of the governing equations evaluated element-wise, 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.

Embedding Methods

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 high-performance 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 high-performance 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 time-dependent Navier-Stokes problem, for example. Fictitious domain methods have been applied to linear and nonlinear 2-D and 3-D elliptic and parabolic problems, to the 2-D Navier-Stokes problem, and to the 3-D Navier-Stokes problem. It is also suitable for flow-related 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 no-slip 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 Navier-Stokes 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 grids---indeed, uniform fixed grids---a 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 time-dependent Navier-Stokes 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 theta-scheme 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 Navier-Stokes equations) show that the Glowinski theta-scheme offers the best compromise for accuracy, stability, robustness and ease of implementation. The theta-scheme 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 time-discretization 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.

Linear System Solvers

Both of the proposed computational schemes require the solution of large systems of linear and nonlinear algebraic equations. For 3-D 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 theta-method 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 3-D 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 right-hand side then they can be used at a much less cost for another right-hand 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 right-hand sides is efficiency in conducting parametric studies.

Preconditioners

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 gradient-type 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 matrix-free 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.

Multiple Right-Hand Sides

The efficient handling of multiple right-hand 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 right-hand 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 |

Solid-Liquid Flows, Computational Methods