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

Authors: A. Khnaijar; R. Benamar
Journal: Advances in Acoustics and Vibration
 (2017)

---

## Abstract

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.


---

## Body

## 1. 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–3]. The Differential Quadrature Method (DQM) is also employed for the solution to similar engineering problems involving beam vibrations and foundations [4–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 two-parameter 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 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 parameterN (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: Section2.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 Section3, 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 Section4, 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.


## 2. 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.


### 2.1. 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 matricesmij and kij. On the other hand, the nonlinear behavior, presented by nonlinearity tensor bijkl 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.


#### 2.1.1. 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 (S,L,H,I,E,ρ, and k) defined in Nomenclature. The total beam strain energy can be written as the sum of the strain energy due to bending denoted as Vls and the strain energy induced by Winkler springs Vll plus the axial strain energy Vnl due to the axial load induced by large deflections. Thus Vls, Vll, Vnl, and the kinetic energy T are as follows [20]:(1)Vls=EI2∫0L∂2W∂x22dx,Vll=EI2∫0LW2dx,Vnl=ES8L∫0L∂W∂x2dx2,T=ρS2∫0L∂W∂t2dx,where W 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](2)Wx,t=qitwixfori,j=1,…,N.Figure 1
A continuous beam supported on Winkler elastic foundation.Figure 2
A multi-degree-of-freedom discrete model made ofN concentrated masses, connected by longitudinal and spiral springs (bar characteristics) and longitudinal springs (soil stiffness).By replacingW (2) in its new form into the energy expressions Vls, Vll, Vnl, and T (1), one gets(3)Vls=ksij¯2qiqj,(4)Vll=klij¯2qiqj,(5)Vnl=bijkl¯2qiqjqkql,(6)T=mij¯2q˙iq˙jin which mij¯, ksij¯, klij¯, and bijkl¯ are defined as follows in the modal basis (MB):(7)mij¯=ρS∫0Lwiwjdxfori,j=1,…,N,ksij¯=EI∫0L∂2wi∂x2∂2wj∂x2dxfori,j=1,…,N,klij¯=k∫0Lwiwjdxfori,j=1,…,N,bijkl¯=ES4L∫0L∂wi∂x∂wj∂xdx∫0L∂wk∂x∂wl∂xdxfori,j,k,l=1,…,N.For a conservative system, the dynamic behavior of the structure may be obtained by Lagrange’s equations [20]:(8)-∂∂t∂T∂q˙r+∂T∂qr-∂V∂qr=0forr=1,…,N,where V=Vls+Vll+Vnl. By introducing the energy expressions, that is, (3) to (6), into (8), one gets(9)mirq¨i+ksirqi+klirqi+2qiqjqkbijkr=0forr=1,…,N.This nonlinear differential system (9) may be expressed in a matrix form as(10)Mq¨+Ksq+Klq+2Bqq=0,where [M], [Ks], [Kl], 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](11)qi=aicos⁡ωnlt,where ai’s are the contribution coefficients to the beam nonlinear mode shapes. By replacing (11) in (10) and applying the harmonic balance method, one has(12)Ks+KlA-ω2MA+32BAA=0.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.


#### 2.1.2. 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 N-degree-of-freedom (DOF) system shown in Figure 2, made of N concentrated masses m1,…,mN, N+1 bars presented by longitudinal linear springs, N+2 torsional springs, and N longitudinal linear soil springs. Cr is the stiffness coefficient of therth spiral spring, for r=1 to N+2. The moment M in the spiral spring is given by M=-CrΔθ; Δθ=θr-θr-1 being the angle between the bars adjacent to the noder, considered as longitudinal springs of length lr and stiffness kbr, r=1, to N. The Winkler foundations are modeled using the longitudinal vertical spring’s distribution, with klr presenting the stiffness coefficient of therth linear spring, for r=1 to N.Using vector notation,d1;d2;…;di;…;dN is the displacement basis, denoted as DB, and y1 to yN are, respectively, the transverse displacements of the masses m1 to mN from the horizontal equilibrium positions (y=0) of the system corresponding to the nondeformed positions of the springs (spiral and linear). The vectors di are defined by(13)d1T=1000⋯0;d2T=0100⋯0;⋮diT=00⋯1⋯0;⋮dNT=0000⋯1,with diT being the unit displacement of theith mass. The displacement vector can be written in DB as(14)y=y1d1+y2d2+⋯+yidi+⋯+yNdN;by considering the modal basis Φ1,Φ2,…,Φi,…,ΦN, obtained by solution of the linear eigenvalue problem and denoted as MB, the displacement vector can be expressed as(15)y=y1¯Φ1+y2¯Φ2+⋯+yi¯Φi+⋯+yN¯ΦN,where Φi is the ith linear mode shape of the system and y1¯,y2¯,…,yi¯,…,yN¯ are the components of the mass displacements expressed in MB. Φi is denoted in DB as(16)ΦiT=φ1iφ2i⋯φ1Nfori=1,…,N.Using the transition matrixΦ from DB to MB, the two sets of coordinates in the two bases are related by(17)y1¯y2¯⋯yN¯=Φy1y2⋯yNT.So, the displacement vector can be written by analogy with the continuous displacement vectorqi, as [9](18)yi=Aidiscrcos⁡ωnldiscrt=ajφijcos⁡ωnldiscrtfori,j=1,…,N,where Aidiscr is the modulus of the displacement yi expressed in DB (or the contribution of the normalized vector di of DB) and aj is the modulus of the displacement yi expressed in MB (or the contribution of the normalized vector Φi of MB). ωnldiscr is the nonlinear frequency of the discrete system associated with the amplitude Aidiscr. Using the usual summation convention for the repeated indices i, j, k, and l, the kinetic, linear, and nonlinear strain energy expressions, obtained by replacement of the expression (18) of yi in the energy expressions (3) to (6), are (19)T=12aiajωnldiscr2mij¯sin2⁡ωnldiscrtfori,j=1,…,N,(20)Vls=12aiajksij¯cos2⁡ωnldiscrtfori,j=1,…,N,(21)Vll=12aiajklij¯cos2⁡ωnldiscrtfori,j=1,…,N,(22)Vnl=12aiajakalbijkl¯cos4⁡ωnldiscrtfori,j,k,l=1,…,N,where mij¯, ksij¯, klij¯, and bijkl¯ 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]:(23)mij¯=φsiφtjmst,ksij¯=φsiφtjksst,klij¯=φsiφtjklst,bijkl¯=φsiφtjφpkφqlbstpq.Replacing in (8) the expressions for T, Vll, Vls, and Vnl given in (19)–(22) and applying the harmonic balance method, in MB, lead to [24, 25](24)K¯s+K¯lA-ωnldiscr2M¯A+32B¯AA=0.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.


### 2.2. The Expressions for the General Terms of the Tensorsmij, ksij, klij, and bijkl

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 numberN. For the linear parameters mij and kij, 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 bijkl (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 N 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.AppendixB gives a summary of the calculation, carried out in [9], for the expression of the parameters mij, ksij, and bijkl. Identically, Sections 2.2.1(1) and 2.2.1(2) give detailed calculations of the new parameter klij.


#### 2.2.1. The Expressions for the General Terms of the Tensorsmij, ksij, klij, and bijkl

The general expression ofmij, ksij, klij, and bijkl for the discrete model are as follows:(i) The Mass tensorM/mij terms are expressed as (Appendix B)(25)mij=ρSLN+1δijfori,j=1,…,N.(ii) The spiral springs tensorKs/ksij terms are expressed as (Appendix B)(26)ksr-2r=N+13EIL3forr=3,…,N,ksr-1r=-4N+13EIL3forr=2,…,N,ksrr=6N+13EIL3forr=1,…,Nand the other values of ksij are obtained by symmetry relations or are equal to zero. The boundaries conditions of the system determine the stiffness kiis for i=1 and N+1. For simply supported ends, the corresponding torsional springs C1 and CN+1 stiffness are equal to zero, since the rotation is free; then kii=5N+13EI/L3 for i=1 and N+1. For clamped ends, the corresponding torsional spring stiffness is infinite, since the rotation is prevented, so kiis=+∞ for i=1 and N+1.(iii) The Winkler springs tensorKl/klij terms are expressed as (details given in Sections 2.2.1(1) and 2.2.1(2))(27)klij=kilδij=kiLN+1δijfori,j=1,…,N.(iv) The nonlinear rigidity tensorB/bijkl terms are expressed as (Appendix B)(28)biiii≅2ES8l3=2N+13ES8L3fori=1,…,N,bii-1i-1i-1=bi-1ii-1i-1=bi-1i-1ii-1=bi-1i-1i-1i≅-N+13ES8L3fori=2,…,N,biii-1i-1=bi-1iii-1=bi-1i-1ii=bii-1i-1i=bii-1ii-1=bi-1ii-1i≅N+13ES8L3fori=2,…,N,biiii-1=bi-1iii=bii-1ii=biii-1i≅-N+13ES8L3fori=2,…,Nand other values of bijkl 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 kli, shown in the discrete model of Figure 3, can be written as [9](29)Vlli=12kliyi2withy0=yN+1=0.Figure 3
Discretization of the linear vertical springs (soil stiffness).On the other hand, the elementary potential energy corresponding to an elementary section of lengthdx in the continuous beam resting on the elastic foundations is given by(30)dVll=12kyx2dxwith k corresponding to the soil stiffness per length unit. By approximation, one can write(31)yxi=yi,ki=ki,dx=li+li+12.By substituting (31) in (30), one obtains(32)dVll=12kili+li+12yi2fori=1,…,N.Since the bars have the same length (30) can be reduced to(33)dVll=12kilyi2fori=1,…,N.The identification between (29) and (33) gives(34)kli=kil=kiLN+1fori=1,…,N.(2) Expressions for the General Terms of the Tensors
k
l
i
j. The relationship existing between the longitudinal elongations Δyi of theith linear spring and the transverse displacements yi is(35)Δyi=yi.The potential energy stored in theith linear spring, of stiffness kli, is given directly in terms of the stiffness kli and the spring length variation yi, by the very well-known relationship, that is, 1/2 kliyi2, which leads, after summation, to(36)Vll=12∑i=1Nkliyi2fori=1toN.The identification between (4) and (30) leads to the following expression for the rigidity tensor associated with the Winkler foundations in DB:(37)klij=δijkli=δijkiLN+1fori,j=1,…,N,where δij is Kronecker’s symbol. Equation (37) can be written in a matrix form as(38)Kl=kilI=kiLN+1Iwith I being the identity matrix.


#### 2.2.2. 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(39)mijm∗ij=mij¯m∗ij¯=ρSLN+1fori,j=1,…,N,(40)klijkl∗ij=klij¯kl∗ij¯=EIN+13L3fori,j=1,…,N,(41)ksijks∗ij=ksij¯ks∗ij¯=EIN+13L3fori,j=1,…,N,(42)bijklb∗ijkl=bijkl¯b∗ijkl¯=ESN+138L3fori,j,k,l=1,…,N,(43)ω∗ω=EIρSL4.The Winkler stiffness coefficient is given by(44)klij=klδij=αiEIN+13L3δijfori,j=1,…,Nα
i and λi being the nondimensional parameter:(45)αi=kiL4EIN+14=λiN+14withλi=kiL4EI.In the case of a constant soil distributionαi=α and λi=λ.The nondimensional amplitudeA∗ is expressed as (46)A∗=ARwithR=IS,where R is the radius of gyration.The expression for the nonlinear frequency [19] is(47)ωnl2=aiajksij¯+klij¯+32aiajakalbijkl¯mij¯aiajfori,j,l,k=1,…,N.As the nondimensional formulation is established, the general expressions form∗ij, k∗sij, k∗lij, and b∗ijkl for the discrete model become as follows:(i) The nondimensional mass tensorM∗/m∗ij is obtained by replacing (25) in (39):(48)m∗ij=δijfori,j=1,…,N.(ii) The nondimensional spiral springs tensorK∗s/k∗sij is obtained by replacing (26) in (40):(49)k∗sr-2r=1forr=3,…,N,k∗sr-1r=-4forr=2,…,N,k∗srr=6forr=1,…,Nand the other values of k∗sij are obtained by symmetry relations or are equal to zero.(iii) The nondimensional Winkler springs tensorK∗l/k∗lij is obtained by replacing (27) in (41):(50)k∗lij=αδijfori,j=1,…,N.(iv) The nondimensional nonlinear rigidity tensorB∗/b∗ijkl is obtained by replacing (28) in (42):(51)b∗iiii=2fori=1,…,N,b∗ii-1i-1i-1=b∗i-1ii-1i-1=b∗i-1i-1ii-1=b∗i-1i-1i-1i=-1fori=2,…,N,b∗iii-1i-1=b∗i-1iii-1=b∗i-1i-1ii=b∗ii-1i-1i=b∗ii-1ii-1=b∗i-1ii-1i=1fori=2,…,N,b∗iiii-1=b∗i-1iii=b∗ii-1ii=b∗iii-1i=-1fori=2,…,Nand other values of b∗ijkl are equal to zero.


### 2.3. 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 εiεja1bijkr are neglected when considering the first nonlinear mode, in (24), rewritten in summation form as(52)klir¯+ksir¯-ω2mir¯ai+32aiajakbijkr¯=0forr=1,…,N.The nonlinear expressionaiajakb∗ijkr¯ is split to terms proportional to a13, terms proportional to a12εi, and terms proportional to a12εiεj which are neglected, leading to(53)aiajakb∗ijkr¯=a13b∗111r¯+a12εib∗11ir¯forr=1,…,N.The 1D ratio of the nonlinear nondimensional frequency established from (47) is given by (54)ω∗nl2=N+14ks∗11¯+kl∗11¯m∗11¯1+316b∗1111ks∗11¯+kl∗11¯¯a1R2.Then (52) can be written in a matrix form as(55)K∗lRI¯+K∗sRI¯-ωnl∗2MRI∗¯A∗RI+316αI∗A∗RI=-316a1R3b∗111r¯with(56)αI∗=a1R2b∗11ir¯,where K∗sRI¯=kij∗s¯, K∗lRI¯=kij∗l¯, and MRI∗¯=mij∗¯ are, respectively, the reduced beam rigidity matrix, the Winkler foundation rigidity matrix, and the mass matrix associated with the first nonlinear mode, obtained by varyingi andj in the set (3,5,…,11). αI∗ is a 5×5 square matrix, depending on a1, whose general term αij∗ is equal to b11ir¯a1/R2, and -3/16a1/R3b111r¯ is a column vector representing the right side of the linear system (55) in which the reduced unknown vector is A∗RIT=ε3/R,…,εi/R,…,εN/R. The modal contributions ε3,…,εi,…,εN can be calculated simply by solving the linear system (55) of five equations and five unknowns [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.


## 2.1. 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 matricesmij and kij. On the other hand, the nonlinear behavior, presented by nonlinearity tensor bijkl 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.


### 2.1.1. 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 (S,L,H,I,E,ρ, and k) defined in Nomenclature. The total beam strain energy can be written as the sum of the strain energy due to bending denoted as Vls and the strain energy induced by Winkler springs Vll plus the axial strain energy Vnl due to the axial load induced by large deflections. Thus Vls, Vll, Vnl, and the kinetic energy T are as follows [20]:(1)Vls=EI2∫0L∂2W∂x22dx,Vll=EI2∫0LW2dx,Vnl=ES8L∫0L∂W∂x2dx2,T=ρS2∫0L∂W∂t2dx,where W 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](2)Wx,t=qitwixfori,j=1,…,N.Figure 1
A continuous beam supported on Winkler elastic foundation.Figure 2
A multi-degree-of-freedom discrete model made ofN concentrated masses, connected by longitudinal and spiral springs (bar characteristics) and longitudinal springs (soil stiffness).By replacingW (2) in its new form into the energy expressions Vls, Vll, Vnl, and T (1), one gets(3)Vls=ksij¯2qiqj,(4)Vll=klij¯2qiqj,(5)Vnl=bijkl¯2qiqjqkql,(6)T=mij¯2q˙iq˙jin which mij¯, ksij¯, klij¯, and bijkl¯ are defined as follows in the modal basis (MB):(7)mij¯=ρS∫0Lwiwjdxfori,j=1,…,N,ksij¯=EI∫0L∂2wi∂x2∂2wj∂x2dxfori,j=1,…,N,klij¯=k∫0Lwiwjdxfori,j=1,…,N,bijkl¯=ES4L∫0L∂wi∂x∂wj∂xdx∫0L∂wk∂x∂wl∂xdxfori,j,k,l=1,…,N.For a conservative system, the dynamic behavior of the structure may be obtained by Lagrange’s equations [20]:(8)-∂∂t∂T∂q˙r+∂T∂qr-∂V∂qr=0forr=1,…,N,where V=Vls+Vll+Vnl. By introducing the energy expressions, that is, (3) to (6), into (8), one gets(9)mirq¨i+ksirqi+klirqi+2qiqjqkbijkr=0forr=1,…,N.This nonlinear differential system (9) may be expressed in a matrix form as(10)Mq¨+Ksq+Klq+2Bqq=0,where [M], [Ks], [Kl], 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](11)qi=aicos⁡ωnlt,where ai’s are the contribution coefficients to the beam nonlinear mode shapes. By replacing (11) in (10) and applying the harmonic balance method, one has(12)Ks+KlA-ω2MA+32BAA=0.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.


### 2.1.2. 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 N-degree-of-freedom (DOF) system shown in Figure 2, made of N concentrated masses m1,…,mN, N+1 bars presented by longitudinal linear springs, N+2 torsional springs, and N longitudinal linear soil springs. Cr is the stiffness coefficient of therth spiral spring, for r=1 to N+2. The moment M in the spiral spring is given by M=-CrΔθ; Δθ=θr-θr-1 being the angle between the bars adjacent to the noder, considered as longitudinal springs of length lr and stiffness kbr, r=1, to N. The Winkler foundations are modeled using the longitudinal vertical spring’s distribution, with klr presenting the stiffness coefficient of therth linear spring, for r=1 to N.Using vector notation,d1;d2;…;di;…;dN is the displacement basis, denoted as DB, and y1 to yN are, respectively, the transverse displacements of the masses m1 to mN from the horizontal equilibrium positions (y=0) of the system corresponding to the nondeformed positions of the springs (spiral and linear). The vectors di are defined by(13)d1T=1000⋯0;d2T=0100⋯0;⋮diT=00⋯1⋯0;⋮dNT=0000⋯1,with diT being the unit displacement of theith mass. The displacement vector can be written in DB as(14)y=y1d1+y2d2+⋯+yidi+⋯+yNdN;by considering the modal basis Φ1,Φ2,…,Φi,…,ΦN, obtained by solution of the linear eigenvalue problem and denoted as MB, the displacement vector can be expressed as(15)y=y1¯Φ1+y2¯Φ2+⋯+yi¯Φi+⋯+yN¯ΦN,where Φi is the ith linear mode shape of the system and y1¯,y2¯,…,yi¯,…,yN¯ are the components of the mass displacements expressed in MB. Φi is denoted in DB as(16)ΦiT=φ1iφ2i⋯φ1Nfori=1,…,N.Using the transition matrixΦ from DB to MB, the two sets of coordinates in the two bases are related by(17)y1¯y2¯⋯yN¯=Φy1y2⋯yNT.So, the displacement vector can be written by analogy with the continuous displacement vectorqi, as [9](18)yi=Aidiscrcos⁡ωnldiscrt=ajφijcos⁡ωnldiscrtfori,j=1,…,N,where Aidiscr is the modulus of the displacement yi expressed in DB (or the contribution of the normalized vector di of DB) and aj is the modulus of the displacement yi expressed in MB (or the contribution of the normalized vector Φi of MB). ωnldiscr is the nonlinear frequency of the discrete system associated with the amplitude Aidiscr. Using the usual summation convention for the repeated indices i, j, k, and l, the kinetic, linear, and nonlinear strain energy expressions, obtained by replacement of the expression (18) of yi in the energy expressions (3) to (6), are (19)T=12aiajωnldiscr2mij¯sin2⁡ωnldiscrtfori,j=1,…,N,(20)Vls=12aiajksij¯cos2⁡ωnldiscrtfori,j=1,…,N,(21)Vll=12aiajklij¯cos2⁡ωnldiscrtfori,j=1,…,N,(22)Vnl=12aiajakalbijkl¯cos4⁡ωnldiscrtfori,j,k,l=1,…,N,where mij¯, ksij¯, klij¯, and bijkl¯ 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]:(23)mij¯=φsiφtjmst,ksij¯=φsiφtjksst,klij¯=φsiφtjklst,bijkl¯=φsiφtjφpkφqlbstpq.Replacing in (8) the expressions for T, Vll, Vls, and Vnl given in (19)–(22) and applying the harmonic balance method, in MB, lead to [24, 25](24)K¯s+K¯lA-ωnldiscr2M¯A+32B¯AA=0.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.


## 2.1.1. 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 (S,L,H,I,E,ρ, and k) defined in Nomenclature. The total beam strain energy can be written as the sum of the strain energy due to bending denoted as Vls and the strain energy induced by Winkler springs Vll plus the axial strain energy Vnl due to the axial load induced by large deflections. Thus Vls, Vll, Vnl, and the kinetic energy T are as follows [20]:(1)Vls=EI2∫0L∂2W∂x22dx,Vll=EI2∫0LW2dx,Vnl=ES8L∫0L∂W∂x2dx2,T=ρS2∫0L∂W∂t2dx,where W 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](2)Wx,t=qitwixfori,j=1,…,N.Figure 1
A continuous beam supported on Winkler elastic foundation.Figure 2
A multi-degree-of-freedom discrete model made ofN concentrated masses, connected by longitudinal and spiral springs (bar characteristics) and longitudinal springs (soil stiffness).By replacingW (2) in its new form into the energy expressions Vls, Vll, Vnl, and T (1), one gets(3)Vls=ksij¯2qiqj,(4)Vll=klij¯2qiqj,(5)Vnl=bijkl¯2qiqjqkql,(6)T=mij¯2q˙iq˙jin which mij¯, ksij¯, klij¯, and bijkl¯ are defined as follows in the modal basis (MB):(7)mij¯=ρS∫0Lwiwjdxfori,j=1,…,N,ksij¯=EI∫0L∂2wi∂x2∂2wj∂x2dxfori,j=1,…,N,klij¯=k∫0Lwiwjdxfori,j=1,…,N,bijkl¯=ES4L∫0L∂wi∂x∂wj∂xdx∫0L∂wk∂x∂wl∂xdxfori,j,k,l=1,…,N.For a conservative system, the dynamic behavior of the structure may be obtained by Lagrange’s equations [20]:(8)-∂∂t∂T∂q˙r+∂T∂qr-∂V∂qr=0forr=1,…,N,where V=Vls+Vll+Vnl. By introducing the energy expressions, that is, (3) to (6), into (8), one gets(9)mirq¨i+ksirqi+klirqi+2qiqjqkbijkr=0forr=1,…,N.This nonlinear differential system (9) may be expressed in a matrix form as(10)Mq¨+Ksq+Klq+2Bqq=0,where [M], [Ks], [Kl], 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](11)qi=aicos⁡ωnlt,where ai’s are the contribution coefficients to the beam nonlinear mode shapes. By replacing (11) in (10) and applying the harmonic balance method, one has(12)Ks+KlA-ω2MA+32BAA=0.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.


## 2.1.2. 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 N-degree-of-freedom (DOF) system shown in Figure 2, made of N concentrated masses m1,…,mN, N+1 bars presented by longitudinal linear springs, N+2 torsional springs, and N longitudinal linear soil springs. Cr is the stiffness coefficient of therth spiral spring, for r=1 to N+2. The moment M in the spiral spring is given by M=-CrΔθ; Δθ=θr-θr-1 being the angle between the bars adjacent to the noder, considered as longitudinal springs of length lr and stiffness kbr, r=1, to N. The Winkler foundations are modeled using the longitudinal vertical spring’s distribution, with klr presenting the stiffness coefficient of therth linear spring, for r=1 to N.Using vector notation,d1;d2;…;di;…;dN is the displacement basis, denoted as DB, and y1 to yN are, respectively, the transverse displacements of the masses m1 to mN from the horizontal equilibrium positions (y=0) of the system corresponding to the nondeformed positions of the springs (spiral and linear). The vectors di are defined by(13)d1T=1000⋯0;d2T=0100⋯0;⋮diT=00⋯1⋯0;⋮dNT=0000⋯1,with diT being the unit displacement of theith mass. The displacement vector can be written in DB as(14)y=y1d1+y2d2+⋯+yidi+⋯+yNdN;by considering the modal basis Φ1,Φ2,…,Φi,…,ΦN, obtained by solution of the linear eigenvalue problem and denoted as MB, the displacement vector can be expressed as(15)y=y1¯Φ1+y2¯Φ2+⋯+yi¯Φi+⋯+yN¯ΦN,where Φi is the ith linear mode shape of the system and y1¯,y2¯,…,yi¯,…,yN¯ are the components of the mass displacements expressed in MB. Φi is denoted in DB as(16)ΦiT=φ1iφ2i⋯φ1Nfori=1,…,N.Using the transition matrixΦ from DB to MB, the two sets of coordinates in the two bases are related by(17)y1¯y2¯⋯yN¯=Φy1y2⋯yNT.So, the displacement vector can be written by analogy with the continuous displacement vectorqi, as [9](18)yi=Aidiscrcos⁡ωnldiscrt=ajφijcos⁡ωnldiscrtfori,j=1,…,N,where Aidiscr is the modulus of the displacement yi expressed in DB (or the contribution of the normalized vector di of DB) and aj is the modulus of the displacement yi expressed in MB (or the contribution of the normalized vector Φi of MB). ωnldiscr is the nonlinear frequency of the discrete system associated with the amplitude Aidiscr. Using the usual summation convention for the repeated indices i, j, k, and l, the kinetic, linear, and nonlinear strain energy expressions, obtained by replacement of the expression (18) of yi in the energy expressions (3) to (6), are (19)T=12aiajωnldiscr2mij¯sin2⁡ωnldiscrtfori,j=1,…,N,(20)Vls=12aiajksij¯cos2⁡ωnldiscrtfori,j=1,…,N,(21)Vll=12aiajklij¯cos2⁡ωnldiscrtfori,j=1,…,N,(22)Vnl=12aiajakalbijkl¯cos4⁡ωnldiscrtfori,j,k,l=1,…,N,where mij¯, ksij¯, klij¯, and bijkl¯ 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]:(23)mij¯=φsiφtjmst,ksij¯=φsiφtjksst,klij¯=φsiφtjklst,bijkl¯=φsiφtjφpkφqlbstpq.Replacing in (8) the expressions for T, Vll, Vls, and Vnl given in (19)–(22) and applying the harmonic balance method, in MB, lead to [24, 25](24)K¯s+K¯lA-ωnldiscr2M¯A+32B¯AA=0.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.


## 2.2. The Expressions for the General Terms of the Tensorsmij, ksij, klij, and bijkl

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 numberN. For the linear parameters mij and kij, 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 bijkl (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 N 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.AppendixB gives a summary of the calculation, carried out in [9], for the expression of the parameters mij, ksij, and bijkl. Identically, Sections 2.2.1(1) and 2.2.1(2) give detailed calculations of the new parameter klij.


### 2.2.1. The Expressions for the General Terms of the Tensorsmij, ksij, klij, and bijkl

The general expression ofmij, ksij, klij, and bijkl for the discrete model are as follows:(i) The Mass tensorM/mij terms are expressed as (Appendix B)(25)mij=ρSLN+1δijfori,j=1,…,N.(ii) The spiral springs tensorKs/ksij terms are expressed as (Appendix B)(26)ksr-2r=N+13EIL3forr=3,…,N,ksr-1r=-4N+13EIL3forr=2,…,N,ksrr=6N+13EIL3forr=1,…,Nand the other values of ksij are obtained by symmetry relations or are equal to zero. The boundaries conditions of the system determine the stiffness kiis for i=1 and N+1. For simply supported ends, the corresponding torsional springs C1 and CN+1 stiffness are equal to zero, since the rotation is free; then kii=5N+13EI/L3 for i=1 and N+1. For clamped ends, the corresponding torsional spring stiffness is infinite, since the rotation is prevented, so kiis=+∞ for i=1 and N+1.(iii) The Winkler springs tensorKl/klij terms are expressed as (details given in Sections 2.2.1(1) and 2.2.1(2))(27)klij=kilδij=kiLN+1δijfori,j=1,…,N.(iv) The nonlinear rigidity tensorB/bijkl terms are expressed as (Appendix B)(28)biiii≅2ES8l3=2N+13ES8L3fori=1,…,N,bii-1i-1i-1=bi-1ii-1i-1=bi-1i-1ii-1=bi-1i-1i-1i≅-N+13ES8L3fori=2,…,N,biii-1i-1=bi-1iii-1=bi-1i-1ii=bii-1i-1i=bii-1ii-1=bi-1ii-1i≅N+13ES8L3fori=2,…,N,biiii-1=bi-1iii=bii-1ii=biii-1i≅-N+13ES8L3fori=2,…,Nand other values of bijkl 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 kli, shown in the discrete model of Figure 3, can be written as [9](29)Vlli=12kliyi2withy0=yN+1=0.Figure 3
Discretization of the linear vertical springs (soil stiffness).On the other hand, the elementary potential energy corresponding to an elementary section of lengthdx in the continuous beam resting on the elastic foundations is given by(30)dVll=12kyx2dxwith k corresponding to the soil stiffness per length unit. By approximation, one can write(31)yxi=yi,ki=ki,dx=li+li+12.By substituting (31) in (30), one obtains(32)dVll=12kili+li+12yi2fori=1,…,N.Since the bars have the same length (30) can be reduced to(33)dVll=12kilyi2fori=1,…,N.The identification between (29) and (33) gives(34)kli=kil=kiLN+1fori=1,…,N.(2) Expressions for the General Terms of the Tensors
k
l
i
j. The relationship existing between the longitudinal elongations Δyi of theith linear spring and the transverse displacements yi is(35)Δyi=yi.The potential energy stored in theith linear spring, of stiffness kli, is given directly in terms of the stiffness kli and the spring length variation yi, by the very well-known relationship, that is, 1/2 kliyi2, which leads, after summation, to(36)Vll=12∑i=1Nkliyi2fori=1toN.The identification between (4) and (30) leads to the following expression for the rigidity tensor associated with the Winkler foundations in DB:(37)klij=δijkli=δijkiLN+1fori,j=1,…,N,where δij is Kronecker’s symbol. Equation (37) can be written in a matrix form as(38)Kl=kilI=kiLN+1Iwith I being the identity matrix.


### 2.2.2. 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(39)mijm∗ij=mij¯m∗ij¯=ρSLN+1fori,j=1,…,N,(40)klijkl∗ij=klij¯kl∗ij¯=EIN+13L3fori,j=1,…,N,(41)ksijks∗ij=ksij¯ks∗ij¯=EIN+13L3fori,j=1,…,N,(42)bijklb∗ijkl=bijkl¯b∗ijkl¯=ESN+138L3fori,j,k,l=1,…,N,(43)ω∗ω=EIρSL4.The Winkler stiffness coefficient is given by(44)klij=klδij=αiEIN+13L3δijfori,j=1,…,Nα
i and λi being the nondimensional parameter:(45)αi=kiL4EIN+14=λiN+14withλi=kiL4EI.In the case of a constant soil distributionαi=α and λi=λ.The nondimensional amplitudeA∗ is expressed as (46)A∗=ARwithR=IS,where R is the radius of gyration.The expression for the nonlinear frequency [19] is(47)ωnl2=aiajksij¯+klij¯+32aiajakalbijkl¯mij¯aiajfori,j,l,k=1,…,N.As the nondimensional formulation is established, the general expressions form∗ij, k∗sij, k∗lij, and b∗ijkl for the discrete model become as follows:(i) The nondimensional mass tensorM∗/m∗ij is obtained by replacing (25) in (39):(48)m∗ij=δijfori,j=1,…,N.(ii) The nondimensional spiral springs tensorK∗s/k∗sij is obtained by replacing (26) in (40):(49)k∗sr-2r=1forr=3,…,N,k∗sr-1r=-4forr=2,…,N,k∗srr=6forr=1,…,Nand the other values of k∗sij are obtained by symmetry relations or are equal to zero.(iii) The nondimensional Winkler springs tensorK∗l/k∗lij is obtained by replacing (27) in (41):(50)k∗lij=αδijfori,j=1,…,N.(iv) The nondimensional nonlinear rigidity tensorB∗/b∗ijkl is obtained by replacing (28) in (42):(51)b∗iiii=2fori=1,…,N,b∗ii-1i-1i-1=b∗i-1ii-1i-1=b∗i-1i-1ii-1=b∗i-1i-1i-1i=-1fori=2,…,N,b∗iii-1i-1=b∗i-1iii-1=b∗i-1i-1ii=b∗ii-1i-1i=b∗ii-1ii-1=b∗i-1ii-1i=1fori=2,…,N,b∗iiii-1=b∗i-1iii=b∗ii-1ii=b∗iii-1i=-1fori=2,…,Nand other values of b∗ijkl are equal to zero.


## 2.2.1. The Expressions for the General Terms of the Tensorsmij, ksij, klij, and bijkl

The general expression ofmij, ksij, klij, and bijkl for the discrete model are as follows:(i) The Mass tensorM/mij terms are expressed as (Appendix B)(25)mij=ρSLN+1δijfori,j=1,…,N.(ii) The spiral springs tensorKs/ksij terms are expressed as (Appendix B)(26)ksr-2r=N+13EIL3forr=3,…,N,ksr-1r=-4N+13EIL3forr=2,…,N,ksrr=6N+13EIL3forr=1,…,Nand the other values of ksij are obtained by symmetry relations or are equal to zero. The boundaries conditions of the system determine the stiffness kiis for i=1 and N+1. For simply supported ends, the corresponding torsional springs C1 and CN+1 stiffness are equal to zero, since the rotation is free; then kii=5N+13EI/L3 for i=1 and N+1. For clamped ends, the corresponding torsional spring stiffness is infinite, since the rotation is prevented, so kiis=+∞ for i=1 and N+1.(iii) The Winkler springs tensorKl/klij terms are expressed as (details given in Sections 2.2.1(1) and 2.2.1(2))(27)klij=kilδij=kiLN+1δijfori,j=1,…,N.(iv) The nonlinear rigidity tensorB/bijkl terms are expressed as (Appendix B)(28)biiii≅2ES8l3=2N+13ES8L3fori=1,…,N,bii-1i-1i-1=bi-1ii-1i-1=bi-1i-1ii-1=bi-1i-1i-1i≅-N+13ES8L3fori=2,…,N,biii-1i-1=bi-1iii-1=bi-1i-1ii=bii-1i-1i=bii-1ii-1=bi-1ii-1i≅N+13ES8L3fori=2,…,N,biiii-1=bi-1iii=bii-1ii=biii-1i≅-N+13ES8L3fori=2,…,Nand other values of bijkl 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 kli, shown in the discrete model of Figure 3, can be written as [9](29)Vlli=12kliyi2withy0=yN+1=0.Figure 3
Discretization of the linear vertical springs (soil stiffness).On the other hand, the elementary potential energy corresponding to an elementary section of lengthdx in the continuous beam resting on the elastic foundations is given by(30)dVll=12kyx2dxwith k corresponding to the soil stiffness per length unit. By approximation, one can write(31)yxi=yi,ki=ki,dx=li+li+12.By substituting (31) in (30), one obtains(32)dVll=12kili+li+12yi2fori=1,…,N.Since the bars have the same length (30) can be reduced to(33)dVll=12kilyi2fori=1,…,N.The identification between (29) and (33) gives(34)kli=kil=kiLN+1fori=1,…,N.(2) Expressions for the General Terms of the Tensors
k
l
i
j. The relationship existing between the longitudinal elongations Δyi of theith linear spring and the transverse displacements yi is(35)Δyi=yi.The potential energy stored in theith linear spring, of stiffness kli, is given directly in terms of the stiffness kli and the spring length variation yi, by the very well-known relationship, that is, 1/2 kliyi2, which leads, after summation, to(36)Vll=12∑i=1Nkliyi2fori=1toN.The identification between (4) and (30) leads to the following expression for the rigidity tensor associated with the Winkler foundations in DB:(37)klij=δijkli=δijkiLN+1fori,j=1,…,N,where δij is Kronecker’s symbol. Equation (37) can be written in a matrix form as(38)Kl=kilI=kiLN+1Iwith I being the identity matrix.


## 2.2.2. 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(39)mijm∗ij=mij¯m∗ij¯=ρSLN+1fori,j=1,…,N,(40)klijkl∗ij=klij¯kl∗ij¯=EIN+13L3fori,j=1,…,N,(41)ksijks∗ij=ksij¯ks∗ij¯=EIN+13L3fori,j=1,…,N,(42)bijklb∗ijkl=bijkl¯b∗ijkl¯=ESN+138L3fori,j,k,l=1,…,N,(43)ω∗ω=EIρSL4.The Winkler stiffness coefficient is given by(44)klij=klδij=αiEIN+13L3δijfori,j=1,…,Nα
i and λi being the nondimensional parameter:(45)αi=kiL4EIN+14=λiN+14withλi=kiL4EI.In the case of a constant soil distributionαi=α and λi=λ.The nondimensional amplitudeA∗ is expressed as (46)A∗=ARwithR=IS,where R is the radius of gyration.The expression for the nonlinear frequency [19] is(47)ωnl2=aiajksij¯+klij¯+32aiajakalbijkl¯mij¯aiajfori,j,l,k=1,…,N.As the nondimensional formulation is established, the general expressions form∗ij, k∗sij, k∗lij, and b∗ijkl for the discrete model become as follows:(i) The nondimensional mass tensorM∗/m∗ij is obtained by replacing (25) in (39):(48)m∗ij=δijfori,j=1,…,N.(ii) The nondimensional spiral springs tensorK∗s/k∗sij is obtained by replacing (26) in (40):(49)k∗sr-2r=1forr=3,…,N,k∗sr-1r=-4forr=2,…,N,k∗srr=6forr=1,…,Nand the other values of k∗sij are obtained by symmetry relations or are equal to zero.(iii) The nondimensional Winkler springs tensorK∗l/k∗lij is obtained by replacing (27) in (41):(50)k∗lij=αδijfori,j=1,…,N.(iv) The nondimensional nonlinear rigidity tensorB∗/b∗ijkl is obtained by replacing (28) in (42):(51)b∗iiii=2fori=1,…,N,b∗ii-1i-1i-1=b∗i-1ii-1i-1=b∗i-1i-1ii-1=b∗i-1i-1i-1i=-1fori=2,…,N,b∗iii-1i-1=b∗i-1iii-1=b∗i-1i-1ii=b∗ii-1i-1i=b∗ii-1ii-1=b∗i-1ii-1i=1fori=2,…,N,b∗iiii-1=b∗i-1iii=b∗ii-1ii=b∗iii-1i=-1fori=2,…,Nand other values of b∗ijkl are equal to zero.


## 2.3. 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 εiεja1bijkr are neglected when considering the first nonlinear mode, in (24), rewritten in summation form as(52)klir¯+ksir¯-ω2mir¯ai+32aiajakbijkr¯=0forr=1,…,N.The nonlinear expressionaiajakb∗ijkr¯ is split to terms proportional to a13, terms proportional to a12εi, and terms proportional to a12εiεj which are neglected, leading to(53)aiajakb∗ijkr¯=a13b∗111r¯+a12εib∗11ir¯forr=1,…,N.The 1D ratio of the nonlinear nondimensional frequency established from (47) is given by (54)ω∗nl2=N+14ks∗11¯+kl∗11¯m∗11¯1+316b∗1111ks∗11¯+kl∗11¯¯a1R2.Then (52) can be written in a matrix form as(55)K∗lRI¯+K∗sRI¯-ωnl∗2MRI∗¯A∗RI+316αI∗A∗RI=-316a1R3b∗111r¯with(56)αI∗=a1R2b∗11ir¯,where K∗sRI¯=kij∗s¯, K∗lRI¯=kij∗l¯, and MRI∗¯=mij∗¯ are, respectively, the reduced beam rigidity matrix, the Winkler foundation rigidity matrix, and the mass matrix associated with the first nonlinear mode, obtained by varyingi andj in the set (3,5,…,11). αI∗ is a 5×5 square matrix, depending on a1, whose general term αij∗ is equal to b11ir¯a1/R2, and -3/16a1/R3b111r¯ is a column vector representing the right side of the linear system (55) in which the reduced unknown vector is A∗RIT=ε3/R,…,εi/R,…,εN/R. The modal contributions ε3,…,εi,…,εN can be calculated simply by solving the linear system (55) of five equations and five unknowns [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.


## 3. Results Validation

In Section2, 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.


### 3.1. Simply Supported Beams Resting on Elastic Foundations

In Figure2, 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.


#### 3.1.1. Linear Case

One has now to make the form of the matrices introduced in Section2 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 KNSS, 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]:(57)KNSS=N+13EIL35+α-410··0-46+α-410··1-4·····01····0·0···-410·01-46+α-40··01-45+α,KNSS=N+13EIL3KNSS∗.With this new expression for the nondimensional rigidity matrixKNSS∗ (24), corresponding to linear vibrations, can be written as(58)KNSS∗y+ωdiscr2ρSL4EIN+14y=0with y being the vector of mass displacements in the MB. Solving (58) gives the eigenvalues β, which allows the corresponding frequencies to be written as(59)ωlSSdiscr=N+12βiEIρSL4.Table1 gives, for the first 3 modes of the simply supported beam, the nondimensional frequencies, for increasing number of masses (N=10,50,100) and increasing values of the Winkler foundation stiffness. Accordingly, these nondimensional frequencies are calculated from (59), by extracting the EI/ρSL4 term, where the eigenvalues βi are those of the nondimensional stiffness matrix KNSS∗. So, they can be presented as(60)ωlSSdiscr∗=N+12βi.Table 1
Comparison of the nondimensional frequencies corresponding to the first 3 modes and the nondimensional Winkler soil stiffness for a S-S beam.
S-S beam (linear)
λ
Present study(N=10)
Present study(N=50)
Present study(N=100)
Ref. [10]
Ref. [11]
ExactRef. [12]
Average
Error %(N=10)
Error %(N=50)
Error %(N=100)
0
3,131
3,141
3,141
3,150
3,141
3,141
3,144
0,427
0,104
0,092
6,198
6,279
6,2822
6,280
6,280
1,304
0,013
−0,035
9,139
9,411
9,421
9,420
9,420
2,982
0,092
−0,014
100
3,742
3,748
3,7483
3,750
3,748
3,7484
3,749
0,181
0,022
0,016
6,300
6,378
6,3807
6,380
6,380
1,245
0,034
−0,010
9,172
9,441
9,4511
9,450
9,450
2,945
0,093
−0,012
10000
10,024
10,024
10,024
10,020
10,024
10,024
10,023
−0,013
−0,016
−0,016
10,350
10,368
10,368
10,360
10,360
0,095
−0,076
−0,082
11,415
11,558
11,563
11,570
11,570
1,343
0,104
0,057
1000000
31,623
31,623
31,623
31,620
31,623
31,623
31,622
−0,004
−0,004
−0,004
31,634
31,635
31,635
31,630
31,630
−0,016
−0,014
−0,014
31,678
31,685
31,685
31,720
31,720
0,133
0,112
0,111In order to validate the model introduced in Section2, for the simply supported beam, a comparison of the first 3 natural frequencies is carried out for increasing values of N (the number of masses) and λ (the soil stiffness), and results are drawn in Table 1.The natural frequencies, Table1, converge faster by increasing the N or λ. A good approximation is already obtained for (N=10, λ =0) where the overall error ((average − present study (N))/average %) is below three percent limit. On the other hand, by increasing the soil stiffness to the case (N=10, λ = 1000000) the error drops below 0.15%. However, once the discretization number N is higher than 50 the overall error is under the 0.12%. Thus, increasing one of the parameters λ and N, or both, makes the discrete model converge to the reference models [10–12].


#### 3.1.2. 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].Table 2
Comparison of the normalized first mode shape at different positions of the beam in the S-S beam case forN=100.
ξ
=
x
/
L
ExactRef. [13]
Present
FEMRef. [13]
Present
FEMRef. [13]
Present
FEMRef. [13]
Present
FEMRef. [13]
A
∗
=
0
A
∗
=
0
A
∗
=
0
A
∗
=
1
A
∗
=
1
A
∗
=
2
A
∗
=
2
A
∗
=
3
A
∗
=
3
0,00
0,0000
0,0000
0,0000
0,0000
0,0000
0,0000
0,0000
0,0000
0,0000
0,05
0,1564
0,1564
0,1564
0,1544
0,1564
0,1498
0,1564
0,1449
0,1564
0,10
0,3090
0,3090
0,3090
0,3054
0,3090
0,2969
0,3090
0,2876
0,3090
0,15
0,4540
0,4540
0,4540
0,4493
0,4540
0,4381
0,4540
0,4258
0,4540
0,20
0,5878
0,5878
0,5880
0,5827
0,5880
0,5707
0,5880
0,5570
0,5880
0,25
0,7071
0,7071
0,7071
0,7024
0,7071
0,6911
0,7071
0,6780
0,7071
0,30
0,8091
0,8090
0,8090
0,8053
0,8090
0,7962
0,8090
0,7853
0,8090
0,35
0,8910
0,8910
0,8910
0,8886
0,8910
0,8825
0,8910
0,8750
0,8910
0,40
0,9511
0,9511
0,9511
0,9498
0,9511
0,9468
0,9511
0,9430
0,9511
0,45
0,9877
0,9877
0,9877
0,9874
0,9877
0,9865
0,9877
0,9855
0,9877
0,50
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000The results are identical to the exact theoretical mode shape in the linear case. However, in the nonlinear case (high amplitudes), in spite of the result often accepted in the literature, according to which the mode shape in the simply supported case is amplitude independent [18], this discrete physical model leads to mode shapes which seem to be amplitude dependent. As a result, the amplitude dependency may be due to the hypothesis (small triangle angle) taken for the nonlinearity effect which is sensitive to high amplitudes.(2) Comparison of Nondimension Frequency
ω
n
l
s
s
∗
d
i
s
c
r
/
ω
l
s
s
∗
d
i
s
c
r
. For validating the results in the nonlinear case, comparing the ratio of the nonlinear frequency parameter (ωnlss∗discr/ωlss∗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 (N=100), because of the good approximation obtained for this degree of discretization in the linear case (error below 0.1%).Table3 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 (ωnlss∗discr/ωlss∗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.Table 3
Comparison of the first three modes’ nondimensional nonlinear frequencyωnlss∗discr/ωlss∗discr, in S-S beam, for various values of the nondimensional Winkler soil stiffness.
SIMPLY SUPPORTED BEAM
λ
0,1π4
1
π
4
10
π
4
A
∗
N
=
100
Ref. [14]
Error %
N
=
100
Ref. [14]
Error %
N
=
100
Ref. [14]
Error %
0,2
1,0026
1,0034
0,0843
1,0014
1,0019
0,0494
1,0003
1,0003
0,0044
0,4
1,0102
1,0135
0,3278
1,0056
1,0075
0,1876
1,0010
1,0014
0,0377
0,6
1,0228
1,0301
0,7128
1,0126
1,0167
0,4053
1,0023
1,0031
0,0799
0,8
1,0401
1,0529
1,2144
1,0223
1,0295
0,7036
1,0041
1,0054
0,1310
1
1,0620
1,0814
1,7928
1,0346
1,0457
1,0649
1,0064
1,0085
0,2110
2
1,2294
1,2924
4,8729
1,1319
1,1708
3,3188
1,0253
1,0334
0,7885
3
1,4666
1,5795
7,1497
1,2779
1,3519
5,4771
1,0560
1,0736
1,6424


### 3.2. 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 Section2.


#### 3.2.1. 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:(61)KNCC=N+13EIL3+∞-410··0-46+α-410··1-4·····01····0·0···-410·01-46+α-40··01-4+∞.However, once introduced in the calculation, a mathematical error about the conditioning of the matrix arises (ill-conditioned 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 aty=0 as shown in Figure 4. Using this method for taking into account the clamped ends conditions, theN-DOF clamped beam problem is reduced to a (N-2)-DOF problem by removing the first and last row and column of the stiffness matrix; then (61) becomes(62)KN-2CC=N+13EIL36+α-410··0-46+α-410··1-4·····01····0·0···-410·01-46+α-40··01-46+α,KN-2CC=N+13EIL3KN-2CC∗.Figure 4
The model considered for the discrete C-C beam by blocking the rotation of the first and the last bars.With this new expression for the nondimensional rigidity matrixKN-2CC∗, the linear vibration equation (24) can be described by(63)KN-2CC∗y+ωlCCdiscr2ρSL4EIN+14y=0with y being the vector of mass displacements in MB. Solving (63) gives the eigenvalues β so one can write the frequencies as(64)ωlCCdiscr=N+12βiEIρSL4.Table4 gives the first three modes of the clamped beam nondimensional frequencies for increasing number of masses (N=10,50,100) and increasing values of the Winkler foundation stiffness. These nondimensional frequencies are calculated from (64) by extracting the EI/ρSL4 term, where the eigenvalues βi are those of the nondimensional stiffness matrix KN-2CC∗, so they can be presented as(65)ωlCCdiscr∗=N+12βi.Table 4
Comparison of the first three modes’ nondimension frequency and the nondimensional Winkler soil stiffness for a C-C beam.
C-C beam (linear)
λ
Present studyN=10
Present studyN=50
Present studyN=100
Ref. [11]
ExactRef. [12]
Error %N=10
Error %N=50
Error %N=100
0
5,2153
4,8251
4,7775
4,7300
4,7300
−10,261
−2,011
−1,003
8,5727
8,0083
7,9313
7,8500
7,8540
−9,150
−1,965
−0,984
11,7658
11,2052
11,1030
11,0000
10,9960
−7,001
−1,902
−0,973
100
5,3833
5,0338
4,9919
4,9500
4,9500
−8,753
−1,692
−0,846
8,6121
8,0565
7,9809
7,9000
7,9040
−8,959
−1,930
−0,973
11,7811
11,2229
11,1212
11,0100
11,0140
−6,965
−1,897
−0,974
10000
10,1800
10,1328
10,1278
10,1200
10,1230
−0,563
−0,097
−0,047
11,1400
10,8995
10,8692
10,8400
10,8390
−2,777
−0,558
−0,279
13,0681
12,6693
12,5991
12,5300
12,5260
−4,328
−1,144
−0,583
1000000
31,6286
31,6271
31,6269
31,6400
31,6260
−0,008
−0,003
−0,003
31,6654
31,6552
31,6540
31,6700
31,6530
−0,039
−0,007
−0,003
31,7732
31,7467
31,7422
31,7500
31,7380
−0,111
−0,027
−0,013In order to validate the discrete model, for clamped beams, Table4 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 massesN= (10/50/100), and the theoretical values of the natural frequencies of continuous beams obtained in [12].In conclusion, the natural frequencies drawn in Table4 converge slower to the continuous model [12] than the simply supported case. The error percentage ((results [12] − present study (N))/results [12]%), for the case of no soil stiffness applied λ = 0, do not get below the one percent until N=100. However, once the soil stiffness λ is increased, a good approximation is already obtained for the case of N=10. This convergence rate is due to the approximations taken by contraction of the matrices in the clamped beam case.


#### 3.2.2. 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.Table 5
Comparison of the first mode shape at different positions of the beam in the C-C beam case.
C-C beam (nonlinear)
ξ
=
x
/
L
PresentN=100
ASMRef. [15]
FEMRef. [13]
PresentN=100
GFEMRef. [16]
FEMRef. [13]
PresentN=100
GFEMRef. [16]
FEMRef. [13]
A
∗
=
0
A
∗
=
2
A
∗
=
5
0,00
0,0000
0,0000
0,0000
0,0000
0,0000
0,0000
0,0000
0,0000
0,0000
0,05
0,0324
0,0325
0,0329
0,0341
0,0339
0,0341
0,0464
0,0411
0,0413
0,10
0,1196
0,1119
0,1203
0,1242
0,1234
0,1237
0,1555
0,1429
0,1432
0,15
0,2444
0,2435
0,2454
0,2508
0,2505
0,2506
0,2927
0,2795
0,2796
0,20
0,3906
0,3900
0,3922
0,3969
0,3981
0,3982
0,4363
0,4307
0,4308
0,25
0,5435
0,5435
0,5455
0,5480
0,5511
0,5513
0,5757
0,5815
0,5817
0,30
0,6895
0,6901
0,6918
0,6917
0,6960
0,6964
0,7058
0,7200
0,7204
0,35
0,8168
0,8178
0,8189
0,8174
0,8214
0,8219
0,8217
0,8372
0,8376
0,40
0,9156
0,9164
0,9170
0,9155
0,9182
0,9185
0,9159
0,9259
0,9262
0,45
0,9783
0,9787
0,9789
0,9782
0,9792
0,9293
0,9783
0,9812
0,9813
0,50
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000For clamped beams, the mode shapes are shown both theoretically and experimentally to be amplitude dependent at large vibration amplitudes. Table5 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
ω
n
l
c
c
∗
d
i
s
c
r
/
ω
l
c
c
∗
d
i
s
c
r. 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 (ωnlcc∗discr/ωlcc∗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 (N=100), for which a good convergence is reached in the linear case (error below 1%).Table 6
Comparison of first modes nondimensional nonlinear frequencyωnlcc∗discr/ωlcc∗discr, for various values of the nondimensional Winkler soil stiffness for a C-C beam.
Clamped-clamped beam
λ
0,1π4
1
π
4
10
π
4
A
∗
N
=
100
Ref. [14]
Error %
N
=
100
Ref. [14]
Error %
N
=
100
Ref. [14]
Error %
0,2
1,0013
1,0009
−0,0397
1,0011
1,0008
−0,0314
1,0005
1,0003
−0,0161
0,4
1,0052
1,0037
−0,1475
1,0044
1,0032
−0,1244
1,0018
1,0013
−0,0541
0,6
1,0116
1,0082
−0,3391
1,0100
1,0071
−0,2859
1,0041
1,0029
−0,1236
0,8
1,0206
1,0146
−0,5879
1,0177
1,0125
−0,5110
1,0073
1,0052
−0,2137
1
1,0320
1,0227
−0,9046
1,0275
1,0195
−0,7829
1,0115
1,0081
−0,3330
2
1,1224
1,0876
−3,1963
1,1058
1,0756
−2,8113
1,0451
1,0320
−1,2670
3
1,2587
1,1870
−6,0403
1,2253
1,1625
−5,4062
1,0988
1,0705
−2,6458The nonlinear behavior of a clamped beam resting on elastic foundations is presented in Table6 by means of the amplitude dependence of the parameter (ωnlcc∗discr/ωlcc∗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 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 λ.


## 3.1. Simply Supported Beams Resting on Elastic Foundations

In Figure2, 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.


### 3.1.1. Linear Case

One has now to make the form of the matrices introduced in Section2 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 KNSS, 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]:(57)KNSS=N+13EIL35+α-410··0-46+α-410··1-4·····01····0·0···-410·01-46+α-40··01-45+α,KNSS=N+13EIL3KNSS∗.With this new expression for the nondimensional rigidity matrixKNSS∗ (24), corresponding to linear vibrations, can be written as(58)KNSS∗y+ωdiscr2ρSL4EIN+14y=0with y being the vector of mass displacements in the MB. Solving (58) gives the eigenvalues β, which allows the corresponding frequencies to be written as(59)ωlSSdiscr=N+12βiEIρSL4.Table1 gives, for the first 3 modes of the simply supported beam, the nondimensional frequencies, for increasing number of masses (N=10,50,100) and increasing values of the Winkler foundation stiffness. Accordingly, these nondimensional frequencies are calculated from (59), by extracting the EI/ρSL4 term, where the eigenvalues βi are those of the nondimensional stiffness matrix KNSS∗. So, they can be presented as(60)ωlSSdiscr∗=N+12βi.Table 1
Comparison of the nondimensional frequencies corresponding to the first 3 modes and the nondimensional Winkler soil stiffness for a S-S beam.
S-S beam (linear)
λ
Present study(N=10)
Present study(N=50)
Present study(N=100)
Ref. [10]
Ref. [11]
ExactRef. [12]
Average
Error %(N=10)
Error %(N=50)
Error %(N=100)
0
3,131
3,141
3,141
3,150
3,141
3,141
3,144
0,427
0,104
0,092
6,198
6,279
6,2822
6,280
6,280
1,304
0,013
−0,035
9,139
9,411
9,421
9,420
9,420
2,982
0,092
−0,014
100
3,742
3,748
3,7483
3,750
3,748
3,7484
3,749
0,181
0,022
0,016
6,300
6,378
6,3807
6,380
6,380
1,245
0,034
−0,010
9,172
9,441
9,4511
9,450
9,450
2,945
0,093
−0,012
10000
10,024
10,024
10,024
10,020
10,024
10,024
10,023
−0,013
−0,016
−0,016
10,350
10,368
10,368
10,360
10,360
0,095
−0,076
−0,082
11,415
11,558
11,563
11,570
11,570
1,343
0,104
0,057
1000000
31,623
31,623
31,623
31,620
31,623
31,623
31,622
−0,004
−0,004
−0,004
31,634
31,635
31,635
31,630
31,630
−0,016
−0,014
−0,014
31,678
31,685
31,685
31,720
31,720
0,133
0,112
0,111In order to validate the model introduced in Section2, for the simply supported beam, a comparison of the first 3 natural frequencies is carried out for increasing values of N (the number of masses) and λ (the soil stiffness), and results are drawn in Table 1.The natural frequencies, Table1, converge faster by increasing the N or λ. A good approximation is already obtained for (N=10, λ =0) where the overall error ((average − present study (N))/average %) is below three percent limit. On the other hand, by increasing the soil stiffness to the case (N=10, λ = 1000000) the error drops below 0.15%. However, once the discretization number N is higher than 50 the overall error is under the 0.12%. Thus, increasing one of the parameters λ and N, or both, makes the discrete model converge to the reference models [10–12].


### 3.1.2. 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].Table 2
Comparison of the normalized first mode shape at different positions of the beam in the S-S beam case forN=100.
ξ
=
x
/
L
ExactRef. [13]
Present
FEMRef. [13]
Present
FEMRef. [13]
Present
FEMRef. [13]
Present
FEMRef. [13]
A
∗
=
0
A
∗
=
0
A
∗
=
0
A
∗
=
1
A
∗
=
1
A
∗
=
2
A
∗
=
2
A
∗
=
3
A
∗
=
3
0,00
0,0000
0,0000
0,0000
0,0000
0,0000
0,0000
0,0000
0,0000
0,0000
0,05
0,1564
0,1564
0,1564
0,1544
0,1564
0,1498
0,1564
0,1449
0,1564
0,10
0,3090
0,3090
0,3090
0,3054
0,3090
0,2969
0,3090
0,2876
0,3090
0,15
0,4540
0,4540
0,4540
0,4493
0,4540
0,4381
0,4540
0,4258
0,4540
0,20
0,5878
0,5878
0,5880
0,5827
0,5880
0,5707
0,5880
0,5570
0,5880
0,25
0,7071
0,7071
0,7071
0,7024
0,7071
0,6911
0,7071
0,6780
0,7071
0,30
0,8091
0,8090
0,8090
0,8053
0,8090
0,7962
0,8090
0,7853
0,8090
0,35
0,8910
0,8910
0,8910
0,8886
0,8910
0,8825
0,8910
0,8750
0,8910
0,40
0,9511
0,9511
0,9511
0,9498
0,9511
0,9468
0,9511
0,9430
0,9511
0,45
0,9877
0,9877
0,9877
0,9874
0,9877
0,9865
0,9877
0,9855
0,9877
0,50
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000The results are identical to the exact theoretical mode shape in the linear case. However, in the nonlinear case (high amplitudes), in spite of the result often accepted in the literature, according to which the mode shape in the simply supported case is amplitude independent [18], this discrete physical model leads to mode shapes which seem to be amplitude dependent. As a result, the amplitude dependency may be due to the hypothesis (small triangle angle) taken for the nonlinearity effect which is sensitive to high amplitudes.(2) Comparison of Nondimension Frequency
ω
n
l
s
s
∗
d
i
s
c
r
/
ω
l
s
s
∗
d
i
s
c
r
. For validating the results in the nonlinear case, comparing the ratio of the nonlinear frequency parameter (ωnlss∗discr/ωlss∗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 (N=100), because of the good approximation obtained for this degree of discretization in the linear case (error below 0.1%).Table3 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 (ωnlss∗discr/ωlss∗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.Table 3
Comparison of the first three modes’ nondimensional nonlinear frequencyωnlss∗discr/ωlss∗discr, in S-S beam, for various values of the nondimensional Winkler soil stiffness.
SIMPLY SUPPORTED BEAM
λ
0,1π4
1
π
4
10
π
4
A
∗
N
=
100
Ref. [14]
Error %
N
=
100
Ref. [14]
Error %
N
=
100
Ref. [14]
Error %
0,2
1,0026
1,0034
0,0843
1,0014
1,0019
0,0494
1,0003
1,0003
0,0044
0,4
1,0102
1,0135
0,3278
1,0056
1,0075
0,1876
1,0010
1,0014
0,0377
0,6
1,0228
1,0301
0,7128
1,0126
1,0167
0,4053
1,0023
1,0031
0,0799
0,8
1,0401
1,0529
1,2144
1,0223
1,0295
0,7036
1,0041
1,0054
0,1310
1
1,0620
1,0814
1,7928
1,0346
1,0457
1,0649
1,0064
1,0085
0,2110
2
1,2294
1,2924
4,8729
1,1319
1,1708
3,3188
1,0253
1,0334
0,7885
3
1,4666
1,5795
7,1497
1,2779
1,3519
5,4771
1,0560
1,0736
1,6424


## 3.1.1. Linear Case

One has now to make the form of the matrices introduced in Section2 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 KNSS, 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]:(57)KNSS=N+13EIL35+α-410··0-46+α-410··1-4·····01····0·0···-410·01-46+α-40··01-45+α,KNSS=N+13EIL3KNSS∗.With this new expression for the nondimensional rigidity matrixKNSS∗ (24), corresponding to linear vibrations, can be written as(58)KNSS∗y+ωdiscr2ρSL4EIN+14y=0with y being the vector of mass displacements in the MB. Solving (58) gives the eigenvalues β, which allows the corresponding frequencies to be written as(59)ωlSSdiscr=N+12βiEIρSL4.Table1 gives, for the first 3 modes of the simply supported beam, the nondimensional frequencies, for increasing number of masses (N=10,50,100) and increasing values of the Winkler foundation stiffness. Accordingly, these nondimensional frequencies are calculated from (59), by extracting the EI/ρSL4 term, where the eigenvalues βi are those of the nondimensional stiffness matrix KNSS∗. So, they can be presented as(60)ωlSSdiscr∗=N+12βi.Table 1
Comparison of the nondimensional frequencies corresponding to the first 3 modes and the nondimensional Winkler soil stiffness for a S-S beam.
S-S beam (linear)
λ
Present study(N=10)
Present study(N=50)
Present study(N=100)
Ref. [10]
Ref. [11]
ExactRef. [12]
Average
Error %(N=10)
Error %(N=50)
Error %(N=100)
0
3,131
3,141
3,141
3,150
3,141
3,141
3,144
0,427
0,104
0,092
6,198
6,279
6,2822
6,280
6,280
1,304
0,013
−0,035
9,139
9,411
9,421
9,420
9,420
2,982
0,092
−0,014
100
3,742
3,748
3,7483
3,750
3,748
3,7484
3,749
0,181
0,022
0,016
6,300
6,378
6,3807
6,380
6,380
1,245
0,034
−0,010
9,172
9,441
9,4511
9,450
9,450
2,945
0,093
−0,012
10000
10,024
10,024
10,024
10,020
10,024
10,024
10,023
−0,013
−0,016
−0,016
10,350
10,368
10,368
10,360
10,360
0,095
−0,076
−0,082
11,415
11,558
11,563
11,570
11,570
1,343
0,104
0,057
1000000
31,623
31,623
31,623
31,620
31,623
31,623
31,622
−0,004
−0,004
−0,004
31,634
31,635
31,635
31,630
31,630
−0,016
−0,014
−0,014
31,678
31,685
31,685
31,720
31,720
0,133
0,112
0,111In order to validate the model introduced in Section2, for the simply supported beam, a comparison of the first 3 natural frequencies is carried out for increasing values of N (the number of masses) and λ (the soil stiffness), and results are drawn in Table 1.The natural frequencies, Table1, converge faster by increasing the N or λ. A good approximation is already obtained for (N=10, λ =0) where the overall error ((average − present study (N))/average %) is below three percent limit. On the other hand, by increasing the soil stiffness to the case (N=10, λ = 1000000) the error drops below 0.15%. However, once the discretization number N is higher than 50 the overall error is under the 0.12%. Thus, increasing one of the parameters λ and N, or both, makes the discrete model converge to the reference models [10–12].


## 3.1.2. 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].Table 2
Comparison of the normalized first mode shape at different positions of the beam in the S-S beam case forN=100.
ξ
=
x
/
L
ExactRef. [13]
Present
FEMRef. [13]
Present
FEMRef. [13]
Present
FEMRef. [13]
Present
FEMRef. [13]
A
∗
=
0
A
∗
=
0
A
∗
=
0
A
∗
=
1
A
∗
=
1
A
∗
=
2
A
∗
=
2
A
∗
=
3
A
∗
=
3
0,00
0,0000
0,0000
0,0000
0,0000
0,0000
0,0000
0,0000
0,0000
0,0000
0,05
0,1564
0,1564
0,1564
0,1544
0,1564
0,1498
0,1564
0,1449
0,1564
0,10
0,3090
0,3090
0,3090
0,3054
0,3090
0,2969
0,3090
0,2876
0,3090
0,15
0,4540
0,4540
0,4540
0,4493
0,4540
0,4381
0,4540
0,4258
0,4540
0,20
0,5878
0,5878
0,5880
0,5827
0,5880
0,5707
0,5880
0,5570
0,5880
0,25
0,7071
0,7071
0,7071
0,7024
0,7071
0,6911
0,7071
0,6780
0,7071
0,30
0,8091
0,8090
0,8090
0,8053
0,8090
0,7962
0,8090
0,7853
0,8090
0,35
0,8910
0,8910
0,8910
0,8886
0,8910
0,8825
0,8910
0,8750
0,8910
0,40
0,9511
0,9511
0,9511
0,9498
0,9511
0,9468
0,9511
0,9430
0,9511
0,45
0,9877
0,9877
0,9877
0,9874
0,9877
0,9865
0,9877
0,9855
0,9877
0,50
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000The results are identical to the exact theoretical mode shape in the linear case. However, in the nonlinear case (high amplitudes), in spite of the result often accepted in the literature, according to which the mode shape in the simply supported case is amplitude independent [18], this discrete physical model leads to mode shapes which seem to be amplitude dependent. As a result, the amplitude dependency may be due to the hypothesis (small triangle angle) taken for the nonlinearity effect which is sensitive to high amplitudes.(2) Comparison of Nondimension Frequency
ω
n
l
s
s
∗
d
i
s
c
r
/
ω
l
s
s
∗
d
i
s
c
r
. For validating the results in the nonlinear case, comparing the ratio of the nonlinear frequency parameter (ωnlss∗discr/ωlss∗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 (N=100), because of the good approximation obtained for this degree of discretization in the linear case (error below 0.1%).Table3 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 (ωnlss∗discr/ωlss∗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.Table 3
Comparison of the first three modes’ nondimensional nonlinear frequencyωnlss∗discr/ωlss∗discr, in S-S beam, for various values of the nondimensional Winkler soil stiffness.
SIMPLY SUPPORTED BEAM
λ
0,1π4
1
π
4
10
π
4
A
∗
N
=
100
Ref. [14]
Error %
N
=
100
Ref. [14]
Error %
N
=
100
Ref. [14]
Error %
0,2
1,0026
1,0034
0,0843
1,0014
1,0019
0,0494
1,0003
1,0003
0,0044
0,4
1,0102
1,0135
0,3278
1,0056
1,0075
0,1876
1,0010
1,0014
0,0377
0,6
1,0228
1,0301
0,7128
1,0126
1,0167
0,4053
1,0023
1,0031
0,0799
0,8
1,0401
1,0529
1,2144
1,0223
1,0295
0,7036
1,0041
1,0054
0,1310
1
1,0620
1,0814
1,7928
1,0346
1,0457
1,0649
1,0064
1,0085
0,2110
2
1,2294
1,2924
4,8729
1,1319
1,1708
3,3188
1,0253
1,0334
0,7885
3
1,4666
1,5795
7,1497
1,2779
1,3519
5,4771
1,0560
1,0736
1,6424


## 3.2. 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 Section2.


### 3.2.1. 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:(61)KNCC=N+13EIL3+∞-410··0-46+α-410··1-4·····01····0·0···-410·01-46+α-40··01-4+∞.However, once introduced in the calculation, a mathematical error about the conditioning of the matrix arises (ill-conditioned 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 aty=0 as shown in Figure 4. Using this method for taking into account the clamped ends conditions, theN-DOF clamped beam problem is reduced to a (N-2)-DOF problem by removing the first and last row and column of the stiffness matrix; then (61) becomes(62)KN-2CC=N+13EIL36+α-410··0-46+α-410··1-4·····01····0·0···-410·01-46+α-40··01-46+α,KN-2CC=N+13EIL3KN-2CC∗.Figure 4
The model considered for the discrete C-C beam by blocking the rotation of the first and the last bars.With this new expression for the nondimensional rigidity matrixKN-2CC∗, the linear vibration equation (24) can be described by(63)KN-2CC∗y+ωlCCdiscr2ρSL4EIN+14y=0with y being the vector of mass displacements in MB. Solving (63) gives the eigenvalues β so one can write the frequencies as(64)ωlCCdiscr=N+12βiEIρSL4.Table4 gives the first three modes of the clamped beam nondimensional frequencies for increasing number of masses (N=10,50,100) and increasing values of the Winkler foundation stiffness. These nondimensional frequencies are calculated from (64) by extracting the EI/ρSL4 term, where the eigenvalues βi are those of the nondimensional stiffness matrix KN-2CC∗, so they can be presented as(65)ωlCCdiscr∗=N+12βi.Table 4
Comparison of the first three modes’ nondimension frequency and the nondimensional Winkler soil stiffness for a C-C beam.
C-C beam (linear)
λ
Present studyN=10
Present studyN=50
Present studyN=100
Ref. [11]
ExactRef. [12]
Error %N=10
Error %N=50
Error %N=100
0
5,2153
4,8251
4,7775
4,7300
4,7300
−10,261
−2,011
−1,003
8,5727
8,0083
7,9313
7,8500
7,8540
−9,150
−1,965
−0,984
11,7658
11,2052
11,1030
11,0000
10,9960
−7,001
−1,902
−0,973
100
5,3833
5,0338
4,9919
4,9500
4,9500
−8,753
−1,692
−0,846
8,6121
8,0565
7,9809
7,9000
7,9040
−8,959
−1,930
−0,973
11,7811
11,2229
11,1212
11,0100
11,0140
−6,965
−1,897
−0,974
10000
10,1800
10,1328
10,1278
10,1200
10,1230
−0,563
−0,097
−0,047
11,1400
10,8995
10,8692
10,8400
10,8390
−2,777
−0,558
−0,279
13,0681
12,6693
12,5991
12,5300
12,5260
−4,328
−1,144
−0,583
1000000
31,6286
31,6271
31,6269
31,6400
31,6260
−0,008
−0,003
−0,003
31,6654
31,6552
31,6540
31,6700
31,6530
−0,039
−0,007
−0,003
31,7732
31,7467
31,7422
31,7500
31,7380
−0,111
−0,027
−0,013In order to validate the discrete model, for clamped beams, Table4 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 massesN= (10/50/100), and the theoretical values of the natural frequencies of continuous beams obtained in [12].In conclusion, the natural frequencies drawn in Table4 converge slower to the continuous model [12] than the simply supported case. The error percentage ((results [12] − present study (N))/results [12]%), for the case of no soil stiffness applied λ = 0, do not get below the one percent until N=100. However, once the soil stiffness λ is increased, a good approximation is already obtained for the case of N=10. This convergence rate is due to the approximations taken by contraction of the matrices in the clamped beam case.


### 3.2.2. 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.Table 5
Comparison of the first mode shape at different positions of the beam in the C-C beam case.
C-C beam (nonlinear)
ξ
=
x
/
L
PresentN=100
ASMRef. [15]
FEMRef. [13]
PresentN=100
GFEMRef. [16]
FEMRef. [13]
PresentN=100
GFEMRef. [16]
FEMRef. [13]
A
∗
=
0
A
∗
=
2
A
∗
=
5
0,00
0,0000
0,0000
0,0000
0,0000
0,0000
0,0000
0,0000
0,0000
0,0000
0,05
0,0324
0,0325
0,0329
0,0341
0,0339
0,0341
0,0464
0,0411
0,0413
0,10
0,1196
0,1119
0,1203
0,1242
0,1234
0,1237
0,1555
0,1429
0,1432
0,15
0,2444
0,2435
0,2454
0,2508
0,2505
0,2506
0,2927
0,2795
0,2796
0,20
0,3906
0,3900
0,3922
0,3969
0,3981
0,3982
0,4363
0,4307
0,4308
0,25
0,5435
0,5435
0,5455
0,5480
0,5511
0,5513
0,5757
0,5815
0,5817
0,30
0,6895
0,6901
0,6918
0,6917
0,6960
0,6964
0,7058
0,7200
0,7204
0,35
0,8168
0,8178
0,8189
0,8174
0,8214
0,8219
0,8217
0,8372
0,8376
0,40
0,9156
0,9164
0,9170
0,9155
0,9182
0,9185
0,9159
0,9259
0,9262
0,45
0,9783
0,9787
0,9789
0,9782
0,9792
0,9293
0,9783
0,9812
0,9813
0,50
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000For clamped beams, the mode shapes are shown both theoretically and experimentally to be amplitude dependent at large vibration amplitudes. Table5 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
ω
n
l
c
c
∗
d
i
s
c
r
/
ω
l
c
c
∗
d
i
s
c
r. 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 (ωnlcc∗discr/ωlcc∗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 (N=100), for which a good convergence is reached in the linear case (error below 1%).Table 6
Comparison of first modes nondimensional nonlinear frequencyωnlcc∗discr/ωlcc∗discr, for various values of the nondimensional Winkler soil stiffness for a C-C beam.
Clamped-clamped beam
λ
0,1π4
1
π
4
10
π
4
A
∗
N
=
100
Ref. [14]
Error %
N
=
100
Ref. [14]
Error %
N
=
100
Ref. [14]
Error %
0,2
1,0013
1,0009
−0,0397
1,0011
1,0008
−0,0314
1,0005
1,0003
−0,0161
0,4
1,0052
1,0037
−0,1475
1,0044
1,0032
−0,1244
1,0018
1,0013
−0,0541
0,6
1,0116
1,0082
−0,3391
1,0100
1,0071
−0,2859
1,0041
1,0029
−0,1236
0,8
1,0206
1,0146
−0,5879
1,0177
1,0125
−0,5110
1,0073
1,0052
−0,2137
1
1,0320
1,0227
−0,9046
1,0275
1,0195
−0,7829
1,0115
1,0081
−0,3330
2
1,1224
1,0876
−3,1963
1,1058
1,0756
−2,8113
1,0451
1,0320
−1,2670
3
1,2587
1,1870
−6,0403
1,2253
1,1625
−5,4062
1,0988
1,0705
−2,6458The nonlinear behavior of a clamped beam resting on elastic foundations is presented in Table6 by means of the amplitude dependence of the parameter (ωnlcc∗discr/ωlcc∗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 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 λ.


## 3.2.1. 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:(61)KNCC=N+13EIL3+∞-410··0-46+α-410··1-4·····01····0·0···-410·01-46+α-40··01-4+∞.However, once introduced in the calculation, a mathematical error about the conditioning of the matrix arises (ill-conditioned 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 aty=0 as shown in Figure 4. Using this method for taking into account the clamped ends conditions, theN-DOF clamped beam problem is reduced to a (N-2)-DOF problem by removing the first and last row and column of the stiffness matrix; then (61) becomes(62)KN-2CC=N+13EIL36+α-410··0-46+α-410··1-4·····01····0·0···-410·01-46+α-40··01-46+α,KN-2CC=N+13EIL3KN-2CC∗.Figure 4
The model considered for the discrete C-C beam by blocking the rotation of the first and the last bars.With this new expression for the nondimensional rigidity matrixKN-2CC∗, the linear vibration equation (24) can be described by(63)KN-2CC∗y+ωlCCdiscr2ρSL4EIN+14y=0with y being the vector of mass displacements in MB. Solving (63) gives the eigenvalues β so one can write the frequencies as(64)ωlCCdiscr=N+12βiEIρSL4.Table4 gives the first three modes of the clamped beam nondimensional frequencies for increasing number of masses (N=10,50,100) and increasing values of the Winkler foundation stiffness. These nondimensional frequencies are calculated from (64) by extracting the EI/ρSL4 term, where the eigenvalues βi are those of the nondimensional stiffness matrix KN-2CC∗, so they can be presented as(65)ωlCCdiscr∗=N+12βi.Table 4
Comparison of the first three modes’ nondimension frequency and the nondimensional Winkler soil stiffness for a C-C beam.
C-C beam (linear)
λ
Present studyN=10
Present studyN=50
Present studyN=100
Ref. [11]
ExactRef. [12]
Error %N=10
Error %N=50
Error %N=100
0
5,2153
4,8251
4,7775
4,7300
4,7300
−10,261
−2,011
−1,003
8,5727
8,0083
7,9313
7,8500
7,8540
−9,150
−1,965
−0,984
11,7658
11,2052
11,1030
11,0000
10,9960
−7,001
−1,902
−0,973
100
5,3833
5,0338
4,9919
4,9500
4,9500
−8,753
−1,692
−0,846
8,6121
8,0565
7,9809
7,9000
7,9040
−8,959
−1,930
−0,973
11,7811
11,2229
11,1212
11,0100
11,0140
−6,965
−1,897
−0,974
10000
10,1800
10,1328
10,1278
10,1200
10,1230
−0,563
−0,097
−0,047
11,1400
10,8995
10,8692
10,8400
10,8390
−2,777
−0,558
−0,279
13,0681
12,6693
12,5991
12,5300
12,5260
−4,328
−1,144
−0,583
1000000
31,6286
31,6271
31,6269
31,6400
31,6260
−0,008
−0,003
−0,003
31,6654
31,6552
31,6540
31,6700
31,6530
−0,039
−0,007
−0,003
31,7732
31,7467
31,7422
31,7500
31,7380
−0,111
−0,027
−0,013In order to validate the discrete model, for clamped beams, Table4 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 massesN= (10/50/100), and the theoretical values of the natural frequencies of continuous beams obtained in [12].In conclusion, the natural frequencies drawn in Table4 converge slower to the continuous model [12] than the simply supported case. The error percentage ((results [12] − present study (N))/results [12]%), for the case of no soil stiffness applied λ = 0, do not get below the one percent until N=100. However, once the soil stiffness λ is increased, a good approximation is already obtained for the case of N=10. This convergence rate is due to the approximations taken by contraction of the matrices in the clamped beam case.


## 3.2.2. 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.Table 5
Comparison of the first mode shape at different positions of the beam in the C-C beam case.
C-C beam (nonlinear)
ξ
=
x
/
L
PresentN=100
ASMRef. [15]
FEMRef. [13]
PresentN=100
GFEMRef. [16]
FEMRef. [13]
PresentN=100
GFEMRef. [16]
FEMRef. [13]
A
∗
=
0
A
∗
=
2
A
∗
=
5
0,00
0,0000
0,0000
0,0000
0,0000
0,0000
0,0000
0,0000
0,0000
0,0000
0,05
0,0324
0,0325
0,0329
0,0341
0,0339
0,0341
0,0464
0,0411
0,0413
0,10
0,1196
0,1119
0,1203
0,1242
0,1234
0,1237
0,1555
0,1429
0,1432
0,15
0,2444
0,2435
0,2454
0,2508
0,2505
0,2506
0,2927
0,2795
0,2796
0,20
0,3906
0,3900
0,3922
0,3969
0,3981
0,3982
0,4363
0,4307
0,4308
0,25
0,5435
0,5435
0,5455
0,5480
0,5511
0,5513
0,5757
0,5815
0,5817
0,30
0,6895
0,6901
0,6918
0,6917
0,6960
0,6964
0,7058
0,7200
0,7204
0,35
0,8168
0,8178
0,8189
0,8174
0,8214
0,8219
0,8217
0,8372
0,8376
0,40
0,9156
0,9164
0,9170
0,9155
0,9182
0,9185
0,9159
0,9259
0,9262
0,45
0,9783
0,9787
0,9789
0,9782
0,9792
0,9293
0,9783
0,9812
0,9813
0,50
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000
1,0000For clamped beams, the mode shapes are shown both theoretically and experimentally to be amplitude dependent at large vibration amplitudes. Table5 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
ω
n
l
c
c
∗
d
i
s
c
r
/
ω
l
c
c
∗
d
i
s
c
r. 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 (ωnlcc∗discr/ωlcc∗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 (N=100), for which a good convergence is reached in the linear case (error below 1%).Table 6
Comparison of first modes nondimensional nonlinear frequencyωnlcc∗discr/ωlcc∗discr, for various values of the nondimensional Winkler soil stiffness for a C-C beam.
Clamped-clamped beam
λ
0,1π4
1
π
4
10
π
4
A
∗
N
=
100
Ref. [14]
Error %
N
=
100
Ref. [14]
Error %
N
=
100
Ref. [14]
Error %
0,2
1,0013
1,0009
−0,0397
1,0011
1,0008
−0,0314
1,0005
1,0003
−0,0161
0,4
1,0052
1,0037
−0,1475
1,0044
1,0032
−0,1244
1,0018
1,0013
−0,0541
0,6
1,0116
1,0082
−0,3391
1,0100
1,0071
−0,2859
1,0041
1,0029
−0,1236
0,8
1,0206
1,0146
−0,5879
1,0177
1,0125
−0,5110
1,0073
1,0052
−0,2137
1
1,0320
1,0227
−0,9046
1,0275
1,0195
−0,7829
1,0115
1,0081
−0,3330
2
1,1224
1,0876
−3,1963
1,1058
1,0756
−2,8113
1,0451
1,0320
−1,2670
3
1,2587
1,1870
−6,0403
1,2253
1,1625
−5,4062
1,0988
1,0705
−2,6458The nonlinear behavior of a clamped beam resting on elastic foundations is presented in Table6 by means of the amplitude dependence of the parameter (ωnlcc∗discr/ωlcc∗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 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 λ.


## 4. 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.


### 4.1. A Partially Supported Beam

In order to get a good approximation in the nonlinear case, the numberN 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 (λi≠0) and the nonsupported span (λi=0).


#### 4.1.1. 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 Figure5(a) and the discrete model of Figure 5(b).Figure 5
A partially supported continuous S-S beam.
(a)
(b)For a practical convenience, a semilogarithmic curve is drawn, showing the relative stiffness parameterλ on the abscissa and the nondimensional frequency parameter ωlss∗discr on the ordinate. The ratio of the supported span to the total span B [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 B=0, no-intermediate supports, and B=1, fully supported. The results correspond to the case N=100, which also enabled a simplification of B 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.Figure 6
Nondimensional frequencyωlSS∗discrλ for a simply supported beam: (a) first mode; (b) second mode; (c) third mode; (d) forth mode; (e) fifth mode.
(a)
(b)
(c)
(d)
(e)


#### 4.1.2. A Simply Supported Beam (Nonlinear Case)

After validating the physical model in the linear regime, one may investigate the nonlinearity effect on the NDFλωlSS∗discr.It may be concluded in Figure7(b), corresponding to a high nondimensional amplitude (A∗=3), that the difference between the curves is more accentuated than in Figure 7(a) (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 B=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 B=0.75, contrary to the case B=1, the gap between the curves seems to widen by increasing the soil stiffness. In the case of B=0.5, the same conclusion may be made as in the case B=0.75 but the gap seems to remain constant. For B=0.25, the same behavior is obtained as for B=0.5. In the case of B=0, the two curves are parallel.Figure 7
Nondimensional frequencyωlSS∗discr(λ) for a simply supported beam: (a) first mode (A∗=0/A∗=1); (b) first mode (A∗=0/A∗=3).
(a)
(b)


#### 4.1.3. 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 Figure8(a) to the discrete model of Figure 8(b).Figure 8
A partially supported C-C beam: (a) continuous model; (b) discrete model.
(a)
(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λωlCC∗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.Figure 9
Nondimensional frequencyωlCC∗discr(λ) for a clamped beam: (a) first mode; (b) second mode; (c) third mode; (d) forth mode; (e) fifth mode.
(a)
(b)
(c)
(d)
(e)


#### 4.1.4. A Clamped Beam (Nonlinear Case)

After establishing and validating the physical model in the linear regime in Section4.1.3, a comparison is carried out here for the nonlinear case. To show the impact of nonlinearity, the NDF(λ)ωlCC∗discr curves in the linear case (A∗=0) and in the nonlinear case (A∗≠0) are drawn in the same graph Figure 10.Figure 10
Nondimensional frequencyωlCC∗discr(λ) for a clamped beam: (a) first mode (A∗=0/A∗=1); (b) first mode (A∗=0/A∗=3).
(a)
(b)In Figure10(b), which corresponds to a high value of the nondimensional amplitude (A∗=3), the difference between the curves is more expressed in Figure 10(a) (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 B=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 B=0.75, contrary to the case B=1, the gap between the curves seems to widen by increasing the soil stiffness. In the case of B=0.5, the same conclusion may be made as in the case B=0.75 but the gap seems to remain constant. For B=0.25, the same behavior is obtained as for B=0.5. In the case of B=0, the two curves are parallel. However, the gaps between the curves are narrower than in the simply supported case.


### 4.2. 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. (66)klx=λ1-πx,(67)kpx=λ1-τx2.Equations (66) and (67) are introduced in the calculation algorithm and results are obtained.


#### 4.2.1. A Simply Supported Beam (Linear Case)

Simplification of variable implementation is the aim of the present discrete model. Figure11 shows a simply supported beam resting on a variable elastic foundation and Table 7 gives the variation of the frequency parameter Ωi associated with theith mode shape, for i=1 to 8, in accordance with the soil variation for the simply supported beam.Table 7
The first eight modes’ frequencies for the linear and parabolic soil distribution in S-S beam case.
Linearπ=0.2
λ
Ω
1
Ω
2
Ω
3
Ω
4
Ω
5
Ω
6
Ω
7
Ω
8
10
Ref. [17]
3,2118
6,2922
9,4275
12,5675
15,7085
18,8499
21,9914
25,1329
Present work
N
=
10
3,2010
6,2074
9,1420
11,8954
14,4077
16,6270
18,5079
20,0122
N
=
100
3,2116
6,2912
9,4240
12,5594
15,6927
18,8225
21,9479
25,0681
100
Ref. [17]
3,6999
6,3720
9,4515
12,5777
15,7138
18,8529
21,9933
25,1342
Present work
N
=
10
3,6884
6,2895
9,1681
11,9073
14,4144
16,6313
18,5111
20,0147
N
=
100
3,6993
6,3709
9,4481
12,5696
15,6980
18,8256
21,9499
25,0694
1000
Ref. [17]
5,6185
7,0420
9,6828
12,6783
15,7657
18,8831
22,0123
25,1469
Present work
N
=
10
5,6021
6,9748
9,4177
12,0242
14,4808
16,6747
18,5426
20,0396
N
=
100
5,6171
7,0406
9,6794
12,6702
15,7500
18,8558
21,9690
25,0822Figure 11
Model of S-S beam supported on variable soil. (a) Linear distribution; (b) parabolic distribution.
(a)
(b)The results obtained for the simply supported beam case, drawn in Table7 for a linear soil distribution and in Table 8 for a parabolic distribution, for the first eight modes (frequency parameter Ωi), 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 N=10.Table 8
The first eight modes’ frequencies for the linear and parabolic soil distribution in S-S beam case.
Parabolic(τ=0.2)
λ
Ω
1
Ω
2
Ω
3
Ω
4
Ω
5
Ω
6
Ω
7
Ω
8
10
Ref. [17]
3,2150
6,2926
9,4276
12,5675
15,7086
18,8499
21,9914
25,1329
Present work
N
=
10
3,2042
6,2078
9,1421
11,8955
14,4077
16,6270
18,5079
20,0122
N
=
100
3,2148
6,2916
9,4241
12,5594
15,6927
18,8226
21,9480
25,0681
100
Ref. [17]
3,7212
6,3755
9,4526
12,5781
15,7140
18,8530
21,9933
25,1342
Present work
N
=
10
3,7090
6,2928
9,1691
11,9077
14,4146
16,6315
18,5112
20,0148
N
=
100
3,7206
6,3744
9,4491
12,5700
15,6982
18,8257
21,9499
25,0694
1000
Ref. [17]
5,6788
7,0676
9,6923
12,6824
15,7679
18,8843
22,0131
25,1474
Present work
N
=
10
5,6600
6,9987
9,4269
12,0285
14,4832
16,6763
18,5438
20,0406
N
=
100
5,6772
7,0660
9,6888
12,6744
15,7521
18,8571
21,9697
25,0827


#### 4.2.2. A Simply Supported Beam (Nonlinear Case)

Figure12 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.Figure 12
Backbone curve nondimensional frequencyωnlSS∗discr/ωlSS∗discrA∗: (a) (τ=0.2, π=0.2, and λ=10); (b) (τ=0.2, π=0.2, and λ=1000).
(a)
(b)


#### 4.2.3. A Clamped Beam (Linear Case)

Figure13 shows a clamped beam resting on a variable elastic foundation. Table 9 indicates the variation of the frequency parameter Ωi associated with theith mode shape, with i=1 to 3, for the linear and parabolic soil stiffness variation, by increasing both the number of masses from N=10 to N=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 (N=10) is not sufficient to get a good approximation (error about 10%). For N=100, the results become accurate with an error around one percent.Table 9
The 3 first modes’ frequencies for the linear and parabolic soil distribution in C-C beam case.
λ
Linear(π=0.2)
Parabolic(τ=0.2)
Ω
1
Ω
2
Ω
3
Ω
1
Ω
2
Ω
3
1
Ref. [17]
4,7322
7,8537
10,9958
4,7323
7,8537
10,9958
Present work
N
=
10
4,2684
7,0143
9,6267
4,2685
7,0143
9,6267
N
=
100
4,6849
7,7746
10,8833
4,6850
7,7747
10,8833
10
Ref. [17]
4,7512
7,8579
10,9973
4,7522
7,8581
10,9974
Present work
N
=
10
4,2799
7,0169
9,6277
4,2805
7,0170
9,6278
N
=
100
4,7029
7,7786
10,8848
4,7039
7,7788
10,8848
100
Ref. [17]
4,9297
7,8993
11,0125
4,9391
7,9013
11,0132
Present work
N
=
10
4,3897
7,0427
9,6378
4,3953
7,0437
9,6381
N
=
100
4,8730
7,8180
10,8992
4,8820
7,8199
10,8999
1000
Ref. [17]
6,1172
8,2815
11,1611
6,1665
8,2988
11,1677
Present work
N
=
10
5,1964
7,2860
9,7362
5,2301
7,2957
9,7398
N
=
100
6,0167
8,1825
11,0407
6,0643
8,1989
11,0469Figure 13
Model of C-C beam supported by a soil with a variable stiffness: (a) linear distribution; (b) parabolic distribution.
(a)
(b)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.


#### 4.2.4. A Clamped Beam (Nonlinear Case)

The backbone curve presenting the nondimensional frequency is drawn for the case of a clamped beam in Figure14. 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.Figure 14
Backbone curve nondimensional frequencyωnlCC∗discr/ωlCC∗discrA∗: (a) (τ=0.2, π=0.2, and λ=10); (b) (τ=0.2, π=0.2, and λ=1000).
(a)
(b)


## 4.1. A Partially Supported Beam

In order to get a good approximation in the nonlinear case, the numberN 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 (λi≠0) and the nonsupported span (λi=0).


### 4.1.1. 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 Figure5(a) and the discrete model of Figure 5(b).Figure 5
A partially supported continuous S-S beam.
(a)
(b)For a practical convenience, a semilogarithmic curve is drawn, showing the relative stiffness parameterλ on the abscissa and the nondimensional frequency parameter ωlss∗discr on the ordinate. The ratio of the supported span to the total span B [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 B=0, no-intermediate supports, and B=1, fully supported. The results correspond to the case N=100, which also enabled a simplification of B 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.Figure 6
Nondimensional frequencyωlSS∗discrλ for a simply supported beam: (a) first mode; (b) second mode; (c) third mode; (d) forth mode; (e) fifth mode.
(a)
(b)
(c)
(d)
(e)


### 4.1.2. A Simply Supported Beam (Nonlinear Case)

After validating the physical model in the linear regime, one may investigate the nonlinearity effect on the NDFλωlSS∗discr.It may be concluded in Figure7(b), corresponding to a high nondimensional amplitude (A∗=3), that the difference between the curves is more accentuated than in Figure 7(a) (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 B=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 B=0.75, contrary to the case B=1, the gap between the curves seems to widen by increasing the soil stiffness. In the case of B=0.5, the same conclusion may be made as in the case B=0.75 but the gap seems to remain constant. For B=0.25, the same behavior is obtained as for B=0.5. In the case of B=0, the two curves are parallel.Figure 7
Nondimensional frequencyωlSS∗discr(λ) for a simply supported beam: (a) first mode (A∗=0/A∗=1); (b) first mode (A∗=0/A∗=3).
(a)
(b)


### 4.1.3. 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 Figure8(a) to the discrete model of Figure 8(b).Figure 8
A partially supported C-C beam: (a) continuous model; (b) discrete model.
(a)
(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λωlCC∗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.Figure 9
Nondimensional frequencyωlCC∗discr(λ) for a clamped beam: (a) first mode; (b) second mode; (c) third mode; (d) forth mode; (e) fifth mode.
(a)
(b)
(c)
(d)
(e)


### 4.1.4. A Clamped Beam (Nonlinear Case)

After establishing and validating the physical model in the linear regime in Section4.1.3, a comparison is carried out here for the nonlinear case. To show the impact of nonlinearity, the NDF(λ)ωlCC∗discr curves in the linear case (A∗=0) and in the nonlinear case (A∗≠0) are drawn in the same graph Figure 10.Figure 10
Nondimensional frequencyωlCC∗discr(λ) for a clamped beam: (a) first mode (A∗=0/A∗=1); (b) first mode (A∗=0/A∗=3).
(a)
(b)In Figure10(b), which corresponds to a high value of the nondimensional amplitude (A∗=3), the difference between the curves is more expressed in Figure 10(a) (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 B=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 B=0.75, contrary to the case B=1, the gap between the curves seems to widen by increasing the soil stiffness. In the case of B=0.5, the same conclusion may be made as in the case B=0.75 but the gap seems to remain constant. For B=0.25, the same behavior is obtained as for B=0.5. In the case of B=0, the two curves are parallel. However, the gaps between the curves are narrower than in the simply supported case.


## 4.1.1. 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 Figure5(a) and the discrete model of Figure 5(b).Figure 5
A partially supported continuous S-S beam.
(a)
(b)For a practical convenience, a semilogarithmic curve is drawn, showing the relative stiffness parameterλ on the abscissa and the nondimensional frequency parameter ωlss∗discr on the ordinate. The ratio of the supported span to the total span B [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 B=0, no-intermediate supports, and B=1, fully supported. The results correspond to the case N=100, which also enabled a simplification of B 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.Figure 6
Nondimensional frequencyωlSS∗discrλ for a simply supported beam: (a) first mode; (b) second mode; (c) third mode; (d) forth mode; (e) fifth mode.
(a)
(b)
(c)
(d)
(e)


## 4.1.2. A Simply Supported Beam (Nonlinear Case)

After validating the physical model in the linear regime, one may investigate the nonlinearity effect on the NDFλωlSS∗discr.It may be concluded in Figure7(b), corresponding to a high nondimensional amplitude (A∗=3), that the difference between the curves is more accentuated than in Figure 7(a) (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 B=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 B=0.75, contrary to the case B=1, the gap between the curves seems to widen by increasing the soil stiffness. In the case of B=0.5, the same conclusion may be made as in the case B=0.75 but the gap seems to remain constant. For B=0.25, the same behavior is obtained as for B=0.5. In the case of B=0, the two curves are parallel.Figure 7
Nondimensional frequencyωlSS∗discr(λ) for a simply supported beam: (a) first mode (A∗=0/A∗=1); (b) first mode (A∗=0/A∗=3).
(a)
(b)


## 4.1.3. 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 Figure8(a) to the discrete model of Figure 8(b).Figure 8
A partially supported C-C beam: (a) continuous model; (b) discrete model.
(a)
(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λωlCC∗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.Figure 9
Nondimensional frequencyωlCC∗discr(λ) for a clamped beam: (a) first mode; (b) second mode; (c) third mode; (d) forth mode; (e) fifth mode.
(a)
(b)
(c)
(d)
(e)


## 4.1.4. A Clamped Beam (Nonlinear Case)

After establishing and validating the physical model in the linear regime in Section4.1.3, a comparison is carried out here for the nonlinear case. To show the impact of nonlinearity, the NDF(λ)ωlCC∗discr curves in the linear case (A∗=0) and in the nonlinear case (A∗≠0) are drawn in the same graph Figure 10.Figure 10
Nondimensional frequencyωlCC∗discr(λ) for a clamped beam: (a) first mode (A∗=0/A∗=1); (b) first mode (A∗=0/A∗=3).
(a)
(b)In Figure10(b), which corresponds to a high value of the nondimensional amplitude (A∗=3), the difference between the curves is more expressed in Figure 10(a) (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 B=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 B=0.75, contrary to the case B=1, the gap between the curves seems to widen by increasing the soil stiffness. In the case of B=0.5, the same conclusion may be made as in the case B=0.75 but the gap seems to remain constant. For B=0.25, the same behavior is obtained as for B=0.5. In the case of B=0, the two curves are parallel. However, the gaps between the curves are narrower than in the simply supported case.


## 4.2. 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. (66)klx=λ1-πx,(67)kpx=λ1-τx2.Equations (66) and (67) are introduced in the calculation algorithm and results are obtained.


### 4.2.1. A Simply Supported Beam (Linear Case)

Simplification of variable implementation is the aim of the present discrete model. Figure11 shows a simply supported beam resting on a variable elastic foundation and Table 7 gives the variation of the frequency parameter Ωi associated with theith mode shape, for i=1 to 8, in accordance with the soil variation for the simply supported beam.Table 7
The first eight modes’ frequencies for the linear and parabolic soil distribution in S-S beam case.
Linearπ=0.2
λ
Ω
1
Ω
2
Ω
3
Ω
4
Ω
5
Ω
6
Ω
7
Ω
8
10
Ref. [17]
3,2118
6,2922
9,4275
12,5675
15,7085
18,8499
21,9914
25,1329
Present work
N
=
10
3,2010
6,2074
9,1420
11,8954
14,4077
16,6270
18,5079
20,0122
N
=
100
3,2116
6,2912
9,4240
12,5594
15,6927
18,8225
21,9479
25,0681
100
Ref. [17]
3,6999
6,3720
9,4515
12,5777
15,7138
18,8529
21,9933
25,1342
Present work
N
=
10
3,6884
6,2895
9,1681
11,9073
14,4144
16,6313
18,5111
20,0147
N
=
100
3,6993
6,3709
9,4481
12,5696
15,6980
18,8256
21,9499
25,0694
1000
Ref. [17]
5,6185
7,0420
9,6828
12,6783
15,7657
18,8831
22,0123
25,1469
Present work
N
=
10
5,6021
6,9748
9,4177
12,0242
14,4808
16,6747
18,5426
20,0396
N
=
100
5,6171
7,0406
9,6794
12,6702
15,7500
18,8558
21,9690
25,0822Figure 11
Model of S-S beam supported on variable soil. (a) Linear distribution; (b) parabolic distribution.
(a)
(b)The results obtained for the simply supported beam case, drawn in Table7 for a linear soil distribution and in Table 8 for a parabolic distribution, for the first eight modes (frequency parameter Ωi), 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 N=10.Table 8
The first eight modes’ frequencies for the linear and parabolic soil distribution in S-S beam case.
Parabolic(τ=0.2)
λ
Ω
1
Ω
2
Ω
3
Ω
4
Ω
5
Ω
6
Ω
7
Ω
8
10
Ref. [17]
3,2150
6,2926
9,4276
12,5675
15,7086
18,8499
21,9914
25,1329
Present work
N
=
10
3,2042
6,2078
9,1421
11,8955
14,4077
16,6270
18,5079
20,0122
N
=
100
3,2148
6,2916
9,4241
12,5594
15,6927
18,8226
21,9480
25,0681
100
Ref. [17]
3,7212
6,3755
9,4526
12,5781
15,7140
18,8530
21,9933
25,1342
Present work
N
=
10
3,7090
6,2928
9,1691
11,9077
14,4146
16,6315
18,5112
20,0148
N
=
100
3,7206
6,3744
9,4491
12,5700
15,6982
18,8257
21,9499
25,0694
1000
Ref. [17]
5,6788
7,0676
9,6923
12,6824
15,7679
18,8843
22,0131
25,1474
Present work
N
=
10
5,6600
6,9987
9,4269
12,0285
14,4832
16,6763
18,5438
20,0406
N
=
100
5,6772
7,0660
9,6888
12,6744
15,7521
18,8571
21,9697
25,0827


### 4.2.2. A Simply Supported Beam (Nonlinear Case)

Figure12 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.Figure 12
Backbone curve nondimensional frequencyωnlSS∗discr/ωlSS∗discrA∗: (a) (τ=0.2, π=0.2, and λ=10); (b) (τ=0.2, π=0.2, and λ=1000).
(a)
(b)


### 4.2.3. A Clamped Beam (Linear Case)

Figure13 shows a clamped beam resting on a variable elastic foundation. Table 9 indicates the variation of the frequency parameter Ωi associated with theith mode shape, with i=1 to 3, for the linear and parabolic soil stiffness variation, by increasing both the number of masses from N=10 to N=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 (N=10) is not sufficient to get a good approximation (error about 10%). For N=100, the results become accurate with an error around one percent.Table 9
The 3 first modes’ frequencies for the linear and parabolic soil distribution in C-C beam case.
λ
Linear(π=0.2)
Parabolic(τ=0.2)
Ω
1
Ω
2
Ω
3
Ω
1
Ω
2
Ω
3
1
Ref. [17]
4,7322
7,8537
10,9958
4,7323
7,8537
10,9958
Present work
N
=
10
4,2684
7,0143
9,6267
4,2685
7,0143
9,6267
N
=
100
4,6849
7,7746
10,8833
4,6850
7,7747
10,8833
10
Ref. [17]
4,7512
7,8579
10,9973
4,7522
7,8581
10,9974
Present work
N
=
10
4,2799
7,0169
9,6277
4,2805
7,0170
9,6278
N
=
100
4,7029
7,7786
10,8848
4,7039
7,7788
10,8848
100
Ref. [17]
4,9297
7,8993
11,0125
4,9391
7,9013
11,0132
Present work
N
=
10
4,3897
7,0427
9,6378
4,3953
7,0437
9,6381
N
=
100
4,8730
7,8180
10,8992
4,8820
7,8199
10,8999
1000
Ref. [17]
6,1172
8,2815
11,1611
6,1665
8,2988
11,1677
Present work
N
=
10
5,1964
7,2860
9,7362
5,2301
7,2957
9,7398
N
=
100
6,0167
8,1825
11,0407
6,0643
8,1989
11,0469Figure 13
Model of C-C beam supported by a soil with a variable stiffness: (a) linear distribution; (b) parabolic distribution.
(a)
(b)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.


### 4.2.4. A Clamped Beam (Nonlinear Case)

The backbone curve presenting the nondimensional frequency is drawn for the case of a clamped beam in Figure14. 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.Figure 14
Backbone curve nondimensional frequencyωnlCC∗discr/ωlCC∗discrA∗: (a) (τ=0.2, π=0.2, and λ=10); (b) (τ=0.2, π=0.2, and λ=1000).
(a)
(b)


## 4.2.1. A Simply Supported Beam (Linear Case)

Simplification of variable implementation is the aim of the present discrete model. Figure11 shows a simply supported beam resting on a variable elastic foundation and Table 7 gives the variation of the frequency parameter Ωi associated with theith mode shape, for i=1 to 8, in accordance with the soil variation for the simply supported beam.Table 7
The first eight modes’ frequencies for the linear and parabolic soil distribution in S-S beam case.
Linearπ=0.2
λ
Ω
1
Ω
2
Ω
3
Ω
4
Ω
5
Ω
6
Ω
7
Ω
8
10
Ref. [17]
3,2118
6,2922
9,4275
12,5675
15,7085
18,8499
21,9914
25,1329
Present work
N
=
10
3,2010
6,2074
9,1420
11,8954
14,4077
16,6270
18,5079
20,0122
N
=
100
3,2116
6,2912
9,4240
12,5594
15,6927
18,8225
21,9479
25,0681
100
Ref. [17]
3,6999
6,3720
9,4515
12,5777
15,7138
18,8529
21,9933
25,1342
Present work
N
=
10
3,6884
6,2895
9,1681
11,9073
14,4144
16,6313
18,5111
20,0147
N
=
100
3,6993
6,3709
9,4481
12,5696
15,6980
18,8256
21,9499
25,0694
1000
Ref. [17]
5,6185
7,0420
9,6828
12,6783
15,7657
18,8831
22,0123
25,1469
Present work
N
=
10
5,6021
6,9748
9,4177
12,0242
14,4808
16,6747
18,5426
20,0396
N
=
100
5,6171
7,0406
9,6794
12,6702
15,7500
18,8558
21,9690
25,0822Figure 11
Model of S-S beam supported on variable soil. (a) Linear distribution; (b) parabolic distribution.
(a)
(b)The results obtained for the simply supported beam case, drawn in Table7 for a linear soil distribution and in Table 8 for a parabolic distribution, for the first eight modes (frequency parameter Ωi), 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 N=10.Table 8
The first eight modes’ frequencies for the linear and parabolic soil distribution in S-S beam case.
Parabolic(τ=0.2)
λ
Ω
1
Ω
2
Ω
3
Ω
4
Ω
5
Ω
6
Ω
7
Ω
8
10
Ref. [17]
3,2150
6,2926
9,4276
12,5675
15,7086
18,8499
21,9914
25,1329
Present work
N
=
10
3,2042
6,2078
9,1421
11,8955
14,4077
16,6270
18,5079
20,0122
N
=
100
3,2148
6,2916
9,4241
12,5594
15,6927
18,8226
21,9480
25,0681
100
Ref. [17]
3,7212
6,3755
9,4526
12,5781
15,7140
18,8530
21,9933
25,1342
Present work
N
=
10
3,7090
6,2928
9,1691
11,9077
14,4146
16,6315
18,5112
20,0148
N
=
100
3,7206
6,3744
9,4491
12,5700
15,6982
18,8257
21,9499
25,0694
1000
Ref. [17]
5,6788
7,0676
9,6923
12,6824
15,7679
18,8843
22,0131
25,1474
Present work
N
=
10
5,6600
6,9987
9,4269
12,0285
14,4832
16,6763
18,5438
20,0406
N
=
100
5,6772
7,0660
9,6888
12,6744
15,7521
18,8571
21,9697
25,0827


## 4.2.2. A Simply Supported Beam (Nonlinear Case)

Figure12 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.Figure 12
Backbone curve nondimensional frequencyωnlSS∗discr/ωlSS∗discrA∗: (a) (τ=0.2, π=0.2, and λ=10); (b) (τ=0.2, π=0.2, and λ=1000).
(a)
(b)


## 4.2.3. A Clamped Beam (Linear Case)

Figure13 shows a clamped beam resting on a variable elastic foundation. Table 9 indicates the variation of the frequency parameter Ωi associated with theith mode shape, with i=1 to 3, for the linear and parabolic soil stiffness variation, by increasing both the number of masses from N=10 to N=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 (N=10) is not sufficient to get a good approximation (error about 10%). For N=100, the results become accurate with an error around one percent.Table 9
The 3 first modes’ frequencies for the linear and parabolic soil distribution in C-C beam case.
λ
Linear(π=0.2)
Parabolic(τ=0.2)
Ω
1
Ω
2
Ω
3
Ω
1
Ω
2
Ω
3
1
Ref. [17]
4,7322
7,8537
10,9958
4,7323
7,8537
10,9958
Present work
N
=
10
4,2684
7,0143
9,6267
4,2685
7,0143
9,6267
N
=
100
4,6849
7,7746
10,8833
4,6850
7,7747
10,8833
10
Ref. [17]
4,7512
7,8579
10,9973
4,7522
7,8581
10,9974
Present work
N
=
10
4,2799
7,0169
9,6277
4,2805
7,0170
9,6278
N
=
100
4,7029
7,7786
10,8848
4,7039
7,7788
10,8848
100
Ref. [17]
4,9297
7,8993
11,0125
4,9391
7,9013
11,0132
Present work
N
=
10
4,3897
7,0427
9,6378
4,3953
7,0437
9,6381
N
=
100
4,8730
7,8180
10,8992
4,8820
7,8199
10,8999
1000
Ref. [17]
6,1172
8,2815
11,1611
6,1665
8,2988
11,1677
Present work
N
=
10
5,1964
7,2860
9,7362
5,2301
7,2957
9,7398
N
=
100
6,0167
8,1825
11,0407
6,0643
8,1989
11,0469Figure 13
Model of C-C beam supported by a soil with a variable stiffness: (a) linear distribution; (b) parabolic distribution.
(a)
(b)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.


## 4.2.4. A Clamped Beam (Nonlinear Case)

The backbone curve presenting the nondimensional frequency is drawn for the case of a clamped beam in Figure14. 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.Figure 14
Backbone curve nondimensional frequencyωnlCC∗discr/ωlCC∗discrA∗: (a) (τ=0.2, π=0.2, and λ=10); (b) (τ=0.2, π=0.2, and λ=1000).
(a)
(b)


## 5. 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 N=10 to an error less than 0.1% for N=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 forN=10 and dropped to one percent for N=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 (A∗=1 and A∗=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.



---

*Source: 4740851-2017-03-15.xml*