Static and transient dynamic nonlinear stress analysis by the boundary element method with implicit techniques

Engineering Analysis with Boundary Elements 14 (1994) 65-74 © 1994 Elsevier Science Limited Printed in Great Britain. All fights re~rved 0955-7997/94/$07.00

ELSEVIER

Static and transient dynamic nonlinear stress analysis by the boundary element method with implicit techniques J. C. F. Telles & J. A. M. Carter COPPE, Universidade Federal do Rio de Janeiro, Programa de Engenharia Civil, Caixa Postal 68506, CEP: 21945-970, Rio de Janeiro, Brazil

This paper is concerned with the application of implicit techniques to solve static and transient dynamic elastoplastic problems. The current formulation keeps the inertial domain integral, which appears due to the use of the static fundamental solution, leading to the discretization of the entire domain. The time-domain problem is solved by a direct integration method (Houbolt scheme). Illustrative examples are presented at the end of the paper. Key words: Implicit procedures, transient dynamic analysis, elastoplasticity,

integral equation.

extended over the plastic region but the inertial one covers the entire domain. (ii) The time dependent problem is solved by the direct integration method with the use of the Houbolt scheme.

INTRODUCTION Alternative boundary element formulations have already been established to deal with inelastic problems for quite some time) '2 These approaches usually differ only in the way in which the inelastic terms are included in the integral equations and normally are classified as: initial strain, initial stress and fictitious surface and body forces. The solution routines, however, have been nearly always of an explicit nature and solve the inelastic problem by taking the inelastic contribution to the right-hand-side of the equations in a recursive way. More recently, a general approach which provides a selection of implicit solution techniques was developed with the aim of performing visco and inviscid plasticity analyses in a unified manner. 3-5 Later on, alternative versions of these implicit techniques were introduced in a formulation implemented by Carrer and Telles 6'7 to solve transient dynamic elastoplastic problems. Note that due to the nature of these problems, at each time step an elastoplastic analysis is performed under the application of the total load, whereas in the static case, the load is incrementally applied. Two aspects of the formulation can be pointed

Another way of performing transient dynamic elastoplastic analysis using the static fundamental solution has been accomplished by Kontoni and Beskos. 8 In their work, the inertial domain integrals are replaced by boundary integrals by means of a suitable transformation and interpolation of the domain variables (dualreciprocity approach9). The potentialities of such a formulation are demonstrated in the analyses of several examples presented and discussed by them. The use of time-dependent fundamental solutions certainly presents great advances in this area of research and, in spite of the complexity involved, is a matter that will deserve special attention from now on.

At the end of the paper, three examples are presented, allowing for an assessment of the current formulation.

out:

BOUNDARY ELEMENT FORMULATION

(i) Due to the nature of the static (Kelvin) fundamental solution the equations present a further domain integral: the integral related to the inelastic terms (initial stresses) is still only

The basic boundary integral equation (Somigliana's identity for the displacement field at a boundary point) associated with the initial stress formulation is defined 65

J. C. F. Telles, J. A. M. Carrer

66 as follows:

NUMERICAL ANALYSIS

Cij(~)Uj(~; t) : Jp Ui*j(~,X )pj(X; t) d F ( X ) - lr p~j(~; X)uj(X; t) d r ( X )

p J, u,.*A¢;x) j(x; t) dI2(X) +

x) k(x; t) da(X)

(1)

In the above equation, t is time, f~ and F represent the domain and the boundary of the body respectively and co. is the usual free coefficient of elastic analysis. The subscripts vary between 1 and 3 (3-D), and 1 and 2 (2-D). The fundamental elastostatic solution u~j(~;X) represents the displacement component at a field point X in the j direction due to a unit concentrated load applied at a source point ~ in the i direction. The same interpretation remains valid for the fundamental tensors Pi*j(~;X) and C~ki(~;X). Furthermore u:-, pj and ~ represent the displacement, traction and 'initial' (plastic) stress components;//j = 02uj/Ot 2 represents the acceleration components and p is the constant mass density. The stress integral equation, obtained by combining the correct derivatives of eqn (1) with respect to the coordinates of ~ E [2 (cij(~) = 1) according to the procedure described by Telles,1 has the following representation:

In order to perform the numerical analysis, the boundary and the domain are discretized by employing linear elements and linear triangular cells respectively. The displacements at internal points are calculated simultaneously with the boundary unknowns (displacements and tractions). Thus, the application of the discretized version of eqn (1) to all boundary nodes and internal points produces a system of equations with the following matrix representation:

[

I-ibb

ub

__

I(; ~b

[Qbb (3)

in which the superscripts b and d denote boundary and domain associated coefficients and unknowns. In the double superscript, the first one corresponds to the source point and the second to the field point positions. The identity matrix I represents the coefficients cu(~) related to the internal points. The simultaneous evaluation of the stress components at the boundary nodes and internal points is accomplished by the system below:

O'ij(¢; t) = [ U~.jk(¢,X)pk(X , t)dF(X) JP " .

~d = [G,dbJ {pb}- [H,abj{Ub}

Ir e,*jk(~;X)uk(X;t) dr(X) -- p J~ Ui~k(~;X)iik( X; t) d[2(X) +

q- I" 6iJkl(~; X)°~kl(X; t) d~(X) + gij[cr~kt(X; t)]

(2)

where the initial stress domain integral must be evaluated in the principal value sense and the presence of the 'free' term go.(oJ~kt) is due to the correct derivative of the initial stress domain integral in eqn (1). The boundary stress components are directly calculated as functions of the derivative of the interpolated displacements and tractions over 1". The static elastoplastic problem is treated as a particular case of the transient dynamic elastoplastic analysis. This is done by simply setting up in eqns (1) and (2) the mass density p equal to zero. Additional aspects of the static case will be discussed later on in the text. In general, the numerical implementation of eqns (1) and (2) are performed with the use of the Kelvin fundamental solution. 1° But alternative fundamental solutions, which satisfy particular boundary conditions over particular domain geometries, can equally be used.l

Q,db

(Q,dd + E,dd) [ o.p,d

(4)

Equations (3) and (4) can be written in a concise manner as Hu(t) = Gp(t) - Mii(t) + Qo'P(t)

(5)

and or(t) = G'p(t) - H'u(t) - M'ii(t) + (Q' + E')~rP(/)

(6) In eqns (5) and (6), matrices H, G, H ' and G' correspond to the boundary integrals while matrices M and M' correspond to the inertial domain integrals and matrices Q and Q' to the initial stress domain integrals. The free terms in eqn (2) are represented by matrix E'. The reader is referred to Telles 1 for additional details on boundary and initial stress numerical integration and Carrer 1° for inertial domain integration. The time-dependent problem is solved by the direct integration method. The Houbolt scheme was adopted

Static and transient dynamic nonlinear stress analysis taking advantage of the experience of its previous applications presented by Loettter n and Carrer. 1° Since this integration scheme is an implicit one, eqns (5) and (6) are written for time ( t + At) as functions of the previous solution values. The finite difference expansions are: 12 fi(t + At) =,6--~ [1 lu(t + At) - 18u(t)

67

In eqn (15), m(t + At) stands for the dynamic elastic solution for boundary unknowns and internal displacements and, in eqn (16), n ( t + A t ) stands for the dynamic elastic stress components. Furthermore: m(t + At) = A*(-0f*(t + At)

(17)

R = A*(-I)Q *

(18)

n(t + At) = f'(t + At) -- A'm(t + At) + 9u(t - At) - 2u(t - 2At)]

(7) - M'ii(t + At)

and

S = (Q' + E') - A'R.

(19) (20)

ii(t + At) = A--~ [2u(t + At) -- 5u(t) + 4u(t - At) - u(t - 2At)]

(8)

for the velocity and acceleration components, respectively. The substitution of eqn (8) in eqn (5) (provided t is substituted by t + At) gives [2M + At2H]u(t + At) - At2Gp(t + At) = M[Su(t) - 4u(t - At) + u(t - 2At)] + At2QtrP(t + At)

(9)

or, in a compact form H*u(t + At) - G*p(t + At) = h*(t) + Q*crP(t + At) (10) where h*(t) contains the contribution of the displacements previously determined multiplied by matrix M, i.e. h*(t) = h l ( t ) + h 2 ( t - A t ) + h a ( t - 2 A t

).

(11)

After the application of the boundary conditions, eqn (10) can be written as A*y(t + At) = f*(t + At) + Q*~rP(t + At)

(12)

where vector f*(t + At) is given by f*(t + At) = f(t + At) + h*(t).

(13)

In eqn (13), the contribution of the prescribed values at time (t + At) is included in the vector f(t + At). Equation (6) can be written as ~r(t + At) = - A ' y ( t + At) + f'(t + At)

The dynamic response is strongly conditioned by the correct choice of the time-step length. 1°'n33 In this implementation, the time step was considered acceptable when defined by the semi-empirical relation At = l/3CL

(21)

in which l is the smallest boundary element length and c L is the longitudinal wave propagation velocity: c2 = A + 2# P where A=

(22)

Eu (1 + v)(1 - 2v)

(23)

and # = G

E

2(1 +

(24)

represent the usual Lam6 constants 14 (E is Young's modulus and v is Poisson's ratio). As mentioned earlier, static analysis requires the particularizations described in what follows. The corresponding version of the system of equations related to the static Somigliana identity is constituted only by the upper part of the system of eqns (3). The computation of internal point displacements is not necessary or at least optional. The computation of the stress components follows the same pattern indicated in eqn (4). The corresponding versions of eqns (5) and (6) are therefore written as

-- M'ii(t + At) + (Q' + E')o'P(/+ At) (14) in which the contribution of the prescribed values is included in the vector f'(t + At). The solution of eqn (12) for the boundary unknowns (displacements and tractions) and displacements at internal points is y(t + At) = RtrP(t + At) + m(t + At).

(15)

After the substitution of (15) in (14), the stress components are calculated as follows o-(t + At) = So-P(t + At) + n(t + At).

(16)

Hu = Gp + Qtr p

(25)

and tr = G'p - H'u + (Q' + E')tr p.

(26)

Following the procedure described for the dynamic case, after the application of the boundary conditions, eqns (25) and (26) are written as Ay = f + Qo "p o" = - A ' y + f' + (Q' + E')o "p.

(27) (28)

J. C. F. Telles, J. A. M. Carrer

68

The solution of eqn (27) for the boundary unknowns is y = Rtr p + m

In addition, the increments of initial stress can be calculated by the relation 1

(29)

¢

do~ij = da~j - daq = -~Cqmnamnaktdakl

(40)

The substitution of eqn (29) in eqn (28) leads to tlr : S t r

p -+- n

(30)

In eqns (29) and (30), the vectors m and n represent the elastic solution corresponding to boundary unknowns and stress components respectively. The following matrix relations, entirely equivalent to those given by expressions (17)-(20), are valid for static analysis:

Hence, the increments of elastic stress can be computed directly from eqn (2) after gq substitutes for the free c o e f f i c i e n t gij: P P gij(O'kl) = gij(OrPl) "l- O'ij

(41)

For numerical implementation, the increments of elastic stress can be computed directly from eqn (6) if matrix E' is replaced by

m = A(-l)f

(31)

R = A(-1)Q

(32)

where I is the identity matrix. This gives (in incremental

n = f' - A ' m

(33)

form):

S = (Q' + E') - A'R

(34)

= E' + I

Atre(t + At) = SAtrP(t + At) + An(t + At)

(42)

(43)

in which = (Q' + E) - A'R

CONSTITUTIVE EQUA'I'IONS

(44)

and The incremental stress-strain relations for inviscid plasticity problems can be written as follows: ep

(35)

dcrij = C ijkldekl

An(/+

1

C ij kl : Cij kl -- "~ Cijmnanmaop Copkl

(36)

where

2Gu

With respect to the static analysis, eqn (43) is written A O "e = S A o "p + A n .

Ao = (37)

dY

Here, ae and e~p are the equivalent (or effective) stress and plastic strain, respectively, Y is the uniaxial yield stress and 6ij is the Kronecker delta. For the implementation of the initial stress formulation, eqn (35) can be conveniently simplified by introducing a fictitious 'elastic stress', defined by

(47)

(48)

where w is the associated load increment (kept constant during the analysis). The increments of elastic stress are therefore defined as

An = (Ai- Ai_l)n

(49)

An = wAon

(50)

or

(38)

Equation (35) can then be written in the following form: 1

Oe anax

Ai = Ai_l + ~OAo

dd"

d~ej = Cijk I dekl.

Y

and defining further values of the load factor as

,-y' = aqCijktakl + H' = aqdij + H' Hi=

(46)

In the dynamic case, the elastic stress increments are defined as a difference between two consecutive elastic responses (see eqn (45)). In the static case, however, the algorithm begins by scaling down the total elastic solution by a load factor )'o given by

Cijkl = 1 - 21., 6ij6kt + G(6ik6yl + 6il6jk)

OF _ Oae akl = Oakl &rkt

(45)

as

where dekt is the total strain increment and within the context of associated isotropic work hardening theory ep

At) = n(t + At) - n(t).

e

dcYij = do'ej - -~ Cijmnamnakl dO'kl

(39)

which means that true stresses can be computed from the corresponding elastic stresses in incremental form.

SOLUTION TECHNIQUES Implicit procedures have been successfully used to solve static and transient dynamic elastoplastic problems, as demonstrated by Telles and Carrer. 1'5'7 In the dynamic case, an elastoplastic problem is solved at the end of

Static and transient dynamic nonlinear stress analysis

69

Table 2. Number of iterations u d C P U time for w = 40% 3.6

R

A 0.6771 0-8705 1.0000 o

•5 k g f / m m

0 : plastified

nodes/points

J

Fig. 1. Perforated aluminium strip.

CPU

( k - 1 ) a e .~ (k)Acre = s [ ( k - 1 ) a P

-~- (k)DP (k)Ao-¢]

+ n(t) + An(t + At)

(52)

A residual vector can be defined as ¢(t) = StrP(t) + n(t) - tre(t)

(53)

Substituting eqn (53) into eqn (52) and rearranging, yields SPAtre(t + At) = An(t + At) + ¢(t)

(54)

where the nonlinear matrix S p is given by S p = I - SD p

(55)

Equation (54) can now be solved for each time step requiring, in principle, no iterations. In practice, however, if the time-step length is not sufficiently small, a few iterations would have to be performed. In this work, the solution schemes tested are the following:

Table 1. Number of iterations and C P U time for w = 20%

A 0-5804 0"6771 0"7738 0'8705 0"9673 1"0000 CPU

SIS

MNR

NR2

10

I0

3

1

28 24 39 83 65

26 8 33 36 7

6 4 5 5 4

1 1 1 1 1

2'44

1"00

2"50

2"12

NR2

PI

28 33 39

9 3 7

1 1 1

3.17

1-00

3.37

SPAtre(t + At) = An(t + At) + ~(t)

(56)

and for the subsequent iterations, until ~ becomes vanishingly small (i.e. ~b ~ 0 within some predetermined tolerance), the above equation is written as

(51)

in which D p is a quasi-diagonal matrix. At a fixed time (t + At), the substitution of eqn (51) into the accumulated form of eqn (43) allows for the application of an iterative elastoplastic algorithm whose step k can be defined by:

3.48

MNR

(a) Pure incremental (PI): eqn. (54) is applied only once for each time step without iterations. Matrix S p is always updated as a function of the previous solution. (b) Newton-Raphson: matrix S p is updated after each iteration (by definition). For the first iteration of a given time step, eqn (54) is written as

each time interval in which the overall time of analysis is discretized. The application of eqn (40) to all boundary nodes and internal points can be represented in matrix form as A O "p = DPo"e

SIS 28 42 70

SPAtTe(t --1-At) = ~b(t + At)

(57)

In practice, however, only alternative versions of this procedure have been utilized. The most common version and perhaps the most efficient, is that known as the modified Newton-Raphson (MNR). Here, matrix S p is updated only once at the start of each time step (it is re-assembled as a function of the converged previous solution) and kept unaltered for the iterations. Another alternative version, designated by NR2, based on the double updating of S p (in the first and second iterations) was also implemented with the purpose of enriching the discussion concerning the number of iterations versus computational effort. It is worth noting that vector ¢(t) is preserved in eqn (56) only to avoid the effects of a possible error propagation (~/,(t) in the vanishing residual vector of the previous time solution). (c) Standard initial stress (SIS): this case is the explicit one extensively tested by Telles 1 and Venturini, 15 being treated here as a particular case of (b) in which S o is never updated (i.e. S p = I).

Table 3. Stress components at the centre section for ~ = 20%

PI

~y/r x/R

SIS

MNR

NR2

PI

0"00 0"30 0"50 0"70 0"90 1"00

0"3709 1"0281 1"1070 1'1037 1"0748 1"0146

0-3722 1-0361 1"1056 1'1013 1"0720 1"0131

0"3716 1'0345 1"1056 1"1011 1"0714 1'0148

0"3721 1"0331 1"1062 1"1017 1"0750 1"0148

J. C. 1;'. Telles, J. A. M. Carrer

70

Table 4. Stress components at the centre section for w = 40%

o : plastified nodes/points

YI x/R

SIS

MNR

NR2

PI

0.00 0.30 0.50 0.70 0.90 1.00

0.3705 1.0228 1.1078 1-1056 1.0786 1.0144

0.3718 1-0322 1.1056 1-1021 1.0754 1.0149

0.3720 1.0307 1.1060 1.1020 1.0737 1.0148

0-3431 0.8627 1.1147 1.1075 1.0754 1.0248

As indicated by Carrer and TeUes, 7 the use of implicit procedures needs a careful definition of the loading/ unloading situations, at each stress point, in order to predict the appropriate tangent matrix S p, i.e. for the first iteration of a given time step, the columns of S p associated with a stress point that now behaves elasticaUy must be set equal to the corresponding columns of the identity matrix. In the static case, the elastoplastic problem is solved incrementally, with the load factor varying from Ao to one, i.e. Ao < A < 1. For the ith load increment, the equivalent version of eqn (52) becomes (k-1)tye + (k)ACre = s[(k-1)tyP + (k)DP (k)Ao-e] + Ai_ln + A n

(58)

I 1.0 ksi

0.40 ksi

x~ Fig. 3. Boundary element and internal cell discretization.

EXAMPLES The first and second examples refer to static analysis, and the third one deals with elastodynamics. Example 1. Perforated aluminium strip (plane stress)

In this classical plasticity problem, the material is assumed to obey the Von Mises yield criterion and has the following parameters: ] E -- 7000 kgf/mm 2 ao = Y + 0"032Eeep

A residual vector can be defined as

Y = 24"3 kgf/mm 2 ~bi_ 1

=

(59)

~ 1 S o ' P _ I + A i _ l n -- O'i_

v = 0.2 Substituting eqn (59) into eqn (58) and rearranging, yields SPAo "e = An + ~bi_l

(60)

where, as before, the nonlinear matrix S p is given by Sp

=

I - SD p

(61)

The definition of the techniques based on eqn (60) follow the same procedures developed for the dynamic case.

The discretization is depicted in Fig. 1, which also presents the plastic points obtained for the total applied load. Tables 1 and 2 present the number of iterations corresponding to each load factor and the associated relative overall computer time for the load increments w = 20% and w = 40%, respectively (not including the time O//C~ -7.50

c~

1.oo x/R

0 -5.00 0.75

"'"

0 0 0 0 0 SIS w=20~ ooooo Pt w=20~ A,A'AAA PI w=40~ r~

0.50

.

-2.50 zx

.

.

.

.

.

Baker (FEM) (.FEM,) (BEM)

0

0.25

×/R 0.00

ooo

o~o

40, 0.40

o~o

o~o

Fig. 2. Results at the centre section.

1.6o

1.2o

0.00 1.0

2.O

Fig. 4. Final stress components along the horizontal section.

Static and transient dynamic nonlinear stress analysis Y

Table 5. Final stress components

x/R

SIS

MNR

NR2

PI

0.0042 -0.3828 -0.4874 -0.5164 -0.4981 -0.4715 -0.4554 -0.4438 -0-4354 -0.4291

0.0023 -0.3793 -0.5064 -0.5185 -0.4996 -0.4725 -0.4561 -0.4438 -0.4357 -0.4293

-0-9547 -1.8918 -1.7138 -1.3608 -1.1909 -1.1100 -1-0769 -1.0569 -1.0440 -1.0350

-0.9593 -1.8975 -1.7218 -1.3543 -1.8888 -1.1088 -1.0760 -1.0563 -1.0435 -1.0346

ax 1.00 1.25 1.50 1.90 2.40 3.00 3.50 4.00 4.50 5.00

0.0023 -0.3758 -0.4886 -0.5166 -0.4981 -0.4715 -0.4554 -0.4438 -0.4353 -0.4290

0.0035 -0.3838 -0.4876 -0.5162 -0-4979 -0.4715 -0.4554 -0.4438 -0.4353 -0.4290

p=l.0ksi

% 1.00 1.25 1-50 1.90 2-40 3-00 3.50 4.00 4.50 5.00

-0.9554 -1.8888 -1.7161 -1.3611 -1.1911 -1.1101 -1.0769 -1.0570 -1.0440 -1.0351

-0.9567 -1.8914 -1.7137 -1.3605 -1.1909 -1.1100 -1.0769 -1.0569 -1.0440 -1.0350

required for the generation of matrices R, S, m and n given by eqns (31)-(34)). The good agreement of the results at the centre section is shown in Tables 3 and 4. Figure 2 depicts a comparison between the final results obtained with the PI technique and those obtained with the SIS technique along this section. It is interesting to note that the PI technique is faster than the others: the computer time required by this technique with w = 20% is 71% of the one required by the SIS technique with ~o = 40%. Example 2. Deep circular tunnel

This plane strain problem was analysed with the assumption of the infinite domain being initially subjected to a uniform stress field of 1.0ksi vertical and 0.4ksi in both horizontal directions (K0 = 0.4). External loads corresponding to the relaxation of this in situ stress field were applied on the surface of the opening. Table 6. Number of iterations and CPU time

A 0"20 0"30 0'40 0"50 0'60 0'70 0-80 0"90 1'00 CPU

SIS 1 16 19 27 39 56 117 87 92 2.29

MNR

NR2

1 16 14 17 27 36 85 42 50

1 8 3 10 7 8 46 11 38

2.43

71

2"57

PI 1 1 1 1 1 1 1 1 1 1"00

~t

Fig. 5. Loading of the cavity. The anlaysis was carried out under the D r u c k e r - P r a g e r simulation of the M o h r - C o u l o m b yield criterion) Finite element analyses of this problem were performed by Reyes and Deere 16 and Baker et al. 17 and boundary element analyses by Telles I and Venturini. 15 The material (rock) constants are: E = 500.0 ksi c' = 0.28 ksi v = 0.2 ~b' = 30 °

(cohesion) (angle of internal friction)

The boundary element and internal cell discretization is presented in Fig. 3. The stress components along the horizontal section computed at the end of the relaxation process (SIS technique) are presented in Fig. 4 in conjunction with the corresponding F E M results. The BEM results are entirely equivalent to those presented by Telles. 1 A critical comment on the discrepancy observed in Fig. 4 between the BEM and F E M results can be found in the same reference. In order to demonstrate the accuracy of the implicit techniques, the stress values along the horizontal section, obtained with the use of alternative techniques, are presented in Table 5. The efficiency of the implicit techniques can be observed in Table 6 where a comparison between the number of iterations

I

a(R,O) B (2.02 R,O)

r~

C(3.45 R,O)

-A

-e

-C

x~

Fig. 6. Boundary element and internal cell discretization.

72

J. C. F. Telles, J. A. M. Carrer 0.5

1.50

o r

. .+ . . . . . /

1.00

.........

/

0,0 0.50

-

-0.5

P

,,'

I ~

~

-

Ref. 18 BEM Ref. 18 BEM Ref. 18 BEM

A

point

ooooo point A --point B ° o o o o point B . . . . . . point C **-** point C

I~

0.00

,~

- oooo* ..... ...... ~,oooo

-0.50

000+o

radial e l a s t i c SIS PI circumferential elastic SIS PI

-1.0 -1.00 ct/R

ct/R -1,5

['~ 0.0

5.'o

-1.50

15'o

+0%

0.0

5.~0

I Ol.O

151.0

Fig. 7. Radial stress components at A, B and C.

Fig. 9. Stress components at point A: c' = 0-90ksi.

corresponding to each load factor and the associated relative overall computer time is indicated. Tables 5 and 6 were compiled for ~v = 100%. It should be pointed out that due to the presence of the in situ stress field, the load factor corresponding to first yield (Ao) could not be computed directly (i.e. in linear fashion) as stated in eqn (47). Thus, a convenient starting value was attributed to Ao such that the analysis began in the elastic range. Here Ao was assumed to be equal 0"1. This value is smaller than the correct one and at A = 0.2 the response was still elastic.

corresponding results obtained by Chow and Koertig 18

Example 3. Circular cavity

fundamental solution, showed a good agreement with those presented by Chow and Koenig 18 and Fu 19 (FE analysis). The advantage of not performing domain discretization becomes more evident in an infinite domain problem, though in the formulation presented here the use of double symmetry provided a significant reduction in the domain discretization effort (no elements over symmetry axes). Since the results presented in Figs 7 and 8 are considered acceptable, a plane strain elastoplastic

This problem consists of a uniform internal pressure suddenly applied to the surface of a circular cavity and kept constant in time as seen in Fig. 5. The boundary element and internal cell discretization is depicted in Fig. 6. The time history of the radial and circumferential stress components at points A, B and C displayed in Fig. 6 is presented in Figs 7 and 8 together with the 1.50

%

(method of characteristics). The latter is for plane stress analysis (E = 100.0 ksi; v = 0.3). Therefore, in order to compare the results, a plane strain boundary element analysis was carried out with modified Young's modulus and Poisson's ratio as follows: E = 94.6769 ksi = 0.2308 p = 3-5 slug/ft 3 The results of another elastic BEM analysis, carried out by Mansur 13 with the use of the time dependent

1.50

a

1.00 /

1.00

/ ]

__

point A A point B oo..o point B . . . . . . point C .+p°'nt C

}/

o o o o o point

___

+I

0.50

~ t

~

~

*

_

"

18 BEM Ref. 18 BEId Ref. 18 BEM Ref.

0.50

0.00

--radial oooo= SIS ~8~.

...... -0.50

ooooo o <~,0++ <,

elastic

PI

circumferential elastic SIS PI

0.00

- 1.00

ct/R -0.50

o+o

5.b

io'.o

15:o

Fig. 8. Circumferential stress components at A, B and C.

ct/R - 1.50

0.0

5.1010.01

15~.0

Fig. 10. Stress components at point A: c' = 0.70 ksi.

Static and transient dynamic nonlinear stress analysis 0.50

0.50

ff

73

cr

0.00

0.00

+ ,tz

-0.50

--radial elastic o°o0° SIS .... A PI ...... circumferential ooooo SIS ooooo PI

t ~

--o*oo° +*,~* ...... ooooo

-0.50

elastic

ooooo

radial elastic SIS PI circumferential SIS PI

ct/R

ct/R -1.00

-1.00 0.0

lo'.o

5.()

elastic

15'.o

5.b

0.0

Io'o

+5'.o

Fig. 11. Stress components at point B: c' = 0.90 ksi.

Fig. 13. Stress components at point C: c' = 0.90ksi.

analysis was carried out with the same modified material constants. The Drucker-Prager simulation of the Mohr-Coulomb yield criterion was adopted with two values of the cohesion of the material: c' = 0.90 ksi and c' = 0"70 ksi. In both cases, ~b' = 30 ° was maintained. The time history of the radial and circumferential elastoplastic stress components at points A, B and C is presented in Figs 9-14. As can be observed from the above results, for c ' = 0.70ksi all solution techniques predict instability of the cavity after a certain level of plastification, leading to unlimited large displacement values over the boundary, as depicted in Fig. 15. The same did not occur with the use of c' -- 0.90 ksi allowing for the complete analysis up to a later time. Note that for this higher c' value, points B and C are still outside the plastic region (Figs 11 and 13). It is important to note that while in the static example the domain discretization can be restricted to the plastic region, in the dynamic example the domain discretization is conditioned by the overall time of

analysis (i.e. the discretization should always reach the undisturbed outer region of the domain within the total time of analysis).

In order to perform static elastoplastic analysis by the boundary element method, only the internal plastic region needs to be discretized. This advantage of the method is not taken into account in the formulation developed here, since it is necessary to discretize the entire finite domain to perform the time-dependent analysis or, in the case of an infinite domain, extend the discretization in such a way that the dynamic response is not affected by undesirable wave reflexions produced by a fictitious outer boundary (see example 3). The simplicity of the kernels in the inertial domain integrals, however, partially counterbalances this disadvantage. Furthermore, former experience from static elastoplastic analysis renders this formulation attractive

o o

0.50

0.00

CONCLUSIONS

0.50

,.%,~

0.00

~L i

¢ ...- ........ ."

a~ /

. ..... ..

__--... -

..

~-

j/i

~ g

-0.50

--,%,¢ i l ~

/

--°°0°° ..... ...... ooooo

g

ooooo

radial elastic SIS PI circumferential SIS PI

--°oo,o * , , +,~ * ......

-0.50 elastic

ooooo

ooooo

radial elastic SIS PI circumferential SiS PI

ct/R

ct/R - 1.00

o.o

5.b

+o'o

15'o

Fig. 12. Stress components at point B: c' = 0.70ksi.

elastic

-

1.00

0.0

5/0

1010

15'.0

Fig. 14. Stress components at point C: c' = 0.70ksi.

74

J. C. F. Telles, J. A. M. Carrer 0.25

u

0.20

0.15 e e o

0.10

o a

a

ii 0.05

ooooo

elastic

oD~oo

elostoplostic elostoplestic

.....

c' c'

= =

0.90 0.70

ksi ksi

e

ct/R 0.00 0.0

5'0

loto

15to

Fig. 15. Elastic and elastoplastic radial displacements at boundary node A. and the so called dual-reciprocity procedure can be alternatively implemented s to avoid domain integrals of the inertial term. With respect to the Houbolt integration scheme, it seems to be an adequate technique, perhaps ideal for dynamic BE solutions. However, it is important not only to use techniques taken directly from the F E M but also to search for alternative ones, which could be more appropriate for the BEM. Even the implementation of the dual-reciprocity procedure calls for alternative time integration schemes. In addition, this very problem could be solved with the use o f time-dependent fundamental solutions; in this case, boundary integral equations are established directly with no inertial domain integrals (the only domain integrals being those related to the inelastic terms) and there is no need for the direct time integration schemes adopted (i.e. time interpolation functions can be directly integrated instead). These alternatives, however, are still the subject of further research.

REFERENCES 1. Telles, J. C. F. The Boundary Element Method Applied to Inelastic Problems, Lecture Notes in Engineering, Vol. 1. Springer-Verlag, Berlin, Heidelberg, 1983. 2. Mukherjee, S. Boundary Element Methods in Creep and Fracture. Applied Science Publishers, London, 1982. 3. Telles, J. C. F. On inelastic analysis algorithms for boundary elements. In Advanced Topics in Boundary Element Analysis, eds T. A. Cruse, A. B. Pifko & H. Armen. AMD-Vol. 72, 1985, pp. 35-44.

4. Telles, J. C. F. & Carrer, J. A. M. Implicit solution techniques for inelastic boundary element analysis. In Boundary Element X, Proceedings of the lOth International Conference on Boundary Element Methods in Engineering, Vol. 3, ed. C. A. Brebbia. Springer-Verlag, Berlin, 1988, pp. 3-15. 5. Telles, J. C. F. & Carrer, J. A. M. Implicit procedures for the solution of elastoplastic problems by the boundary element method. Mathematical Computational Modelling, 1991, 15, 303-11. 6. Carrer, J. A. M. & Telles, J. C. F. Transient dynamic elastoplastic analysis by the boundary element method. In Boundary Element Technology VI, Proceedings of the 6th Int. Conf. on Boundary Element Technology, ed. C. A. Brebbia. Computational Mechanics Publications, Southampton, Elsevier, London, 1991, pp. 265-77. 7. Carrer, J. A. M. & Telles, J. C. F. A boundary element formulation to solve transient dynamic elastoplastic problems. Computers and Structures, 1992, 4, 707-13. 8. Kontoni, D. P, N. & Beskos, D. E. Transient dynamic elastoplastic analysis by the dual reciprocity BEM. Engineering Analysis with Boundary Elements, 1993, 12, 1-16. 9. Nardini, D. & Brebbia, C. A. A new approach to free vibration analysis using boundary elements. In Boundary Element Methods in Engineering, ed. C. A. Brebbia. Springer-Verlag, Berlin, 1982. 10. Carrer, J. A. M. Implicit techniques to solve static and transient dynamic elastoplastic problems by the boundary element method. D Sc thesis, COPPE/UFRJ, Brasil, 1991, (in Portuguese). 11. Loeffler, C. F. An alternative formulation of the boundary element method applied to scalar field problems. D Sc thesis, COPPE/UFRJ, Brasil, 1988, (in Portuguese). 12. Bathe, K. J. Finite Element Procedures in Engineering Analysis. Prentice-Hall, Englewood Cliffs, NJ, 1982. 13. Mansur, W. J. A time-stepping technique to solve wave propagation problems using boundary element method. PhD thesis, University of Southampton, 1983. 14. Achenbach, J. D. Wave Propagation in Elastic Solids. North-Holland, Amsterdam, London, 1973. 15. Venturini, W. S. Boundary Element Method in Geomechanics, Lecture Notes in Engineering, Vol. 4. Springer-Verlag, Berlin, Heidelberg, New York, Tokyo, 1983. 16. Reyes, S. F. & Deere, D. U. Elastic-plastic analysis of underground openings by the finite element method. Proceedings of the 1st Int. Congr. Rock Mechanics, 1966, pp. 477-83. 17. Baker, L. E., Shandu, R. S. & Shieh, W. Y. Application of elastoplastic analysis in rock mechanics by finite element method. Proceedings of the 11th Syrup. Rock Mechanics, ed. Somerton, 1969, pp. 237-51. 18. Chow, P. C. & Koenig, H. A. A unified approach to cylindrical and spherical elastic waves by method of characteristics. Transactions of ASME, Journal of Applied Mechanics, 1966, series E, 33, pp. 159-67. 19. Fu, C. C. A method for the numerical integration of the equations of motion arising from a finite element anlaysis. Transactions of ASME, Journal of Applied Mechanics, 1970, series E, 37, pp. 599-605.