EULER − BERNOULLI BEAM THEORY: FIRST-ORDER ANALYSIS, SECOND-ORDER ANALYSIS, STABILITY, AND VIBRATION ANALYSIS USING THE FINITE DIFFERENCE METHOD

: This paper presents an approach to the Euler − Bernoulli beam theory (EBBT) using the finite difference method (FDM). The EBBT covers the case of small deflections, and shear deformations are not considered. The FDM is an approximate method for solving problems described with differential equations (or partial differential equations). The FDM does not involve solving differential equations; equations are formulated with values at selected points of the structure. The model developed in this paper consists of formulating partial differential equations with finite differences and introducing new points (additional points or imaginary points) at boundaries and positions of discontinuity (concentrated loads or moments, supports, hinges, springs, brutal change of stiffness, etc.). The introduction of additional points permits us to satisfy boundary conditions and continuity conditions. First-order analysis, second-order analysis, and vibration analysis of structures were conducted with this model. Efforts, displacements, stiffness matrices, buckling loads, and vibration frequencies were determined. Tapered beams were analyzed (e.g., element stiffness matrix, second-order analysis). 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. The efforts and displacements could be determined at any time. accuracy increasing


Introduction
The Euler−Bernoulli beam has been widely analyzed in the relevant literature. Several methods have been developed, e.g., the force method, the slope deflection method, and the direct stiffness method, etc. Anley et al. [1] considered a numerical difference approximation for solving two-dimensional Riesz space fractional convection-diffusion problem with source term over a finite domain. Kindelan et al. [2] presented a method to obtain optimal finite difference formulas which maximize their frequency range of validity. Both conventional and staggered equispaced stencils for first and second derivatives were considered. Torabi et al. [3] presented an exact closed-form solution for free vibration 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.
According to the Euler−Bernoulli beam theory (EBBT), the governing equation for the case of a uniform beam loaded with q(x) is as follows: , (1) where EI is the flexural rigidity and k(x) is the stiffness of the elastic Winkler foundation.
The bending moment, shear force, and slope (x) are related to the deflection as follows:  The governing equation (Equation (1)) has a fourth order derivative, the deflection curve is consequently approximated around the point of interest i as a fourth degree polynomial.

FDM Formulation of equations, efforts, and deformations
Thus, the deflection curve can be described with the values of deflections at equidistant grid points: (3a) The shape functions f j (x) (j = i-2; i-1; i; i+1; i+2) can be expressed using the Lagrange polynomials: Thus, a five-point stencil is used to write finite difference approximations to derivatives at grid points. Therefore, the derivatives at i are expressed with values of deflection at points i-2; i-1; i; i+1; i+2. We set W(x) = EIw(x) (5) Substituting Equations (4a) and (5) into Equation (1) yields (6) At point i, the bending moment, shear force, and slope are formulated with finite differences using Equations (2a−c), (4b−d), and (5).
Let us determine here the FDM value q i (Equation (6)) in the case of a varying distributed load q(x). The distributed load q(x) is related to the shear force V(x) as follows: Considering here a three-point stencil, the following FDM formulations of the first derivative can be made. The position i is considered as the left beam's end, an interior point on the beam, or the right beam's end, respectively: The balance of vertical forces applied to a free body diagrams yields the following: The combination of Equations (8a−f) yields the FDM value q i for the position i being the left beam's end, an interior point on the beam, or the right beam's end.
The application of Equations (8g−8i) shows that in the case of a linearly distributed load, q i is equal q(xi).
At point i, the stiffness of the elastic Winkler foundation ki is calculated similarly to Equations (8g−8i).

Analysis at positions of discontinuity
Positions of discontinuity are positions of application of concentrated external loads (force or moment), supports, hinges, springs, abrupt change of cross section, and change of grid spacing.

Concentrated load or moment
Let us analyze here the case of concentrated loads (force P and moment M * ) applied at point i, as represented in Figure 3. The beam has a uniform cross section. The model developed in this paper consists of realizing an opening of the beam at point i and introducing additional points (fictive points ia, ib, ic, and id) in the opening, as represented in Figure 4a,b.
Opening of the beam and introduction of additional points on the left side (4a) and right side (4b).
Thus, the governing equations at positions il and ir yield: The continuity equations express the continuity of the deflection and slope (Equation (7c)), and the equilibrium of the bending moment (Equation (7a)) and shear force (Equation (7b)): An adjustment of the continuity equations is made in case of a hinge (no continuity of the slope, Mil = Mir = 0), a support (Wil = Wir = 0, no equation (10d)), or a spring.
At the beam's ends, additional points are introduced (as shown in Figure 4a,b) and so governing equations are applied at the beam's ends, as well as the boundary conditions.

Abrupt change of cross section
A reference flexural rigidity EI r is introduced. Therefore, we set The governing equation and FDM formulation of the moment and shear force become The continuity equations are formulated in Equations (10a−d) with consideration of Equations (12b,c).

Change of grid spacing
The discretization of the beam may lead to uniform-grid segments, but the grid spacings being different from one segment to another, as represented in Figure 5.

Non-uniform grid
The grid may be such that every node has a non-constant distance from another, as represented in Figure 6.
Here, the Lagrange interpolation polynomial (Equation (3b)) is used for FDM formulation. The resulting equations are complicated, and consequently the non-uniform grid is not further analyzed in this paper. In fact, this case should not be analyzed as a discontinuity position.

First-order analysis of a tapered beam
Let us analyze here the case of a tapered beam, as shown in Figure 7.

Analysis at positions of discontinuity of a tapered beam
Let us determine the continuity equation for the case of a concentrated load or moment.  The additional points are ia and id, and the unknowns are wia, Mia, wid, and Mid.
The continuity equations express the continuity of the deflection and slope (Equation (15d)), and the equilibrium of the bending moment and shear force (Equation (15c)) as follows: At a point i with a hinge the slopes are not more equal, and the moments Mil and Mir are zero.

44 element stiffness matrix
The sign conventions for bending moments, shear forces, displacements, and slope adopted for use in determining the element stiffness matrix in local coordinates is illustrated in Figure 9.
. Sign conventions for moments, shear forces, displacements, and slopes for stiffness matrix.
Let us define following vectors: The 44 element stiffness matrix in local coordinates of the tapered beam is denoted by K 44 .
The vectors above 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) as shown in Figure 10. Considering the sign conventions adopted for bending moments and shear forces in general (see Figure 1) and for bending moments and shear forces in the element stiffness matrix (see Figure 9), we can set following static compatibility boundary conditions in combination with Equations (2b) and (14b): ; ; ; Considering the sign conventions adopted for the displacements and slope in general (see Figure 1) and for displacements and slope in the element stiffness matrix (see Figure 9), we can set following geometric compatibility boundary conditions in combination with Equations (2c) and (14b): The number of equations is 2(n+1) + 4 +4 = 2n + 10. The number of unknowns is 2(n+3) + 4 = 2n + 10, 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 (18), (22), and (23) yields the element stiffness matrix of the beam. (24a) A general matrix formulation of K 44 is as follows: (24b) In Equation (24b), 0 is a zero matrix with four rows and 2n+6 columns, I is the 4  4 identity matrix. 1 1

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 slope is illustrated in Figure 11. Figure 11. Sign conventions for moments, shear forces, displacements, and slope for stiffness matrix.
The 33 element stiffness matrix in local coordinates of the tapered beam is denoted by K 33 .
The vectors of Equations (17a), (17b), and (18) become The matrix K 33 can be formulated with the values of the matrix K 44 (see Equations (24a−b)).
(26) 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 (18) with the presence of a hinge at position k (Mk = 0), and Equation (25c) yields the matrix K 33 as follows: (27)

First-order element stiffness matrix of a uniform beam
The beam is divided in n parts of equal length h (l =nh) as shown in Figure 12. ; ; 12a) with qi = 0 is applied at any point on the grid (positions 1; 2; …n+1 of Figure 12).
Thus, the number of equations is n+9. The number of unknowns is n+5 + 4 = n + 9, especially n+5 unknowns W (at points on the beam and additional points at beam's ends) and 4 efforts at beam's ends (Vi; Mi; Vk; Mk).
The vector becomes, The use of Equations (22) to (24b) yields the element stiffness matrix of the uniform beam.

Second-order analysis
The equation of static equilibrium can be expressed as follows: The axial force (positive in tension) is denoted by N, and the stiffness of the elastic Winkler foundation by k. Let us also consider an external distributed axial load n(x) positive along the + x axis The combination of Equations (30) and (31) Equation (34) is applied at any point on the grid. At point i, the external distributed axial load ni is calculated similarly to Equations (8g−i).
The analysis at positions of discontinuity is conducted similarly to that of the first-order analysis; however the shear force is replaced by the transverse force.

Second-order analysis of a tapered beam
The FDM formulation of Equation (32) of static equilibrium is as follows: Equations (37)  The slope is calculated similarly to Equation (15d).
The analysis at positions of discontinuity is conducted similarly to that of the first-order analysis; however the shear force is replaced by the transverse force.

44 element stiffness matrix
The sign conventions for bending moments, transversal forces, displacements, and slopes adopted for use in determining the element stiffness matrix in local coordinates is illustrated in Figure 13. The FDM discretization is the same as Figure 10. Equations (17a) becomes,  (37) and (15b) are applied at any point on the grid (the distributed load qi being zero).
The static compatibility boundary conditions are expressed similarly to Equations (19a−d); however in Equations (19a) and (19c), the shear forces are replaced by the transverse forces (Equation (38)). The analysis continues similar to the first-order element stiffness matrix (Equations (20a−24)).

Free vibration analysis
Our focus here is to determine the eigenfrequencies of the beams.

Beam with constant stiffness EI within segments
The governing equation is as follows: where  is the beam's mass per unit volume, A is the cross-sectional area, N(x) is the axial force (positive in tension), n(x) is the external distributed axial load positive along + x axis, and k is the stiffness of the elastic Winkler foundation. A harmonic vibration being assumed, w * (x,t) can be expressed as follows: Here,  is the circular frequency of the beam. Substituting Equation (41)

Effect of a concentrated mass, or a spring
We analyzed the dynamic behavior of a beam carrying a concentrated mass or having a spring, as represented in Figure 14. Effect of a spring−mass system: We analyzed the dynamic behavior of a beam carrying a spring−mass system, as represented in Figure 15. The deflection of the mass is denoted by w iM . Figure 15. Vibration of a beam carrying a spring−mass system.   The derivatives with respect to x are formulated with Equations (14a,b); those with respect to t are formulated with Equations (50a) to (50c).
The FDM formulation of Equations (51) and (46b) are applied at any point on the beam and at any time t using a fivepoint stencil and a three-point stencil, respectively. Additional points are introduced to satisfy the boundary and continuity conditions. The boundary conditions are satisfied using a three-point stencil. Thus, the bending moment M * (x,t) and beam deflection w * (x,t) can be determined with the Cartesian model represented in Figure 17. The transverse force T * (x,t) is calculated using Equation (38). With this model, the assumptions previously made can be verified, namely the separation of variables and the harmonic vibration (Equation (41)).

Extrapolation to approximate the exact result
The analysis with the FDM is an approximation. Generally, the accuracy of the results increases by increasing the number of grid points. This means that when the number of points is infinitely high, the results tend toward the exact result. However, in first-order analysis, the exact result is rapidly obtained since the FDM polynomial approximation can match the exact results.
We assume that the relationship between the results R and the number of grid points on the beam N follows a hyperbolic curve with the constants A, B, and C, as follows: The exact result Re is approximated when the number of grid points N → : (54)

Beam subjected to a uniformly distributed load
We analyzed a uniform fixed−pinned beam subjected to a uniformly distributed load, as shown in Figure 18. The governing equation (Equation (6)) is applied at grid points 1, 2, 3, 4, and 5. The boundary conditions are satisfied using Equations (7a) and (7c).
Details of the analysis and results are presented in the supplementary material "fixed−pinned beam subjected to a uniformly distributed load". Table 1 lists the results obtained with the classical beam theory (CBT) and those obtained in the present study. Furthermore, the results are presented for a three-point stencil (TPS) considered for the slope (Equation (14b)) and bending moment (Equation (14a)) when satisfying the boundary conditions. Finally the results for a number of grid points n = 4, 3, and 2 are presented.  The results of the present study are exact for a uniformly distributed load whatever the discretization, since the exact solution for the deflection curve is here a fourth-order polynomial which corresponds to the FDM approximation.
It is noted that the use of a three-point stencil for the bending moment and slope yields less accurate results since a five-point stencil is used for the governing equations.

Beam subjected to a concentrated load
We analyzed a uniform fixed−pinned beam subjected to a concentrated load, as shown in Figure 19. Figure 19. Uniform fixed−pinned beam subjected to a concentrated load The model showing the grid points, as represented in Figure 4a,b, is considered.
Details of the analysis and results are presented in the supplementary material "fixed−pinned beam subjected to a concentrated load". Table 2 lists the results obtained with the classical beam theory (CBT) and those obtained in the present study (FDM). The results are also presented for a three-point discretization of the beam. The results of the present study are exact for a concentrated load whatever the discretization, since the exact solution for the deflection curve is here a third-order polynomial which is exactly described with the fourth-order polynomial FDM approximation.

Beam subjected to a linearly distributed load
We analyzed here a uniform fixed−pinned beam subjected to a linearly distributed load, as shown in Figure 20. The beam is calculated at one hand with a five-point grid, and at another hand with a six-point grid.
Details of the analysis and results are presented in the supplementary material "fixed−pinned beam subjected to a linearly distributed load". Table 3 lists the results obtained with the classical beam theory (the exact results) and those obtained in the present study (FDM). The results of the present study have a high accuracy. It is noted that the exact results cannot be achieved since the exact solution of w(x) for a linearly distributed loading is a fifth-order polynomial and the FDM approximation is a fourth-order polynomial. However the accuracy increases with increasing number of grid points.

Tapered pinned−fixed beam subjected to a uniformly distributed load
We analyzed a tapered pinned−fixed beam subjected to a uniformly distributed load, as shown in Figure 21.
I1 is the second moment of area at x1 = L1. L = 8.0m and L0 = 2.0m First, the beam is calculated with the force method of the classical beam theory (exact results). Then, the calculation is conducted with the FDM using n = 9, 13, and 17 grid points. The results are extrapolated to obtain those for infinite grid points. Details of the analysis and results are presented in Appendix A and in the supplementary material "tapered pinned −fixed beam subjected to a uniformly distributed load". Table 4 lists the results obtained with the classical beam theory (the exact results) and those obtained in the present study (FDM). The results of the present study have a high accuracy.

Beam subjected to a uniformly distributed load and a compressive force
We analyzed a uniform fixed−free beam subjected to a uniformly distributed load and a compressive force, as shown in Figure 22. Fixed−free beam subjected to a uniformly distributed load and a compressive force.
The governing equation (Equation (6)) is applied at grid points.
Details of the analysis and results are presented in the supplementary material "fixed−free beam subjected to a uniformly distributed load and compressive force". The analysis is conducted with n = 9, 13, and 17 grid points. The results are then extrapolated to obtain those for infinite grid points (Equation (54)). Table 5 lists the results obtained with the classical beam theory (CBT) and those obtained in the present study (FDM).

Buckling load of a fixed−pinned beam
We determined the buckling load of a fixed−pinned beam, as shown in Figure 23.
Values of the buckling factor  are listed in Table 6. The results of the present study have a high accuracy.

Buckling load of a tapered beam
We determined here the buckling loads of tapered beams with various support conditions, as shown in Figure 21.
The analysis is conducted with n = 9, 13, and 17 grid points. The results are then extrapolated to obtain those for infinite grid points. Details of the analysis and the results are listed in the supplementary material "stability of a tapered beam". The buckling load Ncr is defined as follows: The buckling factors  are listed in Table 7

Free vibration analysis of a fixed−fixed beam
We determined here the vibration frequencies of a fixed−fixed beam.
The analysis is conducted with n = 9, 13, and 17 grid points. The results are then extrapolated to obtain those for infinite grid points. The details of the analysis and results are listed in the supplementary file "vibration analysis of a fixed−fixed beam". The vibration frequency is . . .
The coefficients  are listed in Table 8 below. The results of the present study have a high accuracy.

Free vibration analysis of a tapered free−fixed beam
We determined the vibration frequencies of the tapered beam represented in Figure 24.  Table 9 lists the results obtained by Torabi [3] and those obtained in the present study. The results of the present study are identical to those presented by Torabi [3].

Conclusions
The FDM-based model developed in this paper enables, with relative easiness, first-order analysis, second-order analysis, and vibration analysis of Euler−Bernoulli beams. The results showed that the calculations conducted as described in this paper yielded accurate results. First-and second-order element stiffness matrices (the axial force being tensile or compressive) in local coordinates were determined. The analysis of tapered beams was also conducted.
The following aspects were not addressed in this study but could be analyzed with the model in future research: ✓ Analysis of Timoshenko beams.
✓ Analysis of linear structures, such as frames, through the transformation of element stiffness matrices from local coordinates in the global coordinates.

( ) / A x A x L =
In the associated statically determinate system, M0(x) and m(x) are the bending moment due to the distributed load and to the virtual unit moment at the fixed-end, respectively.
Let us introduce the dimensionless ordinate  = x/l and  0 = L 0 /L 1 .
M0(x), m(x), and I(x) can be expressed as follows (A1) The bending moment M1 at the fixed end is the solution of the following equations: