The extended finite element method for rigid particles in Stokes flow

A new method for the simulation of particulate flows, based on the extended finite element method (X‐FEM), is described. In this method, the particle surfaces need not conform to the finite element boundaries, so that moving particles can be simulated without remeshing. The near field form of the fluid flow about each particle is built into the finite element basis using a partition of unity enrichment, allowing the simple enforcement of boundary conditions and improved accuracy over other methods on a coarse mesh. We present a weak form of the equations of motion useful for the simulation of freely moving particles, and solve example problems for particles with prescribed and unknown velocities. Copyright © 2001 John Wiley & Sons, Ltd.


INTRODUCTION
A recurring theme in the study of the rheology of particulate suspensions is the use of the analytical solution for a single particle or small group of particles to aid in the construction of the many-particle solution. For example, for very dilute suspensions of rigid spheres [9], rigid ellipsoids [2], or spherical droplets [3], particle-particle interactions can be ignored, and the properties of the suspension can be derived in terms of the vanishing Reynolds number solution for linear ow past a single particle. For more concentrated solutions, as in Batchelor and Green's analysis of a suspension of rigid spheres [4], the two-sphere analytical solution is used to approximate particle-particle interactions. Analytical solutions for single particles have also been exploited in numerical simulations of particulate ows; the method of Stokesian dynamics [6] utilizes a grand mobility matrix, constructed by assembling solutions for pairs of particles, to compute the velocities of all particles in a suspension.
In the last few years, the state of the art in computers and numerical methods has made the direct numerical simulation of particle-laden ows possible. Recently, ÿnite element methods have been used to simulate particulate ows in great detail. These simulations can be divided into two classes: moving mesh simulations, in which the computational grid moves with the particles, and ÿxed mesh simulations, in which the particles move through a ÿxed, regular grid. For example, Hu [24] used an arbitrary Lagrangian-Eulerian moving ÿnite element mesh technique to simulate a large number of rigid particles. Johnson and Tezduyar [15] also used an unstructured moving mesh; their simulation used a stabilized space-time formulation to simulate up to 100 rigid spherical particles. For a ÿnite element problem of this size, the number of degrees of freedom required is extremely large, and the uid domain must be remeshed frequently to avoid element distortion. For example, Johnson and Tezduyar report mesh sizes of around 1.2 million elements and 240 000 nodes for their simulation of 100 particles, although the mesh size varies signiÿcantly due to remeshing and adaptivity. Fixed mesh methods help alleviate these di culties, but the enforcement of boundary conditions at particle surfaces (which do not necessarily coincide with element surfaces) can be di cult. Glowinski et al. [7] used a ÿctitious domain technique for particulate ows, applying a distributed Lagrange multiplier to enforce rigidity inside the particles. Both moving and ÿxed mesh methods require very ÿne meshes in order to capture ow details near the particles.
The ÿxed mesh methods listed above are related to methods in which a complex boundary, moving or stationary, is captured on a regular grid that does not conform to the boundary. In the immersed boundary method [12], the force and velocity at the uid-surface interface are interpolated using approximated Dirac delta functions. This method has been applied to the ow of suspensions [21; 1]. A similar method is the virtual boundary method [17], in which a feedback force is added at the boundary to constrain the velocity. For problems in solid mechanics, Donning and Liu [8] used Lagrange multipliers to enforce boundary conditions at a non-conforming boundary. In their study, the objective was to remove the locking constraints of plates and beams. They concluded that a very ÿne mesh is required to capture detailed solutions near the boundary.
In this paper we describe a method that utilizes an analytical solution near the particle surfaces, allowing for a solution on a coarser mesh than would otherwise be possible. The technique is based on the extended ÿnite element method (X-FEM), which was originally developed for crack growth problems without remeshing [18; 22]; the method uses a partition of unity enrichment [16] to allow knowledge about the partial di erential equation being solved to be included in the ÿnite element space. X-FEM has also been used for elastic solid problems with holes [5] and inclusions [14] which need not conform to the mesh. These geometrical features can pass arbitrarily through element interiors, because the set of basis functions includes discontinuities across the feature surface.
A similar partition of unity method, the multi-scale ÿnite element method, can be found in Liu et al. [19]. In this approach, the analytical wave and beam solutions were embedded in the numerical solutions by multiplying the linear ÿnite elements with the wave-number banded solutions. As a result, the selected band of wavenumber solutions was obtained. An alternative ÿnite element enrichment procedure with a mesh-free method can be found in Liu et al. [23].
The ÿrst class of uid ow problems to which we have applied this technique is the simulation of rigid circular particles in 2D viscous ow. Currently, the 2D problem is solved, In Section 2, we show how a combined weak form of the particle and uid equations of motion for Stokes ow can be derived. The velocity approximation is constructed using X-FEM to satisfy the boundary conditions at the particle surfaces (Section 3). The method is veriÿed by solving a simple problem with a known solution in Section 4, and examples are shown in Section 5.

X-FEM=FICTITIOUS DOMAIN METHOD FOR PARTICULATE STOKES FLOWS
Glowinski et al. [7] presented a ÿctitious domain method for the simulation of particulate ows which used a distributed Lagrange multiplier to enforce rigid body motion inside the particle domains. In our method, this ÿctitious domain method can be used without the need for any Lagrange multipliers, as X-FEM makes it simple to construct a solution space which automatically satisÿes the rigid body constraint inside the particles. Here we present a ÿctitious domain method similar to that of Glowinski et al, but simpliÿed for the case of Stokes ow using X-FEM.
The problem geometry 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 G about the particle centre.

Fluid motion
The uid motion u(x; y) is governed by the Stokes equations, with rigid body motion enforced at the particle surfaces: Figure 1. Fluid and particle system geometry.
where f is the viscosity of the uid, g is the acceleration due to gravity, A is the uid stress, u (t) is a prescribed outer boundary condition and r is the directed distance from the centre of particle .

Particle motion
Because the ow is in the Stokes regime, particle inertia is ignored and the total force and torque on the particle are zero. Then for each particle : where M is the mass of particle . The hydrodynamic force F and torque T on each particle are given by F is a short-range repulsive force between pairs of particles and between particles and the walls included in order to avoid collisions (see for example, Reference [7]). It is ignored in the examples reported here.

Combined weak form
In order to derive a weak form we deÿne a space of velocities, valid both inside and outside the particles, from which we will seek a solution. The test velocity space involves test function v(x; y; t) for the velocity, along with particle translations V and rotations^ . The ÿeld v(x; y; t) is constrained to match rigid body motion of the particles on the particle surfaces. The trial function space is then The pressure p(x; y) lies in the space Choosing velocity and pressure trial functions (u; U ; G ) ∈ V and p ∈ L 2 0 ( ), and test functions (v; V ;^ ) ∈ V 0 and q ∈ L 2 ( ), the weak form of (1) and (2) is (now using index notation): Integrating (7a) by parts, with use of Equation (3) and the constraints on V and V 0 , gives: Note that the forces and torques on the particle surfaces have cancelled because the force of the uid on the particles is equal and opposite to the force of the particles on the uid. This is an attractive property of this method, as it obviates the need to perform cumbersome surface integrals on the particles to calculate the force. This approach for dealing with uidsolid interaction should also be compared with the immersed boundary method [12], which requires an approximation to transfer the solid forces to the uid and vice versa. No such approximation is necessary here.
The ÿnal step in the derivation is to extend the integrals in (8) to the entire domain ∪ P(t), where P(t) is the total domain of all the particles. With this goal in mind, we note that letting s represent the density of the solid, we can write: In the ÿrst equality, we have used the fact that the integral of r over particle is zero. Substituting into (8), and using the fact that v (i; j) = 0 inside the particles because of the rigid body motion constraint, we arrive at the ÿnal form These equations can be solved to give the velocities everywhere in the domain at a given time step; this solution can then be used to update the particle positions at the next time step.

ENRICHMENT SCHEME
The advantage of X-FEM for this problem is that the combined velocity space V(t) given in (4) can be constructed without Lagrange multipliers to enforce particle rigidity. Instead, the total ÿnite element and enrichment solutions are formulated to satisfy this constraint exactly, while particle translational and rotational velocities remain unknowns which are part of the solution. This provides an advantage over the use of Lagrange multipliers, which increase the problem size and add complexity to the problem formulation and solution.
In our approach, the velocity and pressure solutions spaces are augmented by a partition of unity enrichment method [16; 19]. Supposing that the basis for a scalar variable w(x) is to be enriched with functions (x), the approximation for w(x) is where N I (x) are the standard ÿnite element shape functions, (x) are the enrichment functions, n E is the number of enrichment functions, and w I and a I are scalar coe cients. D is the set of nodes whose supports (domains of non-zero shape function) lie in the region D, which is the union of all subdomains in which the solution is to be enriched.
Using (10), the standard ÿnite element interpolation can be recovered by setting the a Ii coe cients to zero. Furthermore, any enrichment function ÿ (x) can be exactly reproduced in the region D by setting the coe cients a Iÿ to unity and all other coe cients to zero, i.e.
The ÿnite element shape functions N I (x) in the enrichment term in (10) preserve the sparsity of the resulting matrix problem; any a I contributes to w(x) only where N I (x) is non-zero, i.e. within the support of node I .

Velocity enrichment
The enrichment of a vector ÿeld such as velocity is similar to that of a scalar ÿeld (10), but requires special attention. Speciÿcally, we must decide whether to enrich each component of the vector separately with a scalar function (x), or to tie the components together by enriching with a vector function M (x). Enriching separately o ers more freedom in the solution, but leads to additional unknowns.
For this particular problem, we will enrich the velocity with vector functions u (x), since vector ÿelds which satisfy the Stokes equation and the incompressibility equation can be derived (see the appendix). The total velocity ÿeld is thus whereê i is the unit vector in the ith spatial direction. The vectors u * (x) are four fundamental solutions of the Stokes equations. They are derived in the appendix, and quoted here for convenience (along with the associated pressure ÿelds)

Pressure enrichment
In enriching the pressure ÿeld, modes that will lead to a singular matrix problem should be excluded from the enrichment basis. We note that (14c) and (16c) can be rewritten in Cartesian co-ordinates Because these linear ÿelds can already be exactly reproduced by the piecewise linear ÿnite element shape functions, their inclusion as enrichment functions would lead to ambiguity; i.e. more than one choice of ÿnite element and enrichment coe cients would lead to the same pressure solution. Therefore, we enrich only with the functions p * 1 (13c) and p * 3 (15c): The unknown pressure coe cients are p I and b I ; the superscript p on N p I (x) is a reminder that we need not choose the same ÿnite element shape functions for pressure as for velocity.

Enrichment of higher-order shape functions
Care must be taken in choosing ÿnite element spaces for velocity and pressure in problems with incompressibility constraints. These spaces should ideally satisfy the BabÄ uska-Brezzi condition [13]. To satisfy this condition for our 2D problem, we have used quadrilateral elements with nine velocity nodes and four pressure nodes, so that the velocity ÿeld is quadratic and the pressure ÿeld is linear. All examples solved in this work use square elements with regularly spaced nodes.
In selecting nodes to be enriched, we follow Mo es et al. [18], who for an elastic crack problem enrich with a discontinuous displacement ÿeld at all nodes whose support is intersected by the crack. This was done for a linear displacement ÿeld (four-node quad elements), but the extension to higher-order elements is straightforward. All nodes whose supports are intersected by a particle surface are enriched. The support of a node is deÿned as the region in which the node's shape function is non-zero. For an edge node, this region consists of the two elements which share that edge; for a node at the center of an element, the support is ( )-nodes with ÿnite element degrees of freedom set to match rigid body particle motion.
that single element. An example showing which nodes are enriched near the particle surfaces is shown in Figure 2.

Boundary conditions
The enrichment functions in Equations (13)-(16) were chosen to exactly satisfy the u = 0 boundary condition at the particle surface. Thus, for ow past a stationary particle, the zero boundary condition is met very simply: all ÿnite element velocity degrees of freedom which contribute to the solution inside the particle are set to zero. For example, in Figure 2 all boxed nodes are set to zero. For a moving particle, we construct a ÿnite element velocity ÿeld in the enrichment region equal to a rigid body translation and rotation of the particle; the other enrichment ÿelds vanish inside the particle and at the surface, so the total solution is equal to a rigid body translation and rotation inside the particle. To construct the proper ÿnite element ÿeld, we ÿrst express the velocity inside the particle as where r is measured from the centre of the particle and U and ! ! ! are the translational and rotational velocities of the particle, respectively. In the enrichment region D associated with particle , the ÿnite element ÿeld is The enrichment term in (20) was ÿrst introduced in Reference [11] to model the rigid rotation of a disk placed on a regular mesh which need not conform to the disk geometry. Note the introduction of the additional degrees of freedom ! i . In order to enforce the desired boundary condition, we simply set U i and ! i equal to the constant particle velocity and rotation. The beneÿt of decomposing the near-particle velocity ÿeld in this way is that the particle translation and rotation can be ÿxed or left as unknowns in the discrete equations.

Numerical integration
In order to integrate the weak form as accurately as possible, we subdivide elements which are intersected by particle surfaces into smaller integration cells. The algorithm we use is the same as that used by Mo es et al. in Reference [18].

METHOD VERIFICATION
In order to ensure that the method works properly, we would like to be able to compare simulation results for a simple problem with an exact solution. However, there is no analytical solution for pure Stokes ow for the problem of a 2D cylinder in a ow which is uniform at inÿnity. Furthermore, since we need to solve on a ÿnite domain, we must look at a problem with known boundary conditions at the outer boundary of the ow domain, rather than at r → ∞. The problem we choose for code veriÿcation is that of a cylinder of radius R at the centre of a square domain, where the velocity on the outer boundary of the domain is u = u * 1 ; v * 1 where u * 1 and v * 1 are taken from (13). The exact solution to this problem is simply u = u * 1 ; v * 1 on the entire domain, since this solution is known to satisfy the Stokes equations. This problem is attractive for code veriÿcation purposes because the enrichment part of the numerical solution is expected to be able to exactly capture the analytical solution.
Although all four velocity enrichments given in Equations (13)- (16) are solutions of the equations of motion and satisfy the boundary conditions, and are therefore valid enrichments, we ÿnd in practice that enriching with all four enrichment ÿelds leads to a poorly conditioned matrix problem. This problem is alleviated by enriching only with the sets of functions given in (13) and (15). The other two sets grow as r 2 far from the cylinder, and are expected to be unimportant in real ow problems. We ÿnd that enriching with only two sets of functions instead of four leads to very little loss in accuracy, and we have therefore used only these two set of functions for all of the problems solved in this paper.
The problem is solved on a square domain of side length 2, with a cylinder of radius 0:2 centred at the origin. Velocity contours for a standard ÿnite element solution and the X-FEM solution are plotted in Figure 3. The ÿnite element mesh has 1042 triangular elements, while the X-FEM mesh is a 20 × 20 array of square elements.
The L 2 error in the velocity solution both for the standard ÿnite element method with a conforming mesh and for X-FEM is plotted in Figure 4. The error is plotted against N 1=2 where N is the total number of degrees of freedom in the simulation for both velocity and pressure. The X-FEM solution shows the same rate of convergence as the conforming ÿnite element solution, but with a smaller error for a given number of degrees of freedom. We conclude that for this example the extended ÿnite element method can achieve a given level of accuracy more e ciently than standard ÿnite elements, since it requires fewer degrees of freedom.

Particles with prescribed motions
The X-FEM technique described in the previous sections has been applied to several problems of Stokes ow past arrays of ÿxed or moving particles. All solutions reported here are for square domains of side length 2:0, centred at the origin. In all cases, solutions are plotted both for standard ÿnite element meshes that confom to particle boundaries and for an X-FEM mesh consisting of a regular array of square elements. 5.1.1. Uniform ow past a stationary particle. A particle of radius 0:2 is located at the origin, with a boundary condition of u = e x applied at the outer boundary. The X-FEM mesh is an array of 20 × 20 square elements. Velocity contours are shown in Figure 5, and velocities on the cross-section y = 0:21 are shown in Figure 6. Note in Figure 6 that the FEM and X-FEM solutions are indistinguishable except very near the particle surface (x = 0), where the X-FEM solution picks up more detail.
N is the total number of velocity and pressure degrees of freedom.

5.1.3.
Flow ÿeld due to four rotating particles. Four particles of radius R = 0:2 are located at (x; y) = (±0:3; ±0:3). The particles in the ÿrst and third quadrants rotate with angular velocity 1:0, while those in the second and fourth quadrants have angular velocity −1:0. The X-FEM mesh consists of 30 × 30 square elements. Velocity contours are shown in Figure 9, and velocities on the cross-section y = 0:0 are plotted in Figure 10.

Freely moving particles in a gravitational ÿeld
As simple test cases of the simulation of freely moving particles, we consider a pair of related problems for which analytical solutions are available. In the ÿrst of these, an inÿnitely long two-dimensional channel of width 2L is oriented such that gravity acts parallel to the channel walls. A solid circular particle of radius R and density s is located at the centre of the channel (Figure 11a). The uid is assumed to have zero density; this is done in order to match the analytical solution, which is derived for a force acting on the particle but not on the uid. The variable to be solved for is the resulting downward particle velocity U . The asymptotic solution to this problem for R=L small is [20] U = s gR 2 f 1 (R=L) 4 (21a)   The second problem is similar, but with gravity acting perpendicular to the channel walls ( Figure 11b). In this case, the asymptotic solution is [20] U = s gR 2 f 2 (R=L) 4 (22a) It should be stressed that the numerical solutions we have computed use the same enrichment scheme as the previous examples. The enrichment functions are those derived in the appendix, and contain no information from the asymptotic analyses mentioned above. Our numerical solutions are on a computational domain with length of 6L, with boundary condition u = 0 at the channel ends as well as the channel walls. The domain was meshed with a rectangular grid with 10 square elements across the channel width by 30 elements down the length. The computed values of f 1 (R=L) and f 2 (R=L) are shown in Figures 12 and  13, respectively. As can be seen, the computed solutions match the asymptotic solutions very      well, especially where R=L is small, which is where the asymptotic solution is most accurate. Note that the smallest value of R=L computed is 0:02; for this case the particle radius is just one ÿfth the size of the element width for the 10 × 30 mesh used. The method is applicable to the motion of a large array of particles on a relatively coarse mesh. Figure 14 shows the y-velocity in the uid for an array of 50 sedimenting particles in a channel of length 3:0 and width 1:2. The ÿnite element mesh consists of 1440 square elements. The falling particles cause a ow in the positive y-direction between columns of the array due to the conservation of mass; the total ux through any horizontal cross section must be zero.

SUMMARY AND CONCLUSIONS
We have presented a method for the simulation of solid particles in a 2D ow in which the computational mesh need not conform to the particle surfaces. The partition of unity enrichment scheme we use is designed to incorporate the analytical solution for ow past a circular particle in the region near the particle surface. We have also formulated the problem so that the particle translational and rotational velocities are unknown, making the simulation of a suspension of freely moving particles possible. The method remains accurate even for particle sizes on the order of or smaller than the size of an element in the mesh.
The eventual goal of this work is the simulation of large numbers of particles in a suspension ow. This will require the extension of the method presented here in several directions. The simulation of three dimensional problems with this method is straightforward; in fact, there are more analytical solutions in three dimensions than in two available in the literature for use as enrichment functions. One potential complication in three dimensions is one of computational geometry; elements are partitioned across particle surfaces for accurate integration, and the partitioning of a 2D element by a 1D surface as is done here is simpler than the partitioning of a 3D element by a 2D surface. However, methods are available which can be applied to this task.
The method used to apply boundary conditions in this work may require modiÿcation for more concentrated suspensions in which particles collide or nearly collide. Enrichment functions based on solutions to two-sphere problems can be used instead of the analytical solutions given in Section 6. However, a ÿnite element degree of freedom cannot be tied to the rigid body velocities of all particles in the vicinity as in Section 3.4, since this would force all particles in the vicinity of a single node to move with the same velocity. For the treatment of closely spaced particles, one can imagine all particles in the ow would be tied together and move as a block, which is obviously incorrect. A more general formulation is needed, such as one in which the ÿnite element degrees of freedom capture a velocity ÿeld which is averaged in some sense and an enrichment solution which represents the perturbation due to the presence of the particles.
Although the formulation and examples given in this work are limited to slow viscous ow, the method is applicable to ows of higher Reynolds number as well. It should be recognized that the analysis of Section 6 is valid for low particle Reynolds number Re p , i.e. the Reynolds number based on particle size. A ow with a much higher Reynolds number based on ow domain dimensions can also be enriched with the given functions provided that the particles are small enough. Higher Re p of order unity call for enrichment based on asymptotic solutions such as those of Proudman and Pearson [10]. For even higher Re p it may be necessary to simply choose an enrichment function which satisÿes the boundary conditions at the particle surface and rely on the reÿnement of the mesh to pick up details like ow separation and trailing wakes.
None of these obstacles in the path to large-scale particulate ow simulations seem insurmountable. Suspension ow is an excellent example of a problem in which multiple scale phenomena are at work, and it seems natural to apply a technique such as the one outlined which allows the long-range, large-scale e ects to be computed while analytical models are employed to account for the near-ÿeld, small-scale perturbations. In this way, the strengths of both numerical and analytical methods can be exploited.

APPENDIX: ANALYTICAL SOLUTIONS OF STOKES FLOW PAST A CYLINDER
It is well known that the problem of two-dimensional Stokes ow past a single cylinder, with no-slip boundary conditions on the cylinder and uniform ow at r = ∞, has no solution. The di culty lies in the inappropriateness of the fundamental assumption that viscous forces are much larger than inertial forces. Very far from the cylinder, where velocity gradients are small, viscous forces become negligible and the assumption is violated; the problem must then be approached using either an approximation like Oseen's (see for instance, Reference [25]) or by matched asymptotics [10].
In the present case, we simulate only ÿnite domains, for which the Stokes equations do have a solution. Furthermore, since we enrich the ÿnite element method only near the cylinder, we can ignore the loss of validity of the solution far from the cylinder. Rather, we will look for solutions which satisfy the boundary conditions at the cylinder surface. One can think of the method as a matching problem, in which the form of the analytical solution near the surface is used as the inner solution, and the outer solution is obtained by ÿnite elements. Matching is done automatically through the Galerkin weak form.
The equations for an incompressible, very viscous uid are −∇ ∇ ∇p + ∇ 2 u = 0 (A1) ∇ ∇ ∇ · u = 0 (A2) The problem to be solved is that of a cylinder of radius R located at the origin. The ow far from the cylinder has uniform velocity U in the x-direction, so the boundary conditions are where r = x 2 + y 2 and e x is the unit vector in the x-direction. The problem will be solved in terms of a streamfunction (x), deÿned such that u = 1 r @ @Â e r − @ @r e Â (A5) This form of the velocity automatically satisÿes the incompressibility condition (A2). Taking the divergence of (A1) gives the biharmonic equation: with boundary conditions on : at r = R : @ @r = 0; @ @Â = 0 (A7) at r → ∞ : @ @r = U sin Â; @ @Â = Ur cos Â (A8) Based on the boundary conditions as r → ∞ (A8), we will look for solutions of the form (r; Â) = f(r) sin Â (A9) for which (A6) reduces to 1 r @ @r r @ @r − 1 r 2 2 f(r) = 0 (A10) Integration of (A9) gives a solution for (r; Â) in terms of the unknown constants A-D: (r; Â) = Ar 3 + Br ln r + Cr + D r sin Â (A11) Because the velocities are in terms of the derivatives of the streamfunction (A5), (r; Â) is unique only up to an additive constant. We are, therefore, free to satisfy (A7) by requiring that (r = R; Â) = 0, leading to the following two independent solutions for : Unfortunately, no linear combination of 1 and 2 satisÿes the boundary conditions at