Geometrically Nonlinear Inelastic Analysis of Timoshenko Beams on Inelastic Foundation

A.E.Kampitsis and E.J.Sapountzakis

1 Introduction

In design of civil engineering structures(e.g.bridges,wind-turbines,offshore platforms,etc.)the analysis of beam–foundation systems is often encountered.In order to conduct precise analysis,without jeopardizing accuracy and thus safety,the thorough understanding of the mechanics of the beam–foundation system is required.Currently,these systems are designed to behave elastically for every type of loading[EC(2004)],however recent research efforts[Gerolymo et al.(2007);Chiou et al.(2012)]have investigated the beneficial character of permitting plastification to occur at the beam-foundation system.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,yet significant,strength reserves that are mobilized in redundant members after inelastic redistribution takes place.Thus,the cost-effective design of infrastructures requires the realistic estimation of the beam–foundation system response,accounting for all sources of nonlinearities;namely nonlinear stress–strain behavior of the structural member and the soil(material nonlinearity)along with the geometrical nonlinearity.Moreover,the contemporary advancements in material science have facilitated the intensive use of materials having relatively high transverse shear modulus;thereby the error incurred from the ignorance of shear deformation effect may be substantial,particularly in the case of heavy lateral loading.

Over the years,the beam-foundation interaction has been an area of extensive research activity and various methods have been developed in order to study the arising uncertainties[Silva et al.(2010)]and the complex behavior of the system,from the material level to the interaction between structural and foundation elements.These methods can be grouped into three major categories;namely the limit equilibrium[Broms(1964a,b);Broms(1965)],the beam-on-Winkler-foundation[Winkler(1867),Filonenko-Borodich(1940),Hetenyi(1946),Pasternak(1954),Vlasov(1966)]and those based on the continuum mechanics.Among them,the most commonly employed in engineering practice is the beam approach due to the significant advantages over the other methods,such as the simplicity in formulation and modeling together with the high level of accuracy with minor computational cost.

Within this framework,several researches have employed the concept ofelastic beam on nonlinear foundation. In this formulation,the foundation loaddisplacement relation is assumed to follow a nonlinear law while the beam remains elastic throughout the analysis.Sharma and Dasgupta[Sharma and Das-Gupta(1975)]employed an iteration method using Green’s functions for the analysis of uniformly loaded axially constrained hinged beams assuming an exponential load-displacement foundation reaction law.Beaufait and Hoadley(1980)approximated the nonlinear load-displacement relationship of the Winkler 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.(1989)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.El Naggar and Novak(1996)used a Winkler model employing a hyperbolic stress-strain relationship to evaluate the lateral response of piles,while Wang et al.(1998)employed the same method to predict results of centrifuge model tests of single piles in a soft clay soil profile.Lately,Sapountzakis and Kampitsis(2011)studied the nonlinear static analysis of shear deformable beam-columns partially supported on tensionless three parameter foundation,undergoing moderate large deflections under general boundary conditions.

Although the nonlinear behavior of the soil due to high strain level has been studied extensively[Brown and Shie(1991);Laman et al.(1999);Kim and Jeong(2011)]only few studies have encountered the inelastic behavior of both the beam and the foundation elements.According to this,the beam stress-strain and the foundation load-displacement relations are assumed to follow nonlinear inelastic constitutive laws.Consequently,such models are not easily formulated due to the complexity of the problem.To start with,Budek et al.(2000)investigated the inelastic response of a reinforced concrete pile in cohesionless soil while,Ayoub(2003)presented an inelastic finite element formulation 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.Mullapudi and Ayoub(2010)expanded this 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 conditions for the soil medium.Recently,Sapountzakis and Kampitsis(2013)presented a BE method for the inelastic analysis of Euler–Bernoulli beam resting on two-parameter tensionless elastoplastic foundation,illustrating the important influence of both the inelastic and the tensionless character of the foundation.

In this paper a Boundary Element Method is developed for the geometrically nonlinear inelastic analysis of Timoshenko beams of arbitrary doubly symmetric simply or multiply connected constant cross-section,resting on inelastic tensionless Winkler foundation.The beam is subjected to the combined action of arbitrarily distributed or concentrated transverse loading and bending moments in both directions as well as to axial loading,while its edges are subjected to the most general boundary conditions.To account for shear deformations,the concept of shear deformation coefficients is used.A displacement based formulation is developed and inelastic redistribution is modeled through a distributed plasticity(fiber)approach exploiting three-dimensional material constitutive laws and numerical integration over the cross-sections.An incremental–iterative solution strategy along with an efficient iterative process are employed[Ortiz and Simo(1986)],while the arising boundary value problem is solved employing the boundary element method[Katsikadelis(2002)].The essential features and novel aspects of the present formulation compared with previous ones are summarized as follows.

1.The proposed beam model accounts for the geometrical nonlinearity by retaining the square of the slope in the strain–displacement relations,avoiding in this way the inaccuracies arising from a linearized second-order analysis.For that purpose the total Lagrange formulation(intermediate non-linear theory)has been adopted.

2.Shear deformation effect is taken into account on the geometrically nonlinear inelastic analysis of beams on nonlinear foundation(explicit axial-shearflexure interaction)

3.The formulation presented adopts a J2 three-dimensional plasticity law(von Mises)to assess the inelastic beam-foundation system response.

4.The formulation is a displacement based one taking into account inelastic redistribution along the beam axis.

5.A distributed plasticity(fiber)approach has been employed,which has been acknowledged in the literature[Teh and Clarke(1999);Nukala and White(2004);Saritas and Filippou(2009)]to capture more rigorously material nonlinearities than cross-sectional stress resultant approaches[Attalla et al.(1994)]or lumped plasticity idealizations[Orbison et al.(1982);Ngo-Huu(2007)].

6.The inelasticity of the soil medium is taken into account,employing an inelastic Winkler foundation model.

7.The tensionless character of the foundation is also taken into consideration.

8.An incremental-iterative solution strategy is adopted to restore global equilibrium of the system.

9.The shear deformation coefficients are evaluated using an energy approach,instead of Timoshenko’s[Timoshenko and Goodier(1984)]and Cowper’s[Cowper(1966)]definitions,for which several authors[Schramm et al.(1994);(1997)]have pointed out that one obtains unsatisfactory results or definitions given by other researchers[Stephen(1980);Hutchinson(2001)]for which these factors take negative values.

10.The beam is supported by the most general nonlinear boundary conditions.

11.The use of BEM permits the effective computation of derivatives of the field functions(e.g.stresses,stress resultants)which is very important during the nonlinear inelastic response of beam-foundation systems.

12.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.

Numericalexamplesares worked outconfirming the accuracy and the computational efficiency of the proposed beam formulation,as well as the significant influence of the geometrical nonlinearity and the shear deformation effect in the response of a beam-foundations system.

2 Statement of the Problem

2.1　Displacements,Strains,Stresses

Let us consider a prismatic beam of lengthl(Fig.1)with an arbitrarily shaped doubly symmetric constant cross-section,occupying the two-dimensional multiply connected region ? of they,zplane bounded by the Γj(j=1,2,...,K)boundary curves,which are piecewise smooth,i.e.they may have a finite number of corners.In Fig.1,Cyzisthe principalbending 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 elasticityE,shear modulusG,post-yield modulus of elasticityEt,yield stressσY0,and yield strainεY0.The beam is partially supported on inelastic tensionless Winkler type soil of initial stiffnessesky,kzyielding loadsPYy,PYz,and hardening modulikyt,kztaccording toyandzaxes,respectively.The tensionless character of the foundation reaction is taken into consideration through the constitutive law of the nonlinear Winkler springs,prohibiting negative values.The beam is subjected to the combined action of the arbitrarily distributed or concentrated axial loadingpx=px(x),transverse loadingpy=py(x),pz=pz(x)and bending momentsmy=my(x),mz=mz(x)acting alongy,zdirections,respectively(Fig.1).

Under the action of the aforementioned loading,the displacement field of the beam taking into account shear deformation effect is given as

Employing the strain-displacement relations of the three-dimensional elasticity for moderate large displacements[Ramm and Hofmann(1995);Rothertand Gensichen(1987)],the following strain components can be easily obtained

Figure 1:x-z plane of a prismatic beam with an arbitrarily shaped doubly symmetric constant cross-section resting on inelastic foundation under axial-flexural loading.

Considering strains to be small,employing the work conjugate second Piola–Kirchhoff stress tensor[Crisfield(1991)],assuming an isotropic and homogeneous material without exhibiting any damage during its plastification and neglecting theSyy,Szz,Syzcomponents,the stress rates are defined in terms of the strain ones as

whered(·)denotes infinitesimal incremental quantities over time(rates),the ofE?(E?≈E)in beam formulations[Vlasov(1963);Armenakas(2006)].This last consideration has been followed throughout the paper,while any other reasonable expression ofE?could also be used without any difficulty in many beam formulations.

As long as the material remains elastic or elastic unloading occurs

the stress rates are given with respect to the total strain ones following the Hooke’s law(Eqn.(4)),while when plastic flow occurs the stress rates are given with respect to the total and plastic strain ones through Eqns.(4)and(6)as

where the superscriptpldenotes the plastic part of the strain component.The von Mises yield criterion(J2 plasticity),an associated flow rule and an isotropic hardening rule for the material are considered[Crisfield(1991)],permitting the determination of the plastic strain components.The yield condition is described with the expression

Using the aforementioned relation linking the yield stress rate and the proportionality factor,Eqns.(3),(5)-(7)and exploiting the plastic loading condition(d f=0),the stress rates-total strain rates relations are resolved as

whereDelplis the elastoplastic constitutive matrix with

By settingh=0 in the above relations,the constitutive matrix presented by Baba and Kajita(1982)is obtained,while if one of the shear stress components(along with the corresponding strain one)is dropped out,the constitutive relations presented by Chen and Trahair(1992)are also precisely recovered.

Figure 2:Normal stress-strain(a)and yield stress-equivalent plastic strain(b)relationships.

2.2　Equations of Global Equilibrium and Boundary Conditions

To establish the global equilibrium equations and the boundary conditions of the beam-foundation system,the principle of virtual work under a Total Lagrangian formulation neglecting body forces is employed

whereδ(·)denotes virtual quantities,Wintis the stain energy of the beam due to normal and shear stress andWextis the external load work,defined as

whereVis the volume andlis the length of the beam in the undeformed configuration,psy,pszare the foundation reaction according toyandzaxes,respectively,whileNb,Vby,Vbz,MbyandMbzare the externally applied forces and moments at the beam boundaries.In this framework,the stress resultants of the beam are defined as

whereN,Qy,Qzcorrespond to the axial and shear forces andMy,Mzcorrespond to the bending moments according toyandzaxes,respectively.Subsequently,substituting the expressions of the stress components given from Eqn.(7)and exploiting the strain-displacement relations(3),the stress resultants are obtained as

where(′)denotes differentiation with respect tox,Npl,Qplz,Qply,MplzandMplyare the plastic parts of the corresponding stress resultants,Ais the cross-section area,Iy,Izthe moments of inertia with respect to the principle bending axes andGAy,GAzare its shear rigidities of the Timoshenko’s beam theory,where

are the shear areas with respect toy,zaxes,respectively withκy,κzthe shear correction factors anday,azthe shear deformation coefficients.It is worth here noting thatthese stressresultants referto the directionsofthe infinitesimalelements of the cross-section atits deformed configuration,since they have been defined with respect to the second Piola-Kirchhoff stress tensor[Crisfield(1991)].

After substituting Eqns.(3)and(15)into Eqn.(12)and conducting some algebraic manipulations,the global equilibrium equations of the beam-foundation system are obtained as

Furthermore,the application of the principle of virtual work yields the corresponding boundary conditions as

In Eqns.(19b-e)the total vertical reactionsVby,Vbz,and the total bending momentsMby,Mbzare given as

2.3　Shear Deformation Coefficients

The evaluation of the shear deformation coefficientsay,az,corresponding to the principal coordinate systemCyz,is implemented according to the theory of elasticity.These coefficients are established equating the approximate formula of the shear strain energy per unit length[Stephen(1980)]

and are obtained as[Sapountzakis and Mokos(2005)]

and Θ(y,z),Φ(y,z)are stress functions,which are evaluated from the solution of the following Neumann type boundary value problems[Sapountzakis and Mokos(2005)]

wherenis the outward normal vector to the boundary Γ.In the case of negligible shear deformationsaz=ay=0.It is also worth here noting that the boundary conditions(25b),(26b)have been derived from the physical consideration that the traction vector in the direction of the normal vectornvanishes on the free surface of the beam.

3 Numerical Solution

3.1　Integral Representations for the Axial and Transverse Displacements u,v,w and Rotations θy,θz

According to the precedent analysis,the geometrically nonlinear inelastic problem of Timoshenko beams supported on inelastic tensionless soil,reduces to establishing the axial and transverse displacement componentsu(x),v(x),w(x)as well as the rotations due to bendingθy(x),θz(x)having continuous derivatives up to the second order with respect toxand satisfying the boundary value problem described by the governing differential equation(17)along the beam and the boundary conditions(19)at the beam endsx=0,l.

This boundary value problem is solved applying a Boundary Integral Equation method[Sellountos et al.(2010);Sellountos et al.(2012a,b)]and employing the BEM[Papacharalampopoulos et al.(2010);Katsikadelis(2002)],as this is developed in[Sapountzakis(2000);Sapountzakis and Kampitsis(2012);Sapountzakis and Kampitsis(2013)]for the solution of coupled second order differential equations,after modifying it as follows.The motivation to use this particular technique is justified from the intention to retain the advantages of a BEM solution over a domain approach[Providakis(2000)],while using simple fundamental solutions and avoiding finite differences to the solution of the problem.

whereu?is the fundamental solution given as

withr=x?ξ,x,ξpoints of the beam.SinceEA,GAz,GAy,EIyandEIzare independent ofx,Eqns.(27)can be written as

where the kernels Λj(r)= Λj(x,ξ)(j=1,2)are given as

Solving Eqns.(17a-e)with respect toEAu′,GAyv′,GAzw′,EIyθ′′yandEIzθ′′zand substituting the result into Eqns.(29a-e),respectively,the following integral representations are obtained

After carrying out several integrations by parts,Eqns.(31)yield

while by assembling the boundary terms in a more convenient form the integral representations are written as

If shear deformation effects are negligible,thenu5≈u′2andu4≈u′3.In such cases,numericalmethods requiring domain approximation ofunknown quantities,such as FEM,exhibit“l(fā)ocking”effects,when Timoshenko theory is applied to cases where the Euler–Bernoulli theory could also be used[Zienkiewicz and Taylor(2005)].The locking effect phenomenon has been addressed lately by several researchers[Zhu et al.(2010);Cai et al.(2011);Useche(2014);Dong et al.(2014)]and resent developments have been presented.Since domain approximation of unknown quantities is employed in the present numerical technique,locking effects are alleviated by employing the same order of approximation foru4,u5andu′2,u′3.In order to achieve explicit appearance ofu′2,u′3in Eqns.(32b,c),respectively these integral representations are differentiated with respect toξ,yielding

Moreover,noting that plastic parts of the stress resultants depend on the derivatives of the displacement components,it is deduced thatu′1,u′4,u′5must also be computed in order to resolve the total stress resultants(as well as strain components),thus the integral representations(32a,d,e)are differentiated with respect toξ,yielding

Thereafter it is deduced that Eqns.(32d,e),(33a,b)and(34a-c)have been brought into a convenient form to establish a numerical computation of the unknown quantities.Thus,the interval(0,l)is divided intoLelements,on each of which the unknown quantities together with the plastic parts of the stress resultants are 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 7(L+1)algebraic equations is obtained.Six additional algebraic equations are obtained by applying the integral representation(32a-c)at the beam endsξ=0,l,while together the ten boundary conditions(Eqns.(19))yield a linear system of 7L+23 simultaneous algebraic equations

whereK(d)is a generalized elastic(geometrically)nonlinear stiffness matrix,j5i0abt0bis a 7L+23 generalized unknown vector given as

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

3.2　Incremental–Iterative Solution Algorithm

The first step of the incremental-iterative procedure is to determine the external load vector.Thus,load control over the incremental steps is used and load stations are chosen according to load history and convergence requirements.At each load station,the system of nonlinear equations(35)is numerically solved employing an iterative solution strategy.In the framework of this study the modification of Powell’s hybrid algorithm[IMSL,User’s manual(1997)]has been used.This algorithm is a variation of Newton’s method requiring the following quantities.

1.A Jacobian matrix of the system of nonlinear equations which corresponds to the generalized stiffness matrix to the problem at hand.In this study,this matrix is defined explicitly,avoiding this way any possible inaccuracy resulting from the finite differences approximation,while significantly improving the computational time.

2.An initial guess of the solution{dinit}(at each load station).The resolved vectorj5i0abt0bof the previously converged load station is employed in this study{dinit}={dconv}whilej5i0abt0b={0}is used at the first load station.

3.A tolerance parametertolto perform the convergence criterion of the algorithm. In this study this parameter takes values of the rangetol=10?7?10?10.

Thereafter,a number of monitoring cross-sections is defined.It is assumed that the monitoring sections coincide with theL+1 nodal points of the beam interval(Fig.3).

The fiber approach is to be followed for the integrating the section internal forces and moments.Each section is divided into a number of triangular or quadrilateral cells and standard two-dimensional Gauss quadrature rules are employed in each cell to resolve the plastic parts of the stress resultants.If the same number of Gauss points is employed in every cell,thenNdof=Ncells×NGaussholds.Thus,the monitoring stations of each cross-section coincide with the Gauss points of its cells,while exact patch between adjacent cells is not required.

At each load station,the system of nonlinear Eqns.(35)is expressed without explicitly deriving its incremental form which is more extensive due to terms associated with geometrical nonlinearity.This is achieved by exploiting the values of the stressesSxx,Sxy,Sxzthe plastic parts of the strainsεpleq,εplxx,γplxy,γplxzand the kinematic componentsu′1,u′2,u′3,u4,u′4,u5,u′5of the previously converged load station at the current monitoring stations and adhering to the following steps(subscriptcurdenotes the current value of a quantity that is iteratively updated through the algorithm and subscriptconvdenotes the converged value of a quantity from a previous load station).

1.Elastic prediction step:At each monitoring station of the beam,evaluate the trial stress components as

2.Yield criterion:At each monitoring station of the beam the Mises yield criterion is performed,employing Eqn.(8).

3.At each monitoring cross-section of the beam,plastic parts of the stress resultants are evaluate numerically employing the two-dimensional numerical integration scheme.

5.Since convergence is achieved,the parameters are updated and the process described by steps(i)-(iv)is repeated until the foundation convergence criterion is achieved by using a prescribed tolerance oftolfound.

6.The increments of the external loading continue until the goal loading is undertaken or until convergence cannot be satisfied,which means that the last additional increment cannot be undertaken(plastic collapse)

Finally,it is worth noting that the monitoring displacement componentsu,vandwat any interior point of the beam are updated after convergence in each increment by employing the integral representations(32a-c),respectively.A step-by-step algorithmic approach of the nonlinear solution is presented in a flowchart form in Fig.4.

Figure 4:Flowchart of the incremental-iterative solution algorithm.

4 Numerical Examples

On the basis of the analytical and numerical procedures presented in the previous sections,a computer program has been written and representative numerical examples have been studied confirming the accuracy and the computational eff iciency of the proposed beam formulation,as well as the significant influence of the geometrically nonlinear and the shear deformation effects in the response of a beam-foundations system.

Example 1

For comparison purposes,in the first example a cantilever beam of lengthl=2munder concentrated transverse and axial forcesPz(l),Px(l),respectively acting at the tip,has been studied.The beam is made out of aluminum with modulus of elasticityE=69GPa,shear modulusG=26GPaand yielding stressσY0=275MPa,with rectangular cross-section of widthb=0.02m,heighth=0.8mand shear correction factorαz=1.2.The efficiency of the proposed formulation is illustrated through a convergence analysis performed in case of linear elastic response as compared with the exact solution for the tip displacementwexactand rotationθyexactevaluated by the analytical expressions

In Fig.5,the percentage error of the maximum tip displacement and rotation for various internal nodal points’discretization schemes is presented,while in Table 1 the converged values are compared with those obtained from the Reduced Integration Element(RIE)proposed by Reddy(1997).From the obtained results it is concluded that the shear locking has been successfully prevented and satisfactory accuracy is achieved(i.e.error≤1%)with small number of nodal points,while it is noted that in order to achieve adequate accuracy with RIE several elements are required[Reddy(1997),Saritas and Filippou(2009)].

Thereinafter,the geometrically nonlinear inelastic response of the cantilever is investigated taking into account the shear deformation effect(axial-shear-flexural interaction),employing 22 linear longitudinal elements,40 quadrilateral cells and a 2×2 Gauss integration scheme for each cell.The influence of the normalized axial loadingnx=Px/Pulton the nonlinear response of the beam is also investigated.The present example was first studied by Triantafyllou and Koumousis(2011)who presented a hysteretic Timoshenko beam element based on the lumped plasticity assumption,accounting for the interaction between axial,shear and bending,implementing the yield criterion proposed by Simo et al.(1984).

Figure 5:Tip displacement and rotation error for different internal domain discretization schemes.

Table 1:Deflection(m)and rotation(rad)of the tip of the cantilever of Example 1.

In Fig.6,the load-displacement curves at the cantilever’s tip are presented for two axial load cases;namely zero axial force andnx=0.9.The results obtained with the proposed formulation are compared with those from Triantafyllou and Koumousis(2011)and from a 3-D FEM solution[NX Nastran(2007)]by employing 640 solid(brick)elements.Excellent agreement between the results is observed in case of zero axial load while very good convergence is achieved fornx=0.9.Moreover,the ultimate load predicted from the proposed formulation for zero axial force

Finally,in Figs.7a,b the von Mises stress distribution along the cantilever’s length is presented for different load stages showing the spread of plasticity,while in Fig.7c,d the normaland shearstresspro fi le along the cross-section atx=0.1mfrom the fi xed end,are presented assuming either constant or a more accurate parabolic shear stress distribution as presented in[Saritas and Filippou(2009)].From these fi gures,the flexural character of the plasti fi cation becomes apparent while it is evident that the in fluence of the shear stress pro fi le is negligible,in this example.

Figure 6:Load–displacement curves at the tip of the cantilever beam of Example 1.

Example 2

Table 2:Geometric Properties of the I-shaped cross-section,of Example 2.

In Figs.8,9 the load-displacement curves are presented performing either,geometric and material nonlinear(GMNL)analysis or material nonlinear(MNL)analysis ignoring the foundation reaction,for both the boundary condition cases.The results are compared with those obtained from a FEM model[NX Nastran(2007)]implemented by employing 2400 quadrilateral shell elements.Excellent convergence between the results is observed.In the same figures the von Mises stressσvMdistribution is also presented illustrating the plastification of the wed,as well as the non-symmetry of the normal stresses due to the developed axial force.Additionally,the flexure-only response ispresented in these figures.Since the beam yieldsin shear,the Euler-Bernoulli model fails to capture the nonlinear response and overestimates the collapse load for both the clamped and the fixed-pinned boundary conditions.

The main reason for that divergence is the inability of the flexure-only model to predict the exact collapse mechanism,as it ignores the development of the shear stresses.In more detail,Fig.10a,b depicts the stress distribution along the length of the web for geometrically nonlinear and linear analysis,respectively indicating the shearcharacterofthe collapse mechanism.In the same figure the corresponding deformed shell FEM contour representations are also presented verifying the accuracy of the presented model.On the contrary,Fig 10c show the von Mises stress distribution assuming a flexure-only model demonstrating the collapse mechanism due to bending,which require the formation of three plastic hinges instead of two in the axial-shear–flexure coupling model.

Moreover,under the scope of efficiency,it is worth noting that even thought the two approaches have fundamental differences(i.e.22 elements for the proposed

model instead of 2400 elements for the shell one),the difference between the computational time required for the analyses is significant,while the obtained result have the same accuracy.Indicatively,it is mentioned that the refined shell model required approximately 30 min to 1.0h depending to the analysis type and model parameters,while the proposed one required from 10sec to 240sec for the same type of analysis.

Figure 7:von Mises stress(MPa)distribution along length(a)and normal&shear stress(MPa)distribution along cross-section(b)for different load stages.

Finally,in Figs.11,12 the load-displacement curves of the beam-foundation system are presented,performing either geometrically nonlinear or linear inelastic analysis for both cases of boundary conditions,making evident the influence of the geometrical nonlinearity to the response of the system.Additionally,the flexureonly response is presented in these figures,illustrating once again the importance of the shear deformation effect on the behavior of the beam-foundation system.

Example 3

The geometrically linear case with absence of foundation reaction has been studied by Papachristidis et al.(2010),who proposed a force-based(FB)3D fiber beam element formulation accounting for the axial–shear–moment interaction.In Fig.14 the load-displacement curve at the beam’s midpoint is presented as compared with those obtained from[Papachristidis et al.(2010)]assuming both force and displacement based(DB)formulations for numerous integration sections and numerical integration schemes.

The accuracy and efficiency of the proposed formulation are confirmed by the excellent agreement between the converged solution of Papachristidis et al.(2010)obtained by2 FBelements with 8 integration sections and the one obtained from the conducted analysis assuming the same number of integration sections(i.e.16).More specifically,from this figure it is concluded that the conventional displacement based elements of equal length fail to capture accurately the collapse load.This shortcoming can be resolved by employing either more dense mesh or adap-tively spaced elements.Contrary to the conventionalDBelements,theFBare capable of describing the inelastic response of the beam with a single element per member.However the results may differ with respect to the number of integration sections and the numerical integration scheme(i.e.Gauss(G)and Gauss–Lobatto(G-L)).

Figure 8:Load–displacement curve at the midpoint of the clamped beam of Example 2.

Figure 9:Load–displacement curve at the midpoint of the fixed-pinned beam of Example 2.

Figure 10:von Mises stress distribution contour diagrams along the length of the web for geometrically nonlinear(a)&linear(b)analysis as compared with the shell FEM model[NX Nastran(2007)].Flexure-only model(c).

Figure 11:Load–displacement curve at the midpoint of the clamped beamfoundation system of Example 2.

Figure 12:Load–displacement curve at the midpoint of the fixed-pinned beamfoundation system of Example 2.

Figure 13:Inelastic beam-foundation system subjected to monotonically increasing concentrated load.

Figure 14:Load–displacement curve at the midpoint of the beam of Example 3,performing geometrically linear inelastic analysis.

In Fig 15 the load-displacement curves are presented,performing either geometrically nonlinear or linear inelastic analysis for different types of material properties ignoring the foundation reaction.From this figure,it is concluded that large displacements,influence significantly the behavior of the beam since the developed restoring force does not allow the evolvement of the plastic hinges and thus the plastic collapse.This can also be evident from the contrast observed between the von Mises stress distribution contour diagram as presented in Fig.16a,b performing either geometrically nonlinear or linear inelastic analysis for perfectly plastic and strain hardening material,respectively.

Finally,in Fig 17 the load-displacement curves of the beam-foundation system are presented,performing either geometrically nonlinear or linear inelastic analysis for different types of material properties,while in Table 3 the extreme values of the von Mises stresses for the all the conducted analyses are shown illustrating once again the paramount importance of both geometrical and material nonlinearity in the beam-foundation system analysis.

Figure 15:Load–displacement curve at the midpoint of the beam of Example 3.

Figure 16:von Mises stress(MPa)distribution along the beam length,of Example 3.

Table 3:Extreme values ofthe von Mises stress σvM ofthe beam-foundation system of Example 3.

Figure 17:Load–displacement curve at the midpoint of the beam of Example 3,resting on nonlinear foundation.

5 Concluding Remarks

In this paper a Boundary Element Method is developed for the geometrically nonlinear inelastic analysis of Timoshenko beams resting on inelastic tensionless Winkler foundation.To account for shear deformations,the concept of shear deformation coefficients is used.A displacement based formulation is developed and inelastic redistribution is modeled through a distributed plasticity(fiber)approach.An incremental–iterative solution strategy along with an efficient iterative process are employed.The main conclusions that can be drawn from this investigation are

1.The proposed beam formulation is capable of yielding results of high accuracy,as verified by comparing with analytical and FEM results,with minimum computational cost,providing a simple and efficient computational tool for the geometrically nonlinear inelastic analysis of beam-foundation systems.

2.The beam character of the developed formulation confers advantages over more refined approaches in the sense of modelling effort,model handling,results interpretation and isolation of structural phenomena.

3.The proposed model accurately captures both the initial yielding and the ultimate(collapse)load.

4.The influence ofgeometricalnonlinearity isillustrated through the significant discrepancy between the results of the linear and the nonlinear analyses.

5.The proposed model takes into account coupling effects of bending,shear and axial deformations,illustrating the paramount importance of this interaction in the inelastic analysis either under small or large displacement formulation.

6.The significant influence of the inelastic character of the foundation is also demonstrated.

7.A small number of cells(fibers)is required in order to achieve satisfactory convergence.

8.The developed procedure retains most of the advantages of a BEM solution over a FEM approach,requiring a small number of nodal points to achieve high accuracy.

9.The use of BEM enables the accurate calculation of the stress resultants which are very important during both the analysis and the design of beamfoundation systems.

Acknowledgement:The authors would like to gratefully acknowledge the f inancial support provided by the“Bodossakis Foundation”.

Armenakas,A.E.(2006):Advanced mechanics of materials and applied elasticity.Taylor&Francis Group:New York.

Attalla,M.R.;Deierlein,G.G.;McGuire,W.(1994):Spread ofplasticity:quasiplastic-hinge approach.Journal Struct.Engrg.,vol.120,pp.2451-2473.

Ayoub,A.(2003):Mixed formulation of nonlinear beam on foundation elements.Computers and Structures,vol.81,pp.411–421.

Baba,S.;Kajita,T.(1982):Plastic analysis of torsion of a prismatic beam.International Journal for Numerical Methods in Engineering,vol.18,pp.927-944.

Beaufait,F.W.;Hoadley,P.W.(1980):Analysis of elastic beams on non-linear foundations.Computer&Structures,vol.12,pp.669-676.

Broms,B.(1964a):Lateral resistance of piles in cohesive soils.Journal of Soil Mechanics and Foundation Division ASCE,vol.90,no.2,pp.27-63.

Broms,B.(1964b):Lateral resistance of piles in cohesionless soils,Journal of Soil Mechanics and Foundation Division ASCE,vol.90,no.3,pp.123-156.

Broms,B.(1965):Design of laterally loaded piles.Journal of Soil Mechanics and Foundation Division ASCE,vol.91,no.3,pp.77-99.

Brown,D.A.;Shie,C.F.(1991):Some numerical experiments with a threedimensional finite element model of a laterally loaded pile.Computers and Geotechnics,vol.12,pp.149-162.

Budek,A.M.;Priestley,M.J.N.;Benzoni,G.(2000):Inelastic seismic response of bridge drilled-shaft RC pile/columns.Journal of Structural Engineering,vol.126,no.4,pp.510-517.

Cai,Y.C.;Tian,L.G.;Atluri,S.N.(2011):A simple locking-free discrete shear triangular plate element.CMES:Computer Modeling in Engineering&Sciences,vol.77,no.4,pp.221-238.

Chen,G.;Trahair,N.S.(1992):Inelastic nonuniform torsion of steel I-beams.Journal of Constructional Steel Research,vol.23,pp.189-207.

Chiou,J.S.;Tsai,Y.C.;Chen,C.H.(2012):Investigating influencing factors on the ductility capacity of a fixed-head reinforced concrete pile in homogeneous clay.Journal of Mechanics,vol.28,no.3,pp.489-498.

Cowper,G.R.(1966):The shear coefficient in Timoshenko’s beam theory.Journal of Applied Mechanics ASME,vol.33,no.2,pp.335-340.

Crisfield,M.A.(1991):Non-linear finite element analysis of solids and structures.vol.1 Essentials.John Wiley and Sons;New York:USA.

Dong,L.;El-Gizawy,A.S.;Juhany,K.A.;Atluri,S.N.(2014):A simple locking-alleviated 4-node mixed-collocation finite element with over-integration,for homogeneous or functionally-graded or thick-section laminated composite beams.CMC-Computers,Materials&Continua,vol.40,no.1,pp.49-77.

El Naggar,M.H.;Novak,M.(1996):Nonlinear analysis for dynamic lateral pile response.Journal of Soil Dynamics and Earthquake Engineering,vol.14,no.3,pp.141-157.

En 1998-1(2004):Eurocode8—Design of structure for earthquake resistance—Part 2:Bridges,Part 5:Foundations,retaining structures and geotechnical aspects.European Committee for Standardization.

Filonenko-Borodich,M.M.(1940):Some approximate theories of elastic foundation.Uchenyie Zapiski Moskovskogo Gosudarstuennogo.Moscow:Universiteta Mekhanika,vol.46,pp.3–18.

Gerolymos,N.;Drosos,V.;Gazetas,G.(2007):Seismic response of singlecolumn bent on pile:Evidence of beneficial role of pile and soil inelasticity.BEE,pp.547-573.

Hetenyi,M.(1946):Beams on elastic foundations.University of Michigan Press:Ann Arbo.

Hutchinson,J.R.(2001):Shearcoefficients forTimoshenko beam theory.Journal of Applied Mechanics ASME,vol.68,pp.87-92.

IMSL(1997):User’s manual,IMSL MATH/Library:Fortran subroutinesformathematical applications(1997)Visual Numerics Inc.,Houston,USA.

Katsikadelis,J.T.(2002):Boundary Elements:Theory and Applications.Amsterdam-London,United Kingdom:Elsevier.

Kim,Y.;Jeong,S.(2011):Analysis of soil resistance on laterally loaded piles based on 3D soil–pile interaction.Computers and Geotechnics,vol.38,pp.248–257.

Laman,M.;King,G.J.W.;Dickin,E.A.(1999):Three-dimensional finite element studies of the moment-carrying capacity of short pier foundations in cohesionless soil.Computers and Geotechnics,vol.25,pp.141-155.

Lubliner,J.(2008):Plasticity Theory.Dover,New York.

Mullapudi,R.;Ayoub,A.(2010):Nonlinear Finite element modeling of beams on two-parameter foundations.Computers and Geotechnics,vol.37,pp.334–342.Ngo-Huu,C.;Kim,S.;Oh,J.(2007):Nonlinear analysis of space steel frames using fiber plastic hinge concept.Engineering Structures,vol.29,pp.649-657.

Nukala,P.;White,D.(2004):A mixed finite element for three-dimensional nonlinear analysis of steel frames.Computer Methods in Applied Mechanics and Engineering,vol.193,pp.2507-2545.

NX Nastran(2007):User’s Guide,Siemens PLM Software Inc.

Orbison,J.G.;McGuire,W.;Abel,J.F.(1982):Yield surface applications in nonlinear steel frame analysis.Computer Methods in Applied Mechanics and Engineering,vol.33,pp.557-573.

Ortiz,M.;Simo,J.(1986):Analysis of a new class of integration algorithms for elastoplastic constitutive relations.International Journal for Numerical Methods in Engineering23,353-366.

Papacharalampopoulos,A.;Karlis,G.F.;Charalambopoulos,A.;Polyzos,D.(2010):BEM solutions for 2D and 3D dynamic problems in mindlin’s strain gradient theory of elasticity.CMES-Computer Modeling in Engineering&Sciences,vol.58,no.1,pp.45-73.

Papachristidis,A.;Fragiadakis,M.;Papadrakakis,M.(2010):A 3D fibre beam-column element with shear modeling for the inelastic analysis of steel structures.Comput Mech.,vol.45,pp.553–572.

Pasternak,P.L.(1954):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.

Providakis,C.P.(2000):BEM/FEM Comparison Studies for the Inelastic Dynamic Analysis of Thick Plates on Elastic Foundation.CMES-Computer Modeling in Engineering&Sciences,vol.1,no.3,pp.123-130.

Ramm,E.;Hofmann,T.J.(1995):Stabtragwerke,Der Ingenieurbau,Ed.G.Mehlhorn,Band Baustatik/Baudynamik,Ernst&Sohn,Berlin.

Reddy,J.N.(1997):On locking-free sheardeformable beam finite elements.Comput Methods Appl Mech Eng.,vol.149,no.1-4,pp.113-132.

Rothert,H.;Gensichen,V.(1987):Nichtlineare Stabstatik,Springer-Verlag,Berlin.

Sapountzakis,E.J.(2000):Solution of non-uniform torsion of bars by an integral equation method.Computers and Structures,vol.77,pp.659-667.

Sapountzakis EJ;Kampitsis AE(2013)Inelastic analysis of beams on twoparameter tensionless elastoplastic foundation.Engineering Structures48,389–401.

Sapountzakis,E.J.;Kampitsis,A.E.(2012):A BEM Approach for Inelastic Analysis of Beam-Foundation Systems under Cyclic Loading.CMES-Computer Modeling in Engineering and Sciences,vol.87,pp.97-124.

Sapountzakis,E.J.;Kampitsis,A.E.(2011):Nonlinear analysis of shear deformable beam-columns partially supported on tensionless three-parameter foundation.Arch Appl Mech.,vol.81,pp.1833-1851.

Sapountzakis,E.J.;Mokos,V.G.(2005):A BEM solution to transverse shear loading of beams.Computational Mechanics,vol.36,pp.384-397.

Saritas,A.;Filippou,F.(2009):Inelastic axial-flexure–shear coupling in a mixed formulation beam finite element.International Journal of Non-Linear Mechanics,vol.44,pp.913-922.

Saritas,A.;Filippou,F.C.(2009):Frame element for metallic shear-yielding members under cyclic loading.Journal Struct.Engrg.,vol.135,pp.1115-1123.

Saritas,A.;Filippou,F.C.(2009):Inelastic axial-flexure–shear coupling in a mixed formulation beam finite element.International Journal of Non-Linear Mechanics,vol.44,no.8,pp.913-922.

Schramm,U.;Kitis,L.;Kang,W.;Pilkey,W.D.(1994):On the shear deformation coefficient in beam theory.Finite Elements in Analysis and Design,vol.16,pp.141-162.

Schramm,U.;Rubenchik,V.;Pilkey,W.D.(1997):Beam stiffness matrix based on the elasticity equations.International Journal for Numerical Methods in Engineering,vol.40,pp.211-232.

Sellountos,E.J.;Sequeira,A.;Polyzos,D.(2010):Solving Elastic Problems with Local Boundary Integral Equations(LBIE)and Radial Basis Functions(RBF)Cells.CMES-Computer Modelling in Engineering and Sciences.,vol.57,no.2,pp.109-135.

Sellountos,E.J.;Tsinopoulos,S.V.;Polyzos,D.(2012a):A LBIE method for solving gradient elastostatic problems.CMES-Computer Modeling in Engineering and Sciences,vol.86,no.2,pp.145-170.

Sellountos,E.J.;Polyzos,D.;Atluri,S.N.(2012b):A new and simple meshless LBIE-RBF numerical scheme in linear elasticity.CMES-Computer Modelling in Engineering and Sciences.,vol.89,no.6,pp.513-551.

Sharma,S.;DasGupta,S.(1975):The bending problem of axially constrained beams on elastic foundations.Int.J.Solids Structures,vol.11,pp.853-9.

Silva,C.R.A.;Azikri de Deus,H.P.;Mantovani,G.E.;Beck,A.T.(2010):Galerkin Solution of Stochastic Beam Bending on Winkler Foundations.CMESComputer Modeling in Engineering&Sciences,vol.67,no.2,pp.119-150.

Simo,J.C.;Hjelmstad,K.D.;Taylor,R.L.(1984):Numerical formulations of elasto-viscoplastic response of beams accounting for the effect of shear.Comput Methods Appl Mech Eng.,vol.42,no.3,pp.301–330.

Stephen,N.G.(1980):Timoshenko’s shear coefficient from a beam subjected to gravity loading.Journal of Applied Mechanics ASME,vol.47,pp.121-127.

Stephen,N.G.(1980):Timoshenko’s shear coefficient from a beam subjected to gravity loading.Journal of Applied Mechanics ASME,vol.47,pp.121-127.

Teh,L.;Clarke,M.(1999):Plastic-zone analysis of 3D steel frames using beam elements.Journal of Structural Engineering,vol.125,pp.1328-1337.

Timoshenko,S.P.;Goodier,J.N.(1984):Theory of Elasticity.3rd edn.New York McGraw-Hill.

Triantafyllou,S.;Koumousis,V.(2011):An inelastic Timoshenko beam element with axial–shear–flexural interaction.Comput Mech.,vol.48,pp.713–727.

Useche,J.(2014):Boundary Element Analysis of shear deformable shallow shells under harmonic excitation.CMES-Computer Modeling in Engineering&Sciences,vol.1,no.1,pp.1-14.

Vlasov,L.(1966):UN beams,plates,and shells on elastic foundation.Jerusalem:Israel Program for Scientific Translations.

Vlasov,V.(1963):Thin-walled elastic beams.Israel Program for Scientific Translations:Jerusalem.

Wang,S.;Kutter,B.L.;Chacko,J.M.;Wilson,D.W.;Boulanger,R.W.;Abghari,A.(1998):Nonlinear seismic soil-pile-superstructure interaction.Earthquake Spectra,vol.14,no.2,pp.377-396.

Winkler,E.(1867):Theory of elasticity and strength.Czechoslovakia:Dominicus Prague.

Yankelevsky,D.Z.;Eisenberger,M.;Adin,M.A.(1989):Analysis of beams on nonlinear Winkler foundation.Computer and Structures,vol.31,no.2,pp.287-292.

Zhu H.H.,Cai Y.C.,Paik J.K.and S.N.Atluri(2010)Locking-free thick-thin rod/beam elementbased on a von Karman type nonlineartheory in rotated reference frames for large deformation analyses of space-frame structures.CMES-Computer Modeling in Engineering&Sciences,57(2),175-204.

Zienkiewicz,O.C.;Taylor,R.L.(2005):The finite element method for solid and structural mechanics.Elsevier,Amsterdam.