Application of embedded discontinuities for softening solids

Engineering Fracture Mechanics 65 (2000) 263±281

www.elsevier.com/locate/engfracmech

Application of embedded discontinuities for softening solids G.N. Wells a,*, L.J. Sluys b a

Faculty of Aerospace Engineering, Koiter Institute Delft/Delft University of Technology, P.O. Box 5048, 2600 GA Delft, The Netherlands b Faculty of Civil Engineering and Geosciences, Koiter Institute Delft/Delft University of Technology, P.O. Box 5048, 2600 GA Delft, The Netherlands

Abstract Additional, discontinuous functions are added to the displacement ®eld of standard ®nite elements in order to capture highly localised zones of intense straining. By embedding discontinuities within an element it is possible to e€ectively model localisation phenomena (such as fracture in concrete) with a relatively small number of ®nite elements and without resorting to higher-order continuum theories. However, the application of continuum constitutive laws in combination with embedded discontinuities, several unique performance related issues arise. The source of these problems and some remedies are addressed here. Numerical results are compared with smeared type models for brittle fracture. # 2000 Elsevier Science Ltd. All rights reserved. Keywords: Strong discontinuity; Embedded crack; Plasticity; Localisation; Enhanced strain

1. Introduction Recently, attention has been paid to so-called `embedded discontinuity' models as a means to overcome the limitations of smeared and discrete crack models. By capturing a crack within an element, it is e€ectively possible to model a discrete (or highly localised) phenomenon within a continuum framework. Embedded discontinuities o€er a hybrid approach between smeared and discrete modelling. The de®ciency of discrete crack models is the dependence on mesh alignment, when separation can occur only at elements interfaces. To remedy this, special, computer intensive re-meshing techniques are required that attempt to predict the crack propagation direction. The weakness of smeared crack models is their inherent dependency on element size and alignment [1], arising from the ill-posed nature of the classical continuum when applied for softening materials [2]. Attempts to overcome element size * Corresponding author. Fax: +31-15-261-1465. E-mail address: [email protected] (G.N. Wells). 0013-7944/00/$ - see front matter # 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 1 3 - 7 9 4 4 ( 9 9 ) 0 0 1 2 0 - 4

264

G.N. Wells, L.J. Sluys / Engineering Fracture Mechanics 65 (2000) 263±281

dependency through simple fracture energy regularisation are only partly successful, as this method provides no information over crack direction nor is the necessary length scale (so-called crack bandwidth) mesh objective. By embedding a discontinuity, direction information is taken from the local stress (or strain) ®eld and kinematic enhancements make elements more able to properly reproduce the real deformation modes. The application of gradient [3,4] and non-local [5,6] continuum theories has proved successful in overcoming the de®ciency of the classical continuum when analysing softening materials. However, the common feature of higher-order models is the need for a ®ne mesh within the localisation zone in order to capture the high strain gradients. This requires a priori knowledge of where failure will occur and limits the potential of such methods for structural scale problems, especially in three-dimensions. With embedded discontinuities, the width of the failure zone is less than the width of a single element, making the method ideal for three-dimensional and very large scale applications. Embedded discontinuity models can be placed broadly in two classes; weak and strong discontinuities. Weak discontinuities involve the addition of discontinuous modes to the strain ®eld without regard to the displacement ®eld [7,8]. This approach avoids the need to deal with unbounded strains as a result of a displacement jump, but the ability of such models to overcome mesh-alignment dependence [9] and to reproduce full crack separation in mode-I type problems is questionable [10] since elements are not kinematically enhanced. The model developed and applied here incorporates a strong discontinuity. The common feature of strong discontinuities is that a jump is added in the displacement ®eld of a ®nite element [11±15]. Additional, discontinuous shape functions are added to an element to represent the deformations in conjunction with a specially constructed interpolation function that is used in calculating the internal force vector. This method is able to fully represent the kinematics of the problem as well as neatly satisfy traction continuity across the discontinuity. The formulation presented here follows consistent steps from the three-®eld variational statements, extended to incorporate a jump in displacements. It has been shown that under certain conditions, continuum constitutive laws and unbounded strains are compatible [11±13]. This then allows strong discontinuities to be treated within a continuum framework. Through the use of a regularised discontinuity jump, conventional continuum relationships can be applied, avoiding the need to formulate discrete traction-displacement relationships. The implications of applying continuum constitutive laws in discontinuity analysis are examined, with special attention to the performance in mode-I type problems with Rankine plasticity. For a benchmark problem, the embedded discontinuity model is compared with a conventional smeared model using crack bandwidth regularisation, with both models using the same constitutive relationship. The performance of the embedded model with respect to mesh sensitivity is shown to be superior to the crack bandwidth regularised models, although the application of a regularised displacement jump can lead to spurious behaviour unique to embedded discontinuities, namely stress locking in the absence of signi®cant rotation of the stress ®eld.

2. Strong discontinuity formulation To begin, consider a body O that is crossed by a material discontinuity (see Fig. 1). The body is composed of two sub-domains, O‡ and Oÿ , that lie on either side of the discontinuity Gd : In addition, it is useful to de®ne a discontinuity band …Od †, centred on Gd and a sub-domain of O: To formulate a model that contains a discontinuity ®rst the kinematic ®elds are examined, then the variational statements to impose traction continuity across the discontinuity.

G.N. Wells, L.J. Sluys / Engineering Fracture Mechanics 65 (2000) 263±281

265

Fig. 1. Body containing a discontinuity.

2.1. Discontinuous displacement ®eld The displacement ®eld of a body containing a discontinuity can be decomposed into two parts; a Ä Ã t†† and a discontinuous …HGd ‰uŠ…x, t†† component. The jump in displacement is provided continuous …u…x, via the Heaviside function …HGd †, centred on the discontinuity and operating on a continuous function Ä ‰uŠ…x, t†: However, in the context of ®nite element implementation, this presents diculties since both components must be taken into account when imposing essential boundary conditions. A more suitable approach is to specially develop the decomposition such that the contribution of the discontinuous component is zero on the boundaries [11,14]. Hence, only the regular, continuous part of the displacement ®eld must be taken into account.
ÿ à t† ‡ HGd ÿ j…x† ‰uŠ…x, t† …1† u…x, t† ˆ u…x, In Eq. (1) ‰uŠ…x, t† is the magnitude of the discrete displacement jump at the discontinuity and j…x† is a continuous function that satis®es  1 in O‡ %Od …2† j…x† ˆ 0 in Oÿ %Od The one-dimensional displacement decomposition is shown in Fig. 2. The strain ®eld can be easily found by calculating the gradient of Eq. (1).

Fig. 2. Displacement decomposition.

266

G.N. Wells, L.J. Sluys / Engineering Fracture Mechanics 65 (2000) 263±281

e …x, t† ˆ rs uà ‡

ÿÿ



dGd nd ÿ rs j ‰uŠ

…3†

Note the appearance of the Dirac delta distribution …dGd † in the strain expression. From the inclusion of a displacement jump, the strain ®eld at the discontinuity is unbounded. 2.2. Mixed variational formulation To construct an element that incorporates a displacement discontinuity, it is necessary to start from the variational formulation. It is useful to begin from a modi®ed form of the three-®eld variational statements that includes so-called enhanced strains [16]. The strain ®eld is decomposed into a compatible …rs u† and incompatible …~e † part, yielding … …4a† rs Z :s s…rs u ‡ e~ † dO ÿ W ext …Z † ˆ 0 O

… O

… O

t :~e dO ˆ 0

…4b†

g :… ÿ s ‡ s …rs u ‡ e~ †† dO ˆ 0

…4c†

where …ZZ, g , t † 2 …V  E  S † are all variations of displacements, strains and stresses, respectively, with …V  E  S † the spaces of admissible displacement, strain and stress variations. In order to allow for a strong discontinuity, the variational statements in Eq. (4c) must be further modi®ed by decomposing the enhanced strains into continuous …~e c † and discontinuous …~e d † parts. Taking into account that the discontinuous part of the enhanced strain ®eld is non-zero only on Gd and using the divergence theorem, Eq. (4c) can be re-phrased allowing for a displacement jump as … …5a† rs Z :s s…rs u ‡ e~ † dO ÿ W ext …Z † ˆ 0 O

… O

… O

… t :~e c dO ‡

Gd

…t n d †
nd dG ˆ 0

g~ c :… ÿ s ‡ s …rs u ‡ e~ †† dO ‡

…5b†

…

 Gd

 … ÿ s ‡ s …rs u ‡ e~ ††Z d
nd dG ˆ 0

…5c†

where n d is the discrete displacement at the discontinuity and Z d is a variation of the discrete displacement at the discontinuity. From this point, there have been two distinct approaches. The ®rst is the kinematic approach, where extra discontinuous shape functions are added to the element [14]. This method is kinematically representative of the real deformations and allows full separation at the interface, but leads to diculties in imposing the traction continuity conditions in Eq. (5). The alternative is to use enhanced assumed strains (EAS), where the enhanced modes are constructed in the strain space, orthogonal to the stress ®eld, thereby automatically satisfying Eq. (5b) and eliminating the stresses from the unknown ®elds. This approach allows traction continuity across the discontinuity to be easily satis®ed, although it is kinematically lacking with respect to the former approach. The approach followed here is a hybrid between the kinematic and the EAS methods. The kinematic

G.N. Wells, L.J. Sluys / Engineering Fracture Mechanics 65 (2000) 263±281

267

approach is used to enhance the kinematic ®elds, while the EAS approach is used in calculating the internal forces and hence impose the traction continuity [12,13,17]. Essentially, the variational statements in Eq. (5) are followed, but unlike normal Galerkin procedures, the space of the enhanced strain ®eld …~e † is not necessarily the same as the space of admissible variations of the enhanced strains …~g †: This hybrid approach exploits the advantages of both the kinematic and the EAS methods. 2.3. Finite element implementation The discretised form of the displacement and strain ®elds including a discontinuity, for timeindependent models, can be written as follows: u…x† ˆ Na …x†ae ‡ Na …x†aae

…6†

e …x† ˆ B…x†ae ‡ G…x†aae

…7†

where Na and B are the normal displacement and strain interpolation matrices, respectively; Na and G (gradient of Na † collect the interpolation functions for the enhanced displacement and strain ®elds, respectively; ae are the nodal displacements; a e are the displacements associated with the enhanced modes at the discontinuity. Taking into account the form of the normal matrix, which in threedimensions has the form 2 3 nx 0 0 ny 0 n z …8† nTd ˆ 4 0 ny 0 nx nz 0 5 0 0 n z 0 ny nx where nx , ny and nz are the normal components to the discontinuity …Gd †, the matrix G for an element is constructed by examining the enhanced part of the kinematic strain ®eld in Eq. (3) 3 2 @ j…x† 0 0 7 6 d Gd n x ÿ @ x 7 6 7 6 @ j…x† 7 60 0 d Gd n y ÿ 7 6 @y 7 6 7 6 6 @j…x† 7 7 60 n ÿ 0 d Gd z 6 @z 7 7 6 …9† Ge ˆ 6 7 @j…x† @j…x† 7 6 7 6 d Gd n y ÿ d Gd n x ÿ 0 7 6 @y @x 7 6 6 @ j…x† @j…x† 7 7 6 dGd n y ÿ d Gd n z ÿ 7 60 6 @z @y 7 7 6 4 @ j…x† @j…x† 5 0 dGd n x ÿ d Gd n z ÿ @z @x Next, attention must be paid to the traction continuity condition. Eq. (5b) is automatically satis®ed by the orthogonality requirement of EAS, allowing the stress ®eld …s s† to be eliminated from Eq. (5c). Eq. (5c) is then rephrased as … … ÿ
s …rs u ‡ e~ †Z d
nd dG ˆ 0 g~ c :s s…rs u ‡ e~ † dO ‡ …10† O

Gd

268

G.N. Wells, L.J. Sluys / Engineering Fracture Mechanics 65 (2000) 263±281

Considering then that Eq. (5b) must apply for each element gives … … … g~ :tt dO ˆ g~ c :tt dO ‡ …t Z d †
nd dG ˆ 0 Oe O e % Gd Gd, e

…11†

For a straight discontinuity (that is, nd is constant within an element), a piecewise constant strain ®eld and taking into account that in order to pass the patch test [18], Eq. (11) must hold for an arbitrary constant stress ®eld, Eq. (11) can be integrated explicitly. Ae g~ c ‡ le …nd Z d † ˆ 0

…12†

Rearranging Eq. (12), an expression for the strain ®eld outside the discontinuity can be found. g~ c ˆ ÿ

le nd Z d Ae

…13†

For the two-dimensional case, le is the length of the discontinuity through the element and Ae is the area of the element. For the three-dimensional case, le is the area of the discontinuity plane through the element and Ae is the volume of the element. It can be shown that the actual length of the discontinuity plays no role in the ®nal numerical result [19]. For the case of a CST element with a discontinuity aligned with an element boundary, using the actual length of the discontinuity leads to a symmetric sti€ness matrix. Taking the traction continuity condition (10), which must also apply for each element, and now expressing s …rs u ‡ e~ † as s …ee† … … ÿ
g~ c :s s…ee † dO ‡ …14† s …ee †ZZd
nd dG ˆ 0 Oe

Gd

Substituting the continuous part of the enhanced strain ®eld (13) into Eq. (14)   … … ÿ
le ÿ s…ee † dO ‡ nd Z d :s s …ee †ZZd
nd dG ˆ 0 Ae O e % Gd Gd Now using properties of the Dirac-delta distribution for the term integrated over Gd ,   … … ÿ
le nd Z d :s dGd s …ee †ZZd
nd dO ˆ 0 ÿ s…ee † dO ‡ A e O e % Gd Oe

…15†

…16†

Eq. (15) can be rearranged to express the traction force at the discontinuity and in the bulk of the element. … … 1 1 s …ee †nd dO ˆ s …ee †nd dG …17† Ae Oe nGd , e le Gd, e |‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚ ‚{z‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚‚} |‚‚‚‚‚‚‚‚‚‚‚‚‚{z‚‚‚‚‚‚‚‚‚‚‚‚‚} For constant strain elements (three-noded triangle (CST) and four-noded tetrahedron), the constant strain ®eld assumption used to integrate Eq. (13) is exact, so the average tractions (17) are also the exact solution. For higher-order elements, the traction continuity requirement is imposed in an average (weak) sense. For ®nite element implementation, it is necessary to form an interpolation matrix …G † for calculating the internal force vector associated with the degrees of freedom at the discontinuity.

G.N. Wells, L.J. Sluys / Engineering Fracture Mechanics 65 (2000) 263±281

f int ˆ

…  O

BT s GT s



 dO ˆ

 int

fu f int a

Taking variations of Z d in Eq. (16), the interpolation matrix Ge can be formed   le  nd Ge ˆ dGd ÿ Ae From the variational statements (5a) and (5c), the element sti€ness matrix can then be formed.    ( int ) …  T _ a_ BTe DGe Be DBe dO e ˆ f u, e T T a_ e Ge DGe 0 Oe Ge DBe

269

…18†

…19†

…20†

As with EAS methods, the enhanced, incompatible modes can be statically condensed at element level. This makes the method ecient and simple to implement in existing ®nite element codes. Notice that the hybrid approach, in general, results in a non-symmetric element sti€ness matrix, irrespective of the constitutive model applied. This is a direct result of the hybrid kinematic-EAS approach. 2.4. Application of continuum constitutive laws In the previous section, a strain tensor was formed at the discontinuity based on a kinematic decomposition of the displacement ®eld. The formation of a strain tensor implies that continuum constitutive laws can be applied. The application of continuum constitutive laws relies primarily on the physical justi®cation that unbounded strain ®elds are notionally possible and that unbounded stresses are impossible. That is, irrespective of the strain ®eld, the stress ®eld is bounded at all times. From a physical point of view, an unbounded strain ®eld arises naturally in the case of a displacement discontinuity and results in the inclusion of the Dirac-delta distribution in the interpolation matrices (9) and (19). From the natural stress boundedness requirement, the application of continuum constitutive laws is conditional upon the chosen constitutive law exhibiting a softening branch. Any hardening behaviour at the discontinuity would lead to an unbounded stress ®eld. For the implementation examined here, a softening inelastic material model is applied at the discontinuity with linear-elastic behaviour outside the discontinuity. 2.4.1. Plasticity For plasticity models, it is simplest to begin from the inverse of the classical elasto-plastic tangential sti€ness matrix for hardening/softening associative plasticity (hardening modulus (H ) is not equal to zero)   1 ÿ1 e T …21† e_ ˆ …D † ‡ mm s_ H with mˆ

@F @s s

…22†

where De is the elastic tangent matrix and F the yield function. The expression for the strain ®eld in Eq. (3) can be rearranged, grouping the bounded, continuous strains …e† and the unbounded parts. e ˆ e ‡ dGd nd ‰uŠ

…23†

270

G.N. Wells, L.J. Sluys / Engineering Fracture Mechanics 65 (2000) 263±281

The expression for the strain ®eld in Eq. (23) can be di€erentiated with respect to time, then substituted into Eq. (21). …Dÿ1 †T s_ ˆ 1 mmT s_ e_ ‡ dGd nd ‰uÇ Š ˆ |‚‚‚‚{z‚‚‚‚} |{z} |‚‚‚‚{z‚‚‚‚} H |‚‚‚{z‚‚‚}

…24†

To impose the stress boundedness condition, the unbounded term in Eq. (24) containing the Dirac-delta distribution must be cancelled by the hardening modulus. This can be done simply by distinguishing  [11,12] and relating the two between the hardening modulus (H ) and the intrinsic hardening modulus …H† by 1 1 ˆ d Gd H H

…25†

The intrinsic softening parameter can be found easily by integrating the one-dimensional linear stressdisplacement relationship (shown in Fig. 3),   f ‡ H‰uŠ ‰uŠ < ‰uŠd …26† s ˆ t 0 ‰uŠr‰uŠd giving f 2 H ˆ ÿ t 2Gf

…27†

where ft is the uni-axial tensile strength and Gf the fracture energy. For exponential softening,   ft s ˆ ft exp ÿ ‰uŠ Gf

…28†

Eq. (28) can be di€erentiated with respect to u, to give the intrinsic softening parameter in terms of ‰uŠ:   f 2t ft  …29† H…u† ˆ ÿ exp ÿ ‰uŠ Gf Gf E€ectively, if the inverse of the softening modulus is negative and has the distributed form of the Dirac-

Fig. 3. 1D linear stress-displacement relationship.

G.N. Wells, L.J. Sluys / Engineering Fracture Mechanics 65 (2000) 263±281

271

delta distribution, the stress ®eld will always be bounded (at least for the one-dimensional case shown in Fig. 3). 2.4.2. Rankine model At the discontinuity, a multi-surface Rankine plasticity model with isotropic softening is applied in order to examine mode-I failure. The uncoupled yield surfaces of the Rankine model are given in Eq. (30), where si are the principal stresses, s comes from the hardening law and ft is the uni-axial tensile strength. For plane strain, i ˆ 1 to 3. Fi ˆ s i , Fi ÿ si R0

…30†

From Eq. (26) for linear softening, s in terms of the hardening modulus (H ) is given by considering Eq. (25). s ˆ ft ÿ

1 f 2t ‰uŠ dGd 2Gf

For exponential softening, the hardening modulus is   1 f 2t ft exp ÿ ‰uŠ H…u† ˆ ÿ Gf dGd Gf  which can be integrated to give s:   1 ft s ˆ ft exp ÿ ‰uŠ dGd Gf

…31†

…32†

…33†

Note the similarity of Eqs. (31) and (33) to the softening laws used in smeared crack bandwidth formulations. Simply the crack bandwidth measure has been replaced by the inverse of the Dirac-delta distribution. This makes application of existing constitutive models designed for smeared modelling with crack bandwidth regularisation possible without any special modi®cations. For computational eciency, in forming the element sti€ness matrix a consistent tangent can be used. 2.5. Displacement jump regularisation The ®nal step to apply the method numerically is to deal with the Dirac-delta distribution that appears in the hardening modulus and in the enhanced strain matrices. Since unbounded strains do not ®t with numerical implementation, the Dirac-delta distribution is regularised with a function that satis®es the necessary integration properties of the distribution [12] dG 1

1 k

…34†

As k 40, the approximation (34) approaches the exact solution. For crack bandwidth-type constitutive laws, the measure of bandwidth is simply replaced by k.

3. Application and performance In the following sections the embedded discontinuity model is examined for the particular case of

272

G.N. Wells, L.J. Sluys / Engineering Fracture Mechanics 65 (2000) 263±281

mode-I type failure using Rankine plasticity. In this context, the terms crack, embedded crack are appropriate and are used interchangeably with embedded discontinuity. 3.1. Crack path determination Cracks within an element are approximated as ¯at planes in three dimensions and as straight lines in two dimensions, with the direction of the plane determined by the stress or strain ®eld at the integration point and the type of material being modelled. For the Rankine model, the discontinuity is inserted when the yield surface is ®rst violated. The discontinuity is then aligned normal to the major principal tensile stress. For other constitutive laws, the normal to the discontinuity can be determined by examining the acoustic tensor at the onset of inelastic deformation [8,20,21]. In this model no crack path continuity across element boundaries is not enforced since the actual length of discontinuity plays no part in the numerical result. For other types of localisation problems (such as shear band modelling with Von Mises plasticity), analysis of the acoustic tensor often leads to multiple roots. In this case, crack path continuity can be used to assist in choosing the appropriate discontinuity direction. 3.2. Energy dissipation and mesh objectivity Ð comparison with crack bandwidth models The similarity between the softening models for smeared crack bandwidth type formulations and strong discontinuity models (in Eqs. (31) and (33)) implies that in some cases the two approaches would yield the same result. For this reason the two approaches are compared for simple tension bar type problems to assess the objectivity of the approaches with respect to element size. The smeared model uses pa simple crack bandwidth (w ) calculated based on the element area. For the CST element, w ˆ 2Ae 3.2.1. Tension bar In Fig. 4 two bars of 20 CST elements are shown. One bar has a height of 1 mm and the second bar a height of 0.2 mm and the shaded elements are slightly weakened to induce localised softening. For both the smeared and embedded models Rankine plasticity with exponential softening is used (the softening law is given by Eq. (33)). The material properties adopted are shown in Table 1. The stress-displacement responses (stress is normalised using the tensile yield strength) for the four analysed cases are shown in Fig. 5, where F is the applied force, ft the tensile yield strength and A the cross-sectional area of the bar. It can be seen that the embedded model produces identical curves for

Fig. 4. Tension bar top (a) bottom (b).

G.N. Wells, L.J. Sluys / Engineering Fracture Mechanics 65 (2000) 263±281

273

Table 1 Material properties for tension bar E n fct Gf k

100 MPa 0.1 1.0 MPa 0.2 Nmm/mm2 1:0  10ÿ4 mm (for discontinuity model)

both bars while the smeared model does not. The energy dissipated in the failure of the thin …h ˆ 0:2 mm) bar is overestimated by the smeared model. This is due to the ambiguity in calculation of the crack bandwidth. Without taking into account the crack orientation within an element, the simple area-based regularisation procedure cannot produce unconditionally objective results. It is seen that calculation of a correct crack bandwidth (in the case of the bar with height = 1 mm) that the two methods yield exactly the same result. It can be concluded that it is not possible to calculate an objective regularisation parameter based simply on element size. The regularisation procedure lacks information about crack direction through the element, so can result in under- or over-estimation of the energy dissipated in failure, depending on the crack direction. The method is particularly not suitable for elements with large aspect ratios since the crack bandwidth becomes increasingly dependent on the direction of the crack through the element, relative to the element sides. The strong discontinuity model contains information over crack direction plus a prede®ned localisation zone width, giving energy dissipation independent of element size, aspect ratio and orientation. Remark. To avoid the erroneous result with respect to energy dissipation of the smeared model, more advanced crack bandwidth (characteristic length) calculation techniques can be applied [22]. However, more advanced crack bandwidth calculations do not o€er any improvements to the underlying kinematics of the ®nite element. As a result, the mesh alignment dependency of smeared models is still not addressed and

Fig. 5. Stress-displacement response of tension bar.

274

G.N. Wells, L.J. Sluys / Engineering Fracture Mechanics 65 (2000) 263±281

the constitutive model is not independent of the spatial discretisation.

3.3. Speci®c performance related issues In Section 2.4 it was shown that by specially constructing the softening parameter, continuum constitutive laws can be applied. The approximation of the Dirac-delta distribution in Eq. (34) allows the numerical application of continuum constitutive laws, where k can be considered as the width of an inner softening band. The strain ®eld in the softening band is given in terms of displacements at the discontinuity as ed ˆ

1 …nd ‰uŠ† k

…35†

The e€ect of the approximation across the discontinuity is here further investigated. In the following sections the e€ect of the displacement jump regularisation in combination with the Rankine plasticity model is examined. 3.3.1. Stress locking For modelling materials where mode-I type failure is dominant (such as concrete), Rankine plasticity has been applied in smeared crack modelling [23,24] and has been shown to largely overcome the spurious stress locking inherent in ®xed crack models [1]. The multi-surface plasticity model performs well in tension type problems where the principal stress axes rotate since the yield surface is constructed in the principal stress space, allowing rotation of the yield surface in the global coordinate system. However, for embedded models the introduction of a discontinuity immediately gives some directional bias which can result in stress locking. This e€ect for a Rankine model has been examined [17,25], with attention focused on stress locking as a consequence of an element being unable to kinematically produce the real deformation modes. The source of this problem and solution strategies are similar to those for smeared ®xed crack models [1], that is multiple and rotating discontinuities. The problem of the kinematic stress locking in embedded models with ®xed discontinuities is not unique to models using continuum constitutive laws. However, the application of a regularised discontinuity in combination with Rankine plasticity introduces a second source of stress locking that is dependent on the regularisation function, arising from the ®ctional nature of the strain ®eld at the discontinuity. The strain ®eld at the discontinuity is proportional to the displacements at the discontinuity (the objective measure of deformation) and inversely proportional to the width of the discontinuity band k (Eq. (35)). So, as k decreases (approaching the Dirac-delta distribution), the strains increase. This e€ect is then accounted for by the adjustment of the softening modulus. However, if the chosen constitutive law does not allow softening in all directions, the very large strains that develop at the discontinuity can result in spurious stress development. This is in contradiction to the condition at the beginning of Section 2.4 that the constitutive law must exhibit a softening branch. For Rankine plasticity, the response of a square patch in pure shear is examined (see Fig. 6) with the material parameters given in Table 2 and linear softening. Since the Rankine yield surface is in the principal stress space, there is no bound on the shear stresses that can develop. The shear stressdisplacement response in Fig. 6 shows that no shear softening occurs, hence the shear stresses that can develop are unbounded, with the response insensitive to the fracture energy parameter. In smeared analysis of tension dominated problems, the shear strains are usually small, hence large shear stresses do not develop and stress locking is not signi®cant.

G.N. Wells, L.J. Sluys / Engineering Fracture Mechanics 65 (2000) 263±281

275

Fig. 6. Shear panel and shear-displacement response for varying Gf :

As a result of the inability of the Rankine model to exhibit softening in shear, for decreasing k, very small sliding motions at the discontinuity result in very large shear stresses. In realistic type problems, the spurious shear stresses that develop at the discontinuity can be of the same order of magnitude as the inverse of k. This build up of very large spurious stresses can then lead to unpredictable results and a loss of robustness. A simple Mohr circle in Fig. 7 shows how the ®ctional strain state at the discontinuity can lead to the development of both large shear and compressive stresses for tension problems in the fully softened state. 3.3.2. Fixed/rotating discontinuity There are two possible solution methods for k-type stress locking with a regularised discontinuity. The ®rst, also common to kinematic type stress locking, is an improved constitutive model that allows both separation and sliding type deformations. The addition of a lower bound to the Rankine yield surfaces would limit the shear stress that could develop, although since it is still formed in the principal stress space this does not directly address the mode-II de®ciency of the Rankine model. Unlike for kinematic type stress locking, allowing multiple cracks is not a feasible solution for k-type stress locking. With decreasing k and Rankine plasticity, severe stress locking can occur with only very minor rotation of the principal axes. For a multiple crack model, as k decreases the threshold angle for the introduction of new discontinues to avoid stress locking approaches zero. For this reason, allowing rotation of discontinuities is examined. To compare ®xed and rotating discontinues using Rankine plasticity, a non-symmetric notched specimen, shown in Fig. 8 and previously examined by Berends et al. [26], is analysed. The material parameters adopted are shown in Table 2. The notched specimen is ®rst analysed using conventional continuum CST elements with the fracture energy regularised Rankine model for two di€erent meshes. The meshes are shown (in magni®ed Table 2 Material properties for notched specimen E n fct Gf k

10000 MPa 0.2 1.0 MPa 0.02 Nmm/mm2 1:0  10ÿ7 mm(for discontinuity model)

276

G.N. Wells, L.J. Sluys / Engineering Fracture Mechanics 65 (2000) 263±281

Fig. 7. Increasing compressive and shear stresses under tension loading in the fully softened state.

deformed con®guration) in Fig. 9. It can be seen from the load displacement response in Fig. 10 that the crack bandwidth model fails to produce mesh objective results when used with unstructured meshes. Further, the same two meshes are analysed using the ®xed embedded discontinuity model. Fig. 11 shows that the load-displacement response is not mesh objective and su€ers from severe stress locking. In extreme cases, hardening behaviour can be observed signi®cantly past the peak load. For the rotating discontinuity, the discontinuity direction is determined by the maximum principal stress direction at the beginning of the each load increment. By allowing rotation, at the beginning of each increment, no shear stresses develop on the discontinuity surface. Fig. 12 shows that by eliminating shear stresses on the discontinuity surface, the load displacement response is mesh objective. The discontinuity directions for the ®xed and rotating models at the end of the analysis are shown in Fig. 13 for mesh (a) in Fig. 9. In most elements the discontinuity directions are indistinguishable between the ®xed and rotating models, although it is evident from the load-displacement response that a signi®cant load is transfered across the ®xed discontinuities. The minor di€erence in discontinuity directions for the ®xed and rotating models excludes the possibility of alleviating the k-type stress locking by allowing multiple cracks per element.

Fig. 8. Notched specimen and deformed mesh.

G.N. Wells, L.J. Sluys / Engineering Fracture Mechanics 65 (2000) 263±281

277

Fig. 9. Deformed specimens for rotating model (magni®ed 400 x-times).

To further investigate the source of the k-type stress locking in the ®xed models, the sensitivity of the ®xed and rotating models to the regularisation function is examined. For mesh (a) in Fig. 9, the loaddisplacement response for the ®xed discontinuity is examined for varying k in Fig. 14. As k approaches zero, thereby approaching the Dirac-delta distribution, the stress locking e€ect increases. For small k, the failure of the constitutive law to allow shear softening becomes increasingly apparent. Allowing rotation of the discontinuity, the load-displacement response in Fig. 15 shows that the model is completely insensitive to the regularisation function (conditional on k being smaller than the element). Allowing small rotations of the discontinuity alleviates a de®ciency of the Rankine plasticity model when coupled with a regularised discontinuity. However, it is not a reasonable solution to kinematic type stress locking for the ®nite element model used here since large rotations of the discontinuity are necessary, making the model unstable as the discontinuity passes from one side of an element node to the other since the kinematically constructed shape functions change suddenly [17,10]. For applications in problems where signi®cant rotation of the principal stresses occurs, a combined rotating-multiple discontinuity approach is required to fully alleviate both k-type and kinematic stress locking.

Fig. 10. Load-displacement response for smeared model.

278

G.N. Wells, L.J. Sluys / Engineering Fracture Mechanics 65 (2000) 263±281

Fig. 11. Load-displacement response for ®xed discontinuity model.

4. Conclusions The embedment of discontinuities can be successfully applied to overcome energy dissipation mesh dependence of the classical continuum in softening problems and signi®cantly improve the objectiveness with respect to failure pattern. It has been shown to overcome the de®ciencies of crack bandwidth type models, particularly with respect to unstructured meshes and elements with large aspect ratios. By taking into account crack direction and using a pre-de®ned softening bandwidth, no mesh information is used in the constitutive softening law. Since the discontinuous modes are incompatible and continuum

Fig. 12. Load-displacement response for rotating discontinuity model.

G.N. Wells, L.J. Sluys / Engineering Fracture Mechanics 65 (2000) 263±281

279

Fig. 13. Embedded cracks for ®xed (a) and rotating (b) models for k ˆ 1:0  10ÿ7 :

constitutive laws are applied, the method can be easily implemented in existing ®nite element codes. The additional modes can be statically condensed, providing the magnitude of the displacement jump. No special modi®cations of strain softening classical continuum laws are required, allowing existing models based on fracture energy regularisation to be applied directly. However, care must still be taken when applying continuum relationships that spurious stress states to do not develop, as shown by the Rankine plasticity examples. The application of continuum constitutive laws in strong discontinuity analysis is based on the regularisation of the displacement jump, resulting in a bounded, although ®ctional, strain ®eld. In this respect embedded models are lacking when compared to higher-order models where the strain ®eld in the localised zone is an objective measure of deformation. However, the application of embedded discontinuities is intended for cases where the typical element size is signi®cantly larger than the identi®able internal length scale of the material. The examples using Rankine plasticity illustrate that a negative hardening modulus does not guarantee satisfaction of the softening branch condition for all constitutive laws, a problem highlighted by the

Fig. 14. Sensitivity to k for ®xed discontinuity.

280

G.N. Wells, L.J. Sluys / Engineering Fracture Mechanics 65 (2000) 263±281

Fig. 15. Sensitivity to k for rotating discontinuity.

®ctional strain ®eld at the discontinuity. The performance of di€erent constitutive laws in discontinuity analysis cannot be inferred directly from their behaviour in conventional continuum analysis and must be examined independently for the limiting case of a true strong discontinuity. Acknowledgements This research is supported by the Technology Foundation STW, applied science division of NWO and the technology program of the Ministry of Economic A€airs, The Netherlands. References [1] Rots JG. Computational modeling of concrete fracture. PhD thesis, Delft University of Technology, 1988. [2] Giessen E van der, Borst R de. Introduction to material instabilities in solids. In: de Borst R, van der Giessen E, editors. Material instabilities in solids. Delft, The Netherlands: Wiley, 1998. p. 1±13. [3] Aifantis EC. On the microstructural origin of certain inelastic models. Journal of Engineering Materials Technology 1984;106:326±34. [4] Borst RuÈde, Sluys LJ, MuÈhlhaus HB, Pamin J. Fundamental issues in ®nite element analyses of localization and deformation. Engineering Computations 1993;10:99±121. [5] Pijaudier-Cabot G, BazÏant Z. Nonlocal damage theory. ASCE Journal of Engineering Mechanics 1987;113(10):1512±33. [6] BazÏant Z, Pijaudier-Cabot G. Nonlocal continuum damage, localization instability and convergence. ASME Journal of Applied Mechanics 1988;55:287±93. [7] Belytschko T, Fish J, Engelmann BE. A ®nite element with embedded localization zones. Computer Methods in Applied Mechanics and Engineering 1988;70:59±89. [8] Sluys LJ, Berends AH. Discontinuous failure analysis for mode-I and mode-II localization problems. International Journal of Solids and Structures 1998;35(31-32):4257±74. [9] Sluys LJ. Discontinuous modelling of shear banding. In: Owen DRJ, OnÄate E, Hinton E, editors. Computational plasticity, fundamentals and applications. Barcelona, Spain: CIMNE, 1997. p. 735±44. [10] JiraÂsek M. Finite elements with embedded cracks. Technical Report LSC Internal Report 98/01, cole Polytechnique Fdrale de Lausanne, Switzerland, 1998.

G.N. Wells, L.J. Sluys / Engineering Fracture Mechanics 65 (2000) 263±281

281

[11] Simo JC, Oliver X, Armero F. An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids. Journal of Computational Mechanics 1993;12:277±96. [12] Oliver J, Simo JC. Modelling strong discontinuities by means of strain softening constitutive equations. In: Mang H, BicÂanic N, de Borst R, editors. EURO-C 1994 computer modelling of concrete structures. Innsbruck, Austria: Pineridge Press, 1994. p. 363±72. [13] Armero F, Garikipati K. Recent advances in the analysis and numerical simulation of strain localization in inelastic solids. In: Owen DRJ, OnÄate E, Hinton E, editors. Computational plasticity, fundamentals and applications. Barcelona, Spain: Pineridge Press, 1995. p. 547±61. [14] Lot® H, Shing PB. Embedded representation of fracture in concrete with mixed ®nite elements. International Journal for Numerical Methods in Engineering 1995;38:1307±25. [15] Larsson R, Runesson K. Element-embedded localization band based on regularized displacement discontinuity. ASCE Journal of Engineering Mechanics 1996;122(5):402±11. [16] Simo JC, Rifai MS. A class of mixed assumed strain methods and the method of incompatible modes. International Journal for Numerical Methods in Engineering 1990;29:1595±638. [17] Wells GN. Strong discontinuity analysis Ð the classical continuum. Technical Report LR-TM-99-01, Delft University of Technology, Delft, The Netherlands, 1999. [18] Taylor RL, Simo JC, Zienkiewicz OC, Chan CH. The patch test Ð a condition for assessing FEM convergence. International Journal for Numerical Methods in Engineering 1986;22:39±62. [19] Armero F, Garikipati K. An analysis of strong discontinuities in multiplicative ®nite strain plasticity and their relation with the numerical simulation of strain localization. International Journal of Solids and Structures 1996;33(20±22):2863±85. [20] Oliver J. Modelling strong discontinuities in solid mechanics via strain softening constitutive equations. Part I: fundamentals. International Journal for Numerical Methods in Engineering 1996;39:3575±600. [21] Ortiz M, Leroy Y, Needleman A. A ®nite element method for localized failure analysis. Computer Methods in Applied Mechanics and Engineering 1987;61:189±214. [22] Oliver J. A consistent characteristic length for smeared cracking models. International Journal for Numerical Methods in Engineering 1989;28:461±74. [23] Feenstra PH. Computational aspects of biaxial stress in plain and reinforced concrete. PhD thesis, Delft University of Technology, 1993. [24] Pearce CJ, BicÂanic N. One multi-surface plasticity and Rankine model. In: Owen DRJ, OnÄate E, editors. Computational plasticity, fundamentals and applications. Barcelona, Spain: Pineridge Press, 1997. p. 957±64. [25] Tano R, Klisinski M, Olofsson T. Stress locking in the inner softening band: A study of the origin and how to reduce the e€ects. In: de Borst R, BicÂanic N, Mang H, Meshke G, editors. EURO-C 1998 Computer Modelling of Concrete Structures. Badgasstein, Austria: Balkema, 1998. p. 363±72. [26] Berends AH, Sluys LJ, Borst R de. Discontinuous modelling of mode-I failure. In: Hendrix MA, Jongedijk H, Rots JG, Spanje WJE, editors. Finite elements in engineering and science. Amsterdam, The Netherlands: Balkema, 1997. p. 351±61.