Chapter 2
Direct method

Starting with only one element beam which is subject to bending and shear forces. There are 4 nodal degrees of freedom. Rotation at the left and right nodes of the beam and transverse displacements at the left and right nodes. The following diagram shows the sign convention used for external forces. Moments are always positive when anti-clockwise direction and vertical forces are positive when in the positive \(y\) direction.

The two nodes are numbered \(1\) and \(2\) from left to right. \(M_{1}\) is the moment at the left node (node 1), \(M_{2}\) is the moment at the right node (node 2). \(V_{1}\) is the vertical force at the left node and \(V_{2}\) is the vertical force at the right node.

The above diagram shows the signs used for the applied forces direction when acting in the positive sense. Since this is a one dimensional problem, the displacement field (the unknown being solved for) will be a function of one independent variable which is the \(x\) coordinate. The displacement field in the vertical direction is called \(v\left ( x\right ) \). This is the vertical displacement of point \(x\) on the beam from the original \(x-axis\). The following diagram shows the notation used for the coordinates.

Angular displacement at distance \(x\) on the beam is found using \(\theta \left ( x\right ) =\frac {dv\left ( x\right ) }{dx}.\) At the left node, the degrees of freedom or the displacements, are called \(v_{1},\theta _{1}\) and at the right node they are called \(v_{2},\theta _{2}\). At an arbitrary location \(x\) in the beam, the vertical displacement is \(v\left ( x\right ) \) and the rotation at that location is \(\theta \left ( x\right ) \).

The following diagram shows the displacement field \(v\left ( x\right ) \)

In the direct method of finding the stiffness matrix, the forces at the ends of the beam are found directly by the use of beam theory. In beam theory the signs are different from what is given in the first diagram above. Therefore, the moment and shear forces obtained using beam theory (\(M_{B}\) and \(V_{B}\) in the diagram below) will have different signs when compared to the external forces. The signs are then adjusted to reflect the convention as shown in the diagram above using \(M\) and \(V\).

For an example, the external moment \(M_{1}\) is opposite in sign to \(M_{B1}\) and the reaction \(V_{2}\) is opposite to \(V_{B2}\). To illustrate this more, a diagram with both sign conventions is given below.

The goal now is to obtain expressions for external loads \(M_{i}\) and \(R_{i}\) in the above diagram as function of the displacements at the nodes \(\left \{ d\right \} =\{v_{1},\theta _{1},v_{2},\theta _{2}\}^{T}\).

In other words, the goal is to obtain an expression of the form \(\left \{ p\right \} =\left [ K\right ] \left \{ d\right \} \) where \(\left [ K\right ] \) is the stiffness matrix where \(\left \{ p\right \} =\left \{ V_{1},M_{1},R_{2},M_{2}\right \} ^{T}\) is the nodal forces or load vector, and \(\left \{ d\right \} \) is the nodal displacement vector.

In this case \(\left [ K\right ] \) will be a \(4\times 4\) matrix and \(\left \{ p\right \} \) is a \(4\times 1\) vector and \(\left \{ d\right \} \) is a \(4\times 1\) vector.

Starting with \(V_{1}\). It is in the same direction as the shear force \(V_{B1}\). Since \(V_{B1}=\frac {dM_{B1}}{dx}\) then

\[ V_{1}=\frac {dM_{B1}}{dx}\]

Since from beam theory \(M_{B1}=-\sigma \left ( x\right ) \frac {I}{y}\), the above becomes

\[ V_{1}=-\frac {I}{y}\frac {d\sigma \left ( x\right ) }{dx}\]

But \(\sigma \left ( x\right ) =E\varepsilon \left ( x\right ) \) and \(\varepsilon \left ( x\right ) =\frac {-y}{\rho }\) where \(\rho \) is radius of curvature, therefore the above becomes

\[ V_{1}=EI\frac {d}{dx}\left ( \frac {1}{\rho }\right ) \]

Since \(\frac {1}{\rho }=\frac {\left ( \frac {d^{2}u}{dx^{2}}\right ) }{\left ( 1+\left ( \frac {du}{dx}\right ) ^{2}\right ) ^{3/2}}\) and for a small angle of deflection \(\frac {du}{dx}\ll 1\) then \(\frac {1}{\rho }=\left ( \frac {d^{2}u}{dx^{2}}\right ) \), and the above now becomes

\[ V_{1}=EI\frac {d^{3}u\left ( x\right ) }{dx^{3}}\]

Before continuing, the following diagram illustrates the above derivation. This comes from beam theory.

Now \(M_{1}\) is obtained. \(M_{1}\) is in the opposite sense of the bending moment \(M_{B1}\) hence a negative sign is added giving \(M_{1}=-M_{B1}\). But \(M_{B1}=-\sigma \left ( x\right ) \frac {I}{y}\) therefore

\begin{align*} M_{1} & =\sigma \left ( x\right ) \frac {I}{y}\\ & =E\varepsilon \left ( x\right ) \frac {I}{y}\\ & =E\left ( \frac {-y}{\rho }\right ) \frac {I}{y}\\ & =-EI\left ( \frac {1}{\rho }\right ) \\ & =-EI\frac {d^{2}w}{dx^{2}}\end{align*}

\(V_{2}\) is now found. It is in the opposite sense of the shear force \(V_{B2}\), hence a negative sign is added giving

\begin{align*} V_{2} & =-V_{B2}\\ & =-\frac {dM_{B2}}{dx}\end{align*}

Since \(M_{B2}=-\sigma \left ( x\right ) \frac {I}{y}\), the above becomes

\[ V_{2}=\frac {I}{y}\frac {d\sigma \left ( x\right ) }{dx}\]

But \(\sigma \left ( x\right ) =E\varepsilon \left ( x\right ) \) and \(\varepsilon \left ( x\right ) =\frac {-y}{\rho }\) where \(\rho \) is radius of curvature. The above becomes

\[ V_{2}=-EI\frac {d}{dx}\left ( \frac {1}{\rho }\right ) \]

But \(\frac {1}{\rho }=\frac {\left ( \frac {d^{2}w}{dx^{2}}\right ) }{\left ( 1+\left ( \frac {dw}{dx}\right ) ^{2}\right ) ^{3/2}}\) and for small angle of deflection \(\frac {dw}{dx}\ll 1\) hence \(\frac {1}{\rho }=\left ( \frac {d^{2}u}{dx^{2}}\right ) \) then the above becomes

\[ V_{2}=-EI\frac {d^{3}u\left ( x\right ) }{dx^{3}}\]

Finally \(M_{2}\) is in the same direction as \(M_{B2}\) so no sign change is needed. \(M_{B2}=-\sigma \left ( x\right ) \frac {I}{y}\), therefore

\begin{align*} M_{2} & =-\sigma \left ( x\right ) \frac {I}{y}\\ & =-E\varepsilon \left ( x\right ) \frac {I}{y}\\ & =-E\left ( \frac {-y}{\rho }\right ) \frac {I}{y}\\ & =EI\left ( \frac {1}{\rho }\right ) \\ & =EI\frac {d^{2}u}{dx^{2}}\end{align*}

The following is a summary of what was found so far. Notice that the above expressions are evaluates at \(x=0\) and at \(x=L\). Accordingly to obtain the nodal end forces vector \(\left \{ p\right \} \)

\begin{equation} \left \{ p\right \} =\begin {Bmatrix} V_{1}\\ M_{1}\\ V_{2}\\ M_{2}\end {Bmatrix} =\begin {Bmatrix} EI\left . \frac {d^{3}u\left ( x\right ) }{dx^{3}}\right \vert _{x=0}\\ -EI\left . \frac {d^{2}u}{dx^{2}}\right \vert _{x=0}\\ -EI\left . \frac {d^{3}u\left ( x\right ) }{dx^{3}}\right \vert _{x=L}\\ EI\left . \frac {d^{2}u}{dx^{2}}\right \vert _{x=L}\end {Bmatrix} \tag {1}\end{equation}

The RHS of the above is now expressed as function of the nodal displacements \(v_{1},\theta _{1},v_{2},\theta _{2}\). To do that, the field displacement \(v\left ( x\right ) \) which is the transverse displacement of the beam is assumed to be a polynomial in \(x\) of degree \(3\) or

\begin{align} v\left ( x\right ) & =a_{0}+a_{1}x+a_{2}x^{2}+a_{3}x^{3}\nonumber \\ \theta \left ( x\right ) & =\frac {dv\left ( x\right ) }{dx}=a_{1}+2a_{2}x+3a_{3}x^{2} \tag {A}\end{align}

Polynomial of degree \(3\) is used since there are \(4\) degrees of freedom, and a minimum of \(4\) free parameters is needed. Hence

\begin{equation} v_{1}=\left . v\left ( x\right ) \right \vert _{x=0}=a_{0} \tag {2}\end{equation}

And

\begin{equation} \theta _{1}=\left . \theta \left ( x\right ) \right \vert _{x=0}=a_{1} \tag {3}\end{equation}

Assuming the length of the beam is \(L\), then

\begin{equation} v_{2}=\left . v\left ( x\right ) \right \vert _{x=L}=a_{0}+a_{1}L+a_{2}L^{2}+a_{3}L^{3} \tag {4}\end{equation}

And

\begin{equation} \theta _{2}=\left . \theta \left ( x\right ) \right \vert _{x=L}=a_{1}+2a_{2}L+3a_{3}L^{2} \tag {5}\end{equation}

Equations (2-5) gives

\[ \left \{ d\right \} =\begin {Bmatrix} v_{1}\\ \theta _{1}\\ v_{2}\\ \theta _{2}\end {Bmatrix} =\begin {Bmatrix} a_{0}\\ a_{1}\\ a_{0}+a_{1}L+a_{2}L^{2}+a_{3}L^{3}\\ a_{1}+2a_{2}L+3a_{3}L^{2}\end {Bmatrix} =\begin {pmatrix} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 1 & L & L^{2} & L^{3}\\ 0 & 1 & 2L & 3L^{2}\end {pmatrix}\begin {Bmatrix} a_{0}\\ a_{1}\\ a_{2}\\ a_{3}\end {Bmatrix} \]

Solving for \(a_{i}\) gives

\begin{align*}\begin {Bmatrix} a_{0}\\ a_{1}\\ a_{2}\\ a_{3}\end {Bmatrix} & =\begin {pmatrix} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ -\frac {3}{L^{2}} & -\frac {2}{L} & \frac {3}{L^{2}} & -\frac {1}{L}\\ \frac {2}{L^{3}} & \frac {1}{L^{2}} & -\frac {2}{L^{3}} & \frac {1}{L^{2}}\end {pmatrix}\begin {Bmatrix} v_{1}\\ \theta _{1}\\ v_{2}\\ \theta _{2}\end {Bmatrix} \\ & =\begin {pmatrix} v_{1}\\ \theta _{1}\\ \frac {3}{L^{2}}v_{2}-\frac {1}{L}\theta _{2}-\frac {3}{L^{2}}v_{1}-\frac {2}{L}\theta _{1}\\ \frac {1}{L^{2}}\theta _{1}+\frac {1}{L^{2}}\theta _{2}+\frac {2}{L^{3}}v_{1}-\frac {2}{L^{3}}v_{2}\end {pmatrix} \end{align*}

\(v\left ( x\right ) \), the field displacement function from Eq. (A) can now be written as a function of the nodal displacements

\begin{align*} v\left ( x\right ) & =a_{0}+a_{1}x+a_{2}x^{2}+a_{3}x^{3}\\ & =v_{1}+\theta _{1}x+\left ( \frac {3}{L^{2}}v_{2}-\frac {1}{L}\theta _{2}-\frac {3}{L^{2}}v_{1}-\frac {2}{L}\theta _{1}\right ) x^{2}+\left ( \frac {1}{L^{2}}\theta _{1}+\frac {1}{L^{2}}\theta _{2}+\frac {2}{L^{3}}v_{1}-\frac {2}{L^{3}}v_{2}\right ) x^{3}\\ & =v_{1}+x\theta _{1}-2\frac {x^{2}}{L}\theta _{1}+\frac {x^{3}}{L^{2}}\theta _{1}-\frac {x^{2}}{L}\theta _{2}+\frac {x^{3}}{L^{2}}\theta _{2}-3\frac {x^{2}}{L^{2}}v_{1}+2\frac {x^{3}}{L^{3}}v_{1}+3\frac {x^{2}}{L^{2}}v_{2}-2\frac {x^{3}}{L^{3}}v_{2}\end{align*}

Or in matrix form

\begin{align*} v\left ( x\right ) & =\begin {pmatrix} 1-3\frac {x^{2}}{L^{2}}+2\frac {x^{3}}{L^{3}} & \ x-2\frac {x^{2}}{L}+\frac {x^{3}}{L^{2}} & \ 3\frac {x^{2}}{L^{2}}-2\frac {x^{3}}{L^{3}} & \ -\frac {x^{2}}{L}+\frac {x^{3}}{L^{2}}\end {pmatrix}\begin {Bmatrix} v_{1}\\ \theta _{1}\\ v_{2}\\ \theta _{2}\end {Bmatrix} \\ & =\begin {pmatrix} \frac {1}{L^{3}}\left ( L^{3}-3Lx^{2}+2x^{3}\right ) \ & \frac {1}{L^{2}}\left ( L^{2}x-2Lx^{2}+x^{3}\right ) & \ \frac {1}{L^{3}}\left ( 3Lx^{2}-2x^{3}\right ) & \ \frac {1}{L^{2}}\left ( -Lx^{2}+x^{3}\right ) \end {pmatrix}\begin {Bmatrix} v_{1}\\ \theta _{1}\\ v_{2}\\ \theta _{2}\end {Bmatrix} \end{align*}

The above can be written as

\begin{align} v\left ( x\right ) & =\begin {pmatrix} N_{1}\left ( x\right ) \ & \ N_{2}\left ( x\right ) & \ N_{3}\left ( x\right ) & \ N_{4}\left ( x\right ) \end {pmatrix}\begin {Bmatrix} v_{1}\\ \theta _{1}\\ v_{2}\\ \theta _{2}\end {Bmatrix} \tag {5A}\\ v\left ( x\right ) & =\left [ N\right ] \left \{ d\right \} \nonumber \end{align}

Where \(N_{i}\) are called the shape functions. The shape functions are

\[\begin {pmatrix} N_{1}\left ( x\right ) \\ N_{2}\left ( x\right ) \\ \ N_{3}\left ( x\right ) \\ \ N_{4}\left ( x\right ) \end {pmatrix} =\begin {pmatrix} \frac {1}{L^{3}}\left ( L^{3}-3Lx^{2}+2x^{3}\right ) \\ \frac {1}{L^{2}}\left ( L^{2}x-2Lx^{2}+x^{3}\right ) \\ \frac {1}{L^{3}}\left ( 3Lx^{2}-2x^{3}\right ) \\ \frac {1}{L^{2}}\left ( -Lx^{2}+x^{3}\right ) \end {pmatrix} \]

We notice that \(N_{1}\left ( 0\right ) =1\) and \(N_{1}\left ( L\right ) =0\) as expected. Also

\[ \left . \frac {dN_{2}\left ( x\right ) }{dx}\right \vert _{x=0}=\left . \frac {1}{L^{2}}\left ( L^{2}-4Lx+3x^{2}\right ) \right \vert _{x=0}=1 \]

And

\[ \left . \frac {dN_{2}\left ( x\right ) }{dx}\right \vert _{x=L}=\left . \frac {1}{L^{2}}\left ( L^{2}-4Lx+3x^{2}\right ) \right \vert _{x=L}=0 \]

Also \(N_{3}\left ( 0\right ) =0\) and \(N_{3}\left ( L\right ) =1\) and

\[ \left . \frac {dN_{4}\left ( x\right ) }{dx}\right \vert _{x=0}=\left . \frac {1}{L^{2}}\left ( -2Lx+3x^{2}\right ) \right \vert _{x=0}=0 \]

and

\[ \left . \frac {dN_{4}\left ( x\right ) }{dx}\right \vert _{x=L}=\left . \frac {1}{L^{2}}\left ( -2Lx+3x^{2}\right ) \right \vert _{x=L}=1 \]

The shape functions have thus been verified. The stiffness matrix is now found by substituting Eq. (5A) into Eq. (1), repeated below

\[ \left \{ p\right \} =\begin {Bmatrix} V_{1}\\ M_{1}\\ V_{2}\\ M_{2}\end {Bmatrix} =\begin {pmatrix} EI\left . \frac {d^{3}v\left ( x\right ) }{dx^{3}}\right \vert _{x=0}\\ -EI\left . \frac {d^{2}v}{dx^{2}}\right \vert _{x=0}\\ -EI\left . \frac {d^{3}v\left ( x\right ) }{dx^{3}}\right \vert _{x=L}\\ EI\left . \frac {d^{2}v}{dx^{2}}\right \vert _{x=L}\end {pmatrix} \]

Hence

\begin{equation} \left \{ p\right \} =\begin {Bmatrix} V_{1}\\ M_{1}\\ V_{2}\\ M_{2}\end {Bmatrix} =\begin {pmatrix} EI\frac {d^{3}}{dx^{3}}\left ( N_{1}v_{1}+N_{2}\theta _{1}+N_{3}v_{2}+N_{4}\theta _{2}\right ) \\ -EI\frac {d^{2}}{dx^{2}}\left ( N_{1}v_{1}+N_{2}\theta _{1}+N_{3}v_{2}+N_{4}\theta _{2}\right ) \\ -EI\frac {d^{3}}{dx^{3}}\left ( N_{1}v_{1}+N_{2}\theta _{1}+N_{3}v_{2}+N_{4}\theta _{2}\right ) \\ EI\frac {d^{2}}{dx^{2}}\left ( N_{1}v_{1}+N_{2}\theta _{1}+N_{3}v_{2}+N_{4}\theta _{2}\right ) \end {pmatrix} \tag {6}\end{equation}

But

\begin{align*} \frac {d^{3}}{dx^{3}}N_{1}\left ( x\right ) & =\frac {1}{L^{3}}\frac {d^{3}}{dx^{3}}\left ( L^{3}-3Lx^{2}+2x^{3}\right ) \\ & =\frac {12}{L^{3}}\end{align*}

And

\begin{align*} \frac {d^{3}}{dx^{3}}N_{2}\left ( x\right ) & =\frac {1}{L^{2}}\frac {d^{3}}{dx^{3}}\left ( L^{2}x-2Lx^{2}+x^{3}\right ) \\ & =\frac {6}{L^{2}}\end{align*}

And

\begin{align*} \frac {d^{3}}{dx^{3}}N_{3}\left ( x\right ) & =\frac {1}{L^{3}}\frac {d^{3}}{dx^{3}}\ \left ( 3Lx^{2}-2x^{3}\right ) \\ & =\frac {-12}{L^{3}}\end{align*}

And

\begin{align*} \frac {d^{3}}{dx^{3}}N_{4}\left ( x\right ) & =\frac {1}{L^{2}}\frac {d^{3}}{dx^{3}}\left ( -Lx^{2}+x^{3}\right ) \\ & =\frac {6}{L^{2}}\end{align*}

For the second derivatives

\begin{align*} \frac {d^{2}}{dx^{2}}N_{1}\left ( x\right ) & =\frac {1}{L^{3}}\frac {d^{2}}{dx^{2}}\left ( L^{3}-3Lx^{2}+2x^{3}\right ) \\ & =\frac {1}{L^{3}}\left ( 12x-6L\right ) \end{align*}

And

\begin{align*} \frac {d^{2}}{dx^{2}}N_{2}\left ( x\right ) & =\frac {1}{L^{2}}\frac {d^{2}}{dx^{2}}\left ( L^{2}x-2Lx^{2}+x^{3}\right ) \\ & =\frac {1}{L^{2}}\left ( 6x-4L\right ) \end{align*}

And

\begin{align*} \frac {d^{2}}{dx^{2}}N_{3}\left ( x\right ) & =\frac {1}{L^{3}}\frac {d^{2}}{dx^{2}}\ \left ( 3Lx^{2}-2x^{3}\right ) \\ & =\frac {1}{L^{3}}\left ( 6L-12x\right ) \end{align*}

And

\begin{align*} \frac {d^{2}}{dx^{2}}N_{4}\left ( x\right ) & =\frac {1}{L^{2}}\frac {d^{2}}{dx^{2}}\left ( -Lx^{2}+x^{3}\right ) \\ & =\frac {1}{L^{2}}\left ( 6x-2L\right ) \end{align*}

Hence Eq. (6) becomes

\begin{align} \left \{ p\right \} & =\begin {Bmatrix} V_{1}\\ M_{1}\\ V_{2}\\ M_{2}\end {Bmatrix} =\begin {pmatrix} EI\left . \frac {d^{3}}{dx^{3}}\left ( N_{1}v_{1}+N_{2}\theta _{1}+N_{3}v_{2}+N_{4}\theta _{2}\right ) \right \vert _{x=0}\\ -EI\left . \frac {d^{2}}{dx^{2}}\left ( N_{1}v_{1}+N_{2}\theta _{1}+N_{3}v_{2}+N_{4}\theta _{2}\right ) \right \vert _{x=0}\\ -EI\left . \frac {d^{3}}{dx^{3}}\left ( N_{1}v_{1}+N_{2}\theta _{1}+N_{3}v_{2}+N_{4}\theta _{2}\right ) \right \vert _{x=L}\\ EI\left . \frac {d^{2}}{dx^{2}}\left ( N_{1}v_{1}+N_{2}\theta _{1}+N_{3}v_{2}+N_{4}\theta _{2}\right ) \right \vert _{x=L}\end {pmatrix} \nonumber \\ & =\begin {pmatrix} EI\left ( \frac {12}{L^{3}}v_{1}+\frac {6}{L^{2}}\theta _{1}-\frac {12}{L^{3}}v_{2}+\frac {6}{L^{2}}\theta _{2}\right ) _{x=0}\\ -EI\left ( \frac {1}{L^{3}}\left ( 12x-6L\right ) v_{1}+\frac {1}{L^{2}}\left ( 6x-4L\right ) \theta _{1}+\frac {1}{L^{3}}\left ( 6L-12x\right ) v_{2}+\frac {1}{L^{2}}\left ( 6x-2L\right ) \theta _{2}\right ) _{x=0}\\ -EI\left ( \frac {12}{L^{3}}v_{1}+\frac {6}{L^{2}}\theta _{1}-\frac {12}{L^{3}}v_{2}+\frac {6}{L^{2}}\theta _{2}\right ) _{x=L}\\ EI\left ( \frac {1}{L^{3}}\left ( 12x-6L\right ) v_{1}+\frac {1}{L^{2}}\left ( 6x-4L\right ) \theta _{1}+\frac {1}{L^{3}}\left ( 6L-12x\right ) v_{2}+\frac {1}{L^{2}}\left ( 6x-2L\right ) \theta _{2}\right ) _{x=L}\end {pmatrix} \nonumber \\ & =\frac {EI}{L^{3}}\begin {pmatrix} 12 & 6L & -12 & 6L\\ -\left ( 12x-6L\right ) _{x=0} & -L\left ( 6x-4L\right ) _{x=0} & -\left ( 6L-12x\right ) _{x=0} & -L\left ( 6x-2L\right ) _{x=0}\\ -12 & -6L & 12 & -6L\\ \left ( 12x-6L\right ) _{x=L} & L\left ( 6x-4L\right ) _{x=L} & \left ( 6L-12x\right ) _{x=L} & L\left ( 6x-2L\right ) _{x=L}\end {pmatrix}\begin {Bmatrix} v_{1}\\ \theta _{1}\\ v_{2}\\ \theta _{2}\end {Bmatrix} \tag {7}\end{align}

Or in matrix form, after evaluating the expressions above for \(x=L\) and \(x=0\) as

\begin{align*}\begin {Bmatrix} V_{1}\\ M_{1}\\ V_{2}\\ M_{2}\end {Bmatrix} & =\frac {EI}{L^{3}}\begin {pmatrix} 12 & 6L & -12 & 6L\\ 6L & 4L^{2} & -6L & 2L^{2}\\ -\frac {12}{L^{3}} & -\frac {6}{L^{2}} & \frac {12}{L^{3}} & -\frac {6}{L^{2}}\\ \left ( 12L-6L\right ) & L\left ( 6L-4L\right ) & \left ( 6L-12L\right ) & L\left ( 6L-2L\right ) \end {pmatrix}\begin {Bmatrix} v_{1}\\ \theta _{1}\\ v_{2}\\ \theta _{2}\end {Bmatrix} \\ & =\frac {EI}{L^{3}}\begin {pmatrix} 12 & 6L & -12 & 6L\\ 6L & 4L^{2} & -6L & 2L^{2}\\ -12 & -6L & 12 & -6L\\ 6L & 2L^{2} & -6L & 4L^{2}\end {pmatrix}\begin {Bmatrix} v_{1}\\ \theta _{1}\\ v_{2}\\ \theta _{2}\end {Bmatrix} \end{align*}

The above now is in the form

\[ \left \{ p\right \} =\left [ K\right ] \left \{ d\right \} \]

Hence the stiffness matrix is

\[ \left [ K\right ] =\frac {EI}{L^{3}}\begin {pmatrix} 12 & 6L & -12 & 6L\\ 6L & 4L^{2} & -6L & 2L^{2}\\ -12 & -6L & 12 & -6L\\ 6L & 2L^{2} & -6L & 4L^{2}\end {pmatrix} \]

Knowing the stiffness matrix means knowing the nodal displacements \(\left \{ d\right \} \) when given the forces at the nodes. The power of the finite element method now comes after all the nodal displacements \(v_{1},\theta _{1},\) \(v_{2},\theta _{2}\) are calculated by solving \(\left \{ p\right \} =\left [ K\right ] \left \{ d\right \} \). This is because the polynomial \(v\left ( x\right ) \) is now completely determined and hence \(v\left ( x\right ) \) and \(\theta \left ( x\right ) \) can now be evaluated for any \(x\) along the beam and not just at its end nodes as the case with finite difference method. Eq. 5A above can now be used to find the displacement \(v\left ( x\right ) \) and \(\theta \left ( x\right ) \) everywhere.

\begin{align*} v\left ( x\right ) & =\left [ N\right ] \left \{ d\right \} \\ v\left ( x\right ) & =\begin {pmatrix} N_{1}\left ( x\right ) \ & \ N_{2}\left ( x\right ) & \ N_{3}\left ( x\right ) & \ N_{4}\left ( x\right ) \end {pmatrix}\begin {Bmatrix} v_{1}\\ \theta _{1}\\ v_{2}\\ \theta _{2}\end {Bmatrix} \end{align*}

To summarise, these are the steps to obtain \(v\left ( x\right ) \)

1. An expression for \(\left [ K\right ] \) is found.
2. \(\left \{ p\right \} =\left [ K\right ] \left \{ d\right \} \) is solved for \(\left \{ d\right \} \)
3. \(v\left ( x\right ) =\left [ N\right ] \left \{ d\right \} \) is calculated by assuming \(v\left ( x\right ) \) is a polynomial. This gives the displacement \(v\left ( x\right ) \) to use to evaluate the transverse displacement anywhere on the beam and not just at the end nodes.
4. \(\theta \left ( x\right ) =\frac {dv\left ( x\right ) }{dx}=\frac {d}{dx}\left [ N\right ] \left \{ d\right \} \) is obtained to evaluate the rotation of the beam any where and not just at the end nodes.
5. The strain \(\epsilon \left ( x\right ) =-y\left [ B\right ] \left \{ d\right \} \) is found, where \(\left [ B\right ] \) is the gradient matrix \(\left [ B\right ] =\frac {d^{2}}{dx^{2}}\left [ N\right ] \).
6. The stress from \(\sigma =E\epsilon =-Ey\left [ B\right ] \left \{ d\right \} \) is found.
7. The bending moment diagram from \(M\left ( x\right ) =EI\left [ B\right ] \left \{ d\right \} \) is found.
8. The shear force diagram from \(V\left ( x\right ) =\frac {d}{dx}M\left ( x\right ) \) is found.

2.1 Examples using the direct beam stiffness matrix

The beam stiffness matrix is now used to solve few beam problems. Starting with simple one span beam

2.1.1 Example 1

A one span beam, a cantilever beam of length \(L\), with point load \(P\) at the free end

The first step is to make the free body diagram and show all moments and forces at the nodes

\(P\) is the given force. \(M_{2}=0\) since there is no external moment at the right end. Hence \(\left \{ p\right \} =\left [ K\right ] \left \{ d\right \} \) for this system is

\[\begin {Bmatrix} R\\ M_{1}\\ -P\\ 0 \end {Bmatrix} =\frac {EI}{L^{3}}\begin {pmatrix} 12 & 6L & -12 & 6L\\ 6L & 4L^{2} & -6L & 2L^{2}\\ -12 & -6L & 12 & -6L\\ 6L & 2L^{2} & -6L & 4L^{2}\end {pmatrix}\begin {Bmatrix} v_{1}\\ \theta _{1}\\ v_{2}\\ \theta _{2}\end {Bmatrix} \]

Now is an important step. The known end displacements from boundary conditions is substituted into \(\left \{ d\right \} \), and the corresponding row and columns from the above system of equations are removed¹. Boundary conditions indicates that there is no rotation on the left end (since fixed). Hence \(v_{1}=0\) and \(\theta _{1}=0\). Hence the only unknowns are \(v_{2}\) and \(\theta _{2}\). Therefore the first and the second rows and columns are removed, giving

\[\begin {Bmatrix} -P\\ 0 \end {Bmatrix} =\frac {EI}{L^{3}}\begin {pmatrix} 12 & -6L\\ -6L & 4L^{2}\end {pmatrix}\begin {Bmatrix} v_{2}\\ \theta _{2}\end {Bmatrix} \]

Now the above is solved for \(\begin {Bmatrix} v_{2}\\ \theta _{2}\end {Bmatrix} \). Let \(E=30\times 10^{6}\) psi (steel), \(I=57\) \(in^{4}\), \(L=144\ in\), and \(P=400\) lb, hence

P=400; 
L=144; 
E=30*10^6; 
I=57.1; 
A=(E*I/L^3)*[12       6*L            -12      6*L; 
6*L      4*L^2          -6*L     2*L^2; 
-12     -6*L            12      -6*L; 
6*L      2*L^2          -6*L     4*L^2 
]; 
load=[P;P*L;-P;0]; 
x=A(3:end,3:end)\load(3:end)
 

which gives

x = 
-0.2324 
-0.0024
 

Therefore the vertical displacement at the right end is \(v_{2}=0.2324\) inches (downwards) and \(\theta _{2}=-0.0024\) radians. Now that all nodal displacements are found, the field displacement function is completely determined.

\[ \left \{ d\right \} =\begin {Bmatrix} v_{1}\\ \theta _{1}\\ v_{2}\\ \theta _{2}\end {Bmatrix} =\begin {Bmatrix} 0\\ 0\\ -0.2324\\ -0.0024 \end {Bmatrix} \]

From Eq. 5A

\begin{align*} v\left ( x\right ) & =\left [ N\right ] \left \{ d\right \} \\ & =\begin {pmatrix} N_{1}\left ( x\right ) \ & \ N_{2}\left ( x\right ) & \ N_{3}\left ( x\right ) & \ N_{4}\left ( x\right ) \end {pmatrix}\begin {Bmatrix} 0\\ 0\\ -0.2324\\ -0.0024 \end {Bmatrix} \\ & =\begin {pmatrix} \frac {1}{L^{3}}\left ( L^{3}-3Lx^{2}+2x^{3}\right ) \ & \ \frac {1}{L^{2}}\left ( L^{2}x-2Lx^{2}+x^{3}\right ) & \ \frac {1}{L^{3}}\left ( 3Lx^{2}-2x^{3}\right ) & \ \frac {1}{L^{2}}\left ( -Lx^{2}+x^{3}\right ) \end {pmatrix}\begin {Bmatrix} 0\\ 0\\ -0.2324\\ -0.0024 \end {Bmatrix} \\ & =\frac {0.002\,4}{L^{2}}\left ( Lx^{2}-x^{3}\right ) +\frac {0.232\,4}{L^{3}}\left ( 2x^{3}-3Lx^{2}\right ) \end{align*}

But \(L=144\) inches, and the above becomes

\[ v\left ( x\right ) =3.992\,0\times 10^{-8}x^{3}-1.695\,6\times 10^{-5}x^{2}\]

To verify, let \(x=144\) in the above

\begin{align*} v\left ( x=144\right ) & =3.992\,0\times 10^{-8}144^{3}-1.695\,6\times 10^{-5}144^{2}\\ & =-0.232\,40 \end{align*}

The following is a plot of the deflection curve for the beam

v=@(x) 3.992*10^-8*x.^3-1.6956*10^-5*x.^2 
x=0:0.1:144; 
plot(x,v(x),'r-','LineWidth',2); 
ylim([-0.8 0.3]); 
title('beam deflection curve'); 
xlabel('x inch'); ylabel('deflection inch'); 
grid
 

Now instead of removing rows/columns for the known boundary conditions, a \(1\) is put on the diagonal. Starting again with

\[\begin {Bmatrix} R\\ M_{1}\\ -P\\ 0 \end {Bmatrix} =\frac {EI}{L^{3}}\begin {pmatrix} 12 & 6L & -12 & 6L\\ 6L & 4L^{2} & -6L & 2L^{2}\\ -12 & -6L & 12 & -6L\\ 6L & 2L^{2} & -6L & 4L^{2}\end {pmatrix}\begin {Bmatrix} v_{1}\\ \theta _{1}\\ v_{2}\\ \theta _{2}\end {Bmatrix} \]

Since \(v_{1}=0\) and \(\theta _{1}=0\), then

\[\begin {Bmatrix} 0\\ 0\\ -P\\ 0 \end {Bmatrix} =\frac {EI}{L^{3}}\begin {pmatrix} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 12 & -6L\\ 0 & 0 & -6L & 4L^{2}\end {pmatrix}\begin {Bmatrix} 0\\ 0\\ v_{2}\\ \theta _{2}\end {Bmatrix} \]

The above system is now solved as before. \(E=30\times 10^{6}\) psi, \(I=57\) \(in^{4}\), \(L=144\ \)in, \(P=400\) lb

P=400; 
L=144; 
E=30*10^6; 
I=57.1; 
A=(E*I/L^3)*[12       6*L            -12      6*L; 
6*L      4*L^2          -6*L     2*L^2; 
-12     -6*L            12      -6*L; 
6*L      2*L^2          -6*L     4*L^2 
]; 
load=[0;0;-P;0];  %put zeros for known B.C. 
A(:,1)=0; A(1,:)=0; A(1,1)=1;  %put 1 on diagonal 
A(:,2)=0; A(2,:)=0; A(2,2)=1;  %put 1 on diagonal 
A
 

Gives

A = 
1.0e+07 * 
0.0000         0         0         0 
0    0.0000         0         0 
0         0    0.0007   -0.0496 
0         0   -0.0496    4.7583
 

Then

sol=A\load  %SOLVE
 

Gives

sol = 
0 
0 
-0.2324 
-0.0024
 

The same solution is obtained as before, but without the need to remove rows/column from the stiffness matrix. This method might be easier for programming than the first method of removing rows/columns.

The rest now is the same as was done earlier and will not be repeated.

2.1.2 Example 2

This is the same example as above, but the vertical load \(P\) is now placed in the middle of the beam

In using stiffness method, all loads must be on the nodes. The vector \(\left \{ p\right \} \) is the nodal forces vector. Hence equivalent nodal loads are found for the load in the middle of the beam. The equivalent loading is the following

Therefore, the problem is as if it was the following problem

Now that equivalent loading is in place, we continue as before. Making a free body diagram showing all loads (including reaction forces)

The stiffness equation is now written as

\begin{align*} \left \{ p\right \} & =\left [ K\right ] \left \{ d\right \} \\\begin {Bmatrix} R-P/2\\ M_{1}-PL/8\\ -P/2\\ PL/8 \end {Bmatrix} & =\frac {EI}{L^{3}}\begin {pmatrix} 12 & 6L & -12 & 6L\\ 6L & 4L^{2} & -6L & 2L^{2}\\ -12 & -6L & 12 & -6L\\ 6L & 2L^{2} & -6L & 4L^{2}\end {pmatrix}\begin {Bmatrix} v_{1}\\ \theta _{1}\\ v_{2}\\ \theta _{2}\end {Bmatrix} \end{align*}

There is no need to determine \(R\) and \(M_{1}\) at this point since these rows will be removed due to boundary conditions \(v_{1}=0\) and \(\theta _{1}=0\) and hence those quantities are not needed to solve the equations. Note that the rows and columns are removed for the known boundary displacements before solving \(\left \{ p\right \} =\left [ K\right ] \left \{ d\right \} \). Hence, after removing the first two rows and columns, the above system simplifies to

\[\begin {Bmatrix} -P/2\\ PL/8 \end {Bmatrix} =\frac {EI}{L^{3}}\begin {pmatrix} 12 & -6L\\ -6L & 4L^{2}\end {pmatrix}\begin {Bmatrix} v_{2}\\ \theta _{2}\end {Bmatrix} \]

The above is now solved for \(\begin {Bmatrix} v_{2}\\ \theta _{2}\end {Bmatrix} \) using the same numerical values for \(P,E,I,L\) as in the first example

P=400; 
L=144; 
E=30*10^6; 
I=57.1; 
A=(E*I/L^3)*[12       6*L            -12      6*L; 
6*L      4*L^2          -6*L     2*L^2; 
-12     -6*L            12      -6*L; 
6*L      2*L^2          -6*L     4*L^2 
]; 
load=[-P/2;P*L/8]; 
x=A(3:end,3:end)\load
 

Gives

x = 
-0.072630472854641 
-0.000605253940455
 

Therefore

\[ \left \{ d\right \} =\begin {Bmatrix} v_{1}\\ \theta _{1}\\ v_{2}\\ \theta _{2}\end {Bmatrix} =\begin {Bmatrix} 0\\ 0\\ -0.072630472854641\\ -0.000605253940455 \end {Bmatrix} \]

This is enough to obtain \(v\left ( x\right ) \) as before. Now the reactions \(R\) and \(M_{1}\) can be determined if needed. Going back to the full \(\left \{ p\right \} =\left [ K\right ] \left \{ d\right \} \), results in

\begin{align*}\begin {Bmatrix} R-P/2\\ M_{1}-PL/8\\ -P/2\\ PL/8 \end {Bmatrix} & =\frac {EI}{L^{3}}\begin {pmatrix} 12 & 6L & -12 & 6L\\ 6L & 4L^{2} & -6L & 2L^{2}\\ -12 & -6L & 12 & -6L\\ 6L & 2L^{2} & -6L & 4L^{2}\end {pmatrix}\begin {Bmatrix} 0\\ 0\\ -0.072630472854641\\ -0.000605253940455 \end {Bmatrix} \\ & =\frac {EI}{L^{3}}\begin {pmatrix} 0.871\,57-3.631\,5\times 10^{-3}L\\ 0.435\,78L-1.210\,5\times 10^{-3}L^{2}\\ 3.631\,5\times 10^{-3}L-0.871\,57\\ 0.435\,78L-2.421\times 10^{-3}L^{2}\end {pmatrix} \end{align*}

Hence the first equation becomes

\[ R-P/2=\frac {EI}{L^{3}}\left ( 0.871\,57-3.631\,5\times 10^{-3}L\right ) \]

and since \(E=30\times 10^{6}\) psi (steel) and \(I=57\) \(in^{4}\)and \(L=144\ \)in and \(P=400\) lb, then

\(R=\frac {\left ( 30\times 10^{6}\right ) 57}{144^{3}}\left ( 0.871\,57-3.631\,5\times 10^{-3}\,\left ( 144\right ) \right ) +400/2\,\). Therefore

\[ R=400\ \text {lb}\]

and \(M_{1}-PL/8=\frac {EI}{L^{3}}\left ( 0.435\,78L-1.210\,5\times 10^{-3}L^{2}\right ) +PL/8\,,\) hence

\[ M_{1}=28\,762\text { lb-ft}\]

Now that all nodal reactions are found, the displacement field is found and the deflection curve can be plotted.

\begin{align*} v\left ( x\right ) & =\left [ N\right ] \left \{ d\right \} \\ v\left ( x\right ) & =\begin {pmatrix} N_{1}\left ( x\right ) \ & \ N_{2}\left ( x\right ) & \ N_{3}\left ( x\right ) & \ N_{4}\left ( x\right ) \end {pmatrix}\begin {Bmatrix} v_{1}\\ \theta _{1}\\ v_{2}\\ \theta _{2}\end {Bmatrix} \\ & =\begin {pmatrix} N_{1}\left ( x\right ) \ & \ N_{2}\left ( x\right ) & \ N_{3}\left ( x\right ) & \ N_{4}\left ( x\right ) \end {pmatrix}\begin {Bmatrix} 0\\ 0\\ -0.072630472854641\\ -0.000605253940455 \end {Bmatrix} \\ & =\begin {pmatrix} \ \frac {1}{L^{3}}\left ( 3Lx^{2}-2x^{3}\right ) & \ \frac {1}{L^{2}}\left ( -Lx^{2}+x^{3}\right ) \end {pmatrix}\begin {pmatrix} -0.072630472854641\\ -0.000605253940455 \end {pmatrix} \\ & =\frac {6.052\,5\times 10^{-4}}{L^{2}}\left ( Lx^{2}-x^{3}\right ) +\frac {0.072\,63}{L^{3}}\left ( 2x^{3}-3Lx^{2}\right ) \end{align*}

Since \(L=144\) inches, the above becomes

\begin{align*} v\left ( x\right ) & =\frac {6.052\,5\times 10^{-4}}{\left ( 144\right ) ^{2}}\left ( \left ( 144\right ) x^{2}-x^{3}\right ) +\frac {0.072\,63}{\left ( 144\right ) ^{3}}\left ( 2x^{3}-3\left ( 144\right ) x^{2}\right ) \\ & =1.945\,9\times 10^{-8}x^{3}-6.304\,7\times 10^{-6}x^{2}\end{align*}

The following is the plot

clear all; close all; 
v=@(x) 1.9459*10^-8*x.^3-6.3047*10^-6*x.^2 
x=0:0.1:144; 
plot(x,v(x),'r-','LineWidth',2); 
ylim([-0.8 0.3]); 
title('beam deflection curve'); 
xlabel('x inch'); ylabel('deflection inch'); 
grid
 

2.1.3 Example 3

Assuming the beam is fixed on the left end as above, but simply supported on the right end, and the vertical load \(P\) now at distance \(a\) from the left end and at distance \(b\) from the right end, and a uniform distributed load of density \(m\) lb/in is on the beam.

Using the following values: \(a=0.625L,b=0.375L,E=30\times 10^{6}\) psi (steel), \(I=57\) \(in^{4}\), \(L=144\ \)in\(,~P=1000\) lb, \(m=200\) lb/in.

In the above, the left end reaction forces are shown as \(R_{1}\) and moment reaction as \(M_{1}\) and the reaction at the right end as \(R_{2}\). Starting by finding equivalent loads for the point load \(P\) and equivalent loads for for the uniform distributed load \(m\). All external loads must be transferred to the nodes for the stiffness method to work. Equivalent load for the above point load is

Equivalent load for the uniform distributed loading is

Using free body diagram, with all the loads on it gives the following diagram (In this diagram \(M\) is the reaction moment and \(R_{1},R_{2}\) are the reaction forces)

Now that all loads are on the nodes, the stiffness equation is applied

\begin{align*} \left \{ p\right \} & =\left [ K\right ] \left \{ d\right \} \\\begin {Bmatrix} R_{1}-\frac {Pb^{2}\left ( L+2a\right ) }{L^{3}}-\frac {mL}{2}\\ M-\frac {Pab^{2}}{L^{2}}-\frac {mL^{2}}{12}\\ R_{2}-\frac {Pa^{2}\left ( L+2b\right ) }{L^{3}}-\frac {mL}{2}\\ \frac {Pa^{2}b}{L^{2}}+\frac {mL^{2}}{12}\end {Bmatrix} & =\frac {EI}{L^{3}}\begin {pmatrix} 12 & 6L & -12 & 6L\\ 6L & 4L^{2} & -6L & 2L^{2}\\ -12 & -6L & 12 & -6L\\ 6L & 2L^{2} & -6L & 4L^{2}\end {pmatrix}\begin {Bmatrix} v_{1}\\ \theta _{1}\\ v_{2}\\ \theta _{2}\end {Bmatrix} \end{align*}

Boundary conditions are now applied. \(v_{1}=0,\theta _{1}=0\),\(v_{2}=0\). Therefore the first, second and third rows/columns are removed giving

\[ \frac {Pa^{2}b}{L^{2}}+\frac {mL^{2}}{12}=\frac {EI}{L^{3}}4L^{2}\theta _{2}\]

Hence

\[ \theta _{2}=\left ( \frac {Pa^{2}b}{L^{2}}+\frac {mL^{2}}{12}\right ) \left ( \frac {L}{4EI}\right ) \]

Substituting numerical values for the above as given at the top of the problem results in

\begin{align*} \theta _{2} & =\left ( \frac {\left ( 1000\right ) \left ( 0.625\left ( 144\right ) \right ) ^{2}\left ( 0.375\left ( 144\right ) \right ) }{\left ( 144\right ) ^{2}}+\frac {\left ( 200\right ) \left ( 144\right ) ^{2}}{12}\right ) \left ( \frac {144}{4\left ( 30\times 10^{6}\right ) \left ( 57\right ) }\right ) \\ & =7.719\,9\times 10^{-3}rad \end{align*}

Hence the field displacement \(u\left ( x\right ) \) is now found

\begin{align*} v\left ( x\right ) & =\left [ N\right ] \left \{ d\right \} \\ v\left ( x\right ) & =\begin {pmatrix} N_{1}\left ( x\right ) \ & \ N_{2}\left ( x\right ) & \ N_{3}\left ( x\right ) & \ N_{4}\left ( x\right ) \end {pmatrix}\begin {Bmatrix} v_{1}\\ \theta _{1}\\ v_{2}\\ \theta _{2}\end {Bmatrix} \\ & =\begin {pmatrix} N_{1}\left ( x\right ) \ & \ N_{2}\left ( x\right ) & \ N_{3}\left ( x\right ) & \ N_{4}\left ( x\right ) \end {pmatrix}\begin {Bmatrix} 0\\ 0\\ 0\\ 7.719\,9\times 10^{-3}\end {Bmatrix} \\ & =\ \frac {1}{L^{2}}\left ( -Lx^{2}+x^{3}\right ) \left ( 7.719\,9\times 10^{-3}\right ) \\ & =3.722\,9\times 10^{-7}x^{3}-5.361\times 10^{-5}\allowbreak x^{2}\end{align*}

And a plot of the deflection curve is

clear all; close all; 
v=@(x) 3.7229*10^-7*x.^3-5.361*10^-5*x.^2 
x=0:0.1:144; 
plot(x,v(x),'r-','LineWidth',2); 
ylim([-0.5 0.3]); 
title('beam deflection curve'); 
xlabel('x inch'); ylabel('deflection inch'); 
grid
 

¹Instead of removing rows/columns for known boundary conditions, we can also just put a \(1\) on the diagonal of the stiffness matrix for that boundary conditions. I will do this example again using this method