The vibration characteristic of a Timoshenko beam resting on non-linear viscoelastic foundation subjected to any number of springs – mass systems (sprung masses) is governed by system of non – linear partial differential equations. The governing differential equations are examined using differential quadrature method to be transformed with boundary conditions into a set of algebraic equations. The non – linear Pasternak foundation is assumed to be cubic. Therefore, the effects of shear deformable beam and the shear deformation of foundations are considered at the same time. The numerical investigations show the dynamic response considering different values for engineering properties for both beam and foundation. Also, the numerical investigations show the efficiency and reliability of using differential quadrature method.

Keywords: beam, viscoelastic, sprung masses, differential quadrature, vlasove

Vibration analysis of beam type structures rested on a non – linear foundation has recently received a remarkable amount of attention due to importance and various applications of the subject. The differential quadrature method is an effective numerical technique for initial and boundary problems; it has not been applied to calculate nonlinear behaviors of Timoshenko beam rested on non-linear viscoelastic foundation. Application of the multiple scales method (MSM), method of Shaw and Pierre, method of normal forms and method of King and Vakakis in free vibration analysis of a simply supported beam rested on non – linear elastic foundation have been summarized by Nayfeh.¹ Ming-Hung Hsu² proposed new version differential quadrature method to obtain the vibration characteristics of rectangular plates resting on elastic foundations and carrying any number of sprung masses. The electrostatic behavior of the fixed-fixed beam type micro-actuators was simulated using the differential quadrature method by Ming-Hung Hsu.³ The vehicle load is one of the most important reasons for road damage. These pavements – vehicle systems can be theoretically modeled as beams supported by foundations subjected to moving forces of load similar to spring – mass. A novel state-space formulation was used by Giuseppe Muscolino & Alessandro Palmeri⁴ to scrutinize the response of beams resting on viscoelastically damped foundation under moving single degree of freedom oscillator.⁴ Li-Qun Chen⁵ paid special attentions to different nonlinear models and the introduction of the material time derivative into the viscoelastic constitutive relations. Iancu-Bogdan Teodoru & Vasile Musat⁶ applied Vlasov approach to beams resting on elastic supports. Davood Younesian et al.,⁷ solve the nonlinear governing differential equation of an elastic beam rested on a nonlinear foundation using Variational Iteration Method (VIM). Numerical solutions based on differential quadrature method were introduced for different structural problems.⁸ EJ Sapountzakis and AE Kampitsis⁹ developed boundary element method for the nonlinear dynamic analysis of beam-columns of an arbitrary doubly symmetric simply or multiply connected constant cross section, partially supported on a nonlinear three-parameter viscoelastic foundation. Galerkin method was used to find the response of a Timoshenko beam supported by a nonlinear foundation by Yan Yang et al.,¹⁰ also the convergence of this method was studied. Hu Ding et al.,¹¹ used the Adomian decomposition method and a perturbation method in conjunction with complex Fourier transformation to get the solution of the governing differential equations for Timoshenko beams with defined length supported by nonlinear viscoelastic foundations subjected to a moving concentrated force.

Consider a beam of length , with cross section of dimensions b x h carrying any number of sprung masses have masses mi and stiffness k_i and resting on viscoelastic Pasternak foundation as shown in the following Figure 1.

Pasternak foundation reaction
The foundation of the considered beam is taken as Pasternak foundation with linear and cubic stiffness and viscous damping:
$P\left(x,t\right)=K_{1}w\left(x,t\right)+K_3{w}^3\left(x,t\right)+\eta \frac{\partial w\left(x,t\right)}{\partial t}-G_P\frac{{\partial }^{2}w\left(x,t\right)}{\partial x^2}$ (1)
Where P(x,t) is the force induced by the foundation per unit length of the beam as a function of the horizontal coordinate x and time t, K₁ and K₃ are the first and third order foundation parameters, respectively. Furthermore G_P and $\eta$ are the shear deformation coefficient and damping coefficient of the foundation respectively. $w$ is the vertical displacement of the beam.

Beam strain energy
By considering Timoshenko beam theory, one can obtain the strain energy per unit length of beam element as:
$U_e=\frac{1}{2}{\displaystyle \underset{0}{\overset{L}{\int }}\left(EI{\left(\frac{\partial \theta }{\partial x}\right)}^2+kAG{\left(\frac{\partial w}{\partial x}-\theta \right)}^2\right)}dx$ (2)
Where $\theta$ is the rotation of the cross section, $\frac{\partial w}{\partial x}$ is the slope of the vertical displacement, E is the modulus of elasticity of the beam material, is the second moment of area, is the shear correction factor, is the cross section area and is the shear section modulus.

Beam Equation of Motion
The total strain energy of a beam resting on Pasternak foundation and carrying any numbers of sprung masses (oscillators) is:
$U=\frac{1}{2}{\displaystyle {\int }_0^L\left(EI{\left(\frac{\partial \theta }{\partial x}\right)}^2+kAG{\left(\frac{\partial w}{\partial x}-\theta \right)}^2+Pw\right)+{\displaystyle {\sum }_{i=1}^s\frac{1}{2}}}{k}_i{\left(y_i\left(t\right)-w\right)}^2$ (3)
Where y(t) is the vertical displacement of the oscillator and S is the number of sprung masses connected to the beam. The kinetic energy of the system can be expressed as:
$T=\frac{1}{2}\rho {\displaystyle {\int }_0^L\left(A{\left(\frac{\partial w}{\partial t}\right)}^2+I{\left(\frac{\partial \theta }{\partial t}\right)}^2+Pw\right)}+{{\displaystyle {\sum }_{i=1}^s\frac{1}{2}{m}_i\left(\frac{\partial y_i}{\partial t}\right)}}^2$ (4)
Where $\rho$ is the density of beam material. By applying Hamilton’s principle:
$\delta {\displaystyle {\int }_{t_1}^{t_2}\left(T-U\right)} dt=0$ (5)
From equations (3), (4) and (5), one can obtain:
$\rho A{w}^{\mathrm{..}}+\eta w^{.}-kAG\left(w^{\backslash \backslash }-{\theta }^{\backslash }\right)+K_{1}w+2K_3{w}^3-G_p{w}^{\backslash \backslash }-{\displaystyle {\sum }_{i=1}^s{k}_i\left(y_i-w\right)=0}$ (6)
$\rho I{\theta }^{\mathrm{..}}-kAG\left(w^{\backslash }-\theta \right)-EI{\theta }^{\backslash \backslash }=0$ 7)
${\displaystyle {\sum }_{i=1}^s{k}_i\left(y_i-w\right)-{\displaystyle {\sum }_{i=1}^s{m}_i{y}_i^{\mathrm{..}}}}=0$  (8)
Let $w\left(x,t\right)=W\left(x\right)e^{\lambda t},\theta \left(x,t\right)=\vartheta \left(x\right)e^{\lambda t}$  and $y_i\left(t\right)=Y_i{e}^{\lambda t}$
Then, equations (6), (7) and (8) yield:
$\left(\rho A{\lambda }^2+\eta \lambda +K_1+2K_3{W}^2{e}^{2\lambda t}\right)W+kAG\left({\vartheta }^{\backslash }\right)-\left(kAG+G_p\right)W^{\backslash \backslash }-{\displaystyle {\sum }_{i=1}^s{k}_i\left(Y_i-W\right)}=0$  (9)
$\rho I{\lambda }^2\vartheta -kAG\left(W^{\backslash }-\vartheta \right)-EI{\vartheta }^{\backslash \backslash }=0$  (10)
${\displaystyle {\sum }_{i=1}^s{k}_i}\left(Y_i-W\right)-{\displaystyle {\sum }_{i=1}^s{m}_i{\lambda }^2{Y}_i=0}$  (11)
For the following non-dimensional variables:
$w^{*}=\frac{w}{L}, y^{*}{}_i=\frac{Y_i}{L}and x^{*}=\frac{x}{L}$  (12)
${\alpha }_1=\frac{\left(\rho A{\lambda }^2+\eta \lambda +K_1+{\displaystyle {\sum }_{i=1}^s{k}_i}\right)L}{EA}$  (13)
${\alpha }_2=\frac{2K_3{L}^3}{EA}$ (14)
${\alpha }_3={\alpha }_4+\frac{G_P}{EA{L}^2}$  (15)
${\alpha }_4=\frac{kG}{EL}$  (16)
${\beta }_1=\frac{\rho {\lambda }^2{}_{L^2}}{E}$  (17)
${\beta }_2=\frac{kAG{L}^2}{EI}$  (18)
${\gamma }_1={\gamma }_2-\frac{m_i{\lambda }^{2}L}{EA}$  (19)
${\gamma }_2=\frac{k_{i}L}{EA}$ (20)

Equations (9), (10) and (11) can be rewritten as:
${\alpha }_1{w}^{*}+{\alpha }_2{e}^{2\lambda t}{w}^{*3}-{\alpha }_3{w}^{*\backslash \backslash }+{\alpha }_4{\vartheta }^{\backslash }-{\displaystyle {\sum }_{i=1}^s{\gamma }_2{y}^{*}{}_i=0}$  (21)
${\beta }_1\vartheta -{\beta }_2\left(W^{*\backslash }-\vartheta \right)-{\vartheta }^{\backslash \backslash }=0$ (22)
${\displaystyle {\sum }_{i=1}^s{\gamma }_1{y}^{*}{}_i-}{\displaystyle {\sum }_{i=1}^s{\gamma }_2{w}^{*}=0}$  (23)

Boundary conditions
The simply supported end conditions can be expressed as:
$w^{*}\left(0\right)=w^{*}\left(1\right)=0$ (24)
$w^{*}\left(0\right)=w^{*\backslash \backslash }\left(L\right)=0$  (25)

Differential quadrature technique
The method of DQ assumes that the function derivatives can be expressed as linear sum of the weighting coefficient times function value at all discrete points in the domain of the concerned variable, and then the function derivative can be written as:
$\frac{{\partial }^{m}f\left(x_i\right)}{\partial x^m}={\displaystyle \sum _{j=1}^N{C}_{ij}^{\left(m\right)}}f\left(x_j\right)$ (26)

Where:
$f\left(x_j\right)$ is the value of a function at a grid point x_j.

$C_{ij}^{\left(m\right)}$ is a weighting coefficient for the derivative of order (m).

By determining the weighting coefficients, the link between the derivatives and the functional values can be established.

By supposing that $f\left(x_j\right)$ is approximated by Fourier series expansion of the form:
$f\left(x_j\right)=c_o+{\displaystyle {\sum }_{k=1}^{N/2}\left(c_k\cos k{x}_j+d_k\sin k{x}_j\right)}$  (27)

Where: N is the number of grid points

By using the above test function, one can obtain explicit formulations to compute weighting coefficients of the first, second and higher order, where the diagonal elements of weighting coefficients are:
$C_{ii}^{\left(1\right)}=a_{ii}=-{\displaystyle {\sum }_{j=0,i\ne j}^N{a}_{ij}}$ (28)
$C_{ii}^{\left(2\right)}=b_{ii}=-{\displaystyle {\sum }_{j=0,i\ne j}^N{b}_{ij}}$ (29)
$C_{ii}^{\left(m\right)}=-{\displaystyle {\sum }_{j=0,i\ne j}^N{C}_{ij}^{\left(m\right)}}$  (30)

Also, the non-diagonal elements of weighting coefficients are:

$C_{ij}^{\left(1\right)}=a_{{}_{ij}=\frac{\alpha q\left({\xi }_i\right)}{2\sin \left(\frac{{\xi }_i-{\xi }_j}{2}\right)q\left({\xi }_j\right)}};i\ne j$  (31)
$C_{ij}^{\left(2\right)}=b_{ij}=a_{ij}\left[2a_{ii}-\cot \left(\frac{x_i-x_j}{2}\right)\right];i\ne j$ (32)
$C_{ij}^{\left(m\right)}=a_{ij}\left(\frac{1}{2}+m{b}_{ii}\right)-\frac{m}{2}{b}_{ij}\cot \frac{x_i-x_j}{2};i\ne j$ (33)

Where:
$C_{ij}^{\left(m\right)}=a_{ij}\left(\frac{1}{2}+m{b}_{ii}\right)-\frac{m}{2}{b}_{ij}\cot \frac{x_i-x_j}{2};i\ne j$ (34)
The above algebraic equations can be applied to periodic problems, i. e $\left(0\le x\le 2\pi \right)$  and non-periodic problems, i.e. $\left(0\le x\le \pi \right)$ . For practical applications the physical domain is not $\left[0,\pi \right]$  or $\left[0,2\pi \right]$ , but rather [a,b]. Then for this case, one can perform coordinates transformation from x – domain to – $\xi$ domain.
$C_{ij}^{\left(1\right)}=a_{ij}=\frac{\alpha q\left({\xi }_i\right)}{2\sin \left(\frac{{\xi }_i-{\xi }_j}{2}\right)q\left({\xi }_j\right)};i\ne j$  (35)
$C_{ij}^{\left(2\right)}=b_{ij}=a_{ij}\left[2a_{ii}-\alpha \cot \left(\frac{{\xi }_i-{\xi }_j}{2}\right)\right];i\ne j$  (36)
$C_{ij}^{\left(m\right)}=a_{ij}\left(\frac{{\alpha }^2}{2}+m{b}_{ii}\right)-\frac{m}{2}\alpha b_{ij}\cot \frac{{\xi }_i-{\xi }_j}{2};i\ne j$ (37)

Where:
$${\xi }_i=2\pi \frac{x_i-a}{b-a}, \alpha =\frac{2\pi }{b-a}$$ (38)

Grid points selection
Chebyshev- Gauss- Lobatto grid points were adopted by Shu and Chen (1999) as the accurate selection of the grid points. The coordinates of the grid points were chosen as:
$$x_i^{*}=\frac{1}{2}\left[1-\cos \left[\frac{i-1}{N-1}\pi \right]\right];i=1,2,\mathrm{....}N$$  (39)

Boundary conditions implementation
The Direct Substitution approach will be applied. The basic idea of this approach is implementing the function condition at the end points, while the derivative condition should be descritized by the DQ method. The descritized Neumman conditions at the two boundaries are then combined to get the $W_2,W_{\left(N-1\right)}$ in terms of W₃,W₄,….,W_(N-2). The dimension of the equation system using this technique is $\left(N-4\right)\times \left(N-4\right)$

For any clamped and simply supported conditions, the descritized end conditions using the DQ method can be expressed as:
${\displaystyle \sum _{k=1}^N{C}_{1k}^{\left(n0\right)}}{W}_k^{*}=0$ (40)
${\displaystyle \sum _{k=1}^N{C}_{1k}^{\left(n1\right)}}{W}_k^{*}=0$  (41)
Where $\left(n0\right),\left(n1\right)$ can be written as 1or 2. By selecting the values of and, one can get the following sets of end conditions:

n0 = 2, n1 = 2 ……….simply supported ----- simply supported

By substitution in equations (40), (41), one can couple these equations together to give

W₂, W_(N-1) as:
$W_2^{*}=\frac{1}{AXN}{\displaystyle \sum _{k=3}^{N-2}AXK1W_k^{*}}$  (42)

$W_{N-1}^{*}=\frac{1}{AXN}{\displaystyle \sum _{k=3}^{N-2}AXKN{W}_k^{*}}$ (43)
where
$AXKN1=C_{1,k}^{\left(n0\right)}.C_{N,N-1}^{\left(n1\right)}-C_{1,N-1}^{\left(n0\right)}.C_{N,k}^{\left(n1\right)}$  (44)
$AXKN=C_{1,2}^{\left(n0\right)}.C_{N,k}^{\left(n1\right)}-C_{1,k}^{\left(n0\right)}.C_{N,k}^{\left(n1\right)}$  (45)
$AXN=C_{N,2}^{\left(n1\right)}.C_{1,N-1}^{\left(n0\right)}-C_{1,2}^{\left(n0\right)}.C_{N,N-1}^{\left(n1\right)}$  (46)
Hence $W_2^{*}$ , $W_{(N-1)}^{*}$ are introduced in terms of $W_3^{*}$ , $W_4^{*}$ ,..., $W_{(N-2)}^{*}$ , to be smoothly inserted into introduced into discretized from of the governing equations (21), (22) and (23) to be applied at (N-4_ grid points, then the matrices of the weighting coefficients can be obtained from

$C_1=C_{i,k}^{\left(2\right)}-\frac{C_{i,2}^{\left(2\right)}AXK1+C_{i,N-1}^{\left(2\right)}.AXKN}{AXN}$ (47)
$C_{m-1}=C_{i,k}^{\left(m\right)}-\frac{C_{i,2}^{\left(m\right)}AXK1+C_{i,N-1}^{\left(m\right)}.AXKN}{AXN}$  (48)
Where:
C₁ is a new weighting coefficient for second order derivative.
Cm_-1 is a new weighting coefficient for m^th order derivative

Numerical results
The introduced problem with differential quadrature solution was verified with the model presented by Y Yang et al.¹⁰ The considered values for geometric and engineering properties of beam, foundation and sprung masses load are shown in the following Table 1-3. The transverse deflection is plotted versus the longitudinal coordinate (x) considering one oscillator as shown in the following Figure 1. Good agreement between proposed solution for 13 Chebyshev- Gauss- Lobatto grid points and the solution presented by Y. Yang et al.¹⁰ for 200-term Galerkin truncation, considering one oscillator, was shown in above Figure 2. The transverse deflection of the beam was investigated for first three modes $\lambda =1$ , $\lambda =2$ and 3, considering 13 sprung masses at time As shown in Figure 3 the trend of the curve is the same due to the considered number of sprung masses and the central deflection increases as the mode number increases. The effect of changing modulus of elasticity of beam material, consequently changing the shear modulus of the material, on the transverse deflection of the beam was carried out. As shown in Figure 4, the central deflection increases as both of the modulus of elasticity and shear modulus increases. The effects of both linear and non-linear foundation parameters and are studied. As shown in (Figure 5) & (Figure 6), as the linear foundation parameter increases the transverse deflection increases but as the nonlinear foundation parameter increases the transverse deflection decreases. Finally, the effect of the Pasternak shear deformation coefficient is investigated. As shown in Figure 7, as the shear deformation increases the transverse deflection increases.

Figure 1 Timoshenko Beam Rested on Viscoelastic Fondation.

Figure 2 Verification of Presented Solution.

Figure 3 Transverse Deflection of the Beam for First Three Modes.

Figure 4 Effects of Both Modulus of Elasticity and Shear Modulus on Transverse Deflection of the Beam.

Figure 5 Effect of Linear Foundation Parameter K₁ on Transverse Deflection of the Beam.

Figure 6 Effect of Nonlinear Foundation Parameter K₃ on Transverse Deflection of the Beam.

Figure 7 Effect of Pasternak Shear Deformation Coefficient (G^p) on Transverse Deflection of the Beam.

Differential quadrature method is an effective numerical technique can be applied to calculate nonlinear behaviors of Timoshenko beam rested on non-linear viscoelastic foundation. Good agreement between differential quadrature technique using 13 grid points and the Galerkin truncation method for 200 terms that reflect efficiency and reliability of differential quadrature method for this non-linear problem. Also differential quadrature gives the availability of considering any number of sprung masses. The numerical investigation shows that both linear and non-linear foundation parameter have more considerable effects on beam transverse deflection than Pasternak shear deformation coefficient.

A is the beam cross section area.
b is the width of beam cross section.
C_ij^(m)is a weighting coefficient for the derivative of order (m).
E is the modulus of elasticity of the beam material.
G is the beam shear modulus.
G_p is the shear deformation coefficient of the foundation.
h is the height of beam cross section.
I is the second moment of area.
K₁ is the linear foundation parameter.
Ka is the non-linear foundation parameter.
k is the shear correction factor.
k_i is the stiffness of sprung masses.
L is the beam length.
m_i is the mass of sprung masses.
N is the number of grid points.
t is the time.
w is the vertical displacement of the beam.
x is the horizontal coordinate.
$\eta$ is the damping coefficient of the foundation.
$\theta$ is the rotation of the cross section.

None.

The authors declare no conflict of interest.

1. Ali Hasan Nayfeh. Introduction to perturbation techniques. A Wiley-Inter science Publication; 1993. p. 1-532.
2. Ming-Hung Hsu. Vibration characteristics of rectangular plates resting on elastic foundations and carrying any number of sprung masses. International Journal of Applied Science and Engineering. 2006;4(1):83-89.
3. Ming-Hung Hsu. Nonlinear deflection analysis of electrostatic micro-actuators with different electrode and beam shapes. Iranian Journal of Electrical and Computer Engineering. 2007;6(1):73-79.
4. Giuseppe Muscolino, Alessandro Palmeri. Response of beams resting on viscoelastically damped foundation to moving oscillators. International Journal of Solids and Structures. 2007;44(5):1317–1336.
5. Li-Qun Chen. Nonlinear dynamics. Croatia: INTECH; 2010. 366 p.
6. Iancu-Bogdan Teodoru and Vasile Musat. The modified vlasov foundation model: an attractive approach for beams resting on elastic supports. Electronic Journal of Geotechnical Engineering. 2010;15:1-13.
7. Davood Younesian, Zia Saadatnia, Hassan Askari. Analytical solutions for free oscillations of beams on nonlinear elastic foundations using the variational iteration method. Journal of Theoretical and Applied Mechanics. 2012;50(2):639-652.
8. Ahmad Salah Edeen Nassef. Structural applications on differential quadrature method. Lap Lambert academic publishing. 2012.
9. EJ Sapountzakis, AE Kampitsis. Nonlinear dynamic analysis of shear deformable beam-columns on nonlinear three-parameter viscoelastic foundation. I: Theory and numerical implementation. Journal of Engineering Mechanics. 2013;139(7):886-896.
10. Yan Yang, Hu Ding, Li-Qun Chen. Dynamic response to a moving load of a timoshenko beam resting on a nonlinear viscoelastic foundation. Acta Mechanica Sinica. 2013;29(5):718–727.
11. Hu Ding, Kang-Li Shi, Li-Qun Chen, et al. Dynamic response of an infinite Timoshenko beam on a nonlinear viscoelastic foundation to a moving load; Nonlinear Dynamics. 2013;73(1-2):285–298.