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.

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

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

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 small

The foundation of the beam is taken as a Pasternak foundation with linear stiffness:where is the reaction of the foundation per unit length, is the beam vertical displacement, is the first order foundation parameter (elastic stiffness), and is the shear deformation coefficient.

Basic energy derivations for the Euler–Bernoulli beam model are

The beam bending energy is given as follows:

The system’s kinetic energy is given as follows:where is the beam density. The work done by the compressive force at the beam’s ends is

Applying Hamilton’s Principle,where L is the Lagrangian. Solving equation (6) results in the following fourth-order differential equation:

To neglect the effect of time, we use the method of separation of variables and assume that

Substituting then simplifying gives

Introducing the dimensionless variables,

Substituting the newly stated variables into equation (9) yields after simplification:where

The following two boundary condition types are considered.

A simply upheld end:

A fixed end:

Assuming a fixed and a simply upheld ends, end conditions can be shown, respectively:

2.2. The Timoshenko Beam Model

From the Timoshenko theory assumptions,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:

According to (16) and (17), the strain components are given as

From Hook’s Law,the beam’s bending energy is calculated as

Substituting (18) and (20) into (21) results in

As a result, the beam bending strain energy resting on a Pasternak foundation is expressed aswhere is the cross-section rotation, 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 is the shear modulus of the cross-section area. The system’s kinetic energy is given as follows:where is the density of the beam material. Applying Hamilton’s Principle,where is the Lagrangian:

Solving the system, the resulting coupled linear diff. equations are

To neglect the effect of time, assume that

Substituting the previous expansions into (27) and (28), respectively, yields

Reintroducing the dimensionless variables,

Substituting in (30),

Substituting in (31),factoring for , equation (34) for a variable cross-section becomeswhere

The following two boundary condition types are considered.

A simply upheld end:

A fixed end:

Assuming fixed and simply upheld ends, end conditions can be shown, respectively, as

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]:where is the m derivative of the function f(x) at the point , is the functional value at point of , and is the weighting coefficient relating the m derivative at to the functional value at . 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 (k = 1: n) aswhere

Differentiating equation (42) and comparing the result with equation (41), the weighting coefficients of the first derivative are obtained as

Applying the chain rule, the weighting coefficients for the m derivative can be related to the weighting coefficients of the (m − 1)-derivative as

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:where ξ is the grid point number and 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

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:

3.2. Implementing the DQM

By using the definition of the DQM and substituting in (11), the governing equation becomesfor i: 3, 4, … (N − 2).

As for the Timoshenko model, substituting (46) in the governing equations (34) and (35) givesfor i: 1, 2, … (N), andfor 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, aswhere and represent the boundary condition at the first end and last end, respectively. and may be taken as either 1 or 2 (1 for simply supported or “upheld” BC, 2 for fixed BC). By choosing and , equations (53) and (54) can give the following sets of boundary conditions:  = 1,  = 1 ⟶ fixed from both sides  = 1,  = 2 ⟶ fixed and upheld  = 2,  = 1 ⟶ upheld and fixed  = 2,  = 2 ⟶ upheld from both sides

The 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 and as

According to (55) and (56), are expressed in terms of 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 . 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 for .

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:where are the coefficients of :

As for the Timoshenko beam, similar to the Euler beam, the discretized conditions can be shown, respectively, aswhere 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), are expressed in terms of 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 2N unknowns and . So, to close the system, the discretized governing system has to be applied at (2N) − 4 grid points. This can be done by applying (51) and (52) at grid points for and for . It is noted that, after substitution, the system has (2N) − 4 equations and (2N) − 4 unknowns, which can be written in an eigenvalue matrix form:where are the coefficients of from equations (61) and (62), respectively:

3.4. Formulation of the Moving Load

In Figure 2, the load is assumed to move with a constant velocity “” 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.

Starting from the basic equation,and from Figure 2,

Therefore, substituting from (64) into (65) gives

Factoring for i, results in the grid-time equation:

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 Table 1, 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).

Using Table 1, 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.

Table 2 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].

The actual physical and geometrical properties of the Euler–Bernoulli beam, foundation, and moving load are listed in Tables 3 and 4.

The 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, , and ), beam (2) (slender, , and ), and beam (3) (nonuniform, , and ). 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.

The results from Table 2 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.

Figure 3 shows the deflection of the beam due to axial loading with a constant applied vertical load at midspan.

Figures 49 show the deflection due to a traversing load (30 N) at midspan at different instances, under different boundary conditions.

Figures 1013 compare the deflection between both beam theories.

Figures 1418 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 49.

Comparing the results from Figure 5 with an exact solution for the deflection of a simply supported beam due to a point load at midspan, it is calculated as [27]. The error is calculated using the following formula:

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 Figure 11, 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. Figures 11 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 Figure 3, between () and (), an axial load “” equivalent to 66,666% of the applied transverse load “F” produces only 25% greater deflection (from 6.6210–5 to 8.3210–5 meters).

Data Availability

The data used to support the findings of this research are included within the article.

Conflicts of Interest

The authors declare that they have no conflicts of interest regarding the publication of this research.