# Finite Element Formulation for Stability and Free Vibration Analysis of Timoshenko Beam

Authors: Abbas Moallemi-Oreh; Mohammad Karkon
Journal: Advances in Acoustics and Vibration
 (2013)

---

## Abstract

A two-node element is suggested for analyzing the stability and free vibration of Timoshenko beam. Cubic displacement polynomial and quadratic rotational fields are selected for this element. Moreover, it is assumed that shear strain of the  element has the constant value. Interpolation functions for displacement field and beam rotation are exactly calculated by employing total beam energy and its stationing to shear strain. By exploiting these interpolation functions, beam elements' stiffness matrix is also examined. Furthermore, geometric stiffness matrix and mass matrix of the proposed element are calculated by writing governing equation on stability and beam free vibration. At last, accuracy and efficiency of proposed element are evaluated through numerical tests. These tests show high accuracy of the element in analyzing beam stability and finding its critical load and free vibration analysis.


---

## Body

## 1. Introduction

Two versions of theories have been developed for analysis of beams. In Euler-Bernoulli theory, the displacement of beam is considered without shear effects. This method gives appropriate and acceptable response in thin beam in which shear effect is insignificant. However, in this approach, by increasing the thickness of beam and shear effect deformation, the error of response is increasing [1]. Correspondingly, the effect of shear transformation is formulated in Timoshenko theory. Therefore, this method has a better result, especially in deep beams in which shear effect is impressive. Although the rotational inertia of thick beams was investigated by Rayleigh for the first time, Timoshenko has developed this theory and formulated shear effect. Due to the complexity of the governing equations of the free vibration and stability of beams in general, numerical methods such as finite element have been developed profoundly. Up to now, many elements have been presented based on Timoshenko theory. These elements are classified into two groups which are simple and higher order elements. Some researchers used simple two-node elements with four degrees of freedom [2–4]. Thomas et al. have examined the elements proposed by other researchers [3]. The first high-order element was proposed by Kapur with eight degrees of freedom [5]. Lees and Thomas formulated a complex element by applying independent polynomial series for displacement and rotation fields [6, 7]. Also, this method has been used by Webster [8]. Rao and Gupta have examined free vibration of rotating beams [9]. In some methods, like isoparametric formulation, displacement and rotation fields are assumed dependently with the same order [10]. Based on Euler-Bernoulli theory, Goncalves et al. have presented frequency equation and vibration modes for classical boundary conditions such as clamped, free, pinned, and sliding supports [11]. Lee and Schultz have considered free vibration of Timoshenko beam through pseudospectral method [12].So far, very little research has been done on buckling of Timoshenko beam with respect to free vibration analysis. Exploiting an approximate method based on finite elements formulation, Wieckowski and Golubiewski examined beam stability of Euler-Bernoulli and Timoshenko theories [13]. Also, Kosmatka has proposed a two-node element for stability and free vibration analysis of Timoshenko beam [14].In this study, a new beam element is proposed for free vibration and stability analysis of Timoshenko beams. In order to compute shape function of the beam element, cubic displacement polynomial and quadratic rotational fields are selected. In addition, shear strain of the element is assumed constant. Then, by exploiting the bending and shear strain energy of the beam and stationary with respect to constant shear strain, interpolation function of this element has been carefully calculated. In the following, by using these shape functions, the stiffness matrix, geometric stiffness matrix, and mass matrix of the proposed element have been exactly determined. Finally, several numerical tests are performed to investigate the robustness of this element for free vibration and stability analysis of beams with different boundary conditions. The findings prove that the suggested element has high level of accuracy and free of shear locking.


## 2. Finite Element Formulation

In the finite element method, displacement and rotation fields of the element are associated with interpolation functions to nodal degrees of freedom. Figure1 shows proposed element with two nodes. The shape functions of the Timoshenko beam are calculated based on Figure 1. In order to compute shape function of the beam in Figure 1, cubic displacement polynomial and quadratic rotational fields are selected. Additionally, it is assumed that shear strain has the constant value of γ0. Based on these assumptions, (1) can be written as follows:
(1)w=wi2(1-s)+wj2(1+s)+β0l(1-s2)+β1ls(1-s2)θ=θi2(1-s)+θj2(1+s)+α0(1-s2)γ=γ0s=2xl-1.
In these equations, β1, β0, α0, and γ0 are unknown parameters. In order to determine their values, at first, the equation of shear strain for Timoshenko beam is established. By utilizing the shear strain value equal to γ0, the subsequent equations will be available:
(2)γ=dwdx-θ=2l·dwds-θγ0=2l(-wi2+wj2-2β0ls+β1l-3β1ls2)-θi(1-s2)-θj(1+s2)-α0(1-s2).
In the present formula, the coefficients of the terms s and s2 are equivalent to zero. Therefore, in the succeeding lines, α0, β1 are determined in terms of the unknown parameter γ0:
(3)Γ=2l(wj-wi)-(θi+θj)β0=18(θi-θj)α0=-32(γ0-12Γ)β1=16α0.
At this stage, there is only one unknown constant γ0, which can be discovered through the condition of minimum strain energy. It should be added that the structural strain energy is the sum of bending and shear strain energy. Bending strain energy is calculated in the following way:
(4)Ub=EI2∫0lκ2dx=EIl4∫-11κ2ds.
In (4), κ represents curvature which is determined as follows:
(5)κ=-2l·dθds=κ0-6sγ0lκ0=1l(θi-θj+3sΓ).
Substituting (5) into (4) leads to the following bending strain energy:
(6)Ub=U0+6EI(-Γγ0l+γ02l)U0=EIl4∫-11κ02ds.
Besides, (7) shows the energy of shear strain:
(7)Us=GA2fs∫01γ2dx=GAl4fs∫-11γ02ds=GAl2fsγ02.
In this equation, fs is shear correction factor which depends on cross section shape of beam. This coefficient for rectangular section is 5/6. By adding the bending and shear strain energy together, total strain energy is found as follows:
(8)U=Ub+Us=U0-6EIΓlγ0+6EIlγ02+GAl2fsγ02.
Implementing ∂U/∂γ0=0 will give the following results:
(9)γ0=6EIΓGAl2/fs+12EI=δΓδ=6λl2+12λ,λ=fsEIGA.
Substituting β1, β0, α0, and γ0 into (1), the succeeding shape functions for Timoshenko beam can be as follows:
(10){w~θ~n}=[N1N5N2N6N3N7N4N8]{wiθniwjθnj}N1=14[2+s3(1-2δ)+s(-3+2δ)]N2=l4[0.5(1-s2)+(s3-s)(0.5-δ)]N3=14[2-s3(1-2δ)-s(-3+2δ)]N4=l4[-0.5(1-s2)+(s3-s)(0.5-δ)]N5=14l[6(1-s2)(-1+2δ)]N6=14[-1+s(-2+3s)+6(1-s2)δ]N7=14l[-6(1-s2)(-1+2δ)]N8=14[-1+s(2+3s)+6(1-s2)δ].
In finite element method, displacements and strain field of the element can be related to the nodal degrees of freedom with interpolation functions. Hence, equations of finite element formulation can be written as follows:
(11){yθ}=[N]{D}E{ε}={dydx-dydx+θ}=[B]{D}E,
where
(12)[B]=[0ddxddx-1][N]
and [N] is the matrix of interpolation function. Also, {D}E is nodal displacement and [B] is strain matrix. Consequently, stiffness matrix of Timoshenko element can be obtained as follows:
(13)[K0]=∫[B]T[Dm][B]dx,
where, [Dm] is the elasticity matrix for Timoshenko beam element:
(14)[Dm]=[EI00GAfs].
By calculating (13), the stiffness matrix of the proposed element is obtained in the succeeding forms:
(15)[K0]=EIl3+12lλ×[126l-126l6l4l2+12λ-6l2l2-12λ-12-6l12-6l6l2l2-12λ-6l4l2+12λ].Figure 1
Timoshenko beam element.


## 3. Mass Matrix

The kinetic energy,T, of an elemental length, l, of a uniform Timoshenko beam is given as follows [1]:
(16)T=12∫-l/2l/2ρAy˙2dx+12∫-l/2l/2ρIθ˙2dx.
In this equation, ρ is the mass density of the material of the beam and I is the second moment of area of cross section. Therefore, the mass matrix of the element has been two parts, one related to translations and the other related to rotations, in the form of
(17)[M]=[M1]+[M2]=l2∫11ρA[Nw]T[Nw]ds+l2∫-l/2l/2ρI[Nθ]T[Nθ]Tds.
The translation mass matrix [M1] is achieved as follows:
(18)[M1]=ρAl5210(12λ+l2)2×[m11m12m13m14m12m15-m14m16m13-m14m11-m12m14m16-m12m15].
In the following, the entries of this matrix are introduced:
(19)m11=6l4(1680λ2+294l2λ+13l4)m12=1l3(1260λ2+231l2λ+11l4)m13=9l4(560λ2+84l2λ+3l4)m14=-12l3(2520λ2+378l2λ+13l4)m15=2l2(126λ2+21l2λ+13l4)m16=-32l2(168λ2+28l2λ+13l4).
In addition, the rotation mass matrix [M2] can be obtained as follows:
(20)[M2]=ρIl230(12λ+l2)2×[m21m22-m21m22m22m23-m22m24-m21-m22m21-m12m22m24-m12m23].
The entries of this matrix are defined in the below form:
(21)m21=36lm22=-3(60λ-l2)m23=4l(360λ2+15l2λ+3l4)m24=1l(720λ2-60l2λ-l4).


## 4. Geometric Stiffness Matrix

The concept of the neutral state of equilibrium is used for buckling analysis of beam. In the mathematical formulation of elastic stability of beam, the neutral equilibrium is used assuming a bifurcation of the deformations. That is, at the critical load, of the possible two paths of deformations (one associated with the stable equilibrium and the other pertinent to the unstable equilibrium condition, as shown in Figure2), the beam always takes the buckled form. In addition to the existence of this bifurcation of equilibrium paths, the elastic stability analysis of plates assumes the validity of Hooke's law. By considering Figure 2 at the buckled form, the axial shortening of beam can be acquired as follows:
(22)ds=(dx2+dy2)≃dx+12(dydx)2dx⟹dδ=12(dydx)2dx.
As a result, the strain energy of axial load P can be expressed as follows:
(23)ΔW=P2∫0l(dydx)2dx.
Based on (23), geometric stiffness matrix of the element can be calculated as follows:
(24)[Kg]=∫0l[dNwdx]TP[dNwdx]dx.
By calculating this equation, the geometric stiffness matrix becomes as follows:
(25)[Kg]=P60l×[kg1kg2kg3kg2kg2kg4-kg2kg5kg3-kg2kg1-kg2kg2kg5-kg2kg4],
where
(26)kg1=60+12βkg2=6βlkg3=-60-12βkg4=5l2+3βl2kg5=3βl2-5l2β=l4(l2+12λ)2.
For stability analysis and determination of the magnitude of a static compressive axial load that will produce beam buckling, the following eigenvalue will be achieved:
(27)det([K])=det([K0]-P[Kg])=0⟹Pcr.
The lowest positive eigenvalue of this equation is the magnitude of buckling load and the corresponding eigenvector is the deformed shape of the buckled beam. The exact solution of beam-column buckling load with shear deformation effect is obtained as succeeding form [15]:
(28)Pcr=π2EILeff2×(11+(π2fsEI)/(Leff2GA)),
where Leff is the effective beam length in which (Leff=L) and (Leff=L/2) are used for pinned-pinned beams and fixed-fixed beams, respectively.Figure 2
Straight and buckled form.


### 4.1. First Example (Vibration Analysis)

The efficiency of proposed element is evaluated by analyzing the free vibration of beam with simply and clamped supports for various lengths to thickness ratio. The Poisson’s ratio of this beam isν=0.3 and shear correction factor is taken 5/6. To compare other researchers’ results, frequency dimensionless parameter λi, defined in (29), has been shown in (Tables 1, 2, 3, 4, 5, and 6) for 10 first frequencies of this beam by exploiting 40 elements:
(29)λi2=ωil2ρAEI
Tables 1, 2, 3, and 4, 5, 6 show dimensionless parameters of natural frequency of beam-free vibration for simply and clamped supports, respectively, considering three ratios of thickness to length. Results have been compared to other researchers’ studies. The mentioned tables demonstrate that the accuracy of the proposed element is very high in analyzing free vibration of beam.Table 1
Dimensionless frequency parameterλi for the simply supported Timoshenko beam.
Mode
Euler theory
h
/
l
=
0.002
Ferreira [10]
Lee and Schultz [12]
Proposed element
1
3.14159
3.1428
3.14158
3.14158
2
6.28319
6.2928
6.28310
6.28310
3
9.42478
9.4573
9.42449
9.42450
4
12.5664
12.6437
12.5657
12.5657
5
15.7080
15.8596
15.7066
15.7068
6
18.8496
19.1127
18.8473
18.8476
7
21.9911
22.4113
21.9875
21.9883
8
25.1327
25.7638
25.1273
25.1288
9
28.2743
29.1793
28.2666
28.2692
10
31.4159
32.6672
31.4053
31.4098Table 2
Dimensionless frequency parameterλi for the simply supported Timoshenko beam.
Mode
Euler theory
h
/
l
=
0.01
Ferreira [10]
Lee and Schultz [12]
Proposed element
1
3.14159
3.1425
3.14133
3.14133
2
6.28319
6.2908
6.28106
6.28106
3
9.42478
9.4503
9.41761
9.41764
4
12.5664
12.6271
12.5494
12.5496
5
15.7080
15.8267
15.6749
15.6754
6
18.8496
19.0552
18.7926
18.7937
7
21.9911
22.3186
21.9011
21.9034
8
25.1327
25.6231
24.9988
25.0034
9
28.2743
28.9749
28.0845
28.0927
10
31.4159
32.3806
31.1568
31.1705Table 3
Dimensionless frequency parameterλi for the simply supported Timoshenko beam.
Mode
Euler theory
h
/
l
=
0.1
Ferreira [10]
Lee and Schultz [12]
Proposed element
1
3.14159
3.1169
3.11568
3.11569
2
6.28319
6.0993
6.09066
6.09094
3
9.42478
8.8668
8.84052
8.84229
4
12.5664
11.3984
11.3431
11.3492
5
15.7080
13.7089
13.6132
13.6282
6
18.8496
15.8266
15.6790
15.7093
7
21.9911
17.7811
17.5705
17.6239
8
25.1327
19.5991
19.3142
19.3997
9
28.2743
21.3030
20.9325
21.0606
10
31.4159
22.9117
22.4441
22.6257Table 4
Dimensionless frequency parameterλi for the clamped supported Timoshenko beam.
Mode
Euler theory
h
/
l
=
0.002
Ferreira [10]
Lee and Schultz [12]
Proposed element
1
4.73004
4.7345
4.72998
4.72998
2
7.85320
7.8736
7.85295
7.85296
3
10.9956
11.0504
10.9950
10.9950
4
14.1372
14.2526
14.1359
14.1360
5
17.2788
17.4888
17.2766
17.2768
6
20.4204
20.7670
20.4168
20.4174
7
23.5619
24.0955
23.5567
23.5578
8
26.7035
27.4833
26.6960
26.6980
9
29.8451
30.9398
29.8348
29.8382
10
32.9867
34.4748
32.9729
32.9786Table 5
Dimensionless frequency parameterλi for the clamped supported Timoshenko beam.
Mode
Euler theory
h
/
l
=
0.01
Ferreira [10]
Lee and Schultz [12]
Proposed element
1
4.73004
4.7330
4.72840
4.72840
2
7.85320
7.8675
7.84690
7.84692
3
10.9956
11.0351
10.9800
10.9801
4
14.1372
14.2218
14.1062
14.1064
5
17.2788
17.4342
17.2246
17.2253
6
20.4204
20.6783
20.3338
20.3355
7
23.5619
23.9600
23.4325
23.4358
8
26.7035
27.2857
26.5192
26.5253
9
29.8451
30.6616
29.5926
29.6032
10
32.9867
34.0944
32.6514
32.6687Table 6
Dimensionless frequency parameterλi for the clamped supported Timoshenko beam.
Mode
Euler theory
h
/
l
=
0.1
Ferreira [10]
Lee and Schultz [12]
Proposed element
1
4.73004
4.5835
4.57955
4.57962
2
7.85320
7.3468
7.33122
7.33193
3
10.9956
9.8924
9.85611
9.85918
4
14.1372
12.2118
12.1454
12.1540
5
17.2788
14.3386
14.2324
14.2513
6
20.4204
16.3046
16.1487
16.1841
7
23.5619
18.1375
17.9215
17.9807
8
26.7035
19.8593
19.5723
19.6641
9
29.8451
21.4875
21.1185
21.2523
10
32.9867
23.0358
22.5735
22.7598


### 4.2. Second Example (Buckling Analysis)

The efficiency of the proposed element is determined in buckling analysis. For this propose, critical load of a simply and clamped supported beam-column is computed. The modulus of elasticity, Poisson’s ratio, and length of this beam-column areE=10e7,  ν=1/3, and l=1, respectively. The buckling loads of simply and clamped supported beam-columns are listed in Table 7. It demonstrates that the proposed element has high accuracy in buckling analysis. Figure 3 shows the first three modes of buckling for simply and clamped supported beam-column, respectively.Table 7
Critical load of the simply and clamped supported beam-column.
l
/
h
Simply supported
Clamped supported
Analytical solution
Ferreira [10]
Proposed element
Analytical solution
Ferreira [10]
Proposed element
10
8013.8
8021.8
8013.86
29766
29877
29770
100
8.223
8.231
8.2225
32.864
32.999
32.864
1000
0.0082
0.0082
0.00822
0.0329
0.0330
0.0329The first three modes of buckling for beam-column, (a) simply supported and (b) clamped supported.
(a)
(b)


## 4.1. First Example (Vibration Analysis)

The efficiency of proposed element is evaluated by analyzing the free vibration of beam with simply and clamped supports for various lengths to thickness ratio. The Poisson’s ratio of this beam isν=0.3 and shear correction factor is taken 5/6. To compare other researchers’ results, frequency dimensionless parameter λi, defined in (29), has been shown in (Tables 1, 2, 3, 4, 5, and 6) for 10 first frequencies of this beam by exploiting 40 elements:
(29)λi2=ωil2ρAEI
Tables 1, 2, 3, and 4, 5, 6 show dimensionless parameters of natural frequency of beam-free vibration for simply and clamped supports, respectively, considering three ratios of thickness to length. Results have been compared to other researchers’ studies. The mentioned tables demonstrate that the accuracy of the proposed element is very high in analyzing free vibration of beam.Table 1
Dimensionless frequency parameterλi for the simply supported Timoshenko beam.
Mode
Euler theory
h
/
l
=
0.002
Ferreira [10]
Lee and Schultz [12]
Proposed element
1
3.14159
3.1428
3.14158
3.14158
2
6.28319
6.2928
6.28310
6.28310
3
9.42478
9.4573
9.42449
9.42450
4
12.5664
12.6437
12.5657
12.5657
5
15.7080
15.8596
15.7066
15.7068
6
18.8496
19.1127
18.8473
18.8476
7
21.9911
22.4113
21.9875
21.9883
8
25.1327
25.7638
25.1273
25.1288
9
28.2743
29.1793
28.2666
28.2692
10
31.4159
32.6672
31.4053
31.4098Table 2
Dimensionless frequency parameterλi for the simply supported Timoshenko beam.
Mode
Euler theory
h
/
l
=
0.01
Ferreira [10]
Lee and Schultz [12]
Proposed element
1
3.14159
3.1425
3.14133
3.14133
2
6.28319
6.2908
6.28106
6.28106
3
9.42478
9.4503
9.41761
9.41764
4
12.5664
12.6271
12.5494
12.5496
5
15.7080
15.8267
15.6749
15.6754
6
18.8496
19.0552
18.7926
18.7937
7
21.9911
22.3186
21.9011
21.9034
8
25.1327
25.6231
24.9988
25.0034
9
28.2743
28.9749
28.0845
28.0927
10
31.4159
32.3806
31.1568
31.1705Table 3
Dimensionless frequency parameterλi for the simply supported Timoshenko beam.
Mode
Euler theory
h
/
l
=
0.1
Ferreira [10]
Lee and Schultz [12]
Proposed element
1
3.14159
3.1169
3.11568
3.11569
2
6.28319
6.0993
6.09066
6.09094
3
9.42478
8.8668
8.84052
8.84229
4
12.5664
11.3984
11.3431
11.3492
5
15.7080
13.7089
13.6132
13.6282
6
18.8496
15.8266
15.6790
15.7093
7
21.9911
17.7811
17.5705
17.6239
8
25.1327
19.5991
19.3142
19.3997
9
28.2743
21.3030
20.9325
21.0606
10
31.4159
22.9117
22.4441
22.6257Table 4
Dimensionless frequency parameterλi for the clamped supported Timoshenko beam.
Mode
Euler theory
h
/
l
=
0.002
Ferreira [10]
Lee and Schultz [12]
Proposed element
1
4.73004
4.7345
4.72998
4.72998
2
7.85320
7.8736
7.85295
7.85296
3
10.9956
11.0504
10.9950
10.9950
4
14.1372
14.2526
14.1359
14.1360
5
17.2788
17.4888
17.2766
17.2768
6
20.4204
20.7670
20.4168
20.4174
7
23.5619
24.0955
23.5567
23.5578
8
26.7035
27.4833
26.6960
26.6980
9
29.8451
30.9398
29.8348
29.8382
10
32.9867
34.4748
32.9729
32.9786Table 5
Dimensionless frequency parameterλi for the clamped supported Timoshenko beam.
Mode
Euler theory
h
/
l
=
0.01
Ferreira [10]
Lee and Schultz [12]
Proposed element
1
4.73004
4.7330
4.72840
4.72840
2
7.85320
7.8675
7.84690
7.84692
3
10.9956
11.0351
10.9800
10.9801
4
14.1372
14.2218
14.1062
14.1064
5
17.2788
17.4342
17.2246
17.2253
6
20.4204
20.6783
20.3338
20.3355
7
23.5619
23.9600
23.4325
23.4358
8
26.7035
27.2857
26.5192
26.5253
9
29.8451
30.6616
29.5926
29.6032
10
32.9867
34.0944
32.6514
32.6687Table 6
Dimensionless frequency parameterλi for the clamped supported Timoshenko beam.
Mode
Euler theory
h
/
l
=
0.1
Ferreira [10]
Lee and Schultz [12]
Proposed element
1
4.73004
4.5835
4.57955
4.57962
2
7.85320
7.3468
7.33122
7.33193
3
10.9956
9.8924
9.85611
9.85918
4
14.1372
12.2118
12.1454
12.1540
5
17.2788
14.3386
14.2324
14.2513
6
20.4204
16.3046
16.1487
16.1841
7
23.5619
18.1375
17.9215
17.9807
8
26.7035
19.8593
19.5723
19.6641
9
29.8451
21.4875
21.1185
21.2523
10
32.9867
23.0358
22.5735
22.7598


## 4.2. Second Example (Buckling Analysis)

The efficiency of the proposed element is determined in buckling analysis. For this propose, critical load of a simply and clamped supported beam-column is computed. The modulus of elasticity, Poisson’s ratio, and length of this beam-column areE=10e7,  ν=1/3, and l=1, respectively. The buckling loads of simply and clamped supported beam-columns are listed in Table 7. It demonstrates that the proposed element has high accuracy in buckling analysis. Figure 3 shows the first three modes of buckling for simply and clamped supported beam-column, respectively.Table 7
Critical load of the simply and clamped supported beam-column.
l
/
h
Simply supported
Clamped supported
Analytical solution
Ferreira [10]
Proposed element
Analytical solution
Ferreira [10]
Proposed element
10
8013.8
8021.8
8013.86
29766
29877
29770
100
8.223
8.231
8.2225
32.864
32.999
32.864
1000
0.0082
0.0082
0.00822
0.0329
0.0330
0.0329The first three modes of buckling for beam-column, (a) simply supported and (b) clamped supported.
(a)
(b)


## 5. Conclusion

This study has proposed a new beam finite element formulation for the stability and free vibration analysis of beams with shear effect deformation. For this purpose, displacement field of the element has been selected from the third degree, rotation field has been selected from the second degree, and shear strain is assumed constant value. By employing the bending and shear strain energy of the element and stationary respect to unknown shear strain, this value is obtained. In the following, using the shear strain, the interpolation functions for displacement, and rotation fields of element has been exactly calculated. Then, these interpolation functions have been used and stiffness matrix, geometric stiffness matrix, and mass matrix of the proposed element have been clearly obtained. Evaluating the efficiency and accuracy of the element for free vibration and stability analysis of beam with simply and clamped supports, desirable results are obtained. The results show high accuracy and efficiency of the proposed element in calculating natural frequencies and critical load of beam with different boundary conditions.



---

*Source: 841215-2013-04-18.xml*