Optimized meshfree methods for acoustics

Comput. Methods Appl. Mech. Engrg. 200 (2011) 2223–2236

Contents lists available at ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

Optimized meshfree methods for acoustics Christina Wenterodt ⇑, Otto von Estorff Institute of Modelling and Computation, Hamburg University of Technology, Denickestr. 17, D-21073 Hamburg, Germany

a r t i c l e

i n f o

Article history: Received 23 November 2010 Received in revised form 10 March 2011 Accepted 16 March 2011 Available online 30 March 2011 Keywords: Meshfree methods EFGM RPIM Dispersion Helmholtz equation Acoustics

a b s t r a c t When the Helmholtz equation is solved by numerical methods as, e.g., the finite element method (FEM), the solution suffers from the so-called pollution effect. The pollution is mainly caused by the dispersion, meaning that the numerical wave number disagrees with the wave number of the exact solution. This leads to inaccurate results, especially for high wave numbers. In order to obtain acceptable results also for higher wave numbers, either a high element resolution or elements of a higher order can be used. For either option the consequence is an increased computation time and memory capacity. Meshfree methods as the element free Galerkin method (EFGM) and the radial point interpolation method (RPIM) are not dispersion-free either, but it has been shown that meshfree methods are able to reduce the dispersion significantly. Both methods offer several parameters, which can be modified in order to obtain optimal results with respect to the dispersion effect. This work presents an exhaustive parameter study on both the EFGM and the RPIM. It is shown, that the methods can be significantly improved if certain parameters as, e.g., weighting functions, shape parameters, size of the influence domain, are chosen appropriately. Ó 2011 Elsevier B.V. All rights reserved.

1. Introduction For acoustic computations in the low-frequency range, the finite element method (FEM) is established as a well-known standard tool. During the last years considerable efforts have been undertaken in order to improve the accuracy and efficiency of the method. Nevertheless, it is proved that the accuracy of the FEM decreases for high wave numbers, which is denoted as pollution effect [10,11]. The pollution is mainly a consequence of the dispersion effect, meaning that the wave number of the numerical solution and the wave number of the exact solution disagree. The dispersion effect was first examined for the wave equation [21] and later also for the Helmholtz equation [4,9]. In order to obtain acceptable solutions for high wave numbers, the FEM has to be applied using a refined mesh (h-FEM) or a higher order of the polynomial approximation (p-FEM). Thereby, the accuracy of the FEM solutions can be improved, but large systems of equations and hence an increased computation time and memory capacity are the consequences, i.e., the efficiency of the FEM decreases. For this reason, a reduction of the dispersion is strongly desirable. Several methods as for example the Galerkin leastsquares (GLS) FEM [27] and the residual-based FEM [22] have been proposed to reach this goal. The partition of unity finite element method (PUFEM) includes a priori knowledge about the solution

⇑ Corresponding author. E-mail address: [email protected] (C. Wenterodt). 0045-7825/$ - see front matter Ó 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2011.03.011

of the problem into the finite elements which can also reduce the dispersion [20]. As a further development of the PUFEM, in [14] a system of plane waves is introduced into the finite elements. A disadvantage of methodologies that are based on the PUFEM is the fact that ill-conditioned system matrices can be obtained and thus the use of classical iterative solvers is difficult [8,23]. Shape functions which reduce the pollution while preserving the positive characteristics of the system matrices are presented in [24]. More recently, an edge-based smoothed finite element method (ES-FEM) was developed [18] and applied to acoustics [7]. Compared to the FEM, the ES-FEM is able to reduce the dispersion effect and very accurate results can be obtained. With respect to the FEM, meshfree methods are relatively new in the field of numerical simulations. These methods do not need a mesh to construct the shape functions. Instead, shape functions are computed based on an arbitrary field point distribution inside the domain and on its boundary. Various meshfree methods, especially the element free Galerkin method (EFGM) and the radial point interpolation method (RPIM), are introduced and compared in [6,15,16]. The EFGM has originally been developed by Belytschko et al. [1]. Meanwhile, it has been applied to numerous engineering problems as e.g. static and dynamic fracture problems [2]. Furthermore, the EFGM has successfully been coupled to meshbased methods as the FEM and the BEM in order to study soil–structure interaction [28]. In the field of acoustic computations, Bouillard et al. [3,25,26] showed that the dispersion effect of the EFGM is significantly reduced in comparison to the FEM.

2224

C. Wenterodt, O. von Estorff / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2223–2236

The RPIM, which has been developed by Wang and Liu [30,31], has some advantages over the EFGM. For example, the RPIM shape functions fulfill the Kronecker delta function property which simplifies the imposition of essential boundary conditions. For computations in the field of structural mechanics the RPIM has shown to perform well [17,19]. In [33] the RPIM has been applied to acoustic problems. Regarding the dispersion effect, it was shown that the RPIM leads to better results than the FEM. In this paper the EFGM and the RPIM are optimized with respect to the dispersion effect by conducting several numerical experiments. Miscellaneous parameters as e.g. the weighting functions of the EFGM, the radial basis functions of the RPIM, the size and form of the influence domains, the integration order, etc. are varied. Section 2 introduces the interior acoustic problem and the associated weak formulation. Section 3 deals with the basics of meshfree methods. Furthermore the shape functions for EFGM and the RPIM are derived. The theory regarding the dispersion effect is outlined in Section 4. In Section 5 the parameter study is presented and optimal parameters for the EFGM as well as for the RPIM are suggested.

integrated over their respective domain. Dirichlet boundary conditions (4) have to be fulfilled by the approximation for p. If Green’s theorem is applied and the surface terms are cancelled out, the weak formulation of the given problem reads

Z

2

ðrp  rw  k pwÞdX þ

Z

ixqA0 pw dS ¼ 

SR

X

n X

ph ðxÞ ¼

ui pi

as an approximation for the sound pressure, with given shape functions ui and coefficients pi to be determined. Substituting (8) into the weak formulation (7), and choosing identical weight and shape functions, results in the linear equation system

Kij ¼

Consider an acoustic fluid domain X. The propagation of pressure waves inside X is governed by the wave equation which is derived using the ideal gas law as well as the balance of mass and momentum. Assuming that the state variables pressure P, velocity v, and density q, experience only small variations, the wave equation

ð9Þ

Z

ruTi ruj dX;

X

ð10Þ

the acoustic damping matrix

Cij ¼ qc

Z SR

A0 ui uj dS;

ð11Þ

the acoustic mass matrix

1 @2P ¼0 c2 @t 2

ð1Þ Mij ¼

is obtained, with c denoting the sound velocity. For time harmonic waves of frequency x the pressure can be expressed as ixt

g

ð2Þ

with p being the complex amplitude of the sound pressure. Introduction of the wave number k = x/c leads to the Helmholtz equation 2

ð8Þ

with the acoustic stiffness matrix

2.1. Helmholtz equation

Pðx; tÞ ¼ RefpðxÞ  e

ð7Þ

i¼1

2

r2 P 

ixqv 0 w dS:

SN

If a numerical solution of the boundary value problem should be obtained, Eq. (7) must be discretized and appropriate shape functions for p and weight functions w have to be introduced. For Galerkin methods, the shape and weight functions are chosen identically. Consider

½K þ ikC  k Mp ¼ ikF 2. Interior acoustics

Z

2

r p þ k p ¼ 0:

ð3Þ

In order to obtain a well-posed boundary value problems, appropriate boundary conditions have to be prescribed. For this reason, the boundary S of the bounded fluid domain X is decomposed into three distinct regions SD, SN, and SR such that S = SD [ SN [ SR. The respective boundary conditions are  Dirichlet boundary conditions:

 on SD : p¼p

ð4Þ

 Neumann boundary conditions:

@p ¼ ixqv 0 @n

on SN :

ð5Þ

 Robin boundary conditions:

@p ¼ ixqA0 p on SR : @n

Z X

Here, v0 represents the normal velocity, n is the normal vector and A0 the admittance.

2.2. Weak formulation The weak formulation of the Helmholtz equation can be obtained using the method of weighted residuals. The residuals of Eqs. (3), (5) and (6) are multiplied by weight functions w and

ð12Þ

and the acoustic excitation vector

Fi ¼ qc

Z SN

v 0 ui dS:

ð13Þ

The vector p contains the unknown coefficients pi. It is important to note that the matrices K, C and M and the vector F are independent of the wave number. For this reason, if a given acoustic problem has to be solved for several different frequencies, it is not necessary to assemble the matrices for each frequency anew. 3. Meshfree methods In comparison to mesh-based methods as, e.g., the FEM, in the case of meshfree methods no mesh is needed for the representation of the considered domain X. In X and on its boundary S, field points xi can be arbitrarily distributed. An approximation of the sound pressure ph at an arbitrary point xQ in X is obtained by an interpolation of the form

ph ðxQ Þ ¼ ð6Þ

ui uj dX;

n X

ui ðxQ Þpi :

ð14Þ

i¼1

Here, pi is the coefficient of the solution vector p belonging to field point xi and ui is the associated shape function which is constructed using all n field points that influence the considered point xQ. In order to obtain an equation system (9) that can be solved for the coefficient vector p, the shape functions are substituted into the weak formulation (7) of the Helmholtz equation, both as shape and as weight functions. Since the resulting equations have to be integrated numerically, a background mesh has to be introduced, which can be chosen independently of the field point distribution and which does not have to represent the domain X (see Section 3.2).

2225

C. Wenterodt, O. von Estorff / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2223–2236

3.1. Influence domain Concerning meshfree methods, a shape function belonging to a field point xi depends on the point xQ where the shape function is evaluated. For this reason, an influence domain is defined for each field point xi. In the two-dimensional case this domain can be a circle or a square with a characteristical length rinfl. In the case of a square, rinfl denotes the half edge length, for a circle it is the radius. Other shapes, e.g. rectangles, are also permitted. Moreover, the influence domain can be different for each field point. The influence domain of a field point xi prescribes in which region of the domain X the respective field point is used for the construction of the shape functions. This also means that two field points xi and xj are associated if their influence domains intersect (cf. Fig. 1). 3.2. Background mesh In order to obtain the system matrices given in Eqs. (10)–(12) an integration over the domain X and the boundary S has to be carried out. Usually a numerical integration is conducted and thus an integration grid is needed. This so-called background mesh does not depend on the field points and can be chosen rather coarse. For example, the matrix K can be obtained by evaluating

Kij ¼

nel Z X el¼1

Xel

ruTi ruj dX

ð15Þ

with nel being the number of elements of the background mesh. Integration over an element of the background mesh is carried out according to

Z Xel

domain. On the other hand, the constructed shape functions lack the Kronecker delta property, which makes the imposition of Dirichlet boundary conditions complicated. In order to obtain the shape functions of the EFGM, the unknown sound pressure ph is interpolated by

ph ðxÞ ¼

ru ruj dX ¼

X

W Q detðJel Þrui ðxQ Þ ruj ðxQ Þ;

ð16Þ

PT ðxÞ ¼ ½1 x y xy    xm ym 

3.3. EFGM

ð18Þ

in the two-dimensional case. The basis can also be enriched with special terms in order to better reproduce special characteristics of the solution, as described in Section 5.3. The coefficient vector a(x) is determined by minimization of a functional of weighted residuals n X

J¼

wi ðxÞ½ph ðxÞ  pi  ¼

i¼1

n X

wi ðxÞ½PT ðxÞaðxÞ  pi 

such that the approximation ph(xi) at the field points yields the best possible agreement with the according coefficients pi of the solution vector. The weighting functions wi(x) can be chosen, e.g., as the cubic spline

8 2 2 3 > < 3  4r þ 4r 4 2 wi ðrÞ ¼ 3  4r þ 4r  43 r 3 > : 0

for r 6 12 ; for

1 2

< r 6 1;

ð20Þ

for r > 1;

the quartic spline

wi ðrÞ ¼

1  6r 2 þ 8r 3  3r 4

for r 6 1;

0

for r > 1;

ð21Þ

or the exponential function

( wi ðrÞ ¼

expðð2rÞ2 Þexpð4Þ 1expð4Þ

for r 6 1;

0

for r > 1;

ð22Þ

with r denoting the distance between point x and field point xi scaled by the influence radius rinfl, i.e.

r¼

kx  xi k : r infl

ð23Þ

Minimization of the functional J leads to

aðxÞ ¼ A1 ðxÞBðxÞp The element free Galerkin method (EFGM) is based on the moving least squares (MLS) approximation. An important feature of the MLS approximation is its capability of producing an approximation with any desired order of consistency. Furthermore, the approximated field function is continuous and smooth in the entire

ð19Þ

i¼1

Q¼1

where nQ denotes the number of Gauss points in the respective element, WQ is the associated Gauss weight and Jel is the element’s Jacobian. Regarding meshfree methods, in most cases, the shape functions u have to be computed individually for each Gauss point. The acoustic mass and stiffness matrix can be obtained accordingly. Especially for the EFGM, the choice of an appropriate background mesh and integration scheme is a delicate issue, which will also be examined in Section 5. An extensive study on the numerical integration procedure in meshfree methods can be found in [15].

ð17Þ

where n is the number of terms in the basis PT(x). In general, a complete polynomial basis of order m is chosen as, e.g.,

( T

Pi ðxÞai ðxÞ ¼ PT ðxÞaðxÞ;

i¼1

nQ

T i

n X

ð24Þ

with

AðxÞ ¼

n X

wi ðxÞPðxi ÞPT ðxi Þ;

ð25Þ

i¼1

BðxÞ ¼ ½w1 ðxÞPðxi Þ; . . . ; wn ðxÞPðxn Þ:

ð26Þ

Substituting Eq. (24) into Eq. (17) yields

p ðxÞ ¼ PT ðxÞA1 ðxÞBðxÞp ¼ UðxÞp h

ð27Þ

with the vector of shape functions

UðxÞ ¼ PT ðxÞA1 ðxÞBðxÞ:

ð28Þ

3.4. RPIM

Fig. 1. Intersecting influence domains of two field points xi and xj: at xQ the shape function is constructed under consideration of both field points.

In comparison to the MLS approximation used for the construction of the EFGM shape functions, point interpolation methods as the RPIM obtain approximations which interpolate the function values at all fieldpoints. For this reason, the shape functions

2226

C. Wenterodt, O. von Estorff / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2223–2236

possess the Kronecker Delta property which makes the imposition of Dirichlet boundary conditions straightforward. RPIM shape functions are constructed under the assumption that an approximation ph of the sound pressure can be represented by h

p ðx; xQ Þ ¼

n X

T

Ri ðxÞai ðxQ Þ ¼ R ðxÞaðxQ Þ:

ð29Þ

and

0

 The multiquadric basis

R1 ðr n Þ R2 ðr n Þ    Rn ðrn Þ

rk ¼

2 q

Ri ðxÞ ¼ ðr þ ðac dc Þ Þ

ð30Þ

with shape parameters ac and q.

!

c

Þ

ð32Þ

Here, r is the distance between the ith field point xi and the point x. The characteristic length dc denotes the average field point spacing of all the field points in the influence domain. Shape parameters ac, q, wc, and ec can be adapted to the problem in order to optimize the results. Investigations regarding the optimal choice of shape parameters have been carried out by Wang and Liu. In [30] they find optimal shape parameters for the simple RPIM with multiquadric and Gaussian basis but also show that those parameters are no longer optimal for the RPIM with polynomial reproduction (see Section 3.5). Furthermore, when applied to acoustic problems, the shape parameters suggested in [30] are not optimal either, which is shown in Section 5. The Wendland functions belong to the compactly supported, strictly positive definite basis functions and have been introduced by Wendland [32]. They can be constructed to have any desired degree of smoothness, e.g. the function given in Eq. (31) is in C6. The Wendland radial basis functions have been applied in collocation methods [5,13]. In conjunction with the RPIM, Wang et al. [29] studied different Wendland functions in comparison to the standard radial basis functions. In Eq. (29), n denotes the number of field points that influence xQ and ai(xQ) is the coefficient for the ith basis function. The components ai of vector a(xQ) = [a1a2   an]T are chosen such that Eq. (29) is satisfied at all n field points that influence xQ. For field point i this yields

for i ¼ 1 . . . n;

ð33Þ

with pi being the coefficient of the solution vector at a field point xi. Since the RPIM shape functions possess the Kronecker delta property, the coefficient pi is equivalent to the sound pressure at the respective field point. Eq. (33) can also be written as a linear system of equations

p ¼ RQ a;

ð34Þ

with

p ¼ ½ p1

is used. It can be shown that the inverse of the n  n-dimensional moment matrix RQ always exists, such that Eq. (34) can be rearranged as

ð38Þ

ph ðxÞ ¼ UðxÞp:

ð39Þ

The matrix U contains the shape functions ui and can be computed according to

with shape parameter ec.

pi ¼ RT ðxi Þa

ð37Þ

Substituting Eq. (38) in (29) yields

ð31Þ

with shape parameter wc  and the Gaussian basis 2 =d

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðxi  xk Þ2 þ ðyi  yk Þ2

a ¼ R1 Q p:

 The Wendland basis

Ri ðxÞ ¼ eðec r

ð36Þ

For the two-dimensional case, in Ri(rk) the distance between the ith and the kth field point

Here, Ri(x) are radial basis functions, e.g.

 8  3  2 r r r r Ri ðxÞ ¼ 1   32 þ 25 þ8 þ1 wc wc wc wc

1

B R1 ðr2 Þ R2 ðr2 Þ    Rn ðr 2 Þ C C B C RQ ¼ B . .. . C: B .. . . @ . . . . A

i¼1

2

R1 ðr1 Þ R2 ðr1 Þ    Rn ðr 1 Þ

UðxÞ ¼ RT ðxÞR1 Q :

ð40Þ

The derivatives of the shape functions that are needed in order to compute the acoustic stiffness matrix K in Eq. (10) are obtained by simple derivation of the basis functions Ri. 3.5. RPIM shape functions with polynomial reproduction A drawback of the RPIM based on the simple shape functions described in Section 3.4 is its lack of consistency. In order to be able to reproduce a linear field exactly, polynomial terms have to be inserted into the basis [15]. Analogously to Eq. (29), the ansatz for the RPIM with polynomial reproduction is given by

ph ðx; xQ Þ ¼

n X

Ri ðxÞai ðxQ Þ þ

i¼1

m X

qj ðxÞbj ðxQ Þ

j¼1

¼ RT ðxÞaðxQ Þ þ qT ðxÞbðxQ Þ:

ð41Þ

The polynomial basis vector can be chosen as for the EFGM in Eq. (18). The coefficients ai and bj of the vectors a and b must be determined such that Eq. (41) is fulfilled at all n field points influencing point xQ. This leads to the system of equations

p ¼ RQ a þ Q m b:

ð42Þ

Moreover, to ensure unique approximation of a function, the condition

Q Tm a ¼ 0

ð43Þ

has to be fulfilled. The moment matrix RQ is constructed according to Eq. (36). The moment matrix Qm resulting from a basis, e.g.,

qT ðxÞ ¼ ½ 1 x y xy 

ð44Þ

reads

0

x1 y 1

1

1 x1

y1

B 1 x2 B Qm ¼ B B .. .. @. .

y2 .. .

x2 y 2 C C C .. C: . A

1 xn

yn

xn y n

ð45Þ

Since RQ is invertible, Eq. (42) can be transformed into

p2

   pn T

ð35Þ

1 a ¼ R1 Q p  R Q Q m b:

ð46Þ

2227

C. Wenterodt, O. von Estorff / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2223–2236

Substituting (46) in Eq. (43) yields

b ¼ Sb p

and

ð47Þ

b ¼ ikF

ð54Þ h

with

h

Sb ¼ Q Tm R1 Q Qm

i1

Q Tm R1 Q :

are introduced. Assuming that the numerical solution p is of the form

ð48Þ

A substitution of b in Eq. (46) by (47) yields

a ¼ Sa p

h

with

ð50Þ

Finally, Eq. (41) can be written as

ph ðxÞ ¼ ½RT ðxÞSa þ qT ðxÞSb p ¼ UðxÞp;

ð51Þ

ð52Þ

4. The dispersion effect The numerical computation of wave propagation problems, as e.g. the Helmholtz equation, is subjected to the so-called pollution effect: with increasing wave number k the accuracy of the numerical solution decreases. The main reason for this phenomenon is the dispersion, which means that the wave number of the numerical solution does not match the wave number of the exact solution. The pollution effect has been studied intensively for the FEM [4,9,10,12]. Concerning meshfree methods, Suleau et al. [3,25,26] have analyzed pollution and dispersion of the EFGM. In [33], the RPIM is investigated with respect to the dispersion effect. The numerical wave number kh can be written as a function of the exact wave number k. In the following, a regular field point distribution is assumed with h being the distance between two adjacent field points. The influence domain of each field point is a rectangle with influence radius rxinfl in x- and ryinfl in y-direction, as shown in Fig. 2. The numerical solution of the acoustic problem (3) with boundary conditions (4)–(6) can be obtained by solving the system of Eq. (9). In the following, the abbreviations 2

A ¼ K þ ikC  k M

k2 ¼ k ðhÞ sin h

ð55Þ h

where h is denoting the propagation angle and k the numerical wave number. Then, for an interior field point xr1 ;r2 ¼ ðxr1 ; yr2 Þ, the equation s1 X

s2 X

Ar1 þj1 ;r2 þj2 pðxr1 þ j1 h; yr2 þ j2 hÞ ¼ 0

ð56Þ

j1 ¼s1 j2 ¼s2

where the matrix of the shape functions U can be obtained from

UðxÞ ¼ ½RT ðxÞSa þ qT ðxÞSb :

h

with k1 ¼ k ðhÞ cos h; ð49Þ

Sa ¼ R1 Q ½1  Q m Sb :

ph ðx; yÞ ¼ eiðk1 xþk2 yÞ

ð53Þ

must hold. Here, Ar1 ;r2 is the element of matrix A that relates to field point xr1 ;r2 . The influence domains of the field points are chosen such that 2s1 + 1 points in x- and 2s2 + 1 points in y-direction are associated (see Fig. 2). For the derivation of the numerical wave number, an interior field point which is not subjected to any boundary conditions is considered. Thus, the damping matrix C in Eq. (53) and the vector b in Eq. (54) vanish, such that the right hand side of Eq. (56) equals zero. Substituting Eq. (55) into (56) yields s1 X

s2 X

h

h

Ar1 þj1 ;r2 þj2 expðik j1 h cos hÞ  expðik j2 h sin hÞ ¼ 0:

ð57Þ

j1 ¼s1 j2 ¼s2

With

Cj ¼



1 for j ¼ 0;

ð58Þ

2 else

and under consideration of the symmetry of A, the expression s1 X s2 X

h

h

Ar1 þj1 ;r2 þj2 C j1 C j2 cosðj1 k h cos hÞ  cosðj2 k h sin hÞ ¼ 0

ð59Þ

j1 ¼0 j2 ¼0

is obtained. In this equation, the only unknown is the numerical wave number kh. With a matrix A computed by any numerical method fulfilling the aforementioned prerequisites, the numerical wave number kh for a given angle of propagation h, wave number k, and field point spacing h can be computed. In order to obtain the numerical solution of Eq. (59), which yields the numerical wave number kh, only one row of the system matrix A belonging to an interior field point has to be assembled. The shape of the domain and the boundary conditions are not taken into account because an interior field point is considered. In the following, the dispersion error e is defined as h

eðh; kÞ ¼

Fig. 2. Influence domains of two associated field points ðxr1 ;r2 Þ and ðxr1 þs1 ;r2 þs2 Þ.

kh  k hðh; kÞ ; kh

ð60Þ

where the wave number k is made dimensionless by multiplication with the field point spacing h. The dimensionless numerical wave number khh only depends on the propagation angle h and the exact wave number k. Thus, for a given method, the dispersion error can be computed for several given wave numbers and propagation angles. As the dispersion error is not dependent on the problem geometry and boundary conditions, it can be applied to assess the employed numerical method objectively. Fig. 3 shows the dispersion error distribution for the standard linear FEM. It should be noted that the distribution is symmetric with respectpto ffiffiffiffiffiffian axis parallel to the kh-axis at h = 45°. Furthermore, if kh 6 12 is not fulfilled, this leads to an evanescent wave [4,9] and thus no reasonable dispersion error can be obtained. Regarding pffiffiffiffiffiffi the linear FEM the limit case of kh ¼ 12 means that one wave is discretized with only two elements, which does not lead to acceptable results. Meshfree methods can yield higher order shape

2228

C. Wenterodt, O. von Estorff / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2223–2236 Table 1 Parameter combinations of the EFGM with polynomial basis. Varied parameter

Parameter values

Number of variations

Weighting function

Cubic spline, Quartic spline, Exponential function Linear, Incomplete quadratic, Quadratic Circle, Square rinfl = 2h, 2.5h, . . . , 4h 4  4,5  5,. . .,10  10

3

Polynomial basis

Fig. 3. Dispersion error for the FEM with respect to the propagation angle h and the dimensionless wave number kh.

functions if the influence domain and polynomial basis are chosen appropriately, but a coarse discretization will still lead to a poor performance. The numerical investigations conducted for this paper showed that kh 6 2 leads to acceptable results. For this reason and in order to keep the number of necessary computations as small as possible, the investigated interval was restricted to kh 6 2. 5. Numerical investigations When the Helmholtz equation is solved by numerical methods as, e.g., the EFGM or the RPIM, the dispersion effect is an appropriate criterion to assess the quality of the respective method. Although the derivation of the numerical wave number given in Section 4 is only valid for a plane wave and a regular field point distribution, it can be assumed that a method that leads to minimal dispersion for this theoretical setup will also perform well in real life applications. For problems which are only marginally dominated by dispersion, other parameters might be better suited, which are, e.g., optimized for certain boundary conditions or problem geometries. Such special cases are not considered here, but might be the topic of future work. In the following, the optimization of the investigated parameters is accomplished by a number of systematic numerical experiments. A general parameter optimization for the considered meshfree methods is not straightforward, since all parameters are nonlinearly dependent on each other and the dependence of the dispersion error on the parameters cannot be determined analytically. A heuristic optimization algorithm as, e.g., a genetic algorithm might be applied. Nevertheless, this kind of optimization method cannot guarantee that an optimal parameter set is found. For this reason, regarding meshfree methods, a parameter optimization by numerical experiments is common practice (cf. [30]). 5.1. Average dispersion error The dispersion error e given in Eq. (60) can be employed to compare different methods with several parameter configurations. For the results presented in this section, countless configurations have been examined, each of which yields an error distribution as shown for the FEM in Fig. 3. As it is impossible to compare all error distributions directly, an appropriate criterion has to be found to directly assess a given error distribution, preferably by just one value. An appropriate measure for the quality of an error distribution is the average dispersion error

e ¼

nkh nh X X 1 jeij j: nh  nkh i¼1 j¼1

Form of influence domain Size of influence domain Number of integration points

3

2 5 7

dispersion error e has been computed for kh 2 [0, 0.05, 0.1, . . . , 2] and h 2 [0°, 1°, 2°, . . . , 45°]. 5.2. Dispersion of the EFGM with polynomial basis In order to find an optimal EFGM in terms of the dispersion error, several parameters are varied. The three weighting functions given in Eqs. (20)–(22) are investigated. In the first step, three different polynomial bases are chosen, namely the linear basis

PT ðxÞ ¼ ½ 1 x y 

ð62Þ

the incomplete quadratic basis

PT ðxÞ ¼ ½ 1 x y xy 

ð63Þ

and the quadratic basis

PT ðxÞ ¼ ½ 1 x y xy x2

y2 :

ð64Þ

Phase-dependent bases are examined in Section 5.3. Furthermore the form and size of the influence domain as well as the number of Gauss integration points per integration cell is varied as listed in Table 1. All parameter variations yield 630 possible combinations. For each parameter combination the distribution of the dispersion error e is computed as described in Section 4. With this error distribution, the average dispersion error e can be obtained from Eq. (61). Thus, each parameter combination yields one average dispersion error and two parameter combinations can be easily compared. As five different parameters have been considered and those parameters appeared to influence each other, only general tendencies of optimal parameters can be obtained. To illustrate this, the 630 parameter combinations are sorted with respect to their average dispersion error e as shown in Fig. 4. Thus, parameter combinations that yield a low error are located on the left hand side while parameter combinations that yield a high error can be found on the right hand side. It can be seen that the error e is in a range of about four orders of magnitude. This kind of representation does not aim to preserve the correlation between a certain parameter

ð61Þ

Here, nh is the number of considered propagation angles and nkh is the number of dimensionless wave numbers that have be taken into account. As mentioned before, values of kh 6 2 lead to reasonable results. Furthermore, for symmetry reasons, the propagation angle has to be considered on the interval 0° 6 h 6 45° only. Thus, for all numerical investigations presented in the following, the average

Fig. 4. Average dispersion error of 630 parameter combinations of the polynomial EFGM.

C. Wenterodt, O. von Estorff / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2223–2236

2229

Fig. 5. Average dispersion error of 630 parameter combinations of the polynomial EFGM; parameter combinations marked according to their polynomial basis and weighting function (same graph, labels and scaling as in Fig. 4).

Fig. 6. Average dispersion error of 280 parameter combinations of the polynomial EFGM (exponential weighting function and quadratic basis excluded).

combination and the respective error. Instead, in the following the graph in Fig. 4 will be employed to mark selected parameter combinations. These combinations can then be easily assessed. In order to find optimal parameters, numerous investigations have been conducted. In the following, only the significant steps which lead to the final parameter recommendations on a direct way are documented. In Fig. 5 the parameter configurations are marked according to their weighting function and polynomial basis. E.g., in the upper left graph all combinations with linear basis and cubic spline weighting function are marked. As the parameter combinations are sorted with respect to the error, the combinations which yield a low error can be found on the left hand side while the combinations leading to a high error are located on the right hand side. In the case of Fig. 5 it can be seen immediately that in most cases the cubic and quartic spline weighting functions lead to better results

Fig. 8. Average dispersion error of the best polynomial EFGM (circular influence domain, linear basis) depending on the number of Gauss integration points.

than the exponential weighting function, as the markers are located on the left side in the graph. Furthermore, although the quadratic basis can lead to reasonable results in some cases, it can also yield a very high average dispersion error. As a consequence, it can be stated that the cubic or quartic spline weighting function and a linear or incomplete quadratic basis should be used. In the following, the exponential function and the quadratic basis are excluded from the further investigations and the average dispersion error of the remaining 280 parameter combinations is shown in Fig. 6. Now, the error is reduced to a range of about three

Fig. 7. Average dispersion error of 280 parameter combinations of the polynomial EFGM; exponential weighting function and quadratic basis excluded; parameter combinations marked according to their influence domain (same graph, labels and scaling as in Fig. 6).

2230

C. Wenterodt, O. von Estorff / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2223–2236

Table 2 Parameter combinations of the EFGM with kh-dependent basis. Varied parameter

Parameter values

Number of variations

Weighting function

Cubic spline, Quartic spline, Exponential function Sine/cosine basis, Extended sine/cosine basis Circle, Square rinfl = 2h, 2.5h, . . . , 4h 4  4, 10  10 ~ h ¼ 0 ; 45

3

kh-Dependent basis Form of influence domain Size of influence domain Number of integration points Angle for basis

2 2 5 2 2

error e is plotted with respect to the number of integration points. In contrast to the aforementioned investigations, parameter combinations with up to 20  20 integration points are added. For the analyzed parameter combinations (cubic and quartic spline weighting function, circular influence domain, rinfl = 3h or 3.5h) the two recommended polynomial basis (linear and incomplete quadratic) did not show any differences. Thus, the results in Fig. 8 are obtained with a linear basis. The cubic spline yields slightly better results than the quartic spline. Although an influence radius of rinfl = 3.5h leads to a lower error than rinfl = 3h, it also causes a stronger dependence of the error on the number of integration points. Surprisingly, the error is not constantly decreasing for a higher number of integration points. An almost constant and thus reliable error can be obtained with at least 8  8 integration points, which does not mean that a lower integration order necessarily leads to worse results. Recommendations that can be drawn from this investigations are summarized in Section 5.4. 5.3. Dispersion of the EFGM with kh-dependent basis It is not possible to completely eliminate the dispersion effect in the two-dimensional case [4]. But, if a kh-dependent basis is chosen as, e.g., sine/cosine basis

PT ðxÞ ¼ ½1 Fig. 9. Average dispersion error of 240 parameter combinations of the khdependent EFGM.

orders of magnitude. In Fig. 7 the parameter configurations which include a certain influence domain (i.e. form and size) are marked. The best results can be obtained with an influence radius of rinfl = 3h or rinfl = 3.5h. A circular influence domain is slightly better than a quadratic one. The remaining parameter configurations can be analyzed with respect to the integration order, i.e. the number of Gauss integration points in each integration cell. In Fig. 8 the average dispersion

cosðkx cos ~h þ ky sin ~hÞ    ~ sinðkx cos h~ þ ky sin hÞ

ð65Þ

or the extended sine/cosine basis

PT ðxÞ ¼ ½1

cosðkx cos ~h þ ky sin ~hÞ . . . ~ ... sinðkx cos h~ þ ky sin hÞ cosðkx sin ~h þ ky cos ~hÞ . . . ~ sinðkx sin h~ þ ky cos hÞ

ð66Þ

the resulting numerical solution is dispersion-free for the given wave number k and the propagation angle ~ h used in the basis. The basis can also be enriched with further propagation angles as

Fig. 10. Average dispersion error of 240 parameter combinations of the kh-dependent EFGM; parameter combinations marked according to their influence domain and weighting function (same graph, labels and scaling as in Fig. 9). Table 3 Optimal parameters for the EFGM. Varied parameter

Weighting function Basis Form of influence domain Size of influence domain Number of integration points ~ Propagation angle h

Optimal parameter values for polynomial EFGM

kh-dependent EFGM

Cubic spline, Quartic spline Linear, Incomplete quadratic Circle rinfl = 3h (3.5h) 5  5, . . . , 8  8 –

Cubic spline Sine/cosine Square rinfl = 4h 44 ~ h ¼ 0

2231

C. Wenterodt, O. von Estorff / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2223–2236

proposed in [26] in order to eliminate the dispersion for this propagation angles as well. In general, this enrichment does not only work for a plane wave but also improves the results of a real life application. Unfortunately, the enriched basis increases the computation time significantly. The investigated 240 parameter variations for the kh-dependent EFGM are listed in Table 2. In addition to the parameters that have been considered for the polynomial EFGM, the propagation angle has to be chosen. Since the numerical wave number is computed as described in Section 4, the shape functions have to be symmetric. Thus, the propagation angle ~ h is chosen to be either 0° or 45°. In Fig. 9 the average dispersion error of the 240 parameter combinations of the kh-dependent EFGM is depicted. As for the polynomial EFGM, the combinations are sorted with respect to the error. If the parameters combinations that include a certain influence domain (size and form) and a certain weighting function are marked, as done in Fig. 10, it can be reasoned that the exponential weighting function is better suited for small influence domains (rinfl = 2h), but fails for larger influence domains which yield much lower errors. For the kh-dependent EFGM it can be stated that the influence domain should be chosen as large as possible. A square influence domain is slightly better suited than a circular one. In most cases, the cubic spline weighting function leads to the lowest error.

Table 4 Parameter combinations of the RPIM. Varied parameter

Parameter values

Number of variations

Basis

No polynomial reproduction Linear Incomplete quadratic Quadratic Sine/cosine (h = 45°) Extended sine/cosine (h = 45°) Square rinfl = 2h, 3h, 4h 44

6

ac = 0.4, 0.8, 1.2. . .10 q = 0.6, 0.7, 0.8

25 3

Wend-RPIM Shape parameter, wc

wc = 0.1, 0.2, 0.3. . .4

40

EXP-RPIM Shape parameter, ec

ec = 0.1, 0.2, 0.3. . .4

40

Form of influence domain Size of influence domain Number of integration points MQ-RPIM Shape parameter, ac Shape parameter, q

1 3 1

Further investigations showed that the kh-dependent EFGM is less dependent on the integration order. In addition, no significant difference between the sine/cosine and the extended sine/cosine basis could be detected.

Fig. 11. Dispersion error distributions of the FEM and three different parameter configurations of the EFGM.

2232

C. Wenterodt, O. von Estorff / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2223–2236

5.4. Optimal EFGM Based on the results from the parameter studies in Sections 5.2 and 5.3, a set of optimal parameters for both the polynomial and the kh-dependent EFGM is given in Table 3. Both methods work well with cubic spline weighting functions. For the polynomial EFGM the quartic spline might be chosen as well. The optimal influence domain of the polynomial EFGM is a circular domain with an influence radius of rinfl = 3h. Although an influence domain with rinfl = 3.5h yields slightly better results, it also leads to an increased computation time and a stronger dependence on the Gauss

Fig. 12. Average dispersion error of 1350 parameter combinations of the MQ-RPIM.

integration order. The polynomial basis should be linear or incomplete quadratic. A complete quadratic basis does not improve the results. For the kh-dependent EFGM a square influence domain is advantageous. Among the investigated influence domains, the influence radius should be chosen as large as possible, i.e. rinfl = 4h. The extended sine/cosine basis does not show any advantages, thus in order to minimize the computation time, the simple sine/cosine basis should applied. The propagation angle ~ h should be adapted to the problem to be solved. Of course, the results will be best if the real propagation angle is similar to the propagation angle used in the basis. With respect to the average dispersion error e, a propagation angle of ~ h ¼ 0 yields slightly better results for the proposed parameter configurations. The dispersion error distributions for the FEM as well as for three parameter configurations of the EFGM and their resulting average dispersion error e are depicted in Fig. 11. It should be noted that the error of the FEM is scaled by a factor of 1000, thus the dispersion error of the polynomial EFGM is only about 0.1% of the FEM dispersion error. The EFGM with sine/cosine basis performs even slightly better. Furthermore, it can be seen that the khdependent EFGM completely eliminates the dispersion if the propagation angle applied in the basis matches the propagation angle of the real problem. On the other hand, the time t that is needed to construct the equation system is significantly increased. A comparison of EFGM, RPIM and FEM with respect to their efficiency is conducted in Section 5.6.

Fig. 13. Average dispersion error of 1350 parameter combinations of the MQ-RPIM; parameter combinations marked according to their influence domain, basis and shape parameter ac (same graph, labels and scaling as in Fig. 12).

C. Wenterodt, O. von Estorff / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2223–2236

2233

Fig. 14. Average dispersion error of the Wend-RPIM with respect to wc for different influence domains and bases.

Fig. 15. Average dispersion error of the EXP-RPIM with respect to ec for different influence domains and bases.

5.5. Dispersion of the RPIM In order to find an optimal RPIM, three different radial basis functions are investigated, namely the multiquadric

(MQ), the Wendland (Wend) and the Gaussian (EXP) basis given in Eqs. (30)–(32), respectively. The results for these radial basis functions are first studied separately and later compared.

2234

C. Wenterodt, O. von Estorff / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2223–2236

Fig. 16. Dispersion error distributions of four different parameter configurations of the RPIM.

As for the EFGM, different polynomial and kh-dependent bases given in Eqs. (62)–(66) are taken into account for the RPIM with polynomial reproduction. Additionally, the simple RPIM without polynomial reproduction is investigated. In contrast to the EFGM, the radial basis functions can be adapted to the problem by different shape parameters. These shape parameters are listed in Table 4. Thus, for the MQ-RPIM 1350 parameter combinations are obtained. For each the Wend-RPIM and the EXP-RPIM, 720 combinations have been examined. For the choice of the shape parameters of the MQ-RPIM and the EXP-RPIM, the recommendations given in [30] had been considered in the first place. Own studies showed that the recommended shape parameters did not lead to optimal results with respect to the dispersion effect. Consequently, the effect of the shape parameters was studied on enlarged intervals as given in Table 4. Further preliminary enquiries showed that a square influence domain is advantageous for the RPIM and that the integration order does not influence the dispersion error. For this reason, the circular influence domain has not been considered in the following and all parameter combinations have been computed with a fixed number of 4  4 Gauss integration points.

5.5.1. MQ-RPIM For the MQ-RPIM the sorted average dispersion errors of all parameter combinations are shown in Fig. 12. In the next step, the parameter combinations belonging to a certain basis, influence domain and shape parameter ac are marked in Fig. 13. Polynomial, as well as kh-dependent bases lead to comparable results, i.e. the performance of the method is more dependent on the choice of the shape parameter ac and the size of the influence domain rinfl than on the choice of the basis. The sine/cosine basis leads only to slightly better results than the polynomial bases. The quadratic basis should be avoided. Even the MQ-RPIM without polynomial reproduction can be considered as a useful method, as the results are in the same range as for the other bases. It is important to note, that reasonable results can only be obtained if the shape parameter ac is chosen appropriately. Furthermore, the optimal ac depends on the chosen size of the influence domain. For rinfl = 3h the shape parameter should be chosen as ac = [3. . .4.5] whereas for rinfl = 4h the optimal interval is ac = [4.5. . .6]. In general, it can be stated that the smaller the influence domain is, the smaller is the optimal ac. Obviously, among the investigated influence domains, the size should be chosen as large as possible i.e. rinfl = 4h. On the other

C. Wenterodt, O. von Estorff / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2223–2236 Table 5 Efficiency of the best meshfree methods in comparison to the FEM and the EFGM with exponential weighting function.

FEM EFGM, linear basis, cubic weighting EFGM, sine cosine basis, h = 0° EFGM, sine cosine basis, h = 45° MQ-RPIM linear basis MQ-RPIM sine cosine basis Wend-RPIM, no polynomial reproduction EXP-RPIM, no polynomial reproduction EFGM, linear basis, exp. weighting [25,26]

e

t (s)

Edisp

3.57  1002 1.56  1005 7.92  1006 1.07  1005 2.53  1005 1.87  1005 4.57  1005 9.86  1006 2.56  1004

0.3 18 24 25 244 243 248 236 15

1.0 38.1 56.3 40.0 1.7 2.4 0.9 4.6 2.8

hand, if the computation time shall be reduced, also a smaller influence domain with rinfl = 3h might be considered. Further investigations showed that the choice of the shape parameter q has a minor influence on the dispersion error. The optimum is obtained for q = [0.6. . .0.8]. 5.5.2. Wend-RPIM For the Wend-RPIM the influence of the shape parameter wc is investigated in Fig. 14 for different bases and influence radii rinfl. On the investigated interval of wc unique minima of the dispersion error can be found, which are marked with circles in Fig. 14. The optimal wc depends mainly on the influence radius and slightly on the basis. As for the MQ-RPIM, the influence radius should be chosen as large as possible. Among the investigated influence radii, rinfl = 4h leads to the best results. A significant advantage of the Wend-RPIM is that very good results can be obtained without polynomial reproduction and thus computation time can be saved. The sine/cosine bases are not superior to the standard polynomial bases. The lowest error of the Wend-RPIM can be obtained if no polynomial reproduction is applied and if the parameters are chosen as rinfl = 4h and wc = 0.7. 5.5.3. EXP-RPIM The EXP-RPIM is investigated in the same way as the WendRPIM. The influence of the shape parameter ec is shown in Fig. 15 for different bases and influence radii rinfl. The minima of the dispersion error are marked with circles. As for the Wend-RPIM, the optimal ec depends mainly on the influence radius and only slightly on the basis. Among the investigated influence radii, rinfl = 4h leads to optimal results. Best results can be obtained if no polynomial reproduction is applied. The quadratic basis should not be used. The sine/cosine bases are not better suited than the standard polynomial bases. The by far lowest error of the EXP-RPIM is obtained if no polynomial reproduction is applied and the parameters are chosen as rinfl = 4h and wc = 2.3. 5.5.4. Comparison of MQ-RPIM, Wend-RPIM and EXP-RPIM The dispersion error distributions of four parameter configurations of the RPIM and their resulting average dispersion error e are depicted in Fig. 16. All results have been obtained using a square influence domain with rinfl = 4h and nqp = 4  4 integration points. Best results can be obtained for the EXP-RPIM without polynomial reproduction. The computation time t for building the equation system is in the same range for all methods, but significantly higher than for the EFGM. A comparison of the efficiency of all investigated meshfree methods and the FEM can be found in Section 5.6.

2235

error of the FEM. However, meshfree methods need an increased computation time to construct the systems of equations. For all investigations, an identical field point distribution has been considered, such that the time for solving the equation system in a real life application should be almost the same for all meshfree methods and the FEM. In order to compare the efficiency of meshfree methods, the average dispersion errors and computation times of the best parameter combinations shown in Figs. 11 and 16 are listed in Table 5. Additionally, an efficiency criterion Edisp is established, which is defined as

Edisp ¼

eFEM  t FEM ; e  t

ð67Þ

where eFEM and tFEM are the dispersion error and the computation time of the FEM, respectively. Thus, the efficiency of the FEM is 1, whereas an Edisp > 1 indicates a method which is more efficient than the FEM. The efficiency of the EFGM with an exponential weighting, which has been used earlier by Suleau et al. [25,26] for acoustic problems is also shown. It gets clear that, except for the WendRPIM, meshfree methods are more efficient than the FEM, although this effect is less significant for the RPIM. Especially the EFGM with sine/cosine basis and h = 0° leads to excellent results in conjunction with a reasonable computation time. The EFGM with exponential weighting given in [25,26] is in the same range as the RPIM but cannot keep up with the optimized EFGMs discovered in this work. As a conclusion, among the investigated methods, the EFGM with cubic spline weighting function and linear or sine/cosine basis is most efficient with respect to the dispersion effect. In comparison to the EFGM used in [25,26], the dispersion error could be reduced by more than one order of magnitude without a significant increase of the computation time. 6. Conclusions Two meshfree methods, namely the element free Galerkin method (EFGM) and the radial point interpolation method (RPIM) have been optimized with respect to the dispersion effect. Numerous parameters have been varied, as e.g. weighting functions, polynomial bases, radial basis functions, support domains, integration orders, and shape parameters. For the polynomial EFGM a cubic spline weighting function, a linear basis, a circular influence domain with influence radius rinfl = 3h and about 5  5 integration points is a good choice. For the kh-dependent EFGM, the cubic spline weighting function, a square influence domain with radius rinfl = 4h and 4  4 integration points are well suited. Regarding the RPIM, the EXP-RPIM without polynomial reproduction performs best. Here, a square influence domain with radius rinfl = 4h and 4  4 integration points should be applied. These optimal meshfree methods yield a significantly lower dispersion error than the FEM but also require an increased computation time. Nevertheless, the efficiency of almost all considered meshfree methods is higher than the efficiency of the standard FEM. It can thus be stated that meshfree methods are a valuable alternative to the FEM, especially in the higher frequency range. In the future, the optimal meshfree methods suggested in this contribution have to be investigated with respect to acoustic problems in order to prove their utilisability for real life applications. References

5.6. Efficiency of meshfree methods For identical discretizations, the dispersion error of meshfree methods is drastically reduced in comparison to the dispersion

[1] T. Belytschko, Y.Y. Lu, L. Gu, Element-free Galerkin methods, Int. J. Numer. Methods Engrg. 37 (2) (1994) 229–256. [2] T. Belytschko, Y.Y. Lu, L. Gu, L. Tabbara, Element-free Galerkin methods for static and dynamic fracture, Int. J. Solids Struct. 32 (17–18) (1995) 2547–2570.

2236

C. Wenterodt, O. von Estorff / Comput. Methods Appl. Mech. Engrg. 200 (2011) 2223–2236

[3] P. Bouillard, S. Suleau, Element-free Galerkin solutions for Helmholtz problems: formulation and numerical assessment of the pollution effect, Comput. Methods Appl. Mech. Engrg. 162 (1) (1998) 317–335. [4] A. Deraemaeker, I. Babuška, P. Bouillard, Dispersion and pollution of the FEM solution for the Helmholtz equation in one, two and three dimensions, Int. J. Numer. Methods Engrg. 46 (4) (1999) 471–499. [5] G.E. Fasshauer, Solving differential equations with radial basis functions: multilevel methods and smoothing, Adv. Comput. Math. 11 (2–3) (1999) 139– 159. [6] Y.T. Gu, Meshfree methods and their comparisons, Int. J. Comput. Methods 2 (4) (2005) 477–515. [7] Z.C. He, G.R. Liu, Z.H. Zhong, S.C. Wu, G.Y. Zhang, A.G. Cheng, An edge-based smoothed finite element method (es-fem) for analyzing three-dimensional acoustic problems, Comput. Methods Appl. Mech. Engrg. 199 (1–4) (2009) 20– 33. [8] T. Huttunen, P. Gamallo, R.J. Astley, Comparison of two wave element methods for the Helmholtz problem, Commun. Numer. Methods Engrg. (2008) (Early View). [9] F. Ihlenburg, I. Babuška, Dispersion analysis and error estimation of Galerkin finite element methods for the Helmholtz equation, Int. J. Numer. Methods Engrg. 38 (22) (1995) 3745–3774. [10] F. Ihlenburg, I. Babuška, Finite element solution of the Helmholtz equation with high wave number Part I: The h-version of the FEM, Comput. Math. Appl. 30 (9) (1995) 9–37. [11] F. Ihlenburg, I. Babuška, Finite element solution of the Helmholtz equation with high wave number Part II: The h-p version of the FEM, SIAM J. Numer. Anal. 34 (1) (1997) 315–358. [12] F. Ihlenburg, I. Babuška, Reliability of finite element methods for the numerical computation of waves, Adv. Engrg. Software 28 (1997) 417–424. [13] Pei-Lin Jiang, Shu-Qing Li, Chi Hou Chan, Analysis of elliptical waveguides by a meshless collocation method with the Wendland radial basis functions, Microwave Opt. Technol. Lett. 32 (2) (2001) 162–165. [14] O. Laghrouche, P. Bettess, Short wave modelling using special finite elements, J. Comput. Acoust. 8 (1) (2000) 189–210. [15] G.R. Liu, Mesh Free Methods: Moving Beyond the Finite Element Method, CRC Press, Boca Raton, 2003. [16] G.R. Liu, Y.T. Gu, An Introduction to Meshfree Methods and their Programming, Springer, Netherlands, 2005. [17] G.R. Liu, Y.T. Gu, K.Y. Dai, Assessment and applications of point interpolation methods for computational mechanics, Int. J. Numer. Methods Engrg. 59 (2004) 1373–1397. [18] G.R. Liu, T. Nguyen-Thoi, K.Y. Lam, An edge-based smoothed finite element method (es-fem) for static, free and forced vibration analyses of solids, J. Sound Vibr. 320 (4–5) (2009) 1100–1130.

[19] G.R. Liu, G.Y. Zhang, Y.T. Gu, Y.Y. Wang, A meshfree radial point interpolation method (RPIM) for three-dimensional solids, Comput. Mech. 36 (2005) 421– 430. [20] J.M. Melenk, I. Babuška, The partition of unity finite element method: basic theory and applications, Comput. Methods Appl. Mech. Engrg. 139 (1–4) (1996) 289–314. [21] R. Mullen, T. Belytschko, Dispersion analysis of finite element semidiscretizations of the two-dimensional wave equation, Int. J. Numer. Methods Engrg. 18 (1982) 11–29. [22] A.A. Oberai, P.M. Pinsky, A residual-based finite element method for the Helmholtz equation, Int. J. Numer. Methods Engrg. 49 (3) (2000) 399– 419. [23] P. Ortiz, E. Sanchez, An improved partition of unity finite element model for diffraction problems, Int. J. Numer. Methods Engrg. 50 (12) (2001) 2727– 2740. [24] S. Petersen, D. Dreyer, O. von Estorff, Assessment of finite and spectral element shape functions for efficient iterative simulations of interior acoustics, Comput. Methods Appl. Mech. Engrg. 195 (44–47) (2006) 6463–6478. [25] S. Suleau, P. Bouillard, One-dimensional dispersion analysis for the elementfree Galerkin method for the Helmholtz equation, Int. J. Numer. Methods Engrg. 47 (6) (2000) 1169–1188. [26] S. Suleau, A. Deraemaeker, P. Bouillard, Dispersion and pollution of meshless solutions for the Helmholtz equation, Comput. Methods Appl. Mech. Engrg. 190 (5) (2000) 639–657. [27] L.L. Thompson, P.M. Pinsky, A Galerkin least-squares finite element method for the two-dimensional Helmholtz equation, Int. J. Numer. Methods Engrg. 38 (3) (1995) 371–397. [28] O. von Estorff, J. Quan, Direct coupling of EFGM–FEM and EFGM–BEM for dynamic soil–structure interactions, Int. J. Comput. Methods 2 (4) (2005) 627– 644. [29] J.G. Wang, Bingyin Zhang, T. Nogami, Wave-induced seabed response analysis by radial point interpolation meshless method, Ocean Engrg. 31 (1) (2004) 21– 42. [30] J.G. Wang, G.R. Liu, On the optimal shape parameters of radial basis functions used for 2-D meshless methods, Comput. Methods Appl. Mech. Engrg. 191 (23–24) (2002) 2611–2630. [31] J.G. Wang, G.R. Liu, A point interpolation meshless method based on radial basis functions, Int. J. Numer. Methods Engrg. 54 (11) (2002) 1623–1648. [32] H. Wendland, Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree, Adv. Comput. Math. 4 (1) (1995) 389–396. [33] C. Wenterodt, O. von Estorff, Dispersion analysis of the meshfree radial point interpolation method for the Helmholtz equation, Int. J. Numer. Methods Engrg. 77 (12) (2009) 1670–1689.