Microstructurally-based modelling of intergranular creep fracture using grain elements

MECHANICS OF MATERIALS ELSEVIER

Mechanics of Materials 26 (1997) 109-126

Microstructurally-based modelling of intergranular creep fracture using grain elements Patrick Onck *, Erik van der Giessen Laboratory for Engineering Mechanics, Delft University of Technology, P.O. Box 5033, 2600 GA Delft, The Netherlands

Received 7 March 1996; received in revised form 10 April 1997

Abstract

A new numerical method is proposed to simulate intergranular creep fracture in large polycrystalline aggregates. The method utilizes so-called grain elements to represent the polycrystal. These grain elements take care of the average elastic and creep deformation of individual grains. Grain boundary processes, like cavitation and sliding, are accounted for by grain boundary elements connecting the grains. Results are compared with full-field finite element calculations. The method is demonstrated to capture the essential features of creep fracture, like creep constrained cavitation and the interlinkage of microcracks. Also the: performance in polycrystals with random variations in microstructure, in terms of grain shape, is shown to be reasonably well. For the size of the unit-cell considered, a factor of around 600 is gained in computer time as compared with the full-field calculations. © 1997 Elsevier Science Ltd.

1. I n t r o d u c t i o n At elevated temperatures and relatively low stresses, failure in polycrystalline metals after years of stationary loading is often intergranular. Observation of the crack surface shows a more or less uniform distribution of dimples, which are the remnants of cavities. These cavities nucleate and grow most rapidly on facets that are normal to the maximum principal stress direction until they coalesce to form grain boundary microcracks (Argon, 1982; Cocks and Ashby, 1982). Final intergranular failure

* Corresponding author. Laboratory for Engineering Mechanics, Melkweg 2, 2628 CD Delft, The Netherlands. Tel.: +31-152782703; fax: +31-15-2782150; e-mail: p.r.onck@wbmt. tudelft.nl.

occurs when these microcracks link-up with the macroscopic crack. In the modelling of creep fracture, three different approaches can be identified. Firstly, at the macroscopic level, classical fracture mechanics approaches are extended to time-dependent behavior in order to model creep crack growth. The idea is to identify the appropriate load parameter (e.g. K I, J, C * ), so that measured crack growth rates can be transferred from the laboratory to the component in service. In this approach it is assumed that the asymptotic near-tip stress and strain fields are unaffected by damage. Secondly, in continuum damage models, cavitation is incorporated in an average, smeared-out manner by means of a damage parameter, which is either micromechanism-based (Hutchinson, 1983; Tvergaard, 1986) or entirely phenomenological (Hayhurst and Leckie, 1984; Hayhurst et al., 1984). Thirdly, micromechanical models have been developed which

0167-6636/97/$17.00 @ 1997 Elsevier Science Ltd. All rights reserved. PII S0167-6636(97),30020-3

110

P. Onck, E. van der Giessen / Mechanics of Materials 26 (1997) 109-126

are directly based on the various physical mechanisms involved. This paper will be concerned exclusively with the third approach. In the micromechanical modelling of creep rupture three length scales can be identified: (i) the scale of cavities (microscopic level), (ii) the scale of grains (mesoscopic level) and (iii) the macroscopic or continuum level. The objective of this kind of modelling is to prescribe the material behavior on a certain level by means of the relevant micromechanical processes of the underlying level. During the last two decades a scale transition has been made from the considerations of a single cavity (Needleman and Rice, 1980; Sham and Needleman, 1983) to that of failure of a polycrystalline aggregate comprising a number of grains (Van der Giessen and Tvergaard, 1994a,b,c). However, the latter computations, using full-field approximations inside the grains and over the grain boundaries, have reached the limits of present computer power. Therefore, in order to reach the macroscopic level, a final micromechanical modelling step is needed. This recently led to the introduction of a latticetype of model: the so-called Delaunay network modelling technique (Van der Burg and Van der Giessen, 1994a). This approach presumes free grain boundary sliding and is based on representing the polycrystalline aggregate by a network of truss-like elements, connecting the centres of neighboring grains. The constitution of such a Delaunay element accounts for the creep response of part of the two grains as well as for cavitation damage on their mutual grain boundary facet. Although this model has been shown capable of showing the main phenomena (Van der Burg and Van der Giessen, 1994a,b), an inherent drawback of this lattice model is that it needs calibration with respect to the detailed analysis (Van der Giessen and Tvergaard, 1994a) and that the dependence of the grain deformations on the applied stress state is not captured well. To circumvent these difficulties a new approximate model is proposed. The idea is to represent the polycrystalline material by an aggregate of so-called grain elements that describe the average elastic and creep deformations of each individual grain. Cavitation and free grain boundary sliding are accounted for by grain boundary elements, connecting the grains. The constitution of the grain elements and the

grain boundary elements is corrected for the strain rate enhancement caused by free grain boundary sliding. With this approach we are geared to describe macroscopic intergranular creep fracture, based on the underlying microstructure and physical mechanisms. This paper is concerned with the formulation of the grain element method for planar problems as well as with the verification by comparison with detailed finite element calculations for a multi-grain unit cell model studied previously by Van der Giessen and Tvergaard (1994a,b,c). After the formulation of that problem, we recapitulate the constitutive equations governing the relevant physical processes in Section 2. Then, we proceed to define the grain elements and grain boundary elements, as well as their numerical implementation (Section 3). Finally, in Section 4, grain element predictions will be shown for a practically relevant range of material conditions and applied stress states, which will be compared with detailed full-field calculations. 2. Problem formulation and constitutive models

A planar polycrystalline material is considered which is represented by hexagonal grains. It is taken to be subjected to a stationary macroscopic stress state specified by the principal true stresses ~1 and ~2 in plane strain (see Fig. 1). In line with the experimental observation that cavitation occurs predominantly on the transverse facets, we take the maximal principal stress direction along the xZ-axis. The material is taken to exhibit a certain periodicity so that only a unit-cell needs to be analyzed. Let 2A 0 and 2B 0 be the dimensions of the unit-cell in the x j and x 2 directions, respectively, and let the cell contain ml × m 2 grains (the cell shown in Fig. 1 is characterized by (m l, m 2) = (6, 5)). It is also assumed that the unit-cell itself exhibits reflection symmetry in the x 1 and x 2 directions, so that only one quadrant needs consideration. Due to symmetries, the four faces of the quarter unit-cell will remain straight and aligned with the coordinate axes, and will not support any shear stress. The analysis is carried in a finite strain, convected coordinate formulation of the governing equations. The contravariant stress components of the Kirchhoff tensor with respect to the deformed base vectors are

P. Onck, E. van der Giessen / Mechanics of Materials 26 (1997) 109-126

111

~2

I ......

m2 grains

£2 Y m 1 grains Fig. 1. Planar polycrystal anit-cell model comprising m I X m 2 grains (in this example, (m j, m 2) = (6, 5). Only a quarter of the unit cell is analyzed.

denoted by r ;j and related to the components of the Cauchy stress tensor:, o-~, by r u = G ~ o-ij with G and g the determinants of the metric tensors Gij and g u in the current and reference configuration, respectively. The strain components conjugate to r s-/ are the covariant components of the Lagrangian strain I tensor, which are expressed as rl~~ = 7(u~,j + uj, i + u~uk, j) in terms of the displacement components u i (with 0,i denoting covariant differentiation in the reference frame). The material inside the grains is assumed to be homogeneous and to deform by power law creep in addition to elasticity. Accordingly, the total Lagrangian strain rate 4-/,.y is written as the sum of an elastic part, i?E, and a creep part, //c, given by #c =

3 sij

(la)

Here, t o and o-0 are reference strain rate and stress quantities and n is the creep exponent. The elastic strains remain small, such that the effective Mises stress ~ = (3sijsiJ/2) 1/2 and the stress deviator s iJ = r i j - GiJr~/3 can be specified directly in terms of the Kirchhoff stresses, by neglect of the elastic volume changes. The elastic response is governed by vii

the hypoelastic relationship r

= "p~jkt~E '~ki, in terms

of the Jaumann stress-rate vii 'r = ,j.ij + ( Gil¢.g jt + Gjk,l.il) i?kl

(2)

and the elastic moduli RiJ kl =

E 1 + v

×{2(GikGJt+GilGJk)+

1 - 2 v GuGkl

'

112

P. Onck, E. van der Giessen / Mechanics of Materials 26 (1997) 109-126

so that the constitutive relations for the elastic-creeping grain material can be written as vu

=

(3)

Experimental evidence shows that the grain boundary cavities responsible for intergranular creep rupture in many metals and alloys remain much smaller than the grain size and that facets often contain a distribution of many cavities. Therefore, it is useful to introduce the density N of cavities along a grain boundary. The nucleation of cavities at the grain boundaries is observed to occur continuously during the creep process, such that the number of cavities is approximately proportional to the overall effective strain. Dyson (1983) explained this by suggesting that cavities nucleate by growing from athermal decohesions which arise when slip bands intersect with grain boundary particles. Acknowledging the fact that cavity nucleation is a complex process and still not understood well, it is incorporated here by the evolution equation (Van der Giessen and Tvergaard, 1990) = G(

co/Z0 ) " ec'

(4)

where tr, is the local normal stress at the grain boundary and F, is a material constant (X0 is a mere normalization factor). The cavities are often observed to have a spherical-caps shape and are specified by their radius a, half-spacing b and cavity tip angle ~b (see Fig. 2a). The cavity volume can b e expressed as V = 4/3~ra3h(~b), where h is the cavity shape parameter, d e f i n e d by h(~b) = [(1 + cos ~b)-i - ½ c o s ~b]/sin ~b. The cavities grow by the diffusional flow of atoms from their surface into the grain boundary layer. If we assume surface diffusion to be much faster than grain boundary diffusion, the cavities retain their quasi-equilibrium spherical-caps shape. Not only do cavities grow by diffusion, they also grow by creep of the surrounding grain material. It seems that creep, additionally, accelerates the plating of atoms from the cavity surface into the grain boundary layer. Needleman and Rice (1980) and Sham and Needleman (1983) analyzed an axisymmetric cavity model, similar to Fig. 2a, and investigated the combined influence of creep and diffusion

2a

,@

---C> ~-

2b

~1

(b)

J (a) Fig. 2. (a) Geometryof a quasi-equilibriumspherical-capsshaped cavity with radius a, half spacing b and cavity tip angle 6. (b) Representation of cavities on a grain boundary in the planar polycrystalmodel. by means of detailed finite element calculations. They summarized their numerical results into an expression for the volumetric cavity growth rate V. This expression was later modified by Tvergaard (1984), resulting in

= v, + v2.

(5)

Here, 91 is the contribution of diffusion, specified through cr,- (1-f)~ 47r'@" I n ( l / f ) - ( 1 / 2 ) ( 3 - f ) ( 1

(6)

-f)

with f=max

-~

,

a+l.5~

'

(7)

and the contribution of creep, I;"2 is given by _-+-2~reeCagh(O) o~n ~

+

.

92=

,

for _ + - -

> 1;

ere 2~r4Ca3h(0)[a, +/3,]" - m %

for am < 1. o,e

(8) The coupling between diffusive and creep contributions to void growth enters in Eq. (7) through the stress and temperature dependent length scale

L= [ ~ / ~ 7 ] '/3,

(9)

P. Onck, E. van der Giessen / Mechanics of Materials 26 (1997) 109-126

introduced by Needleman and Rice (1980). For small values of a/L (say, a/L < 0 . 1 ) cavity growth is dominated by diffusion, while for larger values of alL creep growth becomes more and more important. The range of validity of Eq. (5) was specified to be alL < 10. Furthermore in Eq. (6), .~ is the grain boundary diffusion parameter and % is the sintering stress which will be neglected. The constants a , and ft, in Eq. (8) are given by or. = 3/2n and /3, = (n - 1 ) ( n + 0.4319)/n:". The effective stress tr:, the mean stress trm and the normal stress or, are considered to be stresses remote from each cavity (on the scale of individual cavities). However, on the size scale of grains, with cavities being small compared to the grain facet site, the stresses must be seen as local near the grain boundary. In the axisymmetfic cavity model of Fig. 2a, it is tacitly assumed that the grain boundary facet is built up by identical cylinders, thus assuming a hexagonal close packing of identical cavities. In the planar polycrystal model, the out-of-plane dimension is not incorporated explicitly, and the state of cavitation is represented simply by the in-plane grain boundary parameters a and b, as shown in Fig. 2b. Nevertheless, it should be noted that the cavities are still treated to have the spherical-caps shape.

Conservation of mass requires that any cavity growth needs to be accompanied by a separation of the adjacent grains. This can be expressed in terms of an average separation 6c = V/('n'b2), which is the volume of a cavity smeared out over the grain boundary area associated with the cavity (Rice, 1981). This concept now plays a central role in the polycrystal models considered here for materials where grain facets contain a distribution of many cavities that are small compared to the grain size. Then, following Rice (1981) and Tvergaard (1984), the actual discrete distribution of cavities can be 'smeared out' over the facet, so that the discrete distribution of cavities with radius a and half-spacing b is replaced by a continuous distribution along each facet. Thus, the average separation 6c also varies continuously along the grain boundary facet, and can alternatively be regarded as a grain boundary layer of thickness 6c (see Fig. 3). The rate of separation is given by

/~e

I) ,rrb 2

2V b 7rb 2 b '

(10)

thus expressing that the grains can separate by an increasing cavity volume V and by a decreasing cavity spacing b. The latter is related to the nucleation of new cavities and, to a lesser extent, to finite strain effects associated with deformations of the ligaments between the cavities. Recalling the fact that the cavity density N at a part of the grain boundary with current surface area d A (dA o initially) is defined as dA 1 N = 'n'b-----~ dA ° ,

1 ~5

C

Fig. 3. (a) Individual cavities with varying radius on a grain boundary. (b) 'Smeared..out' model in terms of a layer with thickness 6c.

113

(11)

it can be deduced that b 1 ~ = "~(~, + ~n)

12Q 2N'

(12)

where ~ and ~n are the principal logarithmic strain rates in the plane of the grain boundary, and where /¢ follows from the nucleation law (Eq. (4)). According to Ashby (1972), grain boundaries at elevated temperatures can be modelled as thin layers that slide in a Newtonian viscous manner, so that the

P. Onck, E. van der G iessen / Mechanics of Materials 26 (1997) 109-126

114

relative sliding velocity 3's of adjacent grains due to a shear stress r in the grain boundary is given by T

= w--,

a / b (see Fig. 2). From the definition of the cavity volume V it follows that the cavity growth rate is given by

(13)

r/B

& = (//(4zra2h(qJ)),

where w is the thickness of the boundary. In his model, the grain boundary viscosity r/B is related to the grain boundary diffusivity. Ashby observed that, even though irregularities in the grain boundary, such as ledges or second-phase particles, can increase the boundary viscosity, it still is much lower than the creep resistance of the grain material in situations of sufficiently high temperatures and low stresses. Therefore, following earlier work (e.g. Van der Giessen and Tvergaard, 1991, 1994a,b,c; Rice, 1981; Hsia et al., 1991) the grain boundary viscosity is neglected here all together, 71B/W = 0, SO that the grains can slide freely against each other. It should be noted that the expression for the volumetric growth rate ~', given in Eq. (5), has been derived for cavities on non-sliding grain boundaries. Diffusional growth accompanied with grain boundary sliding gives rise to non-equilibrium cavity shapes (Argon, 1982). However, when the siding rates are low, Eq. (5) may be used as a good approximation for cavity growth on sliding grain boundaries (Chen, 1983). Substitution of Eqs. (4), (5) and (12) in Eq. (10) yields the separation rate t~c as a function of the present state of cavitation and stress. An appropriate way to quantify the damage state is the parameter

(14)

and together with Eq. (12) this describes the damage evolution. When the ratio a / b approaches unity, coalescence of cavities occurs. However, justified by experimental observations, coalescence may occur earlier by ductile tearing or cleavage of the ligament between the cavities (Cocks and Ashby, 1982). As in Van der Giessen and Tvergaard (1994b,c), we assume coalescence to take place at a / b = 0.7.

3. The grain element method

3.1. Numerical formulation In this new numerical method to simulate intergranular creep fracture, each grain is represented by a special finite element, termed a grain element, that accounts for elasticity and creep inside the grain in an average manner. The polycrystalline material is built up directly by connecting these grain elements by way of grain boundary elements which account for grain boundary sliding and cavitation (see Fig. 4). In the approach chosen here, each grain element can be regarded as a six-noded super element, which consists of two quadrilaterals. The quadrilaterals are

Grain,..element /

,

/

?

Grain boundary element Fig. 4. Grain elements representationof a polycrystallineaggregate. Each grain element consists of two quadrilaterals, built up of four constant strain triangles. The grains are connectedby grain boundaryelements.

115 P. Onck, E. van der Giessen/ Mechanics of Materials 26 (1997) 109-126 Vu = 3Sijr /(20"~) together with the Eqs. (la), (lb) each built up of four constant strain triangles that lie along the diagonals of the quadrilateral. The choice for the two quadrilaterals is prompted by the fact that other representations can lead to locking of the grain element mesh. For a more detailed discussion of this numerical artifact, the reader is referred to Appendix A. The governing equations for the creeping grains are basically the same as in earlier publications by Van der Giessen and Tvergaard (1991, 1994a,b,c). The averaging over each grain is carried out by invoking numerical: integration in the virtual work expression. Here, a linear incremental version of the principle of virtual work for a time step A t is used in the form (cf. Van der Giessen and Tvergaard, 1991)

and (3), the modified effective creep rate can be written as a function of ~/kV Finally, upon using this modified effective creep rate in Eqs. (la) and (lb) and combining this with Eq. (2) in Eq. (3), the constitutive equations for the elastic-creeping grain material can be written in the form (Tvergaard, 1984)

with 1 L ijkl = e ijkl - 7 ( T i k G jl + 7"JkG il -F TJlG ik "~- TilGJk)

At fv(+iJrn,s + 7"ikl'~{krUj,i)dV MiJ=

=AtfsT"i~uidS~-[SsTi~uidS-SvTiJ~ijdV] , (15) with 6u~ and 6"Oij being the virtual variations of the displacements and associated strains, respectively. The bracketed terna is an equilibrium correction, included to prevent the solution from drifting off the equilibrium path. The stress rates "i"~j in Eq. (15) are to be eliminated by way of the constitutive equations. Use of the original expressions (Eq. (3)) would lead to rather poor :numerical stability and therefore a forward gradient method is used as suggested by Tvergaard (1984) similar to the one proposed by Peirce et al. (1984). Recognizing the fact that instabilities are mainly caused by the highly nonlinear, power law dependence of the effective strain rate ~c according to Eqs. (la) and (lb), the latter is written as a linear interpolation between the values at time t and t + At: ~c=(1-o)~CO)+O~C('+at',

0~[0,11.

(17)

~.ij = LiJk,i.lkt + ;r~j,

[~

MiJMkl

1+ ~

h

h sit --,

~;C(t) '

oi~c 15j=Oh--At,

0~

~'~J= - - - 1+ ~ Mij' 3 h=---

E

21+v"

Using the constitutive Eq. (17) in Eq. (15) with a value of 0 = 0.9 has been found to improve the numerical stability considerably (Tvergaard, 1984). To account for the thickening of the grain boundary layer due to the nucleation and growth of cavities as well as for the relative sliding of adjacent grains, special-purpose grain boundary elements are used. Each grain boundary facet is modelled by a single grain boundary element, as is shown in Fig. 4. The definition of such a grain boundary element follows the same lines as in Van der Giessen and Tvergaard (1991). The kinematical relations for the grain boundary elements are set up at each step during the incremental procedure, based on the current configuration (see Fig. 5). The grain boundary element is associated with the segment ac'cd, where a' and c' are the orthogonal projections of a and c on the opposite grain facet. The grain boundary nodes p and q are midpoints of the lines connecting a node

(16)

Next, a Taylor series expansion is used to estimate the strain rate at time t + A t, a~c

~c(,+~,) = ~co) +

~at.

0~ Then, by substituting the above in Eq. (16) and using Fig. 5. Definitionof a grain boundaryelement.

116

P. Onck, E. van der Giessen / Mechanics of Materials 26 (1997) 109-126

and its projection and, thus, define the current midplane of the grain boundary facet. A local coordinate ~'~ [0, 1] and a local set of base vectors (e I, e 2) are connected to p q . We can now define a separation rate 8 and shearing rate ~ at each position ff within the grain boundary element by ~ = (b'm- b'n).e 2,

~ = (/)m- b'n).e I,

(18)

where v., and v~ are the associated velocity vectors at the edges of the adjacent grains. Upon discretization, v m and v n are written as linear interpolations between the values in the grain nodes a, b and c, d, respectively. The separation and shearing rates, 6 and ~/, are related to the grain boundary stress state through the constitutive relations, Eqs. (10) and (13). However, the nature of these (nonlinear) viscous-like constitutive equations is different from the elastic-creep constitutive equations for the grains. For numerical convenience, we therefore introduce a fictitious layer of linear elastic springs against opening and shearing, so that the constitutive behaviour of the grain boundary layer is described by the following visco-elastictype relationships: 6-n = k n ( 8 - 6~),

(19)

(20) where k n and k s are the elastic stiffnesses of the normal and shearing spring layers, respectively. The deviations 8 - 8¢ and ~ - 4/s are kept small by using large values for the stiffnesses. Furthermore, in this analysis, the free grain boundary sliding condition is incorporated by using a very small value for ~ ) B / w . The average separation rate 8c is a function of o"n, cre and o-m through the volumetric growth rate relation (Eq. (5)). The normal stress o-n is defined directly in the grain boundary elements; but for and o-m we take the average over the two triangles located on either side of the grain boundary element (see Fig. 5). To increase the stable step size, a forward gradient method is also used for the integration of Eqs. (19) and (20). It is noted that 6c is a function of o"n, and o-m, but we only incorporate the dependence on o-n as this was found to be the main cause for

numerical instabilities (Tvergaard, 1984; Van der Giessen and Tvergaard, 1991). Then, the modified constitutive relations can be written as

withkn* = k n

~'=k[(~-4L),

l+kn0-~At

,

(21)

withk[ =k s

1

,

(22) where again 0 = 0.9 is used in the calculations. The finite element equations for the grain boundary elements are also based on the incremental virtual work equation, A t f7 ( 6-n 88 + "i'ST) dV

= at fj'au, dS + [fsT'au, dS - f ( na8 +

(23)

but with the internal work rate computed in an updated Lagrangean manner based on the current deformed configuration (volume V). This is computationally convenient in view of the grain boundary sliding. With the definition of the grain boundary elements, the volume integral reduces to a line integral, which is evaluated for each element by Gaussian quadrature using two sampling points. Finally, for the unit-cell analysis outlined in Fig. 1, the boundary conditions read ti 1 = 0 ,

T 2=0,

along x l = 0

/~/1 = UI,

T 2 = 0,

along x 1 = A 0

~i2 = 0 ,

T 1=0,

along x 2 = 0

ti 2 = L?n ,

T 1 = 0,

along x 2 = B 0 ,

in terms of the velocity components ti i and the nominal traction components T i in the reference configuration. The uniform velocities 01 and 0 u are determined by a special Rayleigh-Ritz technique (Needleman and Tvergaard, 1984), such that the

P. Onck, E. van der Giessen / Mechanics of Materials 26 (1997) 109-126

macroscopic average true stresses "~l and fined as

~l

=

"~2, de-

/

/

/ ~ - / H - d , 6 X ~ ~ \

117 \

\

//J// /

1 £aOT z Ix' =ao ':Ix 2,

\

Z2 = L [a°r2l~2=B<, d x 1 AJo

'

remain constant. Here A and B are the dimensions of the quarter unit-cell. Details found in (Van der Burg and Van der Giessen, The macroscopic logarithmic strain rates are by

current can be 1994b). defined

\

\ \

\

/ /

(a)

c

(b)

2e

- -

e

1.94

1.67

SEI = l)x/A , 1~2 = (lillB.

1.41 1.14

3.2. A reference to detailed finite element analyses The polycrystal material model, described in Section 2, has been the subject of several investigations by Van der Giessen ~md Tvergaard (1994a,b,c). These have used full-field finite element approximations of the goveming equations, using many 'crossed triangle' quadrilaterals to discretize each grain and similarly many grain boundary elements for each grain boundary facet. A typical finite element discretization is shown in Fig. 6a. With this approach a rather clear picture has been obtained concerning the interaction between the various micromechanical processes responsible fi)r final intergranular fracture. A major difference between such a full-field approach and the grain element representation (see Fig. 4), is that the latter can only account for limited variations in stress inside the grains, as well as variations in cavitation and stress over the grain boundary facets. It will be recalled that the idea behind the grain element approach is to deprive of a lot of details in order to make it feasible to analyze polycrystalline aggregates of macroscopic dimensions. To do so with some confidence in the accuracy of the approach, a reference has to be made to detailed finite element analyses. Therefore, in Section 4, results of the two approaches will be compared in order to verify the accuracy of the grain element method. Detailed finite element calculations will not only serve as a verification tool, they will also be used as a guideline to incorporate physical phenomena that

0.88 0.61 0.35 0.08 Fig. 6. (a) Example of a detailed finite element mesh used in the analyses by Van der Giessen and Tvergaard (1994b,c). The grain boundary elements are not visible. (b) Contour plot of the effective stress o"e, normalized by the applied effective stress Xe, obtained by using the mesh of (a) for a problem with free grain boundary sliding in the absence of cavitation.

act on a smaller size scale. In particular, we shall do so in relation to free grain boundary sliding. In a power-law creeping polycrystalline aggregate with homogeneous grains and no cavitation, grain boundary sliding is accommodated by creep inside the grains, thus resulting in a highly nonuniform creep flow. This nonuniform flow appears to enhance the overall creep rate. Fig. 6b shows the normalized effective stress distribution, obtained with the detailed finite element mesh of Fig. 6a for the problem of a (m 1, m 2) = (2, 1) unit-cell (see Fig. 1) (the complete plot is produced with mirror images of the quarter cell results). Crossman and Ashby (1975) postulate that the overall behavior of a regular polycrystal in plane strain with free sliding still obeys a power-law of the form in Eq. (lb) for macroscopic strain rate /~ as a function of macroscopic effective stress ,~, but modified as /~e = ~0(f* ~ ) " ,

(24)

118

P. Onck, E. van der Giessen / Mechanics o f Materials 26 (1997) 109-126

with f * ( f * > 1) being the so-called stress enhancement factor. The strain rate enhancement depends on the stress exponent n and the geometry of the microstructure. Ghahremani (1980) and Hsia et al. (1991) investigated f * as a function of the creep exponent for an array of regular hexagonal grains in plane strain. For n = 5 they found that f " = 1.2. More recently, Onck and Van der Giessen (1995) also incorporated variations in microstructure, in terms of grain size and shape, into the analysis and found that the strain rate was further enhanced in case of random deviations from the regular hexagonal structure. They fitted their finite element results by the following expression for the stress enhancement factor: 3 f * (n, p) =f,:o + "~( n + 15)( p - p0) 2.

(25)

Here, p is a parameter that characterizes the geometrical randomness in microstructure, f,~0 is the stress enhancement factor of the regular hexagonal geometry for creep exponent n and P0 the parameter for the regular hexagonal geometry (Onck and Van der Giessen, 1995). It is clear that the grain elements cannot pick up the stress concentrations (Fig. 6b) and the associated highly non-uniform creep flow inside the grains that are caused by grain boundary sliding. In order to account for the resulting average enhancement in creep flow over each grain, we simply apply a stress enhancement factor f * in the power-law in Eq. (lb) to be used for that grain, similar to the macroscopic equivalent (Eq. (24)). Strictly speaking, the value of this enhancement factor should be a function of the grain's shape; but, the same macroscopic f * will be attributed to all grains in the aggregate. Thus, in simulating creep fracture in very large polycrystalline aggregates with different microstructures, we use Eq. (25) to determine the stress enhancement factor f * to be used in the constitution of the grain elements. Note that the enhanced creep flow in the grain material also influences the growth of cavities on the grain boundary. This is incorporated by using the free sliding enhanced effective strain rate in the creep part I)2 in Eq. (6), as well as in the diffusion length parameter L in Eq. (9) and in the nucleation law (Eq. (4)).

4. Verification

To show the accuracy of the grain element method we compare the results with detailed finite element calculations, performed by Van der Giessen and Tvergaard (1994a,b,c). The meshes used in these full-field calculations within each element correspond to the one shown in Fig. 6a. We confine attention to a subset of calculations from the above references, which capture the essential features of intergranular creep fracture. First we look at cavitation on an isolated grain boundary facet and the development of damage from an initial microcrack (Sections 4.1 and 4.2, respectively). These calculations are for the regular hexagonal microstructure, shown in Fig. 1. Finally, the effect of random variations in grain shape is briefly considered (Section 4.3). Close examination of the governing equations in Section 2 yields that the whole process of creep rupture can be described by non-dimensional parameters through the relation

a(o tg

k~,-b I'

) R I' N I ' E '

Xe

'

(26) where the subscript I refers to the initial state and 2R is the facet width. The time to failure te is normalized by the reference time t R = ~J(EEC), where X, is the applied effective stress and/~c the corresponding effective strain rate through Eq. (1 b). All calculations are carried out with ~ / E = 0.5 X 10 -3, so that the reference time is the same for all cases. By specifying a value for t R we automatically set the creep parameter ~0/tr0" in Eq. (lb). The number of independent parameters in Eq. (26) is further reduced by taking n = 5, v = 0.3 and ~/,= 75 °. The stress triaxiality Xm/,~e, corresponding to the applied stresses ,~ and X 2, is approximated in plane strain + -~) and .~ through the expressions "~m ~--" ~(-~2 ! = ½V/-31~2 -- "~1[ valid for pure creep. If cavitation is considered to occur on a grain boundary facet it will be assumed to be initially uniform and specified by (a/b) x = 0.01 and (b/R) I = 1. Note that this corresponds to an initial cavity density of N I = 1/(~'R2). Except when specified otherwise, the nucleation pro-

119

P. Onck, E. van der Giessen/ Mechanics of Materials 26 (1997) 109-126

cess is specified through F n / N ! = 100 and the scaling parameter ~0 in Eq. (4) is taken to be equal to "~e. It is assumed that nucleation stops when N = Nma~ = 100N r The initial value for the diffusion length parameter L is specified by substituting the applied effective str,:ss "~e and strain rate /~c in Eq. (9). Then, ( a / L ) I is adjusted through the grain boundary diffusion parameter ..~. Note that with these specifications the set of independent parameters in Eq. (26) reduced to ( a l L ) I a n d ~ m / / , ~ e .

o~ 0.6

0.4

0.3

0.2

.....

Grain

- -

Detailed

element

method

analysis

o.I o.1~

. . . .

r

0.0

(a)

. . . .

i

0.5

. . . .

t

1.0

,

,

1.5

~/~ ul R

4.1. Creep constrained cavitation a/rb

In this subsection we consider the case where only the central facet of the (6, 5) unit-cell is allowed to cavitate, while the surrounding material remains undamaged. When diffusion is much faster than creep, the average separation 8c due to cavitation on the central facet is constrained by creep of the surrounding grains (Dyson, 1976; Rice, 1981; Tvergaard, 1984). This situation corresponds to a relatively small value for the parameter a / L , which determines the relative contribution of creep and diffusion to cavity growth (see Section 2). The evolution of damage on the cavitating facet is monitored through the parameter a / b in the centre of the facet as well as at the triple point (denoted by C and T, respectively, in forthcoming figures). Note that this is only meaningfu]l for the detailed finite element results, since the cavitation described by the grain element method remains uniform over the facet due to symmetry of the unit cell. In the first case, we apply uniaxial tension along the x2-axis ( ~ m / . E e--0.577), while cavitation on the facet is initially uniform and specified by ( a / b ) I = 0.1 and ( b / R ) I = 0.1. No nucleation is accounted for, F n = 0, so that N I = Nmax -- 100/(TrR2). Fig. 7a shows the results when diffusion is much faster than creep, ( a / L ) l = 0.025. The damage evolution is seen to be rather uniform which suits the grain element method well. When the contribution of creep to the cavity growth rate increases ( ( a / L ) t = 0.1, see Fig. 7b), the damage is less uniform. However, the grain elements still give a good description of the damage development averaged over the facet. Next, we consider cases where continuous nucleation of cavities takes place ( F n = 100 NI), with an initial state of lower cavitation damage, ( a / b ) t =

j , ""

0.6

0.5

0.4

0.3

0.2 ment 0.1

- -

o.o

(b)

,

o

,

,

i

lo

. . . .

i

20

,

,

,

Detailed

i

30

,

,

t/tR

method

analysis

,

Fig. 7. The damage parameter a / b versus normalized time t / t R on an isolated cavitating facet in a (6, 5) unit-cell for uniaxial tension, Xm/"~e = 0.577, and no nucleation, F, ~ 0. (a) (a/L) z = 0.025; (b) (a/L) l = 0.1.

0.01, ( b / R ) 1 = 1 (corresponding to N 1 = 1/(TrR2)). Figs. 8a and b correspond to using ( a l L ) I = 0.025 and 0.1, respectively, and again uniaxial tension is applied. Comparison with Fig. 7 shows, as expected, that the introduction of nucleation with little initial damage results in a considerable increase in lifetime. Cavitation is seen to become much more non-uniform, however. For ( a / L ) l = 0.025, the grain elements still give a good description of the average cavitation damage. But, for the case where creep is much more significant, the damage development near the triple point is very much faster than elsewhere, and microcracking towards the center of the facet is much like a tearing process; this cannot be picked up well by the grain element description, which therefore somewhat overestimates the lifetime. Finally, in Fig. 9, we consider the same cases but subject to a higher level of triaxiality, ~ m / / . ~ e = 4, corresponding to the situation in front of a macroscopic plane strain crack tip. The lifetime decreases,

P. Onck, E. van der Giessen / Mechanics of Materials 26 (1997) 109-126

120

as compared with Fig. 8, since the higher level of triaxiality increases the applied mean stress "~m, while ~e remains fixed, so that the cavity growth process on average is accelerated. As a consequence, the creep constraint is more active than in Fig. 8, and the damage evolves in a more uniform manner. However, the grain element description still has some difficulty with picking-up the rapid acceleration of damage in the facet center near the end of the life time. In spite of this, it can be concluded that the grain elements do capture the dependence of the rupture process on the macroscopic stress state through the stress triaxiality ratio ,Ym/Ze.

--)~

a/b 07 06 0s o.4 0.3

02 ol 11o

. . . . . .

~,,,

1

(a)

k . . . . . . .

4

5

/

°'f

/

o.i/

4.2. Damage development from an initial microcrack

t,,,

3

07[ o5F

In this subsection a (6, 5) unit-cell is considered with an initial microcrack located in the center. The surrounding material is virtually undamaged initially but cavitation may develop along all facets in the course of the process. Nucleation is accounted for as

i ....

2

0.3

O. 1 [-/

; ...''""

.'"'"

__

(b)

o

5

lo

7 t/tR

i

/

I"

/ ~/

[' . ~ _ _ ~ - - ~ - -

15

f,,

6

%

---

Grain element method

- -

Detailed analysis

20 l/tR

Fig. 9. The damage parameter a / b vs. normalized time t / t R on an isolated cavitating facet in a (6, 5) unit-cell for a higher triaxiality, 2 m / ' S e = 4, and F, = 100 N I. (a) ( a l L ) l = 0.025; (b) ( a l L ) l = 0.1.

0.7

a/b 0.6 0.5 0.4 0,3 0.2 0.1 0.0

(a)

a/b

,

,

i .... 10

j , , 20

,

k , , 30

O'T/

i 40

J

~tR

"/'""

0.6 0.5

../""

.-'" o.1 [ . - " "

(b)

,

0.0 ",~'-'-'--~. ~ 0 20

j ~ , 40

~ 60

,, 80

..... ---,t,,, 100

Grain element method Detailed analysis ~, 120

t/tR

Fig. 8. The damage parameter a / b vs. normalized time t i t R on an isolated cavitating facet in a (6, 5) unit-cell for uniaxial tension, ~ m / 2 e = 0.577, and F, = 100 N 1. (a) ( a / L ) l = 0.025; (b) ( a l L ) I 0.1. =

before. Uniaxial tension is applied in the x2-direc tion (Em/,~e = 0.577) and ( a / L ) I is specified to be 0.025. Fig. 10 shows snapshots of the damage state at two instants during the lifetime in a quarter of the unit-cell. The time bars denote the elapsed time as a fraction of the ultimate time to failure tf. The state of damage is shown by plotting the value of a / b perpendicular to each facet and with the ordinate along the facet. The facets where microcracking has occurred due to cavity coalescence at a / b = 0.7 are highlighted by a darker gray level. Figs. 10a and b show the results for the detailed finite element calculations and Figs. 10c and d for the grain element method at the same instants during the lifetime. In Fig. 10a we see that damage has developed non-uniformly over the transverse facets and that coalescence has taken place at the triple point nearest to the microcrack and on the upper right facet. Comparison of Figs. 10a and c shows that the grain element method is able to capture the key features of failure propagation. Cavitation along the facets is predicted

P. Onck, E. van der Giessen / Mechanics of Materials 26 (1997) 109-126

121

t = tf

(a)

I

i (b)

t=0

Fig. 10. The state of damage at two instants t / t f in a quarter unit-cell for ( a / L ) I = 0.025 and ,~m/,~e = 0.577. Values of a / b are plotted along, and on either side of the grain boundary facets. Microcracked regions where a / b = 0.7 are indicated by the darker gray scale. (a) Detailed analysis, t / t f :: 0.83; (b) Detailed analysis, t / t f = 1; (c) Grain element method, t / t f = 0.83; (d) Grain element method, t / t f ~ 1.

in an average manner and the distribution of damage over the unit-cell is simulated accurately (note the shielded region above the microcrack). Figs. 10b and d show the damage state at the moment of failure. Three microcracks have developed and opened up due to the sliding off of adjacent grains. Despite the fact that the upper right facet has cracked early, the aggregate fails by a linkage of transverse microcracks from the initial microcrack to the right cell boundary under an angle of 30 ° with the xl-axis. This too is what the grain element characterization predicts. It is worthwhile to note that the gain in computer time is very significant: 20 h (detailed finite element description) versus 2 rain for the grain element description (i.e. a factor 600) on a relatively fast workstation.

To check the influence of the stress triaxiality and the initial value of a/L, calculations have been repeated for ( a / L ) I = 0.1 and .~m/2e = 4. The resuits, in terms of failure times, are summarized in Table 1. It is noted that for all cases analyzed the microcracking process proceeds along the same patTable 1 Normalized times to failure for a (6, 5) unit-cell with an initial microcrack for both the detailed analysis (da) and the grain element method (gem) for two values of ( a / L ) 1 and .~m/~e

(a/L) I

~,,/~

t f / t R (da)

t f / t R (gem)

Deviation(%)

0.025 0.1 0.025 0.1

0.577 0.577 4 4

11.25 37.99 0.66 3.16

11.16 40.91 0.72 3.37

-0.8 7.7 9.1 6.6

122

P. Onck, E. van der Giessen / Mechanics of Materials 26 (1997) 109-126

tern as shown in Fig. 10. For the higher value of ( a l L ) I the cavitation along the facets is more nonuniform (Van der Giessen and Tvergaard, 1994b) and deformation in the grains is much larger. In the higher triaxiality cases, cavitation also develops on the inclined facets, but this does not exceed a / b = 0.2. The resemblance between the damage distribution over the polycrystal according to the grain element method and the detailed analysis is similar as in Fig. 10 and needs not be shown.

(a) ~

~

(b)

4.3. Random variations In this subsection, the influence of variations in grain shape on the rupture process is investigated. Variations are introduced by giving the triple grain junctions of a (6, 5) unit-cell (cf. Fig. 1) a random offset. We consider six different realizations which were analyzed earlier by Van der Giessen and Tvergaard (1994c) using detailed finite element computations. The different microstructures are shown in Fig. 11. Note that special care has been taken for the vertices on the cell boundary in order to maintain symmetry. Uniaxial tension along the x2-axis is applied and the parameter ( a l L ) l is chosen to be 0.025. The material is virtually undamaged initially, ( a / b ) I = 0.01 and ( b / R ) I = 1. Since the initial cavity density N I = 1/rrR~ varies from facet to facet, we specify the nucleation parameter F, to be equal to 100/(~'R 2) in terms of the initial facet radius R 0 in the regular hexagonal microstructure. Fig. 12 compares the evolution of damage for the grain element method and detailed analysis, similar to Fig. 10, for the microstructure of Fig. lid. As discussed by Van der Giessen and Tvergaard (1994c), the results shown in Fig. 12a indicate that cavitation on the transverse facets is roughly inversely proportional to the facet width and apparently the grain element method is able to pick this up accurately (Fig. 12c). Figs. 12b and d show that three microcracks have developed, which means that the lifetime of the specimen is nearly exhausted. In contrast to the regular hexagonal structure, final failure occurs when a percolation of microcracks (transverse and inclined) has developed from the left to the fight hand side of the cell. For all microstructures of Fig. 11 the failure patterns for the grain element method and the detailed analysis were similar, although the

(c)

(e)

(d)

_/

?-j

(f)

Fig. 11. Different realizations of random variations in the geometry of grains inside a quarter of the (6, 5) unit-cell. From (Van der Giessen and Tvergaard, 1994c).

order at which facets cracked in time was sometimes slightly different. Fig. 13 shows the creep curve for both the grain element method and detailed analysis for the microstructure of Fig. 11 d. The grain elements are seen to underestimate the lifetime by about 17%. Note that in the early stages of the creep curve, when there is negligible damage on the facets, the secondary creep rates agree well. This is because we corrected the constitution of the grain elements for the enhancement of creep due to sliding (cf. Section 3.2). It should be mentioned in this respect that here we determined the stress enhancement factor f * for each microstructure numerically and not through Eq. (25) in order to check the performance of the grain elements as accurately as possible. The approximate relation Eq. (25) is used only for very large polycrystalline aggregates that are beyond the grasp of the full-field finite element calculations. Table 2 summarizes the results for the different microstructures shown in Fig. 11. As noted first by Van der Giessen and Tvergaard (1994c), the life-

123

P. Onck, E. van der Giessen / Mechanics of Materials 26 (1997) 109-126

t= tf

(a)

(b) t=O

i

(c)

(d)

I

Fig. 12. The state of damage at two different stages t / t f in a quarter unit-cell with the microstructure of Fig. 1 ld for uniaxial tension and ( a l L ) I = 0.025. (a) Detailed analysis, t / t r = 0.75; (b) Detailed analysis, t / t r = 0.99; (c) Grain element method, t / t t = 0.75; (d) Detailed analysis, t i t e = 0.99.

l

0.08

/

0.0~

t i m e s f o r all m i c r o s t r u c t u r e s a r e s h o r t e r t h a n f o r t h e r e g u l a r h e x a g o n a l g r a i n s t r u c t u r e . A s e x p l a i n e d further by Onck and Van der Giessen (1995), the two

/

0.04

Table 2 Stress enhancement factor f ° and normalized time to failure for the microstructures of Fig. 11 for uniaxial tension and ( a / L ) I = 0.025

------ Grain element method ,~

0.02

I

2

3

4

5

6

7

8 t/tR

Fig. 13. Creep curve for the microstructure of Fig. l i d for uniaxial tension and ( a / L ) t = 0.025.

Figure

f *

t f / t R (da)

t f / t R (gem)

Deviation (%)

Fig. Fig. Fig. Fig. Fig. Fig.

1.224 1.187 1.251 1.296 1.389 1.232

15,76 23.90 15.14 8.71 5.25 18.54

15.87 22.85 17.08 7.19 4.07 16.82

0.7 -4.4 12.8 - 17.5 -22.5 -9.3

ll(a) ll(b) ll(c) ll(d) 11(e) ll(f)

124

P. Onck, E. van der Giessen / Mechanics of Materials 26 (1997) 109-126

main reasons for this are that random variations tend to increase the free sliding creep enhancement, represented by f * , but also increase the normal stresses on small transverse facets. From Fig. 12 and Table 2 it is clear that the grain elements are able to pick up this phenomenon, while the deviation from the fullfield finite element calculations remains within reasonable bounds.

5. Summary and discussion A microstructurally-based numerical method has been presented to simulate intergranular creep fracture. So-called grain elements are used to account for the average elastic and creep deformation of individual grains, while grain boundary elements account for sliding and cavitation. This method allows for the analysis of intergranular creep fracture in large polycrystalline aggregates, based on the physical mechanisms active on a much smaller length scale. Two essential steps in the formulation of the grain elements can be identified. Firstly, the choice of finite elements to represent the grain element is not trivial. One has to account for locking phenomena which can occur in case of creep constrained cavitation (see Appendix A). Secondly, the enhanced creep flow in the grains due to free grain boundary sliding must be incorporated, since it affects the failure process considerably. This is accomplished by substituting a geometry-dependent stress enhancement factor f * in the creep law. To establish the accuracy of the grain element method, results are compared with detailed finite element calculations for a multi-grain unit-cell model. Creep constraint is an essential feature of diffusive cavitation in front of a macroscopic crack as well as in a uniformly stressed material (Dyson, 1976). To investigate the performance of the grain elements for this, cavitation on an isolated facet was analyzed for relevant material parameters and stress triaxialities and the agreement was found to be excellent. Another principal phenomenon of intergranular fracture is the subsequent linking-up of microcracks. We monitored the evolution of damage from an initial microcrack and found that the grain elements were capable of predicting the crack pattern as well as the

ultimate time to failure. Finally, we considered the effect of random variations in microstructure. As compared with the regular hexagonal microstructure this results in higher stresses on small transverse facets and an increase in the strain rate enhancement. Both these aspects result in a decrease of the lifetime, which is shown to be captured by the grain element method with a reasonable accuracy. As a great deal of detail is lost when using the grain elements, a huge gain in computer time can be obtained: a factor of 600 for the present cell size. It must be noted that the accuracy of the grain element method, reported on in this paper, has been checked by means of unit-cell calculations. This assumes a uniform stress state at the macroscopic level, which differs from the high gradient fields near a macroscopic crack tip. However, we have shown that essential features of creep fracture, like creep constrained cavitation and the interlinkage of microcracks, can be captured very accurately at the unit-cell level. Based on our experience with these and the detailed analyses (Van der Giessen and Tvergaard, 1994a,b,c) we are confident about the applicability for creep crack problems; the precise accuracy in those cases, however, cannot be readily established. As compared to the Delaunay network modelling technique of Van der Burg and Van der Giessen (1994a,b), the grain element method does not need any tuning or calibration with respect to the detailed analysis. Furthermore, in contrast to this lattice-type of model, the grain element method is able to pick up the stress state dependence of creep rupture, which is especially important because of the hydrostatic stress effect on void growth. The gap between the microscopic and macroscopic or continuum level is usually filled by employing continuum damage methods, in which damage is modelled in an average, smeared-out manner. This can lead to pathological results when continuum mechanics is applied at length scales comparable to or smaller than the characteristic dimensions of the microstructure. The grain element method, however, can be seen as a mixed continuum/discrete model in which the microstructure is explicitly represented, so that macroscopic creep crack growth can be modelled in a discrete fashion. This will be reported elsewhere.

P. Onck, E. oan der Giessen / Mechanics of Materials 26 (1997) 109-126

AppendixA The driving force for the development of the grain elements is to simulate creep fracture with as little computational effort as possible, but still capturing the essential micromechanical processes. To do so we represent the grains by six-noded super elements to account for elasticity and creep. An obvious choice would be to use six constant strain triangles to subdivide the hexagonal grains and condense out the midnode. However, this would lead to volumetric locking as will be demonstrated by way of the following (academic) problem. Suppose a limited number of transverse grain boundary facets are prone to cavitation and are located in columns above each other (see Fig. 14a). The other facets are allowed to slide freely but will not cavitate. Uniaxial tension is applied in a direction normal to the cavitating facets. Upon loading, the stresses will redistribute due to grain boundary sliding, leaving the inclined facets stress free and concentrating load on the transverse facets. When grain boundary diffusion is much faster than creep (corresponding to a small value for a l L , cf. Section 2) the normal stresses on the cavitating facets relax, thus shedding load to the neighboring columns of non-

I y.

125

cavitating grains. In this situation cavitation is strongly constrained by creep of the surrounding grains, so that a separation of the grains due to cavitation is only possible if the surrounding grains deform. However, the resulting deformation mode of these grains cannot be accomplished by the six triangles and the mesh locks. The reason for the failure of the six triangles is that the grains must elongate in such a way that the angle of the inclined facets with the tensile direction must remain 30 °, while plastic incompressibility demands (elastic strains are small) must be satisfied for the grain and the separate triangles. It appears that the ratio of degrees of freedom to incompressibility constraints is too low in this grain element representation to account for the deformation. To solve this problem we use 'crossed triangle' quadrilaterals, which are known for their good performance in incompressible materials (Nagtegaal et al., 1974). The orientation of the quadrilaterals in the grain is such that their connecting edge is normal to the maximum principal stress direction. Using these elements allows the cavitating facets to fracture and successively open up without numerical artifacts (see Fig. 14b). The underlying reason for the success of the

l y,

(a) Fig. 14. (a) Polycrystal representation of grain elements consisting of six constant strain triangles. Cavitation is assumed to occur on a limited number of transversefacets, arranged in columns. Mesh locking occurs. (b) The same situation as in (a), but with the present grain element representation.

126

P. Onck, E. van der Giessen / Mechanics of Materials 26 (1997) 109-126

c r o s s e d t r i a n g l e q u a d r i l a t e r a l s is that o n l y t h r e e o f t h e f o u r i n c o m p r e s s i b i l i t y c o n s t r a i n t s are i n d e p e n dent, d u e to the c r o s s e d triangle a r r a n g e m e n t ( N a g t e g a a l et al., 1974). H o w e v e r , w h e n large deform a t i o n s take place, the m i d n o d e n o l o n g e r r e m a i n s o n the d i a g o n a l , so that the p e r f o r m a n c e o f the e l e m e n t s d e t e r i o r a t e s . T o fix this, w e s i m p l y r e s t r a i n the m i d n o d e f r o m m o v i n g o f f the d i a g o n a l , thus i n t r o d u c i n g a s m a l l i n c o n s i s t e n c y in the e l e m e n t f o r m u l a t i o n . H o w e v e r , the r e s u l t i n g e r r o r r e m a i n s s m a l l c o m p a r e d to the error m a d e in c o m p a r i s o n w i t h the d e t a i l e d finite e l e m e n t a n a l y s i s (see S e c t i o n 4) a n d m u c h s m a l l e r t h a n that due to locking.

References Argon, A.S., 1982. Mechanisms and mechanics of fracture in creeping alloys. In: Wilshire, B., Owen, D.R.J. (Eds.), Recent Advances in Creep and Fracture of Engineering Materials and Structures. Pineridge press, Swansea, pp. 1-52. Ashby, M.F., 1972. Boundary defects and atomistic aspects of boundary sliding and diffusional creep. Surf. Sci. 31,498-542. Chen, I.-W., 1983. Cavity growth on a sliding grain boundary. Metall. Trans. 14A, 2289-2293. Cocks, A.C.F., Ashby, M.F., 1982. On creep fracture by void growth. Prog. Mater. Sci. 27, 189-244. Crossman, F.W., Ashby, M.F., 1975. The non-uniform flow of polycrystals by grain-boundary sliding accommodated by power-law creep. Acta Metall. 23, 425-440. Dyson, B.F., 1976. Constraints on diffusional cavity growth rates. Metal Sci. 10, 349-353. Dyson, B.F., 1983. Continuous cavity nucleation and creep fracture. Scr. Metall. 17, 31-37. Ghahremani, F., 1980. Effect of grain boundary sliding on steady creep of polycrystals. Int. J. Solids Struct. 16, 847-862. Hsia, K.J., Parks, D.M., Argon, A.S., 1991. Effects of grain boundary sliding on creep-constrained boundary cavitation and creep deformation. Mech. Mater. 11, 43-62. Hayhurst, D.R., Brown, P.R., Morrison, C.J., 1984. The role of continuum damage in creep crack growth. Philos. Trans. R. Soc. London Ser. A: 311, 131-158. Hayhurst, D.R., Leckie, F.A., 1984. Behaviour of materials at high temperatures. In: Carlsson, J., Ohlson, N.G. (Eds.), Mechanical Behaviour of Materials. Proceedings of ICM4, vol. 2. Pergamon Press, Oxford, pp. 1195-1212. Hutchinson, J.W., 1983. Constitutive behavior and crack tip fields

for materials undergoing creep-constrained grain boundary cavitation. Acta Metall. 31, 1079-1088. Nagtegaal, J.C., Parks, D.M., Rice, J.R., 1974. On numerically accurate finite element solutions in the fully plastic range. Comp. Methods Appl. Mech. Eng. 4, 153-177. Needleman, A., Rice, J.R., 1980. Plastic creep flow effects in the diffusive cavitation of grain boundaries. Acta Metall. 28, 1315-1332. Needleman, A., Tvergaard, V., 1984. Finite element analysis of localization in plasticity. In: Oden, J.T., Carey, G.F. (Eds.), Finite Elements, Special Problems in Solid Mechanics. Prentice-Hall, New Jersey, pp. 94-157. Onck, P.R., Van der Giessen, E., 1997. Influence of microstrucrural variations on steady state creep in 2D freely sliding polycrystals. Int. J. Solids Struct. 34, 703-726. Peirce, D., Shih, C.F., Needleman, A., 1984. A tangent modulus method for rate dependent solids. Comp. Struct. 18, 875-887. Rice, J.R., 1981. Constraints on the diffusive cavitation of isolated grain boundary facets in creeping polycrystals. Acta Metall. 29, 675-681. Sham, T.-L., Needleman, A., 1983. Effects of triaxial stressing on creep cavitation of grain boundaries. Acta Metall. 31,919-926. Tvergaard, V., 1984. On the creep constrained diffusive cavitation of grain boundary facets. J. Mech. Phys. Solids 32, 373-393. Tvergaard, V., 1986. Analysis of creep crack growth by grain boundary cavitation. Int. J. Fract. 31, 183-209. Van der Burg, M.W.D., Van der Giessen, E., 1994a. Delaunaynetwork modelling of creep failure in regular polycrystalline aggregates by grain boundary cavitation. Int. J. Dam. Mech. 3, 111-139. Van der Burg, M.W.D., Van der Giessen, E., 1994b. Simulation of microcrack propagation in creeping polycrystals due to diffusive grain boundary cavitation. Appl. Mech. Rev. 47, S122-S131. Van der Giessen, E., Tvergaard, V., 1990. On cavity nucleation effects at sliding grain boundaries in creeping polycrystals. In: Wilshire, B., Evans, R.W. (Eds.), Creep and Fracture of Engineering Materials and Structures. Elsevier, Swansea, pp. 169-178. Van der Giessen, E., Tvergaard, V., 1991. A creep rupture model accounting for cavitation at sliding grain boundaries. Int. J. Fract. 48, 153-178. Van der Giessen, E., Tvergaard, V., 1994a. Interaction of cavitating grain boundary facets in creeping polycrystals. Mech. Mater. 17, 47-69. Van der Giessen, E., Tvergaard, V., 1994b. Development of final creep failure in polycrystalline aggregates. Acta Metall. Mater. 42, 959-973. Van der Giessen, E., Tvergaard, V., 1994c. Effect of random variations in microstructure on the development of final creep failure in polycrystalline aggregates. Modelling Simul. Mater. Sci. Eng. 2, 721-738.