HAL Id: hal-01382196

https://hal.archives-ouvertes.fr/hal-01382196

Submitted on 16 Oct 2016

HAL is a multi-disciplinary open access

archive for the deposit and dissemination of sci-

entic 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 diusion de documents

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

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

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, Wiley, 1983, 19 (1),

pp.73 - 91. �10.1002/nme.1620190109�. �hal-01382196�

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.

1

In a large

number

of these

phenomena,

part

of

the

boundary

is

at

infinity, which

makes the use of integral

equation

methods

particularly attractive, since they

are

capable of

accounting for

the

presence

of

infinite domains. Morse and Feshbach,

2

Eringen

and

Suhubi

3

and

Kupradze

4

are

but

a few of the references

that

contain details

of

the classical work done

in elastodynamics

and

related topics.

Although as

mentioned

above the basic integral equation formulations for wave propagation

problems have

been

known

for

quite

some

time, their adaptation for constructing numerical

algorithms for use in the solution of

boundary

value problems is new. Some of the earliest

such developments were by

Friedman

and

Shaw

5

and

Chen

and

Schweikert,

6

in acoustics,

and

by

Banaugh

and

Goldsmith,

7

in steady-state elastic wave propagation. Since then, many

other

investigators have formulated

and

used numerical implementations known as integral equation

methods (IEM),

boundary

integral

equation

methods

(BIEM)

or

boundary

element methods

(BEM) in various fields of engineering.

In

the field of elastodynamics, Cruse

and

Rizzo

8

and

Cruse

9

derived a

BIEM

approach in conjunction with

the

Laplace transform

and

employed

it

to

solve a half-plane wave propagation problem. A modified version of their approach was

used later

by

Manolis

and

Beskos

10

to

investigate stress concentrations

around

cylindrical

openings of arbitrary

shape

due

to

the

scattering of elastic waves of a general transient nature.

The

determination of

the

steady-state solution

by

BEM

approaches and the reconstitution of

the

transient response by

Fourier

synthesis was

done

by Niwa

et

al.

11

'

12

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.

13

for

the

anti-plane strain case

and

by

Niwa et al.

14

for

the

general two-dimensional

plane

stress/plane

strain

cases.

It

therefore

appears

that

a study comparing

the

various

BEM

approaches

used

to

solve

the

propagation, diffraction

and

scattering

problems

arising in

the

field

of

elastodynamics would

be

of value

to

the

engineering

community.

In

this work,

three

possible

BEM

approaches

are

investigated

and

compared.

All

three

approaches

are used

to

study

the

classical

problem

of

the

scattering

of

a

pressure

wave pulse by a circular cylindrical cavity in a linear elastic,

homogeneous

and

isotropic medium. This is a two-dimensional

problem

of

the

plane

strain

kind.

Analytic-numerical

solutions for

the

displacement, velocity

and

stress fields

at

the

circumference

of

the

cavity

were

obtained

by

Baron

and

Matthews,

15

Baron

and

Parnes

16

and

by

Pao

and

Mow.

17

The

former

two pairs

of

investigators

used

an

integral transform,

separation

on

the

circumferential angle (8)

method.

The

latter

investigators (Pao

and

Mow) used

the

method

of wave function expansion coupled with

the

Fourier

synthesis technique. Miklowitz

18

used a special technique to bring

out

hitherto

undetected

late-time singular, non-decaying

Rayleigh waves

on

the

cavity walls

due

to impulsive loads.

It

is interesting to

note

that

this

classical

problem

has

been

resolved recently

by

other

analytical

means

as, for instance,

the

method

of

matched

asymptotic expansions discussed

by

Datta.

19

The

three

BEM

used

to

resolve this well-documented

problem

are:

1. A time

domain

approach

that

solves

the

problem

as if it

were

a three-dimensional one. In

particular,

the

plane

strain

problem

is

represented

as a circular cylinder of infinite length

and

a numerical integration

is

performed

along

the

axis of

the

cylinder so as

to

collect

the

contributions

to

a

representative

point

from points lying along

the

same

generator.

2.

An

approach

in which

the

integral

equation

formulation is applied

to

the

Laplace transform

of

the

governing

equations

of

motion

for

the

two-dimensional case.

In

this case, a numerical

inverse

transformation

is

required

to bring

the

transformed solution

back

to

the original time

domain.

3.

The

last

approach

is similar

to

case 2, except

that

the

integral transform in question is

the

Fourier

transform

instead

of

the

Laplace transform.

All

of

the

aforementioned

cases

share

some

common

traits. First,

the

direct

BEM

is

employed,

where

the

unknown

boundary

tractions

and

displacements are directly related

to

the known ones.

One

could use

an

indirect

BEM,

as was

done

in Niwa et a!.

11

For

more

details

on

direct

and

indirect

BEM,

reference should

be

made

to

Banerjee

and

Butterfield?

0

Second, all

approaches

utilize a spatial discretization scheme for which

the

boundary

tractions

and

displacements

are

assumed

to

have constant values

over

a given

boundary

segment.

Finally, considerable effort is

directed

towards ensuring

that

all

three

approaches

share

as

many

common

features

as possible, in

order

to

render

the

comparison

study

as fair as possible.

In

the

next few sections,

the

integral

equation

formulation

of

the

general elastodynarnic

problem

along

with

the

numerical

implementation

is discussed. Following that,

the

results

drawn from applying

the

three

aforementioned

BEM

approaches

to

the

same

problem,

which

is shown in

Figure

1,

are

discussed in detail.

STATEMENT

OF

THE

PROBLEM

AND

DEFINITIONS

Let

V

denote

the

area

occupied by a given

body

and

S

denote

its bounding surface. A

co-ordinate

system

(x,

t) is employed,

where

x

denotes

the

cartesian spatial co-ordinates x

1

and

x

2

and

tis

the

time. Associated with

each

surface

point

is

an

outward

pointing unit normal

vector n.

y

Figure

1.

Circular

cavity in

an

infinite

medium

under

the

action

of

an

incoming

pressure

wave

Under

the

assumption

of

small

displacement

theory

and

linear

elastic

isotropic

and

homo~eneous

material

behaviour,

the

equations

of

motion

of

the

body

are

(1)

where

u;(x,

t)

is

the

displacement

vector

and[;

is

the

body

force

vector

per

unit

mass.

The

propagation

velocities

of

the

pressure

and

shear

waves

in

the

body

are

given, respectively, as

ci=(A+2M)/p;

c~=IJ.IP

(2)

where

A

and

IJ.

are

the

Lame

constants

and

p is

the

mass

density.

Furthermore,

commas

indicate

space

differentiation,

dots

indicate

time

differentiation

and

repeated

indices imply

the

summation

convention.

Since in

what

follows

attention

is

restricted

to

two-dimensional

cases,

the

indices

i

and

j

assume

the

values

of

1

and

2, unless

mentioned

otherwise.

The

displacement

vector

uJx,

t)

is

assumed

to

be

twice

differentiable

except

at

possible surfaces

of

discontinuity,

as

discussed in

Eringen

and

Suhubi.

3

The

displacements

and

the

tractions

t; (x,

t)

satisfy

the

boundary

conditions

I;=

U;;n; = p;(X,

t)

ui = q;(x,

t)

where

u;;

is

the

stress

tensor

and

S = S

1

+ Su.

The

displacements

and

velocities also satisfy

the

initial

conditions

U;

(X,

0+) =

U;o(x)

U;(X,

0+) = uiO(x)

XE

V+S

XE

V+S

(3)

(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

cr··=p[(c

2

1-2c2

2

)u

~--+c2

2

(u-

+u··)]

I]

r,,UI]

1,]

],1

(5)

with

8ii

being

the

Kronecker

delta.

The Laplace transform

The

Laplace

transform

of

a function f(x,

t)

is defined as

/(x,

s)

= .P{f(x, 1)} = r f(x,

t)

e

_,

dt

(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:

(

2

2)

- + 2 - +

f-

2 - + . + 0

C1

-c2

Ui,ii

CzUj,ii

j

-s

Uj

UjO

SUjo

=

(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 deriva-

tives,

their

Laplace

transforms

are

very simple

and

become,

respectively,

and

~

=

iJiini

= Pi(x, s)

ui=qi(x,s)

ii"··

= p [(c

2

1-

2c2

2

)u

~--

+ c

2

2(u· ·

+ u ·

·)]

I]

r,,UIJ

1,]

],1

(8)

(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

parameters.

This series

of solutions

must

finally

be

numerically

inverted

back

to

the

original

time

domain.

In

what

follows, a quiescent

state

is

assumed

before

time

t =

o+,

which implies

that

the

initial conditions

are

zero. Also,

the

body

forces

are

taken

equal

to

zero

as well.

These

assumptions

are

made

for

the

sake

of

convenience

and

no

additional conceptual difficulty in

the

solution

procedure

is

encountered

if these assumptions

are

relaxed.

The Fourier transform

The

Fourier

transform

of

a function f(x,

t)

is defined as

ic

x,

w)

= .'F{f(x, t)} = L: t(x,

t)

e -

<w•

dt

(10)

where

w is

the

frequency.

For

the

quiescent

initial conditions

assumed

previously, it

is

elementary

to

see

that

the

Fourier

transform

of

the

governing

equations

of

motion ( 1) will

be

identical

to

equation

(7)

if

the

Laplace

parameter

s

is

set

equal

to

iw in

that

equation.

Therefore,

the

BEM

solution

in the

Fourier

domain

is

identical

to

the

solution

scheme

that

is

to

be

employed for

the

Laplace