Oscillations of a beam on a non-linear elastic foundation under periodic loads

The complexity of the response of a beam resting on a nonlinear elastic foundation makes the design of this structural element rather challenging. Particularly because, apparently, there is no algebraic relation for its load bearing capacity as a function of the problem parameters. Such an algebraic relation would be desirable for design purposes. Our aim is to obtain this relation explicitly. Initially, a mathematical model of a flexible beam resting on a non-linear elastic foundation is presented, and its non-linear vibrations and instabilities are investigated using several numerical methods. At a second stage, a parametric study is carried out, using analytical and semi-analytical perturbation methods. So, the influence of the various physical and geometrical parameters of the mathematical model on the non-linear response of the beam is evaluated, in particular, the relation between the natural frequency and the vibration amplitude and the first period doubling and saddle-node bifurcations. These two instability phenomena are the two basic mechanisms associated with the loss of stability of the beam. Finally Melnikov’s method is used to determine an algebraic expression for the boundary that separates a safe from an unsafe region in the force parameters space. It is shown that this can be used as a basis for a reliable engineering design criterion.


Introduction
Beams on an elastic foundation, or columns and piles supported along their length, typically by soils, are a very old issue, and well known topic in structural mechanics, besides being a very common structural element, with applications in many engineering fields such as mechanical, civil and offshore engineering, in particular, in foundation analysis and design.The study of this structural configuration started with the work of pioneers such as Winkler [1] and has been treated by numerous researchers with diverse theoretical tools, including numerical [2], finite element [3][4][5], analytical [1,6,7] and perturbation methods [8].
The linear behavior of beams on elastic foundations has been extensively studied.Little attention, however, has been given to their behavior in the nonlinear range [3].In many important applications, the elastic foundation is the soil.As an elastic foundation soils are typically highly nonlinear.This nonlinear effect greatly influences the beam behavior by changing the beam's load bearing capacity, and natural frequencies.
In spite of its practical and theoretical importance, studies of the dynamics of beams, or columns, on a nonlinear elastic foundation under dynamic loads are rare.In the present paper the nonlinear dynamics of an axially loaded beam-foundation system is analyzed, considering both the nonlinearity of the beam and foundation.The foundation's non-linearity is of the softening type and has a marked influence on the vibrations of the structural system.The periodic load combined with the elastic foundation's non-linearity makes the beam oscillations rather complex, with supercritical and sub-critical bifurcations, cascade period doubling bifurcations, leading to chaos, and saddle-node bifurcations.As a final outcome, the beam, under certain excitation frequencies, may display a load bearing capacity of only a fraction of its normal capacity when compared to static or step loads.On the other hand, at some frequencies the beam load capacity can be many times superior to its static critical load.
The complexity of the beam's response to a periodic load makes the design of this structural element particularly challenging because there is, apparently, no algebraic relation of the load bearing capacity as a function of the geometrical and elastic parameter of the beam and foundation in design codes or technical literature.The aim of the present work is to study in detail the non-linear dynamic behavior and instabilities of the beam under harmonic forcing and, from the understanding of the overall behavior of the beam, derive upper and lower bound estimates of the dynamic buckling load as a function of the beam, foundation and load parameters

Large deflection beam formulation
To analyze the beam, a nonlinear large-deflection formulation is used together with the Ritz method to obtain a simplified equation of motion, which is in turn solved by perturbation procedures and numerical integration.
Consider the prismatic beam of Fig. 1 with length L and bending stiffness EI, simply supported and loaded by the axial force P, which retains its magnitude and direction as the beam deflects, and a time dependent transversal load q(t).
Additionally consider that the beam lies on an elastic foundation such that the deflection at a particular point is independent of the foundations behavior at all other points.For soils, that assumption, of course, is not strictly true.However, it has been found from experiments that, for the pattern of lateral deflections which can occur in practice, the soil reaction at a point is dependent essentially on the beam lateral deflection at that point and independent of the beam deflections above or below that point [9,10].Thus, for purposes of analysis, the foundation can be replaced by a set of discrete nonlinear soil springs with the same load-deflection characteristics.This relation can be approximated by various mathematical formulas.The most frequently used mathematical functions in soil analysis are parabolas, hyperbolas, splines, and the Ramberg-Osgood functions.Here the foundation nonlinear restoration force is given by the three parameter Ramberg-Osgood function where E ti is the initial tangent modulus, E tf is the final tangent modulus and P u is the ultimate soil resistance.Approximating Eq. ( 1) by a third order odd Taylor expansion one obtains: which can be written in a more concise form as where The transversal uniformly distributed, time dependant periodic load q(t) is given by As the beam deflects each point of the center line positioned at a distance x from the left corner moves to an unspecified horizontal position and a vertical position.Considering that the center line is inextensional, the arc length from the left corner is equal to x, and the deflected form of the beam is totally specified by the transversal displacement w(x) where x ranges from 0 to L.
The curvature χ is, by definition, the rate of change of angle with arc length, so The strain energy stored in a beam element of length dx is dU beam = M χ dx/2, where M is the bending moment given by M = EIχ.So, the total strain energy stored in the beam is Additionally the strain energy stored in a foundation element of length dx is And the total strain energy stored in the foundation is The potential energy, L P , of the conservative load P is given by where e is the end-shortening of the column.Let m be the mass per unit length of the beam, then the kinetic energy, T , can be expressed as: The variation of the work of the transversal and damping forces can be written as The deformed shape of a simply-supported beam under small to moderately large deflections in the vicinity of the lowest natural frequency can be described fairly well by the one term approximation up to moderately large deflections [11,12]: that satisfies all boundary conditions.Substituting (13) into Eqs ( 7) to (12), using the Ritz method to discretize the beam in space and applying Hamilton's principle one obtains the following differential equation of motion: which can be rewritten in a more concise form as where: Equation ( 16) is the well known Duffing equation.

A practical example
Consider a circular cylindrical steel pile with thickness t = 1 cm; external diameter D = 10 cm, length L = 20 m, Young modulus E steel = 210 GPA, mass density, m = 20 kg/m and bending stiffness EI = 5672067 N m 2 .Considering a stiff clay, the parameters used for the Ramberg-Osgood constitutive law are [9]: initial tangent modulus −E ti = 4037 Pa; final tangent modulus −E tf =0; shape parameter −n = 1 and ultimate soil resistance Based on these values, one obtains the following reference parameters: -Lowest natural frequency -17 rad/s; -Euler buckling load −P euler = (π/L) 2 E I = 139953 N For this particular case, the differential equation of motion Eq. ( 16) takes the form where a damping parameter of about 1% of the critical damping value is adopted (μ = 5,34 kg/m/s).
Considering a linear foundation model, Eq. ( 16) reduces to

Non-linear vibration and escape
Special cases of Eq. ( 16) have been extensively studied [13][14][15][16][17][18][19][20][21][22][23][24] and it is known that for the laterally unloaded beam (F = 0) and for axial loads smaller than the buckling load (ω 0 >0), the beam on a softening foundation has a trivial stable equilibrium configuration.This trivial stable equilibrium sits between two non trivial unstable configurations (saddle points).The heteroclinic orbit connecting these two unstable solutions is the boundary of the stable solutions.This is illustrated in Fig. 2. That is, initial conditions within this region lead to bounded solutions (centers).Here we call this region the "conservative basin of attraction" (Fig. 2).For initial conditions outside this region, the solution tends to infinity, indicating that the beam looses its stability or the deflections become larger than the beam can accept, leading as a rule to permanent damage or even the failure of the structure.We call this phenomenon escape from the safe pre-buckling well.Buckling will also occur when the load level makes the basin of attraction too small or zero in size (ω 0 = 0).

Influence of the foundation nonlinearity on the beam's vibration and buckling load
Considering only the static terms in Eqs ( 15) and ( 16) one obtains for the post-buckling path of the column the following equation Now, considering the free vibration model, the frequency-amplitude relation is given by As observed in Fig. 3 the degree of non-linearity of the post-buckling path and of the frequency-amplitude response of a beam without foundation or considering a linear foundation model is rather small.In fact the linear foundation only increases the critical load and natural frequency, without affecting, as expected, the non-linear term in Eq. ( 16).On the other hand, if the softening behavior is considered, the post-buckling path and the frequency-amplitude relation become highly non-linear due to the decreasing stiffness of the beam-foundation system.As a result the beam becomes imperfection sensitive.This softening effect is the dominant factor in the local and global non-linear dynamics of the beam, as shown herein.

The dynamical escape mechanisms
Besides the initial position and velocity conditions, the lateral cyclic load, at certain forcing levels and frequencies, can also lead to escape from this potential well.When small levels of this load are applied, the solutions tend to be periodic with the same period of the applied load.For the same load level F the vibration amplitude varies with the forcing frequency, Ω.It tends to be larger close to the beam's natural frequencies and its super and sub-harmonics.Figure 4 shows a typical resonance curve for the beam on a softening foundation.Note that the peak of the curve bends, as expected, to the left due to the softening character of the cubic nonlinearity.
As the load level is increased, the vibration amplitude increases and the overall response tends to become more complex.Finally escape occurs.Different excitation frequencies lead to different escape loads.The variation of the escape load, F e , as a function of the excitation frequency, Ω, is called the escape boundary.Figure 5 shows the escape boundary obtained for the present example considering a slowly evolving environment.The escape boundary is formed by a sequence of ascending and descending curves.The ascending parts are of fractal nature [25,26].Our investigations have shown that, under a periodic load, there are basically two escape mechanisms [26]: a simple saddle-node, or fold, bifurcation when the stable solution, with the same frequency of the excitation, meets an unstable solution and disappears, and a supercritical period doubling, or flip, bifurcation sometimes cascading to chaos.Figure 6 shows an example of the two different escape mechanisms.Case (a) shows that for F = 0 there is only one trivial stable solution.As the load level is increased an unstable solution of the same period appears and becomes increasingly close to the stable solution.When the two, stable and unstable, solutions meet they cancel out and disappear.This means that the system does not have a stable configuration and escapes.One should note that when the two solutions approach each other, the stable solution's basin of attraction decreases and tends to zero.Thus finite dynamical perturbations to the system can lead to escape.This kind of bifurcation occurs at the escape boundary's descending branches.
In case (b) the solution's behavior prior to escape is conceptually different.At a certain load level the solution doubles its period, but remains stable.As the load increases a succession of period doublings occur at an ever decreasing interval.The solution then becomes chaotic and escape occurs.In this case the basin of attraction's shape may not change as the load increases, but becomes more and more complex.In fact our investigations show that it becomes fractal, such that escape configurations sit very close to stable configurations.This type of bifurcation occurs at the ascending branches of the escape boundary.Besides the main solution, which starts from the trivial solution of the unloaded beam, other stable and unstable solutions may co-exist at the same load level.In fact numerical studies of these secondary paths have shown that they are quite complex.But the knowledge of these two basic escape mechanisms, the fold and the flip, will be the basis for the development of the theory that will enable one to predict the beams load carrying capacity in spite of the underlying complex dynamics.

An approximate solution
In order to be sure that the vibration phenomena described above are not restricted to the parameters used in the present numerical simulations, it is important to obtain analytical expressions based on the system parameters.Due to the impossibility of obtaining an exact analytical solution for the non-linear Eq. ( 16), an approximation using the variant of Lindstedt-Poincar é's method presented by Meirovich [27] was obtained.Assuming that the response, and its frequency can be expanded up to second order as a function of a small parameter ξ, one can write: )  Additionally the other system parameters can be expressed as: As observed here, the force amplitude and the damping are considered to be small of the order of ξ.
The substitution of Eqs ( 22) through (27) into Eq.( 16) generates a sequence of ordinary differential equations in W 0 , W 1 and W 2 .Solving the equations one at a time and adjusting the parameters α 1 , α 2 and φ to eliminate the secular terms, the following relations are obtained: The amplitude parameter a can be obtained by solving Eq. ( 28) and Eq.(29).Once this parameter is obtained, the time response of the beam, Eq. ( 30), is totally defined.Additionally making the damping parameter, c, and the external load, f , zero, in Eq. ( 29) one obtains the following approximation for the frequency-amplitude relation This is the expression of the backbone, that is the nonlinear frequency-amplitude relationship, as shown in Fig. 3.Note that due to the softening effect of the non-linear term the natural frequency decreases with the vibration's amplitude, as obtained numerically in Fig. 4. It should be pointed out that the present approximation compares rather well with the exact solution obtained by numerical integration.

The solution's stability
An important element in the non-linear dynamic analysis of the beam is to consider how it looses stability under external excitation.This can be done by assuming a small perturbation x(t) << 1 on the solution W (t). Substituting this in Eq. ( 16), using Eq. ( 30) as the approximate solution and eliminating the higher order terms of x, the stability equation becomes: where Clearly the trivial solution x ≡ 0 is a solution of Eq. (32) for every set of parameters.Additionally if a = 0 the trivial solution is stable and unique.If, for any a * = a > 0, Eq. (32) admits another solution, then there is a bifurcation at this point.If the secondary solution is periodic with the same period as Eq. ( 34) then the bifurcation is a saddle-node (fold) bifurcation.Thus assuming [22] x Hill's Determinant becomes: When Eq. ( 35) becomes zero a fold bifurcation occurs.
On the other hand, if the secondary solution has twice the period of Eq. ( 30) (half the frequency) a period doubling (flip) bifurcation occurs.Thus assuming Hill's determinant becomes Equation ( 35) together with Eq. (37) allows one to verify if the solution given by Eqs ( 28) and ( 29) are stable.These criteria were used to identify the bifurcation events connected with the stability boundaries shown in Fig. 5.

A lower bound estimate
The use of the fold and the flip bifurcations as a prediction of the escape load appears, initially, to be a safe and reliable criterion for design.Actually it is an unsafe prediction as can be seen in Fig. 5, where black dots show the load carrying capacity of the structure obtained numerically by adding to the lateral periodic load a noise with zero mean and standard deviation of 20% of the maximum load F , and a bandwidth of 20% around the applied frequency [26].From the scatter of results obtained by this numerical simulation, it can be seen that the structure's load capacity is always inferior to the predicted load obtained by the bifurcation criterion.It is expected in practice a similar scatter of dynamic buckling loads.Actual buckling loads can be even lower if the deleterious effect of imperfections and deviations in the value of system parameters is considered in the problem formulation.In essence, prior to the fold bifurcation, the stable solution's basin of attraction can be so small that finite small dynamical perturbations can make the beam escape to infinity.Also the basin boundary can become fractal and, in this case,  the long term behavior of the beam becomes unpredictable [26,28,29].These two possible cases are illustrated in Fig. 7. So, based on these arguments, the escape boundary as predicted above is better defined as an upper bound.Thus, a safer load capacity prediction can be obtained by examining the basin of attraction's boundary.For the laterally unloaded beam this basin of attraction's boundary is the stable manifold of the two saddle points that form the heteroclinic orbits around the stable focus.This is an indication for the use of the Melnikov method to determine the load level that produces the first transversal crossing between the stable and unstable manifolds of the saddle points.When this crossing occurs at a certain point, it occurs at an infinite number of points, making the corresponding basin of attraction fractal.
To apply Melnikov's method one starts with an exact expression of the heteroclinic orbit of the conservative part of Eq. ( 18) [30].This is By this expression the Melnikov function that measures the distance between the stable and unstable manifolds of the saddle point can be written as Simple zeros of the Melnikov function imply the transversal crossing of the stable and unstable manifolds.Therefore a necessary condition for the Melnikov function to be zero is that In terms of the system parameters the lower bound is given by Assuming that the basin of attraction starts to be eroded from f M up, one can use this load level as a safe lower bound estimate for the beam's load carrying capacity.This global change in the nature of the basin of attraction is not another escape mechanism, it only signals the beginning of a process that will ultimately lead to escape by a fold or flip bifurcation.Fig. 5 shows the values of f M for the practical example analyzed in this work.The significant difference between the upper and lower bounds in Fig. 5 is due to the softening character and consequent imperfection-sensitivity of the system which may lead to instability loads in practice much lower than the theoretical critical load.Note the relative simplicity of Eq. (41) which can be easily calculated and used as a design criterion.The scatter of experimental loads of imperfection sensitive systems is well documented in literature [31].This leads usually to design loads as low as a small fraction of the critical load in agreement with the theoretical results of the present research.

Conclusions
Based on the study of the nonlinear dynamic behavior and instabilities of the beam under harmonic forcing and, from the understanding of the overall behavior of the beam, upper Eqs ( 25), ( 26), ( 27), (32) and (34) and lower bound Eq. (38) estimates of the dynamic buckling load as a function of the beam, foundation and load parameters were derived.Specifically, the influence of a non-linear elastic foundation on the non-linear dynamic behavior and stability of slender beams has bean analyzed using a simplified model.Experimental results show that this foundation non-linearity for most soils is of the softening type.For practical values of soil parameters the softening term is particularly large, dominating the response: while beams without foundation or resting on a linear elastic foundation exhibit a stable post-buckling path and a frequency-amplitude relation of the hardening type with a rather small degree of non-linearity, columns resting on a non-linear softening foundation exhibit a strong non-linearity and a high imperfection-sensitivity.This has a profound influence on the load carrying capacity of the beam and on the bifurcations connected with the instability boundaries in control space.To take into account the imperfection sensitivity, the use of the Melnikov function as a base for the derivation of a safe lower-bound of the dynamic buckling loads was proposed in the present work.This can be used as a guidance in pile foundation analysis and design.Additionally, the global loss of engineering integrity is discussed and tools for estimating the degree of safety of a given structure are proposed.

Fig. 2 .
Fig. 2. Basin of attraction of the trivial solution for the laterally unloaded beam.

Fig. 3 .Fig. 4 .
Fig. 3. Influence of the foundation on (a) the buckling and post-buckling behavior of the beam and (b) on the frequency-amplitude relation.

Fig. 5 .
Fig. 5.The escape boundary in force control space.

Fig. 6 .
Fig. 6.The two possible escape mechanisms associated with the instability boundaries.

Fig. 7 .
Fig. 7. Possible configurations of the basin of attraction close to the escape load.