Timoshenko Beam Theory: First-Order Analysis, Second-Order Analysis, Stability, and Vibration Analysis Using the Finite Difference Method

: This paper presents an approach to the Timoshenko beam theory (TBT) using the finite difference method (FDM). The Timoshenko beam theory covers cases associated with small deflections based on shear deformation and rotary inertia considerations. The FDM is an approximate method for solving problems described with differential equations. It does not involve solving differential equations; equations are formulated with values at selected points of the structure. In addition, the boundary conditions and not the governing equations are applied at the beam’s ends. The model developed in this paper consisted of formulating differential equations with finite differences and introducing additional points at the beam’s ends and at positions of discontinuity (concentrated loads or moments, supports, hinges, springs, brutal change of stiffness, spring − mass system, etc.). The introduction of additional points allowed us to apply the governing equations at the beam’s ends. Moreover, grid points with variable spacing were considered, the grid being uniform within beam segments. First-order, second-order, and vibration analyses of structures were conducted with this model. Furthermore, tapered beams were analyzed (element stiffness matrix, second-order analysis, vibration analysis). Finally, a direct time integration method (DTIM) was presented; the FDM-based DTIM enabled the analysis of forced vibration of structures, with damping taken into account. The results obtained in this paper showed good agreement with those of other studies, and the accuracy was increased through a grid refinement. Especially in the first-order analysis of uniform beams, the results were exact for uniformly distributed and concentrated loads regardless of the grid. and those obtained by Hibbitt et al. [17] with the finite element method using ABAQUS software. In this study the analysis was conducted with n = 9, 17, 25, 33, and 41 grid points for different values of the taper ratio ( 1-h r /h l ) and support conditions: h l and h r are the heights at the left and the right beam’s end, respectively. The reference values Ar and Ir are taken at the left beam’s end. Detailed results are listed in the Supplementary Material “Vibration analysis of tapered Timoshenko beams.” The results of this study are compared to those of Soltani [16]

In the classical analysis using the FDM, points outside the beam are not considered. The boundary conditions and not the governing equations are applied at the beam's ends. Consequently, the non-application of the governing equations at the beam's ends has led to inaccurate results, making the FDM less interesting in comparison to other numerical methods such as the finite element method. This paper presented a model based on the FDM; it consisted of formulating the differential equations with finite differences and introducing additional points at the beam's ends and at positions of discontinuity (concentrated loads or moments, supports, hinges, springs, change of grid spacing, brutal change of stiffness, and spring−mass system, etc.). The introduction of additional points allowed us to apply the governing equations at the beam's ends. First-order, second-order, and vibration analyses of structures were conducted using the model. Finally, a direct time integration method (DTIM) was presented; the FDM-based DTIM enabled the analysis of forced vibration of structures, the damping being considered.

Statics
The sign convention adopted for the loads, bending moments, shear forces, and displacements is illustrated in Figure 1. Sign convention for loads, bending moments, shear forces, and displacements Specifically, M(x) is the bending moment in the section, V(x) is the shear force, w(x) is the deflection, and q(x) is the distributed load in the positive downward direction.
In first-order analysis the equations of static equilibrium on an infinitesimal element are as follows: where k(x) is the stiffness of the elastic Winkler foundation. Substituting Equation (1b) into Equation (1a) yields According to Timoshenko beam theory, the bending moment and shear force are related to the deflection and rotation (positive in clockwise) of the cross section (x), as follows: where E is the elastic modulus, I is the second moment of area,  is the shear correction factor, G is the shear modulus, and A is the cross-sectional area.
In the case of a uniform beam, substituting Equations (3) and (4) into Equations (1a) and (1b) yields (5a) (5b) In the case of a tapered beam, substituting Equations (3) and (4) into Equations (1a) and (1b) yields (6a) (6b) Fogang [6] presented the following relationships for a uniform beam and a tapered beam, respectively: Differentiating Equation (6c) twice with respect to x and combining the result with Equation (2) yields the following widely known formula for a uniform beam without Winkler foundation: In the presence of an elastic Winkler foundation, Equation (6e) becomes For a uniform beam, the bending moment, the shear force, and the rotation of the cross section are derived using Equations (6c) and (2), Equation (1b), and Equation (4), respectively, as follows: In summary, a W− finite difference approximation (FDA) using Equations (5a) to (6b), an M−W FDA using Equations (2) and (6c-d), and a W FDA using Equations (6e-i) can be considered.  Thus, curves w(x) and (x) can be described with the deflections values at equidistant grid points:

Fundamentals of FDM
(7a) The shape functions f j (x) (j = i-1, i, i+1) can be expressed using Lagrange polynomials: Thus, a three-point stencil is used to write finite difference approximations to derivatives at grid points. The derivatives (S(x) representing w(x) or (x)) at i are expressed with deflection values at points i-1, i, and i+1.

M−W FDA in first-order analysis of a uniform beam
Equations (2) and (6c) are the governing equations. Applying Equations (8a-b) and (9a-e) in Equations (2)

M−W FDA in first-order analysis of a tapered beam
Equations (2) and (6d) are the governing equations. Substituting Equations (8a-b), (9d), and (14a-d) into Equations (2) and (6d) yields Equation (13d) and the following equation: Thus, Equations (13d) and (16) are the governing equations. The shear force and the rotation of the cross section are calculated using Equations (13f) and (13g), respectively. However,  Vk is replaced by  Vi in Equation (13g).

FDA of q(x) and k(x)
Fogang [1] presented formulas to determine FDA of distributed loads and stiffness of an elastic Winkler foundation. The FDA q i for position i being the left beam's end, an interior point on the beam, or the right beam's end was as follows: The application of Equations (17a) to (18) shows that for a linearly distributed load, q i = q(xi).
At any point i, the stiffness of the elastic Winkler foundation, ki, is calculated similarly to Equations (17a) to (18).

Analysis at positions of discontinuity
Positions of discontinuity are positions of application of concentrated external loads (force or moment), supports, hinges, springs, abrupt change in cross section, positions where EI(x) and GA(x) are not differentiable, and change in grid spacing.

Uniform beam within segments
Concentrated loads (force P and moment M * ) are applied at point i, as shown in Figure 3. The beam has a uniform cross section within segments; at point i, an abrupt change in cross section and a change in grid spacing are assumed.       Figure 5a,b below shows the additional points (fictive points ia, ib, ic, id) introduced in the opening. The unknown at any point is the deflection. In the equations above, il, Mil, and Vil are formulated by adopting for i, i+1, and i+2 the values of il, ia, and ib, respectively. Similarly, ir, Mir, and Vir are formulated by adopting for i, i-1, and i-2 the values of ir, id, and ic, respectively.

M−W FDA of a uniform beam
The additional points of Figure 4a

Tapered beam within segments
As described in Section 2.1.3.1, an opening of the beam is realized at point i and additional points (fictive points ia, id) are introduced in the opening (Figure 4a,b).

W− FDA of a tapered beam
The governing equations (Equations (15a-b)) are applied at any point of the beam: … i-1, il, ir, i+1 … The continuity equations can be expressed through an adjustment of Equations (20a-d), as follows: (22)

M−W FDA of a tapered beam
The governing equations (Equations (13d) and (16)) are applied at any point of the beam: … i-1, il, ir, i+1 … The continuity equations can be expressed using Equations (21a-d), whereby the shear force and the rotation of the cross section are calculated using Equations (13f) and (13g), respectively. However,  Vk is replaced by  Vi in Equation (13g).

Mixed FDA of a tapered beam
Similar to the uniform beam, different approximations (W−, M−W) can be considered on either side of the point of discontinuity. The continuity equations are then formulated with the corresponding formulas.

Non-uniform grid
The grid may be such that every node has a non-constant distance from another ( Figure 6).

Figure 6 Beam with a non-uniform grid
In this case, the Lagrange interpolation polynomial (Equation (7b)) was used for FDM formulations. The resulting equations were complicated, so, the non-uniform grid was not further analyzed. In fact, it should not be analyzed as a discontinuity position.

44 element stiffness matrix
The sign convention for bending moments, shear forces, displacements, and rotations of the cross section adopted to determine the element stiffness matrix in local coordinates is illustrated in Figure 7. Figure 7 Sign convention for moments, shear forces, displacements, and rotations for the stiffness matrix Let us define the following vectors: The 44 element stiffness matrix in local coordinates of the tapered beam is denoted by K 44 .
The vectors defined are related together with the element stiffness matrix K 44 , as follows: Let us divide the beam in n parts of equal length h (l = nh k ), as shown in Figure 8. ; ; ; Considering the sign convention adopted for displacements and rotations of cross sections in general ( Figure 1) and in the element stiffness matrix (Figure 7), the following geometric compatibility boundary conditions can be set: The number of equations is 2(n+1) + 4 + 4 = 2n + 10. The number of unknowns is 2(n+3) + 4 = 2n + 10, especially 2(n+3) unknowns (W; ) at points on the beam and additional points at the beam's ends, and four efforts at the beam's ends (Vi; Mi; Vk; Mk). Let us define the following vector (27) The combination of Equations (15a-b) applied at any point on the grid, Equations (25a−d), and Equations (26a−d) can be expressed with matrix notation as follows, the geometric compatibility boundary conditions (Equations (26a−d)) being at the bottom: The matrix T has 2n+10 rows and 2n+10 columns. The zero vector above has 2n+6 rows.
The matrix Taa has 2n+6 rows and 2n+6 columns, the matrix Tab has 2n+6 rows and 4 columns, the matrix Tba has 4 rows and 2n+6 columns, and the matrix Tbb has 4 rows and 4 columns.
The combination of Equations (24), (28), and (29) yields the element stiffness matrix of the beam. (30a) A general matrix formulation of K 44 is as follows: In Equation (30b), 0 is a zero matrix with 4 rows and 2n+6 columns, and I is the 4  4 identity matrix.

M−W FDA for element stiffness matrix of a tapered beam
Equations (13d) and (16) with qi = 0 and ki = 0 are applied at any point on the grid (nodes 1, 2, …n+1 of Figure 8).

33 element stiffness matrix
Assuming the presence of a hinge at the right end, the sign convention for bending moments, shear forces, displacements, and rotations of the cross section is illustrated in Figure 9. The 33 element stiffness matrix in local coordinates of the tapered beam is denoted by K 33 .
The vectors of Equations (23a-b) and (24) become The matrix K 33 can be formulated with the values of the matrix K 44 (see Equations (30a-b)). (33) The matrix K 44 has 4 rows and 4 columns, the matrix Kaa has 3 rows and 3 columns, the matrix Kab has 3 rows and 1 column, the matrix Kba has 1 row and 3 columns, and the matrix Kbb has 1 row and 1 column (a single value).
The combination of Equation (24) with the presence of a hinge at position k (Mk = 0) and Equation (32c) yields the matrix K 33 , as follows: (34)

Second-order analysis
The equation of static equilibrium can be expressed as follows: The axial force (positive in tension) is denoted by N(x) and the transverse force by T(x). Let us consider an external distributed axial load n(x) positive along the + x axis The transverse force T(x) is related to the shear force V(x), as follows: ; ;

Second-order analysis of a uniform beam within segments
The grid spacing h k , the reference flexural rigidity EI r , the reference shear stiffness GA r , and the parameters  lk ,  Mk ,  Vk , and  r are defined similarly to previous sections.

W− FDA in second-order analysis of a uniform beam
Substituting Equations (4) (3) and (4) The bending moment is calculated using Equation (11a). The analysis at positions of discontinuity is conducted similarly to the first-order analysis; however, the shear force is replaced by the transverse force.

W FDA in second-order analysis of a uniform beam
It is assumed here that the axial force and stiffness of the Winkler foundation are constant along the beam. Substituting Equations (2) (41)  (43b)

M−W FDA in second-order analysis of a uniform beam
Substituting Equations (2)

M−W FDA in second-order analysis of a tapered beam
The governing equations (Equations (16) and (43f)) are applied at any point on the grid. The rotations of the cross sections and the transverse forces are calculated using Equations (13g) and (43g), respectively. However,  Vk is replaced by  Vi in Equation (13g).

Second-order element stiffness matrix of a uniform beam
The beam is divided in n parts of equal length h k , as shown in Figure 10.

Second-order element stiffness matrix of a tapered beam
The M−W FDM approximation is applied here. The W− FDM approximation can also be applied with appropriate formulas developed previously. The sign convention for bending moments, transverse forces, displacements, and rotations of the cross sections adopted to determine the element stiffness matrix in local coordinates is illustrated in Figure 7, the shear forces V i and V k being replaced by the transverse forces T i and T k .

Free vibration analysis
It is focused here on the natural frequencies of the beam. A second-order analysis is conducted; and the first-order analysis can easily be deduced. The equations of dynamic equilibrium on an infinitesimal beam element are as follows: , .
where  is the beam's mass per unit volume, A(x) is the cross-sectional area, N(x) is the axial force (positive in tension), and k(x) is the stiffness of the elastic Winkler foundation. A harmonic vibration being assumed, T * (x,t), M * (x,t), w * (x,t), and  * (x,t) can be expressed as follows (S * (x,t) representing T * (x,t), M * (x,t), w * (x,t), and  * (x,t)): where  is the circular frequency of the beam. Substituting Equation (47)

W− FDA in free vibration analysis of a uniform beam
A reference flexural stiffness EI r , a reference shear stiffness GA r , parameters  Mk ,  Vk , and  r were defined in Equations (9a-c). A reference cross-sectional area A r and a reference length l r are defined and related to the crosssectional area A k and the grid spacing h k in the segment, as follows: The parameter  r (Equation (9c)) is defined with l r instead of l. Substituting Equations (3) and (4)

Effect of a concentrated mass, or a spring
The dynamic behavior of a beam carrying a concentrated mass or having a spring was analyzed, as shown in Figure 11. Effect of a spring−mass system: The dynamic behavior of a beam carrying a spring−mass system was analyzed, as represented in Figure 12. The deflection of the mass is denoted by w iM .
The dynamic behavior of a beam carrying a concentrated mass, a spring, or a spring−mass system was analyzed similarly to the previous section (Equations (55a)-(57b)).

Free vibration analysis of tapered beams
The W− FDM approximation was considered for the vibration analysis of tapered beams.
Equations (61a-b) are applied at any point on the grid. The bending moments and transverse forces are determined using Equations (11a) and (40), respectively,  Mk and  Vk being replaced by  Mi and  Vi .

Effect of a concentrated mass, a spring, or a spring−mass system
The dynamic behavior of a beam carrying a concentrated mass, a spring, or a spring−mass system is analyzed similarly to the previous section (Equations (55a)-(57b). The transverse forces Til and Tir are calculated using Equation (40),  Vk being replaced by  Vi .

Direct time integration method
The direct time integration method (DTIM) developed here describes the dynamic response of a beam as a multi-degreeof-freedom system. Viscosity  and external loading p(x,t) are considered.

DTIM for uniform beams within segments
At initial time t = 0, a three-point forward difference approximation is applied (Equation (18a)): At final time t = T, a three-point backward difference approximation is applied (Equation (18c)): (63c) The governing equations (Equations (62a-b)) can be formulated with the FDM for x = i at time t. The FDM formulations of these equations are applied at any point of the beam at any time t using a five-point stencil. Additional points are introduced to satisfy the boundary and continuity conditions. The boundary conditions are satisfied using a three-point stencil. Thus, beam deflection w * (x,t) and rotation  * (x,t) can be determined with the Cartesian model represented in Figure 13. The bending moment M * (x,t), shear force V * (x,t), and transverse force T * (x,t) are calculated using Equations

DTIM for tapered beams
A similar analysis can be conducted. Thus, Equations (62a-b) become The derivatives with respect to x are formulated using Equations (8a-b), while those with respect to t (time increment is t) are formulated considering a three-point stencil with Equations (63a-c). The FDM formulations of Equations (64a-b) are applied at any point on the beam and at any time t using a five-point stencil. Additional points are introduced to satisfy the boundary and continuity conditions. The boundary conditions are satisfied using a three-point stencil. Thus, beam deflection w * (x,t) and rotation  * (x,t) can be determined with the Cartesian model represented in Figure 13. The bending moment M * (x,t), shear force V * (x,t), and transverse force T * (x,t) are calculated using Equations (11a-b) and (37), respectively,  Vk being replaced by  Vi in Equation (11b) and  Mk by  Mi in Equation (11a).
With this model, the assumptions made previously can be verified, namely the separation of variables and the harmonic vibration (Equation (47)).
The bending shear factor  = EI/GAl² = 0.025 The analysis was conducted with W− FDA, M−W FDA, and W FDA. Details of the analysis and results are presented in Appendix A and in the Supplementary Material "Fixed−pinned beam subjected to a uniformly distributed load." Table 1 lists the exact results obtained with classical beam theory (CBT) and those obtained with W− FDA and M−W FDA. Table 2 lists the results obtained with CBT and those obtained with W FDA.  The results obtained with W FDA are exact for a uniformly distributed load regardless of the grid. An explanation is the fact that the exact solution for the deflection curve is a fourth-order polynomial, which corresponds to the FDM polynomial hypothesis.

Beam subjected to a concentrated load
We analyzed a uniform fixed−pinned beam subjected to a concentrated load, as represented in Figure 15.  Table 3 shows the exact results (bending moments) obtained with CBT and those obtained in this study (W− FDA, M−W FDA, and W FDA). for the deflection curve is a third-order polynomial, which is exactly described with the fourth-order polynomial FDM approximation.

Tapered pinned−fixed beam subjected to a uniformly distributed load
A tapered pinned−fixed beam subjected to a uniformly distributed load, as shown in Figure 16, was analyzed.

Beam subjected to a uniformly distributed load and a compressive force
A uniform fixed−pinned beam subjected to a uniformly distributed load and a compressive force, as shown in Figure 17, was analyzed.  Table 5 lists the results obtained using [6] and those in this study (W− FDA, M−W FDA, and W FDA). The results of all of the approximations converge towards the results of Fogang [6], the accuracy increasing with an increasing number of grid points. W FDA yields here the best results.

Buckling load of a fixed−pinned beam
The buckling load of a uniform fixed−free beam, as represented in Figure 18, was determined. Figure 18 Buckling load of a uniform fixed−free beam The analysis was conducted with n = 9 and 17 grid points. The buckling load Ncr is defined as follows: Hu et al. [5] presented the following closed-form expression of the buckling load of a uniform fixed−free beam: The combination of Equations (9c), (66), and (67a) yields the buckling factor  as follows:  Table 6 shows the results obtained using the formula by Hu et al. [5] and those obtained in this study (W− FDA, M−W FDA, and W FDA) for different values of the bending shear factor  (Equation (9c)). The results of all of the approximations converge towards those using Hu et al. [5], the accuracy increasing with an increasing number of grid points. M−W FDA yields here the best results.

Second-order element stiffness matrix of a uniform beam
The element stiffness matrix of a uniform beam with following characteristics was calculated.  TB  TB  TB  TB   TB  TB  TB  Tbl  TB  TB   TB   T  Q  T The aforementioned characteristics become P = 1.5  EI/L²,  = 1-P/(ksGA) = 1-1.5  0.05 = 0.925, Details of the analysis and results are presented in Appendix D and in the Supplementary Material "Second-order element stiffness matrix of a uniform beam." Table 7 lists the results obtained using Hu et al. [5] and those obtained in this study. The results of both approximations converge towards those using Hu et al. [5], the accuracy increasing with an increasing number of grid points. W FDA yields here the best results.

Free vibration analysis of a fixed−free beam
The natural frequencies of a fixed−free beam were determined, depending on the bending shear factor  and the coefficient of rotary inertia kRI. The analysis was conducted with n = 9, 17, and 25 grid points. Details of the analysis and results are listed in Appendix E and in the Supplementary Material "Vibration analysis of a uniform fixed−free beam." The coefficients  of vibration frequency are defined in Equation (52b). The results obtained in this paper (W− FDA and M−W FDA) are compared to those obtained using Kruszewski [15], and are listed in Table 8.

TIMOSHENKO BEAM THEORY USING THE FINITE DIFFERENCE METHOD
The results of both approximations converge towards those using Kruszewski [15], the accuracy increasing with an increasing number of grid points. M−W FDA yields better results than W− FDA for a given grid.

Free vibration analysis of beams resting on Winkler foundation and subjected to a compression force
The dynamic response of beams resting on an elastic Winkler foundation and subjected to an axial load was determined.
A pinned−pinned beam and a fixed−pinned beam were analyzed.
Ghannadiasl [14] analytically The definition of the stiffness of the Winkler foundation in Ghannadiasl [14] has an error: in the denominator, it should be L 4 instead of L 2 .
In the present study, the analysis was conducted with n = 9, 17, and 33 grid points. Detailed results are listed in the Supplementary Materials "Vibration analysis of a pinned−pinned beam with an axial load" and "Vibration analysis of a fixed−pinned beam with an axial load." Table 9 and Table 10 list the results of Ghannadiasl [14] and those obtained in this study (W− FDA and M−W FDA).  The results of both approximations converge towards those of Ghannadiasl [14], the accuracy increasing with an increasing number of grid points. M−W FDA yields better results than W− FDA for a given grid.

Free vibration analysis of tapered Timoshenko beams
The natural frequencies (coefficients ) of tapered Timoshenko beams were determined. Pinned−pinned, fixed−free, and fixed−fixed beams were considered.
The beams have the following characteristics: Poisson's ratio  = 0.30, Timoshenko shear coefficient  = 5/6, and coefficient of rotary inertia k RI = 0.01. The bending shear factor is therefore calculated as follows: Soltani [16] presented results obtained with the power series method (PSM) and those obtained by Hibbitt et al. [17] with the finite element method using ABAQUS software. In this study the analysis was conducted with n = 9, 17, 25, study are compared to those of Soltani [16] and Hibbitt et al. [17] in Table 11.  The results of both approximations converge towards those of Soltani [16] and Hibbitt et al. [17], whereby the accuracy increases with an increasing number of grid points.

Conclusions
The FDM-based model developed in this paper enabled, with relative easiness, first-order analysis, second-order analysis, and vibration analysis of Timoshenko beams. The results showed that the calculations conducted as described in this paper were accurate, and the accuracy increased with an increasing number of grid points. Especially in first-order analysis of uniform beams, the results were exact for uniformly distributed and concentrated loads regardless of the discretization. First-and second-order element stiffness matrices (tensile or compressive axial force) in local coordinates were determined. In addition, tapered beams were analyzed.
The following aspects were not addressed in this study but could be analyzed with the model in future research: