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

and Applied Analysis 3 Level M+ 1 children Level M boxes 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 x, 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 O(NlogN) 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 Y n [2]. The potential field due to a source point q located at y(ρ, α, β) and calculated at a point x(r, θ, φ) is given by Φ(x, y) = 1 󵄨󵄨󵄨󵄨x − y 󵄨󵄨󵄨󵄨 = ∞


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 ( 2 ) to ( log ) or even to (), and the memory requirements to ( log ) or () [1][2][3], 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 ( log ) 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 ( log ) 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.

Overview
FMM reduces the cost of the product matrix vector from ( 2 ) to ( log ) 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 ( 2 ) for the dense matrix of size  to ( log ) 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 ( log ) algorithm.The FMM was then developed by Greengard and Rokhlin [1][2][3], 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.

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 ( log ) or () in place of ( 2 ) 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  = 3    , 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 th subdivision is obtained from the (−1)th one by subdividing each cube into 8 small cubes (cf. Figure 1).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 (log) 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).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 (log) scheme.

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  ≥ 0, and the second one coincides with solid harmonic    for  ≥ 0.
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.

The Schemes for 3D Stokes Problem
We will develop here the two numerical schemes of the  log() 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  = ∪   =1   ,   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.

𝑂(𝑁logN)
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 ,  , (2)   , () has 9 components and  (1) , () 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.

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  ≥ 1 as follows: where It should be noted that (/  )( 0 0 (  → )) = 0 for all  ∈ {1, 2, 3}.16) at the order , and basing on the error given by (11), the error bound is given by

𝑂(𝑁)
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  = inf(,   ).

Relations Formulas between
which allows a pole shift that will permit us to construct the multipolar to multipolar translation.Moreover, and under the condition |   →  0 | < |   →  0 |, 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 |   →  0 | < |   →  0  1 |, we have a third equation that satisfies a pole relationship that will be used later to construct the local to local translation:

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  > (+1) and  > 1.The integral equation is given by the multipole expansion ( 16): , ()] . ( If the point  is outside the sphere  0 of radius  + , then the integral equation takes the following form: , where ,

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  0 with  > ( + 1) and  > 1.The corresponding multipole expansion (16) converges for any points  inside the sphere  0 of radius  centered at  0 .The velocity is given by , where (28)

Translation of a Local Expansion.
Suppose that V() can be written as in (29), and suppose that  1 ∈  0 ; then for any point  ∈  0 , V() can be written according to the local moment centered at  1 as follows: , where , , , ( 1 ) , , , 5.2.5.Error Bound.With a similar reasoning as in the error bound of the (log) scheme and under the same hypothesis of the representation (29), the evaluation of the velocity V() at an order  generates an error bound expressed as follows: where  = ∫  |()|.

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  =  (10 6 ) 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).

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  > 3 for the (log) scheme and  > 3 for the () scheme.The error given by ( 20) is thus controlled by (3/2)2  (1/3) +1 , 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  = 12 and we worked on a tree of depth equal to 5 which give a plug of error equal to (10 6 ).We obtained the same thing for the error given by (31).

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 10 −6 .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.

𝑂(𝑁logN).
We present here the results obtained with the (log) scheme.Figure 5 that represents the CPU time versus the number   of collocation point at different levels (  = 2( + 1)) shows that the method FMM in the case of the scheme (log) 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  = 3(21).Indeed we can see that the effectiveness of this method appears from a certain number   evaluated by the results between 10 4 and 3 × 10 4 .For   ≤ 9202 the method FMM does not have any influence from the computing times point of view and from   = 27702 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   = 27702 the efficiency of the method starts from level 2 (see Figure 6) and then for   = 351918 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 2.75 × 10 5 and 3 × 10 5 .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.

𝑂(𝑁).
The interpretation of the results found for the () scheme is similar to that obtained for the (log) one.The same remarks given from all figures for the (log()) scheme may be reproduced here from Figures 9, 10

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 ( log ) 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.

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

Figure 8 :Figure 9 :
Figure 8: CPU time versus the level number for different collocation points   .

Figure 12 :Figure 13 :
Figure 12: CPU time versus the level for several values of .

Figure 14 :
Figure 14: CPU time of the direct calculus and the two schemes (log) and () versus the collocation points   for level = 3.

11 CPUFigure 15 :Figure 16 :
Figure 15: CPU time of the direct calculus and the two schemes (log) and () versus the level for several values of   .