Table of Contents Author Guidelines Submit a Manuscript
Abstract and Applied Analysis

Volume 2014, Article ID 563603, 12 pages
Research Article

The FMM-BEM Method for the 3D Particulate Stokes Flow

1Technical and Vocational Training Corporation, Riyadh 11523, Saudi Arabia

2Laboratory of Engineering Mathematics, Polytechnic School of Tunisia, 2078 La Marsa, Tunisia

3Department of Mathematics, College of Sciences, King Saud University, Riyadh 11451, Saudi Arabia

Received 8 April 2014; Accepted 13 May 2014; Published 3 June 2014

Academic Editor: Hassan Eltayeb

Copyright © 2014 Hassib Selmi et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


This work introduces new functions based on the spherical harmonics and the solid harmonics which have been used to construct a multipole development for the 3D Stokes problem in order to reduce the operations costs in the BEM method. We show that the major properties of those functions are inherited from the solid harmonics. The contribution of this paper is the introduction of new formulas that serve to calculate the multipole moments and the transfer functions that are necessary for the schemes of order . Moreover, new translation formulas are introduced to obtain an scheme. The error truncation of the resulting scheme is discussed. In comparison to the BEM that attains a limit storage at , we present here a method based on FMM-BEM that attains a storage at a limit of . The implementation of the method achieves a high accuracy level at a reasonable cost.

1. Introduction

The studies of complex fluids, such as particulate suspensions, emulsions, and sedimentation problem, remain a great challenge in spite of their omnipresence in many physical, chemical, or biomedical processes and industrial applications. This problem represents physically the creeping flow around several solid particles moving in a viscous fluid. Due to the size of the particles, this kind of problems is modelled by the linear Stokes equation. This linearity allows the application of the surface integral equation method or the boundary element method (BEM) which may be more efficient than the boundary value methods such as the finite element method (FEM) or the finite difference method (FDM). Indeed, the advantage of the BEM is that only the surface is discretized, generating fewer elements than a volume discretization and making the system matrix much smaller. However, the integral equation method is a global approach; every point is affected by the points of the entire system, giving a dense matrix. Since there are efficient algorithms, that is, GMRES, for solving the sparse matrix generated by FEM and FDM, it was believed that BEM was more expensive to use than FEM or FDM. The development of the fast multipole method (FMM) with the aim of accelerating the product matrix vector and reducing the complexity of such operations from to or even to , and the memory requirements to or [13], brought new life to the classical BEM. Indeed such method approximates the effect of the far field points of the entire system. Like in FEM, only the terms in a certain band of the matrix are stored in an exact way and all the others are approximated and compressed.

We introduce in this paper new multipole moments and transfer functions for scheme and new formulas for the translations from multipole to multipole, multipole to local, and local to local that are necessary to construct the scheme for solving the 3D particulate Stokes problem. The multipole and local expansions can be useful to develop an efficient method based on the FMM to determine a real time solution of the dynamics of rigid bodies in a 3D infinite or semi-infinite viscous fluid. This solution is obtained by determining the stresses exerted by the fluid on each body and then the velocities of the bodies. These stresses are represented by Stokeslets distributed over the surfaces of the bodies. We begin in the next section by presenting a brief overview of the literature. In Section 3, we present the principle of the fast multipole method. Section 4 presents the two representations used in FMM formulas for the Laplace problem: the first is based on spherical harmonic [2] and the second on the solid harmonic [4]. We present also in this section the new functions that will be used to represent differently the Green function of the Laplace equation. For the sake of conciseness, clarity, and complexity of the algorithm, we choose to work with this new representation with the aim of developing the FMM formulas that permit solving the Stokes problem in an scheme in the first step and in scheme in the final step. In Section 5, we present a detailed description of the FMM formulas founded to be applied to the Green function of the Stokes equation. We also present a discrete algorithm as well as the study of the error bound. Finally, in Section 6, we present, discuss, and interpret the numerical results obtained from a cluster.

2. Overview

FMM reduces the cost of the product matrix vector from to or even by approximating the calculation by handling numerical series. It uses a hierarchical subdivision of space into panels or clusters of sources calculating the multipolar moment that are used in the evaluation of the far field expansion, and it reduces the space of storage from for the dense matrix of size to or even to . From the algorithmic point of view, the FMM works down and up the TREE constructed by the recursive subdivision of the initial domain. FMM was initially introduced by Rokhlin [5] as a fast solution method for integral equations for two-dimensional Laplace equation. Barnes and Hut [6] developed an algorithm. The FMM was then developed by Greengard and Rokhlin [13], as a fast evaluation method for the pairwise force calculation in multibody problems with Coulombic potential. FMM has been applied to problems in various fields such as boundary integral equation method (BIEM) and molecular dynamics (MD). Rokhlin [7] and Coifman et al. [8] used the FMM to solve the Helmholtz equation in 2D. Greengard et al. [9] and Mammoli and Ingber [10, 11] applied the FMM to the Stokes equations in fluid mechanic and Greengard and Helsing [12] for the problem of elasticity. Nishimura et al. [13, 14] studied this topic in order to apply fast multipole boundary integral equation method (FM-BIEM) to practical problems in fracture mechanics and earthquake engineering. Gumerov and Duraiswami used the FMM to solve the biharmonic equation [15]. In this paper, we discuss the applications of the FMM-BIEM to the fundamental boundary value problems in three dimensions given by Stokes equations and show its efficiencies.

3. The FMM Principles for the 3D Stokes Problem

The FMM method replaces the classical product matrix vector by an operation called the multipole product that allows obtaining, by a very fast algorithm, a good approximation of the solution. The complexity is of or in place of for the direct calculus. The basic idea of the multipole method is the separation of the variables in the Green function. This Green function is rewritten differently in a new expression that determines the form of the moments and the transfer functions. To develop the FMM method to the kind of problems of our interest, we will consider the integral equation of the problem to be solved for the unknown values of the intensities of Stokeslets (see Pozrikidis [16]) located on the surfaces of the particles and will consider that this equation has been already discretized by the boundary element method. The problem is then supposed to be reduced to a linear system, of order , of the unknown values of stress at the collocation points of the particles surfaces. is related to the number of the particles by , where is the number of the triangles per particle. The linear system is solved by an iterative method. At each iteration, we start from a great imaginary cube containing the whole collocations points. This cube constitutes level zero of the subdivision. The subdivision is obtained from the one by subdividing each cube into 8 small cubes (cf. Figure 1).

Figure 1: Recursive subdivision of the initial cube.

The level of the constructed tree is chosen after fixing a desired accuracy of computation . The set of collocation points is then distributed into several subdomains. For each one a moment centered at the middle of the subdomain is calculated as well as the transfer function for the far field interaction. To evaluate the integral equation at a given point using the scheme, we have to work in a tree given by the recursive subdivision of the cubic initial domain and evaluate the far field interaction when necessary. The near interactions are evaluated at the finest level. If we have to use the , we must translate the moment of each box at each level to the far field box and convert the information of the far field subdomain at a local moment centered at the middle of the box. Beginning at the coarsest level, these local moments are shifted to the children level (cf. Figure 2).

Figure 2: Transmission of the far field expansion to the boxes in its interaction list and translation to the children.

To evaluate the integral at a given point , we go directly to the box containing this point at the finest level and calculate the far field expansion due to all points outside the neighbors boxes. The near interactions are evaluated directly in a similar manner as in scheme.

4. The Harmonic Functions and Laplace Equation

It is well known that the fundamental solution of the Laplace equation referred to as the Green function of the Laplacian operator can be written using the spherical harmonic [2]. The potential field due to a source point located at and calculated at a point is given by where This potential can be written differently using the solid harmonics and [4]: In this paper, we will use new functions that we define from spherical harmonic as follows: where From (3) and (5), we obtain the following expression for the functions introduced above: Here, the special functions are the associated Legendre functions, defined in [17] as where denotes the Legendre polynomial of degree , defined by Rodrigues formula [17].

The first function coincides with solid harmonic for , and the second one coincides with solid harmonic for .

For the negative values of , we have the following propriety:

4.1. Representation of the Laplacian Kernel

Like in (1) we have the following new representation of the Green function of the Laplace equation: Let us suppose that the point is located inside a sphere of radius with center and is outside the sphere . If we approach by a finite sum obtained from (10) at the order , we obtain an error bound which is the same error obtained by truncating the sums in (1) and (4): The major advantage of this new representation is that we only have to work with positive index , which decreases the complexity of the algorithms in a remarkable way.

5. The Schemes for 3D Stokes Problem

We will develop here the two numerical schemes of the and the order multipole expansion associated with the problem of several solid arbitrarily shaped particles moving in a creeping flow, which are based on a suitable translation formulas of and .

The fundamental solution or the Green function of the free-space Stokes problem [18] is Except for the multiplicative constant, the Green function may be written as follows: and the integral equation for the velocity on a surface of a particle is where , is the number of particles, is the surface of the th particle, and is th component of the surface stress.

The discrete problem associated with (14) is presented as follows: where is the surface element, is the gravity center of a given surface, and denote, respectively, the number of particle and surface element per particle, and is the stress applied on the surface element. Since we are interested here in the far field calculus, the treatment of the singularity problem will not be considered.

5.1. Scheme

Let us suppose that we have uniform distribution points on a certain surface located inside a sphere of radius with center . From (10) and (13), the velocity at a point outside the sphere distant of from is given by the following equation: where One can notice that, for a given pair of and has 9 components and has 3 components, and therefore the number of multipole moments in this formulation is 12. The derivatives that appear in the relations above are calculated in an exact way using the following formulas.

5.1.1. Derivatives Formulas

According to the relation between Cartesian coordinates and polar coordinates and some chain rules [19], the derivatives of can be obtained as follows: In a similar manner, we can obtain the derivatives formulas of for as follows: where . It should be noted that for all .

5.1.2. Error Bound

By truncating (16) at the order , and basing on the error given by (11), the error bound is given by where .

5.2. Scheme

To develop a scheme of order , we introduce three basic translations adequate to the functions and which will be useful to build the formulas of translations multipole to multipole (M2M), multipole to local (M2L), and local to local (L2L). For the rest of the paper, we will take .

5.2.1. Relations Formulas between and

Under the condition , the function satisfies a relationship of the following form: which allows a pole shift that will permit us to construct the multipolar to multipolar translation.

Moreover, and under the condition , the function can be further expressed in the following form that will be used later to construct the conversion multipolar to local moment: Just like the first relation, under the condition , we have a third equation that satisfies a pole relationship that will be used later to construct the local to local translation:

5.2.2. Translation of a Multipole Expansion

Let us suppose that the surface contains uniform distribution points included inside a sphere of radius centered at a point at a distance from the origin with and . The integral equation is given by the multipole expansion (16): If the point is outside the sphere of radius , then the integral equation takes the following form: where

5.2.3. Conversion of a Multipole Expansion into a Local Expansion

Let a surface containing uniform distribution points included inside a sphere of radius centered at a point at a distance from a given point with and . The corresponding multipole expansion (16) converges for any points inside the sphere of radius centered at . The velocity is given by where

5.2.4. Translation of a Local Expansion

Suppose that can be written as in (29), and suppose that ; then for any point , can be written according to the local moment centered at as follows: where

5.2.5. Error Bound

With a similar reasoning as in the error bound of the scheme and under the same hypothesis of the representation (29), the evaluation of the velocity at an order generates an error bound expressed as follows: where .

6. Numerical Results

We present here some numerical results given by applying the FMM-BEM to the Stokes problem. The method uses several parameters on which depend the speed and the precision of the calculation, that is, the size of the cells, the level number of the grids, the order of truncation of the multipole series expansion, and so forth. The optimal values to obtain a good equilibrium between effectiveness and precision are given for the mono- and multilevel approaches. Theoretical complexities are checked numerically. Problems of size are solved in about two hours by iterations on a PC Pentium 4, 3.2GHZ with 4Go-RAM.

6.1. Configuration Test
6.1.1. Sphere Triangulation

The starting point is to generate a grid of many spheres that we put in a cube that represents the initial domain at the initial level of the tree (cf. Figure 3). To generate the collocation points on the spheres, we used a simple technique, presented by Pozrikidis [16] for the discretization of a sphere. Beginning from an octahedron, we developed a triangulation by subdividing recursively each triangle into four subtriangles until obtaining the desired level of discretisation. At each step, the triangle is projected on the surface of the sphere (cf. Figure 4).

Figure 3: Initial configuration.
Figure 4: Sphere triangulation.
6.1.2. Numerical Parameters

For the numerical consideration, we considered that two boxes are far away if they are separated by a third box. In addition, we took for the scheme and for the scheme. The error given by (20) is thus controlled by , where is the level of the tree on which one approximates the far field interactions, and represents the dispatcher of the initial cube. In the results of the tests presented below, we took and we worked on a tree of depth equal to 5 which give a plug of error equal to . We obtained the same thing for the error given by (31).

6.2. Presentation and Interpretation of the Results

It is worth to be pointed out that the maximum value of the error obtained in the different tests is about . These tests permit us to validate the method and show its precision. We remark from the calculation that the error introduced by the FMM compared to the traditional BEM does not have a practical incidence on the quality of the result.


We present here the results obtained with the scheme. Figure 5 that represents the CPU time versus the number of collocation point at different levels shows that the method FMM in the case of the scheme is useful not only to minimize the memory storage of the matrix due to the compression carrying out but also to enable us to gain in terms of computing time. This method appears to be much better than the direct calculation for the large values of where . Indeed we can see that the effectiveness of this method appears from a certain number evaluated by the results between and . For the method FMM does not have any influence from the computing times point of view and from the efficiency of the FMM method becomes to be visible; see Figures 5 and 8. To be more precise, we consider two boxes, shown, respectively, in Figures 6 and 7, that are obtained from Figure 5 by zooming on the regions that seem to be interesting for the interpretation. For the efficiency of the method starts from level 2 (see Figure 6) and then for the efficiency becomes to be visible from level 3 (see Figure 7). We can also see from Figure 7 that level 4 becomes better than level 2 from a certain number evaluated to be between and . Thus, as a general rule about the potential of the FMM method, we can conclude that the more the size of the problem increases, the more we have to work down in the hierarchical tree of the code.

Figure 5: CPU time versus the number of collocation points at different levels.
Figure 6: Zoom 1 of Figure 5.
Figure 7: Zoom 2 of Figure 5.
Figure 8: CPU time versus the level number for different collocation points .

The interpretation of the results found for the scheme is similar to that obtained for the one. The same remarks given from all figures for the scheme may be reproduced here from Figures 9, 10, 11, and 12 which correspond to the . The comparison of the results of the two schemes and is illustrated in Figures 13, 14, 15, and 16. It should be noted that the scheme becomes more efficient than the scheme after a certain number ; see Figure 16. However, one has to work under a certain depth of the tree; otherwise, the becomes more efficient from the CPU time point of view; see Figure 15.

Figure 9: CPU time versus particles number for different levels.
Figure 10: Zoom 1 of Figure 9.
Figure 11: Zoom 2 of Figure 9.
Figure 12: CPU time versus the level for several values of .
Figure 13: CPU time of the two schemes versus the level for .
Figure 14: CPU time of the direct calculus and the two schemes and versus the collocation points for level = 3.
Figure 15: CPU time of the direct calculus and the two schemes and versus the level for several values of .
Figure 16: CPU time of the direct calculus and the two schemes and versus the level for several values of .

7. Concluding Remarks

We presented in this work the fast multipole method applied to the Stokes problem in a 3D infinite viscous fluid in the presence of several rigid bodies. The formulas obtained for both schemes and , where is the size of the matrix obtained from the discretization of the problem, are implemented numerically. The numerical results obtained are very satisfactory and incite us to use the method in the future in a traditional IBEM code that will allow for an efficient resolution of the linear system by iterative methods such as GMRES. This work makes it possible to investigate in future works more complex and realistic configurations such as the sedimentation of dilute suspensions composed by a great number of particles or the movements of deformable particles. Such configurations could enable us to compare between the numerical results and the experimental ones. For the case of deformable bodies, the integral equation must take into account the double potential layer. The development of an adequate FMM-BEM based on formulae given previously constitutes an interesting perspective for a future area of research.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.


This project was supported by King Saud University, Deanship of Scientific Research, College of Sciences Research Center.


  1. L. Greengard and V. Rokhlin, “A fast algorithm for particle simulations,” Journal of Computational Physics, vol. 73, no. 2, pp. 325–348, 1987. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  2. L. Greengard and V. Rokhlin, “On the effficient implementation of the fast multipole algorithm,” Tech. Rep. 602, Department of Computer Science, Yale University, 1988. View at Google Scholar
  3. L. Greengard and V. Rokhlin, “On the evaluation of electrostatic inter-actions in molecular modeling,” Chemica Scripta, vol. 29, pp. 139–144, 1989. View at Google Scholar
  4. H. Y. Wang and R. LeSar, “An efficient fast-multipole algorithm based on an expansion in the solid harmonics,” Journal of Chemical Physics, vol. 104, no. 11, pp. 4173–4179, 1996. View at Google Scholar · View at Scopus
  5. V. Rokhlin, “Rapid solution of integral equations of classical potential theory,” Journal of Computational Physics, vol. 60, no. 2, pp. 187–212, 1985. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  6. J. Barnes and P. Hut, “A hierarchical O(N log N) force-calculation algorithm,” Nature, vol. 324, no. 6096, pp. 446–449, 1986. View at Publisher · View at Google Scholar · View at Scopus
  7. V. Rokhlin, “Rapid solution of integral equations of scattering theory in two dimensions,” Journal of Computational Physics, vol. 86, no. 2, pp. 414–439, 1990. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  8. R. Coifman, V. Rokhlin, and S. Wandzura, “Fast multiple method for the wave equation: a pedestrian prescription,” IEEE Antennas and Propagation Magazine, vol. 35, no. 3, pp. 7–12, 1993. View at Publisher · View at Google Scholar · View at Scopus
  9. L. Greengard, M. C. Kropinski, and A. Mayo, “Integral equation methods for Stokes flow and isotropic elasticity in the plane,” Journal of Computational Physics, vol. 125, no. 2, pp. 403–414, 1996. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  10. A. A. Mammoli and M. S. Ingber, “Stokes flow around cyclinders in a bounded two-dimensional domain using multipole-accelerated boundary element methods,” International Journal for Numerical Methods in Engineering, vol. 44, no. 7, pp. 897–917, 1999. View at Google Scholar · View at MathSciNet
  11. A. A. Mammoli and M. S. Ingber, “Parallel multipole BEM simulation of two-dimensional suspension flows,” Engineering Analysis with Boundary Elements, vol. 24, no. 1, pp. 65–73, 2000. View at Publisher · View at Google Scholar · View at Scopus
  12. L. Greengard and J. Helsing, “On the numerical evaluation of elastostatic fields in locally isotropic two-dimensional composites,” Journal of the Mechanics and Physics of Solids, vol. 46, no. 8, pp. 1441–1462, 1998. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  13. N. Nishimura, K. I. Yoshida, and S. Kobayashi, “A fast multipole boundary integral equation method for crack problems in 3D,” Engineering Analysis with Boundary Elements, vol. 23, no. 1, pp. 97–105, 1999. View at Google Scholar · View at Scopus
  14. N. Nishimura, “Fast multipole accelerated boundary integral equation methods,” Applied Mechanics Reviews, vol. 55, no. 4, pp. 299–324, 2002. View at Publisher · View at Google Scholar · View at Scopus
  15. N. A. Gumerov and R. Duraiswami, “Fast multipole method for the biharmonic equation in three dimensions,” Journal of Computational Physics, vol. 215, no. 1, pp. 363–383, 2006. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet
  16. C. Pozrikidis, Numerical Computation in Science and Engineering, Oxford University Press, New York, NY, USA, 1998. View at MathSciNet
  17. M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, Dover, New York, NY, USA, 1965.
  18. O. A. Ladyzhenskaya, The Mathematical Theory of Viscous Incompressible Flow, Gordon and Breach, Science Publishers, New York, NY, USA, 1969. View at MathSciNet
  19. K. Yoshida, Applications of fast multipole method to boundary integral method [Ph.D. thesis], Department of Global Environmental Engineering, Kyoto University, Kyoto, Japan, 2001.