- Email: [email protected]
A fracture mechanics based seismic analysis of concrete gravity dams using discrete cracks
Engineering Fracture Mechanics
Vol. 35, No. l/2/3, Pp. W-598, 1990
Printed in Great Britain.
0013-7944/90 $3.00 + 0.00 Perpmon Press pk.
A FRACTURE MECHANICS BASED SEISMIC ANALYSIS OF CONCRETE GRAVITY DAMS USING DISCRETE CRACKS MOHAMED
L. AYARIt
and VICTOR E. SAOUMAI
Department of Civil Engineering, University of Colorado, Boulder, Colorado, U.S.A. Abstract-New
models for the efficient simulation of discrete crack closure, and fracture mechanics based model for crack propagation, both under transient dynamic conditions, are developed. Following simple validation problems, those two models are integrated into an interactive graphics program which was used to analyse Koyna dam.
1. INTRODUCTION FOLLOWING over ten years of academic research in concrete fracture mechanics, it is rapidly becoming apparent that dams may well be the most suitable (if not only) application problem for this “maturing” discipline. While it was previously shown by the authors[l] that a fracture mechanics based analysis of cracked concrete dam yielded results which could have been obtained by current and simpler analysis procedures (summation of axial and flexural stresses), no alternative analytical/numerical methods appear to exist for the dynamic case. As such it is speculated by the authors that possibly the most important application of concrete fracture mechanics is in the dynamic analysis of dams. As this critical infrastructure has been receiving increased attention by the concrete fracture community, many problems remain to be addressed especially as they relate to dynamic analysis. In this paper we shall address the transient response of a concrete gravity dam subjected to ground motion with emphasis on its cracking behavior. The Implicit-Explicit time integration scheme of Hughes[2] is adopted for the discrete crack model in which two innovative approaches are introduced: (a) an efficient, yet accurate model for crack closure; and (b) a modified, damaged based, mixed mode crack propagation model for concrete. All those implementations were made on the program SICRAP (Simulation of Crack Propagation) developed by the authors[3,4,5]. Discrete crack, while providing highly “realistic” model for the few and “long” cracks typically encountered in dams, require special modeling when they are subjected to crack closure. While this problem seldom occurs in static analyses, it becomes of paramount importance in dynamic analysis, particularly if discrete cracks are used. It should be noted that this problem is not unique to fracture propagation in concrete gravity dams, but also to joint opening in arch dams. This problem has been partially addressed by models in which either fracture mechanics, or soil, fluid/structure interaction (but not yet both) play the dominant role. Droz[6] has developed a very innovative (and promising) approach to model smeared cracks in the dynamic analysis of dams. Skrikerud and his co-workers in Zurich[7,8] have analysed the Koyna dam using discrete cracks, along with strength approach. Dowling[9, lo], has developed a model for joint opening/closure in arch dams. Finally, Graves and Derucher[l l] used an improved model of the so called interface smeared crack model (ISCM) of Kuo[lZ] in which they model the stress free surface of the crack by eliminating the normal strains to the crack interface. While crack closure problems using primarily smeared cracks have recently been addressed within the context of dynamic analysis of dams, elastostatic analysis of contact problems have also been the subject of much research. In most, if not all, of those studies discrete discontinuities are used, and as discussed by Kikuchi and Oden[ 131and Taylor[ 141the main approaches fall into either one of the following categories: (a) Lagrange Multiplier methods; (b) Penalty methods. Unfortunately, in those methods, just as with joint elements (such as the ones of Goodman[l5] and Ghaboussi[ 161) the analysis requires numerous stiffness matrices decomposition.
tGraduate ResearchAssistant. SAssociate Professor.
588
MOHAME~
L. AYARI and VICTOR E. SAOUMA
Another major problem addressed in this paper is the condition for mixed mode crack propagation in concrete when subjected to dynamic excitation. While concrete response to impact load has received much attention lately[17-191, the fracture response of concrete subjected to random excitation has been addressed only by few researchers. One prominent model is the one of Chappuis[ZO] in which an elasto-plastic damage model which accounts for stress history has been developed and the other one is by Briihwilfer and Wittman[lO] which addresses the effects of preloading on G,, While this model is not strictly a fracture model, it addresses the conditions under which a crack under random excitation would propagate. 2. CO~ACT-I~PA~
MODEL
FO~~LA~O~
Following a formulation based on the one introduced by Kikuchi and Oden[l3] we denote the crack tip as c, and the two surfaces as S’, and S’. Furthermore we define a local coordinate system 0: (c,, 0,) &, t?,) centered at the crack tip, and defined in terms of the undeformed crack geometry (Fig. 1). Oz is taken to be normal to the plane containing the crack, while 8i and 0, are respectively, tangent and normal to the crack surface (@, is oriented along the inward normal of S’, or outward normal of S2). Furthermore, we designate as (I,, Iz, I,) an orthonormal basis associated with 8. We now consider two arbitrary points P(t);, 0:, 0:) and Q(@, S$, 6:) on S’ and S*, respectively, such that under deformation they are in direct contact, The kinematic conditions for incipient contact between P and Q are: Bj+uf
=e;+z$
j=
1,2
e: + V: 2 e: + v:.
(la)
(lb)
where the subject j refers to the local coordinate system 0, and the superscript to the crack face. v is the transformed displacement field given by: v = ru
(2)
where u is the displacement field in the global coordinate system, and r is the appropriate transformation matrix. Noting that the two cracked surfaces can be analytically defined by a function Y defined as: ei,=.!Y~(e:,ei,). The non~netration
j=l,2
(3)
condition of eq. (lb) can be rewritten as: Y’(Bi, 6:) + v: 2 Y2(#, e:> + v:.
Deformed
configuration
Fig. I. Local coordinate system 0.
(4)
Concrete gravity dams
589
Furthermore, for a pointwise contact condition, it can be shown that if we eliminate the first two coordinates ((I:, 0:) of Q from eq. (lb) by rewriting both Y2 and v$ as a first order Taylor series expansion with respect to both 8, and e2 and neglecting higher order terms we obtain (derivation can be found in[5,21]): nrf(v’ - v2) + (9’ - Y2)IBJ < 0
(5)
where n is a unit normal inward vector to Y2 at P. Moreover the normal stress at the interface must be compressive: a,, = San < 0 =o
if there is contact
(6)
otherwise
(7)
where (r is the stress tensor. Combining eq. (7) with (5) we obtain: Au,-B
@>
where: (9)
and Au, = nr(vi - v’). Finally the equations governing the nonfrictional lem) can be summarized as follows: 2 cTij,j
+fi
=
gij=
P
contact-impact
UO)
problem (Signorini’s prob-
3
bji
in G
(11)
Ui= fi/ on FV aijnj = ti on F,
(12)
Au, - 9GO (T,G 0 o,(Au, - Y) = 0.
on Tc
(13)
3. TEST PROBLEMS New algorithms to solve the previously presented contact~impact problem were developed (details of the numerical implementation can be found in[S, 211) and implemented into an existing finite element code written by the authors for the simulation of discrete crack growth[3]. In order to assess the accuracy of the algorithms two test problems (shown in Fig. 2), were analysed and results compared with known exact analytical solutions. Those problems are: 1. Two overlapping cantilevered beams. 2. Unidimensional dynamic Impact Problem. 3.1. Two overlapping ~un~i~e~ered beams In this first example we consider two cantilevered beams of equal rigidity EI and length L such that under zero load the lower surface of the first beam is in contact with the top surface of the other. Assuming frictionless surfaces, a uniform downward distributed load 4 acting on the top
MOHAMED
590
L. AYARI and VICTOR E. SAOUMA
-1.0
1
Fig. 2. (a) Two overlapping cantilevered beams. (b) Impact of a spring mass system against a rigid wall.
beam, can be decomposed into two loadings, a symmetric and an antisymmetric. The symmetric case is made up of two distributed loads of magnitude q/2 acting in opposite direction on the top fibers and the bottom fibers of the upper and lower beams, respectively. Due to symmetry there is no deflection under such a load case. In the second antisymmetric load case we have two distributed loads of magnitude q/2 acting along the same direction on the top and bottom of the two beams. Thus the net deflection would be the same as the one of a single beam under q/2, and is given by y =-&s’(-s2+4s
-6)
where
s = 2 L’
Assuming: Elastic modulus E = 3000 ksi, length L = 1000 in., width B = 10 in., height H = 40 in., and a uniform load of 0.02 kips/in, Fig. 3 shows the lower beam deflection in contrast with the close form solution, Table 1 shows displacement and force convergence norms, while Table 2 shows the end point deflection error. The final deformed mesh is shown in Fig. 4.
17.5
Exact solution -*--*..- Initial solution --Iteration I
r.:
gggos, 3
-t
Length x (in.) EE +31
Fig. 3. Beam deflection (two overlapping cantilevered beams problem).
Concrete gravity dams
591
Fig. 4. Deformed mesh (two overlapping cantilevered beams problem).
3.2. Dynamic impact Consider the uniaxial spring mass system shown in Fig. 2. The mass is initially displaced by x0 and then released. The mass will oscillate freely between x0 and xl were a rigid surface is present. The equation of motion describing the free vibration before impact is written as: x(t)=x,cosot
(14)
where w = &. Denoting by ti the time at which the mass strikes the rigid plane, it can be easily shown that the post-impact response is given by: X(t)
= XI COS fB(t
-
Ti) -
%
sin o(t - Zi)
(15)
where u, is the velocity at impact. For a mass m = 1.O, stiffness k = 1.O, x0 = 1.O, x, = -0.5 y = i and fi = $ (central difference algorithm) the response is compared to the exact solution in Fig. 5 for three different ratios of At/T. The algorithm is found to be very accurate and does not alter the overall energy balance of the system. The post-impact velocities and accelerations are captured correctly. A systematic derivation of the stability and accuracy of the algorithm is presented elsewhere[5].
4. FAILURE MODEL FOR CONCRETE CRACKS UNDER DYNAMIC EXCITATION As mentioned earlier most work on dynamic crack propagation has centered on impact load at very high frequencies. Furthermore those are typically simple impact loads. Under earthquake loading, a crack would typically be subjected to a complex histogram of excitation at frequencies well below those encountered under impact conditions. With the exception of the work of Chappuis in Zurich[20] and of Briihwiller in Lausanne[l9], the authors are not aware of any other related work. As such a new mixed mode quasistatic crack propagation model for concrete is presented.
Table 1. Norm convergence (two overlapping cantilevered beams) Iter. No.
1
2
3
Displacement Norm (%)
100.0
16.2
3.3
Force Norm (%)
100.0
33.9
14.1
4
5
6
7
8
9
10
15
0.70
0.16
0.045
0.017
0.009
0.006
0.004
0.001
9.1
7.8
2.4
18.0
12.0
12.2
11.6
10.4
Table 2. End beam deflection error Iteration No Error (%)
1
2
3
4
5
6
20.37
4.171
0.874
0.188
0.036
0.001
MOHAMED
592
? E 0
L. AYARI and VICTOR E. SAOUMA
0.75
0.50
% E 8
0.25
G i5
0
- 0.25
-0.75
I 2.50
0
Time
Fig. 5. Mass-spring
I 10.0
I 7.50
I 5.00
I 12.5
I!
(z.1
system impact response (y = f and /I = a).
In developing this new empirical model, the following assumptions were made: 1. Linear elastic fracture mechanics model is applicable. 2. A modification of the generalized maximum circumferential tensile stress model Saouma and Ayari[22] is adopted. 3. Excursions beyond the failure surface do not cause immediate crack extensions, they do however cause a localized “damage”. 4. Subcritical crack growth can take place as a result of damage accumulation. This is conceptually analogous to fatigue crack propagation. 5. Damage, as determined by the generalized maximum circumferential tensile stress theory is not localized at a discrete direction 13,)but rather spread around 0, according to a gaussian distribution. 6. Damage is greatly amplified when critical crack growth is to take place. 7. Crack extension takes place: (a) when average normalized damage exceeds unity; and (b) in the direction of maximum accumulated damage. On the basis of the preceding assumptions,
the following semi-empirical
model is adopted:
where: l
~9(0,0,) is an empirical distribution function which controls the damage localization around the angle e,(r). The following function is adopted.
d(e, e,) = =o
1 - t exp(1 - 5”)
if
5<$
if
<>$
(17)
where 5 = 10 - &(r)l/I(/, p is a decay factor which controls damage distribution, the “spread angle” of damage. l
The instantaneous angle of crack extension for the maximum circumferential theory is given by:
and II/ is
tensile stress
(18)
Concrete gravity dams
593
t,, is the starting time since last crack nucleation or extension. K,, is the quasistatic fracture toughness (and not the dynamic fracture toughness). l T is the integration variable. l t is the current time station. l r] is an experimental weight factor analogous to n in Paris Law (da/dN = c(AK)“). l k is a dynamic magnification factor. 0 K,(r), K,,(r) are stress intensity factors at time r. l
l
From this proposed model it can be seen that it is the normalized damage accumulation which controls the crack propagation. 5. CASE STUDY OF SEISMIC CRACKING 5.1. Koyna dam
While Peninsular India has long been considered as non-seismic, frequent tremors in 1962 lead to the installation of seismographs in the Koyna dam. Five years later a major earthquake shook the region and had its epicenter 2 km west and 13 km north of the dam. The depth focus was 12 km. Significant damage were sustained by the dam[23] in which horizontal cracks were observed in both the upstream and downstream faces of the taller nonoverllow monolith. As a result of this earthquake on an instrumented dam, Koyna became the “classical” analysis problem for the seismic analysis of dams. The literature is rich in case studies Koyna, and we shall only mention the early work by Chopra[23] and the recent state of the art paper by Hal1[25] on dynamic analysis of dams. To the best of the authors knowledge the only fracture analysis of Koyna (at least in terms of discrete crack propagation) was undertaken by Skrikerud[7,8] while Chapuis, Rebora and Zimmerman[26] used a fracture mechanics approach in a hybrid smeared-discrete crack model in the analysis of Pine Flat Dam. 5.2. Material parameters
In this Analysis two discretizations of the dam were adopted, a coarse and a fine one (Fig. 6) A summary of the material constants used is found in Table 3. We shall note that with exception to K,, all those values are assumed to be the dynamic ones for concrete.
Fig. 6. Meshes used in the analysis.
Table 3. Koyna analysis, material constants Elastic modulus Poisson’s ratio Tensile strength Compressive strength Mass density Quasistatic fracture toughness
31. 0.2 2.4 24 2640. 1.5
GPa GPa GPa
MOHAMED
L. AYARI and VICTOR E. SAOUMA
Mod! I r, 10.324s
Mode 5 r, ~0.041s
Mode 4 T, -0.062
Mode 3 r, = 0.092s
s
Fig. 7. Mode shapes of Koyna dam.
5.3 Modal analysis At first a modal analysis on the coarse mesh was undertaken using the subspace iteration technique[27]. The first five modes were determined and found to compare very well with the ones obtained by Chopra[24]. The individual mode shapes along with their respective natural periods of vibration are shown in Fig. 7. 5.4. Nonlinear transient dynamic analysis Following this preliminary modal analysis, a non linear analysis was undertaken using the finer mesh. At this point it is essential to indicate that the aim of this study was not to perform a complete comprehensive exact analysis of the dam (which would have had to include soil/structure/fluid interaction) but rather to focus solely on the dam’s fracture response in a hypothetical simulation. Thus this analysis is only a mean of applying models previously developed. The dynamic analysis was thus performed with the modfied version of SICRAP[3] in which dynamic contact-impact models are implemented. A Rayleigh damping model, in which the damping matrix is linearly expressed in terms of the mass and stiffness matrices as: C = aM + bK was used and the parameters a and b were selected
I
I
-0.40 -0.30 -0.20 -0.100
0 0.10 0.20
I
0.30 0.w
.., 0
I’ I,
I
..,
2.3
Time (ssc)
Fig. 8. Horizontal and vertical ground accelerations.
.,
5.0
Time
I,,
7.e
(sex)
,
,
,
I
I6
10.0
I
Ii?.,
-l.no
I.00
I
zoo
I Time
+.oo
I (5)
9.00
I
I 6.00 I.00
I
em
I 9 D
-4.00
-MO0 I
I
.oo
I 2.00
I s.00 !LwJ
t
Time
(I
,‘I , 4.00
Fig. 10. Cracking sequence.
Fig. 9. Horizontal and vertical crest displacements (linear elastic case).
3.00
I
i 1.00
I I.00
I 0.00
9
MOHAMED L. AYARI and VICTOR E. SAOUMA
596
,_._
12.5
10.0
E
too -
7.50 7.50 -
g
5
5.w
8
5.00-
5 f ', %
2.50 2.500
z $
0
-2.50 -500
0
I 1.00
I
I
2.w
5.00
I 4.00
Time
I
I
!mo
6.00
I 7.00
I
I
t
I
I
2.00
5.00
4.00
5.00
I
8.00
-2.501 0
I.00
f 5t
Time
I s.GQ
I 1.00
(I
t s1
Fig. 1I. Crest displacement for the nonlinear case.
such that a 5% and 25% damping
were assumed for the first and fifth fundamental modes, respectively. This led to u = 0.724 and b = 0.00323. The Newmark’s parameters y and /3 used in this analysis were 0.25 and 0.55, respectively. Furthermore all elements were assumed to be implicit as we were seeking the long term response of the structure. The new algorithm produced a solution within ten to fifteen iterations per time step when a crack closes. For this initial analysis, and short of experimental data, failure model parameters had to be assumed. As it is well known that K,,/lu, is much higher than unity {up to IO5 according to{ZSJ, an arbitrary value of 60 was selected with q = 1.1 and p = 0.016. The time step adopted was At = 3.0 ms. Such a small value was selected in order to maintain a good accuracy in the post-impact response taking place at the crack(s) interface(s) and to ensure a moderately small numerical dissipation (less than 10%) after 5 s. 5.5. A~~Iysis rest&s The analysis was performed on an Alliant FX/8 “mini-super” computer. The program was optimized by the compiler to take simultaneous advantage of both concurrency and vectorization. Following each crack nu~leation/extension, results were sent to an Apollo workstation through an Ethernet connection, in which remeshing was performed and results sent back to the Alliant through a restart capability, seven seconds of ground motion excitation were analysed in over 30 h of CPU.
350 -
300-
250 -
2
2Qo-
2
I50 -
Time from
Fig.
12.
cmck
initiation
(s I
No~aljzed stress intensity factors.
Concrete gravity dams
597
This analysis is best summarized through the following figures. Figure 10 for the sequence of crack nucleation/extension, Fig. 11 shows both horizontal and vertical crest displacements. Finally Fig. 12 shows the normalized stress intensity factors of the two cracks. It is interesting to note the alternation in crack closure between the upstream and downstream crack. Also the frequency and duration of crack closure should be studied further in order to come up with an experimental set up that would bring more information about concrete crack behavior in those ranges of excitation. 6. CONCLUSIONS Dynamic fracture of dams is probably the most important application of concrete fracture mechanics. An important amount of work has to be directed to understand crack behavior under seismic excitations, such as modeling the dynamic failure surface and the mixed mode direction of crack extension. In this contribution new crack closure/impact model is validated for the potential applications of joint closure in arch dams and crack closure in gravity dams under both static and dynamic loadings. This model is shown to be both efficient, inexpensive and reliable. Also a new semi-empirical model is proposed for transient crack propagation in concrete, however it requires experimental validation. All these models have been integrated into a single code, which performed a fracture analysis of the Koyna dam. While the analysis yielded realistic results, it must be emphasized that no claim is made regarding its correctness. Such a cautious statement is based on our inability to compare results with similar analysis and the arbitrariness of the fracture parameters selected. Acknowledgemenr-Panial
support from the Electric Power Research Institute for this research is gratefully acknowledged.
REFERENCES [I] V. E. Saouma, M. L. Ayari and H. Boggs, Fracture mechanics of roller compacted concrete gravity dams. Invited Paper, Fracture of Dams Session, Proc. fnr. Conf. on Fracture of Concrere rmd Rocks, Houston (June, 1987). [2] T. J. Hughes and W. K. Liu, Implicit+xplicit finite elements in transient analysis: implementation and numerical examples,J. uppl. Mech. 45, (June, 1978). [3] V. E. Saouma, M. L. Ayari, SICRAP User’s and Theory Manual. University of Colorado, Boulder (1988). [4] V. E. Saouma and M. L. Ayari, Simulation of CRAck Propagation (SICRAP-PC) User’s ond Technical Manual. University of Colorado, Dept. of Civil Eng. (1988). [5] M. L. Ayari. Fracture mechanics of concrete gravity dams. Doctoral Thesis, University of Colorado, Boulder (August, 1988). [6] P. Droz, Modele numerique du comportement non-lineaire d’ouvrages massifs en beton artne. Doctoral Thesis, No. 682, Swiss Polytechnic Institute, Lausanne (1987). (71 P. Skrikerud, Mcdelle und Berechnungsverfahren fur das Rissverhalten von Unarmierten Betonbauten unter Erbebenbeanspruchung (Fracture Model for Concrete Dam under Seismic Loading). Institute for Static & Construction. Swiss Polvtechnic Institute. Zurich IJune. 1983). PI P. Skrikerud, Discrete crack modelling for dynamically loaded, unreinforced concrete structures. Earthquake Engng Structural Dynamics, 14, 297-3 I5 (1986). [91J. F. Hall and M. J. Dowling, Response of jointed arches to earthquake excitation. Earthquake Engng Swucrural Dynamics
13, 779-798
(1985).
WI M. J. Dowling, Nonlinear seismic analysis of arch dams. Doctoral Thesis, Rep. No. EERl 87-03, California Institute of Technology, Pasadena, California (1987). IllI R. H. Graves and K. N. Derucher, Interface smeared crack model analysis of concrete dams in earthquakes. J. engng Me& 113 (November, 1987). WI J. Kuo, Joint opening nonlinear mechanism. Interface smeared crack model, Rep. No. UBC/EERc 82/10 Earthquake Engng Res. Ctr., Berkeley, California (1982). N. KiKuchi and J. T. Gden. Contact problems in elastostatics, in Finite E/emenrs, Special Problems in Solid Mechanics, 1131 (Edited by J. T. Gden and G. F. Carey), Vol. 5. Prentice Hall, New Jersey (1984). (141 R. L. Taylor, Methods for solution of dynamic contact problems. NCEL Contract Report, CR-86.002, Naval Civil Engineering Laboratory, California (December, 1985). (151 R. E. Goodman and C. St John, Numerical methods in geotechnical engineering (Edited by Desai and Christian), pp. 148-175. McGraw-Hill, New York. iI61J. Ghaboussi, E. L. Wilson and J. Henberg, Finite element for rock joints and interfaces. J. Soil Mech. Foundations Div. 833-848 (October, 1973). R. John and S. P. Shah, Effect of high strength and rate of loading on fracture parameters of concrete. Proc. Inr. Conf. 1171 on Fracfure of Concrete and Rocks, Houston (June, 1987).
598
MOHAMED
L. AYARI and VICTOR E. SAOUMA
[ 181 N. Banthia, S. Mindness and A. Bentur, Energy balance in instrumented impact tests on plain concrete beams. Proc. ht. ConJ on Fracture of Concrete and Rocks, Houston (June, 1987). [19] E. Briihwiller and F. H. Wittman, Failure of dam concrete subjected to seismic loading conditions. Engng Fracture Mech. 35, 117-125 (1990). [20] P. Chappuis, Modelisation non-lineaire du comportement du beton sous des sollicitations dynamiques. Doctoral Thesis, Zurich, 1987. Institute for Static ULConstruction, Swiss Polytechnic Institute, Zurich (June, 1983). [21] M. L. Ayari and V. E. Saouma, Static and dynamic contact/impact problems using fictitious interface forces. Internal Report, University of Colorado. [22] V. E. Saouma, M. L. Ayari and D. A. Leavell, Mixed mode crack propagation in homogeneous anisotropic solids. Engng Fracture Mech. 27, 171-184 (1987). [23] N. Pal, Seismic cracking of concrete gravity dams. J. Structural Division, Proc. American Society of Civil Engineers, 102, No. ST9 (1976). [24] A. K. Chopra, Earthquake response of concrete gravity dams. Report No. EERC 70-1, University of California,
Berkely (January, 1970). [25] J. F. Hall, The dynamic and earthquake behavior of concrete dams: review of experimental behaviour and observational evidence. Soil Dynamics Earthquake Engng No. 2 7, 57-121 (1988). [26] J. Chapuis, P. Rebora and Th. Zimmermann, Numerical approach of crack propagation analysis in gravity dams durina earthauakes. Trans. 15th. ICOLE. Vol. 2. DD. 451-473 (1985). [27] K.-J. Bathe, ‘Finite Element Procedures in Enginee&g Analysis.’ Prentice-Hall, New Jersey (1982). [28] M. F. Kanninen and C. H. Popelar, Aduanced Fracture Mechanics. Oxford University Press, Oxford (1985). [29] D. R. J. Owen and E. Hinton, Finite Elements in Plasticity, Theory and Practice. Pineridge Press, Swansea (1980). (Received for publication 16 November 1988)
Vol. 35, No. l/2/3, Pp. W-598, 1990
Printed in Great Britain.
0013-7944/90 $3.00 + 0.00 Perpmon Press pk.
A FRACTURE MECHANICS BASED SEISMIC ANALYSIS OF CONCRETE GRAVITY DAMS USING DISCRETE CRACKS MOHAMED
L. AYARIt
and VICTOR E. SAOUMAI
Department of Civil Engineering, University of Colorado, Boulder, Colorado, U.S.A. Abstract-New
models for the efficient simulation of discrete crack closure, and fracture mechanics based model for crack propagation, both under transient dynamic conditions, are developed. Following simple validation problems, those two models are integrated into an interactive graphics program which was used to analyse Koyna dam.
1. INTRODUCTION FOLLOWING over ten years of academic research in concrete fracture mechanics, it is rapidly becoming apparent that dams may well be the most suitable (if not only) application problem for this “maturing” discipline. While it was previously shown by the authors[l] that a fracture mechanics based analysis of cracked concrete dam yielded results which could have been obtained by current and simpler analysis procedures (summation of axial and flexural stresses), no alternative analytical/numerical methods appear to exist for the dynamic case. As such it is speculated by the authors that possibly the most important application of concrete fracture mechanics is in the dynamic analysis of dams. As this critical infrastructure has been receiving increased attention by the concrete fracture community, many problems remain to be addressed especially as they relate to dynamic analysis. In this paper we shall address the transient response of a concrete gravity dam subjected to ground motion with emphasis on its cracking behavior. The Implicit-Explicit time integration scheme of Hughes[2] is adopted for the discrete crack model in which two innovative approaches are introduced: (a) an efficient, yet accurate model for crack closure; and (b) a modified, damaged based, mixed mode crack propagation model for concrete. All those implementations were made on the program SICRAP (Simulation of Crack Propagation) developed by the authors[3,4,5]. Discrete crack, while providing highly “realistic” model for the few and “long” cracks typically encountered in dams, require special modeling when they are subjected to crack closure. While this problem seldom occurs in static analyses, it becomes of paramount importance in dynamic analysis, particularly if discrete cracks are used. It should be noted that this problem is not unique to fracture propagation in concrete gravity dams, but also to joint opening in arch dams. This problem has been partially addressed by models in which either fracture mechanics, or soil, fluid/structure interaction (but not yet both) play the dominant role. Droz[6] has developed a very innovative (and promising) approach to model smeared cracks in the dynamic analysis of dams. Skrikerud and his co-workers in Zurich[7,8] have analysed the Koyna dam using discrete cracks, along with strength approach. Dowling[9, lo], has developed a model for joint opening/closure in arch dams. Finally, Graves and Derucher[l l] used an improved model of the so called interface smeared crack model (ISCM) of Kuo[lZ] in which they model the stress free surface of the crack by eliminating the normal strains to the crack interface. While crack closure problems using primarily smeared cracks have recently been addressed within the context of dynamic analysis of dams, elastostatic analysis of contact problems have also been the subject of much research. In most, if not all, of those studies discrete discontinuities are used, and as discussed by Kikuchi and Oden[ 131and Taylor[ 141the main approaches fall into either one of the following categories: (a) Lagrange Multiplier methods; (b) Penalty methods. Unfortunately, in those methods, just as with joint elements (such as the ones of Goodman[l5] and Ghaboussi[ 161) the analysis requires numerous stiffness matrices decomposition.
tGraduate ResearchAssistant. SAssociate Professor.
588
MOHAME~
L. AYARI and VICTOR E. SAOUMA
Another major problem addressed in this paper is the condition for mixed mode crack propagation in concrete when subjected to dynamic excitation. While concrete response to impact load has received much attention lately[17-191, the fracture response of concrete subjected to random excitation has been addressed only by few researchers. One prominent model is the one of Chappuis[ZO] in which an elasto-plastic damage model which accounts for stress history has been developed and the other one is by Briihwilfer and Wittman[lO] which addresses the effects of preloading on G,, While this model is not strictly a fracture model, it addresses the conditions under which a crack under random excitation would propagate. 2. CO~ACT-I~PA~
MODEL
FO~~LA~O~
Following a formulation based on the one introduced by Kikuchi and Oden[l3] we denote the crack tip as c, and the two surfaces as S’, and S’. Furthermore we define a local coordinate system 0: (c,, 0,) &, t?,) centered at the crack tip, and defined in terms of the undeformed crack geometry (Fig. 1). Oz is taken to be normal to the plane containing the crack, while 8i and 0, are respectively, tangent and normal to the crack surface (@, is oriented along the inward normal of S’, or outward normal of S2). Furthermore, we designate as (I,, Iz, I,) an orthonormal basis associated with 8. We now consider two arbitrary points P(t);, 0:, 0:) and Q(@, S$, 6:) on S’ and S*, respectively, such that under deformation they are in direct contact, The kinematic conditions for incipient contact between P and Q are: Bj+uf
=e;+z$
j=
1,2
e: + V: 2 e: + v:.
(la)
(lb)
where the subject j refers to the local coordinate system 0, and the superscript to the crack face. v is the transformed displacement field given by: v = ru
(2)
where u is the displacement field in the global coordinate system, and r is the appropriate transformation matrix. Noting that the two cracked surfaces can be analytically defined by a function Y defined as: ei,=.!Y~(e:,ei,). The non~netration
j=l,2
(3)
condition of eq. (lb) can be rewritten as: Y’(Bi, 6:) + v: 2 Y2(#, e:> + v:.
Deformed
configuration
Fig. I. Local coordinate system 0.
(4)
Concrete gravity dams
589
Furthermore, for a pointwise contact condition, it can be shown that if we eliminate the first two coordinates ((I:, 0:) of Q from eq. (lb) by rewriting both Y2 and v$ as a first order Taylor series expansion with respect to both 8, and e2 and neglecting higher order terms we obtain (derivation can be found in[5,21]): nrf(v’ - v2) + (9’ - Y2)IBJ < 0
(5)
where n is a unit normal inward vector to Y2 at P. Moreover the normal stress at the interface must be compressive: a,, = San < 0 =o
if there is contact
(6)
otherwise
(7)
where (r is the stress tensor. Combining eq. (7) with (5) we obtain: Au,-B
@>
where: (9)
and Au, = nr(vi - v’). Finally the equations governing the nonfrictional lem) can be summarized as follows: 2 cTij,j
+fi
=
gij=
P
contact-impact
UO)
problem (Signorini’s prob-
3
bji
in G
(11)
Ui= fi/ on FV aijnj = ti on F,
(12)
Au, - 9GO (T,G 0 o,(Au, - Y) = 0.
on Tc
(13)
3. TEST PROBLEMS New algorithms to solve the previously presented contact~impact problem were developed (details of the numerical implementation can be found in[S, 211) and implemented into an existing finite element code written by the authors for the simulation of discrete crack growth[3]. In order to assess the accuracy of the algorithms two test problems (shown in Fig. 2), were analysed and results compared with known exact analytical solutions. Those problems are: 1. Two overlapping cantilevered beams. 2. Unidimensional dynamic Impact Problem. 3.1. Two overlapping ~un~i~e~ered beams In this first example we consider two cantilevered beams of equal rigidity EI and length L such that under zero load the lower surface of the first beam is in contact with the top surface of the other. Assuming frictionless surfaces, a uniform downward distributed load 4 acting on the top
MOHAMED
590
L. AYARI and VICTOR E. SAOUMA
-1.0
1
Fig. 2. (a) Two overlapping cantilevered beams. (b) Impact of a spring mass system against a rigid wall.
beam, can be decomposed into two loadings, a symmetric and an antisymmetric. The symmetric case is made up of two distributed loads of magnitude q/2 acting in opposite direction on the top fibers and the bottom fibers of the upper and lower beams, respectively. Due to symmetry there is no deflection under such a load case. In the second antisymmetric load case we have two distributed loads of magnitude q/2 acting along the same direction on the top and bottom of the two beams. Thus the net deflection would be the same as the one of a single beam under q/2, and is given by y =-&s’(-s2+4s
-6)
where
s = 2 L’
Assuming: Elastic modulus E = 3000 ksi, length L = 1000 in., width B = 10 in., height H = 40 in., and a uniform load of 0.02 kips/in, Fig. 3 shows the lower beam deflection in contrast with the close form solution, Table 1 shows displacement and force convergence norms, while Table 2 shows the end point deflection error. The final deformed mesh is shown in Fig. 4.
17.5
Exact solution -*--*..- Initial solution --Iteration I
r.:
gggos, 3
-t
Length x (in.) EE +31
Fig. 3. Beam deflection (two overlapping cantilevered beams problem).
Concrete gravity dams
591
Fig. 4. Deformed mesh (two overlapping cantilevered beams problem).
3.2. Dynamic impact Consider the uniaxial spring mass system shown in Fig. 2. The mass is initially displaced by x0 and then released. The mass will oscillate freely between x0 and xl were a rigid surface is present. The equation of motion describing the free vibration before impact is written as: x(t)=x,cosot
(14)
where w = &. Denoting by ti the time at which the mass strikes the rigid plane, it can be easily shown that the post-impact response is given by: X(t)
= XI COS fB(t
-
Ti) -
%
sin o(t - Zi)
(15)
where u, is the velocity at impact. For a mass m = 1.O, stiffness k = 1.O, x0 = 1.O, x, = -0.5 y = i and fi = $ (central difference algorithm) the response is compared to the exact solution in Fig. 5 for three different ratios of At/T. The algorithm is found to be very accurate and does not alter the overall energy balance of the system. The post-impact velocities and accelerations are captured correctly. A systematic derivation of the stability and accuracy of the algorithm is presented elsewhere[5].
4. FAILURE MODEL FOR CONCRETE CRACKS UNDER DYNAMIC EXCITATION As mentioned earlier most work on dynamic crack propagation has centered on impact load at very high frequencies. Furthermore those are typically simple impact loads. Under earthquake loading, a crack would typically be subjected to a complex histogram of excitation at frequencies well below those encountered under impact conditions. With the exception of the work of Chappuis in Zurich[20] and of Briihwiller in Lausanne[l9], the authors are not aware of any other related work. As such a new mixed mode quasistatic crack propagation model for concrete is presented.
Table 1. Norm convergence (two overlapping cantilevered beams) Iter. No.
1
2
3
Displacement Norm (%)
100.0
16.2
3.3
Force Norm (%)
100.0
33.9
14.1
4
5
6
7
8
9
10
15
0.70
0.16
0.045
0.017
0.009
0.006
0.004
0.001
9.1
7.8
2.4
18.0
12.0
12.2
11.6
10.4
Table 2. End beam deflection error Iteration No Error (%)
1
2
3
4
5
6
20.37
4.171
0.874
0.188
0.036
0.001
MOHAMED
592
? E 0
L. AYARI and VICTOR E. SAOUMA
0.75
0.50
% E 8
0.25
G i5
0
- 0.25
-0.75
I 2.50
0
Time
Fig. 5. Mass-spring
I 10.0
I 7.50
I 5.00
I 12.5
I!
(z.1
system impact response (y = f and /I = a).
In developing this new empirical model, the following assumptions were made: 1. Linear elastic fracture mechanics model is applicable. 2. A modification of the generalized maximum circumferential tensile stress model Saouma and Ayari[22] is adopted. 3. Excursions beyond the failure surface do not cause immediate crack extensions, they do however cause a localized “damage”. 4. Subcritical crack growth can take place as a result of damage accumulation. This is conceptually analogous to fatigue crack propagation. 5. Damage, as determined by the generalized maximum circumferential tensile stress theory is not localized at a discrete direction 13,)but rather spread around 0, according to a gaussian distribution. 6. Damage is greatly amplified when critical crack growth is to take place. 7. Crack extension takes place: (a) when average normalized damage exceeds unity; and (b) in the direction of maximum accumulated damage. On the basis of the preceding assumptions,
the following semi-empirical
model is adopted:
where: l
~9(0,0,) is an empirical distribution function which controls the damage localization around the angle e,(r). The following function is adopted.
d(e, e,) = =o
1 - t exp(1 - 5”)
if
5<$
if
<>$
(17)
where 5 = 10 - &(r)l/I(/, p is a decay factor which controls damage distribution, the “spread angle” of damage. l
The instantaneous angle of crack extension for the maximum circumferential theory is given by:
and II/ is
tensile stress
(18)
Concrete gravity dams
593
t,, is the starting time since last crack nucleation or extension. K,, is the quasistatic fracture toughness (and not the dynamic fracture toughness). l T is the integration variable. l t is the current time station. l r] is an experimental weight factor analogous to n in Paris Law (da/dN = c(AK)“). l k is a dynamic magnification factor. 0 K,(r), K,,(r) are stress intensity factors at time r. l
l
From this proposed model it can be seen that it is the normalized damage accumulation which controls the crack propagation. 5. CASE STUDY OF SEISMIC CRACKING 5.1. Koyna dam
While Peninsular India has long been considered as non-seismic, frequent tremors in 1962 lead to the installation of seismographs in the Koyna dam. Five years later a major earthquake shook the region and had its epicenter 2 km west and 13 km north of the dam. The depth focus was 12 km. Significant damage were sustained by the dam[23] in which horizontal cracks were observed in both the upstream and downstream faces of the taller nonoverllow monolith. As a result of this earthquake on an instrumented dam, Koyna became the “classical” analysis problem for the seismic analysis of dams. The literature is rich in case studies Koyna, and we shall only mention the early work by Chopra[23] and the recent state of the art paper by Hal1[25] on dynamic analysis of dams. To the best of the authors knowledge the only fracture analysis of Koyna (at least in terms of discrete crack propagation) was undertaken by Skrikerud[7,8] while Chapuis, Rebora and Zimmerman[26] used a fracture mechanics approach in a hybrid smeared-discrete crack model in the analysis of Pine Flat Dam. 5.2. Material parameters
In this Analysis two discretizations of the dam were adopted, a coarse and a fine one (Fig. 6) A summary of the material constants used is found in Table 3. We shall note that with exception to K,, all those values are assumed to be the dynamic ones for concrete.
Fig. 6. Meshes used in the analysis.
Table 3. Koyna analysis, material constants Elastic modulus Poisson’s ratio Tensile strength Compressive strength Mass density Quasistatic fracture toughness
31. 0.2 2.4 24 2640. 1.5
GPa GPa GPa
MOHAMED
L. AYARI and VICTOR E. SAOUMA
Mod! I r, 10.324s
Mode 5 r, ~0.041s
Mode 4 T, -0.062
Mode 3 r, = 0.092s
s
Fig. 7. Mode shapes of Koyna dam.
5.3 Modal analysis At first a modal analysis on the coarse mesh was undertaken using the subspace iteration technique[27]. The first five modes were determined and found to compare very well with the ones obtained by Chopra[24]. The individual mode shapes along with their respective natural periods of vibration are shown in Fig. 7. 5.4. Nonlinear transient dynamic analysis Following this preliminary modal analysis, a non linear analysis was undertaken using the finer mesh. At this point it is essential to indicate that the aim of this study was not to perform a complete comprehensive exact analysis of the dam (which would have had to include soil/structure/fluid interaction) but rather to focus solely on the dam’s fracture response in a hypothetical simulation. Thus this analysis is only a mean of applying models previously developed. The dynamic analysis was thus performed with the modfied version of SICRAP[3] in which dynamic contact-impact models are implemented. A Rayleigh damping model, in which the damping matrix is linearly expressed in terms of the mass and stiffness matrices as: C = aM + bK was used and the parameters a and b were selected
I
I
-0.40 -0.30 -0.20 -0.100
0 0.10 0.20
I
0.30 0.w
.., 0
I’ I,
I
..,
2.3
Time (ssc)
Fig. 8. Horizontal and vertical ground accelerations.
.,
5.0
Time
I,,
7.e
(sex)
,
,
,
I
I6
10.0
I
Ii?.,
-l.no
I.00
I
zoo
I Time
+.oo
I (5)
9.00
I
I 6.00 I.00
I
em
I 9 D
-4.00
-MO0 I
I
.oo
I 2.00
I s.00 !LwJ
t
Time
(I
,‘I , 4.00
Fig. 10. Cracking sequence.
Fig. 9. Horizontal and vertical crest displacements (linear elastic case).
3.00
I
i 1.00
I I.00
I 0.00
9
MOHAMED L. AYARI and VICTOR E. SAOUMA
596
,_._
12.5
10.0
E
too -
7.50 7.50 -
g
5
5.w
8
5.00-
5 f ', %
2.50 2.500
z $
0
-2.50 -500
0
I 1.00
I
I
2.w
5.00
I 4.00
Time
I
I
!mo
6.00
I 7.00
I
I
t
I
I
2.00
5.00
4.00
5.00
I
8.00
-2.501 0
I.00
f 5t
Time
I s.GQ
I 1.00
(I
t s1
Fig. 1I. Crest displacement for the nonlinear case.
such that a 5% and 25% damping
were assumed for the first and fifth fundamental modes, respectively. This led to u = 0.724 and b = 0.00323. The Newmark’s parameters y and /3 used in this analysis were 0.25 and 0.55, respectively. Furthermore all elements were assumed to be implicit as we were seeking the long term response of the structure. The new algorithm produced a solution within ten to fifteen iterations per time step when a crack closes. For this initial analysis, and short of experimental data, failure model parameters had to be assumed. As it is well known that K,,/lu, is much higher than unity {up to IO5 according to{ZSJ, an arbitrary value of 60 was selected with q = 1.1 and p = 0.016. The time step adopted was At = 3.0 ms. Such a small value was selected in order to maintain a good accuracy in the post-impact response taking place at the crack(s) interface(s) and to ensure a moderately small numerical dissipation (less than 10%) after 5 s. 5.5. A~~Iysis rest&s The analysis was performed on an Alliant FX/8 “mini-super” computer. The program was optimized by the compiler to take simultaneous advantage of both concurrency and vectorization. Following each crack nu~leation/extension, results were sent to an Apollo workstation through an Ethernet connection, in which remeshing was performed and results sent back to the Alliant through a restart capability, seven seconds of ground motion excitation were analysed in over 30 h of CPU.
350 -
300-
250 -
2
2Qo-
2
I50 -
Time from
Fig.
12.
cmck
initiation
(s I
No~aljzed stress intensity factors.
Concrete gravity dams
597
This analysis is best summarized through the following figures. Figure 10 for the sequence of crack nucleation/extension, Fig. 11 shows both horizontal and vertical crest displacements. Finally Fig. 12 shows the normalized stress intensity factors of the two cracks. It is interesting to note the alternation in crack closure between the upstream and downstream crack. Also the frequency and duration of crack closure should be studied further in order to come up with an experimental set up that would bring more information about concrete crack behavior in those ranges of excitation. 6. CONCLUSIONS Dynamic fracture of dams is probably the most important application of concrete fracture mechanics. An important amount of work has to be directed to understand crack behavior under seismic excitations, such as modeling the dynamic failure surface and the mixed mode direction of crack extension. In this contribution new crack closure/impact model is validated for the potential applications of joint closure in arch dams and crack closure in gravity dams under both static and dynamic loadings. This model is shown to be both efficient, inexpensive and reliable. Also a new semi-empirical model is proposed for transient crack propagation in concrete, however it requires experimental validation. All these models have been integrated into a single code, which performed a fracture analysis of the Koyna dam. While the analysis yielded realistic results, it must be emphasized that no claim is made regarding its correctness. Such a cautious statement is based on our inability to compare results with similar analysis and the arbitrariness of the fracture parameters selected. Acknowledgemenr-Panial
support from the Electric Power Research Institute for this research is gratefully acknowledged.
REFERENCES [I] V. E. Saouma, M. L. Ayari and H. Boggs, Fracture mechanics of roller compacted concrete gravity dams. Invited Paper, Fracture of Dams Session, Proc. fnr. Conf. on Fracture of Concrere rmd Rocks, Houston (June, 1987). [2] T. J. Hughes and W. K. Liu, Implicit+xplicit finite elements in transient analysis: implementation and numerical examples,J. uppl. Mech. 45, (June, 1978). [3] V. E. Saouma, M. L. Ayari, SICRAP User’s and Theory Manual. University of Colorado, Boulder (1988). [4] V. E. Saouma and M. L. Ayari, Simulation of CRAck Propagation (SICRAP-PC) User’s ond Technical Manual. University of Colorado, Dept. of Civil Eng. (1988). [5] M. L. Ayari. Fracture mechanics of concrete gravity dams. Doctoral Thesis, University of Colorado, Boulder (August, 1988). [6] P. Droz, Modele numerique du comportement non-lineaire d’ouvrages massifs en beton artne. Doctoral Thesis, No. 682, Swiss Polytechnic Institute, Lausanne (1987). (71 P. Skrikerud, Mcdelle und Berechnungsverfahren fur das Rissverhalten von Unarmierten Betonbauten unter Erbebenbeanspruchung (Fracture Model for Concrete Dam under Seismic Loading). Institute for Static & Construction. Swiss Polvtechnic Institute. Zurich IJune. 1983). PI P. Skrikerud, Discrete crack modelling for dynamically loaded, unreinforced concrete structures. Earthquake Engng Structural Dynamics, 14, 297-3 I5 (1986). [91J. F. Hall and M. J. Dowling, Response of jointed arches to earthquake excitation. Earthquake Engng Swucrural Dynamics
13, 779-798
(1985).
WI M. J. Dowling, Nonlinear seismic analysis of arch dams. Doctoral Thesis, Rep. No. EERl 87-03, California Institute of Technology, Pasadena, California (1987). IllI R. H. Graves and K. N. Derucher, Interface smeared crack model analysis of concrete dams in earthquakes. J. engng Me& 113 (November, 1987). WI J. Kuo, Joint opening nonlinear mechanism. Interface smeared crack model, Rep. No. UBC/EERc 82/10 Earthquake Engng Res. Ctr., Berkeley, California (1982). N. KiKuchi and J. T. Gden. Contact problems in elastostatics, in Finite E/emenrs, Special Problems in Solid Mechanics, 1131 (Edited by J. T. Gden and G. F. Carey), Vol. 5. Prentice Hall, New Jersey (1984). (141 R. L. Taylor, Methods for solution of dynamic contact problems. NCEL Contract Report, CR-86.002, Naval Civil Engineering Laboratory, California (December, 1985). (151 R. E. Goodman and C. St John, Numerical methods in geotechnical engineering (Edited by Desai and Christian), pp. 148-175. McGraw-Hill, New York. iI61J. Ghaboussi, E. L. Wilson and J. Henberg, Finite element for rock joints and interfaces. J. Soil Mech. Foundations Div. 833-848 (October, 1973). R. John and S. P. Shah, Effect of high strength and rate of loading on fracture parameters of concrete. Proc. Inr. Conf. 1171 on Fracfure of Concrete and Rocks, Houston (June, 1987).
598
MOHAMED
L. AYARI and VICTOR E. SAOUMA
[ 181 N. Banthia, S. Mindness and A. Bentur, Energy balance in instrumented impact tests on plain concrete beams. Proc. ht. ConJ on Fracture of Concrete and Rocks, Houston (June, 1987). [19] E. Briihwiller and F. H. Wittman, Failure of dam concrete subjected to seismic loading conditions. Engng Fracture Mech. 35, 117-125 (1990). [20] P. Chappuis, Modelisation non-lineaire du comportement du beton sous des sollicitations dynamiques. Doctoral Thesis, Zurich, 1987. Institute for Static ULConstruction, Swiss Polytechnic Institute, Zurich (June, 1983). [21] M. L. Ayari and V. E. Saouma, Static and dynamic contact/impact problems using fictitious interface forces. Internal Report, University of Colorado. [22] V. E. Saouma, M. L. Ayari and D. A. Leavell, Mixed mode crack propagation in homogeneous anisotropic solids. Engng Fracture Mech. 27, 171-184 (1987). [23] N. Pal, Seismic cracking of concrete gravity dams. J. Structural Division, Proc. American Society of Civil Engineers, 102, No. ST9 (1976). [24] A. K. Chopra, Earthquake response of concrete gravity dams. Report No. EERC 70-1, University of California,
Berkely (January, 1970). [25] J. F. Hall, The dynamic and earthquake behavior of concrete dams: review of experimental behaviour and observational evidence. Soil Dynamics Earthquake Engng No. 2 7, 57-121 (1988). [26] J. Chapuis, P. Rebora and Th. Zimmermann, Numerical approach of crack propagation analysis in gravity dams durina earthauakes. Trans. 15th. ICOLE. Vol. 2. DD. 451-473 (1985). [27] K.-J. Bathe, ‘Finite Element Procedures in Enginee&g Analysis.’ Prentice-Hall, New Jersey (1982). [28] M. F. Kanninen and C. H. Popelar, Aduanced Fracture Mechanics. Oxford University Press, Oxford (1985). [29] D. R. J. Owen and E. Hinton, Finite Elements in Plasticity, Theory and Practice. Pineridge Press, Swansea (1980). (Received for publication 16 November 1988)