Coons-patch macroelements in two-dimensional eigenvalue and scalar wave propagation problems

Computers and Structures 82 (2004) 383–395 www.elsevier.com/locate/compstruc

Coons-patch macroelements in two-dimensional eigenvalue and scalar wave propagation problems Christopher G. Provatidis

*

Department of Mechanical Engineering, National Technical University of Athens, 9 Heroon Polytechniou Avenue, Zografou Campus, 157 73 Athens, Greece Received 1 November 2002; accepted 10 October 2003

Abstract Recently, a global functional set has been proposed for the construction of large finite elements, here called ‘‘Coonspatch macroelements’’, with degrees of freedom at the boundaries only. The background of the method is ‘‘Coons interpolation formula’’; a mathematical formula capable of describing CAD surfaces initially applied to automobile applications. So far, these macroelements were found to be accurate in potential and plane elasticity problems. Based on these encouraging results, this paper continues the investigation on the universality which Coons interpolation formula offers in the development of a class of macroelements including members of the well-known serendipity family. In this context, it is shown in a systematic way that it is not necessary to use Lagrange polynomials to interpolate the potential along each side of the macroelement but one may choose piecewise-linear and/or cubic B-splines interpolants. Moreover, the paper investigates the performance of Coons-patch macroelements in eigenvalue and scalar wave propagation problems by implementing standard time-integration schemes.
2003 Elsevier Ltd. All rights reserved. Keywords: Finite element; Macroelement; Global approximation; Coons interpolation; Scalar wave propagation; Computer-aideddesign (CAD)

1. Introduction Historically, before the application of the conventional finite element method [1,2,43], which may be probably characterized as a local approximation method, early numerical methods used to treat the entire problem domain as a global approximation procedure [3,4]. In the sequence, apart from the use of small-size isoparametric finite elements, the development of ‘‘large’’ elements aiming at reducing the mesh generation workload, the total number of degrees of freedom, as well as the computational effort, has kept researchers busy for a long time [5–7]. According to Zienkiewicz [2, pp. 155– 159], besides the well-known bilinear and quadratic elements, cubic and quartic elements were originally

*

Tel.: +30-210-772-1520; fax: +30-210-772-2347. E-mail address: [email protected] (C.G. Provatidis).

derived by inspection. Progression to yet higher members, which constitute the ÔSerendipity’ family, is difficult and requires some ingenuity. A systematic way of generating the Ôserendipity’ shape functions was first introduced by Zienkiewicz et al. [8] but a simpler formulation is that of Taylor [9]. Moreover, Taig [10] appears to be the first who mentioned the concept of using shape functions for establishing curvilinear co-ordinates in the context of finite element analysis. In his first application basic linear quadrilateral relations were established, while Irons [11] generalized the idea of arbitrarily noded elements. It is interesting that in Ref. [2, pp. 181–182] there is a special remark that ‘‘. . .quite independently the exercises of devising various practical methods of generating curved surfaces for purposes of engineering design led to the establishment of similar definitions by Coons [12], and indeed today the subjects of surface definitions and analysis are drawing closer together due to this

0045-7949/$ - see front matter
2003 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2003.10.012

384

C.G. Provatidis / Computers and Structures 82 (2004) 383–395

activity . . .’’. Despite the latter important remark and apart from some interesting special element families [13– 15], no significant relevant progress seems to appear in the open literature. The Coons’ interpolation method [12], which consists the mathematical background of the proposed global interpolation method, was used initially for the description of regular CAD-surfaces for automotive industry and later for the mesh generation in structured four-sided curvilinear patches [16]. Gordon [13] as well as Gordon and Hall [14] proposed some early theoretical thoughts about the use of Coons interpolation in the construction of (reduced size) transfinite finite elements. Later, Coons’ interpolation was applied to twodimensional problems in conjunction with Lagrange and Hermite interpolation functions by El-Zafrany and Cookson [15] allowing a small number of degrees of freedom per element, as well as in some families of finite element of plates and shells [17]. The above ideas were further extended by applying cubic B-splines and/or piecewise-linear interpolation along each of the four sides of a Coons patch, a procedure that led to large (global) cardinal shape functions within either the whole domain or a significantly large portion of it. Therefore, a series of papers have been recently published on applying large finite elements (with degrees of freedom appearing only at the element boundaries) to potential [18–20] as well as to elasticity problems [21–23]. Clearly, both cubic B-spline and piecewise-linear interpolations have been well estab-

lished before this work [18–23], but not in a systematic way. In addition, a generalization of Coons interpolation to include derivatives has been proposed [24]. This paper consists of two parts. In the first part, it is shown that Coons formula is not only capable of deriving the common serendipity elements (bilinear, quadratic, cubic) but it also offers a systematic way of generating global cardinal shape functions for larger elements with any different number of nodes along each side. In other words, Coons formula achieves a global approximation. In the second part, the theory is applied to a 24-noded Coons-patch macroelement and its performance to hyperbolic problems (two-dimensional eigenvalue and scalar wave propagation) is investigated. For demonstration purposes, usual explicit and implicit time-integration techniques are applied to both Coonspatch macroelements and conventional bilinear finite elements.

2. Development of Coons-patch macroelements 2.1. General theory Two-dimensional Coons-patch macroelements treat the entire problem domain, or a large portion of that, as a four-sided patch ABCD on the (x; y)-plane. This patch is mapped to a reference patch (n; g), where the normalized curvilinear co-ordinates vary between 0 and 1 (0 6 n; g 6 1) as shown in Fig. 1. According to Coons’

η

y

1

η

E1(η)

0 1

(η =1)

D

Eo(η)

C

D

C A (ξ =η = 0)

1

η

P

⇔

P

(ξ = 1)

Eo(ξ)

A

B ξ

0 E1(ξ) B

ξ

1

O

ξ

x

Fig. 1. Reference patch (n; g) surrounded by four sides (AB: g ¼ 0, BC: n ¼ 1, CD: g ¼ 1 and DA: n ¼ 0) with q1 ¼ 9, q2 ¼ 10, q3 ¼ 6 and q4 ¼ 6 nodes, respectively.

C.G. Provatidis / Computers and Structures 82 (2004) 383–395

interpolation formula [12], which can be also found in standard CAD textbooks [25,26], each point xðn; gÞ ¼ fxðn; gÞ; yðn; gÞgT along the patch can be approximated by its boundaries ðxðn; 0Þ; xðn; 1Þ; xð0; gÞ; xð1; gÞÞ as follows: xðn; gÞ ¼ E0 ðnÞxð0; gÞ þ E1 ðnÞxð1; gÞ þ E0 ðgÞxðn; 0Þ

385

with uj ðtÞ denoting time-dependent nodal potentials appearing only at the boundaries of the macroelement and qe the total number of the nodes along that. Following the previously mentioned procedure, one can derive the following expressions: (i) Corner nodes

þ E1 ðgÞxðn; 1Þ  E0 ðnÞE0 ðgÞxð0; 0Þ AB NA ðn; gÞ ¼ E0 ðnÞBDA q4 ðgÞ þ E0 ðgÞB1 ðnÞ  E0 ðnÞE0 ðgÞ

 E1 ðnÞE0 ðgÞxð1; 0Þ  E0 ðnÞE1 ðgÞxð0; 1Þ  E1 ðnÞE1 ðgÞxð1; 1Þ

ð1Þ

where the blending functions can be chosen, for example, to be linear as follows: E0 ðnÞ ¼ 1  n;

E1 ðnÞ ¼ n

E0 ðgÞ ¼ 1  g;

E1 ðgÞ ¼ g

AB NB ðn; gÞ ¼ E1 ðnÞBBC 1 ðgÞ þ E0 ðgÞBq1 ðnÞ  E1 ðnÞE0 ðgÞ CD NC ðn; gÞ ¼ E1 ðnÞBBC q2 ðgÞ þ E1 ðgÞB1 ðnÞ  E1 ðnÞE1 ðgÞ CD ND ðn; gÞ ¼ E0 ðnÞBDA 1 ðgÞ þ E1 ðgÞBq3 ðnÞ  E0 ðnÞE1 ðgÞ

ð7Þ

ð2Þ (ii) Interior nodes to AB (local numbering)

In the sequence, the idea of isoparametric elements is applied to Eq. (1) for the interpolation of the potential uðn; gÞ within the patch, as follows:

Nj ðn; gÞ ¼ E0 ðgÞBAB j ðgÞ;

2 6 j 6 q1  1

ð8Þ

(iii) Interior nodes to BC (local numbering)

uðn; gÞ ¼ E0 ðnÞuð0; gÞ þ E1 ðnÞuð1; gÞ þ E0 ðgÞuðn; 0Þ Nj ðn; gÞ ¼ E1 ðnÞBBC j ðgÞ;

þ E1 ðgÞuðn; 1Þ  E0 ðnÞE0 ðgÞuð0; 0Þ  E1 ðnÞE0 ðgÞuð1; 0Þ  E0 ðnÞE1 ðgÞuð0; 1Þ  E1 ðnÞE1 ðgÞuð1; 1Þ

uðn; 0Þ ¼

q1 X

BAB j ðnÞuðnj ; 0Þ

j¼1

uð1; gÞ ¼

Side BC: Side CD:

uðn; 1Þ ¼

q2 X

BBC j ðgÞuð1; gj Þ

j¼1 q3 X

ð5Þ BCD j ðnÞuðnj ; 1Þ

j¼1

Side DA:

uð0; gÞ ¼

q4 X

BDA j ðgÞuð0; gj Þ

j¼1

and Eq. (3) is then collocated to all boundary nodes of the reference macroelement, then the global cardinal shape functions Ni ðn; gÞ within the Coons-patch can be constructed and the solution uðn; gÞ is approximated by: uðn; gÞ ¼

qe X j¼1

Nj ðn; gÞuj ðtÞ

Nj ðn; gÞ ¼ E1 ðgÞBCD j ðnÞ;

2 6 j 6 q3  1

ð10Þ

(v) Interior nodes to DA (local numbering) Nj ðn; gÞ ¼ E0 ðnÞBDA j ðgÞ;

2 6 j 6 q4  1

ð11Þ

ð4Þ

If the boundary values uðn; 0Þ, uðn; 1Þ, uð0; gÞ and uð1; gÞ ^ in (3) are interpolated by any set of trial functions Bj ðnÞ ^ (n is either n or g; the upper index in Bj below corresponds to the relevant side): Side AB:

ð9Þ

(iv) Interior nodes to CD (local numbering) ð3Þ

Let us assume that the sides AB, BC, CD and DA include q1 , q2 , q3 and q4 nodes, respectively. Then, the total number of nodes along the boundary of the patch becomes: qe ¼ q1 þ q2 þ q3 þ q4  4

2 6 j 6 q2  1

ð6Þ

^ inðnÞ In principle, the set of trial functions BSIDE j volved in Eqs. (7)–(11) can be arbitrarily chosen but they should be linearly independent and complete. Typical cases are Lagrange polynomials [13,14] (see also Ref. [27, pp. 57–59]), piecewise-linear interpolation along a side [28,20,23], cubic B-splines [18], et cetera. Remarks. (a) It is here clarified that it is not necessary to preserve the same interpolation along each of the four sides. Therefore, in principle, the side AB could be approximated using Lagrangian polynomials, the side BC using cubic B-splines, the side CD using piecewise-linear interpolation and the side DA using any of the previous or even a different type of interpolant such as piecewise quadratic interpolation. (b) It is noted that Ref. [23] presents the cubic Bsplines interpolation in the form of Eqs. (7)–(11) and the piecewise-linear interpolation in a different form, which is also found in [20]. However, in this paper it is presented for the first time that the functions Bj appearing in Eqs. (7)–(11) can be of any type, including the particular case of piecewise-linear interpolation (‘‘hat’’-type functions Bj ).

386

C.G. Provatidis / Computers and Structures 82 (2004) 383–395

2.2. Special case: cubic element

uðx; 0Þ ¼ u0 ðxÞ

Recently, it has been clearly shown that in case of linear blending functions and linear (respectively, quadratic variation along each side), Eq. (3) leads to the common four-node bilinear (respectively, eight-node quadratic) finite elements [19]. Here, the procedure is extended to a twelve-node macroelement in which the trial functions are chosen as Lagrange polynomials Lnk giving unity at nk and passing through n ¼ 3 points. It is shown in Appendix A that this macroelement coincides with the well-known cubic element of serendipity family.

vðx; 0Þ ¼

2.3. Arbitrary-noded elements As it was previously mentioned, apart from Lagrange polynomials that tend to produce parasitic phenomena (undesired oscillations), particular choices of trial functions may include cardinal cubic B-splines, piecewiselinear shape functions, and so on. Typical global shape functions for a rectangular macroelement of twenty-four nodes (eight in x- and four in y-direction), using both piecewise-linear and cubic B-splines, are illustrated in Fig. 2. In this figure one can notice that these global shape functions are cardinal ([1-0]-type). Obviously, such shape functions expand the family of conventional serendipity finite elements, in the sense that any type of interpolation functions and any number of nodes along each side of a patch can be used. Apparently, for a large number of nodes it is not necessary to derive closed form expressions of the global shape functions but it is adequate to numerically calculate their values and derivatives at the integration (Gauss) points.

 ouðx; tÞ  ¼ v0 ðtÞ ot t¼0

ð13Þ

for every point x in X, at t ¼ 0; and boundary conditions over the C boundary: u¼ uðtÞ on C1 q¼ qðtÞ on C2

ð14Þ

In Eq. (14), C1 and C2 are parts of the boundary C, where the potentials and the surface tractions ðq ¼ ou=onÞ are prescribed, respectively (n ¼ normal to C unit outward vector). 3.2. Macroelement analysis The usual Galerkin procedure [2,43] is applied to Eq. (12) and finally leads to the following matrix formulation: ½M f€ ug þ ½K fug ¼ ff ðtÞg

ð15Þ

where the Ômass’ ½M and Ôstiffness’ ½K matrices are given in terms of the global shape functions by: Z Z ½M ij ¼ ð1=c2 Þ Ni Nj dX ¼ ð1=c2 Þ Ni Nj dx dy ð16Þ X

½K ij ¼

Z

rNi rNj dX ¼

X

X

Z

rNi rNj dx dy

ð17Þ

X

and the vector of the external force ff g by: ff ðtÞgi ¼

Z

Ni ru dC

ð18Þ

C

3. Application to scalar wave propagation 3.1. Governing equation In the absence of body sources, the differential equation governing scalar wave propagation through an elastic, isotropic, homogeneous body is given by: ð1=c2 Þ

o2 uðx; tÞ  r2 uðx; tÞ ¼ 0 ot2

ð12Þ

where u denotes the potential (e.g. longitudinal displacement), x the position vector, t the time, c the ve-ffi pffiffiffiffiffiffiffiffi locity of the wave propagation given by E=q (E ¼ elastic modulus and q ¼ mass density) and r the Nabla operator. In order to find a particular solution of Eq. (12) within the two-dimensional domain XðxÞ, it is necessary to specify initial conditions:

3.2.1. Eigenvalue analysis In the context of the proposed Coons-patch macroelements, the eigenvalue analysis in the body under the boundary conditions of Eq. (14), is performed as usual, that is by determining the roots of the determinant: k½K  x2 ½M k ¼ 0

ð19Þ

where ½K and ½M are the global matrices given by Eqs. (16) and (17). The above procedure is purely algebraic and this fact emphasizes the superiority of the proposed macroelements in contrast to the non-algebraic boundary element solutions due to their frequency-dependent fundamental solutions (kernels) [29,30]. Of course, the Dual Reciprocity Method (DRM) [31] is an exception to the previous remark, but this technique suffers, at least in this

C.G. Provatidis / Computers and Structures 82 (2004) 383–395

387

Fig. 2. Global shape functions for a twenty-four noded rectangular Coons-patch macroelement, of dimensions 2 · 1, for a corner-node using (a) piecewise-linear interpolation and (b) cubic B-splines interpolation.

specific class of problems, where high frequencies are involved. Clearly, DRM requires internal degrees of freedom [32], a topic that has been previously discussed by several authors [33–38].

3.2.2. Time response As usual, Eq. (15) is solved using explicit (e.g. central difference) and implicit (e.g. Houbolt, h-Wilson, Newmark) techniques. The mass and stiffness matrices are calculated once and the sequential time-integration procedure is applied blindly. This fact differentiates the proposed method in contrast to the BEM that requires complicated time-dependent fundamental solutions and causality conditions [39,40].

3.3. Alternative formulations Coons-patch macroelements were applied using two alternative trial functions Bj :

• Model 1: Piecewise-linear global shape functions (Figs. 2a and 3a). • Model 2: Cubic B-splines (Figs. 2b and 3a). In addition to the above formulations, a third model was used in order to implement a conventional FEM solution with four-node bilinear finite elements. The corresponding mesh is shown in Fig. 3(b) and was generated so as to have the same number of boundary nodes with those of the Coons-patch macroelement shown in Fig. 3(a). With respect to the above models, it is strongly emphasized that Model 1 and the conventional bilinear FEM model are completely different. This difference is due to the fact that a nodal piecewise-linear global shape function is active within a large area while a nodal shape function of a conventional finite element is active only within a significantly smaller area. In fact, the corner node (x ¼ 0; y ¼ 0) shown in Fig. 2a carries a piecewiselinear function (‘‘hat’’-type) that is active within the entire area of two strips (0 6 x 6 0:25, 0 6 y 6 1:00) and (0 6 x 6 2:00, 0 6 y 6 0:25), equivalently (0 6 n 6 0:125,

388

C.G. Provatidis / Computers and Structures 82 (2004) 383–395 y q=0

D

G

C

q = Eqo

H

(a)

b

E A

F q=0

x B

a

(b)

Fig. 3. Wave propagation in elastic rod: Domain discretization using (a) one Coons-patch macroelement (24 nodes) and (b) thirty-two conventional bilinear finite elements (45 nodes).

0 6 g 6 1:00) and (0 6 n 6 1:00, 0 6 g 6 0:25). On the contrary, a conventional finite element at this corner would be active only in the area (0 6 x 6 0:25, 0 6 y 6 0:25). Moreover, according to Eq. (8) applied to Model 1, for any node ‘‘i’’ between the corners A and B (g ¼ 0), the corresponding global shape function Ni is active within a strip that is perpendicular to the n-axis and it is defined by ni1 6 ni 6 niþ1 . On the contrary, as far as Model 2 is concerned the global shape function carried by any nodal point is active within the entire patch. 3.4. Estimation of matrices With respect to the domain integrals in Eqs. (16) and (17), numerous numerical experiments have shown that instead of the usual Gauss integration within the entire domain an alternative scheme should be adopted. Generally, it is suggested to apply Gaussian quadrature using a few points (e.g. 2 · 2 integration points) in more integration cells, which are defined in the reference macroelement by n ¼ constant and g ¼ constant families of the introduced DOFs. Obviously, the so-produced matrices (½M ; ½K ) are both symmetric and fully occupied. In more details, let us assume that the discretization of the quadrilateral structure consists of Nn and Ng subdivisions in the directions n and g, respectively. Then, the Coons-patch macroelement includes NMACRO ¼ 2ðNn þ Ng Þ nodes while the corresponding conventional FEM solution (with the same number of boundary nodes) includes NFEM ¼ ðNn þ 1ÞðNg þ 1Þ nodes. For each partic-

ular case, the number of the integration points involved in the estimation of any matrix is given below: (1) The conventional FEM uses (Nn Ng ) elements and 2 · 2 ¼ 4 Gauss points per element. So, a total number of Nintegration;points ¼ 4Nn Ng integration points are required. When ignoring the symmetry of the matrix under consideration (for simplicity purposes), 4 · 4 ¼ 16 matrix terms for each finite element using 4 Gauss points should be calculated (64 numerical operations in total). So, the total number of numerical operations for all fiFEM nite elements involved becomes Noperations ¼ 64Nn Ng . (2) The cubic B-splines based Coons-patch macroelement (Model 2) could use Ncell ¼ Nn Ng integration cells and, for example, 2 · 2 ¼ 4 Gauss points per cell (totally: 4Ncell integration points). In other words, for the previous choice the total number of integration points coincides with those of the corresponding FEM solution. However, in this case each integration point is related to all nodes of the macroelement. So, since each global 2 matrix consists of NMACRO matrix elements, the total MACRO;2 number of numerical operations becomes Noperations ¼ ðN þN Þ2

g n 2 NMACRO 4Ncell ¼ 16ðNn þ Ng Þ2 Nn Ng , which is 4 times larger than that of the FEM. Of course, it is possible to choose integration cells of double size so that the number of numerical operations reduces to onefourth. (3) The piecewise-linear based Coons-patch macroelement (Model 1) requires less numerical operations than the above mentioned Model 2 because, as explained earlier at the end of Section 3.3, each boundary node influences only an adjacent strip being perpendicular to the side where this node belongs to. Here, at least two integration schemes are possible. The first scheme is to preserve the previously mentioned integrations cells and perform a 2 · 2 quadrature only on those cells where non-vanishing products of global shape functions exist. For the particular case of the numerical example in the next section (Nn ¼ 8, Ng ¼ 4), the number of numerical operations reduces to more than one-fifth than that of Model 2. The second scheme takes advantage from the fact that the global shape function is piecewise-linear. For example, by applying Eq. (8) for a node Ôi’ along AB, the global shape function becomes Ni ¼ ð1 

nni gÞð1  niþ1 Þ, ni 6 n 6 niþ1 . So, only a 2 · 2 quadrature is ni

sufficient for each of the four cells appearing in the overlapping area of the entire strip (ni1 6 ni 6 niþ1 ) with the strip gj1 6 gj 6 gjþ1 that corresponds, for example, to the Ôj’th node along BC. Instead, when considering the node Ôj’ to be along AB with ji  jj > 1, then the matrix elements kij and mij vanish while, when ji  jj ¼ 1, only one cell is required (instead of the aforementioned four). By taking a mean average of two integration cells (each of 2 · 2 quadrature: totally 8 points) throughout, the number of numerical operations MACRO;1 2 is approximated as Noperations ffi 8 NMACRO ¼ 32ðNn þ

C.G. Provatidis / Computers and Structures 82 (2004) 383–395

Ng Þ2 . For the previously mentioned case (Nn ¼ 8, Ng ¼ 4), one obtains that the proposed Coons-patch MACRO;1 macroelement (Model 1) requires Noperations ffi 4608 (the exact value was found to be 4640) numerical operations; that is more than twice compared to those required in FEM the conventional FEM (Noperations ¼ 2048). Despite the above shortcoming of increased computer effort concerning the numerical integration for the estimation of the global stiffness and mass matrices, it should be noted that the total computational cost is reduced in three ways. First, the data preparation cost is minimal because only the geometrical data of the boundary should be given. Second, an additional gain is due to the smaller size of the matrices when compared to the FEM (NMACRO < NFEM ). Third, as it will be shown in the subsequent numerical study, the time step can be chosen to be larger than that of the conventional FEM solution. Finally, for the particular case of rectangular Coons-patch macroelements, analytical expressions can be also derived for both stiffness and mass matrices. In this case the efficiency of the proposed Coons-patch macroelement approach dramatically improves.

4. Numerical results In order to demonstrate the efficiency of the proposed Coons-patch macroelement approach, a study was carried out concerning the propagation of longitudinal waves in an elastic rod fixed at one of its extremities ðx ¼ 0Þ and subjected to a Heaviside type loading q ¼ Eq0 H ðt  0Þ [N/m2 ] on the other one ðx ¼ aÞ. The load, the geometrical characteristics and boundary conditions of this problem are shown in Fig. 3. This problem was chosen because of its degree of difficulty. In this problem, which is actually one-dimensional, all of the odd nodes contribute to potentials (displacements), and tractions are strongly influenced by higher modes [41]. In order to elucidate the proposed method, accuracy tests were performed on the calculation of the lowest eigenvalues before the time response analysis. Then, explicit (central difference) and implicit (h-Wilson, Newmark) schemes were applied.

389

4.1. Mode shapes and eigenvalues The exact solution of the full two-dimensional problem is characterized by mode shapes [42]:  npx   mpy  ^ un;m ¼ sin cos ð20Þ 2a b while the related eigenvalues are given by:  2 n m2 þ ; x2n;m ¼ p2 c2 4a2 b2 n ¼ 1; 3; 5; . . . m ¼ 0; 1; 2; . . .

ð21Þ

Numerical results were obtained using (a) a mesh of 32 four-noded bilinear elements (totally 45 nodes) and (b) a 24-noded macroelement shown in Fig. 3. For both cases the number and the location of the boundary nodes is the same (twenty-four). In the beginning, in order to give insight on the accuracy of the proposed macroelements, the lowest eigenvalues of the domain were calculated using a consistent mass matrix. It is clearly shown in Table 1 that there is no practical difference between the results obtained using the Coons-patch macroelements (Model 1) and the conventional four-noded finite elements. On the contrary, cubic B-splines interpolation (Model 2) is superior. Generally speaking, this conclusion is not really a new finding; a similar superiority has been observed on a free-free acoustical cavity of different dimensions [18]. Finally, it was found that when the domain is divided into two and four equal sub-regions (Fig. 4) the two first eigenvalues do not change at all while the third one improves by only 0.2%. Now, a lumped mass is used for the 24-node Coonspatch macroelement. As usual, it was calculated by adding the elements along each line of the consistent mass matrix. It was also found that only the first calculated eigenvalue is adequately accurate while the rest ones are smaller than the exact values. In addition, two negative eigenvalues appear! This is due to the fact that the lumped mass (at all four corner nodes) equals to )1.25, while the total mass is 2.00. All other lumped masses along the long and the short side equal to 0.125 and 0.25, respectively [Total mass ¼ 4ð0:3125Þ

Table 1 Calculated eigenvalues using consistent mass formulation for (a) one Coons-patch macroelement (two different interpolations: Model 1 and Model 2) and (b) thirty-two conventional four-node finite elements Mode shapes n, m

1, 0 3, 0 0, 1

Exact eigenvalues (s2 )

0.6168 5.5516 10.4865

Errors (in %) of calculated eigenvalues Coons-patch macroelement

Conventional FEM

B-splines Interpolation (Model 2)

Piecewise-linear interpolation (Model 1)

0.02 0.08 0.99

0.32 2.93 5.14

0.32 2.93 4.95

390

C.G. Provatidis / Computers and Structures 82 (2004) 383–395

(a)

(b)

Fig. 4. Subdivision of the problem area into (a) two and (b) four equal sub-regions with 27 and 33 nodes, respectively.

þ14 0:125 þ 6 0:25 ¼ 2:00]. It was also found that when the domain is divided into two and four equal subregions (Fig. 4) the results significantly improve, as shown in Table 2. As a result of the above findings, only the consistent mass formulation was adopted and applied in all sequential time-response results.

4.2. Time response In the sequence, the investigation will focus on conventional time-integration schemes including explicit (central difference method) and implicit (h-Wilson, Newmark) techniques, well known from the common FEM praxis. All relevant algorithms were based on the textbook of Bathe [43] (Chapter 9: Tables 9.1, 9.3 and 9.4, respectively), without any modification. With respect to Fig. 3, the computational results of this section refer to the boundary nodes Bðx ¼ a; y ¼ 0Þ, Fðx ¼ a=2; y ¼ 0Þ and Eðx ¼ a=4; y ¼ 0Þ, where the analytical solution is given by [44]: Point Bðx ¼ a; y ¼ 0Þ:

( uða; tÞ ¼

q0 ct;

0 < t < 2ac

q0 ct  2q0 cðt  2a=cÞ;

2a c

< t < 4ac etc: ð22aÞ

Point Fðx ¼ a=2; y ¼ 0Þ: 8 0; > > > < q0 cðt  a=2cÞ; uða=2; tÞ ¼ > q a; > > 0  : q0 a  q0 c t  5a ; 2c

0 < t < 2ca a 2c 3a 2c 5a 2c

< t < 3a 2c < t < 5a 2c < t < 7a etc: 2c ð22bÞ

Point Eðx ¼ a=4; y ¼ 0Þ: 8 0; > > > < q cðt  3a=4cÞ; 0 uða=4; tÞ ¼ > q 0 a=2; > >  : ; q0 a=2  q0 c t  11a 4c

0 < t < 3a 4c 3a < t < 5a 4c 4c 5a 11a < t < 4c 4c 11a < t < 13a 4c 4c

etc: ð22cÞ

In the beginning, in order to establish the critical time step that ensures numerical stability, the maximum cir-

Table 2 Calculated eigenvalues using lumped mass formulation for (a) one, (b) two and (c) four Coons-patch macroelements (piecewise-linear interpolation: Model 1) Exact eigenvalues (s2 )

One sub-region (a)

Two sub-regions (b)

Four sub-regions (c)

0.6168 5.5516 10.4865

0.6144 3.6202 4.6660

0.6157 4.8051 5.7835

0.6158 5.2929 9.5131

Table 3 Maximum calculated eigenvalues using consistent mass matrices

Maximum calculated circular frequency (xmax ) (s1 ) Minimum calculated period Tmin (s) Critical Dtcr ¼ 2=xmax (s) (Eq. (21))

B-splines (Model 2) (n ¼ 19 DOF)

Piecewise-linear (Model 1) (n ¼ 19 DOF)

FEM (n ¼ 40 DOF)

12.94

14.09

19.46

0.49 0.155

0.45 0.142

0.32 0.103

C.G. Provatidis / Computers and Structures 82 (2004) 383–395

ð23Þ

where Tn ¼ 2p=xn is the smallest period of the finite element assemblage with n degrees of freedom. An estimate for Dtcr , unknown Ôa priori ’, is the time it takes the wave to travel across the smallest element of the mesh, that is: Dlmin D^tcr ¼ c

ð24Þ

For the particular discretization of Fig. 3 (Dlmin ¼ 0:25 m, c ¼ 1 m/s), Eq. (24) gives D^tcr ¼ 0:25 s. In the following, details are given for the three aforementioned time-integration schemes applied to the three alternative element formulations (models) described in Section 3.3.

axial displ, u/(qoa)

Tn 2 Dtcr ¼ ¼ p xn

2.5

F-Model2

E-Model2

B-Analytical

F-Analytical

E-Analytical

4.0 ct/a

6.0

B 1.5 F

1.0

E

0.5

2.0

(a)

2.5

B-Model2

F-Model2

E-Model2

B-Analytical

F-Analytical

E-Analytical

8.0

2.0 B 1.5 F 1.0 E 0.5 0.0 0.0

(b)

4.2.1. Central difference scheme It was further found that: (i) Even for 50% of the D^tcr , i.e. Dt ¼ 0:125 s, the Central Difference Method applied to the conventional FEM (Fig. 3b) diverges, a finding consistent to the (small) time step related to the (large) maximum calculated eigenvalue (DtFEM ¼ 0:103 s, in Table 3). (ii) On the contrary, for the same time step (Dt ¼ 0:125 s), Model 2 is stable probably because it is smaller than the critical time step related to the maximum calculated eigenvalue (0.155 s, in Table 3). Results of time response are presented in Fig. 5a and they are of good quality. (iii) Moreover, as it is shown in Fig. 5b, the above Model 2 converges for even a larger time step, Dt ¼ 0:15 s, an event again consistent to Table 3 (0:15 < 0:155 s). (iv) According to Table 3 and Eq. (23), the conventional FEM becomes stable when Dt < 0:103 s. Now, it is chosen that Dt ¼ 0:10 s and consequently the FEM becomes stable. Moreover, for this particular case, the results obtained through Model 1 (shown in Fig. 6a) were found to coincide with those obtained through the conventional FEM solution, although they concern two different element formulations. For the same time step, Model 2 is of superior quality as shown in Fig. 6b. (v) It is remarkable that, despite that Model 1 and conventional FEM are two different element formulations, for a broad tested spectrum of smaller time-steps (Dt < 0:10 s) in the test case of this paper Model 1 led to identical results to those obtained through conventional FEM solution (four-node finite elements). (vi) For the above mentioned small time steps, Model 2 leads to superior results than Model 1 (here numeri-

B-Model2

2.0

0.0 0.0

axial displ, u/(qoa)

cular frequencies xmax were calculated and are shown in Table 3. According to the usual recipe (see for example Ref. [2, p. 587]) the critical time step in the use of the successful implementation of the central difference method is given by:

391

2.0

4.0

6.0

8.0

ct/a

Fig. 5. Axial displacements at points Bða; 0Þ, Eða=4; 0Þ and Fða=2; 0Þ using Coons-patch macroelement: Model 2 (cubic Bsplines: 24 DOF) in conjunction with the Central Difference Method (l ¼ a=8, segment length): (a) Time step: Dt ¼ 0:50 l=c and (b) time step: Dt ¼ 0:60 l=c.

cally identical to FEM) as it is clearly shown in Fig. 7, selectively for Dt ¼ 0:05 s. (vii) As a remark, in order to ensure stability, FEM requires a time step Dt 6 0:10 s, which is smaller than 40% of the Ôa priori’ known value D^tcr ¼ 0:25 s, while Model 2 requires a larger value (Dt 6 0:15 s), which is 60% of D^tcr . 4.2.2. H-Wilson scheme It is noted that as previously occurred, even when using the h-Wilson scheme the Model 1 leads to identical results with those obtained through the conventional FE solution while Model 2 is of similar quality as shown in Fig. 8 for Dt ¼ 0:10 s. It can be also noticed that, in comparison to the corresponding numerical solution using the central-difference method (shown in Fig. 6), the h-Wilson leads to a rather smoother solution. 4.2.3. Newmark scheme As previously, Newmark scheme is applicable in conjunction with the proposed Coons-patch macroelements. It is again noted that as previously occurred, even using this scheme the Model 1 (Fig. 9a) leads to identical results with those obtained through the conventional FE solution while Model 2 (shown in Fig. 9b) is of similar

392

C.G. Provatidis / Computers and Structures 82 (2004) 383–395 B-Model1

F-Model1

E-Model1

B-Analytical

F-Analytical

E-Analytical

B

2.0 1.5

F 1.0 E

0.5

B-Model1 B-Analytical

2.5 axial displ, u/(qoa)

axial displ, u/(qoa)

2.5

F-Model1 F-Analytical

E-Model1 E-Analytical

2.0 B 1.5 F 1.0 E

0.5

0.0 0.0

2.0

4.0 ct/a

2.5

B-Model2 B-Analytical

2.0

B

F-Model2 F-Analytical

6.0

8.0

0.0 0.0

E-Model2 E-Analytical

2.0

F 1.0 E 0.5

2.0

(b)

4.0 ct/a

6.0

4.0

8.0

Fig. 6. Axial displacements at points Bða; 0Þ, Eða=4; 0Þ and Fða=2; 0Þ using a Coons-patch macroelement in conjunction with the Central Difference Method (Time step: Dt ¼ 0:40 l=c): (a) Model 1 and (b) Model 2. Remark: Here, results obtained using Model 1 (24 nodes, Fig. 3a) coincide with those obtained using four-noded bilinear finite elements (45 DOF, Fig. 3b).

quality with Model 1, for Dt ¼ 0:10 s. It can be also noticed that the quality of the solution using the Newmark is lower than that obtained using the h-Wilson scheme. As a general remark, since Model 1 and conventional FEM are two different element formulations (it was explained in detail in Section 3.3), the numerical coincidence between these models in the example of this paper probably will not be valid for any other problem, especially in cases characterized by entirely two-dimensional non-uniform boundary conditions and/or curvilinear domains. Relevant numerical experience exists from static problems [19,23]. Finally, for all time-integration schemes, numerical results were obtained at internal points where the potential was found to be identical with the corresponding at the boundary nodes. For example, the displacement values at the points H and G (shown in Fig. 3a) were found to be numerically identical to those of the nodes F and E, respectively.

5. Conclusions and discussion In this paper it was shown that Coons’ interpolation, primarily applied in the design of automotive surfaces in

6.0

8.0

ct/a 2.5

1.5

0.0 0.0

2.0

(a)

axial displ, u/(qoa)

axial displ, u/(qoa)

(a)

B-Model2

F-Model2

E-Model2

B-Analytical

F-Analytical

E-Analytical

B

1.5 F 1.0 E 0.5 0.0 0.0

(b)

2.0

4.0 ct/a

6.0

8.0

Fig. 7. Axial displacements at points Bða; 0Þ, Eða=4; 0Þ and Fða=2; 0Þ using a Coons-patch macroelement in conjunction with the Central Difference Method (Time step: Dt ¼ 0:20 l=c): (a) Model 1 and (b) Model 2. Remark: Here, results obtained using Model 1 (24 nodes, Fig. 3a) coincide with those obtained using four-noded bilinear finite elements (45 DOF, Fig. 3b).

1960s, offers the background for the systematic development of large finite elements, closely related to the serendipity family, with any number of nodes along any side of a four-sided patch. Apart from the well-known Lagrange polynomials, this paper offers a unified formulation and proposes the combination of Coons’ interpolation formula with either piecewise-linear or cubic B-splines interpolation along each of the sides of the patch under consideration. Although the numerical example of this paper refers to a regular domain, the proposed Coons-patch macroelements are generally applicable to arbitrary-shaped domains, and the degrees of freedom appear only at the boundaries of the entire field. It is also possible to divide the domain into large sub-regions either for the purpose of reducing the bandwidth or in order to deal with inhomogeneous problems characterized by different material properties. Although previous studies had shown that the so-far behaviour of the so-called Coons-patch macroelements in static and eigenvalue problems was satisfactory, these elements were applied for the first time

C.G. Provatidis / Computers and Structures 82 (2004) 383–395 B-Model1

F-Model1

E-Model1

B-Analytical

F-Analytical

E-Analytical

2.5 axial displ, u/(qoa)

axial displ, u/(qoa)

2.5

B

2.0 1.5

F 1.0 E

0.5 0.0

B-Model1

F-Model1

E-Model1

B-Analytical

F-Analytical

E-Analytical

B

1.5 F 1.0 E

0.5 0.0

0.0

2.0

4.0 ct/a

(a)

B 2.0

6.0

2.0

4.0 ct/a

6.0

B-Model2

F-Model2

E-Model2

B-Analytical

F-Analytical

E-Analytical

(a)

B-Model2

F-Model2

E-Model2

B-Analytical

F-Analytical

E-Analytical

2.5

B

1.5 F 1.0 E

0.5

0.0

8.0

axial displ, u/(qoa)

2.5 axial displ, u/(qoa)

2.0

393

2.0

8.0

B

1.5 F

1.0

E

0.5 0.0

0.0 0.0

2.0

(b)

4.0 ct/a

6.0

8.0

Fig. 8. Axial displacements at points Bða; 0Þ, Eða=4; 0Þ and Fða=2; 0Þ using a Coons-patch macroelement in conjunction with the h-Wilson method (Time step: Dt ¼ 0:40 l=c): (a) Model 1 and (b) Model 2. Remark: Here, results obtained using Model 1 (24 nodes, Fig. 3a) coincide with those obtained using fournoded bilinear finite elements (45 DOF, Fig. 3b).

in two-dimensional scalar wave-propagation problems and showed once more an excellent behaviour. In this particular class of problems, this study suggests that the proposed Coons-patch macroelements must be applied in combination with the consistent mass matrix only. The applicability of lumped mass is still an open task that requires further investigation.

0.0

2.0

4.0 ct/a

(b)

6.0

8.0

Fig. 9. Axial displacements at points Bða; 0Þ, Eða=4; 0Þ and Fða=2; 0Þ using a Coons-patch macroelement in conjunction with Newmark method (Time step: Dt ¼ 0:40 l=c): (a) Model 1 and (b) Model 2. Remark: Here, results obtained using Model 1 (24 nodes, Fig. 3a) coincide with those obtained using fournoded bilinear finite elements (45 DOF, Fig. 3b).

By substituting (A.1) in Eq. (3), one derives: DA uðn; gÞ ¼ ½E0 ðgÞLAB 1 ðnÞ þ E0 ðnÞL4 ðgÞ  E0 ðnÞE0 ðgÞ u1 AB þ ½E0 ðgÞLAB 2 ðnÞ u2 þ ½E0 ðgÞL3 ðnÞ u3 BC þ ½E0 ðgÞLAB 4 ðnÞ þ E1 ðnÞL1 ðgÞ  E1 ðnÞE0 ðgÞ u4 BC þ ½E1 ðnÞLBC 2 ðgÞ u5 þ ½E1 ðnÞL3 ðgÞ u6 CD þ ½E1 ðnÞLBC 4 ðgÞ þ E1 ðgÞL1 ðnÞ  E1 ðnÞE1 ðgÞ u7 CD þ ½E1 ðgÞLCD 2 ðnÞ u8 þ ½E1 ðgÞL3 ðnÞ u9 DA þ ½E1 ðgÞLCD 4 ðnÞ þ E0 ðnÞL1 ðgÞ  E0 ðnÞE1 ðgÞ u10

Appendix A. Application of Coons formula to derive the cubic element of serendipity family Let us consider a twelve-noded uniform Coons-patch macroelement in reference co-ordinates, with q1 ¼ q2 ¼ q3 ¼ q4 ¼ 4 nodes per side as shown in Fig. 10. We consider that along each side (AB, BC, CD and DA) the potential u is interpolated using Lagrange polynomials (j ¼ 1; . . . ; 4) as follows: LSide j AB AB AB uðn; 0Þ ¼ LAB 1 ðnÞu1 þ L2 ðnÞu2 þ L3 ðnÞu3 þ L4 ðnÞu4

uð1; gÞ ¼

LBC 1 ðgÞu4

uðn; 1Þ ¼

LCD 1 ðnÞu7

þ

LBC 2 ðgÞu5

þ

LCD 2 ðnÞu8

þ

LBC 3 ðgÞu6

þ

LCD 3 ðnÞu9

þ

LBC 1 ðgÞu7

þ LCD 4 ðnÞu10

DA DA DA uð0; gÞ ¼ LDA 1 ðgÞu10 þ L2 ðgÞu11 þ L3 ðgÞu12 þ L1 ðgÞu4

ðA:1Þ

DA þ ½E0 ðnÞLDA 2 ðgÞ u11 þ ½E0 ðnÞL3 ðgÞ u12

ðA:2Þ

where 9 LAB 1 ðnÞ ¼  ðn  1=3Þðn  2=3Þðn  1Þ; 2 9 DA L4 ðgÞ ¼  ðg  1=3Þðg  2=3Þðg  1Þ 2

ðA:3Þ

and so on. Also, we recall that E0 ðnÞ ¼ 1  n and

E0 ðgÞ ¼ 1  g

ðA:4Þ

So, the global shape function related, for example, to the node ‘‘1’’ finally becomes:

394

C.G. Provatidis / Computers and Structures 82 (2004) 383–395

η

s

10

9

8

D

7 C

11

6 r

12

5

A

1

ξ

B

2

3

4

Fig. 10. Cubic finite element of Serendipity family.

 N1 ðn; gÞ ¼ ð1  nÞð1  gÞ

 9 2 ðn  n þ g2  gÞ þ 1 2 ðA:5Þ

Obviously, the relation between the usual natural coordinates (r; s) in the parent element and (n; g) is given by n¼

ð1 þ rÞ ; 2

g¼

ð1 þ sÞ 2

ðA:6Þ

By substituting (A.6) in (A.5) one derives: N1 ðr; sÞ ¼

1 ð1  rÞð1  sÞ½10 þ 9ðr2 þ s2 Þ

32

ðA:7Þ

which is the standard expression for the cubic element [2, p. 157]. In other words, it was shown that Coons interpolation in conjunction with Lagrange polynomials leads to conventional finite elements of the serendipity family. In a similar way, the proof may be extended to the mid-side nodes. References [1] Argyris JH, Kelsey S. Energy theorems and structural analysis. London: Butterworths; 1960. Originally published in a series of articles in Aircraft Engineering, October, November 1954, February–May 1955. [2] Zienkiewicz OC. The finite element method. 3rd edition. London: McGraw-Hill; 1977. € [3] Ritz W. Uber eine neue Methode zur L€ osung gewisser Variationsprobleme der mathematischen Physik. Zeitschrift f€ ur Angewandte Mathematik und Mechanik 1908;135(1):1–61. [4] Trefftz E. Ein Gegenst€ uck zum Ritz’schen Verfahren. In: Proceedings, 2nd International Congress in Applied Mechanics, Zurich, 1926.

[5] Wachpress E. A rational finite element basis. New York: Acedemic Press; 1975. [6] Jirousek J, Theodorescu P. Large finite elements method for the solution of problems in the theory of elasticity. Comput Struct 1982;15(5):575–87. [7] Zielinski AP, Zienkiewicz OC. Generalized finite element analysis with T-complete boundary solution functions. Int J Numer Meth Eng 1985;21:509–28. [8] Zienkiewicz OC, Irons BM, Campbell J, Scott FC. Three dimensional stress analysis. Int Un Th Appl Mech Symp High Speed Comput. Elasticity, Liege, 1970. [9] Taylor RL. On completeness of shape functions for finite element analysis. Int J Numer Meth Eng 1972;4:17–22. [10] Taig IC. Structural analysis by the matrix displacement method. Engl. Electric Aviation Report No. S017, 1961. [11] Irons BM. Engineering application of numerical integration in stiffness method. AIAA J 1966;14(11):2035–7. [12] Coons SA. Surfaces for computer aided design of space form. Project MAC, MIT (1964), revised for MAC-TR-41. Springfield, VA, USA: Available by CFSTI, Sills Building, 5285 Port Royal Road, 1967. [13] Gordon WJ. Blending functions methods of bivariate multivariate interpolation and approximation. SIAM J Numer Anal 1971;8:158–77. [14] Gordon WJ, Hall CA. Transfinite element methods blending function interpolation over arbitrary curved element domains. Numer Math 1973;21:109–12. [15] El-Zafrany A, Cookson RA. Derivation of Lagrangian and Hermitian shape functions for quadrilateral elements. Int J Numer Meth Eng 1986;23:1939–58. [16] Yildir YB, Wexler A. MANDAP-A FEM/BEM preparation package. IEEE Trans Mag 1983;19:2562–5. [17] Zhaobei Z, Zhiqiang X. Coons surface method for formulation of finite element of plate and shells. Comput Struct 1987;27:79–88. [18] Kanarachos A, Deriziotis D. On the solution of Laplace and wave propagation problems using C-elements. Finite Element Anal Des 1989;5:97–109. [19] Provatidis Ch, Kanarachos A. Performance of a macroFEM approach using global interpolation (Coons) functions in axisymmetric potential problems. Comput Struct 2001;79(19):1769–79. [20] Provatidis CG. Coons-patch macroelements in potential Robin problems. Forschung im Ingenieurwesen 2002;67: 19–26. [21] Kanarachos AE, Provatidis CG, Deriziotis DG, Foteas N. A new approach of the fem analysis of two-dimensional elastic structures using global (Coons) interpolation functions. In: Wunderlich J, editor. CD Proceedings European Conference on Computational Mechanics, Munich, Germany, 1999. [22] Provatidis CG, Kanarachos AE. On the use of Coons’ interpolation in CAD/CAE systems. In: Mastorakis N, editor. Systems and control: theory and applications. World Scientific and Engineering Society Press; ISBN: 960-8052-11-4, 2000. p. 343–8, Available from: http://www.worldses.org. [23] Provatidis CG. Analysis of axisymmetric structures using Coons’ interpolation. Finite Elements Anal Des 2003;39: 535–58. [24] Kanarachos A, Grekas D, Provatidis C. Generalized formulation of Coons’s interpolation, in: Kaklis PD, Sapidis NS, editors. Computer aided geometric design: from

C.G. Provatidis / Computers and Structures 82 (2004) 383–395

[25] [26] [27] [28] [29]

[30]

[31]

[32]

[33]

theory to practice. Athens: National Technical University of Athens Press; 1995 (Chapter 7: pp.65–76). Faux ID, Pratt MJ. Computational geometry for design and manufacture. Chichester: Ellis Horwood; 1979. Farin G. Curves and surfaces for computer aided geometric design: a practical guide. Boston: Academic Press; 1990. Carey GF, Oden JT. Finite elements: a second course, vol. II. New Jersey: Prentice-Hall; 1983. Irons B, Ahmad S. Techniques of finite elements. Chichester: Ellis Horwood; 1981 (p. 76). Tai G, Shaw RP. Helmholtz equation eigenvalues and eigenmodes for arbitrary domains. J Acoust Soc Am 1974; 56:796–804. Mey De. Calculation of eigenvalues of the Helmholtz equation by an integral equation. Int J Numer Meth Eng 1976;10:59–66. Nardini D, Brebbia CA. The solution of parabolic and hyperbolic problems using an alternative boundary formulation, in: Brebbia CA, Meier G, editors. Proceedings 7th International Conference on BEM, vol1, Italy, 1985. Berlin: Springer; 1985. p. 387–97. Loeffer CF, Mansur WJ. Analysis of time integration schemes for boundary element applications to transient wave propagation problems, in: Brebbia CA, Venturini WS, editors. Boundary element techniques (BETECH 1987): applications in stress analysis and heat transfer. Southampton: Computational Mechanics Publications; 1987. p. 105–22. Provatidis CG. Application of the boundary element method in the analysis of field and dynamic problems. Doctoral Dissertation, National Technical University of Athens, 1987.

395

[34] Kanarachos AE, Provatidis CG. Performance of mass matrices for the BEM dynamic analysis of wave propagation problems. Comput Meth Appl Mech Eng 1987;63: 155–65. [35] Provatidis CG, Kanarachos AE. Further research on the performance of consistent mass matrices using BEM for symmetric/nonsymmetric formulations. Comput Mech 1995;16:197–207. [36] Ali A, Rajakumar C, Yunus SM. Advances in acoustic eigenvalue analysis using boundary element method. Comput Struct 1995;56:837–47. [37] Coyette JP, Fyfe KR. An improved formulation for acoustic eigenmode extraction from boundary element models. J Vib Acoust Trans ASME 1990;112:392–8. [38] Rashed Y. BEM for dynamic analysis using compact supported radial basis functions. Comput Struct 2002;80: 1351–67. [39] Brebbia CA, editor. Topics in boundary element research. vol. 2. Berlin: Springer-Verlag; 1985. [40] Israil ASM, Banerjee PK. Advanced development of timedomain BEM for two-dimensional scalar wave propagation. Int J Numer Meth Eng 1990;29:1003–20. [41] Clough RW, Penzien J. Dynamics of structures. Tokyo: McGraw-Hill-Kogakusha; 1975. [42] Blevins RD. Formulas for natural frequency and mode shape. New York: Van Nostrand; 1979. [43] Bathe KJ. Finite element procedures in engineering analysis. New Jersey: Prentice-Hall; 1982. [44] Pipes LA, Harvill LR. Applied mathematics for engineers and physicists. 3rd edition. McGraw-Hill International Book Company, International Student Edition; 1981. p. 494–6.