# Analysis of Tapered Timoshenko and Euler–Bernoulli Beams on an Elastic Foundation with Moving Loads

Authors: W. Abbas; Omar K. bakr; M. M. Nassar; Mostafa A. M. Abdeen; M. Shabrawy
Journal: Journal of Mathematics
 (2021)

---

## Abstract

This research studies the vibration analysis of Euler–Bernoulli and Timoshenko beams utilizing the differential quadrature method (DQM) which has wide applications in the field of basic vibration of different components, for example, pillars, plates, round and hollow shells, and tanks. The free vibration of uniform and nonuniform beams laying on elastic Pasternak foundation will be studied under three sets of boundary conditions, that is, mixing between being simply upheld and fixed while utilizing the DQM. The natural frequencies and deflection values were produced through the examination of both beam types. Results show great concurrence with solutions from previous research studies. The impact of the nonuniform cross-section area on the vibration was contemplated. A comparison between the results from both beams is obtained. The focus of this work is on studying the deflection difference between both beam theories at different beam dimensions as well as showing the shape of rotation of the cross section while applying a nodal point load equation to simulate the moving load. The results were discussed and a general contemplation about the theories was developed.


---

## Body

## 1. Introduction

The most basic assumption developed about a bending beam resting on an elastic foundation is that the foundation’s reaction forces are directly proportional to the beam deflection at any arbitrary point. The deformation simulation of the foundation due to loading is perceived as independent, similar, and linear elastic springs. This trivial consideration was presented in 1867 by Emil Winkler [1]. This analysis has been widely accepted and has been developed to include more complicated parameters throughout the years. Oppermann [2] studied the beam on elastic foundation using the simplified continuum approach. Avramidis and Morfidis [3] introduced three-parameter elastic foundation for Kerr-type Timoshenko beam. Nassef et al. [4] studied the effect of using spring loads on a Timoshenko beam supported on a viscoelastic damped spring foundation. The dynamic response of a finite beam resting on a Winkler subjected to a moving point load as well as a moving vehicle on two axles moving with a speed varying with time was studied by Beskou and Muho [5].Euler–Bernoulli beam theory is an introductory method to calculate bending of beams when a load is applied. Bernoulli formulated the differential equation of motion of a vibrating beam. After that, Euler studied elastic beams under different loading conditions’ deforms. The Euler–Bernoulli beam theory is the most applied of other theories because of its simple and sensible numerical approximations in many problems, which stated that the model tends to overestimate the natural frequencies slightly. Shu and Du [6] studied the vibration of an Euler beam while applying the clamped and simply supported boundary conditions on the beam. Rajib et al. [7] studied the dynamic response of an Euler–Bernoulli beam subjected to a moving load supported by the Pasternak foundation. Mutman [8] presented a study about the free vibration of an Euler beam with a variable width using the Homotopy perturbation method. Abumandour [9] studied the vibration of a nonuniform Euler–Bernoulli beam resting on two-layer elastic foundation. Bezerra et al. [10] presented and studied the free vibration of an Euler–Bernoulli beam resting on Pasternak foundation.Timoshenko [11] contributed a theory which takes into account the shear and rotary inertia corrections that are neglected in Euler–Bernoulli’s beam theory. The solution of simply supported beam subjected to moving loads through using the power series expansion was studied by Timoshenko [12]. Leszek [13] studied the vibration of uniform Euler and Timoshenko beam resting on an elastic foundation using the green functions. Azam et al. [14] studied the response of Timoshenko beams subject to moving mass using Hamilton’s principle. Lucas et al. [15] studied the dynamic analysis of a Timoshenko beam resting on a two-parameter Pasternak foundation using finite element. Li and Hu [16] investigated the nonlinear bending and free vibration behaviors of size-dependent nonlinear Euler–Bernoulli and Timoshenko beam through thickness power-law variation of two-constituent functionally graded. The cross-section effect on the bending behaviors of nonlocal beams was studied by Li et al. [17]. Mirzaeiet al. [18] studied, using a semianalytical method, the thermoelastic behavior of temperature dependent and independent functionally graded variable thickness cantilever beam.Researching the development of new numerical techniques for solving engineering problems is on par with the rapidly growing advancement of faster computer machines. The differential quadrature method was first developed by Bellman [19]. The DQM is a numerical solution technique which is noticeably growing to be used in solving the boundary and initial value problems; it is considered a well-suited alternative to the well-known conventional numerical methods such as the finite difference technique. Patil and Kadoli [20] utilized the DQM to study the influence of the Winkler foundation on the free vibration of a functionally graded beam.In this paper, a numerical analysis of Timoshenko beams and Euler–Bernoulli beams is presented, using the differential quadrature method, through which the natural frequencies and deflections due to external loads are calculated. The output results using the DQM are compared with the exact solution as well as the work of other researchers to validate the accuracy of the work. Both beam theories resting on a two-parameter Pasternak foundation are compared to describe the behavior of both theories under the same conditions.The contribution of this work is the study of both beam theories simultaneously while using a variable beam depth which subjected to a moving load and an axial load. The moving load was applied without using the Dirac delta function, instead a relation between the time and grid point was formed which interprets the location of the vertical load by inserting the time parameter; this facilitates the application of the vertical load without using an analytical solution.


## 2. Problem Formulation

Figure1 shows the schematic diagram of a beam of length L, carrying a constant traversing point load f, and a compressive point load P at the beam’s ends. The beam rests on an elastic Pasternak foundation.Figure 1
Physical problem.


### 2.1. The Euler–Bernoulli Model

The assumptions stated by the Euler–Bernoulli beam theory are [21]Beam sections remain plain after deformationDeformed angles of rotation of the neutral axis are smallThe foundation of the beam is taken as a Pasternak foundation with linear stiffness:(1)Rx,t=k1wx,t−GR∂2wx,t∂x2,where Rx,t is the reaction of the foundation per unit length, wx,t is the beam vertical displacement, k1 is the first order foundation parameter (elastic stiffness), and GR is the shear deformation coefficient.Basic energy derivations for the Euler–Bernoulli beam model are(2)M=EI∂2w∂x2.The beam bending energyUb is given as follows:(3)Ub=12∫0LEI∂2w∂x22dx.The system’s kinetic energyT is given as follows:(4)T=12∫0LρA∂w∂t2dx,where ρ is the beam density. The work done by the compressive force P at the beam’s ends is(5)Wp=P2∫0L∂w∂x2dx.Applying Hamilton’s Principle,(6)δ∫t1t2∫0LLdxdt=δ12∫0LρA∂w∂t2dx−12∫0LEI∂2w∂x22dx−12∫0LKdx+P2∫0L∂w∂x2dxdxdt,where L is the Lagrangian. Solving equation (6) results in the following fourth-order differential equation:(7)∂2∂x2EIwxx+Pwxx+ρAwtt+k1w−Gpwxx=fx,t.To neglect the effect of time, we use the method of separation of variables and assume that(8)wx,t=eλt∗ωx.Substituting then simplifying gives(9)∂2∂x2EI∂2ω∂x2+P∂2ω∂x2+ρAωλ2+k1∗ωx−Gp∂2ω∂x2=fx,teλt∗L3EI0.Introducing the dimensionless variables,(10)W=ωL,X=xL⟶ω=WL,x=XL.Substituting the newly stated variables into equation (9) yields after simplification:(11)∂2∂X2SX∂2W∂X2+β1W+β4∂2W∂X2=Fx,t+Ω2W,where(12)β1=k1L4EI0,β2=GpL2EI0,β3=PL2EI0,β4=β3−β2,andΩ2=ρAL4λ2EI0.The following two boundary condition types are considered.A simply upheld end:(13)W=0,∂2W∂X2=0.A fixed end:(14)W=0,∂W∂X=0.Assuming a fixed and a simply upheld ends, end conditions can be shown, respectively:(15)W1=0,∂W∂X1=0,WN=0,∂2W∂XN2=0.


### 2.2. The Timoshenko Beam Model

From the Timoshenko theory assumptions,(16)ux=−zφ,uy=0,uz=w,εx=duxdx=−z∗∂φ∂x,where φ is the rotation of the cross section (which is different from the rotation of the neutral axis in Timoshenko beams). This difference results in additional rotation due to shear:(17)γs=∂w∂x−φ.According to (16) and (17), the strain components are given as(18)ε=−z∂φ∂x∂w∂x−φ.From Hook’s Law,(19)σ=E.ε=EI00KAG.−z∂φ∂x∂w∂x−φ,(20)∴σ=−EIz∂φ∂xKAG∂w∂x−φ,the beam’s bending energy U is calculated as(21)U=12∫VσT.εdV.Substituting (18) and (20) into (21) results in(22)U=12∫0L−EIz∂φ∂xKAG∂w∂x−φT−z∂φ∂x∂w∂x−φ=12∫0LEI∂φ∂x2+KAG∂w∂x−φ2dx.As a result, the beam bending strain energy resting on a Pasternak foundation is expressed as(23)U=12∫0LEI∂φ∂x2+KAG∂w∂x−φ2dx,where φ is the cross-section rotation, ∂w/∂x is the slope of vertical displacement, E is the beam’s modulus of elasticity, I is the moment of inertia, A is the cross-section area, G is the shear correction factor, K is the shear section modulus, and AG is the shear modulus of the cross-section area. The system’s kinetic energy T is given as follows:(24)T=12ρ∫0LA∂w∂t2+I∂φ∂t2dx,where ρ is the density of the beam material. Applying Hamilton’s Principle,(25)δ∫t1t2Ldt=0⟶∫t1t2U−T+Wdt=0,where L is the Lagrangian:(26)δ∫t1t2∫0LLdxd=δ∫t1t2∫0L12EI∂φ∂x2+KAG∂w∂x−φ2−12ρA∂w∂t2+I∂φ∂t2+P2∂w∂x2dxdt=0.Solving the system, the resulting coupled linear diff. equations are(27)−∂∂xEIφx−KAGwx−φ+ρIφtt=0,(28)−KAG∂∂xwx−φ+ρAwtt−pwxx+k1w−Gpwxx=fx,t.To neglect the effect of time, assume that(29)φx,t=eλt∗Фx,wx,t=eλt∗ωx.Substituting the previous expansions into (27) and (28), respectively, yields(30)−∂∂xEI∂Ф∂xeλt−KAG∂ω∂xeλt−eλtФ+ρIλ2eλtФ=0,(31)−KAG∂∂x∂ω∂xeλt−eλtФ+ρAωλ2eλt−p∂2ω∂x2eλt+k1eλtω−Gp∂2ω∂x2∗eλt=fx,t.Reintroducing the dimensionless variables,(32)W=ωL,X=xL.Substituting in (30),(33)−β1∂∂XSX∂Ф∂X−β2∂W∂X−Ф=SXΩ2Ф.Substituting in (31),(34)−β3∂2W∂X2+β3∂Ф∂X+WΩ2−β4∂2W∂X2+β7W−β5∂2W∂X2=F,factoring for Ф,W,andWxx, equation (34) for a variable cross-section becomes(35)β3∂Ф∂X−β6∂2W∂X2+β7W=F+WΩ2,where(36)SX=EIEI0=II0=1+αxL3,A=A01+αxL,α=hL−h0hL,β1=AL2I0,β2=KA2GL4EI02,β3=KAGLEI0,β4=PL2EI0,β5=GpL2EI0,β6=β3+β4+β5,β7=K1L4EI0,Ω2=ρAL4λ2EI0,F=fx,teλtL3EI0.The following two boundary condition types are considered.A simply upheld end:(37)W=0,∂2Ф∂X2=0.A fixed end:(38)W=0,∂Ф∂X=0.Assuming fixed and simply upheld ends, end conditions can be shown, respectively, as(39)W1=0,∂Ф∂X1=0,(40)WN=0,∂2Ф∂XN2=0.


## 2.1. The Euler–Bernoulli Model

The assumptions stated by the Euler–Bernoulli beam theory are [21]Beam sections remain plain after deformationDeformed angles of rotation of the neutral axis are smallThe foundation of the beam is taken as a Pasternak foundation with linear stiffness:(1)Rx,t=k1wx,t−GR∂2wx,t∂x2,where Rx,t is the reaction of the foundation per unit length, wx,t is the beam vertical displacement, k1 is the first order foundation parameter (elastic stiffness), and GR is the shear deformation coefficient.Basic energy derivations for the Euler–Bernoulli beam model are(2)M=EI∂2w∂x2.The beam bending energyUb is given as follows:(3)Ub=12∫0LEI∂2w∂x22dx.The system’s kinetic energyT is given as follows:(4)T=12∫0LρA∂w∂t2dx,where ρ is the beam density. The work done by the compressive force P at the beam’s ends is(5)Wp=P2∫0L∂w∂x2dx.Applying Hamilton’s Principle,(6)δ∫t1t2∫0LLdxdt=δ12∫0LρA∂w∂t2dx−12∫0LEI∂2w∂x22dx−12∫0LKdx+P2∫0L∂w∂x2dxdxdt,where L is the Lagrangian. Solving equation (6) results in the following fourth-order differential equation:(7)∂2∂x2EIwxx+Pwxx+ρAwtt+k1w−Gpwxx=fx,t.To neglect the effect of time, we use the method of separation of variables and assume that(8)wx,t=eλt∗ωx.Substituting then simplifying gives(9)∂2∂x2EI∂2ω∂x2+P∂2ω∂x2+ρAωλ2+k1∗ωx−Gp∂2ω∂x2=fx,teλt∗L3EI0.Introducing the dimensionless variables,(10)W=ωL,X=xL⟶ω=WL,x=XL.Substituting the newly stated variables into equation (9) yields after simplification:(11)∂2∂X2SX∂2W∂X2+β1W+β4∂2W∂X2=Fx,t+Ω2W,where(12)β1=k1L4EI0,β2=GpL2EI0,β3=PL2EI0,β4=β3−β2,andΩ2=ρAL4λ2EI0.The following two boundary condition types are considered.A simply upheld end:(13)W=0,∂2W∂X2=0.A fixed end:(14)W=0,∂W∂X=0.Assuming a fixed and a simply upheld ends, end conditions can be shown, respectively:(15)W1=0,∂W∂X1=0,WN=0,∂2W∂XN2=0.


## 2.2. The Timoshenko Beam Model

From the Timoshenko theory assumptions,(16)ux=−zφ,uy=0,uz=w,εx=duxdx=−z∗∂φ∂x,where φ is the rotation of the cross section (which is different from the rotation of the neutral axis in Timoshenko beams). This difference results in additional rotation due to shear:(17)γs=∂w∂x−φ.According to (16) and (17), the strain components are given as(18)ε=−z∂φ∂x∂w∂x−φ.From Hook’s Law,(19)σ=E.ε=EI00KAG.−z∂φ∂x∂w∂x−φ,(20)∴σ=−EIz∂φ∂xKAG∂w∂x−φ,the beam’s bending energy U is calculated as(21)U=12∫VσT.εdV.Substituting (18) and (20) into (21) results in(22)U=12∫0L−EIz∂φ∂xKAG∂w∂x−φT−z∂φ∂x∂w∂x−φ=12∫0LEI∂φ∂x2+KAG∂w∂x−φ2dx.As a result, the beam bending strain energy resting on a Pasternak foundation is expressed as(23)U=12∫0LEI∂φ∂x2+KAG∂w∂x−φ2dx,where φ is the cross-section rotation, ∂w/∂x is the slope of vertical displacement, E is the beam’s modulus of elasticity, I is the moment of inertia, A is the cross-section area, G is the shear correction factor, K is the shear section modulus, and AG is the shear modulus of the cross-section area. The system’s kinetic energy T is given as follows:(24)T=12ρ∫0LA∂w∂t2+I∂φ∂t2dx,where ρ is the density of the beam material. Applying Hamilton’s Principle,(25)δ∫t1t2Ldt=0⟶∫t1t2U−T+Wdt=0,where L is the Lagrangian:(26)δ∫t1t2∫0LLdxd=δ∫t1t2∫0L12EI∂φ∂x2+KAG∂w∂x−φ2−12ρA∂w∂t2+I∂φ∂t2+P2∂w∂x2dxdt=0.Solving the system, the resulting coupled linear diff. equations are(27)−∂∂xEIφx−KAGwx−φ+ρIφtt=0,(28)−KAG∂∂xwx−φ+ρAwtt−pwxx+k1w−Gpwxx=fx,t.To neglect the effect of time, assume that(29)φx,t=eλt∗Фx,wx,t=eλt∗ωx.Substituting the previous expansions into (27) and (28), respectively, yields(30)−∂∂xEI∂Ф∂xeλt−KAG∂ω∂xeλt−eλtФ+ρIλ2eλtФ=0,(31)−KAG∂∂x∂ω∂xeλt−eλtФ+ρAωλ2eλt−p∂2ω∂x2eλt+k1eλtω−Gp∂2ω∂x2∗eλt=fx,t.Reintroducing the dimensionless variables,(32)W=ωL,X=xL.Substituting in (30),(33)−β1∂∂XSX∂Ф∂X−β2∂W∂X−Ф=SXΩ2Ф.Substituting in (31),(34)−β3∂2W∂X2+β3∂Ф∂X+WΩ2−β4∂2W∂X2+β7W−β5∂2W∂X2=F,factoring for Ф,W,andWxx, equation (34) for a variable cross-section becomes(35)β3∂Ф∂X−β6∂2W∂X2+β7W=F+WΩ2,where(36)SX=EIEI0=II0=1+αxL3,A=A01+αxL,α=hL−h0hL,β1=AL2I0,β2=KA2GL4EI02,β3=KAGLEI0,β4=PL2EI0,β5=GpL2EI0,β6=β3+β4+β5,β7=K1L4EI0,Ω2=ρAL4λ2EI0,F=fx,teλtL3EI0.The following two boundary condition types are considered.A simply upheld end:(37)W=0,∂2Ф∂X2=0.A fixed end:(38)W=0,∂Ф∂X=0.Assuming fixed and simply upheld ends, end conditions can be shown, respectively, as(39)W1=0,∂Ф∂X1=0,(40)WN=0,∂2Ф∂XN2=0.


## 3. Differential Quadrature

In the differential quadrature method (DQM), the derivatives of a function at any point are approximated by a weighted linear summation of all the functional values at all other points as follows [22]:(41)fmx=∑j=1nCi,jfxji=1:n,where fmx is the m derivative of the function f(x) at the point x=xi, fxj is the functional value at point of x=xj, and Cmi,j is the weighting coefficient relating the m derivative at x=xi to the functional value at x=xj. To obtain the weighting coefficients, many polynomials with different base functions may be used. Lagrange interpolation formula is the most common one, where the functional value at a point x is approximated by all the functional values fxk (k = 1: n) as(42)fx≈∑k=1nLxx−xkL1xkfxk,where(43)Lx=∏k=1nx−xj,L1xk=∏k=1,k≠inxi−xk,i,k=1:n.Differentiating equation (42) and comparing the result with equation (41), the weighting coefficients of the first derivative are obtained as(44)C1i,j=L1xixi−xjL1xji≠j,i,j=1:n,(45)C1i,j=−∑j=1,j≠inC1i,ji=j,i,j=1:n.Applying the chain rule, the weighting coefficients for them derivative can be related to the weighting coefficients of the (m − 1)-derivative as(46)Cmi,j=∑k=1nC1i,jCm−1i,j,i,k=1:n.After choosing the sampling points which cover the solution domain, the weighting coefficients of the first derivative are calculated first, and then, the weighting coefficients of the higher derivative can be derived. As an example, if a total of 5 grid points are to be taken, the following sample of coefficients are obtained:(47)C11,1=1ξ1−ξ5+1ξ1−ξ4+1ξ1−ξ3+1ξ1−ξ2,C13,2=−ξ3−ξ1ξ3−ξ4ξ3−ξ5ξ2−ξ1ξ3−ξ2ξ2−ξ4ξ2−ξ5,where ξ is the grid point number and C1 is the weighting coefficient of the first derivative. Hence, to solve a boundary value problem, the governing differential equation is transformed at n-sampling points into a system on n-algebraic equations in n-unknown functional values of the solution function. The solution of such algebraic system considering the boundary condition leads to the functional values which satisfy both the governing differential equation and the boundary conditions.


### 3.1. Sampling of the Solution Points

As the DQM is a numerical method, its accuracy is affected by both the number and the distribution of the sampling points. One of the frequently used distributions is the normalized Gauss–Chebyshev–Lobatto distribution which is given as(48)xi=121−cosi−1n−1π,i=1:n.Substituting of the solution points’ coordinates into (44) and (45), the values of the weighting coefficients can be calculated, and then, the governing differential equation may be transformed to the algebraic system:(49)A1n∗mFjn∗1=B1n∗1,j=1:n.


### 3.2. Implementing the DQM

By using the definition of the DQM and substituting in (11), the governing equation becomes(50)∂2∂X2SX∑j=1NC2ij∗Wj+β1Wi+β4∑j=1NC2ij∗Wj=Fix,t+Ω2Wi,for i: 3, 4, … (N − 2).As for the Timoshenko model, substituting (46) in the governing equations (34) and (35) gives(51)β2Φi−β2∑j=1NC1ij∗Wj−SXβ1∑j=1NC2ij∗Φj−β1∂SX∂X∑j=1NC1ij∗Φj=SXΦiΩ2,for i: 1, 2, … (N), and(52)β3∑j=1NC1ij∗Φj−β6∑j=1NC2ij∗Wj+β7Wi=Fi+WiΩ2,for i: 3, 4, … (N − 2).


### 3.3. Shaping the Boundary Conditions

For an Euler beam, utilizing the DQM, the discretized conditions can be shown, respectively, as(53)W1=0,∑k=1NC1,kn0∗Wk=0,(54)WN=0,∑k=1NCN,kn1∗Wk=0,where n0 and n1 represent the boundary condition at the first end and last end, respectively. n0 and n1 may be taken as either 1 or 2 (1 for simply supported or “upheld” BC, 2 for fixed BC). By choosing n0 and n1, equations (53) and (54) can give the following sets of boundary conditions:n0  = 1, n1  = 1 ⟶ fixed from both sidesn0  = 1, n1  = 2 ⟶ fixed and upheldn0  = 2, n1  = 1 ⟶ upheld and fixedn0  = 2, n1  = 2 ⟶ upheld from both sidesThe first equations from (53) and (54) can be substituted easily into the governing equation (50). However, that is not the case for the second equations. Fortunately, they can be coupled to give two solutions W2 and WN−1 as(55)W2=1AXN∗∑k=3N−2AXK1∗Wk,(56)WN−1=1AXN∗∑k=3N−2AXKN∗Wk.According to (55) and (56), W2&WN−1 are expressed in terms of W3,W4,…,WN−2 and can be easily substituted into the governing equation system (50). It is noted that equations (53) and (54) provide four boundary equations. In total, we have N unknowns W1,…,WN. So, to close the system, the discretized governing equation has to be applied at (N −  4) grid points. This can be done by applying equation (50) at grid points X3,X4,…,XN−2 for W.It is noted that, after substitution, the system has (N − 4) equations and (N − 4) unknowns, which can be written in an eigenvalue matrix form:(57)AW=Ω2W,where AandB are the coefficients of W:(58)W=W3,W4,…,WN−2T.As for the Timoshenko beam, similar to the Euler beam, the discretized conditions can be shown, respectively, as(59)W1=0,∑k=1NC1,kn0∗Фk=0,(60)WN=0,∑k=1NCN,kn1∗Фk=0,where n0 and n1 represent the boundary condition at the first end and last end, respectively, as explained before in Section 3.3. According to (59) and (60), W2&WN−1 are expressed in terms of W3,W4,…,WN−2 and can be easily substituted into the governing equation system (51) and (52). It is to be noted that (59) and (60) provide four boundary equations. In total, we have 2∗N unknowns Ф1,…,ФN and W1,…,WN. So, to close the system, the discretized governing system has to be applied at (2∗N) − 4 grid points. This can be done by applying (51) and (52) at grid points X1,X2,…,XN for Ф and X3,X4,…,XN−2 for W. It is noted that, after substitution, the system has (2∗N) − 4 equations and (2∗N) − 4 unknowns, which can be written in an eigenvalue matrix form:(61)AФ+BW=Ω2Ф,(62)CФ+DW=Ω2W,where A,B,C,andD are the coefficients of ФandW from equations (61) and (62), respectively:(63)Ф=Ф1,Ф2,…,ФNT,W=W3,W4,…,WN−2T.


### 3.4. Formulation of the Moving Load

In Figure2, the load is assumed to move with a constant velocity “v” to traverse the beam length “L” through time “T.” A relation between the traversed distance “d” and the corresponding grid point “i” must be developed. A simple representation of how this relation is derived is as follows.Figure 2
Beam divided intoN grid points.Starting from the basic equation,(64)d=v∗t,and from Figure 2,(65)d=LN−1∗i−1.Therefore, substituting from (64) into (65) gives(66)v∗t=LN−1∗i−1.Factoring fori, results in the grid-time equation:(67)i=v∗t∗N−1L+1.


## 3.1. Sampling of the Solution Points

As the DQM is a numerical method, its accuracy is affected by both the number and the distribution of the sampling points. One of the frequently used distributions is the normalized Gauss–Chebyshev–Lobatto distribution which is given as(48)xi=121−cosi−1n−1π,i=1:n.Substituting of the solution points’ coordinates into (44) and (45), the values of the weighting coefficients can be calculated, and then, the governing differential equation may be transformed to the algebraic system:(49)A1n∗mFjn∗1=B1n∗1,j=1:n.


## 3.2. Implementing the DQM

By using the definition of the DQM and substituting in (11), the governing equation becomes(50)∂2∂X2SX∑j=1NC2ij∗Wj+β1Wi+β4∑j=1NC2ij∗Wj=Fix,t+Ω2Wi,for i: 3, 4, … (N − 2).As for the Timoshenko model, substituting (46) in the governing equations (34) and (35) gives(51)β2Φi−β2∑j=1NC1ij∗Wj−SXβ1∑j=1NC2ij∗Φj−β1∂SX∂X∑j=1NC1ij∗Φj=SXΦiΩ2,for i: 1, 2, … (N), and(52)β3∑j=1NC1ij∗Φj−β6∑j=1NC2ij∗Wj+β7Wi=Fi+WiΩ2,for i: 3, 4, … (N − 2).


## 3.3. Shaping the Boundary Conditions

For an Euler beam, utilizing the DQM, the discretized conditions can be shown, respectively, as(53)W1=0,∑k=1NC1,kn0∗Wk=0,(54)WN=0,∑k=1NCN,kn1∗Wk=0,where n0 and n1 represent the boundary condition at the first end and last end, respectively. n0 and n1 may be taken as either 1 or 2 (1 for simply supported or “upheld” BC, 2 for fixed BC). By choosing n0 and n1, equations (53) and (54) can give the following sets of boundary conditions:n0  = 1, n1  = 1 ⟶ fixed from both sidesn0  = 1, n1  = 2 ⟶ fixed and upheldn0  = 2, n1  = 1 ⟶ upheld and fixedn0  = 2, n1  = 2 ⟶ upheld from both sidesThe first equations from (53) and (54) can be substituted easily into the governing equation (50). However, that is not the case for the second equations. Fortunately, they can be coupled to give two solutions W2 and WN−1 as(55)W2=1AXN∗∑k=3N−2AXK1∗Wk,(56)WN−1=1AXN∗∑k=3N−2AXKN∗Wk.According to (55) and (56), W2&WN−1 are expressed in terms of W3,W4,…,WN−2 and can be easily substituted into the governing equation system (50). It is noted that equations (53) and (54) provide four boundary equations. In total, we have N unknowns W1,…,WN. So, to close the system, the discretized governing equation has to be applied at (N −  4) grid points. This can be done by applying equation (50) at grid points X3,X4,…,XN−2 for W.It is noted that, after substitution, the system has (N − 4) equations and (N − 4) unknowns, which can be written in an eigenvalue matrix form:(57)AW=Ω2W,where AandB are the coefficients of W:(58)W=W3,W4,…,WN−2T.As for the Timoshenko beam, similar to the Euler beam, the discretized conditions can be shown, respectively, as(59)W1=0,∑k=1NC1,kn0∗Фk=0,(60)WN=0,∑k=1NCN,kn1∗Фk=0,where n0 and n1 represent the boundary condition at the first end and last end, respectively, as explained before in Section 3.3. According to (59) and (60), W2&WN−1 are expressed in terms of W3,W4,…,WN−2 and can be easily substituted into the governing equation system (51) and (52). It is to be noted that (59) and (60) provide four boundary equations. In total, we have 2∗N unknowns Ф1,…,ФN and W1,…,WN. So, to close the system, the discretized governing system has to be applied at (2∗N) − 4 grid points. This can be done by applying (51) and (52) at grid points X1,X2,…,XN for Ф and X3,X4,…,XN−2 for W. It is noted that, after substitution, the system has (2∗N) − 4 equations and (2∗N) − 4 unknowns, which can be written in an eigenvalue matrix form:(61)AФ+BW=Ω2Ф,(62)CФ+DW=Ω2W,where A,B,C,andD are the coefficients of ФandW from equations (61) and (62), respectively:(63)Ф=Ф1,Ф2,…,ФNT,W=W3,W4,…,WN−2T.


## 3.4. Formulation of the Moving Load

In Figure2, the load is assumed to move with a constant velocity “v” to traverse the beam length “L” through time “T.” A relation between the traversed distance “d” and the corresponding grid point “i” must be developed. A simple representation of how this relation is derived is as follows.Figure 2
Beam divided intoN grid points.Starting from the basic equation,(64)d=v∗t,and from Figure 2,(65)d=LN−1∗i−1.Therefore, substituting from (64) into (65) gives(66)v∗t=LN−1∗i−1.Factoring fori, results in the grid-time equation:(67)i=v∗t∗N−1L+1.


## 4. Results and Verification

All outputs were produced using the open-source software WXMAXIMA. As a starting reference, the temporary physical and geometrical properties of the two models, foundation and the moving load, are listed in Table1, in which all are taken as unity and the resulting effect on the natural frequencies will be studied (due to the beam only, and this means that all external loads will be nullified here).Table 1
Temporary governing properties.
PropertyValueunitBeamYoung’s modulus E1PaShear modulus G• 1• PaMass densityρ1• kg m−3Beam height 0 h01• mBeam height 1 h11• mBeam width b1• mLength L1• mShear coefficient k1• -Α1• -Foundation• K1• 1• -• G• 1• -Moving load• F• 0• NAxial load• P• 0• NUsing Table1, the resulting first four nondimensional frequencies and natural frequencies of a simply supported Euler–Bernoulli beam and Timoshenko beam are free of any loads for n = 13 grid points which have been computed in.Table2 is a verification reference, and this result was also achieved by Taha [23] using the method of recursive differentiation. Another validation is the work of Kumar [24] in which an Euler beam is studied using Green’s function and the Rayleigh–Ritz method. These results have also been confirmed in the work of Ramzy et al. [9, 25].Table 2
Nondimensional and natural frequencies.
Euler–Bernoulli (N = 13)Timoshenko (N = 7)ΩnλnΩnλn19.8693.1419.7933.129239.4796.28334.7375.893388.8669.42661.8097.8614156.33912.5031034.36032.161The actual physical and geometrical properties of the Euler–Bernoulli beam, foundation, and moving load are listed in Tables3 and 4.Table 3
Actual governing properties.
ItemPropertyValueUnitBeamYoung’s modulus E2.1∗1011PaShear modulus G•8.1∗109• PaMass densityρ7860• kg m−3Beam height 0 h0Variable• mBeam height 1 h1Variable• mBeam width b0.02• mLength L1• mShear coefficient k5/6• -ΑVariable• -FoundationK1100• -G80• -Moving loadF30• NAxial loadPVariable• NTable 4
Current values of the natural frequencies.
Beam (1)Beam (2)E-BTimE-BTimλ11178.171169.10442.017443.3λ24712.694146.661767.441567.15λ310608.17378.253978.222781.11The problem introduced is compared with the model presented by Leszek [13]. First, the values of the natural frequency “λ” were compared. It is to be noted that to obtain these values one must first obtain the coefficient matrix and formulate the determinant then equating to zero. The resulting roots of the equation are the required nondimensional frequencies. The computation was done using n = 15 grid points for deflection computations and n = 13 for frequency computations. Four beam dimensions will be considered: beam (1) (deep, b=0.02m, and h=0.08m), beam (2) (slender, b=0.02m, and h=0.03m), and beam (3) (nonuniform, b=0.02m, and h0=0.03m,h1=0.05m). The following table shows the first three natural frequencies of the Euler–Bernoulli model (N = 13 grid points) and the Timoshenko model (N = 7 grid points). Note that the number of grid points taken for the Timoshenko beam (for calculating the natural frequency) is less than that of the Euler–Bernoulli beam, which is due to a limitation in the software where it would not solve the determinant matrix (for the frequency values) with a higher grid point value. However, the deflection analysis was accurately produced using N = 15 grid points since no determinant was required. The resulting natural frequencies to Leszek’s results for beam (1) and beam (2) were compared in Tables 5 and 6.Table 5
First 3 natural frequencies for beam (1).
λCurrent (EB.)Leszek [13] (EB.)Current (Tim.)Leszek [13] (Tim.)λ11178.171178.11169.101160.5λ24712.694712.64146.664447.4λ310608.110603.07378.259386.1Table 6
First 3 natural frequencies for beam (2).
λCurrent (EB.)Leszek [13] (EB.)Current (Tim.)Leszek [13] (Tim.)λ1442.017441.8443.30440.8λ21767.441767.21567.151752.2λ33978.223978.22781.123901.2The results from Table2 were also confirmed by the results achieved by Mohammad et al. [26] in which he studied the nonlocal fully intrinsic equations of Euler–Bernoulli beams with constitutive boundary conditions.Figure3 shows the deflection of the beam due to axial loading with a constant applied vertical load at midspan.Figure 3
S-S EB Beam (2) axial load,F = 30 N.Figures4–9 show the deflection due to a traversing load (30 N) at midspan at different instances, under different boundary conditions.Figure 4
S-S EB Beam (3) deflection,F = 30 N.Figure 5
S-S EB Beam (1) deflection,F = 30 N.Figure 6
C-S EB Beam (2) deflection,F = 30 N.Figure 7
C-C EB Beam (3) deflection,F = 30 N.Figure 8
C-S T Beam (1) deflection,F = 30 N.Figure 9
C-S T Beam (2) deflection,F = 30 N.Figures10–13 compare the deflection between both beam theories.Figure 10
C-S beam (3) deflection,F = 30 N.Figure 11
C-C beam (2) deflection,F = 30 N.Figure 12
C-C Beam (1) deflection,F = 30 N.Figure 13
C-S Beam (2) deflection,F = 30 N.Figures14–18 show the rotation of the cross section of both beams due to the traversing load (30 N). The outputs conform logically with the outputs from Figures 4–9.Figure 14
S-S T Beam (1) cross-section rotation.Figure 15
C-C T Beam (2) cross-section rotation.Figure 16
C-S T Beam (1) cross-section rotation.Figure 17
S-S T Beam (3) cross-section rotation.Figure 18
C-S T Beam (3) cross-section rotation.Comparing the results from Figure5 with an exact solution for the deflection of a simply supported beam due to a point load at midspan, it is calculated as w=FL3/48EI [27]. The error is calculated using the following formula:(68)Error=wapprox−wexactwexact∗100.Table 7
Deflection of S-S beam (1), (F = 30 N).
wapproxmwexactmError%3.4930∗10–63.4877∗10–60.15


## 5. Conclusion

In this work, the method of differential quadrature was applied on two beam models, that is, the Euler–Bernoulli and Timoshenko beams. The natural frequencies and deflection modes were derived from the basic constitutive equations and the numerical outputs of both problems were shown.The aim of this work is to compare the deflection of both beam theories subjected to a moving concentrated load as well as an axial load. The moving load was applied using a different approach other than the Dirac delta function; a relation between the grid point and time parameter was formed so as to allow manipulation of the load location by only changing the time parameter. The models discussed approved the viability of applying the moving load without analytical solutions while dealing with nodal equations. Summarizing the results based on the figures,(1)
From Figure11, it is noted that when dealing with slender beams the difference between the Timoshenko beam theory and Euler–Bernoulli theory tends to vanish. When dealing with deep beams, however, the difference is very noticeable and cannot be neglected.(2)
Euler beams produce smaller deflections than Timoshenko beams, which are explained by the presence of shear effect in Timoshenko beams, which allows for a greater deflection. Figures11 and 12 showcase this property in which Figure 11 (the slender beam) shows nearly identical deflection between both theories, while in Figure 12 (the deep beam), the difference is clearly visible.(3)
The deflection of both beam theories produces a maximum deflection near the center of the beam no matter the location of the moving vertical load. By increasing or decreasing the beam’s cross-sectional inertia, the maximum deflection will shift closer or further, respectively, from the center of the beam.(4)
Axial loading does not affect bending nearly as much as the traversing load. From Figure3, between (P=0N) and (P=20000N), an axial load “P” equivalent to 66,666% of the applied transverse load “F” produces only 25% greater deflection (from 6.62∗10–5 to 8.32∗10–5 meters).



---

*Source: 6616707-2021-04-20.xml*