A stable 3D energetic Galerkin BEM approach for wave propagation interior problems

Engineering Analysis with Boundary Elements 36 (2012) 1756–1765

Contents lists available at SciVerse ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

A stable 3D energetic Galerkin BEM approach for wave propagation interior problems A. Aimi a, M. Diligenti a,n, A. Frangi b, C. Guardasoni a a b

´ di Parma, V.le delle Scienze 53/A, 43124 Parma, Italy Department of Mathematics, Universita Department of Structural Engineering, Politecnico di Milano, Piazza Leonardo da Vinci, 32, 23133 Milano, Italy

a r t i c l e i n f o

a b s t r a c t

Article history: Received 19 March 2012 Accepted 17 June 2012 Available online 21 July 2012

We consider 3D interior wave propagation problems with vanishing initial and mixed boundary conditions, reformulated as a system of two boundary integral equations with retarded potentials. These latter are then set in a weak form, based on a natural energy identity satisfied by the solution of the differential problem, and discretized by the energetic Galerkin boundary element method. Numerical results are presented and discussed in order to show the stability and accuracy of the proposed technique. & 2012 Elsevier Ltd. All rights reserved.

Keywords: Wave propagation Boundary integral equation Energetic Galerkin boundary element method

1. Introduction Time-dependent problems that are frequently modelled by hyperbolic partial differential equations can be dealt with the boundary integral equation (BIE) method in which the transformation of the problem to a BIE follows the same well-known method as for elliptic boundary value problems. Although the basic integral equation formulations for time-dependent problems were known for many years their systematic use for obtaining numerical solutions is a rather recent event. The boundary value problem is formulated directly in terms of boundary values and the solution at interior points does not need to be considered, although it may be evaluated at these points directly from the integral representation, once the BIE has been solved. Boundary element methods (BEMs) have been successfully applied at the discretization level both in the frequencydomain and in the time-domain (see, e.g., the reviews by Beskos [8] and Costabel [9]). Focusing on the latter case, the representation formula in terms of single layer and double layer potentials employs the fundamental solution of the hyperbolic partial differential equation and jump relations, giving rise to the so-called ‘‘retarded’’ BIEs. Usual numerical discretization procedures include the collocation technique with explicit evaluation of the time-convolution, which is very appealing thanks to its simplicity and limited demand in terms of computing effort. However, this approach has been always plagued with limited robustness. Indeed, the time

n

Corresponding author. E-mail address: [email protected] (M. Diligenti).

0955-7997/$ - see front matter & 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.enganabound.2012.06.003

marching scheme of classical space–time collocation methods often turns unstable in an erratic way upon variation of the time step adopted in the analysis (see, e.g., [15,16]). Even if this issue can be somehow cured in BEM only approaches by means of a trial and error choice of the time step, improving the basic formulation seems imperative, e.g., when a coupling with other numerical methods like the FEM is envisaged. It has been shown (see, e.g., [26]), indeed, that in staggered iterative coupling procedure between BEM and FEM the stability and accuracy of the solution in the two subdomains impose requirements which may be contradictory. It is worth mentioning that an alternative to the standard approach for the direct explicit evaluation of the time-convolution is represented by the convolution quadrature method which has been developed in [19] and then extended to different applications in [13,22–24]. It provides a straightforward way to obtain an efficient time stepping scheme in terms of the Laplace transform of the kernel function, which is crucial for many timedependent problems where only the transformed fundamental solutions are available. Although stability and convergence are assured only under strong regularity assumptions on problem data, it is empirically found that the stability of the algorithm is improved compared to the standard approach. Also the time-stepping method, based on the finite difference approximation of time derivatives, is not yet well developed. This is partly because accuracy and stability of the solution are affected by the time step size. Properties of the numerical solution computed by this method have been already studied in [5,11]. It is however only with variational techniques that the most promising results concerning stability have been obtained, also from the theoretical point of view. The application of Galerkin

A. Aimi et al. / Engineering Analysis with Boundary Elements 36 (2012) 1756–1765

BEM in both space and time has been proposed by several authors but in this direction only the weak formulation due to Ha Duong [6,7,18] yields genuine convergence results. The only drawback of the method is that it has stability constants growing exponentially in time, as stated in [9]. Recently, in [4], we have considered 2D interior problems for a temporally homogeneous scalar wave equation in the time interval ½0,T, equipped with mixed boundary and vanishing initial conditions, reformulated as a system of two BIEs with retarded potentials. Special attention has been devoted to a natural energy identity related to the differential problem, that leads to a space– time weak formulation for the BIE, having, under suitable constraints, precise continuity and coerciveness properties. Consequently it can be discretized by unconditionally stable schemes with well-behaved stability constants even for large times, as also confirmed by the numerous numerical benchmarks developed [2,3]. In the present paper we focus on the 3D extension of the above wave propagation analysis, with the aim of testing the stability properties of the so-called energetic Galerkin BEM in this fully general context. It is indeed well known even in the 2D case that interior problems, as opposed to exterior ones, represent a severe test since continuous reflection of waves in a bounded domain can give rise, at the discretization level, to significant noise and instabilities. First we reformulate the differential problem as a system of two BIEs with retarded potentials and provide explicit expressions of the involved boundary integral operators. The integral problem is then set in the energetic space–time weak form. While in previous 2D works, the hypersingular part of the energetic bilinear form was discretized directly with its hypersingular kernel (the interest reader is referred to [10] for a review article on hypersingularity), here we consider an equivalent bilinear form having only weakly singular kernels, obtained by means of the regularization approach introduced in [14] for elastodynamics problems. Galerkin BEM is used in the discretization phase, performing analytically the required double integration in time variables; then, the resulting weakly singular double integrals in space variables, whose accurate evaluation, needed to well capture the wave propagation front, is much more challenging than in 2D case, are evaluated by inner analytical and outer numerical integrations, with high precision. The excellent stability and accuracy properties of the energetic approach are shown by means of several numerical results.

2. Model problem and its energetic weak formulation We will analyze a mixed boundary value problem for the wave equation with homogeneous initial conditions, in a bounded domain O  R3 , with a (piecewise) smooth boundary G referred to a Cartesian reference system ðx1 ,x2 ,x3 Þ> ¼ x and partitioned into two non-intersecting subsets Gu and Gp such that Gu [ Gp ¼ G: utt nu ¼ 0, uðx,0Þ ¼ 0,

x A O, t A ð0,TÞ, ut ðx,0Þ ¼ 0,

uðx,tÞ ¼ uðx,tÞ, pðx,tÞ :¼

x A O,

ðx,tÞ A SuT :¼ Gu  ½0,T,

@u ðx,tÞ ¼ pðx,tÞ, @nx

ðx,tÞ A SpT :¼ Gp  ½0,T,

ð2:1Þ ð2:2Þ ð2:3Þ

ð2:4Þ

where nx is the unit outward normal vector in xA G and u,p are given boundary data of Dirichlet and Neumann type, respectively.

1757

Let us consider the boundary integral representation of the solution of (2.1)–(2.4), for xA O,t A ð0,TÞ Z Z t Z Z t @G Gðr,ttÞpðn, tÞ dt dgn  ðr,ttÞuðn, tÞ dt dgn , uðx,tÞ ¼ G 0 G 0 @nn ð2:5Þ where r ¼ JrJ2 ¼ JxnJ2 and Gðr,ttÞ ¼

1 d½ttr 4pr

ð2:6Þ

is the forward fundamental solution of the 3D wave operator, with d½ the Dirac distribution. With a limiting process for x tending to G we obtain the space–time BIE (see [9]) Z Z t Z Z t 1 @G uðx,tÞ ¼ Gðr,ttÞpðn, tÞ dt dgn  ðr,ttÞuðn, tÞ dt dgn , 2 G 0 G 0 @nn ð2:7Þ which can be written, with obvious meaning of notation, in the compact form: 1 uðx,tÞ ¼ ðVpÞðx,tÞðKuÞðx,tÞ: 2

ð2:8Þ

Note that the Cauchy singular operator K can be rewritten as   Z t Z @r uðn, tÞ dt dgn : Gðr,ttÞ ut ðn, tÞ þ ð2:9Þ Kuðx,tÞ ¼  r G @nn 0 Expression (2.9) can be obtained starting from the definition of the double layer operator K and observing that @ @ d½ttr ¼ d½ttr: @r @t

ð2:10Þ

Indeed, using (2.10) we have:   @G @G @r 1 @ 1 @r ðr,ttÞ ðr,ttÞ ¼ ¼ d½ttr @nn @r @nn 4p @r r @nn   1 1 1 @ @r ¼  2 d½ttr þ d½ttr : 4p r @t @nn r

ð2:11Þ

Now, inserting (2.11) in the definition of K, integrating in the sense of distributions the term containing the derivative with respect to t, one gets, up to the factor 1=ð4pÞ,  Z t Z @r d½ttr d½ttr @ uðn, tÞ þ uðn, tÞ dt dgn , ð2:12Þ 2 r @t r G @nn 0 from which (2.9) can be immediately deduced. The BIE (2.8) is generally used to solve Dirichlet problems. In the case of mixed problems, however, one can consider a second space– time BIE, obtained from (2.5) by taking the normal derivative with respect to nx and operating a limiting process for x tending to G: Z Z t 1 @G pðx,tÞ ¼ ðr,ttÞpðn, tÞ dt dgn 2 G 0 @nx Z Z t @2 G  ðr,ttÞuðn, tÞ dt dgn , ð2:13Þ @n x @nn G 0 which can be written in the compact form 1 pðx,tÞ ¼ ðK 0 pÞðx,tÞðDuÞðx,tÞ: 2 Note that also the operator K 0 can be rewritten as   Z Z t @r pðn, tÞ dt dgn , Gðr,ttÞ pt ðn, tÞ þ K 0 pðx,tÞ ¼  r G @nx 0

ð2:14Þ

ð2:15Þ

further, considering the derivative with respect to nx of (2.9) and operating with the same arguments as before, after a cumbersome but easy calculation the hypersingular integral operator D in (2.14) can be equivalently expressed in the following way:   Z Z t @2 r uðn, tÞ Duðx,tÞ ¼  dt dgn Gðr,ttÞ ut ðn, tÞ þ r G @nx @nn 0

1758

A. Aimi et al. / Engineering Analysis with Boundary Elements 36 (2012) 1756–1765

Z þ

@r @r G @nx @nn

Z

t 0

  2ut ðn, tÞ 2uðn, tÞ þ Gðr,ttÞ utt ðn, tÞ þ dt dgn : r r2 ð2:16Þ

Taking into account the presence of the Dirac distribution in the kernels of the above boundary integral operators V,K,K 0 ,D, the time integration can be performed analytically yielding the following final simplified expressions: Z 1 Vpðx,tÞ ¼ H½trpðn,trÞ dgn , ð2:17Þ G 4pr Z Kuðx,tÞ ¼ 

  1 @r uðn,trÞ dgn , H½tr ut ðn,trÞ þ @nn r G 4pr

ð2:18Þ

  1 @r pðn,trÞ dgn , H½tr pt ðn,trÞ þ @nx r G 4pr

ð2:19Þ

Z K pðx,tÞ ¼  0

 2   1 @ r uðn,trÞ H½tr ut ðn,trÞ þ @nx @nn r G 4pr   @r @r 2ut ðn,trÞ 2uðn,trÞ þ þ utt ðn,trÞ þ dgn , 2 @nx @nn r r

Duðx,tÞ ¼ 

Z

ð2:20Þ where H½ is the Heaviside function. At this stage, using the boundary conditions (2.3) and (2.4), a mixed boundary value wave propagation problem can be rewritten as a system of two BIEs in the boundary unknowns the functions pðx,tÞ and uðx,tÞ on Gu , Gp , respectively: " #" # " # 1 V p Vu K p pðx,tÞ 2 I þ Ku ¼ K 0u Dp uðx,tÞ  12 I þ K 0p Du " # ðx,tÞ A SuT pðx,tÞ ,  , ð2:21Þ ðx,tÞ A SpT uðx,tÞ where the boundary integral operator subscripts b ¼ u,p define their restriction to SbT . Note that in order to assure the correctness of the integral operator splitting, the Dirichlet datum is supposed to be locally extended to zero by continuity outside Gu over Gp . Now, we remark that the solution of (2.1)–(2.4) satisfies the following energy identity: Eðu,TÞ : ¼

1 2

Z O

2

½u2t ðx,TÞ þ9ruðx,TÞ9 dx dt ¼

Z Z G

T 0

ut ðx,tÞpðx,tÞ dt dgx ,

which can be obtained multiplying Eq. (2.1) by ut and integrating by parts over O  ½0,T. Then, also taking into account Eqs. (2.8)–(2.14), the energetic weak formulation of the system (2.21) is finally defined as follows: 1=2 Find u A H1 ð½0,T; H0 ðGp ÞÞ and p A L2 ð½0,T; H1=2 ðGu ÞÞ such that 8 u < /ðV u pÞt , cSL2 ðSu Þ /ðK p uÞt , cSL2 ðSu Þ ¼ /f t , cSL2 ðSu Þ T T T , ð2:22Þ p : /K 0u p, Zt SL2 ðSp Þ þ /Dp u, Zt SL2 ðSp Þ ¼ /f , Zt SL2 ðSp Þ T

T

T

Note that the involved scalar products are represented by a space–time integral; hence, remembering (2.17)–(2.20) we will have to deal with still one integration in time and a double integration in space variables.

3. Galerkin BEM discretization For time discretization we consider a uniform decomposition of the time interval ½0,T with time step Dt ¼ T=NDt , N Dt A N þ , generated by the N Dt þ 1 time-knots t k ¼ kDt,

and we choose temporally piecewise constant shape functions for the approximation of p and piecewise linear shape functions for the approximation of u, although higher degree shape functions can be used. Note that, for this particular choice, temporal shape functions, for k ¼ 0, . . . ,NDt 1, will be defined as vpk ðtÞ ¼ H½tt k H½tt k þ 1 

ð3:1Þ

for the approximation of p, or as vuk ðtÞ ¼ Rðtt k ÞRðtt k þ 1 Þ,

ð3:2Þ

for the approximation of u, where Rðtt k Þ ¼ ððtt k Þ=DtÞH½tt k  is the ramp function. Note that vpk ðtÞ and vuk ðtÞ will give contribution, respectively, to the first time derivative in (2.19) and to the second time derivative in (2.20) in terms of Dirac distributions. For the space discretization, we employ a Galerkin boundary element method. We consider O (suitably approximated by a domain) of polyhedral type and a boundary mesh on Gu constituted by Mu flat triangles feu1 , . . . ,euMu g, with diamðeui Þ r Dxu , eui \ S u u euj ¼ | if i a j and such that M i ¼ 1 e i ¼ G u . The same is done for the Neumann part of the boundary Gp with obvious change in notation. Let Dx ¼ maxfDxu , Dxp g. The functional background compels one to choose spatially shape functions belonging to L2 ðGu Þ for the approximation of p and to H10 ðGp Þ for the approximation of u. Having defined P di the space of algebraic polynomials of degree di on the element ei of G, we consider, respectively, the space of piecewise polynomial functions X 1, Dx :¼ fwp A L2 ðGu Þ : wp9eu A P di 8eui  Gu g,

ð3:3Þ

i

and the space of continuous piecewise polynomial functions X 0, Dx :¼ fwu A C 0 ðGp Þ : wu9ep A P dj 8epj  Gp g:

ð3:4Þ

j

Hence, denoted with M pDx , M uDx the number of unknowns on Gu and Gp , respectively, and having introduced the standard piecewise polynomial boundary element basis functions wpj ðxÞ,j ¼ 1, . . . ,M pDx , in X 1, Dx and wuj ðxÞ,j ¼ 1, . . . ,M uDx in X 0, Dx , the approximate solutions of the problem at hand will be expressed as ~ pðx,tÞ :¼

MpDx NX Dt 1 X

aðkÞ wpj ðxÞvpk ðtÞ, pj

k¼0 j¼1

where   u f t ¼ ðV p pÞt þð 12I þ K u Þu t ,

k ¼ 0, . . . ,NDt ,

p

f ¼ ð12I þK 0p ÞpDu u

and cðx,tÞ, Zðx,tÞ are suitable test functions, belonging to the same functional space of pðx,tÞ, uðx,tÞ, respectively (see also [4], where the energetic weak formulation was introduced for the problem (2.1)–(2.4), but with O  R2 ). In particular, the first equation in (2.21) has been differentiated with respect to time and projected with the L2 ðSuT Þ scalar product by means of functions belonging to L2 ð½0,T; H1=2 ðGu ÞÞ, while the second equation in (2.21) has been projected with the L2 ðSpT Þ scalar product by means of functions 1=2 belonging to H1 ð½0,T; H0 ðGp ÞÞ, derived with respect to time.

~ uðx,tÞ :¼

MuDx NX Dt 1 X

aðkÞ wuj ðxÞvuk ðtÞ: uj

k¼0 j¼1

The Galerkin BEM discretization coming from energetic weak formulation (2.22) produces the linear system

Ea ¼ b,

ð3:5Þ

where matrix E has a block lower triangular Toeplitz structure, since its elements depend on the difference t h t k and in particular they vanish if t h rt k . Each block has dimension MDx :¼ M pDx þ M uDx . If we denote by Eð‘Þ the block obtained when t h t k ¼ ð‘ þ1ÞDt,

A. Aimi et al. / Engineering Analysis with Boundary Elements 36 (2012) 1756–1765

‘ ¼ 0, . . . ,NDt 1, the linear system can be written as 0 B B B B B B @

ð0Þ

E Eð1Þ Eð2Þ

0

0

...

Eð0Þ Eð1Þ

0

...

Eð0Þ

^

^

^

... ^

EðNDt 1Þ

EðNDt 2Þ

...

Eð1Þ

0

b

B B B ¼B B B @

1

ð0Þ

CB 0 CB CB B 0 C CB CB 0 A@

Eð0Þ

ð0Þ

a að1Þ að2Þ ^

b

ð2Þ

b ^

where

1

Dðr, Dh þ a,k þ b Þ ¼

C C C C C C A

aðNDt 1Þ

C C C C, C C A

ð1Þ

b

0

10

ð3:6Þ

ðN Dt 1Þ

where

að‘Þ ¼ ðajð‘Þ Þ, bð‘Þ ¼ ðbjð‘Þ Þ, ‘ ¼ 0, . . . ,NDt 1; j ¼ 1, . . . ,MDx : Note that each block has a 2  2 block sub-structure of the type 2 ð‘Þ 3 ð‘Þ

Eð‘Þ ¼ 4

Euu

Eup

Eðpu‘Þ

ð‘Þ Epp

5,

ð3:7Þ

where symmetric diagonal sub-blocks have dimensions, respectively, M pDx ,M uDx and, since basis functions (3.1) and (3.2) satisfy the property ‘Þ ð‘Þ > ðd=dtÞvuk ðtÞ ¼ vpk ðtÞ, Eðpu ¼ ðEup Þ as proven in Appendix A. Further, the unknowns are organized as follows: ð‘Þ ð‘Þ ð‘Þ ð‘Þ að‘Þ ¼ ðap1 , . . . , apM Þ> : p , au1 , . . . , a uM u Dx

Dx

The solution of (3.6) is obtained with a block forward substitution, i.e. at every time instant t ‘ ¼ ð‘ þ 1ÞDt, ‘ ¼ 0, . . . ,N Dt 1, one computes ð‘Þ

zð‘Þ ¼ b 

‘ X

EðjÞ að‘jÞ

j¼1

and then solves the reduced linear system:

Eð0Þ að‘Þ ¼ zð‘Þ :

1759

ð3:8Þ

Procedure (3.8) is a time-marching technique, where the only matrix to be inverted is the non-singular block Eð0Þ ; therefore the LU factorization needs to be performed only once and saved. At each time step the solution of (3.8) requires only a forward and backward substitution phases, while all the other blocks are used to update at every time step the right-hand side. Owing to this procedure we can construct and store only the blocks Eð0Þ , . . . , EðNDt 1Þ with a considerable reduction in computational cost and memory requirement. Having set Dhk ¼ t h t k , the matrix elements in blocks of the ‘Þ ð‘Þ type Eðuu , Eup , after an analytic integration in the time variable, are of the form, respectively Z Z 1 X ð1Þa þ b 1 H½Dh þ a,k þ b rwpj ðnÞ dgn dgx , wpi ðxÞ ð3:9Þ 4p Gu Gu r a, b ¼ 0

1 n

nn  nx wuj ðnÞwui ðxÞ ðDtÞ2  1 þ ðDh þ a,k þ b rÞ2 ½nn  rwuj ðnÞ  ½nx  rwui ðxÞ : 2 ð3:12Þ

We observe in (3.9)–(3.11) space singularities of type Oðr 1 Þ as r-0, which are typical of weakly singular kernels related to 3D elliptic problems, and also the presence of the Heaviside function which analytically represents the wave front propagation. The double integration in space variables has been developed with the element by element technique over the triangulation of G, following the strategy proposed in [20]. The outer integration on the source triangle has been carried out numerically, using standard Gauss–Hammer rules [12,17]. As a consequence, we can now focus on a fixed point x over the outer (source) triangle that shall be actually a node of the exterior numerical integration. The integration over the inner (field) triangle has been performed as follows. The source point x is projected onto the plane p containing the field triangle so that the distance between the source point and a field point n becomes r 2 ¼ z2 þ r2 with z equal to the distance between the source point and the plane p and r equal to the distance between the projection point Px and the field point, as represented in Fig. 1. Since z is constant for the inner integral, we fix p as a working plane and we select a 2D Cartesian reference system centered at Px . The presence of the Heaviside functions in (3.9)–(3.11) implies that the integration over the inner triangle will be in general limited to the portion enclosed in a spherical surface of radius r ¼ Dh þ a,k þ b centered at x or, equivalently, in a circle (hereafter called wave front) of radius qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi r ¼ D2h þ a,k þ b z2 centered at Px and defined on p. In order to establish a simple and general procedure to account for the presence of the wave front in an exact way, the inner integration over the field triangle ei is expressed as an algebraic sum of integrations over three triangles eðkÞ ,k ¼ 1; 2,3, having Px as one of their vertices, as illui strated for instance in the three pictures of Fig. 2. Each of the three sub-triangles eðkÞ is now addressed sepai rately. For a given eðkÞ a circular wave front induces, in the most i general case, the partition depicted in Fig. 3 for the wave front of radius r , represented by (the sum of) a triangle and two circular sectors. In specific cases, one or more of these three regions might vanish. The integrals over the three regions can then be computed either by numerical or by analytical techniques generally

Z Z 1 X ð1Þa þ b Dh þ a,k þ b 1 H½Dh þ a,k þ b r wpi ðxÞ 4pDt Gu Gp r a, b ¼ 0 

r  nn u wj ðnÞ dgn dgx : r2

ð3:10Þ

‘Þ For what concerns matrix elements in block Eðpp , after a regularizing procedure of the corresponding hypersingular bilinear form oDp u, Zt 4 L2 ðSp Þ (see (2.22)), done as described in [14] for elastoT dynamics, and an analytic integration in the time variable, one obtains Z Z 1 X ð1Þa þ b 1 H½Dh þ a,k þ b rDðr, Dh þ a,k þ b Þ dgn dgx , 4p Gp Gp r a, b ¼ 0

ð3:11Þ

Fig. 1. Projection of the source point x onto the plane of the inner (field) triangle of integration.

1760

A. Aimi et al. / Engineering Analysis with Boundary Elements 36 (2012) 1756–1765

Fig. 2. Decomposition of the inner triangle (upper figure) into three triangles having Px as vertex (lower figure).

Now, in order to test the procedure at element level, let us consider the double integral Z Z 1 wj ðxÞ H½Dhk r wi ðnÞdgn dgx , ð3:13Þ r e e

Fig. 3. General partition of one triangle eðkÞ due to the wave front. i

Table 1 Relative errors obtained evaluating the double integral (3.13) over the reference right-angled triangle e and varying the parameter Dhk . Number of Gauss–Hammer nodes:

1 3 7 12 19 27

over the reference triangle e ¼ fðx1 ,x2 Þ : 0 o x1 o 1; 0 o x2 o 1x1 g, with wj ðxÞ ¼ wi ðnÞ ¼ 1, for different values of Dhk ¼ 10; 1,0:1: Table 1 collects the relative errors (with respect to the analytical reference value) obtained using an increasing number of Gauss– Hammer nodes for the outer integration. It is worth remarking that the three values of Dhk have been selected in order to investigate the accuracy in three qualitatively different configurations. Indeed, for Dhk ¼ 10 the wave front encloses all the field triangle irrespective of the source point position; for Dhk ¼ 1 only small portions are excluded when the source point is close to the vertices; finally, for Dhk ¼ 0:1 wave front effects become dominant. In this latter case, with one or three Gauss–Hammer points, the wave front lies always completely within the field element never intersecting the outer boundary, which explains the rather surprising fact that the relative error remains constant.

4. Numerical results

Relative errors for

Dhk ¼ 10

Dhk ¼ 1

Dhk ¼ 0:1

1.9993E  001 3.1476E  002 5.3177E  003 2.4513E  003 8.2409E  004 2.7793E  004

2.0297E  001 3.3409E  002 5.3168E  003 2.3290E  003 9.1505E  004 2.9173E  004

1.1680E  001 1.1680E  001 5.4003E  002 8.4610E  003 4.1134E  003 8.3591E  005

resorting to polar coordinates with origin in Px . In the present investigation, analytical formulae based on results obtained in [21] have been employed. This procedure, which can account for wave fronts accurately, can be however conveniently replaced with much faster fully numerical approaches under certain conditions on the distance between source and field triangles. An in depth investigation of these and other major numerical issues will be the topic of a forthcoming paper.

To validate the energetic Galerkin BEM approach, in this section we address three numerical examples endowed with mixed boundary conditions. In the current implementation we employ constant space–shape functions for the approximation of p and linear space– shape functions for the approximation of u. As recalled in the previous section, double time integrals are performed analytically as well as the inner space integral over the field element while, on the contrary, the outer space integral over the source element is performed numerically. When not explicitly otherwise stated, a Gauss–Hammer rule with seven nodes has been employed, even if a faster rule with three nodes typically yields very accurate overall results, as verified in the first example. We first consider the classical benchmark (see for instance [15]) of a bar O of height equal to 10 and square cross section with side equal to 2. On the lower surface the Dirichlet boundary data u ¼ 0 is enforced, while the upper surface is subjected to an uniform Neumann condition p ¼ H½t. On the lateral surface p ¼ 0 is assumed.

A. Aimi et al. / Engineering Analysis with Boundary Elements 36 (2012) 1756–1765

1761

Δ t = 0.25

Δ t = 0.125

0

0

−0.5

−0.5

−1

−1

p

0.5

p

0.5

−1.5

−1.5

−2

−2

−2.5 0

100

200

−2.5 0

300

100

t

200

300

t

Fig. 7. Time history of the numerical solution (continuous line) p compared with the analytical one (dashed line), for Dt ¼ 0:25 on the left and Dt ¼ 0:125 on the right.

Fig. 4. Configuration and domain triangulation for the first benchmark.

1

0.5

u evaluated in point P(x=−1,y=−1,z)

u at height = 10 0

Δ t = 0.0625 Δ t = 0.125 Δ t = 0.25 Δ t = 0.5 analytic

Δ t = 0.0625 0

0

0 −0.5

Δ t = 0.125 analytical

−10

p

−5

−2

−10

200

300

−2

z=2 z=4 z=6 z=8 z = 10 50 t

−15

100

−20 0

t

−3

Δt=1

100

−1

−1

−1.5

−1.5

−2

−2

300

280 t

300

0 −0.5

p

0 −0.5 p

0

200

0.5

0.5

−0.5

100

Fig. 8. On the left, time history of the numerical solution (continuous line) p compared with the analytical one (dashed line), for Dt ¼ 0:0625 and, on the right, highlight of the behavior of p at the end of interval of observation refining only the time step.

Δ t = 0.5

0.5

0

−2.5 260

t

Fig. 5. For Dt ¼ 0:125, on the left, time history of the numerical solution u compared with the analytical one and, on the right, evolution of u at different heights along a vertical edge of the bar.

p

−1 −1.5

−15 −20 0

−1

p

−5

−1 −1.5

−2.5 0

100

200 t

300

−2.5 0

100

200

300

−2

t −2.5

Fig. 6. Time history of the numerical solution (continuous line) p compared with the analytical one (dashed line), for Dt ¼ 1 on the left and Dt ¼ 0:5 on the right.

We start considering, on G, a uniform mesh of 396 rightangled p triangles with cathetus equal to 2=3 as shown in Fig. 4 ffiffiffi (Dx ¼ 2 2=3) and the time interval of analysis ½0; 300 discretized with different time steps. The analytical behavior of the time history of u is well captured for almost every reasonable choice of Dt and is shown here in Fig. 5 only for Dt ¼ 0:125. The picture on the left presents the whole time history computed at one of the points of the upper surface while on the right, the behavior of u at some points placed along a vertical edge of the bar is highlighted only up to time 100 to have a clearer representation. As well expected, the accuracy of the history of traction p on the lower surface is more crucial. In Figs. 6–8, p at the center of the lower surface is plotted for the different values Dt ¼ 1,0:5,0:25,0:125, 0:0625. In particular we remark that the respective ratios

0

100

200

300 t

400

500

600

Fig. 9. Time history of the numerical solution p for Dt ¼ 0:125 extended to final time T ¼600.

b ¼ Dt=Dx become much smaller than the unitary wave propagation velocity of the continuum differential problem (2.1)–(2.4) and would be expected to eventually generate unstable results with any standard approach (see, e.g., [15]). On the contrary, adopting the proposed energetic approach, all the analyses remain stable despite the initial presence of spurious oscillations across the discontinuities. Moreover, even smaller, unrealistic, values of Dt were tested without succeeding in turning the analysis unstable. The oscillations in the plot of p are clearly associated with the jump discontinuities of the analytical solution, but anyway they always tend to decrease with time, as can be remarked in Fig. 9 (Dt ¼ 0:125) in the case of a longer analysis where the computation has been carried on to the final time T¼600. Fig. 8 on the

1762

A. Aimi et al. / Engineering Analysis with Boundary Elements 36 (2012) 1756–1765

0.5

Table 2 Relative errors w.r.t. the analytical solution in norm L2 ð½0,TÞ for u at one point of the mesh on the upper surface and for p in one triangle of the mesh on the lower pffiffiffi surface (Dx ¼ 2 2=3, Dt ¼ 1=2).

0

p

−0.5 −1

Number of Gauss–Hammer nodes

−1.5

Relative errors for u at height ¼ 10

p at height ¼ 0

0.0174 0.0170 – –

0.1579 0.1572 0.1570 –

−2 −2.5 0

10

20

30

40

50 t

60

70

80

90

100

110

120

130

140

150 t

160

170

180

190

200

0.5

3 7 12 19

0

p

−0.5 −1 −1.5 −2 −2.5 100 0.5 0

p

−0.5

Fig. 11. Two different triangulation of the considered portion of the sphere (with 421 triangles and 5100 triangles respectively) for the second numerical example.

−1 −1.5

421 triangles

220 Δ4

230

240 Δ3

250 t Δ2

260 Δ1

270

280

290

300

analytical

Fig. 10. Behavior of the numerical solution p along the whole time interval ½0; 300 refining discretization p both in time and space, but keeping constant the ratio ffiffiffi b ¼ Di ¼ Di t=Di x ¼ 3=ð4 2Þ, i ¼ 1, . . . ,4.

0

0

−2

−2

−4

−4

u

210

2

u

−2.5 200

5100 triangles

2

−2

−6

−6

−8

−8 −10 −12 0

10

20

t

right, where the evolution of p at the end of the time interval of investigation is analysed as a function of the time step Dt, should be compared with Fig. 10 on the bottom, where the same evolution is observed when refining both the time step and the space discretization noting improvements in accuracy. In this pffiffiffi latter case we keep the ratio b ¼ 3=ð4 2Þ constant, that is, instead of refining only the discretization pffiffiffi in time, we uniformly refine both the triangulation (D1 x ¼ 2, D2 x ¼ 2D1 x=3, D3 x ¼ 2D1 x=5, D4 x ¼ 2D1 x=7) and the time discretization (D1 t ¼ 3=4, D2 t ¼ 1=2, D3 t ¼ 3=10, D4 t ¼ 3=14). A convergence towards the analytic solution is observed and the peaks of the spurious oscillations of p remain remarkably smaller than when refining only the time step. We remark also that in Fig. 10 these oscillations seem to appear after the jumps, while in Figs. 6–8 we note that they can appear after or before the jumps depending on the parameter b (as, in [25], it can be observed comparing the behavior of the Convolution Quadrature Method with respect to the different Runge–Kutta methods adopted): having fixed Dx, when Dt ¼ 0:5 peaks are located after the jumps, while on the contrary when Dt ¼ 0:0625 peaks are located in front of them. We further stress that, with the energetic formulation, the use of time steps that are not submultiple of the wave period (as Dt ¼ 3=4; 3=10; 3=14 in Fig. 10) does not result in instabilities as it is often the case with other formulations also in 1D (see [1]). Finally, the sensitivity with respect to the number of Gauss– Hammer quadrature nodes employed for the outer space

Δ t = 0.1 Δ t = 0.05 analytic 30

−10 −12 0

10

20

Δ t = 0.05 Δ t = 0.025 analytic 30

t

Fig. 12. u at a point of the outer surface of the sphere with two different meshes and time steps.

integration has been investigated and results are collected in Table 2 which presents relative errors w.r.t. the analytical solution in L2 ð½0,TÞ norm for u and p at selected points of the mesh. In Table 2 the symbol – means that the corresponding relative error is the same of the former line. It appears that the choice of the integration rule has a minimal impact on global error measures. Indeed one major concern about variational approaches is computing time, which can be considerably mitigated using low order Gauss–Hammer rules without loosing accuracy and stability. As a second benchmark to validate the proposed approach, we consider a solid comprised within two spherical surfaces of inner radius Ri ¼1 and outer radius Ro ¼3. The inner surface is constrained and the Dirichlet boundary data u ¼ 0 is assigned. The outer surface is on the contrary subjected to a uniform load p ¼ H½t. Since the problem has spherical symmetry we analyse a portion delimited by the four planes through the origin with normal vectors ð0; 1,0Þ 7 ð0; 0,1Þ and ð0; 1,0Þ 7 ð1; 0,0Þ, along which the symmetry condition p ¼ 0 is enforced, as shown in Fig. 11. The analytical solution has been obtained as briefly described in Appendix B. Reproducing the complicated traction plot represents a formidable challenge for the proposed numerical methodology. Two different meshes, with 421 and 5100 elements respectively (corresponding to Dx C0:7 and Dx C0:2), have been

A. Aimi et al. / Engineering Analysis with Boundary Elements 36 (2012) 1756–1765

421 triangles

5100 triangles

p

−20 0

10

20

t

u(S3,⋅)

p(e ,⋅) 3

1

0

−10

Δ t = 0.1 −20 Δ t = 0.05 0 30 analytic

2

2

0

−50

−15

−15

p(e ,⋅)

2

u(S ,⋅)

50

−5

−5

1

u(S1,⋅)

100

0

0

p

p(e ,⋅)

5

5

−10

1763

−1

10

20

Δ t = 0.05 Δ t = 0.025 30 analytic

t

Fig. 13. p at a point of the inner surface of the sphere with two different meshes and time steps.

−100 −2

−150

−3

−200 0

500 t

1000

0

500 t

1000

Fig. 16. Left: evolution of u at some points of the plate (see Fig. 14). Right: evolution of p in three triangles of the bottom of the plate.

Fig. 14. Geometry and loading conditions for the plate of the third numerical example and adopted mesh.

In Fig. 15 the contour plots of u at two different time values are presented, showing that the wave front, which is initially straight, is deformed by the hole and later recovers its original shape. Fig. 16, on the left, collects the time-history of u at points of coordinates S1 ð17:3,0,96:8Þ, S2 ð17:3,0; 62Þ, S3 ð17:3,0,27:2Þ (see Fig. 14). As expected, the behavior recalls that of the bar in tension of the first example, with perturbations due to the presence of the hole. Further, the plots in Fig. 16 contain, on the right, the traction time-history in the three triangular elements e1 ,e2 ,e3 of the mesh whose centers have coordinates ð19:5,2:4,0Þ,ð39,2:4,0Þ,ð58:5, 2:4,0Þ. In these cases the hole introduces several local reflections, but the analysis remains stable irrespective of the time step adopted.

5. Conclusions

Fig. 15. Snapshot of u over the plate in two different time instants.

tested also varying the time step. The computed solution u, at a point of the outer surface, and p, at a point of the inner surface, rapidly converge towards the analytical prediction as can be observed in Figs. 12 and 13. No instabilities appear. As expected, when the ratio Dt=Dx deviates considerably from unity, the accuracy of the solution is degraded. This is case for the plots on the left associated to the coarser mesh and Dt ¼ 0:05, yielding Dt=Dx C 0:07. Refining the space discretization (figures on the right) an excellent agreement is recovered. For the final example we have adopted the geometry and loading conditions employed in [27] and in several papers of elastodynamics. This example mimics a more realistic problem for which numerical results are sometimes compared with experimental data. The thin plate with a hole depicted in Fig. 14 is subject to the Heaviside load (p ¼ H½t) on the upper surface and it is fixed at the bottom (u¼0); elsewhere the Neumann boundary condition p ¼ 0 is enforced. The numerical solution has been obtained with the triangulation shown in Fig. 14 on the right and with time step Dt ¼ 0:5.

In this paper, the energy based weak form of the boundary element formulation with retarded potentials is presented for 3D interior wave propagation problems with vanishing initial and mixed boundary conditions. Galerkin BEM is used in the discretization phase with regularization of the hypersingular bilinear form and with double analytical integration in time variable. The resulting weakly singular double integrals in space variables are then evaluated by inner analytical and outer numerical integrations. It is well known that interior problems, unlike exterior ones, represent a severe test since continuous reflection of waves can give rise to significant noise and instabilities. Through several numerical computations for some sample problems, stability and accuracy properties of the energetic approach are shown. Based on these results, adaptive procedures and development of hybrid schemes (BEM–FEM coupling) for dynamic analysis are two research areas which will be investigated in the next future by the energetic approach.

Appendix A. Symmetry of matrix blocks Eð‘Þ The following property can be proved. Proposition A.1. If time basis functions vuk ðtÞ and vpk ðtÞ for the approximation, respectively, of the unknowns u and p are chosen in such a way that ðd=dtÞvuk ðtÞ ¼ vpk ðtÞ, the matrix blocks 2 ð‘Þ 3 ð‘ Þ

Eð‘Þ ¼ 4

Euu

Eup

‘Þ Eðpu

‘Þ Eðpp

are symmetric.

5,

‘ ¼ 0, . . . ,N Dt 1

ðA:1Þ

1764

A. Aimi et al. / Engineering Analysis with Boundary Elements 36 (2012) 1756–1765

ð‘Þ ð‘Þ Proof. The sub-blocks Euu and Epp are obviously symmetric (see (3.9) and (3.11)), while, due to the assumption made in the ð‘Þ ‘Þ > statement of the proposition, we prove that Epu ¼ ðEðup Þ . Consider, for h ¼ 1, . . . ,NDt , k ¼ 0, . . . ,h, the generic element of ð‘Þ , defined by: the sub-block Eup Z TZ @ ‘Þ ½Eðup ij ¼  wpi ðxÞvph ðtÞ ðK p ½wuj vuk Þðx,tÞ dgx dt, @t Gu 0

i ¼ 1, . . . ,M pDx , j ¼ 1, . . . ,M uDx :

þ

fr ðRo ,tÞ Ro

¼ H½t:

ðB:3Þ

In terms of the general solution

fðr,tÞ ¼ f ðrtÞ þ gðr þtÞ

f ðRo tÞ gðRo þ tÞ þ , Ro Ro f ðRi tÞ R2i

ðB:8Þ

ðB:9Þ

0

þ

f ðRi tÞ gðRi þ tÞ g 0 ðRi þ tÞ  þ : Ri Ri R2i

Ro g 0 ðzÞgðzÞ ¼ R2o H½zRo ,

with boundary conditions: R2o

ðB:7Þ

ðB:10Þ

The adopted solution strategy is the following. Consider first 0 ot o DR. Initial conditions impose f ðzÞ ¼ 0, z 4 Ri and gðzÞ ¼ 0 for z o Ro . The Neumann datum is applied on the outer surface and generates an in-going wave. From (B.8):

ðB:1Þ

fðRo ,tÞ

gðzÞ ¼ f ð2Ri zÞ

pðRi ,tÞ ¼ 

ðB:2Þ

pðRo ,tÞ ¼ 

We now focus on the evaluation of g(z) and f(z). It is clear that, due to (B.5), the following condition will always hold:

uðRo ,tÞ ¼

which, setting u ¼ f=r, gives

¼ 0,

ðB:6Þ

After obtaining g and f, u at Ro and p at Ri will be computed as

In spherical coordinates the wave equation becomes:

Ri

0

¼ Ro g 0 ðzÞgðzÞ þ Ro g 0 ðz2DRÞ þgðz2DRÞ ¼ R2o H½zRo :

Appendix B. Analytic solution for second test problem

fðRi ,tÞ

gðRo þtÞ þ Ro g 0 ðRo þ tÞf ðRo tÞ þ Ro f ðRo tÞ ¼ R2o H½t:

0

By hypothesis, ðd=dtÞvuk ðtÞ ¼ vpk ðtÞ and remembering also (2.15), one finally gets: Z TZ Z tZ @ u @r d½ttr ‘Þ ½wj ðnÞvuh ðtÞ ij ¼ ½Eðup 4pr Gp @t Gu @nn 0 0 ! wpi ðxÞvpk ðtÞ @  þ ½wpi ðxÞvpk ðtÞ dgx dt dgn dt @t r Z TZ @ u ð‘Þ ½wj ðnÞvuh ðtÞðK 0u ½wpi vpk Þðn,tÞ dgn dt ¼ ½Epu ji : & ¼ Gp @t 0

uðRi ,tÞ ¼

ðB:5Þ

Ro g 0 ðzÞgðzÞ þRo f ðz þ 2Ro Þf ðz þ 2Ro Þ

Now, integrating in the sense of distributions the inner time integral, then commuting the space integrals, one obtains: Z TZ Z tZ u u @r @ wj ðnÞvk ðtÞ ‘Þ ij ¼ wpi ðxÞvph ðtÞ ½Eðup @nn @t r Gu 0 Gp 0 @ d½ttr dgn dt dgx dt þ ½wuj ðnÞvuk ðtÞ @t 4pr

Z TZ Z tZ @r d½ttr 1 @ ¼ wuj ðnÞvph ðtÞ ½wp ðxÞvuk ðtÞ @nn 4pr r @t i Gp 0 Gu 0 @2 þ 2 ½wpi ðxÞvuk ðtÞ dgx dt dgn dt: @t

ftt frr ¼ 0,

gðRi þ tÞ þ f ðRi tÞ ¼ 0,

so that we only need to compute g from (B.6). Moreover, having set DR :¼ Ro Ri , from (B.6):

Remembering (2.6) and (2.9), one has:  Z TZ Z tZ @ @r d½ttr ‘Þ ½Eðup ij ¼ wpi ðxÞvph ðtÞ 4pr Gu Gp @t @nn 0 0 !# u u wj ðnÞvk ðtÞ @ þ ½wuj ðnÞvuk ðtÞ dgn dt dgx dt  @t r Z TZ Z tZ u u @r wj ðnÞvk ðtÞ ¼ wpi ðxÞvph ðtÞ @nn r Gu 0 Gp 0

@ @ d ½t t r dgn dt dgx dt þ ½wuj ðnÞvuk ðtÞ @t @t 4pr Z TZ Z tZ u u @r wj ðnÞvk ðtÞ wpi ðxÞvph ðtÞ ¼ @nn r Gu 0 Gp 0

  @ @ d½ttr dgn dt dgx dt: þ ½wuj ðnÞvuk ðtÞ  @t @t 4pr

2 utt nu ¼ utt urr  ur ¼ 0 r

the two boundary conditions can be rewritten as

ðB:4Þ

gðzÞ ¼ R2o ½eðz=Ro 1Þ 1,

gðRo Þ ¼ 0

Ro o z o Ro þ DR:

ðB:11Þ ðB:12Þ

The wave travels towards the inner surface and reaches Ri at t ¼ DR ¼ Ro Ri . During the interval DR o t o 2DR, from (B.7) we can get f. The condition on the outer surface is still respected by the same g(z) as before, which can be hence extended: gðzÞ ¼ R2o ½eðz=Ro Þ1 1,

Ro þ DR oz o Ro þ 2DR:

ðB:13Þ

We move then to 2DR ot o 3DR. From (B.8), imposing the continuity of g one gets, for Ro þ 2DR o z oRo þ 3DR: gðzÞ ¼ 2R2o Ro eð12DR=Ro þ z=Ro Þ ½4DR3Ro Ro eð2DR=Ro Þ þ 2z ðB:14Þ and so on. For 3DR o t o4DR the same function still satisfies the boundary conditions. Globally it can be seen that function g can be computed piecewise for any Ro þnDR oz o Ro þ ðn þ 1ÞDR, n Z 0, from (B.8) and the integration can be easily automated. References [1] Aimi A, Diligenti M. A new space-time energetic formulation for wave propagation analysis in layered media by BEMs. Int J Numer Meth Eng 2008;75:1102–1132. [2] Aimi A, Diligenti M, Guardasoni C, Mazzieri I, Panizzi S. An energy approach to space–time Galerkin BEM for wave propagation problems. Int J Numer Meth Eng 2009;80:1196–1240. [3] Aimi A, Diligenti M, Panizzi S. Energetic Galerkin BEM for wave propagation Neumann exterior problems. CMES 2010;58:185–219. [4] Aimi A, Diligenti M, Guardasoni C. On the energetic Galerkin boundary element method applied to interior wave propagation problems. J Comp Appl Math 2011;235:1746–1754. [5] Aliabadi MH, Wrobel LC. The boundary element method. New York: John Wiley & Sons; 2002. [6] Bamberger A, Ha Duong T. Formulation variationnelle espace-temps pour le calcul par potentiel retarde´ de la diffraction d’une onde acoustique. I. Math Meth Appl Sci 1986;8:405–435. [7] Bamberger A, Ha Duong T. Formulation variationnelle pour le calcul de la diffraction d’une onde acoustique par une surface rigide. Math Meth Appl Sci 1986;8:598–608. [8] Beskos DE. Boundary element methods in dynamic analysis: part II (1986– 1996). Appl Mech Rev 1997;50:149–197. [9] Costabel M. Time-dependent problems with the boundary integral method. In: Stein E, de Borst R, Hughes TJR, editors. Encyclopedia of computational mechanics. Wiley; 2004.

A. Aimi et al. / Engineering Analysis with Boundary Elements 36 (2012) 1756–1765

[10] Chen JT, Hong HK. Review of dual boundary element methods with emphasis on hypersingular integrals and divergent series. Appl Mech Rev, ASME 1999;52(1):17–33. [11] Curran DAS, Cross M, Lewis BA. Solution of parabolic differential equations by the boundary element method using discretization in time. Appl Math Model 1980;4:398–400. [12] Dunavant DA. High degree efficient symmetrical gaussian quadrature rules for the triangle. Int J Numer Meth Eng 1985;21:1129–1148. [13] Falletta S, Monegato G, Scuderi L. Space–time BIE methods for non homogeneous exterior wave equation problems. The Dirichlet case. IMA J Numer Anal 2011:1–25. [14] Frangi A. Elastodynamics by BEM: a new direct formulation. Int J Num Meth Eng 1999;45:721–740. [15] Frangi A, Novati G. On the numerical stability of time-domain elastodynamic analyses by BEM. Comput Meth Appl Mech Eng 1999;173:403–417. [16] Frangi A. ‘‘Causal’’ shape functions in the time domain boundary element method. Comp Mech 2000;25:533–541. [17] Hammer PC, Stroud AH. Numerical integration over simplexes. Math Tables other Aids Comput 1956;10(55):137–139. [18] Ha Duong T. On retarded potential boundary integral equations and their discretization. In: Davies P, Duncan D, Martin P, Rynnm B, editors. Topics in computational wave propagation. Direct and inverse problems. SpringerVerlag; 2003. p. 301–336.

1765

[19] Lubich C. On the multistep time discretization of linear initial-boundary value problems and their boundary integral equations. Numer Math 1994;67:365–389. [20] Mansur WJ, Brebbia CA. Further developments on the solution of the transient scalar wave equation. In: Brebbia CA, editor. Topics in boundary elements research, vol. 2. Berlin: Springer-Verlag; 1985. p. 87–123. [21] Milroy J, Hinduja S, Davey K. The elastostatic three-dimensional boundary element method: analytical integration for linear isoparametric triangular elements. Appl Math Model 1997;21:763–782. [22] Monegato G, Scuderi L, Stanic MP. Lubich convolution quadratures and their application to problems described by space–time BIEs. Numer Algor 2011;56(3):405–436. [23] Schanz M, Antes H. Application of operational quadrature methods in time domain boundary element methods. Mechanica 1997;32:179–186. [24] Schanz M. Application of 3-D boundary element formulation to wave propagation in poroelastic solids. Eng Anal Boundary Elem 2001;25:363–376. [25] Banjai L, Messner M, Schanz M, Runge–Kutta convolution quadrature for the boundary element method, Preprints of Institute of Applied Mechanics of Graz University of Technology, 3; 2011. [26] Soares D, von Estorff O, Mansur WJ. Iterative coupling of BEM and FEM for nonlinear dynamic analyses. Comput Mech 2004;34:67–73. [27] Xiao J, Ye W, Cai Y, Zhang J. An efficient frequency-domain approach for transient structural analysis based on the precorrected-FFT accelerated BEM, in IABEM Extended Abstracts; 2011:349–354.