- Email: [email protected]
A simple comparison between two 3D time domain elastodynamic boundary element formulations
0955-7997(95)00088-7
ELSEVIER
Engineering Analysis with Boundary Elements 17 (1996) 33-44 Copyright © 1996. Published by Elsevier Science Limited Printed in Great Britain. All rights reserved 0955-7997/96/$15.00
A simple comparison between two 3D time domain elastodynamic boundary element formulations H. B. Coda & W. S. Venturini* Sao Carlos School of Engineering, University of S~o Paulo, Av. Carlos Botelho 1465, 13560-970 S~o Carlos, Brazil
(Received 9 March 1995; revised version received 9 September 1995; accepted 18 October 1995) The three-dimensional transient boundary element method for elastodynamic analysis is studied here to verify its performance concerning some required features when solving practical problems. The classical formulation was established more than 10 years ago and has already been applied several times, mainly in the context of infinite domain analysis. In this work the classical fundamental solution is compared with an alternative one in which a fundamental solution based on the unit impulse is assumed distributed over a time interval. Some interesting aspects of both solutions, regarding the result confidence for practical analysis, are also discussed in the paper. Key words: Boundary element method, time domain, elastodynamic analysis.
the consistence of the developed algorithms, the amount of memory required by them, the expected accuracy, and the stability ranges achieved for numerical applications.
INTRODUCTION The first application of the time domain boundary element method (BEM) apparently was made by Friedman and Shaw, l in 1962, in which the acoustic anti-plane problem was studied. In the context of two dimensional elastodynamic analysis, the first work was presented by Niwa et al., 2 in 1980; nevertheless, in 1983 only, a well established time stepping technique was presented by Mansur and Brebbia. 3 Karabalis and Beskos, 4 as well as Manolis, 5 in 1984, were the first to examine the general three-dimensional elastodynamic methodology in time domain BEM. In order to have a more complete overview of the elastodynamic BEM formulations the reader is invited to see the works written by Beskos 6 and Mansur. 7 Complete works have been recently published describing the BEM elastodynamic formulations, in details, and their practical applications, s'9 More recently, an alternative fundamental solution for BEM transient formulation was presented by Coda and Venturini. l°al This solution is also derived from the Stokes' state representing the responses due to a unit impulse accumulated along an interval of time, herein named, for convenience, the new fundamental solution. In this work a comparative analysis of the boundary element transient formulations using these two fundamental solutions is carried out. Important aspects of those formulations will be discussed: the coherence and
INTEGRAL REPRESENTATION OF DISPLACEMENTS FOR ELASTODYNAMICS Assuming that the bodies under analysis are isotropic and show linear elastic behaviors, the governing differential equations for the dynamic equilibrium, the so-called Navier-Cauchy equations, are given by: ( C 2 - C2)uj,ji + C~ui,jj + bi/p = ui
(1)
where bi represents the body forces, p is the media density, and C1 and C2 are the longitudinal and shear wave velocities, respectively. A family of possible fundamental solutions, Stokes' state, 4a2 is achieved by assuming, in eqn (1), the following body force term, b~i = f ( r ) a ( q - s)tki
(2)
in which s and q are the load and field points respectively, and f ( r ) gives the load time behaviour. Using this general solution together with the Graffi's theorem, 13a4 one can write the final displacement integral representation, in the same form several times shown in the literature. For a body of domain [2 and boundary F, with a quiescent past, submitted to a loading process, that began at a time r = 0 and is extended to a time t, the displacement integration representation
*To whom correspondence should be addressed. 33
34
H. B. Coda, W. S. Venturini
The new solution also exhibits a variable wave front width given by
is reduced as follows:
;o JoJ*
Cki(Q,s )
=
ui(s, 7 " ) f ( l -
7") dT-
Uki(Q,t-7";s/f)p~(Q,7")drdT"
- J; Jr Ui(Q, 7")p*ki(O, t - 7";s / f ) dr dr
+ IiI uzilq,,- 7";s/ytb,lq, ld dT
(3)
where u*ki and P*ki a r e fundamental values and Cki is an independent free term, similar to those used for elastostatic formulations.
FUNDAMENTAL SOLUTIONS
(4)
Adopting the body force distribution given in eqn (4), the fundamental displacement expressions are obtained, as shown in many works. 13-15 For completeness these expressions are displayed in Appendix 1. The new fundamental solution also represents the response to an infinite body subjected to a unit impulse, now distributed along an interval of time denoted by At, instead of applied at a specified instant. The expression adopted to represent this unit impulse is given by b*ki = [H(r) - H ( r - A t ) ] 6 ( q - S)6ki/At
where the subscript 'n' is used to indicate the new solution. The wave front width difference between both fundamental values, more relevant along the first time steps, is a reason for the improved behaviour observed in the numerical responses given by the BEM algorithm developed with the new solution adopted, in comparison with those computed by using the classical solution. It must be said that the unit impulse can be developed along any time interval which leads to a more general representation of eqn (5), b*ki = I n ( 7 ) -- U ( 7 -- M A t ) ] t S ( s
Both fundamental solutions analysed here are derived from the general Stokes' state, 12 assuming an appropriate time dependent load function f ( r ) to define the body force expression, eqn (2). The classical one represents a response of an infinite body subjected to a concentrated unit impulse applied at an instant r = 0. In this case eqn (2) becomes b*ki = 6(r)~(q - s)loki
(9)
In = (C, - C 2 ) N A t + C 2 A t
- q)~Ski/MAt
(10)
Using this expression, we are able to obtain the wave front widths needed for our purposes. The value M stands for the number of time intervals over which the unit impulse is distributed. It has been taken equal to one for the present analysis. The solution of some numerical examples can demonstrate that this assumption gives very precise results and requires a reasonable amount of computing time. Another reason is that it improves the numerical responses regarding the shape exhibited by the new solution with respect to time. The new fundamental values are represented by rather smoother curves in comparison with those given by the classical formula. This characteristic also explain why some time discretization restrictions, as described in the next section, must be enforced on the formulation based on the classical fundamental solution.
(5)
in which H ( r ) is the Heaviside function. For this case, the fundamental values, derived by Coda and Venturini, l°'1l can be represented as follows: U*ki(q, 7"; S, O) = [a~i(q , 7"; s, O) - A*ki(q, 7" -- A t ; s, 0)]
(6)
NUMERICAL ASPECTS First of all, we must rewrite expression (3) taking into account the types of load functions discussed earlier. Considering the load function adopted to derive the classical solution, eqn (4), we can achieve the following integral representation:
P*ki(q, 7-; s, O) = [B~i(q, r; s, O) - B~i(q , r - At; s, 0)]
(7) for which A*ki and BT,i are also given in Appendix 1. By analysing the classical fundamental solution we can verify that the corresponding wave front exhibits a variable width given by lc~ = (C1 - C 2 ) N A t
(8)
where the subscript 'Cl' is used to represent the classical solution and N is the number of time intervals over which the integral has already been performed.
Ck~(Q,s)ui(s,t) =
o
u*~(Q,t;s,r)pi(Q,r)drd7-
-
u~ (Q, r)p*k~(Q, t; s, r) dr d r
+Ji1
U*ki(q, t; S, r) bi (q, r) df~ dr
(11)
Following the same steps now using the new load function, eqn (5), the integral representation of displacements,
A simple comparison between two 3D time domain elastodynamic boundary element formulations becomes
I
'
Cki(Q,s) t-at
35
and as an average value over the time interval [t - At, t], for the new fundamental solution procedure with constant time approximation.
ui(s,T) dr
At
= Ii Iru*ki(Q,t;s, r)pi(Q, r)dF d'r
TIME CONVOLUTION OF FUNDAMENTAL VALUES
- ~ Ir ui ( Q, r) p~i ( Q , t; s, ~-)dF d~" + Ii In U~i(q't;s' "r)bi(q,r)dadr
(12)
As usual for boundary element formulations, eqn (11) can be modified by introducing appropriate space and time discretizations, as follows:
Cki (s) Ui (s, t) Nr =
(j) t~-I
-
When the time convolutions of eqn (13) are carried out, different expressions are achieved for each time approximation and the selected fundamental values. For the classical solution and adopting linear approximation we have
~°~(q, t,s) = A{B[K°~ - K°l~] + [K2°~- (D + E)K2a2~]}
Uki(a, t;s,T)OJ(Q)~b d r d F P
a
(17)
-,o/~(q,t,S) = I { j[KO~ - Kt0~]+ [K ° ~ - ~C2 KO3 ] Pki 21 J
J,.., J,"_pZ,(Q,t;s, "r)q~J(Q)~b~d'rdr U~Oj3 (13)
Similarly eqn (12) can be transformed into Cki(S) Ii_At~Nflt('t') d'r Uifl(s, t )Nr
u[i(Q,t;s,'r)¢J(Q)~ d~-dVPff fi
=
(j) to_~ -
L;' P*ki (J) o-~
(Q, t; s, T)~bJ(Ql~b~ dT d r
u;: (14)
where 0 ranges from 1 to the integration step number Nt and the summation rule is implied. In eqns (13) and (14), displacements and tractions are approximated by
Uio(Q j, t)
= "t'(JL/'(e) wa w~ r~'r°'J~ io
(15)
Pio(Qj, t)
*'h(J)M'(O)DaJfl = we~ wfl ~t iO
(16)
in which j stands for the spatial element, ~ba represents the space approximation functions, while ~ba gives the time approximation. Although any time approximation may be adopted together with the new fundamental solution, we have only implemented the constant element. This approximation has been found to be very efficient in terms of numerical result accuracy. Such a conclusion could be achieved after testing the formulation performance in many representative examples either with infinite or finite domain. Thus, adopting this simple approximation, eqn (14) is reduced to eqn (13) with /3 = 1 and if31 = 1. The convolution aspects of eqn (13) led us to the fact that the nodal displacement and traction values are considered at the instant t, for the classical formulation,
o[.:
.:1}
.8,
where A, B, C, D, E, I, J, M, N and O are functions written in terms of r and O, while Ki~ are expressions written in terms of time only, as shown in Appendix 2. Similarly, using constant time approximation we have tiff(q, t,s) = A{B[KI°I - K~21+ [K~l - (D + E)K2°z]} (19)
[
ro
:ff(q,t,s) = I J[KI°I - K°2] + K2°2-~21221 N[K~,] - O[K~21}
- O [~-~2K°] }
(20)
(21)
where the term ~-~o is used to compute the backward finite difference. For the new approach the final convoluted values are given by
~O(q, t, s) = A {B [KlOl_ K102]+ CKf, - (D + E)1(202} (22)
H. B. Coda, IV. S. Venturini
36
= ( q,t, sl = I
+u
where /3 varies from 1 to 2 for the classical solution with linear approximation, while it is assumed to be 1 for the new solution with constant approximation; the summation rule is implied. In the same way, now dealing with eqn (25), the final matrix representation for the classical solution with constant approximation is obtained, as follows:
- ICl°21+ L I °2 -
[
C 3 K,1° - N [ K f l +
Kfl]
- O[Kf2 ---~2 KJ~]}
(23)
The value Ks~ (for simplicity/3 = 1 is omitted) for eqns (19)-(23) is given in Appendix 2. The final representations of both formulations, obtained after convoluting eqn (13), become
-,0~(Q,t;s)4~(Q)dFP~ j Cki(S)Ui(s, I )NT = fJr uki j~" (J)
t;s)(bJ(Q) d r UiaoJ~
-
(24) Similarly, when the constant time approximation is adopted, together with the classical approach, we have
Cki(S)Ui(s, t)Nx = f
Jr
(J) f
~O(Q, t; s)4~(Q) dr e~'
-- JrIj/P °i(a' t; + JrIjl VP °(Q' t;
-
-
J
r(j)
oj
d r Uio
dr VJ
_~O_l(a,t;s)(aj(Q)d p Uio_l ~J (25)
The spatial integrals in the above equations ((24) and (25)) are performed by considering flat elements with eight nodes and quadratic shape functions. Singular integrals are carried out by a direct scheme based on the Kutt's quadrature formula, particularly written to compute finite part integrals./6'17 Non-singular integrals are performed by employing the well known element subdivision technique. Outside collocation points are considered by using a time translation technique described in Ref. 18.
ALGEBRAIC SYSTEM OF EQUATIONS After carrying out all spatial integrals, in eqn (24), we achieve an algebraic equation relating nodal displacements and tractions. As usual for boundary element formulations, selecting an appropriate number of collocation points 's' and writing their corresponding algebraic representations, the classical system of linear equations is achieved, as follows: o 0 0 0 H~U~ = G~P;~
(26)
(H 1 + HH1)U 1 :
G1p I
(27a)
(0 = 1)
H°U ° + [HH°U ° - HH0-1U0-1 ] = G°p °
(0 = 2, 3,..., Nt)
(27b)
where I-I[H° and I-II-I°-1 are given by the last integral terms in eqn (25). In order to improve the results due to the previous time stepping algebraic system, eqn (27), Manolis et al) 9 has proposed the following displacement correction U °-I = (U °-I + U°-2)/2
(0 = 2, 3 , . . . ,
Ut)
(28)
This procedure is very poor and could be avoided by taking into account the actual stage of the method. Another aspect that can be pointed out is the number of matrices computed at each time step. As can be seen, the BEM formulation used together with constant time elements generates two matrices per each time step, when employing the new solution, while three matrices are defined by the classical procedure. The linear approximation always generates four matrices at each step. To reduce the amount of computer memory, it is usual, in the classical approach, to adopt a combined time approximation model, linear for displacements and constant for tractions. This procedure avoids the use of finite differences and also gives three matrices at each time step, but requires rewriting eqn (26) in the following form: o o + H2U2 o o = Gopo H1UI
(29)
This scheme has been adopted by many authors and seems to be the best among the three possibilities regarding the use of the classical fundamental solution discussed here to model elastodynamic problems. It is important to point out that one way to improve the result accuracy, for both formulations, is to increase the time approximation in order to give smoother kernels. However, the authors prefer using other more appropriate load time distribution to define the unit concentrated impulse and consequently the fundamental solution, as can be seen in Ref. 20. It is worth mentioning that both fundamental solutions vanish after some period of time over the discretized domain. Thus, the maximum time step to compute new matrices is given by dmax. NtC' = C 2 A t '
u r n _ a~nax + 1 - - A----t C2
(30)
where d~ax is defined as the maximum possible length of
A simple comparison between two 3D time domain elastodynamic boundary element formulations
l
80 cm
l
37
P
E • ].1000 kg/(cm sz ) P
p = 0 . 0 0 2 koJ'cm3
t
Fig. 1. Rod geometry and loading. the discretized solid, the superscript 'Cl' indicates the classical solution procedure and 'n' indicates the new one.
EXAMPLES Two examples were selected to compare the performance of the time domain BEM elastodynamic formulation when the classical and the new fundamental solutions are adopted respectively. The objective of performing this numerical analysis is to make clear the differences between both formulations concerning accuracy and stability. In the first example, the classical prismatic rod problem is analysed when subjected to an impact at its free end. The geometry of the analysed solid is given in Fig. 1, together with the assumed material properties and the characteristics of the loading. Altogether, 34 quadratic elements were taken to discretize the solid boundary, using two different meshes, as illustrated in Fig. 2. It is important to stress that, for the classical approach, we have adopted linear time approximation for displacements and constant time approximation for tractions, while for the new formulation we have employed only constant time elements due to the advantages already reported. As we adopted different time discretizations for each formulation, which means different numbers of matrices, the difference between the observed execution times cannot be taken as an efficiency parameter, so it has been omitted. The time integration interval, At, taken for this analysis, varies from 0.001 s to 0.01 s. Some representative numerical results, obtained for both formulations, were selected and displayed in Figs 3, 4 and 5. Taking into account both presented discretizations and performing many other numerical analyses with
(a)
different time integration intervals, not presented here to save space, we have achieved, for this particular example, the stable time step range for both formulations, as illustrated in Fig. 6. Only the discretizations displayed in Fig. 2 have been adopted because it is enough to give stable results for the new approach. The results computed when using the classical approach and discretization 'a' cannot be shown due to the observed lack of convergence for any time interval adopted for the present analysis. It's important to note that the spatial integrals have been performed with the same numerical scheme for both formulations, using the same computer code. In this way, another code with a better numerical integration scheme2~ may be employed to find more stable time steps for the classical approach, which is not required by the new approach. As can be observed, the convergence At range for the new approach is very large. On the other hand, finding convergence for the classical solution was a very hard task due to its very narrow range. T h e numerical results presented above have been obtained by using singular formulations, 16 i.e. adopting load points along the boundary to write their integral representations and the corresponding algebraic relations. When the non-singular approach (outside collocations) was adopted together with the new solution, no significant changes were verified. Basically all numerical results are the same. Is On the other hand, using the classical formulation with non-singular c0110cations, the convergence range was completely missed. The second example taken to compare the performance of both formulation is the case of a dynamic load applied on a narrow finite strip defined over the free half-space surface. This problem has been previously analysed by Mansur and Brebbia 22 using a two-dimensional BEM approach and also by Tseng et al. 23 employing a finite difference code. T h e problem
-'-
Fig. 2. Discrctizations.
(b)
.r
H. B. Coda, W. S. Venturini
38
. . . . N
E U
z
-1
:. _
,,,, ....
OT-O.O02 (X VALU(
0|,
0.~
~.,o
o.~o
0
o.~o
0.00
0.10
•S Z ,X
0.20
0.30
t(s)
t(s) 2000
- EX. V A L U E
~
-
1
- EX. VALUE
1
G.
O. |
,
o.=
i i w ,
o.,o o.2o t(s)
o.~
o.~
,1
ii
w l l w
o.,o o.20 t(s)
o.~
Fig. 3. Reactions at the rod end; new fundamental solution; discretization 'a'.
2000
- (X.
},
VALUE
2000
~
"
ElL W4.IX
"
,
o.
(~. 0.00
0.10
0.20
0.30
0.~
0.10
t(s)
0.20
0.~
t(s)
Fig. 4. Reactions at the rod end; new fundamental solution; discretization 'b'. --'~.
2000
A
% ¢J
0T-0.00257~ (X. VALU(
~ooo
.1¢ v n 0 0.00
O, t 0
0.20
0.30
t(s)
N
-
2000.
DT-O.O02
- ( X . VALUE
-
2000 ~
J
r'"':
.... tx. 01'-0.003 v~u(
E U
5
1000 ¸
0,. 0
0.~
0
0.111
.2O
O.3O
W
Fig. 5. Reactions at the rod end; classical fundamental solution; discretization 'b'. ---NEW
C~V..
CLAS NEW
CON~-
1 Discretization
abt &t (=) (o)
Qoom~
Discretization
Fig. 6. Time BEM convergence ranges.
D
o~6
&t ( • ) (b)
A simple comparison between two 3D time domain elastodynamic boundary element formulations
PL
S.
1,1- t
•
llOms
t
r/
.
/
/
/
~" /
/
/
/
/--J
/"
I~i_/
39
Tlllg
/
/_/
,'/ /////// / // //.///// / / / / / / / / / / / / / / / / / / / /////// / / / / / / / / / / / /
•
-
xI
Fig. 7. Dynamic load applied on the half-space free surface. ~s
S~;5
,
to be analysed is depicted in Fig. 7. The half-space is initially at rest and then its surface is disturbed by a vertical load which is continuous in both space and time. The size of the narrow strip is given together with the adopted 200 linear element discretization, shown in Fig. 8. No residual stress has been assumed in the solid constituted by elastic material with Young's modulus equal to 200ksi and Poisson's ratio of 0.15. The material density was taken as p = 6-885 lbfs2/in 4. The analysis was conducted for the following time intervals: 3 ms, 5 ms, 10 ms, 15 ms and 20 ms. The results presented here and used for comparison were taken from the cases of 5ms (Aq) and 10ms (At2). Rather smaller values were adopted in Refs 23 and 22, 1ms and 3.65 ms respectively, to run the same example. No significant differences have been verified between the results obtained with both formulations when only singular collocations were considered. For this case,
LOAOEO REGION
I
il
ImWIllllll I ~ag~l I I I I I I I
]llllIIm
750'
IIIil11111 IIIII BIBillIlII 11111 I~!111111 I
Fig. 8. Half-space surface discretization.
"l
-8
-6
i)e O
w
;'1
I)D
N tZ
*4
gu =E
tU U 4[ .J O. (/1
o
At I' Art Ta*~I Mansur
100
ms
0
Q
in 0 2O
....;'/ 4O
-
GO
8O t
Fig. 9. Vertical displacement at point
D(0', 0~,70').
12Q
40
H. B. Coda, W. S. Venturini 6
I
x2~
¢
o°°. 7 0 °
oy
N
4
•
I.-
o°i
F
t
Z U.I
0
ILl ~J / n ¢n
-2
J
•
At;
D
At z Tse~ M onsur
o
in
0 20
40
60
80
100
ms
120
Fig. 10. Vertical displacements at point F (80', 0', 60').
almost the same results have been computed, also the observed convergence ranges were very similar. On the other hand, when only external collocations were adopted, no convergence was achieved for the classical approach. Figures 9 and 10 show the displacements computed at two internal points. In spite of comparing results obtained by different models (two and three dimensional hypotheses) a good agreement can be observed, mainly regarding the solution given by Tseng. It is also important to observe that although adopting very large time intervals, 5 and 10 times larger than the ones adopted by Tseng, the results were not disturbed.
CONCLUSIONS Numerical analyses have been carried out to show the performance of two time domain BEM formulations. The first formulation considered here was the classical one where the fundamental solution is given by assuming a unit concentrated impulse applied at some instant. For the second formulation the unit concentrated impulse adopted to derive the fundamental solution is assumed to be distributed along a certain period of time. In spite of being similar, they give completely different results when applied to some particular problems. The classical formulation gives very unstable results for problems characterized by a rather large ratio between boundary and volume. For these cases, for which the same numerical integration technique has been adopted for both formulation, the algebraic representations obtained after convoluting and discretizing the classical integral representations are not capable of giving stable and precise results. On the other hand, the new fundamental solution, being smoother, originates a much more convenient integral
representation, and consequently more appropriate algebraic relations, which give stable and precise results. It is important to confirm that the classical formulation also gives good results for problems with very small boundary/volume ratios, as for infinite domain cases, if only singular collocations were adopted. The algebraic representations given by the classical formulation also lose accuracy and stability if the collocations are not properly selected along the boundary.
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. & Brebbia, C. A. Transient elastodynamic using a time-stepping technique. In Proc. 5th Int. Sere. BEM Eng., ed. C. A. Brebbia, T. Futagami & M. Tanaka. Springer-Verlag, Berlin, 1983, pp. 677-98. 4. 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. 5. Manolis, G. BEM for soil-structure interactions. In Dynamic Soil Structure Interactions, ed. D. E. Beskos et al. Balkema, 1984, pp. 85-91. 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 Boundary Elements X, Vol. 4, ed. C. A. Brebbia. Computational Mechanics Publications, Southampton, Springer-Verlag, Berlin, 1988, pp. 387-99. 8. Dominguez, J. Boundary Elements in Dynamics. Computational Mechanics Publication, Elsevier Applied Science, London, 1993. 9. Manolis, G. & Davies, T. (eds.) Boundary Element Techniques in Geomechanics. Computational Mechanics
A simple comparison between two 3D time domain elastodynamic boundary element formulations Publications, Southampton, Elsevier Applied Science, London, 1993. I0. Coda, H. B. & Venturini,W. S. Further improvements on
three dimensional transient BEM elastodynamic analysis. Engng Anal Boundary Elements (submitted). 11. Coda, H. B. Three-dimensional transient analysis of stsuctures by a BEM/FEM combination. PhD thesis, University of S~o Paulo (BR), 1993 (in Portuguese). 12. Stokes, G. G. On the dynamical theory of diffraction. Trans. Camb. Phil. $oc., 1849, 9, 1. 13. Dominguez, J. & Abascal, R. Dynamics of foundations. In
Topics in Boundary Element Research: Applications in Geomechanics, Vol. 4, ed. Brebbia, C. A. Springer-Verlag, Berlin, 1987, pp. 27-75. 14. Wheeler, L. T. & Sternberg, E. Some theorems in classical elastodynamics. Arch. RationalMech. 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. Coda, H. B. & Venturini, W. S. Numerical evaluation of flat singular boundary elements in elastostatics and elastodynamics. Boundary Elements Commun. 1995, 6, 6-10. 17. Kutt, H. R. WISK 178: quadrature formulae for finite-part integrals. National Research Institute for Mathematical Sciences, Pretoria, 1975. 18. Coda, H. B. & Venturini, W. S. Non-singular timestepping BEM for transient elastodynamic analysis. Engng Anal. Boundary Elements, 1995, 15, 11-18. 19. Manolis, G. D. et al. Boundary element method implementation form three-dimensional transient elastodynamics. In Developments in Boundary Element Methods - - 4, ed. P. K. Banerjee & J. O. Watson. Elsevier, London, 1986, pp. 29-65. 20. Coda, H. B. & Venturini, W. S. Three dimensional elastodynamic formulation with domain initial conditions. In Boundary Elements XVII, ed. C. A. Brebbia et al. Computational Mechanics Publications, Southampton, 1995, pp. 195-202. 21. Dumont, N. A. A procedure for the semi-analytical evaluation of generally singular integrals that occur in the three-dimensional boundary element analysis. In Boundary Elements XVI1, ed. C. A. Brebbia et al. Computational Mechanics Publications, 1995, pp. 83-90. 22. 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. 23. Tseng, M. N. & Robinson, A. R. A transmitting boundary for finite-differenceanalysis of wave propagation in solids. Project no. NR 064-183, University of Illinois, Urbana, 1975.
+ r'rk[~212 7
(~1)
1
r)}
(
41
~22)]
+~-~2 6 r - ~ 2 2
(AI.1)
From eqn (AI.1), one can derive the corresponding tractions,
* { 2 I rirjrk Pki(q,r;s/f) : ~ -6C~ r5 x ~r H r _ r +2
×
-H
[rirjrk 6 r5
r
6ijrk+6kirj+6jkr(] r3 r)]
c,
[(
rirjrk +2r-V C r - Er
r:,, ( , r3
(g]
r)__~2ffr
If(
~ijrk+C~k_~_irj+c~jkri 1 r3 j
1-
)
C~ ~ ( r - - C rl ) ]
C21)
r3 6 (r --C22)+-~2 r r~r_ 6"iry+~jkri[ ( r C---22)]
(A1.2) where r = Is - q[ is the distance between the load point 's' and the field point 'q' and nj stands for the outward normal direction cosines. (b) E x p r e s s i o n s to c o m p u t e the new
fundamental values
The auxiliary values A~i and B*k used to define this fundamental solution, eqns (8) and (9), are given by 1
A*ki(q, r; s, O) -- 47rpAt APPENDIX 1: FUNDAMENTAL SOLUTIONS
X { ( 3rirk ~, r 3
(a) Classical fundamental values For the case of a concentrated unit impulse applied at an instant r = 0, the displacements are given by
u*ki(q,r;s,O)
=~p ~, r3
r2
r2
~k_.i) l [ ( r 2 )---~1 ~r 2
T2
r H('r - -~ll)
r
+--~-rkri[~-~lH(r--~l)-1---H(r~, C2]1
42
H. B. Coda, W. S. Venturini Bi,k(q,T;s,O)= x
nj
-N
4rrAt
-3C 2 5 r5
r3
t--T--
+
t--T--
[ ( r ) r ( r)]} - O 6 t-T-~-~2 +~22~ t-r--~22 (A2.2)
(,
where A, B, C, D, E, I, J, L, M, N and O are functions written in terms of r and 0 and can be easily recognized in eqns (AI.1) and (A1.2). The same procedure can be followed using the new fundamental values,
[ rirjrk 6ijrk+6ikrj-f-6jkri]
+ 2 6 r5 _
r3
r
C~
r
(t -- T)2 -
Aki : A
riryrk [ ( -)- ~r3 C !C~ 6 ( 7 \- _ r +2 r---~C 267-
)3 -
rk6' J[ (C;~] ,~ l-2kc VJ
(t-7-)=-
1(
+C-~1H
tSkirj-FtSjkri[
r3
H t - "r -
( r ) r( r)]} H 7--C-22 +C226 T--C22
+ EH
H t-r-
~)o, (
~-
H
r)
t - T -
t - ~- - ~
(AZ3)
(A1.4) B;i APPENDIX
2: FUNDAMENTAL
=
I
( [ ( r ~( t 1 2T)) 2 -
-
--
H
( t -~ ~- )-
VALUE
CONVOLUTION
In order to simplify the discussion regarding the convolution performance, let us consider the classical fundamental solution appropriately written in its reduced form, {
u*ki= A ( t - r ) n
[ (
+C6 t - 7 - -
(
+Ee t-7--~
Pki* = I
r/
r
-D6
r
)
+ M 6 t--T--
r
--
6( 1 \
--T
r -
t - 7- -
r)}
t--T) H t--7"--
(
C~H t--T--
-c-7,
r)]
t - r - - ~ l 1 - H t - r --~22 B
{, [(~,)
[(
(
[(r)
+L M ~- ~ - ~
- N H t--T--
-F 6 t--T--
[ ( r ) - O n t-z-~22
r ( r)]} +~226 t - r - ~ 2 2
(h2.1)
-- H
(
+ L 6 t - 7---~-22 - t - c 6! \
(l--T-- ~)1
J
7---
[ ( r ) C3~(t_T + M 6 t-- T---~22 - C--~I
r)]
-~1
r))] E
(A2.4)
After performing all convolutions the final achieved forms, eqns (17)-(23), are given, for each fundamental solution and for the chosen time approximations, in terms of coefficients KiT and Ki0 defined as follows.
A simple comparison between two 3D time domain elastodynamic boundary element formulations (a) For the classical solution and Hnear approximation
Kffl(r,t):Jtto°l(~)~(t-7"--~7) dT"
{t
K~(r,t)=J'i_(t-r)(to~A~)H(t-r--~)dr At"
r
if t - to_l >- ~
to,o_,1
2
>- t - to
r
r
~ < t - to
if-~>t-to_lor
_ ~l [to3_ 3tLlto + 27o_~]
r if t - to 3, -~7
(A2.9)
K3"r(r,t) 02
=
1 r
- ~-~ [ T A 1 3 _ 3to2_1ira + 2to3_~] if t -
c~/
r
if t - to-l > C.y >- t - t o
t - to
to_ 1 > -~7 >-
{1
[to (.r_to_l ~ ~ ( t - r - r_ ~ d r
j,,_, C-S?--/
r
0
r
0
43
if -~ >_t - to-i
if -~7 > t
r
-
or -~7 < t - to
to_ l
(A2.10)
(A2.5)
where 3' (3' = 1,2) specifies the wave velocity. K~(r,t)=Iii_(t-r)(~)H(t-~'--~)
dr
(b) For the classical solution and constant approximation
K~(r,t) = to (t-- T)H t - T -
d'r
to-I
[2t~ - 3ta_l tff + t3_ll
1
t [ TA2+t2-1 -At 2 - 6~
r
if t -
t2_, -- t 2 tAt + ~
to > -~7
TAto_1] t( TA - to-l) -t
[2TA3 - 3to_1TA 2 + t~_l]
if t - to >_ t2_l -- T A 2
r
if t-t0_1 > _ ~ > t - t o
2
r
if t - to-i >_-~7 > _ t - to
r
r
0
if ~
>_ t - to-l
(a2.11)
r
if ~ >_t - to-1 (A2.6)
=
6 t-r-
dr
to-1
K°l(r,t)=Jt~i_l(~-ff-)'(t-7"--~) dT"
1
r
if t - to-l --~ ~ ~-~ t -- to t
r
r
(to-t+-~v ) 0
o
i,o
r
K°~(r't)=
r
if t - to_ 1 <_~
0
te-~t° 16 t - r -
dr=Or
6 t-r-
dr r
if t - to-1 >_~ if t - to-1
r
(A2.13)
or -~7 < t -- to
ta-i
1 r --~t(t--to-l---~)
lor-~7
if t - to_ 1 >_- ~ > t - to
(A2.7)
K~7(r,t)=
(A2.12)
r
0 if-~7>t-to_
r
<--- ~
The velocity terms can be neglected by assuming K3~ = 0. Taking into account this assumption, the numerical results will appear rather dumped. In order to avoid this problem, as recommended in Ref. 19, we could introduce the velocity contribution by a finite difference procedure. It is made by replacing eqn (A2.13) by
>_ t - to K~(r,t) =
r
or - ~ <_ t - to (A2.8)
1 A---tt
0
if t - t o _ l > r
r --~7 > - t - t ° r
if-~7>t-to_lor-~7
0 (A2.14)
44
H.B.
Coda, W. S. Venturini
(c) For the new fundamental solution and constant approximation 1
1
r
r
r2
+(t~ - t3_,)/3]
if t - to _>
K20~=
if t - t o _>
\
At
if t - to_ 1 > C7 >- t - to
]
r
0
r
if ~
>_ t - to_ 1
(A2.16) Kj° =
1
r2 r
- t ( r A ~ - ¢ L , / + IrA ~ - t~_~)/~] J
if t-
/, to_j >_ ~ >_ t - t
if t - to_ l >_ - ~ > t - to o
g
0
if ~
g
> t - to_ 1 or ~
r if -~7 >- t - to_
< t - to
(A2.17) (A2.15)
where 7 = 1,2 a n d T A = t - r/CT.
ELSEVIER
Engineering Analysis with Boundary Elements 17 (1996) 33-44 Copyright © 1996. Published by Elsevier Science Limited Printed in Great Britain. All rights reserved 0955-7997/96/$15.00
A simple comparison between two 3D time domain elastodynamic boundary element formulations H. B. Coda & W. S. Venturini* Sao Carlos School of Engineering, University of S~o Paulo, Av. Carlos Botelho 1465, 13560-970 S~o Carlos, Brazil
(Received 9 March 1995; revised version received 9 September 1995; accepted 18 October 1995) The three-dimensional transient boundary element method for elastodynamic analysis is studied here to verify its performance concerning some required features when solving practical problems. The classical formulation was established more than 10 years ago and has already been applied several times, mainly in the context of infinite domain analysis. In this work the classical fundamental solution is compared with an alternative one in which a fundamental solution based on the unit impulse is assumed distributed over a time interval. Some interesting aspects of both solutions, regarding the result confidence for practical analysis, are also discussed in the paper. Key words: Boundary element method, time domain, elastodynamic analysis.
the consistence of the developed algorithms, the amount of memory required by them, the expected accuracy, and the stability ranges achieved for numerical applications.
INTRODUCTION The first application of the time domain boundary element method (BEM) apparently was made by Friedman and Shaw, l in 1962, in which the acoustic anti-plane problem was studied. In the context of two dimensional elastodynamic analysis, the first work was presented by Niwa et al., 2 in 1980; nevertheless, in 1983 only, a well established time stepping technique was presented by Mansur and Brebbia. 3 Karabalis and Beskos, 4 as well as Manolis, 5 in 1984, were the first to examine the general three-dimensional elastodynamic methodology in time domain BEM. In order to have a more complete overview of the elastodynamic BEM formulations the reader is invited to see the works written by Beskos 6 and Mansur. 7 Complete works have been recently published describing the BEM elastodynamic formulations, in details, and their practical applications, s'9 More recently, an alternative fundamental solution for BEM transient formulation was presented by Coda and Venturini. l°al This solution is also derived from the Stokes' state representing the responses due to a unit impulse accumulated along an interval of time, herein named, for convenience, the new fundamental solution. In this work a comparative analysis of the boundary element transient formulations using these two fundamental solutions is carried out. Important aspects of those formulations will be discussed: the coherence and
INTEGRAL REPRESENTATION OF DISPLACEMENTS FOR ELASTODYNAMICS Assuming that the bodies under analysis are isotropic and show linear elastic behaviors, the governing differential equations for the dynamic equilibrium, the so-called Navier-Cauchy equations, are given by: ( C 2 - C2)uj,ji + C~ui,jj + bi/p = ui
(1)
where bi represents the body forces, p is the media density, and C1 and C2 are the longitudinal and shear wave velocities, respectively. A family of possible fundamental solutions, Stokes' state, 4a2 is achieved by assuming, in eqn (1), the following body force term, b~i = f ( r ) a ( q - s)tki
(2)
in which s and q are the load and field points respectively, and f ( r ) gives the load time behaviour. Using this general solution together with the Graffi's theorem, 13a4 one can write the final displacement integral representation, in the same form several times shown in the literature. For a body of domain [2 and boundary F, with a quiescent past, submitted to a loading process, that began at a time r = 0 and is extended to a time t, the displacement integration representation
*To whom correspondence should be addressed. 33
34
H. B. Coda, W. S. Venturini
The new solution also exhibits a variable wave front width given by
is reduced as follows:
;o JoJ*
Cki(Q,s )
=
ui(s, 7 " ) f ( l -
7") dT-
Uki(Q,t-7";s/f)p~(Q,7")drdT"
- J; Jr Ui(Q, 7")p*ki(O, t - 7";s / f ) dr dr
+ IiI uzilq,,- 7";s/ytb,lq, ld dT
(3)
where u*ki and P*ki a r e fundamental values and Cki is an independent free term, similar to those used for elastostatic formulations.
FUNDAMENTAL SOLUTIONS
(4)
Adopting the body force distribution given in eqn (4), the fundamental displacement expressions are obtained, as shown in many works. 13-15 For completeness these expressions are displayed in Appendix 1. The new fundamental solution also represents the response to an infinite body subjected to a unit impulse, now distributed along an interval of time denoted by At, instead of applied at a specified instant. The expression adopted to represent this unit impulse is given by b*ki = [H(r) - H ( r - A t ) ] 6 ( q - S)6ki/At
where the subscript 'n' is used to indicate the new solution. The wave front width difference between both fundamental values, more relevant along the first time steps, is a reason for the improved behaviour observed in the numerical responses given by the BEM algorithm developed with the new solution adopted, in comparison with those computed by using the classical solution. It must be said that the unit impulse can be developed along any time interval which leads to a more general representation of eqn (5), b*ki = I n ( 7 ) -- U ( 7 -- M A t ) ] t S ( s
Both fundamental solutions analysed here are derived from the general Stokes' state, 12 assuming an appropriate time dependent load function f ( r ) to define the body force expression, eqn (2). The classical one represents a response of an infinite body subjected to a concentrated unit impulse applied at an instant r = 0. In this case eqn (2) becomes b*ki = 6(r)~(q - s)loki
(9)
In = (C, - C 2 ) N A t + C 2 A t
- q)~Ski/MAt
(10)
Using this expression, we are able to obtain the wave front widths needed for our purposes. The value M stands for the number of time intervals over which the unit impulse is distributed. It has been taken equal to one for the present analysis. The solution of some numerical examples can demonstrate that this assumption gives very precise results and requires a reasonable amount of computing time. Another reason is that it improves the numerical responses regarding the shape exhibited by the new solution with respect to time. The new fundamental values are represented by rather smoother curves in comparison with those given by the classical formula. This characteristic also explain why some time discretization restrictions, as described in the next section, must be enforced on the formulation based on the classical fundamental solution.
(5)
in which H ( r ) is the Heaviside function. For this case, the fundamental values, derived by Coda and Venturini, l°'1l can be represented as follows: U*ki(q, 7"; S, O) = [a~i(q , 7"; s, O) - A*ki(q, 7" -- A t ; s, 0)]
(6)
NUMERICAL ASPECTS First of all, we must rewrite expression (3) taking into account the types of load functions discussed earlier. Considering the load function adopted to derive the classical solution, eqn (4), we can achieve the following integral representation:
P*ki(q, 7-; s, O) = [B~i(q, r; s, O) - B~i(q , r - At; s, 0)]
(7) for which A*ki and BT,i are also given in Appendix 1. By analysing the classical fundamental solution we can verify that the corresponding wave front exhibits a variable width given by lc~ = (C1 - C 2 ) N A t
(8)
where the subscript 'Cl' is used to represent the classical solution and N is the number of time intervals over which the integral has already been performed.
Ck~(Q,s)ui(s,t) =
o
u*~(Q,t;s,r)pi(Q,r)drd7-
-
u~ (Q, r)p*k~(Q, t; s, r) dr d r
+Ji1
U*ki(q, t; S, r) bi (q, r) df~ dr
(11)
Following the same steps now using the new load function, eqn (5), the integral representation of displacements,
A simple comparison between two 3D time domain elastodynamic boundary element formulations becomes
I
'
Cki(Q,s) t-at
35
and as an average value over the time interval [t - At, t], for the new fundamental solution procedure with constant time approximation.
ui(s,T) dr
At
= Ii Iru*ki(Q,t;s, r)pi(Q, r)dF d'r
TIME CONVOLUTION OF FUNDAMENTAL VALUES
- ~ Ir ui ( Q, r) p~i ( Q , t; s, ~-)dF d~" + Ii In U~i(q't;s' "r)bi(q,r)dadr
(12)
As usual for boundary element formulations, eqn (11) can be modified by introducing appropriate space and time discretizations, as follows:
Cki (s) Ui (s, t) Nr =
(j) t~-I
-
When the time convolutions of eqn (13) are carried out, different expressions are achieved for each time approximation and the selected fundamental values. For the classical solution and adopting linear approximation we have
~°~(q, t,s) = A{B[K°~ - K°l~] + [K2°~- (D + E)K2a2~]}
Uki(a, t;s,T)OJ(Q)~b d r d F P
a
(17)
-,o/~(q,t,S) = I { j[KO~ - Kt0~]+ [K ° ~ - ~C2 KO3 ] Pki 21 J
J,.., J,"_pZ,(Q,t;s, "r)q~J(Q)~b~d'rdr U~Oj3 (13)
Similarly eqn (12) can be transformed into Cki(S) Ii_At~Nflt('t') d'r Uifl(s, t )Nr
u[i(Q,t;s,'r)¢J(Q)~ d~-dVPff fi
=
(j) to_~ -
L;' P*ki (J) o-~
(Q, t; s, T)~bJ(Ql~b~ dT d r
u;: (14)
where 0 ranges from 1 to the integration step number Nt and the summation rule is implied. In eqns (13) and (14), displacements and tractions are approximated by
Uio(Q j, t)
= "t'(JL/'(e) wa w~ r~'r°'J~ io
(15)
Pio(Qj, t)
*'h(J)M'(O)DaJfl = we~ wfl ~t iO
(16)
in which j stands for the spatial element, ~ba represents the space approximation functions, while ~ba gives the time approximation. Although any time approximation may be adopted together with the new fundamental solution, we have only implemented the constant element. This approximation has been found to be very efficient in terms of numerical result accuracy. Such a conclusion could be achieved after testing the formulation performance in many representative examples either with infinite or finite domain. Thus, adopting this simple approximation, eqn (14) is reduced to eqn (13) with /3 = 1 and if31 = 1. The convolution aspects of eqn (13) led us to the fact that the nodal displacement and traction values are considered at the instant t, for the classical formulation,
o[.:
.:1}
.8,
where A, B, C, D, E, I, J, M, N and O are functions written in terms of r and O, while Ki~ are expressions written in terms of time only, as shown in Appendix 2. Similarly, using constant time approximation we have tiff(q, t,s) = A{B[KI°I - K~21+ [K~l - (D + E)K2°z]} (19)
[
ro
:ff(q,t,s) = I J[KI°I - K°2] + K2°2-~21221 N[K~,] - O[K~21}
- O [~-~2K°] }
(20)
(21)
where the term ~-~o is used to compute the backward finite difference. For the new approach the final convoluted values are given by
~O(q, t, s) = A {B [KlOl_ K102]+ CKf, - (D + E)1(202} (22)
H. B. Coda, IV. S. Venturini
36
= ( q,t, sl = I
+u
where /3 varies from 1 to 2 for the classical solution with linear approximation, while it is assumed to be 1 for the new solution with constant approximation; the summation rule is implied. In the same way, now dealing with eqn (25), the final matrix representation for the classical solution with constant approximation is obtained, as follows:
- ICl°21+ L I °2 -
[
C 3 K,1° - N [ K f l +
Kfl]
- O[Kf2 ---~2 KJ~]}
(23)
The value Ks~ (for simplicity/3 = 1 is omitted) for eqns (19)-(23) is given in Appendix 2. The final representations of both formulations, obtained after convoluting eqn (13), become
-,0~(Q,t;s)4~(Q)dFP~ j Cki(S)Ui(s, I )NT = fJr uki j~" (J)
t;s)(bJ(Q) d r UiaoJ~
-
(24) Similarly, when the constant time approximation is adopted, together with the classical approach, we have
Cki(S)Ui(s, t)Nx = f
Jr
(J) f
~O(Q, t; s)4~(Q) dr e~'
-- JrIj/P °i(a' t; + JrIjl VP °(Q' t;
-
-
J
r(j)
oj
d r Uio
dr VJ
_~O_l(a,t;s)(aj(Q)d p Uio_l ~J (25)
The spatial integrals in the above equations ((24) and (25)) are performed by considering flat elements with eight nodes and quadratic shape functions. Singular integrals are carried out by a direct scheme based on the Kutt's quadrature formula, particularly written to compute finite part integrals./6'17 Non-singular integrals are performed by employing the well known element subdivision technique. Outside collocation points are considered by using a time translation technique described in Ref. 18.
ALGEBRAIC SYSTEM OF EQUATIONS After carrying out all spatial integrals, in eqn (24), we achieve an algebraic equation relating nodal displacements and tractions. As usual for boundary element formulations, selecting an appropriate number of collocation points 's' and writing their corresponding algebraic representations, the classical system of linear equations is achieved, as follows: o 0 0 0 H~U~ = G~P;~
(26)
(H 1 + HH1)U 1 :
G1p I
(27a)
(0 = 1)
H°U ° + [HH°U ° - HH0-1U0-1 ] = G°p °
(0 = 2, 3,..., Nt)
(27b)
where I-I[H° and I-II-I°-1 are given by the last integral terms in eqn (25). In order to improve the results due to the previous time stepping algebraic system, eqn (27), Manolis et al) 9 has proposed the following displacement correction U °-I = (U °-I + U°-2)/2
(0 = 2, 3 , . . . ,
Ut)
(28)
This procedure is very poor and could be avoided by taking into account the actual stage of the method. Another aspect that can be pointed out is the number of matrices computed at each time step. As can be seen, the BEM formulation used together with constant time elements generates two matrices per each time step, when employing the new solution, while three matrices are defined by the classical procedure. The linear approximation always generates four matrices at each step. To reduce the amount of computer memory, it is usual, in the classical approach, to adopt a combined time approximation model, linear for displacements and constant for tractions. This procedure avoids the use of finite differences and also gives three matrices at each time step, but requires rewriting eqn (26) in the following form: o o + H2U2 o o = Gopo H1UI
(29)
This scheme has been adopted by many authors and seems to be the best among the three possibilities regarding the use of the classical fundamental solution discussed here to model elastodynamic problems. It is important to point out that one way to improve the result accuracy, for both formulations, is to increase the time approximation in order to give smoother kernels. However, the authors prefer using other more appropriate load time distribution to define the unit concentrated impulse and consequently the fundamental solution, as can be seen in Ref. 20. It is worth mentioning that both fundamental solutions vanish after some period of time over the discretized domain. Thus, the maximum time step to compute new matrices is given by dmax. NtC' = C 2 A t '
u r n _ a~nax + 1 - - A----t C2
(30)
where d~ax is defined as the maximum possible length of
A simple comparison between two 3D time domain elastodynamic boundary element formulations
l
80 cm
l
37
P
E • ].1000 kg/(cm sz ) P
p = 0 . 0 0 2 koJ'cm3
t
Fig. 1. Rod geometry and loading. the discretized solid, the superscript 'Cl' indicates the classical solution procedure and 'n' indicates the new one.
EXAMPLES Two examples were selected to compare the performance of the time domain BEM elastodynamic formulation when the classical and the new fundamental solutions are adopted respectively. The objective of performing this numerical analysis is to make clear the differences between both formulations concerning accuracy and stability. In the first example, the classical prismatic rod problem is analysed when subjected to an impact at its free end. The geometry of the analysed solid is given in Fig. 1, together with the assumed material properties and the characteristics of the loading. Altogether, 34 quadratic elements were taken to discretize the solid boundary, using two different meshes, as illustrated in Fig. 2. It is important to stress that, for the classical approach, we have adopted linear time approximation for displacements and constant time approximation for tractions, while for the new formulation we have employed only constant time elements due to the advantages already reported. As we adopted different time discretizations for each formulation, which means different numbers of matrices, the difference between the observed execution times cannot be taken as an efficiency parameter, so it has been omitted. The time integration interval, At, taken for this analysis, varies from 0.001 s to 0.01 s. Some representative numerical results, obtained for both formulations, were selected and displayed in Figs 3, 4 and 5. Taking into account both presented discretizations and performing many other numerical analyses with
(a)
different time integration intervals, not presented here to save space, we have achieved, for this particular example, the stable time step range for both formulations, as illustrated in Fig. 6. Only the discretizations displayed in Fig. 2 have been adopted because it is enough to give stable results for the new approach. The results computed when using the classical approach and discretization 'a' cannot be shown due to the observed lack of convergence for any time interval adopted for the present analysis. It's important to note that the spatial integrals have been performed with the same numerical scheme for both formulations, using the same computer code. In this way, another code with a better numerical integration scheme2~ may be employed to find more stable time steps for the classical approach, which is not required by the new approach. As can be observed, the convergence At range for the new approach is very large. On the other hand, finding convergence for the classical solution was a very hard task due to its very narrow range. T h e numerical results presented above have been obtained by using singular formulations, 16 i.e. adopting load points along the boundary to write their integral representations and the corresponding algebraic relations. When the non-singular approach (outside collocations) was adopted together with the new solution, no significant changes were verified. Basically all numerical results are the same. Is On the other hand, using the classical formulation with non-singular c0110cations, the convergence range was completely missed. The second example taken to compare the performance of both formulation is the case of a dynamic load applied on a narrow finite strip defined over the free half-space surface. This problem has been previously analysed by Mansur and Brebbia 22 using a two-dimensional BEM approach and also by Tseng et al. 23 employing a finite difference code. T h e problem
-'-
Fig. 2. Discrctizations.
(b)
.r
H. B. Coda, W. S. Venturini
38
. . . . N
E U
z
-1
:. _
,,,, ....
OT-O.O02 (X VALU(
0|,
0.~
~.,o
o.~o
0
o.~o
0.00
0.10
•S Z ,X
0.20
0.30
t(s)
t(s) 2000
- EX. V A L U E
~
-
1
- EX. VALUE
1
G.
O. |
,
o.=
i i w ,
o.,o o.2o t(s)
o.~
o.~
,1
ii
w l l w
o.,o o.20 t(s)
o.~
Fig. 3. Reactions at the rod end; new fundamental solution; discretization 'a'.
2000
- (X.
},
VALUE
2000
~
"
ElL W4.IX
"
,
o.
(~. 0.00
0.10
0.20
0.30
0.~
0.10
t(s)
0.20
0.~
t(s)
Fig. 4. Reactions at the rod end; new fundamental solution; discretization 'b'. --'~.
2000
A
% ¢J
0T-0.00257~ (X. VALU(
~ooo
.1¢ v n 0 0.00
O, t 0
0.20
0.30
t(s)
N
-
2000.
DT-O.O02
- ( X . VALUE
-
2000 ~
J
r'"':
.... tx. 01'-0.003 v~u(
E U
5
1000 ¸
0,. 0
0.~
0
0.111
.2O
O.3O
W
Fig. 5. Reactions at the rod end; classical fundamental solution; discretization 'b'. ---NEW
C~V..
CLAS NEW
CON~-
1 Discretization
abt &t (=) (o)
Qoom~
Discretization
Fig. 6. Time BEM convergence ranges.
D
o~6
&t ( • ) (b)
A simple comparison between two 3D time domain elastodynamic boundary element formulations
PL
S.
1,1- t
•
llOms
t
r/
.
/
/
/
~" /
/
/
/
/--J
/"
I~i_/
39
Tlllg
/
/_/
,'/ /////// / // //.///// / / / / / / / / / / / / / / / / / / / /////// / / / / / / / / / / / /
•
-
xI
Fig. 7. Dynamic load applied on the half-space free surface. ~s
S~;5
,
to be analysed is depicted in Fig. 7. The half-space is initially at rest and then its surface is disturbed by a vertical load which is continuous in both space and time. The size of the narrow strip is given together with the adopted 200 linear element discretization, shown in Fig. 8. No residual stress has been assumed in the solid constituted by elastic material with Young's modulus equal to 200ksi and Poisson's ratio of 0.15. The material density was taken as p = 6-885 lbfs2/in 4. The analysis was conducted for the following time intervals: 3 ms, 5 ms, 10 ms, 15 ms and 20 ms. The results presented here and used for comparison were taken from the cases of 5ms (Aq) and 10ms (At2). Rather smaller values were adopted in Refs 23 and 22, 1ms and 3.65 ms respectively, to run the same example. No significant differences have been verified between the results obtained with both formulations when only singular collocations were considered. For this case,
LOAOEO REGION
I
il
ImWIllllll I ~ag~l I I I I I I I
]llllIIm
750'
IIIil11111 IIIII BIBillIlII 11111 I~!111111 I
Fig. 8. Half-space surface discretization.
"l
-8
-6
i)e O
w
;'1
I)D
N tZ
*4
gu =E
tU U 4[ .J O. (/1
o
At I' Art Ta*~I Mansur
100
ms
0
Q
in 0 2O
....;'/ 4O
-
GO
8O t
Fig. 9. Vertical displacement at point
D(0', 0~,70').
12Q
40
H. B. Coda, W. S. Venturini 6
I
x2~
¢
o°°. 7 0 °
oy
N
4
•
I.-
o°i
F
t
Z U.I
0
ILl ~J / n ¢n
-2
J
•
At;
D
At z Tse~ M onsur
o
in
0 20
40
60
80
100
ms
120
Fig. 10. Vertical displacements at point F (80', 0', 60').
almost the same results have been computed, also the observed convergence ranges were very similar. On the other hand, when only external collocations were adopted, no convergence was achieved for the classical approach. Figures 9 and 10 show the displacements computed at two internal points. In spite of comparing results obtained by different models (two and three dimensional hypotheses) a good agreement can be observed, mainly regarding the solution given by Tseng. It is also important to observe that although adopting very large time intervals, 5 and 10 times larger than the ones adopted by Tseng, the results were not disturbed.
CONCLUSIONS Numerical analyses have been carried out to show the performance of two time domain BEM formulations. The first formulation considered here was the classical one where the fundamental solution is given by assuming a unit concentrated impulse applied at some instant. For the second formulation the unit concentrated impulse adopted to derive the fundamental solution is assumed to be distributed along a certain period of time. In spite of being similar, they give completely different results when applied to some particular problems. The classical formulation gives very unstable results for problems characterized by a rather large ratio between boundary and volume. For these cases, for which the same numerical integration technique has been adopted for both formulation, the algebraic representations obtained after convoluting and discretizing the classical integral representations are not capable of giving stable and precise results. On the other hand, the new fundamental solution, being smoother, originates a much more convenient integral
representation, and consequently more appropriate algebraic relations, which give stable and precise results. It is important to confirm that the classical formulation also gives good results for problems with very small boundary/volume ratios, as for infinite domain cases, if only singular collocations were adopted. The algebraic representations given by the classical formulation also lose accuracy and stability if the collocations are not properly selected along the boundary.
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. & Brebbia, C. A. Transient elastodynamic using a time-stepping technique. In Proc. 5th Int. Sere. BEM Eng., ed. C. A. Brebbia, T. Futagami & M. Tanaka. Springer-Verlag, Berlin, 1983, pp. 677-98. 4. 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. 5. Manolis, G. BEM for soil-structure interactions. In Dynamic Soil Structure Interactions, ed. D. E. Beskos et al. Balkema, 1984, pp. 85-91. 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 Boundary Elements X, Vol. 4, ed. C. A. Brebbia. Computational Mechanics Publications, Southampton, Springer-Verlag, Berlin, 1988, pp. 387-99. 8. Dominguez, J. Boundary Elements in Dynamics. Computational Mechanics Publication, Elsevier Applied Science, London, 1993. 9. Manolis, G. & Davies, T. (eds.) Boundary Element Techniques in Geomechanics. Computational Mechanics
A simple comparison between two 3D time domain elastodynamic boundary element formulations Publications, Southampton, Elsevier Applied Science, London, 1993. I0. Coda, H. B. & Venturini,W. S. Further improvements on
three dimensional transient BEM elastodynamic analysis. Engng Anal Boundary Elements (submitted). 11. Coda, H. B. Three-dimensional transient analysis of stsuctures by a BEM/FEM combination. PhD thesis, University of S~o Paulo (BR), 1993 (in Portuguese). 12. Stokes, G. G. On the dynamical theory of diffraction. Trans. Camb. Phil. $oc., 1849, 9, 1. 13. Dominguez, J. & Abascal, R. Dynamics of foundations. In
Topics in Boundary Element Research: Applications in Geomechanics, Vol. 4, ed. Brebbia, C. A. Springer-Verlag, Berlin, 1987, pp. 27-75. 14. Wheeler, L. T. & Sternberg, E. Some theorems in classical elastodynamics. Arch. RationalMech. 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. Coda, H. B. & Venturini, W. S. Numerical evaluation of flat singular boundary elements in elastostatics and elastodynamics. Boundary Elements Commun. 1995, 6, 6-10. 17. Kutt, H. R. WISK 178: quadrature formulae for finite-part integrals. National Research Institute for Mathematical Sciences, Pretoria, 1975. 18. Coda, H. B. & Venturini, W. S. Non-singular timestepping BEM for transient elastodynamic analysis. Engng Anal. Boundary Elements, 1995, 15, 11-18. 19. Manolis, G. D. et al. Boundary element method implementation form three-dimensional transient elastodynamics. In Developments in Boundary Element Methods - - 4, ed. P. K. Banerjee & J. O. Watson. Elsevier, London, 1986, pp. 29-65. 20. Coda, H. B. & Venturini, W. S. Three dimensional elastodynamic formulation with domain initial conditions. In Boundary Elements XVII, ed. C. A. Brebbia et al. Computational Mechanics Publications, Southampton, 1995, pp. 195-202. 21. Dumont, N. A. A procedure for the semi-analytical evaluation of generally singular integrals that occur in the three-dimensional boundary element analysis. In Boundary Elements XVI1, ed. C. A. Brebbia et al. Computational Mechanics Publications, 1995, pp. 83-90. 22. 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. 23. Tseng, M. N. & Robinson, A. R. A transmitting boundary for finite-differenceanalysis of wave propagation in solids. Project no. NR 064-183, University of Illinois, Urbana, 1975.
+ r'rk[~212 7
(~1)
1
r)}
(
41
~22)]
+~-~2 6 r - ~ 2 2
(AI.1)
From eqn (AI.1), one can derive the corresponding tractions,
* { 2 I rirjrk Pki(q,r;s/f) : ~ -6C~ r5 x ~r H r _ r +2
×
-H
[rirjrk 6 r5
r
6ijrk+6kirj+6jkr(] r3 r)]
c,
[(
rirjrk +2r-V C r - Er
r:,, ( , r3
(g]
r)__~2ffr
If(
~ijrk+C~k_~_irj+c~jkri 1 r3 j
1-
)
C~ ~ ( r - - C rl ) ]
C21)
r3 6 (r --C22)+-~2 r r~r_ 6"iry+~jkri[ ( r C---22)]
(A1.2) where r = Is - q[ is the distance between the load point 's' and the field point 'q' and nj stands for the outward normal direction cosines. (b) E x p r e s s i o n s to c o m p u t e the new
fundamental values
The auxiliary values A~i and B*k used to define this fundamental solution, eqns (8) and (9), are given by 1
A*ki(q, r; s, O) -- 47rpAt APPENDIX 1: FUNDAMENTAL SOLUTIONS
X { ( 3rirk ~, r 3
(a) Classical fundamental values For the case of a concentrated unit impulse applied at an instant r = 0, the displacements are given by
u*ki(q,r;s,O)
=~p ~, r3
r2
r2
~k_.i) l [ ( r 2 )---~1 ~r 2
T2
r H('r - -~ll)
r
+--~-rkri[~-~lH(r--~l)-1---H(r~, C2]1
42
H. B. Coda, W. S. Venturini Bi,k(q,T;s,O)= x
nj
-N
4rrAt
-3C 2 5 r5
r3
t--T--
+
t--T--
[ ( r ) r ( r)]} - O 6 t-T-~-~2 +~22~ t-r--~22 (A2.2)
(,
where A, B, C, D, E, I, J, L, M, N and O are functions written in terms of r and 0 and can be easily recognized in eqns (AI.1) and (A1.2). The same procedure can be followed using the new fundamental values,
[ rirjrk 6ijrk+6ikrj-f-6jkri]
+ 2 6 r5 _
r3
r
C~
r
(t -- T)2 -
Aki : A
riryrk [ ( -)- ~r3 C !C~ 6 ( 7 \- _ r +2 r---~C 267-
)3 -
rk6' J[ (C;~] ,~ l-2kc VJ
(t-7-)=-
1(
+C-~1H
tSkirj-FtSjkri[
r3
H t - "r -
( r ) r( r)]} H 7--C-22 +C226 T--C22
+ EH
H t-r-
~)o, (
~-
H
r)
t - T -
t - ~- - ~
(AZ3)
(A1.4) B;i APPENDIX
2: FUNDAMENTAL
=
I
( [ ( r ~( t 1 2T)) 2 -
-
--
H
( t -~ ~- )-
VALUE
CONVOLUTION
In order to simplify the discussion regarding the convolution performance, let us consider the classical fundamental solution appropriately written in its reduced form, {
u*ki= A ( t - r ) n
[ (
+C6 t - 7 - -
(
+Ee t-7--~
Pki* = I
r/
r
-D6
r
)
+ M 6 t--T--
r
--
6( 1 \
--T
r -
t - 7- -
r)}
t--T) H t--7"--
(
C~H t--T--
-c-7,
r)]
t - r - - ~ l 1 - H t - r --~22 B
{, [(~,)
[(
(
[(r)
+L M ~- ~ - ~
- N H t--T--
-F 6 t--T--
[ ( r ) - O n t-z-~22
r ( r)]} +~226 t - r - ~ 2 2
(h2.1)
-- H
(
+ L 6 t - 7---~-22 - t - c 6! \
(l--T-- ~)1
J
7---
[ ( r ) C3~(t_T + M 6 t-- T---~22 - C--~I
r)]
-~1
r))] E
(A2.4)
After performing all convolutions the final achieved forms, eqns (17)-(23), are given, for each fundamental solution and for the chosen time approximations, in terms of coefficients KiT and Ki0 defined as follows.
A simple comparison between two 3D time domain elastodynamic boundary element formulations (a) For the classical solution and Hnear approximation
Kffl(r,t):Jtto°l(~)~(t-7"--~7) dT"
{t
K~(r,t)=J'i_(t-r)(to~A~)H(t-r--~)dr At"
r
if t - to_l >- ~
to,o_,1
2
>- t - to
r
r
~ < t - to
if-~>t-to_lor
_ ~l [to3_ 3tLlto + 27o_~]
r if t - to 3, -~7
(A2.9)
K3"r(r,t) 02
=
1 r
- ~-~ [ T A 1 3 _ 3to2_1ira + 2to3_~] if t -
c~/
r
if t - to-l > C.y >- t - t o
t - to
to_ 1 > -~7 >-
{1
[to (.r_to_l ~ ~ ( t - r - r_ ~ d r
j,,_, C-S?--/
r
0
r
0
43
if -~ >_t - to-i
if -~7 > t
r
-
or -~7 < t - to
to_ l
(A2.10)
(A2.5)
where 3' (3' = 1,2) specifies the wave velocity. K~(r,t)=Iii_(t-r)(~)H(t-~'--~)
dr
(b) For the classical solution and constant approximation
K~(r,t) = to (t-- T)H t - T -
d'r
to-I
[2t~ - 3ta_l tff + t3_ll
1
t [ TA2+t2-1 -At 2 - 6~
r
if t -
t2_, -- t 2 tAt + ~
to > -~7
TAto_1] t( TA - to-l) -t
[2TA3 - 3to_1TA 2 + t~_l]
if t - to >_ t2_l -- T A 2
r
if t-t0_1 > _ ~ > t - t o
2
r
if t - to-i >_-~7 > _ t - to
r
r
0
if ~
>_ t - to-l
(a2.11)
r
if ~ >_t - to-1 (A2.6)
=
6 t-r-
dr
to-1
K°l(r,t)=Jt~i_l(~-ff-)'(t-7"--~) dT"
1
r
if t - to-l --~ ~ ~-~ t -- to t
r
r
(to-t+-~v ) 0
o
i,o
r
K°~(r't)=
r
if t - to_ 1 <_~
0
te-~t° 16 t - r -
dr=Or
6 t-r-
dr r
if t - to-1 >_~ if t - to-1
r
(A2.13)
or -~7 < t -- to
ta-i
1 r --~t(t--to-l---~)
lor-~7
if t - to_ 1 >_- ~ > t - to
(A2.7)
K~7(r,t)=
(A2.12)
r
0 if-~7>t-to_
r
<--- ~
The velocity terms can be neglected by assuming K3~ = 0. Taking into account this assumption, the numerical results will appear rather dumped. In order to avoid this problem, as recommended in Ref. 19, we could introduce the velocity contribution by a finite difference procedure. It is made by replacing eqn (A2.13) by
>_ t - to K~(r,t) =
r
or - ~ <_ t - to (A2.8)
1 A---tt
0
if t - t o _ l > r
r --~7 > - t - t ° r
if-~7>t-to_lor-~7
0 (A2.14)
44
H.B.
Coda, W. S. Venturini
(c) For the new fundamental solution and constant approximation 1
1
r
r
r2
+(t~ - t3_,)/3]
if t - to _>
K20~=
if t - t o _>
\
At
if t - to_ 1 > C7 >- t - to
]
r
0
r
if ~
>_ t - to_ 1
(A2.16) Kj° =
1
r2 r
- t ( r A ~ - ¢ L , / + IrA ~ - t~_~)/~] J
if t-
/, to_j >_ ~ >_ t - t
if t - to_ l >_ - ~ > t - to o
g
0
if ~
g
> t - to_ 1 or ~
r if -~7 >- t - to_
< t - to
(A2.17) (A2.15)
where 7 = 1,2 a n d T A = t - r/CT.