- Email: [email protected]
The time domain boundary element method for elastodynamic problems
‘Math1 Comput.
hfodelling
Vol. 15, No. 3-5, pp. 119-129,
1991
Printed in Great Britain. All rights reserved
THE
TIME DOMAIN BOUNDARY FOR ELASTODYNAMIC
0895-7177/91 $3.90 + 0.00 Copyright@ 1991 Pergamon Press plc
ELEMENT PROBLEMS
METHOD
Jo& DOMINGUEZ AND RAFAEL GALLEGO Escuela Tkcnica Superior de Ingenieros Industriales Universidad de Sevilla, Sevilla, Spain
INTRODUCTION The number of papers published on the time domain formulation and implementation of the Boundary Element Method (BEM) is rather small. Cole e2 al. [l] were the first to present the formulation for two-dimensional antiplane problems. Plane strain problems were studied by Niwa et al. [2] and Manolis [3] using a three-dimensional approach. More recently Spyrakos [4] and Spyrakos and Beskos [5] d eveloped a general two-dimensional plane stress/plane strain formulation using a piecewise constant approximation in space and time of all the boundary variables. They used the response of the elastic plane to a unit rectangular impulse force of finite duration as fundamental solution. The formulation presented by Mansur problems under zero initial conditions is [6] and Mansur and Brebbia [7] for wave propagation more general. This formulation is based on weighted residual as are other BEM formulations for different kinds of problems [8]. The fundamental solution corresponds to the singular unit impulse point force. The interpolation function for the boundary variables can be of any order in space and time. Antes [9] generalized the formulation to transient elastodynamic problems with arbitrary initial conditions. In the present paper the formulation according to Mansur [6] and Antes [9] is summarized as shown by Gallego first. Spyrakos approach [4] can be included in this general formulation and Dominguez [lo]. In fact, Spyrakos assumed a piecewise constant time interpolation function for both tractions and displacements whereas Mansur and Antes developed the case of piecewise constant time variation for the tractions and piecewise linear for the displacements. Results obtained using both approximations in time and assuming constant space interpolation functions are compared in the present paper. Quadratic space interpolation functions are also used, for the fist time to the authors’ knowledge, in combination with the transient elastodynamic BEM formulation. A comparison is made between results obtained with constant space iuterpolation functions, used in previous research, and those obtained with quadratic space interpolation functions. The use of the latter functions is shown to be necessary for an accurate solution of some problems. An analysis of the range of time increments for which the results are accurate is done for The effect of boundary discretizations constant and quadratic space interpolation functions. containing elements of variable sizes combined with that of the length of the time increment is also studied.
The authors would like to express their gratitude to the Spanish Con&i& Interministerial
de Ciencia y Tecnologia
for supporting this work under the research grant No. PB8G0139. Typeset by AM-?&Ii
119 MCM
15:3/5-I
120
J. DOMINGUEZ,
BASIC
FL.GALLEGO
.
EQUATIONS
The displacement along the direction “j” of an infinite medium due to a unit impluse in the direction ‘5” located at point t and acting at time r (bj = 6ijb(t - r) 6(~ - 0) is given by [ll] 1 li(clt’ - r) 2c:P2 - r2 = Jr,i r,j -RI 6ij 2np 1 cl r2 [ 1 H(C# - r) 2c2tt2 - r2 2 -r,ir,j (:I+ $) bij]) r2 C2 R2 [
urj = ut(+;[)
=
H(ql’
-
+J!?)
‘3
+
H(c2t’
-
&?’
‘I
=
“Ifj(l)
+
$P
‘J
(1)
’
where cl and c2 are the dilational and shear wave velocities, p is the density, R,
1
t’
=
t--r;
r =
1x-t N
1; u
= (Cit’2 - r2)li2 and H is the Heaviside function. The integral representation for the displacement 2 of a point on the boundary at time t may
be written as [ll] CijUj
=
.(Utj
* pj
-
Pi’j
* Uj)
dl?
J
+
J(U:j
* bj)&+
PJn(~j~UL,-ZliaUC*)dR,
(2)
a
where uj and pi stand for the “j” component of the displacement and traction, respectively, “*” denotes time convolution, pt is the traction tensor obtained by differentiation of U~j and the last integral contains the initial conditions. From now on zero initial conditions and zero body forces bj will be assumed. No generality will be lost because of these assumptions. In this case equation (2) b ecomes
*pj -Pt J (u~j
*uj)dI’. (3) r The computation of pzj requires a transformation of the spatial derivatives because, as pointed out by Mansur [G] and by Antes [9], the spatial derivatives of equation (1) give rise to the product of two singular functions that are not suitable for numerical solution. Nevertheless, the spatial derivatives can be related to time derivatives [G,9] and, after a process of integration by parts, an integral representation which can be numerically treated without special difficulty is obtained. Thus, equation (3) becomes CijUj =
(4) where
=
w;i(l) +
w;p.
The upper limit t+ 1s . used to avoid ending the integration at the peak of a delta function [9]. The above equation (4) is the same as that given by Antes [9] except for the initial conditions and the body force terms which have been assumed to be zero in the present paper. After the transformation of the spatial derivatives, the discretization and the integration of equation (4) are accomplished. Interpolation functions of any order can be used.
The time domain boundary
BOUNDARY Displacements tions
and tractions
along
element method
121
APPROXIMATION
the boundary
are approximated
using
interpolation
func-
(5)
where uyq and pyq stand for the j displacement and traction, respectively, of the node Q at the time t ,,, = mAt. The functions &‘(r) and @( T) are space interpolation functions and C(r) and vm(r) time interpolation functions. First the spatial interpolation functions @(r) and @( r ) are assumed to be piecewise constant, the time interpolation function pm(r) piecewise constant and V”‘(T) piecewise linear. After the approximation, equation (4) for the boundary node p at the time step n takes the form:
where Q is the total number of boundary nodes. The first time integral on the right-hand side of equation
+ J +m
urj/.imdr Twz-1
=
+
*, (u$) J 7*-1
+ u;j2))/irnd7
(6) can be written
= u;m(l)
as
+ U,Tf.mf2),
U?m(l)and Ulfmc2)are obtained by analytical integration of 2litj given by equation Explicit yxpressions & those functions may be found in references [9]. The second time integral of equation (6) is written following a similar notation as
where
J TA
(&Fj,,
w;j j”)&
-
=
(1).
pi;m(l) + P,;m(2),
Tm-1
pffmcl)and p.?m(2)can be computed by analytical integration and algebraical condensation. Their explicit ‘ixpressions are also given in reference [9]. Notice that the dependence of UC” and Pi;” on n and m is only on the difference between n and m since urj depends on the time difference t’ = t - T. The boundary element equation can now be written as where
4jzlj”p
=
2
5 { [L (~;~(l)
m=l
q=l
which can be written
+ $m(2)) dr]
I)j”’ -
[@;m’l’
+ P;m(2’)dr]
P
uy’)
I
(9)
as :
cfjuy = f: m=l
5 {G;“‘Pqp~q -
&JmPqu~‘J} ,
q=l
where Ginjtnpq= Jrn
jr*;mpq=
(u*;““(l) + qmq dr,
Jr (p;“‘(‘) + p;m(2))
dr.
0
(10)
J. DOMINGIJEZ, R. GALLECO
122
Spyrakos [4] and Spyrakos a.nd Beskos [5] assumed in their formulation a piecewise constant variation for the space and the time interpolation of both tractions and displacements; i.e., @(r), tiq(r),
pm(r)
P iecewise constant.
and ?(r)
The time integral
of equation
U?.“(2)) remains the same as before. The integral of equation (8) (P”” cz also be computed following the above procedure for the case in w&h piecewise constant.
P$”
r,+ (z$?p
=
- +j”)
&- =
(7) (VG” = utyrn(‘) + = Piyrn(‘) + PZyrn(‘)) 77”’ is assumed to be
p.y(l) + p;m@),
(11)
J Tnr-1
where
vrn
= 1,
when
r,,,_i
= 0,
otherwise,
5 r 5 r,,,
and
By use of the definitions
of zd”’ and u$!“) (equation
4) followed by an analytical
integration
one
obtains the Pnm’“) kernels. Their explicit expressions may be found in reference [lo]. The Ucm ij and Pi;” kernels obtained by the authors using piecewise constant time interpolation functions for the displacements and the tractions are different to some extent from those presented by Spyrakos [4]. Nevertheless, when instead of following the above procedure to obtain PG” the authors followed the procedure proposed by Spyrakos [4], i.e., doing the spatial differentiation of the displacement kernel UG” to compute directly the traction kernel PGm, the result was the same as before [lo]. INTEGRATION
ALONG
THE
ELEMENTS
The integration along the boundary elements indicated by equation (9) is done analytically when n = m for the element or elements to which the collocation point belongs. In all the other cases a Gauss quadrature with ten points is used. When the P or S-wave produced by the point load only reaches part of an element, the integration domain for that element is only the part of it that has been reached by the wave. The ten Gauss points are distributed over that part of the element. NUMERICAL
APPLICATIONS
In order to evaluate the relative advantage of different space and time interpolation functions, as well as to analyse the effect of the size of the time step, two elastodynamic problems have been studied using the above formulation. The first problem corresponds to an infinite strip under a uniform traction applied as a step function at time 1 = 0 (Figure 1). The material properties are: Poisson’s ratio v = 0.25, shear modulus /J = 4. lo4 Pa, P-wave velocity C, = 346.41 m/se and Swave velocity C’s = 200 m/se. The second problem refers to the behaviour of a rectangular plate under a flexural load with a triangular time variation (Figure 2). The material properties are the same as in the first problem. Both problems have been studied using the same rectangular geometry and the four space discretizations shown in Figure 3: (u) constant elements of the same size, (b) quadratic elements of the same size, (c) constant elements of variable sizes and (d) quadratic elements of variable sizes. The cases studied for the first problem are summarized in Table 1. The normal tractions at the clamped end (node 2 of mesh a) are shown in Figure 4 for cases 1 to 5. The normal displacement at the free end (node 14) is also shown in the same figure for case 3. The exact values are represented by a dashed line. It can be said that none of the cases 1 to 5 produces accurate values of the tractions at the clamped end. Only when P 2 1 the solution
The time domain boundary
element method
123
--k-f(t)
1
TINE
Figure
1.
Figure 2.
u0 I
t
N14
t
i
%
I4s 4i
1
L
L
$‘ L
5
Mesh
(a)
Mesh
(b)
Mesh
Figure 3.
N
(u)
Mesh
(d)
J. DOMINGUEZ, R. GALLEGO
124
Table 1.
Case
Space
Time
Time
interpolation
interpolation
interpolation
function
function for p
function for
Boundary
p
/3 = f&At/L
mesh
1
constant
constant
constant
a
2
constant
constant
constant
a
0.5
3
constant
constant
constant
a
1
4
constant
constant
constant
a
1.5
5
constant
constant
constant
a
2 0.25
0.25
6
constant
linear
constant
a
7
constant
linear
constant
a
0.5
8
constant
linear
constant
a
1
9
constant
linear
constant
a
1.5
10
constant
linear
constant
a
2 0.25
11
quadratic
linear
constant
b
12
quadratic
linear
constant
b
0.5
13
quadratic
linear
constant
b
1
14
quadratic
linear
constant
b
1.5
15
quadratic
linear
constant
b
2
16
constant
linear
constant
C
0.5
17
constant
linear
constant
C
1
18
constant
linear
constant
C
1.5
19
quadratic
linear
constant
d
0.5 1
20
quadratic
linear
constant
d
21
quadratic
linear
constant
d
1.5
22
constant
linear
constant
a
1
23
quadratic
linear
constant
b
1
L = distance between nodes for meshes a and b L = distance indicated in Figure 3 for meshes c and d.
is stable for some time, but in those cases the numerical solution is not able to follow the jumps of the exact solution. The computed displacements for the free end are more accurate, as could be expected from the fact that the exact solution is continuous for this variable. When the time variation for the displacement is assumed to be piecewise linear (cases 6 to 10) the solution improves (Figure 5). The normal tractions at the clamped end for p = 0.5 (case 7) are rather accurate but the solution becomes unstable after a short period of time. When /3 2 1 (cases 8 to 10) the solution remains stable for all the time range represented. The computed solution has more difficulties to follow the jumps of the exact solution and becomes smoother as /3 grows. The displacements of the free end are well represented for values of p > 1 (cases 9 and 10) and also for p 5 1 (cases 6 to 8) as long as the solution remains stable. The normal displacement of the node 14 is shown in the figure for case 8. As can be seen in Figure 5, the numerical solution of stresses and displacements for cases 9 and 10 (p = 1.5 and p = 2) present a small amount of damping as time goes on. Some improvement of the solution is obtained when space quadratic elements are used. Figure 6 shows the results for cases 11 to 15. The normal tractions at the centre of the clamped end have clearly a better behaviour in cases 11 and 12 (/? = 0.25 and p = 0.5) than in cases 6 and 7. However the improvement is small for cases 13 to 15 (p 1 1) as compared to cases 8 to 10. The solutions in cases 14 and 15 show less damping as time goes on than the solution obtained with constant space interpolation functions (cases 9 and 10). On the other hand, the space quadratic elements present more oscillations in the constant part of the normal stress plot than the space constant elements. The normal displacement of node 15 is shown versus time in Figure 6 for case 13 (p = 1). The non-uniform meshes (c) and (cl) are used for cases 16 to 18 and 19 to 21, respectively. These meshes have been used not because they were necessary for the problem at hand but
The time domain boundary element method
O
?a
125
*
0
0
lb
TIME (microseconds)
TIME (microseconds) *1
II
0
u
*
100
0
a
u
(0
P
16
TIME (microseconds)
TIME (microseconds)
US?,
TIME (microseconds)
TIME (microseconds)
because meshes like those are necessary geometries.
for problems
with more complicated
boundary
conditions
Figure 7 shows the normal traction at node 5 for the cases 16 to 21. The lower and the higher values of ,B used before have been dropped because good results could not be expected for those cases. The effect of the increase in the value of p is similar to that observed for the uniform meshes. It should be noticed that for a given ,# the equivalent values of ,B for the smaller and the larger elements of the mesh are 4p and p/2, respectively. A comparison between cases 16 to 18 and 19 to 21 shows that the results are similar for space constant and space quadratic elements and no need for the more complicated quadratic elements is apparent. It was also shown by comparison of Figures 5 and 6 that when a uniform mesh
126
J. DOMINGUEZ, Ft. GALLEGO
TIME (micrmccon&)
Figure 5.
is used and p 2 1 the space constant elements give good results. The question is whether a spatial interpolation function better than piecewise constant is needed. To answer this question the behaviour of a rectangular plate under flexural loads is studied using the above formulation. Two cases are analysed (‘22 and 23). The uniform boundary discretization with space constant elements (mesh a) is used in case 22 whereas space quadratic boundary elements are used for case 23 (mesh b). A step size p = CpAt/L = 1, piecewise linear time variation for the displacements and piecewise constant time variation for the tractions are assumed in both cases. The tangential displacements of node 15 for cases 22 and 23 are compared in Figure 8 with those obtained using a Finite Element mesh with 440 degrees of freedom [12]. It can be seen from Figure 8 that the boundary element results obtained using the spatial quadratic discretization
The time domain boundary
.
.
.
.
.
.
.
.
I.
element method
127
I
TIME (microseconds)
TIME (microsccon&) OI 0 z
0
a
40
I
P
“-
IP
TIME (microseconds) Figure
6.
are in excellent agreement with the finite element results. On the other hand, the space constant boundary elements solution shows significant diIferences from the other two. CONCLUSIONS From the above analysis the following conclusions and recommendations can be derived: 1. The use of a piecewise linear time variation of the displacements produces important improvements of the results as compared to those obtained using piecewise constant displacements. The computational effort is in both cases the same. 2. The use of time steps ,B close to 1 is recommended for problems with uniform meshes. Very small time steps tend to give unstable solutions and very large time steps give results that
128
.I. DOMINGIJEZ,Ft.GALLEGO
TIME (nlicroseconds) w
I.
t
40
20
IO
m
Irn
TIME(mhosecon&)
TIME
TIME (&rcsecon&)
(microseconds)
Fistwe7.
cannot represent well quick changes of the exact value /3 = 1 should be referred to an intermediate
solution. When size element.
the mesh is not uniform
the
of the body may well be solved using piece3. Problems that do not present flexural deformation wise constant spatial approximation for the displacements and the tractions. of the body require the use of higher order spatial 4. Problems that present flexural deformation interpolation functions. The use of space quadratic interpolation functions has been shown to produce accurate results. This conclusion is consistent with a similar one that was obtained for elastostatics some time ago [$I. To the authors knowledge, quadratic spatial interpolation functions in the context of time boundary elements have been used in the present paper for the first time.
The
t.ime domain
0
boundary
m
40
element
fm
cm
m
method
129
loo
140
lm
w-4
Figure
8.
REFERENCES
1. D.M. Cole, D.D. Kosloff S&m. 2.
Y.
Sec. Am.
Niwa,
T.
Fukui,
two-dimensional 3.
G. Manolis,
A comparative
4. C.C. Spyrakos, University 5.
C.C.
and K. study
Meth.
Dynamic
response
of Minnesota, and D.E.
method,
Inl.
Fuji,
Theor.
J. Num.
Spyrakos
element
S. Kato
A numerical
boundary
integral
on three
Eng.
Beskos,
application Mech.
Dynamic
Meth.
of Tokyo
element
and J. Dominquez, and McGraw-Hill,
response
of rigid
Eng. 23, 1547-1565
R. Gallego
and J. Dominguez,
App.
Num.
11.
A.C.
Eringen
12.
D. Nardini and C.A. Elements
Meth.,
(in press)
Boundary
Elements.
Southampton,
New
A unified 6, 17-25
and E.S. Suhubi,
(Edited
Brebbia,
by C.A.
method
to
(1980).
to problems
BEM-FEM
method,
in elasto-
Ph.D.
Thesis,
strip-foundations problems
formulation
boundary
using the boundary
using a time-stepping Verlag, (1983).
An Introductory
York,
by a time-domain
Course,
technique, Computational
element
In Eovndary Mechanics
(1989).
of two existing
in tw*dimensional
boundary
element
isotropic elastic
approaches,
Comm.
(1990).
Elaslodynamics
Vol. II, Linear
Transient dynamic
Brebbia),
equation
281-290
(1986).
9. H. Antes, A boundary element procedure for transient wave propagations media, Finite E/em. Anal. Des. 1, 313-322 (1985). 10.
integral 28,
approaches
by the time domain
W.J. Mansur and C.A. Brebbia, Transient elastodynamics Elements (Edited by C.A. Brebbia), pp. 667-698, Springer Brebbia
Bull.
(1984).
7.
C.A.
Press.
method
W.J. Mansur, A time-stepping technique to solve wave propagation method, Ph.D. Thesis, University of Southampton, U.K., (1983).
Publications
for elastodynamics,
(1983).
of strip-foundations Minn.,
of the boundary
Univ.
boundary
19, 73-91
Minneapolis,
J. Num.
An
Appl.
6.
8.
method
(1978).
elastodynamics,
Int.
dynamics,
and J.B. Minster,
68, 1331-1357
pp. 719-i30,
Theory,
Academic
analysis by the boundary Springer
Verlag,
Berlin,
Press, NY.,
element method, (1983).
(1975). In Boundary
hfodelling
Vol. 15, No. 3-5, pp. 119-129,
1991
Printed in Great Britain. All rights reserved
THE
TIME DOMAIN BOUNDARY FOR ELASTODYNAMIC
0895-7177/91 $3.90 + 0.00 Copyright@ 1991 Pergamon Press plc
ELEMENT PROBLEMS
METHOD
Jo& DOMINGUEZ AND RAFAEL GALLEGO Escuela Tkcnica Superior de Ingenieros Industriales Universidad de Sevilla, Sevilla, Spain
INTRODUCTION The number of papers published on the time domain formulation and implementation of the Boundary Element Method (BEM) is rather small. Cole e2 al. [l] were the first to present the formulation for two-dimensional antiplane problems. Plane strain problems were studied by Niwa et al. [2] and Manolis [3] using a three-dimensional approach. More recently Spyrakos [4] and Spyrakos and Beskos [5] d eveloped a general two-dimensional plane stress/plane strain formulation using a piecewise constant approximation in space and time of all the boundary variables. They used the response of the elastic plane to a unit rectangular impulse force of finite duration as fundamental solution. The formulation presented by Mansur problems under zero initial conditions is [6] and Mansur and Brebbia [7] for wave propagation more general. This formulation is based on weighted residual as are other BEM formulations for different kinds of problems [8]. The fundamental solution corresponds to the singular unit impulse point force. The interpolation function for the boundary variables can be of any order in space and time. Antes [9] generalized the formulation to transient elastodynamic problems with arbitrary initial conditions. In the present paper the formulation according to Mansur [6] and Antes [9] is summarized as shown by Gallego first. Spyrakos approach [4] can be included in this general formulation and Dominguez [lo]. In fact, Spyrakos assumed a piecewise constant time interpolation function for both tractions and displacements whereas Mansur and Antes developed the case of piecewise constant time variation for the tractions and piecewise linear for the displacements. Results obtained using both approximations in time and assuming constant space interpolation functions are compared in the present paper. Quadratic space interpolation functions are also used, for the fist time to the authors’ knowledge, in combination with the transient elastodynamic BEM formulation. A comparison is made between results obtained with constant space iuterpolation functions, used in previous research, and those obtained with quadratic space interpolation functions. The use of the latter functions is shown to be necessary for an accurate solution of some problems. An analysis of the range of time increments for which the results are accurate is done for The effect of boundary discretizations constant and quadratic space interpolation functions. containing elements of variable sizes combined with that of the length of the time increment is also studied.
The authors would like to express their gratitude to the Spanish Con&i& Interministerial
de Ciencia y Tecnologia
for supporting this work under the research grant No. PB8G0139. Typeset by AM-?&Ii
119 MCM
15:3/5-I
120
J. DOMINGUEZ,
BASIC
FL.GALLEGO
.
EQUATIONS
The displacement along the direction “j” of an infinite medium due to a unit impluse in the direction ‘5” located at point t and acting at time r (bj = 6ijb(t - r) 6(~ - 0) is given by [ll] 1 li(clt’ - r) 2c:P2 - r2 = Jr,i r,j -RI 6ij 2np 1 cl r2 [ 1 H(C# - r) 2c2tt2 - r2 2 -r,ir,j (:I+ $) bij]) r2 C2 R2 [
urj = ut(+;[)
=
H(ql’
-
+J!?)
‘3
+
H(c2t’
-
&?’
‘I
=
“Ifj(l)
+
$P
‘J
(1)
’
where cl and c2 are the dilational and shear wave velocities, p is the density, R,
1
t’
=
t--r;
r =
1x-t N
1; u
= (Cit’2 - r2)li2 and H is the Heaviside function. The integral representation for the displacement 2 of a point on the boundary at time t may
be written as [ll] CijUj
=
.(Utj
* pj
-
Pi’j
* Uj)
dl?
J
+
J(U:j
* bj)&+
PJn(~j~UL,-ZliaUC*)dR,
(2)
a
where uj and pi stand for the “j” component of the displacement and traction, respectively, “*” denotes time convolution, pt is the traction tensor obtained by differentiation of U~j and the last integral contains the initial conditions. From now on zero initial conditions and zero body forces bj will be assumed. No generality will be lost because of these assumptions. In this case equation (2) b ecomes
*pj -Pt J (u~j
*uj)dI’. (3) r The computation of pzj requires a transformation of the spatial derivatives because, as pointed out by Mansur [G] and by Antes [9], the spatial derivatives of equation (1) give rise to the product of two singular functions that are not suitable for numerical solution. Nevertheless, the spatial derivatives can be related to time derivatives [G,9] and, after a process of integration by parts, an integral representation which can be numerically treated without special difficulty is obtained. Thus, equation (3) becomes CijUj =
(4) where
=
w;i(l) +
w;p.
The upper limit t+ 1s . used to avoid ending the integration at the peak of a delta function [9]. The above equation (4) is the same as that given by Antes [9] except for the initial conditions and the body force terms which have been assumed to be zero in the present paper. After the transformation of the spatial derivatives, the discretization and the integration of equation (4) are accomplished. Interpolation functions of any order can be used.
The time domain boundary
BOUNDARY Displacements tions
and tractions
along
element method
121
APPROXIMATION
the boundary
are approximated
using
interpolation
func-
(5)
where uyq and pyq stand for the j displacement and traction, respectively, of the node Q at the time t ,,, = mAt. The functions &‘(r) and @( T) are space interpolation functions and C(r) and vm(r) time interpolation functions. First the spatial interpolation functions @(r) and @( r ) are assumed to be piecewise constant, the time interpolation function pm(r) piecewise constant and V”‘(T) piecewise linear. After the approximation, equation (4) for the boundary node p at the time step n takes the form:
where Q is the total number of boundary nodes. The first time integral on the right-hand side of equation
+ J +m
urj/.imdr Twz-1
=
+
*, (u$) J 7*-1
+ u;j2))/irnd7
(6) can be written
= u;m(l)
as
+ U,Tf.mf2),
U?m(l)and Ulfmc2)are obtained by analytical integration of 2litj given by equation Explicit yxpressions & those functions may be found in references [9]. The second time integral of equation (6) is written following a similar notation as
where
J TA
(&Fj,,
w;j j”)&
-
=
(1).
pi;m(l) + P,;m(2),
Tm-1
pffmcl)and p.?m(2)can be computed by analytical integration and algebraical condensation. Their explicit ‘ixpressions are also given in reference [9]. Notice that the dependence of UC” and Pi;” on n and m is only on the difference between n and m since urj depends on the time difference t’ = t - T. The boundary element equation can now be written as where
4jzlj”p
=
2
5 { [L (~;~(l)
m=l
q=l
which can be written
+ $m(2)) dr]
I)j”’ -
[@;m’l’
+ P;m(2’)dr]
P
uy’)
I
(9)
as :
cfjuy = f: m=l
5 {G;“‘Pqp~q -
&JmPqu~‘J} ,
q=l
where Ginjtnpq= Jrn
jr*;mpq=
(u*;““(l) + qmq dr,
Jr (p;“‘(‘) + p;m(2))
dr.
0
(10)
J. DOMINGIJEZ, R. GALLECO
122
Spyrakos [4] and Spyrakos a.nd Beskos [5] assumed in their formulation a piecewise constant variation for the space and the time interpolation of both tractions and displacements; i.e., @(r), tiq(r),
pm(r)
P iecewise constant.
and ?(r)
The time integral
of equation
U?.“(2)) remains the same as before. The integral of equation (8) (P”” cz also be computed following the above procedure for the case in w&h piecewise constant.
P$”
r,+ (z$?p
=
- +j”)
&- =
(7) (VG” = utyrn(‘) + = Piyrn(‘) + PZyrn(‘)) 77”’ is assumed to be
p.y(l) + p;m@),
(11)
J Tnr-1
where
vrn
= 1,
when
r,,,_i
= 0,
otherwise,
5 r 5 r,,,
and
By use of the definitions
of zd”’ and u$!“) (equation
4) followed by an analytical
integration
one
obtains the Pnm’“) kernels. Their explicit expressions may be found in reference [lo]. The Ucm ij and Pi;” kernels obtained by the authors using piecewise constant time interpolation functions for the displacements and the tractions are different to some extent from those presented by Spyrakos [4]. Nevertheless, when instead of following the above procedure to obtain PG” the authors followed the procedure proposed by Spyrakos [4], i.e., doing the spatial differentiation of the displacement kernel UG” to compute directly the traction kernel PGm, the result was the same as before [lo]. INTEGRATION
ALONG
THE
ELEMENTS
The integration along the boundary elements indicated by equation (9) is done analytically when n = m for the element or elements to which the collocation point belongs. In all the other cases a Gauss quadrature with ten points is used. When the P or S-wave produced by the point load only reaches part of an element, the integration domain for that element is only the part of it that has been reached by the wave. The ten Gauss points are distributed over that part of the element. NUMERICAL
APPLICATIONS
In order to evaluate the relative advantage of different space and time interpolation functions, as well as to analyse the effect of the size of the time step, two elastodynamic problems have been studied using the above formulation. The first problem corresponds to an infinite strip under a uniform traction applied as a step function at time 1 = 0 (Figure 1). The material properties are: Poisson’s ratio v = 0.25, shear modulus /J = 4. lo4 Pa, P-wave velocity C, = 346.41 m/se and Swave velocity C’s = 200 m/se. The second problem refers to the behaviour of a rectangular plate under a flexural load with a triangular time variation (Figure 2). The material properties are the same as in the first problem. Both problems have been studied using the same rectangular geometry and the four space discretizations shown in Figure 3: (u) constant elements of the same size, (b) quadratic elements of the same size, (c) constant elements of variable sizes and (d) quadratic elements of variable sizes. The cases studied for the first problem are summarized in Table 1. The normal tractions at the clamped end (node 2 of mesh a) are shown in Figure 4 for cases 1 to 5. The normal displacement at the free end (node 14) is also shown in the same figure for case 3. The exact values are represented by a dashed line. It can be said that none of the cases 1 to 5 produces accurate values of the tractions at the clamped end. Only when P 2 1 the solution
The time domain boundary
element method
123
--k-f(t)
1
TINE
Figure
1.
Figure 2.
u0 I
t
N14
t
i
%
I4s 4i
1
L
L
$‘ L
5
Mesh
(a)
Mesh
(b)
Mesh
Figure 3.
N
(u)
Mesh
(d)
J. DOMINGUEZ, R. GALLEGO
124
Table 1.
Case
Space
Time
Time
interpolation
interpolation
interpolation
function
function for p
function for
Boundary
p
/3 = f&At/L
mesh
1
constant
constant
constant
a
2
constant
constant
constant
a
0.5
3
constant
constant
constant
a
1
4
constant
constant
constant
a
1.5
5
constant
constant
constant
a
2 0.25
0.25
6
constant
linear
constant
a
7
constant
linear
constant
a
0.5
8
constant
linear
constant
a
1
9
constant
linear
constant
a
1.5
10
constant
linear
constant
a
2 0.25
11
quadratic
linear
constant
b
12
quadratic
linear
constant
b
0.5
13
quadratic
linear
constant
b
1
14
quadratic
linear
constant
b
1.5
15
quadratic
linear
constant
b
2
16
constant
linear
constant
C
0.5
17
constant
linear
constant
C
1
18
constant
linear
constant
C
1.5
19
quadratic
linear
constant
d
0.5 1
20
quadratic
linear
constant
d
21
quadratic
linear
constant
d
1.5
22
constant
linear
constant
a
1
23
quadratic
linear
constant
b
1
L = distance between nodes for meshes a and b L = distance indicated in Figure 3 for meshes c and d.
is stable for some time, but in those cases the numerical solution is not able to follow the jumps of the exact solution. The computed displacements for the free end are more accurate, as could be expected from the fact that the exact solution is continuous for this variable. When the time variation for the displacement is assumed to be piecewise linear (cases 6 to 10) the solution improves (Figure 5). The normal tractions at the clamped end for p = 0.5 (case 7) are rather accurate but the solution becomes unstable after a short period of time. When /3 2 1 (cases 8 to 10) the solution remains stable for all the time range represented. The computed solution has more difficulties to follow the jumps of the exact solution and becomes smoother as /3 grows. The displacements of the free end are well represented for values of p > 1 (cases 9 and 10) and also for p 5 1 (cases 6 to 8) as long as the solution remains stable. The normal displacement of the node 14 is shown in the figure for case 8. As can be seen in Figure 5, the numerical solution of stresses and displacements for cases 9 and 10 (p = 1.5 and p = 2) present a small amount of damping as time goes on. Some improvement of the solution is obtained when space quadratic elements are used. Figure 6 shows the results for cases 11 to 15. The normal tractions at the centre of the clamped end have clearly a better behaviour in cases 11 and 12 (/? = 0.25 and p = 0.5) than in cases 6 and 7. However the improvement is small for cases 13 to 15 (p 1 1) as compared to cases 8 to 10. The solutions in cases 14 and 15 show less damping as time goes on than the solution obtained with constant space interpolation functions (cases 9 and 10). On the other hand, the space quadratic elements present more oscillations in the constant part of the normal stress plot than the space constant elements. The normal displacement of node 15 is shown versus time in Figure 6 for case 13 (p = 1). The non-uniform meshes (c) and (cl) are used for cases 16 to 18 and 19 to 21, respectively. These meshes have been used not because they were necessary for the problem at hand but
The time domain boundary element method
O
?a
125
*
0
0
lb
TIME (microseconds)
TIME (microseconds) *1
II
0
u
*
100
0
a
u
(0
P
16
TIME (microseconds)
TIME (microseconds)
US?,
TIME (microseconds)
TIME (microseconds)
because meshes like those are necessary geometries.
for problems
with more complicated
boundary
conditions
Figure 7 shows the normal traction at node 5 for the cases 16 to 21. The lower and the higher values of ,B used before have been dropped because good results could not be expected for those cases. The effect of the increase in the value of p is similar to that observed for the uniform meshes. It should be noticed that for a given ,# the equivalent values of ,B for the smaller and the larger elements of the mesh are 4p and p/2, respectively. A comparison between cases 16 to 18 and 19 to 21 shows that the results are similar for space constant and space quadratic elements and no need for the more complicated quadratic elements is apparent. It was also shown by comparison of Figures 5 and 6 that when a uniform mesh
126
J. DOMINGUEZ, Ft. GALLEGO
TIME (micrmccon&)
Figure 5.
is used and p 2 1 the space constant elements give good results. The question is whether a spatial interpolation function better than piecewise constant is needed. To answer this question the behaviour of a rectangular plate under flexural loads is studied using the above formulation. Two cases are analysed (‘22 and 23). The uniform boundary discretization with space constant elements (mesh a) is used in case 22 whereas space quadratic boundary elements are used for case 23 (mesh b). A step size p = CpAt/L = 1, piecewise linear time variation for the displacements and piecewise constant time variation for the tractions are assumed in both cases. The tangential displacements of node 15 for cases 22 and 23 are compared in Figure 8 with those obtained using a Finite Element mesh with 440 degrees of freedom [12]. It can be seen from Figure 8 that the boundary element results obtained using the spatial quadratic discretization
The time domain boundary
.
.
.
.
.
.
.
.
I.
element method
127
I
TIME (microseconds)
TIME (microsccon&) OI 0 z
0
a
40
I
P
“-
IP
TIME (microseconds) Figure
6.
are in excellent agreement with the finite element results. On the other hand, the space constant boundary elements solution shows significant diIferences from the other two. CONCLUSIONS From the above analysis the following conclusions and recommendations can be derived: 1. The use of a piecewise linear time variation of the displacements produces important improvements of the results as compared to those obtained using piecewise constant displacements. The computational effort is in both cases the same. 2. The use of time steps ,B close to 1 is recommended for problems with uniform meshes. Very small time steps tend to give unstable solutions and very large time steps give results that
128
.I. DOMINGIJEZ,Ft.GALLEGO
TIME (nlicroseconds) w
I.
t
40
20
IO
m
Irn
TIME(mhosecon&)
TIME
TIME (&rcsecon&)
(microseconds)
Fistwe7.
cannot represent well quick changes of the exact value /3 = 1 should be referred to an intermediate
solution. When size element.
the mesh is not uniform
the
of the body may well be solved using piece3. Problems that do not present flexural deformation wise constant spatial approximation for the displacements and the tractions. of the body require the use of higher order spatial 4. Problems that present flexural deformation interpolation functions. The use of space quadratic interpolation functions has been shown to produce accurate results. This conclusion is consistent with a similar one that was obtained for elastostatics some time ago [$I. To the authors knowledge, quadratic spatial interpolation functions in the context of time boundary elements have been used in the present paper for the first time.
The
t.ime domain
0
boundary
m
40
element
fm
cm
m
method
129
loo
140
lm
w-4
Figure
8.
REFERENCES
1. D.M. Cole, D.D. Kosloff S&m. 2.
Y.
Sec. Am.
Niwa,
T.
Fukui,
two-dimensional 3.
G. Manolis,
A comparative
4. C.C. Spyrakos, University 5.
C.C.
and K. study
Meth.
Dynamic
response
of Minnesota, and D.E.
method,
Inl.
Fuji,
Theor.
J. Num.
Spyrakos
element
S. Kato
A numerical
boundary
integral
on three
Eng.
Beskos,
application Mech.
Dynamic
Meth.
of Tokyo
element
and J. Dominquez, and McGraw-Hill,
response
of rigid
Eng. 23, 1547-1565
R. Gallego
and J. Dominguez,
App.
Num.
11.
A.C.
Eringen
12.
D. Nardini and C.A. Elements
Meth.,
(in press)
Boundary
Elements.
Southampton,
New
A unified 6, 17-25
and E.S. Suhubi,
(Edited
Brebbia,
by C.A.
method
to
(1980).
to problems
BEM-FEM
method,
in elasto-
Ph.D.
Thesis,
strip-foundations problems
formulation
boundary
using the boundary
using a time-stepping Verlag, (1983).
An Introductory
York,
by a time-domain
Course,
technique, Computational
element
In Eovndary Mechanics
(1989).
of two existing
in tw*dimensional
boundary
element
isotropic elastic
approaches,
Comm.
(1990).
Elaslodynamics
Vol. II, Linear
Transient dynamic
Brebbia),
equation
281-290
(1986).
9. H. Antes, A boundary element procedure for transient wave propagations media, Finite E/em. Anal. Des. 1, 313-322 (1985). 10.
integral 28,
approaches
by the time domain
W.J. Mansur and C.A. Brebbia, Transient elastodynamics Elements (Edited by C.A. Brebbia), pp. 667-698, Springer Brebbia
Bull.
(1984).
7.
C.A.
Press.
method
W.J. Mansur, A time-stepping technique to solve wave propagation method, Ph.D. Thesis, University of Southampton, U.K., (1983).
Publications
for elastodynamics,
(1983).
of strip-foundations Minn.,
of the boundary
Univ.
boundary
19, 73-91
Minneapolis,
J. Num.
An
Appl.
6.
8.
method
(1978).
elastodynamics,
Int.
dynamics,
and J.B. Minster,
68, 1331-1357
pp. 719-i30,
Theory,
Academic
analysis by the boundary Springer
Verlag,
Berlin,
Press, NY.,
element method, (1983).
(1975). In Boundary