Influence of shear deformations on the buckling of columns using the Generalized Beam Theory and energy principles

European Journal of Mechanics A/Solids 61 (2017) 216e234

Contents lists available at ScienceDirect

European Journal of Mechanics A/Solids journal homepage: www.elsevier.com/locate/ejmsol

Influence of shear deformations on the buckling of columns using the Generalized Beam Theory and energy principles ~o a, b, * Pedro Dias Sima a b

Institute of Computers and Systems Engineering of Coimbra (INESC-Coimbra), Rua Antero de Quental 199, 3000-033 Coimbra, Portugal Department of Civil Engineering, University of Coimbra, Rua Luís Reis Santos, 3030-788 Coimbra, Portugal

a r t i c l e i n f o

a b s t r a c t

Article history: Received 10 February 2016 Received in revised form 29 September 2016 Accepted 30 September 2016 Available online 3 October 2016

The paper presents a procedure for the stability analysis of columns that are sensitive to shear deformations, the so-called “weak-in-shear” columns that can be found, in the engineering practice, in build-up or composite columns, or in elastomeric bearings. Two distinct formulas are commonly used to compute the critical load for shear sensitive columns: the Engesser and the Haringx formulae, the latter enabling significantly higher loads. They differ on the choice of the cross section's shear stresses resultant, and this duality has been object of much passionate discussion during the last decades. This problem is analysed here under the perspective of the Generalized Beam Theory, and a specific mode for shear deformations was developed using two distinct strategies: i) a linear shear formulation, corresponding to the Timoshenko beam theory for which cross sections remain plane after deformation, and ii) a nonlinear shear formulation, for which shear warping is allowed in order to accomplish with the condition of null shear distortions at the section's contour. A total potential energy is defined, assuming a linear elastic behaviour, and the correspondent functional is rendered discrete by means of the RayleighRitz method. Finally, the traditional stability procedures are applied and the critical loads are computed. The Engesser critical load is derived by applying the stability procedures to the total potential energy associated with the linear shear formulation. A parametric study on the critical behaviour of a shear deformable column and a set of conclusions end the paper. © 2016 Elsevier Masson SAS. All rights reserved.

Keywords: Engesser and Haringx formulae Generalized Beam Theory Nonlinear shear warping displacements

1. Introduction For a long time, the classical stability problem of “weak-inshear” compressed members has been attracting the attention of many researchers worldwide. This problem has significant practical relevance on the buckling of stocky, build-up or composite compressed members with relatively small shear stiffness when compared to the bending's one. In the engineering practice, this is the case, for example, of elastomeric bearings, commonly used in the seismic isolation of buildings (Kelly and Konstantinidis, 2011), and of composite columns with soft shear core (Allen, 1969). Nowadays, these theories are regaining new interest due to their applications in new fields such as the buckling of nanostructures (Wang et al., 2010; Zhang et al., 2006). There are two classical

theories to access the critical load of shear deformable columns. The first theory is due to Engesser (Engesser, 1891) and computes the critical load by:

PEngesser ¼

http://dx.doi.org/10.1016/j.euromechsol.2016.09.015 0997-7538/© 2016 Elsevier Masson SAS. All rights reserved.

L2e

1 $ 2 EI 1 þ L2p$kGA

(1)

e

where Le is the effective buckling length of the column, and A and I are the cross section area and the cross section's moment of inertia, respectively. k is a shear correction factor (Timoshenko, 1921, 1922) that defines the reduced cross section area Ar ¼ k$A, while E and G are respectively the Young and the shear modulus. The second theory is due to Haringx (Haringx, 1949), for which the critical load is equal to:

PHaringx * Department of Civil Engineering, University of Coimbra, Rua Luís Reis Santos, 3030-788 Coimbra, Portugal E-mail address: [email protected].

p2 EI

where

GAr ¼ 2

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! 4P 1 þ Euler  1 GAr

(2)

~o / European Journal of Mechanics A/Solids 61 (2017) 216e234 P.D. Sima

PEuler ¼

p2 EI L2e

(3)

is the Euler buckling load. Although both formulations assume linear elastic constitutive relations between the shear stress resultant VSH and the distortion angle w:

VSH ¼ GAr $w ¼ GAr $w0SH

(4)

and between the bending moment M and the linearized curvature 1: r

1 00 M ¼ EI$ z  EI$wB

r

(5)

Taking into account Fig. 1, wSH is the transverse deflection due only to shear while wB is the transverse displacement due only to bending. The prime 0 denotes differentiation with respect to the longitudinal coordinate x. Formulas (1) and (2) give quite different values for the smaller slenderness coefficients range due to the different choice of the cross section's shear VSH stress resultants. In Appendix A, both formulas are derived, highlighting their different choices for the transverse shear force. Putting apart the analysis of rings and curved bars, that are not of our concern here, and focusing on straight bars only, during the last decades much passionate discussion took place on this subject, some supporting the Engesser approach while others supporting

217

the Haringx strategy. So far, a consensus is not reached yet and, in the following, some works are briefly referred, in chronological order, to illustrate this long debate. Allen (Allen, 1969) studied the buckling behaviour of sandwich €nni, 1971) structures and derived the Engesser formula. N€ anni (Na studied the effect of shear forces on the buckling load of simply supported members with narrow rectangular cross section (plane stress or strain states) and of elliptic cross section (three-dimensional problem), and concluded the Engesser approach was superior. Short time later, Reissner (Reissner, 1972) computed a onedimensional large-strain beam theory for plane beams and also analysed the buckling of circular rings, and noticed that Haringx theory provided better agreement. Ziegler (Ziegler, 1982) formulated a one-dimensional formulation for bars, and concluded that the Engesser approach is superior to Haringx's for columns and can still be improved if axial shortening is taken into account, while Haringx's is more appropriate for helical springs. These conclusions were, shortly after, rebutted by Reissner (Reissner, 1982), who again supported Haringx's formulation, finding that both theories can be obtained by different forms of the one-dimensional stress-strain relations. Simo and Kelly (Simo and Kelly, 1984) performed a twodimensional buckling analysis of multilayer elastomeric bearings, and concluded that, as long as the beam theory assumptions hold, their formulation was in closer agreement with the Haringx theory. Banerjee and Williams (Banerjee and Williams, 1994) analysed shear-deformable uniform columns, derived the Engesser formula and buckling curves for the most common supporting conditions.

Fig. 1. The equilibrium analysis for the shear-deformable column.

218

~o / European Journal of Mechanics A/Solids 61 (2017) 216e234 P.D. Sima

Kardomateas et al (Kardomateas et al., 2002). made theoretical predictions and FEM analyses of buckling loads for sandwich columns with metallic laminated facings and foam or honeycomb cores under uniform axial compression, and noticed that Engesser formula was too conservative and Haringx's was closer to their results. Attard (Attard, 2003) computed an internal strain energy density for isotropic hyperelastic Hookean materials and showed that this formulation conduces to a buckling load formula in closer agreement with Haringx's. This author reinforced these arguments in a subsequent paper (Attard and Hunt, 2008), and later applied both Engesser's and Haringx's formulations for the stability analysis of non-conservative systems (Attard et al., 2008). Ba zant (Ba zant, 2003) performed some theoretical considerations on the buckling behaviour of weak-in shear columns and, despite some correspondence between Engesser and Haringx formulations, concluded that Haringx formula gave better results. In a previous work, this same author derived only the Engesser formula from the theoretical analysis of the problem (Ba zant and Cedolin, 1991). Tsai and Kelly (Tsai and Kelly, 2005) introduced, as far as known for the first time, the shear warping displacements in the buckling analysis of short beams, and noticed their solutions converged to Haringx formula when warping stiffness tends to infinity. Aristizabal-Ochoa (Aristizabal-Ochoa, 2007, 2008) analysed the stability problem of Timoshenko beam-columns under compressive conservative or non-conservative loads with elastic supports, including geometrical imperfections and loading eccentricities, and found that the developed formulations agreed with Haringx formula. These results were rebutted by Blaauwendraad (Blaauwendraad, 2008, 2010), who revisited a wide range of previous works from other authors and presented very strong arguments in favour of the Engesser formulation. More recently, it is worth to mention Challamel et al. (Challamel et al., 2013), where a consistent energy formulation for the stability analysis of shear beam-columns is presented and, among other things, the influence of shear warping on the buckling load is evaluated. Also in this work, the Engesser formula is derived by imposing some simplifications to the energy model. Finally, Forcellini and Kelly (Forcellini and Kelly, 2014) analysed the stability behaviour of elastomeric bearings and adopted the Haringx strategy. The shear mode of deformation is traditionally modelled by means of the Timoshenko kinematic pattern (Timoshenko, 1921, 1922), henceforth denoted by first-order or linear shear theory. This strategy assumes a constant shear distortion along the whole cross section and plane cross sections after deformation. However, it violates the boundary condition of null shear stresses normal to the member's contour, since a constant shear stress is generated along the whole cross section. Hence, the analysis of shear deformable members requires a correction factor k, used to compute the reduced cross section area Ar ¼ k$A, so that the internal strain energy due to shear equals the one arising from the correct shear stress distribution. Cowper (Cowper, 1966) computed consistently the shear correction factor k for several cross section shapes by means of the elasticity theory, while Frank Pai and Schultz (Frank Pai and Schultz, 1999) derived the shear correction factors by matching the exact shear stress resultants with the ones computed by an equivalent first-order theory. Even so, the linear shear theory still constitutes a simplified strategy and is not equivalent to an enhanced beam theory that accounts for the correct shape of shear deformations and stresses along the member. Since the early eighties, several strategies to define the nonlinear shear warping displacements field for prismatic beams and plates were developed. Complete descriptions on these theories can be found in Wang et al (Wang et al., 2000). or Nayfeh and Frank Pai (Naifeh and Frank Pai, 2004). Here, and focusing on beam analysis only, it is worth mentioning the polynomial shapes of Levinson

(Levinson, 1981, 1987), Bickford (Bickford, 1982), Heyliger and Reddy (Heyliger and Reddy, 1988) and Petrolito (Petrolito, 1995), and the sinusoidal shape of Touratier (Touratier, 1991). Dufort et al (Carrera and Giunta, 2010). refer that there is little difference between the polynomial and the sinusoidal shapes, since the cubic functions adopted in the polynomial formulations are close to the power series development of the Touratier's trigonometric function up to order 3. Finally, it is worth to refer the works of Carrera et al. (Carrera and Giunta, 2010; Carrera et al., 2010, 2011), where the displacements fields for generic cross section members are computed through an unified and systematic procedure for prismatic members under bending, shear or torsion. In the present paper and in order to give some contribution to the clarification of the Engesser-Haringx discord, we analyse the buckling behaviour of general two-dimensional shear-deformable prismatic columns under compression by means of the Generalized ~o et al., Beam Theory (GBT). The GBT theory (Schardt, 1989; Sima ~o and Simo ~ es da Silva, 2004; Sima ~o, 2007) was 2002; Sima initially developed by Schardt (Schardt, 1966) and constitutes a whole subject on structural mechanics. It is based on the definition of unitary kinematic deformation patterns along the cross section, denoted by modes of deformation. Thus, the displacements of a generic point of the member is computed by a linear combination of these modes of deformation, where each mode is multiplied by its correspondent amplitude modal function that only depends on the member's longitudinal coordinate x. Recently, this strategy was applied to the buckling and post-buckling analyses of prismatic slender columns with splices. It was validated for the critical col~o umn behaviour, by the derivation of the Euler buckling load (Gira Coelho et al., 2010), and for the initial post-buckling behaviour, by comparison with the Finite Element Method and with theoretical ~o et al., 2012). The GBT strategy is extended here for solutions (Sima the analysis of two-dimensional “weak-in-shear” prismatic members, by the definition of three modes of deformation: to the traditional axial elongation (mode 1) and bending (mode 2) modes, a third mode (mode 3) related to shear deformations is added. For the shear mode of deformation, two distinct kinematic patterns are developed: i) the Timoshenko shear deformation pattern, which corresponds to a linear shear formulation and keeps cross sections plane after deformation, and ii) a nonlinear shear pattern, which consistently models the cross section warping that occurs in prismatic members due to shear and respects the condition of null shear stresses normal to the member's contour. These kinematic deformation patterns allow a consistent calculus of the internal strains, making resource to the Engineering Strains concept (Crisfield, 1991), and the correspondent normal and shear stresses derive from a linear elastic constitutive relation. An internal strain energy is thus computed and, together with the potential of the external work, they define the total potential energy of the system. Having obtained the total potential energy, it is rendered discrete by means of the Rayleigh-Ritz method, and each modal amplitude function is approximated by a linear combination of preestablished coordinate functions. The adopted coordinate functions are polynomials that derive from the relevant modal boundary ~o, 2007; Sima ~o et al., 2003, 2012). Then, the conditions (Sima traditional stability procedures (Thompson and Hunt, 1973; Hunt, 1981) are applied for buckling analysis, giving rise to a nonlinear eigenvalue problem, whose lowest eigenvalue and the correspondent eigenvector describe the critical state. All computations in the present work were performed in the context of the symbolic algebraic manipulator Mathematica (Corp, 2010), and some illustrative examples for the buckling behaviour of columns that can undergo shear deformations are presented, in order to describe the influence of several aspects on the buckling load, such as the modal boundary conditions. The Engesser critical load is derived from the

~o / European Journal of Mechanics A/Solids 61 (2017) 216e234 P.D. Sima

linear shear formulation, and it was also noted that this formula is a lower bound for the critical load and agrees much better with the results arising from linear and nonlinear shear strategies than Haringx's. A set of conclusions ends the paper.

219

displacements pattern 3 uðzÞ for mode 3. Looking at Figs. 2 and 3, that depict the unitary deformation patterns, the displacement of a generic point of the structural member is defined by

2. The GBT energy formulation

iÞ axial displacements : uðx; zÞ ¼

2.1. General assumptions

where:

i) The member is prismatic, with straight longitudinal axis and constant cross section, and all cross section mechanical properties are constant along the member's length. ii) The member has linear elastic constitutive relation in the form:

sx ¼ E$εx txz ¼ G$gxz

i

uðzÞ$i f u ðxÞ

(7)

i¼1

In the present paragraph, the total potential energy (TPE) for the stability analysis of shear-deformable columns is derived in the context of GBT. The following assumptions are considered:



3 X

2

u ¼ ½ 1u

2u

1 3u  ¼ 4 « 1

z1 « zn

3 uðz Þ 1

«

3 5

(8)

3 uðz Þ n

iiÞ transverse displacements : wðx; zÞ ¼

3 X

i

wðzÞ$i f ðxÞ

(9)

i¼1

where:

2

(6)

where sx and εx are the longitudinal normal stresses and extensions, txz and gxz are the shear stresses and distortions, and E and G are the Young and the shear moduli, respectively. iii) The longitudinal extensions εx and the shear distortions gxz in (6) are derived from the Engineering Strains concept (Crisfield, 1991). iv) The loading is applied in such a way that no stress concentrations occur; hence, local stress variation effects can be ignored and the Saint-Venant principle can be considered valid (Dias da Silva, 2006). v) Cross sections keep their shape as the member deforms, so that cross section deformations in their own plane are neglected. vi) The column is considered to be continuously supported for displacements along the major axis plane, so that buckling along the weak axis is prevented and a two-dimensional analysis is applicable. These assumptions are the ones adopted by Engesser (Engesser, 1891) and Haringx (Haringx, 1949), and follow closely the ones ~o Coelho et al., adopted for the analysis of the Euler column in (Gira ~o et al., 2012), with the exception of assumption ii), since 2010; Sima now we assume that shear deformations may occur. 2.2. The nonlinear shear theory 2.2.1. The shear mode of deformation In compressed members with small shear stiffness, the symmetry condition that guided the Bernoulli's assumption of plane cross sections after deformation clearly does not occur (Dias da Silva, 2006). A closer look at the deformed shape of “weak-inshear” prismatic members puts in evidence that plane cross sections do not remain plane after deformation (Wang et al., 2000). However, the Timoshenko's beam formulation neglects this warping effect and keeps the cross sections plane, although they no longer remain perpendicular to the deformed longitudinal axis. Hence, a shear coefficient k, associated with the reduced cross section area Ar (Ar ¼ k$A), is required in order to compute an average shear stiffness that takes into account, indirectly, the variation of shear stresses along the cross section. Several strategies were presented to extend the Timoshenko for the case where cross section shear warping is considered. In the present work, this deformation pattern is modelled by an appropriate shear warping

w ¼ ½ 1w

2w

0 3w  ¼ 4 « 0

1 « 1

3 1 «5 1

(10)

Functions i f u ðxÞ and i f ðxÞ in (27) and (29) are the modal amplitude functions associated with mode i, respectively for the longitudinal and for the transverse displacements. They become the unknowns of the problem and are related between each other by, for any mode of deformation i, by (Dufort et al., 2001; Schardt, 1989): i

0

f u ðxÞ ¼ i f ðxÞ;

cx2½0; L

(11)

In the present paper, we derive the energy formulae explicitly in terms of the modal amplitude functions i f ðxÞ, following the usual procedures in GBT analysis (Schardt, 1989). This strategy unifies the notation for all modes of deformation, including those associated with cross sections deformations in their own plane in more complex problems, for which a physical meaning for the modal ~o, 2007). This preamplitude functions is not recognizable (Sima sentation may render the mathematical expressions less understandable for those less familiarized with GBT. Given expressions (7), (9) and (11) and to overcome this difficulty for the problem under observation, in the forthcoming paragraphs we suggest 0 simply to replace 1 f ðxÞ by u0 ðxÞ (axial displacement of the section's centre of gravity), 2 f ðxÞ by wB ðxÞ and 3 f ðxÞ by wSH ðxÞ. The task is now to define the appropriate modal warping displacements unitary pattern 3 uðzÞ along the cross section, in order to describe the displacements of any point of the structural member. Following most works on this subject, the modal unitary displacements pattern will be defined here for the rectangular cross section only, with depth h and width b, and for two relevant cases in the engineering practice: the homogeneous section and the composite member made from several layers of distinct materials with constant thickness. Generalization to other cross section shapes is left for future developments, and it shall be highlighted that the Timoshenko theory is modelled simply by assuming null values for the longitudinal displacements associated with mode 3 in expression (8).

2.2.1.1. The homogeneous member. Consider the prismatic member depicted in Fig. 1, made from a homogeneous isotropic material with elastic and shear moduli E and G, respectively. Following most works on this subject (Levinson, 1981, 1987; Bickford, 1982; Heyliger and Reddy, 1988; Petrolito, 1995; Shi and Voyiadjis, 2011), at first it is assumed that the cross section warping

~o / European Journal of Mechanics A/Solids 61 (2017) 216e234 P.D. Sima

220

Fig. 2. The unitary kinematic patterns for modes of deformation 1 (axial elongation) and 2 (bending).

Fig. 3. Mode 3: Shear deformations mode.

function 3u(z) is a third degree polynomial in the form 3

uðzÞ ¼ 3 u0 þ 3 u1 $z þ 3 u2 $z2 þ 3 u3 $z3

i) kinematic condition e null linear shear distortion at the top border of the cross section:

(12)

where constants 3 u0 to 3 u3 are yet to be determined. This option allows the displacements field to show a double curvature along the member's depth and to have a first derivative with parabolic form, in accordance to the classical parabolic shear stresses distribution (Dias da Silva, 2006). From expressions (7)e(10), the displacements 3u and 3w are given by:

3

wðx; zÞ ¼ 1$3 f ðxÞ ¼ 3 uðzÞ$3 f u ðxÞ

3 uðx; zÞ

(13)

h 0 3 _ z ¼  03 gxz ¼ 3 uðzÞ$ f u ðxÞ þ 1$3 f ðxÞ ¼ 0; 2

3

gxz ¼

v3 u v3 w 3 0 3 _ þ ¼ uðzÞ$ f u ðxÞ þ 1$3 f ðxÞ vz vx

(14)

where , denotes differentiation with respect to the transverse coordinate z. The warping function 3 uðzÞ is determined by imposing the following four conditions along the member's rectangular cross section, with depth h and width b:

(15)

Since this condition must be respected for all values of x along the whole member's length, it imposes that: 3

_ gxz jz¼h2 ¼ ½3 uðzÞ þ 1jz¼h ¼ 0 2

(16)

and also that 3

and the correspondent linear shear distortion is (Crisfield, 1991; Dias da Silva, 2006):

cx2½0; L

0

f u ðxÞ ¼ 3 f ðxÞ;

cx2½0; L

(17)

Expression (17) relates the amplitude modal functions for warping and transverse displacements and coincides with the Schardt's formula (11) (Schardt, 1989). This coincidence extends the Schardt's formula to shear modes of deformation, since it physically represents the boundary condition of null shear distortions at the contour along the whole member's length. ii) kinematic condition e null shear distortion at the base of the cross section:

~o / European Journal of Mechanics A/Solids 61 (2017) 216e234 P.D. Sima

3

_ z¼h þ 1 ¼ 0 uj

(18)

2

u_ I jz¼z1 þ 1 ¼ 0

Zh 2 =

1

u$3 u dA ¼

A

1$3 u dz ¼ 0

(19)

3

uI jz¼z2 ¼ 3 uII jz¼z2

=

Z

iii) continuity of the longitudinal shear stresses resultant at the interface between layers I and II:

Zh 2 =

2

u$3 u dz ¼

A

ðzÞ$3 u dz ¼ 0

(20)

h 2

  GI $3 gI jz¼z2 ¼ GII $3 gII jz¼z2 0GI $ 3 u_ I jz¼z2 þ 1   ¼ GII $ 3 u_ II jz¼z2 þ 1

=

For the rectangular homogeneous cross section, these four conditions yield the following warping displacements function for mode 3:

uðzÞ ¼

z 5 z3  $ 4 3 h2

(21)

and mode 3 becomes fully defined. Formula (21) is equal to one used by Vo and Thai (Vo and Thai, 2012; Thai, 2012) for the analysis of composite beams and nanobeams. In opposition to Levinson (Levinson, 1981, 1987), Bickford (Bickford, 1982) and Heyliger and Reddy (Heyliger and Reddy, 1988) strategies, the GBT formulation uncouples the longitudinal displacements between bending and shear. This uncoupling brings significant simplifications to the member's internal strain energy and contributes substantially to the efficiency of the analysis. 2.2.1.2. The composite member. The procedure presented just above for a homogeneous rectangular beam is extended now for prismatic composite beams made from several layers of distinct materials with constant thickness. It is assumed here that displacements continuity and equilibria of longitudinal shear forces are assured between any two adjacent layers, so that, for example, no delamination nor slipping phenomena occurs. The procedure is illustrated for the prismatic member with rectangular cross section made from three uniform layers of distinct materials: layer I between coordinates z1 and z2, layer II between z2 and z3, and layer III between z3 and z4. For each layer I, II and III, the materials have Young and shear moduli given by EI and GI, EII and GII, and EIII and GIII, respectively. The displacement of a generic point of the member is then given by:



(25)

h 2

iv) orthogonality between modes 2 and 3:

3

(24)

ii) continuity of the warping function at the interface between layers I and II:

iii) orthogonality between modes 1 and 3:

Z

3

221

3 uðx; zÞ 3 wðx; zÞ

0

¼ 3 ui ðzÞ$3 f ðxÞ; ¼ 1  3 f ðxÞ

i ¼ I; II or III

3

uII jz¼z3 ¼ 3 uIII jz¼z3

(27)

v) continuity of the longitudinal shear stresses resultant at the interface between layers II and III: 3

tII jz¼z3 ¼ 3 tIII jz¼z3 ⇔GII $3 gII jz¼z3

  ¼ GIII $3 gC jz¼z3 0GII $ 3 u_ II jz¼z3 þ 1   ¼ GIII $ 3 u_ III jz¼z3 þ 1

(28)

vi) null distortion at the bottom edge: 3

u_ C jz¼z3 þ 1 ¼ 0

(29)

vii) orthogonality with mode 1:

Z A

E$1 u$3 u dA ¼ EI

Zz2

1  3 uI dz þ EII

z1

Zz3

1  3 uII dz

z2

Zz4 þ EIII

(30)

3

1  uIII dz ¼ 0 z3

viii) orthogonality with mode 2:

Z

Zz2 E$ðzÞ$3 u dA ¼ EI

A

Zz3 ðzÞ$3 uI dz þ EII

z1

ðzÞ  3 uII dz z2

Zz4 ðzÞ  3 uIII dz ¼ 0

þ EIII z3

(23)

These warping functions must comply with the following boundary conditions: i) null distortion at the top edge:

iv) continuity of the warping function at the interface between layers II and III:

(22)

For each layer, a polynomial of third degree is adopted as warping function because it is pretended to allow change of concavity and a parabolic shear stress distribution along any layer of the cross section, as previously for the homogeneous section:

83 3 3 3 2 3 3 > < uI ðzÞ ¼ uI;0 þ uI;1 $z þ uI;2 $z þ uI;3 $z 3 uII ðzÞ ¼ 3 uII;0 þ 3 uII;1 $z þ 3 uII;2 $z2 þ 3 uII;3 $z3 > :3 uIII ðzÞ ¼ 3 uIII;0 þ 3 uIII;1 $z þ 3 uIII;2 $z2 þ 3 uIII;3 $z3

(26)

(31) In physical terms, no more conditions apply to the problem, hence the finding of the twelve unknowns 3 u ; i ¼ I; II or III; j ¼ 1; …; 4 becomes indeterminate. In physical i;j terms, this algebraic indeterminacy is a consequence of the hyperstatic character of the problem. Therefore, it is legitimate to apply the principle of minimum strain energy (Dias da Silva, 2006;

~o / European Journal of Mechanics A/Solids 61 (2017) 216e234 P.D. Sima

222

longitudinal extension is given by:

Przemieniecki, 1968) in the following way: i) At first, we express the internal strain energy d3 U stored along dx in the following form:

εx ¼ ¼

A2 B2  A1 B1 A1 B1 3 X

i

00

u$i f þ

i¼1

2

  vu 1 vw 2 z þ vx 2 vx

3 X 3 0 0 1X i w$j w$i f $j f 2 i¼1 j¼1

(34)

3 Zz3 Zz4 G G G 0 2 2 2 2 I II III 3 3 3 ð u_ I þ 1Þ dz þ ð u_ II þ 1Þ dz þ ð u_ III þ 1Þ dz5$ ½3 f ðxÞ dx d U¼4 2 2 2 z1 z2 z3 |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Zz2

3

(32)

d3 U

ii) Then, among all admissible functions for 3 ui ðzÞ, i ¼ I, II or III in expression (23), we will chose the one that minimizes the unitary internal strain energy:

2 min d U ¼ min4

Zz2

3

z1

Zz4 þ z3

GI 3 ð u_ I þ 1Þ2 dz þ 2 3

Zz3 z2

The shear distortion corresponds to change of the angle £ðBACÞ, from its initial configuration £ðB1 A1 C1 Þ, where it takes the value p2 by definition, to its final position £ðB2 A2 C2 Þ (Crisfield, 1991; Dias da Silva, 2006):

GII 3 ð u_ II þ 1Þ2 dz 2

GIII 3 ð u_ III þ 1Þ2 dz5 2

gxz (33)

submitted to restraints (24)e(31). A simple and consistent validation of this procedure is easily obtained by analysing a rectangular cross section with layers made from “different” materials with similar mechanical properties E and G. The procedure just presented leads to warping functions, for all layers, equivalent to formula (21). 2.2.2. The strains-displacements relations Taking into account Fig. 4, which illustrates the displacements of three generic points of the member A, B and C that move from position 1 to position 2 as the member deforms, and by the Engineering Strains concept (Crisfield, 1991; Dias da Silva, 2006) the

0

1  !  ! $A A B C 2 2 2 2 A ¼ g1 þ g2 ¼  £ðC2 A2 B2 Þ ¼  arccos@ 2 2 A2 B2 $A2 C2

p

p

  vu vw vu vu vw vw z þ þ $ þ $ vz vz vx vz vx |{z} vz ¼0

3 3 X 3 00 X X 0 i _ j i 0 j ¼ ði u_ þ i wÞi f þ u$ u$ f $ f i¼1

where

(35)

i¼1 j¼1

vw vz

¼ 0 is a consequence of assumption v) in paragraph 3.1.

2.2.3. The internal strain energy In paragraph 2.1, it was assumed that the member follows a linear elastic constitutive relation with respect to the Engineering Strains in the form:



sx ¼ E$εx txz ¼ G$gxz

(36)

Therefore, the internal strain energy is obtained simply by (Crisfield, 1991; Dias da Silva, 2006):

U ¼ UL þ USH ¼

1 2

Z

sx $εx dU þ

U

1 2

Z

txz $gxz dU

(37)

U

where U is the member's volume and εx and gxz are derived from expressions (34) and (35), respectively. Due to large extension of the resulting expressions, the detailed formulas for the internal strain energy are presented in Appendix B, together with the internal strain energy related to the Timoshenko formulation. Here, the internal strain energy is presented for the linear terms only, as follows:

U¼

EA 2

Z 

1

00

f

L

þ

G$3 D 2

2

Z

dx þ

EI 2

Z  L

2

00

f

2

dx þ

E$3 C 2

Z 

3

00

f

2

dx

L

0 2

ð3 f Þ dx L

Fig. 4. The displacements of three generic points of the column.

(38)

~o / European Journal of Mechanics A/Solids 61 (2017) 216e234 P.D. Sima

In expression (38), each part depends on a single mode of deformation and no coupling between modes occurs, owing to the orthogonality imposed between all modes. The first part is associated with axial elongation, the second part with bending and the third and fourth parts with the shear mode. 3 C is an additional cross section property that models the warping inertia for mode 3: 3

Z C¼

½3 uðzÞ2 dz

(39)

each modal amplitude function i f ; i ¼ 1; …; 3 is approximated by a linear combination of pre-established coordinate functions i 4j ðxÞ that must comply with the kinematic restraints of the problem, obedience to the static restraints being optional, in some cases even ~o, 2007). So, each modal amplitude function is not advisable (Sima approximated in the form: i

f ðxÞzi a1 $i 41 ðxÞ þ … þ i ani $i 4ni ðxÞ ¼

and plays, for the shear mode, the same role as the cross section area for the axial elongation, as the moment of inertia for bending, and, in the advanced theory of strength of materials, as the warping ~es da Silva, constant for torsion (Schardt, 1989; Sim~ ao and Simo 2004; Sim~ ao, 2007; Vlasov, 1961; Kollbrunner and Basler, 1969). 3 D is the shear stiffness that corresponds to the reduced cross section area Ar in the Timoshenko shear theory and to the torsional constant J in the torsion theory (Schardt, 1989; Vlasov, 1961), and is given by:

Z D¼

2 _ ½1 þ 3 uðzÞ dA

(40)

A

2.2.4. The potential of the external loading and the total potential energy For the structural system depicted in Fig. 1, it is considered that the applied axial load keeps its direction as the member deforms and it is applied to the member's initial perfect configuration. Consequently, for a compressive load applied at the edge sections of ~es the member, as in Fig. 1, its potential is simply (Sim~ ao and Simo ~o, 2007): da Silva, 2004; Sima 0

0

P ¼ P$1 f jx¼0  P$1 f jx¼L ¼ P

Z 1

00

f dx

(41)

L

Finally, the Total Potential Energy V of the system is simply computed by:

V ¼UP

(42)

where V is a functional on the modal amplitude functions i f ðxÞ (i ¼ 1,2,3 and on the load amplitude parameter P that appears in (42) always in a linear form. At the present stage, the traditional stability procedures (Thompson and Hunt, 1973; Hunt, 1981) and ~o et al., the natural discretization procedure (Sim~ ao, 2007; Sima 2003, 2012), in the context of the Rayleigh-Ritz method, are ready to be applied. 3. Application of the Rayleigh-Ritz method 3.1. Introduction The behaviour of a prismatic shear-deformable member under transverse or longitudinal loading in the two-dimensional space is completely described by the functional V of expression (42), which depends on the modal amplitude functions and on the load parameter P in the following form:

Z V¼

ni X

i

aj $i 4j ðxÞ

(44)

j¼1

A

3

223

F



1

 00 00 00 0 0 f ; 2 f ; 2 f ; 3 f ; 3 f ; P dx

(43)

L

The stability analysis of the compressed column member is made here in the context of the Rayleigh-Ritz method, by which

where iaj are the generalized coordinates of the problem, to be determined, and ni is the number of coordinate functions adopted for the modal amplitude function i f ðxÞ. The efficiency and numerical stability of the Rayleigh-Ritz method is strongly dependent of the judicious choice of the coordinate functions (Richards, 1977; Brown and Stone, 1997). Hence, the adopted pre-established modal coordinate functions i 4j ðxÞ are derived from a natural dis~o, 2007; Sima ~o et al., 2003, 2012) that cretization procedure (Sima creates sequentially a complete set of coordinate polynomials by imposing three types of conditions: the boundary conditions, the normality conditions with respect to the member's length and orthogonality conditions between a coordinate function being computed and all the precedent ones, in terms of the functions themselves or between their first or second order derivative in order to enhance the numerical stability of the problem (Brown and Stone, 1997). The adopted boundary conditions are the simply supported, the clamped-clamped and the clamped-sway conditions that are very common in practice, and the adopted coordinate functions for all modes and all cases are presented in Appendix C. After this transformation, the total potential energy becomes a nonlinear polynomial in the generalized coordinates iaj (after the transformation, a global numbering can be adopted for the generalized coordinates and subscript j can be suppressed) and in the load parameter P, which appears in V always in a linear form:

V ¼ Vði a; PÞ;

i ¼ 1; …; nC

(45)

where nC denotes the total number of adopted generalized coordinates. The traditional stability procedures (Thompson and Hunt, 1973; Hunt, 1981) are now ready to be applied, and the computed critical behaviours are described in the following paragraph.

4. The critical behaviour of the shear-deformable column 4.1. Validation: derivation of the Engesser formula by the GBT formulation The GBT-based stability analysis for the Euler-Bernoulli column under axial compression was already validated by deriving the traditional Euler buckling load formula (Gir~ ao Coelho et al., 2010) and by comparing the correspondent post-buckling behaviour with the ones arising from alternative theoretical formulations and ~o et al., 2012). In the following, the Finite Element solutions (Sima stability procedures of Thompson and Hunt (Thompson and Hunt, 1973; Hunt, 1981) are used to analyse the shear-deformable member under axial compression for the pin-ended column case and adopting the kinematical pattern associated with the Timoshenko linear shear theory. In the context of the Rayleigh-Ritz method, sinusoidal approximations for the modal amplitude functions of the bending and shear modes are chosen for now, in line with (Challamel et al., 2013), while for mode 1 the quadratic polynomial given by expression (88a) in Appendix C is adopted, in order to enable a constant normal force resultant along the

~o / European Journal of Mechanics A/Solids 61 (2017) 216e234 P.D. Sima

224

member's length, proportional to the modal amplitude's second derivative (Schardt, 1989). Hence, the modal amplitude functions are approximated by:

pffiffiffi 3 1 f ðxÞz1 a$ x2 2L

H FP ¼

vW vi q vj q

(46a)

 k q¼0

2 EA 3 6 L 6 6 6 ¼6 0 6 6 4 0

0

p



2

2L

Pþ

p2



2L

3

0

p EI

p

4



2L3 P

7 7 7 7 7 7 7 5

2

2L

p2 G Ar 2L

P



p2 2L

P (52)

2

3

f ðxÞz2 a$sin

f ðxÞz3 a$sin

px L

px L

(46b)

detðH FP Þ ¼

L aFP ¼ pffiffiffi P; 2 aFP ¼ 3 aFP ¼ 0 3 EA

(47)

(48)

The W-transformation (Thompson and Hunt, 1973; Hunt, 1981) introduces a new set of sliding coordinates iq, i ¼ 1,2,3 in the form:

8 L 1 1 1 > > P þ 1q > a ¼ aFP þ q ¼ pffiffiffi > 3 EA < 2a ¼ 2a þ 2q ¼ 2q > FP > > > :3 3 a ¼ aFP þ 3 q ¼ 3 q

(49)

and turns the total potential energy into the form:

Wði q; PÞ ¼ Vði aFP þ i q; PÞ

q ¼ 0;

i ¼ 1; 2; 3

Pcrit ¼

(50)

(51)

The equilibrium and the stability conditions pass over unchanged to the new total potential function W (Thompson and Hunt, 1973; Hunt, 1981), in this case along the whole fundamental path because, being linear with respect to the load parameter P, it does not present any limit point along its domain (Hunt, 1981). The Hessian matrix evaluated along the fundamental path HFP is equal to (i, j and k ¼ 1,2,3):

G Ar p2 EI p2 EI 1 ¼ PEngesser ¼ 2 $ 2 2 G Ar L þ p EI L 1 þ 2p2 EI

(54)

L GAr

which is the Engesser critical load given by expression (1). Consequently, the energy formulation is validated and an additional proof is obtained for the Engesser formula. Note that, for the derivation of expression (54), no assumptions were made with respect to the cross section shear stresses resultant, since the internal strain energy was computed in a way consistent with the concepts of the Theory of Elasticity (Crisfield, 1991), by integrating the semiproduct of the Engineering Strains with the correspondent conjugate stresses, a strategy commonly used in structural stability analysis. It was only assumed that the material obeys a linear elastic constitutive relation, in the present work with respect to the Engineering Strains, similarly to Engesser (Engesser, 1891) and to Haringx (Haringx, 1949). Finally, the form of the Hessian matrix in expression (52) proves that the axial elongation mode is always associated with passive coordinates.

4.2. Parametric studies 4.2.1. Introduction Having validated the GBT energy formulation for its simplest form by the derivation of the Engesser formula, now it is worth to analyse the influence of some relevant properties on the critical behaviour of prismatic compressed members that can undergo shear deformations, such as the adopted kinematic deformation patterns for the mode of deformation 3 and the modal boundary conditions. The buckling behaviour is described by plotting the shear-buckling factor CSH (Banerjee and Williams, 1994):

CSH ¼

for which the fundamental path is trivially defined by:

i

(53)

leads to the critical load Pcrit:

The equilibrium system is obtained by rendering stationary polynomial (47) with respect to the generalized coordinates (Thompson and Hunt, 1973). The fundamental path e the equilibrium path that emerges from the initial unloaded state e is simply determined by:

1

   6 3EA p6 GAr EI p EI p4 GAr P ¼0  þ 4 4 2 L 4L 4L 4L

(46c)

and the generalized coordinates ia, i ¼ 1,2,3 become the unknowns of the problem. Retaining only the second order terms of the internal strain energy derived in Appendix B and neglecting all nonlinear terms associated with the shear effects, together with replacing all amplitude functions by the ones in expression (46), leads to the following total potential energy:

pffiffiffi 3 EA1 2 p4 EI 2 2 p2 G Ar 3 2 3 p2 EA 1 2 2 $ 2 a$ a V¼ $ a þ a þ a þ 2 L 4 4L 4L3 L pffiffiffi pffiffiffi 3 p2 EA 1 2 3 3 p2 EA 1 3 2 pffiffiffi $ 2 a$ a$ a þ $ 2 a$ a þ 3$P$1 a þ 2 4 L L

and setting

Pcr PEuler

(55)

(Pcr is the real critical load accounting for shear deformations) against the slenderness ratio l ¼ Lie where i is the cross section radius of gyration (Chen and Lui, 1987). In the examples presented below, a prismatic member with rectangular cross section is adopted from Dufort et al. (Dufort et al., 2001). It is made from a carbon-epoxy material that presents a longitudinal Young modulus E ¼ 115 GPa, and, depending of the orientation of the carbon-fibre plies, the Young-shear moduli ratio E/G is assumed to vary between 2.6 (isotropic case) and 38.2. The reference cross section has width b ¼ 30 mm and high h ¼ 20.27 mm, and the member's length varies between 5 mm and 400 mm, which corresponds, for the simply supported case, to slenderness coefficients between 0.854 and 68.359. The stability procedures (Thompson and Hunt, 1973; Hunt,

~o / European Journal of Mechanics A/Solids 61 (2017) 216e234 P.D. Sima

1981) are applied here for more refined approximations, due to the adoption of several coordinate functions per mode. Instead of solving the Hessian matrix's characteristic equation as in paragraph 4.1, the critical behaviour is better determined by establishing a nonlinear eigenproblem in the form:

i h H FP $fi qg ¼ H 0 þ P$H 1 þ P 2 $H 2 $fi qg ¼ 0 FP

(56)

where matrix H 0 collects all parts of H FP independent of P, and matrices H 1 and H 2 collect the coefficients of P and P2 in H FP , respectively. This nonlinear eigenproblem is solved iteratively (Mirasso and Godoy, 1992), and the lowest positive eigenvalue is retained, since it is the pretended critical load, while its correspondent eigenvector defines the modal participation at the critical state. In all examples below, the amplitude modal function for axial elongation (mode 1) is approximated by a linear combination of coordinate functions (88) in Appendix C, and the compressive load is applied at point B, as depicted in Fig. 1. For the remaining modes of deformation, it was observed that the first five coordinate polynomials per mode from the sequential procedure (Sim~ ao et al., 2003, 2012) are sufficient to achieve convergence. Achieving convergence with the adoption of n coordinate functions means here that similar results are obtained when the adoption of the first n and n þ 2 polynomials of the sequence, for a modal amplitude function, renders similar critical loads. 4.2.2. The simply supported column The simply supported column is the reference case in column stability analysis, and is analysed here for two distinct conditions: i) edge sections attached to thick end plates (WR case), or ii) edge sections free to warp (WF case). The adopted coordinate functions for bending (mode 2) are always the polynomials (90) in Appendix C, while, for the shear mode (mode 3) the adopted functions are the polynomials (92) for case i) and polynomials (90) for case ii). Fig. 5 displays the shear-buckling factor CSH for GE ¼ 38:2 for several models, together with the theoretical formulas of Engesser, Haringx and Euler. We highlight the remarkable convergence between the Timoshenko models and the nonlinear model with edge sections free to warp with the Engesser buckling factor, while Haring's shear-buckling factor is always bigger. A similar coincidence occurs with the nonlinear model with warping prevented at both edge sections, with the exception of the very small lengths range (l < 5), for which the buckling factor becomes bigger than the one corresponding to the Engesser formula. We notice also that the Timoshenko formulation provides very similar buckling factors for cases i) and ii), meaning that the linear shear formulation is not able to

225

distinguish between edge sections with warping restricted or allowed, as a consequence of the neglect of shear warping in expression (8). Fig. 6 presents the influence of the Young e shear moduli ratio on the shear buckling factor, and, as expected, the higher this ratio the higher the degradation of the buckling load due to shear. Finally, Fig. 7 depicts the modal participation at the critical state for the nonlinear shear formulation, and it evidences that the shear mode becomes less relevant as the column's length increases. 4.2.3. Other boundary conditions Clamped-sway and clamped-clamped boundary conditions are very common in structural engineering, namely in elastomeric bearings. In the context of the Rayleigh-Ritz method, these conditions are modelled by choosing, as coordinate functions, polynomials (94) for both modes 2 and 3 for the clamped-sway case, and polynomials (92) for both modes 2 and 3 for the clampedclamped conditions. Fig. 8 compares the shear-buckling factors arising from the nonlinear shear formulation with the classical formulas of Engesser and Haringx and with the pinned-pinned column, and, once again, the computed critical loads agree much better with the Engesser formula than with Haringx's. One aspect deserves special attention: for the smallest lengths range, the shear buckling factor for the simply supported column with shear warping prevented at the edge sections is significantly higher than for the remaining cases. Looking at Fig. 9, where the buckling modes are depicted for all boundary conditions, it is concluded that warping restriction due to thick edge plates is more effective in areas along the member's length were the shear resultant force VSH is higher e making VSH ¼ vM vx , it is in the simply supported case and near the edge sections that higher values of VSH occur, while, for the clamped-sway and for the clamped-clamped columns, near the edge sections the bending moment reaches a local maximum, thus VSH is near zero and shear warping restriction has negligible effect on the critical load. Finally, the analysis for the clamped-clamped case was made for lengths between 10 and 800 mm, in order to cover the same slenderness ratios range of the other cases. 5. Concluding remarks and scope for further work In this study, the GBT energy formulation and the traditional stability procedures were applied to the buckling analysis of shear deformable columns under axial compression that fed much passionate discussion during the last decades. The analysis was performed in the context of a GBT strategy by establishing three modes of deformation (axial elongation, bending and shear

Fig. 5. Shear-buckling factor: comparison between the linear and nonlinear formulations with the formulas of Engesser and Haringx, for the pinned-pinned column with GE ¼ 38:2.

226

~o / European Journal of Mechanics A/Solids 61 (2017) 216e234 P.D. Sima

Fig. 6. Shear-buckling factor: influence of the E= ratio. G

Fig. 7. Modal participation in the critical state for the pinned-pinned column with GE ¼ 38:2 and with warping restriction at the edge sections, for the nonlinear shear formulation.

Fig. 8. Shear-buckling factor: comparison between different boundary conditions for

deformation), each one associated with a linearly independent kinematical deformation pattern defined along the member's cross section. A major attribute of the GBT formulation, in comparison to other strategies such as the Finite Element method, is that the kinematic description of the structural member is made by orthogonal modes of deformation that have a physical meaning: they correspond to the buckling modes (Schardt, 1994). Furthermore, the GBT strategy dispenses the prior choice of the direction of the cross section shear stresses resultant and follows the classical concepts of the Theory of Elasticity. A total potential energy was

E G

¼ 38:2

computed and the Rayleigh-Ritz method was used to render the problem discrete, by adopting coordinate polynomials to approximate the modal amplitudes. Then, the traditional stability procedures were applied to characterize the critical behaviour by plotting the shear-buckling factor against the member's slenderness for several boundary conditions. All computations were performed in the context of the symbolic algebraic manipulator Mathematica, fully exploring its abilities to manipulate very large and highly nonlinear expressions. It was intended, in this work, to use precisely the same

~o / European Journal of Mechanics A/Solids 61 (2017) 216e234 P.D. Sima

Fig. 9. Critical modes for the compressed column with l ¼ 8:545 and

assumptions of Haringx's and Engesser's formulations, and we must bear in mind that, in stability analysis, the critical load is a theoretical concept that shall have a mathematical reasoning. Therefore, the major conclusion is that, under these circumstances, the Engesser load agrees much better with the results arising from the GBT analysis, regardless of being always on the conservative side. It was also observed that the Timoshenko linear shear formulation provides a much poorer description of the member's behaviour than the nonlinear one, due to the neglect of shear warping. Nevertheless, these conclusions do not end the discussion of the influence of shear deformations on the buckling of columns,

E G

227

¼ 38:2.

and much is left to do. A thorough description of the stability behaviour of columns requires, for example, the determination of the post-buckling behaviour, a topic actually under development. The energy formulation developed in paragraph 3 is appropriate for the analysis of the initial post-buckling behaviour and requires little adaptation for imperfection sensitivity studies. In addition, the stabilizing effect of transverse stiffeners, that are common in elastomeric bearings, can be modelled by the adoption of appropriate coordinate functions, while the Poisson's effect can be represented by transverse extensions modes of deformations. Finally, it is worth to refer that the present formulation is adequate to the

~o / European Journal of Mechanics A/Solids 61 (2017) 216e234 P.D. Sima

228

buckling analysis of shear deformable members submitted to tensile loads, a recent subject that has been attracting much attention by many researchers worldwide (Kelly, 2003; Zaccaria et al., 2011; Caddemi et al., 2015, 2013, 2016) and that will be analysed in the near future.



00

M ¼ EI$wB VSH ¼ GAr $w0SH

(60)

The shear force VSH can be derived from the bending moment M by the traditional statics formula (Dias da Silva, 2006):

Acknowledgments ~o para a Cie ^ncia e a This work was supported by the Fundaça Tecnologia (FCT e Portuguese Foundation for Science and Technology) under project grant UID/MULTI/00308/2013. Appendices Appendix A. The derivation of the Engesser and of the Haringx formulae In this section, the stability analysis of a shear-deformable column by the classical equilibrium strategy is analysed, following closely the procedures presented in (Chen and Lui, 1987) for the Euler-Bernoulli column. Consider again the “weak-in-shear” column [AB] depicted in Fig. 1 that, without loss of generality, we can assume to be simply supported e generalization to other boundary conditions is trivial and provides similar results, with little modifications. It is also supposed that the column is perfectly straight and centrally loaded, follows a linear elastic constitutive relation and residual stresses are ignored. Starting by analysing Fig. 1-a), where the column [AB] is presented in its deformed configuration, global moment equilibria at A and B impose that VA ¼ VB ¼ 0, regardless the distance between A and B in the member's deformed configuration. Equilibrium of segment [AC] in its deformed shape imposes that MC ¼ wC $P, where wC is the total deflection of the generic point C due both to bending and to shear deformations, and also that the transverse force resultant QC, along the z-axis direction at C, is null. So, the only resulting force at the generic cross section C is the longitudinal force P that can be decomposed in two components, N and T, respectively normal and parallel to the deformed cross section shape. If q is the rotation angle of cross section C with respect to its original orientation (along the vertical axis Oz), regardless of whether the cross section remains perpendicular to the deformed longitudinal axis (Euler Bernoulli formulation) or not (Timoshenko or nonlinear shear formulations), we obtain:



and to transverse shear resultants VSH in the form:

VSH ¼

0 000 dM ⇔GAr $wSH ¼ EI$wB dx

0



00 EI ð4Þ w ¼ wSH GAr B

(61)

Introducing (60) and (61) in (59) yields the following equilibrium equation that is expressed in terms of the bending transverse displacements wB only, is independent of the column's boundary conditions and is valid for P < GAr: ð4Þ wB

! 00 P 1 þ $ wB ¼ 0 P EI 1  GA r

(62)

This equation follows closely the traditional fourth order equilibrium equation of the Euler column (Chen and Lui, 1987), and has a global solution given by:

wB ¼ a1 sinðm xÞ þ a2 cosðm xÞ þ a3 x þ a4

(63)

where

m¼

rffiffiffiffiffi P 1 $qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi EI 1 P

(64)

GAr

Application of the boundary conditions associated with the simply supported column to the global solution (63) yields the following values for the coefficients ai, i ¼ 1,..,4 (once again, generalization to other boundary conditions is trivial):

( x¼0 0

M ¼ 0 0 a2 ¼ 0

wB ¼ 0 0 a4 ¼ 0 8 a sinðm LÞ > < wB ¼ 0 0 a3 ¼  1 L x¼L 0 > : M ¼ 0 0 ða1 ¼ 0∧a3 ¼ 0Þ∨ðsinðmLÞ ¼ 0∧a3 ¼ 0Þ (65)

N ¼ P$cosðqÞ T ¼ P$sinðqÞ

(57)

Fig. 1-b) describes the equilibrium of an elementary column segment [DE] with infinitesimal length dx. The displacement w(x) includes the total transverse displacement due both to bending and to shear: ðnÞ

ðnÞ

wðxÞ ¼ wB ðxÞ þ wSH ðxÞ0wðnÞ ¼ wB þ wSH

(58)

where wB ðxÞ and wSH ðxÞ denote the transverse displacements due to bending and to shear deformations, respectively, and ðnÞ denotes the nth order derivative with respect to x. Equilibria for moments and for forces parallel to the z-axis are expressed by:

8 0

dM dM > 0 0 > > < Q ¼ dx  P$w ¼ dx  P$ wB þ wSH   00 > 00 00 d2 M > dQ d2 M > : ¼  P$w ¼  P$ wB þ wSH ¼ 0 2 2 dx dx dx

(59)

where Q is always perpendicular to P. The column obeys to linear elastic constitutive relations with respect to bending moments M

For the calculus of a1 , the boundary conditions impose that it can be either null, yielding the trivial solution with no physical relevance:

wB ¼ 0

∧

wSH ¼ 0;

cx2½0; L∧cP > 0

(66)

or non-null when (for n ¼ 1):

mL ¼ n p 0

rffiffiffiffiffi P L $qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ p EI 1 P

(67)

GAr

Solution of equation (67) with relation to P yields the Engesser critical load:

Pcr ¼

GAr $p2 EI p2 EI 1 ¼ 2 $ ¼ PEngesser 2 2 GAr L þ p EI L 1 þ p2 EI2

(68)

GAr L

Similar results are obtained if it is assumed that the shear force VSH is proportional to the first derivative of the total transverse displacement, thus related to the normal load P by (Allen, 1969):

~o / European Journal of Mechanics A/Solids 61 (2017) 216e234 P.D. Sima



VSH ¼ P w0B þ w0 SH ¼ GAr $w0SH

0

w0SH ¼

P G Ar  P

assumes the following general form:

$w0B (69)

Remembering the analysis of Fig. 1-a), expression (69) could also be obtained by assuming VSH ¼ dM and by rendering to zero the dx force Q in (59). This means that the Engesser formulation respects the equilibrium conditions (59). However and in line with all referred works that support the Haringx theory, if it assumed that cross section rotations are due only to bending and that shear deformations occur with null cross section rotations, following the Timoshenko deformations pattern (Timoshenko, 1921, 1922), by expression (60) the shear force is assumed to be tangential to the cross section's plane and equal to T by making w0B ¼ tanðqÞzsinðqÞ, as shown in Fig. 1-a). Taking into account the constitutive relation (60), VSH is now given by:

VSH ¼ TzP$w0B ¼ G Ar $w0SH

0

w0SH ¼

P $w0 G Ar B

  00 P P 1þ wB ¼ 0 þ EI G Ar

(71)

mL ¼ n p 0

(72)

Pcr

U

Therefore, Haringx's formula assumes, for the constitutive relation (60), a shear force VSH distinct of dM and not related to dx equilibrium condition (59). Furthermore, when the bending stiffness E I trends to infinity, the Engesser load trends to G Ar while Haringx's trends to infinity, and both critical formulas trend to the Euler load when the shear stiffness G Ar trends to infinity. The disagreement between both theories is clearly identified and is a consequence of the choice of the shear force: by comparing expressions (69) and (70), it is observed that Haringx theory choses a smaller shear force in the constitutive relation (60), hence it computes a larger critical load.

Z

txz $gxz dU

(75)

U

(76)

where the terms corresponding to a first order analysis are:

E 2

Z  U

Z

EA ¼ 2

vu vx

2

dU

2

ð1 f 00 Þ dx þ

EI 2

Z

2

ð2 f 00 Þ dx þ

E$3 C 2

L

Z

2

ð3 f 00 Þ dx

L

(77) The second order terms are the parts of third degree in the energy polynomial, as follows:

  Z 00 vu vw 2 EA 0 2 1 $ dU ¼ f $ð2 f Þ dx vx vx 2 L U Z Z 00 00 EA 2 2 0 3 0 1 1 þ EA f $ f $ f dx þ f $ð3 f 0 Þ dx 2

UL;2 ¼

E 2

Z

L

(78)

L

while the third order terms correspond to the parts of fourth degree and are:

Z  4 Z Z vw EA EA 0 4 0 3 0 dU ¼ ð2 f Þ dx þ ð2 f Þ $3 f dx vx 8 2 L L U Z Z 3 EA 2 0 2 3 0 2 2 0 3 0 3 þ EA ð f Þ $ð f Þ dx þ f $ð f Þ dx 4 2 L L Z EA 0 4 ð3 f Þ dx þ 8

UL;3 ¼ (74)

1 2

UL ¼ UL;1 þ UL;2 þ UL;3

(73)

Equation (73) has two roots, the positive one has physical relevance and is equal to:

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

GAr 1 GAr GAr L2 þ 4p2 EI þ ¼ 2L 2 0sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 G Ar @ 4p2 EI 1þ  1A ¼ PHaringx ¼ 2 G Ar L2

sx $εx dU þ

L

Application of the boundary conditions as in (65) leads to the following equation for the calculus of the critical load:

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  ffi P P 1þ ¼p L EI GAr

Z

1 2

where U is the member's volume and εx and gxz are derived from expressions (34) and (35), respectively, and their conjugate stresses are obtained by formula (36). Due to the large extension of the resulting expressions, we split the parts associated with the normal stresses from those related to the shear effects. So, for the longitudinal normal stresses and extensions, the internal strain energy is given by:

UL;1 ¼

and its solution is equal to expression (63), with the coefficient m given now by:

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi   P P 1þ m¼ EI G Ar

U ¼ UL þ USH ¼

(70)

The following equilibrium equation is thus obtained: ðivÞ wB

229

E 8

L

(79) The internal strain energy is completed with the parts related to the shear distortions gxz and the correspondent shear stresses txz. Like above and due to the large extension of the resulting formulas, this energy is split in three parts in the form:

USH ¼ USH;1 þ USH;2 þ USH;3

(80)

The first order terms are given by a single part: Appendix B. The internal strain energy for the nonlinear and for the Timoshenko formulations B.1: The internal strain energy for the nonlinear shear formulation As mentioned in paragraph 2.2.3, the internal strain energy

USH;1 ¼

G$3 D 2

Z

0 2

ð3 f Þ dx

L

The second order terms are given by:

(81)

~o / European Journal of Mechanics A/Solids 61 (2017) 216e234 P.D. Sima

230

Z USH;2 ¼ G U

  vu vu vu vw $ $ þ dU vx vz vz vx Z

¼ G$kSH2;1

0 2

00

ð3 f Þ $3 f dx þ G$kSH2;2

L

f $ f $ f dx

Z  0 2 f $ 3f dx  G$kSH2;4

1 00

L

 G$kSH2;5

1 00 2 0 3 0

L

Z

þ G$kSH2;3

Z

B.2. The Timoshenko linear shear theory The Timoshenko linear shear theory can be derived from the nonlinear theory simply by rendering null all warping displacements in expression (8), in the form:

2 00

0 2

2

f $ð3 f Þ dx u ¼ ½ 1u

L

Z

2 00 3 0 3 00

f $ f $ f dx þ G$kSH2;6

L

Z

2 0 2 00 3 0

f $ f $ f dx

L

(82) where the correspondent generalized cross section constants are listed in Table 1, together with their values for the homogeneous rectangular cross section b  h. The third order terms correspond to the parts of fourth degree of the polynomial:

Z  USH;3 ¼ G U

vu vx

þ2G$kSH3;4 G$kSH3;8 þ 2 G$kSH3;12

Z

L

f $ f $ f $ f dxG$kSH3;5

L

Z  L 2

Z

2 0 2 00 3 0 1 00

2

00

f

2

3 0 2

$ð f Þ dx2G$kSH3;9

G$kSH3;13 f $ð f Þ $ f dxþ 2 00

3 0 2 3

00

L

Z

Z L

2 0 3 0

0 2

f $ f $ð1 f Þ dxþ

L

G$kSH3;6 f $ð f Þ $ f dxþ 2

2 00

3 0 2 1 00

L

G$kSH3;3 2

Z

00 2

L

Z

2 0 2

00

00

f $ f $ f $1 f dxþG$kSH3;10

2 00 2

ð f Þ $ð f Þ dxG$kSH3;7

L

Z

2 0

00

3 0 2 3

00 2

f $ð2 f Þ $3 f 0dx

L

Z L

Z  00 2 ð f Þ $ 3f dxG$kSH3;14 2 0 2

0 2

ð1 f Þ $ð3 f Þ dx

L

2 0 3 0 3

(84)

The neglect of the shear warping displacements in expression (84) implies that the internal strain energy has much fewer parts than in the nonlinear shear theory, thus leading to a poorer description of the member's behaviour. Remembering that now a correction factor k is needed, associated with the reduced cross section area Ar ¼ k,A to compensate the constant shear stresses distribution along the cross section, the internal strain energy associated with the longitudinal extensions is given by:

2  2 Z Z G$kSH3;1 vu 0 2 00 2 $ dU ¼ ð2 f Þ $ð1 f Þ dxG$kSH3;2 vz 2 Z

2u

3 0 «5 0

z1 « zn

1 3u  ¼ 4 « 1

00

ð f Þ $ f $1 f dxþ2G$kSH3;11

Z

2 0 2

00

0

00

f $ f $3 f $3 f dx

L

Z  00 2  00 2 G$kSH3;15 0 2 f $ f $ 3f dxþ ð3 f Þ $ 3 f dx 2

2 0 3 0

L

L

(83) where the third-order cross section constants are also listed in Table 1. For the sake of simplicity, the internal strain energy and the correspondent generalized cross section properties were determined here for the homogeneous member. The generalization to composite members formed by nl layers of different materials is now trivial, but it must be noticed that the generalized properties become dependent of the materials elasticity constants Ei and Gi of each layer. Table 1 The generalized cross section properties and their correspondent values for the rectangular cross section b  h. Generalized property 1st order

3C

¼

R

½3 uðzÞ2 dA

A

2

nd

order

kSH2;1 ¼

R

3 u$ð1

_ 3 u_ dA þ 3 uÞ$

Value for the rectangular cross section b  h

Generalized property

bh3 1008

3D

¼

R

½1 þ

A

kSH2;3 ¼

A

3 u$ð1 _

_ dA þ 3 uÞ

R

_ 3 u_ dA kSH2;4 ¼ z$ð1 þ 3 uÞ$

0

A

kSH2;5 ¼ ð1 þ 3

order

3 uÞ$ _ 3u

dA

RA

kSH3;1 ¼ dA ¼ A

R

kSH2;6 ¼ z$ð1 þ

0

3 uÞ _

kSH3;2 ¼

bh

R

3 u_

dA

dA

A

_ 2 dA kSH3;3 ¼ ð3 uÞ A

R

_ 2 dA kSH3;5 ¼ z$ð3 uÞ

R

kSH3;4 ¼ z$3 u_ dA

bh 6

¼ Ar

5 bh 6

¼ Ar

0

kSH3;7 ¼ z2 $3 u_ dA A

kSH3;6 ¼ z2 d A ¼ I

0

A

3 u$3 u_

dA

3

bh 24 0

A

R

_ 3 u dA kSH3;11 ¼ z$3 u$ A

R

kSH3;13 ¼ ð3 uÞ2 dA A

R

_ 2 dA kSH3;15 ¼ ð3 uÞ2 $ð3 uÞ A

0 bh 6 0

A

R

A

R

kSH3;9 ¼

5 bh 6

A

A

R

R

Value for the rectangular cross section b  h

A

R

rd

R

_ kSH2;2 ¼ ð1 þ 3 uÞdA

0

A

R

2 3 uðzÞ _ dA

R

_ 2 dA kSH3;8 ¼ z2 $ð3 uÞ A

kSH3;10 ¼

R

3 u$ð3 uÞ _ 2 dA

bh3 12 5bh3 168

0

A

bh3 336 bh3 1008 71 3 199584 bh

R

_ 2 ,3 u dA kSH3;12 ¼ z$ð3 uÞ A

kSH3;14 ¼

R A

3 u$ð _ 3 uÞ2 dA

3

bh 378 3

bh 2592

~o / European Journal of Mechanics A/Solids 61 (2017) 216e234 P.D. Sima

UL ¼

1 2

Z

sx $εx dU ¼

U

EA 2

Z 

1

00

f

L

dxþ

EI 2

Z 

2

00

f

2

dx

EA 2

1

00

Z

0 2

1

f $ð2 f Þ dxþEA

L

rffiffiffiffiffiffi rffiffiffiffiffiffi 15 x2 10 x3 $ þ $ 2 L 3 L2

1

42 ¼ 

1

43 ¼

L

Z EA 1 00 3 0 2 f $ð f Þ dx 2 L L L Z Z Z EA 0 4 EA 0 3 0 3 0 2 0 2 ð2 f Þ dxþ ð2 f Þ $3 f dxþ EA ð2 f Þ $ð3 f Þ dx þ 8 2 4 L L L Z Z EA 2 0 3 0 3 EA 3 0 4 f $ð f Þ dxþ ð f Þ dx þ 2 8 þ

Z

2

00

0

0

f $2 f $3 f dxþ

rffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffi 105 x2 pffiffiffiffiffiffiffiffiffi x3 105 x4 $  210$ 2 þ $ 2 L 2 L3 L

231

(88b)

(88c)

L

(85)

Fig. 10. The adopted coordinate functions for mode 1 (first derivative), for L ¼ 3.

and the internal strain energy due to the shear distortions is equal:

USH ¼

1 2

Z U

GAr ¼ 2

txz $gxz dU ¼ Z

1 2

3 0 2

Z U

kG$g2xz dU Z

ð f Þ dx  G Ar L

G Ar þ 2

Z 

1

00

f

2

1

00

0

L 2 0 2

$ð f Þ dx þ k

L

0

f $2 f $3 f dx GI 2

Z

0 2

ð2 f Þ $



2

00

f

2

dx

L

C.2. The coordinate functions for modes 2 and 3 In the present paragraph, the coordinate functions that will be used later in the stability analysis for modes 2 and 3 are calculated for the simply supported, clamped-clamped and clamped-sway conditions. These conditions show similar critical buckling length for bending but, if it is assumed that the edge sections are connected to thick edge plates that prevent shear warping, it will be seen that the column may exhibit slightly different critical behaviours when shear deformations are included e in case of thick edge plates, it will be seen that the boundary conditions for mode 3 correspond to the clamped-clamped beam.

(86) We highlight, in expression (85), the absence of any linear term associated with mode 3.

C.2.1 The pinned-pinned boundary conditions. For the pinnedpinned member, only two kinematic boundary conditions apply for bending (mode 2):

Appendix C. The adopted coordinate functions

ðx ¼ 0 ∨ x ¼ LÞ02 f ¼ 0 C.1 The coordinate functions for the axial elongation mode In the context of the analysis of prismatic members using GBT, the axial elongation mode (mode of deformation 1) has distinct properties from the remaining ones (Schardt, 1989), hence it requires specific coordinate functions. Taking into account Fig. 1-a), it is considered that the member is fixed at one edge (point A at x ¼ 0) and free on the other (point B at x ¼ L), where the compressive load is applied, independently of the supporting conditions related to transverse displacements. So, only one boundary condition applies to the problem:

x¼0

0

1

1 0

f u jx¼0 ¼ f jx¼0 ¼ 0

(87)

and the natural discretization procedure enables the following coordinate functions, whose graphics are depicted in Fig. 10: 1

41 ¼

pffiffiffi 3 x2 $ 2 L

(88a)

(89)

and the first five coordinate functions given by the sequential procedure are: 2

41 ¼

pffiffiffiffiffiffi x x2  30  þ 2 L L

(90a)

2

42 ¼

 pffiffiffiffiffiffiffiffiffix x2 x3 210  3 2 þ 2 3 L L L

(90b)

2

43 ¼

 pffiffiffiffiffiffi x x2 x3 x4 70  3 þ 18 2  30 3 þ 15 4 L L L L

(90c)

 pffiffiffiffiffiffiffiffiffix x2 x3 x4 x5 (90d) 154  30 2 þ 90 3  105 4 þ 42 5 L L L L L  2 3 4 5 pffiffiffiffiffiffiffiffiffi x x x x x x6 2 45 ¼ 286  3 þ 45 2  210 3 þ 420 4  378 5 þ 126 6 L L L L L L (90e) 2

44 ¼

Their graphics, for L ¼ 3, are depicted in Fig. 11.

232

~o / European Journal of Mechanics A/Solids 61 (2017) 216e234 P.D. Sima

Fig. 11. Coordinate functions for modes 2 and 3, for L ¼ 3.

~o / European Journal of Mechanics A/Solids 61 (2017) 216e234 P.D. Sima

C.2.2 The clamped-clamped conditions. For modes 2 and 3 and if the structural member is clamped at both edges, the following four boundary conditions apply:

i f ¼0 ; ðx ¼ 0 ∨ x ¼ LÞ0 i 0 f ¼0

i ¼ 2; 3

(91)

These boundary conditions apply to modes 2 and 3 for the clamped-clamped member and also for mode 3 only in the simply supported case when edge cross sections are restrained by thick edge plates that prevent shear warping movements. Under these circumstances, the natural discretization procedure enables the following coordinate functions, the first five of the infinite sequence being given by (i ¼ 2; 3): i

41 ¼

 pffiffiffiffiffiffi x2 x3 x4 70 3 2  6 3 þ 3 4 L L L

(92a)

i

42 ¼

 pffiffiffiffiffiffiffiffiffi x2 x3 x4 x5 770  3 2 þ 12 3  15 4 þ 6 5 L L L L

(92b)

i

43 ¼

 pffiffiffiffiffiffiffiffiffiffiffiffiffiffix2 x3 x4 x5 x6 30030 2  6 3 þ 13 4  12 5 þ 4 6 L L L L L

(92c)

i

44 ¼

i

45 ¼

rffiffiffiffiffiffiffiffiffiffiffiffiffiffi 10010 x2 x3 x4 x5 x6  27 2 þ 240 3  795 4 þ 1242 5  924 6 41 L L L L L  x7 þ 264 7 L (92d) rffiffiffiffiffiffiffiffiffiffiffiffiffiffi 46410 x2 x3 x4 x5 x6 18 2  212 3 þ 975 4  2244 5 þ 2750 6 41 L L L L L  x7 x8  1716 7 þ 429 8 L L (92e)

Their graphics, for L ¼ 3, are depicted also in Fig. 11. C.2.3 The clamped-sway conditions. The clamped-sway supporting conditions are very common in many structural elements such as elastomeric bearings. In mathematical form, they are expressed by:

 x¼0

0

if

¼0 and x ¼ L f ¼0

i 0

0

i 0

f ¼0

(93)

and they are applicable both to mode 2 and to mode 3. The first five coordinate polynomials given by the sequential procedure are: 2

41 ¼

rffiffiffiffiffiffi  35 x3 x2 2 3 3 2 13 L L

2

42 ¼

pffiffiffiffiffiffi3 x4 6 x3 3 x2  70  3 þ 2 L4 L L

2

2

rffiffiffiffiffiffiffiffiffi  231 x5 x4 x3 x2 84 5  210 4 þ 170 3  45 2 43 ¼ 83 L L L L 44 ¼

 pffiffiffiffiffiffiffiffiffiffiffiffiffiffi x6 x5 x4 x3 x2 30030 4 6  12 5 þ 13 4  6 3 þ 2 L L L L L

(94a)

(94b)

(94c)

(94d)

2

pffiffiffiffiffiffiffiffiffi x7 pffiffiffiffiffiffiffiffiffi x6 pffiffiffiffiffiffiffiffiffi x5 45 ¼ 110 858 7  385 858 6 þ 518 858 5 L L L rffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffi x2 429 x4 910 286 x3  665 þ þ 140 858 2 2 L4 3 3 L3 L rffiffiffiffiffiffiffiffiffi 286 x2  455 3 L2

233

(94e)

Their graphics, for L ¼ 3, are depicted in Fig. 11. References Allen, H.G., 1969. Analysis and Design of Structural Sandwich Panels. Pergamon Press, Oxford. Aristizabal-Ochoa, J.D., 2007. Large deflection and post-buckling behaviour of Timoshenko beam-columns with semi-rigid connections including shear and axial effects. Eng. Struct. 29, 991e1003. Aristizabal-Ochoa, J.D., 2008. Slope-deflection equations for stability analysis and second-order analysis of Timoshenko beam-columns with semi-rigid connections. Eng. Struct. 30, pp. 2517e2527 and 3394-3395. Attard, M.M., 2003. Finite strain e beam theory. Int. J. Solids Struct. 40, 4563e4584. Attard, M.M., Hunt, G.W., 2008. Column buckling with shear deformations e a hyperelastic formulation. Int. J. Solids Struct. 4322e4339. Attard, M.M., Lee, J.-S., Kim, M.-Y., 2008. Dynamic stability of shear flexible Beck's column based on Engesser's and Haringx's buckling theories. Comput. Struct. 86, 2042e2055. Banerjee, J.R., Williams, F.W., 1994. The effect of shear deformation on the critical buckling of columns. J. Sound Vib. 174, 607e616. Ba zant, Z.P., 2003. Shear buckling of sandwich fibre composite and lattice columns, bearings and helical springs: paradox resolved. ASME e J. Appl. Mech. 70, 75e83. Ba zant, Z.P., Cedolin, L., 1991. Stability of Structures e Elastic, Inelastic, Fracture and Damage Theories. Oxford University Press, NY. Bickford, W.B., 1982. A consistent higher order beam theory. In: Chung, T.J., Karr, G.R. (Eds.), Developments in Theoretical and Applied Mechanics, vol. XI. Department of Mechanical Engineering e The University of Alabama in Huntsville, Huntsville, pp. 137e150. Blaauwendraad, J., 2008. Timoshenko beam-column buckling. Does Dario stand the test? Eng. Struct. 30, 3389e3393. Blaauwendraad, J., 2010. Shear in structural stability: on the Engesser-Haringx discord. ASME - J. Appl. Mech. 77 (3) (May), paper n. 031005, pp. 1e8. Brown, R.E., Stone, M.A., 1997. On the use of polynomial series with the RayleighRitz method. Compos. Struct. 39 (3e4), 191e197.  , I., Cannizzaro, F., 2013. The influence of multiple cracks on tensile Caddemi, S., Calio and compressive buckling of shear deformable beams. Int. J. Solids Struct. 50 (20e21), 3166e3183.  , I., Cannizzaro, F., 2015. Tensile and compressive buckling of Caddemi, S., Calio columns with shear deformation singularities. Meccanica 50, 707e720. , I., Cannizzaro, F., 2016. On the dynamic instability of shear Caddemi, S., Calio deformable beams under a tensile load. J. Sound Vib. 373, 89e103. Carrera, E., Giunta, G., 2010. Refined beam theories based on a unified formulation. Int. J. Appl. Mech. 2 (1), 117e143. Carrera, E., Giunta, G., Nali, P., Petrolo, M., 2010. Refined beam elements with arbitrary cross-section geometries. Comput. Struct. 88, 283e293. Carrera, E., Giunta, G., Petrolo, M., 2011. Beam Structures e Classical and Advanced Theories. John Wiley & Sons, Chichester. Challamel, N., Mechab, I., Elmeiche, N., Houari, M.S.A., Ameur, M., Atmane, H.A., 2013. Buckling of generic higher order shear beam-columns with elastic connections: local and nonlocal formulation. ASCE e J. Eng. Mech. 139 (8) (August), pp. 1091e1109. Chen, W.F., Lui, E.M., 1987. Structural Stability e Theory and Implementation. Prentice Hall, Upper Saddle River (NJ). Corp, Wolfram, 2010. Mathematica e Version 8. Wolfram Corp, Champaign (IL). Cowper, G.R., 1966. The shear coefficient in Timoshenko's beam theory. ASME e J. Appl. Mech. 33 (2), 335e340 (Jun 01, 1966). Crisfield, M.A., 1991. Non-linear Finite Element Analysis of Solids and Structures e Vol. 1: Essentials. John Wiley & Sons, Chichester. Dias da Silva, V., 2006. Mechanics and Strength of Materials. Springer Verlag, BerlinHeidelberg. diac, M., 2001. Closed-form solution for the cross-section Dufort, L., Drapier, S., Gre warping in short beams under three point bending. Compos. Struct. 52, 233e246. €be. Cent. Bauveraltung 49, Engesser, F., 1891. Die knickfestigkeit gerader sta 483e486. Forcellini, D., Kelly, J.M., 2014. The analysis of the large deformation stability of elastomeric bearings. ASCE e J. Eng. Mech. 140 (6) paper n. 04014036. Frank Pai, P., Schultz, M.J., 1999. Shear correction factors and an energy-consistent beam theory. Int. J. solids Struct. 36, 1523e1540. ~o, P.D., Bijlaard, F.S.K., 2010. Stability design criteria for steel Gir~ ao Coelho, A.M., Sima column splices. J. Constr. Steel Res. 66, 1261e1277. Haringx, J.A., 1949. Elastic stability of helical springs at a compression larger than

234

~o / European Journal of Mechanics A/Solids 61 (2017) 216e234 P.D. Sima

original length. Appl. Sci. Res. 1 (1), 417e434. Heyliger, P.R., Reddy, J.N., 1988. A higher order beam finite element for bending and vibration problems. J. Sound Vib. 126 (2), 309e326. Hunt, G.W., 1981. An algorithm for the nonlinear analysis of compound bifurcation. Philos. Trans. R. Soc. Lond. e Ser. A 300eA1455, 443e471. Kardomateas, G.A., Simitses, G.J., Shen, L., Li, R., 2002. Buckling of sandwich wide columns. Int. J. Non-Linear Mech. 37, 1239e1247. Kelly, J.M., 2003. Tension buckling in multilayer elastomeric bearings. ASCE e J. Eng. Mech. 129 (12), 1363e1368. Kelly, J.M., Konstantinidis, D.A., 2011. Mechanics of Rubber Bearings for Seismic and Vibration Isolation. John Wiley & Sons, Chichester. Kollbrunner, C., Basler, K., 1969. Torsion in Structures. Springer-Verlag, BerlinHeidelberg. Levinson, M., 1981. A new rectangular beam theory. J. Sound Vib. 74 (1), 81e87. Levinson, M., 1987. Consistent and inconsistent higher order beam and plate theories: some surprising comparisons. In: Elishakoff, I., Irretier, H. (Eds.), “Refined Dynamical Beams, Plates and Shells and Their Applications” e Proceedings of the Euromech Colloquium 219. Springer Verlag, Berlin-Heidelberg. Mirasso, A.E., Godoy, L.A., 1992. Iterative techniques for non-linear eigenvalue buckling problems. Commun. Appl. Numer. Methods 8, 311e317. Naifeh, A.H., Frank Pai, P., 2004. Linear and Nonlinear Structural Mechanics. Wiley, Weinheim. €nni, J., 1971. Das Eulerische Knickproblem unter Berücksichtigung der Querkr€ Na afte. ZAMP e Z. Angew. Math. Phys. 22 (1), 156e185. Petrolito, J., 1995. Stiffness analysis of beams using a higher-order theory. Comput. Struct. 55 (1), 33e39. Przemieniecki, J.S., 1968. Theory of Matrix Structural Analysis. McGraw Hill, NY. Reissner, E., 1972. On one-dimensional finite-strain beam theory: the plane problem. ZAMP - J. Appl. Math. Phys. 23, 795e804. Reissner, E., 1982. Some remarks on the problem of column buckling. Ing. Archiv. 52, 115e119. Richards, T.H., 1977. Energy Methods in Stress Analysis. Ellis Horwood Ltd, Chichester. Schardt, R., 1966. Eine Erweiterung der Technischen Biegetheorie zur Berechnung prismatischer Faltwerke. Der Stahlbau 35, 161e171. Schardt, R., 1989. Verallgemeinerte Technische Biegetheorie. Springer-Verlag, Berlin-Heidelberg. Schardt, R., 1994. Generalized Beam Theory e an adequate method for coupled stability problems. Thin-Walled Struct. 19, 161e180. Shi, G., Voyiadjis, G.Z., 2011. A sixth-order theory of shear deformable beams with variational consistent boundary conditions. ASME e J. Appl. Mech. 78 paper n. 021019. ~o, P.D., 2007. Post-buckling Bifurcational Analysis of Thin-walled Prismatic Sima Members in the Context of the Generalized Beam Theory (Ph.D. thesis). University of Coimbra, Coimbra.

~es da Silva, L., 2004. A unified energy formulation for the stability Sim~ ao, P.D., Simo analysis of open and closed thin-walled members in the framework of the Generalized Beam Theory. Thin-Walled Struct. 42 (10), 1495e1517. ~es da Silva, L., 2002. Comparative analysis of the stability of open and Sim~ ao, P., Simo closed thin-walled section members in the framework of Generalised Beam ~ es da Silva, L. (Eds.), Proceedings of the 3rd European Theory. In: Lamas, A., Simo Conference on Steel Structures e EUROSTEEL 2002, Coimbra, Portugal, 19-20 September. CMM, Coimbra, pp. 711e721. ~ es da Silva, L., 2003. Post-buckling behaviour of open cross-section Sim~ ao, P., Simo thin-walled columns in the context of the Generalized Beam Theory. In: Iu, V.P., Lamas, L.N., Pi, Y.-P., Mok, K.M. (Eds.), Computational Methods in Engineering and Science e Proceedings of the 9th International Conference on Enhancement and Promotion of Computational Methods in Engineering and Science e EPMESC IX, 25-28 November. Lisse: A. A. Balkema Publishers, Macao, pp. 891e898. ~o Coelho, A.M., Bijlaard, F.S.K., 2012. Influence of splices on the Sim~ ao, P.D., Gira buckling of columns. Int. J. Non-Linear Mech. 47, 806e822. Simo, J.C., Kelly, J.M., 1984. Finite element analysis of the stability of multilayer elastomeric bearings. Eng. Struct. 6, 162e174. Thai, H.-T., 2012. A nonlocal beam theory for bending, buckling and vibration of nanobeams. Int. J. Eng. Sci. 52, 56e64. Thompson, J.M.T., Hunt, G.W., 1973. A General Theory of Elastic Stability. John Wiley and Sons, London. Timoshenko, S.P., 1921. On the transverse vibrations of bars of uniform cross-section. Philos. Mag. 41, 744e746. Timoshenko, S.P., 1922. On the correction for shear of the differential equation for transverse vibration of prismatic bars. Philos. Mag. 43, 125e131. Touratier, M., 1991. An efficient standard plate theory. Int. J. Eng. Sci. 29 (8), 901e916. Tsai, H.-C., Kelly, J.M., 2005. Buckling of short beams with warping effect included. Int. J. Solids Struct. 42, 239e253. Vlasov, V.Z., 1961. Thin-walled Elastic Beams. The National Science Foundation and the Department of Commerce, Washington. Vo, T.P., Thai, H.-T., 2012. Vibration and buckling of composite beams using refined shear deformation theory. Int. J. Mech. Sci. 62, 67e76. Wang, C.M., Reddy, J.N., Lee, K.H., 2000. Shear Deformable Beams and Plates. Elsevier, Oxford. Wang, C.M., Zhang, Y.Y., Xiang, Y., Reddy, J.N., 2010. Recent studies on buckling of carbon nanotubes. ASME e Appl. Mech. Rev. 63 paper n. 030804 (18 pages). Zaccaria, D., Bigoni, D., Noselli, G., Misseroni, D., 2011. Structures buckling under tensile dead load. Proc. R. Soc. Lond. e Ser. A 467, 1686e1700. Zhang, Y.Y., Wang, C.M., Tan, V.B.C., 2006. Buckling of multiwalled carbon nanotubes using Timoshenko beam theory. ASCE e J. Eng. Mech. 132 (9), 952e958. Ziegler, H., 1982. Arguments for and against Engesser's buckling formulas. Ing. Archiv. 52, 105e113.