Hypotheses:
Strain Small
Displacement Small
Material Linear
Timoshenko beam model.
Beam element: straight, prismatic (constant cross-section), two nodes.
The local reference frame contains the beam centroidal axis, the local x ˉ \bar x x ˉ y ˉ \bar y y ˉ z ˉ \bar z z ˉ
In the local reference frame, these DOFs are three displacements
( u ˉ , v ˉ , w ˉ )
(\bar u,\bar v,\bar w)
( u ˉ , v ˉ , w ˉ ) and three rotations of the cross-section
( ψ x ˉ , ψ y ˉ , ψ z ˉ ) .
(\psi_{\bar x},\psi_{\bar y},\psi_{\bar z}).
( ψ x ˉ , ψ y ˉ , ψ z ˉ ) . For the Timoshenko beam, the nodal rotations are the rotations of the cross-section ψ x ˉ , ψ y ˉ , ψ z ˉ \psi_{\bar x},\psi_{\bar y},\psi_{\bar z} ψ x ˉ , ψ y ˉ , ψ z ˉ φ x ˉ , φ y ˉ , φ z ˉ \varphi_{\bar x},\varphi_{\bar y},\varphi_{\bar z} φ x ˉ , φ y ˉ , φ z ˉ
Figure 1. The 3D Timoshenko beam is described by the following two groups of equations in the local reference frame.
For the x ˉ y ˉ \bar x\bar y x ˉ y ˉ
{ d ψ z ˉ d x ˉ = M z ˉ E I z , β y = T y G A y , φ z ˉ = d v ˉ d x ˉ , φ z ˉ = ψ z ˉ + β y .
\left\{
\begin{array}{l}
\dfrac{d\psi_{\bar z}}{d\bar x}
=
\dfrac{M_{\bar z}}{EI_z},
\\[3mm]
\beta_y=\dfrac{T_y}{GA_y},
\\[3mm]
\varphi_{\bar z}=\dfrac{d\bar v}{d\bar x},
\\[3mm]
\varphi_{\bar z}=\psi_{\bar z}+\beta_y .
\end{array}
\right.
⎩ ⎨ ⎧ d x ˉ d ψ z ˉ = E I z M z ˉ , β y = G A y T y , φ z ˉ = d x ˉ d v ˉ , φ z ˉ = ψ z ˉ + β y . For the x ˉ z ˉ \bar x\bar z x ˉ z ˉ Section 23.1 :
{ d ψ y ˉ d x ˉ = M y ˉ E I y , β z = T z G A z , φ y ˉ = − d w ˉ d x ˉ , φ y ˉ = ψ y ˉ − β z .
\left\{
\begin{array}{l}
\dfrac{d\psi_{\bar y}}{d\bar x}
=
\dfrac{M_{\bar y}}{EI_y},
\\[3mm]
\beta_z=\dfrac{T_z}{GA_z},
\\[3mm]
\varphi_{\bar y}
=
-\dfrac{d\bar w}{d\bar x},
\\[3mm]
\varphi_{\bar y}
=
\psi_{\bar y}-\beta_z .
\end{array}
\right.
⎩ ⎨ ⎧ d x ˉ d ψ y ˉ = E I y M y ˉ , β z = G A z T z , φ y ˉ = − d x ˉ d w ˉ , φ y ˉ = ψ y ˉ − β z . β y \beta_y β y β z \beta_z β z φ y ˉ \varphi_{\bar y} φ y ˉ φ z ˉ \varphi_{\bar z} φ z ˉ
The interpolation functions for u ˉ , v ˉ , w ˉ \bar u,\bar v,\bar w u ˉ , v ˉ , w ˉ Section 23.1 . However, in the Timoshenko beam model the nodal rotational DOFs are the rotations of the cross-section ψ y ˉ , ψ z ˉ \psi_{\bar y},\psi_{\bar z} ψ y ˉ , ψ z ˉ φ y ˉ , φ z ˉ \varphi_{\bar y},\varphi_{\bar z} φ y ˉ , φ z ˉ
Because forces are applied only at the nodes, there are no distributed transverse loads along the element. Therefore the shear forces, and consequently the shear angles, are constant along the beam element:
d β y d x ˉ = 0 , d β z d x ˉ = 0.
\frac{d\beta_y}{d\bar x}=0,
\qquad
\frac{d\beta_z}{d\bar x}=0.
d x ˉ d β y = 0 , d x ˉ d β z = 0. Thus,
d ψ z ˉ d x ˉ = d φ z ˉ d x ˉ , d ψ y ˉ d x ˉ = d φ y ˉ d x ˉ .
\frac{d\psi_{\bar z}}{d\bar x}
=
\frac{d\varphi_{\bar z}}{d\bar x},
\qquad
\frac{d\psi_{\bar y}}{d\bar x}
=
\frac{d\varphi_{\bar y}}{d\bar x}.
d x ˉ d ψ z ˉ = d x ˉ d φ z ˉ , d x ˉ d ψ y ˉ = d x ˉ d φ y ˉ . Remembering that
φ z ˉ = d v ˉ d x ˉ , φ y ˉ = − d w ˉ d x ˉ ,
\varphi_{\bar z}
=
\frac{d\bar v}{d\bar x},
\qquad
\varphi_{\bar y}
=
-\frac{d\bar w}{d\bar x},
φ z ˉ = d x ˉ d v ˉ , φ y ˉ = − d x ˉ d w ˉ , and using the equilibrium equations
T y = − d M z ˉ d x ˉ , T z = d M y ˉ d x ˉ ,
T_y=-\frac{dM_{\bar z}}{d\bar x},
\qquad
T_z=\frac{dM_{\bar y}}{d\bar x},
T y = − d x ˉ d M z ˉ , T z = d x ˉ d M y ˉ , it results:
T y = − E I z d 3 v ˉ d x ˉ 3 ,
T_y
=
-EI_z
\frac{d^3\bar v}{d\bar x^3},
T y = − E I z d x ˉ 3 d 3 v ˉ , T z = − E I y d 3 w ˉ d x ˉ 3 .
T_z
=
-EI_y
\frac{d^3\bar w}{d\bar x^3}.
T z = − E I y d x ˉ 3 d 3 w ˉ . Consequently,
β y = − E I z G A y d 3 v ˉ d x ˉ 3
\boxed{
\beta_y
=
-\frac{EI_z}{GA_y}
\frac{d^3\bar v}{d\bar x^3}
}
β y = − G A y E I z d x ˉ 3 d 3 v ˉ and
β z = − E I y G A z d 3 w ˉ d x ˉ 3 .
\boxed{
\beta_z
=
-\frac{EI_y}{GA_z}
\frac{d^3\bar w}{d\bar x^3}
}.
β z = − G A z E I y d x ˉ 3 d 3 w ˉ . T y , T z T_y,T_z T y , T z A y , A z A_y,A_z A y , A z
A y = α A , A z = α A ,
A_y=\alpha A,
\qquad
A_z=\alpha A,
A y = α A , A z = α A , where, for a rectangular cross-section,
α = 5 6 .
\alpha=\frac56.
α = 6 5 . A A A
The transverse displacements v ˉ \bar v v ˉ w ˉ \bar w w ˉ x ˉ \bar x x ˉ β y , β z \beta_y,\beta_z β y , β z
01
Interpolation functions in the local reference frame The element nodal displacement vector is
u ‾ e l = { u ˉ 1 v ˉ 1 w ˉ 1 ψ x ˉ 1 ψ y ˉ 1 ψ z ˉ 1 u ˉ 2 v ˉ 2 w ˉ 2 ψ x ˉ 2 ψ y ˉ 2 ψ z ˉ 2 } .
\overline{\htmlClass{fea-math-bold-italic}{\boldsymbol{u}}}_{el}
=
\left\{
\begin{array}{c}
\bar u_1\\
\bar v_1\\
\bar w_1\\
\psi_{\bar x1}\\
\psi_{\bar y1}\\
\psi_{\bar z1}\\
\bar u_2\\
\bar v_2\\
\bar w_2\\
\psi_{\bar x2}\\
\psi_{\bar y2}\\
\psi_{\bar z2}
\end{array}
\right\}.
u e l = ⎩ ⎨ ⎧ u ˉ 1 v ˉ 1 w ˉ 1 ψ x ˉ 1 ψ y ˉ 1 ψ z ˉ 1 u ˉ 2 v ˉ 2 w ˉ 2 ψ x ˉ 2 ψ y ˉ 2 ψ z ˉ 2 ⎭ ⎬ ⎫ . As in Section 23.1 , introduce
ξ = x ˉ L ,
\xi=\frac{\bar x}{L},
ξ = L x ˉ , where L L L
H 1 = 2 ξ 3 − 3 ξ 2 + 1 ,
H_1=2\xi^3-3\xi^2+1,
H 1 = 2 ξ 3 − 3 ξ 2 + 1 , H 2 = L ( ξ 3 − 2 ξ 2 + ξ ) ,
H_2=L(\xi^3-2\xi^2+\xi),
H 2 = L ( ξ 3 − 2 ξ 2 + ξ ) , H 3 = − 2 ξ 3 + 3 ξ 2 ,
H_3=-2\xi^3+3\xi^2,
H 3 = − 2 ξ 3 + 3 ξ 2 , H 4 = L ( ξ 3 − ξ 2 ) .
H_4=L(\xi^3-\xi^2).
H 4 = L ( ξ 3 − ξ 2 ) . The axial displacement and the torsional rotation remain linear:
u ˉ = ( 1 − ξ ) u ˉ 1 + ξ u ˉ 2 ,
\bar u
=
(1-\xi)\bar u_1+\xi\bar u_2,
u ˉ = ( 1 − ξ ) u ˉ 1 + ξ u ˉ 2 , ψ x ˉ = ( 1 − ξ ) ψ x ˉ 1 + ξ ψ x ˉ 2 .
\psi_{\bar x}
=
(1-\xi)\psi_{\bar x1}
+
\xi\psi_{\bar x2}.
ψ x ˉ = ( 1 − ξ ) ψ x ˉ 1 + ξ ψ x ˉ 2 . For the x ˉ y ˉ \bar x\bar y x ˉ y ˉ
φ z ˉ = ψ z ˉ + β y ,
\varphi_{\bar z}
=
\psi_{\bar z}+\beta_y,
φ z ˉ = ψ z ˉ + β y , therefore
v ˉ = H 1 v ˉ 1 + H 2 ( ψ z ˉ 1 + β y ) + H 3 v ˉ 2 + H 4 ( ψ z ˉ 2 + β y ) .
\bar v
=
H_1\bar v_1
+
H_2(\psi_{\bar z1}+\beta_y)
+
H_3\bar v_2
+
H_4(\psi_{\bar z2}+\beta_y).
v ˉ = H 1 v ˉ 1 + H 2 ( ψ z ˉ 1 + β y ) + H 3 v ˉ 2 + H 4 ( ψ z ˉ 2 + β y ) . For the x ˉ z ˉ \bar x\bar z x ˉ z ˉ
φ y ˉ = ψ y ˉ − β z ,
\varphi_{\bar y}
=
\psi_{\bar y}-\beta_z,
φ y ˉ = ψ y ˉ − β z , therefore
w ˉ = H 1 w ˉ 1 − H 2 ( ψ y ˉ 1 − β z ) + H 3 w ˉ 2 − H 4 ( ψ y ˉ 2 − β z ) .
\bar w
=
H_1\bar w_1
-
H_2(\psi_{\bar y1}-\beta_z)
+
H_3\bar w_2
-
H_4(\psi_{\bar y2}-\beta_z).
w ˉ = H 1 w ˉ 1 − H 2 ( ψ y ˉ 1 − β z ) + H 3 w ˉ 2 − H 4 ( ψ y ˉ 2 − β z ) . Introduce the shear parameters:
Φ y = 12 E I z G A y L 2
\boxed{
\mathit{\Phi}_y=
\frac{12EI_z}{GA_yL^2}
}
Φ y = G A y L 2 12 E I z and
Φ z = 12 E I y G A z L 2 .
\boxed{
\mathit{\Phi}_z=
\frac{12EI_y}{GA_zL^2}
}.
Φ z = G A z L 2 12 E I y . From the preceding equations, the constant shear angles can be expressed directly in terms of the nodal DOFs:
β y = Φ y 1 + Φ y [ v ˉ 2 − v ˉ 1 L − ψ z ˉ 1 + ψ z ˉ 2 2 ]
\boxed{
\beta_y
=
\frac{\mathit{\Phi}_y}{1+\mathit{\Phi}_y}
\left[
\frac{\bar v_2-\bar v_1}{L}
-
\frac{\psi_{\bar z1}+\psi_{\bar z2}}{2}
\right]
}
β y = 1 + Φ y Φ y [ L v ˉ 2 − v ˉ 1 − 2 ψ z ˉ 1 + ψ z ˉ 2 ] and
β z = Φ z 1 + Φ z [ w ˉ 2 − w ˉ 1 L + ψ y ˉ 1 + ψ y ˉ 2 2 ] .
\boxed{
\beta_z
=
\frac{\mathit{\Phi}_z}{1+\mathit{\Phi}_z}
\left[
\frac{\bar w_2-\bar w_1}{L}
+
\frac{\psi_{\bar y1}+\psi_{\bar y2}}{2}
\right].
}
β z = 1 + Φ z Φ z [ L w ˉ 2 − w ˉ 1 + 2 ψ y ˉ 1 + ψ y ˉ 2 ] . Substituting β y \beta_y β y β z \beta_z β z 6 × 12 6\times12 6 × 12
u ‾ ( x ˉ ) = N ( x ˉ ) u ‾ e l ,
\overline{\htmlClass{fea-math-bold-italic}{\boldsymbol{u}}}(\bar x)
=
\htmlClass{fea-math-bold-italic}{\boldsymbol{N}}(\bar x)\,
\overline{\htmlClass{fea-math-bold-italic}{\boldsymbol{u}}}_{el},
u ( x ˉ ) = N ( x ˉ ) u e l , where
u ‾ = { u ˉ v ˉ w ˉ ψ x ˉ ψ y ˉ ψ z ˉ } .
\overline{\htmlClass{fea-math-bold-italic}{\boldsymbol{u}}}
=
\left\{
\begin{array}{c}
\bar u\\
\bar v\\
\bar w\\
\psi_{\bar x}\\
\psi_{\bar y}\\
\psi_{\bar z}
\end{array}
\right\}.
u = ⎩ ⎨ ⎧ u ˉ v ˉ w ˉ ψ x ˉ ψ y ˉ ψ z ˉ ⎭ ⎬ ⎫ . 02
Generalized strains and matrix B \htmlClass{fea-math-bold-italic}{\boldsymbol{B}} B The generalized strain vector is
ε = { ε 0 β y β z κ x ˉ κ y ˉ κ z ˉ } = B u ‾ e l ,
\htmlClass{fea-math-bold-italic}{\boldsymbol{\varepsilon}}
=
\left\{
\begin{array}{c}
\varepsilon_0\\
\beta_y\\
\beta_z\\
\kappa_{\bar x}\\
\kappa_{\bar y}\\
\kappa_{\bar z}
\end{array}
\right\}
=
\htmlClass{fea-math-bold-italic}{\boldsymbol{B}}\,
\overline{\htmlClass{fea-math-bold-italic}{\boldsymbol{u}}}_{el},
ε = ⎩ ⎨ ⎧ ε 0 β y β z κ x ˉ κ y ˉ κ z ˉ ⎭ ⎬ ⎫ = B u e l , where
ε 0 = d u ˉ d x ˉ , κ x ˉ = d ψ x ˉ d x ˉ ,
\varepsilon_0
=
\frac{d\bar u}{d\bar x},
\qquad
\kappa_{\bar x}
=
\frac{d\psi_{\bar x}}{d\bar x},
ε 0 = d x ˉ d u ˉ , κ x ˉ = d x ˉ d ψ x ˉ , κ y ˉ = d ψ y ˉ d x ˉ , κ z ˉ = d ψ z ˉ d x ˉ .
\kappa_{\bar y}
=
\frac{d\psi_{\bar y}}{d\bar x},
\qquad
\kappa_{\bar z}
=
\frac{d\psi_{\bar z}}{d\bar x}.
κ y ˉ = d x ˉ d ψ y ˉ , κ z ˉ = d x ˉ d ψ z ˉ . Consistently with the convention adopted in Section 23.1 , the matrix B \htmlClass{fea-math-bold-italic}{\boldsymbol{B}} B
B = [ − 1 L 0 0 0 0 0 1 L 0 0 0 0 0 0 − Φ y L ( 1 + Φ y ) 0 0 0 − Φ y 2 ( 1 + Φ y ) 0 Φ y L ( 1 + Φ y ) 0 0 0 − Φ y 2 ( 1 + Φ y ) 0 0 − Φ z L ( 1 + Φ z ) 0 Φ z 2 ( 1 + Φ z ) 0 0 0 Φ z L ( 1 + Φ z ) 0 Φ z 2 ( 1 + Φ z ) 0 0 0 0 − 1 L 0 0 0 0 0 1 L 0 0 0 0 6 ( 1 − 2 ξ ) L 2 ( 1 + Φ z ) 0 − 4 − 6 ξ + Φ z L ( 1 + Φ z ) 0 0 0 − 6 ( 1 − 2 ξ ) L 2 ( 1 + Φ z ) 0 − 2 + 6 ξ + Φ z L ( 1 + Φ z ) 0 0 − 6 ( 1 − 2 ξ ) L 2 ( 1 + Φ y ) 0 0 0 − 4 − 6 ξ + Φ y L ( 1 + Φ y ) 0 6 ( 1 − 2 ξ ) L 2 ( 1 + Φ y ) 0 0 0 − 2 + 6 ξ + Φ y L ( 1 + Φ y ) ] .
\htmlClass{fea-math-bold-italic}{\boldsymbol{B}}=
\begin{bmatrix}
-\dfrac1L&0&0&0&0&0&
\dfrac1L&0&0&0&0&0
\\[3mm]
0&
-\dfrac{\mathit{\Phi}_y}{L(1+\mathit{\Phi}_y)}
&0&0&0&
-\dfrac{\mathit{\Phi}_y}{2(1+\mathit{\Phi}_y)}
&0&
\dfrac{\mathit{\Phi}_y}{L(1+\mathit{\Phi}_y)}
&0&0&0&
-\dfrac{\mathit{\Phi}_y}{2(1+\mathit{\Phi}_y)}
\\[3mm]
0&0&
-\dfrac{\mathit{\Phi}_z}{L(1+\mathit{\Phi}_z)}
&0&
\dfrac{\mathit{\Phi}_z}{2(1+\mathit{\Phi}_z)}
&0&0&0&
\dfrac{\mathit{\Phi}_z}{L(1+\mathit{\Phi}_z)}
&0&
\dfrac{\mathit{\Phi}_z}{2(1+\mathit{\Phi}_z)}
&0
\\[3mm]
0&0&0&
-\dfrac1L&0&0&
0&0&0&
\dfrac1L&0&0
\\[3mm]
0&0&
\dfrac{6(1-2\xi)}
{L^2(1+\mathit{\Phi}_z)}
&0&
-\dfrac{4-6\xi+\mathit{\Phi}_z}
{L(1+\mathit{\Phi}_z)}
&0&
0&0&
-\dfrac{6(1-2\xi)}
{L^2(1+\mathit{\Phi}_z)}
&0&
\dfrac{-2+6\xi+\mathit{\Phi}_z}
{L(1+\mathit{\Phi}_z)}
&0
\\[3mm]
0&
-\dfrac{6(1-2\xi)}
{L^2(1+\mathit{\Phi}_y)}
&0&0&0&
-\dfrac{4-6\xi+\mathit{\Phi}_y}
{L(1+\mathit{\Phi}_y)}
&0&
\dfrac{6(1-2\xi)}
{L^2(1+\mathit{\Phi}_y)}
&0&0&0&
\dfrac{-2+6\xi+\mathit{\Phi}_y}
{L(1+\mathit{\Phi}_y)}
\end{bmatrix}.
B = − L 1 0 0 0 0 0 0 − L ( 1 + Φ y ) Φ y 0 0 0 − L 2 ( 1 + Φ y ) 6 ( 1 − 2 ξ ) 0 0 − L ( 1 + Φ z ) Φ z 0 L 2 ( 1 + Φ z ) 6 ( 1 − 2 ξ ) 0 0 0 0 − L 1 0 0 0 0 2 ( 1 + Φ z ) Φ z 0 − L ( 1 + Φ z ) 4 − 6 ξ + Φ z 0 0 − 2 ( 1 + Φ y ) Φ y 0 0 0 − L ( 1 + Φ y ) 4 − 6 ξ + Φ y L 1 0 0 0 0 0 0 L ( 1 + Φ y ) Φ y 0 0 0 L 2 ( 1 + Φ y ) 6 ( 1 − 2 ξ ) 0 0 L ( 1 + Φ z ) Φ z 0 − L 2 ( 1 + Φ z ) 6 ( 1 − 2 ξ ) 0 0 0 0 L 1 0 0 0 0 2 ( 1 + Φ z ) Φ z 0 L ( 1 + Φ z ) − 2 + 6 ξ + Φ z 0 0 − 2 ( 1 + Φ y ) Φ y 0 0 0 L ( 1 + Φ y ) − 2 + 6 ξ + Φ y . 03
Hooke’s law The generalized stress vector is
σ = { N T y T z M x ˉ M y ˉ M z ˉ } .
\htmlClass{fea-math-bold-italic}{\boldsymbol{\sigma}}
=
\left\{
\begin{array}{c}
N\\
T_y\\
T_z\\
M_{\bar x}\\
M_{\bar y}\\
M_{\bar z}
\end{array}
\right\}.
σ = ⎩ ⎨ ⎧ N T y T z M x ˉ M y ˉ M z ˉ ⎭ ⎬ ⎫ . Hooke’s law in matrix form is
σ = D ε ,
\htmlClass{fea-math-bold-italic}{\boldsymbol{\sigma}}
=
\htmlClass{fea-math-bold-italic}{\boldsymbol{D}}
\htmlClass{fea-math-bold-italic}{\boldsymbol{\varepsilon}},
σ = D ε , with
D = [ E A 0 0 0 0 0 0 G A y 0 0 0 0 0 0 G A z 0 0 0 0 0 0 G I t 0 0 0 0 0 0 E I y 0 0 0 0 0 0 E I z ] .
\htmlClass{fea-math-bold-italic}{\boldsymbol{D}}=
\begin{bmatrix}
EA&0&0&0&0&0\\
0&GA_y&0&0&0&0\\
0&0&GA_z&0&0&0\\
0&0&0&GI_t&0&0\\
0&0&0&0&EI_y&0\\
0&0&0&0&0&EI_z
\end{bmatrix}.
D = E A 0 0 0 0 0 0 G A y 0 0 0 0 0 0 G A z 0 0 0 0 0 0 G I t 0 0 0 0 0 0 E I y 0 0 0 0 0 0 E I z . Thus,
M y ˉ = E I y κ y ˉ , M z ˉ = E I z κ z ˉ ,
M_{\bar y}=EI_y\kappa_{\bar y},
\qquad
M_{\bar z}=EI_z\kappa_{\bar z},
M y ˉ = E I y κ y ˉ , M z ˉ = E I z κ z ˉ , with no additional sign introduced in the constitutive law; the signs are contained in the kinematic definitions.
04
Element stiffness matrix The stiffness matrix in the local reference frame is
k ‾ e l = ∫ L B T D B d x ˉ .
\boxed{
\overline{\htmlClass{fea-math-bold-italic}{\boldsymbol{k}}}_{el}
=
\int_L
\htmlClass{fea-math-bold-italic}{\boldsymbol{B}}^T
\htmlClass{fea-math-bold-italic}{\boldsymbol{D}}
\htmlClass{fea-math-bold-italic}{\boldsymbol{B}}\,d\bar x
}.
k e l = ∫ L B T D B d x ˉ . Performing the integral analytically yields the following result [1] :
k ‾ e l = [ E A L 0 0 0 0 0 − E A L 0 0 0 0 0 0 k 11 z 0 0 0 k 12 z 0 k 13 z 0 0 0 k 14 z 0 0 k 11 y 0 k 12 y 0 0 0 k 13 y 0 k 14 y 0 0 0 0 G I t L 0 0 0 0 0 − G I t L 0 0 0 0 k 12 y 0 k 22 y 0 0 0 k 23 y 0 k 24 y 0 0 k 12 z 0 0 0 k 22 z 0 k 23 z 0 0 0 k 24 z − E A L 0 0 0 0 0 E A L 0 0 0 0 0 0 k 13 z 0 0 0 k 23 z 0 k 33 z 0 0 0 k 34 z 0 0 k 13 y 0 k 23 y 0 0 0 k 33 y 0 k 34 y 0 0 0 0 − G I t L 0 0 0 0 0 G I t L 0 0 0 0 k 14 y 0 k 24 y 0 0 0 k 34 y 0 k 44 y 0 0 k 14 z 0 0 0 k 24 z 0 k 34 z 0 0 0 k 44 z ] .
\overline{\htmlClass{fea-math-bold-italic}{\boldsymbol{k}}}_{el}
=
\begin{bmatrix}
\dfrac{EA}{L}&0&0&0&0&0&-\dfrac{EA}{L}&0&0&0&0&0\\
0&k_{11}^{z}&0&0&0&k_{12}^{z}&0&k_{13}^{z}&0&0&0&k_{14}^{z}\\
0&0&k_{11}^{y}&0&k_{12}^{y}&0&0&0&k_{13}^{y}&0&k_{14}^{y}&0\\
0&0&0&\dfrac{GI_t}{L}&0&0&0&0&0&-\dfrac{GI_t}{L}&0&0\\
0&0&k_{12}^{y}&0&k_{22}^{y}&0&0&0&k_{23}^{y}&0&k_{24}^{y}&0\\
0&k_{12}^{z}&0&0&0&k_{22}^{z}&0&k_{23}^{z}&0&0&0&k_{24}^{z}\\
-\dfrac{EA}{L}&0&0&0&0&0&\dfrac{EA}{L}&0&0&0&0&0\\
0&k_{13}^{z}&0&0&0&k_{23}^{z}&0&k_{33}^{z}&0&0&0&k_{34}^{z}\\
0&0&k_{13}^{y}&0&k_{23}^{y}&0&0&0&k_{33}^{y}&0&k_{34}^{y}&0\\
0&0&0&-\dfrac{GI_t}{L}&0&0&0&0&0&\dfrac{GI_t}{L}&0&0\\
0&0&k_{14}^{y}&0&k_{24}^{y}&0&0&0&k_{34}^{y}&0&k_{44}^{y}&0\\
0&k_{14}^{z}&0&0&0&k_{24}^{z}&0&k_{34}^{z}&0&0&0&k_{44}^{z}
\end{bmatrix}.
k e l = L E A 0 0 0 0 0 − L E A 0 0 0 0 0 0 k 11 z 0 0 0 k 12 z 0 k 13 z 0 0 0 k 14 z 0 0 k 11 y 0 k 12 y 0 0 0 k 13 y 0 k 14 y 0 0 0 0 L G I t 0 0 0 0 0 − L G I t 0 0 0 0 k 12 y 0 k 22 y 0 0 0 k 23 y 0 k 24 y 0 0 k 12 z 0 0 0 k 22 z 0 k 23 z 0 0 0 k 24 z − L E A 0 0 0 0 0 L E A 0 0 0 0 0 0 k 13 z 0 0 0 k 23 z 0 k 33 z 0 0 0 k 34 z 0 0 k 13 y 0 k 23 y 0 0 0 k 33 y 0 k 34 y 0 0 0 0 − L G I t 0 0 0 0 0 L G I t 0 0 0 0 k 14 y 0 k 24 y 0 0 0 k 34 y 0 k 44 y 0 0 k 14 z 0 0 0 k 24 z 0 k 34 z 0 0 0 k 44 z . For bending in the x ˉ y ˉ \bar x\bar y x ˉ y ˉ
[ k 11 z k 12 z k 13 z k 14 z k 22 z k 23 z k 24 z k 33 z k 34 z k 44 z ] = E I z ( 1 + Φ y ) L 3 [ 12 6 L − 12 6 L ( 4 + Φ y ) L 2 − 6 L ( 2 − Φ y ) L 2 12 − 6 L ( 4 + Φ y ) L 2 ] .
\begin{bmatrix}
k_{11}^{z}&k_{12}^{z}&k_{13}^{z}&k_{14}^{z}\\
&k_{22}^{z}&k_{23}^{z}&k_{24}^{z}\\
&&k_{33}^{z}&k_{34}^{z}\\
&&&k_{44}^{z}
\end{bmatrix}
=
\frac{EI_z}
{(1+\mathit{\Phi}_y)L^3}
\begin{bmatrix}
12&6L&-12&6L\\
&(4+\mathit{\Phi}_y)L^2&-6L&(2-\mathit{\Phi}_y)L^2\\
&&12&-6L\\
&&&(4+\mathit{\Phi}_y)L^2
\end{bmatrix}.
k 11 z k 12 z k 22 z k 13 z k 23 z k 33 z k 14 z k 24 z k 34 z k 44 z = ( 1 + Φ y ) L 3 E I z 12 6 L ( 4 + Φ y ) L 2 − 12 − 6 L 12 6 L ( 2 − Φ y ) L 2 − 6 L ( 4 + Φ y ) L 2 . For bending in the x ˉ z ˉ \bar x\bar z x ˉ z ˉ
[ k 11 y k 12 y k 13 y k 14 y k 22 y k 23 y k 24 y k 33 y k 34 y k 44 y ] = E I y ( 1 + Φ z ) L 3 [ 12 − 6 L − 12 − 6 L ( 4 + Φ z ) L 2 6 L ( 2 − Φ z ) L 2 12 6 L ( 4 + Φ z ) L 2 ] .
\begin{bmatrix}
k_{11}^{y}&k_{12}^{y}&k_{13}^{y}&k_{14}^{y}\\
&k_{22}^{y}&k_{23}^{y}&k_{24}^{y}\\
&&k_{33}^{y}&k_{34}^{y}\\
&&&k_{44}^{y}
\end{bmatrix}
=
\frac{EI_y}
{(1+\mathit{\Phi}_z)L^3}
\begin{bmatrix}
12&-6L&-12&-6L\\
&(4+\mathit{\Phi}_z)L^2&6L&(2-\mathit{\Phi}_z)L^2\\
&&12&6L\\
&&&(4+\mathit{\Phi}_z)L^2
\end{bmatrix}.
k 11 y k 12 y k 22 y k 13 y k 23 y k 33 y k 14 y k 24 y k 34 y k 44 y = ( 1 + Φ z ) L 3 E I y 12 − 6 L ( 4 + Φ z ) L 2 − 12 6 L 12 − 6 L ( 2 − Φ z ) L 2 6 L ( 4 + Φ z ) L 2 . The shear parameters are
Φ y = 12 E I z G A y L 2 , Φ z = 12 E I y G A z L 2 .
\mathit{\Phi}_y
=
\frac{12EI_z}{GA_yL^2},
\qquad
\mathit{\Phi}_z
=
\frac{12EI_y}{GA_zL^2}.
Φ y = G A y L 2 12 E I z , Φ z = G A z L 2 12 E I y . 05
Numerical integration As in Section 23.1 , we prefer to use Gauss quadrature with two integration points.
This integration is exact because the entries of B \htmlClass{fea-math-bold-italic}{\boldsymbol{B}} B x ˉ \bar x x ˉ
B T D B
\htmlClass{fea-math-bold-italic}{\boldsymbol{B}}^T\htmlClass{fea-math-bold-italic}{\boldsymbol{D}}\htmlClass{fea-math-bold-italic}{\boldsymbol{B}}
B T D B contains polynomials of at most second degree.
06
Transformation to the global reference frame The same passive rotation convention as in Section 23.1 is used:
u ‾ e l = R u e l ,
\overline{\htmlClass{fea-math-bold-italic}{\boldsymbol{u}}}_{el}
=
\htmlClass{fea-math-bold-italic}{\boldsymbol{R}}\htmlClass{fea-math-bold-italic}{\boldsymbol{u}}_{el},
u e l = R u e l , where R \htmlClass{fea-math-bold-italic}{\boldsymbol{R}} R
Consequently,
k e l = R T k ‾ e l R .
\boxed{
\htmlClass{fea-math-bold-italic}{\boldsymbol{k}}_{el}
=
\htmlClass{fea-math-bold-italic}{\boldsymbol{R}}^T
\overline{\htmlClass{fea-math-bold-italic}{\boldsymbol{k}}}_{el}
\htmlClass{fea-math-bold-italic}{\boldsymbol{R}}
}.
k e l = R T k e l R . Refer to Section 23.1 for details concerning the construction of R \htmlClass{fea-math-bold-italic}{\boldsymbol{R}} R
Within the hypotheses of the corresponding beam models, the finite elements presented in Sections 23.1 and 23.2 reproduce the Euler–Bernoulli and Timoshenko beam behaviour, respectively.
A MATLAB computer code is available for download.
MATLAB 3D beam element — Timoshenko model Download ZIP 07
Example Let us consider again the example from Section 23.1 : a clamped curved helicoidal beam is analysed,
R = 50 m m , H = 80 m m ,
R=50\ {\it mm},
\qquad
H=80\ {\it mm},
R = 50 mm , H = 80 mm , with rectangular cross-section
h × b = 9 × 6 m m .
h\times b=9\times6\ {\it mm}.
h × b = 9 × 6 mm . Refer to the MATLAB subprogram gen for all input data.
Figure 2. Figure 3. The displacements of the loaded beam end are:
Euler–Bernoulli beam model Timoshenko beam model u u u − 9.814 -9.814 − 9.814 − 9.824 -9.824 − 9.824 v v v 10.814 10.814 10.814 10.824 10.824 10.824 w w w 17.105 17.105 17.105 17.123 17.123 17.123 ψ x \psi_x ψ x − 0.074273 -0.074273 − 0.074273 − 0.074273 -0.074273 − 0.074273 ψ y \psi_y ψ y − 0.19325 -0.19325 − 0.19325 − 0.19325 -0.19325 − 0.19325 ψ z \psi_z ψ z 0.068256 0.068256 0.068256 0.068256 0.068256 0.068256
Remark. The displacements are given in the global reference frame. The angles ψ x , ψ y , ψ z \psi_x,\psi_y,\psi_z ψ x , ψ y , ψ z
As can be seen, the linear displacements are very slightly larger in the case of the Timoshenko beam model, because the shear angles β y \beta_y β y β z \beta_z β z
In both models, the angular degrees of freedom are the rotations of the cross-section. In the Euler–Bernoulli beam model, the rotation of the cross-section coincides with the rotation of the normal.
PREVIOUS SECTION ← 23.1. 3D beam element. Euler–Bernoulli model NEXT SECTION 23.3. 3D isoparametric beam element. Two-node finite element. Timoshenko model → 
