A Discrete Model for Nonlinear Vibration of Beams Resting on Various Types of Elastic Foundations

This paper presents a discrete physical model to approach the problem of nonlinear vibrations of beams resting on elastic foundations. The model consists of a beam made of several small bars, evenly spaced. The bending stiffness is modeled by spiral springs, and the Winkler soil stiffness is modeled using linear vertical springs. Concentrated masses, presenting the inertia of the beam, are located at the bar ends. Finally, the nonlinear effect is presented by the axial forces in the bars, assumed to behave as longitudinal springs, due to the change in their length induced by the Pythagorean Theorem. This model has the advantage of simplifying parametric studies, because of its discrete nature, allowing any modification in the mass matrix, the stiffness matrix, and the nonlinearity tensor to be made separately. Therefore, once the model is established, various practical applications may be performed without the need of going through all the formulation again. The study of the nonlinear behavior makes the solution of the movement equation rise in complexity. By considering this discrete model and using the linearization method, one can achieve an idealized approach to this nonlinear problem and obtain quite easily approximate solutions.


Introduction
The analysis of the vibration of beams resting on elastic foundations wears a practical and theoretical interest in many fields such as civil, mechanical, and transportation engineering.Analytical and numerical methods applied to the modeling of such a problem are extensively addressed in the literature.However, the combination between the emergence of high-performance computers and problems complexity made the discrete methods more appealing.
A discrete method such as the Finite Element Method (FEM) is the first to address the problem numerically [1][2][3].The Differential Quadrature Method (DQM) is also employed for the solution to similar engineering problems involving beam vibrations and foundations [4][5][6].Therefore, Malekzadeh and Karami [7] gather the advantages of both previous methods (DQM and FEM) to perform the free vibration and buckling analysis of thick beams on twoparameter elastic foundations.
On the other hand, soil behavior wears a great complexity, because of its nonlinear nature.In that order, literature has investigated multiple soil models starting from the idealized elastic theories modeled by Winkler to the more complex viscoelastic theories of Pasternak and Hetinyi.Similarly, Kerr [8] carries a study showing the evolution of soil theories with their extensions.
The aim of this work is the development of a flexible general discrete model for linear and nonlinear vibrations of beams resting on elastic foundations, via the adaptation and the extension of the approach presented in [9].Accordingly, the process of discretization carried out is intended to allow efficient variation of beam and soil geometrical and mechanical characteristics, in order to make it easy to perform multiple parametric studies.
For nonlinear vibrations, the multimode analysis leads to a coupled nonlinear differential system involving the contribution to various modes.Using the harmonic balance method, Benamar [18] and Benamar et al. [19] reduce the 2 Advances in Acoustics and Vibration nonlinear free vibration problem to a set of nonlinear algebraic equations.Consequently, the objective here is the investigation into linear and nonlinear vibrations, using the discrete physical model developed for various types of soil stiffness, and beam end conditions.
This study establishes a simplification to the nonlinear problem induced by large deflections.However, the nonlinearity tensor obtained becomes heavy in terms of calculation for high values of the discretization parameter  (presenting the number of masses in the physical model in Figure 2).In fact, we addressed this complexity and proposed a simplification of the calculation, leading to an efficient algorithm.
The structure of this paper is as follows: Section 2.1 presents the general formulation, where the continuous and discrete models are detailed, in order to show their theoretical similarities.After the establishment of the discrete model, the following Section 2.2 is concerned with determination of the discrete model parameters and the nondimensional formulation.Finally, Section 2.3 introduces a solution to the multimodal nonlinear problem using a linearization procedure.
The physical discrete model validation is done in Section 3, devoted to details and comments on the results obtained in the cases considered.An examination is made of how far the results deviate from the linear and the nonlinear continuous theories, and a discussion is given about their range of validity under the different assumptions adopted.
The main purpose of this work is achieved in Section 4, in which straight applications of the theory are developed and validated via comparisons with previous references.The first application deals with a partially supported beam on elastic foundations, where the springs presenting the soil reaction are applied to only a limited portion of the beam span.In the second application, the beam examined is resting on a variable elastic foundation, where the distribution of soil spring stiffness is linear or parabolic.

General Theory
A nonlinear discrete model of a beam, made of extensional bars, concentrated masses and rotational springs located at the bar ends, is introduced in [9] to approach the continuous theory using different assumptions.In the present work, a beam resting on an elastic soil is undergoing the same discretization process by extending and completing the previous model.

General Formulation.
This section presents the discrete and continuous models for nonlinear vibrations of beams resting on elastic foundations.Both models are based on the Lagrangian principle, leading to a formally identical description.The linear behavior of the mechanical system examined is determined by the classical mass and rigidity matrices   and   .On the other hand, the nonlinear behavior, presented by nonlinearity tensor   due to the axial forces induced by the large vibration amplitudes in the continuous case, is established by the application of the Pythagorean Theorem in the discrete physical model.

Continuous Model.
Before introducing the extended formulation of the discrete model, a quick review of the theory behind the continuous model [18] is introduced.Let us consider transverse vibrations of a beam, shown in Figure 1, resting on elastic foundations, having the characteristics (, , , , , , and ) defined in Nomenclature.The total beam strain energy can be written as the sum of the strain energy due to bending denoted as    and the strain energy induced by Winkler springs    plus the axial strain energy  nl due to the axial load induced by large deflections.Thus    ,    ,  nl , and the kinetic energy  are as follows [20]: where  is the beam transverse displacement.Using a generalized parameterization and the usual summation convention used in [18], the transverse displacement can be described as [20]  (, ) =   ()   () for ,  = 1, . . ., .
For a conservative system, the dynamic behavior of the structure may be obtained by Lagrange's equations [20]: where  =    +    +  nl .By introducing the energy expressions, that is, (3) to (6), into (8), one gets This nonlinear differential system (9) may be expressed in a matrix form as where [M], [K  ], [K  ], and [B(q)] are, respectively, the mass matrix, the bending rigidity matrix, the Winkler springs rigidity matrix, and the nonlinear rigidity matrix.
According to the conclusion of [9,18,21,22], a harmonic distortion of the nonlinear response of a harmonically exited structure happens at large transverse vibration amplitudes.However, it is shown both experimentally and theoretically that the higher harmonic components remain very small, compared to the response first harmonic component.Consequently, the nonlinear free response of the beam is assumed in the present paper, as in many previous ones, to be drawn as [18,20] where   's are the contribution coefficients to the beam nonlinear mode shapes.By replacing (11) in (10) and applying the harmonic balance method, one has At this point, the vibration of a continuous Euler-Bernoulli beam resting on elastic foundations is put under a matrix form.It is interesting to note that the model presented in this section and summarized in (12) can be considered as a multidimensional form of the Duffing equation, very often encountered in nonlinear vibration analyses of structures with cubic nonlinearities.The solution to this nonlinear equation is discussed in Section 2.3.

Discrete Model.
In order to adapt the discrete model [9] to the case of a beam resting on Winkler elastic foundations and to compare it to the continuous model established in Section 2.1.1,let us consider the -degree-of-freedom (DOF) system shown in Figure 2, made of  concentrated masses  1 , . . .,   ,  + 1 bars presented by longitudinal linear springs,  + 2 torsional springs, and  longitudinal linear soil springs.  is the stiffness coefficient of the rth spiral spring, for  = 1 to  + 2. The moment M in the spiral spring is given by M = −  Δ; Δ =   −  −1 being the angle between the bars adjacent to the node r, considered as longitudinal springs of length   and stiffness    ,  = 1, to .The Winkler foundations are modeled using the longitudinal vertical spring's distribution, with    presenting the stiffness coefficient of the rth linear spring, for  = 1 to .
Using vector notation, {d 1 }; {d 2 }; . . .; {d  }; . . .; {d  } is the displacement basis, denoted as DB, and  1 to   are, respectively, the transverse displacements of the masses  1 to   from the horizontal equilibrium positions ( = 0) of the system corresponding to the nondeformed positions of the springs (spiral and linear).The vectors {d  } are defined by with {d  }  being the unit displacement of the ith mass.The displacement vector can be written in DB as by considering the modal basis {{Φ 1 }, {Φ 2 }, . . ., {Φ  }, . . ., {Φ  }}, obtained by solution of the linear eigenvalue problem and denoted as MB, the displacement vector can be expressed as where {Φ  } is the th linear mode shape of the system and  1 ,  2 , . . .,   , . . .,   are the components of the mass displacements expressed in MB. {Φ  } is denoted in DB as Using the transition matrix [Φ] from DB to MB, the two sets of coordinates in the two bases are related by So, the displacement vector can be written by analogy with the continuous displacement vector   , as [9] for ,  = 1, . . ., , where  discr  is the modulus of the displacement   expressed in DB (or the contribution of the normalized vector {d  } of DB) and   is the modulus of the displacement   expressed in MB (or the contribution of the normalized vector {Φ  } of MB). nl discr is the nonlinear frequency of the discrete system associated with the amplitude  discr  .Using the usual summation convention for the repeated indices , , , and , the kinetic, linear, and nonlinear strain energy expressions, obtained by replacement of the expression (18) of   in the energy expressions (3) to (6) for , , ,  = 1, . . ., , where   ,    ,    , and   are, respectively, the general terms of the mass tensor, the linear rigidity tensors corresponding to the spiral and longitudinal vertical springs, and the nonlinear rigidity tensor in MB.The relationships between the expressions for these tensors in DB and MB can be obtained using the transition matrix [Φ] as follows [19,23]: Replacing in (8) the expressions for ,    ,    , and  nl given in ( 19)- (22) and applying the harmonic balance method, in MB, lead to [24,25] The linearized approximate methods of solution of nonlinear algebraic systems of the type (24) are proposed in [20].In the present paper, to improve the accuracy of the numerical solutions in the range of vibration amplitudes taken into consideration, the so-called second formulation in [20] is combined with the Newton-Raphson method in Section 2.3.

The Expressions for the General Terms of the Tensors 𝑚 𝑖𝑗 , 𝑘 𝑠
,    , and   .One of the major advantages of this model is its ability to explicitly express each tensor and give a clear insight of the parameters chosen, before going through the solution of the dynamic problem.However, in a discrete problem, the calculation time depends on the discretization number .For the linear parameters   and   , increasing the number of masses in the physical discrete model does not heavily affect the calculation time, because filling the matrices and performing the change of bases from the DB to the MB do not require more than a two loops' process.
Nevertheless, for the nonlinear parameter   (4D tensor), it requires only one loop in the filling process but does need eight nested loops for the inversion, in order to transform the tensor from DB to MB.In fact, this limitation made the solution of the nonlinear problem for high values of the discretization parameter  a quite lengthy process.Accordingly, by analyzing the form of the nonlinearity tensor, we noticed a certain number of patterns.Then, we carried a meticulous analysis leading to a simplified procedure as detailed in Appendix A, where the eight nested loops are reduced to only three.
Appendix B gives a summary of the calculation, carried out in [9], for the expression of the parameters   , and the other values of    are obtained by symmetry relations or are equal to zero.The boundaries conditions of the system determine the stiffness    for  = 1 and  + 1.For simply supported ends, the corresponding torsional springs  1 and  +1 stiffness are equal to zero, since the rotation is free; then   = 5( + 1) 3 / 3 for  = 1 and  + 1.For clamped ends, the corresponding torsional spring stiffness is infinite, since the rotation is prevented, so    = +∞ for  = 1 and  + 1.
(iii) The Winkler springs tensor [K  ]/   terms are expressed as (details given in Sections 2.2.1(1) and 2.2.1(2)) for ,  = 1, . . ., . (iv) The nonlinear rigidity tensor [B]/  terms are expressed as (Appendix B) and other values of   are obtained by symmetry relations or are equal to zero.
(1) Calculation of the Equivalent Stiffness for the Winkler Linear Springs.The linear potential energy stored in the Winkler vertical spring    , shown in the discrete model of Figure 3, can be written as [9] On the other hand, the elementary potential energy corresponding to an elementary section of length  in the continuous beam resting on the elastic foundations is given by with  corresponding to the soil stiffness per length unit.By approximation, one can write By substituting (31) in (30), one obtains Since the bars have the same length (30) can be reduced to The identification between ( 29) and (33) gives (2) Expressions for the General Terms of the Tensors    .The relationship existing between the longitudinal elongations Δ  of the ith linear spring and the transverse displacements   is The potential energy stored in the ith linear spring, of stiffness    , is given directly in terms of the stiffness    and the spring length variation   , by the very well-known relationship, that is, 1/2    (  ) 2 , which leads, after summation, to The identification between ( 4) and (30) leads to the following expression for the rigidity tensor associated with the Winkler foundations in DB: for ,  = 1, . . ., , where   is Kronecker's symbol.Equation (37) can be written in a matrix form as with [I] being the identity matrix.

Nondimensionalization.
In order to establish the necessary comparisons with previous publications without further manipulations, the nondimensional formulation of the problem presented is essential.To define the nondimensional parameters, let us put The Winkler stiffness coefficient is given by and   being the nondimensional parameter: In the case of a constant soil distribution   =  and   = .
The nondimensional amplitude  * is expressed as where  is the radius of gyration.
The expression for the nonlinear frequency [19] is As the nondimensional formulation is established, the general expressions for  *  ,  *   ,  *   , and  *  for the discrete model become as follows: (i) The nondimensional mass tensor [M * ]/ *  is obtained by replacing ( 25) in (39): and other values of  *  are equal to zero.

Solution of the Nonlinear Algebraic System: Linearization.
In order to solve the nonlinear amplitude equation and establish a good approximation for large vibration amplitudes of the beam examined, a linearization of the nonlinear algebraic problem is performed, based on the so-called second formulation introduced in [20].For this linearization to be accomplished, only the second-order terms of the type      1   are neglected when considering the first nonlinear mode, in (24), rewritten in summation form as The nonlinear expression        *  is split to terms proportional to  1 3 , terms proportional to Then (52) can be written in a matrix form as with where  [20].
This formulation is proven [20] to give good results when applied to various nonlinear vibration problems.However, to get an even better approximation, it is combined with a Newton-Raphson algorithm, which requires a good initial estimate of the solution (given by the actual linearization method).Accurate results (up to 10 −14 ) are obtained within few numbers of iteration.

Results Validation
In Section 2, a general theory is developed for the physically discrete model, with demonstrations carried out for calculating the new model parameters.In order to confirm the extended theory, this section is devoted to the validation of the results in the linear and nonlinear cases for simply supported and clamped beams.

Simply Supported Beams Resting on Elastic Foundations.
In Figure 2, a discrete physical model, made of masses, bars, and springs, is described.In the following section, a validation of the results obtained for a simply supported beam resting on Winkler elastic foundation is discussed.Namely, calculations are performed in the linear and nonlinear cases using Matlab's software, based on the procedure described in Section 2. , combining the effects of the extensional and spiral springs, obtained by addition of the Winkler soil stiffness matrix and the spiral spring matrix, for the simply supported case, can be presented as follows [9]: With this new expression for the nondimensional rigidity matrix [K SS *  ] (24), corresponding to linear vibrations, can be written as with {y} being the vector of mass displacements in the MB.Solving (58) gives the eigenvalues {}, which allows the corresponding frequencies to be written as Table 1 gives, for the first 3 modes of the simply supported beam, the nondimensional frequencies, for increasing number of masses ( = 10, 50, 100) and increasing values of the Winkler foundation stiffness.Accordingly, these nondimensional frequencies are calculated from (59), by extracting the √/ 4 term, where the eigenvalues   are those of the nondimensional stiffness matrix [K SS *  ].So, they can be presented as In order to validate the model introduced in Section 2, for the simply supported beam, a comparison of the first 3 natural frequencies is carried out for increasing values of  (the number of masses) and  (the soil stiffness), and results are drawn in Table 1.
The natural frequencies, Table 1, converge faster by increasing the  or .A good approximation is already obtained for ( = 10,  =0) where the overall error ((average − present study ())/average %) is below three percent limit.On the other hand, by increasing the soil stiffness to the case ( = 10,  = 1000000) the error drops below 0.15%.However, once the discretization number  is higher than 50 the overall error is under the 0.12%.Thus, increasing one of the parameters  and , or both, makes the discrete model converge to the reference models [10][11][12].

Nonlinear Case
(1) Mode Shapes.In [9], the discussions covered only the validation of nonlinear frequencies.In the present work, the validation is extended to the beam mode shapes.Likewise, Table 2 gives a comparison between the mode shapes calculated using the present discrete physical model with previous results [13].
(2) Comparison of Nondimension Frequency  ss * nl discr / ss *  discr .For validating the results in the nonlinear case, comparing the ratio of the nonlinear frequency parameter ( ss * nl discr /  ss *  discr ) (54) and (60) may be one of the best indicators of the nonlinearity type and acuity.The number of masses considered in the calculation is ( = 100), because of the good approximation obtained for this degree of discretization in the linear case (error below 0.1%).
Table 3 gives the nonlinear frequencies for increasing values of the vibration amplitude and the soil parameter .These numerical results give a comparison of the beam nonlinear behavior by the mean of frequency ratio ( ss * nl discr / ss *  discr ).
For small values of the nondimensional vibration amplitude (<1), the results are close to those of [14] (error ≈ 1%).However, once the nondimensional amplitude rises to values higher than one, the difference becomes more noticeable making the error exceed seven percent.This late conclusion may be due to the hypothesis taken for the nonlinearity effect.Another aspect that one may extract from this comparison is that the soil stiffness affects the convergence between the discrete and the continuous models [14], since the error decreases by increasing the soil stiffness coefficient.

Clamped Beam Resting on Elastic
Foundations.This section is concerned with the validation of the results obtained for the clamped beam case.Also, calculations of the frequencies and mode shapes are performed for the linear and nonlinear cases using Matlab's software, based on the procedure described in Section 2.

Linear Case.
As it is stated in [9], to establish the clamped beam stiffness matrix, an infinite coefficient should be taken for the first and last nodes.In that order the following matrix is built: However, once introduced in the calculation, a mathematical error about the conditioning of the matrix arises (illconditioned matrix).
Since the present physical model considers that the bars are rigid transversally, it may be concluded that, by blocking (nonrotation) the first bar, the second node located at its end would be blocked too at  = 0 as shown in Figure 4. Using this method for taking into account the clamped ends conditions, the N-DOF clamped beam problem is reduced to a ( − 2)-DOF problem by removing the first and last row and column of the stiffness matrix; then (61) becomes With this new expression for the nondimensional rigidity matrix [K CC −2 * ], the linear vibration equation ( 24) can be described by with {y} being the vector of mass displacements in MB.Solving (63) gives the eigenvalues {} so one can write the frequencies as Table 4 gives the first three modes of the clamped beam nondimensional frequencies for increasing number of masses ( = 10, 50, 100) and increasing values of the Winkler foundation stiffness.These nondimensional frequencies are calculated from (64) by extracting the √/ 4 term, where the eigenvalues   are those of the nondimensional stiffness matrix [K CC −2 * ], so they can be presented as In order to validate the discrete model, for clamped beams, Table 4 shows a comparison between the results obtained here for the first three modes of a clamped beam, for various values of the soil stiffness coefficients and increasing number of masses N = (10/50/100), and the theoretical values of the natural frequencies of continuous beams obtained in [12].

Nonlinear Case
(1) Mode Shapes.In [9], the discussions covered only the validation of frequencies.In order to have a more complete validation, Table 5 compares the present work mode shapes with previous works.For clamped beams, the mode shapes are shown both theoretically and experimentally to be amplitude dependent at large vibration amplitudes.Table 5 shows that the discrete model leads to amplitude dependent mode shapes.By comparing these mode shapes to the ones obtained in [13,15,16], it can be seen that these results complies with the results obtained by the two methods FEM [13] and GFEM [16].
(2) Comparison of Nondimension Frequency  cc * nl discr /  cc *  discr .As it is explained in the linear case Section 3.2.1, the same procedure for implementing the beam end conditions is applied for the nonlinear rigidity tensor, which is also subjected to dimension reduction by removing the first and last column and row of the 4D tensor.Once this manipulation done, the results of ( cc * nl discr / cc *  discr ) given by ( 54) and (65) for various values of the vibration amplitude and soil stiffness are summarized in Table 6.The number of masses taken is ( = 100), for which a good convergence is reached in the linear case (error below 1%).
The nonlinear behavior of a clamped beam resting on elastic foundations is presented in Table 6 by means of the amplitude dependence of the parameter ( cc * nl discr / cc *  discr ).For small values of nondimension amplitude (<1), the results are very close to [14] (error <1%).However, once the nondimensional amplitude rises (>1), the difference becomes more significant.The last conclusion may be due as it is stated Ref. [14] Error %  = 100 Ref. [14] Error %  = 100 Ref. [14] Error % 0,2 Figure 5: A partially supported continuous S-S beam.
for the simply supported case to the model taken for the nonlinearity effect.It may also be noticed that the effect of the soil stiffness on the nonlinear frequency affects the error between the discrete and the continuous model where the error decreases for the same nondimensional amplitude for high values of the soil stiffness parameter .

Applications
One of the main purposes in developing the discrete physical model presented here is to allow easy variation of the different parameters, compared with the continuous models, especially in the nonlinear case for which analytical solutions may be very laborious and, quite often, impossible to obtain.Since the key factor in the present work is the soil stiffness, two applications of an important practical interest are examined.The first application is concerned with a beam supported partially on elastic foundations.The second is concerned with a beam supported on a Winkler elastic foundation with a variable stiffness.In both cases, simply supported and clamped beams are examined.

A Partially Supported Beam.
In order to get a good approximation in the nonlinear case, the number  of masses used in the discretization process is taken equal to 100, taking the same algorithm developed according to the steps discussed in Section 2 and varying the soil stiffness parameter  values and the repartition inside the algorithm between the supported span (  ̸ = 0) and the nonsupported span (  = 0).

A Simply Supported Beam (Linear Case).
In this section, a simply supported beam is assumed to be subjected to the effect of partial intermediate supports and comparison is made between the results of the continuous model of Figure 5(a) and the discrete model of Figure 5(b).For a practical convenience, a semilogarithmic curve is drawn, showing the relative stiffness parameter  on the abscissa and the nondimensional frequency parameter (√ ss *  discr ) on the ordinate.The ratio of the supported span to the total span  [26] chosen covers the cases where the beam is supported along 25, 50, and 75 percent of its span; the extreme cases are also taken where  = 0, no-intermediate supports, and  = 1, fully supported.The results correspond to the case  = 100, which also enabled a simplification of  ratios (25% corresponds to 25 nodes supported and so on).
In order to establish a comparison with the results given in [26], the first five linear modes NDF() curves are drawn.Since the results extracted from [26] consisted only of curves, a graphical comparison is conducted showing a very good approximation.The conclusions drawn from Figure 6 are that the first modes are more sensitive to the soil stiffness variation.For the first mode the influence on the frequency is noticeable at an early stage ( = 10).However, the 5th mode seems to be less sensitive to the soil stiffness variation since the frequency jump is noticeable only for  = 1000.

A Simply Supported Beam (Nonlinear Case).
After validating the physical model in the linear regime, one may investigate the nonlinearity effect on the NDF()(√ SS *  discr ).It may be concluded in Figure 7(b), corresponding to a high nondimensional amplitude ( * = 3), that the difference between the curves is more accentuated than in Figure 7(a) ( * = 1), corresponding to a lower nondimensional amplitude and, consequently, a less accentuated nonlinear effect.By deepening the analysis, the gaps between the curves seem to take different paths: In the case of  = 1, the gap is very small, but by taking a zoom, the gap is more pronounced for low soil stiffness and tends to become narrower as the soil stiffness rises.In the case of  = 0.75, contrary to the case  = 1, the gap between the curves seems to widen by increasing the soil stiffness.In the case of  = 0.5, the same conclusion may be made as in the case  = 0.75 but the gap seems to remain constant.For  = 0.25, the same behavior is obtained as for  = 0.5.In the case of  = 0, the two curves are parallel.

A Clamped Beam (Linear Case).
In this section, the case of partially supported beam with clamped ends is investigated by comparing the results of the continuous model of Figure 8(a) to the discrete model of Figure 8(b).
The study of partially supported beams in [26] covered only the simply supported and free ends cases.In [27] and based on the work of [26], the clamped beam case is investigated.
In [27], a study is conducted using continuous beam models and results are given for the 4 first linear modes shapes.By graphically comparing the curves NDF()(√ CC *  discr ) of Figure 9 to [27], a very good approximation is reached by the present physically discrete model.The same conclusion is obtained for the simply supported beam in terms of frequency variation.

A Clamped Beam (Nonlinear Case).
After establishing and validating the physical model in the linear regime in Section 4.1.3,a comparison is carried out here for the nonlinear case.To show the impact of nonlinearity, the NDF()(√ CC *  discr ) curves in the linear case ( * = 0) and in the nonlinear case ( * ̸ = 0) are drawn in the same graph Figure 10.
In Figure 10(b), which corresponds to a high value of the nondimensional amplitude ( * = 3), the difference between the curves is more expressed in Figure 10(a) ( * = 1), corresponding to a lower nondimensional amplitude and, consequently, a less accentuated nonlinear effect.By deepening the analysis, the gaps between the curves seem to take different paths: in the case of  = 1, the gap is very small, but by taking a zoom, the gap is more pronounced for low soil stiffness and tends to become narrower as the soil stiffness rises.In the case of  = 0.75, contrary to the case  = 1, the gap between the curves seems to widen by increasing the soil stiffness.In the case of  = 0.5, the same conclusion may be made as in the case  = 0.75 but the gap seems to remain constant.For  = 0.25, the same behavior is obtained as for  = 0.5.In the case of  = 0, the two curves are parallel.However, the gaps between the curves are narrower than in the simply supported case.

Variable Elastic Foundation.
Foundations are often not as homogenous as one may expect, especially for long beams such as pipes and piers, because soil is made of multiple layers with a variable stiffness.In order to approach this reality, the foundation may be sometimes modeled with a distribution of stiffness with linear (66) or parabolic (67) variations.In [17], this problem is solved using the Differential Transform Method (DTM) and in the following Sections 4.2.1 and 4.2.2, comparison of frequencies is made between the discrete model introduced in this paper and the DTM.
Equations ( 66) and (67) are introduced in the calculation algorithm and results are obtained.

A Simply Supported Beam (Linear Case).
Simplification of variable implementation is the aim of the present discrete model.Figure 11 shows a simply supported beam resting on a variable elastic foundation and Table 7 gives the variation of the frequency parameter Ω  associated with the ith mode shape, for  = 1 to 8, in accordance with the soil variation for the simply supported beam.
The results obtained for the simply supported beam case, drawn in Table 7 for a linear soil distribution and in Table 8 for a parabolic distribution, for the first eight modes (frequency parameter Ω  ), show that the soil distribution has a low impact on the convergence speed.On the other hand, the discrete model converges well to the results of [17] for a low value of the number of masses since the error is already below one percent for  = 10.

A Simply Supported Beam (Nonlinear Case).
Figure 12 shows the backbone curve of a simply supported beam  resting on a variable elastic foundation.In Figure 12(a), corresponding to a low soil stiffness ( = 10), the nonlinear behavior is more expressed since the nonlinear frequency is higher than the one of the high soil stiffness Figure 12(b) ( = 1000).However, for this last case the curves are more distinct.

A Clamped Beam (Linear Case).
Figure 13 shows a clamped beam resting on a variable elastic foundation.
Table 9 indicates the variation of the frequency parameter Ω  associated with the ith mode shape, with  = 1 to 3, for the linear and parabolic soil stiffness variation, by increasing both the number of masses from  = 10 to  = 100 and the soil stiffness from 1 to 1000.It is obvious form the results obtained that the soil distribution (linear and parabolic) has no effect on the convergence.However, a low number of masses ( = 10) is not sufficient to get a good approximation (error about Advances in Acoustics and Vibration    10%).For  = 100, the results become accurate with an error around one percent.
In conclusion, by applying the present model to the problem of vibrations of beams resting on variable elastic foundations, the simply supported beam converges faster than the clamped beam.This convergence may be expected because of the hypothesis taken to solve the clamped beam problem, based on the reduced matrix, and form the behavior of the clamped beam, which is more restrained than the simply supported one.

A Clamped Beam (Nonlinear Case).
The backbone curve presenting the nondimensional frequency is drawn for the case of a clamped beam in Figure 14.In Figure 14(a), with a low soil stiffness ( = 10), the nonlinear behavior is more expressed since the nonlinear frequency is higher than the one of the high soil stiffness Figure 14(b) ( = 1000).On the other hand, the nonlinear frequency seems to be lower than the simply supported case.This frequency difference is less distinguishable for high soil stiffness  = 1000.

Conclusion
A discrete model of a beam made of multiple bars, masses, and spiral springs is introduced in [9].Accordingly, the present paper extends the previous model to the case of beams resting on elastic foundations, and validations are carried out for both the linear and nonlinear cases with various soil stiffness's and beam's end conditions.
Before validation of the results obtained via the extended theory and to complete the results of the previous paper [9], a mode shape analysis using the present model is investigated in the linear and nonlinear cases for simply supported and clamped beams.For the clamped beam, a good approximation is established, since amplitude dependent mode shapes are obtained, in agreement with previous theoretical and experimental works.However, in the simply supported case, the mode shape seemed to be slightly deviating from the original shape for high amplitudes, which is not conformed to the result mentioned in [9,18], according to which it is amplitude independent, namely, because of the hypothesis taken in the discrete model for the nonlinearity effect, which is sensitive to the amplitude variation (small triangle angle approximation).On the other hand, the nonlinear/linear frequency ratios are in a good agreement with the result given in [13], based on the FEM method.For the case of a simply supported beam, the linear theory is validated with an error less than three percent for  = 10 to an error less than 0.1% for  = 100 in all the cases considered.In the nonlinear case, the error increased for high amplitudes to seven percent but decreased to one percent by increasing the soil stiffness.
For the clamped beam, the convergence is relatively slower, which is expected because of the hypothesis taken.In the linear case, the error is about ten percent for  = 10 and dropped to one percent for  = 100.In the nonlinear case, the error increased for high amplitudes to six percent but it decreased to about three percent by increasing the soil stiffness.
The main objective of this paper is to construct a general discrete physical model enabling easy solutions in various configurations of a practical interest.Two applications were conducted.The first application deals with beams partially supported on elastic foundations, and the results obtained complies graphically with [26,27] in the linear case.For the nonlinear case, similar curves corresponding to the first nonlinear mode are plotted for various amplitudes ( * = 1 and  * = 3), to show the shifting from the linear to the nonlinear curves.The second application is concerned with beams resting on elastic foundations with a variable stiffness (linear and parabolic).The comparison of the results obtained with previous results [17] shows the convergence of the present method for different modes, soil stiffness, and beam ends conditions.For this application, the simply supported case is proven faster to converge on the linear case.For the nonlinear case, the backbone curves are drawn with various soil stiffness coefficients ( = 10,  = 1000), showing the deviation introduced by the nonlinearity effect.

B. Parameters Calculation Summary
The following appendix summarized the demonstrations already carried out in [9].B.1.Expressions for the General Terms of the Tensors m ij .The kinetic energy of the N-DOF system is the sum of the kinetic energies of the masses, considered as particles, concentrated at the nodes: The value of the concentrated mass   for a beam divided into ( + 1) bars is equal to   = /( + 1),  = 1, . . ., .
Identifying the terms of ( 6) and (B.1) leads to the following expression for the mass tensor in DB: usually encountered in the previous cases examined by the present method, for example, [20,23], are adopted here as follows: The rest of   are equal to zero.

Nomenclature * :
Nondimensional parameter symbol   : Contribution coefficient of the jth linear mode shape in MB for the discrete system   : The discrete modulus of the displacement   expressed in DB {A}: Displacement amplitudes of the masses  1 ; . . .;   ; . . .;   in DB   : The general term of the nonlinear rigidity tensor in DB   : The general term of the nonlinear rigidity tensor in MB : The ratio of the supported span to the total span [B]: The nonlinear rigidity matrix in DB [B]: The nonlinear rigidity matrix in MB   : Thestiffnesscoefficientoftheth spiral spring DB: Displacement basis Δ  : The longitudinal displacement following the vector   : Kronecker'ssymbol {  }: Theth displacement vector in DB : TheY oungmodulusofthebar' smaterial   : TheY oungmodulusoftheth bar : Thethicknessofthebarinm : The quadratic moment relative to the neutral fibre for the cross section in m 4 [I]: The identity matrix : The soil stiffness spring stiffness    : The general term of the linear rigidity tensor in DB (spiral spring)

Figure 1 :
Figure 1: A continuous beam supported on Winkler elastic foundation.

Figure 2 :
Figure 2: A multi-degree-of-freedom discrete model made of  concentrated masses, connected by longitudinal and spiral springs (bar characteristics) and longitudinal springs (soil stiffness).

Figure 4 :
Figure 4: The model considered for the discrete C-C beam by blocking the rotation of the first and the last bars.

Table 4 :Table 5 :
Comparison of the first three modes' nondimension frequency and the nondimensional Winkler soil stiffness for a C-C beam.C-C beam (linear) Comparison of the first mode shape at different positions of the beam in the C-C beam case.C-C beam (nonlinear)

Figure 13 :Figure 14 :
Figure 13: Model of C-C beam supported by a soil with a variable stiffness: (a) linear distribution; (b) parabolic distribution.

Table 1 :
Comparison of the nondimensional frequencies corresponding to the first 3 modes and the nondimensional Winkler soil stiffness for a S-S beam.
3.1.1.Linear Case.One has now to make the form of the matrices introduced in Section 2 explicit, knowing that the calculations are developed in the MB.The nondimensional mass matrix is reduced to the unity matrix [I], and the stiffness matrix [K SS  ]

Table 3 :
Comparison of the first three modes' nondimensional nonlinear frequency  ss * nl discr / ss *  discr , in S-S beam, for various values of the nondimensional Winkler soil stiffness.

Table 6 :
Comparison of first modes nondimensional nonlinear frequency  cc * nl discr / cc *  discr , for various values of the nondimensional Winkler soil stiffness for a C-C beam.

Table 7 :
The first eight modes' frequencies for the linear and parabolic soil distribution in S-S beam case.

Table 8 :
The first eight modes' frequencies for the linear and parabolic soil distribution in S-S beam case.

Table 9 :
The 3 first modes' frequencies for the linear and parabolic soil distribution in C-C beam case.