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

George D. Manolis

► To cite this version:

George D. Manolis. A comparative study on three boundary element method approaches to problems in elastodynamics. International Journal for Numerical Methods in Engineering, 1983, 19 (1), pp.73 - 91. 10.1002/nme.1620190109. hal-01382196

HAL Id: hal-01382196

https://hal.science/hal-01382196

Submitted on 16 Oct 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A COMPARATIVE STUDY ON THREE BOUNDARY ELEMENT METHOD APPROACHES TO PROBLEMS IN ELASTODYNAMICS

GEORGE D. MANOLIS

Department of Civil Engineering, State University of New York at Buffalo, Buffalo, New York, U.S.A.

## SUMMARY

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.¹ 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,² Eringen and Suhubi³ and Kupradze⁴ 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⁵ and Chen and Schweikert,⁶ in acoustics, and by Banaugh and Goldsmith,⁷ 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⁸ and Cruse⁹ 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¹⁰ 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.¹¹,¹² in the course of their investigations of the transient stresses produced around cavities during the passage of travelling
waves. A time domain formulation and solution for elastic wave scattering problems was done by Cole *et al.*¹³ for the anti-plane strain case and by Niwa *et al.*¹⁴ 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,¹⁵ Baron and Parnes¹⁶ and by Pao and Mow.¹⁷ The former two pairs of investigators used an integral transform, separation on the circumferential angle (θ) method. The latter investigators (Pao and Mow) used the method of wave function expansion coupled with the Fourier synthesis technique. Miklowitz¹⁸ 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.¹⁹

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 al.*¹¹ For more details on direct and indirect BEM, reference should be made to Banerjee and Butterfield.²⁰ 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 elastodynamic 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 t is the time. Associated with each surface point is an outward pointing unit normal vector **n**.
<figure><img src="image_2.png" /><figcaption>Figure 1. Circular cavity in an infinite medium under the action of an incoming pressure wave</figcaption></figure>

Under the assumption of small displacement theory and linear elastic isotropic and homogeneous material behaviour, the equations of motion of the body are

$$ (c_1^2 - c_2^2)u_{i,ij} + c_2^2 u_{j,ii} + f_j = \ddot{u}_i \quad (1) $$

where $u_i(\mathbf{x}, t)$ is the displacement vector and $f_j$ is the body force vector per unit mass. The propagation velocities of the pressure and shear waves in the body are given, respectively, as

$$ c_1^2 = (\lambda + 2\mu)/\rho; \quad c_2^2 = \mu/\rho \quad (2) $$

where $\lambda$ and $\mu$ are the Lame constants and $\rho$ 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 $u_i(\mathbf{x}, t)$ is assumed to be twice differentiable except at possible surfaces of discontinuity, as discussed in Eringen and Suhubi.³

The displacements and the tractions $t_i(\mathbf{x}, t)$ satisfy the boundary conditions

$$ \begin{aligned} t_i &= \sigma_{ij} n_j = p_i(\mathbf{x}, t) \quad \mathbf{x} \in S_t \\ u_i &= q_i(\mathbf{x}, t) \quad \mathbf{x} \in S_u \end{aligned} \quad (3) $$

where $\sigma_{ij}$ is the stress tensor and $S = S_t + S_u$.

The displacements and velocities also satisfy the initial conditions

$$ \begin{aligned} u_i(\mathbf{x}, 0^+) &= u_{i0}(\mathbf{x}) \quad \mathbf{x} \in V + S \\ \dot{u}_i(\mathbf{x}, 0^+) &= \dot{u}_{i0}(\mathbf{x}) \quad \mathbf{x} \in V + S \end{aligned} \quad (4) $$
In addition, the Sommerfeld radiation condition must be satisfied if the body in question is of infinite extent.

Finally, the constitutive equations are of the form

$$ \sigma_{ij} = \rho [(c_1^2 - 2c_2^2)u_{r,i} \delta_{ij} + c_2^2(u_{i,j} + u_{j,i})] \quad (5) $$

with $\delta_{ij}$ being the Kronecker delta.

## The Laplace transform

The Laplace transform of a function $f(x, t)$ is defined as

$$ \tilde{f}(\mathbf{x}, s) = \mathcal{L}\{f(\mathbf{x}, t)\} = \int_{0}^{\infty} f(\mathbf{x}, t) e^{-st} dt \quad (6) $$

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:

$$ (c_1^2 - c_2^2)\bar{u}_{i,ij} + c_2^2\bar{u}_{j,ii} + \bar{f}_j - s^2\bar{u}_j + \dot{u}_{j0} + su_{j0} = 0 \quad (7) $$

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,

$$ \bar{t}_i = \bar{\sigma}_{ij}n_j = \bar{p}_i(\mathbf{x}, s) \quad \mathbf{x} \in S_t \quad (8) $$

$$ \bar{u}_i = \bar{q}_i(\mathbf{x}, s) \quad \mathbf{x} \in S_u $$

and

$$ \bar{\sigma}_{ij} = \rho [(c_1^2 - 2c_2^2)\bar{u}_{r,i} \delta_{ij} + c_2^2(\bar{u}_{i,j} + \bar{u}_{j,i})] \quad (9) $$

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 parameter $s$. 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 = 0^+$, 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

$$ \tilde{f}(\mathbf{x}, w) = \mathcal{F}\{f(\mathbf{x}, t)\} = \int_{-\infty}^{\infty} f(\mathbf{x}, t) e^{-iwt} dt \quad (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,²¹ Wheeler and Sternberg²² and Kupradze⁴) 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.²³

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.²⁰

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

$$c'_{ij}u_i(\xi', t) = \int_{S'} \{G'_{ij}(\mathbf{x}', t; \xi') * t_i(\mathbf{x}', t) - F'_{ij}(\mathbf{x}', t; \xi') * u_i(\mathbf{x}', t)\} dS'(\mathbf{x}') \quad (11)$$

where $\xi'$ is the source point and $\mathbf{x}'$ 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,

$$c'_{ij} = \begin{cases} \delta_{ij} & \text{for } \xi' \in V' \\ 0.5 \delta_{ij} & \text{for } \xi' \in S' \\ 0 & \text{for } \xi' \notin V' + S' \end{cases} \quad (12)$$

If the boundary $S'$ is not smooth at $\xi'$, $c'_{ij}$ becomes a function of the angle $\theta$ at point $\xi'$. The operator $*$ denotes a convolution in time, that is

$$f(\mathbf{x}', t) * g(\mathbf{x}', t) = \begin{cases} 0 & \text{for } t \le 0 \\ \int_0^t f(\mathbf{x}', \tau)g(\mathbf{x}', t-\tau) \, d\tau & \text{for } t > 0 \end{cases} \quad (13)$$

The Green's function $G'_{ij}$ represents the displacement vector at $\mathbf{x}'$ and at time $t$ due to a unit impulse of the form $\delta(t-\tau) \cdot \delta(\mathbf{x}'-\xi')$ applied in the three principal directions at $\xi'$ and time $\tau$, with $\delta$ being the Dirac delta function. Usually $\tau$ is set equal to zero. Finally, $F'_{ij}$ represents the resulting tractions due to the same unit impulse and can be obtained from $G'_{ij}$ through the use of the constitutive equations and surface normals.
The kernels $G'_{ij}$ and $F'_{ij}$ 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.

<figure><img src="image_3.png" /><figcaption>Figure 2. Circular cylindrical cavity of infinite length—region affected by wave propagating from point $\xi$</figcaption></figure>

There is an alternative path for solving two-dimensional problems in the time domain by integral equations. This simply requires the integration of the tensors $G'_{ij}$ and $F'_{ij}$ along the longitudinal axis so that their true two-dimensional forms $G_{ij}$ and $F_{ij}$ are obtained, as shown in Appendix II. Equation (11) then changes to

$$c_{ij}u_i(\xi, t) = \int_0^t \int_S \{G_{ij}(\mathbf{x}, t; \xi, t_0)t_i(\mathbf{x}, t_0) - F_{ij}(\mathbf{x}, t; \xi, t_0)u_i(\mathbf{x}, t_0)\} dS(\mathbf{x}) dt_0 \quad (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 $u_i$ and $t_i$ remain constant over the time step $[t_0 + \Delta t, t_0]$.

### Laplace transform domain integral equation formulation

In the Laplace transform domain, the Green's function can be found in Cruse and Rizzo⁸ and Nowacki²⁴ and the singular integral equations are of the form

$$c_{ij}\bar{u}_i(\xi, s) = \int_S \{\bar{G}_{ji}(\mathbf{x}, \xi, s)\bar{t}_i(\mathbf{x}, s) - \bar{F}_{ji}(\mathbf{x}, \xi, s)\bar{u}_i(\mathbf{x}, s)\} dS(\mathbf{x}) \quad (15)$$

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 $\bar{G}_{ij}$ and $\bar{F}_{ij}$ 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¹⁰ 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 $\bar{u}_j$ and $\bar{t}_j$ 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

$$f(\mathbf{x}, t) = \frac{1}{2\pi i} \int_{\beta - i\infty}^{\beta + i\infty} \bar{f}(\mathbf{x}, s) e^{st} ds \quad i = \sqrt{-1} \quad (16)$$

where $\beta > 0$ is arbitrary but greater than the real part of all the singularities of $\bar{f}(\mathbf{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 parameter $s$, that is

$$s = \beta + i(2\pi/T)n \quad (17)$$

In equation (17), $T$ is the total time interval of interest and $n$ ranges from 1 to $N$. Good accuracy of the results cannot be obtained by numerical Laplace inversion algorithms based on real values of $s$ only, as shown in Reference 10. The numerical Laplace inversion algorithm used, based on a trapezoidal quadrature scheme for the integral of equation (16), is that of Durbin²⁵ which employs the fast Fourier transform concept of Cooley and Tukey²⁶ for expediting the calculations involved.

### Fourier transform domain integral equation formulation

As far as the solution in the Fourier transform domain is concerned, equation (15) can still be used under the assumptions made earlier if $s = iw$. The only difference is in the numerical inversion of the transform solution, since the inversion integral is of a form different than equation (16):

$$f(\mathbf{x}, t) = \int_{-\infty}^{\infty} \bar{f}(\mathbf{x}, w) e^{i2\pi w t} dw \quad (18)$$

For this case, the fast Fourier transform algorithm of Cooley and Tukey,²⁶ 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 $\Omega$ (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.²⁷ 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²⁸ 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.²⁹

## 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.

<figure><img src="image_4.png" /><figcaption>Figure 3. Circular cavity at different transform parameter (*s*) values</figcaption></figure>
In the three-dimensional approach of the time domain formulation, the information from Appendix I can be collected and used to rewrite the singular integral equations (11) as follows:

$$0.5\delta_{ij}u_i(\xi, t) = \sum_{K=1}^{N} \int_{S} \{U_{ij}(\mathbf{x}, K\Delta t; \xi)t_i(\mathbf{x}, (N-K+1)\Delta t) \\ - n_m[(T_{ijm}(\mathbf{x}, K\Delta t; \xi) + Q_{ijm}(\mathbf{x}, K\Delta t; \xi))u_i(\mathbf{x}, (N-K+1)\Delta t) \\ - Q_{ijm}(\mathbf{x}, K\Delta t; \xi)u_i(\mathbf{x}, (N-K)\Delta t)]\} dS(\mathbf{x}) \quad (19)$$

The time axis is now discretized into L segments so that the total time of interest is $T = L\Delta 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:

$$0.5\mathbf{u}_j^n = \sum_m \sum_i \mathbf{DG}_{ij}^{nm} \mathbf{t}_i^m - \sum_m \sum_i \mathbf{DF}_{ij}^{nm} \mathbf{u}_i^m \quad (20)$$

with $i, j = 1, 2, \dots, J$, $n = 1, 2, \dots, L$ and $m = 1, 2, \dots, n$. The matrices $\mathbf{DG}_{ij}^{nm}$ and $\mathbf{DF}_{ij}^{nm}$ are the spatial integrals of the kernels $U_{ij}$, $T_{ijm}$ and $Q_{ijm}$, respectively. In the latter case, use was made of the finite difference scheme, equation (38), which allowed the velocity to be expressed terms of the displacement.

For the two-dimensional approach of the time domain formulation, equation (14) can also be recast as equation (20), where $\mathbf{DG}_{ij}^{nm}$ and $\mathbf{DF}_{ij}^{nm}$ are now given by

$$\mathbf{DG}_{ij}^{nm} = \int_{(m-1)\Delta t}^{m\Delta t} \int_{\Delta S} G_{ij} \, dS \, dt_0 \\ \mathbf{DF}_{ij}^{nm} = \int_{(m-1)\Delta t}^{m\Delta t} \int_{\Delta S} F_{ij} \, dS \, dt_0 \quad (21)$$

and

with $\Delta 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.¹³ If a two-dimensional approach were to be used, 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.5\bar{\mathbf{u}}_j = \sum_i \mathbf{D}\bar{\mathbf{G}}_{ji}\bar{\mathbf{t}}_i - \sum_i \mathbf{D}\bar{\mathbf{F}}_{ji}\bar{\mathbf{u}}_i \quad (22)$$

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 $\ln r$ for the $\bar{\mathbf{G}}_{ij}$ tensor and of the form $r, i n_j / r$ for the $\bar{\mathbf{F}}_{ij}$ tensor. The behaviour of these singularities as $r \to 0$ is well documented.³⁰
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 \to 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'_{ij}$ and $F'_{ij}$. This property implies that at time step $n\Delta t$ only the matrices $DG_{ij}^{nm}$ and $DF_{ij}^{nm}$ need be computed, while the remaining matrices $DG_{ij}^{nm-1}$, $DG_{ij}^{nm-2}$, ..., $DG_{ij}^{n1}$ (similarly for $DF_{ij}^{nm}$) 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\Delta t/L \le 1 \cdot 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,¹⁵ Pao and Mow¹⁷) 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 $\sigma_x = S_0H(t-t_n)$, $\sigma_y = (\nu/(1-\nu))S_0H(t-t_n)$ and $\tau_{xy} = 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, $\nu$ is Poisson's ratio and $t_n$ 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:

$$ 
\begin{aligned}
\lambda &= \mu = 3,597,000 \text{ lb/in}^2 \\
\rho &= 0.00025 \text{ lb-sec}^2/\text{in}^4 \\
a &= 212 \text{ in.}
\end{aligned} 
\qquad (23) 
$$

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 \text{ in/sec}, \quad c_2 = 120,000 \text{ in/sec} \qquad (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 $2\tau$, is equal to $2a/c_1 = 0.002$ sec. The intensity of the applied load $S_0$ is conveniently taken as equal to $-1.0$ lb/in.²

In all cases, the boundary of the circular hole is discretized into 20 equal segments. Figures 4-6 depict the time history of the ratio $\sigma_{\theta\theta}/S_0$ at the boundary of the cavity for values of the
<figure><img src="image_5.png" /><figcaption>Figure 4. Circumferential stress time history at $\theta = 90^\circ$</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5. Circumferential stress time history at $\theta = 0^\circ$</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6. Circumferential stress time history at $\theta = 180^\circ$</figcaption></figure>

polar angle $\theta$ of 90°, 0° and 180°, respectively. The total time interval of interest $T$ is taken as equal to 0.010 sec or 10 half-transit times $\tau$. In the aforementioned figures the results obtained by the time domain and the two integral transform BEMs are plotted along with the analytic-numerical results of Baron and Matthews¹⁵ and Pao and Mow¹⁷. The results of the static equivalent of this dynamic problem are also plotted.

In the time stepping algorithm used in conjunction with the time domain BEM, the time step $\Delta t$ is taken as equal to 0.0003 sec and 23 time steps are employed in order to get a solution accurate for 0.00069 sec. A total of $N = 10$ values of the transform parameter $s$, $s = 6/T + i(2\pi/T)n$, $n = 1, 2, ..., N$, are considered in the Laplace transform BEM. Finally, $N = 16$ values of the transform parameter $w$ are used in the Fourier transform BEM, where $w = (\Omega/N)n$, $n = 1, 2, ..., N$, with $\Omega = 3100$ cycles/sec along with an additional 16 values of 'trailing' zeros added. It should be noted that both integral transform results are accurate for the first 80 per cent of the total time interval $T$.

At this point it should be mentioned that the time domain formulation employing the two-dimensional kernels $G_{ij}$ and $F_{ij}$ 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 $\Delta t = 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 = 5\tau$. The maximum value for the stress concentration $\sigma_{\theta\theta}/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 = 4\tau$ and $t = 6\tau$.

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 $\theta = 0^\circ$ act as a rigid wall and as a consequence the applied stress $\bar{\sigma}_y$ 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 $\theta = 0^\circ$. The reason the transformed BEM results (and Pao's as well) do not exhibit this sudden rise in stress at $t = 0^+$ 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 $\Delta t$. As a consequence, the sudden rise in $\sigma_{\theta\theta}$ does not materialize at $t = 0^+$. 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.

<table><caption>Table I. Computer time and memory requirements for the three BEM approaches</caption><thead><tr><th>BEM</th><th>CP execution time (sec)</th><th>Total CPU time (sec)</th><th>Field length during execution (octal words)</th><th>Cost ($)</th></tr></thead><tbody><tr><td>Time domain</td><td>953.2</td><td>962.2</td><td>74,200</td><td>51.9</td></tr><tr><td>Laplace domain</td><td>201.2</td><td>212.6</td><td>44,700</td><td>10.2</td></tr><tr><td>Fourier domain</td><td>358.9</td><td>369.3</td><td>46,100</td><td>17.5</td></tr></tbody></table>

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 $λ$ and $μ$ 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.

ACKNOWLEDGEMENTS

The author would like to thank a number of individuals for useful discussions and encouragement during the preparation of this work. These are Professors D. E. Beskos, P. K. Banerjee, R. P. Shaw and T. G. Davies. The receipt of a summer (1981) grant from the Research
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'_{ij}$ and $F'_{ij}$ can be found in Eringen and Suhubi.³ In particular,

$$G'_{ij}(\mathbf{x}', t; \xi') = \frac{1}{4\pi\rho} \left\{ \left( \frac{3r_i r_j}{r^3} - \frac{\delta_{ij}}{r} \right) \int_{c_1^{-1}}^{c_2^{-1}} \lambda \delta(t - \lambda r) \, d\lambda \right. \\ \left. + \frac{r_i r_j}{r^3} \left[ \frac{1}{c_1^2} \delta\left(t - \frac{r}{c_1}\right) - \frac{1}{c_2^2} \delta\left(t - \frac{r}{c_2}\right) \right] + \frac{\delta_{ij}}{rc_2^2} \delta\left(t - \frac{r}{c_2}\right) \right\} \quad (25)$$

and

$$F'_{ik} = H'_{ijk} n_j = \rho [(c_1^2 - 2c_2^2)G'_{mk,m}\delta_{ij} + c^2(G'_{ik,j} + G'_{jk,i})]n_j \quad (26)$$

where

$$H'_{ijk}(\mathbf{x}', t; \xi') = \frac{1}{4\pi} \left\{ -6c_2^2 \left[ 5 \frac{r_i r_j r_k}{r^5} - \frac{\delta_{ij} r_k + \delta_{ik} r_j + \delta_{jk} r_i}{r^3} \right] \cdot \int_{c_1^{-1}}^{c_2^{-1}} \lambda \delta(t - \lambda r) \, d\lambda \right. \\ + 2 \left[ 6 \frac{r_i r_j r_k}{r^5} - \frac{\delta_{ij} r_k + \delta_{ik} r_j + \delta_{jk} r_i}{r^3} \right] \left[ \delta\left(t - \frac{r}{c_2}\right) - \frac{c_2^2}{c_1^2} \delta\left(t - \frac{r}{c_1}\right) \right] \\ + 2 \frac{r_i r_j r_k}{r^4 c_2} \left[ \dot{\delta}\left(t - \frac{r}{c_2}\right) - \frac{c_2^3}{c_1^3} \dot{\delta}\left(t - \frac{r}{c_1}\right) \right] - \frac{r_k \delta_{ij}}{r^3} \left( 1 - 2 \frac{c_2^2}{c_1^2} \right) \left[ \delta\left(t - \frac{r}{c_1}\right) \right. \\ \left. + \frac{r}{c_1} \dot{\delta}\left(t - \frac{r}{c_1}\right) \right] - \frac{\delta_{ik} r_j + \delta_{jk} r_i}{r^3} \left[ \delta\left(t - \frac{r}{c_2}\right) + \frac{r}{c_2} \dot{\delta}\left(t - \frac{r}{c_2}\right) \right] \} \quad (27)$$

In the above expressions, $r_i = (x'_i - \xi'_i)$, $r = |\mathbf{x}' - \xi'|$ and the subscripts range from 1 to 3. The Green's function $G'_{ij}(\mathbf{x}', t; \xi', \tau)$ obeys the causality condition

$$G'_{ij}(\mathbf{x}', t; \xi', \tau) = 0 \quad \text{if } c_1(t - \tau) < r \quad (28)$$

and has the following time translation property:

$$G'_{ij}(\mathbf{x}', t + t_0; \xi', \tau + t_0) = G'_{ij}(\mathbf{x}', t; \xi', \tau) \quad (29)$$

The following property of the derivatives of the Dirac delta function will be used shortly:

$$\int_a^b \delta^{(j)}(e-d)f(e) \, de = (-1)^{j} f^{(j)}(d) \quad a \le d \le b \quad (30)$$

In addition, it should be noted that $r_i = r_i/r$.

The tensors $G'_{ij}$ and $F'_{ij}$ 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 $\xi' = \xi$ envelops a spherical region of radius $R = ct$, where $c$ is the appropriate propagation velocity. Notice that $\xi'$ 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) becomes

$$ \sum_{K=1}^{N} \int_{S} \int_{Z_{K-1}}^{Z_{K}} G_{ij}''(\mathbf{x}', K \Delta t; \xi) t_i(\mathbf{x}, (N-K+1)\Delta t) \, dz \, dS(\mathbf{x}) \quad (31) $$

where the kernel function $G_{ij}''$ is the same as $G_{ij}'$ 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\Delta t$ with $\Delta t$ being an appropriate time increment and $\mathbf{x}$ is the projection of $\mathbf{x}'$ on the XY plane. Furthermore, $Z_k$ is the vertical distance $x_3$ between $\mathbf{x}'$ and $\mathbf{x}$ at time $K\Delta t$, that is

$$ Z_K = \sqrt{[(cK \Delta t)^2 - (\mathbf{x} - \xi)^2]} \quad (32) $$

Therefore, the integration in (31) takes place over area patches whose horizontal ordinates are along $S$ and whose vertical ordinates are the distances $Z_K$ that mark the location of the spherical wave originating from $\xi$ at time increments $K\Delta t, K = 1, 2, \dots, N$.

Assuming that the tractions remain constant during the time interval $[K\Delta t, (K-1)\Delta t]$, equation (31) becomes

$$ \sum_{K=1}^{N} \int_{S} U_{ij}(\mathbf{x}, K \Delta t; \xi) t_i(\mathbf{x}, (N-K+1)\Delta t) \, dS(\mathbf{x}) \quad (33) $$

where

$$ U_{ij} = (2/4\pi\rho) \left[ \int_{c_2}^{c_1} \frac{dc}{c^3} \int_{Z_{K-1,c}}^{Z_{K,c}} \frac{1}{r} (3r_{,i}r_{,j} - \delta_{ij}) \, dz + \frac{1}{c_1^2} \int_{Z_{K-1,c_1}}^{Z_{K,c_1}} \frac{r_{,i}r_{,j}}{r} \, dz + \frac{1}{c_2^2} \int_{Z_{K-1,c_2}}^{Z_{K,c_2}} \frac{1}{r} (\delta_{ij} - r_{,i}r_{,j}) \, dz \right] \quad (34) $$

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 $\mathbf{x}'$ and $\xi$ and that indices $i$ and $j$ assume the values 1 and 2 only. The above equation can also be found in Niwa *et al.*¹⁴

Similarly, the second integral of equation (11) can be written as

$$ \sum_{K=1}^{N} \int_{S} n_{m} \{ T_{ijm}(\mathbf{x}, K \Delta t; \xi) u_{i}(\mathbf{x}, (N-K+1)\Delta t) + Q_{ijm}(\mathbf{x}, K \Delta t; \xi) \dot{u}_{i}(\mathbf{x}, (N-K+1)\Delta t) \} dS(\mathbf{x}) \quad (35) $$

where

$$ T_{ijm} = \frac{2}{4\pi} \left[ -6c_2^2 \int_{c_2}^{c_1} \frac{dc}{c^3} \int_{Z_{K-1,c}}^{Z_{K,c}} \frac{1}{r^2} \left\{ 5r_{,i}r_{,j}r_{,m} - (\delta_{im}r_{,j} + \delta_{ij}r_{,m} + \delta_{jm}r_{,i}) \right\} dz - \int_{Z_{K-1,c_1}}^{Z_{K,c_1}} \frac{1}{r^2} \left\{ 2 \left(\frac{c_2}{c_1}\right)^2 \left[ 6r_{,i}r_{,j}r_{,m} - (\delta_{im}r_{,j} + \delta_{ij}r_{,m} + \delta_{jm}r_{,i}) \right] + \left[ 1 - 2 \left(\frac{c_2}{c_1}\right)^2 \right] \delta_{im}r_{,j} \right\} dz + \int_{Z_{K-1,c_2}}^{Z_{K,c_2}} \frac{1}{r^2} \left\{ 12r_{,i}r_{,j}r_{,m} - (2\delta_{im}r_{,j} + \delta_{ij}r_{,m} + \delta_{jm}r_{,i}) \right\} dz \right] \quad (36) $$
and

$$Q_{ijm} = \frac{2}{4\pi} \left[ -\frac{1}{c_1} \int_{Z_{K-1,c_1}}^{Z_{K,c_1}} \frac{1}{r} \left\{ 2 \left( \frac{c_2}{c_1} \right)^2 r_{,i} r_{,j} r_{,m} + \left[ 1 - 2 \left( \frac{c_2}{c_1} \right)^2 \right] \cdot \delta_{im} r_{,j} \right\} dz + \frac{1}{c_2} \int_{Z_{K-1,c_2}}^{Z_{K,c_2}} \frac{1}{r} \left\{ 2r_{,i} r_{,j} r_{,m} - (\delta_{ij} r_{,m} + \delta_{jm} r_{,i}) \right\} dz \right] \quad (37)$$

Furthermore, the velocity $\dot{u}_i$ in equation (35) can be replaced by a backward finite difference in time:

$$\dot{u}_i(\mathbf{x}, (N-K+1)\Delta t) \approx (u_i(\mathbf{x}, (N-K+1)\Delta t) - u_i(\mathbf{x}, (N-K)\Delta t))/\Delta t \quad (38)$$

Finally, all that remains to be done is the spatial integration of $U_{ij}$, $T_{ijm}$ and $Q_{ijm}$ over the bounding surface S of a representative two-dimensional slice of the cylindrical body. Then, the discrete form of equation (11) can be used at each time step to solve for the unknown boundary quantities in terms of the known ones.

## APPENDIX II

If the general three-dimensional kernel $G'_{ij}$ is spatially integrated over the $x_3$ co-ordinate, then the two-dimensional form of this kernel is obtained³ as follows:

$$G_{ij} = \int_{-\infty}^{\infty} G'_{ij}(\mathbf{x}, t; \xi') \, dx_3 = \frac{1}{2\pi\rho} \left\{ \begin{aligned} & \frac{[2t^2 - r^2/c_1^2]H(t-r/c_1)}{[t^2 - r^2/c_1^2]^{1/2}} \\ & - \frac{[2t^2 - r^2/c_2^2]H(t-r/c_2)}{[t^2 - r^2/c_2^2]^{1/2}} \left( \frac{r_i r_j}{r^4} \right) \\ & - [H(t-r/c_1)(t^2-r^2/c_1^2)^{1/2} - H(t-r/c_2)(t^2-r^2/c_2^2)^{1/2}] \\ & \cdot \frac{\delta_{ij}}{r^2} + \frac{H(t-r/c_2^2)}{c_2^2[t^2-r^2/c_2^2]^{1/2}} \delta_{ij} \end{aligned} \right\} \quad (39)$$

The square brackets in equation (39) denote retarded time values, i.e. if the argument inside the brackets becomes non-positive the whole expression is set equal to zero and $r$ is the distance between $\xi$ and $\mathbf{x}$.

An expression for the two-dimensional form $F_{ij}$ of $F'_{ij}$ can easily be derived in view of equation (26).

## APPENDIX III

The kernel $\bar{G}_{ij}$ appearing in equation (15) has the following form:

$$\bar{G}_{ij} = (1/2\pi\mu) \{\bar{\Psi}\delta_{ij} - \bar{X} r_{,i} r_{,j}\} \quad (40)$$

with

$$\bar{\Psi} = K_0\left(\frac{sr}{c_2}\right) + \frac{c_2}{sr}\left[K_1\left(\frac{sr}{c_2}\right) - \frac{c_2}{c_1}K_1\left(\frac{sr}{c_1}\right)\right] \quad (41)$$

and

$$\bar{X} = K_2 \left( \frac{sr}{c_2} \right) - \frac{c_2^2}{c_1^2} K_2 \left( \frac{sr}{c_1} \right) \quad (42)$$
In the above expressions, *r* is the radial distance between the source point *ξ* (the origin) and the receiver point *x*, while *K*₀, *K*₁ and *K*₂ are the modified Bessel functions of the second kind and order zero, one and two, respectively.

The kernel $\bar{F}_{ij}$, derived from $\bar{G}_{ij}$ through the use of the constitutive relations, is as follows:

$$ \bar{F}_{ij} = \rho [(c_1^2 - 2c_2^2)\bar{G}_{im,m}n_j + c_2^2(\bar{G}_{ij,m} + \bar{G}_{im,j})n_m] \quad (43) $$

The resulting expression for $\bar{F}_{ij}$ may be condensed to

$$ \bar{F}_{ij} = (1/2\pi) \left\{ \left( \frac{d\bar{\Psi}}{dr} - \frac{1}{r}\bar{X} \right) \left( \delta_{ij} \frac{\partial r}{\partial n} + r_i n_j \right) - \frac{2}{r}\bar{X} \left( n_j r_i - 2r_i r_j \frac{\partial r}{\partial n} \right) \right. \\ \left. - 2 \frac{d\bar{X}}{dr} r_i r_j \frac{\partial r}{\partial n} + \left( \frac{c_2^2}{c_1^2} - 2 \right) \left( \frac{d\bar{\Psi}}{dr} - \frac{d\bar{X}}{dr} - \frac{1}{r}\bar{X} \right) r_i n_j \right\} \quad (44) $$

where use is made of the following relation:

$$ \frac{\partial r}{\partial n} = \frac{\partial r}{\partial x} \cdot \frac{\partial x}{\partial n} = r_1 n_1 + r_2 n_2 \quad (45) $$

## REFERENCES

1. B. B. Baker and E. T. Copson, *The Mathematical Theory of Huygens' Principle*, Oxford Univ. Press, London, 1939.

2. P. M. Morse and H. Feshbach, *Methods of Theoretical Physics*, McGraw-Hill, New York, 1953.

3. A. C. Eringen and E. S. Suhubi, *Elastodynamics*, Academic Press, New York, 1975.

4. V. D. Kupradze, Editor, *Three-Dimensional Problems of the Mathematical Theory of Elasticity and Thermoelasticity*, North-Holland, Amsterdam, 1979.

5. M. B. Friedman and R. P. Shaw, 'Diffraction of pulses by cylindrical obstacles of arbitrary cross-section', *J. Appl. Mech.*, **29**, 40-46 (1962).

6. L. H. Chen and J. Schweikert, 'Sound radiation from an arbitrary body', *J. Acoust. Soc. Amer.*, **35**, 1626-1632 (1963).

7. R. P. Banaugh and W. Goldsmith, 'Diffraction of steady elastic waves by surfaces of arbitrary shape', *J. Appl. Mech.*, **30**, 589-597 (1963).

8. T. A. Cruse and F. J. Rizzo, 'A direct formulation and numerical solution of the general transient elastodynamic problem, I', *J. Math. Anal. Appl.*, **22**, 244-259 (1968).

9. T. A. Cruse, 'A direct formulation and numerical solution of the general transient elastodynamic problem, II', *J. Math. Anal. Appl.*, **22**, 341-355 (1968).

10. G. D. Manolis and D. E. Beskos, 'Dynamic stress concentration studies by boundary integrals and Laplace transform', *Int. J. num. Meth. Engng*, **17**, 573-599 (1981).

11. Y. Niwa, S. Kobayashi and N. Azuma, 'An analysis of transient stresses produced around cavities of arbitrary shape during the passage of travelling waves', *Memoirs Fac. Engng. Kyoto Univ.*, **37**, pt 2, 28-46 (1975).

12. Y. Niwa, S. Kobayashi and T. Fukui, 'Applications of integral equation methods to some geomechanical problems', in *Numerical Methods in Geomechanics* (Ed. C. S. Desai), A.S.C.E., New York, 1976, pp.120-131.

13. D. M. Cole, D. D. Kosloff and J. B. Minster, 'A numerical boundary integral equation method for elastodynamics, I', *Bul. Seism. Soc. Am.*, **68** (5), 1331-1357 (1978).

14. Y. Niwa, T. Fukui, S. Kato and K. Fujiki, 'An application of the integral equation method to two-dimensional elastodynamics', *Theor. Appl. Mech.*, **28**, Univ. of Tokyo Press, 281-290 (1980).

15. M. L. Baron and A. T. Matthews, 'Diffraction of a pressure wave by a cylindrical cavity in an elastic medium', *J. Appl. Mech.*, **28**, 347-354 (1961).

16. M. L. Baron and R. Parnes, 'Displacements and velocities produced by the diffraction of a pressure wave by a cylindrical cavity in an elastic medium', *J. Appl. Mech.*, **29**, 385-395 (1962).

17. Y. H. Pao and C. C. Mow, *Diffraction of Elastic Waves and Dynamic Stress Concentrations*, Crane, Russak, New York, 1972.

18. J. Miklowitz, 'Scattering of a plane elastic compressional pulse by a cylindrical cavity', *Proc. 11th Congr. Appl. Mech.*, Munich, 469-483 (1964).

19. S. K. Datta, 'Scattering of elastic waves', in *Mechanics Today Vol 4* (Ed. S. Nemat-Nasser), Pergamon Press, Oxford, 1978, pp.149-208.

20. P. K. Banerjee and R. Butterfield, *Boundary Element Methods in Engineering Science*, McGraw-Hill, New York, 1981.
21. A. Y. DeHoop, 'Representation theorems for the displacement in an elastic solid and their application to elasto-dynamic diffraction theory', *Doctoral dissert.*, Tech. Hogeschole, Delft (1958).

22. L. T. Wheeler and E. Sternberg, 'Some theorems in classical elastodynamics', *Arch. Rat. Mech. Anal.*, **31**, 51–90 (1968).

23. E. Sternberg and R. A. Eubanks, 'On the concept of concentrated loads and an extension of the uniqueness theorem in the linear theory of elasticity', *J. Rat. Mech. Anal.*, **4**, 135–145 (1955).

24. W. Nowacki, *Dynamics of Elastic Systems*, Wiley, New York, 1963.

25. F. Durbin, 'Numerical inversion of Laplace transforms: an efficient improvement to Dubner and Abate's method', *Computer J.*, **17**, 371–376 (1974).

26. J. W. Cooley and J. W. Tukey, 'An algorithm for the machine calculation of complex Fourier series', *Math. Comp.*, **19**, 297–301 (1965).

27. H. A. Schenk, 'Improved integral formulation to acoustic radiation problems', *J. Acoust. Soc. Am.*, **44**(1), 41–58 (1968).

28. G. R. C. Tai and R. P. Shaw, 'Helmholtz equation eigenvalues and eigenmodes for arbitrary domains', *J. Acoust. Soc. Am.*, **56**(3), 796–804 (1974).

29. S. Kobayashi and N. Nishimura, 'Transient stress analysis of tunnels and caverns of arbitrary shape due to travelling waves', in *Developments in Boundary Element Methods*, Vol. 2 (Eds P. K. Banerjee and R. P. Shaw), Applied Science, London, 1981.

30. S. G. Mikhlin, *Multidimensional Singular Integrals and Integral Equations*, Pergamon Press, London, 1965.