INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING
Int. J. Numer. Meth. Engng 2003; 56:1261–1289 (DOI: 10.1002/nme.608)
Particulate
ow simulations using lubrication
theory solution enrichment
G. J. Wagner
∗;†
, S. Ghosal and W. K. Liu
Department of Mechanical Engineering; Northwestern University; 2145 Sheridan Road;
Evanston; IL 60208; U.S.A.
SUMMARY
A technique for the numerical simulation of suspensions of particles in
uid based on the extended
nite element method (X-FEM) is developed. In this method, the particle surfaces need not conform to
the nite element boundaries, so that moving particles can be simulated without remeshing. The nite
element basis is enriched with the Stokes
ow solution for
ow past a single particle and the lubrication
theory solution for
ow between particles. The latter enrichment allows the simulation of particles that
come arbitrarily close together without rening the mesh in the gap between them. Example problems
illustrating both types of enrichment are shown, along with a study of a 50% solution in channel
ow.
Copyright ? 2003 John Wiley & Sons, Ltd.
KEY WORDS: suspension
ow; enrichment methods; partition of unity; lubrication theory
1. INTRODUCTION
In a previous paper [1], we presented a method for the simulation of particulate
ows that
allows moving particles to be simulated on a xed grid without remeshing. Like the ctitious
domain method of Glowinski et al. [2], our method utilizes a global momentum equation
that is valid in both the
uid and solid domains, allowing the particle surfaces to cut through
elements. This approach also obviates the need to explicitly compute the hydrodynamic forces
and torques on each particle in order to solve for its velocity and rotation. In order to apply
boundary conditions on the surface of each particle so that the
uid and solid velocities are
equal at the interface, we utilize the extended nite element method (X-FEM) [3, 4] to enrich
the solution with the known Stokes
ow solutions for
ow past a particle. The incorporation
of these solutions, together with the ability of the linear nite element shape functions to
∗
Correspondence to: G. J. Wagner, Department of Mechanical Engineering, Northwestern University, 2145 Sheridan
Road, Evanston, IL 60208, U.S.A.
†
E-mail: g-wagner@northwestern.edu
Contract=grant sponsor: DOD National Defence Science
Contract=grant sponsor: National Science Foundation
Received 26 September 2001
Revised 30 January 2002
Copyright
?
2003 John Wiley & Sons, Ltd. Accepted 16 May 2002

1262 G. J. WAGNER, S. GHOSAL AND W. K. LIU
represent rigid body motion, allow the construction of a trial function basis that satises the
proper boundary conditions on the particle surfaces.
Because the Stokes
ow velocity solutions used for enrichment in Reference [1] are designed
to go to zero on the surface of a single particle, boundary conditions cannot be properly
enforced in cases where a node has its domain of in
uence cut by more than one particle.
Particle collisions (or near collisions) will occur in any real suspension, and so a successful
simulation method must be able to handle the case of two particles approaching each other
closely. In the current work, we greatly improve the method presented in Reference [1]
by enriching the basis with lubrication theory solutions, valid for small sizes of the gap
between particles, in regions between pairs of particles and between particles and walls.
These regions are characterized by high-velocity gradients and large pressure, and use of
a standard computational technique like nite elements requires substantial local renement
of the mesh in order to accurately capture the solution [5]. Because particles can approach
each other arbitrarily closely in a real solution, there is no lower limit to the element size
in such a scheme unless an articial repelling force is added between particles. Since the
macroscopic rheological properties of a suspension are dominated by the forces between pairs
of particles in close proximity (see for example Reference [6]), failure to solve for these
localized forces may lead to large errors in the full solution. Our method has the advantage
of giving accurate solutions in these regions even when the size of the gap between particles
is much smaller than the size of an element in the mesh. Mesh renement is unnecessary,
and particles can approach each other arbitrarily closely. In addition, we use the fact that the
lubrication solution for velocity and pressure in the region between particles can be completely
described through just seven degrees of freedom in the 2D case, greatly reducing the total
number of degrees of freedom necessary for the standard nite element approach.
In Section 2, we outline the derivation of the combined
uid=particle weak form. In Sec-
tions 3 and 4, the method of enrichment using X-FEM is presented. After discussing important
implementation issues in Section 5, we demonstrate the method through several example prob-
lems in Section 6, including the study of a 50% solution of circular particles in channel
ow.
We conclude with a discussion in Section 7.
2. COMBINED FLUID=PARTICLE WEAK FORM
In Reference [1], a weak form that combines the equations of motion in both the
uid and
the solid was derived. The problem geometry, comprising a set of suspended particles and
a surrounding
uid, is shown in Figure 1. The
uid domain  is bounded externally by
surface . Each particle
has interior domain P


(t), bounded externally by surface S


(t). The
unit normal at the particle surface pointing into the
uid is denoted n
∗
, and that pointing
inward by n. At a given time t, the particle
has translational velocity U


(measured at the
particle centre) and rotation



about the particle centre.
The
uid motion u(x; y) is governed by the Stokes equations, with rigid body motion
enforced at the particle surfaces:

f
g+∇ · = 0 in  (1a)
∇ · u=0 in  (1b)
Copyright ? 2003 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng 2003; 56:1261–1289