A comparative study on three boundary element method approaches to problems in elastodynamics

In this study, three boundary element method approaches are compared by solving a typical problem in linear elastodynamics. The first approach formulates and solves the problem in the real time domain in conjunction with a time-stepping algorithm. The other two approaches formulate and solve the problem in the Laplace and the Fourier transformed domains, which necessitates a numerical inversion of the results to the time domain. The results obtained compare favourably with available analytic solutions. A detailed tabulation of the computer time and memory requirements of each approach is presented.


INTRODUCTION
Perhaps surprisingly, the use of integral equation formulations in the analysis of transient phenomena in solids and fluids has a long history, going back over one hundred years to the Helmholtz-Kirchhoff integral formula, which is the mathematical interpretation of Huygens' principle. 1 In a large number of these phenomena, part of the boundary is at infinity, which makes the use of integral equation methods particularly attractive, since they are capable of accounting for the presence of infinite domains. Morse and Feshbach, 2 Eringen and Suhubi 3 and Kupradze 4 are but a few of the references that contain details of the classical work done in elastodynamics and related topics.
Although as mentioned above the basic integral equation formulations for wave propagation problems have been known for quite some time, their adaptation for constructing numerical algorithms for use in the solution of boundary value problems is new. Some of the earliest such developments were by Friedman and Shaw 5 and Chen and Schweikert, 6 in acoustics, and by Banaugh and Goldsmith,7 in steady-state elastic wave propagation. Since then, many other investigators have formulated and used numerical implementations known as integral equation methods (IEM), boundary integral equation methods (BIEM) or boundary element methods (BEM) in various fields of engineering. In the field of elastodynamics, Cruse and Rizzo 8 and Cruse 9 derived a BIEM approach in conjunction with the Laplace transform and employed it to solve a half-plane wave propagation problem. A modified version of their approach was used later by Manolis and Beskos 10 to investigate stress concentrations around cylindrical openings of arbitrary shape due to the scattering of elastic waves of a general transient nature. The determination of the steady-state solution by BEM approaches and the reconstitution of the transient response by Fourier synthesis was done by Niwa et al. 11 waves. A time domain formulation and solution for elastic wave scattering problems was done by Cole et al. 13 for the anti-plane strain case and by Niwa et al. 14 for the general two-dimensional plane stress/plane strain cases.
It therefore appears that a study comparing the various BEM approaches used to solve the propagation, diffraction and scattering problems arising in the field of elastodynamics would be of value to the engineering community. In this work, three possible BEM approaches are investigated and compared. All three approaches are used to study the classical problem of the scattering of a pressure wave pulse by a circular cylindrical cavity in a linear elastic, homogeneous and isotropic medium. This is a two-dimensional problem of the plane strain kind. Analytic-numerical solutions for the displacement, velocity and stress fields at the circumference of the cavity were obtained by Baron and Matthews, 15 Baron and Parnes 16 and by Pao and Mow. 17 The former two pairs of investigators used an integral transform, separation on the circumferential angle (8) method. The latter investigators (Pao and Mow) used the method of wave function expansion coupled with the Fourier synthesis technique. Miklowitz 18 used a special technique to bring out hitherto undetected late-time singular, non-decaying Rayleigh waves on the cavity walls due to impulsive loads. It is interesting to note that this classical problem has been resolved recently by other analytical means as, for instance, the method of matched asymptotic expansions discussed by Datta. 19 The three BEM used to resolve this well-documented problem are: 1. A time domain approach that solves the problem as if it were a three-dimensional one. In particular, the plane strain problem is represented as a circular cylinder of infinite length and a numerical integration is performed along the axis of the cylinder so as to collect the contributions to a representative point from points lying along the same generator.
2. An approach in which the integral equation formulation is applied to the Laplace transform of the governing equations of motion for the two-dimensional case. In this case, a numerical inverse transformation is required to bring the transformed solution back to the original time domain.
3. The last approach is similar to case 2, except that the integral transform in question is the Fourier transform instead of the Laplace transform.
All of the aforementioned cases share some common traits. First, the direct BEM is employed, where the unknown boundary tractions and displacements are directly related to the known ones. One could use an indirect BEM, as was done in Niwa et a!. 11 For more details on direct and indirect BEM, reference should be made to Banerjee and Butterfield? 0 Second, all approaches utilize a spatial discretization scheme for which the boundary tractions and displacements are assumed to have constant values over a given boundary segment. Finally, considerable effort is directed towards ensuring that all three approaches share as many common features as possible, in order to render the comparison study as fair as possible.
In the next few sections, the integral equation formulation of the general elastodynarnic problem along with the numerical implementation is discussed. Following that, the results drawn from applying the three aforementioned BEM approaches to the same problem, which is shown in Figure 1, are discussed in detail.

STATEMENT OF THE PROBLEM AND DEFINITIONS
Let V denote the area occupied by a given body and S denote its bounding surface. A co-ordinate system (x, t) is employed, where x denotes the cartesian spatial co-ordinates x 1 and x 2 and tis the time. Associated with each surface point is an outward pointing unit normal vector n. Under the assumption of small displacement theory and linear elastic isotropic and homo~eneous material behaviour, the equations of motion of the body are (1) where u;(x, t) is the displacement vector and[; is the body force vector per unit mass. The propagation velocities of the pressure and shear waves in the body are given, respectively, as ci=(A+2M)/p; c~=IJ.IP (2) where A and IJ. are the Lame constants and p is the mass density. Furthermore, commas indicate space differentiation, dots indicate time differentiation and repeated indices imply the summation convention. Since in what follows attention is restricted to two-dimensional cases, the indices i and j assume the values of 1 and 2, unless mentioned otherwise. The displacement vector uJx, t) is assumed to be twice differentiable except at possible surfaces of discontinuity, as discussed in Eringen and Suhubi. 3 The displacements and the tractions t; (x, t) satisfy the boundary conditions I;= U;;n; = p;(X, t) where u;; is the stress tensor and S = S 1 + Su.
The displacements and velocities also satisfy the initial conditions In addition, the Sommerfeld radiation condition must be satisfied if the body in question is of infinite extent.  (5) with 8ii being the Kronecker delta.

The Laplace transform
The Laplace transform of a function f(x, t) is defined as where s is the Laplace transform parameter and f(x, t) must be at least piecewise continuous in time. A Laplace transform of the equations of motion (1), which are hyperbolic partial differential equations, results in the following system of elliptic partial differential equations: which is more amenable to numerical solutions. It should be noted in passing that the Navier equations of equilibrium are also of elliptic nature, i.e. they are similar in form to equation (7). Since the boundary conditions and the constitutive equations do not involve time derivatives, their Laplace transforms are very simple and become, respectively, and ui=qi(x,s) ii"·· = p [(c 2 1-2c2 2 )u ~--+ c 2 2(u· · + u · ·)] I] r,, UIJ 1,] ],1 Essentially, the solution to a transient elastodynamic problem when integral transforms are employed, such as the Laplace or Fourier transform, consists of a series of solutions to a static-like problem for a number of discrete values of the transformed parameters. This series of solutions must finally be numerically inverted back to the original time domain.
In what follows, a quiescent state is assumed before time t = o+, which implies that the initial conditions are zero. Also, the body forces are taken equal to zero as well. These assumptions are made for the sake of convenience and no additional conceptual difficulty in the solution procedure is encountered if these assumptions are relaxed.

The Fourier transform
The Fourier transform of a function f(x, t) is defined as ic x, w) = .'F{f(x, t)} = L: t(x, t) e -<w• dt (10) where w is the frequency.
For the quiescent initial conditions assumed previously, it is elementary to see that the Fourier transform of the governing equations of motion ( 1) will be identical to equation (7) if the Laplace parameter s is set equal to iw in that equation. Therefore, the BEM solution in the Fourier domain is identical to the solution scheme that is to be employed for the Laplace transform case. The difference comes in the numerical inversion from the Fourier transform domain to the time domain, and this subject is elaborated in the next section.

INTEGRAL EQUATION FORMULATION
It is well known (DeHoop/ 1 Wheeler and Sternberg 22 and Kupradze 4 ) that a system of partial differential equations along with the appropriate boundary and initial conditions may be cast in an integral equation form. This holds true for either the system of equations (1)- (5) or the system of equations (7)-(9), depending on whether one wishes to establish an integral equation formulation in either the time domain or in a transformed domain. Various mathematical questions regarding the validity of such representations in elastodynamics and the existence and uniqueness of solutions that may be obtained can be found in the aforementioned references as well as in Sternberg and Eubanks. 23 The basis for an integral equation formulation in elastodynamics is the dynamic extension of Betti's reciprocal theorem. This theorem is derived from virtual work considerations and essentially relates two different dynamic states for the same elastic body. The obvious choice for one of the two elastodynamic states is the unknown state we are seeking to find. A judicious choice for the other state is the appropriate Green's function that satisfies the governing equations of motion in the region V. For more details, an appropriate reference is Banerjee and Butterfield. 20 After some manipulations, the dynamic equivalent to Somigliana's identity can be obtained and the limiting form of this identity results in the singular integral equations that form the basis of the BEM method.

Time domain integral equation formulations
In the time domain, the Green's function is the fundamental point force solution of Stokes and the singular integral equations are of the form where ~~ is the source point and X 1 is the receiver point. In (11), the prime symbol used as a superscript denotes three-dimensional quantities with the pertinent subscripts ranging from 1 to 3. Also, If the boundary S 1 is not smooth at~~, c;i becomes a function of the angle e at point~~. The operator* denotes a convolution in time, that is The Green's function o;i represents the displacement vector at X 1 and at time t due to a unit impulse of the form 8 (t-r ) . 8 (x 1 -~' ) applied in the three principal directions at ~, and time r, with 8 being the Dirac delta function. Usually r is set equal to zero. Finally, F ;i represents the resulting tractions due to the same unit impulse and can be obtained from G~i through the use of the constitutive equations and surface normals.
The kernels G ;i and F;i are given in Appendix I, along with details of the solution procedure.
The approach delineated there is essentially a three-dimensional formulation applied to an infinitely long cylindrical body whose longitudinal axis coincides with the Z axis of the cartesian co-ordinate system ( Figure 2). The problem, of course, remains a plane strain one. The infinitely long cylindrical body may be viewed as a cavity if an interior problem is under consideration.  (14) where S is now the circumference of the two-dimensional body in question. As far as the integration over time is concerned, it appears satisfactory to assume that Ui and ti remain constant over the time step [to+ ilt, t 0 ].

Laplace transform domain integral e'quation formulation
In the Laplace transform domain, the Green's function can be found in Cruse and Rizzo 8 and Nowacki 24 and the singular integral equations are of the form The explanatory notes following equation (11) apply here as well. It should be reminded at this point that bars over a quantity denote the Laplace transform of that quantity. The tensors G;; and F;; for the case of plane strain or plane stress are given in Appendix III. The solution procedure in the Laplace transform domain is elaborated in detail in Manolis and Beskos 10 and only a few remarks will be made here. The problem is essentially static-like for a fixed value of the parameter s. The solution procedure consists of solving for the unknown boundary values of u; and 0 in terms of the prescribed boundary conditions with the aid of equation (15), for a series of values of the parameter s. The final step is to numerically invert the transformed solution to the time domain. The pertinent inversion integral is where {3 > 0 is arbitrary but greater than the real part of all the singularities of f(x, s ).
Equation (16) is an integral over the complex plane, which implies that for accuracy purposes equation (15) must be solved for complex values of the parameters, that is In equation (17)  f(x, 1) = t: f(x, W) e' 2 ww< dw (18) For this case, the fast Fourier transform algorithm of Cooley and Tukey, 26 which has been extensively used by numerous investigators in the past is employed. Since the upper and lower limits in (18) are replaced by 0 and !1 (the maximum frequency of interest), the well-known problem associated with the finite Fourier transform, namely the approximation a non-periodic motion by a periodic one, arises. This problem is circumvented by adding a sufficient number of "trailing' zeros to the transformed solution.
Finally, one should be aware of a difficulty associated with the Fourier transform solution. Equation (15) fails to render a solution for an exterior (or interior) problem for a particular set of frequencies w corresponding to the eigenvalues of the associated interior (or exterior) problem. The original boundary value problem possesses, of course, a unique solution at those frequencies. This problem was noticed by workers in the field of acoustics and a number of methods were introduced to combat it, such as the technique of Schenk. 27 In the particular example used here, this numerical difficulty was not encountered. It is possible, however, for problems in elastodynamics to isolate the eigenvalues of the associated interior (or exterior) problem by methods such as the one described in Tai and Shaw 28 and then to modify the solution procedure of the original exterior (or interior) problem in the spirit of the techniques developed in acoustics so as to avoid the aforementioned difficulty. Some details on this subject can be found in Kobayashi and Nishimura. 29

NUMERICAL IMPLEMENTATION
It is obvious that equations (11), (14) and (15) cannot, in general, be solved analytically and therefore resort must be made to numerical methods of solution. The bounding surface S (circumference) of the two-dimensional body is therefore discretized into J segments with a node defined at the midpoint of a given segment, as shown in Figure 3. The spatial variation of the displacements and tractions over a given segment is assumed to be constant, although it is possible to define linear or higher order variations.
The time axis is now discretized into L segments so that the total time of interest is T = L6.t and a constant temporal value of the displacements and tractions over a time interval is adopted. Equation (19) can be recast in the following convenient for numerical computations form: with ~S being the length of segment i. Essentially, the difference between the two aforementioned time domain approaches is that in the former the integration is done over area segments that move with time, while in the latter the integration is done over a line segment and over a time step. The advantage of using a three-dimensional approach is that singularities in the kernels are encountered only once in the time marching scheme, when n = m = 1 and for i = j. These singularities are discussed in Cole et al. 13 If a two-dimensional approach were to be tlsed, then a singularity would be encountered at every time step m when i = j. This comment becomes obvious when reference is made to the radius r in Figures 2 and 3.
In the integral transform formulations, as described in Reference 10, the static-like system of equations 0·5iij = l:DGjii -l:DF;;ii; (22) i i is solved for a number of values of the transformed parameter s. The additional price that is paid is that a numerical inverse transformation is required to bring the tractions and the displacements back to the time domain. The singularities encountered in these formulations, when i = j, are of the form In r for the G;; tensor and of the form r, ;nJI r for the Fi; tensor. The behaviour of these singularities as r ~ 0 is well documented. 30 In all the aforementioned formulations, 4-or 6-point Gaussian quadrature formulae are used for the non-singular space and time integrations. In the case of singularities, the interval of integration of the singular segment is broken into two parts with the singularity (r = 0) excluded and Gaussian quadrature is employed, while the principal value of the integrand is found analytically as r-+ 0. The matrix inversions are done by Gauss elimination with Crout's algorithm for real cases and by Cholesky's decomposition for complex cases.
As a final note, it should be mentioned that considerable economy in the solution scheme of equation (20) can be achieved if advantage is taken of the time translation property, equation (29), of the integrands G ~i and F:i . . This property implies that at time step n llt only the matrices DGijm and DFijm need be computed, while the remaining matrices DGijm -\ DG ':t --2 , ••. , DG ij 1 (similarly for DF~t) need not be recomputed. The time-stepping algorithm for solving equation (20) becomes an explicit algorithm if the time step is selected from the relation c 1 tl.t/ L ~ 1·0, where L is a typical segment length.

NUMERICAL EXAMPLE
This section describes the detailed solution of a numerical example that serves as the vehicle for a comparison study of the three BEM approaches.
Consider an infinitely extending linear elastic medium with a circular cylindrical cavity of radius a and of infinite length under the action of a compressional plane shock (P) wave whose front is parallel to the axis of the cavity, as shown in Figure 1. For this wave scattering problem, analytic-numerical solutions exist (Baron and Matthews, 15 Pao and Mow 17 ) and the state of stress around the cavity wall is determined numerically by the three proposed approaches for comparison purposes. The problem is of the plane strain kind and the P-wave results in an applied load tensor Ux = S 0 H(ttn), Uy = (v/(1 -v ))SoH(t-t") and 1'.-cy = 0 as it starts enveloping the cavity at the generic time t = 0. In the expression for the load tensor, H is the Heaviside function, v is Poisson's ratio and tn is the time required for the wave moving with velocity c 1 and starting at point (x 1 = a, x 2 == 0) to reach a point n on the circular boundary.
The resulting stress distribution in the medium is obtained by superimposing the stress field produced by the P-wave in the medium without the hole to the stress field produced by the application of corrective tractions on the boundary of the cavity in order to render the cavity surface traction-free. It is obvious that the first problem is trivial, so that attention is confined to the second problem, which is solved by the BEM. The following numerical values for the constants of the problem are used: These material properties correspond physically to granite and lead with the aid of equation (2) to the following values for the propagation velocities: c 1 = 208,000 in/sec, c2 = 120,000 in/sec (24) In addition, Poisson's ratio becomes 0·25 and the time required for the wave to travel the cavity diameter, the full transit time 27', is equal to 2aj c 1 = 0·002 sec. The intensity of the applied load So is conveniently taken as equal to -1·0 lb/in. 2 In all cases, the boundary of the circular hole is discretized into 20 equal segments. At this point it should be mentioned that the time domain formulation employing the two-dimensional kernels Gi; and Fi; given in Appendix II was also used but was abandoned for three reasons: 1. The difficulties encountered in the analytic determination of the singularities for every time step m. 2. The inefficiency of numerically integrating the kernels because of the presence of the Heaviside function. Note how easy the time integration is in the three-dimensional case in Appendix I because of the presence of the delta function. For instance, for a time step llt = 0·00015 sec and N = 6 time steps, the time domain formulation with the two-dimensional kernels required 208·6 sec of total central processor unit (CPU) time and cost $10·8, while the formulation with the three-dimensional kernels required 75·5 CPU sec and cost $4·2. 3. The numerical results obtained were not good.
It can be observed from Figures 4 and 6 that the results obtained from the three BEM formulations agree reasonably well with the results of analytic-numerical approaches. What may be considered unusual is that the two analytic-numerical solutions do not coincide, but this is not surprising since in both a large amount of numerical work, accompanied by some assumptions, is involved. In particular, in Figure 4 the BEM results appear to be bounded by Pao's solution from above and by Baron and Matthews' solution from below until time t =ST. The maximum value for the stress concentration a 66 / S 0 occurs at this location and the analytic-numerical solutions predict a 10 per cent overshoot over the static value of 2·67, while the BEM solutions predict a 15 per cent overshoot. These maximum values occur between t = 4T and t = 6T.
In Figure 5, only the results from the integral transform BEM show reasonably good agreement with the analytic-numerical predictions, while the time domain BEM results appear to be somewhat low. The true behaviour at this location is probably better represented by the results obtained by Baron and Matthews, because when the P-wave pulse first reaches the cavity, the points around(}= oo act as a rigid wall and as a consequence the applied stress iiy should double. At later times, it is reasonable to expect the solution to drop to very low values, since the static stress concentration factor is zero at(}= 0°. The reason the transformed BEM results (and Pao's as well) do not exhibit this sudden rise in stress at t = o+ is because the application of the Laplace and Fourier transforms smoothens the vertical front of the P-wave pulse. The numerical implementation of the time domain BEM also tends to smoothen the sharp wave front to a ramp-like rise over a time interval of a few llt. As a consequence, the sudden rise in a 88 does not materialize at t = o+. The possibility, however, of treating the sudden jump in stress carried by the P-wave pulse in a more satisfactory manner in the time domain BEM formulation exists. This would require the addition of a displacement term to equation (11) resulting from an analytic integration of the time derivatives of the kernels appearing in (11) along the path the wave front traces in time. An interesting discussion of this problem in relation to the acoustic case appears in Reference 5. The results of such an attempt will be communicated at a later date. An additional approximation introduced in the BEM formulations is that the stresses at the boundary involve a simple central finite difference scheme for the displacements that is bound to introduce a small error. Finally, Table I presents the memory and time requirements of the three BEM approaches. It should be noted there that the difference between the execution time and the total time is the time required for the programs to compile. For comparable accuracy, it appears that the time domain program is more expensive by a factor of 5 and 3 than the Laplace and Fourier transform domain programs, respectively. Also, it should be mentioned that all computations were performed on a CDC Cyber 176 computing machine.

CONCLUSIONS
On the basis of the preceding presentation, the following conclusions can be drawn: 1. Three general BEM formulations for solving problems in transient elastodynamics involving bodies of arbitrary shape under conditions of plane stress or plane strain are presented. One of the formulations is in the time domain, while the remaining two involve the Laplace and Fourier transforms with the additional task of the numerical inversion of the transformed solution to the time domain. 3. All three formulations are viable alternatives for the solution of general elastodynamic problems. From the present study it may be concluded that, for comparable accuracy, the Laplace transform domain BEM formulation is more economical than either the Fourier transform domain or the time domain formulations by a factor of 1·7 and 5·0, respectively. However, the time domain formulation is useful because it gives a better picture for very early times up to 1 full transit time. Integral transform formulations tend to reach a plateau corresponding to the solution of the static equivalent of the problem, while time domain formulations tend to diverge after a very large number of time steps. 3. It should be noted that all three BEM approaches can result in better estimates with more refined boundary discretization and more refined time or transform parameter sampling. Also, the use of isoparametric shape functions for the description of the both geometry and boundary data, as well as of a linear time variation of the boundary quantities in the time domain BEM, will certainly result in improvements. 4. Generally speaking, the problem with time domain BEM is that the time step must be kept a small fraction of the total time interval of interest, which necessitates the use of a relatively large number of time steps. Unfortunately, the solution algorithm becomes more involved at later time steps because the information from all earlier time steps must be used. The basic problem with integral transform BEM is that the solution procedure is performed in complex arithmetic, which requires at least twice the number of operations than a comparable procedure done in real arithmetic. 5. It will be observed that the BEM solutions are not inexpensive. To circumvent this economic difficulty, one may employ, in the beginning stages of an analysis/design procedure, simplified methods of analysis and resort to the powerful BEM approaches at the end of the procedure. It should be noted that at present, and at least as far as the field of elastostatics is concerned where the BEM has reached maturity, numerical integral equation methods are competitive with finite element and finite difference methodologies. 6. Some further advantages of the integral transform BEM are that the steady-state (harmonic problem) becomes a special case of the Fourier domain BEM not requiring a numerical inversion from the transformed domain and that viscoelastic material behaviour can be recovered from the Laplace domain BEM by a simple change in the elastic constants A and 1-L according to the correspondence principle. 7. Some further advantages of the time domain BEM are that it is more suitable for extension to nonlinear material behaviour through the use of incremental and iteration schemes and that the extension to three-dimensional problems is straightforward.
Foundation of the State University of New York is acknowledged, as are the University Computing Services of the SUNY at Buffalo for making their facilities available to the author.

APPENDIX I
The general three-dimensional form of the tensors G ;i and F;i can be found in Eringen and Suhubi. 3 In particular, and (26) where (27) In the above expressions, r; = (x! -~D, r = lx'-~'I and the subscripts range from 1 to 3. The Green's function a;i(x', t; ~', T) obeys the causality condition (28) and has the following time translation property: (29) The following property of the derivatives of the Dirac delta function will be used shortly: b L 8<n(e-d)f(e) de = ( -l)uy(i)(d) (30) In addition, it should be noted that r,; = rj r.
The tensors o;i and F;i are now cast into a form appropriate for the case of plane strain.
Consider a cylindrical body V' + S' of infinite extent along the Z direction, as shown in Figure   2. At time t a wave that emanated from source point ~' = ~ envelops a spherical region of radius R = ct, where c is the appropriate propagation velocity. Notice that f is taken to lie on the XY plane for convenience, since both displacement and velocity vectors are now independent of the Z co-ordinate. Therefore, at time t the first of the two integrals appearing in equation (11) where the kernel function G~j is the same as G~i except that the delta functions in (25) have been used in accordance with (30) to determine the limits of integration along the Z axis. In (31), the current time t = N D..t with ~t being an appropriate time increment and x is the projection of x' on the XY plane. Furthermore, Zk is the vertical distance x 3 between x' and x at time KAt, that is (32) Therefore, the integration in (31) takes place over area patches whose horizontal ordinates are along S and whose vertical ordinates are the distances ZK that mark the location of the spherical wave originating from~ at time increments K !:it, K = 1, 2, ... , N.
Assuming In equation (34 ), the factor 2 comes from the fact that the spherical wave propagates both above and below the XY plane. Additionally, it should be noted that r is the distance between x' and~ and that indices i and j assume the values 1 and 2 only. The above equation can also be found in Niwa et al.