...

Document 2097929

by user

on
Category: Documents
4

views

Report

Comments

Transcript

Document 2097929
This thesis is accepted for the degree of Doctor of Philosophy at the Faculty of Engineering and Science,
Aalborg University, Denmark. Defended publicly at Aalborg University, October 29, 1999.
Supervisor:
Søren R. K. Nielsen, dr. techn., Aalborg University, Denmark
Assesment Committee:
Francisco Jose Sánchez-Sesma, Instituto de Ingeniera, UNAM, Coyacan, Mexico
Steen Krenk, Professor, Technical University of Denmark, Denmark
Poul Henning Kirkegaard, Associate Professor, Aalborg University, Denmark
Printed at Aalborg University
i
Acknowledgement
The present thesis "Stress Wave Propagation in Soils Modelled by the Boundary Element
Method" has been prepared in connection with a Ph.D. study carried out in the period
August 1996 to August 1999 at the Department of Building Technology and Structural
Engineering, Aalborg University, Denmark.
Special thanks to my supervisor Reader, dr. techn, Sren R.K. Nielsen for the encouragement and guidance oered during my study. Furthermore, I would like to thank Professor
Francisco Jose Sanchez-Sesma for giving me new inspiration during my 6 month visit to
the Instituto de Ingenieria at Universidad Nacional Autonoma de Mexico and Associate
Professor, Ph.D. Poul Henning Kirkegaard for his assistance with the FEM calculations
for a moving force.
Some of the drawings included in the thesis have been prepared by draughtsperson Norma
Hornung and the proofreading of the manuscript has been done by Senior Secretary
Kirsten Aakjr and their careful and professional work is greatly appreciated.
Finally, I would like to thank colleagues and friends at the department for moral support
and fruitful discussions over the years.
Aalborg, 20. August 1999
Klaus M. Rasmussen
ii
iii
Summary in English
Stress wave propagation in soils has attracted a signicant amount of attention in recent
years. The main reason for this is that this subject is of great importance in numerous elds of application such as earthquake engineering, traÆc induced vibrations and
machinery induced vibrations. Within all of these engineering problems soil-structure
interaction plays an important role. In Chapter 1 the motivation for the thesis will be
further expanded and a short review of the Boundary Element Method (BEM) is given.
Further, the scope of the thesis is dened and a short outline of the chapters is presented.
In Chapter 2 the basic boundary element formulations are derived. Emphasis is placed on
the derivations of Somigliana's identity which is the basic equation for both the direct and
indirect BEM. In the literature simultaneous derivations of both the direct and indirect
BEM are seldom seen, and therefore formulations with both methods will be derived and
practical examples shown.
One of the main diÆculties with the BEM is numerical integration of the kernel functions
with strong singularities. For a 3D indirect boundary element formulation a new analytical integration method for the singularity of Green's function for the displacement eld
is derived in Chapter 3. The analytical solution is used to solve the classical problem
of a Heaviside type loading on an innite 3D isotropic, homogeneous and linear elastic
halfspace. Further improvement of the method is made by a new integration method
for 3D BEM. Usually, integration in 3D BEM is made by ordinary Gauss quadrature in
2D. Subdivision of the elements is used to enhance accuracy of the numerical integration.
The involved parameters are transformed to local polar coordinates, due to the inherent
polar behaviour of the problem. A new wave front based method is proposed where the
integration is directly performed in polar coordinates. This leads to an integration routine
where the causal behaviour of Green's functions can be used in intelligent subdivision of
the integration area of interest.
Coupling of the FEM and BEM gives advantages for e.g. soil-structure interaction problems. A coupled BEM and FEM formulation is presented in Chapter 4 and again both
direct and indirect BEM are used. The principles from the new integration routine for
3D BEM are translated into 2D and implemented in both a direct and indirect 2D BEM
formulation. The singularity problem is solved by a numerical method, and detailed tests
of both the new wave front based integration method and the numerical method for the
singularity problem are performed. The 2D BEM formulation is implemented in the coupled BEM and FEM formulation and tests are performed. Finally it is shown how the
coupled BEM and FEM formulation can easily be used for non-linear problems.
Due to the diÆculties in deriving the kernel functions, the BEM it is not as exible as
e.g. the FEM. But in some cases reformulation of the method can lead to very elegant
solutions for special problems. A new formulation is proposed in Chapter 5 for a moving
load. The method involves redenition of Green's functions, while the boundary element
equations remain the same. The method is tested by the 2D BEM formulation used in
Chapter 4 and the results are compared with a similar solution obtained by the FEM in
convected coordinates with absorbing boundaries. In Chapter 6 the main conclusions of
the thesis are summarised and new contributions are mentioned.
iv
v
Summary in Danish
Udbredelse af spndingsblger i jord er et emne, som i de sidste ar har tiltrukket en del
interesse. Spndingsblgers udbredelse har interesse inden for mange omrader som f.eks.
jordsklv, trak genererede vibrationer eller vibrationer fra maskineri. For alle disse emner er interaktion mellem jord og konstruktion meget vigtigt. I Kapitel 1 motiveres valget
af emne yderligere og en kort litteraturgennemgang prsenteres for randelementmetoden
(BEM). Der gives yderligere en afgrnsning af afhandlingen og en kort disposition.
I Kapitel 2 udledes de grundlggende ligninger for BEM. Vgten er lagt pa udledningen af Somiglianas identitet, som er den grundlggende ligning for bade den direkte og
indirekte BEM. Simultane udledelser af den direkte og indirekte BEM ses sjldent i litteraturen og derfor udledes disse. For begge metoder vises eksempler pa formuleringer af
konkrete problemer.
Et af problemerne med BEM er den numeriske integration af kernefunktionerne, som har
strke singulariteter. En ny analytisk integrationsmetode udledes i Kapitel 3 for Greens
funktion for ytningsfeltet til brug med en 3D indirekte BEM formulering. Den analytiske
metode bruges til at lse det klassiske problem med en Heaviside last pa et uendeligt 3D
isotropisk, homogent og linert elastisk halvrum.
Metoden er yderligere forbedret med en ny integrationsmetode for BEM i 3D. Normalt
udfres integrationerne i BEM ved hjlp af Gauss integration i 2D. Underinddeling af elementerne bruges for at forbedre njagtigheden af den numeriske integration. Problemet
transformeres til polre koordinater pa grund af den grundlggende radire opfrsel. En
ny blgefrontbaseret integrationsmetode prsenteres, hvor integrationen udfres direkte i
polre koordinater. Dette frer til en formulering af problemet, hvor den radire opfrsel
af Greens funktioner udnyttes til intelligent underinddeling af integrationsomraderne.
Kobling af BEM og FEM giver fordele f.eks. i forbindelse med interaktions problemer
mellem jord og konstruktion. En formulering med kobling af BEM og FEM prsenteres
i Kapitel 4, og igen benyttes bade direkte og indirekte BEM. Den nye integrationsrutine fra Kapitel 3 formuleres i 2D og implementeres i bade en direkte og indirekte 2D
BEM formulering. Singularitetsproblemet lses ved hjlp af en numerisk metode, og
detaljerede tests for bade den nye integrationsrutine og den numeriske lsning af singularitetsproblemet udfres. 2D BEM formuleringen implementeres i den koblede BEM og
FEM formulering, og implementeringen testes. Til sidst vises, hvorledes den koblede BEM
og FEM formulering nemt kan bruges for ikke-linere problemer.
Pa grund af problemer med at udlede kernefunktionerne er BEM ikke sa eksibel som
f.eks. FEM, men i nogle tilflde kan en reformulering af metoden fre til meget elegante lsninger for specielle problemer. En ny formulering foreslas i Kapitel 5 for en
bevgelig last. Metoden involverer en ny denition af Greens funktioner, mens selve
BEM ligningerne forbliver de samme. Metoden testes med BEM formuleringen i 2D fra
Kapitel 4, og resultaterne sammenlignes med en lignende lsning fundet ved hjlp af
FEM. Den brugte FEM bruger konvektive koordinater og absorberende randeelementer.
I Kapitel 6 opsummeres de vigtigste konklusioner, og nye bidrag i afhandlingen nvnes.
vi
Contents
1
Introduction
2
Boundary Element Method Formulations
3
BEM Formulation in 3D
1.1 Motivation . . . . . . . . . . . . . . . .
1.2 Review of BEM . . . . . . . . . . . . .
1.2.1 Direct BEM . . . . . . . . . . .
1.2.2 Indirect BEM . . . . . . . . . .
1.2.3 Boundary Methods . . . . . . .
1.2.4 T-matrix Methods . . . . . . .
1.2.5 Dual Reciprocity BEM . . . . .
1.2.6 Symmetric BEM (Variational) .
1.2.7 Coupled BEM/FEM . . . . . .
1.2.8 Applications . . . . . . . . . . .
1.3 Scope of Thesis . . . . . . . . . . . . .
1.4 Thesis Outline . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Basic BEM Time Domain Formulation . . . . . . . . . . . . . . . . . .
2.3 Discretization in Time and Space . . . . . . . . . . . . . . . . . . . . .
2.3.1 Space Discretization . . . . . . . . . . . . . . . . . . . . . . . .
2.4 Numerical Integration . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.1 Singularity for Green's Functions . . . . . . . . . . . . . . . . .
2.5 Time Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.6 Direct Boundary Element Method . . . . . . . . . . . . . . . . . . . . .
Example 2.1 Direct BEM Formulation for Alluvial Valley . . . . . . . .
2.7 Indirect Boundary Element Methods . . . . . . . . . . . . . . . . . . .
Example 2.2 Single-Layer Potential Formulation for Elastic Halfspace .
2.8 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Approximation of Single-Layer Potential . . . . . . . . . . . . . . . . . . .
3.2.1 Singularity Problems of Green's Function for the Displacement Field
3.2.2 Analytical Solution of Singularity Problem . . . . . . . . . . . . . .
3.2.3 Singularity Problem of Green's Function for the Traction Field . . .
vii
1
1
3
3
4
4
4
5
5
5
6
7
7
9
9
9
14
14
16
17
17
18
19
22
23
24
25
25
25
26
27
29
viii
Contents
Example 3.1 Results for 3D Problem Subjected to a Concentrated Heaviside Load . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 A new Integration Scheme for Surface Problems . . . . . . . . . . . . . . .
Example 3.2 Results with new Integration Scheme . . . . . . . . . . . . .
3.4 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
29
31
34
36
4
Coupled BEM and FEM Solutions
39
5
Boundary Element Solution for a Moving Force
6
Conclusions
91
Bibliography
95
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Coupling of BEM and FEM Formulations . . . . . . .
4.3 2D Time Domain Boundary Element Formulation . . .
4.3.1 Space Discretization . . . . . . . . . . . . . . .
4.3.2 Numerical Space Integration . . . . . . . . . . .
4.3.3 Singularity Problem . . . . . . . . . . . . . . .
4.3.4 Enclosing Elements . . . . . . . . . . . . . . . .
Example 4.1 Tests of Singularity Procedure . . . . . .
Example 4.2 Tests of Numerical Space Integration . . .
4.3.5 Time Discretization . . . . . . . . . . . . . . . .
Example 4.3 Results from 2D FEM Implementation . .
4.4 Results from 2D BEM Implementation . . . . . . . . .
Example 4.4 Heaviside Load on Elastic Halfspace . . .
Example 4.5 Lamb's Problem in 2D . . . . . . . . . .
Example 4.6 Direct BEM Compared to Indirect BEM .
4.5 Results from Coupled BEM/FEM Implementation . . .
Example 4.7 Heaviside Load on Elastic Halfspace . . .
Example 4.8 Heaviside Load on Layered Media . . . .
4.6 Non-linear FEM Formulation . . . . . . . . . . . . . .
4.7 Concluding Remarks . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 Formulation of Boundary Element Equations for a Moving Force .
5.3 Time Approximation . . . . . . . . . . . . . . . . . . . . . . . . .
Example 5.1 Results from 2D BEM Solution for Moving Force . .
Example 5.2 Comparison with Corresponding FEM Solutions . .
5.4 Concluding Remarks . . . . . . . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
39
41
45
46
47
47
48
50
53
57
58
60
60
63
64
65
65
68
73
75
77
77
78
79
81
85
89
6.1 Summary of Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
6.2 Overall Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
Appendices
A Integral of Normal Vectors over Semicircular Surface
103
Contents
ix
B Element Types and Gauss Quadrature
B.1 Element Types in Rectangular Coordinates . . . . . . . . . .
B.1.1 2-Node Line Element . . . . . . . . . . . . . . . . . .
B.1.2 3-Node Triangular Element . . . . . . . . . . . . . .
B.1.3 8-Node Quadrilateral Element . . . . . . . . . . . . .
B.1.4 6-Node Triangular Element . . . . . . . . . . . . . .
B.2 Element Types in Cylindrical Coordinates . . . . . . . . . .
B.2.1 8-Node Quadrilateral Element . . . . . . . . . . . . .
B.2.2 6-Node Triangular Element . . . . . . . . . . . . . .
B.3 Gauss Quadrature . . . . . . . . . . . . . . . . . . . . . . . .
B.3.1 3rd Order Gauss Quadrature . . . . . . . . . . . . . .
B.3.2 3rd Order Gauss Quadrature on Triangular Element .
C Analytical Solutions
C.1 Lamb's Problem in 2D . . . . . . . . .
C.1.1 Solution for Distributed Loads .
C.2 Lamb's Problem in 3D . . . . . . . . .
C.2.1 Vertical Displacement . . . . .
C.2.2 Horizontal Displacement . . . .
C.2.3 Solution for Distributed Load .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
105
. 105
. 105
. 105
. 107
. 108
. 109
. 109
. 110
. 111
. 111
. 112
113
. 113
. 115
. 117
. 117
. 118
. 119
D Newmark Integration Scheme
121
E Sparse Formulation
123
F
127
Example E.1 3 3 Matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . 124
Example E.2 FEM Stiness Matrix . . . . . . . . . . . . . . . . . . . . . . 124
E.1 Sparse Solver . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
Integration of Functions for Moving Force
x
Contents
List of Figures
1.1 Train induced vibrations have gained a signicant amount of attention in
recent years due to increased velocity. . . . . . . . . . . . . . . . . . . . . .
1.2 The probably most classical soil-structure interaction problems are earthquake problems. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.1 Cartesian coordinate system. . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Boundary around singular point. . . . . . . . . . . . . . . . . . . . . . . .
2.3 Mapping from parent to global domain. . . . . . . . . . . . . . . . . . . . .
2.4 Element with subdivision. . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.5 Shape functions for time step l. . . . . . . . . . . . . . . . . . . . . . . . .
2.6 Alluvial valley. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.7 Elastic halfspace with coordinate system. . . . . . . . . . . . . . . . . . . . . . .
3.1 Circular plane element mesh for halfspace. . . . . . . . . . . . . . . . . . .
3.2 a) Green's function for the traction eld for component 1,2,1. b) Section
through x2 x3 plane. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Vertical displacement uz at (r; ; z)=(0.90 m, 0 m, 0 m). . . . . . . . . . . . . . . .
3.4 Horizontal displacement ur at (r; ; z)=(0.90 m, 0 m, 0 m). . . . . . . . . . . . . .
3.5 Horizontal displacement perpendicular to the wave front u at (r; ; z)=(0.90 m, 0 m, 0
m). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.6 Vertical displacement uz at (r; ; z)=(0.90 m, 0 m, 0 m). . . . . . . . . . . . . . . .
3.7 Horizontal displacement ur at (r; ; z)=(0.90 m, 0 m, 0 m). . . . . . . . . . . . . .
3.8 A wave front emerging from source point x. . . . . . . . . . . . . . . . . .
3.9 Section of the wave front over which the integration is performed. rref =
(l 1)tc2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.10 The wave front on a single element x. . . . . . . . . . . . . . . . . . . . . .
3.11 Wave front based integration. . . . . . . . . . . . . . . . . . . . . . . . . .
3.12 Vertical displacement uz at (r; ; z)=(0.90 m, 0 m, 0 m) with the new integration scheme.
3.13 Horizontal displacement ur at (r; ; z)=(0.90 m, 0 m, 0 m) with the new integration
scheme. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.14 Horizontal displacement perpendicular to the wave front u at (r; ; z)=(0.90 m, 0 m, 0
m) with the new integration scheme. . . . . . . . . . . . . . . . . . . . . . . . .
3.15 Vertical displacement uz at (r; ; z)=(0.90 m, 0 m, 0 m). . . . . . . . . . . . . . . .
3.16 Horizontal displacement ur at (r; ; z)=(0.90 m, 0 m, 0 m). . . . . . . . . . . . . .
xi
1
2
9
11
15
16
18
19
23
25
29
30
30
30
31
31
32
33
33
34
35
35
35
36
36
xii
List of Figures
4.1 Examples of BEM/FEM coupling. a) Non-linear alluvial valley interacting
with linear elastic media for an incoming stress wave eld. b) Same as a)
but with impulse load applied within the medium. c) Interaction between
linear elastic soil with incoming stress wave eld and structure with impulse
loads. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Denition of local axis . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 Halfspace model with enclosing boundary. . . . . . . . . . . . . . . . . . .
4.4 Model boundary for test example. . . . . . . . . . . . . . . . . . . . . . . . . .
4.5 Halfspace model. A) Circular enclosing boundary without line segments B) Circular
4.6
4.7
4.8
4.9
4.10
4.11
4.12
4.13
4.14
4.15
4.16
4.17
4.18
4.19
enclosing boundary with line segments C ) Square enclosing boundary with line segments.
a) Normalized average errors for 10 integration subdivsions. b) Normalized average
errors with 11 nodes on the enclosing boundary.
Normalized average errors for 10 integration subdivsions. The y-axis is logarithmically
scaled.
t1;2 multiplied by the shape function for second node for x at node 6, step 2. Numbers
within circles are element numbers and dashed lines signify element boundaries. Wave
boundaries are shown as asterisks.
Element mesh for !n = 140 rad/s ) h = 3000 in generated with automatic mesh
generation tool. Length is given in inches.
u1 found from FEM analysis. [!n] =rad/s. A : (x1 ; x2 ) =(0 in, 0 in); B : (x1 ; x2 )=(6000
in, 0 in); C : (x1 ; x2 )=(12000 in, 0 in).
u2 found from FEM analysis. [!n] =rad/s. A : (x1 ; x2 ) =(0 in, 0 in); B : (x1 ; x2 )=(6000
in, 0 in); C : (x1 ; x2 )=(12000 in, 0 in).
u1 found from BEM analysis. B : (x1 ; x2 )=(6000 in, 0 in); C : (x1 ; x2 )=(12000 in, 0 in).
u2 found from BEM analysis. A : (x1 ; x2 ) =(0 in, 0 in); B : (x1 ; x2 )=(6000 in, 0 in); C :
(x1 ; x2 )=(12000 in, 0 in).
u1 found from BEM and FEM analysis compared with the results obtained by Antes.
B : (x1 ; x2 )=(6000 in, 0 in); C : (x1 ; x2 )=(12000 in, 0 in).
u2 found from BEM and FEM analysis compared with the results obtained by Mansur
and Antes and the analytical solution. A : (x1 ; x2 ) =(0 in, 0 in); B : (x1 ; x2 )=(6000 in,
0 in); C : (x1 ; x2 )=(12000 in, 0 in).
Numerical and analytical solution for u1 at point (x1 ; x2 )=(12000 in, 0 in).
Numerical and analytical solution for u2 at point (x1 ; x2 )=(12000 in, 0 in).
u1 for Heaviside type load on elastic halfspace. B : (x1 ; x2 )=(6000 in, 0 in); C :
(x1 ; x2 )=(12000 in, 0 in).
u2 for Heaviside type load on elastic halfspace. A : (x1 ; x2 ) =(0 in, 0 in); B : (x1 ; x2 )=(6000
in, 0 in); C : (x1 ; x2 )=(12000 in, 0 in).
40
46
49
50
50
. . . . . . . . . . . . . . . . . . . 52
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
. . . . . . . . . . . . . . . . . . . . . . . . . 55
. . . . . . . . . . . . . . . . . . . . . . 58
. . . . . . . . . . . . . . . . . . . . . . . . 59
. . . . . . . . . . . . . . . . . . . . . . . . 60
61
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
. . . . . . . . . . . . . . 62
. . . . . . . . . . . . . . . . . . . . . . . . . 62
. . . . . . 63
. . . . . . 64
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
........................
4.20 BEM/FEM coupling for alluvial valley, examples. . . . . . . . . . . . . . .
4.21 Example of coupled element mesh. a) Coupled BEM and FEM mesh. b) FEM mesh
alone. c) BEM mesh alone. Length in inches. . . . . . . . . . . . . . . . . . . . .
4.22 u found from BEM/FEM analysis compared with the results obtained by Antes. B :
(x ; x )=(6000 in, 0 in); C : (x ; x )=(12000 in, 0 in). . . . . . . . . . . . . . . . .
4.23 u found from BEM/FEM analysis compared with analytical results. A : (x ; x ) =(0
in, 0 in); B : (x ; x )=(6000 in, 0 in); C : (x ; x )=(12000 in, 0 in). . . . . . . . . . .
1
1
2
1
2
2
1
1
2
1
2
2
65
66
66
67
67
List of Figures
4.24
4.25
4.26
4.27
4.28
4.29
4.30
4.31
4.32
4.33
4.34
4.35
xiii
+
+
+
for surface points, CASE 1, BEM/FEM analysis. +
= .
1 = 1 = 2 = 2 = Thick lines indicate the edges of the alluvial valley.
+
+
+
u2 for surface points, CASE 1, BEM/FEM analysis. +
= .
1 = 1 = 2 = 2 = Thick lines indicate the edges of the alluvial valley.
u1 for points on the boundary between layers 1 and 2, CASE 1, BEM/FEM analysis.
+
+
+
+
= .
1 = 1 = 2 = 2 = u2 for points on the boundary between layers 1 and 2, CASE 1, BEM/FEM analysis.
+
+
+
+
= .
1 = 1 = 2 = 2 = +
+
+
u1 for surface points, CASE 2, BEM/FEM analysis. +
1 = 1 = 2 = 2 = 0:1 and = . Thick lines indicate the edges of the alluvial valley.
+
+
+
u2 for surface points, CASE 2, BEM/FEM analysis. +
1 = 1 = 2 = 2 == 0:1 and = . Thick lines indicate the edges of the alluvial valley.
u1 for points on the boundary between layers 1 and 2, CASE 2, BEM/FEM analysis.
+
+
+
+
and = .
1 = 1 = 2 = 2 == 0:1 u2 for points on the boundary between layers 1 and 2, CASE 2, BEM/FEM analysis.
+
+
+
+
and = .
1 = 1 = 2 = 2 = 0:1 +
+
+
u1 for surface points, CASE 3, BEM/FEM analysis. +
1 = 1 = 0:1 , 2 = 2 =
0:5 and = . Thick lines indicate the edges of the alluvial valley.
+
+
+
u2 for surface points, CASE 3, BEM/FEM analysis. +
1 = 1 = 0:1 , 2 = 2 =
0:5 and = . Thick lines indicate the edges of the alluvial valley.
u1 for points on the boundary between layers 1 and 2, CASE 3, BEM/FEM analysis.
+
+
+
+
and = .
1 = 1 = 0:1 , 2 = 2 = 0:5 u2 for points on the boundary between layers 1 and 2, CASE 3, BEM/FEM analysis.
+
+
+
+
and = .
1 = 1 = 0:1 , 2 = 2 = 0:5 u1
. . . . . . . . . . . . . . . . . . 68
. . . . . . . . . . . . . . . . . . 68
. . . . . . . . . . . . . . . . . . . . . . . . . . 69
. . . . . . . . . . . . . . . . . . . . . . . . . . 69
. . . . . . . . . . 70
. . . . . . . . . . 70
. . . . . . . . . . . . . . . . . . . 71
. . . . . . . . . . . . . . . . . . . . 71
. . . . . . . 72
. . . . . . . 72
. . . . . . . . . . . . . . . . 73
................
5.1 Moving force in xed and moving coordinate systems. . . . . . . . . . . . .
5.2 u from Heaviside load. v = 0:01 cR and v = 0 . B : (~x ; x~ )=(6000 in, 0 in);
C : (~
x ;x
~ )=(12000 in, 0 in). . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3 u from Heaviside load. v = 0:01 cR and v = 0 . A : (~x ; x~ ) =(0 in, 0 in); B :
(~x ; x~ )=(6000 in, 0 in); C : (~x ; x~ )=(12000 in, 0 in). . . . . . . . . . . . . . . . .
5.4 u from Heaviside load. v = 0:01 cR , v = 0:2 cR and v = 0:5 cR . B :
(~x ; x~ )=(-6000 in, 0 in); C : (~x ; x~ )=(-12000 in, 0 in). . . . . . . . . . . . . . . .
5.5 u from Heaviside load. v = 0:01 cR , v = 0:2 cR and v = 0:5 cR . A :
(~x ; x~ ) =(0 in, 0 in); B : (~x ; x~ )=(-6000 in, 0 in); C : (~x ; x~ )=(-12000 in, 0 in). . . .
5.6 u from Heaviside load. v = 0:01 cR , v = 0:2 cR and v = 0:5 cR . B :
(~x ; x~ )=(6000 in, 0 in); C : (~x ; x~ )=(12000 in, 0 in). . . . . . . . . . . . . . . . .
5.7 u from Heaviside load. v = 0:01 cR , v = 0:2 cR and v = 0:5 cR . A :
(~x ; x~ ) =(0 in, 0 in); B : (~x ; x~ )=(6000 in, 0 in); C : (~x ; x~ )=(12000 in, 0 in). . . . .
5.8 a) BEM and b) FEM meshes. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.9 Horizontal displacements at (~x ; x~ )=(+20 m, 0 m) obtained using BEM (top) and FEM
1
1
1
2
1
1
1
2
1
5.11
83
1
1
83
1
84
1
84
85
2
1
1
2
1
1
1
2
1
2
82
1
1
1
2
2
82
1
2
1
5.10
1
2
1
2
2
2
1
2
1
1
2
1
1
1
2
2
1
2
1
1
1
1
73
78
2
(bottom) for v1 = 0 c2 (BEM v1 = 0:001 c2 ) , v1 = 0:2 c2
and v1 = 0:5 c2 .
Vertical displacements at (~x1 ; x~2 )=(+20 m, 0 m) obtained using BEM (top) and FEM
(bottom) for v1 = 0 c2 (BEM v1 = 0:001 c2 ) , v1 = 0:2 c2
and v1 = 0:5 c2 .
Horizontal displacements at (~x1 ; x~2 )=(-20 m, 0 m) obtained using BEM (top) and FEM
(bottom) for v1 = 0 c2 , v1 = 0:2 c2
and v1 = 0:5 c2 .
86
86
. . . . . . . . . . 87
xiv
5.12
5.13
5.14
List of Figures
Vertical displacements at (~x1 ; x~2 )=(-20 m, 0 m) obtained using BEM (top) and FEM
(bottom) for v1 = 0 c2 , v1 = 0:2 c2
and v1 = 0:5 c2 .
Horizontal displacements at (~x1 ; x~2 )=(+4 m, 0 m) obtained using BEM (top) and FEM
(bottom) for v1 = 0:5 c2 , v1 = 0:7 c2
and v1 = 0:9 c2 .
Vertical displacements at (~x1 ; x~2 )=(+4 m, 0 m) obtained using BEM (top) and FEM
(bottom) for v1 = 0:5 c2 , v1 = 0:7 c2
and v1 = 0:9 c2 .
. . . . . . . . . . 87
. . . . . . . . . 88
. . . . . . . . . 88
A.1 Hemisphere with coordinate system and normal vector. . . . . . . . . . . . 103
B.1 2-node line element. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
B.2 3-node triangular element. . . . . . . . . . . . . . . . . . . . . . . . . . . . 106
B.3 8-node quadrilateral element. . . . . . . . . . . . . . . . . . . . . . . . . . 107
B.4 6-node triangular element. . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
B.5 8-node quadrilateral element in cylindrical coordinates. . . . . . . . . . . . 109
B.6 6-node triangular element in cylindrical coordinates. . . . . . . . . . . . . . 110
B.7 Gauss points for 3rd order Gauss quadrature on quadrilateral element. . . 111
B.8 Gauss points for triangular element. . . . . . . . . . . . . . . . . . . . . . . 112
C.1 u1 for Lamb's problem in 2D, impulse. . . . . . . . . . . . . . . . . . . . . 114
C.2 u2 for Lamb's problem in 2D, impulse. . . . . . . . . . . . . . . . . . . . . 114
C.3 u1 for Lamb's problem in 2D, Heaviside load. . . . . . . . . . . . . . . . . . 114
C.4 u2 for Lamb's problem in 2D, Heaviside load. . . . . . . . . . . . . . . . . . 115
C.5 u1 for Lamb's problem in 2D for a distributed load. . . . . . . . . . . . . . 115
C.6 u2 for Lamb's problem in 2D for a distributed load. . . . . . . . . . . . . . 116
C.7 u1 for Lamb's problem in 2D with a distributed Heaviside load. . . . . . . 116
C.8 u2 for Lamb's problem in 2D with a distributed Heaviside load. . . . . . . 116
C.9 Vertical displacement uz at the surface of a Poisson material from a downward directed unit concentrated force H (t). . . . . . . . . . . . . . . . . . 118
C.10 Horizontal displacement ur at the surface from a downward directed unit
concentrated force H (t). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
C.11 Vertical displacement uz at the surface from a downward directed unit force
H (t) applied over a nite area. . . . . . . . . . . . . . . . . . . . . . . . . . 119
C.12 Horizontal displacement ur at the surface from a downward directed unit
force H (t) applied over a nite area. . . . . . . . . . . . . . . . . . . . . . 120
List of Tables
4.1
4.2
4.3
4.4
4.5
4.6
4.7
4.8
4.9
4.10
B.1
E.1
Average errors of diagonal terms for type A enclosing boundary if the last node at each
end of the element mesh is not included in the average errors. The error is dened as
jval 0:5j , where val is the numerical value.
error =
0:5
Diagonal terms for type A with 22 nodes on the enclosing boundary and 25 integration
subdivisions.
Average relative errors of diagonal terms for type B enclosing boundary. f = 41(12 ) .
The error is dened as error = jval0:50:5j , where val is the numerical value.
Average relative errors of diagonal terms for type C enclosing boundary. f = 41(12 ) .
The error is dened as error = jval0:50:5j , where val is the numerical value.
g1;1 multiplied by shape function for second node integrated over elements for x at node
6, step 1.
t1;2 multiplied by shape function for second node integrated over elements for x at node
6, step 1.
g1;1 multiplied by shape function for second node integrated over elements for x at node
6, step 2.
t1;2 multiplied by shape function for second node integrated over elements for x at node
6, step 2.
g1;1 multiplied by shape function for second node integrated over elements for x at node
6, step 3.
t1;2 multiplied by shape function for second node integrated over elements for x at node
6, step 3.
. . . . . . . . . . . . . . . . . 51
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
. . . . . . . 52
...
.................................
.................................
.................................
.................................
.................................
.................................
Gauss points for 3rd order Gauss quadrature on triangular element.
Memory use for sparse stiness matrix from FEM program. . . . . . . . . .
xv
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
. 53
. 54
. 54
. 54
. 55
. 56
. 56
. 112
. 124
xvi
List of Tables
Chapter 1
Introduction
1.1
Motivation
Stress wave propagation in soils has attracted a signicant amount of attention in recent
years. The main reason for this is that this subject is of great importance in numerous elds of application such as earthquake engineering, traÆc induced vibrations and
machinery induced vibrations. Within all of these engineering problems soil-structure
interaction plays an important role.
Figure 1.1: Train induced vibrations have gained a signicant amount of attention in
recent years due to increased velocity.
The main diÆculty with wave propagation in soils and soil-structure interaction problems
is the complexity of the problems. First of all the inhomogeneous composition of the soil
provides major diÆculties. Secondly, the fact that soil in practical cases must be viewed
as an unbounded medium provides problems for most numerical methods. The problem
1
2
1. Introduction
is mainly visible for soil-structure interaction problems, where a highly bounded structure
is connected to an unbounded medium, i.e. the soil.
Due to the complexity of the problems analytical solution methods can only be obtained
for very simplied and not practicable problems. Therefore, the main solution methods for these types of problems are numerical methods or semi-analytical methods. The
main numerical method used within structural engineering is the Finite Element Method
(FEM). But for wave propagation in soils the FEM has problems due to the unbounded
property of the soil. Waves propagating in a bounded FEM mesh will be reected at
the boundary. There is a variety of solutions with absorbing boundary conditions for the
FEM but they all represent numerical approximations to the problem and only ensure
accuracy under special conditions. Therefore, the Boundary Element Method (BEM) is
an attractive alternative to the FEM. The BEM has the advantage that it reduces the
dimension of the problem by one, e.g. from 3 to 2. Furthermore, and most importantly,
the BEM is capable of producing accurate results for unbounded media due to its inherent
ability to satisfy the radiation condition. But the BEM has the disadvantage that it requires extensive calculations in the time domain, and hardly can be applied to non-linear,
non-homogeneous or non-isotropic media. Therefore, an obvious solution to soil-structure
interaction problems is to use a coupled FEM/BEM.
Figure 1.2: The probably most classical soil-structure interaction problems are earthquake
problems.
1.2. Review of BEM
1.2
3
Review of BEM
oundary Integral Equation Methods (BIEM) have been known since the start of this
century but only when computers became available in the sixties did they really take o as
Boundary Element Methods (BEM). The BEM was initiated in the early sixties and the
paper by Jaswon from 1963 [32] is considered the rst explicit and general introduction
of the direct BEM. Some of the early work with the BEM envolving elastodynamics was
done by Cruse & Rizzo in the late sixties [16].
B
Since the subject of this thesis is wave propagation in soils and some of the the problems
which are of interest involve transients it has been chosen to use time domain formulations.
Thus, this review will focus on time domain formulations but brief sections on frequency
domain formulations will be included where appropriate. The main focus on this brief
review will be of the direct and indirect BEM but other methods will be mentioned, and
examples of a few applications are given.
1.2.1
Direct BEM
The direct BEM applied to elastodynamics starts with an integral representation of the
displacement eld, which can be obtained through the use of the Betti reciprocal theorem.
This integral representation of the displacements and surface tractions with the aid of the
fundamental solutions, Green's functions, is moved to the boundary through a limiting
process. The nal boundary integral equation, often denoted the Somigliana identity, is
solved numerically for the unknown displacements and tractions. The basic derivations of
the elastodynamic equations in the time and frequency domains can be found in numerous
references, see e.g. Eringen and Suhubi [17].
The numerical formulation of the direct BEM has followed the same basic principles since
the mid-eighties, see e.g. Banerjee et al. [10]. The BEM involves numerical integration
of Green's functions times the chosen shape functions. The integrals which must be evaluated are by some authors denoted hypersingular, strongly singular or weakly singular if
the integrand includes terms which are O(1=r3), O(1=r2) or O(1=r) in 3D and O(1=r2),
O(1=r) or O(ln r) in 2D. The weak integrals can be solved using normal numerical techniques, but the strongly singular and hypersingular integrals need special methods. The
evaluation of the singular integrals is either done directly interpreting them as Cauchy
principal values, see Giuggiani [25], or using regularization to convert them to weak singular integrals, e.g. Krishnasamy et al. [41]. The regularization techniques include the
socalled rigid body motion principle, where the static Green function is used to convert
the integration of Green's function for traction eld to a weak integral, see e.g. Manolis
& Beskos [49]. The rigid body motion principle requires closed surfaces and thus when
applied to halfspaces the socalled enclosing elements should be employed, see Banerjee &
Ahmad [10]. Several other methods used to solve the singularity problems exist, see e.g.
the review by Krishnasamy et al. [41].
There is a very large number of papers dealing with time domain BEM formulations for
elastodynamics. A few of these papers, which will be cited later in the thesis, are Ahmad
4
1. Introduction
& Banerjee [3], Antes [5], Beskos [14], Israil & Banerjee [31], Kobayashi [36] and Mansur
[50].
The basic, direct BEM formulation uses the full space Green functions but special methods have been developed using halfspace Green's functions. Various versions of halfspace
Green's functions have been presented but they all have the disadvantage of a very high
complexity which means that very large computational eort is necessary, see e.g. Song
& Wolf [72].
1.2.2
Indirect BEM
The direct BEM is called direct because it deals with quantities of physical signicance,
such as displacements and tractions and because it represents a direct discretization of the
Somigliana identity. The indirect BEM also represents a discretization of the Somigliana
identity but it formulates the problem in terms of some ctitious quantities denoted potentials. The indirect BEM is not very popular compared to the direct version despite the
fact that it has a longer history than the direct version and that it is closely related to
the classical work on integral equation methods. This approach has not been used very
much for time domain elastodynamic problems, but some work exists on the subject, and
Antes et al. [6] made a overview paper on the subject in 1991. The indirect BEM has
mainly been used in the frequency domain where a special, indirect BEM has been developed for the solution of wave scattering problems under conditions of plane or anti-plane
strain harmonic motion. The elastodynamic problems of diraction of SH-, SV-, P- and
Rayleigh waves by surface or underground irregularities have been studied. Some of the
pioneer work on this method was done by Sanchez-Sesma & Esquivel [68] and Wong [85]
in the late seventies. Later the method has been further developed in both 2D and 3D,
see e.g. Sanchez-Sesma & Luzon [69] and Sanchez-Sesma et al. [70, 71].
1.2.3
Boundary Methods
A class of boundary methods exists where the solution of an elastodynamic problem in the
frequency domain is expressed as a linear combination of functions. The functions satisfy
the governing dierential equation of the problem and must form a C-complete system.
The problem is then solved by letting these functions satisfy the boundary conditions.
Some of the contributions to this method were given by Herrera [29], Sanchez-Sesma [67]
and Eshraghi & Davinski. [18].
1.2.4
T
-matrix Methods
The T-matrix method was initially developed for acoustic problems but has later been
applied to a number of elastodynamic problems. A known incident eld and an unknown
scattered eld are both expanded in terms of a set of basis functions (e.g. spherical harmonics) and the transition or T-matrix connecting the unknown expansion coeÆcients is
1.2. Review of BEM
5
constructed by an integral representation using Green's function also expanded in terms
of the basis functions. A description of the method can be found in Pao [59]. To the
author's knowledge, these methods have been abandoned by most people since the early
nineties.
1.2.5
Dual Reciprocity BEM
The dual reciprocity BEM, or DR/BEM as it is often denoted, was rst introduced in
1982 in a paper by Nardini & Brebbia [52]. The method represents a way of dealing
with equations in BEM for which a complete fundamental solution is either unavailable
or inconvenient. In the case of elastodynamics the method uses the much simpler elastostatic fundamental solution which creates an inertial volume integral. The reciprocal
theorem is then applied for the second time which converts the volume integral into a
boundary integral. The problem with the method is the lack of accuracy, but in recent
years a great deal of eort has been put into improvements of this method. One of the
main research areas for this method is the socalled radial basis approximation functions
which is well known by mathematicians but is a new concept for engineers. The choice of
these functions aects the accuracy of the method but depends on the specic problem
at hand, see e.g. Partridge et al. [60].
1.2.6
Symmetric BEM (Variational)
All the direct and indirect BEMs give non-symmetric matrices which increase the computational work of solving the systems. When using variational methods it is possible to
construct BEMs with symmetric matrices. The problem is that these methods involve
double spatial integrals instead of the usual single spatial integral. The computational
work of the extra spatial integration often exceeds the gain of solving symmetric matrices
making the method unnecessary. There are several variational formulations applied to
elastodynamics, see e.g. Maier et al. [47] and Gaul et al. [24].
1.2.7
Coupled BEM/FEM
In some problems it is avantageous to combine the BEM and FEM in coupled BEM/FEM
formulations. When coupling BEM and FEM formulations the advantages of both methods are combined and some of the disadvantages are eliminated. The coupling is usually
done through equilibrium and compatibility at the interface of the two domains modelled by the two methods. In the context of this thesis, time domain BEM/FEM coupled
methods used for soil-structure interaction problems are mainly of interest. A number of
coupled formulations exist for time domain soil-structure interaction problems in 2D and
3D. A few of the main contributions are due to Abe & Yoshida [1] who looked at alluvial
valleys, Belytschko & Lu [12], who used a variational formulation creating symmetric
matrices. The early contributions of the principle may be attributed to Fukuri & Ishida
6
1. Introduction
[21] and the works of Von Estor [78] and Von Estor & Kausel [81].
1.2.8
Applications
Time domain BEM formulations have been applied to most elastodynamic problems, e.g.:
Halfspace problems.
Cavities in full or halfspaces.
Surface or subsurface topography diraction.
Diraction by cracks.
Beams and membranes.
Plates.
Shallow shells.
Foundations and above ground structures.
Piles and pile groups.
Underground structures.
Fluid-structure interaction.
Viscoelasticity.
Inhomogeneity.
Anisotropy.
Contact analysis.
Optimum design and control.
Review of papers including these problems will exceed the scope of this review, see e.g.
the reviews by Beskos [13, 14].
1.3. Scope of Thesis
1.3
7
Scope of Thesis
As mentioned in Section 1.1, wave propagation in soils requires modelling of innite
domains. The FEM has problems with the modelling of absorbing boundaries whereas
the BEM has an inherent ability to solve these kinds of problems due to the fullment
of the radiation conditions. Problems involving wave propagation induced by traÆc are
intrinsic in nature and the obvious solution is therefore to apply time domain solutions.
It has been chosen in this thesis to investigate the use of time domain BEM soil-structure
interaction problems, where a moving or stationary structure is placed on the surface.
Either the structure is responsible for the generation of stress waves in the soil or stress
waves propagating in the soil coincide with the structure. In the thesis it has been chosen
to look closer at a number of dierent aspects of these kinds of problems. The following
items will be included in the thesis
Derivation of the time domain BEM equations for both the direct and indirect
formulation.
Formulation of the BEM in 2D and 3D.
Analytical and numerical methods to solve the singularity problems.
Coupling of BEM and FEM for soil-structure interaction or non-linear problems.
Development of BEM dealing with moving loads.
Implementation of the above and testing of the methods.
1.4
Thesis Outline
Here a very short outline of the thesis is presented.
Chapter 2
The basic direct and indirect BEM equations are derived in 3D and general principles of the numerical procedures are presented. For both the direct and indirect
BEM an example is shown for a specic stress wave propagation problem in soil.
Chapter 3
A 3D BEM formulation in polar coordinates is presented. The singularity
problem for Green's function for the displacement eld is solved by an analytical solution,
and a new integration principle is presented. Finally, the method is tested on a classic
problem from elastodynamics.
Chapter 4
In this chapter a coupled BEM/FEM formulation is shown and 2D BEM
direct and indirect formulations are applied to the formulation. Extensive tests of both
the 2D BEM formulations and the coupled formulation are shown. Further, the new integration principle is extended to 2D and tested.
Chapter 5
The 2D BEM formulation from Chapter 4 is used with a new formulation
using Green's functions for a moving load. The results are compared to results obtained
8
1. Introduction
by a similar FEM formulation.
Chapter 6
Here the conclusions of the thesis are summarised and the overall conclusions are presented. Further, new contributions in the thesis are mentioned.
Appendices
In the appendices various numerical and analytical expressions, results
and methods are presented.
Chapter 2
Boundary Element Method
Formulations
2.1
Introduction
In this chapter the basic boundary element method formulation in the time domain is
derived. The basic formulation is further developed into the direct and the indirect BEM
formulations. For both formulations an example of a formulation for a specic problem
is shown.
The example shown for the direct formulation is the case of an alluvial valley subjected
to an incoming wave. In the indirect case the example is the formulation for an elastic
halfspace subjected to an external load on the surface of the halfspace. Furthermore,
numerical approximations in space and time are shown.
In the rest of the thesis the following Cartesian coordinate system will be used
Figure 2.1: Cartesian coordinate system.
2.2
Basic BEM Time Domain Formulation
Let Uij (x; t; ) signify the displacement eld in coordinate direction i from a concentrated
time dependent force f (t) applied at position in direction j in a homogeneous, isotropic
and linear elastic media with Lame constants and and mass density . Assuming the
media at rest at the time t = 0 the displacement is given as
Uij (x; t; ) =
Z t
0
gij (x; t; ; )f ( )d
(2.1)
9
10
2. Boundary Element Method Formulations
The upper limit t of the time integral signies that impulses up to, but not including
the time t, are considered. Since an impulse at the time t will cause discontinuities in the
velocity eld @[email protected] Uij (x; t; ), but not in the displacement eld Uij (x; t; ), this distinction
is immaterial for the present case. The limit will only be indicated where it may be of
importance. gij (x; t; ; ) is Green's function for the displacement eld, which has the
following properties, Aki & Richards [4]
1) Causality :
gij (x; t; ; )q= 0 for c1 (t q
) < jx j _ c2 (t ) > jx j
+2
where c1 = and c2 = denotes the velocities of P-
and S-waves. The statement expresses that the contribution
is null if the fastest wave has not yet reached the receiver
point x, or if the slowest wave has passed x.
2) Reciprocity : gij (x; t; ; )= gij (; ; x; t)
3) Time translation : gij (x; t; ; )= gij (x; t + t1 ; ; + t1 )= gij (x; t ; ; 0)
4) Symmetry :
gij (x; t; ; )= gji(; t; x; )
gij (x; t; ; 0) is given as, Aki & Richards [4], Banerjee et al.
(
1
3rirj
gij (x; t; ; 0) =
4 r3
ri rj 1 r 1 Æ t
Æ t
r3
c21
c1
c22
[10]
!
Æij Z c2 1
Æ (t r)d +
r c1 1
!
)
r
Æij r + rc2 Æ t c
c2
2
2
(2.2)
= jx j, ri = xi i, where xi and i signify the components of x and . Æ(t) is
the Dirac delta function. ijk (x; t; ; 0) signies Green's function for the Cauchy stress
tensor eld, i.e. the stress tensor components at positions x and the time t due to an unit
impulse applied at at the time t = 0 in the coordinate direction k. This becomes, Aki
& Richards [4], Banerjee et al. [10]
r
ijk (x; t; ; 0) = gmk;m (x; t; ; 0)Æij + (gik;j (x; t; ; 0) + gjk;i(x; t; ; 0)) =
1 ( 6c2 5 rirj rk Æij rk + Æik rj + Æjk ri ! Z c2 1 Æ(t r)d +
2
4
r5
r3
c1 1
!
!
ri rj rk Æij rk + Æik rj + Æjk ri
r c22 r 2 6 r5
Æ t c
Æ t
+
r3
c21
c1
2
!
ri rj rk _ r c32 _ r 2 r4 c Æ t c c3 Æ t c
2
2
1
1
! 2
rk Æij
c2
r
r _ r 1 2 c2 Æ t c c Æ t c
r3
1
1
1
1
)
Æik rj + Æjk ri
r
r_
r
Æ t
Æ t
3
r
c2
c2
c2
(2.3)
2.2. Basic BEM Time Domain Formulation
11
The interior and exterior of the elastic medium are denoted V + and V and the surface
is denoted S . Using the Betti reciprocal theorem and Green's functions (2.2) and (2.3)
leads to the representation theorem for any x 2 V +, Kobayashi [36]
u+(x; t) =
i
Z Z t
S
0
Z Z t
S
0
gij (x; t ; ; 0)t+j (; ; n+ ())d dS( )
ikj (x; t ; ; 0)n+k (x)u+j (; )d dS( )
(2.4)
where the body forces are assumed to be zero, and n+i () is the component of the outward
directed unit normal vector to the volume V +. Further,
the symmetry properties
gij (x; t
+
+
; ; 0) = gji (; t ; x; 0) and ikj (x; t ; ; 0)nk (x) = jki (; t ; x; 0)nk () have been
used. The surface integration for the traction (last part of (2.4)) requires special attention
for x 2 S . A cut is made on the positive side of the boundary, as shown in gure 2.2,
with the subboundaries @S;1 and @S;2 completely within V +, where @S;1 is a disk and
@S;2 is a hemisphere of radius . The last part of (2.4) can hereby be obtained as
Z Z t
S
0
S
0
Z Z t
lim
!0
Z
ikj (x; t ; ; 0)n+k (x)u+j (; )d dS( ) =
ikj (x; t ; ; 0)n+k (x)u+j (; )d dS( ) +
Z t
@S;1
0
ikj (x; t ; ; 0)n+k (x)u+j (; )d dS( )
(2.5)
where R ()dS() denotes the principal value integral in Cauchy's sense.
Figure 2.2: Boundary around singular point.
The rst part of (2.5) is an ordinary principal value integral in Cauchy's sense, but the
last part requires special attention. Since there are no concentrated forces within the
small hemisphere bounded by @S;1 and @S;2 it follows that
Z
(x; t
lim
!0 @S 1 ikj
Z
lim
(x; t
!0 @S 2 ikj
;
;
; ; 0)n+k ()dS()
; ; 0)n+k ()dS() = 0
(2.6)
12
2. Boundary Element Method Formulations
(2.6) is merely an equilibrium equation expressing that the tensions applied at the surface
of the hemisphere and the disk must be in equilibrium, since the body and inertial loads
within the hemisphere become ignorable at the limit. From (2.6) it then follows that
Z
lim
(x; t
!0 @S 1 ikj
Z
lim
(x; t
!0 @S 2 ikj
; ; 0)n+k ()dS() =
;
; ; 0)n+k ()dS()
;
(2.7)
where n+k () in both cases
signies
the normal vector in the outward direction from V +.
+
+
Along @S;1 one has nk () = nk (x). Hence
lim
!0
lim
!0
Z
where
@S;1
Z
Z t+
0
Z t
@S;2
0
ikj (x; t ; ; 0)n+k (x)u+j (; )d dS( ) =
0
ikj (x; t ; ; 0)n+k ()u+j (; )d dS( ) =
Z t
u+j (x; )sij (x; t )d
sij (x; t ) = lim
!0
Z
@S;2
(2.8)
ikj (x; t ; ; 0)n+k ()dS()
(2.9)
In (2.8) it has been used,+ that the eld u+j (; ) is a continuous function of , which is
constant and equal to uj (x; ) on the surface @S;2 as ! 0. By (2.8) the singularity
at = x in the integrand on the left hand side of, (2.8) has been omitted. Instead,
one needs to evaluate the surface integral
(2.9) along the subsurface @S;2. Along @S;2,
1
ijk (x; t ; ; 0) is of the order O 2 . Expressions in ijk (x; t ; ; 0) of the order
O 1 therefore vanish in the limit as ! 0. Using (2.3) with r = , r = n+i (), (2.9)
can be written
i
sij (x; t ) =
(
!
!
Z
n+i ()n+j () 2n+i ()n+j () + Æij
1
c22
lim
3 5 2
1 c2 Æ(t )
!0 4 @S 2
2
1
!
!
2
n+i ()n+j () 2n+i ()n+j () + Æij
c2
2 6 2
1
Æ (t ) +
2
c21
!
)
n+i ()n+j ()
n+i ()n+j () + Æij
c22
1 2 Æ(t ) +
Æ (t ) dS d =
;
c21
2
1 Z
lim
3 1
!0 42 @S 2
;
2
!
!
c22
c22 +
+
n ()nj () + Æij 2 Æ (t )dS() d
c21 i
c1
( )
(2.10)
2.2. Basic BEM Time Domain Formulation
13
where the limits
lim
Æ t
= Æ(t)
!0
c
i
1
c1 2 )
2
have been inserted. In Appendix A it is shown that [email protected] 2 n+i ()n+j ()dS() = 23 2 Æij . Use
of this result leads to
1
sij (x; t ) = Æ (t )Æij
(2.11)
2
lim
!0
Z c 1
2
c1 1
Æ (t )d = Æ (t)(c2 2
;
Hence it follows from (2.8) that
lim
!0
Z
1
ikj (x; t ; ; 0)n+k (x)u+j (; )d dS( ) = ui (x; t)
Z t
(2.12)
2
It should be noted that in this derivation it has been assumed that the boundary
is
smooth. For non-smooth points on the boundary the results will be dierent from 21 .
With the result in (2.12), (2.5) can be written as
@S;1
Z Z t
S
0
0
ikj (x; t ; ; 0)n+k (x)u+j (; )d dS( ) =
c+ (x)u+i (x; t) +
where
c+ (x) =
(
Z Z t
0
;
1 ;
2
0
S
x
x
ikj (x; t ; ; 0)n+k (x)u+j (; )d dS( )
2 V+
2 S
(2.13)
(2.14)
With the result in (2.13), (2.4) can be reformulated to
C + (x)u+(x; t) =
i
Z Z t
where
S
0
Z Z t
S
0
gij (x; t ; ; 0)t+j (; ; n+ ())d dS( )
ikj (x; t ; ; 0)n+k (x)u+j (; )d dS( )
1
C + (x) = > 0
: 1
8
>
<
;
;
2 ;
x
x
x
2 V+
2 V
2 S
(2.15)
(2.16)
(2.15) is often denoted Somigliana's identity, and is the basis for the boundary element
method, Banerjee et al. [10].
14
2. Boundary Element Method Formulations
2.3
Discretization in Time and Space
In order to solve equation (2.15) discretization in time and+ space is +necessary. In the
following the simplied notations for the surface traction tj (x; t) = tj (x; t; n+(x)) and
the contraction+ of the Green's function for the Cauchy stress tensor eld tij (x; t; ; 0) =
ikj (x; t; ; 0)nk (x) , denoted the Green's function for the traction eld, are used.
2.3.1
Space Discretization
(2.15) must be discretizised in space. Due to the complexity of the following equations
the last part of (2.15) is not included in the expressions, but the discretization is similar
to the one used for the rst part
Z Z t
C + (x)u+(x; t) =
i
Z
S
gij (x; t ; ; 0)t+j (; )d dS( ) =
0
+
gij (x; t; ; 0) t (; t)dS
S
(2.17)
( )
j
where the following designation has been introduced for the time convolution integral
( R
t
+
0
+
gij (x; t; ; 0) tj (; t) = 0 gij (x; t ;0; 0)tj (; )d ;; tt (2.18)
<0
The boundary surface integral in (2.17) is divided into M elements
u+i (x; t) =
M Z
X
m=1
Sm
gij (x; t; ; 0) t+j (; t)dS()
(2.19)
Occasionally, isoparametric elements with N node points are used. A mapping from
parent domain to global domain is dened as
i () =
N
X
n=1
N (n) ( )x(in) ;
i = 1; 2; 3
(2.20)
where x(in) is the ith component in the global domain of the node n and where the shape
functions N (n) () full
N (n) ( (m) ) = Æmn
(2.21)
(m) denotes the parent coordinates for node m.
The integration in (2.19) is performed in the parent domain and then transformed to
the global domain. The line element d in the global domain is the mapping of the
corresponding line element in the parent domain, see gure 2.3. These are related as
N
(n)
X
d = x(n) @N () d ; = 1; 2
(2.22)
n=1
@
2.3. Discretization in Time and Space
15
Figure 2.3: Mapping from parent to global domain.
where the summation convention for is suspended. (2.22) leads to the following area
transformation formula
dS() = jd1 d2j =
N X
N
X
(n)
x
n=1 m=1
x(m) @N
(n) @N (m) @1
@2
dS() = det(J())dS() (2.23)
where signies the vectorial product and J() is the Jacobian of the mapping, given as
det(J()) =
N X
N
X
(n)
x
n=1 m=1
x(m) @N
(n) @N (m) @1
@2
(2.24)
The traction is approximated by the same shape functions that are used in the coordinate
mapping
t+i ((); t) =
N
X
n=1
N (n) ( )t(in)+ (t)
(2.25)
where t(in)+ (t) = t+i (x(n) ; t) is the value of the surface traction at the node n.
Using (2.20), (2.23) and (2.25), (2.19) can be written as
C + (x)u+i (x; t) =
M X
N Z
X
m=1 n=1 S
m
gij (x; t; ( ); 0) N (m;n) ()t(jm;n)+ (t)det(J(m) ())dS( )
(2.26)
where N (m;n) (), det(J(m) ()) and t(jm;n)+ (t) signify the shape functions, the Jacobian
and the nodal values belonging to element m. The integration in (2.26) is performed by
Gauss quadrature which is described in Section 2.4.
Element types, shape functions and summary of Gauss quadratures are presented in
Appendix B.
16
2. Boundary Element Method Formulations
2.4
Numerical Integration
Evaluation of the integral in (2.26) is done by Gauss quadrature as explained in Appendix
B. In order to improve accuracy, element subdivision is used. Each element is subdivided
into a number of elements, and Gauss quadrature is used to calculate the integral of each
element. The number of subelements, i.e. the side lenghts of the subelements, depend
on the distance between the source point and the nearest point of the element. The
subdivision also depends on the order of the Gauss quadrature, which inuences the
magnitude of the error.
Manolis et al. [48] suggested a side length criterion for the subelements to be used with
3D transient elastodynamic BEM. The side length criterion uses a tolerance criterion
established by Lachat & Watson [43]. The subelement side length criterion by Manolis et
al. [48] is given as
(2Kj + 1)
Lj 2Kj
2R < cj LmaxLmin
(2.27)
where R is the minimum distance from the source point and to the element, Lj is the side
length of the subelement (j = 1; 2), Kj is the number of Gauss points, Lmax and Lmin are
the maximum and minimum side lenghts of the element and cj is a constant which must
be chosen from experience.
cj is chosen so the distance R from the centre of the element mesh to the outer elements
gives a level of subdivision that is known to produce reasonable results.
Integration of a regular quadratic element with subdivision is shown in (2.28) and an
example of subdivision of a single element is shown in gure 2.4
Figure 2.4: Element with subdivision.
Z
1 X
2 X
K1 X
K2
1Z1
1 X
f ()d wi;j (k)wi;j (l)f (~1 ; ~2 )
1 2 i=1 j =1 k=1 l=1
1 1
where i is the number of subelements in direction i and ~1, ~2 are given as
~1 = 1 + 11 (1 2i + 1(i;j ) (rk ; rl ))
~2 = 1 + 12 (1 2j + 2(i;j ) (rk ; rl ))
(2.28)
(2.29)
2.5. Time Discretization
17
where 1(i;j) and 2(i;j) are the local coordinates in subelement i; j (see gure 2.4) and rk
are the Gauss points.
In the case of triangular elements, a similar subdivision is used, but the element is divided
into a number of quadratic elements and some triangular elements.
2.4.1
Singularity for Green's Functions
When the source and the reciever x points coincide, the Green function for both the
displacement eld and the traction eld is singular. Numerical integration in this case
requires heavy calculations and accurate results are diÆcult to achieve. In Chapters 3
and 4 two dierent approaches are developed.
2.5
Time Discretization
(2.15) is reformulated as
C + (x)u+(x; t) =
Z Z t
S
0
i
Z Z t
gij (x; ; ; 0)t+j (; t )d dS( )
0
+
tij (x; ; ; 0)u (; t
S
)d dS( )
j
(2.30)
The time interval [0; t] is divided into L subintervals, each of the length t = Lt . Thereby
(2.30) can be written as
L Z Z tl
X
C + (x)u+i (x; tL ) =
l=1
L Z Z tl
X
l=1 S tl 1
S tl 1
gij (x; ; ; 0)t+j (; tL
tij (x; ; ; 0)u+j (; tL
)d dS( )
)d dS( )
(2.31)
where tL = t and t0 = 0. A linear variation of the eld variables is assumed in each time
step, and (2.31) becomes
L Z Z tl
X
C + (x)u+i (x; tL ) =
l=1
i
t+ (; tL tl )l d dS
j
i
S tl 1
L Z Z tl
X
( )
tl )l d dS( )
u+j (; tL
where l = t t 1 .
l
h
gij (x; ; ; 0) t+j (; tL
l=1 S tl 1
tl 1 )(1 l ) +
h
tij (x; ; ; 0) u+j (; tL
tl 1 )(1 l ) +
(2.32)
18
2. Boundary Element Method Formulations
Figure 2.5: Shape functions for time step l.
2.6
Direct Boundary Element Method
The initial conditions ui(; 0) = 0 are used. After introducing the space discretization
(2.26), (2.32) can be represented in the following matrix format
1 Iu+ = G t+ T u+ + g+
(2.33)
2 L 1;1 L 1;1 L L
+=
gL
LX1
(GL+1 l;1 + GL l;2) t+l
l=1
LX1
l=1
(TL+1 l;1 + TL l;2) u+l
(2.34)
where uL = u(tL) and tL = t(tL) are assembly of all displacements and tractions in the
boundary element mesh. I is an identity matrix and the boundary element matrices Gl;1,
Gl;2 , Tl;1 and Tl;2 are dened by
=
Gl;1
=
Gl;2
Tl;1
Tl;2
=
=
"Z
Z tl
Sm tl 1
"Z
Z t
l
Sm tl 1
"Z
Z t
l
Sm tl 1
"Z
Z t
l
Sm tl 1
gij (x; ; n ; 0)(1 l )d dS( )
gij (x; ; n ; 0)l d dS( )
(2.35)
#
tij (x; ; n ; 0)(1 l )d dS( )
tij (x; ; n ; 0)l d dS( )
#
(2.36)
#
(2.37)
#
(2.38)
The brackets symbolize ordering in matrix form of the components relating contributions
to the nth node from the mth element, during the lth time interval. The time integrals
in (2.35)-(2.38) may all be evaluated analytically by the use of the following results
Z c 1
1
c2 1
Z tl
tl 1
t Æ (t r)d = 2 H t
r
tl
t
H r
c
8
>
>
>
<
d = >>
>
:
r
H t
c1
1 t3
6t l
1 3
6t tl
0
r c2
(2.39)
3tl t2l 1+ 2t3l 1 3tl rc 2 + 2 rc 3
; tl 1 rc
; tl rc tl 1
; rc tl
(2.40)
2.6. Direct Boundary Element Method
Z tl
tl 1
tl 1 t H 19
r
d =
c
1 2t3 3t2 tl 1 + t3 l
l 1
6t l
1 2t3 3t2 tl 1 + 3tl 1 r 2
l
6t l
c
8
>
>
>
<
>
>
>
:
2
0
where H (t) is the Heaviside unit step function.
Z t
l
tl 1
Z t
l
tl 1
Z tl
tl 1
tl Æ
t
tl 1
t Æ
r
c
r
c
_ tl Æ
t
r
c
1 tl
t
r
c
(
1 r
t c
tl 1
d = 0
(
d = 0
(
1
d = 0t
3 r
c
; tl 1 rc
; tl rc tl 1 (2.41)
; rc tl
;
;
tl 1 rc tl
r t _r t
l c
l 1
c
(2.42)
;
;
tl 1 rc tl
r t _r t
l c
l 1
c
(2.43)
;
;
tl 1 rc tl
r t _r t
l c
l 1
c
(2.44)
(
1
;
tl 1 rc tl
d
= 0 t
(2.45)
r t _r t
;
t 1
l c
l 1
c
where R 11 (t)Æ_ (t t0 )dt = _ (t0) have been used. If boundary conditions for a specic
problem are introduced then (2.33) can be solved. On each part of the boundary either
displacements or surface tractions much be specied. With boundary conditions applied
to (2.33) the +equations can be rearranged
and solved with respect to the unknown displacements uL and surface tractions t+L .
Z tl
l
_
tl 1
t Æ
Example 2.1
r
c
Direct BEM Formulation for Alluvial Valley
Figure 2.6: Alluvial valley.
20
2. Boundary Element Method Formulations
In this example a sedimentary basin also called an alluvial valley is subjected to an earthquake some
distance away. Figure 2.6 shows an alluvial valley V + submerged into an innite medium V . Both
volumes are linear elastic, homogeneous, isotropic and are dened by the Lame constants + ; + and
; and the mass densities + and . S1 and S2 signify the free surfaces of the valley and the elastic
media, whereas the interface has been denoted S12 . S0 signies an articial boundary of the elastic
medium at innity. The system has been loaded by a stress wave eld at innity with particle motion
u0(x; t), which may contain components
from both P-, SV- and SH-waves. The far eld is decomposed
into a refracted component u+ (x; t) in the alluvial valley and a diracted component ud (x; t), which
will dissappear on the articial boundary. It is assumed, that the refracted eld u+ (x; t) does not cause
non-linear responses in the alluvial valley.
Let u (t) and t (t) be an assembly of nodal point displacements and discretized surface tractions related
to the nodes on the surface S2 [ S12 of the medium V . These are given by
u (t) = u (t) + ud(t) (2.46)
t (t) = t (t) + td(t)
where u (t) and ud (t) are the nodal point displacements from the far eld and the continuous eld ud (x; t)
on S [ S , and t (t) and td (t) are the corresponding surface tractions derived from these elds.
Next the eld ud (x; t) is considered. As indicated in gure 2.6, this will disappear at innity on the
0
0
0
2
12
0
articial boundary S0 . Extending the surface integrals merely over S2 [ S12 , the following equations for
the diracted eld at the time tL are obtained from the direct BEM, cf. (2.33) and (2.34).
1
I + T1;1
2
ud;L = G ; td;L + gL
LX1 gL =
l=1
(2.47)
11
GL
l;1 +
+1
GL
l;2
td;l
LX1 l=1
TL
T
l;1 + L l;2
+1
ud;l
(2.48)
where GL; has been dened in (2.35) and (2.36) and TL; in (2.37) and (2.38). ud;l and td;l signify the
displacements and surface tractions at the nodes on S2 [ S12 at the time tl . These are assumed to be
known for l = 0; 1; ; L 1, so ud;L and td;L are the unknows of (2.47). ud and td are eliminated from
(2.47) and (2.48) in favour of u (t) and t (t) by means of (2.46), and the resulting equations are solved
for uL . The results read
uL = D tL + fL
D
(2.49)
1
=
I + T1;1
2
1
0
LX1 l=1
TL
+1
(2.50)
11
fL = u ;L D t ;L +
0
G;
1
I + T1;1
2
T
u
l;1 + L l;2 ( l
1
!
LX1 l=1
GL
+1
l;1 + GL l;2 (tl
u ;l)
0
t ;l )
0
(2.51)
In a similar manner a relationship for u+ (t) and t+ (t) can be established as
uL = D tL + fL
+
D
+
+ +
(2.52)
+
1
=
I + T+1;1
2
1
G;
+
11
(2.53)
2.6. Direct Boundary Element Method
fL = 12 I + T+1;1
LX1 TL
l=1
LX1 1
l=1
T
l;1 + L l;2
+
+1
GL
ul
+
l;1 + GL l;2
+
+1
21
tl
+
!
(2.54)
+
(2.49) and (2.52) are partitioned into parts representing interior and boundary nodes, designated supscripts i and b, respectively
ui;L
ub;L
ui;L
ub;L
+
+
=
=
Dii Dib
Dbi Dbb
ti;L
tb;L
Dii Dib
Dbi Dbb
ti;L
tb;L
+
+
+
+
+
+
+
+
fi;L
fb;L
fi;L
fb;L
+
(2.55)
(2.56)
+
On the boundary S12 the following conditions are required
ub;L = ub;L
(2.57)
tb;L = tb;L
(2.58)
+
+
For the system shown in gure 2.6 the surface has been assumed to be unloaded, so ti;L = t+
i;L = 0. In
what follows a formulation with prescribed surface tractions ti;L and t+
,
which
may
be
dierent
from
i;L
+
+
zero, is worked out. Due to (2.57) and (2.58), the unknowns of the system are ui;L , ui;L , ub;L and t+
b;L .
The second equations of (2.55) and (2.56) are solved for tb;L and t+
respectively.
Using
(2.57)
and
b;L
(2.58), the following results are obtained
tb;L = (Dbb )
+
1
tb;L = (Dbb )
+
+
Dbi ti;L + fb;L ub;L
1
+
Dbi ti;L + fb;L ub;L
+ +
+
+
(2.59)
(2.60)
(2.59) and (2.60) the provide
ub;L = (Dbb )
+
+
1
+ (Dbb )
1
1
(Dbb )
1
Dbi ti;L + (Dbb ) Dbi ti;L + fb;L + fb;L
+
1
+ +
+
(2.61)
All quantities on the right hand side of (2.61) are known at the time tL , so that u+
b;L can be calculated.
+
+
Next, tb;L can be calculated from (2.59). Finally, ui;L and ui;L can be calculated from the rst equations
of (2.55) and (2.56).
22
2. Boundary Element Method Formulations
2.7
Indirect Boundary Element Methods
Alternatively to Section 2.6 the integral equations can be formulated from a single-layer
potential theory, as explained below, Kobayashi [36]. The domain is divided into ctitious
interior and exterior domains, see gure 2.2.
If both the interior volume V + and the exterior volume V are lled with the same elastic
material, then the same Green functions apply to both volumes. But the boundary conditions for V along S , symbolised by ui (; ) and ti (; ; n ()) = ti (; ; n+()) =
ti (; ; n+ ()) may not be the same as for V + . From the assumptions, (2.15) applies
for the volume V with the same Green functions
C
(x)ui (x; t) =
Z Z t
where
S
C
Z Z t
0
S
0
gij (x; t ; ; 0)tj (; ; n+ ())d dS( ) +
ikj (x; t ; ; 0)n+k (x)uj (; )d dS( )
0
(x) = >: 11
8
>
<
;
;
2 ;
x
x
x
(2.62)
2 V+
2 V
2 S
(2.63)
Adding (2.15) and (2.62) leads to the displacement eld for the entire volume V + [ V [ S
C + (x)u+i (x; t) + C
Z Z t
S
0
S
0
Z Z t
gij (x; t
(x)ui (x; t) =
h
; ; 0) t+j (; ; n+ ())
i
tj (; ; n+ ()) d dS( )
h
i
ikj (x; t ; ; 0)n+k (x) u+j (; ) uj (; ) d dS( )
(2.64)
In what follows a socalled single-layer potential formulation is+applied, i.e. on the boundary S it is required that the displacements are continuous, ui (; ) = ui (+; ), but the
tractions may not be equal, requiring single-layer potential i(; ) satisfying ti (; ; n+()) =
ti (; ; n+ ()) + i (; ) to be introduced. Using these conditions, (2.64) can be rearranged, giving the displacement u+i (; ) as an integral equation of the densities i(; )
ui (x; t) = u+i (x; t) =
Z Z t
S
0
gij (x; t ; ; 0)j (; )d dS( )
(2.65)
Since gij (x; t ; ; 0) is of the order O 1r and j (; ) is continuous, the integral (2.65)
exists even for x 2 S . For x 2 V + the surface traction follows directly from (2.65) by
dierentiation as an integral equation of i(; ). Performing the dierentiation of (2.65)
by using the result in (2.13), the following is obtained
t+i (x; t; n+ (x)) = c+ (x)i (x; t) +
Z Z t
S
0
ikj (x; t ; ; 0)n+k (x)j (; )d dS( ) (2.66)
2.7. Indirect Boundary Element Methods
where
c+ (x) =
(
01
;
2 ;
23
2 V+
2 S
x
x
(2.67)
(2.65) and (2.66) are often referred to as the single-layer potential integral equations.
The surface traction ti (x; t; n+(x)) on +the exterior volume V can now be established
from the relationship ti (x; t; n+(x)) = ti (x; t; n+(x)) i(x; t) which gives
ti (x; t; n+ (x)) = c
where
(
c
(x) = 01
;
2 ;
x
x
Z Z t
(x)i(x; t)+
S
0
ikj (x; t ; ; 0)n+k (x)j (; )d dS( ) (2.68)
2 V
2 S
(2.69)
Due to the indirect determination of the physical quantities the formulation by potentials
is called the indirect BEM formulation.
Example 2.2
Single-Layer Potential Formulation for Elastic Halfspace
Figure 2.7: Elastic halfspace with coordinate system.
In this example an isotropic, homogenous, linear elastic halfspace is subjected to either precribed tractions,
displacements or both on the surface. S = S1 [ S2 signies the surface of the halfspace, see gure 2.7.
0
Assume that surface tractions t+
i (x; t) = ti (x; t) are prescribed on a part S1 of the surface and the
+
0
displacements ui (x; t) = ui (x; t) on the remaining part S2 = S=S1 . Then i (x; t) throughout S = S1 [ S2
can be found from the following integral equations provided by (2.65) and (2.66)
Z
Z t
S2
0
x
gij ( ; t
Z
1
(x; t) +
2 i
S1
; ; 0)j ( ; )d dS( )
Z t
0
x
ikj ( ; t
= u0i (x; t)
x2S
;
x)j (; )d dS ; ; 0)n+
k(
(2.70)
2
( )
= t0i (x; t) ;
x2S
1
(2.71)
After solving (2.70) and (2.71) with respect to i (x; t) the displacement can be determined from (2.65)
in V + [ S1 and the surface tractions from (2.66) throughout V + [ S2 .
Formulation of (2.70) and (2.71) with matrix notation as for (2.33) leads to
uS2;L = G ; S2;L + gS2;L
+
11
+
(2.72)
24
2. Boundary Element Method Formulations
tS1;L =
+
1
I + T1;1
2
S1;L + hS1;L
+
(2.73)
where
gS2;L =
+
hS1;L =
+
2.8
LX1
l=1
LX1
l=1
(GL+1 l;1 + GL l;2 ) S2 ;l
(2.74)
(TL+1 l;1 + TL l;2 ) S1 ;l
(2.75)
Concluding Remarks
In the present chapter the basic equations for the BEM are derived. It is shown how the
Betti reciprocal theorem with Green functions is transformed into the Somigliana identity.
The basic principles for discretization of Somigliana's identity in time and space are shown
and the principle of integration by subdivision is mentioned. The Somigliana identity is
directly discretized to obtain the direct BEM matrix equations and it is shown for a
linear time variation discretization. Further, an example of a direct BEM formulation for
the specic problem of an alluvial valley subjected to an incoming earthquake is shown.
Similarly the derivation of the indirect BEM from the Somigliana identity is presented.
The indirect BEM formulation is shown with the example of an elastic halfspace subjected
to tractions or prescribed displacements on the surface.
In the literature, simultaneous derivations of the direct and indirect BEM are seldom seen
and therefore their mutual background is often missed. As shown in this chapter, both
methods represent ways of discretisizing the Somigliana identity, and the same numerical
methods applies for both formulations.
Chapter 3
BEM Formulation in 3D
3.1
Introduction
The indirect BEM formulation described in Chapter 2 is used to model wave propagation
in an isotropic homogeneous halfspace in 3D. The applied circular plane element mesh on
the free surface is shown in gure 3.1.
Figure 3.1: Circular plane element mesh for halfspace.
The discretization from Section 2.3 is modied due to the applied circular plane mesh. An
analytical solution of the singularity problem is presented and the numerical integration
scheme is outlined. An analytical solution for Lamb's problem is compared with numerical
results from the boundary element program. Finally, a substantial improvement of the
method is discussed and implemented.
3.2
Approximation of Single-Layer Potential
The division of the boundary (2.19) into M elements is still valid. However, the surface
integrals are now evaluated in polar coordinates, i.e.
u+(x; t) =
i
M Z
X
m=1 S m
gij (x; t; r; ; 0) j (r; ; t)rdrd
25
(3.1)
26
3. BEM Formulation in 3D
The single-layer potential is approximated by shape functions N (n) (r; ) in the following
way
i (r; ; t) =
N
X
n=1
N (n) (r; )(in) (t)
(3.2)
where (in) (t) = i(r; ; t) is the value of the single-layer potential at node n. The shape
functions N (n) (r; ) full
N (n) (rm ; m ) = Æmn
(3.3)
where rm; m denotes the polar coordinates for node m. With (3.2), (3.1) can be written
as
u+i (x; t) =
M X
N Z
X
m=1 n=1
Sm
gij (x; t; r; ; 0) N (m;n) (r; )(jm;n) (t)rdrd
(3.4)
where N (m;n) (r; ) and (jm;n) (t) signify the shape functions and the nodal values belonging
to element m. Shape functions in cylindrical coordinates have been derived and can be
seen in Appendix B.
3.2.1
Singularity Problems of Green's Function for the Displacement Field
When source and reciever x points coincide, the Green's function for both the displacement eld and the traction eld is singular. For this case a special analytical solution is
devised in order to avoid troublesome numerical quadratures.
The analytical solution is obtained by making a circular integration in polar coordinates
around the singular node. If no additional error is to be introduced the size of the time
steps must be so small that the wave front only moves the distance equal to the smallest
length from any node and to the boundaries of the surrounding elements. This requirement also implies that the singularity problem only arises in the rst time step. Hence
the following bounding at the time step must be fullled.
t Rcmin
1
; R = tc1
(3.5)
where Rmin is the minimum distance from any node and to the boundaries of the surrounding elements and R is the radius where the analytical integration is performed.
3.2. Approximation of Single-Layer Potential
3.2.2
27
Analytical Solution of Singularity Problem
Green's function for the displacement eld can be written as, cf. (2.2) and (2.39)
gij (x; t; ; 0) = gijÆ + gijH
(3.6)
where
(
!
!)
1
3
ri rj Æij
r
r
= 4 r5 r3 t H t c H t c
1
2
(
)
!
1
ri rj 1
r
1
r
Æij r Æ
gij =
(3.7)
4 r3 c21 Æ t c1 c22 Æ t c2 + rc22 Æ t c2
With the requirement in (3.5), the quadrature that must be solved is, cf. (2.33)-(2.38)
Z Z t t gij (x; ; ; 0)
(3.8)
t + gij (x; ; ; 0) t d dS()
S 0
where SR is the circular area with radius R. The following integrals become useful
gijH
R
t
Z
Z
0
r
c
2
r
c
1
Z
t
r
r H 1
d =
c1
c2 !
t
!
r2 1 1
r3 1 1
r
d
=
;
2
2
3
3
t
2 c2 c1 3t c2 c1 c2
H 1
r 1
c
1
r
r
d
=1
; t
t
ct c
0
Besides, the following nominations will be used
Æ t ; cr t (3.9)
(3.10)
t ; gÆ2 = gÆ ij
ij t
t
t ; gH 2 = gH gijH 1 = gijH
ij
ij t
t
(3.8) and (3.10) are combined into the following result for gijÆ1
gijÆ1 = gijÆ
Z
Z
SR
0
(3.11)
t Æ1
gij d dS =
(3.12)
ri rj 1 1 r 1 1 r
Æij 1 r
t c c2 t t c + rc2 t t c dS =
r3 c21 t
SR 4
1
2
2
2
2
(
!Z
!Z
)
Z
1 1 1
1
1
Æij
Æij Z
ri rj
1
ri rj
1
4 c21 c22 SR r3 dS t c31 c32 SR r2 dS + c22 SR r dS c32 t SR dS
For gijÆ2, gijH 1 and gijH 2 the following expressions have been found in a similar manner
Z
1
(
!
)
28
3. BEM Formulation in 3D
Z
Z
0
SR
t Æ2
gij d dS =
1 ( 1
4t c31
Z
t
R
gijH 1
2
Z
t
Z
SR
0
gijH 2
1
c32
(
ri rj
2 dS
R r
R
)
ri rj
Æij Z
dS + c3 S dS
SR r 2
2 R
! Z
1
1
d dS = 24t 2 c3 c3 3 S
S 0
1
2
! Z
)
Z
1
ri rj
1
1
+3t c2 c2 3 S r3 dS Æij S r dS
Z
1
1 !Z
Æij
Z
SR
(3.13)
dS
(3.14)
R
1
d dS = 12t
1
c32
1 ! 3 Z
c31
ri rj
dS
SR r 2
Æij
Z
SR
dS
(3.15)
From (3.12)-(3.15) it can be seen that the following four integrals must be solved
2
3
1
0
0
Z
ri rj
d
S = R 64 0 0 0 75
(3.16)
3
S r
0 0 1
2
3
1
0
0
Z
ri rj
d
S = R2 64 0 0 0 75
(3.17)
2
2 0 0 1
S r
Z
1 dS = 2R
(3.18)
S r
R
R
R
Z
SR
dS = R2
(3.8) then follows by combination of (3.12)-(3.19).
(3.19)
3.2. Approximation of Single-Layer Potential
3.2.3
29
Singularity Problem of Green's Function for the Traction
Field
Figure 3.2: a) Green's function for the traction eld for component 1,2,1. b) Section
through x2 x3 plane.
Integration near the singularity for Green's function for the traction eld is very diÆcult
to handle both numerically and analytically. But when the area of interest is a plane
homogeneous halfspace, as in this case, it is not necessary to make any calculations to
nd a reasonable result. Sanchez-Sesma & Luzon [69] showed that Green's function for
the traction eld on a plane surface posseses planes of antisymmetry. As an example,
the Green's function for the traction eld around the singular point for the component
121 (x; t; ; 0) is shown in gure 3.2. As seen the function is antisymmetrical around the
x1 x2 planes.
When the integration over an element mesh is performed for a specic location of the
source point, the nal result will be equal to zero if the condition (3.5) is met.
Then, the surface integral vanishes in equation (2.13) for the singular nodes, i.e. the nodes
where source and reciever point are close, which again means that for the rst time step
ti (x0 ; t; n+ (x)) = c(x)i (x; t)
(3.20)
Example 3.1
Results for 3D Problem Subjected to a Concentrated Heaviside Load
The 3D BEM formulation described will be tested using a classical wave progration problem. An isotropic,
homogeneous, linear elastic 3D halfspace is subjected to a concentrated downward directed unit load H (t)
normal to the surface z = 0 at the origin of the cylindrical (r; ; z ) coordinate system at the time t = 0.
The z -axis coincides with the x2 -axis. This problem was rst introduced by Lamb [44], and the solution
of the displacement eld was obtained for a Poisson material by Pekeris [63], see Appendix C. Based on
this solution, the result for a uniformly distributed load over a circular area can be obtained by numerical
integration see Appendix C.2.3, which can be compared with the numerical solutions in a meaningful
sense.
The material parameters are chosen as = = 3:5 108 N/m2 and = 2000 kg/m3 which corresponds
to c1 = 725 m/s and c2 = 418 m/s. A circular symmetric mesh with 4 concentric circles is used. The
radius of the 4 circles increases from 0.05 m to 0.045 m from the centre and outwards. An applied load
of magnitude 1 MN is distributed over the inner circular area.
In gures 3.3 - 3.5 the results from the BEM analysis are shown for 8, 12 and 16 elements, respectively,
in each circular area. Numerical tests with dierent widths and number of concentric circles in the radius
30
3. BEM Formulation in 3D
direction have been performed. The tests showed that increasing element sizes from the centre of the
mesh, where the load is applied, and outwards could be used without loss of accuracy. Therefore, the
element meshes for which the results are presented are made with concentric circles with diameters of
0.05 m, 0.15 m, 0.25 m and 0.45 m, respectively. Further tests have showed that 4 concentric circles are
suÆcient.
5
15πµruz
0
−5
8 elements
12 elements
16 elements
−10
−15
0
0.5
1
1.5
2
2.5
c t/r
Figure 3.3: Vertical displacement uz at (r; ; z)=(0.90 m, 0 m, 0 m).
2
5
15πµrur
0
−5
8 elements
12 elements
16 elements
−10
−15
0
0.5
1
1.5
2
2.5
c2t/r
Figure 3.4: Horizontal displacement ur at (r; ; z)=(0.90 m, 0 m, 0 m).
−14
15πµruθ
10
x 10
8 elements
12 elements
16 elements
5
0
−5
0
0.5
1
1.5
2
2.5
c2t/r
Figure 3.5: Horizontal displacement perpendicular to the wave front u at (r; ; z)=(0.90 m, 0 m, 0 m).
It is apparent from gure 3.5 that the displacements perpendicular to the wave front u are negligible
compared to uz and ur , which was expected.
3.3. A new Integration Scheme for Surface Problems
31
In gures 3.3 and 3.4 the results obtained with 8 elements in each concentric ring dier from the other
results, especially in the later part of the plots where the displacements are fully developed. The results
with 12 and 16 elements in each concentric ring are in better agreement with each other. But all the
results show signs of instability and after the development of absolutely maximum displacements the
exitation does not calm down as expected. In gures 3.6 and 3.7 below the results from gures 3.3 and
3.4 with 12 elements in each concentric circle are compared with analytical results.
5
Analytical result
12 elements
15πµruz
0
−5
−10
−15
0
0.5
1
1.5
2
2.5
c2t/r
Figure 3.6: Vertical displacement uz at (r; ; z)=(0.90 m, 0 m, 0 m).
5
15πµrur
0
−5
−10
−15
0
Analytical result
12 elements
0.5
1
1.5
2
2.5
c2t/r
Figure 3.7: Horizontal displacement ur at (r; ; z)=(0.90 m, 0 m, 0 m).
It is clear from gures 3.6 and 3.7 that the numerical results agree qualitatively with the analytical results
up to approximately c2 t=r=1.2, where the absolutely maximum displacements are reached. From this
point the numerical results show tendencies of instability. Further analysis of the results and tests of the
integration procedures applied has shown that the problem can be solved by improving the numerical
integration. The improved numerical integration scheme will be introduced in the next section.
3.3
A new Integration Scheme for Surface Problems
In Chapter 2 and in Section 3.2 integration schemes based on the ordinary boundary
element descriptions are applied. As shown previously, element subdivision was used
in order to achieve reasonably accurate results. Evaluation of the integral equations in
the dynamic BEM provides diÆculties due to two properties of Green's functions. The
rst problem is the singularity when ! x. The second problem comes from numerical
problems with the spatial integrations in (3.4). The rst problem has been addressed in
32
3. BEM Formulation in 3D
Sections 3.2.1 and 3.2.3 where an analytical approach was introduced. Whereas the rst
problem arises in both the static and the dynamic BEM, the second problem only appears
in the dynamic BEM. Green's functions after the analytical time integration (see Section
2.3) are discontinuous in space as seen from (2.40)-(2.45). This discontinuity is diÆcult to
handle with suÆcient accuracy with ordinary numerical element by element integration.
The new integration scheme is a way of tracking the discontinuities better and thereby
achieving higher accuracy at fewer element subdivisions.
According to the formulation (3.4) the basic quadrature previously described is to be
solved. Normally, x is seen as the reciever point and ( expressed in cylindrical coordinates r; ) is seen as the source point. However, since the important factor for Green's
functions is the length of the vector between the two points, r = x we can choose
which one shall be the source or reciever. Therefore, if is chosen as the reciever and x
as the source, the result is that for a given position x at a given time t, the integration
in (3.4) is an integration of a wave front emerging from x. In gure 3.8 an illustration of
the problem is shown. A wave front emerging from x is shown on an element mesh.
Figure 3.8: A wave front emerging from source point x.
Because of the mentioned duality the hatched areas in gure 3.8 can also be seen as the
source areas contributing to the receiver point x at various instants of time. Outside
these areas the contribution to the surface integral in (3.4) is 0. Hence the integration
should concentrate on these hatched areas. The ordinary way of solving (3.4) as described
previously involves subdividing the element and then using numerical integration. What
must be realized is that the shape functions N (m;n) (r; ) are continuous functions over the
element, however, this does not apply to the integrand in (3.4). A close look at (2.40)(2.45) reveals that the integrand (the hatched ring in gure 3.8) has 4 discontinuous
circular lines. It is zero up till r = (l 1)tc2 and after r = ltc1 and has discontinuities
at r = ltc2 and r = (l 1)tc1 (l=current time step). A section of the integrand
in radial direction is shown in gure 3.9. These discontinuities give poor results for the
ordinary numerical integration over a single element, which can easily be seen if the wave
front is drawn for dierent time steps on a single element, see gure 3.10.
3.3. A new Integration Scheme for Surface Problems
Figure 3.9: Section of the wave front over which the integration is performed.
(l 1)tc2 .
33
rref
=
Figure 3.10: The wave front on a single element x.
It is important to remember that the spatial integration is performed on an element
basis. This means that for all the time steps there will be discontinuities in the area of
integration, and for the main part of the time steps the area of interest (the wave front)
will only ll a small part of the element.
In order to improve the accuracy of the integration a wave front based integration is
introduced instead of the approach used in Example 3.1. Since the discontinuities of the
wave front are well dened it is easy to discretize the wave front into elements and perform
the numerical integration on this basis. In gure 3.11 the principle is shown for a single
element.
34
3. BEM Formulation in 3D
Figure 3.11: Wave front based integration.
As seen in gure 3.11 the integration is performed in a local cylindrical coordinate system
with the source point x as the origin. It is clear that there will be discontinuities within
the new integration area as was the case within the old area. However, contrary to the old
approach, the new formulation will imply that the main part of the integration area will
consist of the area between the wave fronts. Furthermore, the two discontinuities within
the wave front will not provide any errors since the discretization/subdivision of the wave
front will follow the discontinuity lines.
In the new formulation (3.4) attains the form
u+i (x; t) =
M X
N Z
X
~
gij (x; t; r~; ; 0) N
m
m=1 n=1 S~
(m;n) (~r; ~)(m;n) (t)~rd~rd~
j
(3.21)
where (~r; ~) are the cylindrical coordinates for a coordinate system with origin at x, and
S~m is the part of the wave front which is within element m as illustrated in gure 3.11.
Example 3.2
Results with new Integration Scheme
Example 3.1 is repeated using the new integration scheme, and the improvement of the results in comparison to the previous results is demonstrated.
In gures 3.12 - 3.14 the results from the BEM analysis with the new integration scheme are shown for
8, 12 and 16 elements, respectively, in each concentric ring.
3.3. A new Integration Scheme for Surface Problems
35
5
15πµruz
0
−5
8 elements
12 elements
16 elements
−10
−15
0
0.5
1
1.5
2
2.5
c2t/r
Figure 3.12: Vertical displacement uz at (r; ; z)=(0.90 m, 0 m, 0 m) with the new integration scheme.
5
15πµrur
0
−5
8 elements
12 elements
16 elements
−10
−15
0
Figure 3.13:
scheme.
0.5
1
1.5
2
2.5
c t/r
2
Horizontal displacement ur at (r; ; z )=(0.90 m, 0 m, 0 m) with the new integration
−14
1
x 10
8 elements
12 elements
16 elements
15πµruθ
0
−1
−2
−3
0
Figure 3.14:
0.5
1
1.5
2
2.5
c t/r
2
Horizontal displacement perpendicular to the wave front u at (r; ; z )=(0.90 m, 0 m, 0
m) with the new integration scheme.
The displacements perpendicular to the wave front u in gure 3.14 are, as in gure 3.5, negligible
compared to uz and ur . Comparison of gures 3.5 and 3.14 reveals that the absolute maximum u
obtained by the new integration scheme is a factor 3 smaller than the absolute maximum u obtained by
the old integration scheme.
The results for uz and ur in gures 3.12 and 3.13 for 8, 12 and 16 elements, respectively, in each concentric
ring show good agreement with each other. The eect of using more and smaller elements is signicantly
reduced by the new integration scheme, indicating that the new scheme is much more stable for larger
elements. The instability of the results after development of the absolutely maximum displacements has
36
3. BEM Formulation in 3D
now completely disappeared.
In gures 3.15 and 3.16 below, the results from gures 3.12 and 3.13 with 12 elements in each concentric
ring are compared with the analytical results and the previous results from gures 3.3 and 3.4 with 12
elements in each concentric ring.
5
Analytical result
New integration
Old integration
15πµruz
0
−5
−10
−15
0
0.5
1
1.5
2
2.5
c t/r
Figure 3.15: Vertical displacement uz at (r; ; z)=(0.90 m, 0 m, 0 m).
2
5
15πµrur
0
−5
Analytical result
New integration
Old integration
−10
−15
0
0.5
1
1.5
2
2.5
c2t/r
Figure 3.16: Horizontal displacement ur at (r; ; z)=(0.90 m, 0 m, 0 m).
As seen from gures 3.15 and 3.16, the new integration scheme provides a fair agreement with the
analytical results both before and after the development of the absolutely maximum displacements. It
should be noted that the analytical results for a distributed load represent a combination of analytical
and numerical results, see Appendix C. However, the indicated result labelled as analytical is suÆciently
accurate to be representative.
3.4
Concluding Remarks
In this chapter a BEM formulation for a specic stress wave propagation problem is addressed. First it is shown how the BEM in 3D can be reformulated in polar coordinates,
making it possible to construct an analytical solution to the singularity problem. Numerical results are presented with the derived formulation. The results show good agreement
but also signs of instability. Instability is found to arise due to inaccuracy of the numerical
integration. Therefore an alternative integration scheme has been suggested.
The new integration scheme, which is based on integration by wave fronts, is explained
and the example is tested using this method. The results show signicant improvements
3.4. Concluding Remarks
37
and the instabilities have disappeared.
The conclusion of this chapter is that numerical wave front based integration is a possibility of ensuring accurate results for 3D BEM formulations in the time domain. In the
next chapters the new integration scheme will be applied to 2D problems.
38
3. BEM Formulation in 3D
Chapter 4
Coupled BEM and FEM Solutions
4.1
Introduction
The boundary element method is superior to the FEM in solving elastodynamic problems
in innite domains due to its inherent ability to satisfy the radiation conditions exactly.
In contrast, any impulse propagating in an FEM mesh will be reected in the solution
domain, when the impulse reaches the boundary of the mesh. However, a class of special
absorbing FEM elements and some special absorbing boundary conditions exists. But
common for all these are that they only give accurate results under special conditions and
they all represent numerical approximations to the radiation conditions, see e.g. Krenk
et al. [39], Krenk & Kirkegaard [40], Takemiya & Kellezi [75] and Takemiya et al. [76].
On the other hand, the FEM is more favourable for modelling conned domains with
non-homogeneous, non-isotropic or non-linear material properties. A large class of problems exists, such as soil-structure and uid-structure interaction problems, where a conned linear or non-linear structural system is interacting with an unlimited surrounding
continuum with linear homogeneous mechanical properties. In such cases it would be
advantageous, if the best of the BEM and the best of the FEM could be combined, so the
surrounding continuum is modelled using the BEM and the structure using the FEM. In
this chapter the coupling of such methods is investigated.
Dynamic BEM and FEM coupling schemes have previously been applied to soil-structure
interaction problems by Kobayashi & Mori [38], Kobayashi & Kawakami [37] and Mita
& Luco [51]. These formulations were performed in the frequency domain, which is appropriate, since both the soil and the structure were assumed to be linear. Formulations
of the same problem in the time domain have been performed by Fukuri [21], Fukuri &
Ishida [22], Spyrakos & Beskos [73], Spyrakos et al. [74], Karabalis [33], Karabalis &
Beskos [34], Karabalis & Gaitanaros [35], Von Estor [78, 79], Von Estor & Kausel [81],
Von Estor & Prabucki [82], Von Estor & Antes [80], Belytschko & Lu [12] and Abe &
Yoshida [1]. Time domain formulations become mandatory, whenever non-linearities are
present in the structural system, since frequency response analysis is prohibited in this
case.
In this chapter the BEM/FEM coupling scheme will be specied for systems as those
shown in gure 4.1.
Figure 4.1a shows the problem from Example 2.1. However, the alluvial valley V + consist
39
40
4. Coupled BEM and FEM Solutions
Figure 4.1: Examples of BEM/FEM coupling. a) Non-linear alluvial valley interacting
with linear elastic media for an incoming stress wave eld. b) Same as a) but with impulse
load applied within the medium. c) Interaction between linear elastic soil with incoming
stress wave eld and structure with impulse loads.
4.2. Coupling of BEM and FEM Formulations
41
in this case not of a linear elastic material but of a material with non linear properties.
The same notation as used in Example 2.1 is adopted here.
Figure 4.1b shows an original linear elastic homogeneous, isotropic innite continuum with
the
Lame constants ; and the mass density . At the time t = 0 impulse forces
fi+(t) occur in the medium, e.g. due to an explosion, and a stress wave u+ (x; t) starts
propagating. The stress wave is supposed to cause plastic deformation close to its origin.
Due to geometrical and material damping, the stress wave will eventually attenuate in
intensity, so a borderline S12 exists, which separates the domains V + and V , where
material are supposed to show non-linear and linear behaviour, respectively. The nonlinear domain V + is analysed by means of FEM and the linear domain V by means of
BEM. The elastic eld ud (x; t) and associated surface tractions td(x; t) dissappear on the
articial boundary S0 at innity.
Figure 4.1c shows a linear or non-linear structure V + connected to a surrounding linearelastic homogeneous isotropic innite media V with properties as dened for the system
in gure 4.1a. The interface nodes S12 may now be separated by parts of the surface S2.
The load on the system is caused partly by
a propagating far eld u0 (x; t) as in gure
4.1a and partly by the nodal point loads fi+(t) on the structure from wind, waves, traÆc
etc.
Below, BEM/FEM coupling schemes will be formulated, making it possible to handle
problems of the specied type. Since the FEM formulation may be non-linear the schemes
will be formulated in the time domain. Schemes applying both the direct and the indirect
BEM formulation will be described.
4.2
Coupling of BEM and FEM Formulations
Below, quantities belonging to the alluvial valley/plastic domain/structural system and
to the elastic media will be dened by superscripts + and , respectively. Quantities
belonging to the interior nodes of the nite element mesh and to the boundary S2 are
designated subscript i ("interior") whereas quantities related to the coupling nodes of the
interface S12 are designated subscript b ("boundary").
The FEM can then be written
M+u + + C+u_ + + R+(t) = f + (t)
(4.1)
where u+ is the nodal displacements from the undeformed state, M+ is the mass matrix
and C+ is the linear viscous damping matrix. f +(t) is the external nodal load vector
conjugated to u(t) and R+(t) is the reaction forces from the internal stresses on the nodal
points. The time axis is discretized in n nite time steps of length t. (4.1) must apply
for the beginning of the Lth time step, denoted tL, but also for the end of the Lth time
step tL+1 = tL + t, where
M+u + (tL) + C+u_ +(tL ) + R+(tL ) = f +(tL)
(4.2)
M+u + (tL+1) + C+u_ +(tL+1 ) + R+(tL+1) = f + (tL+1)
(4.3)
42
4. Coupled BEM and FEM Solutions
Subtracting (4.2) from (4.3) yields
M+u+ (tL) + C+u_ + (tL) + R+(tL) = f +(tL )
(4.4)
where e.g. R+(tL) = R+(tL+1) R+(tL ) and likewise for the other contributions in
(4.4). R+(tL) is assumed to depend linearly on the displacement increment u+
R+(tL) = K+(tL)u+
(4.5)
where K+(tL) in case of a non-linear problem expresses the tangential characteristics at
the beginning of the time step, see Section 4.6. (4.5) is inserted into (4.4)
M+u+ + C+u_ + + K+(tL )u+ = f +(tL )
(4.6)
u+ is an assemby of all incremental nodal point displacements related to the nite
+ is
element mesh, and f + are the conjugated incremental
nodal
point
loads.
Next,
u
partitioned into interior degrees of freedom u+i , and boundary degrees of freedom u+b .
The corresponding partitioning of f + , M+, C+ and K+(tL) have been indicated in (4.7)
and (4.8).
"
"
+#
+#
u
f
+
+
i
i
u = u+ ; f = f +
(4.7)
b
b
M+ii+ M++ib ; C+ = C+ii+ C++ib
M"bi Mbb
Cbi Cbb
#
+ (tL ) K+ (tL )
ii
ib
K+(tL) = K
K+bi (tL) K+bb(tL )
M+ =
"
#
"
#
;
(4.8)
It is assumed that u+L = u+(tL ) and u_ +L = u_ +(tL) are known at the time tL = L t.
Based on these values the tangent stiness matrix K+(tL) can then be+calculated,
and
(4.6) can be integrated to obtain the solution u+, +u_ +, which gives uL+1, u_ +L+1 at the
time tL+1 , if only the incremental coupling forces fb from the surrounding continuum
were known.
First, the coupling problem will be formulated in case of a direct formulation of the
BEM. The direct formulation follows the one used in Example 2.1 until (2.49). (2.49) is
formulated at the time tL+1 and solved for tL+1.
tL+1 = K uL+1 + pL+1
(4.9)
K = (D ) 1 = (G1;1 ) 1 21 I + T1;1
(4.10)
where pL+1 = K fL+1 follows from (2.51) at the time tL+1 and with u0;l = t0;l = 0
4.2. Coupling of BEM and FEM Formulations
pL+1 = (G1;1 )
1
L X
l=1
L X
l=1
TL+2 l;1 + TL+1
l;2
GL+2 l;1 + GL+1
ul
l;2
43
tl
!
(4.11)
As was the case for the FEM analysis, (4.9) is next partitioned into parts representing
interior and boundary nodes
"
#
"
#"
#
"
#
ti;L+1 = Kii Kib ui;L+1 + pi;L+1
(4.12)
tb;L+1
Kbi Kbb ub;L+1
pb;L+1
ti;L+1 on the surface S2 is prescribed. For the systems in gure 4.1 one has ti;L+1 = 0,
since the surface is unloaded. From the rst equations of (4.12) it follows that
ui;L+1 = (Kii ) 1 ti;L+1 pi;L+1 Kib ub;L+1
(4.13)
Using (4.13), ui;L+1 can next be eliminated from the last equations of (4.12) giving the
condensed stiness relationships for the interface degrees of freedom
tb;L+1 = K~ bbub;L+1 + ~fb;L+1
(4.14)
K~ bb = Kbb Kbi (Kii ) 1Kib
(4.15)
~fb;L+1 = Kbi (Kii ) 1 ti;L+1 pi;L+1 + pb;L+1
(4.16)
The coupling of the BEM and FEM is now performed using the following boundary
conditions
u+b;L+1 = ub;L+1
(4.17)
+ = At
~
~
fb;L
(4.18)
+1
b;L+1 = AKbb ub;L+1 Afb;L+1
where A is a diagonal matrix with the area attached to each node in its diagonal, transforming the surface traction into nodal forces. This corresponds to a constant shape
function for the surface traction. Improvement can be made by assuming linear
variation
+
of the surface traction. Since the coupling forces at the previous time step fb;L are assumed
to be known, the incremental coupling forces can now be calculated as
+
+
+ =
~ bb u+b;L+1 A~fb;L+1 fb;L
fb+ = fb;L
+1 fb;L = AK
+ =
AK~ bb u+b + u+b;L A~fb;L+1 fb;L
+
AK~ bbu+b + AK~ bbu+b;L A~fb;L+1 fb;L
(4.19)
44
4. Coupled BEM and FEM Solutions
(4.6), (4.7) and (4.8) can then be written in the following form at the time t = tL+1
M+ii+ M++ib u+i+ + C+ii+ C++ib u_ +i+ +
Mbi Mbb ub
Cbi Cbb # " u_ b
"
#"
#
fi+
K+ii (tL)
K+ib (tL)
u+i+ =
+
AK~ bbu+b;L A~fb;L+1 fb;L
K+bi (tL) K+bb(tL) + AK~ bb ub
"
#"
#
"
#"
#
(4.20)
The algorithm now proceeds as follows. + (4.20) is +solved using numerical integration
+ and u_ + , which gives u
_ L+1. From the boundary condition
to obtain u
L+1 and u
+
ub;L+1 = ub;L+1 , ui;L+1 is calculated from (4.13). As integrator the Newmark scheme is
used, see Appendix D.
Next, the indirect BEM formulation is used. From (2.72), (2.73), (2.74) and (2.75) it
follows that
uL+1 = G1;1L+1 + gL+1
(4.21)
tL+1 = 21 I + T1;1 L+1 + hL+1
(4.22)
gL+1 =
L X
hL+1 =
L X
l=1
l=1
GL+2 l;1 + GL+1
TL+2 l;1 + TL+1
l
(4.23)
l
(4.24)
l;2
l;2
From (4.21) it follows that
L+1 = (G1;1 ) 1 uL+1 gL+1
(4.25)
Inserting (4.25) into (4.22) and an equation similar to (4.9) is obtained but with K and
fL+1 dened as
K = 21 I + T1;1 (G1;1)
1
(4.26)
fL+1 = hL+1 (G1;1) 1 gL+1
(4.27)
The remaining analysis proceeds as for the direct method.
4.3. 2D Time Domain Boundary Element Formulation
4.3
45
2D Time Domain Boundary Element Formulation
Green's functions for the displacement eld in 3D, gij (x; t; ; 0) is given in (2.2), but in
the 2D BEM formulation a 2D solution is needed. This is established by assuming plane
strain and by an integration of (2.2) along the third spatial dimension. The integration
can be seen in Eringen & Suhubi [17] and gives the result
8
2
1 < 1 H (c t r) 4 q2c21t2 r2 rr g (x; t; ; 0) =
2 : c1 1
c21 t2 r2 r4
2
1 H (c t r) 4 q2c22t2 r2 rr + Æ qc2 t2 r2 + q
2
c2
c22 t2
r4
r2
r2 5 +
39
Æ 5=
; ; = 1; 2
c22 t2 r2 ;
2
r2
3
Æ q 2 2
ct
r2 1
(4.28)
where x and are 2D vectors with components x and and r = x .
Green's function for the traction eld in 2D t (x; t; ; 0) is found from (4.28), see Israil
& Banerjee [30]
t (x;8t; ; 0) =
1 H c1 t
>
2r >: c1 r
>
>
<
1H
c2
2
c2 t
r
1
6
6
4
2
1
(
6
6
4
(
c2 t 2
r
c1 t 2
r
1
)
)
1
3=2 0 A(1) 1
@
3=2 0 A(3) 1
@
A
r
A
r
2
+ r
2
+ r
c2 t
r
c2 t
r
2
2
2
c1 t
r
c1 t
r
[email protected]
1
0
2
[email protected]
1
0
3
1
2A(2) 7
A7
r
39
1
2A(2) 7>>=
A7
r
5
5>
>
;
; ; = 1; 2
(4.29)
where
r r r1 n1 + r2 n2 + 2 r2
r
r
r
r
r r r1 n1 + r2 n2 (2)
A = n + n + Æ 4 2
r
r r r
r
r
r
n
+
r
n
r
1
1
2
2
A(3)
Æ
n
= 2 r 2
r
r
r
A(1)
= n
(4.30)
Based on identical derivations as presented in the 3D case, Somigliana's identity in 2D
can then be written as
46
4. Coupled BEM and FEM Solutions
C + (x)u+ (x; t) =
Z Z
S
where
t
0
Z Z
S
t
0
g (x; t ; ; 0)t+ (; ; n+ ())d dS( )
t (x; t ; ; 0)u+ (; )d dS( )
(4.31)
1
+
C (x) = > 0
:
x 2 A+
x 2 A
(4.32)
x 2 S
As for the 3D case discretization in time as well as space is necessary in order to solve
(4.31).
8
>
<
4.3.1
;
;
1 ;
2
Space Discretization
The spatial discretization in 2D is made from the same integration scheme as used for the
3D case explained in Section 2.3.1. This leads to
C + (x)u+ (x; t) =
M X
N Z
X
g
Sm
m=1 n=1
M X
N Z
X
m=1 n=1
(x; t; (); 0) N (m;n) ()t(m;n)+(t)d
t
Sm
(4.33)
(x; t; (); 0) N (m;n) ()u(m;n)+(t)d
where M is the number of elements into which the surface have been divided, N is the
number of nodal points within each element, denotes the convolution integral as shown in
(2.18), is the local axis over which the integration is performed see gure 4.2, N (m;n) (),
t(m;n)+ (t) and u(m;n)+ (t) are the shape functions and the nodal values belonging to element
m as dened in (2.21) and (2.25).
Figure 4.2: Denition of local axis .
The integration in (4.33) is performed by the numerical scheme explained in 4.3.2.
4.3. 2D Time Domain Boundary Element Formulation
4.3.2
47
Numerical Space Integration
The numerical scheme used to evaluate the integrations in (4.33) closely follows the
method described in Section 3.3. In 2D a wave front can be represented as a single point
on the element over which the integration is performed, and not as the 3D case where
each wave front is a circular curve on the element. Hence the disection of the dierent
wave parts becomes much simpler. The integration over each wave part is performed by
simple 3-point Gauss quadrature, see Appendix B, and subdivision of the element over
which integration is performed is done as explained in Section 2.4.
The integration is further enhanced by using directional subdivision for cases where weak
singularities occur within a wave part. Directional subdivision is performed by decreasing
the size of the subelements in the direction of the weak singularity, which gives higher
accuracy with fewer subelements than subdivision where a constant subelement size is
used. Element types and shape functions can be seen in Appendix B.
4.3.3
Singularity Problem
As explained in Section 2.4, Green's functions are singular when the source and the
reciever x points coinside. In the 2D case the singularity of Green's function for the
displacement eld is weak (proportional to ln(jx j)) and can therefore be evaluated
numerically by Gauss quadrature. But
for Green's function for the traction eld, where
the singularity is proportional to jx 1 j as for the displacement eld in the 3D-case, a
special procedure is implemented.
The singularity only arises for the rst time step, and if zero initial conditions are assumed,
only singular integrals of the following type needs to be evaluated
T2D
1;1 =
"Z
t
Z
Sp
0
t (x(p;q) ; t; ( ); 0)(1
1 ( ))N (p;q) ( )d d
#
(4.34)
where RSp signies integration over element p and x(p;q) denotes the coordinates for node
q on element p.
An analytical solution of the singularity in the 3D case was possible over a circular area,
but in the 2D case a dierent approach is necessary. An analytical approach should
integrate the singularity over a line segment instead of a circular area. But it has been
chosen to use a numerical approach instead. The chosen method is the well-known rigid
body motion principle, see e.g. Banerjee & Ahmad [10].
An expression similar to (4.34) can be written for the static case
T2D
static =
Z
S
tstatic (x(p;q) ; ( ))N (p;q) ( )d
p The static solution tstatic
on [27]
(x; ) is given as, see e.g. Gomez-Lera & Alarc
(4.35)
48
4. Coupled BEM and FEM Solutions
tstatic
(x; ) =
(1 2 )
4(1
r
r
n
1
)r
r n
r
r1 n1 + r2 n2 (1
r
2 )Æ + 2 rrr2 (4.36)
where is the Poisson ratio and E is the elasticity modulus. Combining (4.34) and (4.35)
leads to
1 I + T2D = 1 I + T2D + T2D T2D =
1;1
static
static
2
2 1;1
1 I + T2D + T2D T2D (4.37)
static
1;1
static
2
where the last part of the right hand side is non-singular, duec tto the fact that the singui
larities in t (x; t; (); 0) and tstatic
(x; ( )) are identical as r ! 1 and r ! 0. The rst
part of the right hand side in (4.37) is evaluated using the rigid body motion principle for
the static case. If an unit translational ridig body motion is applied to the static version
of (4.33) and the fact that t = 0 is used, the following result is obtained
#
N X
N Z
1I = " X
static
(
p;q
)
(
m;n
)
t (x ; ( ))N
()d
(4.38)
m 2
m=1 n=1 S
which leads to
1 I + Z tstatic (x(p;q); ())N (p;q)()d = 1 I + T2D =
static
2 2 Sp 2
N Z
X
(p;q) ; ( ))N (p;n) ( )d +
4
tstatic
(x
Sp
n=1;n6=q
N
N Z
X
X
m=1;m6=p n=1
Sm
3
(p;q) ; ( ))N (m;n) ( )d 5
tstatic
(x
(4.39)
Note that for halfspace problems a special procedure instead of (4.39) is needed, see
Section 4.3.4. With (4.37) and (4.39) a full calculation of (4.34) is now possible.
In Example 4.1 numerical tests of the described procedure are performed and compared
to analytical results.
4.3.4
Enclosing Elements
The procedure explained in Section 4.3.3 will only work for problems with a closed boundary, due to the use of a rigid body motion. For problems with innite boundaries the
summations in (4.39) will not give the correct results. If the unit translatorial stibody
motion used in Section 4.3.3 is applied to the static version of (4.33), for a problem with
4.3. 2D Time Domain Boundary Element Formulation
49
an innite boundary, the result will not be zero traction for all nodes and (4.38) and
(4.39) do not apply. Therefore a special procedure is needed in those cases.
For problems with innite boundaries the element mesh is extended with a ctitious enclosing boundary solely for the evaluation of (4.39), as suggested by Ahmad & Banerjee
[3]. Since the enclosing elements are only used for the evaluation of the singularity for
Green's function for the traction eld, they will have no inuence on the other integrations.
If L enclosing elements are added to an element mesh with N elements, (4.39) can be
reformulated as
2
1 I + T2D =
static
2
N
X
4
N Z
X
m=1;m6=p n=1
LX
+N LX
+N
Sm
m=N +1 n=N +1
N
X
n=1;n6=q
Z
static (x(p;q) ; ( ))N (p;n) ( )d +
t
Sp
(p;q) ; ( ))N (m;n) ( )d +
tstatic
(x
Z
Sm
3
(p;q) ; ( ))N (m;n) ( )d 5
tstatic
(x
(4.40)
Figure 4.3 shows an example with a halfspace and an enclosing boundary.
Figure 4.3: Halfspace model with enclosing boundary.
The choice of enclosing boundary cannot be made freely. The chosen boundary must
ensure accurate results. The prime demand is that the distance between the enclosing
elements and the model elements is suÆciently large, so that the enclosing elements will
not inuence the results for the model elements. This distance depend on the wave
velocities and in Example 4.1 a non-dimensional parameter f for the distance is found.
Through the material parameters f depends on the wave velocities.
Furthermore, sharp corners between the model elements and the enclosing elements must
be avoided. If the model elements and enclosing elements meet at sharp corners, the jump
term will change for the model nodes adjoined to the enclosing boundary. The jump term
is equal to 2 , where is the angle between two elements, so e.g. at a node where
two elements meet at 90Æ the jump term will change to 0.25. If this happens at the node
between the model elements and the enclosing elements, the boundary element model will
react as if the elastic halfspace outside the model boundary is at 90Æ compared to the
model boundary.
50
4. Coupled BEM and FEM Solutions
In order to full the above requirements the optimal enclosing boundary must have two
plane line segments at the end of the modelled boundary. In gure 4.3 an example of such
an enclosing element mesh is shown.
In Example 4.1 dierent examples of enclosing boundaries are shown and compared to
analytical results.
Example 4.1
Tests of Singularity Procedure
The procedures for calculation of the singular integrals for Green's function for the traction eld as
explained in Sections 4.3.3 and 4.3.4 are tested using a simple example where the solution is known.
A plane 2D isotropic, homogeneous, linear elastic halfspace consisting of 10 elements with 11 nodes is
used as test example. In gure 4.4 the surface model nodes and elements are shown and in gure 4.5
the enclosing boundaries are sketched. The test is performed using the Courant numbers for the P- and
S-waves respectively C = htc1 =0.906, C = htc2 =0.523 and = .
h
1
2
1
3
2
4
3
5
4
6
5
7
6
8
7
9
8
10
9
11
10
Figure 4.4: Model boundary for test example.
Figure 4.5: Halfspace model. A) Circular enclosing boundary without line segments B) Circular
enclosing boundary with line segments C ) Square enclosing boundary with line segments.
For plane elements, as the ones used in the model elements, the integration of t(x; t; ; 0) gives identical
zero, see (4.29). Therefore (4.40) must give exactly 0.5 on the diagonal.
Each of the three enclosing boundaries shown in gure 4.5 is tested using 5 dierent levels of integrational
subdivision: 2 , 5, 10, 15 and 25 subdivisions (see Section 4.3.2). A is tested using 5 dierent discretizations of the enclosing boundary. B is tested using 5 dierent lengths b1 of line segments and C is tested
with 5 dierent heights b2, see gure 4.5.
First A is tested using 5, 8, 11, 16 and 22 nodes on the enclosing boundary, the resulting average errors
are shown in table 4.1.
The averaged errors in table 4.1 are very small and it is clear that with more integrational subdivisions,
e.g. more precise integration, the error drops signicantly. The same is the case with increasing numbers
of nodes on the enclosing boundary. In gure 4.6 a graphic representation of the results can be seen and
in table 4.2 the diagonal terms for the case with 22 nodes on the enclosing boundary and 25 integration
subdivisions are shown.
4.3. 2D Time Domain Boundary Element Formulation
Number of nodes
on enclosing boundary
5
8
11
16
22
2
2:32 10
1:26 10
2:62 10
1:06 10
7:80 10
4
5
7
7
9
Number of subdivisions
5
10
15
7
9
8:07 10
4:48 10 3:66 10
8
2:81 10 3:19 10 10 2:69 10
3:50 10 9 4:74 10 11 4:07 10
2:75 10 10 3:97 10 12 3:44 10
1:78 10 11 2:59 10 13 2:23 10
51
10
11
12
13
14
25
1:65 10
1:23 10
1:88 10
1:58 10
1:09 10
11
12
13
14
15
Table 4.1: Average errors of diagonal terms for type A enclosing boundary if the last node
at each end
jval 0:5j
of the element mesh is not included in the average errors. The error is dened as error =
val is the numerical value.
Node Direction Value
1
1
0.2464
1
2
0.2319
2
1
0.5000
2
2
0.4999
3
1
0.5
9
10
10
11
11
2
1
2
1
2
0:5
, where
0.5
0.5000
0.4999
0.2464
0.2319
Table 4.2: Diagonal terms for type A with 22 nodes on the enclosing boundary and 25 integration
subdivisions.
In table 4.2 it is clear that the jump term will change if the model elements and enclosing elements meet
in sharp corners. The jump term for nodes 1 and 11, e.g. the nodes at each end of the element mesh, are
close to 0.25 since the angle of the intersections is close to 90Æ.
In gure 4.6 a) it is clear that more than 11 nodes on the enclosing boundary give no further signicant
gain in accuracy. Instead of using more nodes, further accuracy can be achieved through more precise
integration. But gure 4.6 b) shows that 10 integration subdivisions are the highest reasonable choice.
Next the type B enclosing boundary is tested using 11 nodes on the enclosing boundary and 2 nodes on
each line segment.
52
4. Coupled BEM and FEM Solutions
b)
1
1
0.9
0.9
0.8
0.8
Value normalized with respect to first value
Value normalized with respect to first value
a)
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0
0
5
10
15
20
Number of nodes on enclosing boundary
5
10
15
20
Integration subdivision
25
Figure 4.6: a) Normalized average errors for 10 integration subdivsions. b) Normalized average errors
with 11 nodes on the enclosing boundary.
Number of subdivisions
l
2
5
10
15
25
4
5
7
2f
0.0047 2:00 10
4:25 10
4:79 10
6:63 10 8
10 f 7:46 10 5 1:30 10 7 1:17 10 9 9:48 10 11 4:26 10 12
15 f 5:51 10 6 1:15 10 8 1:41 10 10 1:19 10 11 5:47 10 13
20 f 2:68 10 7 2:62 10 9 3:62 10 11 3:11 10 12 1:43 10 13
30 f 1:62 10 7 3:30 10 10 4:74 10 12 4:09 10 13 1:87 10 14
Table 4.3: Average relative errors of diagonal terms for type B enclosing boundary. f = 41(12) . The
error is dened as error = jval0:50:5j , where val is the numerical value.
The parameter f which is chosen as the length factor for the line segments, is given as
1 2
f=
(4.41)
4(1 )
f is the deciding factor for the length of the line segments. is the equivalent plane strain found when
modifying the 3D parameters to 2D. If the problems with the end nodes, as shown in table 4.2, is to be
avoided, the model nodes may not inuence the enclosing boundary nodes at the sharp corner. Since
tstatic
(x; t; ; 0) is the kernel function from which the jump term is calculated, it is this function that
must be without signicant inuence at the corner nodes. tstatic
(x; t; ; 0) behaves proportionally with
1 2 1 = f and therefore, f becomes the deciding factor. Note that f depends on the wave velocities
4(1 ) r r
since the wave velocities are functions of the Lame constants which again are related to and E .
It is clear from table 4.3 that 10 f ensures accurate results, which is not unexpected due to the hyperbolic
behaviour of tstatic
(x; t; ; 0).
Finally, the type C enclosing boundary is tested using line segments of the length 10 f , 2 nodes on each
line segment, 2 nodes on each vertical part of the enclosing boundary and 11 nodes on the horizontal part
of the enclosing boundary.
As was the case for the length of the line segments in the type B boundary, it is clear from table 4.4
that the gain in accuracy is negligible above 10 f . Even from 8 f and to 10 f the error changes only
slightly. In gure 4.7 the averaged errors for 10 integration subdivisions have been shown graphically.
4.3. 2D Time Domain Boundary Element Formulation
Number of subdivisions
5
10
15
4
0:00341 1:32 10
1:71 10 6
2:85 10 5 1:63 10 7 1:25 10 10
5:81 10 8 1:41 10 10 1:22 10 11
4:06 10 9 8:08 10 11 7:01 10 12
9:80 10 10 1:45 10 11 1:26 10 12
2
0:0927
0:00684
3:01 10 4
9:22 10 5
3:83 10 6
h
2f
4f
8f
10 f
15 f
53
25
2:65 10 9
1:94 10 12
5:66 10 13
3:25 10 13
5:83 10 14
Table 4.4: Average relative errors of diagonal terms for type C enclosing boundary. f = 41(12) . The
error is dened as error = jval0:50:5j , where val is the numerical value.
0
10
−1
Value normalized with respect to first value
10
−2
10
−3
10
−4
10
−5
10
−6
10
2
4
6
8
10
12
Length of vertical part on enclosing boundary
14
Figure 4.7: Normalized average errors for 10 integration subdivsions. The y-axis is logarithmically
scaled.
From gure 4.7 it becomes even clearer that the gain in accuracy declines dramatically above 8 f .
Therefore, to ensure accurate results when using an enclosing boundary for innite domain problems,
line segments of 8 to 10 times the size of f are required. And if a square enclosing boundary is used, the
height of the vertical parts must be at least 8 to 10 times the size of f .
Example 4.2
Tests of Numerical Space Integration
In this example the numerical integration scheme is tested using a simple example.
A plane halfspace consisting of 10 elements with 11 nodes as shown in gure 4.4 is used as a test example.
The test is performed using the Courant numbers for the P- and S-waves, C =0.906, C =0.523 and = ,
respectively.
Since all elements have equal size and are plane the integration of g (x; t; ; 0) and t (x; t; ; 0) is the
same for all locations of x. Therefore, results for the integration are in the following only shown for x
located at node 6. The integration results shown are all for Green's functions multiplied by the shape
function for the second node. The second node for an element is the right hand node. Again, since
the elements are all plane and of the same size, integration of Green's functions multiplied by the shape
function for the rst node would give the numerically same results, but for the opposite elements.
54
4. Coupled BEM and FEM Solutions
Element number
Number of 5 6 10 11
subdivisions 10 11
2
1.2849378 6.0742377
5
1.2842802 6.1748658
10
1.2842665 6.1855266
15
1.2842654 6.1858546
20
1.2842654 6.1858650
25
1.2842654 6.1858658
Table 4.5: g1;1 multiplied by shape function for second node integrated over elements for x at node 6,
step 1.
Element number
Number of
5
6
subdivisions 10 2
10 2
2
-1.0265921 -3.1594396
5
-1.0217577 -3.1643247
10
-1.0212882 -3.1644831
15
-1.0212738 -3.1644841
20
-1.0212733 -3.1644841
25
-1.0212719 -3.1644841
Table 4.6: t1;2 multiplied by shape function for second node integrated over elements for x at node 6,
step 1.
From tables 4.5 and 4.6 it is clear that the accuracy is excellent for both g1;1(x; t; ; 0) and t1;2(x; t; ; 0)
at the rst time step. Accurate results with 6 decimals are achieved from 15 subdivions and onwards.
Number of
subdivisions
2
5
10
15
20
25
4 10 12
3.3980781
3.3986854
3.3987161
3.3987236
3.3987258
3.3987267
Element number
5 6 10 11
10 11
5.7452389 7.5746247
5.7449037 7.5746145
5.7449003 7.5746139
5.7448998 7.5746138
5.7448997 7.5746138
5.7448996 7.5746138
7 10 11
1.2809548
1.2809701
1.2809711
1.2809713
1.2809713
1.2809713
Table 4.7: g1;1 multiplied by shape function for second node integrated over elements for x at node 6,
step 2.
4.3. 2D Time Domain Boundary Element Formulation
Number of
subdivisions
2
5
10
15
20
25
4 10 3
5.0970140
5.0973609
5.0974694
5.0974834
5.0974884
5.0974907
Element number
5 6 2
10
10 2
-5.0711486 4.8220696
-5.0774486 4.8223557
-5.0774499 4.8223561
-5.0774501 4.8223561
-5.0774501 4.8223561
-5.0774501 4.8223561
55
7 10 2
-2.0352030
-2.0355476
-2.0356061
-2.0356205
-2.0356265
-2.0356297
Table 4.8: t1;2 multiplied by shape function for second node integrated over elements for x at node 6,
step 2.
Table 4.7 with the results for step 2 are equivalent with the tables 4.5 and 4.6. Excellent accuracy is
achieved with 15 subdivisions. But the integration results for t1;2(x; t; ; 0) at step 2 in table 4.8 show
less accuracy for especially element 7. t1;2(x; t; ; 0) multiplied by the shape function for the second node
is drawn over elements 4-7 for step 2 in gure 4.8.
0.5
0.4
0.3
0.2
7
5
4
0.1
6
0
−0.1
−0.2
6
−0.3
−0.4
Figure 4.8:
−0.5
−0.8
−0.6
−0.4
−0.2
0
0.2
0.4
0.6
0.8
t1;2 multiplied by the shape function for second node for x at node 6, step 2. Numbers
within circles are element numbers and dashed lines signify element boundaries. Wave boundaries are
shown as asterisks.
Since the singularity in element 7 is weak compared to the singularities in elements 5 and 6, it has
been chosen not to use directional subdivision for element 7, but it has been used for elements 5 and 6.
Therefore, convergence is a little slower for element 7.
56
4. Coupled BEM and FEM Solutions
Number of
subdivisions
2
5
10
15
20
25
3 10 12
1.5751564
1.5775893
1.5776236
1.5776282
1.5776295
1.5776301
4 10 11
2.8929884
2.8929836
2.8929831
2.8929830
2.8929830
2.8929830
Element number
5 6 10 11
10 11
3.8776810 3.3864729
3.8777031 3.3864780
3.8777049 3.3864783
3.8777051 3.3864783
3.8777051 3.3864783
3.8777052 3.3864783
7 10 11
3.9765854
3.9765832
3.9765831
3.9765830
3.9765830
3.9765830
8 10 12
6.4883319
6.4884339
6.4884453
6.4884470
6.4884475
6.4884477
Table 4.9: g1;1 multiplied by shape function for second node integrated over elements for x at node 6,
step 3.
Number of
subdivisions
2
5
10
15
20
25
3
10 3
1.2850657
1.2852953
1.2853170
1.2853199
1.2853208
1.2853212
4
10 2
1.0042238
1.0059655
1.0063460
1.0064407
1.0064804
1.0065013
Element number
5
6
10 3
10 2
-7.7916867 -1.0787177
-7.7167173 -1.0856133
-7.6999132 -1.0864872
-7.6956923 -1.0868639
-7.6939199 -1.0870231
-7.6929841 -1.0871074
7
10 2
4.8081990
4.8069494
4.8066666
4.8065956
4.8065657
4.8065500
8
10 3
-5.7724101
-5.7724883
-5.7724956
-5.7724967
-5.7724970
-5.7724971
Table 4.10: t1;2 multiplied by shape function for second node integrated over elements for x at node 6,
step 3.
Again, the integration results for step 3 listed in tables 4.9 and 4.10 show execellent convergence up to 6
decimals accuracy. It can hereby be concluded that the integration routines work satisfactorily.
4.3. 2D Time Domain Boundary Element Formulation
4.3.5
57
Time Discretization
The time interval [0; t] is divided into L subintervals, each of the length t = Lt . For the
direct boundary element method (2.33)-(2.38) remain valid, if the 3D Green's functions
gij (x; t; ; 0) and tij (x; t; ; 0) are replaced by the 2D versus (4.28) and (4.29).
The integrals entering the matrix elements on the right hand sides of (2.35)-(2.38) have
all been calculated analytically, Israil & Banerjee [30]
"
(
"
H (1 1) tl Æ 1
1 (1 )
g (x; ; ; 0)(1 l )d =
cosh
2
2
c1
t 2
tl
"
#)
r Æ 3 r r 2 3
r r
1 1 +
+
+
2
r
c1 t 3 1
r2 3 1 1
(
"
#
H (2 1) tl Æ r r
1
c22
t 2 cosh (2) + 22 r2 22 +
"
#)#z =tl
r
Æ 3
r r 2 3
c t
3 (2 + 32) + r2 3 2 + 2
Z
tl
1
2
z =tl
+
(4.42)
1
"
(
"
1
H (1 1) tl Æ g (x; ; ; 0)l d =
1 t 2 cosh 1 (1)
2
c21
tl
"
#)
r r
r Æ 3 r r 2 3
+
+
r2 1 1 c1 t 3 1
r2 3 1 1
(
"
#
H (2 1) tl Æ r r
1
1 t 2 cosh (2) + 22 r2 22
c22
"
#)#z =tl
r
Æ 3
r r 2 3
c2 t
3 (2 + 32) + r2 3 2 + 2 z=tl
Z
1 1
tl
1
1 1
+
(4.43)
1
H (1
t (x; ; ; 0)(1 l )d =
2r
c21
tl 1
"
#)
r 1 (1)
2
H (2
(2)
3
A 2A 1 + 1
+
c1 t 1
3
c22
"
#)#
r
1 A(3) + 2A(2) 2 3 + z=tl
2
3 2
c2 t 2 z =tl 1
Z
"
tl
1) ( tl " 1 A(1) + 2 A(2) # +
t " 1 1 1 #
(
1) tl 2 A(3) 2 A(2) +
t 2 2 2 (4.44)
tl H (1 1) 1
t (x; ; ; 0)l d =
2r #) c21
t
tl 1
(
"
r 1 (1)
2 3 + + H (2 1) 1
A 2A(2)
c1 t 1
31 1
c22
"
#)#
r
1 A(3) + 2A(2) 2 3 + z=tl
2
3 2
c2 t 2 z =tl 1
Z
tl
"
(
"
#
1 (1)
A + 21 1 A(2)
1 #
"
tl 2 (3)
(2)
t 2 A 222 A
(4.45)
58
4. Coupled BEM and FEM Solutions
q
(2) (3)
where = cr z , = 2 1 and A(1)
, A , A are dened in (4.30). From (2.33) and
(2.34) it can be seen that for each time interval the only new system matrices (eg. matrices
that have not been evaluated for an earlier time interval) that need to be evaluated are
2D
2D
2D
(G2D
l;1 +Gl 1;2 ) and (Tl;1 +Tl 1;2 ). Therefore, alternatively to (2.35)-(2.38) the combined
2D
2D
2D
expressions for (G2D
l;1 + Gl 1;2 ) and (Tl;1 + Tl 1;2 ) are used, whereby the computational
eort will be greatly reduced.
Hereafter the BEM formulations from Sections 2.6 and 2.7 or Section 4.2 can be used
directly.
Example 4.3
Results from 2D FEM Implementation
Before the BEM/FEM formulation is tested, separate tests of the BEM and FEM formulations are
performed. The formulations are tested using a distributed constant downward directed vertical Heaviside
type uniform line load on an isotropic, homogeneous, linear elastic 2D plane strain halfspace. This
example has been chosen to make comparison with Mansur [50] and Antes [5] possible. The distributed
load per unit length p(x; t) is applied as
p(x; t) = p0 [H (x1 + b)
H (x1
b)]H (t)
(4.46)
The analysis is performed using the following parameters : Intensity of line
load p0 = 104 ps/in (6:896 104
2
KN/m),
b = 3000 in (1.146.4 m), E = 2:5 106 ps/in (1:724 107 KN/m ), = 2:891 10 3 ps2 =in4 (3:151
3
kg/m ) and = 0:25 (modied for 2D
plane strain). On the assumption of plane strain this gives
= = 106 ps/in (6:895 106 KN/m2 ) and the wave velocities c1 = 3:27 104 ips (inches per second)
(8:306 102 m/s) and c2 = 1:86 104 ips (4:742 102 m/s).
For the FEM analysis a semicircular region with radius R = 1:5 104 in is discretisized, see gure 4.9.
The elements are 3 node triangle constant strain elements with linear interpolation, see Appendix B.1.2.
When evaluating the results reection at the boundaries must be taken into account.
R
b
h
B
C
0
A
−5000
−10000
−15000
−15000
−10000
−5000
0
5000
10000
15000
Figure 4.9: Element mesh for !n = 140 rad/s ) h = 3000 in generated with automatic mesh generation
tool. Length is given in inches.
Before the analysis can be performed an investigation of the load is necessary. Based on a Fourier transform F (!) of the spectrum a limiting cut-o frequency !0 with the signicant circular frequencies in
the oscillation has to2 be chosen. For the Heaviside unit step function the Fourier spectrum decreases
relatively slowly as ! , so a rather large value of !0 needs to be chosen. The minimum period associated
with !0 becomes T0 = 2!0 . In the numerical integration the period T0 is covered in n time steps, so
4.3. 2D Time Domain Boundary Element Formulation
59
t = Tn0 = n!20 . The element length h is adjusted so the Courant number for the P-waves becomes,
c1 . ! n is varied through the values ! n=140 rad/s, 280 rad/s, 400
C = c1ht = 12 ) h = 2c1 t = 4n!
0
0
0
rad/s and 560 rad/s, corresponding
to
the
time
intervals
t=0.04488 s, 0.02244 s, 0.01571 s and 0.01122
s. With c1 = 3:3422 104 ips the corresponding elements lengths become h =3000 in, 1500 in, 1000 in
and 750 in. In the following !0 will simply be denoted !.
Discretization of the region in gure 4.9 is made by an automatic mesh generation tool which uses the
constrained Delaunay triangulation method to generate 2D unstructured meshes with 3 node triangular
elements, Rebay [65].
In gure 4.9 the element mesh for !n = 140 rad/s ) h = 3000 in is shown.
In gure 4.10 and 4.11 the results from the four cases are compared.
1
A
0
B
−1
C
−2
inches
−3
−4
−5
−6
−7
ωn=140
ωn=280
ωn=400
ωn=560
−8
−9
−10
0
0.1
0.2
0.3
0.4
0.5
s
0.6
0.7
0.8
0.9
1
Figure 4.10: u1 found from FEM analysis. [!n] =rad/s. A : (x1 ; x2) =(0 in, 0 in); B : (x1 ; x2)=(6000
in, 0 in); C : (x1 ; x2 )=(12000 in, 0 in).
From gures 4.10 and 4.11 it is evident that !n = 400 rad/s give satisfactory results for this problem.
Comparison between these results and the results obtained by Antes, Mansur and the 2D BEM implementation can be seen in Section 4.4.
In gure 4.11 for point A the eect of reection clearly arises at t = 0:9 s. The reected waves have
to travel between 27000 in and 33000 in due to the distribution of the load, this correspond to a travel
time for the P-wave between 1.0 s and 0.83 s which explains the observed reection at t = 0:9 s. Further
it is clear from gure 4.10 that the boundary conditions give very small horizontal displacements at the
surface.
For point B the reection inuences the displacement at t = 0:9 s in gure 4.11. In gure 4.10 the
increase in the displacement takes place at both t = 0:6 s and t = 0:8 which is in agreement with the
theoretical travel time for the P-wave at point B, which is between t = 0:64 s and t = 0:83 s.
60
4. Coupled BEM and FEM Solutions
5
0
C
A
−5
B
−10
inches
−15
−20
−25
−30
ωn=140
ωn=280
ωn=400
ωn=560
−35
−40
−45
0
0.1
0.2
0.3
0.4
0.5
s
0.6
0.7
0.8
0.9
1
Figure 4.11: u2 found from FEM analysis. [!n] =rad/s. A : (x1 ; x2) =(0 in, 0 in); B : (x1 ; x2)=(6000
in, 0 in); C : (x1 ; x2 )=(12000 in, 0 in).
4.4
Results from 2D BEM Implementation
The direct 2D BEM implementation described in Section 4.3 is tested on two classic
problems from soil dynamics. First the problem of a Heaviside type uniform line load on
an elastic halfspace that was presented in the previous example is analysed. The results
from this problem will be compared to the FEM results and to numerical and analytical
results from the literature.
The second problem is very similar to the rst, but the Heaviside type load is replaced
by an impulse. This is Lamb's problem in 2D and the results are compared to the results
obtained by Lamb.
Finally, results from the direct BEM formulation are compared to similar results obtained
by the indirect BEM.
Example 4.4
Heaviside Load on Elastic Halfspace
The problem of a distributed constant downward directed vertical Heaviside type uniform line load on
an isotropic, homogeneous, linear elastic 2D plane strain halfspace which was described in Example 4.3
is analysed using the direct BEM formulation from Section 2.6.
The same parameters as in Example 4.3 are used. Additionally, the element lengths and size of time
steps have been chosen similar to the FEM analysis. The elements are 2 node line elements with linear
interpolation, see Appendix B.1.1, corresponding to the elements with linear interpolation used for the
FEM analysis.
The results are shown in gures 4.12 and 4.13.
4.4. Results from 2D BEM Implementation
61
2
0
−2
inches
B
C
−4
−6
h=3000 ∆t=0.04488
h=1500 ∆t=0.02244
h=1000 ∆t=0.01571
h=750 ∆t=0.01122
−8
−10
0
0.1
0.2
0.3
0.4
0.5
s
0.6
0.7
0.8
0.9
1
0.9
1
Figure 4.12: u1 found from BEM analysis. B : (x1 ; x2 )=(6000 in, 0 in); C : (x1 ; x2 )=(12000 in, 0 in).
5
0
C
A
−5
B
−10
inches
−15
−20
−25
−30
−35
h=3000 ∆t=0.04488
h=1500 ∆t=0.02244
h=1000 ∆t=0.01571
h=750 ∆t=0.01122
−40
−45
−50
0
0.1
0.2
0.3
0.4
0.5
s
0.6
0.7
0.8
Figure 4.13: u2 found from BEM analysis. A : (x1 ; x2) =(0 in, 0 in); B : (x1 ; x2)=(6000 in, 0 in); C :
(x1 ; x2 )=(12000 in, 0 in).
In gures 4.12 and 4.13 it is clear that the eects of reection that arose for the FEM analysis are not
present for the BEM analsysis. Since h = 1000 in and t = 0:01571 s give satisfactory results for both
the FEM and BEM analysis, these results will be compared to previous results obtained by an early BEM
method by Mansur [50] and Antes [5] and to analytical results from Appendix C.
62
4. Coupled BEM and FEM Solutions
2
0
−2
inches
B
C
−4
−6
BEM h=1000 ∆t=0.01571
FEM h=1000 ∆t=0.01571
Antes
−8
−10
0
0.1
0.2
0.3
0.4
0.5
s
0.6
0.7
0.8
0.9
1
0.9
1
Figure 4.14: u1 found from BEM and FEM analysis compared with the results obtained by Antes. B :
(x1 ; x2 )=(6000 in, 0 in); C : (x1 ; x2)=(12000 in, 0 in).
5
0
C
A
−5
B
−10
inches
−15
−20
−25
−30
BEM h=1000 ∆t=0.01571
FEM h=1000 ∆t=0.01571
Antes
Mansur
Analytical solution
−35
−40
−45
−50
0
0.1
0.2
0.3
0.4
0.5
s
0.6
0.7
0.8
Figure 4.15: u2 found from BEM and FEM analysis compared with the results obtained by Mansur
and Antes and the analytical solution. A : (x1 ; x2 ) =(0 in, 0 in); B : (x1 ; x2 )=(6000 in, 0 in); C :
(x1 ; x2 )=(12000 in, 0 in).
Figures 4.14 and 4.15 clearly show good agreement with the result from Mansur [50] and Antes [5],
whereas the results obtained by the FEM diverge due to reection eects.
Comparison with the analytical results is not shown for u1 since the analytical solution for this component
does not include the Rayleigh wave component. This leads to very divergent results for the BEM and
FEM solutions compared to the analytical solution. The BEM results compared to the analytical solution
for point B in gure 4.15 show very good agreement, whereas the result diverges slightly for point C .
4.4. Results from 2D BEM Implementation
63
Analytical results for point A cannot be obtained due to the singularities when the analytical solution is
integrated numerically, see Appendix C.
Example 4.5
Lamb's Problem in 2D
A downward directed vertical concentrated unit line impulse is acting at the origin on the boundary of
an isotropic, homogeneous, linear elastic 2D plane strain halfspace at the time t = 0. In Appendix C the
analytical solution to this problem can be seen. The BEM analysis is performed with the load applied
over a nite area, therefore the analytical solution is numerically integrated over a nite area before
comparison with the BEM results, see Appendix C.
The BEM analysis for this problem is made using the same material parameters and element mesh as in
Example 4.4. The load is applied over a single element and results are shown for point (x1 ; x2 )=(12000
in, 0 in).
1
0.5
πµxu1/c2
0
−0.5
−1
−1.5
−2
−2.5
0
2D BEM
Analytical solution
0.5
1
1.5
2
2.5
3
3.5
c1t/x1
Figure 4.16: Numerical and analytical solution for u1 at point (x1 ; x2)=(12000 in, 0 in).
Figure 4.16 shows good agreement between the BEM result and the analytical solution for the rst part
of the time history. The last part of the BEM result does not agree with the analytical solution. The
reason is that the analytical solution for the horizontal displacement is only given until the arrival of
the S-wave, which means that the Rayleigh wave component is missing. So the dience between the
numerical result and the analytical solution is an expression of the inuence of the Rayleigh wave.
64
4. Coupled BEM and FEM Solutions
3
πµxu2/c2
2
1
0
−1
−2
−3
−4
2D BEM
Analytical solution
−5
0
0.5
1
1.5
2
2.5
c1t/x1
3
3.5
4
4.5
5
Figure 4.17: Numerical and analytical solution for u2 at point (x1 ; x2)=(12000 in, 0 in).
Figure 4.17 on the other hand shows good agreement between the BEM result and the analytical solution
for the entire time history.
It should be noted that this problem is known to be very diÆcult to handle with numerical methods.
The problem is the singularity of the impulse, which imposes problems for most numerical methods like
the FEM.
Example 4.6
Direct BEM Compared to Indirect BEM
The problem from Example 4.4 is analysed using both the direct BEM formulation from Section 2.6 and
the indirect BEM formulation from Section 2.7.
In gure 4.18 and 4.19 the results are compared.
1
0
−1
inches
−2
B
−3
C
−4
−5
Direct BEM
Indirect BEM
−6
Figure 4.18:
0
0.1
0.2
0.3
0.4
0.5
s
0.6
0.7
0.8
u1 for Heaviside type load on elastic halfspace.
(x1 ; x2 )=(12000 in, 0 in).
0.9
B
1
: (x1 ; x2 )=(6000 in, 0 in);
C
:
4.5. Results from Coupled BEM/FEM Implementation
65
5
0
C
A
−5
B
−10
inches
−15
−20
−25
−30
−35
Direct BEM
Indirect BEM
−40
−45
0
0.1
0.2
0.3
0.4
0.5
s
0.6
0.7
0.8
0.9
1
Figure 4.19: u2 for Heaviside type load on elastic halfspace. A : (x1 ; x2) =(0 in, 0 in); B : (x1 ; x2)=(6000
in, 0 in); C : (x1 ; x2 )=(12000 in, 0 in).
It is very clear from gures 4.18 and 4.19 that the direct and indirect formulations give identical results.
A closer look at the results reveals that the dierences are less than 0.1%.
The indirect formulation gives a reduction in calculation time compared to the direct formulation for this
problem. For the 4 dierent element meshes used in Example 4.4, reductions of approximate 4%, 6%, 8%
and 11% are obtained. From this comparison it seems as if the indirect BEM formulation is preferable to
the direct BEM. For problems, however, where mixed boundary conditions are necessary, the direct BEM
will have an advantage due to its inherent ability of solving displacements and tractions simultaneously.
On the other hand, this ability requires more storage space than the indirect BEM and therefore, the
choice between the direct and indirect BEM will depend on the specic problem.
Some prefer the indirect formulations over the direct formulation since the equations are simpler and
lead to intuitive visualizations of problems. However, this analysis shows that the choice cannot be made
based on accuracy or calculation speed, but merely on the `look' of the formulations.
4.5
Results from Coupled BEM/FEM Implementation
The coupled BEM/FEM formulation described in Section 4.2 is used to analyse the problem of a distributed constant Heaviside type uniform line load on an isotropic, homogeneous, linear elastic 2D plane strain halfspace. The problem is similar to the one analysed
in Example 4.3 and 4.4, except that in this case the elastic halfspace is occupied by an
alluvial valley.
A circular alluvial+valley
with+the +radius R is considered, made up of two layers with the
+
Lame constants 1 , 1 and 2 , 2 , respectively. In the centre of the valley a distributed
Heaviside load of the width 2b and intensity p0 is applied, see (4.46). The problem is
shown in gure 4.20.
Example 4.7
Heaviside Load on Elastic Halfspace
In+order to test the accuracy of the method an elastic halfspace is rst considered, i.e. +1 = +1 = +2 =
2 = = . Here an articial boundary for the FEM is introduced. The same material properties
and loading parametes as indicated in Example 4.3 and Example 4.4 have been chosen.
66
4. Coupled BEM and FEM Solutions
Figure 4.20: BEM/FEM coupling for alluvial valley, examples.
The time step has been chosen as t = 0:01571 s and for both the FEM and BEM meshes an average
element length of h = 1000 in, has been used. This results in 960 FEM elements and 74 BEM elements.
5000
a)
0
−5000
−10000
−15000
−40000
5000
−20000
0
20000
40000
−20000
0
20000
40000
−20000
0
20000
40000
b)
0
−5000
−10000
−15000
−40000
5000
c)
0
−5000
−10000
−15000
−40000
Figure 4.21: Example of coupled element mesh. a) Coupled BEM and FEM mesh. b) FEM mesh
alone. c) BEM mesh alone. Length in inches.
Figure 4.21 corresponds to the mesh used to generate the results in gures 4.22 and 4.23, except that the
number of nodes and elements in gure 4.21 have been reduced for visual clearness.
4.5. Results from Coupled BEM/FEM Implementation
67
2
0
−2
inches
B
C
−4
−6
−8
−10
Figure 4.22:
BEM/FEM
BEM
Antes
0
0.1
0.2
0.3
0.4
0.5
s
0.6
0.7
0.8
0.9
1
u1 found from BEM/FEM analysis compared with the results obtained by Antes.
(x1 ; x2 )=(6000 in, 0 in); C : (x1 ; x2)=(12000 in, 0 in).
B
:
5
0
C
A
−5
B
−10
inches
−15
−20
−25
−30
−35
−40
BEM/FEM
BEM
Analytical solution
−45
−50
0
0.1
0.2
0.3
0.4
0.5
s
0.6
0.7
0.8
0.9
1
Figure 4.23: u2 found from BEM/FEM analysis compared with analytical results. A : (x1 ; x2) =(0 in,
0 in); B : (x1 ; x2 )=(6000 in, 0 in); C : (x1 ; x2 )=(12000 in, 0 in).
The results for the coupled BEM/FEM are obtained by the direct BEM. However, comparisons between
the coupled BEM/FEM with direct and indirect BEM show that the dierence between the results is less
than 0.5%, which was expected due to the results from Example 4.6.
In both gures 4.22 and 4.23 the agreement between the results obtained by the BEM and the results
from the coupled BEM/FEM is very good. The results from Antes [5] also agree well with the BEM/FEM
and BEM results for point C . At point C the BEM/FEM and BEM results do not agree completely with
neither the analytical solution nor the results from Antes [5]. As was the case at point C , the analytical
solution at point B shows some dierences from the BEM/FEM and BEM results. The dierence between
68
4. Coupled BEM and FEM Solutions
the numerical results and the analytical results can be caused by the numerical integration performed on
the analytical solution in order to obtain a solution for a Heaviside type load for a distributed load.
Example 4.8
Heaviside Load on Layered Media
The same parameters and discretization as in Example 4.7 are used in this example, but with the Lame
parameters chosen dierently. 3 dierent cases are examined. First, results from the previous example
are presented for comparison
(CASE
1). +Next the Lame parameters in layers 1 and 2 within the alluvial
valley+ are changed
to
+
=
+
=
+
=
2 = 0:1 (CASE 2). Finally the Lame parameters are chosen
1
1
2
as 1 = +1 = 0:1 and +2 = +2 = 0:5 (CASE 3). In all 3 cases the Lame parameters for the
elastic halfspace are chosen as = = 106 ps/in, see gure 4.20.
Notice that the u2 results are presented upside down. This has been chosen due to visuality.
Figure 4.24:
u1
for surface points, CASE 1, BEM/FEM analysis.
Thick lines indicate the edges of the alluvial valley.
+
1
= +1 = +2 = +2 = = .
Figure 4.25:
u2
+
1
= +1 = +2 = +2 = = .
for surface points, CASE 1, BEM/FEM analysis.
Thick lines indicate the edges of the alluvial valley.
The surface displacements from CASE 1 shown in gures 4.24 and 4.25, visualise the initial development
of the Rayleigh surface waves which propagate from the centre of the mesh, where the load is applied, to
the outer boundaries of the mesh. It is clearly seen in both gures 4.24 and 4.25 that no reection arises.
4.5. Results from Coupled BEM/FEM Implementation
69
Notice especially that the displacements around the alluvial valley boundaries, where the FEM and BEM
meshes intersect, are smooth and without reection. From approximately 2 s and forward the waves have
left the mesh and only the nontransient displacements from the Heaviside type load remain.
Figure
4.26: u for points on the boundary between layers 1 and 2, CASE 1, BEM/FEM analysis.
+ + + 1+
1
= 1 = 2 = 2 = = .
Figure
4.27: u for points on the boundary between layers 1 and 2, CASE 1, BEM/FEM analysis.
+ + + 2+
1
= 1 = 2 = 2 = = .
Figures 4.26 and 4.27 show displacements at the articial boundary within the alluvial valley. Since the
entire elastic halfspace and the alluvial valley in this case consist of the same material the displacements
correspond to the halfspace solution. The main displacements arise due to the Rayleigh surface waves
which decreases exponentially from the surface and downwards. The horizontal displacements from the
Rayleigh surface waves teoretically change sign at a depth of approximately 0.2 times the Rayleigh wave
length. When comparing gures 4.26 and 4.24 it is noticed that this sign change has already taken
place at the boundary between layers 1 and 2. Similary, it is seen that the amplitudes in gure 4.27 are
approximately 60% of the amplitudes in gure 4.25, which corresponds to approximately 0.5 times the
Rayleigh wave length. The exact wave length is diÆcult to establish since the Heaviside load generate a
multitude of Rayleigh waves with dierent wave periods, see Example 4.3. Furthermore, the displacements
do not solely arise due to Rayleigh waves but also the P- and S-waves contribute.
70
4. Coupled BEM and FEM Solutions
Figure 4.28: u1 for surface points, CASE 2, BEM/FEM analysis. +1 = +1 = +2 = +2 = 0:1 and
= . Thick lines indicate the edges of the alluvial valley.
Figure 4.29:
and for surface points, CASE 2, BEM/FEM analysis.
= . Thick lines indicate the edges of the alluvial valley.
u2
+
1
= +1 = +2 = +2 == 0:1 Figures 4.28 and 4.29 show the surface displacements from CASE 2. Comparison of these results with
gures 4.24 and 4.25 from CASE 1 shows that in CASE 2 the main displacements are conned within the
alluvial valley. Since the Lame parameters of the alluvial valley are reduced by 90% in proportion to those
of the halfspace the boundaries of the alluvial valley reect most of the wave energy which is trapped
within the valley. In gures 4.28 and 4.29, but mainly in gure 4.29, the reected wave components
can be identied. In both gures 4.28 and 4.29 the intersection of the Rayleigh surface wave with the
boundary can be seen at approximately 2.7 s. Shortly hereafter the rst reection of the P-wave give the
small decline at the centre of the mesh at 3 s. As the P-wave reection moves from the centre again the
rst S- and Rayleigh wave reections arive at the centre at approximately 5 s and 5.5 s. Finally another
P-wave reection arrive after about 6-7 s.
4.5. Results from Coupled BEM/FEM Implementation
71
Figure
4.30: u for points on the boundary between layers 1 and 2, CASE 2, BEM/FEM analysis.
+ + + 1+
1
= 1 = 2 = 2 == 0:1 and = .
Figure
4.31: u for points on the boundary between layers 1 and 2, CASE 2, BEM/FEM analysis.
+ + + 2+
1
= 1 = 2 = 2 = 0:1 and = .
Displacements at the articial boundary between the two layers within the alluvial valley show the same
standing waves as for the surface displacements, see gures 4.30 and 4.31. The amplitudes in gure 4.31
are approximately half of the amplitudes in gure 4.29, which corresponds to a depth of 0.7 times the
Rayleigh wave length. At this depth the amplitudes of the Rayleigh waves in the horizontal direction are
approximately 0.2 times the amplitudes at the surface. Therefore the horizontal displacements in gure
4.30 are mainly caused by the displacements from the P- and S- waves and their reections.
72
4. Coupled BEM and FEM Solutions
Figure 4.32: u1 for surface points, CASE 3, BEM/FEM analysis. +1 = +1 = 0:1 , +2 = +2 = 0:5
and = . Thick lines indicate the edges of the alluvial valley.
Figure 4.33: u2 for surface points, CASE 3, BEM/FEM analysis. +1 = +1 = 0:1 , +2 = +2 = 0:5
and = . Thick lines indicate the edges of the alluvial valley.
Figures 4.32 and 4.33 show the components of the surface displacements from CASE 3. In this case the
Lame parameters of the lower layer in the alluvial valley have only been reduced by 50% in proportion
to the elastic halfspace values. The standing waves within the alluvial valley, which were observed in
CASE 2, are seen again. As in the previous case most of the wave energy is trapped within the soft layer,
but since layer 2 is not as sti as the elastic halfspace a relatively larger part of the wave energy will
propagate into this layer and is not reected. Therefore, identication of the wave components and their
reections is more diÆcult, but at approximately 3 s in gure 4.33 the rst P-wave reection from the
alluvial valley boundary is noticed. The P-wave reection from the boundary between the 2 layers should
arive about 1.3 s and since this is not seen it indicates that reection of the P-waves is not signicant
at the layer boundary. At approximately 1.3 s in gure 4.32 something is noticed, which is probably the
reection from the layer boundary of the S-waves. In gure 4.33 the intersection with and reection from
the alluvial valley boundaries for the Rayleigh surface waves are seen at 2.7 s and 5.5 s, respectively.
Finally the second P-wave reections from the alluvial valley boundaries are seen at approximately 6.5 s.
4.6. Non-linear FEM Formulation
73
Figure
4.34:
+ +
1
u1 for points on the boundary between layers 1 and 2, CASE 3, BEM/FEM analysis.
= 1 = 0:1 , +2 = +2 = 0:5 and = .
Figure
4.35:
+ +
1
u2 for points on the boundary between layers 1 and 2, CASE 3, BEM/FEM analysis.
= 1 = 0:1 , +2 = +2 = 0:5 and = .
Again, the displacements at the internal boundary of the alluvial valley show the characteristics of the
surface displacements. But in this case the amplitudes in gures 4.34 and 4.35 are signicantly reduced
compared to gures 4.30 and 4.31 due to the increased stiness of the lower layer. Therefore the Rayleigh
wave components become clearer in this case, especially in gure 4.34 where they can be identied at
approximately 2.7 s and 5.2 s.
4.6
Non-linear FEM Formulation
In Section 4.2 it was mentioned that K+(tL) requires special attention for non-linear problems. In this section the method in case of non-linear material properties for the FEM
domain will be outlined.
Notice that the superscript + signifying that the quantities belonging to the alluvial valley
will be omitted in this section.
(4.6) gives the basic time domain FEM formulation
Mu + Cu_ + K(tL)u = f(tL )
(4.47)
If the material which is discretized by the FEM, is non-linear the stiness matrix K is
no longer constant, but depends on the actual loads and displacements and is exchanged
74
4. Coupled BEM and FEM Solutions
with the secant stiness matrix at the time tL , Kt(tL). Thus the non-linear time domain
can be expressed by the FEM formulation
Mu + Cu_ + Kt(tL)u = f(tL )
(4.48)
As for the linear case (4.48) is solved by the Newmark iteration scheme. But instead of
solving directly for u (see Appendix D) the solution is found by a Newton-Raphson iteration scheme, i.e. rst a solution using the Newmark iteration scheme for the rst correction to u, Æu(1) is found using Kt;(0) (tL )= Kt(tL 1 ) and f (0) (tL )=f(tL 1 ). Hereafter
the i'th correction to u, Æu(i) is found using Kt;(i 1) (tL ) and f (i 1) (tL). u is found
as
u(i) = u(i 1) + Æu(i)
(4.49)
The iterations continue until Æu(i) fulls a convergence criterion, e.g.
(i) Æ u ju(tL )j
(i) (i)
Æ u = max(Æ u )
(4.50)
ju(tL)j = max(u(tL))
where is some specic fraction.
Notice that Kt;(i) (tL) depends on the current stresses and strains and therefore must
be recalculated for each iteration. Alternatively to this approach the Modied NewtonRaphson scheme can be used. In this scheme Kt;(i) (tL) is replaced by Kt(tL 1). The
obvious advantage here is that Kt(tL 1) does not need recalculation for each iteration but
only once for each time step. The procedure outlined here is performed for each time step
and the time stepping procedure then follows exactly as in Section 4.2.
The secant stiness matrix Kt is calculated as usual for the FEM
Z
t
K = S BT DepBdV ; B = r~ N
(4.51)
where N is the element shape functions and r~ is the generalized dierential operator.
Dep is the tangential elasticity matrix also called the elasto plastic constitutive matrix.
Dep = dd
(4.52)
(4.52) is known as the incremental linear elasto plastic constitutive relation. It is assumed
that the stresses and strains can be separated in an elastic part and a plastic part
= e + p
= e + p
(4.53)
The elastic part is described by the usual elastic constitutive relation
e = De
(4.54)
4.7. Concluding Remarks
75
The plastic part is described by the yield surface f = f (e; ), the ow rule given by the
plastic potential g = g(e; ) and the hardening parameter . If f = g the material is
known as an associated plastic material.
Using (4.53) and (4.54) it can be shown that the elasto plastic constitutive matrix can be
found as
Dep = D D
"
@g
@
#"
@f
@
#T
2
D 4H +
"
@f
@
#T
D
"
@g
@
#3
5
1
(4.55)
where for `work' hardening
H=
@f @g
@ @
(4.56)
If the material is assumed to be ideally elasto plastic, H = 0.
For a specic elasto plastic model the dierentials in (4.55) and (4.56) need to be evaluated.
4.7
Concluding Remarks
In the present chapter coupled BEM and FEM schemes in the time domain for both the
direct and indirect BEM formulations are presented. The schemes are used in 2D formulations. The 2D BEM formulation is formulated using the principles from Chapter 2 and
the numerical integration is done by the new integration principle described in Section
3.3. The singularity problem in 2D is solved by the ridgid body motion principle with
enclosing elements. The numerical procedures are tested and evaluated in Example 4.1
and Example 4.2, showing that the procedures work very well and give accurate results
even with few subdivisions.
Prior to these examples the BEM and FEM implementations are tested separately in
Example 4.3 and Example 4.4. The FEM results very clearly show why absorbing boundaries are necessary and furthermore they give an indication of the time step which would
be appropriate for the problem. The 2D BEM implementation gives very accurate results
compared to an analytical solution and is further tested for another problem in Example
4.5. Both the direct and indirect 2D BEM are implemented, and in Example 4.6 results
from the two formulations are compared and show that the dierences are negligible.
The coupled formulation is applied to two problems. First the problem of a Heaviside
load on an elastic halfspace is applied (Example 4.7), and hereafter the same problem
but for a layered medium with dierent ratios of the material parameters is evaluated
(Example 4.8). The rst example shows good agreement with both the 2D BEM results
and the analytical solution and the second problem shows some of the advantages of a
coupled formulation.
Finally it is shown how the coupled formulation can be used in connection with a nonlinear FEM formulation.
76
4. Coupled BEM and FEM Solutions
Chapter 5
Boundary Element Solution for a
Moving Force
5.1
Introduction
TraÆc induced vibrations are a known source of discomfort and damage problems in
nearby buildings, either directly in form of vibrations of the buildings or through structure
born acoustic noise, Okumura & Kuno [57], Trochides [77]. During the last couple of
decades the general amount of traÆc has increased considerable. The increased road traÆc
along with heavier trucks have emphazized the problems of traÆc induced vibrations. But
also train induced vibrations have been given some attention in recent years, mainly due
to the dramatic increase in speed, Krylov [42]. The increased speed and the still greater
number of underground railway systems in urban areas have resulted in increasingly larger
vibration problems. An considerable eort has been done to minimize or soften the
eects of the trac induced vibrations, Chouw & Schmid [15], Nelson [53]. With respect
to analysis, semi-empirical methods for estimating the vibrations have been introduced,
Madshus et al. [46], and some numerical methods used to estimate the vibrations have
been introduced, but mainly in the frequency domain.
Simulation of traÆc induced vibrations in an elastic media can be achieved as a sum of
moving concentrated forces due to the superposition principle. Thus, the basic problem
that must be solved to simulate traÆc induced vibrations is the problem of a moving
time dependent concentrated force on the surface of an elastic half-space. When using
the BEM (or FEM) to solve this problem an important problem arises. With an element
mesh of nite size the moving force will soon move beyond the boundary of the mesh.
Therefore it will be preferable to formulate the problem in a coordinate system moving
along with the force.
FEM formulations in convected coordinates where the mesh is following the force are well
known, see e.g. Zienkiewicz & Taylor [87]. If this is coupled with transmitting boundaries
an eÆcient method for modelling transient problems of moving forces on innite media
arises, see e.g. Krenk et al. [39].
In this chapter Green's functions for the time domain in 2D are formulated in a coordinate
system moving along with the force and thereby a BEM formulation with a moving
coordinate system is established similar to the method of coupling convected coordinates
77
78
5. Boundary Element Solution for a Moving Force
with absorbing boundaries for the FEM.
5.2
Formulation of Boundary Element Equations for
a Moving Force
The displacement eld Uik (x; t; ) in coordinate direction i from a concentrated time
dependent unit force f ( ) applied at position in direction k in a homogeneous, isotropic
and linear elastic media is given by (2.1)
Uik (x; t; ) =
Z
t
0
gik (x; t; ; )f ( )d
(5.1)
where Green's function for the displacement eld gik (x; t; ; ) is the solution to the differential equations
@ijk
@xj
@ 2 gik
+ Æik Æ(x
@t2
)Æ(t
) = 0
(5.2)
ijk (x; t; ; ) is the stress
tensor derived from the displacement eld gik (x; t; ; ). Further the summation convention for Cartesian tensors has been used. Next consider a
concentrated unit force moving with a constant velocity vi in a linear elastic media, see
gure 5.1.
Figure 5.1: Moving force in xed and moving coordinate systems.
Next a (~x1 ; x~2 ; x~3)-coordinate system is introduced, which follows the moving force. ~i
signies the coordinates of the force, and g~ik (~x; t; ~; ) is the Green's function for the
displacement eld u~i(~x; t) in the moving coordinate system. The two coordinate systems
coalesce at the time t = which corresponds to the Galilean transformation
xi = x~i + vi (t )
(5.3)
The equations of motion in the moving coordinate system are expressed by introducing
the partial dierentiation operators
9
@ = @
>
>
@xj
@ x~j >
>
=
@
@ = @ v
j
(5.4)
@t x j
@t x~j
@ x~j
>
>
>
@
@
2vj @[email protected] @t x~j + vj vk @[email protected] @x~k >;
@t xj = @t x~j
2
2
2
2
2
2
5.3. Time Approximation
79
Hereby the partial dierential equations for the Green's functions in the moving coordinate
system reads
@ ~ijk
@ x~j
0
@ 2 g~ @ 2ik @t x~j
1
2
@ 2 g~ik 2vj @ x~ @t + vj vl @@x~ [email protected]~ A + Æik Æ(~x
j
j l
x~j
~)Æ(t
) = 0
(5.5)
Introducing the coordinate transformation (5.3) in (5.5) gives the following equations for
g~ij as a function of the xed coordinates xi
@ ~ijk
@xj
@ 2 g~ 2ik + Æik Æ (x
@t x~j
v(t
)
~)Æ(t
) = 0
(5.6)
Since Æ(x v(t ) ~)Æ(t ) = Æ(x ~)Æ(t ) the solution of (5.6) becomes
gij (x; t; ~; ). The solution to (5.5) can then be expressed in terms of the Green's function
for the displacement eld in the xed coordinate system in the following way
g~ik (~x; t; ~; ) = gik (x; t ; ~; 0) = gik (~x + v(t ); t ; ~; 0)
(5.7)
The Green's function for Cauchys stress tensor eld ~ijk (~x; t; ~; ) is obtained upon insertion of the solution of (5.7) into the constitutive equations for the considered linear
elastic material. Since these only involves spatial dierential operations with respect to
x, the result is seen to be
~ijk (~x; t; ~; ) = ijk (~x + v(t ); t ; ~; 0)
(5.8)
Hence a BEM formulation in the moving coordinate system, that calculates the displacement eld following the force, can be established by using the formal Green's functions
g~ik (~x; t; ~; ) and ~ijk (~x; t; ~; ) as given by (5.7) and (5.8) instead of the original ones.
Notice that this result is valid generally even for non-homogeneous or anisotropic linear
materials.
5.3
Time Approximation
As mentioned implementation of the BEM formulation in the moving coordinate system
is done similar to the procedure used for the non moving BEM formulations described in
Section 4.3.
The space integration is done by the numerical integration scheme introduced in Section
4.3.2 and the time discretization is done by analytical integration as in Section 4.3.5.
But in this case the Green's functions are exchanged with the formal Green's functions
g~ik (~x; t; ~; ) and ~ijk (~x; t; ~; ).
The formal Green's functions in a moving coordinate system in 2D are found from (4.28)
and (4.29) and can be written as
80
5. Boundary Element Solution for a Moving Force
1
1 H (c t r~(t)) hg~(1)(t; r~(t); c ) Æ g~(2) (t; r~(t); c )i +
~
g~ (~x; t; ; 0) =
1
1
2 c1 1
1 H (c t r~(t)) h g~(1) (t; r~(t); c ) + Æ g~(2) (t; r~(t); c ) + Æ g~(3) (t; r~(t); c )i (5.9)
2
2
2
c2 2
; = 1; 2
where
2c2t2 r~(t)2 r~ (t)~r (t)
(1) (t; r~(t); c) = q
g~
r~(t)4
c2 t2 r~(t)2
q
1
(2)
g~ (t; r~(t); c) =
c2 t2 r~(t)2
r~(t)2
(3) (t; r~(t); c) = q 1
g~
c2 t2 r~(t)2
!
(5.10)
1 H (c t r~(t)) ht~(1) (t; r~(t); c ) + t~(2) (t; r~(t); c )i
~
t~ (~x; t; ; 0) =
1
1
2 c1 1
1 H (c t r~(t)) ht~(3) (t; r~(t); c ) + t~(2) (t; r~(t); c )i ; ; = 1; 2
2
2
2
c2
(5.11)
3=2
!!
r~ (t)~r (t) r~1 (t)n1 + r~2 (t)n2
r~ (t)
1 < ct !2 1=
n
+ 2 r~(t)2
;
r~(t)2 : r~(t)
r~(t)
r~(t)
2
!
ct
2
1
2
r
~
(
t
)
r
~
(
t
)
r
~
(
t
)~
r
(
t
)
r
~
(
t
)
(2)
r
t~ (t; r~(t); c) =
n
+ n r~(t) + Æ 4 r~(t)2 r~(t)2 ct 2 1 r~(t)
r~(t)
!!
r~1 (t)n1 + r~2 (t)n2
(5.12)
r~(t)
8
9 3=2
!2
!
!
<
=
1
r
~
(
t
)
n
+
r
~
(
t
)
n
ct
r
~
(
t
)~
r
(
t
)
1
1
2
2
(3)
t~ (t; r~(t); c) =
1;
2 r~(t)2 Æ
r~(t)2 : r~(t)
r~(t)
!
r~ (t)
n r~(t)
t~(1) (t; r~(t); c) =
8
9
r is given as
q
r~(t) = r~12 + r~22 ;
r~1 (t) = x~1 + v1 t ~1 ;
r~2 (t) = x~2
~2
(5.13)
5.3. Time Approximation
81
where it has been assumed that the force moves along the x1 -axis. The following integrals
need to be evaluated analytically
Z t
Z
(2) = t g~ (~x; ; ~; 0) d
~
G(1)
=
g
~
(~
x
;
;
;
0)(1
)d
G
(5.14)
l
l
1
tl
1
tl
1
1
where l is dened in gure 2.5. Similar for Green's function for the traction eld for a
moving force.
In both (5.9) and (5.11) it is seen that integration of H (ct r~(t)) times the functions in
(5.10) and (5.12) is necessary. Integration of H (ct r~(t)) times a function of t can be
done as follows
8 Rt
;
h(0) tl 1
>
Z t
< Rtl f ( )d
t
f ( )H (c r~( ))d = > h(0) f ( )d ; t1 h(0) tl 1
(5.15)
tl
:
0
;
h(0) tl 1
where f ( ) is an abitrary function of .
h(0) is found as
q
2(~1 x~1 )v1 + D(0)
h(0) =
(5.16)
2(v12 c2 )
D(0) = 4(~x1 ~1 )2 v12 4(v12 c2 )~r(0)2
(5.17)
With (5.15) the integration of (5.9) and (5.11) is strait forward and is done by integration
of each function in (5.10) and (5.12) separately. Due to the complexity of the results only
results for x~2 = ~2 = r~2 = 0 are presented. The results can be seen in Appendix F.
1
1
1
1
1
Example 5.1
Results from 2D BEM Solution for Moving Force
The problem of a distributed constant downward directed Heaviside type uniform line load moving along
the x~1 -axis with the constant velocity v1 on the surface of an isotropic, homogeneous, linear elastic 2D
plane strain halfspace is rst analysed. The problem is similar to the problem considered in Example
4.3-4.5 except that the load in this case is moving.
The distributed moving load pr. unit length p(~x; t) is dened as
p(~x; t) = p0 [H (~x1 + b) H (~x1 b)]H (t)
(5.18)
The analysis is performed with the following parameters : Intensity of line load
p0 = 104 ps/in (6:896 104
2
6
7
KN/m), b =3 3000 in (1.146.4 m), E = 2:5 10 ps/in (1:724 10 KN/m ), = 2:891 10 3 ps2 =in4
(3:151 kg/m ) and = 0:25 (modied2 for 2D plane strain). With assumption of plane strain this gives
= = 106 ps/in (6:895 106 KN/m ) and the wave velocities for the P-, S- and Rayleigh waves respectively c1 =2 3:27 104 ips (8:306 102 m/s), c2 = 1:86 104 ips (4:742 102 m/s) and cR = 1:71 104 ips
(4:360 10 m/s).
The elements are 2 node line elements with linear interpolation, see Appendix B.1.1.
The problem is rst analysed for v1 = 0:01 cR. Dierent element meshes have been tested and the results
have shown that h=1000 in and t=0.01571 s give results which are satisfactory, similar to the results
for the problem for a stationary load, see Section 4.4. The results for h=1000 in and t=0.01571 s
82
5. Boundary Element Solution for a Moving Force
are compared to the stationary results from Section 4.4 in gures 5.2 and 5.3.
2
1
0
inches
−1
−2
B
−3
C
−4
−5
−6
−7
0
0.1
0.2
0.3
0.4
0.5
s
0.6
0.7
Figure 5.2: u1 from Heaviside load. v1 = 0:01 cR and v1 = 0
(~x1 ; x~2 )=(12000 in, 0 in).
0.8
0.9
1
0.9
1
. B : (~x1 ; x~2)=(6000 in, 0 in); C :
5
0
C
A
−5
B
−10
inches
−15
−20
−25
−30
−35
−40
−45
Figure 5.3:
0
0.1
0.2
0.3
0.4
0.5
s
0.6
0.7
u2 from Heaviside load. v1 = 0:01 cR and v1 = 0
(~x1 ; x~2 )=(6000 in, 0 in); C : (~x1 ; x~2)=(12000 in, 0 in).
0.8
. A : (~x1 ; x~2 ) =(0 in, 0 in); B :
The velocity of the load have been chosen to v1 = 0:01 cR due to problems with stability of the solution.
When v1 ! 0 the moving force solution becomes instabil due to numerical problems with Green's functions
for a moving force, see the expressions in Appendix F.
In gures 5.2 and 5.3 it is clear that the results obtained with the moving force solution are less smooth
than the results from the stationary solution. This eect will decrease when the velocity is increased and
the numerical instability decreases, see gures 5.4-5.7. Eventhough the moving force solution have been
obtained with a velocity dierent from 0 the agreement with the stationary solution is execellent.
Next the problem is analysed with v1 = 0:2 cR and v1 = 0:5 cR and the results are compared with the
results for v1 = 0:01 cR
5.3. Time Approximation
83
6
5
inches
4
3
B
2
C
1
0
−1
0
0.1
0.2
0.3
0.4
Figure 5.4:
0.5
s
0.6
0.7
0.8
0.9
u1 from Heaviside load. v1 = 0:01 cR , v1 = 0:2 cR
(~x1 ; x~2 )=(-6000 in, 0 in); C : (~x1 ; x~2 )=(-12000 in, 0 in).
1
and v1 = 0:5 cR .
B
:
A
:
5
0
A
−5
B
C
−10
inches
−15
−20
−25
−30
−35
−40
−45
−50
Figure 5.5:
0
0.1
0.2
0.3
0.4
0.5
s
0.6
0.7
0.8
0.9
1
u2 from Heaviside load. v1 = 0:01 cR , v1 = 0:2 cR
and v1 = 0:5 cR .
(~x1 ; x~2 ) =(0 in, 0 in); B : (~x1 ; x~2 )=(-6000 in, 0 in); C : (~x1 ; x~2)=(-12000 in, 0 in).
In gures 5.4 and 5.5 the results `behind the force' are shown. It is seen as expected that the results
become smoother with increasing velocity. The wave fronts arrive earlier with increasing velocity. This is
due to the motion of the force. The wave velocity observed at a point following the force is the velocity
of the force subtracted from the wave velocity (c v) and therefore prior arrival of the of the wave fronts
are observed in gures 5.4 and 5.5. For a point `in front of the force' the observed wave velocity is the
velocity of the force added to the wave velocity (c + v) and later arrival of the wave fronts are expected,
see gures 5.6 and 5.7.
84
5. Boundary Element Solution for a Moving Force
2
0
−2
inches
B
C
−4
−6
−8
−10
0
0.1
0.2
0.3
0.4
Figure 5.6:
0.5
s
0.6
0.7
0.8
0.9
u1 from Heaviside load. v1 = 0:01 cR , v1 = 0:2 cR
(~x1 ; x~2 )=(6000 in, 0 in); C : (~x1 ; x~2)=(12000 in, 0 in).
1
and v1 = 0:5 cR .
B
:
A
:
5
0
C
A
−5
B
−10
inches
−15
−20
−25
−30
−35
−40
−45
−50
Figure 5.7:
0
0.1
0.2
0.3
0.4
0.5
s
0.6
0.7
0.8
0.9
1
u2 from Heaviside load. v1 = 0:01 cR , v1 = 0:2 cR
and v1 = 0:5 cR .
(~x1 ; x~2 ) =(0 in, 0 in); B : (~x1 ; x~2 )=(6000 in, 0 in); C : (~x1 ; x~2 )=(12000 in, 0 in).
Comparing gures 5.5 and 5.7 reveals that `in front of the force' the vertical displacements increases
whereas `behind the force' the vertical displacements decreases compared to the results with v1 = 0:01 cR.
The force is continuosly moving away from the points behind it and thereby increasing the distance
between the waves which gives a reduction of the displacements. 'In front of the force' on the contrary
the force moves in the same direction as the waves which decreases the distance between the waves which
again increases the displacements. For a stationary oberserver this will look like a small lump of soil
moving in front of the force.
5.3. Time Approximation
Example 5.2
85
Comparison with Corresponding FEM Solutions
The BEM solution for a moving force is next compared with results found by a FEM analysis with
convected coordinates using elastic wave impedance conditions as transmitting boundaries, Krenk et al.
[39]. In this example an isotropic, homogeneous, linear elastic 2D plane strain halfspace isR 1subjected to a
moving time dependent downward directed uniform line load P (t) with zero net impulse 1 P (t)dt = 0.
The time variation of the impulse is given as
P (t) = P0 (1 2 )2
;
1< <1 (5.19)
= 2t=T 1
P0 is the amplitude of the load, which appears almost sinusonidal with period T . However, the load has
vanishing derivatives at t = 0 and t = T .
The example
is performed
with3 the following parameters : Intensity of load P0 = 1 MN, = = 100 106
2
3
N/m and = 2:0 10 kg/m . This gives the wave velocities for the P-, S- and Rayleigh waves respectively c1 = 387 m/s, c2 = 224 m/s and cR = 206 m/s.
160 m of the surface is discretized and the load is applied at the centre. The applied FEM mesh is rectangular with depth 80 m. The BEM mesh consists of 4 m 2 node line elements with linear interpolation
and similarly the FEM mesh consists of 4 m by 4 m, 4 node quadrilateral elements, see gure 5.8. At the
articial boundary of the modelled domain socalled absorbing boundary conditions of the impendance
type are prescreibed,
Wolf & Song [84]. These are unable to support static loadings and loadings with
a net impulse R 11 P (t)dt 6= 0 will cause a permanent sti body oset of the modelled domain. The
defect has been taken into consideration at the specication of the loading of the present problem. Krenk
& Kirkegaard [40] have derived impedance boundary conditions, which are able to absorb plane P- and
S-waves, whereas other types of waves (Rayleigh, Stonely, Love waves and non-plane P- and S-waves)
will all be more or less reected back into the modelled domain at the articial boundary.
Figure 5.8: a) BEM and b) FEM meshes.
The horizontal, u1 and vertical, u2 displacement histories at (~x1 ; x~2) = (+20 m, 0 m), `in front of the
force' and at (~x1 ; x~2) = (-20 m, 0 m), `behind the force' are presented in the following gures.
86
5. Boundary Element Solution for a Moving Force
2
u [mm]
1
1
0
−1
−2
−3
2
u [mm]
1
1
0
−1
−2
−3
0
0.1
0.2
0.3
0.4
0.5
t [s]
0.6
0.7
0.8
0.9
1
Figure 5.9: Horizontal displacements at (~x1 ; x~2 )=(+20 m, 0 m) obtained using BEM (top) and FEM
(bottom) for v1 = 0 c2 (BEM v1 = 0:001 c2) , v1 = 0:2 c2
and v1 = 0:5 c2 .
2
u [mm]
1
2
0
−1
−2
−3
2
0
2
u [mm]
1
−1
−2
−3
0
0.1
0.2
0.3
0.4
0.5
t [s]
0.6
0.7
0.8
0.9
1
Figure 5.10: Vertical displacements at (~x1 ; x~2 )=(+20 m, 0 m) obtained using BEM (top) and FEM
(bottom) for v1 = 0 c2 (BEM v1 = 0:001 c2) , v1 = 0:2 c2
and v1 = 0:5 c2 .
The convection
terms in the FEM scheme tend to destabilize the numerical results at large Mach-numbers
(M = vc21 ). This was cured by a kind of upwind streamline nite dierence scheme, Krenk et al. [39].
The indicated stabilization is tantamount to the introduction of articial numerical damping, which may
change the magnitude of the observed response slightly. This defect should be kept in mind in the
following comparison with the BEM solution.
Comparison between gures 5.9 and 5.10 shows that both the minimum and maximum displacements
decreases with greater speed. However, for v1 = 0 c2 in gure 5.9 the BEM results for the minimum
displacements are signicantly smaller than the equivalent minima from the FEM results. As mentioned
in the previous example the solution of the time integrals of the Green's functions become instable when
v1 approaches null. Therefore, the results for very small velocities obtained by the BEM are not reliable.
5.3. Time Approximation
87
2
u1 [mm]
1
0
−1
−2
−3
2
u1 [mm]
1
0
−1
−2
−3
0
0.1
0.2
0.3
0.4
0.5
t [s]
0.6
0.7
0.8
0.9
1
0.8
0.9
1
Figure 5.11: Horizontal displacements at (~x1 ; x~2 )=(-20 m, 0 m) obtained using BEM (top) and FEM
(bottom) for v1 = 0 c2 , v1 = 0:2 c2
and v1 = 0:5 c2 .
2
u2 [mm]
1
0
−1
−2
−3
2
u2 [mm]
1
0
−1
−2
−3
0
0.1
0.2
0.3
0.4
0.5
t [s]
0.6
0.7
Figure 5.12: Vertical displacements at (~x1 ; x~2 )=(-20 m, 0 m) obtained using BEM (top) and FEM
(bottom) for v1 = 0 c2 , v1 = 0:2 c2
and v1 = 0:5 c2 .
In all the gures 5.9-5.12, but especially in gure 5.12, the FEM results show tendencies to reect Rayleigh
waves from the boundaries. However, the BEM results do not show any signs of reection in any of the
gures. The calculation time for the FEM solution was 62 seconds on a 200 MHz PC while the calculation
time for the BEM solution was several hours on a SiliconGraphics Indy R4400.
88
5. Boundary Element Solution for a Moving Force
In gures 5.13-5.14 the FEM and BEM results for velocities approaching the S-wave velocity are shown
'in front of the force'.
2
u [mm]
1
1
0
−1
−2
−3
2
u [mm]
1
1
0
−1
−2
−3
0
0.1
0.2
0.3
0.4
0.5
t [s]
0.6
0.7
0.8
0.9
1
0.9
1
Figure 5.13: Horizontal displacements at (~x1 ; x~2 )=(+4 m, 0 m) obtained using BEM (top) and FEM
(bottom) for v1 = 0:5 c2 , v1 = 0:7 c2
and v1 = 0:9 c2 .
2
0
−1
2
u [mm]
1
−2
−3
−4
2
0
−1
2
u [mm]
1
−2
−3
−4
0
0.1
0.2
0.3
0.4
0.5
t [s]
0.6
0.7
0.8
Figure 5.14: Vertical displacements at (~x1 ; x~2 )=(+4 m, 0 m) obtained using BEM (top) and FEM
(bottom) for v1 = 0:5 c2 , v1 = 0:7 c2
and v1 = 0:9 c2 .
The Rayleigh wave velocity for Poisson materials ( = ) is approximately 0:92 c2. The results for
velocities approaching the Rayleigh wave velocity are shown for a point only 4 m in front of force. The
reason for this change is that the wavefront 'in front of the force' at the surface do not move very far
from the force when the force moves almost as fast as the Rayleigh wave. Krylov [42] assumes that
the displacements increases dramatically when the velocity approach the Rayleigh wave velocity. The
results in gures 5.13-5.14 obviously disagrees with that statement. Actually, the displacements decreases
compared to the results obtained with lower velocities shown in gures 5.9-5.12. This result could be
anticipated since the load is a limited function and the energi within the system must hereby also be
limited.
5.4. Concluding Remarks
5.4
89
Concluding Remarks
A new BEM formulation in a moving coordinate system is given. The formulation introduces Green's functions for a moving force and shows how the normal BEM formulation
for a xed coordinate system can be modied to a formulation in the moving coordinate
system by using the Green's functions for a moving force.
The BEM results are compared to the results obtained by the stationary BEM for the
problem of a Heaviside type load on the surface of an elastic halfspace. The results shows
that for small velocities the moving force solution becomes instable due to numerical
problems with Green's functions for a moving force.
The problem of a moving time dependent force of zero net impulse have next been tested
with the method and the results compared to results obtained by a FEM formulation
using an impedance boundary condition as absorbing boundary. The results shows good
agreement and it is demonstrated that the reection of Rayleigh waves is omitted in the
BEM analysis. Further it is shown that the displacements do not increase signicantly
when the velocities approaches the Rayleigh wave velocity.
The problem has been formulated to circumvent some of the shortcomings of the FEM
modelling, making a fair comparison possible. The results show good agreement. As
expected the partial reection of various wave types such as Rayleigh wave propagation,
which are accounted for in the considered absorption conditions, are not present in the
BEM analysis. This result has been achieved at the expense of much larger computational
eorts.
90
5. Boundary Element Solution for a Moving Force
Chapter 6
Conclusions
The present chapter is divided into two parts. Section 6.1 includes a summarized description of the contents of the thesis and some conclusions. General conclusions and a
description of the new contributions provided by this thesis are given in Section 6.2.
6.1
Summary of Thesis
Chapter 1 contains the introduction to the present thesis. Various
boundary element formulations used to solve stress wave propagation problems in soils
are briey described and relevant references are listed. The scope of the thesis is dened
and a brief outline of the coming chapters is given.
Simultaneous derivations of the direct and indirect time domain BEM
are given in this chapter. It is shown how both methods represent ways of discretising
the Somigliana identity and how the same numerical methods apply to both formulations.
Both methods are formulated for simple examples involving stress wave propagation in
soils.
In this chapter the importance of using accurate numerical integration in
the BEM is shown. The problem of a 3D elastic halfspace with a at surface subjected to
a concentrated load is formulated using polar elements. The formulation is tested using
a Heaviside load and the results show signs of instability which is found to arise from
inaccurate numerical integration. A new integration scheme using a wave front based integration is suggested and implemented. The results show signicant improvements and
the conclusion is that a wave front based numerical integration should be used to ensure
accurate results for 3D BEM formulations in the time domain. An analytical solution for
the singularity problem has been applied to Green's functions for the displacement eld.
Chapter 1.
Chapter 2.
Chapter 3.
91
92
6. Conclusions
Chapter 4 shows how the advantages of the FEM and BEM can be
combined into a coupled method. A coupling scheme is presented for both the direct and
indirect BEM. Both schemes are formulated in 2D time domain using the principles from
Chapter 2. In Chapter 4 extensive tests of the wave front based integrations scheme are
shown and the results show that the scheme is working well and producing accurate results
even with few subdivisions. The singularity problem for Green's function for the traction
eld in the 2D time domain BEM has been solved by the ridgid body motion principle.
Prior to examples of the coupled methods an example is applied to the FEM and BEM
schemes separately, which clearly shows why absorbing boundaries or coupled methods
are necessary for the FEM. The 2D BEM implementations of both the direct and indirect
BEM show good agreement with analytical results and it is shown that similar direct and
indirect BEM implementations produce almost identical results. The coupled schemes are
tested in two examples : The problem of a Heaviside load on an elastic halfspace and the
same problem but with dierent ratios between the material parameters in the interior
(FEM) and exterior (BEM) domains. Finally it is briey shown how the coupled schemes
can be used with a non-linear FEM formulation.
A formulation of the time domain BEM for a moving force is presented
in Chapter 5. Green's functions for a moving force are introduced and the method is
implemented. The method uses the BEM formulation shown in Chapters 2 and 4, but
with Green's functions modied for a moving force. Results of the method are compared
to stationary results obtained earlier in the thesis and to results obtained by an FEM
formulation in convected coordinates using absorbing boundary conditions. The results
show that the moving force formulation with the BEM becomes instable when the velocity of the force approaches zero due to numerical problems with Green's function for
a moving force. The FEM formulation reveals partial reection of the Rayleigh waves,
which is not present for the BEM results due to its inherent satisfaction of the radiation
conditions. Further, it is shown that when the velocity of the load is approaching the
Rayleigh wave velocity the displacements do not increase signicantly.
The appendices address various analytical and numerical aspects. In
Appendix A the product of two normal vectors integrated over a hemisphere is derived
analytically to be used in Chapter 2. In Appendix B the shape functions and Gauss
quadratures used in the BEM and FEM implementations are shown. Appendix C contains
the rules by which the sparse formulation used in the BEM and FEM implementations
are dened and a short description of the sparse solver is given. The time domain FEM
implementations use the Newmark integration scheme, which is described in Appendix
D. In Appendix E the analytical solutions from the literature, which have been used for
comparison with the numerical results, are presented. Appendix F contains the analytical
expressions of Green's functions for a moving force after time integration.
Chapter 4.
Chapter 5.
Appendices.
6.2. Overall Conclusions
6.2
93
Overall Conclusions
The conclusions of the present thesis can be summarized in the following main points
1. The accuracy of the BEM is superior to most other numerical methods when dealing
with stress wave propagation problems in soil due to its inherent ability to satisfy
the radiation conditions exactly.
2. Simultaneous derivations of the indirect and direct time domain BEM reveal that
both methods represent numerical discretizations of the Somigliana identity and
hereby that the same numerical techniques apply to both methods.
3. For stress wave propagation problems in soil in 2D or 3D, a wave front based numerical integration scheme is a way of ensuring acceptable accuracy of the numerical
integration, which is crucial to the accuracy of the method.
4. For 3D time domain problems an analytical solution of the singularity problem for
Green's function for the displacement eld has been derived and the results for a
simple 3D problem proved the accuracy.
5. 2D time domain formulations of the direct and indirect BEM give almost identical
results when implemented using the same numerical techniques. This further substantiates the fact that both formulations represent numerical discretizations of the
Somigliana identity.
6. Green's functions for a moving force have been indicated and the obtained results are
compared with results from a similar FEM formulation. It is shown that the BEM
does not give reection of the Rayleigh waves, whereas the applied FEM formulation
does. This result is achieved at the expense of a much larger computational time.
New contributions in this thesis are, to the author's knowledge, the following:
1. Simultaneous derivations of the direct and indirect time domain BEM and implementation of those in 2D with the same numerical techniques.
2. Analytical solution to the singularity problem for Green's function for the displacement eld in 3D time domain BEM.
94
6. Conclusions
3. Formulation of 3D time domain BEM in polar coordinates for a concentrated load
on an elastic halfspace.
4. Simultaneous formulation of coupled BEM/FEM for direct and indirect 2D time
domain BEM.
5. Derivation and implementation of wave front based numerical integration for both
2D and 3D time domain BEM.
6. Formulation of Green's functions and their implementation in the 2D time domain
BEM for a moving force.
Bibliography
[1] Abe, K. & Yoshida, Y., Seismic Response Analysis of Alluvial Valley with two Dipping Layers by a Boundary Element-Finite Element Coupling Method. International
Journal for Numerical Methods in Engineering. Vol. 35, pp. 1753-1769, 1992.
[2] Abramowitz, M. & Stegun, I.A., Handbook of Mathematical Functions. National Bureau of Standards, USA, Applied Mathematics Serie 55, 1970.
[3] Ahmad, S. & Banerjee, P.K., Multi-Domain BEM for Two-Dimensional Problems of
Elastodynamics. International Journal for Numerical Methods in Engineering. Vol.
26, pp. 891-911, 1988.
[4] Aki, K. & Richards, P.G., Quantative Seismology, Volume I, Theory and Methods.
W.H. Freeman and Company, 1980.
[5] Antes, H., A Boundary Element Procedure for Transient Wave Propagation in TwoDimensional Isotropic Elastic Media. Finite Elements in Analysis and Design. Vol.
1, pp. 313-322, 1985.
[6] Antes, H., Steinfeld B. & Trondle, G., Recent Developments in Dynamic Stress Analysis by Time Domain BEM. Engineering Analysis with Boundary Elements. Vol. 8,
pp. 179-191, 1991.
[7] Argyris, J. & Chan, A.S.L., Applications of Finite Elements in Time and Space.
Ingenieur Archiv. Vol. 41, pp. 235-257, 1972.
[8] Argyris, J. & Mlejrek, H.P., Texts on Computational Mechanics. Volume V. Dynamics of Structures. North-Holland, Amsterdam, 1991.
[9] Argyris, J. & Scharpf, D., Finite Elements in Time and Space. The Aeronautical
Journal of the Royal Aeronautical Society. Vol. 73, pp. 1041-1044, 1969.
[10] Banerjee, P.K., Ahmad, S. & Manolis, G.D., Advanced Elastodynamic Analysis. In
Computational Methods in Mechanics, Volume III. Boundary Element Methods in
Mechanics, Ed. Beskos. D.E., Elsevier Science Publishers B.V., 1987.
[11] Bathe, K.J., Finite Element Procedures in Engineering Analysis. Prentice Hall, 1982.
[12] Belytschko, T. & Lu, Y.Y., A Variationally Coupled FE-BE Method for Transient
Problems. International Journal for Numerical Methods in Engineering. Vol. 37, pp.
91-105, 1994.
95
96
Bibliography
[13] Beskos, D.E., Boundary Element Methods in Dynamic Analysis. Applied Mechanics
Reviews. Vol. 40, No. 1, pp. 1-23, 1987.
[14] Beskos, D.E., Boundary Element Methods in Dynamic Analysis: Part II. Applied
Mechanics Reviews. Vol. 50, No. 3, pp. 149-197, 1997.
[15] Chouw, N., Le, R. & Schmid, G., Reduction of Train-induced Vibrations and Waves.
Ruhr-University Bochum, Germany, 1993.
[16] Cruse, T.A. & Rizzo, F.J., A Direct Formulation and Numerical Solution of the
General Transient Elastodynamic Problem, I. Journal of Matematical Analysis and
Applications. Vol. 22, pp. 244-259, 1968.
[17] Eringen, A.C. & Suhubi, E.S., Elastodynamics, Volume II Linear Theory. Academic
Press, 1975.
[18] Eshraghi, H. & Dravinski, M., Transient Scattering of Elastic Waves by Three Dimensional Non-Axisymmetric Dipping Layers. International Journal for Numerical
Methods in Engineering. Vol. 31, pp. 1009-1026, 1991.
[19] Fletcher, R., in Numerical Analysis, Ed. Dold, A. & Eckmann. B. Springer-Verlag,
Berlin, pp. 73-89, 1975.
[20] Fried, I., Finite Element Analysis of Time-Dependent Phenomena. AIAA Journal.
Vol. 7, pp. 1170-1172, 1969.
[21] Fukuri, T., Time Marching BE-FE Method in 2-D Elastodynamic Problems. International Conference on BEM IX. Stuttgart, 1987.
[22] Fukuri, T. & Ishida, Y., Time Marching BE-FE Method in Wave Problem. In Theory
and Applications of Boundary Element Methods. Ed. Tanaka, M. & Du, Q., pp. 95106, Pergamon Press, Oxford, 1988.
[23] Gakenheimer, D.C. & Miklowitz, J., Transient Excitation of an Elastic Half Space
by a Point Load Traveling on the Surface. Journal of Applied Mechanics. Vol. 36, pp.
505-515, 1969.
[24] Gaul, L., Fiedler, C. & Ricoeur, A., A New Hybrid Symmetric Boundary Element
Method in Elastodynamics. In Computational Mechanics '95. Ed. Atluri, S.N., Yagawa, G. & Cruse, T.A., Springer-Verlag, Berlin, Vol. 2, pp. 3074-3079, 1995.
[25] Giuggiani, M., Computing Principal Value Integrals in 3D BEM for Time-Harmonic
Elastodynamic - a Direct Approach. Communication in Applied Numerical Mechanics. Vol. 8, pp. 141-149, 1992.
[26] Golub, G.H. & Van Loan, C.F., Matrix Computations 2nd ed. Johns Hopkins University Press, Baltimore, 1989.
Bibliography
97
[27] Gomez-Lera, M.S. & Alarcon, E., Elastostatics. In Computational Methods in Mechanics, Volume III. Boundary Element Methods in Mechanics. Ed. Beskos, D.E.
Elsevier Science Publishers B.V. 1987.
[28] Groen, A.E., Three-Dimensional Elasto-Plastic Analysis of Soils. Ph.D. Thesis, Technische Universiteit Delft, 1997.
[29] Herrera, I., Boundary Method: An Algebraic Theory. Pitman, Boston, 1984.
[30] Israil, A.S.M. & Banerjee, P.K., Advanced Time-Domain Formulation of BEM for
Two-Dimensional Transient Elastodynamics. International Journal for Numerical
Methods in Engineering. Vol.29, pp. 1421-1440, 1990.
[31] Israil, A.S.M. & Banerjee, P.K., Two-Dimensional Transient Wave-Propagation Problems by Time-Domain BEM. International Journal for Solids and Structures. Vol.
26, No. 8, pp. 851-864, 1990.
[32] Jaswon, M.A., Integral Equation Methods in Potential Theory: I. Proceedings of the
Royal Society of London. Vol. 275 A, pp. 23-32, 1963.
[33] Karabalis, D.L., Formulation of 3-D Dynamic SSI Analysis Involving Contact Nonlinearities by Time Domain BEM-FEM. Engineering Analysis with Boundary Elements.
Vol. 11, pp. 277-284, 1993.
[34] Karabalis, D.L. & Beskos, D.E., Dynamic Response of 3-D Flexible Foundations by
Time Domain BEM and FEM. Soil Dynamics and Earthquake Engineering. Vol. 4,
pp 91-101, 1985.
[35] Karabalis, D.L. & Gaitanaros, A.P., Contact Non-Linearities in Dynamic Soil Structure Interaction Problems by a Time Domain BEM-FEM Approach. In Boundary
Element Methods in Applied Mechanics. Ed. Tanaka, M. & Cruse, T.A., pp. 501-510,
Pergamon Press, Oxford, 1988.
[36] Kobayashi, S., Elastodynamics. Computational Methods in Mechanics, Volume III,
Boundary Element Methods in Mechanics. Ed. Beskos, D.E. Elsevier Science Publishers B.V. 1987.
[37] Kobayashi, S. & Kawakami, T., Application of BE-FE combined Method to Analysis
of Dynamic Interactions Between Structures and Viscoelastic Soil. Proceedings of the
International Conference on Boundary Elements VII. Como, Italy, pp. 6.3-6.12, 1985.
[38] Kobayashi, S. & Mori, K., Three-dimensional Dynamic Analysis of Soil-structure
Interaction by Boundary Integral Equation-nite Element Combined Method. Innovative Numerical Methods in Engineering. Ed. Shaw, R.P. et al. Springer, Berlin,
1996.
[39] Krenk, S, Kellezi, L., Nielsen, S.R.K & Kirkegaard, P.H., Finite Element and Transmitting Boundary Conditions for Moving Loads. Proceedings of the 4th European
Conference on Structural Dynamics, Eurodyn'99, Praha. pp. 447-452, 1999.
98
Bibliography
[40] Krenk & Kirkegaard, P.H., Radiation conditions for elastic waves. Technical University of Denmark, 1999. (to be published).
[41] Krishnasamy, G., Rizzo, F.J. & Liu, Y., Scattering of Acoustic and Elastic Waves by
Crack-Like Objects: the Role of Hypersingular Integrals. In Review of Progress in
Qantitative Nondescructive Evaluation, Vol. 11. Ed. Thompson, D.O. & Chimenti,
D.E., Plenum Press, New York, pp. 25-32, 1992.
[42] Krylov, V., Generation of Ground Vibrations by Superfast Trains. Applied Acoustics.
Vol. 44, pp. 149-164, 1995.
[43] Lachat, J.C. & Watson, J.O., Eective Numerical Treatment of Boundary Integral
Equations: A Formulation for Three-Dimensional Elastostatics. International Journal for Numerical Methods in Engineering. Vol. 10, pp. 991-1005, 1976.
[44] Lamb, H., On the Propagation of Tremors over the Surface of an Elastic Solid.
Philosocal Transactions Royal Society London Serie. A 203, pp. 1-42, 1904.
[45] Leung, K.L & Zavareh, P.B., 2-D Elastostatic Analysis by a Symmetric BEM/FEM
scheme. Engineering Analysis with Boundary Elements. Vol. 15, pp. 67-78, 1995.
[46] Madshus, C., Bessason, B. & Harvik, L., Prediction Model for Low Frequency Vibration from High Speed Railways on Soft Ground. Journal of Sound and Vibration.
Vol. 193, No. 1, pp. 195-203, 1996.
[47] Maier, G., Diligenti, M. & Carini, A., Variational Approach to Boundary Element
Elastodynamic Analysis and Extension to Multidomain Problems. Computational
Methods Applied to Mechanical Engineering. Vol. 92, pp. 193-213, 1991.
[48] Manolis, G.D., Ahmad, S. & Banerjee, P.K., Boundary Element Method Implementation for Three-Dimensional Transient Elastodynamics. In Developments in Boundary
Element Methods - 4. Ed. Banerjee, P.K. & Watson, J.O., Elsevier Applied Science
Publishers LTD, 1986.
[49] Manolis, G.D. & Beskos, D.E., Boundary Element Methods in Elastodynamics.
Unwin-Hyman, London, 1988.
[50] Mansur, W.J., A Time-Stepping Technique to Solve Wave Propagation Problems
Using the Boundary Element Method. Ph.D. Thesis, University of Southampton,
1983.
[51] Mita, A. & Luco, J.E., Dynamic Response of Embedded Foundation: a Hybrid Approach. Computational Methods Applied to Mechanical Engineering. Vol. 63, pp. 233259, 1987.
[52] Nardini, D. & Brebbia, C.A., A New Approach for Free Vibration Analysis using
Boundary Elements. Boundary Element Methods in Engineering. Ed. Brebbia C.A.,
Springer-Verlag, Berlin and New York, pp. 312-326, 1982.
Bibliography
99
[53] Nelson, J.T., Recent Developments in Ground-borne Noise and Vibration Control.
Journal of Sound and Vibration. Vol. 193, No. 1, pp. 367-376, 1996.
[54] Newmark, N.M., A Method of Computation for Structural Dynamics. Journal of the
Engineering Mechanics Division. ASCE 85, EM3, pp. 67-94, 1959.
[55] Nickel, R.E. & Sackman, J.J., Approximate Solutions in Linear Coupled Thermoelasticity. Journal of Applied Mechanics. Ser. E, Vol. 35, pp. 255-266, 1968.
[56] Oden, J.T., A General Theory of Finite Elements, part 2 - applications. International
Journal of Numerical Methods in Engineering. Vol. 1, pp. 247-259, 1969.
[57] Okumura, Y. & Kuno, K., Statistical Analysis of Field Data of Railway Noise and
Vibration Collected in an Urban Area. Applied Acoustics. Vol. 33, pp. 263-280, 1991.
[58] Ottosen, N. & Petersson H., Introduction to the Finite Element Method. Prentice
Hall International Ltd, 1992.
[59] Pao, Y.H., The Transition Matrix for the Scattering of Acoustic Waves and for
Elastic Waves. In Modern Problems in Elastic Wave Propagation. Ed. Miklowitz, J.,
& Achenbach, J.D. John Wiley & Sons, New York, pp. 123-144, 1978.
[60] Partridge, P.W., Brebbia, C.A. & Wrobel, L.C., The Dual Reciprocity Boundary
Element Method. Computational Mechanics Publications and Elsevier, Southampton
and Boston, 1992.
[61] Pavlatos, G.D. & Beskos, D.E., Dynamic Elastoplastic Analysis by BEM/FEM. Engineering Analysis with Boundary Elements. Vol. 14, pp. 51-63, 1994.
[62] Payton, R.G., An Application of the Dynamic Betti-Rayleigh Reciprocal Theorem
to Moving Point Loads in Elastic Media. Quarterly of Applied Mathematics. Vol. 21,
pp. 299-313, 1964.
[63] Pekeris, C.L., The seismic surface pulse. Proceedings of the National Academy of
Sciences of the United States of America. Vol. 41, pp. 469-480, 1955.
[64] Press, W.H., Flannery, B.P., Teukolsky, S.A. & Vetterling, W.T., Numerical Recipes
in Fortran. Cambridge U Press, Cambridge, 1994.
[65] Rebay, S., EÆcient Unstructured Mesh Generation by Means of Delaunay Triangulation and Bowyer-Watson Algorithm. Journal of Computational Physics. Vol. 105,
pp. 125-138, 1993.
[66] Richter, C., A Green's Function Time-Domain BEM of Elastodynamics. Computational Mechanics Publications. Southampton UK and Boston USA, 1997.
[67] Sanchez-Sesma, F.J., Diraction of Elastic Waves by Three-Dimensional Surface Irregularities. Bulletin of the Seismological Society of America. Vol. 73, pp. 1621-1636,
1983.
100
Bibliography
[68] Sanchez-Sesma, F.J. & Esquivel, J.A., Ground Motion on Alluvial Velleys under
Incident Plane SH Waves. Bulletin of the Seismological Society of America. Vol. 69,
pp. 1107-1120, 1979.
[69] Sanchez-Sesma, F.J. & Luzon, F., Seismic Response of Three-Dimensional Alluvial
Valleys for Incident P, S, and Rayleigh Waves. Bulletin of the Seismological Society
of America. Vol. 85, No. 1, pp. 269-284, February 1995.
[70] Sanchez-Sesma, F.J., Ramos-Matinez, J. & Campillo, M., An Indirect Boundary
Element Method Applied to Simulate the Seismic Response of Alluvial Valleys for
Incident P, S and Rayleigh Waves. Earthquake Engineering and Structural Dynamics.
Vol. 22, pp. 279-295, 1993.
[71] Sanchez-Sesma, F.J., Vai, R., Castillo, J.M. & Luzon, F., Seismic Response of Irregular Stratied Media Using the Indirect Boundary Element Method. In Soil Dynamics
and Earthquake Engineering VII. Ed. Cakmak, A.S. & Brebbia, C.A. pp. 225-233,
1995.
[72] Song, C. & Wolf, J.P., Unit-Impulse Response Matrix of Unbounded Medium by
Finite-Element Based Forecasting. International Journal for Numerical Methods in
Engineering. Vol. 38, pp. 1073-1086, 1995.
[73] Spyrakos, C.C. & Beskos, D.E., Dynamic Response of Flexible Strip-foundations by
Boundary and Finite Elements. Soil Dynamics and Earthquake Engineering. Vol. 5,
pp. 84-96, 1986.
[74] Spyrakos, C.C., Patel, P. & Beskos, D.E., Dynamic Analysis of Flexible Embedded Foundations: Plane Strain Case. In Computational Methods and Experimental
Measurements. Ed. Keramidas, G.A. & Brebbia, C.A. pp. 615-624, Springer-Verlag,
Berlin, 1986.
[75] Takemiya, H. & Kellezi, L., Paraseismic Behavior of Wave Impeding Block (WIB)
Measured for Ground Vibration Reduction. The 10th Japan Earthquake Engineering
Symposium. Yokohama, Japan, 1998
[76] Takemiya, H., Kellezi, L. & Imamura, T., Measurement of Vibration Propagation by
Pile Driving and Computer Simulation. Experimentation et Calcul en Genie Civil.
1997.
[77] Trochides, A., Ground-Borne Vibrations in Buildings Near Subways. Applied Acoustics. Vol. 32, pp. 289-296, 1991.
[78] Von Estor, O., Coupling of BEM and FEM in the Time Domain: Some Remarks
on its Applicability and EÆciency. Computers and Structures. Vol. 44, pp. 325-337,
1992.
[79] Von Estor, O., Dynamic Response of Elastic Blocks by Time Domain BEM and
FEM. Computers and Structures. Vol. 38, pp. 289-300, 1991.
Bibliography
101
[80] Von Estor, O. & Antes, H., Vibration Isolation Analysis in the case of Transient
Excitations by a Coupled BE/FE procedure. In Computational Mechanics '91. Ed.
Atluri, S.N., Beskos, D.E., Jones, R. & Yagawa, G. pp. 1148-1153, ICES Publications,
Atlanta, GA, 1991.
[81] Von Estor, O. & Kausel, E., Coupling of Boundary and Finite Elements for SoilStructure Interaction Problems. Earthquake Engineering Structural Dynamics. Vol.
18, pp. 1065-1075, 1989.
[82] Von Estor, O. & Prabucki, M.J., Dynamic Response in the Time Domain by Coupled Boundary and Finite Elements. Computational Mechanics. Vol. 6, pp. 35-46,
1990.
[83] Walther, R., Direct Integration using the Newmark Method, in ASKA Part II Linear Dynamic analysis. ASKA UM 228, Stuttgart, 1980.
[84] Wolf, J.P. & Chongmin, S., Finite-Element Modelling of Unbounded Media. John
Wiley & Sons, New York, 1996.
[85] Wong, H.L., Diraction of P, SV and Rayleigh Waves by Surface Topographies. USC,
Report No CE 79-05, Department og Civil Engineering, Los Angeles, USA, 1979.
[86] Zienkiewicz, O.C., The Finite Element Method third edition. McGraw-Hill Bokk Company Limited, 1977.
[87] Zienkiewicz, O.C. & Taylor, R.L., The Finite Element Method, Vol. 2. McGraw-Hill,
London, 1991.
102
Bibliography
Appendix A
Integral of Normal Vectors over
Semicircular Surface
The product of two normal vectors must be integrated over a hemisphere S .
Figure A.1: Hemisphere with coordinate system and normal vector.
The integral is
Iij =
Z
S
ni (x)nj (x)dS
(A.1)
Using a spherical coordinate system leads to the following expressions for the normal
vectors
n1
n2
n3
= sin cos = sin sin = cos (A.2)
The dierential areal dS is in spherical coordinates given as
dS = 2 sin dd
(A.3)
103
104
A. Integral of Normal Vectors over Semicircular Surface
where is radius in the semicircle. (A.2) and (A.3) is inserted into (A.1)
Z
2 Z =2
ni nj 2 sin dd
0 0
(A.4)
Each of the 9 combinations of the directions in (A.4) is now considered
3 d = 2 2 = I22
sin
3
0
0
2
=2
2 sin d = 2 2 1
d
cos
3
0
0
2
=2
cos sin d 0 sin3 d = 0
0
I11
=
2
Z
I33
=
2
Z
I12
=
2
Z
I13
=
2
Z
0
2
cos2 d
Z
=2
Z
Z
2
cos d
Z
0
=2
sin2 cos d = 0 = I23
(A.5)
From (A.5) it is seen that
2
Iij = 2 Æij
(A.6)
3
In another coordinate system (X1 ; X2; X3), the coordinates Ni to the normal vector are
given by the original coordinates ni as
Ni = Tij nj ; ni = Tij 1 Nj = Tji Nj
(A.7)
where Tij is the transformation tensor. Insertion of (A.7) into (A.1) gives the integral in
the (X1; X2; X3 ) coordinate system
Iij =
Z
S
Ni (x)Nj (x)dS = Tik
Z
S
nk (x)nl (x)dS Tjl =
2 2Æ T T = 2 2 T T = 2 2Æ
(A.8)
3 kl ik jl 3 il jl 3 ij
From (A.8) it is seen that (A.6) contains an arbitrary Cartesian coordinate system.
Appendix B
Element Types and Gauss
Quadrature
B.1
Element Types in Rectangular Coordinates
2-node line, 3-node triangular, 8-node quadrilateral and 6-node triangular elements are
used.
B.1.1
2-Node Line Element
Figure B.1: 2-node line element.
The shape functions for the 2-node line element shown in gure B.1 are
N (1) ( ) = 12 1 + 1(1)
N (2) ( ) = 12 1 + 1(2)
9
=
;
(B.1)
where i(n) are the parent domain coordinates for node number n.
B.1.2
3-Node Triangular Element
The shape functions for the 3-node triangular element shown in gure B.2 are in the
standard form, Zienkiewicz [86]
= Ni Nj Nm
(B.2)
where are a 2 2 identity matrix, and
h
N
i
I
I
I
I
105
106
B. Element Types and Gauss Quadrature
Figure B.2: 3-node triangular element.
where
1 (ai + bi x1 + ci x2 )
Ni = 2
1 (aj + bj x1 + cj x2 )
Nj = 2
1 (am + bm x1 + cm x2 )
Nm = 2
(B.3)
ai = xj1 xm2 xm1 xj2
bi = xj2 xm2
ci = xm1 xj1
(B.4)
with the other coeÆcients obtained by a cyclic permutation of subscripts in the order i,
j , m, and
1 xi1j xi2j
2 = det 1 x1 x2 = 2 (area of triangle ijm)
(B.5)
1 xm1 xm2
First order dierentiation of (B.3) yields
@
1 bi b j bm
=
(B.6)
@x
2 ci cj cm
"
N
#
B.1. Element Types in Rectangular Coordinates
B.1.3
107
8-Node Quadrilateral Element
Figure B.3: 8-node quadrilateral element.
Introduce i(0) = i i(n), i = 1; 2, n = 1; :::; 8. Then the shape functions and their gradients
for the corner nodes shown in gure B.3 become, Ottosen & Petersson [58]
N (n) ()
= 41 (1 + 1(0) )(1 + 2(0) )(1(0) + 2(0) 1)
(0)
(0)
(0) (n)
@
1
(n)
; n = 1; 2; 3; 4:
(B.7)
@ N ( ) = 4 (1 + 2 )(21 + 2 )1
(0)
(0)
(0) (n)
@
1
(n)
@ N ( ) = 4 (1 + 1 )(22 + 1 )2
The shape functions for the midside nodes are
9
>
>
>
>
=
1
>
>
>
>
;
2
N (n) ()
@
(n)
@1 N ( )
@
(n)
@2 N ( )
N (n) ()
@ N (n) ( )
=
=
=
=
=
@
@
(n)
@ N ( ) =
1
2
(0)
1
2
2 (1 (1 ) )(1 + 2 )
1 (1 + 2(0) )
1 (1 (1 )2 ) (n)
2
2
(0)
1
2
2 (1 + 1 )(1 (2 ) )
1 (1 (2 )2 ) (n)
1
2
2 (1 + 1(0) )
9
>
>
>
>
=
>
>
>
>
;
9
>
>
>
>
=
>
>
>
>
;
; n = 5; 7
(B.8)
; n = 6; 8
(B.9)
108
B.1.4
B. Element Types and Gauss Quadrature
6-Node Triangular Element
Figure B.4: 6-node triangular element.
The shape functions with node denition as shown in gure B.4 are given by, Banerjee et
al. [10]
N (1) ( ) = (0) (2 (0) 1)
N (2) ( ) = 1 (21 1)
N (3) ( ) = 2 (22 1)
(B.10)
N (4) ( ) = 41 (0)
N (5) ( ) = 41 2
N (6) ( ) = 42 (0)
9
>
>
>
>
>
>
>
>
>
>
>
>
>
=
>
>
>
>
>
>
>
>
>
>
>
>
>
;
1 4(0) 41 1 0 4((0) 1) 42
42
(0)
(0)
1 4
0 42 1
41 41 4( 1)
where (0) = 1 1 2.
@ (n)
N ( ) =
@
"
#
(B.11)
B.2. Element Types in Cylindrical Coordinates
B.2
109
Element Types in Cylindrical Coordinates
8-node quadrilateral and 6-node triangular elements are used in cylindrical coordinates.
B.2.1
8-Node Quadrilateral Element
Figure B.5: 8-node quadrilateral element in cylindrical coordinates.
Assume that the following conditions are fullled (B.12):
ri > 0, i 2 [0; 2 ], r1 = r4 = r8 , r2 = r3 = r6 , r5 = r7 , 1 = 2 = 5 , 3 = 4 = 7 and
6 = 8 . The shape functions are then given by the following expressions
N (1) (r; )
= (r
1
N (2) (r; ) =
(r
2
N (3) (r; )
N (4) (r; )
N (5) (r; )
N (6) (r; )
N (7) (r; )
N (8) (r; )
= (r
2
= (r
1
= (r(r
5
= (r(r
2
= ((rr
5
= (r(r
1
(r
r2 )( 3 )
(r
r2 )(1 3 )(6 1 )
"
(r r1 )( 3 )
(r
r1 )(1 3 )(6 1 )
"
(r r1 )( 1 )
(r
r1 )(3 1 )(6 3 )
"
(r r2 )( 1 )
(r
r2 )(3 1 )(6 3 )
r1 )(r r2 )( 3 )
r1 )(r5 r2 )(1 3 )
r1 )( 1 )( 3 )
r1 )(6 1 )(6 3 )
r1 )(r r2 )( 1 )
r1 )(r5 2 )(3 1 )
r2 )( 3 )( 1 )
r2 )(6 3 )(6 1 )
"
r5 ) 6
r1
r5 ) 6
r2
r2 ) 3
r5
r1 ) 3
r5
1
+ (1
r5
1
+ (1
r5
6
+ (6
r2
6
+ (6
r1
where ri; i, i = 1; :::; 8 are the polar coordinates for node number i.
#
)
#
)
#
)
#
)
(B.12)
110
B. Element Types and Gauss Quadrature
Figure B.6: 6-node triangular element in cylindrical coordinates.
B.2.2
6-Node Triangular Element
Assume that the following conditions are fullled (B.13):
ri > 0 for i 6= 1, i 2 [0; 2 ], r1 = 0, r2 = r3 = r5 , r4 = r6 , 1 = 0, 2 = 4 and 3 = 6 .
The shape functions is then given by the following expressions
N (1) (r; )
= (r
r2 )(r r4 )
r2 r4
r(r r4 )( 3 )(
r2 (r2 r4 )(2 3 )(2
r(r r4 )( 5 )(
r3 (r3 r4 )(3 5 )(3
r(r r2 )( 3 )
r4 (r4 r2 )(4 3 )
r(r r4 )( 2 )(
r5 (r5 r4 )(5 2 )(5
r(r r2 )( 2 )
r6 (r6 r2 )(6 2 )
N (2) (r; )
=
N (3) (r; )
=
N (4) (r; )
=
N (5) (r; )
=
N (6) (r; )
=
5 )
5 )
2 )
2 )
3 )
3 )
(B.13)
B.3. Gauss Quadrature
B.3
111
Gauss Quadrature
Gauss quadrature is used in the evaluation of the integrals in (2.26), (3.4) and (4.33). For
all elements, 3rd order Gauss integration is used.
B.3.1
3rd Order Gauss Quadrature
Figure B.7: Gauss points for 3rd order Gauss quadrature on quadrilateral element.
The Gauss points in a certain coordinate direction are given by, Ottosen & Petersson [58]
2
6
4
3
r
0
r
7
5
; r=
s
3
5
The corresponding Gauss weights become, Ottosen & Petersson [58]
5 1
w= 8
5 9
The quadrature then becomes
2
3
6
4
7
5
(B.14)
(B.15)
3
X
I1D = m f ( )d ' w(i)f (ri)
S
i=1
(B.16)
Z
3 X
3
X
I2D = m f (1 ; 2 )d1 d2 '
w(i)w(j )f (ri; rj )
S
i=1 j =1
(B.17)
Z
112
B. Element Types and Gauss Quadrature
B.3.2
3rd Order Gauss Quadrature on Triangular Element
Figure B.8: Gauss points for triangular element.
The Gauss points and weights are shown in table B.1, Bathe [11].
Points
Weights
L1 L2 L3
a 1/3 1/3 1/3
0.2250
1 =0.0597158717898
b 1 1 1
c 1 1 1 0.1323941527885 1 =0.4701420641051
d 1 1 1
e 2 2 2
2 =0.7974269853531
f 2 2 2 0.1259391805448
g 2 2 2
2 =0.1012865073235
Table B.1: Gauss points for 3rd order Gauss quadrature on triangular element.
The quadrature is evaluated as
1 7 w(i)f ( ; )
I=
(B.18)
1;i 2;i
2 i=1
where 1;i and 2;i are the coordinates of the Gauss points. These are related to the area
coordinates L1 , L2 and L3 as follows
X
1
2
=
=
L1 1;1 + L2 1;2 + L3 1;3
L1 2;1 + L2 2;2 + L3 2;3
(B.19)
Appendix C
Analytical Solutions
In what follows classic solutions to elastodynamic problems in 2D and 3D are presented.
Lambs problems all consist of concentrated unit loads applied at the surface of an isotropic,
homogeneous, linear elastic halfspace. In 2D the load is applied as a concentrated unit
line impulse, whereas the 3D solutions apply to a concentrated unit Heaviside type load.
Numerical integration is used on the solutions to obtain results for loads applied over
nite areas.
C.1
Lamb's Problem in 2D
A downward directed vertical concentrated line impulse of magnitude 1 is acting at the
origin on the boundary of an isotropic, homogeneous, linear elastic halfspace at the time
t = 0. A complete derivation for an arbitrary angle of the impulse can be found in Eringen
& Suhubi [17].
With the impulse placed at the origin the Green functions, u1(x1; t) and u2(x1 ; t) at surface
point (x1 ; x2) = (x1 ; 0) become
u1 (x1 ; t) =
u2 (x1 ; t) =
2
U ( )
c22 jx1 j 1
1
U ( )
c22 jx1 j 2
where = xt and
p
1
8
<
U1 ( ) = :
8
>
>
>
>
<
U2 ( ) = >
>
>
>
:
0
(C.1)
p
(2 2 c2 2 ) 2 c1 2 c2 2 2
(c2 2 2 2 )4 +16 4 ( 2 c1 2 )(c2 2 2 )
p
(c2 2
(c2 2 2 2 )2 2 c1 2
2 2 )4 +16
p4 (2 c1 2)(c2 2
(c2 2 2 )2
2
p
p
2)
2 c1 2
4 2 2 c1 2 2 c2 2
0
u1 and u2 are plotted in gures C.1 and C.2.
113
; c1 1 < j j < c2 1
;
otherwise
(C.2)
; c1 1 < j j < c2 1
;
j j > c2 1
;
otherwise
(C.3)
114
C. Analytical Solutions
1.4
1.2
1
πµx1u1/c2
0.8
0.6
0.4
0.2
0
−0.2
−0.4
0.8
1
1.2
1.4
1.6
1.8
2
2.2
Figure C.1: u1 for Lamb's problem in 2D, impulse.
c1t/x1
10
πµx1u2/c2
5
0
−5
1
1.5
2
2.5
3
Figure C.2: u2 for Lamb's problem in 2D, impulse.
Using numerical integration (C.1) is integrated from 0 to t to obtain the solution for a
Heaviside load, see gures C.3 and C.4.
c1t/x1
−7
12
x 10
10
πµx1u1/c2
8
6
4
2
0
0
0.05
0.1
0.15
0.2
0.25
Figure C.3: u1 for Lamb's problem in 2D, Heaviside load.
c1t/x1
C.1. Lamb's Problem in 2D
115
−6
x 10
12
10
πµx1u2/c2
8
6
4
2
0
−2
−4
0.1
0.15
0.2
0.25
0.3
0.35
c1t/x1
Figure C.4: u2 for Lamb's problem in 2D, Heaviside load.
u1 is only given until the arrival of the S-wave, see gure C.1, which means that the
Rayleigh
wave component is missing. u2 in gure C.4 becomes constant for large values
c
t
of x .
1
1
C.1.1
Solution for Distributed Loads
Both the analytical solution for the impulse and the numerical result for the Heaviside load
are integrated numerically over a nite area in space to obtain results comparable with
the BEM and FEM results. In gures C.5-C.8 examples of the solutions for distributed
loads are shown.
0.7
0.6
0.5
πµx1u1/c2
0.4
0.3
0.2
0.1
0
−0.1
−0.2
−0.3
0
0.5
1
1.5
c1t/x1
2
2.5
Figure C.5: u1 for Lamb's problem in 2D for a distributed load.
116
C. Analytical Solutions
2.5
2
1.5
1
πµx1u2/c2
0.5
0
−0.5
−1
−1.5
−2
−2.5
−3
0.5
1
1.5
2
2.5
3
Figure C.6: u2 for Lamb's problem in 2D for a distributed load.
c1t/x1
−3
16
x 10
14
12
πµx1u1/c2
10
8
6
4
2
0
−2
0
0.5
1
1.5
2
2.5
Figure C.7: u1 for Lamb's problem in 2D with a distributed Heaviside load.
c1t/x1
0.05
0
πµx1u2/c2
−0.05
−0.1
−0.15
−0.2
0.5
1
1.5
2
2.5
Figure C.8: u2 for Lamb's problem in 2D with a distributed Heaviside load.
c t/x
1
1
C.2. Lamb's Problem in 3D
C.2
117
Lamb's Problem in 3D
An isotropic, homogeneous, linear elastic halfspace is subjected to a unit concentrated
downward directed force H (t) normal to the surface z = 0 (cylindrical coordinates (r; ; z))
at the origin of the coordinate system at t = 0. The z-axis coincide with the x2 -axis. This
problem was rst introduced by Lamb [44] but the solution of the displacement eld was
obtained for Poisson materials by Pekeris [63]. A complete deriviation of the general
solution can be found in Eringen & Suhubi [17].
C.2.1
Vertical Displacement
The force is placed at the point (1,0,3). The vertical displacement at the surface point
(x1 ,0,x3 ) at the distance R = (x1 1)2 + (x3 3)2 from the force is then given by
q
uz = 0 ; < c22 =c1
2 Ki (1 2v 2 )2 [(c2 =c2 ) v 2 ]1=2
i
2 1
i
+
+
2R 4(1 c22 =c21) i=1
( 2 vi2)1=2
K3 (1 2v32 )2 (v32 c22 =c21 )1=2
; c2 =c1 < < 1
(v32 2 )1=2
1
K3 (1 2v32 )2 (v32 c22 =c21 )1=2
1
+
H ( )
uz =
R 4(1 c2 =c2 )
(v2 2 )1=2
uz =
1
4
1
X
#
"
(C.4)
#
2 1
3
; >1
where = c2t=R and Ki are given by
1
(C.5)
Ki =
16 [1 (c22=c21)] ( i2 j2)( i2 k2 ) ; i 6= j 6= k
i are the Rayleigh poles of which one is the true Rayleigh pole = c2 =cR and is in the
interval 1 < 3 < 1. 1 < 2 < 3 are found as solutions to the following equation
D( ) = 16 6 1 (c22 =c21 ) + 8 4 3 2(c22 =c21 ) 8 2 + 1
(C.6)
v1 , v2 are found as the real solutions to the complex equation
(1 2v2)2 + 4v2(v2 1)1=2(v2 c22 =c21)1=2 = 0
(C.7)
and v3 = > 1 is the true Rayleigh root found as the only solution in the interval [0; 1]
to the Rayleigh equation
2 8 = 0
3 8 2 + 8
(C.8)
1 1 For Poisson materials ( = ) one has
p
p
3
3
c22 1 2 1 2 3
3
+
2
=
; v1 = ; v2 =
; v3 =
(C.9)
2
c1 3
4
4
4
h
i
h
i
118
C. Analytical Solutions
6
4
15πµruz
2
0
−2
−4
−6
0.5
1
1.5
c2t/r
Figure C.9: Vertical displacement uz at the surface of a Poisson material from a downward
directed unit concentrated force H (t).
C.2.2
Horizontal Displacement
ur = 0 ; < c2 =c1
"
2
1 K (k ) +
ur = 2
2
2
r(1 c2 =c1 )1=2 8(1 c22 =c21 )
#
3
X
2
2
2
Ki (1 2vi )(1 vi )(ni k ; k) ; c2 =c1 < < 1
i=1
"
2
1 K (1=k) +
ur = 2
2
2
1
=
2
rk(1 c2 =c1 )
8(1 c#22=c21)
3
X
AF Ki (1 2vi2 )(1 vi2 )(ni ; 1=k) +
2r( 2 2 )1=2 H ( ) ; > 1
i=1
where K (k) and (n; k) are complete elliptic integrals of the rst and third kind
Z =2
d
K (k) =
2
0 (1 k sin2 )1=2
Z =2
(n; k) = 0 (1 + n sin2 )(1d k2 sin2 )1=2
k and n are given as
2 c22 =c21
1 c22=c21
k2 ( ) =
;
n
=
i
1 c22=c21
c22 =c21 vi2
and A is given by
i 1
h
1
A = (2 2 1)3 8(1 c22 =c21 ) 6 4 2 + 1
2
for Poisson materials (C.14) reduces to A = 41
(C.10)
(C.11)
(C.12)
(C.13)
(C.14)
The elliptic integrals (C.11) and (C.12) can be solved by numerical integration or they can
be expressed by incomplete elliptic integrals that have been tabulated, see e.g. Abramowitz
& Stegun [2].
C.2. Lamb's Problem in 3D
4
2
0
−2
−4
−6
−8
0.5
1
r
u
1.5
c2t/r
2
2.5
119
at the surface from a downward directed unit
Solution for Distributed Load
5
0
−5
−15
0
1
2
c t/r
1.5
2
2.5
Figure C.11: Vertical displacement uz at the surface from a downward directed unit force
H (t) applied over a nite area.
0.5
The exact solutions (C.4) and (C.10) for concentrated forces are integrated over a nite
area in order to nd a result comparable with the numerical results. Special attention is
given to the singular behaviour of the solutions. In gures C.11 and C.12 examples of the
results are shown.
C.2.3
Figure C.10: Horizontal displacement
concentrated force H (t).
−12
−10
15πµrur
−10
15πµruz
120
5
0
−5
0.5
Figure C.12: Horizontal displacement
force H (t) applied over a nite area.
−15
0
−10
15πµrur
r
u
1
c2t/r
1.5
2.5
C. Analytical Solutions
2
at the surface from a downward directed unit
Appendix D
Newmark Integration Scheme
As mentioned in Section 4.2 the dynamic system of equations must be solved with respect
to time by numerical integration. The procedure that has been chosen for the numerical
integration is the implicit integration scheme known as the Newmark scheme, Newmark
[54], Walther [83], Argyris & Mlejrek [8].
A non-dimensional time is introduced
= t ttL
(D.1)
The displacement and velocity and the beginning of a time step are assumed known and
it is assumed that the acceleration  ( ) can be interpolated between the start and end
values  L and  L . From this the displacement and velocity ( ), _ ( ) can be calculated
within a time step
u
u
u +1
u
_ ( ) = _ L + t
u
u
Z
u
( ) =
uL
0
 (r)dr ; 0 1
(D.2)
u
+ t _ L + (t)
2
u
ZZr
0
0
u
 (s)dsdr ; 0 1
u
(D.3)
Applying the reverse order of integration in (D.3) and using = 1 give the following
expressions for the increments of displacement and velocity
=
u
uL+1
_ = _L
u
u +1
uL
= t _ L + (t)
_ = t
u
uL
Z
0
1
2
Z
1
0
(1
)u( )d
 ( )d
(D.4)
(D.5)
u
Within each time step a variation of the acceleration  ( ) is now assumed analogous to
the one used in space
u
 ( ) = N ( ) L + (1
u
u

N ( ))uL+1
(D.6)
121
122
D. Newmark Integration Scheme
where N ( ) is the interpolating function which determines the variation over the time
step. N ( ) must satisfy the conditions N (0) = 1 and N (1) = 0. Using (D.6) in (D.4) and
(D.5) gives
"
# h
Z
i
N
(
)


_
= L
=
t
+
(
t
)
(1
)
d
(D.7)
L
L
L
L
(1 N ( ))
u
u +1
u
1
2
u
Z 1"
u
0
N ( )
(1 N ( ))
#
h
i
_ = _ L _ L = t
d  L  L
The following denitions are introduced, Newmark [54]
Z
(1 )N ( )d = 21 u
u +1
u
0
u
u +1
(D.8)
u +1
(D.9)
1
Z
0
=1 (D.10)
Now from (D.7), (D.8), (D.9) and (D.10)
= t _ L + 21 (t)  L + (t) 
(D.11)
_ = t  L + t 
(D.12)
In (D.11) and (D.12) the unknowns are , _ and  so an extra equation is required.
This is found from (4.6), which is solved for 
 =
[ (tL ) _ ]
(D.13)
Now (D.13) is substituted into (D.11) and (D.12) and a linear system of equations which
gives the sought solution is established
"
#"
# "
#
+ t t = t _ L + t  L + t (tL ) (D.14)
t  L + t (tL )
t
+ t
_
Alternatively to (D.13) and (D.14), (D.11) and (D.12) can be substituted into (D.13) and
rearranged
1
0
N ( )d
u
u
u
u
2
2
u
u
u
u
u
u
u
u
1
M
2
M
f
2
K
K

C
u
K
u
C
M
u
C
h
Mu
2
2
Mu
Mu
u
i
1
2
f
f
= + t + t ( (tL ) [ _ L + (1 )t  L ]
h
i
(D.15)
L + t _ L + t (1=2 ) L
Rearranging (D.11) and (D.12) gives the displacement and velocity to use with (D.15)
(D.16)
L = L + t _ L + t (1=2 ) L + t  L
_ L = _ L + (1 )t  L + t  L
(D.17)
The main advantage of the Newmark integration scheme is that with = 0:25 and = 0:5
it can be shown to be unconditionally stable.
In the implementation of the time domain FEM a sparse formulation of the system matrices has been used. Solving of the sparse system was done by an iterative biconjugate
gradient approach, Press et al. [64]. Description of the sparse formulation, the solver and
the preconditioner can be seen in Appendix E.
uL+1
M
K u
2
u
u +1
u
u +1
u
2
C
u
K
1
f
C u
+1
u
2
u
u
u +1
2
u +1
u
Appendix E
Sparse Formulation
The sparse formulation used closely follows the one described by Press et al. [64].
The purpose of the sparse formulation is to reduce the memory required for the storage
of a non-symmetric matrix of the dimension N N . The sparse storage consists of two
vectors D and I
aDi ; i = 1; M containing double precision values
aIi ; i = 1; M containing integer values
where M is the maximum size allowed for the sparse arrays.
The sparse formulation follows a set of 7 rules:
A
a
a
In the rst N -locations of aD all the values from the main diagonal of A are
stored. This includes zero values.
In the rst N -locations of I the indices from D containing the rst non-diagonal
terms of the corresponding row in are stored. If a row do not have any nondiagonal terms, one plus the index of the last row with non-diagonal terms is
stored.
Location 1 of I is always N + 2.
Location N + 1 of I is one plus the last index in D for the last non-diagonal
term for the last row. I.e. I (N + 1) gives one plus the number of locations used
in D and I .
Locations in D N + 2 are used to store the non-diagonal terms row by row,
and within each row the non-diagonal terms are stored in column order.
Locations in I N + 2 are used to store the column indices corresponding to
the non-diagonal terms in D .
Location N + 1 of D is always 0.
a
a
A
a
a
a
a
a
a
a
a
a
a
123
124
E. Sparse Formulation
Example E.1
33
Matrix
A small example to illustrate the sparse formulation
2
3
1:00 0 3:17
A = 4 1:51 2:40 8:20 5
0 7:33 0:35
(E.1)
This gives the following sparse formulation
i
Di
I
ai
a
1
1.00
5
2
2.40
6
3
0.35
8
4
0
9
5
3.17
3
6
1.51
1
7
8.20
3
8
7.33
2
The gain in this example is negative (meaning that the sparse formulation uses more memory than the
full formulation), but for larger sparse matrices the gain is signicant (see Example E.2).
Example E.2 FEM Stiness Matrix
Table E.1 shows the memory used for the storage of stiness matrices generated by the FEM program
used in the FEM/BEM coupling, see Chapter 4.
Dof
48
478
726
1150
2506
Elements
31
420
654
1060
2373
Full matrix [Mb]
0.0180
1.79
4.12
10.3
49.06
Sparse matrix [Mb]
0.00518
0.0608
0.0936
0.150
0.332
Table E.1: Memory use for sparse stiness matrix from FEM program.
In table E.1 it can be seen that the memory used for the sparse formulation is reduced to 29 %, 3.4 %,
2.3 %, 1.5 % and 0.68 % of the memory used for the full matrix storage.
E.1. Sparse Solver
E.1
125
Sparse Solver
(E.2)
The N N linear system in (E.2), where is sparse, is solved by an iterative biconjugate
gradient approach, Press et al. [64].
The idea behind the conjugate gradient method is to minimize the function
1
f( ) = T
(E.3)
2
And this function is minimized when its gradient rf is zero. The gradient of (E.3) is
rf =
(E.4)
which is equivalent to the problem of solving (E.2).
The biconjugate method, on the other hand, is not simply connected with function minimization. For further reading, see e.g. Press et al. [64], Fletcher [19] or Golub & Van
Loan [26]. The preconditioner used with the solver is simply the diagonal part of .
Ax
=
b
A
x
x Ax
Ax
bx
b
A
126
E. Sparse Formulation
Appendix F
Integration of Functions for Moving
Force
The assumptions for the integration results can be seen in Section 5.3.
2q
4
1 (c v ) + 2~r v r~
r~ v arctan(
s)
q
g~ (; r~(t); c)(1 l )d =
t
c v
a
(c v ) c + v
q
c (c v ) + 2~r v r~
c r~ arctan(
s)
r~ cqarctan(s )
q
+
2
2
4
+
v (c v )
v (c v ) c + v
v
c +v
q
t
(c v )s + 2~r c vs + s
l
c tq
tq
l arctan(s )
l arctan(s )
2
2
vs
v
c +v
( c +v )
q
3
r~ (c v )s + 2~r c vs + s
r~ c tl arctanh(s )
r~ c arctanh(s ) 5
2 v ps
+2
+ 4 v ps
(F.1)
vs
Zb
(1)
2
2
1
2
2
2
1
2
1
2
1
where
2
1
3
3
1
2 1
2
2
1
2
2
1
2
1 1
2
1
2
1
3
1 1
2
2
1
2
1
3
2
1
2
2
1
2
2
1
2
4
3
1
1
2 1
1
2
2 2
1
4
4
1
2
r~1
v1
s1 = s2 =
2
1
2
2
1
2
1
1 1
2
1
3
2
3
2
1
1 1
1
2
1
3
2
1
2
2
1
2
2
2
2
1
1 1
2
1
2
r~12 c2
v12q
s3 = q
c2 + v12 + cv21 r~v112
r~12 + 2~r1 v1 v12 2
2s2 + 2~r1c2 s1
s4 = p q 2 2 2 v1 2 s1
2 s2 (c v1 )s1 + 2~r1 c v1 + s2
c2 2
127
(F.2)
2
1
128
F. Integration of Functions for Moving Force
2
6
6
t4
2s
p
tl r~ c arctanh p p
tl c s 1
g~ (; r~(t); c)(1 l )d =
+ r~ (c v )
p
v s
a
p
p
p
p
v tl s v tl s v tl s
v tl s
tl c arctan(
s
)
q
+
2 r~ (c v ) + r~ (c v )
+
(c v ) c + v r~ c (c v ) r~ c (c v )
s rc 2~r c arctanh ps ps
v tl arctan(
s)
q
+ v r~s c s 2~r cqarctan(s ) +
p
v s
(c v ) c + v
v
c +v
p
ps p
2
v s
v
c s
r~ c arctan(
s
)
q
+
+
v r~ (c v ) r~ (c v ) v (c v ) c + v r~ c (c v )
p
t
v r~ arctan(
s
)
t
t
l s
l c arctan(s )
ls
q
+ q
+
+
(c v ) c + v r~ c s # v c + v v r~
p
v s
2ps ps
(F.3)
c (c v )
v
c v
Zb
(2)
2
1
2
1
1
2
2
1
2
1
2
1
2
1 1
2
1
2
1
2
1
2
where
2 2
1
1
2
1
2
2
s2 = (c2
2
1
2
2
1
3
1
2
1
2
1
2
v12 )s21 + 2~r1 c2
3
2
2
1
s1 2 c2
+ r~
v1 1 v12
(3)
g~
(; r~(t); c)(1 l )d
ps p s
2
2
2
1
r~1
r~1 c2
+
v1 v1 (c2 v12 )
q
s
c2 + v12 p 3
s4 =
s2
2
c
s5 = r~12 2
v1
a
2
5
1 1
s3 = Zb
2
1
2
4
2
1
2
2
s
2 5 +2~1 2 v1
1
2
5
2
4
1
4
2
2
2
1
r~1
v1
s1 = 2
1
2 2
1
4
2
1
3
1
2
1
2
2
2
2
2
1
2
2
2
2
2
1
2
1
3
2
4
2
1
2
1
2 2
1
2
1
2
3
1
2
2
1
2
1 1
1
2
2
1
1
2
2
1
5
4
1
2
2
2
1
4
2
1
2
3
2
2
2
1 1
2
1
4
2
1
2
3
1
2
2
2
1
r1 c2 vs11
5 +2~
2 s5 s2
2
r~1 v1 arctan
p
s2 s2
p
(F.4)
= 1t
pss21s4
2
tl
4
3
5
p
arctan
p
s
2
pss21s4 s2
s2
(F.5)
where
s1 = c2 2 r~12 + 2~r1 v1 v12 2
s2 = c2 v12
s3 = s2 2 + 2~r1 v1 r~12
v
s4 = + r~1 1
s2
(F.6)
129
Zb
(1)
t~
(; r~(t); c)(1 l )d
= 0 ; (; ) = [(1; 2); (1; 1); (2; 2)]
Zb
1 c + r~ 2~r v + v t~ (; r~(t); c)(1 l )d =
t
a
!
!
c + v v r~
s s t c + r~ s n s t v 2v s r~ + v r~ arctan
a
(1)
2
2
1
l
1
1
1
2
1
2
2
1
1 1 1
( r~ + v ) r~ (
; (; ) = (2; 1)
2
6
2
1
1
2
1
2
1
2
1 1
2
2
1 l
s1 s2
! 3
c2 2 + r~12 2~r1 v1 + v12 2 2
( r~1 + v1 )2
1 1
c2 + v12 )
1 1
3
2
2
1 1
(F.7)
q
s1 =
c2 + v12
q
s2 = c2 2 r~12 + 2~r1 v1 Zb
(F.8)
v12 2
(2)
t~
(; r~(t); c)(1 l )d
= 0 ; (; ) = [(1; 1); (2; 2)]
Zb
1
t~ (; r~(t); c)(1 l )d = a
! t
!
"
c v + r~ v
r~ c v r~ s c +
2n s 2s arctanh v s s r~ c 2 arctan
ss
!
c v + r~ v
2s v s r~ s + arctan
v r~ s 3s r~ v s s +
a
(2)
2 2
1
3
2
2
3 1 2 1 1
2
1 1 2
2
!
2
4 2
1
2
1
s3 s2
2
1
1 1
3 2
1 1
3 3
1 1 1
!
3
2
1 1 1
3
3 1 1 1 2
2
2
4s3arctanh vr~1sc s r~13c2 v1 2s3 arctanh vr~1sc s r~12c2v12 2 + s3v14t1 s1s2 +
1 1 2
1 1 2
!
!
2
2
c2 v12 + r~1 v1 3
c v1 + r~1 v1 2 2 2
4 arctan
v1 r~1 s1 c 2 arctan
v1 r~1 s1 2 c2
s3 s2
s3 s2
!
!
#
c2 v12 + r~1 v1 4 2
c2 v12 + r~1 v1 5
2
2 arctan
v1 r~1 s1 + arctan
v1 r~1 s1 s3 s2
s3 s2
! 12
2 2
2
2 2
6
c
+
r
~
2~
r
1 v1 + v1 1
4
1
(~r1 v1 ) 2 v1 r~1
(s3 s1) 1
(~r1 v1 )2
; (; ) = [(1; 2); (2; 1)]
(F.9)
where
s
c2
v2
q 1
s2 = c2 2 r~12 + 2~r1 v1 q
s3 =
c2 + v12
s1 = r~12
v12 2
(F.10)
130
F. Integration of Functions for Moving Force
Zb
(3)
t~
(; r~(t); c)(1 l )d
= 0 ; (; ) = [(2; 1); (1; 1); (2; 2)]
Zb
1 t~ (; r~(t); c)(1 l )d = n c + r~ 2~r v + v t
a
!
!
c + r~ v v s t v + v r~ arctan
s + 2v s r~ r~ s + s t c a
(3)
2
1 1
(~r
v1 ) r~1
;
s1 =
( c +v )
(; ) = (1; 2)
q
q
1 1
1 1
6
2
1
where
2
2
1
1
2
2
1
3
2
2
2
2
1
s1 s2
c2 2 + r~12
(~r1
2
1
2
2~r v + v v )
1 1
1
2
r~12 + 2~r1 v1 v12 2
The results for the functions multiplied by l are similar.
2
2
1 1
1 1 1
c2 + v12
s2 = c2 2
2
1
1 1
2
1
2
!
1
l
2
3
2
(F.11)
(F.12)
Fly UP