- Email: [email protected]
Further improvements on three dimensional transient BEM elastodynamic analysis
Engineering Analysis with Boundary Elements 17 (1996) 231-243
PII:
ELSEVIER
S0955-7997(96)00019-7
Copyright © 1996. Published by Elsevier Science Ltd Printed in Great Britain. All rights reserved 0955-7997/96/$15.00
Further improvements on three dimensional transient BEM elastodynamic analysis Humberto Breves Coda & Wilson Sergio Venturini* S~o Carlos School of Engineering, University of S~o Paulo, Av. Carlos Botelho 1465, 1356-970 S~o Carlos, Brazil
(Received 14 January 1994; revised version received 6 March 1996; accepted 20 May 1996) In this work the transient boundary element formulation is studied together with a modified time dependent fundamental solution. This modification permits one to obtain more stable results during numerical analysis. It also reduces the amount of matrices usually required to be stored. Furthermore, some tricks, usually needed when the classical fundamental solution is adopted, are now avoided, improving the result quality and eliminating some discretization restrictions. Copyright © 1996. Published by Elsevier Science Ltd Key words: Boundary element method, elastodynamics, time domain, transient
problem. INTRODUCTION
applied at a single instant. In this work, another assumption regarding the load choice is made: a unit impulse distributed along an interval of time was taken to derive the necessary fundamental values. The numerical algorithms developed from the standard fundamental solution described above, which has been used so far, required the use of some inconvenient tricks to deal with terms containing Dirac's delta derivatives. For instance, we can point out the introduction of techniques, usually adopted to compute velocity terms (backward finite differences), to allow the use of constant time approximation for displacements. 11 In order to improve the procedure, many authors have adopted linear approximation for displacements and constant approximation for tractions. These necessary tricks also increase the number of parameters to be computed for each time integration step. The alternative fundamental solution proposed here has been demonstrated to produce accurate results when analyzing practical examples. The use of backward finite differences has been avoided. The computing time is also reduced due to the smaller amount of parameters to be evaluated in comparison with the classical approach.
The first time domain formulation of the boundary element method (BEM) was proposed by Friedman and Shaw, 1 in 1962, in which the acoustic anti-plane problem was treated. In the context of two dimensional elastodynamic analysis, the first work was presented by Niwa et al., 2 in 1980; nevertheless it was only in 1983 that a well established time stepping technique was presented by Mansur and Brebbia. 3,4 Karabalis and Beskos, 5 in 1984, were the first to develop a general three-dimensional elastodynamic methodology in time domain BEM. Interesting overviews on BEM formulations to treat elastic transient problems can be seen in Beskos 6 and Mansur. 7 The BEM formulations, employed to study any physical problem, are derived from some well established integral equations. These equations are usually achieved by using fundamental solutions, which are obtained by solving the corresponding differential equation of the problem, for an infinite domain, where a convenient unit load is applied. For the case of elastic bodies, the fundamental solution comes from the Navier-Cauchy equilibrium equations. 8'9 For elastodynamic problems, the general singular solution is due to Stokes, 1° representing the behavior of an infinite body subjected to a time dependent concentrated load. The classical fundamental solution, used in the time domain BEM approach, is a particular case of the Stokes' state, assuming that the load is a unit concentrated impulse
FUNDAMENTAL SOLUTION Assuming that the bodies under analysis are isotropic and show linear elastic behavior, the governing differential equations for the dynamic equilibrium, the so-called
*To whom correspondence should be addressed. 231
H. B. Coda, W. S. Venturini
232 Navier-Cauchy equations, are given by
(C~ -- C2)u],ji + C2bli,j] + hi~ p : Hi
(1
Bik(q,v-;s,O)=
3C 2 5 r5
r3
j
where b i represents body forces, p is the media density, and Cl and C2 are longitudinal and shear wave propagation velocities, respectively. As usual, we can represent a time-dependent concentrated load in the k-direction, acting at any point's' of the infinite domain, by the following expression:
b*ki =.f(r)5(q -- s)6kg
+216F
(2)
rirjr k7
where q is a field point, ¢5ki is the Kronecker delta, 5 ( s - q ) the Dirac delta and f ( r ) a time-dependent load function. Replacing the body force term in eqn (1) according to eqn (2) gives, (C~
2 • 2 . . . - C2)uj,ji + C2ug,jj +f(r)~5(q -. S)6ki/P = ui
(3)
in which the symbol * indicates fundamental values. By integrating eqn (3) one obtains the so-called Stokes' fundamental solution,m easily found in the literature.12'13 The usual fundamental solutions, employed for all previous works mentioned here, are the particular case of the Stokes' state in which f ( r ) = 5(r) is assumed. This load represents a unit impulse at time r = 0. The u*kg and P*kg expressions, for this particular case, are found in several well known works. 12,14 Trying to obtain an alternative which could improve the quality of the numerical responses for BEM transient analysis, we decided to use another load function together with the Stokes' state, given by the following expression,
b*~i = [H(r) - H(r - At)]6(q - s)6ki/At
(4)
in which H( ) is the well known Heaviside function and At a constant time interval. This load no longer represents a unit impulse at an instant, as in the classical approach, but a unit impulse developed over a time interval At. Displacement and traction values, derived from this particular time dependent load, are given by,
u*ki(q, "r;s, O) = [A*ki(q, r ;s, O) - A*ki(q, r - At; s, 0)] (5) P*ki(q, 7-; s, O) = [B*ki(q,r; s, O) - B~i(q , 7- - At; s, 0)] (6) where Akg and B~g are:
A~i(q,r;s,o)--
l { (3rgrk 4rrp/Xt \ r 3
~5~g) 1 2r 2
rirjrk [ ~
+ 2r-CdTC
r)
C~5(r
+~6
7-r
r)]
\ -U,
rk~50[ 2 fC2"~1 [ H ( r r3 1 - ~C~] J
~-2) -
r3 r (5
r
Equations (7) and (8), when compared with the corresponding classical ones, show a smoother time behavior, i.e. the new solution does not exhibit the undesirable Dirac's delta derivatives. As a consequence, no velocity terms will appear when time convolutions are performed. As for the classical fundamental solution, the new one goes to zero when r>>.r/C2+At, which makes it possible to compute only a finite amount of coefficients.
INTEGRAL REPRESENTATIONS
The Graffi reciprocity theorem, for dynamic problems, can be written for any two elastodynamic states, here represented by their displacement fields u and u', as follows 14-t6
/'J [ui(q,r)b~(q,t
" H r 2 _ rC2
tSijrk + tSikr j r3 + 6jkri]
!
0
~
' t - r)bi(q, r ) ] d a d r - r) - ui(q,
f f [ui(q, r ) b : ( q , t - r ) - u i ( q , t - r)bi(q, r ) ] d ~ d r 3o J ~2 rkr"[~{H(r--~l + r3
) -- I~H C2
= I l l [ui(Q,t-r)pi(Q,r)-ui(Q, ' 3o3 I'
, - r)]dFdr r)pi(Q,t
+ p [ [tJ'i(q, t)ui(q, O) + u'i(q, t)ui(q, 0)]df~ •J~2
(9)
Further improvements on three dimensional transient B E M elastodynamic analysis where 9t is a bounded or unbounded region and P its contour. Without lost of generality, only the case of quiescent past elastodynamic states will be studied in this work. Thus, eqn (9) becomes, [t [ [ui(q, "r)bti(q, t _ T ) - u l ( q , t -
"r)bi(q, T ) ] d ~ d T
JO d f~
233
and P*ki and ~:i are the Dirac's delta fundamental solution. When constant time approximation is adopted, for both displacement and traction, the explicit terms P~,i and U~i become similar to those derived for the new solution. Nevertheless, one can easily realize that they are applied in different equations, (13) and (11) respectively, therefore they have different meanings.
= I'o Jr[u,(Q, t-r)pi(Q, r)-ui(Q, T)pi(Q, ' t-T)]dPdT (10) Assuming the new fundamental solution as one of those elastodynamic states, eqn (10) gives,
ki[~,S)Jt_At- u,ls ~
l d r = Jol* rUki(Q,t;s,T)pi(Q,T)dFd7 +
Jol ui(Q,r)pki(Q, *
t; s,r)drd~-
NUMERICAL ASPECTS As usual for boundary element formulations, eqn (11) can be treated numerically be adopting appropriate space and time discretizations. For the particular case of the time discretization, constant elements have been adopted. Thus, eqn (11) can be written as follows,
Cki(s)UiN,(s,t) = [
[to U*ki(Q,t;s,7)¢J(Q)~°dTdFP~
- [
[ P*k~(Q,t;s,r)¢~(Q)¢°drdrU~d
Jr(/) Jto 1
Uo
U*ki(q,t; s, T)bi(q, T)dfMT
to (ll)
where Cki is the classical free term of displacement integral representations computed as in Hartmann. 17 In eqn (11), one has to use U*kiand p~i, as defined in (5) and (6), taking into account the following property of the Stokes state, U*ki(Q, t - r ; s, O) = U*ki(Q, t; s, r )
(12a)
P*ki(Q, t - r ; s, O) = P*ki(Q, t; s, y)
(12b)
It is worth mentioning, that, when the classical fundamental solution (Dirac's delta load function) is used, as in Karabalis and Beskos, 5'18 eqn (10) reduces to Love's identity and the time convolution leads to functional operators belonging to the Stokes' state, l°J2J4 This expression can be written as:
- Ir P*ki(Q, t; s I ui)dFdT
where e*ki and U~¢i
(13)
given by,
are
j u*ki(q,t, s, r)pi(q, T)dT * Pki(Q,t;s I ui) = 0 fi*ki(q, t, s, T)ui(q, I' t
U~i(Q, t; s I pn) =
o
7")d~-
(14)
Jto_1
(15) where 0 ranges from 1 to the integration step number, Nt, and the summation rule is implied. For the new approach, the value of Uis,(S, t), in eqn (15), means an average value computed along the time step Nt. It is not a single value taken at an instant t = NtAt, as for the classical procedure. The space approximations of displacements and tractions are given by,
Uio(Q j, t) = ¢~)¢(0)V;~oj
(16a)
Pio(Q j, t) = A(JL/,(°) j ~c~ ~ P?iO
(16b)
in which i j 0 ¢, ~0
Cki(Q, S)Ui(S~ T) = [ U~i(Q, t; s I pi)drdT JF
+J;~ P*ki(q,t;slbi)dfMT
JF(/)
= = = = =
displacement or traction directions space elements time step shape function for the space element time shape function (¢e = 1)
Flat quadrangular elements with linear and quadratic approximations have been adopted for the space discretization. The singular integrals have been computed with a special algorithm 19 based on Kutt's quadrature, For non-singular integrals, the standard Gauss quadrature was adopted together with a well known element subdivision technique) l In order to integrate the time dependent kernels in eqn (15), the time shape function ¢ = 1 was adopted,
234
H. B. Coda, W. S. Venturini
giving,
"E
J k5 to_~
r21/,- C,')
-
r
- - -
dr
if t - to >>- r
J
-- ~[(12 r~_~(TA_IoI)_I(TA~[o.I)~_(TA3[;_I)/31
(17)
if t - to_ 1 ~ -~, ~ t -- to
0
,o , A t
C~,]
if ~r / > t -- to-I
i
dr=
t-r-
dXt
) if,
0
Jii , 1 6 (
t-r-
r~'~dr
={~T
C~J
0
if t - t o _ j >
i f ~ >~t - to
ift-t° if~>t-to
and 7 is equal to 1 or 2. Subtracting from the values given by eqns (17)-(19), the corresponding ones computed at the next time step, [to, to+l], the appropriate time coefficients of the new fundamental solution, eqns (5) and (6), are obtained. As for any boundary element formulation, the displacement integral representation is transformed into an algebraic equation, after carrying out all space and time integrals over the elements defined by both discretizations. Selecting an appropriate number of collocation points 's' to write their algebraic representations, a system of linear equations is achieved as follows: (21 )
where 0 = 1,2, 3 . . . . , Nt.
(19)
i Ol'?7
(20)
H°U ° = G ° P °
=
,-
where TA = t - r/CA
(18)
)t-to
Diagonal sub-matrices of h/~) are obtained by the following expression, h~) = h Nt ,s + C
(22)
In eqn (22) I]sNt is obtained by carrying out singular integrals only at the later time step, i.e. exactly where we are computing the unknowns. Analogously for the classical transient BEM techniques employed by other authors] ~'ls the maximum time step recommended to compute new matrices is given by N, -
dmax C2At
+ 1
(23)
where dmax is defined as the maximum possible length of the discretized solid. We have to point out that when the classical solution is used together with constant/constant (displacement/traction) time approximations, the achieved
Om
(,
f E = 1,1 x 10 5N/m 2 p = 2 kg/m 2 v=0
t
Fig. 1. Prismatic beam: geometry and loading.
235
Further improvements on three dimensional transient B E M elastodynamic analysis
/
.
employed to verify numerical models, consists of analyzing a prismatic rod subjected to an impact at the free end. The results obtained, by using the proposed formulation, are compared with the classical numerical solution obtained by Banerjee et al., 20 using Dirac's delta fundamental solution and constant time approximation. The geometry and the material properties, as well as the load time behavior, are given in Fig. 1. A total of 18 quadratic elements, as illustrated in Fig. 2, were used to discretize the boundary of this solid. The relaxed boundary conditions, 0-31= 0 and 0-32= 0 , were also assumed. The analysis has been performed for two different time step intervals: 0.00356s, as used in Ref. 20, and At = 0.00712 s. The results obtained, in terms of stresses at the center of this beam, are shown in Fig. 3. As can be seen (Fig. 3), for both selected time steps, the numerical values obtained follow the analytical solution very well. Even for a rather large time step, good quality numerical results are achieved. In this case, the solution does not show abrupt variation of stresses due to the fundamental solution characteristics, which gives averaged values for each time step. For the smaller time step, the sharp jumps in stresses, for t I = 0.017s and t2 = 0.01 s, are very well captured. In the second example, a three dimensional halfspace, when subjected to a strip load, is analyzed. The half-space and the load intensity, as well as its time behavior, are shown in Fig. 4. The results obtained by adopting the proposed formulation are compared with those taken from Mansur and Brebbia, 16 computed by solving the corresponding two dimensional problem, also using a time domain BEM formulation. The finite difference solution, taken from Tseng and Robinson, El is also used for this comparison. The following values, related to material properties, have been adopted to solve the problem: E = 2 0 0 k s i , v--0-15 and p = 6.885 gbfsE/in4.
/
Fig. 2. Prismatic beam: discretization.
system of equations, equivalent to (21) is represented by, (H l + HHI)U 1 - - G I p 1
(0-- 1)
H°U ° + [HH°U° - n n 0 - 1 U 0 - 1
G°p °
] =
(0 = 2, 3,... N,)
(24)
In the same way, if linear/constant time approximation were adopted, one obtains, o 0+ HIU1
o 0
H2U2
Gop0
:
(25)
In eqn (24), HH ° contains the coefficients related with the backward difference scheme, as shown by Manolis et al. I1 On the other hand, eqn (25) exhibits two matrices, H°l and H °, achieved when linear approximation is used for displacements. 15 Thus, the new approach avoids dealing with finite differences and reduces the global amount of computing.
EXAMPLES Four examples have been selected to show the performance of this new fundamental solution, when used appropriately together with the BEM formulation. The first example, a classical problem usually
1
ANALYTICAl.
o BANERJEE + THIS WORK • THIS WORK
P ( N/m z ) .,...4,,222,,,.m
2000-
| At-O.OO356s) ( At, OOO366s| ( At,O.OOT12 s)
°
•
0 0
I000-
t
~
w
"
"
w ~r 0
-
4" v 0
t l 4. ÷
!
0
. Z
I
T 4
,
, 6
I
l II
! I;21 10
l I l , l , 14 IQ 111
I , ZO
l I l Im |4
I
I l , I ~ l 26 all
Fig. 3. Normal stress at beam middle point.
,
32
•
,
N
t ( 3.56 K)'ss )
236
H. B. Coda, W. S. Venturini
rt"
lksi
,
~,.0 m i
,,
/
///////
TIME
,/'/'//
/ /
11
/ / / / / / / / / / / / / / / / / / /// //// // / / / / / / / / // Fig. 4. Half space and loading description.
The time intervals adopted for this analysis are: At = 5ms and At z = 10ms, while in the works taken for comparison the time steps selected were: At = l m s , in Tseng and Robinson, 21 and A t = 3-65ms in Mansur and Brebbia.16 In Fig. 5(a, b and c), the spatial discretizations adopted in Refs 21 and 16, as well as the one adopted here, are defined. The numerical values of displacements in the x2 direction, at a single point underneath the loaded area, are shown in Fig. 6. As can be seen, in spite of using greater time steps than those adopted by the other mentioned authors, the results computed with the proposed formulation seem to exhibit very good accuracy. It is important to say that the maximum difference found between our results (At~ = 5ms and At z = 10ms) is less than 1%. The third example, taken from a work written by Karabalis and Beskos, 18 consists of computing the horizontal and vertical displacements and rotations at the center of a rigid square footing placed on the elastic half-space. Their results were also obtained by using a three dimensional time domain BEM formulation, already reported in this work, based on the Dirac's delta fundamental solution and employing constant time approximation. The numerical results to be compared are due to external rectangular impulses, applied in the corresponding motion directions. For simplicity, the load distributions along the time are given together with the results obtained. The foundation material is assumed to be elastic with linear behavior defined by the following parameters: Modulus of elasticity, E = 2-5898 10 9 lbf/ft 2 Material density p = 10-368 lbfs2/ft 4 Poisson ratio u = 0.33 Relaxed boundary conditions are assumed for all
generalized displacement components. The contact shear stresses, cq3 and c~23, are taken equal to zero when vertical displacements and rotations are analyzed, while cr33 = c~23 = 0 is assumed to compute horizontal displacements. Figure 7(a,b) gives the discretizations adopted to run the problem with the proposed formulation, using only one linear boundary element to discretize the contact footing area. For the first discretization a reasonable amount of the free surface around the footing is discretized together with the footing area, while for the second one, only this area is taken into consideration. The four dots shown inside the footing area, Fig. 7(a,b), represent the selected collocation points. Figure 7(c) gives the discretization originally adopted by Karabalis and Beskos, is with 64 constant boundary elements. As expected, due to the nature of this problem, the same results are obtained for both discretizations considered here (Figs 7(a) and 7(b)). Vertical and horizontal displacements, due to the corresponding impulses, are given in Figs 8 and 9, respectively, while rotations around axis x 2 are presented in Fig. 10. These results are illustrated together with the values taken for comparison. Although they approximately describe the same behavior, one can see the smoothness of the curve achieved with the new fundamental solution. Some significant differences may be noted between our solution for vertical displacements and the one obtained by Karabalis and Beskos.18 One possible reason that justifies this difference could be the poor discretization adopted here, which enforces constant reaction underneath the footing (the code used to run this problem is prepared to analyze building structure-soil interaction where footings are taken as rigid elements with linear traction distributions). Although the reason pointed out above could be acceptable, we perform a similar analysis for which is possible to obtain analytical results from the
Further improvements on three dimensional transient BEM elastodynamic analys& Stokes general solution. The same rectangular impulse was assumed to be applied at a single point o f the infinite space, then the displacements have been computed at several points in the vicinity under the load
point. All these computed values exhibit a behavior similar to the one verified for the footing case. One can see, in Fig. 11, values obtained along some period o f time for a point placed at 2 ft below the load.
IO' IL
i
II
i
!
I
,
~ '
I a ) TSENG
[21.]
t
1,80'
z| b)
MANSUR
[16]
x1 r
:
"
;
;
:
;
:
"
:
:
;
',
:
237
'
'
'
'
'
:
'
;
:
SZO,
LOAOEO
REGION
c)THIS
WORK
110 °
x1
7~9' Fig. 5. Half space discretizations.
1
238
H. B. Coda, W. S. Venturini -6
I
k
oO.? O
/
X1
-4
•
F
7 LJJ
0
LIJ J
0
o#
-2
in
r
020
..j
•
&t~'
(3
At2 Tleng
o
Manl~r
I
40
60
80
I00
ms
ZO
t
Fig. 6. Vertical displacements at point D(0[ 0 I, -70/),
ColPoints ocotion~
,
",4 \ "-..., "'--..., ,
(o)
\
,
i
5
i
5
i
5
mmmmmmmm
mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmam mmmmmmmm
mmmmmmmm (b)
Fig. 7. Footing and half space discretizations.
(c)
Further improvements on three dimensional transient BEM elastodynamic analysis 0.5
lo-~
1
ft
I
1
~3(t) I
P3 (t)
0.4
~
l I I
I .~,
100 --
. 0.3
/_
239
~"~'x t
___~_~
~
0123
~
---
[1.el
I
0.2
~\\\x J
'
--~ ....
/ ""
i
~x
..~
\
0.I
~
0
0
5
~0
20
15
~ " -
Z5
30
~ "~" ~" ~
35 -4 40 x 0.18x10 s
Fig. 8. Footing vertical displacements.
Although Karabalis and Beskos 18 have proposed obtaining the response of that foundation, for any time varying force function from the impulse responses given in Figs 8-10, applying a convolutional approach,
in our work we compute these values employing directly the presented BEM formulation, to illustrate the generality of this approach. As can be seen, in Fig. 12, applying to this footing another load function, the
0,75,
[
ft 0.60
!
I
PX (t), ZOO--
i •
t
kips
~l
0.45
III
-...s
;~,~
I
1
I
t
-
01,23 x 0.18 x lO-#s
[18] This Work
I
,,
0.30
\
I
I
L
I
0,15
0
0
5
10
15
20
25
Fig. 9. Footing horizontal displacements.
30 38 40 x 0.18 x 10-4s
H. B. Coda, W. S. Venturini
240
;':', t,,
I.!
I
0.15
1
0.10
i
\i i :-i :`.-or.,
-
0.18x10
s
-0.0~ 0
S
~
lS
!
80
N
3O
Q m xO.18x10 -4 s
Fig. 10. Footing rotations.
o
8 -
llO'° tt )
uu~upuTuuluufuuu 0
/
ip|
IVllllrllllllllrll|
1
2
3
| 1~O-4s )
Fig. I 1. Vertical displacements due to a rectangular impulse in an infinite domain.
Further improvements on three dimensional transient B E M elastodynamic analysis
241
0.6 10-5
I
;o.
I/i:_
rk
0
O' -0.4
V
0
~0
40
IO
xO.18x10 -4 s
W
Fig. 12. Vertical displacements for a load given by P(t) = 180sin(13 000t)kips.
computed vertical displacements agree very well with the values given in Refs 5 and 13. The convolutions of both solutions, presented in Fig. 8, lead to very close results when the frequency of the load is not very high, as for this particular case. In the fourth example, the analysis of a short beam submitted to an impact of tangential forces is shown. The following material properties are taken to run this example: E = 52 x 103N/m 2, u = 0.30, p = 1.8 × 103kg/m 3. The assumed time interval is At = 0.72 x 10-2/s. The distributed tangential force has constant
value p -= 1000 N/m 2 and is assumed to be proportional to the Heaviside function in time. The geometry and discretization are shown in Fig. 13. The computed transversal displacement of the central point at the free end, compared with the static value of the shear deflection, is presented in Fig. 14. In Fig. 15, the longitudinal displacement due to a similar horizontal impact is plotted against the transversal one. The frequency behavior of these movements confirms the relation between the longitudinal (Cl = 180 m/s) and transversal (Ct = 100 m/s) wave velocities.
h:2m
:-'! •"
iI o°i
.~.: . . . . .
o
~-4m
t Fig. 13. Geometry and discretization.
i : T •
:,: .......
! Fr.. ! (..L.'I u ~
oor :~ ......
242
H. B. Coda, W. S. Venturini I
E
--=
n
Transversal displacement
numerical stability, when employing constant time approximation for displacements, are no longer required.
4.
J stotic deft. E 0
2-
£ r~
t (0.0072s)
Fig. 14. Transversal displacement.
One can observe that the displacement curves obtained when running this example are rather smooth. This is due to the poor discretization employed. It can be pointed out that the convergence of this formulation is reached when the approximate relation (Cl × At > 0.9A/) is followed, Al being the smaller element length. In this example relaxed conditions have been assumed to compute shear displacements only.
CONCLUSIONS A simple modification in the classical time dependent fundamental solution for three-dimensional elastic analysis was proposed. This new solution, when adopted together with the BEM formulation, reduces the number of coefficients to be stored and also improves the numerical results. Moreover, more flexibility is given to the method regarding space and time discretizations. Differentiations based on numerical differences, usually adopted by the classical formulation to achieve
---------
5-
Longitudinal d i s p l a c e m e n t Transversal displacement • ~°
E o
4-
,j
0 c Q E e
•
g 2.
,(
\
.!
Q
• r =K
o- O
0-~
t (0 0072s)
Fig. 15. Transversal versus longitudinal displacements.
REFERENCES
1. Friedman, M. B. & Shaw, R. P. Diffraction of pulses by cylindrical obstacles of arbitrary cross-section. J. Appl. Mech., 1962, 29, 40-6. 2. Niwa, Y., Fukui, T., Kato, S. & Fujiki, K. An application of the integral equation method to two-dimensional elastodynamics. Theor. Appl. Mech., 1980, 28, 281 90. 3. Mansur, W. J. A time-stepping technique to solve wave propagation problems using the boundary element method. Ph D thesis, Southampton University, 1983. 4. Mansur, W. J. & Brebbia, C. A. Transient elastodynamics using a time-stepping technique. In (ed.) C. A. Brebbia, T. Futagami & M. Tanaka, Proc. 5th. Int. Sem. BEM Eng. Springer-Verlag, Berlin, 1983, pp. 677-98. 5. Karabalis, D. L. & Beskos, D. E. Dynamics response of 3-D rigid surface foundations by time domain boundary element method. Earthqu. Engng Struct. Dynam., 1984, 12, 73-93. 6. Beskos, D. E. Boundary element methods in dynamic analysis. Appl. Mech. Rev., 1988, 40, 1-23. 7. Mansur, W. J. Boundary element method applications in two-dimensional transient elastodynamics. In Boundao' Elements X, Vol. 4, ed. C. A. Brebbia. CMP, Southampton, Springer-Verlag, Berlin, 1988, pp. 387 99. 8. Navier, ?. M~m. Acad. Sciences, Paris, 1827, 7; Love, A. E. H. A Treatise on the Mathematical Theory of Elasticity, 4th ed. Dover, New York, 1944. 9. Cauchy, A. L. Exercices de math6matique, sur l'bquilibre et le mouvement d'un systemb de points mat6riels par des forces d'attraction ou de r+pulsion mutuelle, 1828; Love, A. E. H. A Treatise on the Mathematical Theory of Elasticity, Vol. 4. Dover, New York, 1944. 10. Stokes, G. G. On the dynamical theory of diffraction. Trans. Camb. Phil. Soc., 1849, 9, 1. 11. Manolis, G. D. et al. Boundary element method implementation form three-dimensional transient elastodynamics. In P. K. Banerjee & J. O. Watson, ed. Developments in Boundary Element Methods-4, Elsevier, London, 1986, pp. 29-65. 12. Eringen, A. C. & Suhubi, E. S. Elastodynamics: Linear Theoo,, Vol. 2. Academic Press, New York/London, 1974. 13. Dominguez, J. & Aascal, R. Dynamics of foundations. In ed. C. A. Brebbia. Topics in Boundary Element Research. Applications in Geomechanics, Vol. 4. Springer-Verlag, Berlin, 1987, pp. 27-75. 14. Wheeler, L. T. & Sternberg, E. Some theorems in classical elastodynamics. Arch. Rational Mech. Anal., 1968, 31(1), 87 90. 15. Kobayashi, S. Elastodynamics. In Boundary Element Methods in Mechanics, Vol. 3, ed. D. E. Beskos. North Holland, Amsterdam, 1987, pp. 191-255. 16. Mansur, W. J. & Brebbia, C. A. Transient elastodynamics. In Topics in Boundary Element Research-2, ed. C. A. Brebbia. Springer-Verlag, Berlin, 1985, pp. 124-55. 17. Hartmann, F. Computing the C-matrix on nonsmooth boundary points. In ed. C. A. Brebbia, New Developments in Boundary Element Methods. CMP, Southampton, 1980. 18. Karabalis, D. L. & Beskos, D. E. Three-dimensional soilstructure interaction by boundary element methods. In Topics in Boundary Element Research: Applications in
Further improvements on three dimensional transient B E M elastodynamic analysis Geomechanics, Vol. 4, ed. C. A. Brebbia. Springer-Verlag, Berlin, 1987, pp. 1-26. 19. Coda, H. B. Three-dimensional transient analysis of structures by a BEM/FEM combination. Ph D thesis, University of S~o Paulo, S~o Carlos (Br), 1993 (in Portuguese). 20. Banerjee, P. K., Ahmad, S. & Manolis, G. D. Advanced
243
elastodynamic analysis. In Boundary Element Methods in Mechanics, Vol. 3, ed. D. E. Beskos. North Holland, Amsterdam, 1987, pp 257-284. 21. Tseng, M. N. & Robinson, A. R. A transmitting boundary for finite-difference analysis of wave propagation in solids. Project n. NR 064-183, University of Illinois, Urbana, 1975.
PII:
ELSEVIER
S0955-7997(96)00019-7
Copyright © 1996. Published by Elsevier Science Ltd Printed in Great Britain. All rights reserved 0955-7997/96/$15.00
Further improvements on three dimensional transient BEM elastodynamic analysis Humberto Breves Coda & Wilson Sergio Venturini* S~o Carlos School of Engineering, University of S~o Paulo, Av. Carlos Botelho 1465, 1356-970 S~o Carlos, Brazil
(Received 14 January 1994; revised version received 6 March 1996; accepted 20 May 1996) In this work the transient boundary element formulation is studied together with a modified time dependent fundamental solution. This modification permits one to obtain more stable results during numerical analysis. It also reduces the amount of matrices usually required to be stored. Furthermore, some tricks, usually needed when the classical fundamental solution is adopted, are now avoided, improving the result quality and eliminating some discretization restrictions. Copyright © 1996. Published by Elsevier Science Ltd Key words: Boundary element method, elastodynamics, time domain, transient
problem. INTRODUCTION
applied at a single instant. In this work, another assumption regarding the load choice is made: a unit impulse distributed along an interval of time was taken to derive the necessary fundamental values. The numerical algorithms developed from the standard fundamental solution described above, which has been used so far, required the use of some inconvenient tricks to deal with terms containing Dirac's delta derivatives. For instance, we can point out the introduction of techniques, usually adopted to compute velocity terms (backward finite differences), to allow the use of constant time approximation for displacements. 11 In order to improve the procedure, many authors have adopted linear approximation for displacements and constant approximation for tractions. These necessary tricks also increase the number of parameters to be computed for each time integration step. The alternative fundamental solution proposed here has been demonstrated to produce accurate results when analyzing practical examples. The use of backward finite differences has been avoided. The computing time is also reduced due to the smaller amount of parameters to be evaluated in comparison with the classical approach.
The first time domain formulation of the boundary element method (BEM) was proposed by Friedman and Shaw, 1 in 1962, in which the acoustic anti-plane problem was treated. In the context of two dimensional elastodynamic analysis, the first work was presented by Niwa et al., 2 in 1980; nevertheless it was only in 1983 that a well established time stepping technique was presented by Mansur and Brebbia. 3,4 Karabalis and Beskos, 5 in 1984, were the first to develop a general three-dimensional elastodynamic methodology in time domain BEM. Interesting overviews on BEM formulations to treat elastic transient problems can be seen in Beskos 6 and Mansur. 7 The BEM formulations, employed to study any physical problem, are derived from some well established integral equations. These equations are usually achieved by using fundamental solutions, which are obtained by solving the corresponding differential equation of the problem, for an infinite domain, where a convenient unit load is applied. For the case of elastic bodies, the fundamental solution comes from the Navier-Cauchy equilibrium equations. 8'9 For elastodynamic problems, the general singular solution is due to Stokes, 1° representing the behavior of an infinite body subjected to a time dependent concentrated load. The classical fundamental solution, used in the time domain BEM approach, is a particular case of the Stokes' state, assuming that the load is a unit concentrated impulse
FUNDAMENTAL SOLUTION Assuming that the bodies under analysis are isotropic and show linear elastic behavior, the governing differential equations for the dynamic equilibrium, the so-called
*To whom correspondence should be addressed. 231
H. B. Coda, W. S. Venturini
232 Navier-Cauchy equations, are given by
(C~ -- C2)u],ji + C2bli,j] + hi~ p : Hi
(1
Bik(q,v-;s,O)=
3C 2 5 r5
r3
j
where b i represents body forces, p is the media density, and Cl and C2 are longitudinal and shear wave propagation velocities, respectively. As usual, we can represent a time-dependent concentrated load in the k-direction, acting at any point's' of the infinite domain, by the following expression:
b*ki =.f(r)5(q -- s)6kg
+216F
(2)
rirjr k7
where q is a field point, ¢5ki is the Kronecker delta, 5 ( s - q ) the Dirac delta and f ( r ) a time-dependent load function. Replacing the body force term in eqn (1) according to eqn (2) gives, (C~
2 • 2 . . . - C2)uj,ji + C2ug,jj +f(r)~5(q -. S)6ki/P = ui
(3)
in which the symbol * indicates fundamental values. By integrating eqn (3) one obtains the so-called Stokes' fundamental solution,m easily found in the literature.12'13 The usual fundamental solutions, employed for all previous works mentioned here, are the particular case of the Stokes' state in which f ( r ) = 5(r) is assumed. This load represents a unit impulse at time r = 0. The u*kg and P*kg expressions, for this particular case, are found in several well known works. 12,14 Trying to obtain an alternative which could improve the quality of the numerical responses for BEM transient analysis, we decided to use another load function together with the Stokes' state, given by the following expression,
b*~i = [H(r) - H(r - At)]6(q - s)6ki/At
(4)
in which H( ) is the well known Heaviside function and At a constant time interval. This load no longer represents a unit impulse at an instant, as in the classical approach, but a unit impulse developed over a time interval At. Displacement and traction values, derived from this particular time dependent load, are given by,
u*ki(q, "r;s, O) = [A*ki(q, r ;s, O) - A*ki(q, r - At; s, 0)] (5) P*ki(q, 7-; s, O) = [B*ki(q,r; s, O) - B~i(q , 7- - At; s, 0)] (6) where Akg and B~g are:
A~i(q,r;s,o)--
l { (3rgrk 4rrp/Xt \ r 3
~5~g) 1 2r 2
rirjrk [ ~
+ 2r-CdTC
r)
C~5(r
+~6
7-r
r)]
\ -U,
rk~50[ 2 fC2"~1 [ H ( r r3 1 - ~C~] J
~-2) -
r3 r (5
r
Equations (7) and (8), when compared with the corresponding classical ones, show a smoother time behavior, i.e. the new solution does not exhibit the undesirable Dirac's delta derivatives. As a consequence, no velocity terms will appear when time convolutions are performed. As for the classical fundamental solution, the new one goes to zero when r>>.r/C2+At, which makes it possible to compute only a finite amount of coefficients.
INTEGRAL REPRESENTATIONS
The Graffi reciprocity theorem, for dynamic problems, can be written for any two elastodynamic states, here represented by their displacement fields u and u', as follows 14-t6
/'J [ui(q,r)b~(q,t
" H r 2 _ rC2
tSijrk + tSikr j r3 + 6jkri]
!
0
~
' t - r)bi(q, r ) ] d a d r - r) - ui(q,
f f [ui(q, r ) b : ( q , t - r ) - u i ( q , t - r)bi(q, r ) ] d ~ d r 3o J ~2 rkr"[~{H(r--~l + r3
) -- I~H C2
= I l l [ui(Q,t-r)pi(Q,r)-ui(Q, ' 3o3 I'
, - r)]dFdr r)pi(Q,t
+ p [ [tJ'i(q, t)ui(q, O) + u'i(q, t)ui(q, 0)]df~ •J~2
(9)
Further improvements on three dimensional transient B E M elastodynamic analysis where 9t is a bounded or unbounded region and P its contour. Without lost of generality, only the case of quiescent past elastodynamic states will be studied in this work. Thus, eqn (9) becomes, [t [ [ui(q, "r)bti(q, t _ T ) - u l ( q , t -
"r)bi(q, T ) ] d ~ d T
JO d f~
233
and P*ki and ~:i are the Dirac's delta fundamental solution. When constant time approximation is adopted, for both displacement and traction, the explicit terms P~,i and U~i become similar to those derived for the new solution. Nevertheless, one can easily realize that they are applied in different equations, (13) and (11) respectively, therefore they have different meanings.
= I'o Jr[u,(Q, t-r)pi(Q, r)-ui(Q, T)pi(Q, ' t-T)]dPdT (10) Assuming the new fundamental solution as one of those elastodynamic states, eqn (10) gives,
ki[~,S)Jt_At- u,ls ~
l d r = Jol* rUki(Q,t;s,T)pi(Q,T)dFd7 +
Jol ui(Q,r)pki(Q, *
t; s,r)drd~-
NUMERICAL ASPECTS As usual for boundary element formulations, eqn (11) can be treated numerically be adopting appropriate space and time discretizations. For the particular case of the time discretization, constant elements have been adopted. Thus, eqn (11) can be written as follows,
Cki(s)UiN,(s,t) = [
[to U*ki(Q,t;s,7)¢J(Q)~°dTdFP~
- [
[ P*k~(Q,t;s,r)¢~(Q)¢°drdrU~d
Jr(/) Jto 1
Uo
U*ki(q,t; s, T)bi(q, T)dfMT
to (ll)
where Cki is the classical free term of displacement integral representations computed as in Hartmann. 17 In eqn (11), one has to use U*kiand p~i, as defined in (5) and (6), taking into account the following property of the Stokes state, U*ki(Q, t - r ; s, O) = U*ki(Q, t; s, r )
(12a)
P*ki(Q, t - r ; s, O) = P*ki(Q, t; s, y)
(12b)
It is worth mentioning, that, when the classical fundamental solution (Dirac's delta load function) is used, as in Karabalis and Beskos, 5'18 eqn (10) reduces to Love's identity and the time convolution leads to functional operators belonging to the Stokes' state, l°J2J4 This expression can be written as:
- Ir P*ki(Q, t; s I ui)dFdT
where e*ki and U~¢i
(13)
given by,
are
j u*ki(q,t, s, r)pi(q, T)dT * Pki(Q,t;s I ui) = 0 fi*ki(q, t, s, T)ui(q, I' t
U~i(Q, t; s I pn) =
o
7")d~-
(14)
Jto_1
(15) where 0 ranges from 1 to the integration step number, Nt, and the summation rule is implied. For the new approach, the value of Uis,(S, t), in eqn (15), means an average value computed along the time step Nt. It is not a single value taken at an instant t = NtAt, as for the classical procedure. The space approximations of displacements and tractions are given by,
Uio(Q j, t) = ¢~)¢(0)V;~oj
(16a)
Pio(Q j, t) = A(JL/,(°) j ~c~ ~ P?iO
(16b)
in which i j 0 ¢, ~0
Cki(Q, S)Ui(S~ T) = [ U~i(Q, t; s I pi)drdT JF
+J;~ P*ki(q,t;slbi)dfMT
JF(/)
= = = = =
displacement or traction directions space elements time step shape function for the space element time shape function (¢e = 1)
Flat quadrangular elements with linear and quadratic approximations have been adopted for the space discretization. The singular integrals have been computed with a special algorithm 19 based on Kutt's quadrature, For non-singular integrals, the standard Gauss quadrature was adopted together with a well known element subdivision technique) l In order to integrate the time dependent kernels in eqn (15), the time shape function ¢ = 1 was adopted,
234
H. B. Coda, W. S. Venturini
giving,
"E
J k5 to_~
r21/,- C,')
-
r
- - -
dr
if t - to >>- r
J
-- ~[(12 r~_~(TA_IoI)_I(TA~[o.I)~_(TA3[;_I)/31
(17)
if t - to_ 1 ~ -~, ~ t -- to
0
,o , A t
C~,]
if ~r / > t -- to-I
i
dr=
t-r-
dXt
) if,
0
Jii , 1 6 (
t-r-
r~'~dr
={~T
C~J
0
if t - t o _ j >
i f ~ >~t - to
ift-t° if~>t-to
and 7 is equal to 1 or 2. Subtracting from the values given by eqns (17)-(19), the corresponding ones computed at the next time step, [to, to+l], the appropriate time coefficients of the new fundamental solution, eqns (5) and (6), are obtained. As for any boundary element formulation, the displacement integral representation is transformed into an algebraic equation, after carrying out all space and time integrals over the elements defined by both discretizations. Selecting an appropriate number of collocation points 's' to write their algebraic representations, a system of linear equations is achieved as follows: (21 )
where 0 = 1,2, 3 . . . . , Nt.
(19)
i Ol'?7
(20)
H°U ° = G ° P °
=
,-
where TA = t - r/CA
(18)
)t-to
Diagonal sub-matrices of h/~) are obtained by the following expression, h~) = h Nt ,s + C
(22)
In eqn (22) I]sNt is obtained by carrying out singular integrals only at the later time step, i.e. exactly where we are computing the unknowns. Analogously for the classical transient BEM techniques employed by other authors] ~'ls the maximum time step recommended to compute new matrices is given by N, -
dmax C2At
+ 1
(23)
where dmax is defined as the maximum possible length of the discretized solid. We have to point out that when the classical solution is used together with constant/constant (displacement/traction) time approximations, the achieved
Om
(,
f E = 1,1 x 10 5N/m 2 p = 2 kg/m 2 v=0
t
Fig. 1. Prismatic beam: geometry and loading.
235
Further improvements on three dimensional transient B E M elastodynamic analysis
/
.
employed to verify numerical models, consists of analyzing a prismatic rod subjected to an impact at the free end. The results obtained, by using the proposed formulation, are compared with the classical numerical solution obtained by Banerjee et al., 20 using Dirac's delta fundamental solution and constant time approximation. The geometry and the material properties, as well as the load time behavior, are given in Fig. 1. A total of 18 quadratic elements, as illustrated in Fig. 2, were used to discretize the boundary of this solid. The relaxed boundary conditions, 0-31= 0 and 0-32= 0 , were also assumed. The analysis has been performed for two different time step intervals: 0.00356s, as used in Ref. 20, and At = 0.00712 s. The results obtained, in terms of stresses at the center of this beam, are shown in Fig. 3. As can be seen (Fig. 3), for both selected time steps, the numerical values obtained follow the analytical solution very well. Even for a rather large time step, good quality numerical results are achieved. In this case, the solution does not show abrupt variation of stresses due to the fundamental solution characteristics, which gives averaged values for each time step. For the smaller time step, the sharp jumps in stresses, for t I = 0.017s and t2 = 0.01 s, are very well captured. In the second example, a three dimensional halfspace, when subjected to a strip load, is analyzed. The half-space and the load intensity, as well as its time behavior, are shown in Fig. 4. The results obtained by adopting the proposed formulation are compared with those taken from Mansur and Brebbia, 16 computed by solving the corresponding two dimensional problem, also using a time domain BEM formulation. The finite difference solution, taken from Tseng and Robinson, El is also used for this comparison. The following values, related to material properties, have been adopted to solve the problem: E = 2 0 0 k s i , v--0-15 and p = 6.885 gbfsE/in4.
/
Fig. 2. Prismatic beam: discretization.
system of equations, equivalent to (21) is represented by, (H l + HHI)U 1 - - G I p 1
(0-- 1)
H°U ° + [HH°U° - n n 0 - 1 U 0 - 1
G°p °
] =
(0 = 2, 3,... N,)
(24)
In the same way, if linear/constant time approximation were adopted, one obtains, o 0+ HIU1
o 0
H2U2
Gop0
:
(25)
In eqn (24), HH ° contains the coefficients related with the backward difference scheme, as shown by Manolis et al. I1 On the other hand, eqn (25) exhibits two matrices, H°l and H °, achieved when linear approximation is used for displacements. 15 Thus, the new approach avoids dealing with finite differences and reduces the global amount of computing.
EXAMPLES Four examples have been selected to show the performance of this new fundamental solution, when used appropriately together with the BEM formulation. The first example, a classical problem usually
1
ANALYTICAl.
o BANERJEE + THIS WORK • THIS WORK
P ( N/m z ) .,...4,,222,,,.m
2000-
| At-O.OO356s) ( At, OOO366s| ( At,O.OOT12 s)
°
•
0 0
I000-
t
~
w
"
"
w ~r 0
-
4" v 0
t l 4. ÷
!
0
. Z
I
T 4
,
, 6
I
l II
! I;21 10
l I l , l , 14 IQ 111
I , ZO
l I l Im |4
I
I l , I ~ l 26 all
Fig. 3. Normal stress at beam middle point.
,
32
•
,
N
t ( 3.56 K)'ss )
236
H. B. Coda, W. S. Venturini
rt"
lksi
,
~,.0 m i
,,
/
///////
TIME
,/'/'//
/ /
11
/ / / / / / / / / / / / / / / / / / /// //// // / / / / / / / / // Fig. 4. Half space and loading description.
The time intervals adopted for this analysis are: At = 5ms and At z = 10ms, while in the works taken for comparison the time steps selected were: At = l m s , in Tseng and Robinson, 21 and A t = 3-65ms in Mansur and Brebbia.16 In Fig. 5(a, b and c), the spatial discretizations adopted in Refs 21 and 16, as well as the one adopted here, are defined. The numerical values of displacements in the x2 direction, at a single point underneath the loaded area, are shown in Fig. 6. As can be seen, in spite of using greater time steps than those adopted by the other mentioned authors, the results computed with the proposed formulation seem to exhibit very good accuracy. It is important to say that the maximum difference found between our results (At~ = 5ms and At z = 10ms) is less than 1%. The third example, taken from a work written by Karabalis and Beskos, 18 consists of computing the horizontal and vertical displacements and rotations at the center of a rigid square footing placed on the elastic half-space. Their results were also obtained by using a three dimensional time domain BEM formulation, already reported in this work, based on the Dirac's delta fundamental solution and employing constant time approximation. The numerical results to be compared are due to external rectangular impulses, applied in the corresponding motion directions. For simplicity, the load distributions along the time are given together with the results obtained. The foundation material is assumed to be elastic with linear behavior defined by the following parameters: Modulus of elasticity, E = 2-5898 10 9 lbf/ft 2 Material density p = 10-368 lbfs2/ft 4 Poisson ratio u = 0.33 Relaxed boundary conditions are assumed for all
generalized displacement components. The contact shear stresses, cq3 and c~23, are taken equal to zero when vertical displacements and rotations are analyzed, while cr33 = c~23 = 0 is assumed to compute horizontal displacements. Figure 7(a,b) gives the discretizations adopted to run the problem with the proposed formulation, using only one linear boundary element to discretize the contact footing area. For the first discretization a reasonable amount of the free surface around the footing is discretized together with the footing area, while for the second one, only this area is taken into consideration. The four dots shown inside the footing area, Fig. 7(a,b), represent the selected collocation points. Figure 7(c) gives the discretization originally adopted by Karabalis and Beskos, is with 64 constant boundary elements. As expected, due to the nature of this problem, the same results are obtained for both discretizations considered here (Figs 7(a) and 7(b)). Vertical and horizontal displacements, due to the corresponding impulses, are given in Figs 8 and 9, respectively, while rotations around axis x 2 are presented in Fig. 10. These results are illustrated together with the values taken for comparison. Although they approximately describe the same behavior, one can see the smoothness of the curve achieved with the new fundamental solution. Some significant differences may be noted between our solution for vertical displacements and the one obtained by Karabalis and Beskos.18 One possible reason that justifies this difference could be the poor discretization adopted here, which enforces constant reaction underneath the footing (the code used to run this problem is prepared to analyze building structure-soil interaction where footings are taken as rigid elements with linear traction distributions). Although the reason pointed out above could be acceptable, we perform a similar analysis for which is possible to obtain analytical results from the
Further improvements on three dimensional transient BEM elastodynamic analys& Stokes general solution. The same rectangular impulse was assumed to be applied at a single point o f the infinite space, then the displacements have been computed at several points in the vicinity under the load
point. All these computed values exhibit a behavior similar to the one verified for the footing case. One can see, in Fig. 11, values obtained along some period o f time for a point placed at 2 ft below the load.
IO' IL
i
II
i
!
I
,
~ '
I a ) TSENG
[21.]
t
1,80'
z| b)
MANSUR
[16]
x1 r
:
"
;
;
:
;
:
"
:
:
;
',
:
237
'
'
'
'
'
:
'
;
:
SZO,
LOAOEO
REGION
c)THIS
WORK
110 °
x1
7~9' Fig. 5. Half space discretizations.
1
238
H. B. Coda, W. S. Venturini -6
I
k
oO.? O
/
X1
-4
•
F
7 LJJ
0
LIJ J
0
o#
-2
in
r
020
..j
•
&t~'
(3
At2 Tleng
o
Manl~r
I
40
60
80
I00
ms
ZO
t
Fig. 6. Vertical displacements at point D(0[ 0 I, -70/),
ColPoints ocotion~
,
",4 \ "-..., "'--..., ,
(o)
\
,
i
5
i
5
i
5
mmmmmmmm
mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmmm mmmmmmam mmmmmmmm
mmmmmmmm (b)
Fig. 7. Footing and half space discretizations.
(c)
Further improvements on three dimensional transient BEM elastodynamic analysis 0.5
lo-~
1
ft
I
1
~3(t) I
P3 (t)
0.4
~
l I I
I .~,
100 --
. 0.3
/_
239
~"~'x t
___~_~
~
0123
~
---
[1.el
I
0.2
~\\\x J
'
--~ ....
/ ""
i
~x
..~
\
0.I
~
0
0
5
~0
20
15
~ " -
Z5
30
~ "~" ~" ~
35 -4 40 x 0.18x10 s
Fig. 8. Footing vertical displacements.
Although Karabalis and Beskos 18 have proposed obtaining the response of that foundation, for any time varying force function from the impulse responses given in Figs 8-10, applying a convolutional approach,
in our work we compute these values employing directly the presented BEM formulation, to illustrate the generality of this approach. As can be seen, in Fig. 12, applying to this footing another load function, the
0,75,
[
ft 0.60
!
I
PX (t), ZOO--
i •
t
kips
~l
0.45
III
-...s
;~,~
I
1
I
t
-
01,23 x 0.18 x lO-#s
[18] This Work
I
,,
0.30
\
I
I
L
I
0,15
0
0
5
10
15
20
25
Fig. 9. Footing horizontal displacements.
30 38 40 x 0.18 x 10-4s
H. B. Coda, W. S. Venturini
240
;':', t,,
I.!
I
0.15
1
0.10
i
\i i :-i :`.-or.,
-
0.18x10
s
-0.0~ 0
S
~
lS
!
80
N
3O
Q m xO.18x10 -4 s
Fig. 10. Footing rotations.
o
8 -
llO'° tt )
uu~upuTuuluufuuu 0
/
ip|
IVllllrllllllllrll|
1
2
3
| 1~O-4s )
Fig. I 1. Vertical displacements due to a rectangular impulse in an infinite domain.
Further improvements on three dimensional transient B E M elastodynamic analysis
241
0.6 10-5
I
;o.
I/i:_
rk
0
O' -0.4
V
0
~0
40
IO
xO.18x10 -4 s
W
Fig. 12. Vertical displacements for a load given by P(t) = 180sin(13 000t)kips.
computed vertical displacements agree very well with the values given in Refs 5 and 13. The convolutions of both solutions, presented in Fig. 8, lead to very close results when the frequency of the load is not very high, as for this particular case. In the fourth example, the analysis of a short beam submitted to an impact of tangential forces is shown. The following material properties are taken to run this example: E = 52 x 103N/m 2, u = 0.30, p = 1.8 × 103kg/m 3. The assumed time interval is At = 0.72 x 10-2/s. The distributed tangential force has constant
value p -= 1000 N/m 2 and is assumed to be proportional to the Heaviside function in time. The geometry and discretization are shown in Fig. 13. The computed transversal displacement of the central point at the free end, compared with the static value of the shear deflection, is presented in Fig. 14. In Fig. 15, the longitudinal displacement due to a similar horizontal impact is plotted against the transversal one. The frequency behavior of these movements confirms the relation between the longitudinal (Cl = 180 m/s) and transversal (Ct = 100 m/s) wave velocities.
h:2m
:-'! •"
iI o°i
.~.: . . . . .
o
~-4m
t Fig. 13. Geometry and discretization.
i : T •
:,: .......
! Fr.. ! (..L.'I u ~
oor :~ ......
242
H. B. Coda, W. S. Venturini I
E
--=
n
Transversal displacement
numerical stability, when employing constant time approximation for displacements, are no longer required.
4.
J stotic deft. E 0
2-
£ r~
t (0.0072s)
Fig. 14. Transversal displacement.
One can observe that the displacement curves obtained when running this example are rather smooth. This is due to the poor discretization employed. It can be pointed out that the convergence of this formulation is reached when the approximate relation (Cl × At > 0.9A/) is followed, Al being the smaller element length. In this example relaxed conditions have been assumed to compute shear displacements only.
CONCLUSIONS A simple modification in the classical time dependent fundamental solution for three-dimensional elastic analysis was proposed. This new solution, when adopted together with the BEM formulation, reduces the number of coefficients to be stored and also improves the numerical results. Moreover, more flexibility is given to the method regarding space and time discretizations. Differentiations based on numerical differences, usually adopted by the classical formulation to achieve
---------
5-
Longitudinal d i s p l a c e m e n t Transversal displacement • ~°
E o
4-
,j
0 c Q E e
•
g 2.
,(
\
.!
Q
• r =K
o- O
0-~
t (0 0072s)
Fig. 15. Transversal versus longitudinal displacements.
REFERENCES
1. Friedman, M. B. & Shaw, R. P. Diffraction of pulses by cylindrical obstacles of arbitrary cross-section. J. Appl. Mech., 1962, 29, 40-6. 2. Niwa, Y., Fukui, T., Kato, S. & Fujiki, K. An application of the integral equation method to two-dimensional elastodynamics. Theor. Appl. Mech., 1980, 28, 281 90. 3. Mansur, W. J. A time-stepping technique to solve wave propagation problems using the boundary element method. Ph D thesis, Southampton University, 1983. 4. Mansur, W. J. & Brebbia, C. A. Transient elastodynamics using a time-stepping technique. In (ed.) C. A. Brebbia, T. Futagami & M. Tanaka, Proc. 5th. Int. Sem. BEM Eng. Springer-Verlag, Berlin, 1983, pp. 677-98. 5. Karabalis, D. L. & Beskos, D. E. Dynamics response of 3-D rigid surface foundations by time domain boundary element method. Earthqu. Engng Struct. Dynam., 1984, 12, 73-93. 6. Beskos, D. E. Boundary element methods in dynamic analysis. Appl. Mech. Rev., 1988, 40, 1-23. 7. Mansur, W. J. Boundary element method applications in two-dimensional transient elastodynamics. In Boundao' Elements X, Vol. 4, ed. C. A. Brebbia. CMP, Southampton, Springer-Verlag, Berlin, 1988, pp. 387 99. 8. Navier, ?. M~m. Acad. Sciences, Paris, 1827, 7; Love, A. E. H. A Treatise on the Mathematical Theory of Elasticity, 4th ed. Dover, New York, 1944. 9. Cauchy, A. L. Exercices de math6matique, sur l'bquilibre et le mouvement d'un systemb de points mat6riels par des forces d'attraction ou de r+pulsion mutuelle, 1828; Love, A. E. H. A Treatise on the Mathematical Theory of Elasticity, Vol. 4. Dover, New York, 1944. 10. Stokes, G. G. On the dynamical theory of diffraction. Trans. Camb. Phil. Soc., 1849, 9, 1. 11. Manolis, G. D. et al. Boundary element method implementation form three-dimensional transient elastodynamics. In P. K. Banerjee & J. O. Watson, ed. Developments in Boundary Element Methods-4, Elsevier, London, 1986, pp. 29-65. 12. Eringen, A. C. & Suhubi, E. S. Elastodynamics: Linear Theoo,, Vol. 2. Academic Press, New York/London, 1974. 13. Dominguez, J. & Aascal, R. Dynamics of foundations. In ed. C. A. Brebbia. Topics in Boundary Element Research. Applications in Geomechanics, Vol. 4. Springer-Verlag, Berlin, 1987, pp. 27-75. 14. Wheeler, L. T. & Sternberg, E. Some theorems in classical elastodynamics. Arch. Rational Mech. Anal., 1968, 31(1), 87 90. 15. Kobayashi, S. Elastodynamics. In Boundary Element Methods in Mechanics, Vol. 3, ed. D. E. Beskos. North Holland, Amsterdam, 1987, pp. 191-255. 16. Mansur, W. J. & Brebbia, C. A. Transient elastodynamics. In Topics in Boundary Element Research-2, ed. C. A. Brebbia. Springer-Verlag, Berlin, 1985, pp. 124-55. 17. Hartmann, F. Computing the C-matrix on nonsmooth boundary points. In ed. C. A. Brebbia, New Developments in Boundary Element Methods. CMP, Southampton, 1980. 18. Karabalis, D. L. & Beskos, D. E. Three-dimensional soilstructure interaction by boundary element methods. In Topics in Boundary Element Research: Applications in
Further improvements on three dimensional transient B E M elastodynamic analysis Geomechanics, Vol. 4, ed. C. A. Brebbia. Springer-Verlag, Berlin, 1987, pp. 1-26. 19. Coda, H. B. Three-dimensional transient analysis of structures by a BEM/FEM combination. Ph D thesis, University of S~o Paulo, S~o Carlos (Br), 1993 (in Portuguese). 20. Banerjee, P. K., Ahmad, S. & Manolis, G. D. Advanced
243
elastodynamic analysis. In Boundary Element Methods in Mechanics, Vol. 3, ed. D. E. Beskos. North Holland, Amsterdam, 1987, pp 257-284. 21. Tseng, M. N. & Robinson, A. R. A transmitting boundary for finite-difference analysis of wave propagation in solids. Project n. NR 064-183, University of Illinois, Urbana, 1975.