- Email: [email protected]
Inelastic mixed fiber beam finite element for steel cyclic behavior
Engineering Structures 106 (2016) 399–409
Contents lists available at ScienceDirect
Engineering Structures journal homepage: www.elsevier.com/locate/engstruct
Inelastic mixed fiber beam finite element for steel cyclic behavior I.A. Gkimousis, V.K. Koumousis ⇑ Institute of Structural Analysis & Aseismic Research, National Technical University of Athens, Zografou Campus, 15780 Athens, Greece
a r t i c l e
i n f o
Article history: Received 25 April 2015 Revised 26 August 2015 Accepted 12 October 2015 Available online 11 November 2015 Keywords: Steel cyclic behavior Hysteretic model Hellinger–Reissner variational principle Mixed beam formulation Nonlinear frame analysis
a b s t r a c t In this work a new hysteretic uniaxial steel model is determined to describe steel cyclic behavior, which is further implemented to derive a fiber beam–column element on the basis of Hellinger–Reissner principle. The proposed model maintains full memory of the loading path and evolves following a single nonlinear differential equation expressing the entire hysteresis. The element is capable of addressing the overshooting problem of the existing models which occurs during short reversals. The state determination of the proposed element is investigated numerically following the linearization of the derived equations. Numerical results are presented that validate the proposed approach and demonstrate its computational efficiency. Ó 2015 Elsevier Ltd. All rights reserved.
1. Introduction The state of a deformable body subjected to body forces, tractions and kinematic boundary conditions is considered fully defined when the displacements, stresses and deformations are determined at any point of the body. In particular, for earthquake engineering and structural analysis of skeletal structures, beam elements usually based on the Euler–Bernoulli theory assumption that plane sections remain plane and perpendicular to the deformed axis, are typically considered. This facilitates considerably state determination of such elements. Displacement based beam elements were initially used following the classical stiffness method in which displacements were the only considered independent field [1]. When cubic and linear shape functions are employed for the transverse and axial displacements respectively, the resulting displacement field leads to constant axial deformation and linear curvature, which however is not appropriate when plastic deformations occur. To address this deficiency a structural member should be discretized in more than one element at the expense of increasing computational cost. Also, equilibrium equations are only fully satisfied at element nodes, while within the element they are satisfied in weak form as they are not valid for all possible displacement fields that satisfy essential boundary conditions. To resolve this problem, force based models were proposed that interpolate nodal forces within the element maintaining ⇑ Corresponding author. Tel.: +30 2107721657; fax: +30 2107721651. E-mail addresses: [email protected] (I.A. Gkimousis), [email protected]. gr (V.K. Koumousis). http://dx.doi.org/10.1016/j.engstruct.2015.10.023 0141-0296/Ó 2015 Elsevier Ltd. All rights reserved.
equilibrium. These models were implemented in the framework of the stiffness method of structural analysis and in that respect they are considered ‘‘mixed” as they use both force and displacement fields as independent ones. One of the first consistent and general force based beam model was proposed by Spacone et al. [2] and was later simplified numerically by Neuenhofer and Filippou [3]. Although the force based method proved very efficient and is currently widely used, there were some concerns about its variational consistency that were resolved by Hjelmstad and Taciroglu [4]. Moreover, the same authors in [5] showed that it is possible to provide non-variationally consistent force-based elements within the ‘‘nonlinear flexibility” framework by enforcing equilibrium directly. Also, various local and global solution strategies originate from the variational structure of the mixed beam elements as described in [6,7]. Thereafter, mixed methods seem to dominate the research field of nonlinear beam problems and corresponding numerical procedures as they proved more efficient following also the work of Hjelmstad and Taciroglu [8], Taylor et al. [9], Alemdar and White [10], Alsafadie et al. [11] and Correia et al. [12]. In this context elastoplastic material models were represented in classical form relying on the notions of yield surface, flow rule and hardening parameters, while their incorporation in the state determination process in linearized form gave rise to the returnmapping algorithm [13]. However, cyclic behavior can also be modeled using hysteretic evolution differential equations. Following this approach, Simeonov et al. [14] developed a force based element where material constitutive relations are considered in rate form and are solved simultaneously with the global differential equations of motion in state-space form. Also, Jafari et al. [15] extended
400
I.A. Gkimousis, V.K. Koumousis / Engineering Structures 106 (2016) 399–409
this formulation in large displacement analysis following a displacement based formulation. In addition Triantafyllou and Koumousis [16,17] proposed a finite element procedure where material nonlinearity is treated constitutively at the element level through proper implementation of the Bouc–Wen hysteretic rule. In the case of elastoplastic steel cyclic behavior various models have been proposed among which the widely used Menegotto– Pinto model [18] originating from the generalization of Ramberg–Osgood model [19]. Moreover, for efficient modeling of the inelastic buckling of reinforcing steel bars under cyclic behavior, the Monti–Nuti model [20] is commonly used as an enhancement of the Menegotto–Pinto model, which accounts for four different hardening rules. Furthermore, inelastic buckling of reinforcing bars has been investigated in several other studies, i.e. [21–24], while also modeling of high strength structural steel is addressed in Refs. [25,26]. Although the core of several steel models is the Menegotto– Pinto model, it is characterized by an overshooting of the reloading branch after short reversals [27]. This feature was considered of minor importance and originates from an effort to reduce computational cost by truncating model’s memory. Although attempts were made to tackle this aspect by distinguishing the reversals as major and minor [28] or complete and incomplete [29] the problem can be also addressed from a different perspective. In this work a new small displacement fiber beam–column element is proposed, which incorporates a uniaxial steel model for cyclic loading that resolves the overshooting behavior after short reversals, being on the same time computationally efficient. To accomplish this, a hysteretic model is developed that maintains full memory of the loading path and evolves following a single differential equation expressing the entire hysteresis. This constitutive behavior is implemented in the general framework of the two-field Hellinger–Reissner formulation [30] resulting in a wellestablished state determination algorithm. This is further investigated numerically following the linearization of the equilibrium and compatibility element equations. The rest of the paper is organized as follows: In the first section the basic structure of the proposed model is derived from classical plasticity considerations incorporating kinematic hardening. Then, the modifications for modeling steel cyclic behavior are implemented and the comparison with the Menegotto–Pinto model is performed underlining the features of the proposed formulation. Cross-sectional constitutive equations are derived in the sequel from fiber integration and the two-field Hellinger–Reissner principle is used to derive the element equilibrium and compatibility equations. The standard linearization method is implemented for solving the element equations and the respective state determination process is discussed. Finally, three numerical examples are presented that verify the proposed beam element and demonstrate the performance of the hysteretic constitutive relations embedded into a two-field variational formulation for the inelastic analysis of steel frames. 2. Inelastic Euler–Bernoulli beam theory 2.1. Fiber plasticity A hysteretic model incorporates the entire inelastic loading path of a deformable body, namely elastic loading, yielding, hardening and unloading in a single nonlinear differential equation that embodies both the yield surface and hardening rule. Mathematically this addresses the entire evolution process without the need of incremental considerations [31]. In this context Sivaselvan and Reinhorn [32] based on Bouc–Wen model [33] proposed a hysteretic model in stress resultant terms derived explicitly from classical plasticity theory.
In this work the stress–strain constitutive law of an elastoplastic fiber-rod subjected to uniaxial tension is determined according to the strain decomposition rule in rate form as [34]:
r_ ¼ Eðe_ e_ p Þ
ð1Þ
where r_ is the rate of normal stress, e_ is the rate of total strain and e_ p is the rate of plastic strain. This relation indicates that stresses evolve proportionally to the evolution of the elastic strains. For the same fiber a yield function with linear kinematic hardening and the same initial yield stress in tension and compression rþy0 ¼ ry0 ¼ ry0 is expressed in the following form:
U r; b; ryo ¼ jr bj ry0 6 0 with b being the back stress. Plastic strain rate according to the flow rule:
e_ p ¼ k_
@U _ ¼ k sgnðr bÞ @r
ð2Þ _p
e is determined ð3Þ
where k_ P 0 is the plastic multiplier which is actually the magnitude of the strain rate, with the signum function defining its sense. Axial stress and plastic multiplier are restricted by unilateral constrains representing restrictions that signify whether the material has yielded or not, resulting from Karush–Kuhn–Tucker (KKT) optimality conditions. These are expressed in the following form:
k_ Uðr; b; ry0 Þ ¼ 0
ð4Þ
During plastic response ðk_ > 0; Uðr; b; ry0 Þ ¼ 0Þ the consistency condition is derived by differentiating Eq. (4):
_ ðr; b; ry0 Þ ¼ 0 ) U _ ðr; b; ry0 Þ ¼ 0 ) r_ ¼ b_ k_ U
ð5Þ
Also, when the fiber deforms elastically ðk_ ¼ 0; Uðr; b; ry0 Þ < 0Þ the sign of the rate of the yield function indicates loading or unloading _ < 0 respectively. _ > 0 and U with U Taking into account Eq. (5) meaning that the rate of the back stress during plastic loading is equal to the rate of the axial stress, the following relation is considered:
b_ ¼ H e_ p ¼ H k_ sgnðr bÞ
ð6Þ
where H is the hardening ratio, i.e. the slope of the stress–strain curve in plastic strain terms (Fig. 1). Substituting Eq. (1), (3) and (6) in Eq. (5) the following relation is obtained:
k_ ¼ sgnðr bÞ
E e_ EþH
ð7Þ
Furthermore, introducing a ¼ Et =E as the ratio of the post-yield tangent modulus over the elastic modulus, the relation between the hardening ratio H and tangent modulus Et is established as:
H¼
Et a E ¼ 1 Et =E 1 a
ð8Þ
Then Eq. (3), using relations (7) and (8) obtains the following simple form:
e_ p ¼ ð1 aÞ e_
ð9Þ
Eq. (9) holds only for the plastic deformation phase. During elastic loading or unloading plastic strain is not induced as k_ ¼ 0. This is controlled by the compact form of KKT conditions (4). To include all loading phases in a single relation, Eq. (9) is extended by considering two Heaviside type functions acting as switches as follows:
e_ p ¼ H1 H2 ð1 aÞ e_
ð10Þ
where H1 controls yielding and H2 controls loading/unloading states. Consequently the following distinct phases can be described in a single equation:
401
I.A. Gkimousis, V.K. Koumousis / Engineering Structures 106 (2016) 399–409
Fig. 1. Uniaxial bilinear constitutive law in terms of (a) total strain and (b) plastic strain.
H1 ¼ H2 ¼
0; elastic deformation ðe_ p ¼ 0Þ 1; plastic flow ðe_ p ¼ H2 ð1 aÞ e_ Þ
ð11Þ
0; unloading ðe_ p ¼ 0Þ 1; loading ðe_ p ¼ H1 ð1 aÞ e_ Þ
ð12Þ
Function H1 emerges from the yield function (2) by properly expressing the back stress b as a function of the calculated plastic strain according to the following relation [35]:
b¼He ¼ p
a
1a
E
ee
el
¼
a
1a
ðE e rÞ
ð13Þ
Taking into account Eq. (13), the yield function can be expressed as follows:
raEe 61 ð1 aÞ r
U ¼
ð14Þ
y0
which, in the limit case of yielding
r ¼ ry0 ; e ¼ ey0 sets U ¼ 1. If
the above relation is raised to the nth power the absolute value approximates the step function with discrete values [0, 1] as n ! 1. Consequently the following relation holds:
r a E e n H1 ð1 aÞ ry0
ð15Þ
Parameter n controls the form of transition from the elastic to plastic branch of the stress–strain law, which is smoother for lower values, approaching a bilinear behavior for higher ones ðn P 8Þ. Also the expression ðrh ¼ r a E eÞ refers to the hysteretic part of the normal stress. Moreover, function H2 aims at controlling loading and unloading following the sign of the yield function rate with the positive sign indicating loading and negative sign unloading. Furthermore, smoothening of the Heaviside function to yield zero value for unloading and unity for loading is handled as follows: _ 1 þ sgn U
_ 1 þ sgnðr bÞ sgnðr_ bÞ 1 þ sgnðr bÞ sgnðe_ Þ ¼ ) ¼ 2 2 2 H2 ¼ 0:5 þ 0:5 sgn r aðE e ry0 Þ e_ 0:5 þ 0:5 sgnððr a E eÞ e_ Þ H2 ¼
ð16Þ
In the above relation, b_ ¼ 0 for elastic response and the sign of the stress rate r_ is replaced by the sign of the strain rate e_ . After defining functions H1 and H2 , Eq. (10) is substituted in Eq. (1) and the final form of the constitutive relation is obtained:
r_ ¼ ½1 ð1 aÞ H1 H2 E e_ ¼ Et e_
ð17Þ
where Et is the tangent modulus. 2.2. Steel cyclic behavior Steel stress–strain relationship under cyclic loading primarily aims at expressing the linear elastic and hardening branches in the first excursion, together with the Bauschinger effect for subsequent cycles. Yield plateau [36], plasticity-damage coupling [37] and low cycle fatigue [38] are not considered in this work. To account for steel cyclic behavior Eq. (17) along with relations (15) and (16) for kinematic hardening are implemented. The initial bilinear constitutive behavior is obtained using a high value for the exponent parameter n ðn1 P 8Þ. The bilinear behavior is incorporated just for the initial yield point, while for subsequent yielding in the opposite direction gradual transition in the plastic regime is manifested. In the proposed hysteretic model this is accomplished by subsequently reducing parameter n to lower values n2 ðn2 < 2Þ. Hysteretic models have been criticized in the past for their incompatibility with Drucker’s or Ilyushin’s postulates of plasticity for partial unloading–reloading in the nonlinear path [39]. Based on previous work [40] an accurate solution was recently proposed that eliminates the problem by properly modifying the reloading path of the hysteretic spring for the Bouc–Wen model [41] or the Sivaselvan–Reinhorn model [42]. In this work the same approach is utilized, which leads to the modification of the switch function H2 in the following form:
2 ¼ H2 ð1 RÞ H
ð18Þ
where Rðrh ; e; ry0 Þ is a stiffening factor that controls stiffness recovery in the reloading branch, which is expressed as:
R ¼ 0:5 ½1 þ sgnðrh rr Þ sgnðrh Þ
er ec er e
ð19Þ
402
I.A. Gkimousis, V.K. Koumousis / Engineering Structures 106 (2016) 399–409
where rh is the hysteretic part of the normal stress ðrh ¼ r a E eÞ; rr and er are the stress and strain respectively corresponding to the last observed reversal point and ec is the modified strain in the unloading path that corresponds to the given hysteretic stress rh (Fig. 2):
ec ¼
rh rr ey0 þ er sgnðrh Þ ry0
ð20Þ
where also appropriate tracking of the reversal point is needed to accurately handle different types of short reversals as indicated in [41,42]. 2.2.1. Comparison with experimental data Experiments on steel rebars under cyclic response can be found in the work of Thompson and Park [43]. One of these experiments is presented in the following Fig. 3 along with the predicted response of the proposed model. Steel grade is S275, while post to pre-yield ratio is considered as a ¼ 0:01. For model parameters, initial and final exponents were selected as n1 ¼ 8 and n2 ¼ 0:9 respectively to achieve satisfactory agreement with observed experimental behavior. 2.2.2. Comparison with Menegotto–Pinto model To illustrate the performance of the proposed model, which incorporates the modifications made for modeling steel cyclic behavior, a direct comparison with the well-established Menegotto–Pinto model [18] is performed. For both models an external strain history presented in Table 1 is applied and model parameters are selected in such way as to induce similar outcome. More specifically, material modulus is E = 210 GPa, yield stress ry0 ¼ 235 MPa, post-yield to pre-yield ratio a ¼ 0:02, while for the proposed model the initial exponent is n1 ¼ 10 and the final exponent is n2 ¼ 0:2. The parameters of Menegotto–Pinto model are considered as: R0 ¼ 20, cR1 ¼ 18:5, cR2 ¼ 0:1. In Fig. 4 it is observed that the proposed model overcomes the overshooting of the Menegotto–Pinto model after the partial unloading branch. This overshooting originates from the restricted memory of the Menegotto–Pinto model as underlined by Filippou et al. [27]. To address this issue, Dodd and Restrepo–Posada [28] distinguished reversals as major, minor and simple, while Balan et al. [29] separated them into complete and incomplete, with a value of 2f y being the threshold parameter of their distinction. However, these modifications are perplexing the original model increasing bookkeeping of data and thus computational cost. On the other hand the proposed model imposes no restriction on memory returning exactly on the previous loading path as it is evident in Fig. 4.
Fig. 3. Comparison with experiment [43].
Table 1 Strain history values. Time (s)
Strain (–)
0 1 2 2.5 2.7 3.5
0 0.04 0.005 0.009 0.0085 0.04
Fig. 4. Comparison between the proposed and Menegotto–Pinto model.
3. Cross-sectional relations Cross-sectional constitutive relations are derived from integration of the local constitutive equations of all unit length fibers using midpoint integration. Especially for planar analysis and symmetrical cross sections, fibers can be treated in layers (Fig. 5). Following the Euler–Bernoulli beam theory the normal stress rðx; yÞ; the axial strain eðx; yÞ and the axial and transverse displacement ux ðx; yÞ; uy ðx; yÞ are considered. By considering the kinematic assumption that plane sections before deformation remain plane and normal to the elastic line after deformation, the axial displacement ux ðx; yÞ of any fiber can be expressed in terms of the displacement of the beam reference axis x as:
ux ðx; yÞ ¼ uðxÞ y
Fig. 2. Treatment of short reversals in the proposed formulation.
@wðxÞ @x
ð21Þ
where uðxÞ and wðxÞ are the axial and transverse displacements of the cross section at the reference axis. Thus the following relation holds for the axial strain:
I.A. Gkimousis, V.K. Koumousis / Engineering Structures 106 (2016) 399–409
403
the tangent stiffness matrix at every point of the nonlinear path during dynamic loading. Moreover, at the boundaries the three fields are manifested as nodal values. The standard six nodal displacements ¼ fu1 ; u2 ; r1 ; u4 ; u5 ; r2 gT , and rotations for the 2D beam element u are related with deformations q ¼ fd; h1 ; h2 gT i.e. axial deformation d along the reference axis and the chord rotations at both ends. In addition the stress resultants Q ¼ fN; M 1 ; M 2 gT are introduced as presented in Fig. 6. The relation between the nodal displacements and member deformations is established with the transformation matrix T, which removes rigid body motion as follows (Fig. 7):
2
; q¼T u
Element equilibrium under distributed applied
T bðxÞ ¼ bx ðxÞ; by ðxÞ in both directions x, y requires:
Fig. 5. Cross-sectional local axis and layer definition.
T
eðx; yÞ ¼ e0 ðxÞ y uðxÞ ¼ ½1 y dðxÞ ¼ l dðxÞ
ð22Þ T
where the deformation field dðxÞ ¼ fe0 ðxÞ; uðxÞg consists of the axial strain e0 ðxÞ and curvature uðxÞ at the reference axis. Euler–Bernoulli beam theory considers only the effects of axial stress r, thus by applying equilibrium conditions within the cross section, the stress resultants are evaluated as follows:
Z
Z
NðxÞ ¼ A
rðx; yÞdA; MðxÞ ¼
A
y rðx; yÞdA
ð23Þ
The same relations can be casted in matrix form as:
^ DðxÞ ¼
Z
A
½1 yT rðx; yÞdA
ð24Þ
^ where DðxÞ ¼ fN ðxÞ; MðxÞgT . By replacing Eq. (17) along with Eq. (22) in the rate form of Eq. (24) the ‘‘cross sectional constitutive relation” is derived:
_ ^_ DðxÞ ¼ kðxÞ dðxÞ with
Z
kðxÞ ¼
½1 yT Et ½1 ydA
ð25Þ
ð26Þ
A
Finally, the flexibility matrix of the cross section is set as the inverse of the cross-sectional stiffness.
f ðxÞ ¼ kðxÞ1
3 1 0 0 1 0 0 6 7 T ¼ 4 0 1=L 1 0 1=L 0 5 0 1=L 0 0 1=L 1
ð27Þ
4. Two field beam element variational formulation Mixed methods appear to dominate the classical displacement based formulation as they satisfy element equilibrium in strong form avoiding the need for further element discretization. The equilibrium error is of major importance, as internal forces are overestimated leading to inaccurate structural response. Two or three field methods based on Hellinger–Reissner or Hu–Washizu principles respectively are able to diminish this error successfully. In this work the two-field Hellinger–Reissner principle is employed, disregarding the inherent error in the deformation field which is considered herein of secondary importance. Having the inelastic constitutive equation for a beam fiber in a cross section via relations (17) and (18), the next step is to define the overall nonlinear response of the element. This refers to the displacement, deformation and stress fields across the beam length and
@ 2 MðxÞ by ðxÞ ¼ 0 @x2 @NðxÞ bx ðxÞ ¼ 0 @x
ð28Þ
loads
ð29Þ
A mixed formulation is presented in the sequel based on the two-field Hellinger–Reissner potential [30], where both displacement and stress fields are considered as independent. Classical stiffness analysis of frame structures is based on the continuity of the displacement field between elements. This is accomplished explicitly as the cubic interpolation functions satisfy the kinematic boundary conditions at the end nodes. In the pure force formulation based on the principle of virtual forces, only the stress field is independent and interpolation functions based on element equilibrium relate internal and nodal forces. For a single element analysis, force based method is straightforward as interpolation base functions that satisfy essential boundary conditions and external applied loading can easily be defined. However, in a redundant structure a set of redundant forces must be introduced [4] the selection of which though is difficult to automate. For these reasons stiffness formulation is usually the efficient way to express equilibrium at the structural level. In addition the advantages of the flexibility approach to impose element equilibrium in strong form and derive an exact stiffness matrix in the nonlinear range are explored at the element level. Hence, flexibility approach for the element state determination requires the independency of the stress field, while stiffness based formulation at the structural level suggests the displacement field as independent. Considering the two fields as independent the Hellinger–Reissner energy principle is employed, hence the cross sectional deformations are expressed in two ways; in terms of displacements dðuÞ and as functions of stress resultants dðDÞ. For a fiber beam under uniaxial deformation with stress r and strain e the Hellinger–Reissner functional [30] can be stated in terms of the two independent sets of variables, i.e. the normal stress r and axial displacement u as follows:
PHR ðr; uÞ ¼
Z
V
feðuÞ r vðrÞgdV Pext ðuÞ
ð30Þ
where vðrÞ is the complementary energy density function from which strains can be calculated in terms of stresses, Pext ðuÞ is the functional of external loading and integration is performed at the undeformed volume V of the element. Taking in account Eqs. (22) and (24) the above functional can be written also in stress resultant and deformation terms omitting for simplicity the argument in the integrals i.e. ðDðxÞ ¼ DÞ.
404
I.A. Gkimousis, V.K. Koumousis / Engineering Structures 106 (2016) 399–409
Fig. 6. Basic forces, displacements, deformations and distributed loads on the beam element.
nodal values. For the displacement field the classical cubic interpolation functions N u are considered while deformation interpolation functions N d result from the derivatives of N u (1st Xderivative for the axial strain and second derivative for curvature terms).
8 9 u1 > > > > > > > > > u2 > > > > > > > < ux ðx; yÞ r1 = ) uðxÞ ¼ N u ðxÞ u ¼ N u ðxÞ > uy ðx; yÞ > u4 > > > > > > > > > > u5 > > > : > ; r2 ¼ N d ðxÞ q ¼ N d ðxÞ T u dðuðxÞÞ ¼ @ ðN u ðxÞÞ u
ð35Þ
where
"
Nu ¼
1 xL
0
0
2 xL3 3 xL2 þ 1
x L
0
3
2
x3 L2
0
2
2
#
0 3
2 xL þ x 0 3 xL2 2 xL3
x3 L2
2
xL
ð36Þ
Fig. 7. Displacements and deformations of beam element i.
PHR ðD; uÞ ¼
Z
L
n
0
o dðuÞT D vðDÞ dx Pext ðuÞ
ð31Þ
Complementary energy functional can now be expressed as:
dðDÞ ¼
@ vðDÞ @D
ð32Þ
The functional of external loading for nodal loads F n and distributed
T loads bðxÞ ¼ bx ðxÞ; by ðxÞ obtains the form:
Fn þ Pext ¼ u
Z
L
uT b dx
ð33Þ
0
In order to calculate the state of the element, stationarity of the Hellinger–Reissner functional is imposed by setting its first variation with respect to the two independent fields equal to zero. Then for all variations duðxÞ and dDðxÞ the state defined by ðuðxÞ; DðxÞÞ satisfies the classical equilibrium and strain–displacement compatibility equations.
dPHR ¼
Z
L
Z ddðuÞT D dx þ
0
0
Z
L
T Fn dDT ðdðuÞ dðDÞÞ dx du
L
duT b dx ¼ 0
ð34Þ
0
Then, in the above relation the discretized independent fields are inserted as linear combinations of shape functions and corresponding
For the stress field expressed in stress resultant terms, cross sectional internal forces are calculated from nodal internal forces based on equilibrium considerations. More specifically, equilibrium Eq. (29) are directly integrated and by applying as essential boundary conditions the internal nodal forces, the following interpolation relation is derived.
"
DðxÞ ¼ N D ðxÞ Q þ Dp ;
N D ðxÞ ¼
1 0
0 0 x x 1 L L
#
ð37Þ
where Dp is the particular solution of the direct integration and can be completely determined by the applied loads. Also, the interpolation of the variation of the stress field, since Dp is a known function, has the following form:
dDðxÞ ¼ N D ðxÞ dQ
ð38Þ
Substituting Eq. (35), (37) and (38) into the functional (34) the following discrete form is derived:
Z T TT dPHR ¼ du
L
0 L
Z T
þ dQ 0
N Td N D Q þ Dp dx
N TD
Z 0
L
N Tu b dx F n
dðDÞÞ dx ¼ 0 ðN d T u
ð39Þ
Since Eq. (39) is valid for any arbitrary variation of the two independent fields, both expressions inside the brackets are set equal to zero. Moreover, the shape functions used are orthogonal i.e. RL T N d N D dx ¼ I as geometric nonlinear effects are ignored. Thus 0 the following Euler–Lagrange equations are deduced that describe the state of the element under external loading incorporating the nonlinear material behavior.
I.A. Gkimousis, V.K. Koumousis / Engineering Structures 106 (2016) 399–409
T T Q þ Q p P ext ¼ 0 T u
Z
ð40Þ
L
N TD dðDÞ dx ¼ 0
ð41Þ
where the vector of equivalent nodal external actions has the form:
Z P ext ¼ F n þ
L 0
N Tu b dx
ð42Þ
and Q p is the element vector resulting from the additional cross sectional internal forces due to distributed loading:
Qp ¼
L 0
N Td Dp dx
ð43Þ
Eq. (40) stands for the equilibrium equation and (41) for the compatibility equation. Interelement continuity is imposed by mapping the structure’s global nodal displacements to the element local dis . The only source of nonlinearity in the above equaplacements u tions concerns the derivation of the cross sectional deformations dðDÞ as functions of stress resultants. In the following section the classical linearization technique of the nonlinear constitutive equation is presented, accompanied with a standard Newton–Raphson solution procedure for the solution of the nonlinear system of equations. 4.1. Solution via linearization Various techniques have been proposed for the linearization of the equations that describe the nonlinear problem. Nukala and White [7] have proposed four different solution procedures investigating all different combinations of linearization. Also, Alemdar and White [10] have linearized the constitutive law along with equilibrium and compatibility equations in terms of the independent variables for both geometric and material nonlinearity. The proposed scheme extends the one presented by Hjelmstad and Taciroglu [5] by implementing the rate Eq. (17) describing the constitutive behavior. This is further linearized with respect to the current configuration i following the forward Euler method thus leading to the current tangent fiber modulus expressed as:
h
i
r ¼ ri þ 1 ð1 aÞ Hi1 H i2 E Dei ¼ ri þ Eit Dei
ð44Þ
The current cross sectional flexibility matrix is determined as the inverse of the stiffness matrix using relation (26). Since cross sectional constitutive equation is unique the interpolated forces DðxÞ ^ should be equal to the stress resultants DðxÞ derived from material constitutive law. This equation is linearized as follows:
^ iþ1 ) Ddi ¼ f i Diþ1 D ^i Diþ1 ¼ D
ð45Þ
Consequently, updated cross sectional deformations can be calculated from the following relation:
d
iþ1
L
0
0
Z
Z qi ¼
i ^i ¼ d þ f i Diþ1 D
ð46Þ
Z
L
ci ¼ 0
405
i NTD d Di dx
ð50Þ
^ i þ DDi dx NTD f i Di D p
ð51Þ
i and iþ1 ¼ u i þ Du Considering the expressions Q iþ1 ¼ Q i þ DQ i , u eliminating the increment of nodal forces DQ i , the increment of i can be determined, provided that the elenodal displacements Du ment flexibility matrix F i is invertible. i K i Dui ¼ P iþ1 ext Q int
ð52Þ
h i 1 Q iint ¼ K i ui þ T T Q i ðF i Þ qi þ c i þ Q iþ1 p
ð53Þ
1
where K i ¼ T T ðF i Þ T is the element’s tangent stiffness matrix. Eq. (52) has the same form with the standard incremental equations of the displacement based method the difference being that internal element forces and stiffness matrix are exact as they satisfy equilibrium equations in strong form. On the other hand compatibility is satisfied in weak form as nodal deformations are calculated as the weighted average of the cross sectional deformations. Considering the computational aspects of this formulation, cross sectional residuals are directed to the element’s residuals and then to structure’s residuals. Hence this algorithm is classified as an L–L algorithm according to Nukala and White [7] characterization. Following a Newton–Raphson procedure for the global equilibrium equations, the procedure of calculating updated global stiffness matrix and internal nodal forces proceeds as follows: (a) After solving for an incremental step for the whole structure’s equilibrium equation, the vector of incremental nodal displacements is evaluated and nodal total displacements are updated. Mapping the global structure’s displacements iþ1 of each element to every element, the displacements u are defined and the element state determination begins. (b) Element incremental nodal forces DQ i are calculated from Eq. (48) and total nodal forces are updated Q iþ1 ¼ Q i þ DQ i . (c) Cross-sectional forces are calculated from interpolation: Diþ1 ¼ Di þ N D DQ i þ DDp , where DDp is the increment of cross sectional stress resultants due to distributed loads. (d) Cross-sectional deformations are evaluated from the linearized constitutive Eq. (46) and updated fiber strains eiþ1 are derived from Eq. (22). (e) The linearized material Eq. (44) is solved and updated tangent fiber modulus Eit and fiber stresses riþ1 are evaluated. Then integration over the cross-section evaluates the updated tangent stiffness kðxÞiþ1 which gives the cross sec-
Substituting Eq. (46) in (41) results in the linearized incremental equations:
T T Q iþ1 þ Q piþ1 P iþ1 ext ¼ 0
tional flexibility f ðxÞiþ1 by its inverse. (f) Integration over the length using Eqs. (49) and (50) provides
ð47Þ
element flexibility matrix F iþ1 and element’s nodal deformations qiþ1 .
iþ1 F i DQ i qi c i ¼ 0 T u
ð48Þ
(g) Eq. (53) determines the nodal internal forces Q iþ1 for the next step, while by assembling every element’s contribution
where the element’s flexibility, nodal deformations and element’s deformations residual are derived integrating over the beam length as:
Z
L
Fi ¼ 0
N TD f i N D dx
ð49Þ
the global structure’s stiffness matrix K iþ1 and nodal internal S forces P iþ1 int are derived. (h) Every element’s residual is directed to the structural level
iþ1
P iþ1 and tolerance is checked. P ext int 6 tol.
406
I.A. Gkimousis, V.K. Koumousis / Engineering Structures 106 (2016) 399–409
5. Examples 5.1. Cantilever column In this example the proposed element formulation is validated solving the problem presented in Ref. [6]. In this reference, different numerical procedures and element discretization schemes of a mixed three-field Hu–Washizu element, first proposed in [9] are compared. Also, the proposed steel constitutive law (17) is implemented following the classical displacement based formulation and the results are presented for comparison reasons. The problem consists of a non-dimensional cantilever column with height of 120 units and rectangular cross section of 15 20 units. Material is considered bilinear with elastic modulus of elasticity Eel ¼ 29; 000, yield stress f y ¼ 50 and hardening post-yield to pre-yield ratio a ¼ 0:01. The cross-section is discretized in 5 layers, each one participating in the overall cross-section response according to mid-point integration. Element behavior is integrated along its length using Gauss–Lobatto integration with 5 control points. Implementing the classical displacement based formulation requires several elements in the discretization to achieve a satisfactory and objective response [44]. For this reason two different discretization schemes were tested consisting of one (1DB) and four elements (4DB) respectively. The results are presented in Fig. 8 relating base shear and top displacement, with the proposed hysteretic formulation proven very efficient in generating the accurate response. Also, the proposed two-field mixed element (Proposed-1HR) is compared for the same example against the three-field force based element of Ref. [6] (reference-1FB) and the results are quite satisfactory as presented in Fig. 9.
Fig. 9. Base shear vs top displacement of a cantilever column, comparison with mixed formulations.
Fig. 10. Continuous beam under vertical dynamic excitation.
5.2. Continuous beam under cyclic vertical load In this example the proposed formulation is employed for the solution of a 2-span continuous beam (Fig. 10) under a pseudostatic cyclic external loading consisting of a varying concentrated vertical load at the middle of the first span following a piecewise linear pattern presented in Fig. 11 and Table 2. Also, cross section is an IPE 300 European steel cross section which is discretized in 50 layers. The comparison is twofold, first the displacement based approach is compared with the proposed mixed formulation, and then the proposed hysteretic steel model is compared with the Menegotto–Pinto model. Material and model properties are given in Table 3. In Fig. 12 the two element formulations are compared. For the displacement based formulation two different discretization
Fig. 11. Vertical load input.
Table 2 External loading values. Time (s)
External force (KN)
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5
0 400 300 200 150 250 50 300 0 350 50 380
Table 3 Steel properties and model parameters. Steel properties
Fig. 8. Base shear vs top displacement of a cantilever column, comparison with displacement based formulation.
Model parameters
Young modulus E Yield stress fy
210 GPa 235 MPa
Kinematic hardening ratio a
0.015
Proposed
Menegotto–Pinto
n1 = 10 n2 = 1
R0 = 21 cR1 = 17.5 cR2 = 0.15
I.A. Gkimousis, V.K. Koumousis / Engineering Structures 106 (2016) 399–409
407
Fig. 12. Comparison of vertical midspan displacement vs external loading history. Fig. 15. Northridge_TAR090 ground motion excitation.
a single element captures directly the real cyclic response with no need for further refinement. In Fig. 13 the proposed steel hysteretic model is compare with the Menegotto–Pinto model. The advantage of the proposed model to avoid overshooting after short reversals is reflected in this case also in the overall behavior of the structure. More specifically, after a series of plastic short reversals, which occur often in real seismic excitations, Menegotto–Pinto model underestimates structural displacements. Hence, the proposed model produces results, which are more accurate and on the safe side from the engineering practice perspective. 5.3. Steel frame under seismic excitation Fig. 13. Load–displacement history, comparison of the proposed model with Menegotto–Pinto model.
schemes with 3 control sections are utilized, while for the mixed formulation a single element discretization with 15 control sections is sufficient to yield accurate results. A dense discretization consisting of 40 equal beam elements with 3 Gauss points per element for the total length of the continuous beam is used to overcome the deficiencies of the classical displacement based stiffness method. Moreover, steel behavior is modeled according to the proposed hysteretic model. It is evident that the discretization with 3 elements following the stiffness method (DB 3 elements) delivers inaccurate results and is not capable of capturing the real behavior. This inefficiency is magnified for intensive external loading, as overestimation of the maximum strength results into more than 2.5 times decreased plastic deformations compared to a more refined discretization. On the other hand the proposed two-field mixed element (HR 3 elements) satisfies the equilibrium equations in strong form and
In this example a 2 story–3 bay steel plane frame (Fig. 14) is excited with the Northridge_TAR090 seismic strong motion record (Fig. 15). The frame is designed according to Eurocode 3 specifications and consists of HEA-450 columns and IPE-220 beams discretized in 50 layers, while steel material class is S355 with 1.5% kinematic hardening. Also, distributed load of G = 15.0 KN/m and the associated mass is considered on every beam including selfweight. Lumped Mass and no viscous damping are considered, while direct time integration was performed according to Newmark numerical scheme with beta = 0.25, gamma = 0.5 and time step dt = 0.01. This example is used to compare the results of the proposed hysteretic model with the Menegotto–Pinto steel model for the response of the structure in terms of accuracy and computing time. The proposed mixed element is used in both simulations and the model parameters for the two steel models are those presented in Table 3. First, the overall response of both models is compared by monitoring the top story lateral displacement for the first 10 s of the
Fig. 14. 2 story–3 bay steel plane frame.
408
I.A. Gkimousis, V.K. Koumousis / Engineering Structures 106 (2016) 399–409
20% as compared to Menegotto–Pinto model. The comparison is performed using a general program for frame structural analysis developed in Matlab which accounts for the two different steel models. The time step of the analysis is dt = 0.01 s, while similar deviations were obtained also for different time steps, indicating that both models demonstrate the same tendency for convergence.
6. Conclusions
Fig. 16. Top lateral displacement history for the proposed Hysteretic and Menegotto–Pinto models.
Fig. 17. Moment–curvature diagrams of beam section A, using the proposed Hysteretic and Menegotto–Pinto models.
In this work a new hysteretic uniaxial steel model is proposed describing steel cyclic behavior, which is implemented within a fiber beam–column element based on Hellinger–Reissner mixed variational principle. The hysteretic model is developed on the basis of classical plasticity considerations and has the advantage of embodying all nonlinear phases in a single nonlinear differential equation. The model is further modified to account for steel kinematic hardening, Bauschinger effect and compliance with plasticity postulates after short reversals. The model avoids overshooting after short reversals, being on the same time computationally effective and in that respect more advantageous compared to the well-established Menegotto–Pinto model. Incorporating the constitutive equations into the general formulation of mixed beam finite elements, element equilibrium and compatibility equations are derived. Linearization of the differential constitutive equation is employed, which leads to the incremental equilibrium and compatibility equations, which are solved using a Newton’s type iterative numerical scheme. The main feature and computational advantage of the proposed scheme relies on the explicit and direct calculation of stresses, with no need for predictor–corrector schemes by solving the hysteretic constitutive equations in every time step. Verification of the proposed model was performed on the basis of displacement based or mixed methods using examples from the literature demonstrating the computational efficiency of the proposed formulation.
References
Fig. 18. Total computational time for hysteretic and Menegotto–Pinto models.
excitation as presented in Fig. 16. Minor differences between the two models are observed in this case. This happens besides the fact that certain elements demonstrate a more pronounced difference as observed in Fig. 17, which at the structural level seem to compensate. Also, as for the vertical displacements in example 2, curvatures in this case are underestimated by the Menegotto–Pinto model. After several cycles of excitation, the critical section is oscillating elastically at larger negative curvatures according to the proposed model compared with the Menegotto–Pinto model. The main advantage of the proposed formulation, which avoids overshooting, relies on its computational efficiency. The total computing time for four different schemes with 3, 5, 7, 9 control sections per element is recorded and the results are plotted in Fig. 18. The reduction of computational time is of the order of
[1] Bathe KJ. Finite element procedures. New York: Prentice Hall Engineering, Science, Mathematics; 2007. [2] Spacone E, Filippou FC, Taucer FF. Fibre beam column model for nonlinear analysis of R/C frames I: formulation. Earthq Eng Struct Dyn 1996;25 (7):711–25. [3] Neuenhofer A, Filippou FC. Evaluation of nonlinear frame finite-element models. J Struct Eng 1997;123(7):958–66. [4] Hjelmstad KD, Taciroglu E. Mixed methods and flexibility approaches for nonlinear frame analysis. J Constr Steel Res 2002;58:967–93. [5] Hjelmstad KD, Taciroglu E. Variational basis of nonlinear flexibility methods for structural analysis of frames. J Eng Mech 2005;131(11):1157–69. [6] Saritas A, Soydas O. Variational base and solution strategies for the non-linear force-based beam finite elements. Int J Non-linear Mech 2012;47:54–64. [7] Nukala PK, White DW. Variationally consistent state determination algorithms for nonlinear mixed beam finite elements. Comput Methods Appl Mech Eng 2004;193:3647–66. [8] Hjelmstad KD, Taciroglu E. Mixed variational methods for finite element analysis of geometrically non-linear, inelastic Bernoulli–Euler beams. Commun Numer Methods Eng 2003;19(10):809–32. [9] Taylor RL, Filippou FC, Saritas A, Auricchio F. A mixed finite element method for beam and frame problems. Comput Mech 2003;31:192–203. [10] Alemdar BN, White DW. Displacement, flexibility and mixed beam–column finite element formulations for distributed plasticity analysis. J Struct Eng 2005;131:12. 1811. [11] Alsafadie R, Hjiaj M, Somja H, Battini J-M. A comparative study of displacement and mixed-based corotational finite element formulations for elasto-plastic, three-dimensional beam analysis. Eng Comput: Int J Comput-Aided Eng Softw 2011;28(7):939–82. [12] Correia AA, Almeida JP, Pinho R. Force-based higher-order beam element with flexural-shear-torsional interaction in 3D frames. Part I: theory. Eng Struct 2015;89:204–17. [13] Simo JC, Taylor RL. Consistent tangent operators for rate – independent elastoplasticity. Comput Methods Appl Mech Eng 1985;48:101–18. [14] Simeonov VK, Sivaselvan MV, Reinhorn AM. Nonlinear analysis of structural frame systems by the state-space approach. Comput Aided Civ Infrastruct Eng 2000;15:76–89.
I.A. Gkimousis, V.K. Koumousis / Engineering Structures 106 (2016) 399–409 [15] Jafari V, Rahimian M, Vahdan SH. Improved formulation in mixed-based statespace approach for large displacement inelastic analysis of frames. J Eng Mech 2011;137(3):196–204. [16] Triantafyllou S, Koumousis V. Small and large displacement dynamic analysis of frame structures based on hysteretic beam element. J Eng Mech 2011;138 (1):36–49. [17] Triantafyllou SP, Koumousis VK. Hysteretic finite elements for the nonlinear static and dynamic analysis of structures. J Eng Mech 2013;140(6). [18] Menegotto M, Pinto PE. Method of analysis for cyclically loaded R.C. plane frames including changes in geometry and non-elastic behaviour of elements under combined normal force and bending. In: Symposium on the resistance and ultimate deformability of structures acted on by well defined repeated loads. International association for bridge and structural engineering (IABSE) Lisbon, Portugal; 1973. p. 15–22. [19] Ramberg R, Osgood WR. Description of stress–strain curve by three parameters. Technical note 902. NACA; 1943. [20] Monti G, Nuti C. Nonlinear behaviour of reinforcing bars including buckling. J Struct Eng 1992;118(12):3268–84. [21] Kashani MM, Crewe AJ, Alexander NA. Nonlinear stress–strain behavior of corrosion-damaged reinforcing bars including inelastic buckling. Eng Struct 2013;48:417–29. [22] Kunnath S, Heo Y, Mohle J. Nonlinear uniaxial material model for reinforcing steel bars. J Struct Eng 2009;135(4):335–43. [23] Zong Z, Kunnath S, Monti G. Material model incorporating buckling of reinforcing bars in RC columns. J Struct Eng 2014;140(1). [24] Kashani MM, Crewe AJ, Alexander NA. Nonlinear stress–strain behaviour of corrosion-damaged reinforcing bars including inelastic buckling. Eng Struct 2013;48:417–29. [25] Shi G, Wang M, Bai Y, Wang F, Shi YJ, Wang YQ. Experimental and modeling study of high-strength structural steel under cyclic loading. Eng Struct 2012;37:1–13. [26] Wang YB, Li GQ, Cui W, Chen SW. Experimental investigation and modeling of cyclic behavior of high strength steel. J Constr Steel Res 2015;104:37–48. [27] Filippou FC, Popov EP, Bertero VV. Effects of bond deterioration on hysteretic behavior of reinforced concrete joints. Report. UCB/EERC-83/19. Berkeley: Earthquake Engineering Research Centre, University of California.
409
[28] Dodd LL, Restrepo-Posada JI. Model for predicting cyclic behaviour of reinforcing steel. J Struct Eng ASCE 1995;121(3):433–45. [29] Balan TA, Filippou FC, Popov EP. Hysteretic model of ordinary and highstrength reinforcing steel. J Struct Eng 1998;124(3):288–97. [30] Reissner E. On a variational theorem in elasticity. J Math Phys 1950;29:90–5. [31] Valanis KC, lee CF. Endochronic theory of cyclic plasticity with applications. J Appl Mech 1984;51(2):367–74. [32] Sivaselvan MV, Reinhorn AM. Nonlinear Structural analysis towards collapse simulation – a dynamical systems approach. Technical report, Buffalo, New York; 2003. [33] Wen YK. Method of random vibration of hysteretic systems. J Eng Mech Div 1976;102:249–63. [34] Lubliner J. Plasticity theory. New York: Dover Publications, Inc.; 2008. [35] Borja RJ. Plasticity, modeling & computation. Berlin, Heidelberg: SpringerVerlag; 2013. [36] Ucak A, Tsopelas P. Accurate modeling of the cyclic response of structural components constructed of steel with yield plateau. Eng Struct 2012;35:272–80. [37] Ayhan B, Jehel P, Brancherie D, Ibrahimbegovic A. Coupled damage-plasticity model for cyclic loading: theoretical formulation and numerical implementation. Eng Struct 2013;50:30–42. [38] Mendes LAM, Castro LMSS. A simplified reinforcing steel model suitable for cyclic loading including ultra-low-cycle fatigue effects. Eng Struct 2014;68:155–64. [39] Thyagarajan RS. Modeling and analysis of hysteretic structural behavior. Technical report. Caltech EERL, EERL-89-03; 1989. [40] Wang CH, Foliente GC. Discussion on hysteretic models for deteriorating inelastic structures. J Eng Mech ASCE 2001;127(11):1200–2. [41] Charalampakis AE, Koumousis VK. A Bouc–Wen model compatible with plasticity postulates. J Sound Vib 2009;322(4-5):954–68. [42] Kottari AK, Charalampakis AE, Koumousis VK. A consistent degrading Bouc– Wen model. Eng Struct 2014;60:235–40. [43] Thompson KJ, Park R. Stress–strain model for grade 275 reinforcing steel with cyclic loading. Bull NZ Natil Soc Earthq Eng 1978;11(2). [44] Calabrese A, Almeida JP, Pinho R. Numerical issues in distributed inelasticity modeling of RC frame elements for seismic analysis. J Earthq Eng 2010;14:38–68.
Contents lists available at ScienceDirect
Engineering Structures journal homepage: www.elsevier.com/locate/engstruct
Inelastic mixed fiber beam finite element for steel cyclic behavior I.A. Gkimousis, V.K. Koumousis ⇑ Institute of Structural Analysis & Aseismic Research, National Technical University of Athens, Zografou Campus, 15780 Athens, Greece
a r t i c l e
i n f o
Article history: Received 25 April 2015 Revised 26 August 2015 Accepted 12 October 2015 Available online 11 November 2015 Keywords: Steel cyclic behavior Hysteretic model Hellinger–Reissner variational principle Mixed beam formulation Nonlinear frame analysis
a b s t r a c t In this work a new hysteretic uniaxial steel model is determined to describe steel cyclic behavior, which is further implemented to derive a fiber beam–column element on the basis of Hellinger–Reissner principle. The proposed model maintains full memory of the loading path and evolves following a single nonlinear differential equation expressing the entire hysteresis. The element is capable of addressing the overshooting problem of the existing models which occurs during short reversals. The state determination of the proposed element is investigated numerically following the linearization of the derived equations. Numerical results are presented that validate the proposed approach and demonstrate its computational efficiency. Ó 2015 Elsevier Ltd. All rights reserved.
1. Introduction The state of a deformable body subjected to body forces, tractions and kinematic boundary conditions is considered fully defined when the displacements, stresses and deformations are determined at any point of the body. In particular, for earthquake engineering and structural analysis of skeletal structures, beam elements usually based on the Euler–Bernoulli theory assumption that plane sections remain plane and perpendicular to the deformed axis, are typically considered. This facilitates considerably state determination of such elements. Displacement based beam elements were initially used following the classical stiffness method in which displacements were the only considered independent field [1]. When cubic and linear shape functions are employed for the transverse and axial displacements respectively, the resulting displacement field leads to constant axial deformation and linear curvature, which however is not appropriate when plastic deformations occur. To address this deficiency a structural member should be discretized in more than one element at the expense of increasing computational cost. Also, equilibrium equations are only fully satisfied at element nodes, while within the element they are satisfied in weak form as they are not valid for all possible displacement fields that satisfy essential boundary conditions. To resolve this problem, force based models were proposed that interpolate nodal forces within the element maintaining ⇑ Corresponding author. Tel.: +30 2107721657; fax: +30 2107721651. E-mail addresses: [email protected] (I.A. Gkimousis), [email protected]. gr (V.K. Koumousis). http://dx.doi.org/10.1016/j.engstruct.2015.10.023 0141-0296/Ó 2015 Elsevier Ltd. All rights reserved.
equilibrium. These models were implemented in the framework of the stiffness method of structural analysis and in that respect they are considered ‘‘mixed” as they use both force and displacement fields as independent ones. One of the first consistent and general force based beam model was proposed by Spacone et al. [2] and was later simplified numerically by Neuenhofer and Filippou [3]. Although the force based method proved very efficient and is currently widely used, there were some concerns about its variational consistency that were resolved by Hjelmstad and Taciroglu [4]. Moreover, the same authors in [5] showed that it is possible to provide non-variationally consistent force-based elements within the ‘‘nonlinear flexibility” framework by enforcing equilibrium directly. Also, various local and global solution strategies originate from the variational structure of the mixed beam elements as described in [6,7]. Thereafter, mixed methods seem to dominate the research field of nonlinear beam problems and corresponding numerical procedures as they proved more efficient following also the work of Hjelmstad and Taciroglu [8], Taylor et al. [9], Alemdar and White [10], Alsafadie et al. [11] and Correia et al. [12]. In this context elastoplastic material models were represented in classical form relying on the notions of yield surface, flow rule and hardening parameters, while their incorporation in the state determination process in linearized form gave rise to the returnmapping algorithm [13]. However, cyclic behavior can also be modeled using hysteretic evolution differential equations. Following this approach, Simeonov et al. [14] developed a force based element where material constitutive relations are considered in rate form and are solved simultaneously with the global differential equations of motion in state-space form. Also, Jafari et al. [15] extended
400
I.A. Gkimousis, V.K. Koumousis / Engineering Structures 106 (2016) 399–409
this formulation in large displacement analysis following a displacement based formulation. In addition Triantafyllou and Koumousis [16,17] proposed a finite element procedure where material nonlinearity is treated constitutively at the element level through proper implementation of the Bouc–Wen hysteretic rule. In the case of elastoplastic steel cyclic behavior various models have been proposed among which the widely used Menegotto– Pinto model [18] originating from the generalization of Ramberg–Osgood model [19]. Moreover, for efficient modeling of the inelastic buckling of reinforcing steel bars under cyclic behavior, the Monti–Nuti model [20] is commonly used as an enhancement of the Menegotto–Pinto model, which accounts for four different hardening rules. Furthermore, inelastic buckling of reinforcing bars has been investigated in several other studies, i.e. [21–24], while also modeling of high strength structural steel is addressed in Refs. [25,26]. Although the core of several steel models is the Menegotto– Pinto model, it is characterized by an overshooting of the reloading branch after short reversals [27]. This feature was considered of minor importance and originates from an effort to reduce computational cost by truncating model’s memory. Although attempts were made to tackle this aspect by distinguishing the reversals as major and minor [28] or complete and incomplete [29] the problem can be also addressed from a different perspective. In this work a new small displacement fiber beam–column element is proposed, which incorporates a uniaxial steel model for cyclic loading that resolves the overshooting behavior after short reversals, being on the same time computationally efficient. To accomplish this, a hysteretic model is developed that maintains full memory of the loading path and evolves following a single differential equation expressing the entire hysteresis. This constitutive behavior is implemented in the general framework of the two-field Hellinger–Reissner formulation [30] resulting in a wellestablished state determination algorithm. This is further investigated numerically following the linearization of the equilibrium and compatibility element equations. The rest of the paper is organized as follows: In the first section the basic structure of the proposed model is derived from classical plasticity considerations incorporating kinematic hardening. Then, the modifications for modeling steel cyclic behavior are implemented and the comparison with the Menegotto–Pinto model is performed underlining the features of the proposed formulation. Cross-sectional constitutive equations are derived in the sequel from fiber integration and the two-field Hellinger–Reissner principle is used to derive the element equilibrium and compatibility equations. The standard linearization method is implemented for solving the element equations and the respective state determination process is discussed. Finally, three numerical examples are presented that verify the proposed beam element and demonstrate the performance of the hysteretic constitutive relations embedded into a two-field variational formulation for the inelastic analysis of steel frames. 2. Inelastic Euler–Bernoulli beam theory 2.1. Fiber plasticity A hysteretic model incorporates the entire inelastic loading path of a deformable body, namely elastic loading, yielding, hardening and unloading in a single nonlinear differential equation that embodies both the yield surface and hardening rule. Mathematically this addresses the entire evolution process without the need of incremental considerations [31]. In this context Sivaselvan and Reinhorn [32] based on Bouc–Wen model [33] proposed a hysteretic model in stress resultant terms derived explicitly from classical plasticity theory.
In this work the stress–strain constitutive law of an elastoplastic fiber-rod subjected to uniaxial tension is determined according to the strain decomposition rule in rate form as [34]:
r_ ¼ Eðe_ e_ p Þ
ð1Þ
where r_ is the rate of normal stress, e_ is the rate of total strain and e_ p is the rate of plastic strain. This relation indicates that stresses evolve proportionally to the evolution of the elastic strains. For the same fiber a yield function with linear kinematic hardening and the same initial yield stress in tension and compression rþy0 ¼ ry0 ¼ ry0 is expressed in the following form:
U r; b; ryo ¼ jr bj ry0 6 0 with b being the back stress. Plastic strain rate according to the flow rule:
e_ p ¼ k_
@U _ ¼ k sgnðr bÞ @r
ð2Þ _p
e is determined ð3Þ
where k_ P 0 is the plastic multiplier which is actually the magnitude of the strain rate, with the signum function defining its sense. Axial stress and plastic multiplier are restricted by unilateral constrains representing restrictions that signify whether the material has yielded or not, resulting from Karush–Kuhn–Tucker (KKT) optimality conditions. These are expressed in the following form:
k_ Uðr; b; ry0 Þ ¼ 0
ð4Þ
During plastic response ðk_ > 0; Uðr; b; ry0 Þ ¼ 0Þ the consistency condition is derived by differentiating Eq. (4):
_ ðr; b; ry0 Þ ¼ 0 ) U _ ðr; b; ry0 Þ ¼ 0 ) r_ ¼ b_ k_ U
ð5Þ
Also, when the fiber deforms elastically ðk_ ¼ 0; Uðr; b; ry0 Þ < 0Þ the sign of the rate of the yield function indicates loading or unloading _ < 0 respectively. _ > 0 and U with U Taking into account Eq. (5) meaning that the rate of the back stress during plastic loading is equal to the rate of the axial stress, the following relation is considered:
b_ ¼ H e_ p ¼ H k_ sgnðr bÞ
ð6Þ
where H is the hardening ratio, i.e. the slope of the stress–strain curve in plastic strain terms (Fig. 1). Substituting Eq. (1), (3) and (6) in Eq. (5) the following relation is obtained:
k_ ¼ sgnðr bÞ
E e_ EþH
ð7Þ
Furthermore, introducing a ¼ Et =E as the ratio of the post-yield tangent modulus over the elastic modulus, the relation between the hardening ratio H and tangent modulus Et is established as:
H¼
Et a E ¼ 1 Et =E 1 a
ð8Þ
Then Eq. (3), using relations (7) and (8) obtains the following simple form:
e_ p ¼ ð1 aÞ e_
ð9Þ
Eq. (9) holds only for the plastic deformation phase. During elastic loading or unloading plastic strain is not induced as k_ ¼ 0. This is controlled by the compact form of KKT conditions (4). To include all loading phases in a single relation, Eq. (9) is extended by considering two Heaviside type functions acting as switches as follows:
e_ p ¼ H1 H2 ð1 aÞ e_
ð10Þ
where H1 controls yielding and H2 controls loading/unloading states. Consequently the following distinct phases can be described in a single equation:
401
I.A. Gkimousis, V.K. Koumousis / Engineering Structures 106 (2016) 399–409
Fig. 1. Uniaxial bilinear constitutive law in terms of (a) total strain and (b) plastic strain.
H1 ¼ H2 ¼
0; elastic deformation ðe_ p ¼ 0Þ 1; plastic flow ðe_ p ¼ H2 ð1 aÞ e_ Þ
ð11Þ
0; unloading ðe_ p ¼ 0Þ 1; loading ðe_ p ¼ H1 ð1 aÞ e_ Þ
ð12Þ
Function H1 emerges from the yield function (2) by properly expressing the back stress b as a function of the calculated plastic strain according to the following relation [35]:
b¼He ¼ p
a
1a
E
ee
el
¼
a
1a
ðE e rÞ
ð13Þ
Taking into account Eq. (13), the yield function can be expressed as follows:
raEe 61 ð1 aÞ r
U ¼
ð14Þ
y0
which, in the limit case of yielding
r ¼ ry0 ; e ¼ ey0 sets U ¼ 1. If
the above relation is raised to the nth power the absolute value approximates the step function with discrete values [0, 1] as n ! 1. Consequently the following relation holds:
r a E e n H1 ð1 aÞ ry0
ð15Þ
Parameter n controls the form of transition from the elastic to plastic branch of the stress–strain law, which is smoother for lower values, approaching a bilinear behavior for higher ones ðn P 8Þ. Also the expression ðrh ¼ r a E eÞ refers to the hysteretic part of the normal stress. Moreover, function H2 aims at controlling loading and unloading following the sign of the yield function rate with the positive sign indicating loading and negative sign unloading. Furthermore, smoothening of the Heaviside function to yield zero value for unloading and unity for loading is handled as follows: _ 1 þ sgn U
_ 1 þ sgnðr bÞ sgnðr_ bÞ 1 þ sgnðr bÞ sgnðe_ Þ ¼ ) ¼ 2 2 2 H2 ¼ 0:5 þ 0:5 sgn r aðE e ry0 Þ e_ 0:5 þ 0:5 sgnððr a E eÞ e_ Þ H2 ¼
ð16Þ
In the above relation, b_ ¼ 0 for elastic response and the sign of the stress rate r_ is replaced by the sign of the strain rate e_ . After defining functions H1 and H2 , Eq. (10) is substituted in Eq. (1) and the final form of the constitutive relation is obtained:
r_ ¼ ½1 ð1 aÞ H1 H2 E e_ ¼ Et e_
ð17Þ
where Et is the tangent modulus. 2.2. Steel cyclic behavior Steel stress–strain relationship under cyclic loading primarily aims at expressing the linear elastic and hardening branches in the first excursion, together with the Bauschinger effect for subsequent cycles. Yield plateau [36], plasticity-damage coupling [37] and low cycle fatigue [38] are not considered in this work. To account for steel cyclic behavior Eq. (17) along with relations (15) and (16) for kinematic hardening are implemented. The initial bilinear constitutive behavior is obtained using a high value for the exponent parameter n ðn1 P 8Þ. The bilinear behavior is incorporated just for the initial yield point, while for subsequent yielding in the opposite direction gradual transition in the plastic regime is manifested. In the proposed hysteretic model this is accomplished by subsequently reducing parameter n to lower values n2 ðn2 < 2Þ. Hysteretic models have been criticized in the past for their incompatibility with Drucker’s or Ilyushin’s postulates of plasticity for partial unloading–reloading in the nonlinear path [39]. Based on previous work [40] an accurate solution was recently proposed that eliminates the problem by properly modifying the reloading path of the hysteretic spring for the Bouc–Wen model [41] or the Sivaselvan–Reinhorn model [42]. In this work the same approach is utilized, which leads to the modification of the switch function H2 in the following form:
2 ¼ H2 ð1 RÞ H
ð18Þ
where Rðrh ; e; ry0 Þ is a stiffening factor that controls stiffness recovery in the reloading branch, which is expressed as:
R ¼ 0:5 ½1 þ sgnðrh rr Þ sgnðrh Þ
er ec er e
ð19Þ
402
I.A. Gkimousis, V.K. Koumousis / Engineering Structures 106 (2016) 399–409
where rh is the hysteretic part of the normal stress ðrh ¼ r a E eÞ; rr and er are the stress and strain respectively corresponding to the last observed reversal point and ec is the modified strain in the unloading path that corresponds to the given hysteretic stress rh (Fig. 2):
ec ¼
rh rr ey0 þ er sgnðrh Þ ry0
ð20Þ
where also appropriate tracking of the reversal point is needed to accurately handle different types of short reversals as indicated in [41,42]. 2.2.1. Comparison with experimental data Experiments on steel rebars under cyclic response can be found in the work of Thompson and Park [43]. One of these experiments is presented in the following Fig. 3 along with the predicted response of the proposed model. Steel grade is S275, while post to pre-yield ratio is considered as a ¼ 0:01. For model parameters, initial and final exponents were selected as n1 ¼ 8 and n2 ¼ 0:9 respectively to achieve satisfactory agreement with observed experimental behavior. 2.2.2. Comparison with Menegotto–Pinto model To illustrate the performance of the proposed model, which incorporates the modifications made for modeling steel cyclic behavior, a direct comparison with the well-established Menegotto–Pinto model [18] is performed. For both models an external strain history presented in Table 1 is applied and model parameters are selected in such way as to induce similar outcome. More specifically, material modulus is E = 210 GPa, yield stress ry0 ¼ 235 MPa, post-yield to pre-yield ratio a ¼ 0:02, while for the proposed model the initial exponent is n1 ¼ 10 and the final exponent is n2 ¼ 0:2. The parameters of Menegotto–Pinto model are considered as: R0 ¼ 20, cR1 ¼ 18:5, cR2 ¼ 0:1. In Fig. 4 it is observed that the proposed model overcomes the overshooting of the Menegotto–Pinto model after the partial unloading branch. This overshooting originates from the restricted memory of the Menegotto–Pinto model as underlined by Filippou et al. [27]. To address this issue, Dodd and Restrepo–Posada [28] distinguished reversals as major, minor and simple, while Balan et al. [29] separated them into complete and incomplete, with a value of 2f y being the threshold parameter of their distinction. However, these modifications are perplexing the original model increasing bookkeeping of data and thus computational cost. On the other hand the proposed model imposes no restriction on memory returning exactly on the previous loading path as it is evident in Fig. 4.
Fig. 3. Comparison with experiment [43].
Table 1 Strain history values. Time (s)
Strain (–)
0 1 2 2.5 2.7 3.5
0 0.04 0.005 0.009 0.0085 0.04
Fig. 4. Comparison between the proposed and Menegotto–Pinto model.
3. Cross-sectional relations Cross-sectional constitutive relations are derived from integration of the local constitutive equations of all unit length fibers using midpoint integration. Especially for planar analysis and symmetrical cross sections, fibers can be treated in layers (Fig. 5). Following the Euler–Bernoulli beam theory the normal stress rðx; yÞ; the axial strain eðx; yÞ and the axial and transverse displacement ux ðx; yÞ; uy ðx; yÞ are considered. By considering the kinematic assumption that plane sections before deformation remain plane and normal to the elastic line after deformation, the axial displacement ux ðx; yÞ of any fiber can be expressed in terms of the displacement of the beam reference axis x as:
ux ðx; yÞ ¼ uðxÞ y
Fig. 2. Treatment of short reversals in the proposed formulation.
@wðxÞ @x
ð21Þ
where uðxÞ and wðxÞ are the axial and transverse displacements of the cross section at the reference axis. Thus the following relation holds for the axial strain:
I.A. Gkimousis, V.K. Koumousis / Engineering Structures 106 (2016) 399–409
403
the tangent stiffness matrix at every point of the nonlinear path during dynamic loading. Moreover, at the boundaries the three fields are manifested as nodal values. The standard six nodal displacements ¼ fu1 ; u2 ; r1 ; u4 ; u5 ; r2 gT , and rotations for the 2D beam element u are related with deformations q ¼ fd; h1 ; h2 gT i.e. axial deformation d along the reference axis and the chord rotations at both ends. In addition the stress resultants Q ¼ fN; M 1 ; M 2 gT are introduced as presented in Fig. 6. The relation between the nodal displacements and member deformations is established with the transformation matrix T, which removes rigid body motion as follows (Fig. 7):
2
; q¼T u
Element equilibrium under distributed applied
T bðxÞ ¼ bx ðxÞ; by ðxÞ in both directions x, y requires:
Fig. 5. Cross-sectional local axis and layer definition.
T
eðx; yÞ ¼ e0 ðxÞ y uðxÞ ¼ ½1 y dðxÞ ¼ l dðxÞ
ð22Þ T
where the deformation field dðxÞ ¼ fe0 ðxÞ; uðxÞg consists of the axial strain e0 ðxÞ and curvature uðxÞ at the reference axis. Euler–Bernoulli beam theory considers only the effects of axial stress r, thus by applying equilibrium conditions within the cross section, the stress resultants are evaluated as follows:
Z
Z
NðxÞ ¼ A
rðx; yÞdA; MðxÞ ¼
A
y rðx; yÞdA
ð23Þ
The same relations can be casted in matrix form as:
^ DðxÞ ¼
Z
A
½1 yT rðx; yÞdA
ð24Þ
^ where DðxÞ ¼ fN ðxÞ; MðxÞgT . By replacing Eq. (17) along with Eq. (22) in the rate form of Eq. (24) the ‘‘cross sectional constitutive relation” is derived:
_ ^_ DðxÞ ¼ kðxÞ dðxÞ with
Z
kðxÞ ¼
½1 yT Et ½1 ydA
ð25Þ
ð26Þ
A
Finally, the flexibility matrix of the cross section is set as the inverse of the cross-sectional stiffness.
f ðxÞ ¼ kðxÞ1
3 1 0 0 1 0 0 6 7 T ¼ 4 0 1=L 1 0 1=L 0 5 0 1=L 0 0 1=L 1
ð27Þ
4. Two field beam element variational formulation Mixed methods appear to dominate the classical displacement based formulation as they satisfy element equilibrium in strong form avoiding the need for further element discretization. The equilibrium error is of major importance, as internal forces are overestimated leading to inaccurate structural response. Two or three field methods based on Hellinger–Reissner or Hu–Washizu principles respectively are able to diminish this error successfully. In this work the two-field Hellinger–Reissner principle is employed, disregarding the inherent error in the deformation field which is considered herein of secondary importance. Having the inelastic constitutive equation for a beam fiber in a cross section via relations (17) and (18), the next step is to define the overall nonlinear response of the element. This refers to the displacement, deformation and stress fields across the beam length and
@ 2 MðxÞ by ðxÞ ¼ 0 @x2 @NðxÞ bx ðxÞ ¼ 0 @x
ð28Þ
loads
ð29Þ
A mixed formulation is presented in the sequel based on the two-field Hellinger–Reissner potential [30], where both displacement and stress fields are considered as independent. Classical stiffness analysis of frame structures is based on the continuity of the displacement field between elements. This is accomplished explicitly as the cubic interpolation functions satisfy the kinematic boundary conditions at the end nodes. In the pure force formulation based on the principle of virtual forces, only the stress field is independent and interpolation functions based on element equilibrium relate internal and nodal forces. For a single element analysis, force based method is straightforward as interpolation base functions that satisfy essential boundary conditions and external applied loading can easily be defined. However, in a redundant structure a set of redundant forces must be introduced [4] the selection of which though is difficult to automate. For these reasons stiffness formulation is usually the efficient way to express equilibrium at the structural level. In addition the advantages of the flexibility approach to impose element equilibrium in strong form and derive an exact stiffness matrix in the nonlinear range are explored at the element level. Hence, flexibility approach for the element state determination requires the independency of the stress field, while stiffness based formulation at the structural level suggests the displacement field as independent. Considering the two fields as independent the Hellinger–Reissner energy principle is employed, hence the cross sectional deformations are expressed in two ways; in terms of displacements dðuÞ and as functions of stress resultants dðDÞ. For a fiber beam under uniaxial deformation with stress r and strain e the Hellinger–Reissner functional [30] can be stated in terms of the two independent sets of variables, i.e. the normal stress r and axial displacement u as follows:
PHR ðr; uÞ ¼
Z
V
feðuÞ r vðrÞgdV Pext ðuÞ
ð30Þ
where vðrÞ is the complementary energy density function from which strains can be calculated in terms of stresses, Pext ðuÞ is the functional of external loading and integration is performed at the undeformed volume V of the element. Taking in account Eqs. (22) and (24) the above functional can be written also in stress resultant and deformation terms omitting for simplicity the argument in the integrals i.e. ðDðxÞ ¼ DÞ.
404
I.A. Gkimousis, V.K. Koumousis / Engineering Structures 106 (2016) 399–409
Fig. 6. Basic forces, displacements, deformations and distributed loads on the beam element.
nodal values. For the displacement field the classical cubic interpolation functions N u are considered while deformation interpolation functions N d result from the derivatives of N u (1st Xderivative for the axial strain and second derivative for curvature terms).
8 9 u1 > > > > > > > > > u2 > > > > > > > < ux ðx; yÞ r1 = ) uðxÞ ¼ N u ðxÞ u ¼ N u ðxÞ > uy ðx; yÞ > u4 > > > > > > > > > > u5 > > > : > ; r2 ¼ N d ðxÞ q ¼ N d ðxÞ T u dðuðxÞÞ ¼ @ ðN u ðxÞÞ u
ð35Þ
where
"
Nu ¼
1 xL
0
0
2 xL3 3 xL2 þ 1
x L
0
3
2
x3 L2
0
2
2
#
0 3
2 xL þ x 0 3 xL2 2 xL3
x3 L2
2
xL
ð36Þ
Fig. 7. Displacements and deformations of beam element i.
PHR ðD; uÞ ¼
Z
L
n
0
o dðuÞT D vðDÞ dx Pext ðuÞ
ð31Þ
Complementary energy functional can now be expressed as:
dðDÞ ¼
@ vðDÞ @D
ð32Þ
The functional of external loading for nodal loads F n and distributed
T loads bðxÞ ¼ bx ðxÞ; by ðxÞ obtains the form:
Fn þ Pext ¼ u
Z
L
uT b dx
ð33Þ
0
In order to calculate the state of the element, stationarity of the Hellinger–Reissner functional is imposed by setting its first variation with respect to the two independent fields equal to zero. Then for all variations duðxÞ and dDðxÞ the state defined by ðuðxÞ; DðxÞÞ satisfies the classical equilibrium and strain–displacement compatibility equations.
dPHR ¼
Z
L
Z ddðuÞT D dx þ
0
0
Z
L
T Fn dDT ðdðuÞ dðDÞÞ dx du
L
duT b dx ¼ 0
ð34Þ
0
Then, in the above relation the discretized independent fields are inserted as linear combinations of shape functions and corresponding
For the stress field expressed in stress resultant terms, cross sectional internal forces are calculated from nodal internal forces based on equilibrium considerations. More specifically, equilibrium Eq. (29) are directly integrated and by applying as essential boundary conditions the internal nodal forces, the following interpolation relation is derived.
"
DðxÞ ¼ N D ðxÞ Q þ Dp ;
N D ðxÞ ¼
1 0
0 0 x x 1 L L
#
ð37Þ
where Dp is the particular solution of the direct integration and can be completely determined by the applied loads. Also, the interpolation of the variation of the stress field, since Dp is a known function, has the following form:
dDðxÞ ¼ N D ðxÞ dQ
ð38Þ
Substituting Eq. (35), (37) and (38) into the functional (34) the following discrete form is derived:
Z T TT dPHR ¼ du
L
0 L
Z T
þ dQ 0
N Td N D Q þ Dp dx
N TD
Z 0
L
N Tu b dx F n
dðDÞÞ dx ¼ 0 ðN d T u
ð39Þ
Since Eq. (39) is valid for any arbitrary variation of the two independent fields, both expressions inside the brackets are set equal to zero. Moreover, the shape functions used are orthogonal i.e. RL T N d N D dx ¼ I as geometric nonlinear effects are ignored. Thus 0 the following Euler–Lagrange equations are deduced that describe the state of the element under external loading incorporating the nonlinear material behavior.
I.A. Gkimousis, V.K. Koumousis / Engineering Structures 106 (2016) 399–409
T T Q þ Q p P ext ¼ 0 T u
Z
ð40Þ
L
N TD dðDÞ dx ¼ 0
ð41Þ
where the vector of equivalent nodal external actions has the form:
Z P ext ¼ F n þ
L 0
N Tu b dx
ð42Þ
and Q p is the element vector resulting from the additional cross sectional internal forces due to distributed loading:
Qp ¼
L 0
N Td Dp dx
ð43Þ
Eq. (40) stands for the equilibrium equation and (41) for the compatibility equation. Interelement continuity is imposed by mapping the structure’s global nodal displacements to the element local dis . The only source of nonlinearity in the above equaplacements u tions concerns the derivation of the cross sectional deformations dðDÞ as functions of stress resultants. In the following section the classical linearization technique of the nonlinear constitutive equation is presented, accompanied with a standard Newton–Raphson solution procedure for the solution of the nonlinear system of equations. 4.1. Solution via linearization Various techniques have been proposed for the linearization of the equations that describe the nonlinear problem. Nukala and White [7] have proposed four different solution procedures investigating all different combinations of linearization. Also, Alemdar and White [10] have linearized the constitutive law along with equilibrium and compatibility equations in terms of the independent variables for both geometric and material nonlinearity. The proposed scheme extends the one presented by Hjelmstad and Taciroglu [5] by implementing the rate Eq. (17) describing the constitutive behavior. This is further linearized with respect to the current configuration i following the forward Euler method thus leading to the current tangent fiber modulus expressed as:
h
i
r ¼ ri þ 1 ð1 aÞ Hi1 H i2 E Dei ¼ ri þ Eit Dei
ð44Þ
The current cross sectional flexibility matrix is determined as the inverse of the stiffness matrix using relation (26). Since cross sectional constitutive equation is unique the interpolated forces DðxÞ ^ should be equal to the stress resultants DðxÞ derived from material constitutive law. This equation is linearized as follows:
^ iþ1 ) Ddi ¼ f i Diþ1 D ^i Diþ1 ¼ D
ð45Þ
Consequently, updated cross sectional deformations can be calculated from the following relation:
d
iþ1
L
0
0
Z
Z qi ¼
i ^i ¼ d þ f i Diþ1 D
ð46Þ
Z
L
ci ¼ 0
405
i NTD d Di dx
ð50Þ
^ i þ DDi dx NTD f i Di D p
ð51Þ
i and iþ1 ¼ u i þ Du Considering the expressions Q iþ1 ¼ Q i þ DQ i , u eliminating the increment of nodal forces DQ i , the increment of i can be determined, provided that the elenodal displacements Du ment flexibility matrix F i is invertible. i K i Dui ¼ P iþ1 ext Q int
ð52Þ
h i 1 Q iint ¼ K i ui þ T T Q i ðF i Þ qi þ c i þ Q iþ1 p
ð53Þ
1
where K i ¼ T T ðF i Þ T is the element’s tangent stiffness matrix. Eq. (52) has the same form with the standard incremental equations of the displacement based method the difference being that internal element forces and stiffness matrix are exact as they satisfy equilibrium equations in strong form. On the other hand compatibility is satisfied in weak form as nodal deformations are calculated as the weighted average of the cross sectional deformations. Considering the computational aspects of this formulation, cross sectional residuals are directed to the element’s residuals and then to structure’s residuals. Hence this algorithm is classified as an L–L algorithm according to Nukala and White [7] characterization. Following a Newton–Raphson procedure for the global equilibrium equations, the procedure of calculating updated global stiffness matrix and internal nodal forces proceeds as follows: (a) After solving for an incremental step for the whole structure’s equilibrium equation, the vector of incremental nodal displacements is evaluated and nodal total displacements are updated. Mapping the global structure’s displacements iþ1 of each element to every element, the displacements u are defined and the element state determination begins. (b) Element incremental nodal forces DQ i are calculated from Eq. (48) and total nodal forces are updated Q iþ1 ¼ Q i þ DQ i . (c) Cross-sectional forces are calculated from interpolation: Diþ1 ¼ Di þ N D DQ i þ DDp , where DDp is the increment of cross sectional stress resultants due to distributed loads. (d) Cross-sectional deformations are evaluated from the linearized constitutive Eq. (46) and updated fiber strains eiþ1 are derived from Eq. (22). (e) The linearized material Eq. (44) is solved and updated tangent fiber modulus Eit and fiber stresses riþ1 are evaluated. Then integration over the cross-section evaluates the updated tangent stiffness kðxÞiþ1 which gives the cross sec-
Substituting Eq. (46) in (41) results in the linearized incremental equations:
T T Q iþ1 þ Q piþ1 P iþ1 ext ¼ 0
tional flexibility f ðxÞiþ1 by its inverse. (f) Integration over the length using Eqs. (49) and (50) provides
ð47Þ
element flexibility matrix F iþ1 and element’s nodal deformations qiþ1 .
iþ1 F i DQ i qi c i ¼ 0 T u
ð48Þ
(g) Eq. (53) determines the nodal internal forces Q iþ1 for the next step, while by assembling every element’s contribution
where the element’s flexibility, nodal deformations and element’s deformations residual are derived integrating over the beam length as:
Z
L
Fi ¼ 0
N TD f i N D dx
ð49Þ
the global structure’s stiffness matrix K iþ1 and nodal internal S forces P iþ1 int are derived. (h) Every element’s residual is directed to the structural level
iþ1
P iþ1 and tolerance is checked. P ext int 6 tol.
406
I.A. Gkimousis, V.K. Koumousis / Engineering Structures 106 (2016) 399–409
5. Examples 5.1. Cantilever column In this example the proposed element formulation is validated solving the problem presented in Ref. [6]. In this reference, different numerical procedures and element discretization schemes of a mixed three-field Hu–Washizu element, first proposed in [9] are compared. Also, the proposed steel constitutive law (17) is implemented following the classical displacement based formulation and the results are presented for comparison reasons. The problem consists of a non-dimensional cantilever column with height of 120 units and rectangular cross section of 15 20 units. Material is considered bilinear with elastic modulus of elasticity Eel ¼ 29; 000, yield stress f y ¼ 50 and hardening post-yield to pre-yield ratio a ¼ 0:01. The cross-section is discretized in 5 layers, each one participating in the overall cross-section response according to mid-point integration. Element behavior is integrated along its length using Gauss–Lobatto integration with 5 control points. Implementing the classical displacement based formulation requires several elements in the discretization to achieve a satisfactory and objective response [44]. For this reason two different discretization schemes were tested consisting of one (1DB) and four elements (4DB) respectively. The results are presented in Fig. 8 relating base shear and top displacement, with the proposed hysteretic formulation proven very efficient in generating the accurate response. Also, the proposed two-field mixed element (Proposed-1HR) is compared for the same example against the three-field force based element of Ref. [6] (reference-1FB) and the results are quite satisfactory as presented in Fig. 9.
Fig. 9. Base shear vs top displacement of a cantilever column, comparison with mixed formulations.
Fig. 10. Continuous beam under vertical dynamic excitation.
5.2. Continuous beam under cyclic vertical load In this example the proposed formulation is employed for the solution of a 2-span continuous beam (Fig. 10) under a pseudostatic cyclic external loading consisting of a varying concentrated vertical load at the middle of the first span following a piecewise linear pattern presented in Fig. 11 and Table 2. Also, cross section is an IPE 300 European steel cross section which is discretized in 50 layers. The comparison is twofold, first the displacement based approach is compared with the proposed mixed formulation, and then the proposed hysteretic steel model is compared with the Menegotto–Pinto model. Material and model properties are given in Table 3. In Fig. 12 the two element formulations are compared. For the displacement based formulation two different discretization
Fig. 11. Vertical load input.
Table 2 External loading values. Time (s)
External force (KN)
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 5.5
0 400 300 200 150 250 50 300 0 350 50 380
Table 3 Steel properties and model parameters. Steel properties
Fig. 8. Base shear vs top displacement of a cantilever column, comparison with displacement based formulation.
Model parameters
Young modulus E Yield stress fy
210 GPa 235 MPa
Kinematic hardening ratio a
0.015
Proposed
Menegotto–Pinto
n1 = 10 n2 = 1
R0 = 21 cR1 = 17.5 cR2 = 0.15
I.A. Gkimousis, V.K. Koumousis / Engineering Structures 106 (2016) 399–409
407
Fig. 12. Comparison of vertical midspan displacement vs external loading history. Fig. 15. Northridge_TAR090 ground motion excitation.
a single element captures directly the real cyclic response with no need for further refinement. In Fig. 13 the proposed steel hysteretic model is compare with the Menegotto–Pinto model. The advantage of the proposed model to avoid overshooting after short reversals is reflected in this case also in the overall behavior of the structure. More specifically, after a series of plastic short reversals, which occur often in real seismic excitations, Menegotto–Pinto model underestimates structural displacements. Hence, the proposed model produces results, which are more accurate and on the safe side from the engineering practice perspective. 5.3. Steel frame under seismic excitation Fig. 13. Load–displacement history, comparison of the proposed model with Menegotto–Pinto model.
schemes with 3 control sections are utilized, while for the mixed formulation a single element discretization with 15 control sections is sufficient to yield accurate results. A dense discretization consisting of 40 equal beam elements with 3 Gauss points per element for the total length of the continuous beam is used to overcome the deficiencies of the classical displacement based stiffness method. Moreover, steel behavior is modeled according to the proposed hysteretic model. It is evident that the discretization with 3 elements following the stiffness method (DB 3 elements) delivers inaccurate results and is not capable of capturing the real behavior. This inefficiency is magnified for intensive external loading, as overestimation of the maximum strength results into more than 2.5 times decreased plastic deformations compared to a more refined discretization. On the other hand the proposed two-field mixed element (HR 3 elements) satisfies the equilibrium equations in strong form and
In this example a 2 story–3 bay steel plane frame (Fig. 14) is excited with the Northridge_TAR090 seismic strong motion record (Fig. 15). The frame is designed according to Eurocode 3 specifications and consists of HEA-450 columns and IPE-220 beams discretized in 50 layers, while steel material class is S355 with 1.5% kinematic hardening. Also, distributed load of G = 15.0 KN/m and the associated mass is considered on every beam including selfweight. Lumped Mass and no viscous damping are considered, while direct time integration was performed according to Newmark numerical scheme with beta = 0.25, gamma = 0.5 and time step dt = 0.01. This example is used to compare the results of the proposed hysteretic model with the Menegotto–Pinto steel model for the response of the structure in terms of accuracy and computing time. The proposed mixed element is used in both simulations and the model parameters for the two steel models are those presented in Table 3. First, the overall response of both models is compared by monitoring the top story lateral displacement for the first 10 s of the
Fig. 14. 2 story–3 bay steel plane frame.
408
I.A. Gkimousis, V.K. Koumousis / Engineering Structures 106 (2016) 399–409
20% as compared to Menegotto–Pinto model. The comparison is performed using a general program for frame structural analysis developed in Matlab which accounts for the two different steel models. The time step of the analysis is dt = 0.01 s, while similar deviations were obtained also for different time steps, indicating that both models demonstrate the same tendency for convergence.
6. Conclusions
Fig. 16. Top lateral displacement history for the proposed Hysteretic and Menegotto–Pinto models.
Fig. 17. Moment–curvature diagrams of beam section A, using the proposed Hysteretic and Menegotto–Pinto models.
In this work a new hysteretic uniaxial steel model is proposed describing steel cyclic behavior, which is implemented within a fiber beam–column element based on Hellinger–Reissner mixed variational principle. The hysteretic model is developed on the basis of classical plasticity considerations and has the advantage of embodying all nonlinear phases in a single nonlinear differential equation. The model is further modified to account for steel kinematic hardening, Bauschinger effect and compliance with plasticity postulates after short reversals. The model avoids overshooting after short reversals, being on the same time computationally effective and in that respect more advantageous compared to the well-established Menegotto–Pinto model. Incorporating the constitutive equations into the general formulation of mixed beam finite elements, element equilibrium and compatibility equations are derived. Linearization of the differential constitutive equation is employed, which leads to the incremental equilibrium and compatibility equations, which are solved using a Newton’s type iterative numerical scheme. The main feature and computational advantage of the proposed scheme relies on the explicit and direct calculation of stresses, with no need for predictor–corrector schemes by solving the hysteretic constitutive equations in every time step. Verification of the proposed model was performed on the basis of displacement based or mixed methods using examples from the literature demonstrating the computational efficiency of the proposed formulation.
References
Fig. 18. Total computational time for hysteretic and Menegotto–Pinto models.
excitation as presented in Fig. 16. Minor differences between the two models are observed in this case. This happens besides the fact that certain elements demonstrate a more pronounced difference as observed in Fig. 17, which at the structural level seem to compensate. Also, as for the vertical displacements in example 2, curvatures in this case are underestimated by the Menegotto–Pinto model. After several cycles of excitation, the critical section is oscillating elastically at larger negative curvatures according to the proposed model compared with the Menegotto–Pinto model. The main advantage of the proposed formulation, which avoids overshooting, relies on its computational efficiency. The total computing time for four different schemes with 3, 5, 7, 9 control sections per element is recorded and the results are plotted in Fig. 18. The reduction of computational time is of the order of
[1] Bathe KJ. Finite element procedures. New York: Prentice Hall Engineering, Science, Mathematics; 2007. [2] Spacone E, Filippou FC, Taucer FF. Fibre beam column model for nonlinear analysis of R/C frames I: formulation. Earthq Eng Struct Dyn 1996;25 (7):711–25. [3] Neuenhofer A, Filippou FC. Evaluation of nonlinear frame finite-element models. J Struct Eng 1997;123(7):958–66. [4] Hjelmstad KD, Taciroglu E. Mixed methods and flexibility approaches for nonlinear frame analysis. J Constr Steel Res 2002;58:967–93. [5] Hjelmstad KD, Taciroglu E. Variational basis of nonlinear flexibility methods for structural analysis of frames. J Eng Mech 2005;131(11):1157–69. [6] Saritas A, Soydas O. Variational base and solution strategies for the non-linear force-based beam finite elements. Int J Non-linear Mech 2012;47:54–64. [7] Nukala PK, White DW. Variationally consistent state determination algorithms for nonlinear mixed beam finite elements. Comput Methods Appl Mech Eng 2004;193:3647–66. [8] Hjelmstad KD, Taciroglu E. Mixed variational methods for finite element analysis of geometrically non-linear, inelastic Bernoulli–Euler beams. Commun Numer Methods Eng 2003;19(10):809–32. [9] Taylor RL, Filippou FC, Saritas A, Auricchio F. A mixed finite element method for beam and frame problems. Comput Mech 2003;31:192–203. [10] Alemdar BN, White DW. Displacement, flexibility and mixed beam–column finite element formulations for distributed plasticity analysis. J Struct Eng 2005;131:12. 1811. [11] Alsafadie R, Hjiaj M, Somja H, Battini J-M. A comparative study of displacement and mixed-based corotational finite element formulations for elasto-plastic, three-dimensional beam analysis. Eng Comput: Int J Comput-Aided Eng Softw 2011;28(7):939–82. [12] Correia AA, Almeida JP, Pinho R. Force-based higher-order beam element with flexural-shear-torsional interaction in 3D frames. Part I: theory. Eng Struct 2015;89:204–17. [13] Simo JC, Taylor RL. Consistent tangent operators for rate – independent elastoplasticity. Comput Methods Appl Mech Eng 1985;48:101–18. [14] Simeonov VK, Sivaselvan MV, Reinhorn AM. Nonlinear analysis of structural frame systems by the state-space approach. Comput Aided Civ Infrastruct Eng 2000;15:76–89.
I.A. Gkimousis, V.K. Koumousis / Engineering Structures 106 (2016) 399–409 [15] Jafari V, Rahimian M, Vahdan SH. Improved formulation in mixed-based statespace approach for large displacement inelastic analysis of frames. J Eng Mech 2011;137(3):196–204. [16] Triantafyllou S, Koumousis V. Small and large displacement dynamic analysis of frame structures based on hysteretic beam element. J Eng Mech 2011;138 (1):36–49. [17] Triantafyllou SP, Koumousis VK. Hysteretic finite elements for the nonlinear static and dynamic analysis of structures. J Eng Mech 2013;140(6). [18] Menegotto M, Pinto PE. Method of analysis for cyclically loaded R.C. plane frames including changes in geometry and non-elastic behaviour of elements under combined normal force and bending. In: Symposium on the resistance and ultimate deformability of structures acted on by well defined repeated loads. International association for bridge and structural engineering (IABSE) Lisbon, Portugal; 1973. p. 15–22. [19] Ramberg R, Osgood WR. Description of stress–strain curve by three parameters. Technical note 902. NACA; 1943. [20] Monti G, Nuti C. Nonlinear behaviour of reinforcing bars including buckling. J Struct Eng 1992;118(12):3268–84. [21] Kashani MM, Crewe AJ, Alexander NA. Nonlinear stress–strain behavior of corrosion-damaged reinforcing bars including inelastic buckling. Eng Struct 2013;48:417–29. [22] Kunnath S, Heo Y, Mohle J. Nonlinear uniaxial material model for reinforcing steel bars. J Struct Eng 2009;135(4):335–43. [23] Zong Z, Kunnath S, Monti G. Material model incorporating buckling of reinforcing bars in RC columns. J Struct Eng 2014;140(1). [24] Kashani MM, Crewe AJ, Alexander NA. Nonlinear stress–strain behaviour of corrosion-damaged reinforcing bars including inelastic buckling. Eng Struct 2013;48:417–29. [25] Shi G, Wang M, Bai Y, Wang F, Shi YJ, Wang YQ. Experimental and modeling study of high-strength structural steel under cyclic loading. Eng Struct 2012;37:1–13. [26] Wang YB, Li GQ, Cui W, Chen SW. Experimental investigation and modeling of cyclic behavior of high strength steel. J Constr Steel Res 2015;104:37–48. [27] Filippou FC, Popov EP, Bertero VV. Effects of bond deterioration on hysteretic behavior of reinforced concrete joints. Report. UCB/EERC-83/19. Berkeley: Earthquake Engineering Research Centre, University of California.
409
[28] Dodd LL, Restrepo-Posada JI. Model for predicting cyclic behaviour of reinforcing steel. J Struct Eng ASCE 1995;121(3):433–45. [29] Balan TA, Filippou FC, Popov EP. Hysteretic model of ordinary and highstrength reinforcing steel. J Struct Eng 1998;124(3):288–97. [30] Reissner E. On a variational theorem in elasticity. J Math Phys 1950;29:90–5. [31] Valanis KC, lee CF. Endochronic theory of cyclic plasticity with applications. J Appl Mech 1984;51(2):367–74. [32] Sivaselvan MV, Reinhorn AM. Nonlinear Structural analysis towards collapse simulation – a dynamical systems approach. Technical report, Buffalo, New York; 2003. [33] Wen YK. Method of random vibration of hysteretic systems. J Eng Mech Div 1976;102:249–63. [34] Lubliner J. Plasticity theory. New York: Dover Publications, Inc.; 2008. [35] Borja RJ. Plasticity, modeling & computation. Berlin, Heidelberg: SpringerVerlag; 2013. [36] Ucak A, Tsopelas P. Accurate modeling of the cyclic response of structural components constructed of steel with yield plateau. Eng Struct 2012;35:272–80. [37] Ayhan B, Jehel P, Brancherie D, Ibrahimbegovic A. Coupled damage-plasticity model for cyclic loading: theoretical formulation and numerical implementation. Eng Struct 2013;50:30–42. [38] Mendes LAM, Castro LMSS. A simplified reinforcing steel model suitable for cyclic loading including ultra-low-cycle fatigue effects. Eng Struct 2014;68:155–64. [39] Thyagarajan RS. Modeling and analysis of hysteretic structural behavior. Technical report. Caltech EERL, EERL-89-03; 1989. [40] Wang CH, Foliente GC. Discussion on hysteretic models for deteriorating inelastic structures. J Eng Mech ASCE 2001;127(11):1200–2. [41] Charalampakis AE, Koumousis VK. A Bouc–Wen model compatible with plasticity postulates. J Sound Vib 2009;322(4-5):954–68. [42] Kottari AK, Charalampakis AE, Koumousis VK. A consistent degrading Bouc– Wen model. Eng Struct 2014;60:235–40. [43] Thompson KJ, Park R. Stress–strain model for grade 275 reinforcing steel with cyclic loading. Bull NZ Natil Soc Earthq Eng 1978;11(2). [44] Calabrese A, Almeida JP, Pinho R. Numerical issues in distributed inelasticity modeling of RC frame elements for seismic analysis. J Earthq Eng 2010;14:38–68.