In this paper, a massively parallel implementation of the boundary element method to study particle transport in Stokes flow is discussed. The numerical algorithm couples the quasistatic Stokes equations for the fluid with kinematic and equilibrium equations for the particles. The formation and assembly of the discretized boundary element equations is based on the torus-wrap mapping as opposed to the more traditional row- or column-wrap mappings. The equation set is solved using a block Jacobi iteration method. Results are shown for an example application problem, which requires solving a dense system of 6240 equations more than 1200 times

Ingber, M. S.

Mondy, L. A.

Womble, D. E.

UNT Digital Library

Simulations of Hydrodynamic Interactions Among Immersed Par- ticles in Stokes Flow Using a Massively Parallel Computer M. S. Ingber*, D. E. Wamble**, L. A. Mondy*** *Department of Mechanical Engineering, University of New Mexico, Albu- querque, New Mexico 87131 **Applied and Numerical Mathematics Department, Sandia National Labo- ratories, Albuquerque, NM 87185-1110 ***Energetic and Multiphase Processes Department, Sandia National Labo- ratories, Albuquerque, NM 871 85-083.4 ABSTRACT In this paper, a massively parallel implementation of the boundary element method to study particle transport in Stokes flow is discussed. The numer- ical algorithm couples the quasistatic Stokes equations for the fluid with kinematic and equilibrium equations for the particles. The formation and assembly of the discretized boundary element equations is based on the torus-wrap mapping as opposed to the more traditional row- or column- wrap mappings. The equation set is solved using a block Jacobi iteration method. Results are shown for an example application problem, which re- quires solving a dense system of 6240 equations more than 1200 times. INTRODUCTION The ability to model the hydrodynamic interactions among immersed par- ticles in Stokes flow could impact many areas of scientific and engineer- ing interest, including pipeline transport of slurries, petroleum recovery, composite materials processing, blood flow, sediment transport and others. However, the macroscopic properties of suspensions, such as sedimentation rates, self-diffusion coefficients, and rheological properties, have been diffi- cult to predict because of the complex way in which the microstructural me- chanics influences macroscopic properties of the suspension. Furthermore, the complexity of discretizing a domain that includes moving particles has typically precluded the solution of the many-bodied problem using stan- dard domain techniques. Recently, the boundary element method (BEM) has been established as an effective means for the analysis of this class of problems [l-31. However, as the number of suspended particles included DISCLAIMER This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, make any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof. DISCLAIMER Portions of this document may be illegible in electronic image products. Images are produced from the best available original document. I in a numerical simulation is increased, so are the computational demands. Massively parallel computer architectures provide the resources necessary to accommodate these demands. Several BEM formulations on Multiple In- struction Multiple Data (MIMD) computer architectures have been success- fully implemented [4-61. Most of these implementations have been based on either distributing collocation nodes or boundary elements among the pr& cessors so that either rows or columns in the'discretized BEM equations are assembled in parallel. Recently, Ingber, Womble, and Mondy [7] developed a parallel algorithm based on the torus-wrap mapping in which rectangular blocks of the discretized BEM equations were generated in parallel. They were able to show that the torus-wrap mapping could substantially reduce the time required to perform an LU decomposition during the solution phase of the algorithm as compared to row-wrapped or column-wrapped mappings because of the reduction in communication volume. Nevertheless, for large problems, the solution of the discretized equations dominated the overall CPU effort. Several iterative equation solvers have been applied to discretized bound- ary element equations [8-91. A difficulty arising in most boundary element formulations is that the resulting system of equations is usually nonsym- metric. Although the systems may not be positive definite, they often have relatively large diagonal elements because of the singular nature of the kernel functions, regardless of whether the associated integral equations are of the first or second kind. Some iterative linear equation solvers that have been successfully implemented with boundary element methods include block Ja- cobi iteration, generalized minimal residual (GMRES), conjugate gradient squared (CGS), and stabilized bi-conjugate gradient (Bi-CGSTAB) meth- ods. In this paper, we extend our previous parallel boundary element formu- lations [6,7] in order to apply them to time-evolving particulate systems. In some cases, we have included several thousand time steps in a numerical simulation. We consider a block Jacobi iteration scheme to reduce the CPU effort involved in solving the discretized boundary element equations. This method is similar to the one discussed by Phan-Thien and Tullock [lo] in that it essentially makes use of a domain decomposition. However, unlike Phan-Thien and Tullock's solver, each subdomain consists of both particles and walls and also includes kinematic and equilibrium equations. The &- ciency of the parallel boundary element formulation is demonstrated in an example application. NUMERICAL ALGORITHM The numerical algorithm couples the quasistatic Stokes equations for the fluid with kinematic and equilibrium equations for the particles. The qua- sistatic Stokes equations are transformed into velocity boundary integral equations using standard techniques. The kinematic’ equations relate parti- cle surface velocities to particle centroidal linear and angular velocities. The equilibrium equations relate the particle surface tractions to the imposed body forces or moments. The discretized set of equations corresponds to a grand resistance matrix formulation [113 in that the particle velocities are determined solely as a function of the problem geometry. Once the particle velocities are determined, the particles are repositioned in space using a variable-time-step Runge-Kutta method. Details of this procedure can be found in Dingman et al. [12]. A parallel version of this algorithm for implementation on MIMD com- puter architectures was developed by the authors based on the torus-wrap mapping [7]. In the torus-wrap mapping, each processor is assigned dense rectangular submatrices to assemble. The actual dimensions of the rectan- gular blocks can be chosen as input to the computer program. In particular, row-wrap and column-wrap schemes in which collocation nodes or bound- ary elements are mapped to processors can be recovered as special cases. The torus-wrap mapping was shown capable of minimizing communication volume during the direct equation solving portion of the algorithm by choos- ing the proper dimensions of the rectangular block. Nevertheless, the LU factorization and back solve of the linear equations constituted the largest portion of the overall CPU effort for large problems. When solving a time-dependent problem, the solution at a time t typi- cally provides a good estimate for the solution at time t + At. The problem with a solution algorithm based on a direct method such as LU factor- ization is that it cannot take advantage of a good initial estimate. In an iterative method, on the other hand, a good initial guess of the solution can considerably reduce the solution time. One of the problems in using iterative methods for boundary element methods is that the matrices are dense. This means that a single matrix- vector multiplication requires O(n2) operations. Hence to be effective, an iterative method must require relatively few iterations to compute a so- lution. Another problem that arises in our boundary element application is that the matrices are typically not diagonally dominant and are often indefinite. These problems may be alleviated somewhat, but not eliminated, by reordering and scaling the matrix entries. If the desired ordering is not known a priori, the matrix must be reordered after it is generated. This is a very expensive operation on a distributed memory massively parallel computer. On the other hand, once we know a reordering, we can design an effective iterative method based on matrix splittings (e.g., block Jacobi, block SOR) or an effective preconditioner for Krylov subspace methods (e.g., GMRES, Bi-CGSTAB, CGS). The most effective reordering of the matrix that we have found is to collect in a diagonal block the entries corresponding to one particle’s self- Figure 1: Grayscale representation of the matrix corresponding to four par- ticles moving in a container. Darker entries are small, lighter entries are large. Entries to be collected into a single diagonal block are outlined. Other diagonal blocks are collected in a similar pattern. interactions, that particle’s interactions with a subset of the wall elements, and that particle’s rigid body interactions. As an example, Figure 1 shows a grayscale image of the values (dark entries are small and light entries are large) in the matrix corresponding to four particles moving in a container, when the natural ordering of equations and unknowns is used. Note that there are large off-diagonal entries and some diagonal elements are zero (cor- responding to the particles’ rigid body interactions). The matrix subblocks that must be collected into a diagonal block are outlined in white. We note then that each diagonal block captures the structure of the full matrix. Once the matrix is reordered as shown (or generated in this order), the LU factors of the diagonal blocks can be computed and used to give a convergent block Jacobi iteration or a block Jacobi preconditioner for a Krylov method. To date we have implemented only the block Jacobi iteration. Finally, we remark here that the diagonal blocks described above are relatively small (approximately 50 x SO), +d can easily be collected onto one processor. This leads to an efficient preconditioner. Most of the remaining computations lie in the matrix-vector multiplication subroutine that is the heart of iterative methods. This too can be parallelized efficiently (although not as efficiently as LU factorization). The block Jacobi iteration was first tested on a problem consisting of 4 spheres inside a closed cylinder of diameter 10 and height 10 (top at z = 5, bottom at z = -5). The walls of the cylinder rotated with velocity ug = L?TZ, simulating a parallel plate viscometer in which the top and bottom plates rotate at equal but opposite rates. Each Jacobi iteration was stopped when the residual was less than 1 x The time integrator was given by a variable-time-step second-order Runge-Kutta method. The average number of iterations required for convergence was approximately 15. The simulatiolri was run for three revolutions of the top and bottom plate. After the three revolutions, the solutions obtained, with the block Jacobi iteration and with the direct LU solver were identical in terms of the locations of the spheres to 3 significant figures. f i l l  optimization of the block Jacobi method has not been completed to date; however, preliminary results indicate that it will improve the performance of our boundary element method. APPLICATION PROBLEM We consider spinning ball rheometry as an example numerical appli- cation problem to show the capabilities of this parallel boundary element implementation. Our spinning ball rheometer consists of a suspension of particles placed inside a cylinder with a rotating sphere placed at the center of the cylinder. The effective viscosity of the suspension can be inferred from the torque required to maintain a constant angular velocity of the spinning ball. In application, the particles tend to  migrate away from the spinning ball because of gradients in the local shear rate, and this results in a reduction in the measured torque over time. Modeling of such migration of particles has recently become a topic of considerable interest [13-151. In this example problem, we consider 89 neutrally buoyant spheres placed in a monolayer about a central ball spinning within a stationary cylinder. The original locations of the 89 spheres, illustrated on the left side of Fig- ure 2, were determined using a random number generator. The locations of the 89 spheres after the spinning ball rotated (counterclockwise) 100 revolu- tions are shown on the right side of Figure 2. The spheres have aggregated somewhat and, on average, have tended to migrate away from the spinning ball. For example, the average of (the magnitude of) the radial positions T of the inner 10 spheres increases from about 1.2 initially to about 1.3 after 100 revolutions, where r=O represents the cylinder axis and center of the spinning ball. The torque on the spinning ball is plotted as a function of the number of revolutions in Figure 3. The torque fluctuates as the spin- ning ball interacts with suspended particles. Fluctuations are also seen in actual experiments, as is the general trend for the torque to diminish as the dynamic system evolves [16]. CONCLUSIONS A parallel, boundary element method to study the hydrodynamic interac- tions of immersed particles in Stokes flow has been presented. The formation Fi&e 2: Initial (left) and final (right) configuration spinning ball rheometer. 0.660 0.655 0.650 e c :: 5 0.640 0.635 0.650 0.625 of the spheres in the V J t 1 I 10 20 30 4 0  50 60  70 80 90 ' Revolutions 0 Figure 3: Torque on the spinning ball as a function of the number of revo- lutions. and assembly of the discretized boundary element equations is based on the torus-wrap mapping as opposed to the more traditional row- or column- wrap mappings. The equation set which is’ comprised of the discretized boundary element equations coupled with kinematic and equilibrium equa- tions is solved using a block Jacobi iteration method. The results shown in the application problem were generated by solving a dense system of 6240 equations over 1200 times. The efficiencies provided by this parallel bound- ary element implementation allowed this computation to be performed. ACKNOWLEDGMENTS This work was sponsored by the U. S. Department of Energy at Sandia National Laboratories under Contract DE-AC0494AL85000. This work was funded in part by the U. S. Department of Energy, Office of Energy Research, Division of Engineering and Geosciences. The authors are grateful for the computational time provided by Sandia National Laboratories on the 1840-node Paragon. REFERENCES 1. Ban-Cong, T. and Phan-Thien, N., ‘Stokes problems of multiparticle systems: a numerical method for arbitraxy flows,’ Phys. Fluids A,  Vol. 1(3), pp. 453-461, 1989. 2. Karrila, S. J. and Kim, S., ‘Integral equation of the second kind for Stokes flow: Direct solution for physical variables and removal of in- herent accuracy limitation,’ Chem. Engrg. Comm., Vol. 82, pp. 123-161, 1989. 3. Ingber, M. S., ‘Dynamic simulation of the hydrodynamic interaction among immersed particles in Stokes flow,’ Int. J. Num. Meth. FZuids, Vol. 10, pp. 791-809, 1990. 4. Karrila, S. J., F’uentes, Y. O., and Kim, S., ‘Parallel computational strategies for hydrodynamic interactions between rigid particles of ar- bitrary shape in a viscous fluid,’ J. Rheol., Vol. 33, pp. 913-947, 1989. 5. Kreienmeyer, M. and Stein, E., ‘Parallel implementation of the bound- ary element method for linear elastic problems on a MIMD parallel computer,’ Comp. Mech., Vol. 15, pp. 342-349, 1995. 6. Ingber, M. S. and Womble, D. E., ‘Hindered settling computations using a parallel boundary element method,’ in Parallel Computing in Multiphase Flow Systems, ASME FED-Vol. 199, Eds. Kim, S., Karniadakis, G., and Vernon, M. D., ASME, New York, pp. 53-60, 1994. 7. Ingber, M. S., Womble, D. E., and Mondy, L. .A., ‘A parallel bound- ary element formulation for determining effective properties of het- erogeneous media,’ Int. J. Num. Meths. Engrg., Vol. 37(22), pp. 3905-3920, 1994. 8. Atkinson, K. A. and Graham, I. G., ‘Iterative solution of linear sys= tems arising from the boundary integral method,’ SIAM J. Sci. Stat. Cornput., Vol 13, pp. 694722, 1992.. 9. Kane, J. H., Prasad, K. G., and Keyes, D. E., ‘Preconditioned Krylov equation solvers in elastoplastic boundary element analysis,’ Eng. Anal. BE, Vol. 14(1), pp. 3-14, 1994. 10. Phan-Thien, N. and Tullock, D., ‘Completed double layer boundary element method in elasticity and Stokes flow: Distributed computing through PVM,’ Comp. Mech., Vol. 14, pp. 370-383, 1994. 11. Brenner, H. and O’Neill, M. E., ‘On the Stokes resistance of multi- particle systems in a linear shear field,’ Chem. Eng. Sci., Vol. 27, pp. 1421-1439, 1972. 12. Dingman, S., hgber, M. S., Mondy L. A., Abbott, J., and Brenner, H., ‘Particle tracking in three-dimensional Stokes flow,’ J. Rheol. Vol. 36, pp. 413-440, 1992. 13. Phillips, R. J., Armstrong, R. C., Brown, R. A., Graham, A. L., and Abbott, J. R., ‘A constitutive equation for concentrated suspensions that accounts for shear-induced particle migration,’ Phys. FZuids A,  Vol. 4(1), pp. 30-40, 1992. 14. Koh, C. J., Hookham, P., and Leal, L. G., ‘An experimental investi- gation of concentrated suspension flows in a rectangular channel,’ J. Fluid Mech., Vol. 266, pp. 1-32, 1994. 15. Nott, P. R. and Brady, J. F., ‘Pressure-driven flow of suspensions: simulation and theory,’ J. FZuid Mech., Vol. 275, pp. 157-199, 1994. 16. Mondy, L. A., Graham, A. L., and Brenner, H., ‘Spinning Ball Rheom- etry in Concentrated Suspensions,’ presented at the 66th Annual Meet- ing of the Society of Rheology, Philadelphia, PA, 1994. 

Simulations of hydrodynamic interactions among immersed particles in stokes flow using a massively parallel computer

https://digital.library.unt.edu/ark:/67531/metadc708265/m2/1/high_res_d/69125.pdf

Simulations of hydrodynamic interactions among immersed particles in stokes flow using a massively parallel computer

Abstract

Similar works

Full text

Available Versions

UNT Digital Library