Finite element variational formulation for beams with discontinuities

Finite Elements in Analysis and Design 54 (2012) 37–47

Contents lists available at SciVerse ScienceDirect

Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel

Finite element variational formulation for beams with discontinuities G. Juarez a,n, A.G. Ayala b,1 a b

´noma Metopolitana, Av. San Pablo 180, Col. Reynosa Tamaulipas, 02200, Me´xico D.F., Me´xico Departamento de Materiales, Universidad Auto ´noma de Me´xico, Ciudad Universitaria, 04510, Apdo. 70-642, Mexico D.F., Mexico Institute of Engineering, Universidad Nacional Auto

a r t i c l e i n f o

abstract

Article history: Received 29 December 2010 Received in revised form 31 December 2011 Accepted 2 January 2012 Available online 24 February 2012

This paper presents a variational formulation of the mechanical behaviour of beams with strong discontinuities, enhanced to simulate the strain localization process. The considered strain localization zones represent the formation of dislocations and hinges in beams. The presented general formulation applies to thick beams, which takes into account the internal strain energy due to bending and shear, and also a simpler formulation which takes into account only bending induced strains which applies to thin beams. It is shown that the developed energy functional for the beams with discontinuities has as stationarity conditions the strong formulation of the associated boundary value problem. As illustration, the energy functionals for Timoshenko and Euler–Bernoulli beams with embedded discontinuities are approximated by finite elements with embedded discontinuities. The development of a local material failure (leading to hinge-like strain localization zones) is in terms of continuum constitutive models furnished with strain softening capabilities. To show the validity of this formulation and its consistency with its continuum counterpart, representative numerical examples illustrating the performance of the proposed formulation are presented. & 2012 Elsevier B.V. All rights reserved.

Keywords: Bending elements Strain localization Embedded discontinuities Hinge development Dislocations

1. Introduction The process of material failure in solids is often preceded by a strain-localization phenomenon, characterized by the formation of strain localization zones in which damage and other inelastic effects concentrate. These zones gradually turn into physical discontinuities in the medium which, depending on the type of material, physically occur as: fissures in concrete, fractures in rocks and shear lines in soils and metals [1]. The idea of incorporating displacement or strain discontinuities into standard finite element interpolations to model regions with high localization of strains, called Embedded Discontinuities Model, has triggered the development of powerful techniques such as the Strong Discontinuity Approach which simulates material failure in the continuum during deformation represented as jumps in the displacement field. The solution of a problem with the Strong Discontinuity Approach, originally presented by Simo et al. [2] and Simo and Oliver [3], captures the jumps in the displacement field across a surface with zero bandwidth using standard solid mechanics models with continuum constitutive equations. This idea has induced the study and development of finite elements with

n

Corresponding author. Tel./fax:þ 55 5318 9085. E-mail addresses: [email protected], [email protected] (G. Juarez), [email protected] (A.G. Ayala). URL: http://www.azc.uam.mx/ (G. Juarez). 1 Tel.: þ55 5623 3508; fax: þ 55 5616 1514. 0168-874X/$ - see front matter & 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.finel.2012.01.004

embedded discontinuities (FEEDs) which capture the jump of the displacement field by additional degrees of freedom in the bulk of the element, nevertheless, some of these finite elements present problems such as mesh dependence, and stress locking [4,5]. Beam element is a particular case of a three-dimensional solid reduced, through a behaviour assumption, to a prismatic element along an axis. The study of beam elements is mainly divided into thick and thin beam theories, i.e., Timoshenko and Euler– Bernoulli beam theories. Thick elements consider the contribution of shear and bending internal strain energy whereas thin beams ignore the contribution of the shear strains [6]. For frame structures, formed with beams, the damage evolution process has been modelled [7] with the location of plastic hinges; this approximation, however, cannot simulate the process to collapse through softening. There are other approximations which allow the simulation of damage and the evolution to collapse by means of softening hinges of fixed length [8,18,9–11,12,20]. To study the failure process in beam elements, Ehrlich and Armero [13] recently developed a formulation for the analysis of localized failures, as softening dislocations and hinges, in Timoshenko beam elements. These dislocations and hinges are modelled as embedded discontinuities considering the jumps in the transverse displacement and in the rotation fields. These authors developed also a formulation for modelling softening hinges in thin Euler–Bernoulli beams [14] and an extended formulation for the analysis of softening hinge lines in inelastic thick plates [15]. Their mathematical model is formulated using a

38

G. Juarez, A.G. Ayala / Finite Elements in Analysis and Design 54 (2012) 37–47

displacement based energy functional, which satisfies equilibrium in a weak form, and a strong equation for the inner moment and shear continuity. The finite element approximation of this mathematical model includes an enhanced strain operator to avoid locking. This formulation shows an acceptable performance in the numerical simulation of the formation of hinges in plates and beams. Nevertheless, the fact that in this approximation the inner moment and shear continuities are imposed in a strong form leads to non-symmetric stiffness matrices. The objective of this work is to develop a consistent variational formulation of beam structural members with embedded discontinuities (leading to dislocation and hinge-like strain localization zones) to simulate the material failure problem in beams by finite element approximations. Subsequently, the variational formulations for Timoshenko and Euler–Bernoulli beams are developed as particular cases of the general variational formulation. Even though a complete hierarchy of finite element formulations is derived from the presented general formulation, to demonstrate its validity only the Timoshenko and Euler–Bernoulli finite element approximations for beams with discontinuities are implemented in this work. The approximations of this functional by the Finite Element Method (FEM) lead naturally to a FEEDs formulation able to capture the discontinuities and to dissipate the energy due to damage in a consistent way. The resulting finite element matrices of this formulation are symmetric; the stability and convergence of the numerical solutions is guaranteed. The outline of the rest of this paper is as follows. Section 2 presents the kinematics and the boundary value problem (BVP) of beam elements with discontinuities. Section 3 provides the constitutive models to describe the material behaviour of the material in the continuum, at the localization zone and a damage model for the development of hinges. Section 4 presents the development of a general energy functional of beams with strong discontinuities and a consequential formulation of a hierarchy of energy functionals. Section 5 shows the FEM approximation of Timoshenko and Bernoulli variational formulation with strong discontinuities. Some numerical examples of bending elements with discontinuities which validate the proposed formulation are presented in Section 6. Finally, in Section 7 some conclusions derived from the work and some suggestions about future research are given.

2. Beam members theory A beam is a structural member bounded by two ends, called edge or boundary. The beam theory considers a prismatic onedimensional body represented by its neutral axis, with an open bounded domain, O A R1 , material points, x, and boundary, G (Fig. 1). The transverse load per unit length, q, is given on O, the natural and essential boundary conditions are the transverse shear force, V n , on GV , the bending moment, Mn , on GM , the transverse n displacement, wn , on Gw , and the rotation, y , on, Gy , respectively, such that GV [ GM ¼ Gw [ Gy ¼ G and GV \ Gw ¼ GM \ Gy ¼ |.

Fig. 2. Kinematics of a beam with a strain localization zone in S: (a) graphic description, (b) rotation, (c) transverse displacement, (d) curvature and (e) shear strain.

2.1. Kinematics of discontinuous displacement fields 2.1.1. Thick beams with discontinuities Consider the beam shown in Fig. 2a loaded until it undergoes a transverse displacement jump, ½9w9, and/or a rotation jump, ½9y9, across a localization zone, S, where the domain splits into subdomains, O ¼ O þ O þ , and their corresponding ends G ¼ G þ G þ . The localization zone, S, is characterized by the concentration of inelastic strains in a point, which starts with the formation of voids that progressively turn into a macroscopic discontinuities, while the surrounding material undergoes unloading. To describe the kinematics of a beam, which presents rotations and transverse displacements discontinuities, consider that the inelastic strains concentrate into a localization zone of zero width, S, and that the discontinuous rotation, yðxÞ, and transverse displacement, wðxÞ, fields illustrated in Fig. 2b and c are given by

yðxÞ ¼ y þ Hs ½9y9

ð1Þ

wðxÞ ¼ w þHs ½9w9

ð2Þ

where y and w are the continuous part of the rotation and the transverse displacement, respectively, with the jumps ½9y9 and ½9w9 at a given (material) point S, inducing unbounded curvature, ky , and shear strain, gw, y , fields. Both fields can be expressed as (Fig. 2d and e)

ky ¼

@y @y ¼ þ ds ½9y9 @x |{z} @x |fflfflffl{zfflfflffl} ½9y9 ky

gw, y ¼

k~

    @w @w y ¼ y þ ds ½9w9Hs ½9y9 @x @x |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl} ½9w9,½9y9 g w ,y

Fig. 1. Beam with boundary conditions on G.

ð3Þ

ð4Þ

g~

where HS is the Heaviside function defined on S (HS ðxÞ ¼ 0 8x A O and HS ðxÞ ¼ 1 8x A O þ ) and ds is the Dirac delta function. The continuous curvature, k y , and the continuous shear strain, g w, y , are defined in O\S, whereas the curvature, k~ ½9y9 , and the shear strain, g~ ½9w9,½9y9 , are defined on O. In general, the development of hinges and/or dislocations in thick beam elements, may be

G. Juarez, A.G. Ayala / Finite Elements in Analysis and Design 54 (2012) 37–47

39

Eqs. (7b) correspond to the constitutive compatibility between the moment–curvature and the shear force–shear strain relations. The internal equilibrium equations, Eq. (7c), are defined by the following relations, V¼

2

dM , dx

dV d M q ¼ q ¼ 0 2 dx dx

In Eqs. (7c), qz includes the contribution of the transverse load, p, and the body force bz .

Fig. 3. Jumps into a discontinuity: (a) hinge, (b) dislocation and (c) both.

qz ¼ ðpz þ tbz Þ modelled as: a rotation jump, ½9h9, a transverse displacement jump, ½9w9, or both, see Fig. 3. 2.1.2. Thin beam with discontinuities In thin beams, Euler–Bernoulli elements, the evolution of hinges is simulated only as rotation jumps, the moment field is continuous and bounded in the domain of the beam, whereas the curvature is unbounded in the localization zone and bounded in the rest of the element. When a hinge is completely developed, the element is unable to transmit moments in the localization zone. In this approach, the rotation field with a jump, ½9y9, in the localization zone, S, is defined as

yðxÞ ¼ y ðxÞ þ Hs ðxÞ½9y9 ¼ y^ ðxÞ þ cS ðxÞ½9y9 ¼

^ dw þ cS ðxÞ½9y9 dx

ð5Þ

where y and y^ are the continuous and regular rotation fields, respectively. The curvature according to Eqs. (3) and (5) is given by

k¼

^ @jðxÞ½9y9 @y @2 w @2 w ¼ þ ds ½9y9 ¼ 2  þ ds ½9y9 2 @x |ffl{zffl} @x ffl} |fflfflffl{zfflfflffl} @x @x |fflfflffl{zfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl ½9y9 ½9y9 kw

k~

kw

ð8Þ

ð9Þ

3. Constitutive models in bending dominated elements 3.1. Continuum models The constitutive models describe the physical properties of a given material. The constitutive compatibility of beam elements needs two laws, the moment–curvature and the shear force–shear strain relations as shown in Fig. 4. In these constitutive models, the area under the moment function, MðkÞ, corresponds to the moment free energy density CM ðkÞ, given by Z k CM ðkÞ ¼ MðkÞ dk ð10Þ 0

and the area above the curvature function, kðMÞ, to the complementary moment free energy density Z M CM kðMÞ dM ð11Þ c ðMÞ ¼ 0

ð6Þ

k~ b

The area under the shear function, VðgÞ, is defined as the shear free energy density Z g CV ðgÞ ¼ VðgÞ dg ð12Þ 0

2.2. Boundary value problems The BVP for a thick beam with discontinuities is defined by the following equations and boundary conditions:

0

k y ðx,tÞk ðx,tÞ ¼ 0 g w, y ðx,tÞg ðx,tÞ ¼ 0 M k ðx,tÞMðx,tÞ ¼ 0 V g ðx,tÞVðx,tÞ ¼ 0 Mðx,tÞ, x Vðx,tÞ ¼ 0 Vðx,tÞx qz ðx,tÞ ¼ 0

in O=S kinematical compatibility

ð7aÞ

The complementary moment energy density and complementary shear energy are defined respectively by the Legendre transformation as

in O=S constitutive compatibility

ð7bÞ

M CM c ðMÞ ¼ M  kC ðkÞ

ð14Þ

CVc ðVÞ ¼ V  gCV ðgÞ

ð15Þ

in O=S internal equilibrium

Mðx,tÞ  nM n ðx,tÞ ¼ 0

on GM

Vðx,tÞ  nV n ðx,tÞ ¼ 0

on GV

Mðx,tÞ  nMe ðx,tÞ ¼ 0

on Gy

Vðx,tÞ  nV e ðx,tÞ ¼ 0

on Gw

yðx,tÞ ¼ yn ðx,tÞ wðx,tÞ ¼ wn ðx,tÞ

and the area above the shear strain, gðVÞ, to the complementary shear free energy density Z V CVc ðVÞ ¼ gðVÞ dV ð13Þ

on Gy on Gw

ðV O V S Þ ¼ ðV S V O þ Þ ¼ 0

The moment is derived from Eq. (10) as Mk ðkÞ ¼

external equilibrium

essential boundary condition

ðM O M S Þ ¼ ðMS M O þ Þ ¼ 0

ð7cÞ

on S inner continuity

@CM ðkÞ @k

ð7dÞ

ð7eÞ

ð7fÞ

Eqs. (7a) correspond to the kinematical compatibility between the curvature, ky , and the shear strain, cw, y , with the rotations, y, and the transverse displacement, w, given respectively in Eqs. (3) and (4).

Fig. 4. Constitutive models: (a) Mk and (b) V g.

ð16Þ

40

G. Juarez, A.G. Ayala / Finite Elements in Analysis and Design 54 (2012) 37–47

and the curvature from Eq. (11):

kM ðMÞ ¼

@CM c ðMÞ @M

ð17Þ

In the same way, the shear force is derived from Eq. (12) as V g ðgÞ ¼

@CV ðgÞ @g

ð18Þ

and the shear strain from Eq. (13):

kV ðgÞ ¼

@CVc ðgÞ @V

ð19Þ

In any variational formulation, involving the strain energy density, once k and g are computed, the constitutive model provides in an elastic domain, the moment M, the shear force V, and the tangent constitutive operator C T from a tangent constitutive equation in terms of rates: _ ¼ C T : k_ M M

ð20Þ

V_ ¼ C TV : g_

ð21Þ

Fig. 6. Force damage model: (a) continuum and (b) discrete.

In the same way, the shear force is derived from Eq. (25) as V S ð½9w9Þ ¼

@CV S ð½9w9Þ @½9w9

ð27Þ

3.3. Damage model _ is the rate of moment, V_ of the shear force, k_ of the where M curvature and g_ of the shear strain. Otherwise, in a variational formulation, involving the complementary strain energy density, once the M and V are computed, the constitutive model provides in an elastic domain, the curvature k, the shear strain g, and the tangent compliance constitutive tensor DT from a tangent constitutive equation in terms of rates:

k_ ¼ DTM : M_

ð22Þ

g_ ¼ DTV : V_

ð23Þ

In the Embedded Discontinuities Model, the constitutive model is one of the important ingredients, because it describes the material behaviour. In a continuum approach, the behaviour of the strain localization zone is modelled by a standard continuum constitutive equation (Fig. 6a), a moment/transverse shear forces, F , and curvature/shear strains, w, relationship. However, in a Discrete Approach (Fig. 6b) the behaviour of the discontinuity is described by a moment/shear force, F , and rotation/transverse displacement jump, ½9w9, relationship. In this paper, the discrete damage model, defined in Eq. (28), is used. The equations describing this model are

jð½9w9, a Þ ¼ ð1oÞj0 ða Þ, fj0 ð½9w9Þg ¼ 12 ½9w9  C  ½9w9 jð½9w9, a Þ ¼ ð1oÞC  ½9w9 Constitutive equation F ¼ @½9w9 qða Þ Damage variable o ¼ 1 ; o A ½1,1 a __ @ Evolution law a ¼ l ¼ ða Þ, a A ½0,1 Free energy

3.2. Discrete models In addition to the continuum models, the beam elements with discontinuities need two laws in the localization zone, the moment–rotation jump and the shear force–transverse displacement jump relations as shown in Fig. 5. In these constitutive models, the area under the moment function, MS ð½9y9Þ, corresponds to the moment free energy density CMS ð½9y9Þ, given by Z ½9y9 CMS ð½9y9Þ ¼ M S ð½9y9Þ dk ð24Þ 0

The area under the shear function, V S ð½9w9Þ, is defined as the shear free energy density Z ½9w9 CV S ðgÞ ¼ V S ð½9w9Þ d½9w9 ð25Þ 0

The moment is derived from Eq. (24) as M S ð½9y9Þ ¼

@CMS ð½9y9Þ @½9y9

ð26Þ

@t

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi F  C1  F qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi t½9w9 ¼ J½9w9J ¼ ½9w9  C  ½9w9

Damage criterion f ðF ,qÞ ¼ tT q; gð½9w9, a Þ ¼ t½9w9 a ; __

__

Hardening rule qða Þ ¼ H a;

tT ¼ JF JQ e1 ¼

H ¼ q 0 ða Þ r0

Loading2unloading conditions f r 0; l f_ ¼ 0 ðconsistencyÞ

l Z 0;

l f ¼ 0; ð28Þ

where j is the discrete free energy, C is the elastic tensor, F is the force vector containing the moment, M S , and transverse shear force, V S . The damage variable o is defined in terms of a __ hardening/softening variable q, which is dependent on the hardening/softening parameter H: The damage multiplier, l , defines the loading–unloading conditions, the function f ðF , qÞ bounds the elastic domain. The value, r 0 , is the threshold that determines the initial inelastic behaviour. The vector, ½9w9, contains the rotation and transverse displacement jumps, ½9y9 and ½9w9. From the models defined in Eq. (28), it is possible to obtain the tangent constitutive equation in terms of rates of the variables as F_ ¼ CT  ½9w_ 9

ð29Þ

where the tangent constitutive operator for the nonlinear loading range is defined by Fig. 5. Constitutive models: (a) MS ½9y9 and (b) V S ½9w9.

CT ¼ ð1oÞC

qH a

a3

ðC  ½9w9  ½9w9  CÞ

ð30Þ

G. Juarez, A.G. Ayala / Finite Elements in Analysis and Design 54 (2012) 37–47

_ ¼ 0Þ: and for the elastic loading and unloading range ðo T

C ¼ ð1oÞC

ð31Þ

41

displacement–strain energy functional with discontinuities, Z Pðw, y, k, g,½9w9,½9y9Þ ¼ ½M k : k y þ V g  g w, y CM ðk ÞCV ðg Þqz  w dO O\S

þ CMS ð½9y9Þ þ CV S ð½9w9ÞM n  y9G V n  w9G M k  ðyy Þ9G V  ðwwn Þ9G n

4. Variational formulations of beam with discontinuities 4.1. General energy functional Multiplying the BVP given in Eq. (7) by their corresponding Lagrange multipliers, lY , and integrating over an interval: Z " O\S



@M V þ ðV VÞ  lg þ @x



þ CMS ð½9y9Þ þ CV S ð½9w9ÞM n  y9G V n  w9G n

M  ðyy Þ9G V  ðwwn Þ9G

4.3. Total potential energy functional

O\S

þ ðV g VÞ  dg þ M  dk y þV  dg w, y qz  dw dO þ ðMM n Þ  dy þ ðVV n Þ  dwþ ðMMe Þ  dy þ ðVV e Þ  dw ½9y9

þ ðM S

½9w9

MÞ  d½9y9 þðV S

VÞ  d½9w9 ¼ 0

Considering in Eq. (33) the divergence theorem Z Z @M  dy dO þ M  dyM  d½9y9 M  dk y dO ¼ O\S O\S @x Z V  dg w, y ¼ 

V  dy dO þ O\S

Z

ð33Þ

ð34Þ

@V  dw dO þ V  dwV  d½9w9 O\S @x

Substituting Eqs. (34) and (35) in Eq. (33) Z  dP ¼ ðk y k Þ  dM þ ðg w, y g Þ  dV þðM k MÞ  dk O\S      @M @V V  dy þ qz  dw dO þ ðV g VÞ  dg þ @x @x 

 d½9w9 ¼ 0

Z

ð40Þ

The second functional is derived neglecting shear deformations, the continuous rotations, y , are dependent on the continuous transverse displacement, w, by the strong relationship (5). Then, from Eq. (6), the curvature, k , is now dependent on the continuous transverse displacement, w, @2 w @x2

ð41Þ

Thus, the energy functional defined in Eq. (40), with four independent variables, reduces to two independent variables, corresponding to the Euler–Bernoulli theory with discontinuities.

Pðw,½9h9Þ ¼ ð36Þ

Eq. (36) corresponds to the first variation of a general energy functional, with eight independent fields given by

Pðw, y,M,V, k, g,½9w9,½9y9Þ ¼

þ CMS ð½9y9Þ þ CV S ð½9w9ÞMn  y9G V n  w9G

kw ¼

þ ðMM n Þ  dy þðVV n Þ  dw þðMM e Þ  dy þ ðVV e Þ  dw ½9w9 d½9y9 þ ðV S VÞ

Other two energy functionals for beams with discontinuities may be derived from the general energy functional. The first is derived by assuming that the kinematical compatibility in Eq. (7a) is a priori satisfied in the energy functional defined in Eq. (37). As result, the corresponding version of Timoshenko functional for beams with discontinuities yields: Z Pðw, y,½9w9,½9y9Þ ¼ ½CM ðk y Þ þ CV ðg w, y Þqz  w dO O\S

ð35Þ

½9y9 þ ðMS MÞ

ð39Þ

ð32Þ

The meaning of internal energy to each term of Eq. (32) suggests identifying lY with the variation of the field Y, then Z dP ¼ ½ðk y k Þ  dM þðg w, y g Þ  dV þ ðM k MÞ  dk

O\S

O\S

   @V  ly þ qz lw dO @x

þ ðMM n Þ  ly þðVV n Þ  lw þðMMÞ  ly þ ðVVÞ  lw þ ðM S MÞ  l½9y9 þ ðV S VÞ  l½9w9 ¼ 0

Z

The Hellinger–Reissner energy functional for beams with discontinuities may also be derived as a particular case of the general energy functional by substituting Eqs. (14) and (15) into Eq. (37), yielding: Z V Pðw, y,M,V,½9w9,½9y9Þ ¼ ½M : k y þ V  gw, y CM c ðMÞCc ðVÞqz  w dO

ðk y k Þ  lM þðg w, y g Þ  lV þ ðM k MÞ  lk g

ð38Þ

½M  ðk y k Þ þ V  ðg w, y g Þ þ CM ðk Þ þ CV ðg Þ dO

Z O\S

½CM ðk w Þqz  w dO þ CMS ð½9y9ÞM n  y9G V n  w9G

ð42Þ

5. Finite element approximation of beams with discontinuities

O\S

Z  O\S

qz  w dO þ CMS ð½9y9Þ þ CV S ð½9w9Þ

M n  h9G V n  w9G M  ðyy Þ9G V  ðwwn Þ9G n

ð37Þ The forces M S and V S are the moment and shear on the localization zone S; these forces, dependent on the jumps ½9y9 and ½9w9, are calculated from a discrete damage model defined in Eq. (28). 4.2. Mixed variational formulation From the general energy functional of beams with discontinuities a hierarchy of energy functionals may be derived, including the mixed and potential energy formulations. Suppose that constitutive compatibility in Eq. (7b) is a priori satisfied. Then, the energy functional defined in Eq. (37) reduces to an energy functional with six independent variables, corresponding to a

In this section, the variational formulations of the Timoshenko and the Euler–Bernoulli beams with discontinuities, developed in the previous section, are approximated by the FEM. 5.1. Finite element approximation of beams 5.1.1. Discretization Let a beam under loading that has a localization zone, S (Fig. 7a), be discretized by beam finite elements, as shown in Fig. 7b. The discontinuities in the localization zone, S, are simulated into the finite element mesh as a dislocation, a hinge line or both by introducing an internal node whose degrees of freedom incorporate the jumps ½9y9 and ½9w9 (Fig. 3). 5.1.1.1. Regularization of the displacement kinematics. Considering that it is not possible to prescribe the boundary conditions, wn n and y , on the transverse displacement and rotation fields in only

42

G. Juarez, A.G. Ayala / Finite Elements in Analysis and Design 54 (2012) 37–47

y ðx,tÞ ¼ y^ ðx,tÞjðxÞ½9y9ðx,tÞ

ð50Þ

Substituting Eqs. (49) and (50) into the definition of continuous curvatures, k y and k w , and the continuous shear strain, g w, y given in Eqs. (3), (4) and (41) respectively, leads to Fig. 7. Beam with a localization zone S: (a) geometric description and boundary conditions and (b) discretization with finite elements.

k y ðx,tÞ ¼

@y^ ðx,tÞ @jðxÞ  ½9y9ðx,tÞ @x @x

ð51Þ

k w ðx,tÞ ¼

^ @2 wðx,tÞ @2 jðxÞ  ½9y9ðx,tÞ 2 @x2 @x

ð52Þ

g w, y ¼

^ @wðx,tÞ @jðxÞ  ½9w9ðx,tÞy^ ðx,tÞ þ jðxÞ½9y9ðx,tÞ @x @x

ð53Þ

5.1.1.2. Approximation of transverse displacement, rotation, curvature and shear strain fields. The regular transverse displacement and rotation fields are approximated by ^ ¼ Nw wd w

ð54Þ

y^ ¼ N h yd

ð55Þ

where Nw and Nh are respectively the standard vector of interpolation functions of the nodal transverse displacement, wd , and rotation, yd , vectors. The function, cS ðxÞ, is defined in the finite element approach as

ceS ðxÞ ¼ HeS ðxÞje

ð56Þ

e

where j is constructed by

je ¼ Ni þ

Fig. 8. Graphic representation of (a) transverse displacement, (b) rotation and (c) function M S .

one part of the field as defined in Eqs. (1) and (2), i.e., in the continuous part, w and y , or in the jumps, ½9y9 and ½9w9, these fields may be described by the following equations (Fig. 8a and b): ^ wðx,tÞ ¼ wðx,tÞ þ cS ½9w9ðx,tÞ

ð43Þ

yðx,tÞ ¼ y^ ðx,tÞ þ cS ½9y9ðx,tÞ

ð44Þ

Then, the curvature and the shear strains, dependent on the regular transverse displacement and rotation, are defined by

ky ðx,tÞ ¼

@½9y9ðx,tÞ @yðx,tÞ @y^ ðx,tÞ @cS ¼ þ ½9y9ðx,tÞ þ cS @x @x @x @x

ð45Þ

^ @wðx,tÞ @cS þ ½9w9ðx,tÞy^ ðx,tÞcS ½9y9ðx,tÞ ð46Þ @x @x ^ where wðx,tÞ and y^ ðx,tÞ are the regular transverse displacement and rotation fields, and cS ðxÞ is a function given by

gw, y ¼

cS ðxÞ ¼ HS ðxÞjðxÞ

e

wðx,tÞ ¼ Nw wd þ cS ½9w9

ð58Þ

yðx,tÞ ¼ Nh yd þ cS ½9y9

ð59Þ

The continuous curvature and shear strain fields in Eqs. (51)–(53), dependent on the transverse displacement and rotation, are approximated as

ky ¼

@Nh @je ½9y9 y  @x d @x

ð60Þ

kw ¼

@2 N w @2 je ðxÞ wd  ½9y9 @x2 @x2

ð61Þ

g w, y ¼

@N w @ je w  ½9w9Nh yd þ je ½9y9 @x d @x

ð62Þ

5.2. Finite element approximation of beams 5.2.1. Timoshenko beam with discontinuities The finite element approximation of the regular transverse displacement and rotation is given by

jðxÞ ¼ 0 8x A O ð48Þ

The function, cS , has two properties: cS ðxÞ ¼ 1 8x A S and cS ðxÞ ¼ 0 8x A O [ O þ (Fig. 8c). The continuous transverse displacement and rotation field of Eqs. (1) and (2) are ^ wðx,tÞ ¼ wðx,tÞ jðxÞ½9w9ðx,tÞ

where Ni þ are the interpolation functions of the nodes placed on O þ of the finite element which contains the discontinuity, in agreement with the definition of j in Eq. (48). The transverse displacement and rotation fields defined in Eqs. (43) and (44), dependent on the regular parts and jumps, are given by

ð47Þ

where jðxÞ is a continuous function such that

jðxÞ ¼ 1 8x A O þ

ð57Þ

ð49Þ

y ðxÞ ¼ Ny yj½9y9 wðxÞ ¼ Nw wj½9w9

ð63Þ

Assuming a linear interpolation in both fields Nw ¼ N y ¼ ½N1 N2 

with

N1 ¼ Lx L N2 ¼ xL

ð64Þ

G. Juarez, A.G. Ayala / Finite Elements in Analysis and Design 54 (2012) 37–47

43

and j given by

j¼

x L

ð65Þ

Substituting Eq. (63) into the functional given in Eq. (32), setting its derivatives with respect to each variable to zero, and after its linearization using Taylor series gives  2 3 a a a  aL  a2   a2 2 L 6 La 7  a a a 6 L 7  a2   aL L 2 2 6 7  6 a 7 a EI a L EI a L a EI a L L þ 6  6 2 7 2 2 L þ 3 L  6 6 7  6  a  a  EI þ aL 7 EI a L  a EI a L þ  6 2 7 L 6 L 3  L 3 2 2 6 7 6 7  a a a a a þ IC T a 6 7   2 s L L 2 2 L 6 7  4 5 a a EI a L EI a L T a 1 a 2 2 þ þ IC L  6 L  3  2 b L 3



8 > > > > > > > > > > <

_1 w _2 w y_ 1

9 > > > > > > > > > > =

8 F1 > > > > > F2 > > > > < M1

9 > > > > > > > > > =

y_ 2 > ¼ > M2 > > > > > > > > > > > > > > > > > _ > > > > > ½9 w9 > > > > 0 > > > > > > > > > ; : 0 ; : ½9y_ 9 >

ð66Þ

5.2.2. Bernoulli–Euler beam with discontinuities For Euler–Bernoulli beam, the regular transverse displacement is approximated as ^ wðxÞ ¼ Nd

ð67Þ

where Nd from the classical beam theory stands for the approximated transverse displacement, w, by the following equation: wðxÞ ¼ N 1 w1 þ N2 y1 þ N3 w2 þN 4 y2 ¼ Nd

N2 ¼

3

L

1 L3

ð2x3 3Lx2 þ L3 Þ,

ðLx3 2L2 x2 þ L3 xÞ:

12EI L3

6 6 6EI 6 2 6 L 6 6 12EI 6 3 L 6 6 4 6EI L2

6EI 2

L   4 E þ EIDT L EIL b



 6EI L2



2þ E þEIC T L b

 12EI L3  6EI L2 12EI L3

EI L

 6EI L2



3

6EI L2

 78 _ 9 8 9 7> w 1 > > F 1 > > > > > 2þ E þ EIDT L EIL 7 > > > > _ = < < 7 M1 = b 7 y1 ¼ 7 _ > > F > w 7>  6EI > L2 > > > > 2> > 2 >   7 7: y_ 2 ; : M 2 ; 5 4 E þEIC T L EIL b

N3 ¼

1 3

L

N4 ¼

1

ðx3 LL2 x2 Þ

which corresponds to the stiffness matrix for a beam element with a hinge along the span [16], Fig. 9. Thus, this formulation represents accurately the onset of a strain localization zone until the completed formation of a hinge, where the value of the bending moment is zero. 6. Numerical examples To show the validity and performance of the formulation presented in this paper, this section presents some representative numerical examples.

ð2x3 þ 3Lx2 Þ

L3

Suppose that the element with the previous stiffness matrix has been completely damaged, so the tangential discrete material stiffness, C Tb , is null, i.e., there is no transmission of bending moment in the localization zone. It follows that the stiffness matrix defined by Eq. (71) changes and the equation becomes 2 3 12EI 6EI 6EI 9 8 _ 9 8  12EI L3 L2 L3 L2 1> > 6 7> w > F1 > > > 3EI 6EI 3EI 7> > > > 6 6EI _ = < y 1 = < M1 >  L2 6 L2 L L 7 6 7 ¼ ð72Þ 12EI 6EI 7> w 6  12EI  6EI _ F > >  2> > > 2 > 6 L3 L2 L3 L2 7> > ; > ; : _ > : 5> 4 6EI 3EI 3EI y2 M2  6EI L L L2 L2

ð68Þ

with the interpolation functions defined by 1

2

ð71Þ

As it may be seen, in this equation the stiffness matrix of a Timoshenko beam finite element with embedded discontinuities is symmetric. The additional degrees of freedom in this equation can be statically condensed, facilitating its implementation.

N1 ¼

Fig. 9. Beam element with a hinge.

ð69Þ

Substituting Eqs. (67) and (65) into the functional of Eq. (42), setting the variations with respect to each variable are zero, and after its linearization with Taylor series also yields the following stiffness matrix, in incremental form:  2 12EI 3 6EI 6EI   12EI 8 _ 9 8 9 0 L3 L2 L3 L2  w1 > F1 > 6 6EI 7> > > > 4EI 6EI 2EI  EI 6 7 > > > >  L2 _ > > > >  L L > > > 6 L2 7 y1 > L > > > M1 >  > > > > 6 12EI 7 = = < < 6EI 12EI 6EI  6 3  2 7 0  _ 6 7 w2 F L L L3 L2  2 ¼ ð70Þ EI 6 6EI 7 > > > > y_ > >  L 7> 2EI 4EI  6  6EI 2 > > > M2 >  L L > > > > 6 L2 7> L2 > > > > > : > > > 6 7> ;  ; 4 5: ½9y_ 9 > 0 T EI 0 EIL 0  EIL  þIC b L It is of interest to mention that this variational formulation also yields symmetric stiffness matrices for the Euler–Bernoulli beam. To reduce the additional degree of freedom, the resulting stiffness matrix can be statically condensed as follows:

6.1. Simple supported beam Consider the simply supported beam shown in Fig. 10, loaded at the centre of the span, with length, L ¼ 1, and elastic properties, EI ¼ 1. In the simulation of this example, Euler–Bernoulli beam finite elements were used. The behaviour of the material is initially linear elastic until the bending moment reaches its yield value M y ¼ 1. At this instant, a hinge-like strain localization zone develops, and a rotation jump, ½9y9, takes place at the centre of the beam. The bending moment softens following a linear softening law, H ¼ 0:8: Fig. 11 shows the variation of load, P, vs. transverse displacement, w, and variation of the moment, M, vs. rotation jump, ½9y9, at the centre of the span, L. As shown, the load, P, and the moment, M, decrease after it reaches its yield value, M y , and the beam evolves into collapse. Integrating both curves shows that the energy supplied to the system by the load, 0.625 N m, is equal to the energy dissipated in the discontinuity. This shows that the performance of the elements to simulate this problem is correct.

44

G. Juarez, A.G. Ayala / Finite Elements in Analysis and Design 54 (2012) 37–47

Fig. 10. Beam loaded at the centre: (a) geometric description and (b) with a localization zone (hinge).

Fig. 11. Variation: (a) load, P, vs. transverse displacement, w, and (b) moment, M, vs. jump, ½9y9.

Fig. 12. Cantilever beam: (a) geometric description and (b) with a localization zone (hinge).

Fig. 13. Variation: (a) load, Ry2 , vs. imposed displacement, w2, (b) shear, V S , vs. transverse displacement jump, ½9w9 and (c) moment, M S , vs. rotation jump, ½9h9.

6.2. Cantilever beam Consider the cantilever beam with length, L¼1, shown in Fig. 12a subjected to an imposed transversal displacement, w2, at the free end, the elastic properties are EI¼1 and GAs¼0.641. To study this problem Timoshenko beam finite elements were used. In this example, the behaviour of the material is initially linear elastic until the shear force reaches its yield value, V y ¼ 1, and consequently, the bending moment reaches its maximum value at the fixed end of the beam, where the dislocation-like strain localization zone occurs. The bending moment softens following a linear softening law, H ¼ 0:1. Fig. 13 shows the variation of the reaction, Ry2 , vs. the imposed displacement, w 2 , at the free end of the beam and the variation of the shear, V S , vs. the transverse displacement jump, ½9w9, at the end. As shown, the reaction, Ry2 , and the moment, M S , decrease after the yield value is reached, M y , and the beam evolves into

collapse. Integrating the area under the three curves shows that the energy supplied to the system in Fig. 13a, 5.0 N m, is equal to the energy dissipated during the development of the discontinuity in Fig. 13b and c. In this particular case, the kind of discontinuity developed is a dislocation, so the transverse shear force, V S , presents softening and the moment, MS , unloads in an elastic way. 6.3. Fixed beam Consider the beam with fixed ends shown in Fig. 14a, loaded at the centre of the span, with length, L¼2 m, and elastic properties, E¼40,000 MPa. The beam is 0.20 m thick and 0.10 m wide. The behaviour of the material is initially linear elastic until the bending moment reaches a yield value, My ¼0.8 MN m. In this case, the localization zone (hinge) takes place at the centre of the span (Fig. 14b). The discrete softening modulus is H¼  3.0. In the simulation of this example, Euler–Bernoulli beam finite elements were used.

G. Juarez, A.G. Ayala / Finite Elements in Analysis and Design 54 (2012) 37–47

45

Fig. 14. Fixed beam: (a) geometric description and (b) with a localization zone (hinge).

Fig. 15. Variation: (a) load, P, vs.displacement, w, and (b) moment, M, vs. rotation jumps, ½9y9.

Fig. 16. Simple supported beam: (a) geometric description, (b) idealization with plane FE, (c) idealization with beam FE and (d) with a localization zone (hinge).

Fig. 15 shows the variation of the load, P, vs. the displacement, w, at the span centre and the variation of the moment, M, vs. the rotation jump, ½9h9. As can be seen, once the yield value of the bending moment is reached, the maximum displacement takes place, w¼0.05 m; at this instant, the development of a hinge initiates, while the surrounding material of the beam undergoes unloading. When the development of the hinge is completed, there is a maximum displacement at the centre of the span, w¼0.025 m. The hatched area under the curves in Fig. 15 corresponds to the energy spent in the development of the hinge.

6.4. Simple supported beam with plane stress approximation A simply supported beam with a depth of 1 mm is loaded at the centre of the span, see Fig. 16. This example was analysed using two-dimensional six-node plane stress triangular finite elements by Wells and Sluys [17]. The following material properties were used: Young’s modulus, E¼100 MPa, Poisson’s ratio, n ¼ 0:0, elastic limit moment, Mu ¼ 2:5 N mm, and fracture energy

Gf ¼ 0:1 J mm2 . To study this example Euler–Bernoulli and Timoshenko beam finite elements were used. Fig. 17 shows the variation of the load, P, vs. transverse displacement, u, at the centre of the span, L. The differences in the slopes of ascending part of the curves is because in the beam approximation, the beginning of damage occurs when the value of Mu is reached, whereas in the plane stress approximation, damage occurs before that value is reached, fact that reduces its slope. Nevertheless, integrating the load–displacement response for the three curves in Fig. 17 shows that the energy dissipated is equal to 0.3080 N mm for plane elements [17] and 0.30 N mm for beam elements with discontinuities. This agrees well with the fracture energy multiplied by the depth of the beam that shows the performance or the elements to simulate this problem. 6.5. Darvall–Mendis frame Fig. 18a shows the same clamped portal frame with a vertical load applied at a distance of 0.55L to the top left corner used as example by [18,14]. The one level frame is 3 m high and 3 m wide.

46

G. Juarez, A.G. Ayala / Finite Elements in Analysis and Design 54 (2012) 37–47

The columns and the beam have a cross-section area, A¼ 0.10 m2, and moment of inertia, I ¼ 1  103 m4 . The material properties of the elements were: Young’s modulus, E¼20.7 GPa, yield moments for the columns, M ycol ¼ 158 kN m, and for the beam, M ybe ¼ 169 kN m. Four values of the discrete softening modulus were considered, H¼ 0.0 (elastoplastic behaviour), H¼  0.053EI, H¼ 0.081EI and H¼  0.098EI MN m/rad for softening. An elastic response under axial force was assumed for all elements. To reproduce the exact analytical solution of this example, the frame was discretized using four Euler–Bernoulli beam finite elements, Fig. 18b. A monotonically increasing displacement, d, was imposed at node 3 until it reaches a magnitude of d ¼ 4:75 mm, i.e., the first hinge was formed at this node. As expected, this hinge occurred in all cases considered at the same d value. For the elastoplastic case, i.e., H¼0.0, when the imposed displacement reached a value d ¼ 11:35 a new hinge was formed at the upper part of the right column. The simulation continues until a new hinge at the upper part of the left column occurred, i.e., an imposed displacement d ¼ 13:42. For the three cases with softening behaviour when the second hinge was formed at the upper node of the right column, the load reached its maximum value. Table 1 shows the ultimate load, Pu , and its corresponding imposed displacement, d, for different values of the discrete softening modulus. Fig. 19 shows the variation of the moment, M ybe , and the rotation, y, of the hinge formed at node 3. As seen in Fig. 20 snap-back responses occurred for all cases with softening hinges, this behaviour is consistent with the theoretical description given by [19]. It is important to mention that the results obtained with this formulation up to the ultimate

load match exactly, without any modification, the results of [18] and that if the analysis is continued considering displacement control, this formulation simulates the progression to collapse correctly. For comparison sake, in the paper by [14] to get the same results up to the maximum load, it is necessary to increase the number of elements or enrich the definition of the curvature of the element, both somehow inconsistent with the needed finite element approximation of the problem.

Fig. 17. Variation load, P, vs. displacement, u.

Fig. 19. Variation of moment, M ybe , vs. rotation, y.

7. Conclusions This paper presented the development and application of a general variational formulation for beams with strong discontinuities adapted to represent strain localization in bending-dominated problems in structural mechanics. The main conclusions are

 It is shown that it is possible to develop, in a natural way, a general variational formulation for beam members with discontinuities based on energy sense of the Lagrange multipliers. Table 1 Displacement and ultimate load. H (kN m/rad)

d (mm)

P u (kN)

0  1.089  1.681  2.017

13.4 11.8 12.2 12.8

434.32 388.02 358.04 337.10

Fig. 18. Frame: (a) geometric description, (b) discretization with four elements and (c) discretization with six elements.

G. Juarez, A.G. Ayala / Finite Elements in Analysis and Design 54 (2012) 37–47

47

Fig. 20. Variation load, P, vs. displacement, u.

 The developed general variational formulation leads to a









hierarchy of particular formulations associated to mixed, displacement and other formulations encountered within the framework of finite element approximations. The mathematical model of each beam structural member was developed in variational form, using an energy functional, whose stationarity conditions provide the field equations which govern the problem (strong formulation). However, the use of a variational form allows a deeper mathematical treatment of questions of existence, stability and convergence of its numerical approximations. The advantage of this formulation is that the explored finite element matrices, Timoshenko and Euler–Bernoulli beams with discontinuities, are symmetric, do not present shear locking problems despite of the use of piecewise constant jumps. The numerical implementation of this formulation permits the correct simulation of structural engineering problems involving beam structural members with discontinuities until collapse is reached. It is evident that the application of this methodology to real problems incorporating shear deformation needs a failure criterion involving moments and the shear forces. The development of this criterion falls out of the scope of this paper.

Acknowledgements The first author acknowledges the support given by the Universidad Auto´noma Metropolitana (UAM). The authors acknowledge the support given by the Institute of Engineering of the Universidad Nacional Auto´noma de Mexico (UNAM). References [1] Z.P. Bazant, J. Planas, Fracture and Size Effect in Concrete and Other Quasibrittle Materials, CRC Press, 1998.

[2] J.C. Simo, J. Oliver, F. Armero, An analysis of strong discontinuities induced by strain-softening in rate-independent inelastic solids, Comput. Mech. 12 (1993) 277–296. [3] J.C. Simo, J. Oliver, A new approach to the analysis and simulation of strain softening in solids, in: Z.P. Bazant, et al. (Eds.), Fracture and Damage in Quasibrittle Structures, E. and F.N. Spon, London, 1994, pp. 25–39. [4] M. Jira´sek, Comparative study on finite elements with embedded discontinuities, Comput. Methods Appl. Mech. Eng. 188 (2000) 307–330. [5] J. Oliver, A.E. Huespe, E. Samaniego, A study on finite elements for capturing strong discontinuities, Int. J. Numer. Methods Eng. 56 (2003) 2135–2161. [6] C.A. Felippa, The Amusing History of Shear Flexible Beam Elements, Report CU-CAS-05-1, University of Colorado, Boulder, CO, 2005. [7] J.F. Baker, J. Hyman, Plastic Design of Frames, vols. 1 and 2, Cambridge University Press, 1969. [8] G. Maier, A. Zavelani, J.C. Dotreppe, Equilibrium branching due to flexural softening, ASCE J. Eng. Mech. 4 (1973) 897–901. [9] Z.P. Bazant, M.T. Kazemi, Localization of softening damage in frames and implications for earthquake resistance, in: Proceedings of the 5th US National Conference on Earthquake Engineering, vol. 1, Earthquake Engineering Research Institute, Oakland, CA, 1994, pp. 313–322. [10] J.G. Sanjayan, P.L. Darvall, Singularities in RC beam elements with finitelength hinges, ASCE J. Struct. Eng. 121 (1995) 39–47. [11] Z.P. Bazant, M. Jirasek, Softening-induced dynamic localization instability: seismic damage in frames, ASCE J. Eng. Mech. 122 (1996) 1149–1158. [12] M. Jira´sek, Analytical and numerical solutions for frames with softening hinges, ASCE J. Eng. Mech. 123 (1997) 8–14. [13] D. Ehrlich, F. Armero, Finite element methods for the analysis of softening plastic hinges in beams and frames, Comput. Mech. 35 (2004) 237–264. [14] F. Armero, D. Ehrlich, Numerical modeling of softening hinges in thin Euler– Bernoulli beams, Comput. Struct. 84 (2006) 641–656. [15] F. Armero, D. Ehrlich, Finite element methods for the multi-scale modelling of softening hinge lines in plates at failure, Comput. Methods Appl. Mech. Eng. 195 (2006) 1283–1324. [16] C.A. Felippa, Introduction to Finite Element Methods, Course Notes, /http:// caswww.colorado.edu/Felippa.d/FelippaHome.d/Home.htmlS, 2004. [17] G.N. Wells, L.J. Sluys, A new method for modelling cohesive cracks using finite elements, Int. J. Numer. Methods Eng. 50 (2001) 2667–2682. [18] P.L. Darvall, P.A. Mendis, Elastic–plastic-softening analysis of plane fames, ASCE J. Struct. Eng. 111 (1984) 871–888. [19] Z.P. Bazant, Asymptotic matching analysis of scaling structural failure due to softening hinges. I: theory. II: implications, ASCE J. Eng. Mech. 129 (2003). 641–650, 651–654. [20] M.E. Marante, R. Pico´n, J. Florez-Lopez, Analysis of localization in frame members with plastic hinges, Int. J. Solids Struct. 41 (2004) 3961–3975.