Inelastic analysis of beams on two-parameter tensionless elastoplastic foundation

Engineering Structures 48 (2013) 389–401

Contents lists available at SciVerse ScienceDirect

Engineering Structures journal homepage: www.elsevier.com/locate/engstruct

Inelastic analysis of beams on two-parameter tensionless elastoplastic foundation E.J. Sapountzakis ⇑, A.E. Kampitsis School of Civil Engineering, National Technical University, Zografou Campus, GR–157 80 Athens, Greece

a r t i c l e

i n f o

Article history: Received 21 May 2012 Revised 20 August 2012 Accepted 13 September 2012 Available online 27 November 2012 Keywords: Inelastic analysis Beam on foundation Two-parameter foundation Winkler model Pasternak model Distributed plasticity Boundary element method

a b s t r a c t In this paper a boundary element method is developed for the inelastic analysis of Euler–Bernoulli beams of simply or multiply connected constant cross-section having at least one axis of symmetry, resting on two-parameter tensionless elastoplastic foundation. The beam is subjected to arbitrarily distributed or concentrated vertical loading along its length, while its edges are subjected to the most general boundary conditions. A displacement based formulation is developed and inelastic redistribution is modeled through a distributed plasticity model exploiting material constitutive laws and numerical integration over the cross-sections. An incremental–iterative solution strategy along with an efficient iterative process are employed, while the arising boundary value problem is solved employing the boundary element method. Numerical results are worked out to illustrate the method, demonstrate its efficiency and wherever possible its accuracy. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction In engineering practice we often come across the analysis of beams on/in soil medium. The beam-foundation design is commonly used in modeling piles, pile-columns and pile groups embedded in soil medium as well as beam-columns and railway tracks resting on soil half space. Moreover, design of beams and engineering structures based on elastic analysis are most likely to be extremely conservative not only due to significant difference between initial yield and full plastification in a cross-section, but also due to the unaccounted for yet significant reserves of strength that are not mobilized in redundant members until after inelastic redistribution takes place. Thus, material nonlinearity is important for investigating the ultimate strength of a beam that resists bending loading, while distributed plasticity models are acknowledged in the literature [1–3] to capture more rigorously material nonlinearities than cross-sectional stress resultant approaches [4] or lumped plasticity idealizations [5,6]. According to the modeling of the mechanical behavior of the soil and the beam-foundation interaction, the earliest and probably the simplest mathematical model adopted is the single-parameter Winkler elastic foundation one [7]. In this approach the supporting soil is modeled as a single layer and its behavior is approximated

⇑ Corresponding author. Fax: +30 2107721720. E-mail addresses: [email protected] (E.J. Sapountzakis), cvakamb@ gmail.com (A.E. Kampitsis). 0141-0296/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.engstruct.2012.09.012

by a series of closely spaced, mutually independent, linear elastic vertical spring elements, providing resistance in direct proportion to the deflection of the beam. Although this model serves as a rather rough approximation of the true mechanical behavior of the soil reaction without being able to take into account the continuity or cohesion of the soil (interaction dependency between adjacent springs), it has been extensively used to solve many beam–foundation interaction problems. To overcome this weakness, several modified, multiple-parameter approaches have been proposed as the ones developed by Filonenko–Borodich [8], Hetenyi [9], Pasternak [10] and Vlasov and Leontiev [11] where a second parameter is introduced to account for the interaction among the adjacent springs. During past decades extensive research efforts have been presented concerning elastic analysis of beams resting on single or multi-layer elastic bilateral foundation models. The majority of the adopted models employed numerical methods such as the finite element method [12–17] while analytical formulations [18,19] have also been employed. However, elastic analysis of beams on elastic foundation, taking into account the realistic tensionless character of the subgrade reaction has received limited amount of literature. To begin with, Weitsman [20] studied the contact area developed between a beam under concentrated load and an elastic half-space. Kaschiev and Mikhajlov [21] presented a formulation based on the use of Newton’s method and the finite element method for beams resting on tensionless Winkler foundation, while for the same problem Zhang and Murphy [22]

390

E.J. Sapountzakis, A.E. Kampitsis / Engineering Structures 48 (2013) 389–401

presented an analytical/numerical solution making no assumption about either the contact area or the kinematics associated with the transverse deflection of the beam subjected to lateral point load assuming either free or pinned ends. Silveira et al. [23] presented a nonlinear semi-analytical methodology using a Ritz type approach, for the elastic equilibrium and instability analysis of beams, columns and arches resting on a tensionless Winkler elastic foundation, employing Newton-Raphson method together with arc-length iteration procedure in order to solve the nonlinear equations. Zhang [24] analyzed a beam resting on a tensionless Reissner foundation and demonstrated the improvements of the Reissner foundation model compared with the Winkler one, while Ma et al. [25] used the transfer displacement function method (TDFM) to present the response of an infinite beam resting on a tensionless Pasternak foundation subjected to linearly varying distributed loads. Finally, Sapountzakis and Kampitsis [26] developed a boundary element method for the nonlinear dynamic analysis of beamcolumns of arbitrary doubly symmetric simply or multiply connected constant cross-section, partially supported on tensionless Winkler foundation undergoing moderate large deflections under general boundary conditions, taking into account the effects of shear deformation and rotary inertia. To account for the nonlinear nature of the soil, several researches employed the concept of elastic beam on nonlinear foundation. In that formulation, the foundation load–displacement relation is assumed to follow a nonlinear law while the beam remains elastic throughout the analysis. Sharma and Dasgupta [27] employed an iteration method using Green’s functions for the analysis of uniformly loaded axially constrained Bernoulli beams hinged at both ends assuming an exponential build-up of the foundation reaction with the displacement increment. Beaufait and Hoadley [28] approximated the nonlinear load–displacement relationship of the Winkle foundation with a bilinear curve and utilized the midpoint difference method to analyze the beam coupled with the weighted averages scheme to estimate the spring stiffness for each iteration, followed by Yankelevsky et al. [29] who presented an iterative procedure based on the exact stiffness matrix for the beam on Winkler foundation by approximating the load–displacement curve by three to five regions rather than two. Kaliszky and Logo [30] adopted the extremum principle to analyze a nonlinear elastic beam on nonlinear elastic foundation. Both the beam and the Winkler springs were assumed to follow a bilinear material model while the beam was subdivided into series of rigid bars and the deformation was concentrated in the hinges and spring elements. Sapountzakis and Kampitsis [31] studied the nonlinear static analysis of shear deformable beam-columns of arbitrary doubly symmetric simply or multiply connected constant cross-section, partially supported on tensionless three parameter foundation, undergoing moderate large deflections under general boundary conditions, while in [32] the same authors confirmed the importance of the tensionless nonlinear three-parameter viscoelastic foundation in the nonlinear dynamic analysis of beams under the combined action of arbitrarily distributed or concentrated transverse moving loading. Contrary to the large amount of research concerning the elastic analysis of beams on either linear or nonlinear elastic foundation, only few studies have taken into account the inelastic behavior of both the beam and the foundation. According to this, the beam stress-strain and the foundation load–displacement relations are assumed to follow nonlinear inelastic constitutive laws. Consequently, such beam-foundation models are not commonly used due to the complexity of the problem. Ayoub [33] presented an inelastic finite element formulation for that is capable of capturing the nonlinear behavior of both the beam and the foundation. The element is derived from a two-field mixed formulation with independent approximation of forces and displacements and compared

with the displacement based formulation. Finally, Mullapudi and Ayoub [34] expanded the research in inelastic analysis of beams resting on two-parameter foundation where the values for the parameters are derived through an iterative technique that is based on an assumption of plane strain for the soil medium. In this paper, a boundary element method is developed for the inelastic analysis of Euler–Bernoulli beams of simply or multiply connected constant cross-section having at least one axis of symmetry, resting on two-parameter tensionless inelastic foundation. The beam is subjected to arbitrarily distributed or concentrated vertical loading along its length, while its edges are subjected to the most general boundary conditions. A displacement based formulation is developed and inelastic redistribution is modeled through a distributed plasticity model exploiting material constitutive laws and numerical integration over the cross-sections. An incremental–iterative solution strategy is adopted to resolve both the plastic part of stress resultants and the foundation reaction, while an efficient iterative process is employed to integrate the inelastic rate equations [35]. The arising boundary value problem is solved employing the boundary element method [36]. The essential features and novel aspects of the present formulation compared with previous ones are summarized as follows. i. The formulation is a displacement based one taking into account inelastic redistribution along the beam axis by exploiting material constitutive laws and numerical integration over the cross-sections (distributed plasticity approach). ii. The inelasticity of the soil medium is taken into account, employing a two-parameter foundation model. iii. The tensionless character of the foundation is also taken into consideration. iv. An incremental–iterative solution strategy is adopted to restore global equilibrium of the beam. v. The beam is supported by the most general nonlinear boundary conditions including elastic support or restrain, while its cross-section is an arbitrarily one having at least one axis of symmetry (z-axis). vi. To the authors’ knowledge, a BEM approach has not yet been used for the solution of the aforementioned problem, while the developed procedure retains most of the advantages of a BEM solution even though domain discretization is required. Numerical results are worked out to illustrate the method, demonstrate its efficiency and wherever possible its accuracy. 2. Statement of the problem 2.1. Displacements, strains, stresses Let us consider a prismatic beam of length l (Fig. 1a) of arbitrary constant cross-section having at least one axis of symmetry (zaxis), occupying the two dimensional multiply connected region X of the y, z plane bounded by the Cj (j = 1, 2, . . ., K) boundary curves, which are piecewise smooth, i.e. they may have a finite number of corners. In Fig. 1b Cyz is the principal bending coordinate system through the cross-section’s centroid. The normal stress-strain relationship for the material is assumed to be elastic-plastic-strain hardening with initial modulus of elasticity and yield stress E and rY0, respectively. The beam is resting on homogeneous tensionless elastoplastic two-parameter foundation. According to Pasternak hypothesis, the foundation reaction is expressed as

pf ¼ pw  p0p

ð1Þ

391

E.J. Sapountzakis, A.E. Kampitsis / Engineering Structures 48 (2013) 389–401

taken then E ¼ 1Em2 holds [37], while E is frequently considered instead of E⁄ (E⁄  E) in beam formulations [37,38]. This last consideration has been followed throughout the paper, while any other reasonable expression of E⁄ could also be used without any difficulty in many beam formulations. As long as the material remains elastic or elastic unloading oc  curs dexx ¼ deelxx the stress rate is given with respect to the strain one from Eq. (5). If plastic flow occurs then dexx ¼ deelxx þ depl xx , where the superscript pl denotes the plastic part of the strain component. A simplified Von Mises yielding criterion is considered ignoring the influence of shear stresses and the yield condition is satisfied when the normal stress is equated with the yield stress of the material, that is

(a)

f ¼ rxx  rY

(b) Fig. 1. Prismatic beam resting on an inelastic two-parameter foundation subjected to bending loading (a) with an arbitrary cross-section occupying the two dimensional region X (b).

with Winkler and Pasternak reaction related to the transverse displacement and rotation as

( pw ¼

kw w if pf > 0 0

if pf 6 0

( and pp ¼

if pf > 0 kp dw dx 0

if pf 6 0

ð3aÞ

2

d w

ð4aÞ

2

dx

cxz ¼ 0 ) hy ¼ 

dw dx

ð4bÞ

where (d/dx, d2/dx2) are the first and the second derivatives with respect to x. Considering strains to be small, employing the Cauchy stress tensor and assuming an isotropic and homogeneous material without exhibiting any damage during its plastification, the normal stress rate is defined in terms of the corresponding strain one as

drxx ¼ E deelxx

ð6Þ

where rY is the yield stress of the material and epl eq is the equivalent plastic strain, the rate of which is defined in [39] and is given as depl eq ¼ dk (dk is the proportionality factor [39]). Moreover, the plastic modulus h is defined as h ¼ drY =depl eq or drY = hdk and can be estimated from a tension test as h = EtE/(E  Et) (Fig. 2). The stress rate is given with respect to the total strain one through Eq. (4) and the strain components as

drxx ¼ Edexx  Edepl xx

ð7Þ

To establish global equilibrium equations, the principle of virtual work neglecting body forces is employed, that is

Z

ðrxx dexx ÞdV ¼

ð5Þ

where d() denotes infinitesimal incremental quantities over time (rates), the superscript el denotes the elastic part of the strain commÞ ponent and E ¼ ð1þEð1 mÞð12mÞ. If the plane stress hypothesis is under-

Z

V

l

ðpz dw  my dw0 Þdx 

Z l

½ðpw Þdw

þ ðpP Þdw0 dx

ð8Þ

where the integral quantities represent the strain energy, the external load and foundation reaction work while d() denotes virtual quantities, V is the volume and l is the length of the beam. The stress resultant corresponding to the internal bending moment of the beam is defined as

SMy ¼

Z

rxx zdX

ð9Þ

X

ð3bÞ

 and w  are the axial and transverse beam displacement where u components with respect to the Cyz system of axes; w(x) is the corresponding component of the centroid C and hy(x, t) is the angle of rotation due to bending of the cross-section with respect to its centroid. Employing the strain–displacement relations considering small deflections and adopting the Euler–Bernoulli assumption the following strain components are obtained

exx ¼ z



epleq ¼ 0

2.2. Equations of global equilibrium

ð2a; bÞ

where (0 ) denotes with respect to x, kw = kw(w, wy)    differentiation dw are the Winkler and Pasternak inelastic nonand kp ¼ kp dw ; dx dx y linear functions, respectively depending on the yielding quantities and the current ones. The beam is subjected to the combined action of arbitrarily distributed or concentrated transverse loading pz = pz(x) and bending moment my = my(x) acting in the x direction (Fig. 1a). Under the action of the aforementioned loading and having in mind that no axial loading is applied, the axial displacement vanishes and the displacement field of the beam is given as

 ðx; zÞ ¼ zhy u  wðxÞ ¼ wðxÞ



After substituting Eq. (9) into Eq. (8) and conducting some algebraic manipulations, the global equilibrium equations of the beam is obtained as 2



d SM y 2

dx

þ pf ðxÞ ¼ pz ðxÞ þ

dmy ðxÞ dx

ð10Þ

along with its corresponding boundary conditions

dSM y þ a2 w ¼ a3 dx dw b1 SM y þ b2 ¼ b3 dx

a1

ð11aÞ ð11bÞ

where pf is the foundation reaction defined from Eq. (1), while ai, bi (i = 1, 2, 3) are functions specified at the beam ends. The boundary conditions (11) are the most general ones for the problem at hand, including also the elastic support. It is apparent that all types of the conventional boundary conditions (clamped, simply supported, free or guided edge) may be derived from Eq. (11) by specifying appropriately the functions ai and bi (e.g. for a clamped edge it is a2 = b2 = 1, a1 = a3 = b1 = b3 = 0). Since an incremental–iterative approach is adopted for the problem at hand, the incremental version of Eqs. (10) and (11) is firstly written down as

392

E.J. Sapountzakis, A.E. Kampitsis / Engineering Structures 48 (2013) 389–401

(a)

(b)

Fig. 2. Normal stress–strain (a) and yield stress–equivalent plastic strain (b) relationships.

2



d DSM y 2

dx

þ Dpf ðxÞ ¼ Dpz ðxÞ þ

dDmy ðxÞ dx

ð12Þ

uðnÞ ¼

0

where D() denotes incremental quantities (over time), while the incremental stress resultant is given by virtue of Eqs. (9) and (7) as

DSMy ¼ EIy Dw00  DSMpl y

ð13Þ

where Iy is the moment of inertia with respect to the principle bending axis y, defined as

Iy ¼

Z

z2 dX

Depl zdX

Dw b1 DSMy þ b2 ddx ¼ Db3

#l 2 @u d u @ 2 u du @ 3 u  þ u @x dx2 @x2 dx @x3 0

ð18Þ

u ¼

1 3 ðjrj3  3ljrj2 þ 2l Þ 12

ð19Þ

Z

4

l

EIy

0

"

d u 4

dx

K4 ðrÞdx

 EIy K4 ðrÞ

3

d u 3

dx

2

 K3 ðrÞ

d u 2

dx

þ K2 ðrÞ

du  K1 ðrÞu dx

#l ð20Þ 0

where the kernels Kj(r) (j = 1, 2, 3, 4) are given as

ð16Þ

) at the beam ends x ¼ 0; l

3

dx



EIy uðnÞ ¼

inside the beam

þ a2 Dw ¼ Da3

3

d u

ð15Þ

2

dDSM a1 dx y

"

u dx

where u⁄ is the fundamental solution given as [36]

By substituting Eq. (13) in Eq. (12) and forming the incremental version of the boundary conditions (Eq. (11)), the following boundary value problem is obtained pl dDmy ðxÞ d DSM y  2 dx dx

4

dx

 u

X

EIy Dw0000 þ Dðpw Þ  Dðpp Þ0 ¼ Dpz ðxÞ þ

d u

with r = x  n, x, n points of the beam. Since EIy is independent of x, Eq. (18) can be written as

and SMpl y is the plastic quantity defined as

Z

4

l

ð14Þ

X

DSMpl y ¼ E

Z

ð17a; bÞ

By dropping the plastic quantities of the above equations, the boundary value problem of the examined problem under elastic conditions is formulated. 3. Integral representations – numerical solution

1 1 K1 ðrÞ ¼ sgnr K2 ðrÞ ¼ ðjrj  lÞ 2 2

ð21a; bÞ

1 1 3 ðjrj3  3ljrj2 þ 2l Þ K3 ðrÞ ¼ jrjðjrj  2lÞsgnr K4 ðrÞ ¼ 4 12

ð21c; dÞ

Solving Eq. (16) with respect to EIyDw0000 and substituting the result in Eq. (20), the following integral representation is obtained

! 2 pl dDmy ðxÞ d DSM y 0  EIy uðnÞ ¼ Dpz ðxÞ þ  Dðpw Þ þ Dðpp Þ K4 ðrÞdx 2 dx 0 dx " #l 3 2 d u d u du  EIy K4 ðrÞ 3  K3 ðrÞ 2 þ K2 ðrÞ  K1 ðrÞu dx dx dx 0 Z

l

ð22Þ

3.1. Integral representations for the displacement w According to the precedent analysis, the inelastic problem of beams resting on two-parameter tensionless elastoplastic foundation reduces to establishing the displacement component Dw(x) having continuous derivatives up to the fourth order with respect to x and satisfying the boundary value problem described by the governing differential Eq. (16) along the beam and the boundary conditions (17) at the beam ends x = 0, l. This boundary value problem (Eqs. (16) and (17)) is solved employing the BEM [36], as this is developed in [40] for the solution of a fourth order differential equation with constant coefficients. According to this method, let u(x) = Dw(x) be the sought solution of the problem. The solution of the fourth order differential equation d4u/dx4 = Dw0000 is given in integral form as [36]

After carrying out several integrations by parts and employing equations (2a,b), Eq. (22) yields

EIy uðnÞ ¼

Z 0



l

Dpz ðxÞK4 ðrÞdx  Z

l

0

þ

Z

Z

l

Dmy ðxÞK3 ðrÞdx

0

DSMpl y K2 ðrÞdx 

Z

l

Dðkw wÞK4 ðrÞdx

0

l

Dðkp wÞK2 ðrÞdx þ ½Dðkp w0 ÞK4 ðrÞ

0

 Dðkp wÞK3 ðrÞl0  l dDSMy du þ K4 ðrÞ  K3 ðrÞDSM y þ EIy K2 ðrÞ  EIy K1 ðrÞu dx dx 0 ð23Þ

393

E.J. Sapountzakis, A.E. Kampitsis / Engineering Structures 48 (2013) 389–401

Having in mind that the first derivative of u is contained in Eq. (17), Eq. (23) is differentiated once with respect to n yielding

EIy

duðnÞ ¼ dn þ

Z

l

Dpz ðxÞK3 ðrÞdx þ

0

Z

l

0

þ

Z

DSMpl y

K1 ðrÞdx 

Z Z

derivative is obtained from Eq. (24), while the second one from the following integral representation 2

l

EIy

Dmy ðxÞK2 ðrÞdx

0

d uðnÞ ¼ dn2

l

Z

l

Dpz ðxÞK2 ðrÞdx  0



Dðkw wÞK3 ðrÞdx

Z



0

l

0

Dðkp wÞK1 ðrÞdx  ½Dðkp w ÞK3 ðrÞ

l

Dmy ðxÞK1 ðrÞdx  DSM pl y

l

Dðkw wÞK2 ðrÞdx þ Dðkp wÞ

0

þ K2 ðrÞ

0

Z

dDSM y  K1 ðrÞDSM y þ Dðkp w0 ÞK2 ðrÞ  Dðkp wÞK1 ðrÞ dx

0

ð24Þ

Observing Eqs. (23) and (24), it is deduced that they have been brought into a convenient form to establish a numerical computation of the unknown quantity Dw. Thus, the interval (0,l) is divided into L elements, on each of which Dw is assumed to vary according to a certain law (constant, linear, parabolic etc). The linear element assumption is employed here (Fig. 3) as the numerical implementation is simple and the obtained results are very good. It is worth here noting that this technique does not require either differentiation of shape functions or finite differences application. Employing the aforementioned procedure and a collocation technique, a set of L + 1 algebraic equations is obtained with respect to L + 7 unknowns, namely the values of (Dw)i, (i = 2, 3, . . ., L) at the  nodal  points and the boundary values of  DLw 1 internal dDSMy ðDwÞj ; ddx ðDSMy Þj ; (j = 1 and L + 1) at the beam ends dx j j n1 = 0, nL+1 = l (Fig. 3). Two additional algebraic equations are obtained by applying the integral representation (24) at the beam ends n = 0, l. These L + 3 equations along with the four boundary conditions (Eq. (17)) yield a linear system of L + 7 simultaneous algebraic equations

ð25Þ

where [K] is a known generalized stiffness matrix, {Dd} is a generalized incremental unknown vector given as



 d Dw fDdg ¼ ðDwÞ1 ðDwÞ2 . . . ðDwÞLþ1 dx



1 dDSM y dDSM y DSM y Þ1 ðDSM y ÞLþ1 dx dx 1 Lþ1 T

0

ð27Þ

 Dðkp wÞK2 ðrÞl0  l dDSMy du  K3 ðrÞ  K2 ðrÞDSM y þ EIy K1 ðrÞ dx 0 dx

½KfDdg ¼ fDbext g þ fDbpl g

l



dDw dx Lþ1 ð26Þ

while {Dbext}, {Dbpl} are known vectors representing all the terms related to the incremental externally applied loading and incremental plastic quantities, respectively. After solving the system of Eq. (25), a post-processing step is re2 quired to obtain the derivatives du ; ddx2u at any nodal point ni (i = 1, 2, dx . . ., L + 1) which is essential to the solution algorithm. The first

which is obtained after differentiating Eq. (23) twice with respect to n. The presented integral representations may also be employed to compute u(x) and its derivatives at any interior point of the beam other than ni (i = 1, 2, . . ., L + 1). 3.2. Incremental–iterative solution algorithm Load control over the incremental steps is used and load stations are chosen according to load history and convergence requirements. Incremental stress resultants are decomposed into elastic and plastic part (Eq. (13)). They are computed through an iterative procedure since usually, changes between the plastic part of incremental stress resultants of two successive iterations are not negligible. Thus, using the subscript m to denote the load step, the superscript l to denote the iterative cycle and the symbol D() to denote incremental quantities, the lth iteration of the mth load step of the incremental–iterative solution algorithm will be described. A step-by-step algorithmic approach of the nonlinear solution is presented in a flowchart form in Fig. 4 i. Evaluation of the generalized iterative unknown vector fDdglm from the solution of the linear system of equations (Eq. (25))

½KfDdglm ¼ fDbext gm þ fDbpl gl1 m

ð28Þ

If m = 1 and l = 1, it is assumed that fDbpl g01 ¼ f0g. If m > 1 and l = 1, it is assumed that fDbpl g0m ¼ f0g or fDbpl g0m ¼ fDbpl gnm1 , where n is the total number of iterations performed in the previous increment m  1. ii. Evaluation of the incremental unknown derivatives ½Dw0 ðni Þlm ; ½Dw00 ðni Þlm ði ¼ 1; 2; . . . ; L þ 1Þ

by

introducing

fDdglm into Eqs. (24) and (27), respectively. iii. Elastic prediction step: For each monitoring station k of the ith cross-section of the beam (k = 1, 2, . . ., Ndof, i = 1, 2, . . ., L + 1): evaluation of the trial stress components as



l

rTrxx ðni ; zk Þ

m

 l ¼ ðrxx ðni ; zk ÞÞ0m þ DrTr xx ðni ; zk Þ m

Fig. 3. Discretization of the beam interval into linear elements, distribution of the nodal points and approximation of several quantities.

ð29Þ

394

E.J. Sapountzakis, A.E. Kampitsis / Engineering Structures 48 (2013) 389–401

iv. Perform the yield criterion at each monitoring station k of the ith cross-section of the beam (k = 1, 2, . . ., Ndof, i = 1, 2, . . ., L + 1) employing Eq. (6).  If fTr 6 0 then the trial state is the final state and the incremental plastic strain components along with the equivalent plastic strain are updated as

Start Beam/Foundation Inelastic Data



Increment m

Depl xx

l m

¼0



epleq

l

¼

m



epleq

0

ð31a; bÞ

m

 If fTr > 0 then plastic flow occurs and return must be made to yield surface (plastic correction step). The plastic flow rule and the loading/unloading conditions can easily be formulated [41,42]. v. For each monitoring cross-section of the beam i: Evaluation h il of the plastic quantities DSM pl ði ¼ 1; 2; . . . ; L þ 1Þ by y ðni Þ

Loop over Foundation Convergence Iteration l

m

Global Convergence

Next Cross Section i=i+1

Plastic flow rule loading/unloading

No

Next Gauss Points k=k+1

Von Mises Criterion

Next Increment m=m+1

Gauss Points k

Next Iteration l=l+1

Cross Section i

Next Foundation Loop

Generalized iterative unknown vector [Δwi,Δw'0,l,ΔM0,l,ΔQ0,l]

employing a two-dimensional numerical integration scheme to approximate the domain integrals of Eq. (15). In this study, the beam’s monitoring cross-sections are divided into a number of triangular or quadrilateral cells and standard two-dimensional Gauss quadrature rules are employed in each cell. Thus, the monitoring stations of each cross-section coincide with the Gauss points of its cells. If the same number of Gauss points is employed in every cell, then Ndof = Ncells  NGauss holds. Exact patch between adjacent cells is not required due to the combined use of BEM to compute Iy and the fact that only domain integrals are approximated without performing any structural discretization.h i vi. Employ the obtained plastic quantities DSMpl y ðni Þ ði ¼ 1; 2; . . . ; L þ 1Þ to evaluate (a) the vector fDbpl glm representing all the terms of Eq. (23), related to plastic quantities and (b) terms of Eqs. (24) and (27) related to plastic quantities required to perform step (ii) for the next iteration l + 1. Apart from elementary computations, the current step requires the computation of line integrals of the form Rl Kj ðrÞDSMpl y dx ðj ¼ 1; 2Þ (Eqs. (23) and (24)). A numerical 0 integration scheme must be employed to resolve these integrals since plastic quantities are not known in the whole beam interval (0,l). A semi-analytical scheme has been l m

No

implemented, according to which DSMpl y vary on an element k (k = 1, 2, . . ., L) (Fig. 2) of the beam interval following the same law that is used to approximate Dw (see Section 3.1). This leads to the integration of kernels being products of functions Ki(r) and two-node linear shape functions, thus it is performed analytically without any difficulty.

h il

pl vii. Check convergence. Convergence occurs if

DSM y ðni Þ m  h il1 DSMpl k is sufficiently small. The iterations continue y ðni Þ

Total displacement Stress resultants Foundation reaction

Foundation Convergence

No

m

until the specified accuracy is reached. If convergence is achieved after n iterations then: 1. For each monitoring station k of the ith cross-section of the beam (k = 1, 2, . . ., Ndof, i = 1, 2, . . ., L + 1): Initialize the stress components along with the equivalent plastic strain for the next increment m + 1 as

Final Total Results End

ðrxx Þ0mþ1 ¼ ðrxx Þnm

Fig. 4. Flowchart of the incremental–iterative solution algorithm.

where the incremental trial stress components are obtained by employing Eqs. (4) and (5) as

" #l 2  Tr l d Dwðni Þ Drxx ðni ; zk Þ m ¼ E  ðzÞk 2 dx m

ð30Þ



epleq

0 mþ1

¼



epleq

n m

ð32a; bÞ

2. Resolve (a) the vector {Dbext}m+1 representing all the terms of Eqs. (17) and (23) related to externally applied loading and (b) terms of Eqs. (24) and (27) related to externally applied loading required to perform step (ii) for the next increment. Apart from elementary computations, the current step requires the computation of line integrals of the Rl Rl form 0 Ki ðrÞDpz dx; 0 Ki ðrÞDmy dx ði ¼ 1; 2; 3; 4Þ. Since the distributions of pz and my are usually prescribed in codes

395

E.J. Sapountzakis, A.E. Kampitsis / Engineering Structures 48 (2013) 389–401

Elastic Beam

My (l/2) 1.0

x kp

z

kw

0.4

l=5m Fig. 5. Prismatic elastic beam on elastic two parameter foundation subjected to a bending moment at its midpoint.

Table 1 Deflection (mm) and rotation (rad) of the beam of example 1.

Deflection w(l) Rotation w0 (l/ 2)  103

Pasternak

Present analysis

Mullapudi and Ayoub [34]

Present analysis

Mullapudi and Ayoub [34]

3.88

3.90

1.32

1.30

1.599

1.600

0.5849

0.5850

Bending Moment M (kNm)

Winkler

180

140

100

60

Present Study Mullapudi & Ayoub [34]

30

Bending Moment M y (kNm)

20 20

0 10

0.004

0.008

0.012

Rotation w'(l/2) (rad) Fig. 7. Moment vs. rotation at the midpoint of the beam of example 1, for tensionless Winkler foundation.

0 -10 -20 -30 0

1

2

3

4

5

Length Along the Beam (m) Present Study-Winkler Present Study-Pasternak Mullapudi & Ayoub [34]-Winkler Mullapudi & Ayoub [34]-Pasternak Fig. 6. Bending moment distribution My(x) along the beam of example 1.

and regulations with simple analytical relations, these integrals are evaluated analytically, demonstrating the efficiency of the developed numerical procedure (e.g. concentrated loads may be treated using the Dirac function, without adhering to any simplifications). viii. Since convergence is achieved then the foundation reaction is computed employing a simplified Von Mises yielding criterion. The parameters are updated and the process described by steps (i)-(vii) is repeated until the foundation convergence criterion is achieved by using a prescribed tolerance of tolfound = 1010. ix. The increments of the external loading continue till total loading is undertaken or till convergence cannot be satisfied, which means that the last additional increment cannot be undertaken (plastic collapse). 4. Numerical examples

Table 2 Deflection (cm) at both ends of the beam of example 1, for tensionless Winkler foundation. M (kNm)

w(0) w(l)

Present analysis 20

50

100

150

168

0.555 0.710

0.245 1.020

0.180 1.425

1.510 2.154

2.633 2.570

Mullapudi and Ayoub [34] – M = 168 kNm w(0) w(l)

2.490 2.450

On the basis of the analytical and numerical procedures presented in the previous sections, a computer program has been written and representative examples have been studied to demonstrate the efficiency, the range of applications and wherever possible the accuracy of the developed method. 4.1. Example 1 For comparison purposes, in the first example the elastic analysis of a free-free beam of length l = 5 m resting on Pasternak elastic

396

E.J. Sapountzakis, A.E. Kampitsis / Engineering Structures 48 (2013) 389–401

Table 3 Maximum values of the deflection wmax(l/2) of the beam of example 2 and divergence values for various discretization schemes. Nodal points

6

11

17

20

Ayoub [33]

Max deflection (cm) Divergence (%)

3.894 2.03

3.906 1.73

3.937 0.95

3.925 0.01

3.975 –

Displacement w (in)

160

Load Pz (l/2)

0

1

200

2

3

120

4

80

5 0

Present Study Ayoub [33]-Converged Solution Ayoub [33]-Mixed Model_Cubic Moment Ayoub [33]-Mixed Model_Linear Moment Ayoub [33]-Displ. Model_5 th order Displ.

40

Ayoub [33]-Displ. Model_Cubic Displ. 0 0

1

2

3

4

5

6

Displacement w (l/2) Fig. 8. Load–displacement curve at the midpoint of the beam of example 2.

foundation subjected to a concentrated bending moment My = 50 kNm acting at its midpoint, as shown in Fig. 5 is examined. The beam is made out of timber with elastic modulus E = 10,500 MPa, Poison ratio v = 0.25and has rectangular cross-section of width

3000

2000

Bending Moment My (k.in)

Present Study Ayoub [33]-Converged Solution Ayoub [33]-Mixed Model_Cubic Moment Ayoub [33]-Displ. Model_Cubic Displ.

1000

0

-1000

-2000 0

100

200

300

Length (in) Present Study Ayoub [33]-Converged Solution Ayoub [33]-Mixed Model_Cubic Moment Ayoub [33]-Displ. Model_Cubic Displ. Fig. 9. Bending moment distribution My(x) along the beam of example 2.

100

200

300

Length (in) Fig. 10. Displacement w(x) along the beam of example 2.

b = 0.4 m and depth h = 1.0 m. The elastic foundation is sandy clay with modulus of elasticity Es = 45.5 MPa and Poisson ratio vs = 0.21. The values of Winkler and Pasternak foundation parameters are evaluated as kw = 3.081 MPa and kp = 12,449 kN, respectively according to both Zhaohua and Cook [13] and Mullapudi and Ayoub [34] considering vs = 0.25 and c = 1.0. The present example was first studied by Shirima and Ginger [15] who presented a complete solution of the stiffness matrix and nodal action vectors for a Timoshenko beam element resting on a two-parameter elastic foundation employing the displacement method. Later Mullapudi and Ayoub [34] solved the same problem deriving the values of the foundation parameters through an iterative technique that is based on the plain strain assumption for the soil medium, while the beam is discretized into four mixed elements with cubic moment interpolation functions, assuming five integration points for each finite element. In Table 1, the evaluated deflection at the beam’s right end w(l) and the midpoint rotation w0 (l/2) for both Winkler and Pasternak formulations are presented as compared with those obtained from the literature (the compared values have been extracted from a graph). Moreover, in Fig. 6 the bending moment distribution along the beam length is also presented as compared with available results from the literature for the aforementioned cases, demonstrating the accuracy of the proposed method in elastic analysis and noting that the bending moment is slightly underestimated if ignoring the Pasternak effects. In order to demonstrate the effect of the unilateral character of the soil the same free-free elastic beam resting on tensionless Winkler foundation is analyzed, being subjected to a concentrated vertical force of P(l/2) = 100 kN and to a concentrated bending moment M(l/2) acting at its midpoint. In Table 2 the deflection of both ends for various values of the applied moment and in Fig. 7 the moment-rotation curve at the beam’s midpoint are presented as compared (wherever possible) with those obtained from Mullapudi and Ayoub [34]. A very good agreement is once again verified.

4.2. Example 2 As a second example, also for comparison purposes, a simply supported beam of length l = 300 in. and square cross-section of

397

E.J. Sapountzakis, A.E. Kampitsis / Engineering Structures 48 (2013) 389–401 Table 4 Curvature (1/in) and foundation reaction (k/in) of the beam of example 2.

300

Present Ayoub [33] analysis Converged Mixed model Displ. model solution Cubic Linear 5th Cubic moment moment order 9.98

10.66

11.85

0.601

0.60

0.60

9.95

4.91

2.26

0.611

0.619 0.689

200

Load (kN/m)

Curvature w00 (l/2)  103 (1/in.) Foundation reaction pf(l/6) (k/in.)

100

Present Study - Perfectly Plastic Present Study - Strain Hardening FEM Solid Model [43] - Perfectly Plastic FEM Solid Model [43] - Strain Hardening FEM Beam Model [43] - Perfectly Plastic FEM Beam Model [43] - Strain Hardening

1000

0

800

0

0.005

0.01

0.015

0.02

0.025

Load Pz (2l/3)

Displacement w (l/2) 600

Fig. 12. Load–displacement curve at the midpoint of the beam of example 3.

400

Present Study - Perfectly Plastic Present Study - Strain Hardening FEM Solid Model [43] - Perfectly Plastic FEM Solid Model [43] - Strain Hardening FEM Beam Model [43] - Perfectly Plastic FEM Beam Model [43] - Strain Hardening Step by Step

200

0 0

0.004

0.008

0.012

0.016

Displacement w (2l/3) Fig. 11. Load–displacement curve at the loading point of the beam of example 3.

side d = 6.26 in. subjected to a monotonically increasing concentrated vertical load P at its midpoint has been studied. The material of the beam is assumed to follow an elastic–plastic behavior with modulus of elasticity E = 29,000 ksi, yielding stress rY0 = 30 ksi and a strain hardening slope of 1.4% (tangent modulus Et = 406 ksi), while the Winkler foundation load–displacement relation is also considered to be elastic–plastic with initial stiffness equals to kw = 0.5 kip/in.2, yielding force PwY = 1.0 k/in. and hardening slope of 1.0% (tangent stiffness kwt = 5  103 kip/in.2). For the longitudinal discretization 20 linear elements have been employed, while the cross-section has been discretized into 36 quadrilateral cells (6 fibers) and a 2  2 Gauss integration scheme has been used for each cell. The present example has been studied by Ayoub [33], developing both displacement and mixed-based finite element formulation capable of capturing the nonlinear behavior of both the beam and the foundation. The beam’s section has been discretized into 16 fibers, while for both the displacement and mixed models two different order of interpolation functions were used, employing a six element discretization. The results were compared with the converged solution obtained by a displacement-based model with fifth order polynomial and a mesh consisting of 32 elements. In Table 3, the maximum deflections of the midpoint of the beam for various internal nodal points’ discretization schemes are presented, illustrating that convergence is achieved for a small number of nodal points. The obtained values are compared for the same load stage Pz = 160 kN (the compared values have been

extracted from a graph), with the one presented in Ayoub [33]. In Fig. 8, the load–displacement curve at the beam’s midpoint is presented, as compared with those obtained from a FEM solution [33] demonstrating a very good agreement. More specifically, the obtained curve is almost identical with the results of the mixed model with cubic moment function and the converged solution, which is assumed to describe the exact behavior of the beam foundation system, while the other models present a perceptible amount of error in the inelastic region. In Figs. 9 and 10 the bending moment My(x) and the displacement w(x) curves along the beam length are presented as compared with those from [33] for the load stage producing deflection at the midspan equal to 5. As it can be observed, the corresponding curves of the present study capture the exact behavior rather accurately and agree with the converged solution and the results of the mixed model with cubic moment interpolation function, while differ from the curves of the displacement model with cubic displacement function. Finally, in Table 4 the midpoint curvature w00 (l/2) and the foundation reaction ps at x = l/6 are also presented and compared for the same load stage (the compared values have been extracted from a graph), verifying that the maximum values of the curvature and the soil reaction are accurately represented. 4.3. Example 3 In order to demonstrate the range of applications of the developed method, in the third example a rectangular cross-section (h = 0.60 m, b = 0.30 m) beam of length l = 6.0 m, clamped at both ends has been studied, employing 20 linear longitudinal elements, 400 boundary elements, 72 quadrilateral cells (12 fibers) and a 3  3 Gauss integration scheme for each cell (cross-sectional discretization). Two material cases have been analyzed, namely an elastic–perfectly plastic with E = 32318.4 MPa, rY0 = 20 MN/m2 and Et = 0 and an elastoplastic-strain hardening with Et = 650 MPa, while two load cases have been examined namely a concentrated load at 2l/3 and a uniformly distributed one, both monotonically increasing. In Figs. 11 and 12 the load–displacement curves at the load position are presented for the two load cases, as compared with a FEM solution [43] obtained by employing 60 beam elements and a 3-D FEM solution [43] obtained by employing 8250 solid (brick)

398

E.J. Sapountzakis, A.E. Kampitsis / Engineering Structures 48 (2013) 389–401

Table 5 Normal stress distribution along the beam length for different load stages for concentrated (a) or uniformly distributed (b) load.

elements. In the case of elastic-perfectly plastic material an additional curve is presented in Fig. 11 applying the Step by Step concentrated plasticity method. A good agreement between the results of the present method and the 3-D FEM solution is observed, especially in the elastoplastic-strain hardening case. Moreover, in Table 5 the normal stress distribution along the beam length is presented for different load stages for both load cases (P(2l/3) = 610 kN, 780 kN and 810 kN corresponding to the occurrence of the plastic hinges).

As a variant of this example, the same beam subjected to the same load cases resting on a tensionless Pasternak foundation is examined. The nonlinear load–displacement foundation reaction is characterized by the perfectly plastic Winkler part with initial stiffness kw = 20 MPa and yielding forcePwY = 60 kN/m and the Pasternak part with stiffness kp = 5000 kN. In Figs. 13 and 14 the load– displacement curves are presented for different types of beam and soil material properties and for both load cases, verifying the significant influence of the inelastic analysis to the soil-beam system

399

E.J. Sapountzakis, A.E. Kampitsis / Engineering Structures 48 (2013) 389–401

600

1600

500

Load Pz (kN) at x=3m

Load Pz (2l/3)

1200

800

Elastic Solution Perfectly Plastic - Elastic Winkler Perfectly Plastic - Plastic Winkler Strain Hardening - Elastic Winkler Strain Hardening - Plastic Winkler

400

400

300

50

150 250 σxx(MPa)

200

Present Study - Perfectly Plastic Present Study - Strain Hardening FEM Shell Model [43] - Perfectly Plastic FEM Shell Model [43] - Strain Hardening

100

0

0 0

0.004

0.008

0.012

0

0.016

0.02

0.04

0.06

0.08

0.1

Displacement w (m) at x=3m

Displacement w (2l/3) Fig. 13. Load–displacement curve at the loading point of the beam of example 3 resting on nonlinear foundation.

Fig. 15. Load–displacement curve at the loading point of the beam of example 4.

120 700

100

80

Load (kN/m)

Load p z (kN/m)

500

300

Elastic Solution Perfectly Plastic - Elastic Winkler Perfectly Plastic - Plastic Winkler Strain Hardening - Elastic Winkler Strain Hardening - Plastic Winkler 0

0.01

0.02

150

250

σxx(MPa)

Present Study - Perfectly Plastic Present Study - Strain Hardening FEM Shell Model [43] - Perfectly Plastic FEM Shell Model [43] - Strain Hardening

20

0 0

0.02

0.04

0.06

0.08

0.1

Displacement w (l/2)

0.03

Displacement w (l/2)

Fig. 16. Load–displacement curve at the midpoint of the beam of example 4.

Fig. 14. Load–displacement curve at the midpoint of the beam of example 3 resting on nonlinear foundation.

Table 6 Deflection (mm) of the midpoint of the beam of example 3, subjected to uniformly distributed loading (load step 2 kN/m).

100 250 300 600

50 40

100

w(l/2)

60

for both types of foundation modeling for various load stages taking into account or ignoring the beam’s material strain hardening slope. From this figure and table, it is easily concluded that the inelastic analysis and the soil nonlinearity are of paramount importance.

Strain hardening Et = 650 MN/m2

Perfectly plastic Et = 0 Elastic Winkler

Elastic Pasternak

Plastic Pasternak

Elastic Winkler

Elastic Pasternak

Plastic Pasternak

4.4. Example 4

1.49 4.31 6.16 24.01

1.46 4.16 5.94 22.30

1.46 4.78 10.40 –

1.49 4.15 5.63 20.39

1.46 4.03 5.46 19.08

1.46 4.39 7.10 –

As a final numerical application, an I-shaped cross-section (total height h = 0.3 m, total width b = 0.3 m, flange width tf = 0.02 m, web width tw = 0.01 m) fixed-pinned beam (E = 213,400 MPa, rY0 = 285 MPa) of length l = 8 m resting on an elastic-plastic Winkler foundation (kw = 25 MPa, PwY = 100 kN/m, kwt = 1.25 MPa) has been studied, employing 32 linear longitudinal elements, 400 boundary elements, 43 quadrilateral cells (15 fibers) and a 3  3 Gauss integration scheme for each cell (cross-sectional discretization). The beam is subjected to either a concentrated load at

response and the importance to the deflections of the subgrade modeling. Finally, in Table 6 the deflection of the midpoint of the beam w(l/2) is presented in the case of a uniformly distributed load

400

E.J. Sapountzakis, A.E. Kampitsis / Engineering Structures 48 (2013) 389–401

From these figures, the significant influence of the inelastic analysis to the soil-beam system response and the accuracy of the proposed formulation are verified.

1400

5. Concluding remarks

Load Pz (kN) at x=3m

1000

In this paper a boundary element method is developed for the inelastic analysis of beams of simply or multiply connected constant cross-section resting on two-parameter tensionless elastoplastic foundation. The main conclusions that can be drawn from this investigation are

600

Present Study

FEM Shell Model [43]

Elastic Solution Perfectly Plastic - Plastic Winkler Perfectly Plastic - Hardening Winkler Strain Hardening - Plastic Winkler Strain Hardening - Hardening Winkler

200

0

0.04

0.08

P P P P 0.12

Displacement w (m) at x=3m Fig. 17. Load–displacement curve at the loading point of the beam of example 4 resting on nonlinear foundation.

350

Load p z (kN/m)

250

150

a. The numerical technique presented in this investigation is well suited for computer aided analysis of prismatic beams of arbitrary simply or multiply connected cross-section having at least one axis of symmetry, supported by the most general boundary conditions and subjected to the action of arbitrarily distributed or concentrated vertical loading. b. Accurate results are obtained using a relatively small number of nodal points across the longitudinal axis. c. A small number of cells (fibers) is required in order to achieve satisfactory convergence. d. The influence of both the inelastic and the tensionless character of the foundation is confirmed. More specifically, the unilateral character of the foundation leads to nonlinear behavior of the soil-beam system even if both the beam and the springs are considered to be elastic, while in all the examined cases the foundation’s perfectly plastic behavior leads to significantly lower collapse load (load– displacement curves) compared to those employing elastic and strain hardening constitutive laws. e. The developed procedure retains most of the advantages of a BEM solution even though domain discretization is required.

Acknowledgements Present Study

FEM Shell Model [43]

Elastic Solution

P P P P

Perfectly Plastic - Plastic Winkler

Perfectly Plastic - Hardening Winkler Strain Hardening - Plastic Winkler

50

Strain Hardening - Hardening Winkler 0

0.04

0.08

0.12

Displacement w (l/2) Fig. 18. Load–displacement curve at the midpoint of the beam of example 4 resting on nonlinear foundation.

position x = 3 m from the fixed end or to a uniformly distributed one, both monotonically increasing. In Figs. 15 and 16 the load–displacement curves are presented for different types of beam material properties ignoring the foundation reaction and for both load cases, as compared with a FEM solution [43] obtained by employing 2882 quadrilateral shell elements. Excellent agreement between the results is observed. In the same figures the normal stress rxx distribution is also presented for several inelastic load levels illustrating the spread of plasticity along the cross-section. Finally, in Figs. 17 and 18 the load–displacement curves are presented for different types of beam and soil material properties and for both load cases, as compared with a FEM solution [43] obtained by employing 2882 quadrilateral shell elements for the I-shaped beam and assuming 161 nonlinear springs following the elastic-plastic law given above.

The work of this paper was conducted from the ‘‘DARE’’ project, financially supported by a European Research Council (ERC) Advanced Grant under the ‘‘Ideas’’ Programme in Support of Frontier Research [Grant Agreement 228254]. References [1] Nukala P, White D. A mixed finite element for three-dimensional nonlinear analysis of steel frames. Comput Methods Appl Mech Eng 2004;193:2507–45. [2] Teh L, Clarke M. Plastic-zone analysis of 3D steel frames using beam elements. J Struct Eng 1999;125:1328–37. [3] Saritas A, Filippou FC. Frame element for metallic shear-yielding members under cyclic loading. J Struct Eng 2009;135:1115–23. [4] Attalla MR, Deierlein GG, McGuire W. Spread of plasticity: quasi-plastic-hinge approach. J Struct Eng 1994;120:2451–73. [5] Orbison JG, McGuire W, Abel JF. Yield surface applications in nonlinear steel frame analysis. Comput Methods Appl Mech Eng 1982;33:557–73. [6] Ngo-Huu C, Kim S, Oh J. Nonlinear analysis of space steel frames using fiber plastic hinge concept. Eng Struct 2007;29:649–57. [7] Winkler E. Theory of elasticity and strength. Czechoslovakia: Dominicus Prague; 1867. [8] Filonenko-Borodich MM. Some approximate theories of elastic foundation. Uchenyie Zapiski Moskovskogo Gosudarstuennogo, Universiteta Mekhanika, Moscow, vol. 46; 1940. p. 3–18 [in Russian]. [9] Hetenyi M. Beams on elastic foundations. Ann Arbo: University of Michigan Press; 1946. [10] Pasternak PL. Fundamentals of a new method of analyzing structures on an elastic foundation by means of two foundation moduli. Gosudarstvennoe Izdatelstro Liberaturi po Stroitelstvui Arkhitekture, Moscow; 1954 [in Russian]. [11] Vlasov L. UN beams, plates, and shells on elastic foundation. Israel Program for Scientific Translations, Jerusalem; 1966. [12] DasGupta S. Axially constrained beams on elastic foundation. Int J Mech Sci 1974;16:305–10.

E.J. Sapountzakis, A.E. Kampitsis / Engineering Structures 48 (2013) 389–401 [13] Zhaohua F, Cook RD. Beam elements on two parameter elastic foundations. J Eng Mech ASCE 1983;109(6):1390–402. [14] Chiwanga M, Valsangkar AJ. Generalized beam element on two-parameter elastic foundation. J Struct Eng 1988;114(6):1414–27. [15] Shirm LM, Giger MW. Timoshenko beam element resting on two-parameter elastic foundation. J Eng Mech 1990;118(2):280–95. [16] Onu G. Shear effect in beam finite element on two-parameter elastic foundation. J Struct Eng 2000;126:1104–7. [17] Kim NI, Shin DK. Stiffness matrices of thin-walled composite beam with monosymmetric i-and channel-sections on two-parameter elastic foundation. Arch Appl Mech 2010;80(7):747–70. [18] Niyogi AK. Bending of axially constrained beams on elastic foundation. Int J Mech Sci 1973;15:781–7. [19] Avramidis IE, Morfidis K. Bending of beams on three-parameter elastic foundation. Int J Solids Struct 2006;43:357–75. [20] Weitsman Y. A tensionless contact between a beam and an elastic half-space. Int J Eng Sci 1972;10:73–81. [21] Kaschiev MS, Mikhajlov K. A beam resting on a tensionless Winkler foundation. Comput Struct 1995;55(2):261–4. [22] Zhang Y, Murphy KD. Response of a finite beam in contact with a tensionless foundation under symmetric and asymmetric loading. Int J Solids Struct 2004;41:6745–58. [23] Silveira RAM, Pereira WLA, Goncalves PB. Nonlinear analysis of structural elements under unilateral contact constraints by a Ritz type approach. Int J Solids Struct 2008;45:2629–50. [24] Zhang Y. Tensionless contact of a finite beam resting on Reissner foundation. Int J Mech Sci 2008;50:1035–41. [25] Ma X, Butterworth JW, Clifton GC. Static analysis of an infinite beam resting on a tensionless Pasternak foundation. Eur J Mech A/Solids 2009;28:697–703. [26] Sapountzakis EJ, Kampitsis AE. Nonlinear dynamic analysis of Timoshenko beam-columns partially supported on tensionless Winkler foundation. Comput Struct 2010;88:1206–19. [27] Sharma S, DasGupta S. The bending problem of axially constrained beams on elastic foundations. Int J Solids Struct 1975;11:853–9.

401

[28] Beaufait FW, Hoadley PW. Analysis of elastic beams on non-linear foundations. Comput Struct 1980;12:669–76. [29] Yankelevsky DZ, Eisenberger M, Adin MA. Analysis of beams on nonlinear Winkler foundation. Comput Struct 1989;31(2):287–92. [30] Kaliszky S, Logo J. Analysis of nonlinear beams on nonlinear foundation by the use of mixed extremum principles. Theocaris PS, Kounadis AN, editors. National Technical University of Athens; 1994. p. 65–80. [31] Sapountzakis EJ, Kampitsis AE. Nonlinear analysis of shear deformable beamcolumns partially supported on tensionless three-parameter foundation. Arch Appl Mech 2011;81:1833–51. [32] Sapountzakis EJ, Kampitsis AE. Nonlinear response of shear deformable beams on tensionless nonlinear viscoelastic foundation under moving loads. J Sound Vib 2011;330:5410–26. [33] Ayoub A. Mixed formulation of nonlinear beam on foundation elements. Comput Struct 2003;81:411–21. [34] Mullapudi R, Ayoub A. Nonlinear finite element modeling of beams on twoparameter foundations. Comput Geotech 2010;37:334–42. [35] Ortiz M, Simo J. Analysis of a new class of integration algorithms for elastoplastic constitutive relations. Int J Numer Methods Eng 1986;23:353–66. [36] Katsikadelis JT. Boundary elements: theory and applications. AmsterdamLondon, United Kingdom: Elsevier; 2002. [37] Vlasov V. Thin-walled elastic beams. Israel Program for Scientific Translations, Jerusalem; 1963. [38] Armenakas AE. Advanced mechanics of materials and applied elasticity. New York: Taylor & Francis Group; 2006. [39] Crisfield MA. Non-linear finite element analysis of solids and structures. Essentials, vol. 1. New York, USA: John Wiley and Sons; 1991. [40] Sapountzakis EJ. Solution of non-uniform torsion of bars by an integral equation method. Comput Struct 2000;77:659–67. [41] Simo JC, Hughes TJR. Computational inelasticity. New York: Springer; 1998. [42] EA de Souza Neto, Peric D, Owen DRJ. Computational methods for plasticity theory and applications. Wiley; 2008. [43] NX Nastran User’s Guide, Siemens PLM Software Inc.; 2007.