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

Abstract

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 𝑂(𝑁log𝑁). 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 𝑂(10⁴), we present here a method based on FMM-BEM that attains a storage at a limit of 𝑂(10⁶). 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 𝑂(𝑁log𝑁) or even to 𝑂(𝑁), and the memory requirements to 𝑂(𝑁log𝑁) or 𝑂(𝑁) [1–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.

2. Overview

FMM reduces the cost of the product matrix vector from 𝑂(𝑁²) 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 𝑂(𝑁²) 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–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.

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

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

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 𝑂(𝑁log𝑁) 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 Φ=1𝑥,𝑦=𝑥−𝑦^∞𝑛_{𝑛=0𝑚=−𝑛}𝑀^𝑚_𝑛𝑟^𝑛+1𝑌^𝑚_𝑛𝜃,𝜙, (1) where 𝑀^𝑚_𝑛=𝑞𝜌^𝑛𝑌_𝑛^−𝑚𝛼,𝛽,𝑌^𝑚_𝑛𝜃,𝜙=(𝑛−|𝑚|)!(|𝑚!𝑃𝑛+|)^𝑚𝑛(cos(𝜃))𝑒^𝑖𝑚𝜙. (2)(3) This potential can be written differently using the solid harmonics 𝑅^𝑚_𝑛 and 𝑆^𝑚_𝑛 [4]: Φ𝑥,𝑦=^∞𝑛_{𝑛=0𝑚=−𝑛}𝑅^𝑚_𝑛(𝑜^→𝑦)𝑆^𝑚_𝑛^→=𝑜𝑥_{∞𝑛𝑛=0𝑚=−𝑛}𝑅𝑜^𝑚_𝑛^→𝑦𝑆𝑜^𝑚_𝑛^→𝑥. (4) In this paper, we will use new functions that we define from spherical harmonic as follows: 𝑆𝑦^𝑚_𝑛=𝛼^𝑚_𝑛1𝑟^𝑛+1𝑌^𝑚_𝑛𝜃,𝜙,𝑅𝑦^𝑚_𝑛=1𝛼^𝑚_𝑛𝜌^𝑛𝑌^𝑚_𝑛𝛼,𝛽, (5) where 𝛼^𝑚_𝑛=(𝑛+|𝑚|)!(𝑛−|𝑚|)!. (6) From (3) and (5), we obtain the following expression for the functions introduced above: 𝑆𝑦𝑜^𝑚_𝑛^→𝑥=(𝑛−|𝑚|)!𝑃_𝑛^|𝑚|(cos𝜃)𝑒^𝑖𝑚𝜃1𝑟^𝑛+1,𝑅𝑦^𝑚_𝑛^→=1𝑜𝑦(|𝑚!𝑃𝑛+|)^𝑛|𝑚|(cos𝛼)𝑒^𝑖𝑚𝛽𝜌^𝑛. (7) Here, the special functions 𝑃^𝑚_𝑛 are the associated Legendre functions, defined in [17] as 𝑃^𝑚_𝑛(𝑥)=(−1)^𝑚1−𝑥^2𝑚/2𝑑^𝑚𝑑𝑥^𝑚𝑃_𝑛(𝑥), (8) 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: 𝑆𝑜_𝑛^−𝑚→𝑥=(−1)^𝑚𝑆^𝑚_𝑛(^→𝑅𝑜𝑥);_{𝑛−𝑚→}𝑜𝑥=(−1)^𝑚𝑅𝑜^𝑚_𝑛^→𝑥,𝑆𝑦𝑜_𝑛^−𝑚→𝑥=𝑆𝑦^𝑚_𝑛(𝑜^→𝑥);𝑅𝑦𝑜_𝑛^−𝑚→𝑥=𝑅𝑦^𝑚_𝑛(𝑜^→𝑥)𝑆𝑦𝑜_𝑛^−𝑚→𝑥=(−1)^𝑚𝑆𝑜_𝑛^−𝑚→𝑥;𝑅𝑦𝑜_𝑛^−𝑚→𝑦=(−1)^𝑚𝑅𝑜_𝑛^−𝑚→𝑦;𝑚>0. (9)

4.1. Representation of the Laplacian Kernel

Like in (1) we have the following new representation of the Green function of the Laplace equation: Φ𝑥,𝑦=^∞𝑛_{𝑛=0𝑚=−𝑛}𝑅𝑦𝑜^𝑚_𝑛^→𝑦𝑆𝑦^𝑚_𝑛^→=𝑜𝑥_{∞𝑛𝑛=0𝑚=−𝑛}𝑅𝑦𝑜^𝑚_𝑛^→𝑦𝑆𝑦𝑜^𝑚_𝑛^→𝑥. (10) 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): 𝐸_𝑝=1𝑎𝑟−𝑎𝑟^𝑝+1. (11) 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 𝑁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 𝐺^𝑗_𝑘1𝑥,𝑦=−𝛿8𝜋𝜈^𝑗𝑘+𝑥𝑥−𝑦^𝑗−𝑦_𝑗𝑥_𝑘−𝑦_𝑘𝑥−𝑦³. (12) Except for the multiplicative constant, the Green function may be written as follows: 𝐺^𝑗_𝑘𝑥,𝑦=𝛿^𝑗_𝑘1−𝑥−𝑦^→𝑜𝑥_𝑘𝜕𝜕𝑥_𝑗1+𝑥−𝑦_→𝑜𝑦_𝑘𝜕𝜕𝑥_𝑗1,𝑥−𝑦 (13) and the integral equation for the velocity on a surface 𝑆_𝑙 of a particle 𝑙 is 𝑣_𝑘(𝑥)=_𝑆𝐺^𝑗_𝑘𝑥,𝑦𝑞_𝑗𝑦𝑑𝑠;𝑥∈𝑆_𝑙, (14) 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: 𝑣_𝑘𝑥^{𝑙′,𝑠′}_𝐺=^𝑁𝑝𝑁_𝑙=1^𝑡_{𝑠=1𝑇𝑙𝑠}𝐺^𝑗_𝑘𝑥^{𝑙′,𝑠′}_𝐺,𝑦𝑑𝑠𝑦𝑞_𝑗^𝑠,𝑙, (15) 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. 𝑂(𝑁log𝑁) 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: 𝑣_𝑘(𝑥)=_𝑆𝐺^𝑗_𝑘𝑥,𝑦𝑞_𝑗𝑦=𝑑𝑠_{∞𝑛𝑛=0𝑚=−𝑛}𝐹𝑦𝑜^𝑗𝑘_𝑛𝑚^→𝑥𝑀^𝑚(1)_𝑛,𝑗(𝑜)+𝐻𝑦𝑜^𝑗_𝑛𝑚^→𝑥𝑀^𝑚,𝑘(2)_𝑛,𝑗(𝑜), (16) where 𝑀^𝑚,(1)_𝑛,𝑗(𝑜)=_𝑠𝑅𝑦𝑜^𝑚_𝑛^→𝑦𝑞_𝑗𝑦𝑀𝑑𝑠,_{𝑚,𝑘(2)𝑛,𝑗}(𝑜)=𝑜_𝑠^→𝑦_𝑘𝑅𝑦𝑜^𝑚_𝑛^→𝑦𝑞_𝑗𝑦𝑑𝑠,𝐹𝑦𝑜^𝑗𝑘_𝑛𝑚^→𝑥=𝛿^𝑗_𝑘𝑆𝑦𝑜^𝑚_𝑛^→𝑥−𝑜^→𝑥_𝑘𝜕𝜕𝑥_𝑗𝑆𝑦𝑜^𝑚_𝑛^→𝑥,𝐻𝑦^𝑗→_𝑛,𝑚𝜕𝑜𝑥=𝜕𝑥_𝑗𝑆𝑦𝑜^𝑚_𝑛^→𝑥. (17) 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.

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: 𝜕𝜕𝑥₁𝑆𝑦^𝑚_𝑛^→=1𝑜𝑥2𝑆𝑦𝑜^𝑚−1_𝑛+1^→𝑥−𝑆𝑦𝑜^𝑚+1_𝑛+1^→𝑥,if𝑚>0;−Real𝑆𝑦𝑜¹_𝑛+1^→𝑥,if𝑚=0;𝜕𝜕𝑥₁𝑆𝑦_𝑛^−𝑚→𝜕𝑜𝑥,if𝑚<0.𝜕𝑥₂𝑆𝑦^𝑚_𝑛^→=𝑖𝑜𝑥2𝑆𝑦𝑜^𝑚−1_𝑛+1^→𝑥+𝑆𝑦𝑜^𝑚+1_𝑛+1^→𝑥,if𝑚>0;𝑖Imag𝑆𝑦𝑜¹_𝑛+1^→𝑥,if𝑚=0;𝜕𝜕𝑥₂𝑆𝑦_𝑛^−𝑚→𝜕𝑜𝑥,if𝑚<0.𝜕𝑥₃𝑆𝑦𝑜^𝑚_𝑛^→𝑥=−𝑆𝑦𝑜^𝑚_𝑛+1^→𝑥,if𝑚>0;𝜕𝜕𝑥₃𝑆𝑦𝑜_𝑛^−𝑚→𝑥,else. (18) In a similar manner, we can obtain the derivatives formulas of 𝑅𝑦^𝑚_𝑛 for 𝑛≥1 as follows: 𝜕𝜕𝑥₁𝑅𝑦^𝑚_𝑛^→=1𝑜𝑥2𝑅𝑦𝑜^𝑚−1_𝑛+1^→𝑥−𝑅𝑦𝑜^𝑚+1_𝑛+1^→𝑥,if𝑚>0;−Real𝑅𝑦𝑜¹_𝑛+1^→𝑥,if𝑚=0;𝜕𝜕𝑥₁𝑅𝑦_𝑛^−𝑚→𝜕𝑜𝑥,if𝑚<0,𝜕𝑥₂𝑅𝑦^𝑚_𝑛^→=𝑖𝑜𝑥2𝑅𝑦𝑜^𝑚−1_𝑛+1^→𝑥+𝑅𝑦𝑜^𝑚+1_𝑛+1^→𝑥,if𝑚>0;𝑖Imag𝑅𝑦𝑜¹_𝑛+1^→𝑥,if𝑚=0;𝜕𝜕𝑥₂𝑅𝑦_𝑛^−𝑚→𝜕𝑜𝑥,if𝑚<0,𝜕𝑥₃𝑅𝑦𝑜^𝑚_𝑛^→𝑥=−𝑅𝑦𝑜^𝑚_𝑛+1^→𝑥,if𝑚>0;𝜕𝜕𝑥₃𝑅𝑦𝑜_𝑛^−𝑚→𝑥,else, (19) where 𝑖²=−1. It should be noted that (𝜕/𝜕𝑥_𝑗)(𝑅𝑦⁰₀(𝑜^→𝑥))=0 for all 𝑗∈{1,2,3}.

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 𝐸^′_𝑝𝛿=𝐶^𝑘𝑗+1𝑟−𝑎+𝑎𝑟−𝑎𝑟²𝑎−𝑟𝑎𝑟^𝑝+1, (20) 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 𝑘=inf(𝑚,𝑛^′).

5.2.1. Relations Formulas between 𝑆𝑦^𝑚_𝑛 and 𝑅𝑦^𝑚_𝑛

Under the condition |𝑜^→𝑦|<|𝑜^→𝑜^′|, the function 𝑅𝑦^𝑚_𝑛 satisfies a relationship of the following form: 𝑅𝑦𝑜^𝑚_𝑛^→′𝑦=^𝑛_𝑛′=0^𝑘_𝑚′=0𝑅𝑦𝑜^𝑚′_𝑛′^→′𝑜𝑅𝑦^{𝑚−𝑚′}_{𝑛−𝑛′}^→+𝑜𝑦_{𝑛′𝑚′=𝑘+1}(−1)^{𝑚′−𝑚}𝑅𝑦𝑜^𝑚′_𝑛′^→′𝑜𝑅𝑦^{𝑚′−𝑚}_{𝑛−𝑛′}(^→+𝑜𝑦)_{𝑛′𝑚′=0}(−1)^𝑚′𝑅𝑦𝑜^𝑚′_𝑛′^→′𝑜𝑅𝑦𝑜^𝑚+𝑚′_{𝑛−𝑛′}^→𝑦 (21) 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: 𝑆𝑦^𝑚_𝑛^→=𝑜𝑥_{∞𝑛′=0𝑛′𝑚′=0}(−1)^𝑛′𝑅𝑦𝑥^𝑚′_𝑛′^→₀𝑥𝑆𝑦𝑜^𝑚+𝑚′_𝑛+𝑛′^→𝑥₀+^𝑘_𝑚′=0(−1)^{𝑛′+𝑚′}𝑅𝑦𝑥^𝑚′_𝑛′^→₀𝑥𝑆𝑦𝑜^{𝑚−𝑚′}_𝑛+𝑛′^→𝑥₀+^𝑛′_{𝑚′=𝑘+1}(−1)^{𝑛+2𝑚′−𝑚}𝑅𝑦𝑥^𝑚′_𝑛′^→₀𝑥𝑆𝑦𝑜^{𝑚′−𝑚}_𝑛+𝑛′^→𝑥₀. (22) 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: 𝑅𝑦𝑥^𝑚_𝑛^→₀𝑥=^∞_𝑛′=𝑛^𝑘_𝑚′=0𝑅𝑦𝑥^𝑚′_𝑛′^→₁𝑥𝑅𝑦𝑥^{𝑚−𝑚′}_{𝑛−𝑛′}^→₀𝑥₁+^𝑛′_{𝑚′=𝑘+1}(−1)^{𝑚′−𝑚}𝑅𝑦𝑥^𝑚′_𝑛′^→₁𝑥𝑅𝑦𝑥^{𝑚′−𝑚}_{𝑛−𝑛′}^→₀𝑥₁+^𝑛′_𝑚′=0(−1)^𝑚′𝑅𝑦𝑥^𝑚′_𝑛′^→₁𝑥𝑅𝑦𝑥^𝑚+𝑚′_{𝑛−𝑛′}^→₀𝑥₁. (23)

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 𝜌>(𝑐+1)𝑎 and 𝑐>1. The integral equation is given by the multipole expansion (16): 𝑣_𝑘(𝑥)=^∞𝑛_{𝑛=0𝑚=−𝑛}𝐹𝑦𝑜^𝑗𝑘_𝑛𝑚^→𝑥𝑀^𝑚(1)_𝑛,𝑗(𝑜)+𝐻𝑦𝑜^𝑗_𝑛𝑚^→𝑥𝑀^𝑚,𝑘(2)_𝑗,𝑚(𝑜). (24) If the point 𝑥 is outside the sphere 𝐷₀ of radius 𝑎+𝜌, then the integral equation takes the following form: 𝑣_𝑘(𝑥)=^∞𝑛_{𝑛=0𝑚=−𝑛}𝐹𝑦𝑜^𝑗𝑘_𝑛𝑚^→′𝑥𝑀^𝑚(1)_𝑛,𝑗𝑜^′+𝐻𝑦𝑜^𝑗_𝑛𝑚^→′𝑥𝑀^𝑚,𝑘(2)_𝑗,𝑚𝑜^′, (25) where 𝑀^𝑚(1)_𝑛,𝑗𝑜^′=^𝑛_𝑛′=0^𝑘_𝑚′=0𝑅𝑦𝑜^𝑚′_𝑛′^→′𝑜𝑀^{𝑚−𝑚′(1)}_{𝑛−𝑛′,𝑗}(𝑜)+^𝑛′_{𝑚′=𝑘+1}(−1)^{𝑚′−𝑚}𝑅𝑦𝑜^𝑚′_𝑛′^→′𝑜𝑀^{𝑚′−𝑚(1)}_{𝑛−𝑛′,𝑗}(𝑜)+^𝑛′_𝑚′=1(−1)^𝑚′𝑅𝑦𝑜^𝑚′_𝑛′^→′𝑜𝑀^{𝑚+𝑚′(1)}_{𝑛−𝑛′,𝑗}(𝑜),𝑀^𝑚,𝑘(2)_𝑛,𝑗𝑜^′=𝑜^→′𝑜_𝑘𝑀^𝑚(1)_𝑛,𝑗𝑜^′+^𝑛_𝑛′=0^𝑘_𝑚′=0𝑅𝑦𝑜^𝑚′_𝑛′^→′𝑜𝑀^{𝑚−𝑚′,𝑘(2)}_{𝑛−𝑛′,𝑗}(𝑜)+^𝑛′_{𝑚′=𝑘+1}(−1)^{𝑚′−𝑚}𝑅𝑦𝑜^𝑚′_𝑛′^→′𝑜𝑀^{𝑚′−𝑚,𝑘(2)}_{𝑛−𝑛′,𝑗}(𝑜)+^𝑛′_𝑚′=1(−1)^𝑚′𝑅𝑦𝑜^𝑚′_𝑛′^→′𝑜𝑀^{𝑚+𝑚′,𝑘(2)}_{𝑛−𝑛′,𝑗}(𝑜). (26)

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 𝜌>(𝑐+1)𝑎 and 𝑐>1. The corresponding multipole expansion (16) converges for any points 𝑥 inside the sphere 𝐷₀ of radius 𝑎 centered at 𝑥₀. The velocity is given by 𝑣_𝑘(𝑥)=_𝑆𝐺^𝑗_𝑘𝑥,𝑦𝑞_𝑗𝑦=𝑑𝑠_{∞𝑛𝑛=0𝑚=−𝑛}𝐹𝑦𝑅𝑥^𝑗𝑘_𝑛𝑚^→₀𝑥𝐿^𝑚(1)_𝑛,𝑗𝑥₀+𝐻𝑦𝑅𝑥^𝑗_𝑛𝑚^→₀𝑥𝐿^𝑚,𝑘(2)_𝑛,𝑗𝑥₀, (27) where 𝐿^𝑚(1)_𝑛,𝑗𝑥₀=^∞_𝑛′=0^𝑛′_𝑚′=0(−1)^𝑛𝑀^𝑚′_𝑛⁽¹⁾_′,𝑗(𝑜)𝑆𝑦𝑜^𝑚+𝑚′_𝑛+𝑛′^→𝑥₀+^𝑘_𝑚′=1(−1)^𝑛+𝑚′𝑀^𝑚′_𝑛⁽¹⁾_′,𝑗(𝑜)𝑆𝑦𝑜^{𝑚−𝑚′}_𝑛+𝑛′^→𝑥₀+^𝑛′_{𝑚′=𝑘+1}(−1)^𝑛−𝑚𝑀^𝑚′_𝑛⁽¹⁾_′,𝑗(𝑜)𝑆𝑦𝑜^{𝑚′−𝑚}_𝑛+𝑛′^→𝑥₀,𝐿^𝑚,𝑘(2)_𝑛,𝑗𝑥₀=−𝑜^→𝑥_0𝑘𝐿^𝑚(1)_𝑛,𝑗𝑥₀+^∞_𝑛′=0^𝑛′_𝑚′=0(−1)^𝑛𝑀^𝑚′_𝑛^,𝑘(2)_′,𝑗(𝑜)𝑆𝑦𝑜^𝑚+𝑚′_𝑛+𝑛′^→𝑥₀+^𝑘_𝑚′=1(−1)^𝑛+𝑚′𝑀^𝑚′_𝑛^,𝑘(2)_′,𝑗(𝑜)𝑆𝑦𝑜^{𝑚−𝑚′}_𝑛+𝑛′^→𝑥₀+^𝑛′_{𝑚′=𝑘+1}(−1)^𝑛−𝑚𝑀^𝑚′_𝑛^,𝑘(2)_′,𝑗(𝑜)𝑆𝑦𝑜^{𝑚′−𝑚}_𝑛+𝑛′^→𝑥₀,𝐹_𝑦𝑅𝑥^𝑗𝑘_𝑛𝑚^→₀𝑥=𝛿^𝑗_𝑘𝑅𝑦𝑥^𝑚_𝑛^→₀𝑥−𝑥^→₀𝑥𝜕𝜕𝑥_𝑗𝑅𝑦𝑥^𝑚_𝑛^→₀𝑥,𝐻𝑦𝑅𝑥^𝑗_𝑛𝑚^→₀𝑥=𝜕𝜕𝑥_𝑗𝑅𝑦𝑥^𝑚_𝑛^→₀𝑥. (28)

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: 𝑣_𝑘(𝑥)=_𝑆𝐺^𝑗_𝑘𝑥,𝑦𝑞_𝑗𝑦=𝑑𝑠_{∞𝑛𝑛=0𝑚=−𝑛}𝐹𝑦𝑅𝑥^𝑗𝑘_𝑛𝑚^→₁𝑥𝐿^𝑚(1)_𝑛,𝑗𝑥₁+𝐻𝑦𝑅𝑥^𝑗_𝑛𝑚^→₁𝑥𝐿^𝑚,𝑘(2)_𝑛,𝑗𝑥₁, (29) where 𝐿^𝑚,𝑘(1)_𝑛,𝑗𝑥₁=^∞_𝑛′=0^𝑘_𝑚′=0(−1)^{𝑚−𝑚′}𝐿^𝑚′_𝑛⁽¹⁾_′,𝑗𝑥₀𝑅𝑦𝑥^{𝑚−𝑚′}_{𝑛−𝑛′}^→₀𝑥₁+^𝑛′_{𝑚′=𝑘+1}𝐿^𝑚′_𝑛⁽¹⁾_′,𝑗𝑥₀𝑅𝑦𝑥^{𝑚′−𝑚}_{𝑛−𝑛′}^→₀𝑥₁+^𝑛′_𝑚′=1(−1)^𝑚𝐿^𝑚′_𝑛⁽¹⁾_′,𝑗𝑥₀𝑅𝑦𝑥^𝑚+𝑚′_{𝑛−𝑛′}^→₀𝑥₁,𝐿^𝑚,𝑘(2)_𝑛,𝑗𝑥₁=−𝑥^→₀𝑥_1𝑘𝐿^𝑚(1)_𝑛,𝑗𝑥₁+^∞_𝑛′=0^𝑘_𝑚′=0(−1)^{𝑚−𝑚′}𝐿^𝑚′_𝑛^,𝑘(2)_′,𝑗𝑥₀𝑅𝑦𝑥^{𝑚−𝑚′}_{𝑛−𝑛′}^→₀𝑥₁+^𝑛′_{𝑚′=𝑘+1}𝐿^𝑚′_𝑛^,𝑘(2)_′,𝑗𝑥₀𝑅𝑦𝑥^{𝑚′−𝑚}_{𝑛−𝑛′}^→₀𝑥₁+^𝑛′_𝑚′=1(−1)^𝑚𝐿^𝑚′_𝑛^,𝑘(2)_′,𝑗𝑥₀𝑅𝑦𝑥^𝑚+𝑚′_{𝑛−𝑛′}^→₀𝑥₁. (30)

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 𝑣(𝑥) at an order 𝑝 generates an error bound expressed as follows: 𝐸_𝑝^′′𝛿=𝐶^𝑘𝑗+1𝜌−𝑟+𝑟𝜌−𝑟𝜌²𝑟−𝜌𝑟𝜌^𝑝+1, (31) 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 𝑁=𝑂(10⁶) 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).

(a)

(b)

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 𝑟>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⁶). 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 10⁻⁶. 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.

6.2.1. 𝑂(𝑁log𝑁)

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𝑁𝑝(𝑁2𝑐1). Indeed we can see that the effectiveness of this method appears from a certain number 𝑁_𝑐 evaluated by the results between 10⁴ and 3×10⁴. 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⁵ and 3×10⁵. 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.

(a)

(b)

Figure 8 

CPU time versus the level number for different collocation points 𝑁𝑐.

6.2.2. 𝑂(𝑁)

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, 11, and 12 which correspond to the 𝑂(𝑁). The comparison of the results of the two schemes 𝑂(𝑁log𝑁) and 𝑂(𝑁) is illustrated in Figures 13, 14, 15, and 16. It should be noted that the 𝑂(𝑁) scheme becomes more efficient than the 𝑂(𝑁log𝑁) scheme after a certain number 𝑁; see Figure 16. However, one has to work under a certain depth of the tree; otherwise, the 𝑂(𝑁log𝑁) 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.

(a)

(b)

Figure 12 

CPU time versus the level for several values of 𝑁.

Figure 13 

CPU time of the two schemes versus the level for 𝑁𝑐=351918.

Figure 14 

CPU time of the direct calculus and the two schemes 𝑂(𝑁log𝑁) and 𝑂(𝑁) versus the collocation points 𝑁𝑐 for level = 3.

Figure 15 

CPU time of the direct calculus and the two schemes 𝑂(𝑁log𝑁) and 𝑂(𝑁) versus the level for several values of 𝑁𝑐.

Figure 16 

CPU time of the direct calculus and the two schemes 𝑂(𝑁log𝑁) 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 𝑂(𝑁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.

Conflict of Interests

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

Acknowledgment

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