Numerical modeling of softening hinges in thin Euler–Bernoulli beams

Computers and Structures 84 (2006) 641–656 www.elsevier.com/locate/compstruc

Numerical modeling of softening hinges in thin Euler–Bernoulli beams F. Armero *, D. Ehrlich Structural Engineering, Mechanics and Materials, Department of Civil and Environmental Engineering, University of California, Berkeley, CA 94720-1710, USA Received 21 December 2004; accepted 21 November 2005

Abstract This paper presents new finite elements for thin Euler–Bernoulli beams that incorporate the softening hinges observed at failure. The proposed methods rely crucially on the identification of the classical notion of inelastic hinge with strong discontinuities of the generalized displacements describing the beam’s deformation. The development of a multi-scale framework that effectively incorporates the localized dissipative mechanisms associated with these discontinuous solutions into the large-scale problem of the beam, or general frame system, defines a crucial step undertaken here. This framework defines the setting of its numerical implementation by finite elements enhanced with the singular strains corresponding to the discontinuities. A general procedure is presented that leads, in particular, to finite elements free of stress-locking and that resolve exactly the kinematics of the hinge. Several numerical simulations are presented to illustrate the performance of the proposed methods.
2005 Elsevier Ltd. All rights reserved. Keywords: Euler–Bernoulli beams; Localized failures of frames; Softening hinges; Enhanced finite elements

1. Introduction The failure of beams and frames is usually accompanied by the localization of the yielding and damage in particular critical cross-sections. This situation leads naturally to the concept of plastic hinges where, in the case of bending, very high values of the curvature of the beam’s axis are observed in a very small region, leading in the limit to a relative rotation between the two segments of the beam. The consideration of these failure mechanisms is the basis of the classical limit analysis of structures, where particular sections are assumed to have reached its limit capacity, kept constant while other plastic hinges develop. A good and brief overview of these ideas can be found in [20]. Classical limit analysis does not account, however, for the post-peak structural response, which eventually results in a decreasing load capacity with increasing deformation. The need of considering softening hinges to capture this overall softening structural response was soon realized *

Corresponding author. E-mail address: [email protected] (F. Armero).

0045-7949/$ - see front matter
2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2005.11.010

for some important cases: namely, prestressed and overreinforced concrete beams due to the compressive failure of the concrete, as well as steel beams when local buckling occurs. We can refer to [21,22,4,10,32,7,16], among many of the early works with an analytical description of the phenomenon. A crucial aspect added by the consideration of softening hinges is the size effect that appears in the prediction of failure. This situation is also characteristic of other localized failures in materials (e.g. classical linear elastic fracture). We refer to [5,6] for a complete discussion on all of these aspects. The introduction of softening hinges leads to an involved load redistribution in redundant structures (see e.g. [6,8]), requiring then the development of numerical techniques, like finite element methods, for its resolution. In this setting, the hinge is often modeled with a short segment of the beam exhibiting large localized strains (e.g. curvature) with a relation between these strains and the conjugate stress resultants (e.g. bending moment–curvature relation) exhibiting strain softening. Representative examples of this approach are the work of Darvall and Mendis [11] based on a fixed hinge length, usually given in terms of

642

F. Armero, D. Ehrlich / Computers and Structures 84 (2006) 641–656

the beam’s depth, or the regularization based on a characteristic size of the finite element mesh used in the analysis, as in [9], generalizing similar classical developments for continuum problems. We note that the two approaches accomplish the same objective by the proper scaling of the moment curvature law since, even in the case that a material characteristic length exists, this length is much smaller that the scales resolved by the considered meshes. Analyses of finite elements with finite length hinges can be found in [31]. However, the analysis presented recently in [3] clearly indicates that the elastoplastic problem of a Timoshenko beam, accounting for all the different couplings between axial and shear forces and bending moments in the inelastic range, is mathematically ill-posed in the presence of strain softening if no viscous effects are present. This situation was also corroborated in this reference by the finding of the exact solution of the particular problem of wave propagation in a simply supported beam at failure and the realization of its unphysical character if a strain softening moment–curvature is considered: no energy is dissipated despite the failure of the beam. The fact that these difficulties appear even in a Timoshenko beam model, which possesses a length scale (the beam depth), is significant. The same difficulties are in fact observed in an Euler–Bernoulli beam model, the case of interest here, where that length scale is not present. On the other hand, the consideration of a softening hinge involving a discontinuous rotation, with a localized softening law between the bending moment and rotation jump, was shown in this last reference to lead to the regularization of the mathematical problem and to an exact solution for that particular problem that recovers its physical significance. Importantly, the plastic hinges are understood as strong discontinuities in the beam problem, that is, as solutions involving discontinuities of the generalized displacements describing the beam’s deformation (axial displacement, transversal deflection and cross-section rotation). The analysis and numerical simulation of discontinuities in solids have received a large amount of attention in the recent literature. We refer to [34,26,27,24,38], for representative references, and [29], and references therein, for the related approach based on cohesive elements. Here we follow the so-called strong discontinuity approach along the lines presented in [1,2] for continuum problems. This particular treatment defines a multi-scale framework for the incorporation of the dissipative effects associated to the localized failures into the large-scale problem with the usual regularity requirement on its solutions. The strong discontinuity becomes a tool for the modeling of the localized dissipative effects associated with these failures in the small scales of the material, incorporating the localized dissipation through an enhanced local strain field depending on the discontinuity jumps. These jumps are eventually expressed in terms of the large-scale (smooth) displacements in the large-scale limit of vanishing small scales,

leading then to an enhanced large-scale model of the structure capturing efficiently and objectively the dissipation associated with the localized failures. These considerations lead to a natural and efficient numerical implementation in the context of enhanced finite element methods. We develop these ideas for Euler–Bernoulli thin beams and frames in this paper. Finite elements with plastic hinges lumped at the element’s ends can be found in [18,23], among others. The strong discontinuity approach allows, however, to develop finite elements that fully incorporate the kinematics of the softening hinge at any position. It is precisely the main goal of this paper to develop new such enhanced finite elements for the case of thin beams, that is, in the framework of Euler–Bernoulli beam finite elements exhibiting C1 continuity. The different nature of the kinematics in this setting when compared with the kinematics of a Timoshenko beam considered in [12] motivates the separate treatment presented in this paper. Note that the strain measures in an Euler–Bernoulli beam involve second derivatives of the generalized displacements, a situation never encountered in previous problems approached through the strong discontinuity approach. In particular, we arrive to new finite elements involving Hermite type interpolations in the large scale, that exactly capture the kinematics of the hinges at the element level and that avoid the so-called stress locking (that is, a numerically spurious transfer of stresses through fully softened hinges). An outline of the rest of the paper is as follows. Section 2 develops the aforementioned multi-scale framework for the modeling of strong discontinuities (hinges) in Euler– Bernoulli thin beams. The discussion is relatively brief in order to focus the attention on the development of the enhanced finite elements for thin beams incorporating the softening hinges. This is carried out in Section 3 where a general methodology is developed for the design of these elements, including some specific model finite elements. Section 4 presents several representative numerical simulations that illustrate the ability of the proposed framework of capturing the localized failures of beam and frame structural systems, evaluating at the same time the performance of the newly proposed finite elements. Finally, Section 5 concludes with some final remarks. 2. Multi-scale modeling of softening hinges in thin-beams We present in this section the proposed multi-scale framework of inelastic Euler–Bernoulli beams at failure. To this purpose, the basic well-known equations describing the kinematics and equilibrium of such thin beam are summarized in Section 2.1. These equations define the so-called large-scale problem of interest. Section 2.2 considers the introduction of strong discontinuities in this setting, leading to the classical notion of plastic hinges. These solutions are used in the modeling of the dissipative effects associated to the small scales and resulting in the localized failure of the beams.

F. Armero, D. Ehrlich / Computers and Structures 84 (2006) 641–656

643

2.1. The Euler–Bernoulli beam problem

ance of the localized dissipative mechanism has been detected. We consider then the local displacement fields

We consider the case of a straight Euler–Bernoulli beam, modeled by its axis B ¼ ½0; L
 R. Its infinitesimal plane deformations are described by the axis transversal deflection w : B ! R and the longitudinal displacement u : B ! R, with u :¼ {u, w}. We denote by x 2 [0, L] the local coordinate system along the beam’s axis, so the cross-section rotation is given by w 0 :¼ dw/dx (the Kirchhoff assumption). As shown by the numerical examples presented in Section 4, the case of interest may involve frames consisting of beams with different orientations. The standard considerations defining this local coordinate system x apply then in the assembly of the different equations corresponding to each member; details are omitted for brevity. In the infinitesimal range of interest, the characteristic strain measures are given by the axial strain ea = u 0 and the curvature j = w00 . We use the short-hand notation e(u) :¼ {ea(u), j(w)}. Similarly, the conjugate stress resultants are denoted by r = {N, M} for the axial force N and bending moment M. The equilibrium of the beam is then characterized by the weak statement (principle of virtual work) Z Z rT eðduÞ dx ¼ ½
n du þ
q dw
dx þ N dujoN B

ul ¼ u þ nu W u

B

B

þ V dwjoV B þ M dw0 joM B

ð2:1Þ

for all admissible variations du (that is, du = 0 on ou B, dw = 0 on ow B and dw 0 = 0 on ow0 B, corresponding to the parts of the boundary oB with imposed axial displace
and rotation w0 ¼ w
0 , respecment u ¼
u, deflection w ¼ w tively). In (2.1), the external distributed loading is denoted by
n and
q, with boundary axial N and transversal V loads and moment M. Eq. (2.1) defines the large-scale problem in the multiscale framework considered in this work. In this way, it is posed in terms of the large-scale displacements u = {u, w} with the standard regularity conditions, that is, u 2 H 1 ðBÞ  H 2 ðBÞ involving, in particular, continuous longitudinal and transversal displacements and rotation, and thus regular distributions for the strains e (as opposed to singular distributions; see e.g. [36]). This regularity is characteristic of elastic beams but it breaks down when one considers plastic models with strain softening (i.e. softening r(e) relation) since this leads to an ill-posed problem; see [3]. The interest in this work is to develop the proper framework that incorporates discontinuities in the displacements and the corresponding singular strains with the purpose to model the localized dissipative mechanisms associated with the small scale response of the material at failure in the classical form of softening hinges. The resulting small-scale problem is developed next. 2.2. Softening hinges as strong discontinuities Following the discussion above, we consider a local neighborhood Bxd  B at a point xd 2 B where the appear-

and

wl ¼ w þ nh Ww

in Bxd

ð2:2Þ

for two parameters n :¼ {nu, nh} and two local functions Wu ; Ww : Bxd ! R. These local displacement distributions allow to model the localized dissipative effects of the hinge developing at xd. In particular, we consider local displacement functions satisfying the jump conditions sWu t ¼ 1

and

sW0w t ¼ 1

at xd 2 B

ð2:3Þ

so we identify the local parameters n with the longitudinal displacement jump nu = sulb and the rotation jump nh ¼ sw0l t at the hinge. In this way, the strains el :¼ fu0l ; w00l g associated with the small-scale displacements (2.2) (that is, the small-scale strains) can be written as el ¼ eðuÞ þ Gn þndxd |fflfflfflfflfflffl{zfflfflfflfflfflffl}

ð2:4Þ


e

for a local strain distribution G and the Dirac delta function dxd associated with the discontinuities (2.3). We observe that nu defines the axial opening of the hinge, whereas nh describes the relative rotation of the two parts of the beam about the hinge in the thin-beam limit of interest in this paper. We refer to [12] for a discussion of the kinematics of hinges in a Timoshenko beam, including also a jump in the transversal deflection. This would account for a localized shear strain, but such strain is not present in the Euler–Bernoulli beam considered here. As noted in Section 1, the goal of the approach developed here is to capture the energy dissipation associated with the localized dissipative mechanism of the strong discontinuity in Bxd by the large-scale problem in B and in the large-scale limit of vanishing small scales, that is, hxd ¼ measureðBxd Þ ! 0. This is accomplished by assuring that the stress power in the large scale (that is, in terms of the large-scale displacements u) is the same as the stress power in the small scale (that is, for the small-scale displacements ul). To impose this condition, we first note the relation "Z # Z Z T T T T _ r dx þ n_ eðuÞ G r dx þ rðxd Þ e_ r dx ¼ l

Bx d

Bx d

Bxd

ð2:5Þ for a stress resultant r distribution in Bxd and its value r(xd) at xd, after using the properties of the Dirac delta function in (2.4). Therefore, we conclude that the bracket in Eq. (2.5) must vanish, hence obtaining the equation Z G T r dx þ rðxd Þ ¼ 0 ð2:6Þ Bx d

in the large-scale limit of vanishing small scales hxd ! 0. The consideration of the generalized displacement jumps n allows to develop rigorously constitutive models of the hinge response in the form of rd(n) for the associated

644

F. Armero, D. Ehrlich / Computers and Structures 84 (2006) 641–656

driving stress resultant rd :¼ {Nd, Md}, as described in the next section. We refer to the resulting constitutive model of the hinge as the localized model. It remains to connect the driving stress resultant rd with the original stress resultant r in Bxd . The latter is given by the regular part of the strain
e (2.4) in the bulk Bxdnxd , through the proper constitutive relation, as also discussed in the next section. Indeed this connection must supply enough equations to solve for the newly introduced parameters n. A natural choice of this definition is the identification of the driving stress rd with r(xd), that is, the actual value of the stress resultant at the hinge. Noting the relation (2.6), we can write Z G T rðu; nÞ dx þ rd ðnÞ ¼ 0; ð2:7Þ Bx d

which defines the sought equation. We refer to this relation as ‘‘local equilibrium’’. Eq. (2.7) allows to solve for the local parameters n in terms of the large-scale displacements u in the large-scale limit hxd ! 0, thus eliminating these local parameters from the large-scale problem (2.1) and, hence, arriving to an ‘‘enhanced large-scale problem’’ (i.e. formulated entirely in terms of the regular large-scale displacement u) that captures the localized dissipation associated to the small scales. We refer to [1,2] for a complete discussion in the context of continuum problems. We emphasize that all these considerations are understood in the limit of vanishing small scales hxd ! 0. The consideration of this limit allows, in fact, the modeling of the localized failures of interest (i.e. hinges) very efficiently without a detailed definition of the small scales Bxd and the local distribution functions Wu and Ww. We only observe the requirement Z G dx ¼ 1 as hxd ! 0 ð2:8Þ Bx d

for the identity matrix 1, a relation obtained easily from (2.5). We observe that (2.8) implies that G  1=hxd þ oð1Þ as hxd ! 0. We also note that the operators G in the kinematic relation (2.4) defining the small-scale strains and in the local equilibrium relation (2.7) can be in principle different if they comply with this requirement in the limit hxd ! 0. In fact, this is exploited in the numerical implementation of these ideas developed in Section 3. 2.3. The constitutive models The above relations need to be complemented by the proper constitutive models for the bulk and localized responses of the beam and the hinge, respectively. The setting considered above fits well the so-called plastic models in stress resultants for the constitutive modeling of the bulk response of the beam. In these models, the stress resultants r = {N, M} are given directly in terms of the strain measures e = {ea, j}. We refer to [25,37,28,30,20,13–15], among others, for complete details on a number of such models of steel, reinforced concrete and composite beams, together

with their detailed calibration and validation. Of particular interest is the calibration of these models through the so-called fiber models considering plasticity developing through the cross-section, as presented in [13,15]. See also the recent review article [35] for applications to steel–concrete composite structures. As shown in these references, the models in stress resultants can indeed capture partially plastified sections and their evolution through the proper choice of yield surfaces and hardening laws, including involved two-surface kinematic/isotropic hardening laws; see [14] for a discussion and comparison of different alternatives proposed to date. For the purpose of this paper (namely the novel kinematic treatment of softening hinges developed in the previous section in its theoretical part), it suffices to consider an elastoplastic model in stress resultants given by 9 8 o/ > > p > > = < e_ ¼ k or p r ¼ Cðe  e Þ with o/ > > > > ; : a_ ¼ k oq and   / 6 0; k P 0 ð2:9Þ k/ ¼ 0; k/_ ¼ 0 for a yield surface (or interaction diagram) /(r, q) in the stress resultants space following a hardening law q(a). As usual, ðÞ_ denotes time rate. The general relations (2.9) follow the classical structure of plastic model where the strain are additively decomposed in an elastic and plastic part; see [33]. A linear elastic relation has been assumed in (2.9) which is typically given by C ¼ DIAG½EA; EI
for the Young modulus E of the material and the beam’s crosssection area A and inertia I. The analysis presented in [3] identifies the ill-posedness of the mathematical problems based on inviscid stress resultant models like (2.9) when combined with strain softening laws (that is, H :¼ oq=oa < 0), requiring in this case the introduction of the strong discontinuities considered in the previous sections. This situation even happens in the context of Timoshenko beams where the constitutive model incorporates a length scale (the beam depth). In this way, we consider this softening response as triggering the appearance of the strong discontinuities, which follow then a localized inelastic model between the generalized displacement jumps n and the driving stress resultants rd. The consideration of this pair of quantities allow to develop the localized models in a fully thermodynamic consistent framework. As an example, we can have the relations similar to the bulk model (2.9), that is, ~ o/ n_ ¼ ~k ord

and

~ o/ a~_ ¼ ~k o~q

ð2:10Þ

~ d; q ~Þ based on an yield surface (or interaction diagram) /ðr function of the stress resultant rd driving the hinge and the

F. Armero, D. Ehrlich / Computers and Structures 84 (2006) 641–656

localized internal variables ~ q and ~ a, related by a given soft~ ening cohesive law ~ q¼^ qð~ aÞ. The localized plastic multiplier ~ k is given by the Kuhn–Tucker loading/unloading and consistency conditions ~ 6 0; k ~/ ~ ¼ 0; and k ~/ ~_ ¼ 0 k~ P 0; / ð2:11Þ resulting in rigid-plastic localized relation between the jump n and the driving stress resultants rd. The notation ~ has been employed to distinguish between the quantities ðÞ in the localized model (2.10) and the constitutive model in the bulk (2.9). We refer again to [1,2] for a complete discussion of localized models in the continuum, including the consideration of elastic, damage and finite deformation effects. The goal here is to develop finite element methods that incorporate the considerations developed above, namely, finite elements that exactly capture the kinematics of the hinges identified in Section 2.2 for thin beams and frames. 3. Enhanced finite elements A clear advantage of the constructive style used in the development of the multi-scale framework presented in the previous section, going through the consideration of the small scales in the form of the local neighborhoods Bxd , is its direct translation in the final finite element formulation. In this context, we identify the finite elements with these neighborhoods Be ¼ Bxd , with the strains (2.4) defining a local enhanced strain field at the element level. The large-scale limit he ¼ measureðBe Þ ! 0 has then a clear significance, obtained it in the limit of finite element mesh refinement. However, the actual computations are carried out with a finite mesh he > 0 thus requiring the construction of the displacement functions (2.3) or, equivalently, the local strain operator G in (2.4). The construction of the appropriate enhanced operators is developed in Section 3.3 after presenting in Section 3.1 general considerations about the finite element formulation. 3.1. The finite element formulation WeS consider a spatial discretization of the beam’s axis nelem B ¼ e¼1 Be in nelem finite elements Be ¼ ½xe1 ; xe2
of size he = xe2  xe1 > 0, not necessarily equal. The interpolated large-scale displacements uh = {uh, wh} are given in terms of a standard conforming interpolation uh ¼ Nd ) eðuh Þ ¼ Bd

ð3:1Þ

for a standard set of shape functions N and associated nodal degrees of freedom d, with the corresponding strain operator B defining the large-scale strains. As usual in the Euler–Bernoulli framework under consideration, a C1 interpolation is required for the transversal deflection. The local enhanced strains (2.4) are written then as ehl ¼ Bd þ G c nh þnh dxd |fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl}
eh

ð3:2Þ

645

for a local (element) set of enhanced parameters nhe and a local enhanced strain operator Gc (‘‘c’’ for compatibility as defining the enhanced strains). No regularization of the Dirac delta function is considered in the formulation proposed here. Standard arguments lead to the finite element discretization of the large-scale equation (2.1) in B and the smallscale equation (2.7) locally in each element Beloc where a discontinuity has been activated. This results in the set of non-linear algebraic equations 9 nelem R h > > R ¼ f ext  A 0 e B T rðd e ; ne Þ ds ¼ 0; = e¼1 Z he ð3:3Þ > ; G Te rðd e ; ne Þ ds þ rd ðne Þ ¼ 0 in Be;loc > reenh ¼ 0

in the global degrees of freedom d and the local enhanced parameters ne, where A denotes the standard finite element assembly and fext the external nodal forces. Here, and for future use, we have employed the local element coordinate s = x  xe1, with sd = xd  xe1 representing the point in the element Beloc where a hinge has formed. In Eq. (3.3)2 we have used an enhanced operator Ge (‘‘e’’ for equilibrium) different in general than the strain operator Gc defining the enhanced strains in (3.2). As noted in Section 2.2, this is possible for a finite size he > 0 if the limit requirement (2.8) is maintained. A simple choice is given by the constant enhanced operators G ðconstÞ ¼ G ðconstÞ ¼ e c

1 1. he

ð3:4Þ

However, the separate and non-constant treatment of the Gc and Ge leads to an improved resolution of the kinematics and the equilibrium enforcement for a finite size mesh he > 0. We note that the use of different enhanced operators may lead to a non-symmetric linearized problem in general. Different but constant enhanced strain operators has been successfully employed in the modeling of strong discontinuities in continuum problems; see [1,2] and references therein. As shown by the developments presented below, a better choice is required for the considered beam problem given its more intricate kinematics, involving second derivatives in the strain definition (i.e. the curvature) and discontinuities in the first derivative of the displacement (i.e. the rotation). 3.2. The enhanced strain equilibrium operator The improved performance can be obtained, in particular, by considering an equilibrium enhanced operator, denoted by G ðexactÞ , such that e Z he T G ðexactÞ ðsÞrðsÞ ds ¼ rðsd Þ ð3:5Þ e 0

for the assumed interpolation of the stress resultants r(s) in the finite element. This is accomplished by a diagonal operator

646

G ðexactÞ e

F. Armero, D. Ehrlich / Computers and Structures 84 (2006) 641–656

" 1 ge;
a ðsÞ ¼ he 0

0

# ð3:6Þ

ge;j ðsÞ

for scalar polynomial functions, denoted generically by ge(s), that satisfy Z he 1 g ðsÞrðsÞ ds ¼ rðsd Þ ð3:7Þ he 0 e for the particular distribution r(s) of each stress resultant in the finite element. For example, for a linear distribution, like the bending moment of the two-node element considered in Section 3.3.1, we have    2sd 2s ð1Þ ge ðsÞ ¼ 1 þ 3 1  1 ð3:8Þ he he after noting that Z he 1 gð1Þ ðsÞ ds ¼ 1 he 0 e

and

1 he

Z

This observation provides a clear criterion for the design of the enhanced strains. We first identify the nodal degrees of freedom d^ hinge corresponding to the hinge mode. This can be written in general as d^ hinge ¼ d^ rigid þ Dhinge ^ne

for the nodal displacements corresponding to a general rigid deformation of the beam element d^ rigid and a geometric operator Dhinge defining the kinematics of the hinge mode in terms of the opening jumps ^ne . The operator Dhinge depends in general on the geometry of the particular base finite element under consideration and the local position for the hinge sd . Fig. 1 illustrates these ideas for a two-node Hermite finite element. The requirement that the regular part of the strain must vanish for the hinge mode leads to the relation ^ne ¼ 0
eh ¼ Bd^ hinge þ G ðexactÞ c

he

sgeð1Þ ðsÞ ds ¼ sd

ð3:9Þ

0

as a straightforward calculation shows. For a constant distribution, like the axial force in that same finite element, the use of (3.8) reduces to the constant function gð0Þ e ðsÞ ¼ 1 as in (3.4). Given (3.9)1, the enhanced strain operator G ðexactÞ defined by (3.8) trivially satisfies the e asymptotic requirement (2.8). Furthermore, we observe that with the choice (3.5), Eq. (3.3)2 enforces exactly the local equilibrium in the sense that the driving stress resultant on the hinge is now given exactly by the stress resultant at the hinge (i.e. rd = r(sd)) for any mesh size he > 0. 3.3. The enhanced strain compatibility operator The correct design of the enhanced strain compatibility operator Gc requires a full understanding of the hinge kinematics. In particular, the finite element should allow for a zero regular strain
eh in (3.2) for any value of the jumps ne in order to model the final state of a fully softened hinge (rd = 0). We refer to the final deformation of the beam as the ‘‘hinge mode’’. Failure to do so may result in spurious transfer of stress resultants across the hinge, a phenomenon usually referred to as stress-locking; see [12] for a discussion in the context of Timoshenko beams.

ð3:10Þ

8^ne ) G cðexactÞ ¼ BDhinge ð3:11Þ

after noting that no strains occur for the rigid part of the deformation, that is, B d^ rigid ¼ 0 for the strain operator B of the base finite element. We observe that the resulting enhanced strain operator is not constant in general. The reason for the use of the name ‘‘(exact)’’ will become apparent in the next section. 3.3.1. An enhanced two-node Hermite finite element We particularize the above developments to the case of a two-node finite element as a model example. The longitudinal displacement is defined by the piece-wise linear interpolation   s s uh ðsÞ ¼ 1  ð3:12Þ u1 þ u2 he he for the nodal axial displacements du :¼ {u1, u2}. The standard axial strain operator    1 1 Bea ¼   ð3:13Þ he he so eha ¼ B ea d u , is recovered. We consider the C1 interpolation of the transversal deflection

Fig. 1. Sketch of the ‘‘hinge mode’’ for a two node element. The finite element has to be able to represent an opening of the hinge without a strain in the rest of the element, that is, the regular part of the strain must vanish for these nodal values.

F. Armero, D. Ehrlich / Computers and Structures 84 (2006) 641–656

^ 1 N w1 ðgÞ þ w ^ 01 N w2 ðgÞ þ w ^ 2 N w3 ðgÞ þ w ^ 02 N w4 ðgÞ wh ðgÞ ¼ w

ð3:14Þ

based on the cubic Hermitian polynomials N w1 ðgÞ ¼ 2g3  3g2 þ 1;

N w2 ðgÞ ¼ he ðg3  2g2 þ gÞ;

N w3 ðgÞ ¼ 2g3 þ 3g2 ;

N w4 ðgÞ ¼ he ðg3  g2 Þ ð3:15Þ

in terms of the natural coordinate g = s/he and the nodal degrees of freedom d w :¼ fw1 ; w01 ; w2 ; w02 g corresponding to the nodal deflections and rotations. The associated curvature strain operator is then given by "

       # 6 2s  2 3s  6 2s  2 3s   Bj ¼  2 1  2 1 1 he  he he  h2e he  he he he

ð3:16Þ

so jh = Bjd w. The kinematics of the hinge mode for the considered two-node finite element is depicted in Fig. 1. We observe the decoupling of the longitudinal and transversal components in the assumed infinitesimal range, leading to the relations " # " # " # ^u1 ^ u1 0 ¼ þ nu ^ ^ u2 u1 1 |fflffl{zfflffl} |fflffl{zfflffl} |ffl{zffl} u d^ hinge

u d^ rigid

Duhinge

and 3 2 3 2 3 2 ^1 ^1 w w 0 7 6 7 6 07 6 ^ 7 6 ^ 01 w 7 6 0 7 6w 7þ6 7 nh 6 17¼6 7 6^ þw 7 6w 6 ^ 01 he 7 5 4 he  s d 5 4 ^2 5 4 w 1 ^ 02 w |fflfflffl{zfflffl ffl}

^ 01 w |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl}

w d^ hinge

w d^ rigid

ð3:17Þ

Dw hinge

ð3:18Þ

recovering the constant enhanced strain operator (3.4), as it can be expected for the underlying constant axial strain element. On the other hand, for the curvature component relation (3.11) combined with (3.16) and (3.17)2 leads to the enhanced strain operator G ðexactÞ c;j

¼

Bj Dwhinge

¼

gð1Þ e ðsÞ he

where we have introduced the Heaviside function H sd ðsÞ ¼ hs  sd i= j s  sd j for the Macaulay brackets h Æ i (i.e. hgi = g for g P 0, hgi = 0 for g 6 0). However, the functions (3.20) and (3.21) are not required in the final numerical integration, which involves only the enhanced strain operators (3.6), (3.18) and (3.19). 3.4. Additional remarks on the numerical implementation (i) The discrete linearized equations. The non-linear set of algebraic equations (3.3) is solved through a Newton– Raphson iterative scheme, resulting in the linear algebraic system   " ðiÞ # " ðiÞ # nelem K e;dd K e;dn Dd nþ1 Rnþ1 A ¼ ðiÞ ðiÞ ð3:22Þ e¼1 K e;nd K e;nn rnþ1 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Dnnþ1 ðiÞ

for the axial and transversal displacements, respectively. Using relation (3.11) with the axial strain operator (3.13) and the hinge matrix (3.17)1, we obtain the enhanced strain operator for the axial strain component ¼ Bea Duhinge G ðexactÞ c;ea

Remark 3.1. The interpolation functions Wu(s) and Ww(s) in (2.2) can be obtained by the integration of the relation (3.11) given (2.3) and (2.4). For the enhanced two-node Hermite finite element considered in this section, with the enhanced strain operators (3.18) and (3.19), this calculation leads to the local interpolation functions s Whu ðsÞ ¼  H sd ðsÞ; ð3:20Þ he       s2 2sd 2s 1þ3 1 Whw ðsÞ ¼ hs  sd i  1 ð3:21Þ he he he

K

1 |fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl}

1 gð0Þ ¼ ¼ e he he

647

    1 2sd 2s ¼ 1þ3 1 1 he he he ð3:19Þ

recovering, remarkably, the equilibrium enhanced operator ¼ G ðexactÞ . This situation has the important result G ðexactÞ e c that the final formulation leads to symmetric linearized systems of equations as described in the next section.

ðiÞ

on the incremental values Dd nþ1 and Dnnþ1 for iteration (i) ðiÞ ðiÞ of increment tn+1, with the residuals Rnþ1 and rnþ1 given by (3.3). The stiffness matrices in (3.22) are obtained by the consistent linearization of the residuals (3.3) leading to the elemental expressions ) Rh Rh K e;dd ¼ 0 e BT CB ds; K e;dn ¼ 0 e BT CG c ds; Rh Rh ed K e;nd ¼ e G T CB ds; K e;nn ¼ e G T CG c ds þ C 0

e

0

e

ð3:23Þ for the material tangents e d n_ e r_ ¼ C
e_ h and r_ d ¼ C e

ð3:24Þ

of the bulk and localized models, respectively. The increment and iteration indices have been omitted when referring to the different stiffness matrices to simplify the notation. We note that, given the limit values Ge, Gc  1/he, the matrix Ke,nn is invertible for small enough he despite the negative definiteness of the localized softening e d in (3.23); see [2] for a complete discussion of tangent C this issue in the context of continuum problems. The element character of the enhanced equation (3.3)2 allows the elimination the enhanced parameters at the element level as h i ðiÞ ðiÞ ðiþ1Þ Dne;nþ1 ¼ K 1 ð3:25Þ in Be;loc e;nn re;nþ1 þ K e;nd Dd e;nþ1 which combined with (3.22) results in the statically condensed system

648

K



F. Armero, D. Ehrlich / Computers and Structures 84 (2006) 641–656

ðiÞ Dd nþ1

¼

ðiÞ Rnþ1

8 i nelem h > K < K  ¼ A K e;dd  K e;dn K 1 e;nd e;nn e¼1 for h i n elem ðiÞ > : R ¼ A RðiÞ  K 1 rðiÞ e;nþ1 nþ1 e;nn e;nþ1 e¼1

ð3:26Þ ðiÞ Dd nþ1

involving the nodal degrees of freedom only. That is, the implementation reduces to an enhanced large-scale problem capturing efficiently and objectively the localized effects of the softening hinge. (ii) Eigenvalue analysis. We gain an additional insight on the improvement obtained with the use of the enhanced strain operator G ðexactÞ developed in Section 3.3 by performc ing an eigenvalue analysis of the uncondensed stiffness matrix K in (3.22). We consider only the bending contributions as the axial part coincides with the original constant operator (3.18). In this case, and considering a fully softened hinge (i.e. the hinge mode characterized by Md = 0 _ d ¼ 0), we find that the enhanced element based on and M the non-constant enhanced strain operator (3.19) has three zero-energy modes spanning its kernel 3 2 ð Dw1 Þ 1 0 0 7 6 Dw01 Þ 60 1 0 7ð 7 6 7 6 Dw2 Þ kern½K G ðexactÞ
¼ span6 1 he he  sd 7 ð ð3:27Þ c 7 6 7 60 1 0 Dw2 Þ 1 5ð 4 0

0

ð

1

Dnhe Þ

corresponding to the original rigid body modes of the base finite element (a transversal rigid translation and infinitesimal rotation), and an additional zero-energy mode corresponding to a relative rotation about the hinge at s = sd (the hinge mode). On the other hand, an analysis of the stiffness matrix of the finite element based on the constant enhanced operator G ðconstÞ in (3.4) shows a different kernel, spanned by the c eigenvectors 3 2 ð Dw1 Þ 1 0 0 7 6 Dw01 Þ 60 1 0 7ð 7 6 7 6 Dw2 Þ kern½K G ðconstÞ
¼ span6 1 he he =2 7 ð ð3:28Þ c 7 6 7 60 1 Dw02 Þ 1 5ð 4 0

0

1

ð

Dnhe Þ

The original rigid body modes are also present, and still an extra zero-energy mode so the element should not exhibit stress-locking. However, this additional mode corresponds to a relative rotation about s = he/2, instead of the true physical location of the hinge at s = sd. Note that the stiffness matrices of both elements coincide in the limit he ! 0, showing the consistency of both approaches. However, the element based on the enhanced operator G ðexactÞ is able to c reproduce exactly the hinge mode for finite size meshes, improving considerably on the accuracy of the solutions in these meshes. These observations are confirmed by the numerical examples presented in Section 4.

(iii) Detection of the hinge. Following the discussion in Section 2.3, the plastic hinge is activated when the bulk model rð
eÞ detects the start of perfect plasticity/damage. The check for the hinge activation is done at the quadrature points employed in the evaluation of the different integrals above. In particular, for the two-node Hermite finite element described in Section 3.3.1, we use a three-point Gauss–Lobatto quadrature with two quadrature points located at the end nodes and a third point at the elements center (see e.g. [17, p. 440]). The hinge is activated in the quadrature point where it is detected. In this way, and focusing on the bending response, the hinge is first activated at either end of the beam element due to the linear distribution of the bending moment or, if the moment distribution is constant, so the yielding is observed at all quadrature points simultaneously, the hinge is activated at the center of the element only. (iv) Numerical integration of the hinge model. We refer to [1] and references therein for details on the numerical integration of localized models in the context of continuum problems. This corresponds to the evolution equations (2.10) and (2.11) for the softening hinges of interest here. A particularly interesting approach is obtained by considering an elasto-plastic regularization of these rigid-plastic e e , for a penalty equations (i.e. n = ne + np with rd ¼ Cn e and with n_ p given by (2.10)1). In this form, tangent C, the integration scheme follows the standard closest-pointprojection return mapping algorithms as applied to the elastoplastic models like (2.9) in the continuum; see e.g. [33]. We also note that if Eq. (2.10)1 defines a constant flow ~ for the particular problem/model at vector n/ :¼ ord / hand, only one component of the enhanced parameters n needs to be activated, that is, nn/ ¼ n  n/ . The governing equation (3.3)2 is then also projected in this direction defining rdn/ ¼ n/  rd . This is, for example, the case for a beam developing softening hinges in bending only, like in some of the numerical simulations presented in Section 4. The finite elements involve only an enhanced curvature in bending through the rotation jump nh in this particular case. 4. Representative numerical simulations We present in this section several numerical simulations that illustrate the performance of the new enhanced finite elements developed in this paper. We consider the twonode Hermite element of Section 3.3.1 enhanced with the original constant enhancement (3.4), denoting the final element by EB-(constant), and with the enhancement defined by (3.6), (3.18) and (3.19), referred to as the EB-(exact) element. Section 4.1 presents a numerical evaluation of the two finite elements for a cantilever beam, with the different resolution of the kinematics of the plastic hinge becoming apparent in this simple setting. Section 4.2 presents a comparison of the results obtained with the proposed finite elements with some existing experimental results. Finally, we

F. Armero, D. Ehrlich / Computers and Structures 84 (2006) 641–656

consider in Section 4.3 the collapse of a frame first considered in [11]. 4.1. Numerical evaluation We evaluate in this section how the two enhanced finite elements developed in the previous section resolve the kinematics of a softening hinge. To this purpose we consider a cantilever beam of length L = 1, with a fixed end on the left
and with imposed axial
u and transversal displacements w at the right end. Given the linear distribution of the resulting bending moment M, with a maximum at the fixed end, and constant axial force N along the beam, we expect the formation of a softening hinge at the fixed end. We consider a linear elastic response in the bulk, characterized by a axial stiffness EA = 350 and a bending stiffness EI = 1000, with yielding detected according to the interaction diagram  2   M  N /ðN ; MÞ ¼ ð4:1Þ þ    ð1  qÞ 6 0 Ny My for the limit moment My > 0 and axial force Ny > 0. Perfect plasticity is assumed (that is, q  0 in (4.1)) or, in other words, that a hinge is activated when / = 0 according to the discussion in Section 2.3. The response of the hinge follows the localized plastic model (2.10) in terms of the localized interaction diagram  2   M d  Nd ~ þ    ð1  ~ qÞ 6 0 ð4:2Þ /ðN d ; M d Þ ¼ Nu Mu in terms of the driving moment Md and driving axial force Nd at the hinge, for Nu = Ny and Mu = My in this particular case (that is, localized softening is activated upon yielding). The simulations presented here were obtained with the normalized values Nu = 1 and Mu = 1. We consider the associated evolution equations (2.10) and (2.11), with the linear localized softening law f ag ~ qð~ aÞ ¼ minf1;  H~

ð4:3Þ

f The value H f¼ for the localized softening modulus H. 100 is considered in the simulations. The interaction diagram (4.1) corresponds to the ultimate yield locus for an homogeneous rectangular cross-section, with a symmetric response of the material in tension and compression; see [20, p. 338] for its derivation and calibration. See also formulas (4.6) below. In this example we have neglected the hardening phase from an initial yield surface to the ultimate diagram (4.1); we refer to the example presented in the next section for an experimental evaluation accounting for this effect. We also refer to the discussion in the beginning of Section 2.3, which includes several references focused entirely in the derivation and validation of yield surfaces for specific structural members of steel, reinforced concrete and composite cross-sections, and their evolution in stress space (i.e. hardening) to model the gradual plastification over the cross-section. We note

649

again that the interest here is the purely numerical evaluation of how the finite elements resolve the kinematics of the hinges once they form. The different conclusions obtained here regarding this numerical resolution apply entirely to other situations, with alternative evolving surfaces through general hardening laws. In particular, the interaction diagram (4.2) results in the coupling of the axial and bending response in the hinges, allowing the evaluation of the numerical treatment of the coupled kinematics by the finite elements developed in this work. We also observe the non-smooth character of the interaction diagrams (4.1) and (4.2) at M = 0. The associated plastic evolution equations (2.9) and (2.10), and their numerical integration, are generalized to this case by following the standard arguments behind the framework of multi-surface plasticity presented in [33] for continuum plasticity. Details are omitted here. We first consider the case of no imposed axial displacement
u ¼ 0 with a monotonically increasing imposed
. No axial force develops in this transversal displacement w case. Following the comments in Section 3.4, the finite element capturing the hinge activates only one enhanced parameter corresponding to this rotation jump. Fig. 2 depicts the solution obtained in this case for both the EB-(constant) and EB-(exact) finite elements and different spatial discretizations. We consider meshes with 1, 3, 9 and 27 equally spaced elements. The response is initially linear elastic, where both elements coincide (no enhancements active yet) and with all the meshes giving the exact solution in this problem involving a linear bending moment for the considered Hermite element. This can be observed in the included plots depicting the evolution of the computed transversal reaction Ry for the imposed
at the beams right end. This transversal displacement w situation changes, however, when the softening hinge is formed at the fixed end after the detection of yielding at this point according to Eq. (4.1). The first observation is the different solutions obtained by the EB-(constant) element for different level of discretization in the softening hinge range. We have also included in Fig. 2 the spatial distributions of the bending moment at
¼ 6  103 , before the plastic the imposed displacement w hinge is fully softened for both elements, showing again a completely different solution for different meshes. In contrast, the EB-(exact) element continues providing the exact solution in this range with any finite element mesh. Clearly no mesh dependence occurs for this element. We note that the EB-(constant) does not suffer of the pathological mesh-size dependence characteristic of a standard Galerkin finite element in combination of a bulk plastic model like (2.9) with strain softening. Both enhanced elements capture objectively the localized dissipation associated with the softening hinge. The EB-(constant) element is not able to capture the exact solution for coarse meshes, but the numerical solutions converge eventually to the exact solution as the mesh is refined showing the consistency of the choice (3.4) and its objectivity.

650

F. Armero, D. Ehrlich / Computers and Structures 84 (2006) 641–656


only. Fig. 2. Cantilever beam: solution obtained with the EB-(constant) and EB-(exact) finite elements for an imposed transversal displacement w

No stress-locking is observed for either element, since the fully softened hinge (i.e. vanishing reaction force for increasing deformation) is captured in both cases. These results are in perfect agreement with the observations made from the eigenvalue analysis presented in Section 3.4. The EB-(exact) element reproduces exactly the kinematics of the softening hinge at its position, in contrast with the EB-(constant) element which cannot reproduce this mode exactly as observed when comparing Eqs. (3.27) and (3.28). The superiority of the EB-(exact) element is apparent. The same conclusions hold when we consider the element performance in a problem involving not only bending moment but also an axial force and their couplings as governed by the interaction diagram (4.2). To illustrate this point we consider the case where an axial displacement
u ¼ 2:5  103 is applied first at the right end and kept con-


is increased stant while the transversal displacement w monotonically as before. The response during the elastic range coincides with the one obtained in the previous case, involving now a constant value of the axial force as the transversal displacement is applied. This situation changes completely when yielding is detected, again at the fixed end, and a softening hinge forms following the coupled interaction diagram (4.2). Obviously, by equilibrium, the spatial distributions of the bending moment and axial force remain linear and constant, respectively. However, their values evolve for the increasing transversal displacement as the plastic hinge at the fixed end softens and both a displacement and rotation jumps, nu and nh respectively, develop. Figs. 3 and 4 depict the solutions obtained with the EB(constant) and the EB-(exact) elements respectively. We have included the evolution of the transversal reaction Ry, the axial reaction Rx, the axial displacement jump nu

F. Armero, D. Ehrlich / Computers and Structures 84 (2006) 641–656

651


at the free end while the axial Fig. 3. Cantilever beam: solution obtained with the EB-(constant) finite element for imposed transversal displacement w displacement
u is kept fixed.

and the rotation jump nh at the hinge, all versus the
. The reader can also imposed transversal displacement w

find a plot showing the combined evolution of the reacting moment and axial force at the fixed end (that is, at the

652

F. Armero, D. Ehrlich / Computers and Structures 84 (2006) 641–656


at the free end while the axial Fig. 4. Cantilever beam: solution obtained with the EB-(exact) finite element for imposed transversal displacement w displacement
u is kept fixed.

F. Armero, D. Ehrlich / Computers and Structures 84 (2006) 641–656

hinge). We have also depicted in that plot the original elastic domain given by (4.2). We observe that upon reaching this original surface, after the moment has increased for a constant axial, the coupling of the two effects becomes clear with the elastic domain softening according to the linear softening law (4.3). The stress resultant pair (Nd, Md) follows the normality rule (2.10) which combined with the localized softening leads to a reduction of the axial force with an eventual reduction of the driving moment after some initial stages where it increases. The hinge is finally fully softened at Nd = 0 and Md = 0 with the rest of the beam, fully unloaded now, rotating as an (infinitesimal) rigid body about the opened hinge. This is observed clearly in the included plot showing the evolution of the displacement jumps at the hinge versus the imposed transversal displacement. The increasing transversal displacement is taken completely in this final stages by the hinge through its rotation, while the hinge takes all the imposed axial dis
¼ 1:5  103 in these final stages after placement nu ¼ w the hinge has softened completely its load capacity. We observe that the EB-(exact) element is able to reproduce exactly these different aspects of the solution for all the different meshes. The EB-(constant) element in Fig. 3 only captures this solution in the limit case of fine meshes. It provides again a consistent and objective solution for coarse meshes, but with a lesser performance when compared with the EB-(exact) element due to the crude approximation that the constant enhanced strain fields (3.4) provide. This situation is to be contrasted with similar elements that have been proposed in the context of continuum problems involving the very same type of constant enhancement; see e.g. [1,2] and references therein. This situation is to be traced to the more complex kinematics of the beam when compared to the continuum, involving in the thin limit under consideration here the second derivative of the displacements (the curvature) and the need for an enhancement of the first derivative (the rotation). 4.2. Experimental evaluation We assess in this section the proposed formulation for modeling localized failures of beams through the consideration of experimental results reporting the post-limit softening response of the beam. The reinforced concrete experiments performed by Lane [19], as reported by Schreyer and Chen [32], are considered to this purpose. Several simply supported beams of span 3.82 m were loaded to failure through the application of imposed displacements at distances of La to the support, for La equal to 0.96 m, 1.30 m and 1.60 m. The inset in Fig. 5 illustrates the geometric configuration of the specimens. The reported experimental results are depicted in Fig. 5 through the plot of the vertical reaction force versus the displacement at the center of the beam span. Results were obtained for different loading geometries characterized by different distances from the supports La of the points of application of the vertical dis
. placement w

653

Fig. 5. Experimental evaluation. Vertical reaction force versus displacement at the center of the beam for different displacement application points. Solid lines represent experimental results.

We consider the following model in the numerical simulations. The bulk response follows the elastoplastic model (2.9) based on the yield surface (4.1) which in the considered case of vanishing axial reduces to   M  ð4:4Þ /ðM; qÞ ¼    ð1  qÞ with qðaÞ ¼ Ha M y

for the maximum elastic moment My > 0 and the linear hardening modulus H > 0. This is followed by a localized softening response once the bending moment reaches the ultimate capacity Mu > My, characterized through the ultimate surface (4.2) for the driving bending moment Md at the plastic hinge. The localized softening response of the hinge follows the parabolic softening law ~a2 g ~qð~aÞ ¼ minf1; @~

ð4:5Þ

for the localized softening parabolic coefficient @~ < 0. The numerical simulations that follow assumed the following material parameters: linear elastic stiffness EI = 2.04 · 107 N m2, elastic limit My = 2.1 N m, ultimate bending moment capacity Mu = 2.7 N m, strain softening hardening parameter H ¼ 3:5 N1 and localized softening parabolic coefficient @~ ¼ 60 N2 m2 . The focus is again on the bending response of the beam. The particular choice of the model and material parameters follow the lines discussed in [32]. Due to the symmetry of the problem, we consider simulations for only half of the beam, with symmetry boundary conditions applied at the right-hand end. The half-beam was modeled with 6 and 12 element meshes, discretized as follows. For the La = 0.96 m, Lb = 0.95 m case, we discretized both La and Lb into 3 elements in the 6-element mesh and into 6 elements in the 12-element mesh, due to their similar span. In the La = 1.30 m and Lb = 0.61 m case we discretize La into 4 elements and Lb into 2 elements in the 6-element case, while discretizing La into 8 elements and Lb into 4 elements in the 12-element mesh case. For

654

F. Armero, D. Ehrlich / Computers and Structures 84 (2006) 641–656

the La = 1.60 m and Lb = 0.31 m loading geometry we discretize La into 5 elements and Lb into a single element in the 6-element mesh and in the 12-element mesh we discretize La and Lb into 10 and 2 elements, respectively. The above discretizations result in meshes comprised of elements of similar sizes, regardless of the loading geometry. Due to the constant moment distribution along the Lb span, an imperfection is introduced at the extreme righthand element, with a slightly lower ultimate moment capacity Mu. Fig. 5 depicts the numerical solutions through the plot of the vertical reaction force versus the displacement at the center of the beam span. We observe the ability of the proposed formulation to adequately model the experimental data. After the initial elastic range, the central Lb region enters the hardening regime upon reaching the limit moment My. A plastic hinge is eventually triggered by the imperfection at the extreme right-hand element when reaching Mu. All quadrature points reach the ultimate capacity simultaneously, resulting in a plastic hinge located at the middle of the element, according to the procedure we outlined in Section 3.4. Since the G ðexactÞ enhanced operator c collapses to the G ðconstÞ operator for a hinge located at the c middle of the element, the same results are obtained with the EB-(exact) and EB-(constant) elements in this case. We emphasize again that the main goal of this numerical example is the evaluation of the proposed formulation to model softening hinges in beams. The simulation continues with the moment capacity at the plastic hinge decreasing with added imposed displacement, while the rest of the beam unloads elastically. We note that similar results are obtained with 6- and 12-element meshes, with the slope of the load–deflection response independent of the mesh discretization, confirming again the objectivity of the formulation. Note that, in contrast with the simpler example considered in the previous section, the different meshes in the current example need to resolve the propagation of the hardening plastic response in the bulk. Once again we note the good agreement of the solutions obtained with different mesh discretizations.

Fig. 6. Darvall–Mendis portal frame: geometry and loading.

numerical simulations reported here consider the load applied at fixed increments of DP = 22.2 N. We discretize the frame with a total eight finite elements, with two elements of equal size for each column and four elements for the beam (two elements of 0.275L and two elements of 0.225L). The softening response follows the initial elastic range, that is, the hinges are activated when the limit moment is reached. The analysis is performed for a number of differf in order to compare the ent localized softening moduli H influence of the softening hinge on the predicted ultimate load of the frame. As discussed in Section 1, the appearance of softening hinges leads to significant size effects in the response of the structures at failure. For example, for a rectangular cross-section of depth D and width b, the maximum elastic moment and softening modulus are given, in non-dimensional form, 450 Perf. plast. H = –1.089 H = –1.681 H = –2.017

440 430

4.3. Darvall–Mendis frame Vertical load, kN

We consider next the clamped portal frame under vertical loading first studied in [11]. The geometry of the frame is shown in Fig. 6 with the numerical data as follows: the length of both columns and the beam is L = 3 m, their cross-section area A is 0.10 m2 and the moment of inertia I is 1 · 103 m4. The elastic material response is modeled through the Young’s modulus E = 20.7 GPa and moment limits M y col ¼ 158 kN m for the columns and M y col ¼ 169 kN m for the beam. Elastic response is assumed in axial (i.e. Ny, Nu ! 1), focusing again the attention on the bending problem with a localized model involving linear softening in the driving moment. The frame is loaded through a vertical load applied at a distance of 0.55L to the top left corner; see Fig. 6. The

420 410 400 390 380 370 360 350 340 4

5

6

7

8

9

10

11

12

13

14

Vertical deflection at load point, mm

Fig. 7. Darvall–Mendis portal frame: vertical load versus vertical deflecf The markers tion at the load point for different softening moduli H. indicate structural collapse under load increase.

F. Armero, D. Ehrlich / Computers and Structures 84 (2006) 641–656

450

450 8 elements 16 elements 32 elements 64 elements

440 430

8 elements 16 elements 32 elements 64 elements

440 430

420

420

Vertical load, kN

Vertical load, kN

655

410 400 390 380

410 400 390 380

370

370

360

360

350

350

340

340 4

5

6

7

8

9

10

11

12

13

14

Vertical deflection at load point, mm

4

5

6

7

8

9

10

11

12

13

14

Vertical deflection at load point, mm

Fig. 8. Darvall–Mendis portal frame: comparison of EB-(constant) and EB-(exact) vertical load versus vertical deflection at load point for different discretizations. Red markers indicate structural collapse under load increase.

My fr ¼ 2 6 EbD

and

f fr2 H ¼  D 72EGf EbD2

ð4:6Þ

for the strength fr and the fracture energy Gf of the material; see [5,6]. Larger structures (larger D) lead then to a f ¼ 0; 1:089; 1:681 more brittle response. The values H and 2.017 kN m/rad are considered. This example is particularly relevant since it shows how the ultimate load of a structure varies with the localized softening in the hinges. We compute the evolution of the vertical displacement under the load up to the maximum value attained by the load, which is the only goal in this example. The post-peak response can be traced, as in the previous examples, by imposing the displacement instead or, more generally, by standard arc-length continuation methods to account for possible snap-back responses. See [5,6] for a study on the conditions for these responses to occur. In all cases, hinges are formed first on both elements sharing node 5, where the point-load is applied, when its magnitude reaches 337.10 kN. The hinges are formed next to this f ¼ 0Þ, subsecommon node. In the perfect plastic hinge ð H quent loading causes redistribution of the moments on the frame, causing a new hinge to be formed in element 7, at node 7, when the load reaches 428.23 kN. As the point load reaches 434.32 kN, a hinge is formed at node 3 of element 2, and the structure collapses under load control (that is, the load would eventually be reduced for an increasing deformation). Fig. 7 includes the load–displacement laws obtained in the simulation, depicting clearly these transition points. For the case of softening hinges, the structure collapses as soon as the second hinge, at node 7 on element 7, is formed. The final collapse load reduces drastically as the softening modulus is reduced (increased in absolute value). f ¼ 1:089 kN m/rad the structure is not able to For H support an additional load at the 388.02 kN level. For f ¼ 1:681 kN m/rad, the collapse happens at 358.04 H f ¼ 2:017 kN m/rad, it happens almost kN, while for H

as soon as the first hinge is formed. The reduction by 22.4% of the collapse load shows how important is the correct modeling of the softening hinges for practical purposes. Fig. 8 presents a comparison of the performance of the EB-(exact) and EB-(constant) elements when estimating the collapse load of the structure for a localized softening f ¼ 1:089 kN m/rad. Clearly, the EB-(conmodulus of H stant) element overestimates the ultimate load for coarse meshes. Meaningful results are attained, however, after the discretization is refined, as discussed in Section 4.1. On the other hand, the EB-(exact) element shows again to be more suitable for the analysis of the failure of the structural system. This element is able to obtain the correct response even for small number of elements per member for this problem involving linear bending and constant axial force distributions, and the assumed elastic/localized softening model. 5. Concluding remarks We have presented in this paper a general framework for the incorporation of softening hinges in thin beam finite elements. A careful examination of the local kinematics associated with the hinges, understood here as discontinuities of the generalized displacements characterizing the deformation of the beams, allowed the design of new enhanced strain finite elements that capture exactly the kinematics of these discontinuous solutions. The focus in this paper has been in the finite element development, but the general framework proposed here can easily incorporate additional different effects in the response of the localized model of the hinge. We have presented several numerical simulations that show the added numerical accuracy gained by the proposed approach. In particular, it leads to a very efficient method to capture objectively the

656

F. Armero, D. Ehrlich / Computers and Structures 84 (2006) 641–656

localized dissipative mechanisms characteristic of the localized failures of beams and general frames. We are currently extending these results to additional cases of interest. The incorporation of finite deformation in the beam response is one particular extension that can be easily envisioned. Especially interesting is the use of the ideas and techniques described in this paper for the modeling of hinge lines in plates and shells. We intend to present these results in forthcoming publications. Acknowledgements Financial support for this research has been provided by the ONR under grant no. N00014-00-1-0306 with UC Berkeley. This support is gratefully acknowledged. References [1] Armero F. Large-scale modeling of localized dissipative mechanisms in a local continuum: applications to the numerical simulation of strain localization in rate-dependent inelastic solids. Mech CohesiveFrictional Mat 1999;4:101–31. [2] Armero F. On the characterization of localized solutions in inelastic solids: an analysis of wave propagation in a softening bar. Comp Meth Appl Mech Eng 2001;191:181–213. [3] Armero F, Ehrlich D. An analysis of strain localization and wave propagation in plastic models of beams at failure. Comp Meth Appl Mech Eng 2003;193:3129–71. [4] Bazant ZP. Instability, ductility, and size effect in strain-softening concrete. ASCE J Eng Mech 1976;102:331–44. [5] Bazant ZP. Scaling of structural strength. New York: Hermes Penton Ltd; 2002. [6] Bazant ZP. Asymptotic matching analysis of scaling structural failure due to softening hinges. I: Theory. II: Implications. ASCE J Eng Mech 2003;129:641–50. 651–4. [7] Bazant ZP, Pan J, Pijaudier-Cabot G. Softening in reinforced concrete beams and frames. ASCE J Struct Eng 1987;113:2333–47. [8] Cocchetti G, Maier G. Elastic–plastic and limit-state analyses of frames with softening plasti-hinge models by mathematical programming. Int J Solids Struct 2003;40:7219–44. [9] Coleman J, Spacone E. Localization issues in force-based frame elements. ASCE J Struct Eng 2000;127:1257–65. [10] Darvall PL. Critical softening of hinges in portal frames. ASCE J Struct Eng 1984;110:157–62. [11] Darvall PL, Mendis PA. Elastic–plastic-softening analysis of plane fames. ASCE J Struct Eng 1984;111:871–88. [12] Ehrlich D, Armero F. Finite element methods for the analysis of plastic hinges in beams and frames. Comput Mech 2005;35:237–64. [13] El-Tawil S, Deierlein G. Stress-resultant plasticity for frame structures. ASCE J Eng Mech 1998;124:1360–70. [14] El-Tawil S, Deierlein G. Nonlinear analysis of mixed steel–concrete frames. I: Element formulation. ASCE J Struct Eng 2001;127:647–55. [15] El-Tawil S, Deierlein G. Nonlinear analysis of mixed steel–concrete frames. II: Implementation and verification. ASCE J Struct Eng 2001;127:656–65.

[16] Hillerborg A. Fracture mechanics concepts applied to moment capacity and rotational capacity of reinforced concrete beams. Eng Fract Mech 1990;35:233–40. [17] Hughes TJR. The finite element method. Englewood Cliffs, NJ: Prentice-Hall; 1987. [18] Jirasek M. Analytical and numerical solutions for frames with softening hinges. ASCE J Eng Mech 1997;123:8–14. [19] Lane GE. Behavior of reinforced concrete beams under combined axial and lateral loading. Air Force Weapons Laboratory. Kirtland Air Force Base. New Mexico, AFWL-TR-76-130, May 1972. [20] Lubliner J. Plasticity theory. New York: MacMilan Publisher; 1990. [21] Maier G. Instability due to strain softening. Stability of continuous systems, Proc. IUTAM symp. Berlin: Springer-Verlag; 1971. p. 411–7. [22] Maier G, Zavelani A, Dotreppe JC. Equilibrium branching due to flexural softening. ASCE J Eng Mech 1973;99:897–901. [23] Marante ME, Picon R, Florez-Lopez J. Analysis of localization in frame members with plastic hinges. Int J Solids Struct 2004;41: 3961–75. [24] Moe¨s N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. Int J Numer Meth Eng 1999;46: 131–50. [25] Nigam NC. Yielding in framed structures under dynamic loads. ASCE J Eng Mech 1970;96:687–709. [26] Oliver J. Modelling strong discontinuities in solid mechanics via strain softening constitutive equations. Part 1: Fundamentals. Int J Numer Meth Eng 1996;39:3575–600. [27] Oliver J. Modelling strong discontinuities in solid mechanics via strain softening constitutive equations. Part 2: Numerical simulation. Int J Numer Meth Eng 1996;39:3600–23. [28] Orbison JG, McGuire W, Abel JF. Yield surface applications in nonlinear steel frame analysis. Comp Meth Appl Mech Eng 1982;33:557–73. [29] Ortiz M, Pandolfi A. Finite deformation irreversible cohesive elements for the three-dimensional crack-propagation analysis. Int J Numer Meth Eng 1999;44:1267–82. [30] Powell GH, Chen PF. 3D beam-column element with generalized plastic hinges. ASCE J Eng Mech 1986;112:627–41. [31] Sanjayan JG, Darvall PL. Singularities in RC beam elements with finite-length hinges. ASCE J Struct Eng 1995;121:39–47. [32] Schreyer HL, Chen Z. The effect of localization on the softening behavior of structural members. In: Willam KJ, editor. Constitutive equations: macro and computational aspects, ASME winter annual meeting, New Orleans, LA. New York: ASME; 1984. p. 193–203. [33] Simo JC, Hughes TJR. Computational inelasticity. Springer-Verlag; 1998. [34] Simo JC, Oliver J, Armero F. An analysis of strong discontinuities induced by softening solutions in rate independent solids. J Comput Mech 1993;12:277–96. [35] Spacone E, El-Tawil S. Nonlinear analysis of steel–concrete composite structures: state of the art. ASCE J Struct Eng 2004;130:159–68. [36] Stakgold I. Green’s functions and boundary value problems. 2nd ed. New York: John Wiley & Sons; 1998. [37] Takizawa H, Aoyama H. Biaxial effects in modeling earthquake response of RC structures. Earthquake Eng Struct Dyn 1976;4: 523–52. [38] Wells GN, Sluys LJ. A new method for modeling cohesive cracks using finite elements. Int J Numer Meth Eng 2001;50:2667–82.