``Causal'' shape functions in the time domain boundary element method A. Frangi Abstract Since failing to respect the causality condition has been identi®ed as one of the main sources of inaccu- racies in the time domain boundary element method for elastodynamics and scalar wave propagation problems, in this contribution new shape functions are investigated, which permit a more accurate simulation of the continuous propagation of wave fronts. The performance of these shape functions in 2D scalar wave propagation problems is tested both for the potential (displacement) and for the time gradient (velocity) equations. Analytical time inte- grations are developed and numerical results are presented. 1 Introduction In the literature on the time domain boundary element method (BEM) for elastodynamics and scalar wave prop- agation problems (see the survey article by Beskos, 1997) great attention has been devoted to the choice and im- plementation of different interpolation schemes for dis- placements (or potential functions) and tractions (or normal derivatives of potential functions): combinations of constant, linear and quadratic functions (e.g. Domin- guez, 1993; Mansur et al., 1999) or more advanced solutions adopting B-splines (Rizos and Karabalis, 1994) are but few examples. However, to the author's best knowledge, all these options have a common feature: space±time interpolations are expressed by means of ``product'' shape functions, e.g. in 2D, the generic dis- placement component us; t is approximated as: us; t h m t g n sU mn , where t is time and s is the curvilinear coordinate running along the boundary C of the body X. Probably, the main advantage of this class of interpo- lation functions is simplicity stemming from the intrinsic decoupling of space and time coordinates which makes analytical integration with respect to time (within the convolution integrals) particularly straightforward. How- ever, each function h m t  just modulates the space dis- tribution g n sU mn , without modifying its domain of de®nition (e.g., for linear hat-shaped functions, the two elements adjacent to the space node). If a wave propagates along the boundary (the wave front having a continuous motion), this interpolation proves inaccurate and does not respect, in general, the causality condition, forcing the disturbance to propagate in a non-physical manner. This aspect is of crucial importance since, as reported in Mansur et al. (1998) ``the most important global charac- teristic to be well represented in a numerical time domain analysis is causality''. Similar conclusions have been re- cently presented in Frangi and Novati (1999) and Frangi (1999). In order to visualize a situation in which a suitable simulation of wave front propagation is strongly advo- cated, let us consider a simple benchmark which has been extensively investigated in the literature on time domain BEMs (e.g. Mansur, 1983; Dominguez, 1993): a strip of height L  1 in antiplane-shear conditions (Fig. 1), ex- tending inde®nitely in the x 1 (and x 3 ) direction, is sub- jected on the upper surface to a uniform shearing traction s x 2 x 3  P  Ht, where H is the Heaviside step function. In the right side of Fig. 1 are presented the exact (``cyclic'') non-dimensional displacements u a  ul=PL (u  u 3 , l shear modulus) at B and C and tractions p a  p=P (p  s x 2 x 3 ) at A (l is the shear modulus). The problem is essentially one-dimensional in space, since all the variables involved depend only on the x 2 coordinate. A plane elastic wave keeps travelling (and thus re¯ecting) inde®nitely between the upper and lower surfaces with constant speed c (c 2  l=q, q material density); accelera- tions are concentrated along the wave front where they behave like a Dirac's delta distribution. Let us consider a time knot  t such that 0 <  t < c=L, i.e. before the elastic wave reaches the constrained lower side for the ®rst time. The displacement time history along each of the vertical sides (dashed lines in Fig. 1a) in the interval considered is graphically represented as a plane in Fig. 2 (s is the cur- vilinear coordinate running along these sides; s  0 at the free edge of the strip). For each s the intersection of the plane with t  s yields displacements at s along the ver- tical sides. Moreover, its intersection with the u  0 plane denotes the wave front propagating along the vertical sides as time progresses. This latter line has constant slope c. It is worth stressing that no linear combination of product shape functions could represent the plane in Fig. 2 exactly. The ®nite strip delimited by the dashed lines in Fig. 1 can be analysed as a two-dimensional body employing the customary BE formulation for 2D wave propagation which can be written as (Dominguez, 1993): Cyuy; t Z C U y; x; t  ? px; t  T y; x; t  ? ux; t dC 1 Computational Mechanics 25 (2000) 533±541 Ó Springer-Verlag 2000 533 A. Frangi Department of Structural Engineering, Politecnico di Milano P.za Leonardo da Vinci 32, 20133 Milano, Italy