ACCELERATION OF MODAL ANALYSIS BY FMM BASED ON DRBEM

VIB24

August 12-15, 2012, Chicago, Illinois, USA

DETC2012-70596

DRAFT: ACCELERATION OF MODAL ANALYSIS BY FMM BASED ON DRBEM

Yixiong Wei Qifu Wang* Yingjun Wang Yunbao Huang Linchi Zhang

National CAD Support Software Engineering Research Center

Huazhong University of Science and Technology

Wuhan, Hubei, China 430074

employed in various fields, such as acoustics, electromagnetics, ABSTRACT

This paper proposes a novel algorithm to accelerate the elastostatics, etc. process of modal analysis in 3D elastodynamic problems in Nevertheless, the BEM forms non-symmetric and full BEM (boundary element method) with high accuracy. Because matrix, which results in huge workload during computation in of low efficiency and high cost, conventional BEM is rarely elastodynamics. While introducing the time dimension, the used for solving 3D elastodynamics problems in engineering integral equation will need another discretization for time. problems. With applying the DRBEM (dual reciprocity Suppose that one solves an elastodynamics boundary value boundary element method) to form new integral equations of problem for a cube using 𝑂(𝑛) nodes and 𝑂(𝑚) time steps. 3D elastodynamics problems to reduce time complexity by 𝑁 is the scale of coefficient matrix. In that case, the using reciprocity method twice, we introduce modified FMM computational complexity of BEM is 𝑂(𝑛4)∗𝑂(𝑚4) since (fast multipole method) to simplify the computation process 𝑁=𝑛2∗𝑚2, while that for FDM (Finite difference method) or and improve the efficiency from 𝑂(𝑛2)to 𝑂(𝑛)in matrix FEM is 𝑂(𝑛3)∗𝑂(𝑚) since 𝑁=𝑛3. Because of the non-multiplication. The main features in this method are: (1) banded, non-symmetric coefficient matrices and time Position Location (PL) algorithm is used to eliminate one layer dimension, the computational complexity will be bigger as the of nested loops in conventional FMM, and which achieve a scale of time reference problems improves. Many efforts have good performance in efficiency; (2) time dimension been made to accelerate this process. In Frequency domain, the integrations in the element of matrices are canceled for high papers by Chen et al. [2], Fukui and Inoue K. [3], Fujiwara [4] efficiency; (3) instead of the interaction between points, we in 2D and Yoshida et al. [6], Fujiwara [7] in 3D applied FMM apply point to element interaction method for saving plenty of into elastodynamics to improve efficiency, however the reverse the CPU cost in modified FMM; (4) it does not need to transformation to time domain is high cost but necessary in compute complex dynamic fundamental solutions which are solving engineering problems. In time domain, Takahashi et al. necessary in conventional BEM. In this algorithm, the [5], Lu et al. [8] also applied FMM to accelerate computation corresponding eigenvalue problem is solved by Hessenberg process, but the computation is still with low efficiency when matrix and QR reduction algorithm iteratively. We have tested the fundamental solutions have reference with time. In addition, our method in numerical examples during last section, and have the frequency variable is implicit in the matrices during observed significant optimal results in efficiency and accuracy. analysis, which is very hard to extract in traditional way.

Nardini and Brebbia introduced the Dual Reciprocity 1. INTRODUCTION Boundary Element Method [9], in which the inertial volume

The finite element method has some drawbacks such as integral is transformed into a surface integral with the aid of the numerical cost in 3D elastodynamics [1] and in element reciprocal theorem applied to the second time and after the meshing. Comparing with FEM(Finite element method), the approximation of the displacement field by a finite series boundary element method is more fit for fast pretreatment, fast involving radial basis functions (RBF selected 1+r in this paper, calculation, and self-adaptive structure analysis engineering r is defined as the distance between source and field points). software, due to the advantages such as dimensionality [10] This method combines the advantage of dimensionality reduction, high precision and so on. At present, BEM has been reduction with simply elastostatic foundation solutions. Ahmad

*

Corresponding Author Email: hero471@gmail.com

1 Copyright © 2012 by ASME

and Banerjee [11] proposed a particular integral BEM approach for conversion from domain to surface ones, which is proved to be mathematically equivalent to DRBEM [10]. With respect to modal analysis, the governing equation in this method takes the frequency variable explicitly, also the formulation is similar to that of the finite element with full and non-symmetric matrices.

FMM was introduced by Rokhlin [12] as an 𝑂(𝑁)numerical method for solving an integral equation with 2D Laplace’s equation. This method becomes well noted with the introduction into multi-body problems by Greengard [13][14] for solving potential and particle simulation problems. The influence of FMM to science and engineering was so profound that this method is ranked into the top ten algorithms of the 20th century along with Dantzig’s simplex method for linear programming, Krylov subspace iteration, QR algorithm, FFT, etc [15].

Two aspects limit the application development of DRBEM in elastodynamics. First, different from the banded, sparse and symmetric matrix formed by FEM, plenty of non-symmetric and full matrix multiplications in BEM consumes huge memory space and costs long CPU time; second, during the process of forming coefficient matrix by conventional BEM, time integral computation needs another discretization which leads to huge time cost. The present paper applies DRBEM with static fundamental solution to avoid time discretization computation while forming coefficient matrices. For improving the efficiency of matrix multiplication process, the modified FMM is introduced. Comparing with conventional method, improvements are made in two aspects: (1) calculating the interaction of point to element instead of the interaction between points; (2) using PL (position location) method to omit unnecessary traversal of all cell data.

This paper deals with the acceleration of DRBEM when it is applied into modal analysis. The first section briefly introduces the DRBEM and conventional FMM theories, the second section describes the advantages of our algorithm and introduces the improvements of modified FMM in detail. We have tested our method in numerical examples during last section. This paper lists the efficiency comparison with traditional DRBEM in the CPU time cost, and the accuracy of natural frequency comparison with FEM. We have observed significant optimal results, and the proposed method is concluded to outperform the traditional DRBEM.

2. DRBEM IN ELASTODYNAMICS

A brief review of DRBEM for 3D elastodynamics is presented here for completeness.

It can simplify the process of computation greatly by using static fundamental functions instead of dynamic ones. But the inertial item will appear in the equation such that the domain field needs discretization, and which make the technique lose the attraction of its “boundary only” character. The DRBEM is essentially a generalized way of constructing particular solutions that can be used to represent internal source distribution.

For the homogenous medium, the motion of linearly elastic body of volume Ω and surface 𝑆 is described by Navier-Cauchy partial differential equation as

(𝜆+2𝜇)∇∇∙𝒖(𝒙,𝑡)−𝜇∇×∇×𝒖(𝒙,𝑡)+𝜌𝑏

=𝜌𝒖(𝒙,𝑡) (1) in which 𝒖(𝒙,𝑡) is the displacement vector at point 𝒙 and time 𝑡, and 𝜆 and 𝜇 are Lame elastic constants and 𝜌 is density of body. Assuming zero body forces and initial conditions, as apply traditional BEM one can obtain an integral representation of the Eq. (1) in the form 𝒄(𝑥)𝑢(𝑥,𝑡)=∫,𝑼∗(𝑥,𝜉)𝑝(𝜉,𝑡)

𝛤

−𝑾∗(𝑥,𝜉)𝑢(𝜉,𝑡)-𝑑𝛤(𝜉) (2)

in which 𝑼∗(x,ξ) and 𝑾∗(x,ξ) is the dynamic fundamental solutions as displacement and traction tensor in respect, 𝑝(𝜉,𝑡) is the traction vector at point x and time t.𝑐(𝑥) is the jump tensor and overdots indicate differentiation with respect to time. The fundamental solution is very complicated for computation. Because of the time variables in integral equation, the complexity will be 𝑂(𝑛4)∗𝑂(𝑚4) compare to the 𝑂(𝑛3)∗𝑂(𝑚)in FEM with 𝑂(𝑛) nodes and 𝑂(𝑚) time steps, that is very inefficiency. For simplify the computation process, one can use static functional solutions instead of dynamic ones, such as

𝑐(𝑥)𝒖(𝑥,𝑡)=∫,𝒖∗(𝑥,𝜉)𝑝(𝜉,𝑡)−𝒑∗(𝑥,𝜉)𝑢(𝜉,𝑡)-𝑑𝛤(𝜉)

𝛤

−∫𝒖∗(𝑥,𝜉)𝜌𝑢(𝜉,𝑡)𝑑𝛺(𝜉) (3)

𝛺

in which 𝒖∗(𝑥,𝜉) and 𝒑∗(𝑥,𝜉) is static fundamental solutions respectively. Eq. (3) has obtained advantage of much simpler elastostatic fundamental solution pair and also avoiding the time convolutions in conventional time domain. However, the presence of the inertial volume integral in Eq. (3) indicates the discretization in volume domain is necessary, which would eliminate the biggest advantage of BEM. Nardini and Brebbia [9] transformed this volume integral into a boundary one, thereby creating an all-boundary integral formulation that leading to DRBEM.

The start point of DRBEM is to express the unknown 𝒖(𝒙,𝑡) inside Ω as a series of unknown time dependent coefficients 𝛼𝑖𝑚(𝑡) and known basis function 𝑓𝑚(𝑥)

𝑀

𝑢𝑖(𝑥,𝑡)=∑𝛼𝑖𝑚(𝑡)𝑓𝑚(𝑥),𝑥∈Ω (4)

𝑚=1

in which 𝑀=𝑁+𝐿 with 𝑁 and 𝐿 being the number of boundary and internal collocation points in respect, and L could be zero. Using the reciprocity theorem again one succeeds in transforming Eq. (3) into a boundary integral of the form

2 Copyright © 2012 by ASME

−∫𝒖∗(𝒙,𝜉)𝜌𝒖(𝜉,𝑡)𝑑𝛀(𝜉)

Ω

𝑀

=𝜌∑𝛼̈𝑚𝑛*𝑐𝑖𝑗(𝒙)𝜅𝑗𝑛

𝑚(𝒙)𝑚=1

+∫𝑝∗(𝒙,𝜉)𝜅Γ

𝑖𝑗

𝑗𝑛𝑚

(𝒙)𝑑𝚪(𝜉)−∫𝑢∗(𝒙,𝜉)𝜁𝑗𝑛𝑚(𝒙)𝑑𝚪(𝜉)+ (5) Γ

𝑖𝑗

in which 𝜅𝑗𝑛

𝑚(𝒙) and 𝜁𝑗𝑛𝑚

(𝒙) is the particular solution for displacement and traction field, respectively. Then discretization of the boundary 𝚪 into number of triangle elements with total number of 𝑁 nodes, Eq. (3) could form the matrix equation

,𝑴-*𝒖++,𝑯-*𝒖+=,𝑮-*𝒑+ (6)

where

,𝑴-=𝜌(,𝑮-,𝑷-−,𝑯-,𝑾-),𝑭-−1 (7)

in which ,𝑯- and ,𝑮- are the integral coefficient matrices, and ,𝑷- and ,𝑾- are matrices containing submatrices of

particular solution 𝜿𝑗𝑚

and 𝜻𝑚𝑗 each column of which corresponds to the th𝑚-order radial function and each row to the j nodal point. When introducing the collocation points to improve accuracy, the matrices could be wrote as

,𝑴-=𝜌([𝑮

𝑮𝑩𝑩𝑫𝑩

],𝑷𝑩𝑩𝑷𝑩𝑫-−[𝑯𝑯𝑩𝑩𝑯𝑩𝑫𝑫𝑩𝑯][𝑾𝑩𝑩𝑾𝑩𝑫]),𝑭-−1 (8) in which the subscript 𝐵 and 𝑫𝑫𝑾𝐷 represent the boundary nodes 𝑫𝑩𝟎

and interior collocation points, respectively.

In modal analysis, considering time harmonic dependence for the boundary displacement and traction vectors appearing in Eq. (6). The frequency domain equation form

(−𝜔2,𝑴-+,𝑯-)*𝒖𝟎+=,𝑮-*𝒑𝟎+ (8)

in which 𝜔 is the circular frequency of the harmonic excitation of 𝑢 and 𝑝vectors with amplitude 𝑢0 and 𝑝0, respectively. Just setting the external disturbances equal to zero, one can obtain natural modes and frequencies of vibration.

,𝑨-*𝒙+=𝜔2,𝑴∗-*𝒙+ (9)

in which ,𝐴- is the BEM influence matrix referring to all unknown boundary variables contained in *𝑥+ and ,𝑀∗- is obtained by setting zeros in ,𝑀- where sub-columns refer to specified displacements. From the work of J.P. Agnantiaris [10] for 3D elastodynamics, the effect of augmentation in the linear radial basis function for accuracy is very small, and for non-axisymmetric 3D structures polynomial 1+𝑟 is simplest and high accuracy RBF. So, the corresponding particular solution equation for displacement and traction can be obtained in respect.[16]

3𝜅𝑚

=14𝜇(1−𝜈)*(3−4𝜈)(𝑟26+𝑟12)+𝑟210+𝑟3

𝑖𝑗

18+𝛿𝑖𝑗

−14𝜇(1−𝜈)(2𝑟15+12)𝑟𝑖𝑟𝑗 (10) 𝜍𝑚1𝑖𝑗=

2(1−𝜈)[(1−2𝜈)(1+𝑟)+1+𝑟

](𝑟𝑖𝑛𝑗+𝑟𝑘𝑛𝑘𝛿−1345𝑖𝑗)

2(1[(1−2𝜈)(16+𝑟)−1−𝑟]𝑟𝑗−1−𝜈)𝑛𝑖134562(1−𝜈)∙12𝑟𝑟𝑖𝑟𝑗𝑟𝑘𝑛𝑘

(11) in which 𝜇 is the shear modulus, 𝜈 the Poisson’s ratio, 𝑟 connecting any two points of 𝑖 is the components of the vector 𝑟boundary nodes and interior points, and 𝑛𝑖 is the components of the normal outward vector n at the point where the particular solution is evaluated, and 𝛿𝑖𝑗 is the Dirac function.

3. FAST MULTIPOLE METHOD IN DRBEM

The fast multipole method can be considered providing a method to compute the multiplication between matrix and vector. The conventional algorithm needs 𝑂(𝑁2) complexity to this end, which is too expensive to large problems. On the contrast, the FMM provides the matrix-trial vector product 𝑂(𝑁) operations, which is a huge improvement.

In this section, the expansions of kernel function, such as M2M, M2L and L2L, are brief reviewed at first, more detail information can be found in N Nishimura’s work [15]. Then the new version of FMM algorithm is featured up in detail, that applies new cell structure (node-element cell structure) and objection concepts for improving searching efficiency in procedure of computing the expansion of kernel function.

3.1 Mathematical Equations

In 3D elastodynamics, the expansion of kernel function is different with 2D, which is more complicated as form [18] ∞

𝑛

𝑈1

𝑖𝑗(𝑥,𝑦)=8𝜋𝜇∑∑[𝐹𝑖𝑗,𝑛,𝑚(𝑥−𝑦𝑐)̅𝑅̅̅𝑛,𝑚̅̅̅(𝑦−𝑦𝑐)𝑛=0+𝐺𝑚=−𝑛

𝑖,𝑛,𝑚(𝑥−𝑦𝑐)(𝑦−𝑦𝑐)̅𝑅̅̅𝑛,𝑚̅̅̅(𝑦−𝑦𝑐)] (12) where

𝐹𝜆+3𝜇

𝑖𝑗,𝑛,𝑚(𝑥−𝑦𝑐)=

𝜆+2𝜇𝛿𝑖𝑗𝑆𝑛,𝑚

(𝑥−𝑦𝑐)

−𝜆+𝜇𝜆+2𝜇(𝑥−𝑦𝜕

𝑐)𝑗𝜕𝑥𝑆(𝑥−𝑦𝑐) (13)

𝑖𝑛,𝑚𝐺=

𝜆+𝜇𝜕

𝑖,𝑛,𝑚(𝑥−𝑦𝑐)𝜆+2𝜇𝜕𝑥𝑆(𝑥−𝑦𝑖𝑛,𝑚

𝑐) (14)

Using the expansion of kernel into the integral equation, one could obtain the multipole expansion with substitution point 𝑦𝑐 or 𝑦𝑐

3 Copyright © 2012 by ASME

𝑀𝑗,𝑛,𝑚(𝑦𝑐)=∫𝑅𝑛,𝑚(𝑦−𝑦𝑐)𝑡𝑗(𝑦)𝑑𝑆(𝑦) (15)

𝑆

𝑀𝑛,𝑚(𝑦𝑐)=∫(𝑦−𝑦𝑐)𝑗𝑅𝑛,𝑚(𝑦−𝑦𝑐)𝑡𝑗(𝑦)𝑑𝑆(𝑦) (16)

𝑆

𝑛

𝑛

𝑀𝑗,𝑛,𝑚(𝑦𝑐 )=∑∑𝑅𝑛 ,𝑚 (𝑦𝑐

𝑛 =0−𝑚 𝑦=−𝑛

𝑐)𝑀𝑗,𝑛−𝑛 ,𝑚−𝑚 (𝑦𝑐) (17) 𝑛

𝑛

𝑀𝑛,𝑚(𝑦𝑐 )=∑

∑𝑅𝑛 ,𝑚 (𝑦𝑐 −𝑦𝑐),𝑀𝑛−𝑛 ,𝑚−𝑚 (𝑦𝑐)𝑛 =0+𝑚 (=−𝑛

𝑦𝑐 −𝑦𝑐)𝑗𝑀𝑗,𝑛−𝑛 ,𝑚−𝑚 (𝑦𝑐) (18)

in which 𝑀𝑗,𝑛,𝑚 and𝑀𝑛,𝑚 are multipole expansion of first integral for displacement and traction in DRBEM equation, respectively. The 𝑀𝑛,𝑚(𝑦𝑐 ) and 𝑀𝑗,𝑛,𝑚(𝑦𝑐 ) are M2M expansion when point moved from 𝑦𝑐 to 𝑦move the node point from𝑐 .

FMM also 𝑥 to𝑥𝐿 or 𝑥𝐿 to improve the computational efficiency, which called local expansion. The local expansion coefficients are given by following M2L translations

∞

𝑛

𝐿𝑗,𝑛,𝑚(𝑥𝐿)=(−1)𝑛∑

∑̅̅̅̅̅̅̅̅̅̅̅̅̅̅𝑆𝑛+𝑛 ,𝑚+𝑚 (𝑥𝐿

−𝑦𝑛 =0𝑚 =−𝑛

𝑐)𝑀𝑗,𝑛 ,𝑚 (𝑦𝑐) (19)

∞

𝑛

𝐿𝑛,𝑚(𝑥𝐿)=(−1)𝑛∑

∑̅̅̅̅̅̅̅̅̅̅̅̅̅̅𝑆𝑛+𝑛 ,𝑚+𝑚 (𝑥𝐿

−𝑛 𝑦=0𝑚 =−𝑛

𝑐),𝑀𝑗,𝑛 ,𝑚 (𝑦𝑐)−(𝑥𝐿

−𝑦𝑐)𝑗𝑀𝑗,𝑛 ,𝑚 (𝑦𝑐)- (20)

When the local expansion point is moved from 𝑥L2L translations

𝐿 to 𝑥𝐿 , one have the following ∞

𝑛

𝐿𝑗,𝑛,𝑚(𝑥𝐿 )=(−1)𝑛∑∑𝑅𝑛 −𝑛,𝑚 −𝑚(𝑥𝐿

−𝑥𝑛 =𝑛𝑚 =−𝑛

𝐿)𝐿𝑗,𝑛 ,𝑚 (𝑥𝐿) (21)

∞

𝑛

𝐿𝑛,𝑚(𝑥𝐿 )=(−1)𝑛∑∑𝑅𝑛 −𝑛,𝑚 −𝑚(𝑥𝐿

−𝑛 𝑥=𝑛𝑚 =−𝑛

𝐿),𝐿𝑛 ,𝑚 (𝑥𝐿)

−(𝑥𝐿 −𝑥𝐿)𝑗𝐿𝑗,𝑛 ,𝑚 (𝑥𝐿)- (22)

These formulas are derived from results in Epton MA and Dembart B [17] in terms of the Wigner-3j symbols and concluded by Ken-ichi Yoshida [18].

3.2 Algorithm

We note that Eq. (7) contains two multiplications between matrices as ,𝑮-∗,𝑷- and,𝑯-∗,𝑾-, in which ,𝑮- and ,𝑷- are integral coefficient, ,𝑷- and ,𝑾- are particular solution matrices. And the matrices like ,𝑮- and ,𝑷- is formed through discretization of integration equations, so we could apply FMM into,𝑮-∗,𝑷- and,𝑯-∗,𝑾-. The infinite sum of numerical computation (approximation for integral computation) will be truncated by p, which is a certain number defined in advance.

a

a

b

Figure 1 Form hierarchical oct-tree structure (a,b) Comparing with the conventional FMM algorithm, some improvements have been taken in this paper. First, instead of interaction between points, we apply the point to element interaction method to cancel plenty of repeated traversal and computation for saving lots of CPU time. So there will be two sets of cell data, one stores the node point data and the other stores the element data. Second, during the realization of the algorithm, this paper presented position location method to avoid travelling all the cells in direct computational, which is very efficient in large problems. For realizing this destination, we apply the STL (standard template library) to form cell object for storing data and operation function, instead of using simply c language structure.

In traditional way, each cell only stores the node point data. It must find all elements adjacent to the field node point and compute the interaction between source point and all the elements. (As Figure 2) The interaction computation is result expressed as 𝜂 and 𝑤 express the weight of correspond interaction computation. So the value of correspond position in the coefficient matrix is achieved by the following equation:

𝜂1∗𝑤1+𝜂2∗𝑤2+∙∙∙+𝜂6∗𝑤6

It is easy to find that this method consumes plenty of time in computing repeated interaction results. We apply two sets of

4 Copyright © 2012 by ASME

cell data to avoid repeated computation that is named node-element structure. Firstly, forming all the cell objects according to the node point position in conventional FMM method, then all the elements will be put into the cell objects according to the shape center position of each element.

Figure 2 Point to point interaction

Just as the Figure 3, if the node points of one element (E1, E2, E3, E4, E5) locate in different cells, the shape center decides which cell the element belongs. (E1, E5, E3 belong to Cell1; E2, E4 belong to Cell2)As the result, all the integral during element for each source point just need to be computed once.

This FMM algorithm applied oct-tree structure for complicated situation instead of quad-tree in 2D problems, more details goes as follows. Note that, the number of cells in level l is start from left to right, front to back, downside to upside in turn.

Figure 3 Point to element interaction method  Discretize boundary S into EN elements and N nodes.  Obtain hierarchical tree structure of cells for nodes and

elements. Namely, take a cube box contain all the body of the modal (See Figure 1 (a)), which called cell 0 of level 0. Then we take a cell (a parent cell named icell) of level l to divide into eight sub cubes. If the number of nodes contained in the sub cube does not equal to zero, which will be called child cell of icell in level l+1. A cell of

level l+1 is further divide into eight sub cubes if the number of nodes contained in the cell is more than limitnumber (determined beforehand). (See Figure 1 (b) and Figure 4) Otherwise one terminates this subdivision and this icell would be called leaf.

After the node oct-tree structure, it will be necessary to form another element oct-tree structure for computation in future. (Notice that the element cell must be node cell, the inversion does not work.) Before this operation, we need compute the shape center point of elements (See red points in Figure 4), and then mark every element to corresponding cell (through coordinate location). During the process of forming cell by subdivision, we should store the position (i) of each cell in current level and renumber the node and element number for efficiency of code. Detail pseudo-code is given in Function1, 2.

Function 1 Form oct-tree structure

1: Form root cell 0 in level 0 2: divide the cell into 8 sub cubes 3: for i=0 to 8 do 4: Build Recursion computation for sub cube 5: cell 0 in level 0 6: end for Function 2 Recursion computation (level l)

1: if node numbers of sub cube greater than 2: limit number then 3: form node cell icell 4: icell levell

5: if element numbers of icell not equal to 6: zero then

7: icell also element cell

8: divide the cell into 8 sub cubes 9: for i=0 to 8 do

10: Build Recursion computation in 11: Level l+1 12: end for

13: else icell only node cell 14: end if

15: else if node numbers of child blocks less than 16: or equal to limit number then

17: icellleaf

18: else this block not cell 19: end if

Upward pass. Applied Eq. (15-16) to compute moments in each leaf element cell (compute with element cell data). For non-leaf cells, add child cell’s moments to current icell (parent cell) for obtaining the current cell’s moment, this is called M2M and realized by Eq (21-22). This process is actually transforming interaction source point from child cell center to parent cell center. In this manner, we just need tracing the element tree structure from lowest level and terminate in level 2.

5 Copyright © 2012 by ASME



Figure 4 Tree structure



Downward pass. Firstly, tracing the node tree cells (icell) for computing the interaction with all the elements to form coefficient matrix.

Through the position i which has been recorded when forming oct-tree structure, one could find neighbor cells by Eq. (21-24) instead of tracing all cells and computing the distance applied in conventional FMM. (See Figure 2) That is position location algorithm applied in the modified FMM. During the method, we do not need to travel all the cells for detecting which cell should be computed in directly. Whatever the magnitude of the problem is, the computational complexity for searching neighbor cells is certain, the efficiency will be more obvious as the scale of problem is bigger.

𝑃𝑜𝑠𝑖𝑡𝑖𝑜𝑛=(𝑧+𝑚)∗2𝐷𝑖𝑣+(𝑦+𝑛)∗𝐷𝑖𝑣+𝑥+𝑙,0<(𝑧+𝑚)<(𝐷𝑖𝑣⁄2−1)

{0<(𝑦+𝑛)<(𝐷𝑖𝑣⁄2−1) (21) 0<(𝑥+𝑙)<(𝐷𝑖𝑣⁄2−1)

in which

𝑧=𝑖/2𝐷𝑖𝑣 (22) 𝑦=(𝑖%2𝐷𝑖𝑣)/𝐷𝑖𝑣 (23) 𝑦=(𝑖%2𝐷𝑖𝑣)%𝐷𝑖𝑣 (24)

In which 𝐷𝑖𝑣 is the numbers of division in level l, and m, n, l is neighbor cell position (front or back, left or right, up or down) with respect to icell, with 1, 0 and 1. N is the number of cells we obtained in algorithm, there will need N computer time in conventional way for determining when we need directly integration. And now, 27 determinations are needed in mostly to locate neighbor cells position number by position location method. For example, the position of orange cell in Figure 5 is 63 (i=63, start from 0), and the neighbor cell’s position number is 42, 43, 46, 47, 58, 59, 62, respectively (the silvery cells in Figure 3). The interaction with neighbor cell will be computed directly by using conventional BEM if icell or ccell is leaf cell.

Figure 5 Obtain position number of neighbor cell

If ccell is not adjacent with icell and the parents of them are neighborhood, this relationship is called interaction. We start from level 2 for tracing all the node tree structure to obtain icell, and obtaining interactionlist which contain ccell by tracing element tree structure. The local expansion will be computed by M2L that is transforming filed point from current cell center to child cell center by using Eq. (19-20) in icell and ccell. And the result will be recorded in the node tree structure cell icell. This manner will be terminated when the level l equal to lowest level. Computing the L2L by Eq. (21-22) while the ccell is not interaction or neighbor with icell during the procession of tracing level from 2 to lowest level.



Finally. As tracing to the leaf, applying local expansion values stored in each computation by above step into Eq (12), we can obtain the multiplication result between matrix and vector with respect to each node on the boundary.

4. NUMERICAL EXAMPLE

We can find the complexity reduces from𝑂(𝑛4) to 𝑂(𝑛3) during calculation and the complexity 2)of coefficient matrices formulation reduces from 𝑂(𝑛 to 𝑂(𝑛) through this algorithm. The theory is demonstrated by numerical example in this section.

One 3D elastodynamic problem dealing with the free vibration of an elastic bearing support is solved by DRBEM with modified FMM, which is used to examine the efficiency and accuracy of the algorithm presented in this paper. The time consumed in the method is listed in Table 1 Comparison of time needed with conventional, and the natural frequencies are presented in Figure 7. The material parameters 6for the above problem are: Shear modulus 𝜇=10𝑃𝑎, Poisson’s ratio 𝜈=0.3, and mass density 𝜌=7400𝑘𝑔/𝑚3.

The first issue in concerned is the efficiency of modified FMM applied in DRBEM with respect to conventional method. The time consumed for multiplication between matrix and vector was recorded while the element numbers are different (this example enumerates five different numbers of element, 1276, 5076, 10336, 29466, 52684, respectively. See Figure 6 ). All the results are presented in Table 1 Comparison of time needed with conventional and are compared with those data obtained by conventional way. (For the limit of memory we cannot obtain the result while the number of elements

6 Copyright © 2012 by ASME

increasing more than 6000 in conventional method, so the second row of Table 1 Comparison of time needed with conventional just presents two columns.) As the evident of this table, there is a good improvement in efficiency by introducing the conventional FMM, and we also notice that the modified FMM have a very nice performance comparing with traditional. The effect of modified algorithm is more and more obvious as the number of elements increase, which could be found from Table 1 Comparison of time needed with conventional.

During the first issue, we find that the value of limitnumber has great influence for results. The cost of subdivision will be high when the value is too small. On the other hand, the error will be unacceptable if the value selected is too large. Thus, this parameter is very important in analysis.

Table 1 Comparison of time needed with conventional Element 1276 5076 10336 29466 52684 numbers (20) (100) (200) (400) (850) (limit number)

Conventional 1.0201 20.116 / / / BEM (s) DRBEM 0.9398 2.0017 6.483 36.1429 124.902 +FMM(s) DRBEM+ 0.9221 1.7631 4.3227 29.7522 91.5834 Modified FMM (s)

The second issue is the accuracy of results computed by the presented algorithm. High accuracy is an important character of BEM with respect to FEM and FDM. However, due to the approximation of unknown 𝒖(𝒙,𝑡) (See Eq. (4)), the accuracy of DRBEM has been affected, but which could be fixed by introducing collocation points during computation. After experiment we find that the result of DRBEM agrees with FEM, and the accuracy is high while taking appreciated number of collocation points. In this example 10 order natural frequencies is extracted in the presented algorithm and using in turn 0, 101, 165 and 241 uniformly distributed interior collocation points. Figure 7 is the result comparison with FEM (computed in Ansys and Hyperworks software).

The graphic information of Figure 7 (a-d) conveys that the results achieved by the novel algorithm are reliable for which have good agreement with FEM. Especially, during low frequencies, the agreement of results is significant. However, just as mentioned above, the accuracy could improve as introducing the interior collocation points, which could be concluded by the comparison during Figure 7 (b, d). But this kind of improvement is not endless and even could lead to horrible results, which could be observed by Figure 7 (c). Through this graphic, we could find that the results have great errors in high frequencies. From order 5 to higher frequency, the differences with FMM (computed by Ansys and Hyperworks) are huge and unacceptable. In Figure 7 (b, d), this phenomenon is disappearing for different number of interior

points are counted, and the accuracy has improved with respect to Figure 7 (a). Thus, it should be cautious for the choice of interior points when we want to achieve more accurate result. More details about selection of the collocation points can refer to J.P. Agnantiaris’s work [10][16].

Figure 6 Different element numbers during mesh There is still another drawback: it will consume lots of computer resource with introducing plenty of collocation points, and the accuracy would not improve with the amount of points in some time.

5. CONCLUSION

The character of reducing dimension in BEM could simplify the preprocessing such as mesh procedure and the degeneration of time dimension in DRBEM makes the complexity from 𝑂(𝑛4)∗𝑂(𝑚4) to 𝑂(𝑛4). DRBEM-FMM for 3D elastodynamic modal analysis has been proposed in this paper which is efficient and has satisfactorily accurate. Absorbing the superiorities such as static fundamental solutions and similar formulation with FEM, this algorithm avoids the defect things like high cost of matrix multiplication by introducing FMM. At the same time, we have made some improvements in traditional FMM. This algorithm saving plenty of repeated computation time by using direct interaction between point and element; the PL method reduces one layer loop during nested loops.

In the numerical example, we tested the efficiency of FMM for DRBEM and observed a nice performance. The accuracy of the novel algorithm has been tested through comparison with FEM. Using DRBEM with fast multipole method would reduce the accuracy of result through the analysis of equations, however, the forth graphics show that the agreement of results is totally great with FEM, the error could be acceptable. At the end of the numerical example, we found that although interior collocation points could improve the accuracy of results in some condition, but it also need additional CPU and memory cost. So we should balance the efficiency and accuracy when to apply this algorithm in practical problems.

7 Copyright © 2012 by ASME

(a)

(b)

8 Copyright © 2012 by ASME

(c)

(d)

Figure 7 Comparison of results with FMM (Ansys and Hyperworks)

9

Copyright © 2012 by ASME

The present algorithm not only can be applied in modal analysis, but also for transient response problems, harmonic response analyses, and so on. Although the efficiency of the algorithm is significant, some large magnitude problem still cannot be solved in this way for the limitation of memory space. We need to finger out some methods for solving mass storage problem and make this algorithm suiting for large-scale problems.

ACKNOWLEDGMENTS

This research has been supported by the National Natural Science Foundation of China (50975107).

REFERENCES

[1] Patrick Dangla et al., 2005. A simple and efficient

regularization method for 3D BEM: Application to frequency-Domain elastodynamics. Bulletin of the Seismological Society of America, pp. 1916-1927

[2] Chen YH, Chew WC and Zeroug S., 1997. Fast multiple

method as an efficient solver for 2D elastic wave surface integral equations. Comp Mech, pp. 495–506.

[3] Fukui T and Inoue K., 1998. Fast multipole boundary

element method in 2Delastodynamics. Journal of Applied Mechanics JSCE, pp. 373–80.

[4] Fujiwara H., 1998. The fast multipole method for integral

equations of seismic scattering problems. Geophys J Int, pp. 773–82.

[5] Takahashi T, Nishimura N and Kobayashi S., 2001. Fast

boundary integral equation method for elastodynamic problems in 2D in time domain. Trans JSME, A67-661, pp. 1409–16.

[6] Yoshida K, Nishimura N and Kobayashi S., 2000. Analysis

of three dimensional scattering of elastic waves by crack with fast multipole boundary integral equation method. Journal of Applied Mechanics JSCE, pp. 143–50.

[7] Fujiwara H., 2000. The fast multiple method for solving

integral equations of three-dimensional topography and basin problems. Geophys J Int, pp. 198–210.

[8] Lu M, Wang J, Ergin AA and Michielssen E., 2000. Fast

evaluation of two dimensional transient wave fields. J Comput Phys, pp. 161–85.

[9] Nardini D and Brebbia CA., 1982. A new approach to free

vibration analysis using boundary elements. In: Boundary Element Methods in Engineering. Springer, Berlin, pp. 313-325.

[10] J.P. Agnantiaris, D. Polyzos and D.E. Beskos., 2001. Free

vibration analysis of non-axisymmetric and axisymmetric structure by the dual reciprocity BEM. Elsevier, pp. 713-723.

[11] Ahmad S and Banerjee PK., 1986. Free vibration analysis

by BEM using particular integrals. Journal of Engineering Mechanics ASCE, pp. 682-95.

[12] Rokhlin V., 1985. Rapid solution of integral equations of

classical potential theory, J. Comput. Phys, pp. 187–207. [13] Greengard L., 1988. The Rapid Evaluation of Potential

Fields in Particle Systems, MIT Press, Cambridge, MA. [14] Greengard L and Rokhlin V., 1987. A fast algorithm for

particle simulations, J. Comput. Phys, pp. 325–348.

[15] N Nishimura., 2002. Fast multipole accelerated boundary

integral equation methods. Appl. Mech. Rev, pp. 299-325. [16] J.P. Agnantiaris et al., 1998. Three-dimensional structural

vibration analysis by Dual Reciprocity BEM. Computational Mechanics, pp. 372-381.

[17] Epton MA and Dembart B., 1995. Multipole translation

theory for the three-dimensional Laplace and Helmholtz equations, SIAM J. Sci. Comput. (USA), pp. 865–897.

[18] Ken-ichi Yoshida., 2001. Applications of Fast Multipole

Method to Boundary Integral Equation Method. PHD Thesis.

10 Copyright © 2012 by ASME