A consistent Timoshenko hysteretic beam finite element model

Accepted Manuscript A consistent Timoshenko hysteretic beam finite element model M. Amir, K.G. Papakonstantinou, G.P. Warn

PII: DOI: Reference:

S0020-7462(19)30083-6 https://doi.org/10.1016/j.ijnonlinmec.2019.07.003 NLM 3218

To appear in:

International Journal of Non-Linear Mechanics

Received date : 4 February 2019 Revised date : 3 June 2019 Accepted date : 7 July 2019 Please cite this article as: M. Amir, K.G. Papakonstantinou and G.P. Warn, A consistent Timoshenko hysteretic beam finite element model, International Journal of Non-Linear Mechanics (2019), https://doi.org/10.1016/j.ijnonlinmec.2019.07.003 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

A consistent Timoshenko hysteretic beam finite element model M. Amira,∗, K.G. Papakonstantinoub , G.P. Warnc a Graduate

Student Researcher, Dept. of Civil and Environmental Engineering, The Pennsylvania State University, University Park, PA 16802, USA. b Assistant Professor, Dept. of Civil and Environmental Engineering, The Pennsylvania State University, University Park, PA 16802, USA. c Associate Professor, Dept. of Civil and Environmental Engineering, The Pennsylvania State University, University Park, PA 16802, USA.

Abstract A parametrized Timoshenko hysteretic beam finite element model is developed, using consistent two-node elements, to efficiently simulate the nonlinear behavior of structures. New displacement interpolation functions are derived that satisfy both the exact equilibrium and kinematic conditions. Therefore, strain field continuity along the beam element is achieved, resulting in a consistent formulation without shear locking effects. Multiaxial interactions are considered through yield/capacity functions and distributed plasticity is accounted for by appropriate hysteretic interpolation functions. The suggested interpolation functions yield constant element matrices that do not require updating throughout the analysis, and the nonlinearities are captured in terms of hysteretic curvatures, axial and shear deformations, which are set to evolve through ordinary differential equations. A computationally efficient solution scheme is also suggested, without requiring any linearization, ∗ Corresponding

author Email addresses: [email protected] (M. Amir), [email protected] (K.G. Papakonstantinou), [email protected] (G.P. Warn)

Preprint submitted to International Journal of Non-Linear Mechanics

July 9, 2019

to straightforwardly solve the resulting system of equations for quasi-static problems. The consistency, efficiency, versatility and validity of the suggested model is explained in detail and demonstrated through several numerical examples and comparisons with experimental data from available tests. Keywords: Hysteretic finite element, multiaxial interactions, distributed plasticity, ODEs/DAEs. 1

1. Introduction

2

In this paper, a consistent and efficient Timoshenko hysteretic beam finite element

3

model is developed. The presented approach deviates from traditional nonlinear

4

structural modeling techniques [1], including discrete plasticity uniaxial nonlinear

5

spring models [2], fiber-based approaches [3, 4], and solid finite elements mod-

6

eling, and suggests an alternative hysteretic framework, conveniently based on

7

spatially distributed, multiaxial, interacting stress resultants. The formulation is

8

largely inspired by the notable works of Triantafyllou and Koumousis [5, 6], yet this

9

paper suggests a newly developed formulation, mostly related to the consistency

10

and efficiency of the element.

11

The phenomenon of hysteresis appears in many scientific areas and a large

12

number of theoretical and mathematical models can be found in the literature,

13

e.g. [7–9]. Among the various available hysteretic models for structural systems,

14

e.g. [10–12], a considerably extended multiaxial variant of the smooth Bouc-Wen

15

model [13, 14] is emphasized and utilized herein, although the element formulation

16

idea and representation are general and alternative hysteretic models could also be

17

used. As shown in [15–17] however, the Bouc-Wen model expressions have the

18

significant attributes that they can be also derived following the physics of classical 2

19

multiaxial plasticity [18], and very compactly incorporate and satisfy the loading

20

rate, the yield/capacity criterion and the flow rule.

21

The original hysteretic Bouc-Wen model was initially developed by Bouc [13]

22

and later modified by Wen [14], and is described by first order nonlinear ordinary

23

differential equations (ODEs), often referred to as evolution equations. Over the

24

years, several extensions, modifications and variations have been suggested, e.g.

25

[19–25], and the model and its many variants have been extensively used in various

26

applications to simulate the response of diverse structural components/devices, as

27

for example described in [26]. Further extensions are also possible, for example,

28

relating hysteretic damping to viscoelasticity as suggested by [27–29]. Although

29

the original Bouc-Wen model and the vast majority of works in the literature only

30

consider a uniaxial formulation, in Casciati [15] and Sivaselvan and Reinhorn [16]

31

the model is elegantly formulated in a multiaxial plasticity context.

32

More recently, Triantafyllou and Koumousis [5, 6] suggested further multiaxial

33

plasticity enhancements and interactions, and developed an integrated hysteretic

34

beam-column element formulation, able to characterize both Timoshenko and

35

Euler-Bernoulli beam elements. In the considered beam element formulation,

36

hysteretic degrees of freedom are introduced, in terms of hysteretic curvatures, ax-

37

ial and shear deformations, and are set to evolve according to multiaxial Bouc-Wen

38

type evolution equations. The described element in Triantafyllou and Koumousis

39

[6] however, only satisfies weak equilibrium conditions for the Timoshenko case

40

and thus mesh convergence and universal/general performance for both Euler-

41

Bernoulli and Timoshenko beam types are importantly negatively affected. In this

42

paper, the exact equilibrium and kinematic conditions are simultaneously satisfied

43

due to the introduction of newly derived displacement field interpolation functions, 3

44

for a nonlinear Timoshenko beam element subjected to nodal forces. Therefore,

45

displacement and strain field continuity along the beam element are achieved,

46

resulting in a consistent formulation without any shear locking effects, regardless

47

of the element discretization, boundary conditions and beam thickness [30]. Dis-

48

tributed loads along the length of the element can be treated using the standard

49

finite element approach of equivalent nodal loads [31]. In its most comprehensive

50

form, of an inelastic Timoshenko beam formulation, the herein suggested element

51

can now also seamlessly provide compatible and accordant results to relevant

52

slender Euler-Bernoulli cases, and thus this most extensive Timoshenko format is

53

mainly described in this paper. Nonetheless, due to its parametrized nature, the

54

developed element is also directly transformable, by simple parametric changes, to

55

elastic Timoshenko and Euler-Bernoulli beam elements, exhibiting either elastic

56

or inelastic material behavior, and therefore can also be considered as a general

57

beam element model. Additionally, the advantages of the described element in

58

Triantafyllou and Koumousis [6] are also preserved in our suggested approach,

59

full interaction between axial and shear forces and moments can be considered, if

60

needed, through yield/capacity functions, and distributed plasticity is accounted

61

for by appropriate hysteretic interpolation functions.

62

One of the major advantages of the formulation in [6], and in this paper, is its

63

computational efficiency, since the time-varying element tangent stiffness matrix

64

is replaced by a constant elastic stiffness matrix and a hysteretic matrix, which both

65

need to be evaluated only once at the beginning of the analysis and do not need to be

66

updated, given that inelasticity is treated through the localized hysteretic evolution

67

equations. The whole formulation can be thus conveniently presented in state-

68

space form for dynamic cases, as shown in [5, 6], and be efficiently solved by any 4

69

generic ODE solver, avoiding any type of linearizations. In quasi-static problems,

70

however, the solution process is not as straightforward due to the presence of

71

algebraic equations now, in addition to the ODEs, resulting in a set of complex

72

and hard to solve differential algebraic equations (DAEs). Yet, a new and efficient

73

solution scheme for such problems is presented in this paper, where instead of

74

DAEs only a set of ODEs are again required to be solved, similar to the dynamic

75

cases. Through the suggested approach, the number of unknowns to be updated

76

at each pseudo time step of the resulting ODEs are also significantly reduced,

77

in comparison to both DAE systems and the typically used Newton’s method,

78

reducing computational complexity and increasing efficiency even further.

79

To summarize, the presented formulation is able to simulate the nonlinear

80

behavior of beam elements, using distributed plasticity and interaction effects, with

81

sufficient accuracy and efficient computational effort. As previously mentioned,

82

accuracy is achieved through the formulation of exact interpolation fields for the

83

nonlinear element, while efficiency is achieved through the direct use of stress

84

resultants at the section level, constant matrices and the newly devised solution

85

procedure. Thus, the suggested modeling scheme can be importantly employed

86

for the nonlinear static/dynamic analysis of large-scale frame structures, as well as

87

for system identification techniques, among many other applications dictating and

88

demanding computational efficiency.

89

In the remainder of this paper, the generality, consistency, efficiency, versa-

90

tility and validity of the suggested model is explained in detail and demonstrated

91

through several numerical examples and comparisons with experimental data from

92

available tests.

5

93

2. Hysteretic model

94

2.1. Uniaxial hysteretic model

95

Among the various hysteretic modeling approaches, a Bouc-Wen type model, based

96

on stress resultants, is adopted in this work, as already explained and justified. In

97

its simplest form, the hysteretic behavior is modeled through a set of two parallel

98

components; a linear one, contributing to the post-elastic kinematic hardening,

99

and a nonlinear, describing the hysteretic behavior according to a first order ODE,

100

otherwise known as a hysteretic evolution equation. As such, a uniaxial force

101

(P)-displacement (u) relation or any other work conjugate pair can be expressed

102

as:

n  z Û uÛ P = αku + (1 − α)k z ; zÛ = 1 − h (β + γsgn(zu)) z 

(1)

c

103

where k is the elastic loading stiffness, zch is the hysteretic displacement capacity,

104

α is the ratio of post-elastic to elastic stiffness, z is a hysteretic variable, β, γ,

105

n are model parameters, sgn(.) is the signum function and overdot indicates

106

differentiation with respect to time. Parameter n controls the smoothness of the

107

transition from elastic to inelastic regimes, whereas parameters β and γ control the

108

shape of the hysteretic loop. For the specific case in which the elastic unloading

109

stiffness is equal to the elastic loading stiffness, it can be proven that both β and γ

110

are equal to 0.5. Furthermore, to be thermodynamically admissible, the following

111

constraints should be imposed on the parameters β and γ [32]: −γ ≤ β ≤ γ

;

β+γ =1

(2)

112

In this paper, a generalized law is adopted that equivalently to Eq. (1) maps

113

generalized strains to their corresponding actions, in particular, centerline axial

114

strain (εu ) to axial force (N), curvature (εφ ) to moment (M), and shear strain (εγ ) 6

115

to shear force (Q) [6]. Hence, at any distance x along a beam element of length L,

116

the local M-εφ relation is given by: M(x) = αφ E Iεφ(x) + (1 − αφ )E I zφ(x) h n    M h zÛφ(x) = 1 − h β + γsgn(M εÛφ ) εÛφ(x) M

(3)

c

117

where

118

moment, Mch = (1 − αφ )Mp is the hysteretic bending capacity, Mp is an estimated

Mh

= (1 − αφ )E I zφ represents the hysteretic component of the bending

119

plastic moment capacity of the section, E is the elastic modulus, I is the moment

120

of inertia, zφ is the hysteretic bending deformation, and αφ is the post-elastic to

121

elastic stiffness ratio for bending. Similar equations are also available for N-εu ,

122

and Q-εγ , which are omitted here for brevity. Combining the expressions for axial

123

force, bending moment and shear force, results in the following: e h P(x) = P(x) + P(x) = De ε(x) + HD z(x)

124

(4)

P = {N Q M }T ; ε = {εu εγ εφ }T ; z = {zu zγ zφ }T D=

"

EA

G As

EI

#

;

α=

"

αu

αγ αφ

#

;

De = αD

(5)

HD = (I − α)D

125

where I is the identity matrix, De is the elastic rigidity matrix, HD is the hysteretic

126

rigidity matrix, zu and zγ are the hysteretic axial and shear deformations, respec-

127

tively, A is the cross-section area, As is the effective shear area, G is the shear

128

modulus, αu and αγ are the axial and shear hardening parameters, respectively.

129

2.2. Multiaxial model with interaction

130

In the previous subsection, the axial force, moment and shear force have been

131

assumed to evolve following independent uniaxial plasticity laws, such that the 7

132

individual capacities are not constrained by a multiaxial yield/capacity surface.

133

If needed, the evolution laws can be modified to incorporate axial-moment-shear

134

interactions as follows [6]:

n zÛ (x) = (I − H1 H2 R) εÛ (x) ; H1 = Φ(Ph ) + 1   Ph = {N h Q h M h }T ; H2 = β + γsgn (Ph )T εÛ " T   # −1 "   T # ∂Φ ∂Φ ∂Φ ∂Φ D D R= ∂Ph ∂Ph ∂Ph ∂Ph

(6)

135

where H1 is a smooth function ranging in [0,1] corresponding to the elastic and

136

plastic regions respectively, H2 is a Heaviside function, Ph is the hysteretic com-

137

ponent of the axial force, shear force, bending moment, and Φ(Ph ) is the specified

138

yield/capacity function. The expression for Φ(Ph ) can be specified from pre-

139

defined yield/capacity functions, depending on the material and cross-sectional

140

properties of the beam element, e.g. [33] among many others. It should be also

141

noted, however, that an expression for Φ(Ph ) can also be numerically derived using

142

fiber analysis, [34] for example, and the function H1 and the matrix R are then

143

evaluated accordingly. Similarly, axial-moment interaction can be incorporated by

144

modifying the evolution equations as follows:       εÛu  

      zÛu  

= (I2×2 − H1 H2 R)      εÛφ   zÛφ    (x)   (x) h n    Q h zÛγ(x) = 1 − h β + γsgn(Q εÛγ ) εÛγ(x) Qc

(7)

145

where functions H1 , H2 and the interaction matrix R are dependent on N h and M h

146

now, whereas Q h evolves independently based on the hysteretic shear capacity Q ch .

8

147

3. Suggested hysteretic finite element formulation

148

Friedman and Kosmatka [35] derived transverse displacement and rotation inter-

149

polation functions for linear Timoshenko beam elements that exactly satisfy the

150

differential equations associated with Timoshenko’s beam theory. These interpo-

151

lation functions, also used in Triantafyllou and Koumousis [6], are not entirely

152

compatible with the hysteretic formulations due to the additional hysteretic defor-

153

mation components z. Therefore in the present work, new consistent interpolation

154

functions are explicitly derived for a two-noded beam element, that comply with

155

the equilibrium equations of homogeneous nonlinear Timoshenko beam elements.

156

Six additional DOF are considered to account for hysteretic axial, bending and

157

shear deformations, both at the start (Node 1) and the end nodes (Node 2) of the

158

beam element. Hence, the nodal displacement and hysteretic DOF vectors can be

159

expressed as: d = {u1 w1 θ 1 u2 w2 θ 2 }T ; z = {zu1 zγ1 zφ1 zu2 zγ2 zφ2 }T

(8)

160

where d is the nodal displacement vector in local coordinates, consisting of longi-

161

tudinal displacement u, transverse displacement w, and rotation θ for both nodes,

162

and z represents the hysteretic deformation vector or else the six additional DOF,

163

each corresponding to their respective strains at the two nodes. To account for

164

distributed plasticity over the length of the beam element, compatible interpolation

165

functions are also suggested for z. A detailed description is presented in the sub-

166

sequent subsections, starting from the derivation of the interpolation functions for

167

displacement and hysteretic deformation, to the final formulation of the element

168

matrices using the variational principle of virtual work.

9

169

3.1. Beam kinematics and equilibrium equations

170

Timoshenko beam theory assumes that in the deformed configuration the cross

171

sections of a beam element remain plane but not necessarily perpendicular to the

172

longitudinal axis. Hence, for a 2D homogeneous beam element of length L in a

173

x − y plane, as in Fig. 1, the kinematic relations are expressed as: εu(x) =

∂u(x) ∂w(x) ∂θ (x) ; εγ(x) = − θ (x) ; εφ(x) = ∂x ∂x ∂x

(9)

174

In this formulation, only nodal loading is considered for the beam element, resulting

175

in the following equilibrium equations: ∂N(x) =0 ∂x

176

177

where N(x) = αu E Aεu(x) + (1 − αu )E Azu(x)

∂Q(x) = 0 where Q(x) = αγ G As εγ(x) + (1 − αγ )G As zγ(x) ∂x ∂ M(x) + Q(x) = 0 where M(x) = αφ E Iεφ(x) + (1 − αφ )E I zφ(x) ∂x

(10) (11) (12)

178

Eqs. (9) to (12) represent the kinematic conditions and the fundamental equilibrium

179

equations of the hysteretic Timoshenko beam element, without any nonlinear

180

geometric effects.

181

3.2. Consistent interpolation functions with distributed nonlinearity

182

In order to satisfy the equilibrium conditions in Eqs. (10) to (12), εu(x) , εγ(x) , zu(x)

183

and zγ(x) are assumed to be constant, whereas εφ(x) and zφ(x) are considered to

184

vary linearly along the length of the element. Therefore, based on the hysteretic

185

boundary conditions in Eq. (8), the following hysteretic axial, shear and bending

10

∂w ∂x εγ θ

y

y

∂w ∂x

w Q,w

M,θ N,u

x

x

(a)

(b)

Fig. and (b) (b) deformed deformed configuration. configuration. Fig.1.1. Beam Beamelement: element: (a) (a) undeformed undeformed configuration, configuration, and 104 186

187

188

3.2. Consistentareinterpolation functions with distributed nonlinearity deformations obtained at any distance x along the element:     In order to satisfy the equilibrium conditions 1 1in Eqs. (10) to (12), εu(x) , εγ(x) , zu(x) zu(x) = zu1 + zu2 2 whereas 2 εφ(x) and zφ(x) are considered to and zγ(x) are assumed to be constant,     1 1 (13) zγ1 + zγ2 zγ(x) of = the element. vary linearly along the length Therefore, based on the hysteretic 2 2  following x boundary conditions in Eq. (8), the hysteretic axial, shear and bending x zφ(x) = 1 − zφ1 + zφ2 L x alongLthe element: deformations are obtained at any distance By substituting Eqs. (9) and (13) inEqs. (10)  to (12), the following fundamental 1 1 zu(x)hysteretic = zu1Timoshenko + zu2 beam element are obtained: differential equations for the 2 2     ! ! 1∂ 2 w 1 (13) ∂ 2 u(x) ∂θ (x) zγ1 + zγ2(x) zγ(x) = = 0 ; α G A − = 0 αu E A 2 2 γ s ∂x ∂ x2 ∂ x 2  x x + z  zφ2 φ1 − ∂ 2 θ (x)zφ(x) = 1 − L  zzφ2 φ1 (14) L + αφ E I + (1 − αφ )E I 2 L ∂x    ∂w(x) zγ2 + zγ1 αγ G As − θ (x) + (1 − αγ )G As =0 ∂x 2

11

10

189

Based on the typical beam element assumptions, the displacement fields along the

190

element are assumed as: u(x) = B0 + B1 x w(x) = C0 + C1 x + C2 x 2 + C3 x 3

(15)

θ (x) = D0 + D1 x + D2 x 2 191

where the coefficients Bi (i ∈ {0, 1}), C j ( j ∈ {0, 1, 2, 3}) and D k (k ∈ {0, 1, 2})

192

are derived using appropriate displacement boundary conditions in Eq. (8) and

193

satisfying Eq. (14). Therefore, the overall displacement and rotational fields can

194

be expressed as:         N1 0 0 N2 0 u 0              =  0 N3 N4 0 N5 N6  d w             θ  0 N7 N8 0 N9 N10     (x)    0 0 0 0 0 0      + 0 N11 N12 0 N13 N14  z     0 N15 N16 0 N17 N18              Nu  HN u                  = Nw d + HN w z                 Nθ   HN θ      

12

(16)

195

196

where the derived displacement interpolation functions, Ni (i ∈ {1, 2, ... , 18} ), are provided in Eq. (17):



N1 = 1 − Lx ; N2 = Lx   3α 12α λ 2α N3 = L 3γ x 3 − L 2γ x 2 − Lφ x µ0 + 1   α 2(α +3α λ) N4 = Lγ2 x 3 − γ L φ x 2 + (αγ + 6αφ λ)x µ0   3α 12α λ 2α N5 = − L 3γ x 3 + L 2γ x 2 + Lφ x µ0   α (α −6α λ) N6 = Lγ2 x 3 − γ L φ x 2 − 6αφ λx µ0   6 2 6 N7 = L 3 x − L 2 x αγ µ0   3α 4(α +3α λ) N8 = L 2γ x 2 − γ L φ x µ0 + 1   6 2 6 N9 = − L 3 x + L 2 x αγ µ0   3α 2(α −6α λ) N10 = L 2γ x 2 − γ L φ x µ0   3 2 x − 12 x (1 − αγ )µ0 N11 = N13 = − L12 x 3 + 2L   N12 = L2 x 3 − 3x 2 + L x (1 − αφ )µ0λ   N14 = − L2 x 3 + 3x 2 − L x (1 − αφ )µ0λ   3 2 3 N15 = N17 = − L 2 x + L x (1 − αγ )µ0   N16 = L6 x 2 − 6x (1 − αφ )µ0λ   6 2 N18 = − L x + 6x (1 − αφ )µ0λ

and µ0 = 197

1 αγ +12αφ λ

; λ=

(17)

EI G As L 2

Similar to the linear Timoshenko beam element case, functions N1 − N10 map the

198

conventional nodal displacement and rotational DOF to the displacement fields.

199

To satisfy the equilibrium conditions and obtain consistent beam elements, the

200

additional functions, N11 − N18 , are also needed, connecting the displacement

201

202

fields to the introduced nodal hysteretic DOF as well.

Although the presented interpolation functions in Eq. (17) are particularly 13

203

derived for a fully nonlinear Timoshenko element, they are sufficiently general,

204

having the notable property of being entirely compatible to all other simpler cases

205

by an appropriate and simple choice of parameter values. For example, for the case

206

of linear beam elements, when αu = αγ = αφ = 1, the developed expressions for

207

N1 − N10 revert to the interpolation functions suggested by Friedman and Kosmatka

208

[35], and the hysteretic components of the interpolation functions, N11 − N18 ,

209

become equal to zero, nullifying the hysteretic DOF contributions. Similarly,

210

for the case of inelastic Euler-Bernoulli beam element, by specifying λ = 0 and

211

αγ = 1, the suggested interpolation functions coincide with the typically used ones

212

for these elements, as also proposed by Triantafyllou and Koumousis [5].

213

214

Similar to the displacement fields, the hysteretic deformation fields from Eq. (13), can be expressed in matrix form as: z(x)

 1 0   z   2  u     = zγ =  0 12      zφ      (x)  0 0

1  0 2 0 0  0 0 12 0  z = Nz(x) z 
1 − Lx 0 0 Lx 

(18)

215

where Nz is the hysteretic interpolation matrix, with each element representing a

216

hysteretic interpolation function. From the derived displacement and hysteretic

217

interpolation functions, as shown in Eqs. (16) to (18), generalized strain fields

218

can be obtained at any distance x along the beam element, using the kinematic

219

relations of Eq. (9). Substituting Eq. (16) to Eq. (9), results in the following

14

220

strain-displacement relations:

ε(x)

221

222

where B(x)

∂u        ε       u     ∂x   ∂w = εγ = ∂x − θ         εφ    ∂θ    (x)  ∂ x  ∂HNu ∂Nu              ∂N ∂ x       ∂H ∂ x  Nw w = ∂ x − Nθ d + − H N θ z = B(x) d + HB(x) z ∂ x        ∂HNθ  ∂Nθ     ∂x ∂x    

(19)

T − 1 0 0   L   6αγ µ0 (L−2x) 12αφ λµ0   0 − L − 3 L    0 −6α λµ0 6αγ µ0 x − 4µ0(αγ +3λαφ )    φ 2 L L B(x) =  1 (20)    L 0 0   6αγ µ0 (L−2x)   0 12αφ λµ0   L L 30 0    0 −6αφ λµ0 6αγ µ2 x − 2µ (αγL−6λαφ )  L   0 T 0 0    (1−αγ )µ0 3µ0 (1−αγ )(L−2x)  0  − 2 L2   0 0 (1 − α )λµ0 L − 6λµ (1−αφ )(L−2x)    φ L (21) HB(x) =   0  0 0   (1−α )µ0 3µ0 (1−αγ )(L−2x)  0 − 2γ   2 L  6λµ0 (1−αφ )(L−2x)  0 0 −(1 − αφ )λµ L  L   and HB(x) are the linear and hysteretic strain matrices, respectively.

223

Similar to the discussion about the displacement interpolation fields, the strain

224

fields are also generally applicable. For example, the second row elements of

225

matrices B(x) and HB(x) in Eqs. (20) and (21), representing shear strain, become

226

zero for an Euler-Bernoulli beam case (αγ = 1 and λ = 0).

227

Substituting the strain and hysteretic deformation fields expressions from

228

Eqs. (18) and (19) to Eq. (4), the internal force vector of the beam element

15

229

can be reformulated as: P(x) = De ε(x) + HD z(x) = De (B(x) d + HB(x) z) + HD (Nz(x) z) ∴ P(x) = (De B(x) )d + (De HB(x) + HD Nz(x) )z

(22)

230

3.3. Derivation of element matrices

231

The stiffness and hysteretic matrices of the beam element are formulated in this

232

section by means of the virtual work principle [31], which can be written as: ∫ L
T δd F = N(x) δεu(x) + Q(x) δεγ(x) + M(x) δεφ(x) dx (23) 0

233

where δd is the virtual nodal displacement vector, F is external nodal load, and

234

(δεu(x), δεγ(x), δεφ(x) ) are the relevant internal virtual generalized strains at any

235

distance x. By substituting the kinematic conditions of Eq. (9) to Eq. (23), the

236

237

following relations are obtained:     ∫ L ∂w ∂u T + Q(x) δ −θ δd F = N(x) δ ∂x ∂x 0   ∂θ + M(x) δ dx ∂x     ∫ L ∂(δw) ∂(δu) T δd F = + Q(x) N(x) ∂x ∂x 0    ∂(δθ) − Q(x) δθ dx + M(x) ∂x Integration by parts, results in:    ∫ L L ∂N(x) T δd F = N(x) (δu) 0 − δu dx ∂x 0    ∫ L L ∂Q(x) + Q(x) (δw) 0 − δw dx ∂x 0    ∫ L ∫ L L ∂ M(x)
+ M(x) (δθ) 0 − δθ dx − Q(x) δθ dx ∂x 0 0 16

(24)

(25)

238

By imposing the integration limits at nodes 1 and 2, and considering Eqs. (10)

239

to (12), Eq. (25) becomes: δdT F = −δu1 N1 − δw1 Q1 − δθ 1 M1 + δu2 N2 + δw2 Q2 + δθ 2 M2

(26)

240

Hence, equilibrium is satisfied in a strong sense, such that the internal nodal forces

241

are exactly equal to the external ones, since the present formulation is based on

242

the exact equilibrium equations of the nonlinear beam element. Eq. (26) can be

243

further expressed as: δd F = δd T

T

    −P1       P2    

where

P1 = {N1 Q1 M1 }T P2 = {N2 Q2 M2 }T

(27)

244

Substituting the internal force vector expression obtained in Eq. (22) to Eq. (27)

245

and considering the nodal boundary conditions, results in:       −P1  

  −De B(x=0)  d+ =  F=     De B(x=L)   P2      = Kd + Hz

  −De HB (x=0) − HD Nz(x=0)   z    De HB (x=L) + HD Nz(x=L)   

(28)

246

Eq. (28) represents the final concise force-displacement expression in local co-

247

ordinates for the hysteretic Timoshenko beam formulation. Consistent element

248

stiffness matrix K and hysteretic matrix H are obtained by substituting in Eq. (28)

249

the rigidity matrices De and HD from Eq. (5), the strain matrices B and HB from

250

Eqs. (20) and (21), and the hysteretic interpolation matrix Nz from Eq. (18),

251

resulting in:

17

 αu AE  0 0 − αuLAE 0 0   L   6ψ1 6ψ1  12ψ1 1  0 0 − 12ψ L3 L2 L3 L2    0 6ψ1 4ψ2 2ψ3  1 0 − 6ψ  2 L L  L2 K =  α AE L  αu AE − uL 0 0 0 0  L    0 − 12ψ1 − 6ψ1 0 12ψ1 6ψ1  −  3 2 3 2 L L L L    2ψ 6ψ 4ψ 6ψ 3 1 2  1  0 0 − L2 2 L L L   0 where ψ1 = αφ αγ µ E I ;

(29)

ψ2 = αφ (αγ + 3αφ λ)µ0 E I ; ψ3 = αφ (αγ − 6αφ λ)µ0 E I −h 0 0  u  h φ αγ 6hγ  0 − L2 − L   0 − 3hγ −h (α + 6α λ) φ γ φ  L H=   hu 0 0  h φ αγ  0 6hγ  L L2   0 − 3hLγ 6hφ αφ λ  where hu = (1−α2u )AE ;

−hu

0 6hγ L2 3hγ − L

0 −

0

hu 0 0

0

6hγ L2 3h − Lγ

   h φ αγ  L   −6hφ αφ λ    0  h φ αγ  − L   hφ (αγ + 6αφ λ)  0

(30)

hφ = (1 − αφ )µ0 E I ; hγ = (1 − αγ )µ0 E Iαφ

252

3.4. Generality of the beam element model

253

In relation to the previous discussion about the generality of the suggested interpo-

254

lation functions, the derived stiffness and hysteretic element matrices in Eqs. (29)

255

and (30) are applicable for both Euler-Bernoulli and Timoshenko beam element

256

formulations, and for both linear and nonlinear behaviors. For a pure elastic case

257

(αu = αγ = αφ = 1) for example, the hysteretic matrix H in Eq. (30) becomes

258

0, similar to what was mentioned earlier for the hysteretic interpolation functions,

259

and the elastic stiffness matrix K in Eq. (29) reverts to the typical stiffness matrix

260

of an elastic Timoshenko beam [35]. In addition, for λ = 0 and αγ = 1, the beam

261

matrices take the proper form for a nonlinear Euler-Bernoulli beam element [5], 18

262

and by combining the two cases, i.e. for λ = 0 and αu = αγ = αφ = 1, the conven-

263

tional stiffness matrix of a linear Euler-Bernoulli beam element is obtained [31].

264

The resulting stiffness and hysteretic matrices for each of these common cases are

265

shown in the Appendix. Evidently, other cases, combining linear and nonlinear

266

components, are immediately and readily possible with this parametrized form,

267

by merely changing the parameter values, and without altering anything in the

268

modeling of the problem or its solution process.

269

3.5. Numerical solution process

270

Eq. (28) can be interpreted as the conventional stiffness based formulation, aside

271

from the additional hysteretic DOF vector, z. Matrices K and H, and the vectors

272

d and F in Eq. (28) correspond to local element coordinates. Applying coordinate

273

transformation, the global force vector, Fg , can be obtained as: Fg = ΛT KΛdg + ΛT Hz

(Since d = Λdg ; Fg = ΛT F)

(31)

274

where Λ is the typical transformation matrix from global to local coordinate

275

system, and dg corresponds to the element nodal displacement DOF vector in

276

global coordinates [31]. Eq. (31) is thus expressed in terms of global matrices as: Fg = Kg dg + Hg z ; Kg = ΛT KΛ ;

Hg = ΛT H

(32)

277

where Kg and Hg are the element stiffness and hysteretic matrices in global

278

coordinates, with the former comprising the symmetric and positive definite system

279

stiffness matrix, Ks , by the typical direct stiffness matrix assembly techniques [31],

280

and the latter forming the system hysteretic matrix, Hs , by consistently appending

281

the Hg matrices for all the elements. The overall system force-displacement

282

expression can be then expressed as: 19

283

284

285

     z(el=1)    . . Fs = Ks ds + Hs zs ; zs = (33) .   z   (el=Nel )  where Fs and ds are the system nodal force and displacement DOF vectors respectively, z(el= j) is the hysteretic DOF vector of size 6 × 1 for the j th element

( j = {1, 2, · · · , Nel }), where Nel is the total number of elements in the system, and

286

zs is a vector consisting of the hysteretic DOF for all the elements. For example, in

287

a 2-dimensional problem of a cantilever beam, subjected to a concentrated external

288

force at its free end and discretized into 2-elements (Nel = 2), the size of known

289

external force vector, Fs , is (6 × 1), while the unknown displacement, ds , and

290

hysteretic DOF, zs , vectors are of size (6 × 1) and (12 × 1) respectively. Based

291

on Eq. (6), the evolution equation for the j th beam element, of length L j , can be

292

expressed as: zÛ (el= j) = Ω(el= j) εÛ (el= j) where, Ω(el= j)

  (I − H1 H2 R)(x=0)  0  =    (I − H1 H2 R)(x=L j )  0   (el= j)

(34)

z(el= j) = {zu1 zγ1 zφ1 zu2 zγ2 zφ1 }T(el= j)

ε (el= j) = {εu1 εγ1 εφ1 εu2 εγ2 εφ2 }T(el= j)

293

and 0 is a square null matrix of size 3, function H1 and matrix R are shown in

294

Eq. (6). For numerical implementation purposes in quasi-static problems, function

295

H2 is now given as H2 ≈ β + γsgn((Ph )T (ε − ε0 )), which is slightly different from

296

Eq. (6), where ε0 is the initial strain at the beginning of each solution step. Eqs. (33)

297

and (34) have to be solved simultaneously to obtain the solutions, which results in 20

298

a system of DAEs.

299

The formulated DAEs can be solved using a Newton solution scheme, as

300

explained in more detail subsequently. However, in this paper a new, alternative

301

and more efficient approach is suggested for solving the resulting DAEs. The main

302

idea is to substitute the displacement expression from Eq. (33) to the evolution

303

equation of Eq. (34), so that the resulting expressions will only and simply be an

304

explicit first order ODE system with respect to the vector zs . Therefore, Eq. (33)

305

is reformulated as: −1 ds = K−1 s Fs − Ks Hs zs

(35)

306

and based on Eq. (19) the vector of strains for all the elements of the system, εs ,

307

referred to herein as system strain vector, is expressed in terms of the system DOF

308

as follows:    (Bd + HB z)(el=1)      ε(el=1)         . .. .. εs = = .    ε      (el=Nel )   (Bd + HB z)(el=Nel )   (HB Lz )(el=1) zs       (BΛLd )(el=1) ds +  . .. =     (BΛL ) d (el=Nel ) ds + (HB Lz )(el=Nel ) zs  

(36)

= Bs ds + HBs zs

309

where Ld and Lz are the incidence (or connectivity) matrices for transforming

310

dg and z vectors, for each element, to ds and zs vectors respectively, resulting in

311

d = Λdg = ΛLd ds and z = Lz zs . The elements in the incidence matrices are 1,

312

corresponding to the same DOF in both row and column, and 0 otherwise [31].

313

Substituting now Eq. (35) to Eq. (36), the system strain vector is obtained as: −1 εs = (Bs K−1 s )Fs + (HBs − Bs Ks Hs )zs

21

(37)

314

where all the resulting matrices in the parentheses are constant. As seen in Eq. (37),

315

the strain is now a function of the unknown zs vector and therefore H2 , in Eq. (34),

316

also becomes a function of hysteretic DOF, similar to H1 and R. Consequently,

317

the whole matrix Ω for each element is now a function of hysteretic DOF. Based

318

on Eqs. (34) and (36), the overall evolution equation for the system is expressed as

319

follows:

320

 Ω(el=1)  Û   ε(el=1)      . .   .. . zÛ s =   .    Ω(el=Nel )     εÛ (el=Nel )  −1 Û = Ωs εÛ s = Ωs (Bs K−1 zs s )Fs + Ωs (HBs − Bs Ks Hs )Û

(38)

Further simplification of Eq. (38) results in the following ODEs:

−1 zÛ s = (I − Ωs A)−1 Ωs EFÛ s ; A = HBs − Bs K−1 s Hs ; E = Bs Ks

(39)

321

where A and E are constant system matrices and Ωs is a function of zs , as already

322

explained. The ODEs in Eq. (39) are now the only equations needed to solve

323

the problem, having also a reduced number of unknowns, since the displacement

324

DOF have been eliminated from the evolution equations. Eq. (39) can thus be

325

straightforwardly solved now to obtain the hysteretic DOF, using any customary

326

numerical ODE solver, e.g. based on Runge-Kutta methods. As a post-processing

327

step, the resulting zs , at each or from all time steps, can be substituted to Eq. (35),

328

where Ks and Hs are constant, to evaluate the nodal displacement DOF.

329

Alternatively, the suggested DAEs in Eqs. (33) and (34) can be solved using

330

a Newton scheme, by approximating instead the ODE in Eq. (34) in an algebraic

331

form, resulting in the following expressions at the i th step, based also on Eq. (19):

332

Fs (i) = Ks ds (i) + Hs zs (i) z(el= j) (i) = z(el= j) (i − 1) + Ω(el= j) ∆ε(el= j) ; 22

(40) (41)

333

where ∆ε(el= j) = (B(d(i) − d(i − 1)) + HB (z(i) − z(i − 1)))(el= j)

(42)

334

Therefore, by solving Eqs. (40) and (41) by a typical iterative Newton scheme,

335

the unknown nodal displacements, ds , and hysteretic DOF, zs , at step (i) can be

336

obtained.

337

In the present work, both Newton and Runge-Kutta approaches have been

338

implemented, with both methods yielding the same results. The suggested ODEs

339

approach of Eq. (39) using a Runge-Kutta method has been however significantly

340

more computationally efficient in all the examined problems. Finally, for dynamic

341

cases, as also mentioned earlier, the whole formulation can be favorably and

342

straightforwardly described in a state-space ODE form, as in [6], and be again

343

efficiently solved by any generic ODE solver.

344

4. Numerical studies

345

The generality, consistency, versatility and efficiency of the proposed formulation is

346

demonstrated through numerical examples of a cantilever beam model subjected

347

to a concentrated external force at its free end, as shown in Fig. 2. For all

348

the examples, the beam is a W30 × 148 steel beam with yield strength of 380

349

MPa, elastic modulus of 200 GPa and shear modulus of 77 GPa. Based on the

350

cross section properties and yield strength, the moment and shear capacities are

351

calculated to be 3,100 kNm (Mp ) and 2,800 kN (Q p ), respectively. In addition,

352

other considered model parameters are αγ = αφ = 0.01, β = γ = 0.5, and n = 2,

353

for the Timoshenko case. The Euler-Bernoulli formulation is simply achieved by

354

setting the model parameters λ and αγ to 0 and 1, respectively, without any other

355

changes. The formulation is implemented in MATLAB and the solution schemes 23

described in Section 3.5 are employed. L x

W30x148

F

Fig. 2. Cantilever beam subjected to tip force. 356

357

4.1. Shear locking elimination

358

To assess the performance of the suggested formulation in relation to shear locking,

359

the responses of the cantilever beam over a broad range of scenarios are obtained

360

and compared. Since the shear deformation of a long slender beam is negligible

361

as compared to its flexural deformation, the solution obtained from a Timoshenko

362

formulation should converge to that obtained from an Euler-Bernoulli formula-

363

tion. However, this is not often achieved, mostly due to inconsistent interpolation

364

functions, providing reasonable results for short deep beams but overly stiff results

365

for long slender beams, a phenomenon known as shear locking. The suggested

366

formulation provides consistent interpolation functions and results naturally con-

367

verge toward the expected behavior of an Euler-Bernoulli formulation, even with

368

one element and without the need to resort to any approximate methods.

369

To demonstrate this point, the response of the cantilever beam shown in Fig. 2

370

is compared using both Timoshenko and Euler-Bernoulli formulations for different

371

length-to-depth (L/d) ratios, achieved by varying the length of the cantilever beam,

372

ranging from 0.4 m to 10.0 m. For all scenarios, the tip displacement response is

373

obtained when the beam is subjected to a static force of 100 kN at the free end. The 24

374

results are shown in Fig. 3, where y-axis shows the normalized tip displacement

375

with respect to the corresponding Euler-Bernoulli deflections, and x-axis reports

376

the L/d beam ratio. As seen in the figure, for small L/d ratios, the displacements

377

obtained from the Timoshenko formulation are substantially larger than the Euler-

378

Bernoulli case, suggesting significant shear deformations. However, for large L/d

379

values, the Timoshenko response gradually converges to the Euler-Bernoulli one. 12 Euler Beam Timoshenko Beam

10

w/wEuler

8 6 4 2 0 0

5

10

15

L/d Fig. 3. Difference between Euler-Bernoulli and Timoshenko beam formulation.

380

4.2. Timoshenko vs Euler-Bernoulli response investigation

381

In this subsection, the nonlinear global response is evaluated when the free end

382

of the cantilever beam is subjected to a monotonically increasing force, for three

383

different cases: (a) shear dominant beam, (b) intermediate beam, and (c) flexure

384

dominant beam. For these cases, different beam lengths and external forces

385

are considered, as shown in Table 1. Both Timoshenko and Euler-Bernoulli 25

386

formulations are tested and the results are compared, as shown in Fig. 4, where

387

x-axis represents the beam rotation, obtained by normalizing the displacement

388

of the free end with respect to the beam length, and y-axis reports the fixed end

389

bending moment of the beam.

390

Fig. 4(a) shows the results for the shear dominant beam with length 0.7 m. The

391

beam is subjected to a maximum tip force of 3,000 kN, which exceeds the plastic

392

shear capacity, while the corresponding maximum bending moment of 2,100 kNm

393

is within the elastic range. As expected, the Euler-Bernoulli formulation is not able

394

to capture the global inelastic behavior of the beam element, since the maximum

395

bending moment is below the plastic capacity. In contrast, the nonlinear behavior

396

is evident in the Timoshenko beam formulation, with an accurate transition from

397

the elastic to the inelastic regimes, corresponding to the shear capacity of the

398

section.

399

Fig. 4(b) shows the moment-rotation response for the intermediate beam with

400

length 1.1 m subjected to a maximum force of 3,000 kN, which corresponds

401

to a moment of 3,300 kNm. From these results, it is evident that the beam is

402

yielding both in shear and bending, hence the term intermediate beam. However,

403

the significant shear deformation of the beam correctly results in reduced stiff-

404

ness and increased rotation for the Timoshenko formulation, as compared to the

405

corresponding Euler-Bernoulli case.

406

Lastly, Fig. 4(c) shows the moment-rotation response of a 5 m long beam

407

subjected to a maximum tip force of 700 kN, which is less than the shear capacity

408

of the section but produces a moment of 3,500 kNm, which exceeds the flexural

409

capacity. Due to this beam geometry, there is negligible shear deformation and

410

the Euler-Bernoulli formulation agrees sufficiently well with the more general 26

Table 1. Different beam elements with Q p =2,800kN and Mp =3,100kNm Type

411

Beam length (L)

Max applied force (F)

Max moment (F L)

Shear dominant beam

0.7 m

3,000 kN (>Q p )

2,100 kNm (
Intermediate beam

1.1 m

3,000 kN (>Q p )

3,300 kNm (>M p )

Flexure dominant beam

5.0 m

700 kN (
3,500 kNm (>M p )

Timoshenko one, that proved reliable in all cases. (a) 0.7 m long beam

3000 2500

Moment (kNm)

Moment (kNm)

(b) 1.1 m long beam

Timoshenko formulation Euler formulation

2000 1500

(c) 5 m long beam

3500

3500

3000

3000

2500

2500

Moment (kNm)

3500

2000 1500

2000 1500

1000

1000

1000

500

500

500

0

0 0

0.005

0.01

0.015

0.02

Tip rotation

0 0

0.01

0.02

0.03

Tip rotation

0.04

0

0.04

0.08

0.12

0.16

Tip rotation

Fig. 4. Comparison of Timoshenko and Euler-Bernoulli beam formulations for: (a) shear dominant, (b) intermediate, and (c) flexure dominant cantilever beams.

412

4.3. Element consistency

413

The Timoshenko formulation is here employed for the same shear dominant,

414

intermediate and flexure dominant beams, now discretized into five elements, to

415

analyze the local beam behavior including shear strain and curvature distributions

416

within the beam. All the relevant example values are shown in Table 1, as before.

417

Figs. 5(a-c) give the response of the shear dominant beam. Fig. 5(a) shows

418

the global moment-rotation response, where points A and B represent indicative

419

points in the elastic and plastic domains, respectively. The local responses of the 27

420

beam element are presented in Figs. 5(b) and (c). The curvature for both points is

421

shown in Fig. 5(b) and the shear strain is reported in Fig. 5(c). Figs. 5(b-c) show a

422

slight curvature change but a large shear strain variation between points A and B,

423

verifying the shear dominant behavior at the local level.

424

Similarly, Fig. 5(d) shows the global moment-rotation response of the interme-

425

diate beam, whereas Figs. 5(e) and (f) show the local curvature and shear strain

426

plots, respectively, for points C and D. Significant changes in both curvature and

427

shear strain responses are observed from C to D, as expected, since the bending mo-

428

ment and shear force are both exceeding their yield capacities for the intermediate

429

beam.

430

Lastly, Figs. 5(g-i) correspond to the flexure dominant beam, such that Fig. 5(g)

431

displays the global response and Figs. 5(h-i) the local responses at points E and F.

432

As expected, these figures indicate slight shear strain change and large curvature

433

variation from E to F, correctly suggesting a flexure dominant response.

434

From all the curvature and shear strain plots, namely Figs. 5(b-c), (e-f) and

435

(h-i), model consistency in terms of nodal strain field continuity is observed for

436

both linear and nonlinear regimes. Thus, the suggested formulation is compatible

437

with all theoretical assumptions and produces consistent results at both global and

438

local levels for the entire range of beam behaviors.

439

4.4. Discretization and convergence

440

In this subsection, to further investigate the beam response both at the global and

441

local levels, a convergence study for the Timoshenko formulation is conducted

442

by discretizing the beam into different number of elements. In this illustrative

443

example, only the intermediate beam is considered and the beam is discretized

444

into 1, 5, 10 and 30 elements. Fig. 6(a) shows the global response in terms of 28

4000

0.08

(a-c): 0.7 m long beam response (b)

0.06 B

2000

1000

0 (c)

A B

A B

-0.01

Shear strain

3000

Curvature

Moment (kNm)

(a)

0.04

0.02

-0.02 -0.03 -0.04

A 0

0 0

0.005

0.01

0.015

0.02

0.025

-0.05 1

2

3

4

5

1

2

Nodes

Tip rotation 4000

0.08

(e)

1000

5

(f)

6

C D

-0.01

Shear strain

Curvature

0.06

2000

4

0

C D

D 3000

3

Nodes

(d-f): 1.1 m long beam response

(d)

Moment (kNm)

6

0.04

0.02

-0.02 -0.03 -0.04

C 0

0 0

0.01

0.02

0.03

-0.05 1

2

4000

5

6

1

2

0.06

2000

1000

4

5

6

0 (i)

E F

E F

-0.01

Shear strain

3000

3

Nodes

(h)

F

Curvature

Moment (kNm)

4

(g-i): 5 m long beam response

0.08 (g)

3

Nodes

Tip rotation

0.04

0.02

-0.02 -0.03 -0.04

E 0

0 0

0.02

0.04

Tip rotation

0.06

-0.05 1

2

3

4

Nodes

5

6

1

2

3

4

5

6

Nodes

Fig. 5. Strain continuity at the nodes considering 5 elements for: (a-c) shear dominant, (d-f) intermediate, and (g-i) flexure dominant cantilever beams.

29

445

support moment and beam rotation, and Fig. 6(b) represents the local curvature

446

response of the beam element corresponding to the peak moment of 3,300 kNm.

447

For both the global and local responses, convergence is nearly achieved with

448

10 elements. The local curvature response also indicates that by increasing the

449

number of elements a proper continuous and smooth curvature distribution along

450

the length of the beam element is obtained. As can be seen in Fig. 6(b), the

451

estimated plastic hinge length, i.e. the beam length with curvature values greater

452

than the yield curvature, has converged with 10 elements. For this example, the

453

yield curvature is given by My /E I, with My = 2, 700 kNm. Hence, the developed

454

formulation is able to provide consistent, continuous, converged and smooth strain

455

fields. (a)

3500

(b)

0.05 Nel=1 Nel=5 Nel=10 Nel=30 Yield

3000 0.04

Curvature (1/m)

Moment (kNm)

2500 2000 1500 1000

Nel=1 Nel=5 Nel=10 Nel=30

500

0.03

0.02

0.01

0

0 0

0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04

0

Tip rotation

0.2 L

0.4 L

0.6 L

0.8 L

L

Beam length

Fig. 6. Element discretization and convergence analysis for intermediate beam.

456

4.5. Combined axial-flexural loading

457

The effect of interaction is shown in this subsection by subjecting the beam element

458

to combined axial and transverse forces. As seen in the previous subsection, 30

459

convergence is achieved for 10 elements; therefore, only the converged case of the

460

intermediate beam with the Timoshenko formulation is presented here, with the

461

same monotonically increasing transverse force, as in the previous subsections,

462

and an additional constant axial force, applied at the free end of the beam element.

463

Four cases are considered, based on the magnitude of the applied axial force, i.e.

464

np = 0, 0.2, 0.3, 0.4, where np is the ratio of the applied axial force, N, to the

465

axial capacity of the cross-section, Np = 10, 600 kN. Axial-moment interaction is

466

considered using the Orbison criterion [33] and the results are presented in Fig. 7.

467

Fig. 7(a) presents the global moment-rotation response, whereas Fig. 7(b) shows

468

the local moment-curvature behavior at the fixed end. A reduction in only flexural

469

capacity is observed in these figures with the increasing magnitude of axial force

470

due to the axial-moment interaction, while the shear capacity remains unchanged for all the cases. (a) Global response

3000

3000

2500

2500

2000 1500

np= 0 1000

2000 1500 1000

np= 0.2 np= 0.3

500

(b) Local response

3500

Moment (kNm)

Moment (kNm)

3500

500

np= 0.4

0

0 0

0.01

0.02

0.03

0.04

0.05

0.06

0

Tip rotation

0.05

0.1

0.15

0.2

Curvature (1/m)

Fig. 7. Effect of axial-moment interaction for intermediate beam with Timoshenko formulation. 471

472

To further discuss and connect all the presented analyses in this section to31

473

gether, the L/d ratios of the analyzed shear dominant, intermediate and flexure

474

dominant beams in Fig. 3 are obtained as 0.9, 1.4 and 6.4, respectively. As ex-

475

pected, Fig. 3 shows a large variation between Euler-Bernoulli and Timoshenko

476

formulation responses for the L/d ratio of 0.9, deviation decreases for the L/d

477

ratio of 1.4, whereas identical responses are observed for the 6.4 L/d ratio. In

478

Fig. 4, only these three representative beam cases are considered with monotoni-

479

cally increasing load, in contrast to Fig. 3, to compare the overall global response

480

of Timoshenko and Euler-Bernoulli formulations for both elastic and inelastic

481

regimes. The global behavior of Figs. 4(a-c) correspond to Fig. 5(a), Fig. 5(d)

482

and Fig. 5(g), respectively, and due to different element discretizations, global

483

response variations are observed between the two figures. Note that for the shear

484

dominant beam, the global responses in Figs. 4 and 5 are alike, the deviation is

485

increasing for the intermediate beam, and a large variation is observed for the

486

flexure dominant beam. This is due to the fact that the constant shear strain field in

487

this case is equally captured by either one or five elements, however, as the bending

488

contribution increases, the curvature response is more accurately computed with

489

larger number of elements. To further analyze this, in Fig. 6, the intermediate beam

490

is discretized with different number of elements and the responses are compared

491

for both the global moment-rotation and the local curvature responses. Lastly, to

492

emphasize the effect of interaction, the converged case of Fig. 6 is employed and

493

presented in Fig. 7, with combined axial-flexural loading.

494

5. Experimental validation

495

To validate the suggested formulation, the simulated responses are compared

496

with the reported experimental beam test results from Berman and Bruneau [36]. 32

497

Berman and Bruneau [36] carried out extensive testing on twelve beam specimens,

498

consisting of four different lengths with three cross sections, observing varied

499

responses, ranging from shear to flexure dominant behaviors. In the present study,

500

two such beams are simulated for validation purposes: a shear dominant beam

501

with the shortest length, where the total beam rotation is largely governed by the

502

shear strain, and a flexure dominant beam with the longest length, such that the

503

beam deformation is predominantly due to flexure. The two ends of the beams are

504

restrained in the experimental setup, except for the transverse displacement at one

505

end, where the beam is subjected to a typical cyclic loading protocol with increasing

506

amplitude. To account for some observed flexibility at the end restraints of the

507

beams, linear rotational springs are considered at the beam ends in the analysis.

508

Further details about the experimental setup, cross sectional geometry, loading

509

protocol, etc. can be seen in [36].

510

5.1. Shear dominant beam

511

The shear dominant beam response is simulated using both the suggested Timo-

512

shenko and Euler-Bernoulli formulations. Specimen X3L1.2 in [36] is considered

513

with a length of 578 mm, web yield strength of 430 MPa, flange yield strength

514

of 370 MPa, elastic modulus of 200 GPa and shear modulus of 77 GPa. From

515

the beam mechanical properties, the moment and shear capacities are specified to

516

be as 310 kNm (Mp ) and 1,000 kN (Q p ), respectively. Other considered model

517

parameters, without any inverse method utilization, are αφ = 0.01, αγ = 0.01,

518

β = γ = 0.5, and n = 1, for the Timoshenko case. Fig. 8(a) presents the global

519

force-rotation response, where rotation is the ratio of transverse displacement at

520

the free end to the effective length of the member. Figs. 8(b-c) show the local

521

moment-curvature and shear force-shear strain responses at the fixed end, respec33

522

tively. Experimental curvature is calculated by taking the ratio of the difference

523

between the top flange and bottom flange strains over the total depth of the section,

524

while experimental shear strains are obtained from the strain gauges installed at

525

the mid-section and spanning at 45 degrees from the longitudinal axis. Similarly,

526

Figs. 8(d-f) present the Euler-Bernoulli response compared to the same experi-

527

mental data, where Fig. 8(d) shows the global behavior and Figs. 8(e-f) the local

528

cross section responses.

529

As seen in Figs. 8(a-c), the response of the beam is well represented by the

530

Timoshenko formulation. In contrast, using the same parameters, the obtained

531

responses from the Euler-Bernoulli formulation show significant differences with

532

the experimental data in Figs. 8(d-f). It is however interesting to note that the global

533

behavior in Fig. 8(d) approximates the response to some level of accuracy, but due

534

to the limitation of the Euler-Bernoulli formulation to capture shear deformations,

535

large curvature deviation is obtained, as compared to the experimental values.

536

This example, hence, demonstrates the suggested Timoshenko formulation

537

ability to effectively simulate the beam behavior both at the global and local levels

538

with increased accuracy. In the next validation, the same process is repeated for

539

the flexure dominant beam to demonstrate the overall model robustness.

540

5.2. Flexure dominant beam

541

In this example, the X1L3.0 flexure dominant beam specimen in [36] is simulated,

542

with a length of 1,441 mm, web yield strength of 480 MPa, flange yield strength

543

of 371 MPa, elastic modulus of 200 GPa and shear modulus of 77 GPa. Similar

544

to the previous example, model parameters for bending and shear capacities are

545

calculated as Mp = 302 kNm and Q p = 936 kN, respectively. Other used param-

546

eters, without any inverse method utilization again, are αφ = 0.02, αγ = 0.02, 34

Experiment

500

Moment (kNm)

Force (kN)

750

0

-750

-1500 -0.1

250

0

-250

0

0.05

0.1

-0.2

(d)Euler:Global

-750

250

0

-500 0.05

-750

-0.05

0.1

Total rotation (rad)

-0.7

0

0

0.05

0.1

(f)Euler:Local

1500

-250

0

0

Shear strain

Shear force (kN)

0

-0.05

750

-1500 -0.1

0.2

(e)Euler:Local

500

Moment (kNm)

Force (kN)

0

(c)Timoshenko:Local

Curvature (1/m)

750

-1500 -0.1

1500

-500 -0.05

Total rotation (rad)

1500

Hysteretic finite element

(b)Timoshenko:Local

Shear force (kN)

1500

(a)Timoshenko:Global

0.7

750

0

-750

-1500 -0.1

-0.05

Curvature (1/m)

0

0.05

0.1

Shear strain

Fig. 8. Comparison of experimental and simulated response for shear dominant beam.

547

β = γ = 0.5, n = 1, for the Timoshenko formulation. For the Euler-Bernoulli

548

formulation all the parameters again remain the same, except for λ and αγ , which

549

are set to 0 and 1, respectively.

550

In Fig. 9, the Euler-Bernoulli and Timoshenko formulation responses agree

551

well, both at the global and local levels. It is also interesting to note that the

552

response obtained from the suggested Timoshenko formulation is able to capture

553

the small linear shear strain in the element without any shear locking effects. In

554

this scenario, the Timoshenko formulation is again capable to capture the response,

555

however, the Euler-Bernoulli formulation is also sufficient here.

35

Experiment

(c)Timoshenko:Local

600

600

300

300

300

0

-600 -0.06

Shear force (kN)

600

-300

0

-300

-600 -0.03

0

0.03

0.06

-0.5

Total rotation (rad)

0

-300

-600 -0.25

0

0.25

0.5

-0.01

(e)Euler:Local

300

300

-600

Shear force (kN)

300

Moment (kNm)

600

-300

0

-300

-600 -0.03

0

0.03

Total rotation (rad)

0.06

-0.5

0.005

0.01

(f)Euler:Local

600

0

0

Shear strain

600

-0.06

-0.005

Curvature (1/m)

(d)Euler:Global

Force (kN)

Hysteretic finite element

(b)Timoshenko:Local

Moment (kNm)

Force (kN)

(a)Timoshenko:Global

0

-300

-600 -0.25

0

0.25

Curvature (1/m)

0.5

-0.01

-0.005

0

0.005

0.01

Shear strain

Fig. 9. Comparison of experimental and simulated response for flexure dominant beam.

556

6. Conclusions

557

In this paper, a consistent and computationally efficient hysteretic beam finite ele-

558

ment model is presented. The hysteretic modeling component, expressed in terms

559

of stress resultants and featuring distributed multiaxial plasticity with interactions,

560

is fully integrated into the element formulation, enabling numerous important

561

features. Nodal displacement and strain fields continuity along the element is

562

achieved, due to development of new displacement field interpolation functions.

563

These functions take into account the introduced hysteretic degrees of freedom

564

and satisfy both the exact equilibrium equations and kinematic conditions, leading

565

to a consistent element formulation. Although it is shown in this work that the 36

566

suggested Timoshenko beam element can seamlessly model all possible scenarios,

567

the parametrized nature of the element allows also a direct transformation to an

568

Euler-Bernoulli beam element, by simple parametric changes. Additionally, the

569

suggested formulation produces constant element stiffness and hysteretic matrices

570

that do not need to be updated throughout the analysis, given that inelasticity is

571

treated through localized hysteretic evolution equations. Finally, a new compu-

572

tationally efficient solution scheme is suggested for quasi-static problems, where

573

the entire framework is presented as an ODE system and can be efficiently solved

574

without the need for any type of linearization. Several numerical examples are

575

included in this work in order to showcase the characteristics and capabilities of

576

the element, and comparisons with available experimental data are also provided.

577

Overall, the presented element possesses significant, favorable properties and can

578

be easily modeled and even incorporated in existing finite element codes. Several

579

further extensions are also available and future work will include, among others,

580

incorporation of large displacements and damage-plasticity attributes, concertedly

581

preserving its computational efficiency.

582

7. Acknowledgements

583

The authors would like to acknowledge the support of the U.S. National Science

584

Foundation, which partially supported this research under Grant No. CMMI-

585

1351591, and to thank Dr. Jeffrey W. Berman and Dr. Michel Bruneau for

586

providing the detailed experimental data.

37

587

Appendix. Stiffness and hysteretic matrices for common cases

588

CASE 1. Nonlinear Euler-Bernoulli beam element (λ = 0 and αγ = 1):

589

590

591

 αu AE  0 0 − αuLAE 0 0  L     0 12αLφ3E I 6αLφ 2E I 0 − 12αLφ3E I 6αLφ 2E I    6α E I 2αφ E I  4αφ E I 6αφ E I  0 0 − φ2 2 L L L L  K =  α AE  αu AE 0 0 0 0 − uL  L  12αφ E I 6αφ E I 12αφ E I 6αφ E I   0 − L3 − L2  0 − L3 L2   6αφ E I 4αφ E I  2αφ E I 6αφ E I  0 0 − L L L2 L2  

− (1−αu ) AE 0  0 − (1−αu2 ) AE 0 0 2   (1−αφ )E I (1−αφ )E I   0 0 − 0 0 L L    0  0 −(1 − αφ )E I 0 0 0   H =  (1−αu ) AE (1−αu ) AE  0 0 0 0 2 2   (1−α )E I (1−α )E I φ φ  0  0 0 − 0 L L    0 0 0 0 0 (1 − αφ )E I    CASE 2. Linear Timoshenko beam element (αu = αγ = αφ = 1):  AE 0  L   0 12µLE3 I   0 6µ E I K =  AE L 2 − L 0   0 − 12µLE3 I   0 6µ E2 I L 

0

6µ E I L2 4(1+3λ)µ E I L

0 0

0

12µ E I L3 6µ E I − 2 L

−

0

AE L

0

6µ E I L2 2(1−6λ)µ E I L

0

12µ E I L3 6µ E I − 2 L

−

H = 0 ; where, µ =

592

− AE L

0

   6µ E I  L2  2(1−6λ)µ E I  L   0   6µ E I  − 2 L  4(1+3λ)µ E I  L  0

1 EI ;λ = 1 + 12λ G As L 2

CASE 3. Linear Euler-Bernoulli beam element (λ = 0 and αu = αγ = αφ = 1):  AE 0 0 − AE 0 0  L  L    I 6E I 0 − 12E3 I 6E2I   0 12E 3 2 L L L L    0 6E2I 4EL I 0 − 6E2I 2EL I  L L   ;H = 0 K =  AE  AE 0 0 0 0 −  L  L   I 6E I 12E I 6E I   0 − 12E − 0 − L3 L2 L3 L2    0 6E2I 2EL I 0 − 6E2I 4EL I  L L   38

593

References

594

[1] G. G. Deierlein, A. M. Reinhorn, M. R. Willford, Nonlinear structural analysis

595

for seismic design, NEHRP Seismic Design Technical Brief No. 4 (2010) 1–

596

36.

597

[2] D. G. Lignos, H. Krawinkler, Deterioration modeling of steel components

598

in support of collapse prediction of steel moment frames under earthquake

599

loading, Journal of Structural Engineering 137 (11) (2010) 1291–1302.

600

[3] E. Spacone, F. C. Filippou, F. F. Taucer, Fibre beam–column model for non-

601

linear analysis of R/C frames: Part I. Formulation, Earthquake Engineering

602

& Structural Dynamics 25 (7) (1996) 711–725.

603

[4] M. H. Scott, G. L. Fenves, F. McKenna, F. C. Filippou, Software patterns for

604

nonlinear beam-column models, Journal of Structural Engineering 134 (4)

605

(2008) 562–571.

606

[5] S. P. Triantafyllou, V. K. Koumousis, Small and large displacement dynamic

607

analysis of frame structures based on hysteretic beam elements, Journal of

608

Engineering Mechanics 138 (1) (2011) 36–49.

609

[6] S. P. Triantafyllou, V. K. Koumousis, An inelastic Timoshenko beam ele-

610

ment with axial–shear–flexural interaction, Computational Mechanics 48 (6)

611

(2011) 713–727.

612

613

[7] J. W. Macki, P. Nistri, P. Zecca, Mathematical models for hysteresis, SIAM Review 35 (1) (1993) 94–123.

39

614

615

616

617

618

619

[8] A. Visintin, Mathematical models of hysteresis, in: The Science of Hysteresis, vol. 1, Elsevier, 2005. [9] M. Dimian, P. Andrei, Noise-driven phenomena in hysteretic systems, vol. 218, Springer, 2014. [10] W. D. Iwan, A distributed-element model for hysteresis and its steady-state dynamic response, Journal of Applied Mechanics 33 (4) (1966) 893–900.

620

[11] L. F. Ibarra, R. A. Medina, H. Krawinkler, Hysteretic models that incorporate

621

strength and stiffness deterioration, Earthquake Engineering & Structural

622

Dynamics 34 (12) (2005) 1489–1511.

623

[12] S. Erlicher, Hysteretic degrading models for the low-cycle fatigue behaviour

624

of structural elements: Theory, numerical aspects and applications, Ph.D.

625

thesis, University of Trento, 2003.

626

627

628

629

630

631

[13] R. Bouc, Forced vibrations of mechanical systems with hysteresis, in: Proc. of the 4th Conference on Nonlinear Oscillations, Prague, 1967. [14] Y.-K. Wen, Method for random vibration of hysteretic systems, Journal of the Engineering Mechanics Division 102 (2) (1976) 249–263. [15] F. Casciati, Stochastic dynamics of hysteretic media, Structural Safety 6 (2-4) (1989) 259–269.

632

[16] M. Sivaselvan, A. Reinhorn, Nonlinear structural analysis towards collapse

633

simulation: A dynamical systems approach, Technical Report MCEER-04-

634

0005, 2004.

40

635

[17] S. P. Triantafyllou, Hysteretic finite elements and macro-elements for non-

636

linear dynamic analysis of structures, Ph.D. thesis, National Technical Uni-

637

versity of Athens, 2011.

638

[18] J. Lubliner, Plasticity theory, Dover Publications, 2008.

639

[19] T. T. Baber, Y.-K. Wen, Random vibration hysteretic, degrading systems,

640

641

642

643

644

645

646

Journal of the Engineering Mechanics Division 107 (6) (1981) 1069–1087. [20] T. T. Baber, M. N. Noori, Random vibration of degrading, pinching systems, Journal of Engineering Mechanics 111 (8) (1985) 1010–1026. [21] G. C. Foliente, Hysteresis modeling of wood joints and structural systems, Journal of Structural Engineering 121 (6) (1995) 1013–1022. [22] M. V. Sivaselvan, A. M. Reinhorn, Hysteretic models for deteriorating inelastic structures, Journal of Engineering Mechanics 126 (6) (2000) 633–640.

647

[23] J. Song, A. Der Kiureghian, Generalized Bouc–Wen model for highly asym-

648

metric hysteresis, Journal of Engineering Mechanics 132 (6) (2006) 610–618.

649

[24] K. G. Papakonstantinou, P. C. Dimizas, V. K. Koumousis, An inelastic beam

650

element with hysteretic damping, Shock and Vibration 15 (3, 4) (2008) 273–

651

290.

652

[25] P. S. Harvey Jr, H. P. Gavin, Truly isotropic biaxial hysteresis with arbitrary

653

knee sharpness, Earthquake Engineering & Structural Dynamics 43 (13)

654

(2014) 2051–2057.

655

[26] M. Ismail, F. Ikhouane, J. Rodellar, The hysteresis Bouc-Wen model, a survey,

656

Archives of Computational Methods in Engineering 16 (2) (2009) 161–188. 41

657

[27] M. Amabili, Derivation of nonlinear damping from viscoelasticity in case of

658

nonlinear vibrations, Nonlinear Dynamics (2018) doi:10.1007/s11071–018–

659

4312–0.

660

[28] M. Amabili, Nonlinear damping in nonlinear vibrations of rectangular plates:

661

Derivation from viscoelasticity and experimental validation, Journal of the

662

Mechanics and Physics of Solids 118 (2018) 275–292.

663

664

665

666

[29] M. Amabili, Nonlinear damping in large-amplitude vibrations: modelling and experiments, Nonlinear Dynamics 93 (2018) 5–18. [30] J. Rakowski, The interpretation of the shear locking in beam elements, Computers & Structures 37 (5) (1990) 769–776.

667

[31] K. J. Bathe, Finite element procedures, Prentice Hall, 1996.

668

[32] S. Erlicher, N. Point, Thermodynamic admissibility of Bouc–Wen type hys-

669

teresis models, Comptes Rendus Mecanique 332 (1) (2004) 51–57.

670

[33] J. G. Orbison, W. McGuire, J. F. Abel, Yield surface applications in non-

671

linear steel frame analysis, Computer Methods in Applied Mechanics and

672

Engineering 33 (1-3) (1982) 557–573.

673

[34] A. Charalampakis, V. Koumousis, Ultimate strength analysis of composite

674

sections under biaxial bending and axial load, Advances in Engineering

675

Software 39 (11) (2008) 923–936.

676

677

[35] Z. Friedman, J. B. Kosmatka, An improved two-node Timoshenko beam finite element, Computers & Structures 47 (3) (1993) 473–481.

42

678

[36] J. Berman, M. Bruneau, Further development of tubular eccentrically braced

679

frame links for the seismic retrofit of braced steel truss bridge piers, Technical

680

Report MCEER-06-0006, 2006.

43

Highlights:

• • • • •

A general consistent hysteretic beam finite element model is developed New displacement interpolation functions are introduced to satisfy exact equilibrium conditions Strain continuity is achieved without any shear-locking effects A new computationally efficient solution scheme is suggested Axial-moment-shear interactions and distributed plasticity are considered