INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, VOL.

26, 1709-1728 (1988)

TIME-DOMAIN TRANSIENT ELASTODYNAMIC ANALYSIS

OF 3-D SOLIDS BY BEM

S. AHMAD' AND P. K. BANERJEE'

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

SUMMARY
The numerical implementation of the Direct Boundary Element formulation for time-domain transient
analysis of three-dimensional solids is presented in a most general and complete manner. The present
formulation employs the space and time dependent fundamental solution (Stokes' solution) and Graffis
dynamic reciprocal theorem to derive the boundary integral equations in the time domain. A time-stepping
scheme is then used to solve the boundary initial value problem by marching forward in time. Higher order
shape functions are used to approximate the field quantities in space as well as in time, and a combination of
analytical (time-integration) and numerical (spatial-integration) integration is carried out to form a system
of linear equations. At the end of each time step, these equations are solved to obtain the unknown field
quantities at that time.
Finally, the accuracy and reliability of this algorithm is demonstrated by solving a number of example
problems and comparing the results against the available analytical and numerical solution.

INTRODUCTION
The ability to predict the dynamic response of solid bodies subjected to time and space dependent
loads and boundary conditions has gained considerable importance in all engineering fields such
as structural dynamics, machine foundation design, aircraft structure design, seismology and
soil-structure interaction analysis. For a realistic problem from any of the above engineering
fields, it is extremely difficult to obtain a closed-form solution. Therefore, numerical methods such
as the Finite Difference Method (FDM), the Finite Element Method (FEM) or the Boundary
Element Method (BEM) need to be used. Of the numerical methods, the most versatile and widely
used method is the FEM. It has, however, two major deficiencies:(a) it cannot model an infinite or
semi-infinite medium properly; and (b) the computational cost (both CPU and data preparation)
involved in analysing three-dimensional transient dynamic problems by the FEM is so enormous
that only a few researchers can afford it. Similarly, the FDM has not been used frequently,
primarily because of the difficulties associated with handling complicated geometries and
boundary conditions.
In contrast, it is convincingly demonstrated that accurate and efficient solutions to dynamic
problems can be easily obtained by using the BEM'-3 because an infinite or semi-infinite
medium can be modelled easily, and for linear problems only the surface of the geometry needs to
be discretized. Thus, a considerable reduction in the size of the problem can be achieved for a very
wide class of engineering problems.

* Assistant Professor
+ Professor

0029-598 1/88/08 1709-20s 10.00 Received 11 September 1987

0 1988 by John Wiley & Sons, Ltd. Revised 23 November 1987
 10970207, 1988, 8, Downloaded from [Link] by Ben Gurion University, Wiley Online Library on [04/07/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
1710 S. AHMAD AND P. K. BANERJEE

The boundary element formulation for the time-domain transient dynamic analysis with
constant boundary elements has been implemented by Cole et aL4 for the anti-plane strain case
(i.e. scalar problem), by Niwa et al.' for the two-dimensional wave scattering problem, by Rice
and Sadd6 for the anti-plane strain wave scattering problem, by Spyrakos and Beskos' and
Antes* for the strip-footing problem and by Mansur and Brebbiag-l o for the two-dimensional
problem.
Three-dimensional problems of transient dynamic by time-domain BEM were not attempted
until recently, principally because of the formidable task of numerical implementation and
enormous computing requirements. In order to reduce the computation and complications
involved, simplifications of the BEM formulation dictated by the nature of the problem to be
'
solved have been developed by Karabalis and Beskos.' They have simplified the BEM formul-
ation for the special case of transient surface loading on rectangular foundations. A more general
numerical implementation of the BEM formulation for time-domain transient dynamic analysis
has been attempted by Banerjee et aL2 and Banerjee and Ahmad' using constant temporal and
quadratic spatial variations of the field variables.
Most of the above mentioned work other than that of References 1 and 2 suffers from one or
more of the following: lack of generality, crude assumption of constant variation of the field
variables in space and time, inadequate treatment of singular integrals and unacceptable level of
accuracy. For example, Cole et aL4 found the transient dynamic formulation to be unstable,
leading to a building up of errors as the time-stepping progresses: Rice and Sadd6 found that
dominant errors in the method arise from integrating the Green's function over the singularity
and their time-domain formulation when applied to time harmonic problems reveals a large build
up of errors; and Niwa et aL5 find considerable errors in their numerical results and suggest that
use of higher order interpolation functions for approximating the temporal and spatial variations
of field variables may improve the accuracy and the stability of their method. Mansur and
BrebbiagV1Ore-implemented the scalar dynamic formulation of Cole et al.," as well as the
elastodynamic formulation outlined in Niwa et al.,' both for two-dimensional problems of
dynamics. They did not appear to have encountered any serious errors in their solutions.
It is not really surprising that a number of the above mentioned authors report some problems
of accuracy. Unlike the FEM method which has an energy based error minimization principle
embodied in the discretized system, the BEM system in the discretized state does not satisfy any
energy principle. Only the original boundary integral equation satisfies the energy principle
(namely, the reciprocal work theorem) exactly. Such a system therefore needs to be implemented
as accurately as possible.
In References 1 and 2 three-dimensional elastodynamic formulations were implemented using
quadratic isoparametric shape functions and a very sophisticated error control in integration. For
the transient dynamic analyses, however, only a constant temporal variation was used. While
these analyses eliminated most of the inaccuracies in the boundary solution, the interior stresses
were found to be somewhat inaccurate. This is primarily due to the presence of a boundary
acceleration term in the stress kernel which cannot be satisfactorily handled when the boundary
displacements are assumed to be constant over a time step.
Almost all of the above mentioned problems have been eliminated in the present work by using
a quadratic shape function in space and a linear shape function in time, taking care of the singular
integral in an accurate and elegant manner, using superior and sophisticated integration
techniques and implementing the BEM formulation in a complete and general manner. The time-
domain transient algorithm developed in this work is capable of producing accurate results for
general three-dimensional elastodynamic problems. Furthermore, the use of a linear temporal
 10970207, 1988, 8, Downloaded from [Link] by Ben Gurion University, Wiley Online Library on [04/07/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
TIME-DOMAIN TRANSIENT ELASTODYNAMIC ANALYSIS 1711

shape function makes it possible to extend the present algorithm to non-linear dynamic analysis,
which will be the subject of a forthcoming paper by the authors.

TRANSIENT BOUNDARY INTEGRAL FORMULATION

The direct boundary integral formulation for a general, transient, elastodynamic problem can be
constructed by combining the fundamental point-force solution (Stoke’s solution) of the govern-
ing equation of motion (Navier-Cauchy equations) with Graffis12 dynamic reciprocal theorem.
Details of this construction can be found in Banerjee et a1.2 For zero initial conditions and zero
body forces, the boundary integral formulation for transient elastodynamics reduces to
r

where
Gij*ti= [:Gij(x, r 5, z)t,(x, T)dr
Fij*ui=j:Fij(x. T; 5, z)ui(x,r)dr

are Riemann convolution integrals, and 5 and x are the space positions of the receiver (field point)
and the source (source point). The fundamental solutions G , and Fij are the displacements ui and
tractions ti at a point x at time Tdue to a unit force vector applied at a point 5 at a preceding time
T. These functions are listed in a compact and elegant form in Appendix I. The tensor C i j is the
discontinuity (or jump) tensor, arising from the singularity of the Fij kernel.
Equation (1) represents an exact formulation involving integration over the surface as well as
the time history. It should also be noted that this is,an implicit time-domain formulation because
the response at time T is calculated by taking into account the history of surface tractions and
displacements up to and including the time T. Furthermore, equation (1) is valid for both regular
and unbounded domains.
Once the boundary solution is obtained, the stresses at the boundary can be calculated by
combining the constitutive equations, the directional derivatives of the displacement vector and
the values of field variables in an accurate matrix formulation. For calculating displacements at
interior points equation (1) can be used with cij=6,, and the interior stresses can be obtained from

Ojk(5, TI=
I
CG;k(x, 5, T ) * t i ( x , T)-FGk(x, 5, T ) * u i ( x , T ) I dS(x)
The functions GZkand FGk of the above equation are listed in Appendix 11. These functions are
(3)

presented for the first time; they do not exist in any published literature. It can be seen that these
functions contain first and second order derivatives with respect to time. Therefore, use of
constant temporal variation of the field variables, as done in all previous work, produces
unsatisfactory results for interior stress.

NUMERICAL IMPLEMENTATION
Details of the numerical schemes for surface discretization, spatial integration, assembly and
solution of system equations can be found in References 1, 2 and 3. In this paper, only the
temporal integration and the time-marching scheme used in the present work are presented in
 10970207, 1988, 8, Downloaded from [Link] by Ben Gurion University, Wiley Online Library on [04/07/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
1712 S. AHMAD A N D P. K. BANERJEE

detail. However, for completeness a brief description of the rest of the numerical schemes is also
provided.
The surface modelling is done using either six noded triangles or eight or nine noded
quadrilaterals. The functional variation over them can be either linear or quadratic or a mixture
of linear and quadratic. The spatial integration is carried out to a pre-selected 3 to 5 digit
precision. Facilities to include all possible local and global boundary conditions have been
provided. The entire analysis can be carried out for an assemblage of up to 15 substructured
regions. At the interfaces between sub-regions, facilities to include sliding, bonded or spring
loaded connections have been included.
This implementation is a part of BEST3D (Boundary Element Solution Technique, 3-
Dimensional) system.

Time-marching scheme
In order to obtain the transient response at a time TN, the time axis is discretized into N equal
time intervals, i.e.
N
TN= 1 nAT
n= 1
(4)

where AT is the time step.

Utilizing equations (4) and (2), equation (1) can be written as

cijui(k,TN)- jTN
TN- I '::j
j S [ G i j t t - F i j u i ]dSdz= js[Gijti-Fijui] dSdz (5)

were the integral on the right hand side is the contribution due to past dynamic history.
It is of interest that equation (5), like equation (l), still remains an exact formulation of the
problem since no approximation has yet been introduced. However, in order to solve equation (5),
one has to approximate the time variation of the field quantities in addition to the usual
approximation of spatial variation. For this purpose a linear interpolation function is used which
is described with the resulting time-stepping algorithms as follows.
The field variables (i.e. displacements, tractions or stresses) are assumed to vary linearly during
a time step, i.e.

U i ( X , 7)= c [M;uy-'(x)+M",y(x)]
N

n= 1
(64
N
ti(X, T)= 1 [M;t;-'(x)+M"Fy(x)]
n= 1

where
N is the total number of time steps;
ul(x) and tl(x) represent the spatial variation of ui and t i , respectively, at time T,;
M,and MF are the temporal interpolation functions related to local time nodes I and F (see
Figure l), and are of the form
 10970207, 1988, 8, Downloaded from [Link] by Ben Gurion University, Wiley Online Library on [04/07/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
TIME-DOMAIN TRANSIENT ELASTODYNAMIC ANALYSIS 1713

Function
f(t)

Node 'F"

Tn-l T" time (1)

Figure 1. Definition of temporal interpolation nodes

where
4,(7)= 1 for (n- l ) A T < t S n A T , and
=O otherwise; and
or
4 " ( ~=)[H{ 7 - (n - l)A T }- H { T - nA T } ]
H being the Heaviside function.
For illustration purposes, consider the boundary integral equation ( 5 ) for the first time step, i.e.

CijUi(t9 TI)- j l l [ s [ G i j t i - F i j ~ i dSdr=O

] (8)

The time integration in equation (8) after using equation (6) is done analytically and the spatial
integration is performed numerically (see References 1, 2 and 3). This analytical temporal
integration will be discussed in the following section.
After the integrations and the usual assembly process, the resulting system equation has the
form
c41 {X')- CBkH yl>+ { 4 1{XO)-CB:l{ YO}= (0) (9)
where
[A] and [B] are the matrices related to the unknown and known field quantities, respectively;
{ X} and { Y } are the vectors of unknown and known field quantities, respectively;
for {X}and { Y } the superscript refers to the time step;
for [A] and [B] the superscript denotes the time step at which they are calculated, and the
subscript denotes the local time nodes (I or F) during that time-stepping interval.
 10970207, 1988, 8, Downloaded from [Link] by Ben Gurion University, Wiley Online Library on [04/07/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
1714 [Link] AND P. K. BANERJEE

Since all the unknowns at time T=O are assumed to be zero, equation (9) reduces to
C4l { X 1 }= CBkI { y '1 + D:l{ y O) (10)
Now consider the boundary integral equation for the second time step, i.e.

C i j U i ( S 9 T2)- 11;Is[Gijti-F, ui] dSdz = [I: ] s [ G i j 4 - F i j ~ i dSdt

] (11)

If the time interval ( T, - Tl) is the same as ( T , - To)the resulting coefficient matrices of the left
hand sides of equations (8) and (11) become identical. This is due to the time translation
properties of the fundamental solutions G , and F,, which contain time functions with arguments
(T-z), and therefore the convoluted integral corresponding to the interval T1<z < Tt with T
= T2is identical to that of the interval To< t < TI with T = Tl.
The right hand side of equation (11) is evaluated at time T = T2with the time integration over
the interval To to T, and thus provides the effects of the dynamic history of the first time interval
on the current time node (i.e. T2).
Now, the resulting system equation for this time node ( T2)is of the form
[&I {X'} -[BbI{ Y ' } + CAtI {XI} -CB:I { Y1}= - { C A : I { X ' } - C G I {Yl}
+ [A:] { X O }- [GI { YO)} (12)
Equation (12) can be rearranged such that
+
[Ak] {X'} = [ B k ] { Y '} - [A: +A:] {XI} [B: +B:] { Y '} + [B:] { Y O}
(13)
In the above equation, all the quantities on the right hand side are known. Therefore the
unknown vector { X 2 } at time Tz can be obtained by solving the above equation.
Thus, for the Nth time step, the boundary integral equation (5) can be written in a discretized
form as
N T 1

or
CAkl{xN)=CB;l{ Y N > + { R N 1 (15)
where vector { R"} is the effect of the past dynamic history on the current time node.
The above equation can be solved to find the unknown vector { X " ) at time [Link] may appear
at first glance that a prodigious calculation of coefficients is involved. However, a closer
examination will reveal that:
(i) if the time step size is constant, the [A;] and [B;]matrices do not change from time step to
time step;
(ii) for each time step, a new { R N }vector needs to be formed. This involves the evaluation of a
new set of coefficient matrices [A;], [B;], [A;] and [B;] involving the effects of the
dynamic history of the first time interval on the current time node. Eventually, however,
this contribution to (RN}reduces to zero and from that point onwards no new coefficients
need to be evaluated.
Finally, it is of interest to note that if the linear temporal shape functions M ; and M; are
replaced by Mf = MF = 0.5 $"(z), the present time-stepping scheme reduces to a constant tem-
 10970207, 1988, 8, Downloaded from [Link] by Ben Gurion University, Wiley Online Library on [04/07/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
TIME-DOMAIN TRANSIENT ELASTODYNAMIC ANALYSIS 1715

poral variation scheme. In this constant temporal variation scheme, the averaged value of a field
variable between the two local time nodes (I and F) is taken as the representative value for the
field variable during a time step. This averaging yields more accurate results'. than that obtained '
without averaging. Nevertheless, the latter approach has been invariably used by most of the
researchers in past.

Analytical temporal integration

For linear time interpolation, the field variables are expressed as

f = c W!f
N

n= 1
I- '(x) + M)f"(x)l (16)

where fn(x) represents the spatial variation of the field variablef(x, t)at time Tn(= nAT), and M,
and M, are the temporal shape functions expressed by equation (7).
The transient dynamic kernels listed in Appendices I and I1 have one or more of the following
time functions embedded in them.

Time function (I): S( T-t-r/c)

Time function (2): Ah( T-r-r/c)dt

Time function (3): 6( T-t-r/c) (19)

Time function (4): 8( T - t - r l c ) (20)

where 'c' is either the pressure wave velocity c1 or the shear wave velocity c2, and 6 is the Dirac
delta function.
Using equation (16), the time integrals related to any of the above listed time functions
(equations (17), (18X (19) and (20)) can be expressed as

I:g( T- z - r / c ) f ( x ,t)d t =

=
1
n=l r (n-1AT

[fn-'(x)!
g( T - T - ~ / C ) ~ ( X
nAT
t)dz
,

M!(t)g( T-.r-r/c)dT
n= 1 ( n - 1)AT

+ f "(x) yT
(n - 1)AT
Mt(r)g( T- t - r/c)d t
1
An important characteristic of the transient dynamic kernels is the time translation pro-
perty.' -3Because of this characteristic,at each time step, only the effects of the dynamic history of
the first time interval on the current time node need to be evaluated: i.e. at each time step the
analytical time integration has to be done only for n = 1. Thus, equation (21) reduces to
AT
g( T-r-r/c)f(x,t)dt=f'(x) 41(t)g( T-t-r/c)dr

+f '(x) 10
AT(t/AT)&(T)g(T- T -r/c)dr
 10970207, 1988, 8, Downloaded from [Link] by Ben Gurion University, Wiley Online Library on [04/07/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
1716 S. AHMAD AND P. [Link]

The time integrals on the right hand side of equation (22)can be rearranged into the following two
types of integrals.
Type A:

Type B: loAT r4 ,(r)g( T - t-r/c)dz

The convoluted time integrals of types A and B with time functions ( 1 ) and (2) are evaluated
analytically as follows.

Time function 1:

Type A J +,(t)s( T-r-r/c)dt=4,( T-r/c)=H( T-r/c)-H( T-AT-r/c)

0

Type B: J:T tq5,(t)6 ( T - t - r/c)=( T - r/c)& ( T - r/c)

where
[(A2/2)H(T-lr)]$t =O if T< r/cl

- 1 1
if T>r/c2
2c: 2c:
The second term on the right hand side of equation (27)can be obtained in a similar manner by
replacing ‘T’ by ‘T-AT’ in equation (28).

Type B: IOAT
1”” l/Cl
AS( T-r-lr)t4,(t)dldt=
[Link] A( T - l r ) 4 , ( T - b ) d l

= [:I:TA+l( T-Ar)dl-

= TC(A2/2)4,(7.-vr)l;::
i:::: rA24,(T-Ar)dA

- [(n3/3)4 ( T - wi :;:: (29)

The terms on the right hand side of equation (29)are evaluated in a similar manner as that of
equation (27).
For the time functions (3) and (4), the time integrals of equations (21)are evaluated as follows.
 10970207, 1988, 8, Downloaded from [Link] by Ben Gurion University, Wiley Online Library on [04/07/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
TIME-DOMAIN TRANSIENT ELASTODYNAMIC ANALYSIS 1717

Time function 3

Timefunction 4. The temporal integration involving the time function (4)(i.e. 8( T - z - r l c ) ) is

approximated by using a backward finite difference scheme as follows:

f' (a) Mesh I 1

3 six-noded triangular
elements per octan

i' (b) Mesh 12

3 nine-noded quadrilateral
elements per octan

i' ( c ) Mesh 13
3 eight-noded quadrilateral
elements per octan

x
2

Figure 2. Boundary element models for radial expansion of a cavity

 10970207, 1988, 8, Downloaded from [Link] by Ben Gurion University, Wiley Online Library on [04/07/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
1718 S. AHMAD A N D P. K. BANERJEE

EXAMPLES O F APPLICATION

(a) Spherical cavity

A spherical cavity is embedded in an infinitely extending medium with E = 8-993x lo6, v =025
and p = 2 3 x in self-consistent units. The radius of the cavity is a=212. Three different
meshes for its surface discretization are shown in Figure 2. Using the built in symmetry
capabilities, this problem is modelled by one octant only. The characteristic times required for the

'r/(u, ) AM
&St
:Ilpo
1.4

.e

:1 *

Loadin; Curve
--

a
o
Mesh #I

Mesh 13
Mesh 1 2

.e t
0

.8 -LA,
.0 1.0 2.0 3.0 4.8 5.0 6.0 7.0 .0
-t tc
- 1

( a ) Time History of Radial Displacement at Cavity Wall

.9
I IP(t)

%/Po

/!il ---
a Mesh 13

o Mesh 12

.8

-.I .0 1.0 2.8 3.0 4.0 5.8 6.0 7.0

1.0

a
(b) Time Histocy of Tangential Stress at Cavity H a l l
Figure 3. (a) Time history of radial displacement at cavity wall. (b) Time history of tangential stress at cavity Wall
 10970207, 1988, 8, Downloaded from [Link] by Ben Gurion University, Wiley Online Library on [04/07/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
TIME-DOMAIN TRANSIENT ELASTODYNAMIC ANALYSIS 1719

pressure and shear waves to travel a cavity radius are 0*00102and 0.00177 secs, respectively. Four
cases are considered.

(i) Spherical cavity under sudden radial expansion. A radial pressure P = lo00 is suddenly
applied and maintained at the cavity surface. Figure 3(a) shows the time variation of radial
displacement at the cavity surface obtained by the time-domain algorithm. Concurrently plotted
are the analytical results of Tirno~henko.'~ The dynamic displacements are normalized by the
analytical static solution. Figure 3(b) shows the normalized hoop stress at the surface of the cavity

OGIO.

LOAOING CURVE
0008.

c
z
w
z
0 0006
<
0 TIME
Y
0
2
<
0
a
0-004

0002
.
- NUMERICAL (A*=-)
NUMERICAL IAi-00006)

OQO

-9002 TIME

Figure 4. Radial expansion of a cavity under triangular pulse

LOADING CURVE

TIME

- NUMERICAL t&r00003 SEC)

NUMERY:AL~ASDWOO~SEC)

TIME

Figure 5. Radial expansion of a cavity under rectangular pulse

 10970207, 1988, 8, Downloaded from [Link] by Ben Gurion University, Wiley Online Library on [04/07/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
1720 S. AHMAD A N D P. K. BANERJEE

compared with the analytical solution. The numerical results of displacements for all three meshes
agree well with the analytical solution. The stresses, however, show some oscillations at later
times. This is principally due to the abruptness of the applied loading, which results in a shock
that cannot be modelled in a numerical scheme. The time step used in this analysis is equal to
0.33t*, where t* is the time required for the pressure wave to travel the cavity radius.

(ii) Spherical cavity subjected to a triangular pulse of radial pressure. A triangular pulse of
radial pressure, as shown in Figure 4, is applied at the cavity surface. This example is solved by
using the present analysis with two different time steps. The radial displacements at the cavity
1.4

1.2 - / *gk6 xe e
/ 1 f,a
* a
fir,
1.n
- I @ fi-8. 8-g 8-1 9-1 -g s-8
.e - I - Analytical [Link] Step
Loadi
A Mesh 8,
Linear Var.
r'
(ur
St
)Anal.
.6 -
'
I 4 P(t)
0
X Mesh 11, Constant Var.
Mesh 13, Linear Var.

I .t
P(t)
X
1.n - X

.e - I\ X

-
2t*
I x'Q a x 3
-s
-I ,
\ @

* d_ItN/
a x
-9 5-B 4-w- m-8 H @-

.4
-I X B S Analytical [Link] Step

I a E%,Linear Var.
Mesh 11,constant Var.
.P x Mesh 13,Linear Var.
O Mesh 12,Linear Var.

Figure 6. (a) Time history of radial displacement at a distance R=0.2a from the cavity wall. (b) Time history of radial
stress at R=0.2a
 10970207, 1988, 8, Downloaded from [Link] by Ben Gurion University, Wiley Online Library on [04/07/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
TIME-DOMAINTRANSIENT ELASTODYNAMIC ANALYSIS 1721

surface are plotted in Figure 4. The numerical results from both the time steps are almost
identical. Thus, this example demonstrates the stability of the present algorithm.

(iii) Spherical cavity subjected to a rectangular pulse of radial pressure. A rectangular pulse of
radial pressure, as shown in Figure 5, is applied at the cavity surface. This example is also solved
by using two different time-increments. Figure 5 shows the time history of the radial displacement
of the cavity. By comparing these results with those due to the triangular pulse (i.e. Figure 4), it
can be seen that, in general, the displacement at any time interval due to the rectangular pulse is
twice of that due to the triangular pulse. This is because the response depends upon the total
impulse, and the total impulse due to the rectangular pulses is double that due to the triangular
one. In addition, since the duration of the energy input is the same in both problems, the response
curves for both cases have the same shape.

(iu) Spherical caoity subjected to ramp loading. It was found from the analysis of previous cases
that a suddenly applied loading results in a shock which could not be modelled very effectively.
This is so because the fundamental solution due to Stoke which has been used in the development
of the BEM formulation does not satisfy the kinetic energy requirement at the shock front. The
boundary integral representation also does not admit any such discontinuity.
It was therefore decided to consider the case of a ramp loading history, as shown in Figure
qa, b). The time t* is defined as the time for pressure wave to travel the radius of the cavity. Figure
6(a)shows the normalized displacement responses from a number of analyses compared with the
analytical solution for the suddenly applied loading (since the exact solution is not known for the
present case). The results clearly show a delayed response for the ramp loading. All results
(obtained by using time step 0*33t*),including the one for a constant time variation for mesh 1,
seem satisfactory. Figure qb) shows the normalized radial stress at a distance of 0 2 a (a = radius of
the cavity) outside the cavity surface obtained from a series of analyses compared with that of an
abruptly loaded cavity. Although the results for all meshes using the linear variation are similar,
those obtained from a constant time variation are clearly not satisfactory. This is because of the
presence of the boundary acceleration term in the interior stress calculations which cannot be
satisfactorily represented in a constant time interpolation scheme. In Figure qb), for constant

-5

.4 -
.3 - / O 0 0 0
F-@T. n m H II
9 0 3 -.
-2 - k 0 - Analytical [Link] Step
Loading
- Mesh #I, Linear Var.
0.3

.1 - I
A

0
Mesh #l, Constant Var.
Mesh t3, Linear Var.

Loading Curve
-.z - \'
-
2t*
3
-.3
 10970207, 1988, 8, Downloaded from [Link] by Ben Gurion University, Wiley Online Library on [04/07/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
1722 S. AHMAD A N D P. K. BANERJEE

time interpolation, the radial stress obtained at time t = t* is greater than the applied pressure,
which is physically impossible.
Figure 7 shows the time history of the hoop stress at a distance 0.2~ from the surface of the
cavity, where, once again, all linear time variation analyses give consistent results but the constant
time averaging does not, particularly at earlier times. The analytical result shown here once again
is for the case of suddenly applied loading.

(b) Circular loaded area on a halfspace

In this example, a circular area of radius R is subjected to suddenly applied vertical loading
of Po= 1OOO. The mesh for this problem is shown in Figure 8(a). The necessary details of

,
'RI 9R I
(a) Boundary Element Mesh for a Circular Footing on an Elastic Half-space

I&a.

888.

688.

488.

288.

0.
.0

( b ) Time History of Vertical Displacements at the Center and at

the Edge of a Circular Footing

Figure 8. (a) Boundary element mesh for a circular footing on an elastic half-space. (b) Time history of vertical
displacements at the centre and at the edge of a circular footing
 10970207, 1988, 8, Downloaded from [Link] by Ben Gurion University, Wiley Online Library on [04/07/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
TIME-DOMAIN TRANSIENT ELASTODYNAMIC ANALYSIS 1723

geometry and material properties are radius R = 1.0, v=0-333, density= 1.0, pressure and shear
wave velocities of 2.0 and 1.0, respectively, in self-consistent units.
The vertical displacement response (pRu,), where p is the shear modulus, is shown in Figure
8(b) plotted against non-dimensional time t = tc,/R. As expected, the responses of the centre and
the edge are different. Whereas the centre responds very quickly and deflects to the full static
value, the edge takes nearly twice as long to reach its static value. It is also of interest to note the
use of different time steps AT=0-1, 0.15 and 0-2 leads essentially to the same results.
(c) A Jexible square plate foundation on an elastic half space
Karabalis and Beskos,' who described an approximate BEM implementation for half-space
problems, examined the transient dynamic response of a 60 x 60 inch square, 11.517 in thick

(a) Footing Mesh

(b) Half-space Mesh

Figure 9. Boundary element model Tor a flexible square footing on an elastic half-space (Mesh # 1)
 10970207, 1988, 8, Downloaded from [Link] by Ben Gurion University, Wiley Online Library on [04/07/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
1724 S. AHMAD AND P. K. BANERJEE

flexible foundation resting on the surface of an elastic half space. The elastic constants for the
plate were modulus of elasticity E , = 30.004 x lo6 lbjin', Poisson's ratio vp = 0.3 and mass
density pp=7.34 x lb-s2/in4. Those of the half-space were E, =844 x lo6lb/in2, v s = 0 3 ,
p,=2.82 x lb-sz/in.4 The above choice of elastic constants led Karabalis and Beskos" to
obtain a relative plate stiffness ratio K=0.004, where K is defined by
~ , h 3 ( 1 - v,)
K=
12(1 - v ; ) ~ , b
where h and b are the plate thickness and width, respectively, and ,usis the shear modulus of the
half space.

I
(a) Footing Mesh

I 1 I
(b) Half-space Mesh

Figure 10. Boundary element model for a flexible square footing on an elastic, half-space (Mesh #2)
 10970207, 1988, 8, Downloaded from [Link] by Ben Gurion University, Wiley Online Library on [04/07/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
TIME-DOMAIN TRANSIENT ELASTODYNAMIC ANALYSIS 1725

Karabalis and Beskos” used a finite element idealization (4x4mesh) for the plate and a
boundary element idealization (5 x 5 mesh) for the half-space. The boundary elements used were
of constant variation over both space and time. This was obviously a very courageous attempt to
solve such a difficult transient dynamic problem. Those authors probably gained the necessary
confidence by comparing with known static solutions which must have been satisfactory.
In the present attempt at the analyses of the same problem, both the plate and the surface of the
half-space were modelled by BEM. The mesh used for one quarter of the plate region is shown in
Figure 9(a) for the coarse mesh and in Figure l q a ) for the fine mesh. Similarly two sets of meshes
(coarse and fine) for the half-space region are shown in Figures 9(b) and lqb). In all cases,
quadratic elements with linear time variations were used.

.* 30.

-
K-0.004

(a) Loading Curve and the Dimensions of the Footing

.0 4.0 9.0 12.0 16.0 20.0

time (sec. x 4.12E-053
(b) Time History of Vertical Displacement at t h e Center of
the Footing-

Figure 11. (a) Loading curve and the dimensions of the footing. (b) Time history of vertical displacement at the center of
the footing
 10970207, 1988, 8, Downloaded from [Link] by Ben Gurion University, Wiley Online Library on [04/07/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
1726 S. AHMAD AND P. K. BANERJEE

The results of the present analysis for a dynamic pressure pulse shown in Figure ll(a) are
plotted in Figure ll(b), where the results of Karabalis and Beskos” are also shown for
comparison. The results of Karabalis and Beskos” are obviously incorrect, while both sets of
BEM results from the present analysis give almost identical results. It should be noted that in the
present model a plate of finite thickness has resulted in a slight delay in the wave front reaching
the half-space. This, however, does not account for the very large discrepancies in the response
times.
CONCLUDING REMARKS
An advanced algorithm based on the direct boundary element formulation for time-dependent
elastodynamic analysis of three-dimensional solids has been presented. The algorithm is an
unconditionally stable implicit time-marching scheme and is capable of producing very accurate
results. However, for better accuracy, it is recommended that the time step should remain smaller
than L/c,, L being the smallest distance measured along the surface between two corner nodes of
an element. This algorithm is a viable alternative to that based on the finite element methodology,
particularly for soil-structure interaction problems.

ACKNOWLEDGEMENT

The work presented in this paper was made possible by the NASA contract NAS3 23697. The
authors are indebted to Dr C. Chamis, the NASA Project Manager, and Dr E. Todd, the Pratt
and Whitney Project Manager, for their support and encouragement,and to Dr R. B. Wilson and
N. Miller of Pratt and Whitney for valuable discussions.

APPENDIX I

Boundary kernels for transient dynamics

The tensors G i j and F,, are of the form

Gij(x, T ; 6,r )= A6(u - Ar) d l + aij{(l/c2)6(u- r/cl)

A6(u-Ar)dA+(12aij-2bij){6(u-r/cz)

- ( C ~ / C ~ ) ~ ~ ( U+ -2 r~a/i jC/ c~z {)6}’ ( u - r / c 2 ) - ( c z / c l ) 3 6 ’ ( u - r / ~ 1 ) }

-cij(l - 2c~/c~){6(u-r/c,)+(r/cl)~’(u-r/cl)}
- d i j { 6 ( u -r/c2)

+(r/cz)Wo - I/.,)>
1 (A2)
 10970207, 1988, 8, Downloaded from [Link] by Ben Gurion University, Wiley Online Library on [04/07/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
TIME-DOMAIN TRANSIENT ELASTODYNAMIC ANALYSIS 1727

where
v=T-t
dij = yi yj ~mnrn/r'
cij= yjni/r3
dij=(y,nj+Gijymnm)/r3
bIJ. . = c IJ. . + d li
.

APPENDIX I1

Interior stress kernels for transient dynamics

The tensors G t k and FGk are of the form

-(12aijk - 2bijk)[6(u-r/c2) - ( C 2 / C 1 ) 2 S ( U -r / c l ) }
-2raijk/c2{6'(u- r/c2)-( c ~ / c , ) ~ ~r/c,)f
(u-
+ cijk(1 - +
2ci/c:) {S(u - r/cl) (r/cl)d'(u- r/c ,)}

+ d i j k {S ( v - r / ~ 2 )+( ~ / C ~ P ' ( L J1- ~ / C ~ )

1 643)

.2c:(35aijk - 5bijk +
-4ci(45aijk-6b, + cijk) (6 (u -r/c2)-(cz/cl)2S(u -r / c l ) }
-4~,r(lOa~~~-b~~~)(~'(u-r/c~)-(c~/c~)~d'(u-r/c,)}
+ 2c31 -2c~/c~)(3dijk-2eijk){6(u-r/cl)+(r/cl)b'(u-r/cl)}
+ c:(3fjk -2 g i j k ) { -r / c 2 ) +(r/cZ a'(u - r / c 2 ) >
-4aijkr2{S"(u-r/c2)-(c2/cl)"S"(u-r / c l ) )
+ r2(1-2c2,/c:){2c$dijk/c:
+(I - 2 c ~ / c ~ ) e i j k } ~ " ( u - r / crl ~) +j k 6 " ( u - r / c Z )
1 (A41
 10970207, 1988, 8, Downloaded from [Link] by Ben Gurion University, Wiley Online Library on [04/07/2025]. See the Terms and Conditions ([Link] on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
1728 S. A H M A D A N D P. K. BANERJEE

REFERENCES

1. P. K. Banerjee and S. Ahmad, ‘Advanced three-dimensional dynamic analysis by boundary element’, Proc. ASME
Conf: on Advanced Topics in Boundary Element Analysis, Florida, Nov. 1985, AMD, Vol. 12, ASME, 1985, pp. 65-81.
2. P. K. Banerjee, S. Ahmad and G. D. Manolis, ‘Transient elastodynamic analysis of 3-D problems by boundary
element method’, Earthquake eng. s t r u t . dyn., 14,933-949 (1986).
3. S. Ahmad, ‘Linear and nonlinear dynamic analysis by boundary element method’, Ph.D. Thesis, Department of Civil
Engineering, State University of New York at Buffalo, 1986.
4. D. M. Cole, D. D. Koslof and J. B. Minster, ‘A numerical boundary integral equation method for elastodynamics I’,
Bull. seism. SOC. Am., 68, 1331-1357 (1978).
5. [Link], T. Fukui, S. Kato and K. Fujiki, ‘An application of the integral equation method to two-dimensional
elastodynamics’, Theor. Appl. Mech. Tokyo Univ., 28, 281-290 (1980).
6. J. M. Rice and M. H. Sadd, ‘Propagation and scattering of SH-waves in semi-infinite domains using a time-dependent
boundary element method, J . Appl. Mech. ASME, 51, 641-645 (1984).
7. C. C. Spyrakos and D. E. Beskos, ‘Dynamic response of rigid strip foundations by time domain boundary element
method, Int. j . numer. methods eng., 23, 1547-1565 (1986).
8. H. Antes, ‘A boundary element procedure for transient wave propagation in two-dimensional isotropic elastic media’,
Finite Elem. Anal. Des., I, 313-322 (1985).
9. W. J. Mansur and C. A. Brebbia, ‘Transient elastodynamics,’ in C. A. Brebbia (ed.), Topics in Boundary Element
Research, Vol. 2, Springer-Verlag, Berlin, 1985, Chapter 5.
10. W. J. Mansur and C. A. Brebbia, ‘Numerical implementation of the boundary element method for two dimensional
transient scalar wave propagation problems’, Appl. Math. Model, 6, 299-306 (1982).
11. D. L. Karabalis and D. E. Beskos, ‘Dynamic response of 3-D foundations by time domain boundary element method‘,
Final Report Part A, NSF Grant No. CEE-8024725, Department of Civil and Mineral Engineering, University of
Minnesota, Minneapolis, Minnesota, Jan. 1984.
12. D. Graffi, ‘Memorie d e h Academia Scienze, Bologna Series 10, Vol. 4, 1947, pp. 103-109.
13. S. P. Timoshenko and J. N. Goodier, Theory ojElasticity, 3rd edn, McGraw-Hill, New York, 1979, pp. 51C513.