Efficient Beam-Column Element

with Variable Inelastic End Zones

Chin-Long Lee1 and Filip C. Filippou, M.ASCE2
Downloaded from ascelibrary.org by INDIAN INSTITUTE OF TECHNOLOGY on 06/05/14. Copyright ASCE. For personal use only; all rights reserved.

Abstract: The paper presents a new beam-column element that combines computational efficiency with accuracy. The element uses only
one monitoring section in each end inelastic zone of the structural member, but accounts for the spread of inelastic deformations under
strain hardening response. It is, therefore, a variable inelastic zone model that combines the benefit of integrating the section moment-
curvature response with the computational efficiency of concentrated hinge models for beams and columns. The element is more accurate
than current distributed inelasticity models in simulating the monotonic and cyclic inelastic response of beams and columns under typical
curvature distributions, while being much more cost effective. The element is also suitable for softening response, as long as a relation is
available between the length of the softening zone and element response parameters.
DOI: 10.1061/共ASCE兲ST.1943-541X.0000064
CE Database subject headings: Beam columns; Plasticity; Inelasticity; Plastic hinges; Parameters.

Introduction moment by means of resultant plasticity theory 共Hilmy and Abel

1985; Powell and Chen 1986兲, the inclusion of the effect of the
Beam-column elements are in common use in structural engineer- pullout of reinforcing steel 共Filippou and Issa 1988兲, the effect of
ing practice, because of their computational efficiency and ease of composite action in concrete-filled tubes 共Hajjar and Gourley
response interpretation. In performance-based earthquake engi- 1997兲 and other refinements.
neering such elements are often used in nonlinear push-over In spite of these efforts a serious limitation of the concentrated
analyses and in dynamic response studies. Under the action of hinge model remains: the inelastic flexural deformations are con-
high lateral forces or intense base excitations large inelastic de- centrated in a plastic hinge whose moment-rotation characteristics
formations take place at the end of the structural members. The need to include not only material parameters, but also the effect of
zone of inelastic response grows with the response history and the inelastic zone length, which depends on boundary conditions
affects the determination of local response quantities such as plas- and element response parameters like the moment distribution.
tic strains. Thus, the calibration of the plastic hinge properties becomes im-
Because inelastic deformations concentrate at the end of the portant and is fraught with misinterpretation 共Mahin and Bertero
structural members, the most economical approach for modeling 1976; Anagnostopoulos 1981兲.
the cyclic inelastic behavior of beams and columns is an element To overcome the limitation of the concentrated hinge model to
with concentrated inelasticity in the form of a hinge to represent explicitly account for the spread of inelastic deformations into the
the inelastic zone. The earliest such model is the two-component element a couple of studies proposed a spreading inelastic zone
parallel model 共Clough et al. 1965兲. To overcome some of its beam element 共Soleimani et al. 1979; Meyer et al. 1983兲. This
serious limitations the one-component model was proposed a few element retains the computational efficiency of the concentrated
years later, which consists of a linear elastic frame element that is hinge model by monitoring the inelastic response at a single sec-
connected in series with two rigid-hardening springs, one spring tion. For determining the element deformations the element
for each inelastic zone 共Giberson 1967兲. The versatility of this makes assumptions about the curvature distribution in the inelas-
approach lies in the ability to insert additional inelastic springs in tic zone and is governed by several ad hoc rules during loading
series to account for different deformation mechanisms, such as and reloading of the element. Another proposal in the same vein
semirigid connectors, pullout of reinforcing steel from the joint, is the quasi-plastic-hinge model by Attalla et al. 共1994兲. In this
etc. In the years following its introduction many extensions of the model a closed-form expression of the element force-deformation
original model have been proposed. Developments include the relationship is derived by exact integration of the section defor-
description of the interaction between axial force and bending mations. A rather complex expression results under monotonic
loading whose extension to cyclic load conditions is not clear.
1
Research Assistant, Dept. of Civil and Environmental Engineering, The spreading inelastic zone elements proposed to date lack a
Univ. of California, Berkeley, CA 94720-1710. consistent numerical integration framework to assess their accu-
2
Professor, Dept. of Civil and Environmental Engineering, Univ. of racy under general load conditions.
California, Berkeley, CA 94720-1710 共corresponding author兲. E-mail: In recent years distributed inelasticity elements based on the
filippou@ce.berkeley.edu
force formulation have gained attention, because this approach
Note. This manuscript was submitted on May 2, 2008; approved on
April 3, 2009; published online on April 23, 2009. Discussion period
offers excellent accuracy of inelastic response with a single ele-
open until April 1, 2010; separate discussions must be submitted for ment discretization of the structural member. The consistent for-
individual papers. This paper is part of the Journal of Structural Engi- mulation of such an element originates with the work of Ciampi
neering, Vol. 135, No. 11, November 1, 2009. ©ASCE, ISSN 0733-9445/ and Carlesimo 共1986兲 with subsequent elaborations and refine-
2009/11-1310–1319/$25.00. ments by Spacone et al. 共1996b兲, Petrangeli and Ciampi 共1997兲,

1310 / JOURNAL OF STRUCTURAL ENGINEERING © ASCE / NOVEMBER 2009

J. Struct. Eng. 2009.135:1310-1319.

 and Neuenhofer and Filippou 共1997兲. With the implementation in y v2 v1 q3
the widely used simulation platform OpenSees 共http://opensees. x
q2 q1
berkeley.edu/index.php兲, such elements are now available for use
in performance-based earthquake engineering. While this type of (a) L v3
element is an excellent choice for the response simulation of y
structural members that may experience inelastic deformations at
several locations along the span, e.g., under distributed element M
loading, they are not computationally efficient for members with z
εa N
inelastic deformations only at the ends. In such a case, a distrib- κ
uted inelasticity element behaves like an element with an inelastic
Downloaded from ascelibrary.org by INDIAN INSTITUTE OF TECHNOLOGY on 06/05/14. Copyright ASCE. For personal use only; all rights reserved.

zone of fixed length at each end, if only the integration point (b) Cross−Section Strain Diagram Section Forces
closest to the end of the member experiences yielding. The length
of the end inelastic zone is equal to the integration weight of the Fig. 1. 2d distributed inelasticity frame element
monitoring point and varies with the numerical integration
scheme and the number of integration points. Increasing the num-
ber of integration points to account for the spread of inelastic consist of the axial strain at the reference axis ␧a and the curva-
deformations is computationally wasteful, unless the element is ture ␬. The corresponding section forces s are the axial force N
subdivided into an inelastic region at each end, and an elastic and the bending moment M.
region in between. Even so, several integration points are required In the force formulation the section forces s are derived from
in each inelastic region. Addessi and Ciampi 共2007兲 propose the element forces q by satisfying the element equilibrium in the
using three integration points in each inelastic zone for a total of undeformed configuration
seven integration points in the element. Even so, it is possible that
s共q,x兲 = b共x兲q + sw共x兲 共1兲
only the integration point closest to the end of the member expe-
riences plastic deformations over a significant range of element where b = force interpolation functions and sw = section forces
deformations, particularly for small strain hardening ratios. under the element loads in the basic system without rigid body
Under strain-softening response the distributed inelasticity el- modes. The force interpolation functions b consist of a constant
ement suffers from lack of objectivity unless measures are taken polynomial for the axial force and of linear polynomials for the
to normalize the material or response to a reference value of bending moment
fracture energy 共Coleman and Spacone 2001兲. This problem led

冤 冥
Scott and Fenves 共2006兲 to the formulation of an efficient beam 01 0
element with fixed inelastic zone length for softening response. b共x兲 = x x 共2兲
0 −1
This study is the first to address the issue of computational effi- L L
ciency within the framework of a numerically consistent integra-
tion scheme: a single monitoring point is used at the end of each The compatibility relation between the element deformations v
inelastic zone whose integration weight is derived from the and the section deformations e is linear under small deformations
Gauss-Radau integration scheme. The study, however, does not
冕
L
address strain hardening response. v= bT共x兲e共x兲dx 共3兲
The element proposed in this paper aims to address the limi- 0

tations of distributed inelasticity elements, while retaining the The section deformations e and the section forces s are related by
computational efficiency of a fixed inelastic zone length model by the section constitutive law, which can be expressed in terms of
monitoring the inelastic section response at only one location of the independent fields of the problem in the following form
the end inelastic zone of the member. A numerically consistent
integration scheme is proposed for the transformation of the sec- ŝ关e共x兲兴 − s共q,x兲 = 0 共4兲
tion to element deformations, while the inelastic zone length 共and
Various aspects of the implementation of the force formulation
thus the integration weight of the inelastic section deformations兲
in a general purpose finite element analysis program are discussed
is variable during the response time history. The element is based
by Spacone et al. 共1996a兲, Neuenhofer and Filippou 共1997兲, Tay-
on the force formulation which furnishes the framework for its
lor et al. 共2003兲, and Lee and Filippou 共2009兲.
consistent derivation and numerical implementation.
An important aspect of the efficiency of the computer imple-
mentation is the numerical evaluation of the integral in Eq. 共3兲.
This is expressed as
Force Formulation of Beam Element
Np

The discussion in this paper is limited to two dimensional Euler- v= 兺

i=1
共bTe兲兩x=␰ wi i
共5兲
Bernoulli beams with bilinear moment-curvature relation. The ex-
tension of the concepts to beams with integration of the material where N p = number of integration points; ␰i = integration locations;
response over the cross section is possible with further refine- and wi = weights of the numerical scheme. The challenge is to
ments, but is beyond the scope of this paper. The basic forces q select a consistent and accurate numerical integration scheme
and deformations v of a two-dimensional 共2d兲 frame element with the fewest number of integration points.
without rigid body modes are shown in Fig. 1共a兲. The section If little is known a priori about the curvature distribution in
forces and deformations are shown in Fig. 1共b兲. The element Eq. 共5兲, except that the end points of the element are of particular
length is denoted with L. importance, then the Gauss-Lobatto integration scheme offers the
Under the assumption that plane sections remain plane and highest accuracy for the smallest number of integration points. If
normal to the axis after deformation, the section deformations e information about the curvature distribution is available, it is pos-

JOURNAL OF STRUCTURAL ENGINEERING © ASCE / NOVEMBER 2009 / 1311

J. Struct. Eng. 2009.135:1310-1319.

 Node 1 Node 2 e = ee + ep 共6兲
Inelastic Zone Elastic Zone Inelastic Zone e
The linear elastic section deformations e can be readily obtained
LP,1 LP,2
from the section forces s with the help of the elastic section flex-
L
ibility matrix fse according to
Fig. 2. Both elastic and inelastic zones of element
ee = fses 共7兲
Substituting Eqs. 共6兲 and 共7兲 and making use of Eq. 共1兲 without
sible to take advantage of its characteristics in devising an opti- the element load term 共sw = 0兲 gives

冕 冕
mum integration scheme, as will be shown in the following.
Downloaded from ascelibrary.org by INDIAN INSTITUTE OF TECHNOLOGY on 06/05/14. Copyright ASCE. For personal use only; all rights reserved.

L L
v= bT共x兲fse共x兲b共x兲dxq + bT共x兲ep共x兲dx 共8兲
0 0
Element Response with Spread of Plasticity
With the elastic element flexibility fe defined according to

冕
In moment resisting frames that are subjected to high lateral L
forces, inelastic deformations concentrate at the ends of the mem- fe = bT共x兲fse共x兲b共x兲dx 共9兲
bers, as shown in Fig. 2. In metallic structures the inelastic be- 0
havior is characterized by strain hardening and the gradual spread
of inelasticity into the member, as long as local buckling or frac- the element deformations are written as
ture is prevented. The same is true in well designed reinforced

冕
L
concrete members as long as the axial compression remains
below the balance point. Strain-softening response is common in v = f eq + bT共x兲ep共x兲dx 共10兲
0
older RC elements with poor reinforcement detailing and in col-
umns under very high levels of axial force. Under strain-softening Identifying the first term in Eq. 共10兲 as the elastic element defor-
response, a damage zone forms at the ends of the member. mations ve
In determining the element deformations the length of the in-
elastic zone plays an important role. This length is denoted with v e = f eq 共11兲
L p. Under strain-softening conditions the length of the zone may
and the second term in Eq. 共10兲 as the plastic element deforma-
be derived from considerations of objectivity of the analytical
tions vp
response, as proposed by Coleman and Spacone 共2001兲 and Scott
and Fenves 共2006兲. Alternatively, empirical formulas from experi-
冕
L
mental evidence have been proposed, e.g., by Paulay and Priest- vp = bT共x兲ep共x兲dx 共12兲
ley 共1992兲. For strain hardening response the spread of the 0
inelastic zone is related to the increase of the internal forces at the
ends of the member. More generally, the length of the inelastic gives the decomposition of the element deformations v into the
zone can be expressed in terms of response variables and is, elastic contribution ve and the plastic contribution vp
therefore, time dependent, i.e., L p共t兲.
v = ve + vp 共13兲
In the presence of element loads, an additional term 兰L0 bTfseswdx
Spreading Inelastic Zone Element „SIZE… appears on the right hand side of Eq. 共13兲. Since this term is
constant, it does not change the generality of the derivation. Con-
To address the limitations of concentrated plasticity frame ele- sequently, it is left out of the following discussion for the sake of
ments or frame elements with fixed inelastic zone lengths a new brevity in the resulting expressions.
element is proposed in this paper. The basic idea of the element is Splitting off the linear elastic section contribution from the
to make optimum use of only two integration points at the ends of element deformations has the advantage of uncoupling the evalu-
the frame element while allowing the inelastic zone to spread ation of plastic element deformations, which arise exclusively at
during the response time history. The element’s name is accord- the element ends from the linear elastic response of the frame
ingly SIZE, which stands for spreading inelastic zone at element element. The latter can be evaluated independently with a suitable
ends. The integration scheme is designed to be numerically con- integration scheme, or may be available in closed form for special
sistent and to achieve great accuracy for hardening element re- cases. All that is required is that an elastic contribution can be
sponse, while being as efficient as any concentrated inelasticity identified in the section response. Splitting off the linear elastic
element with two monitoring sections. The element also improves section contribution is also advantageous during the elastic un-
the accuracy and consistency of the response under softening con- loading and reloading response to a general load history.
ditions. With the elastic contribution to the element deformations ac-
counted for the focus of developing an accurate integration
scheme for the section deformations is on the plastic curvatures,
Determination of Element Deformations
which arise at the element ends. A typical plastic curvature dis-
A key idea in the development of the new element is to split off tribution for a beam with bilinear moment-curvature relation
the elastic section deformations similar to the additive strain de- under monotonic loading is shown in Fig. 3. The integration of
composition of small strain plasticity theory. Thus, the section the plastic curvatures according to Eq. 共12兲 is then the focus of
deformations e are decomposed into the linear elastic contribution the following discussion.
ee and the plastic contribution ep so that The integral in Eq. 共12兲 is split into two separate integrals

1312 / JOURNAL OF STRUCTURAL ENGINEERING © ASCE / NOVEMBER 2009

J. Struct. Eng. 2009.135:1310-1319.

 0.6
Plastic curvature at t = 1.68 sec
κp1 0.5 Plastic curvature at t = 1.7 sec
Lp,2 Plastic curvature increment between t = [1.68,1.7] sec
0.4

0.3

max
Lp,1

p
or ∆κ /∆κ
κp2 0.2

p
0.1
Fig. 3. Plastic curvature ␬ p extracted from total curvature ␬

max
0

p
κ /κ
p
−0.1

冕 冕
L p,1共t兲 L
Downloaded from ascelibrary.org by INDIAN INSTITUTE OF TECHNOLOGY on 06/05/14. Copyright ASCE. For personal use only; all rights reserved.

−0.2
vp = bT共x兲ep共x兲dx + bT共x兲ep共x兲dx 共14兲 −0.3
0 L−L p,2共t兲
−0.4
0 0.2 0.4 0.6 0.8 1
where the subscripts 1 and 2 refer to the inelastic zone at the left x/L
and right element end, respectively. This split highlights the fact
that the plastic element deformations arise from two separate seg- Fig. 4. Plastic curvature increment distribution of a beam element
ments that can be integrated without regard for the intervening with bilinear moment-curvature
segment of zero plastic deformations.
Under general cyclic load conditions, the distribution of the
plastic section deformations ep in each end inelastic zone may be
simple one point integration suffices to give results of very good
rather irregular and may not follow in general the smooth trian-
accuracy for the smooth distribution of section deformation incre-
gular distribution of Fig. 3. Such an irregular distribution is im-
ments over the PDI zone at each element end.
possible to evaluate accurately with only one integration point.
With only one integration point in the PDI zone at each ele-
On the other hand, the distribution of plastic deformation rates at
ment end, Eq. 共17兲 becomes
the element ends is much smoother and offers an efficient solution
of the problem. Fig. 4 shows the plastic curvature distribution 2
during a two time step sequence of a frame element under a
general load history. The figure shows clearly that the distribution
v̇p = 兺
i=1
bTi ėpi wi共t兲 共18兲
of plastic curvature increments is rather smooth, even though the
distribution of plastic curvatures is quite irregular. Moreover, the where the subscripts 1 and 2 refer to the PDI zone in the left and
figure reveals an additional length parameter that plays an impor- right end inelastic zone of the element, respectively. Since the
tant role in the proposed element: the length of the zone of plastic integration interval at each end inelastic zone is time dependent,
deformation increments 共PDIs兲, which is the portion of the end the integration weight in Eq. 共18兲 is time dependent. By contrast,
inelastic zone with nonzero section deformation rates ėp. The the monitoring or integration point location remains the same
length of the PDI zone is denoted with l p. In the proposed element over the response time history, since it is essential that the mate-
the PDI zone is always assumed to stretch to the element end. In rial history be traced at the same physical location over the entire
contrast to L p, which is monotonically increasing with the ele- response time history.
ment response, the length l p varies between zero, when the ele- Using the backward Euler method to numerically integrate Eq.
ment undergoes elastic loading and unloading, and L, when 共18兲 in the time domain gives
plastic section deformations spread over the entire element.
2
To take advantage of the relatively smooth distribution of PDIs
over the PDI zone of the element, the integral of the plastic ele- ⌬vp = 兺 bTi ⌬epi wn+1,i 共19兲
ment deformations vp in Eq. 共14兲 is converted to an integral of i=1

deformation rates over the space domain where ⌬共 · 兲 = 共 · 兲n+1 − 共 · 兲n and the subscripts n and n + 1 denote the

冕
L variable values at time steps tn and tn+1, respectively.
v̇p共t兲 = bT共x兲ėp共x,t兲dx 共15兲 There now remains to select the integration point location and
0 the integration weight of the PDI zone for maximum accuracy.
This selection is governed by the following considerations: 共1兲 the
Which is subsequently integrated over the time domain to give plastic deformation history needs to be monitored at the ends of
the plastic element deformations vp according to the frame element where the bending moments are largest in mo-

冕 ment resisting frames subjected to high lateral forces; 共2兲 since

t
vp共t兲 = vp共t0兲 + v̇p共␶兲d␶ 共16兲 the PDIs exhibit a relatively smooth distribution over the PDI
t0 zone and can be approximated by a uniform distribution, as
with t0 denoting the starting point of time integration. According shown in Fig. 5, the integration scheme needs to be exact for a
to the earlier discussion for the plastic element deformations the uniform distribution over a zone of variable length. The second
integral in Eq. 共15兲 is split into separate integrals requirement suggests that the midpoint integration method is most
appropriate, which however brings up the issue of the location of

冕 l1p共t兲
冕
L
the monitoring point: it is not possible to select the midpoint of
v̇p = b共x兲Tėp共x兲dx + b共x兲Tėp共x兲dx 共17兲 the variable PDI zone, since it violates the first requirement. With
0 L−l2p共t兲
the recognition that the distribution of the curvature increments is
with each integral extending over the PDI zone at the correspond- uniform over the PDI Zone a new midend point integration
ing element end. The significant advantage of recasting the deter- 共MEPI兲 scheme is therefore selected that satisfies both of the
mination of element deformations in the form of Eq. 共17兲 is that a aforementioned requirements. With the selection of this custom

JOURNAL OF STRUCTURAL ENGINEERING © ASCE / NOVEMBER 2009 / 1313

J. Struct. Eng. 2009.135:1310-1319.

 ∆κp To conclude the determination of the element deformation in-
∆κp crements, the selection of a suitable reduction factor ␤ needs to be
Plastic Curvature Increment Diagram
lpn+1,1 lpn+1,2
discussed. With the typical distribution of PDIs varying between
β∆κ 1
p linear and trapezoidal, as shown in Fig. 5, the value of ␤ lies
between 1/2 and 1. In numerical simulations of structural re-
β∆κ 2
p
Approximate Plastic Curvature Increment Diagram sponse to time dependent excitations, the time step is usually
small enough so that the change of the PDI zone length is a small
Fig. 5. Approximate distribution of plastic curvature increment ⌬␬ p fraction of the PDI zone length. Consequently, the distribution of
the deformation increments is nearly uniform and the value of ␤
is closer to 1 than 1/2. In the validation studies of the proposed
Downloaded from ascelibrary.org by INDIAN INSTITUTE OF TECHNOLOGY on 06/05/14. Copyright ASCE. For personal use only; all rights reserved.

integration scheme for the proposed element the evaluation of the element a value of 7/8 gave the most satisfactory results and is
plastic element deformation increments in Eq. 共19兲 consists of the used throughout the numerical studies.
following steps: With the deformation increments according to Eq. 共21兲 the
1. The plastic curvature increments in the PDI zone are ap- total element deformations are determined by the simple expres-
proximated by a uniform distribution with the introduction of sion
a reduction factor ␤ that multiplies the end curvature values,
as shown in Fig. 5; with this approximation Eq. 共19兲 be- 2
comes vn+1 = vn + f ⌬q + ␤
e
兺
i=1
T
bmi ⌬epi ln+1,i
p
共23兲
2

⌬vp = ␤ 兺
i=1
bTi ⌬epi wn+1,i 共20兲 One limitation of the proposed element was already noted: the
PDI zone is assumed to stretch to the element end in each inelas-
2. The plastic curvature increments in each PDI zone are moni- tic end zone, since the plastic section deformation increments are
tored at the end point of the corresponding inelastic zone. monitored at the end of the element. Thus, a PDI zone forming
p
3. The integration weight is equal to the current length ln+1,i of inside the inelastic end zone of the element is not possible for
the PDI zone at the corresponding end. lack of a monitoring point. The validation studies, however, re-
4. The interpolation functions bi for the transformation of the veal that this limitation is not serious, since the case either rarely
PDI zone deformation increment ⌬epi ln+1,i
p
to the element de- arises in the response history of frame elements, or the error from
formation increments are evaluated at the midpoint of the this assumption is sufficiently small.
PDI zone.
With these steps the element deformation increments in Eq. Determination of PDI Zone Length
共20兲 become
The determination of the PDI zone length depends on whether the
2
section response exhibits hardening or softening with increasing
⌬vp = ␤ 兺
i=1
T
bmi ⌬epi ln+1,i
p
共21兲 deformation.

where bm⫽force interpolation matrix b at the midpoint of the Deformation Softening Response
corresponding PDI zone which is denoted with subscript i = 1 or 2 For softening response with increasing deformation the PDI zone
for the left or the right end of the element, respectively. length l p is presently assumed to be fixed and time independent.
However, the model can be extended to accommodate the evolu-

冉 冊 冤 冥
p
1 0 0 tion of the inelastic zone as function of response parameters with
ln+1,1 p p the use of one or more damage variables.
bm1 = b x = = ln+1,1 ln+1,1
2 0 −1 Under the assumption of time independence the PDI zone
2L 2L
length can be related to the inelastic zone length L p. The latter can
be obtained either from empirical relations, or from the expres-

冉 冊 冤 冥
p
1 0 0 sion
ln+1,2
bm2 = b x = L − = l p
lp 共22兲
2 0 − n+1,2 1 − n+1,2
2L 2L L p = 0.08L + 0.022f ydb 共kN,mm兲 共24兲
It is worth noting that the formulation for the element defor-
mation increments in Eq. 共21兲 leads to the dependence of the recommended for reinforced concrete elements by Paulay and
plastic rotation increment at each end from the plastic curvature Priestley 共1992兲, where f y and db are the yield strength and diam-
increments at both ends; this is due to the evaluation of the inter- eter of the longitudinal reinforcing bars, respectively.
polation functions bmi at the midpoint of the PDI zone at each The relation between plastic deformations and time indepen-
end. Without this coupling term the accuracy of the determination dent PDI zone length l p can be obtained from Eq. 共21兲 after drop-
of the element deformation increments ⌬vp decreases apprecia- ping the subscript n + 1 for the time step
bly, particularly when the plastic curvature is zero at one end, as
is the case in cantilever bridge piers. Inelastic zone elements with
2
end point evaluation, such as the modified Gauss-Radau quadra-
ture of the inelastic zone element by Scott and Fenves 共2006兲, vp = ␤ 兺
i=1
T p p
bmi ei li 共25兲
suffer from this limitation. This aspect of the proposed element is
relevant in the consistent evaluation of softening response, as will
be discussed in the validation studies. which, in detail, gives

1314 / JOURNAL OF STRUCTURAL ENGINEERING © ASCE / NOVEMBER 2009

J. Struct. Eng. 2009.135:1310-1319.

 −
M1 Mp
MP
Bending Moment Diagram Inelastic Zone
LP,2 lp1 −
Mp
Node 1 Node 2
Elastic Zone +
Mp
LP,1
Inelastic Zone MP +
LP,1 Mp
M2 LP,2
L Mp lines
Bending Moment Diagram
Fig. 6. Calculation of L p by comparing bending moment diagram
Downloaded from ascelibrary.org by INDIAN INSTITUTE OF TECHNOLOGY on 06/05/14. Copyright ASCE. For personal use only; all rights reserved.

with section capacity for monotonic loadings Fig. 7. Calculation of l p by comparing bending moment diagram
with section capacity for cyclic loading

冤 冥
␧a,1
p p
l1 + ␧a,2
p p
l2 Validation Studies for SIZE Element

冉 冊 冉 冊
冤冥
v1p l1p l2p
␬1pl1p − 1 − ␬2pl2p The proposed element is validated by studying the inelastic re-
v2p =␤ 2L 2L 共26兲

冉 冊 冉 冊
sponse of a beam under characteristic monotonic and cyclic load-
v3p l1p l2p ing conditions. The beam is simply supported and represents the
␬1pl1p + ␬2pl2p 1 −
2L 2L basic system without rigid body modes. The nonlinear geometric
transformation of the node variables according to the corotational
Eq. 共26兲 shows the dependency of the end rotations on the plastic formulation introduces the effect of nonlinear geometry under
curvature at both element ends, which is missing from the ele- large displacements independently from the basic element formu-
ment with fixed plastic hinge length 共Scott and Fenves 2006兲 lation, as discussed in Filippou and Fenves 共2004兲.
based on the modified Gauss-Radau integration method. The beam is assumed to have linear elastic axial behavior with
stiffness EA. The axial behavior is assumed to be uncoupled from
Deformation Hardening Response the flexural behavior. A bilinear moment-curvature relationship
For hardening response with increasing curvature the PDI zone with initial stiffness EI governs the flexural response. The post-
length l p is time dependent. It can be related to element response yield stiffness is set at 3% of the initial stiffness for strain hard-
variables. The most straightforward approach determines the PDI ening and strain-softening response. Only kinematic hardening is
zone length l p from the element forces q and the plastic capacity included in the analyses without limitation of the generality of the
of the cross section. Fig. 6 shows the relation between the end results.
forces and the bending moment distribution in the absence of For the strain hardening case the response of the beam with the
element loads. By comparing the bending moment with the plas- SIZE element is compared with the exact analytical response
tic moment capacity M p of the section the inelastic zone length under monotonic loading. The latter is obtained by integration of
can be readily established, as earlier spreading inelastic zone the bilinear moment-curvature relation. Under cyclic loading con-
models have also proposed 共Soleimani et al. 1979; Meyer et al. ditions, the response of the beam with one SIZE element is com-
1983兲. Under monotonic response the PDI zone length is equal to pared with the response of the beam using 10 distributed
the inelastic zone length, and it can be readily concluded from inelasticity elements with five Gauss-Lobatto integration points in
Fig. 6 that each. The distributed inelasticity element is also based on the
force formulation 共Spacone et al. 1996a兲 and is subsequently

冉 冊
identified as FBDI5 element. The response of the beam with one
1 SIZE element is also compared with the response of the beam
lip = L p,i = max 共qi+1 − sign共V兲M p兲,0 , V ⫽ 0, i = 1,2
V using a single FBDI5 element, so as to reveal the limitations of
共27兲 the latter approach, which amounts to using a fixed inelastic zone
length for the common case that inelasticity only arises at the end
where V = element shear force given by V = 共q2 + q3兲 / L. The sim- monitoring points.
plicity of the diagram in Fig. 6 does not account for cyclic loading For the strain-softening case, the SIZE element reproduces the
conditions. Under cyclic response the plastic moment capacity exact response for a given inelastic zone length. The latter is
needs to be adjusted in the inelastic zone to account for isotropic obtained from empirical formulas or experimental calibration and
and kinematic hardening of the section response. Instead of a line presently remains fixed during the response time history. The ef-
parallel to the element axis a piecewise linear capacity line is fect of a time dependent damage zone will be considered in future
introduced that is continuously updated during the response time studies. To demonstrate the benefit of the proposed MEPI integra-
history 共see Fig. 7兲. Moreover, the diagram in Fig. 6 does not tion method, the response of the element under softening condi-
account for the effect of a variable axial force on the moment tions is compared with the inelastic zone element by Scott-
capacity of the section. The proposed element is, therefore, cur- Fenves.
rently limited to inelastic members with a small variation of axial In the following studies the reduction factor ␤ of the SIZE
force during the response time history, as is the case in girders element is consistently selected equal to 7/8, since this value was
and columns of low to midrise moment resisting frames or bridge found to give the best agreement between the model and the exact
piers. The extension of the proposed element to account for the solution for a range of curvature distributions under monotonic
interaction between axial force and bending moment during the and cyclic load conditions. These studies reveal that the reduction
response time history is under progress. factor is not sensitive to the loading history.

JOURNAL OF STRUCTURAL ENGINEERING © ASCE / NOVEMBER 2009 / 1315

J. Struct. Eng. 2009.135:1310-1319.

 1.5 1.5
Exact α =1
v
SIZE 0.25
1.4 α =1 1
v
0.25 0
1.3 0
0.5
−0.25

p
p

q /M
q /M 1.2

2
3

0
1.1 −0.5 −0.25

−0.75
αv = −1 −0.5
1
Downloaded from ascelibrary.org by INDIAN INSTITUTE OF TECHNOLOGY on 06/05/14. Copyright ASCE. For personal use only; all rights reserved.

−0.5
−0.75 Exact
0.9 −1 α = −1 SIZE
v

0 1 2 3 4 5 6 7 −5 0 5
(a) v3 / vy (c) v2 / vy

1.5 1.5
Exact α =1
v
FBDI5 0.25
1.4 αv = 1 1
0.25 0
1.3 0
0.5
−0.25

p
q3 / Mp

q /M
1.2

2
0
1.1 −0.5 −0.25

−0.75
αv = −1 −0.5
1
−0.5
−0.75 Exact
0.9 −1 α = −1 FBDI5
v

0 1 2 3 4 5 6 7 −5 0 5
(b) v /v (d) v /v
3 y 2 y

Fig. 8. Comparison of moment-rotation relationships of the SIZE element, the FBDI5 element, and the exact solution for hardening ratio= 3%
and fixed rotation ratios

Strain-Hardening Response lent for Node 2, which yields first for all rotation ratios. The
accuracy of the SIZE element is less satisfactory for single cur-
Monotonic Loading vature conditions with rotation ratios from ⫺0.5 to ⫺1 in Fig.
The accuracy of the proposed SIZE element is validated by com- 8共a兲. The inability of the SIZE element to represent the strength
paring the inelastic response of a beam with the exact solution for softening for increasing rotation values is noteworthy. This can be
a range of rotation ratios ␣v = v2 / v3 and moment ratios ␣m attributed to the formation of the PDI zone in the interior of the
= q2 / q3 from ⫺1 to 1 in increments of 0.25 units, where v2 and q2 inelastic zone of the element. Since the SIZE element assumes
denote the rotation and moment at the left end, and v3 and q3 the that the PDI zone always stretches to the element end, this soft-
rotation and moment at the right end. The ratio of ⫺1 corresponds ening behavior cannot be captured. Fig. 8共c兲 shows that the accu-
to uniform curvature, while the ratio of 1 corresponds to antisym- racy is also quite remarkable for Node 1, which trails Node 2 in
metric double curvature. Only for the case of uniform curvature
its inelastic response and experiences much smaller plastic curva-
the PDI zone length of the SIZE element is fixed at l p = L / 2,
tures.
where L is the length of the beam. Since there is uniform yielding
In contrast to the SIZE element the beam with a single FBDI5
over the entire element, the value of the reduction factor ␤ should
element is unable to capture the gradual spread of plastic defor-
be set equal to 1. To maintain consistency with the other curvature
distribution cases, the following results use the factor of 7/8, since mations into the beam under double curvature conditions. This is
it is expected that the beam will experience different curvature evident from the constant postyield slope of the moment-rotation
distributions during a dynamic response analysis. relation in Fig. 8共b兲: such response is typical of an element with
Fig. 8 compares the exact moment-rotation relation at the ends fixed inelastic zone length, as is the case with the FBDI5 element
of a beam element with that of the SIZE element and with that of with plastic deformations taking place only at the end sections of
a single FBDI5 element under increasing end rotations with fixed the element. The accuracy of the FBDI5 element is significantly
ratio. In this case the point of inflection shifts during loading. The better for single curvature conditions, since in this case several
end rotation values in Fig. 8 are normalized to the yield rotation integration points experience plastic deformations and the inelas-
under antisymmetric double curvature vy = M pL / 6EI. Fig. 8 tic zone length is no longer constant with response history. In this
shows the remarkable accuracy of the SIZE element in represent- case the FBDI5 element shows better agreement with the exact
ing the moment-rotation response of the beam under spreading solution than the SIZE element. It is, however, important to note
inelastic deformations. Fig. 8共a兲 shows that the agreement be- that single curvature conditions with end rotation ratios from
tween the proposed SIZE element and the exact solution is excel- ⫺0.5 up to the uniform curvature case with a rotation ratio of ⫺1

1316 / JOURNAL OF STRUCTURAL ENGINEERING © ASCE / NOVEMBER 2009

J. Struct. Eng. 2009.135:1310-1319.

 1.5
Exact Exact αv = 1
1.8 SIZE
SIZE
1.4 αm = 1 1.7 0.25
0.75 1.6 0
1.3 0.5
0.25 1.5
−0.25

q3 / Mp
0
p

1.4
q /M

1.2
−0.25
−0.5
3

1.3 −0.5
1.1 −0.75
1.2
αm = −1 −0.75
1.1
1 αv = −1
Downloaded from ascelibrary.org by INDIAN INSTITUTE OF TECHNOLOGY on 06/05/14. Copyright ASCE. For personal use only; all rights reserved.

1
0.9 0.9
0 1 2 3 4 5 6 7 0 1 2 3 4 5 6 7
v /v v3 / vy
(a) 3 y

1.5 Fig. 10. Comparison of moment-rotation relationships of Node 2 of

α =1
m the SIZE element and the exact solution for hardening ratio= 10%
1 0.75

0.5 Cyclic Loading

0.5
The validation of the element under cyclic loading conditions is
0.25
q2 / Mp

shown in Fig. 12. In this case, the simply supported beam is

0 0 subjected to cycles of imposed end rotation values with a fixed
−0.25 rotation ratio ␣v = 0.5. The rotation history is such that the refer-
−0.5 ence rotation is first increased monotonically to a value equal to
−0.5
4.8 times the yield rotation of the beam and is subsequently re-
−0.75
Exact
−1 α = −1 SIZE
m

−5 0 5
v /v 16 Exact
(b) 2 y
SIZE
14 α =1
v
Fig. 9. Comparison of moment-rotation relationships of the SIZE
0.25
element and the exact solution for hardening ratio= 3% and fixed 12
0
moment ratios
10
κ / κy

8 −0.25

6
are rather unusual in moment resisting frame elements under the −0.5
effect of gravity loads and high lateral forces. Thus, in most prac- 4
tical situations the SIZE element is much more accurate than the −0.75
2
αv = −1
FBDI5 element, while also being significantly more economical
on account of the response evaluation at only two integration 0
0 1 2 3 4 5 6 7
points. (a) v3 / vy
Fig. 9 compares the exact moment-rotation relation at the ends
of a beam element with that of the SIZE element under increasing 0
0 Exact
end moments with fixed ratio in which case the point of inflection 2 −2 SIZE
αv = −1
remains fixed during loading. The results show that the agreement
between the proposed SIZE element and the exact solution is −4
1.5
excellent for both ends of the element and for all moment ratios. −6 0.25
To show that the reduction factor ␤ is not sensitive to the ratio
y

−0.75 −8
κ/κ

of the postyield stiffness to the initial stiffness of the moment- 1

curvature relationship, Fig. 10 compares the moment-rotation re- −0.5
−10
lationship of Node 2 for a hardening ratio equal to 10% and fixed 0.5
rotation ratios. The results show that the accuracy of the SIZE −12
αv = 1
element is as good as that for the smaller hardening ratio in Figs. −0.25 −14
0
8 and 9. These results are also obtained with a value of ␤ equal to
7/8, which is consequently the recommended value for the use of −16

the SIZE element in the case of bilinear moment-curvature re- −8 −6 −4 −2 0 0 2 4 6 8

(b) v2 / vy v /v
sponse. 2 y

Fig. 11 compares the curvature of the end sections of the beam

for the SIZE element and the exact solution under fixed rotation Fig. 11. Comparison of curvature-rotation relationships of the SIZE
ratios. The results show that the accuracy of the SIZE element is element and the exact solution for hardening ratio= 3% and fixed
excellent in representing the local response of the beam. rotation ratios

JOURNAL OF STRUCTURAL ENGINEERING © ASCE / NOVEMBER 2009 / 1317

J. Struct. Eng. 2009.135:1310-1319.

 1.5 Scott−Fenves
0.25 −0.25−0.5−0.75 αv = −1 SIZE
1
α =1
1 v

0.5
0.95

q2 / Mp
q3 / Mp

−0.5
0.9

−1
Downloaded from ascelibrary.org by INDIAN INSTITUTE OF TECHNOLOGY on 06/05/14. Copyright ASCE. For personal use only; all rights reserved.

10 FBDI5
FBDI5
−1.5 SIZE
0.85
−8 −6 −4 −2 0 2 4 6 8 0 1 2 3 4 5 6
v3 / vy (a) v /v
(a) 2 y

1.5 1
αv = 1
1 0.25
0.5
0.5 −0.25

q3 / Mp
p
q /M

0 0
2

−0.5
−0.5
−1 −0.5
10 FBDI5
FBDI5
αv = −1 Scott−Fenves
SIZE −1 −0.75 ↑ SIZE
−1.5
−8 −6 −4 −2 0 2 4 6 8 −8 −6 −4 −2 0 2 4 6 8
v /v (b) v /v
(b) 2 y 3 y

Fig. 12. Comparison of moment-rotation relationships of the SIZE, Fig. 13. Comparison of moment-rotation relationships of the SIZE
FBDI5, and 10 FBDI5 elements for ␣v = 0.5 and hardening ratio and the Scott-Fenves elements for a strain-softening beam subjected
= 3% to monotonic loadings with fixed rotation ratios

duced to a value of three times the yield rotation. This is followed restricts the effect of the plastic curvature at one end of the beam
by three symmetric cycles with peak rotation values of 6, 3, and element on the plastic beam rotation at that end only. This fact
7.2 times the yield rotation. With this load history the beam ex- needs to be taken into account when calibrating the length of the
periences the spread of plastic deformations, elastic unloading inelastic zone of the element with experimental data.
and reloading, and plastic deformations without increase of the
inelastic zone at the element ends. Fig. 12 shows that the pro-
posed SIZE element exhibits remarkable accuracy when com- Conclusions
pared to the reference solution obtained with 10 FBDI5 elements
for the beam. On the contrary, the solution obtained with one The paper presents a new beam-column element that combines
FBDI5 element exhibits relatively poorer accuracy. computational efficiency with accuracy. The element uses only
one monitoring section in each end inelastic zone of the structural
member, but accounts for the spread of inelastic deformations
Strain-Softening Response
under strain hardening response. It is, therefore, a spreading in-
For the validation of the proposed SIZE element under strain- elastic end zone model called SIZE that combines the benefit of
softening conditions the same monotonic and cyclic load history integrating the section moment-curvature response with the com-
is used. The negative postyield stiffness is equal to 3% of the putational efficiency of concentrated hinge models for beams and
initial stiffness of the moment-curvature relation. In this case the columns. Under strain-softening response the element uses a fixed
inelastic zone length of the element is constant. It is set equal to plastic hinge length, but accounts for the dependence of the end
20% of the element length. Since there is no exact solution in this rotations on the plastic curvatures at both element ends, which is
case, the inelastic response of the element is compared in Fig. 13 missing from a recently proposed element with fixed plastic hinge
with the recently proposed fixed inelastic zone element by Scott- length.
Fenves that is based on a modified Gauss-Radau integration of The accuracy of the element derives from the integration of the
inelastic deformations. This comparison intends to demonstrate section deformation rate over a PDI zone. With the use of the
the slight discrepancy between the two models in the representa- deformation rate instead of the total deformations a novel integra-
tion of the softening response of the beam element at either end. tion method called MEPI that combines the benefits of the mid-
It stems from the dependence of the plastic element rotations in point method with the end point monitoring of the plastic
the SIZE element on the plastic curvatures at both element ends. response is introduced with excellent accuracy.
By contrast, the integration scheme of the Scott-Fenves element The accuracy of the element is confirmed under monotonic

1318 / JOURNAL OF STRUCTURAL ENGINEERING © ASCE / NOVEMBER 2009

J. Struct. Eng. 2009.135:1310-1319.

 and cyclic loading conditions with different end moment and end and V. V. Bertero, eds., CRC, Boca Raton, Fla.
rotation ratios. It is exceptional for double curvature conditions, Giberson, M. 共1967兲. “The response of nonlinear multistory structures
but slightly less accurate for single curvature conditions. The lat- subjected to earthquake excitation.” Ph.D. dissertation, California In-
ter are, however, less relevant in dynamic response applications stitute of Technology 共CalTech兲, Pasadena, Calif.
of structural systems under earthquake loading. Hajjar, J. F., and Gourley, B. C. 共1997兲. “A cyclic nonlinear model for
The presentation of the SIZE element in this paper is limited to concrete-filled tubes. 1. Formulation.” J. Struct. Eng., 123共6兲, 736–
744.
bilinear moment-curvature response under fixed axial force. This
Hilmy, S. I., and Abel, J. F. 共1985兲. “Material and geometric nonlinear
is not a serious limitation for beams and columns under axial
dynamic analysis of steel frames using computer-graphics.” Comput.
force with small variation during the dynamic response history. Struct., 21共4兲, 825–840.
The extension of the proposed element to fiber section response Lee, C.-L., and Filippou, F. 共2009兲. “Frame elements with mixed formu-
Downloaded from ascelibrary.org by INDIAN INSTITUTE OF TECHNOLOGY on 06/05/14. Copyright ASCE. For personal use only; all rights reserved.

that will directly account for the interaction between axial force lation for singular section response and bifurcation.” Int. J. Numer.
and bending moment is currently in progress. Methods Eng., 78共11兲, 1320–1344.
Mahin, S., and Bertero, V. 共1976兲. “Problems in establishing and predict-
ing ductility in seismic design.” Proc., Int. Symp. on Earthquake
References Structural Engineering, EERI, St. Louis, Mo., 613–628.
Meyer, C., Roufaiel, M., and Arzoumanidis, S. 共1983兲. “Analysis of dam-
Addessi, D., and Ciampi, V. 共2007兲. “A regularized force-based beam aged concrete frames for cyclic loads.” Earthquake Eng. Struct. Dyn.,
element with a damage-plastic section constitutive law.” Int. J. 11共2兲, 207–228.
Numer. Methods Eng., 70共5兲, 610–629. Neuenhofer, A., and Filippou, F. C. 共1997兲. “Evaluation of nonlinear
Anagnostopoulos, S. 共1981兲. “Inelastic beams for seismic analysis of frame finite-element models.” J. Struct. Eng., 123共7兲, 958–966.
structures.” J. Struct. Eng., 107共7兲, 1297–1311. Paulay, T., and Priestley, M. J. N. 共1992兲. Seismic design of reinforced
Attalla, M. R., Deierlein, G. G., and McGuire, W. 共1994兲. “Spread of concrete and masonry buildings, Wiley, New York.
plasticity—Quasi-plastic-hinge approach.” J. Struct. Eng., 120共8兲, Petrangeli, M., and Ciampi, V. 共1997兲. “Equilibrium based iterative solu-
2451–2473. tions for the non-linear beam problem.” Int. J. Numer. Methods Eng.,
Ciampi, V., and Carlesimo, L. 共1986兲. “A nonlinear beam element for 40, 423–437.
seismic analysis of structures.” 8th European Conf. on Earthquake Powell, G. H., and Chen, P. F. S. 共1986兲. “3D beam-column element with
Engineering, Laboratorio Nacional de Engenharia Civil, Lisbon, 6.3/ generalized plastic hinges.” J. Eng. Mech., 112共7兲, 627–641.
73–6.3/80. Scott, M. H., and Fenves, G. 共2006兲. “Plastic hinge integration methods
Clough, R., Benuska, K., and Wilson, E. 共1965兲. “Inelastic earthquake for force-based beam column elements.” J. Struct. Eng., 132共2兲, 244–
response of tall buildings.” Proc., 3rd World Conf. on Earthquake 252.
Engineering, Vol. 2, New Zealand Society for Earthquake Engineer- Soleimani, D., Popov, E., and Bertero, V. 共1979兲. “Nonlinear beam model
ing, Auckland, New Zealand, 68–69. for r/c frame analysis.” Proc., Seventh Conf. on Electronic Computa-
Coleman, J., and Spacone, E. 共2001兲. “Localization issues in force-based tion, ASCE, New York, 483–509.
frame elements.” J. Struct. Eng., 127共11兲, 1257–1265. Spacone, E., Ciampi, V., and Filippou, F. C. 共1996a兲. “Mixed formulation
Filippou, F., and Issa, A. 共1988兲. “Nonlinear analysis of reinforced con- of nonlinear beam finite-element.” Comput. Struct., 58共1兲, 71–83.
crete frames under cyclic load reversals.” Rep. No. UCB/EERC-88/12, Spacone, E., Filippou, F., and Taucer, F. 共1996b兲. “Fiber beam-column
Earthquake Engineering Research Center, Univ. of California, Berke- model for nonlinear analysis of rc frames: I: Formulation.” Earth-
ley, Calif. quake Eng. Struct. Dyn., 25共7兲, 711–725.
Filippou, F. C., and Fenves, G. L. 共2004兲. “Methods of analysis for Taylor, R. L., Filippou, F. C., Saritas, A., and Auricchio, F. 共2003兲.
earthquake-resistant structures.” Earthquake engineering from engi- “Mixed finite element method for beam and frame problems.” Com-
neering seismology to performance-based engineering, Y. Bozorgnia put. Mech., 31共1–2兲, 192–203.

JOURNAL OF STRUCTURAL ENGINEERING © ASCE / NOVEMBER 2009 / 1319

J. Struct. Eng. 2009.135:1310-1319.