home

PDF (letter size)

Derivation Of The Beam Stiffness Matrix


July 7, 2016   Compiled on January 30, 2024 at 2:57am
1 Introduction
2 Direct method
2.1 Examples using the direct beam stiffness matrix
2.1.1 Example 1
2.1.2 Example 2
2.1.3 Example 3
3 Finite elements (adding more elements)
3.1 Example 3 redone with 2 elements
4 Generating shear and bending moments diagrams
5 Finding the stiffness matrix using methods other than direct method
5.1 Virtual work method for derivation of the stiffness matrix
5.2 Potential energy (minimize a functional) method to derive the stiffness matrix
6 References

Chapter 1
Introduction

A short review for solving the beam problem in 2D is given. The deflection curve, bending moment and shear force diagrams are calculated for a beam subject to bending moment and shear force using direct stiffness method and then using finite elements method by adding more elements. The problem is solved first by finding the stiffness matrix using the direct method and then using the virtual work method.

Chapter 2
Direct method

2.1 Examples using the direct beam stiffness matrix
2.1.1 Example 1
2.1.2 Example 2
2.1.3 Example 3

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

2.1.1 Example 1
2.1.2 Example 2
2.1.3 Example 3

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 removed1. 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
 

Chapter 3
Finite elements (adding more elements)

3.1 Example 3 redone with 2 elements

Finite elements generated displacements are smaller in value than the actual analytical values. To improve the accuracy, more elements are added. To add more elements, the beam is divided into 2,3,4 and more beam elements. To show how this works, example 3 above is solved again using two elements. It is found that displacement field \(v\left ( x\right ) \) becomes more accurate (By comparing the the result with the exact solution based on using the beam 4th order differential equation. It is found to be almost the same with only 2 elements)

3.1 Example 3 redone with 2 elements

The first step is to divide the beam into two elements and number the degrees of freedom and global nodes as follows

There is 6 total degrees of freedom. two at each node. Hence the stiffness matrix for the whole beam (including both elements) will be 6 by 6. For each element however, the same stiffness matrix will be used as above and that will remain as before 4 by 4.

The stiffness matrix for each element is found and then the global stiffness matrix is assembled. Then \(\left \{ p_{global}\right \} =\left [ K_{global}\right ] \left \{ d_{global}\right \} \) is solved as before.  The first step is to move all loads to the nodes as was done before. This is done for each element. The formulas for equivalent loads remain the same, but now \(L\) becomes \(L/2\). The following diagram show the equivalent loading for \(P\)

The equivalent loading for distributed load \(m\) is

Now the above two diagrams are put together to show all equivalent loads with the original reaction forces to obtain the following diagram

\(\left \{ p\right \} =\left [ K\right ] \left \{ d\right \} \) for each element is now constructed. Starting with the first element\[\begin {pmatrix} R_{1}-\frac {m\left ( \frac {L}{2}\right ) }{2}\\ M_{1}-\frac {m\left ( \frac {L}{2}\right ) ^{2}}{12}\\ -\frac {Pb^{2}\left ( \frac {L}{2}+2a\right ) }{\left ( L/2\right ) ^{3}}-\frac {m\left ( \frac {L}{2}\right ) }{2}\\ \frac {m\left ( \frac {L}{2}\right ) ^{2}}{12}-\frac {m\left ( \frac {L}{2}\right ) ^{2}}{12}-\frac {Pab^{2}}{\left ( \frac {L}{2}\right ) ^{2}}\end {pmatrix} =\frac {EI}{\left ( \frac {L}{2}\right ) ^{3}}\begin {pmatrix} 12 & 6\left ( \frac {L}{2}\right ) & -12 & 6\left ( \frac {L}{2}\right ) \\ 6\left ( \frac {L}{2}\right ) & 4\left ( \frac {L}{2}\right ) ^{2} & -6\left ( \frac {L}{2}\right ) & 2\left ( \frac {L}{2}\right ) ^{2}\\ -12 & -6\left ( \frac {L}{2}\right ) & 12 & -6\left ( \frac {L}{2}\right ) \\ 6\left ( \frac {L}{2}\right ) & 2\left ( \frac {L}{2}\right ) ^{2} & -6\left ( \frac {L}{2}\right ) & 4\left ( \frac {L}{2}\right ) ^{2}\end {pmatrix}\begin {Bmatrix} v_{1}\\ \theta _{1}\\ v_{2}\\ \theta _{2}\end {Bmatrix} \] And for the second element\[\begin {pmatrix} -\frac {Pb^{2}\left ( \frac {L}{2}+2a\right ) }{\left ( L/2\right ) ^{3}}-\frac {m\left ( \frac {L}{2}\right ) }{2}\\ \frac {m\left ( \frac {L}{2}\right ) ^{2}}{12}-\frac {m\left ( \frac {L}{2}\right ) ^{2}}{12}-\frac {Pab^{2}}{\left ( \frac {L}{2}\right ) ^{2}}\\ R_{3}-\frac {Pa^{2}\left ( \frac {L}{2}+2b\right ) }{\left ( \frac {L}{2}\right ) ^{3}}-\frac {m\left ( \frac {L}{2}\right ) }{2}\\ \frac {m\left ( \frac {L}{2}\right ) ^{2}}{12}+\frac {Pa^{2}b}{\left ( \frac {L}{2}\right ) ^{2}}\end {pmatrix} =\frac {EI}{\left ( \frac {L}{2}\right ) ^{3}}\begin {pmatrix} 12 & 6\left ( \frac {L}{2}\right ) & -12 & 6\left ( \frac {L}{2}\right ) \\ 6\left ( \frac {L}{2}\right ) & 4\left ( \frac {L}{2}\right ) ^{2} & -6\left ( \frac {L}{2}\right ) & 2\left ( \frac {L}{2}\right ) ^{2}\\ -12 & -6\left ( \frac {L}{2}\right ) & 12 & -6\left ( \frac {L}{2}\right ) \\ 6\left ( \frac {L}{2}\right ) & 2\left ( \frac {L}{2}\right ) ^{2} & -6\left ( \frac {L}{2}\right ) & 4\left ( \frac {L}{2}\right ) ^{2}\end {pmatrix}\begin {Bmatrix} v_{2}\\ \theta _{2}\\ v_{3}\\ \theta _{3}\end {Bmatrix} \] The 2 systems above are assembled to obtain the global stiffness matrix equation giving\[\begin {pmatrix} R_{1}-\frac {m\left ( \frac {L}{2}\right ) }{2}\\ M_{1}-\frac {m\left ( \frac {L}{2}\right ) ^{2}}{12}\\ -2\frac {Pb^{2}\left ( \frac {L}{2}+2a\right ) }{\left ( \frac {L}{2}\right ) ^{3}}-2\frac {m\left ( \frac {L}{2}\right ) }{2}\\ -2\frac {Pab^{2}}{\left ( \frac {L}{2}\right ) ^{2}}\\ R_{3}-\frac {Pa^{2}\left ( \frac {L}{2}+2b\right ) }{\left ( \frac {L}{2}\right ) ^{3}}-\frac {m\left ( \frac {L}{2}\right ) }{2}\\ \frac {m\left ( \frac {L}{2}\right ) ^{2}}{12}+\frac {Pa^{2}b}{\left ( \frac {L}{2}\right ) ^{2}}\end {pmatrix} =\frac {EI}{\left ( \frac {L}{2}\right ) ^{3}}\begin {pmatrix} 12 & 6\left ( \frac {L}{2}\right ) & -12 & 6\left ( \frac {L}{2}\right ) & 0 & 0\\ 6\left ( \frac {L}{2}\right ) & 4\left ( \frac {L}{2}\right ) ^{2} & -6\left ( \frac {L}{2}\right ) & 2\left ( \frac {L}{2}\right ) ^{2} & 0 & 0\\ -12 & -6\left ( \frac {L}{2}\right ) & 12+12 & -6\left ( \frac {L}{2}\right ) +6\left ( \frac {L}{2}\right ) & -12 & 6\left ( \frac {L}{2}\right ) \\ 6\left ( \frac {L}{2}\right ) & 2\left ( \frac {L}{2}\right ) ^{2} & -6\left ( \frac {L}{2}\right ) +6\left ( \frac {L}{2}\right ) & 4\left ( \frac {L}{2}\right ) ^{2}+4\left ( \frac {L}{2}\right ) ^{2} & -6\left ( \frac {L}{2}\right ) & 2\left ( \frac {L}{2}\right ) ^{2}\\ 0 & 0 & -12 & -6\left ( \frac {L}{2}\right ) & 12 & -6\left ( \frac {L}{2}\right ) \\ 0 & 0 & 6\left ( \frac {L}{2}\right ) & 2\left ( \frac {L}{2}\right ) ^{2} & -6\left ( \frac {L}{2}\right ) & 4\left ( \frac {L}{2}\right ) ^{2}\end {pmatrix}\begin {Bmatrix} v_{1}\\ \theta _{1}\\ v_{2}\\ \theta _{2}\\ v_{3}\\ \theta _{3}\end {Bmatrix} \] The boundary conditions are now applied giving \(v_{1}=0,\theta _{1}=0\) since the first node is fixed, and \(v_{3}=0\). And putting 1 on the diagonal of the stiffness matrix corresponding to these known boundary conditions results in\[\begin {pmatrix} 0\\ 0\\ -2\frac {Pb^{2}\left ( \frac {L}{2}+2a\right ) }{\left ( L/2\right ) ^{3}}-2\frac {m\left ( L/2\right ) }{2}\\ -2\frac {Pab^{2}}{\left ( L/2\right ) ^{2}}\\ 0\\ \frac {m\left ( \frac {L}{2}\right ) ^{2}}{12}+\frac {Pa^{2}b}{\left ( \frac {L}{2}\right ) ^{2}}\end {pmatrix} =\frac {EI}{\left ( \frac {L}{2}\right ) ^{3}}\begin {pmatrix} 1 & 0 & 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 & 0 & 0\\ 0 & 0 & 24 & 0 & 0 & 6\left ( \frac {L}{2}\right ) \\ 0 & 0 & 0 & 8\left ( \frac {L}{2}\right ) ^{2} & 0 & 2\left ( \frac {L}{2}\right ) ^{2}\\ 0 & 0 & 0 & 0 & 1 & 0\\ 0 & 0 & 6\left ( \frac {L}{2}\right ) & 2\left ( \frac {L}{2}\right ) ^{2} & 0 & 4\left ( \frac {L}{2}\right ) ^{2}\end {pmatrix}\begin {pmatrix} 0\\ 0\\ v_{2}\\ \theta _{2}\\ 0\\ \theta _{3}\end {pmatrix} \] As was mentioned earlier, another method would be to remove the rows/columns which results in\[\begin {pmatrix} -2\frac {Pb^{2}\left ( \frac {L}{2}+2a\right ) }{\left ( L/2\right ) ^{3}}-2\frac {m\left ( \frac {L}{2}\right ) }{2}\\ -2\frac {Pab^{2}}{\left ( \frac {L}{2}\right ) ^{2}}\\ \frac {m\left ( \frac {L}{2}\right ) ^{2}}{12}+\frac {Pa^{2}b}{\left ( \frac {L}{2}\right ) ^{2}}\end {pmatrix} =\frac {EI}{\left ( \frac {L}{2}\right ) ^{3}}\begin {pmatrix} 24 & 0 & 6\left ( \frac {L}{2}\right ) \\ 0 & 8\left ( \frac {L}{2}\right ) ^{2} & 2\left ( \frac {L}{2}\right ) ^{2}\\ 6\left ( \frac {L}{2}\right ) & 2\left ( \frac {L}{2}\right ) ^{2} & 4\left ( \frac {L}{2}\right ) ^{2}\end {pmatrix}\begin {pmatrix} v_{2}\\ \theta _{2}\\ \theta _{3}\end {pmatrix} \] Giving the same solution. There are 3 unknowns to solve for. Once these unknowns are solved for, \(v\left ( x\right ) \) for the first element and for the second element are fully determined. The following code displays the deflection curve for the above beam

clear all; close all; 
P=400; 
L=144; 
E=30*10^6; 
I=57.1; 
m=200; 
a=0.125*L; 
b=0.375*L; 
A=E*I/(L/2)^3*[1       0        0       0        0     0; 
0       1        0       0        0     0; 
0       0        24      0        0    6*L/2; 
0       0        0   8*(L/2)^2    0 2*(L/2)^2; 
0       0        0       0        1     0; 
0       0      6*L/2  2*(L/2)^2   0  4*(L/2)^2]; 
A 
load = [0; 
0; 
-(m*L/2) - 2*P*b^2*(L/2+2*a)/(L/2)^3; 
-2*P*a*b^2/(L/2)^2; 
0; 
P*a^2*b/(L/2)^2+(m*(L/2)^2)/12] 
sol=A\load
 

Gives

A = 
1.0e+08 * 
0.0000         0         0         0         0         0 
0    0.0000         0         0         0         0 
0         0    0.0011         0         0    0.0198 
0         0         0    1.9033         0    0.4758 
0         0         0         0    0.0000         0 
0         0    0.0198    0.4758         0    0.9517 
load = 
0 
0 
-15075 
-8100 
0 
87750 
sol = 
0 
0 
-0.2735 
-0.0019 
0 
0.0076
 

The above solution gives \(v_{2}=-0.2735\) in (downwards displacement) and \(\theta _{2}=-0.0019\) radians and \(\theta _{3}=0.0076\) radians. \(v\left ( x\right ) \) polynomial is now found for each element\begin {align*} v_{elem1}\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\\ v_{2}\\ \theta _{2}\end {Bmatrix} \\ & =\begin {pmatrix} \ \frac {1}{\left ( \frac {L}{2}\right ) ^{3}}\left ( 3\left ( \frac {L}{2}\right ) x^{2}-2x^{3}\right ) & \ \frac {1}{\left ( \frac {L}{2}\right ) ^{2}}\left ( -\left ( \frac {L}{2}\right ) x^{2}+x^{3}\right ) \end {pmatrix}\begin {pmatrix} -0.2735\\ -0.0019 \end {pmatrix} \\ & =\frac {2.188}{L^{3}}\left ( 2x^{3}-\frac {3}{2}Lx^{2}\right ) -\frac {0.007\,6}{L^{2}}\left ( x^{3}-\frac {1}{2}Lx^{2}\right ) \allowbreak \end {align*}

The above polynomial is the transverse deflection of the beam for the region \(0\leq x\) \(\leq L/2\). \(v\left ( x\right ) \) for the second element is found in similar way\begin {align*} v_{elem2}\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_{2}\\ \theta _{2}\\ 0\\ \theta _{3}\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 {pmatrix} -0.2735\\ -0.0019\\ 0\\ 0.0076 \end {pmatrix} \\ & =\begin {pmatrix} N_{1}\left ( x\right ) \ & \ N_{2}\left ( x\right ) & \ N_{4}\left ( x\right ) \end {pmatrix}\begin {pmatrix} -0.2735\\ -0.0019\\ 0.0076 \end {pmatrix} \\ & =\begin {pmatrix} \frac {1}{\left ( \frac {L}{2}\right ) ^{3}}\left ( \left ( \frac {L}{2}\right ) ^{3}-3\left ( \frac {L}{2}\right ) x^{2}+2x^{3}\right ) \ & \ \frac {1}{\left ( \frac {L}{2}\right ) ^{2}}\left ( \left ( \frac {L}{2}\right ) ^{2}x-2\left ( \frac {L}{2}\right ) x^{2}+x^{3}\right ) & \ \frac {1}{\left ( \frac {L}{2}\right ) ^{2}}\left ( -\left ( \frac {L}{2}\right ) x^{2}+x^{3}\right ) \end {pmatrix}\begin {pmatrix} -0.2735\\ -0.0019\\ 0.0076 \end {pmatrix} \end {align*}

Hence \[ v_{elem2}\left ( x\right ) =\frac {0.030\,4}{L^{2}}\left ( x^{3}-\frac {1}{2}Lx^{2}\right ) -\frac {0.007\,6}{L^{2}}\left ( \frac {1}{4}L^{2}x-Lx^{2}+x^{3}\right ) -\frac {2.188}{L^{3}}\left ( \frac {1}{8}L^{3}-\frac {3}{2}Lx^{2}+2x^{3}\right ) \allowbreak \] Which is valid for \(\allowbreak L/2\leq x\) \(\leq L\). The following is a plot of the deflection curve using the above 2 equations

When the above plot is compared to the case with one element, the deflection is seen to be larger now. Comparing the above to the analytical solution shows that the deflection now is almost exactly the same as the analytical solution. Hence by using only two elements instead of one element, the solution has become more accurate and almost agrees with the analytical solution.

The following diagram shows the deflection curve of problem three above when using one element and two elements on the same plot to help illustrate the difference in the result more clearly.

The analytical deflection for the beam in problem three above (fixed on the left and simply supported at the right end) when there is uniformly loaded with \(w\) lbs per unit length is given by \[ v\left ( x\right ) =-\frac {wx^{2}}{48EI}\left ( 3L^{2}-5Lx+2x^{2}\right ) \] While the analytical deflection for the same beam but when there is a point load P at distance \(a\) from the left end is given by \[ v\left ( x\right ) =-P(\left \langle L-a\right \rangle ^{3}(3L-x)x^{2}+L^{2}(3(L-a)(L-x)x^{2}+2L\left \langle x-a\right \rangle ^{3})) \] Therefore, the analytical expression for deflection is given by the sum of the above expressions, giving\[ v\left ( x\right ) =-\frac {wx^{2}}{48EI}\left ( 3L^{2}-5Lx+2x^{2}\right ) -P(\left \langle L-a\right \rangle ^{3}(3L-x)x^{2}+L^{2}(3(L-a)(L-x)x^{2}+2L\left \langle x-a\right \rangle ^{3})) \] Where \(\left \langle x-a\right \rangle \) means it is zero when \(x-a\) is negative. In other words \(\left \langle x-a\right \rangle =\left ( x-a\right ) UnitStep\left ( x-a\right ) \)

The following diagram is a plot of the analytical deflection with the two elements deflection calculated using Finite elements above.

In the above, the blue dashed curve is the analytical solution, and the red curve is the finite elements solution using 2 elements. It can be seen that the finite element solution for the deflection is now in a very good agreement with the analytical solution.

Chapter 4
Generating shear and bending moments diagrams

After solving the problem using finite elements and obtaining the field displacement function \(v\left ( x\right ) \) as was shown in the above examples, the shear force and bending moments along the beam can be calculated. Since the bending moment is given by \(M\left ( x\right ) =-EI\frac {d^{2}v\left ( x\right ) }{dx^{2}}\) and shear force is given by \(V\left ( x\right ) =\frac {dM}{dx}=-EI\frac {d^{3}v\left ( x\right ) }{dx^{3}}\) then these diagrams are now readily plotted as shown below for example three above using the result from the finite elements with 2 elements. Recalling from above that \[ v\left ( x\right ) =\left . \begin {array} [c]{c}\frac {2.188}{L^{3}}\left ( 2x^{3}-\frac {3}{2}Lx^{2}\right ) -\frac {0.007\,6}{L^{2}}\left ( x^{3}-\frac {1}{2}Lx^{2}\right ) \\ \frac {0.030\,4}{L^{2}}\left ( x^{3}-\frac {1}{2}Lx^{2}\right ) -\frac {0.007\,6}{L^{2}}\left ( \frac {1}{4}L^{2}x-Lx^{2}+x^{3}\right ) -\frac {2.188}{L^{3}}\left ( \frac {1}{8}L^{3}-\frac {3}{2}Lx^{2}+2x^{3}\right ) \end {array} \right \} \begin {array} [c]{c}0\leq x\leq L/2\\ L/2\leq x\leq L \end {array} \] Hence\[ M\left ( x\right ) =-EI\left . \begin {array} [c]{c}\frac {-0.0076\left ( 6x-L\right ) }{L^{2}}+\frac {2.188\left ( 12x-3L\right ) }{L^{3}}\\ -\frac {0.0076\left ( 6x-2L\right ) }{L^{2}}+\frac {0.0304\left ( 6x-L\right ) }{L^{2}}-\frac {2.188\left ( 12x-3L\right ) }{L^{3}}\end {array} \right \} \begin {array} [c]{c}0\leq x\leq L/2\\ L/2\leq x\leq L \end {array} \] using \(E=30\times 10^{6}\) psi and \(I=57\) in\(^{4}\) and \(L=144\) in, the bending moment diagram plot is

The bending moment diagram clearly does not agree with the bending moment diagram that can be generated from the analytical solution given below (generated using my other program which solves this problem analytically)

The reason for this is because the solution \(v\left ( x\right ) \) obtained using the finite elements method is a third degree polynomial and after differentiating twice to obtain the bending moment (\(M\left ( x\right ) =-EI\frac {d^{2}v}{dx^{2}}\)) the result becomes a linear function in \(x\) while in the analytical solution case, when the load is distributed, the solution \(v\left ( x\right ) \) is a fourth degree polynomial. Hence the bending moment will be quadratic function in \(x\) in the analytical case.

Therefore, in order to obtain good approximation for the bending moment and shear force diagrams using finite elements, more elements will be needed.

Chapter 5
Finding the stiffness matrix using methods other than direct method

5.1 Virtual work method for derivation of the stiffness matrix
5.2 Potential energy (minimize a functional) method to derive the stiffness matrix

There are three main methods to obtain the stiffness matrix

  1. Variational method (minimizing a functional). This functional is the potential energy of the structure and loads.
  2. Weighted residual. Requires the differential equation as a starting point. Approximated in weighted average. Galerkin weighted residual method is the most common method for implementation.
  3. Virtual work method. Making the virtual work zero for an arbitrary allowed displacement.

5.1 Virtual work method for derivation of the stiffness matrix

In virtual work method, a small displacement is assumed to occur. Looking at small volume element, the amount of work done by external loads to cause the small displacement is set equal to amount of increased internal strain energy. Assuming the field of displacement is given by \(\mathbf {u}=\left \{ u,v,w\right \} \) and assuming the external loads are given by \(\left \{ p\right \} \) acting on the nodes, hence these point loads will do work given by \(\left \{ \delta d\right \} ^{T}\left \{ p\right \} \) on that unit volume where \(\left \{ d\right \} \) is the nodal displacements. In all these derivations, only loads acting directly on the nodes are considered for now. In other words, body forces and traction forces are not considered in order to simplify the derivations.

The increase of strain energy is \(\left \{ \delta \varepsilon \right \} ^{T}\left \{ \sigma \right \} \) in that same unit volume.

Hence, for a unit volume \[ \left \{ \delta d\right \} ^{T}\left \{ p\right \} =\left \{ \delta \varepsilon \right \} ^{T}\left \{ \sigma \right \} \] And for the whole volume\begin {equation} \left \{ \delta d\right \} ^{T}\left \{ p\right \} =\int _{V}\left \{ \delta \varepsilon \right \} ^{T}\left \{ \sigma \right \} dV \tag {1} \end {equation} Assuming that displacement can be written as a function of the nodal displacements of the element results in\[ \mathbf {u=}\left [ N\right ] \left \{ d\right \} \] Therefore \begin {align} \delta \mathbf {u} & \mathbf {=}\left [ N\right ] \left \{ \delta d\right \} \nonumber \\ \left \{ \delta \mathbf {u}\right \} ^{T} & \mathbf {=}\left \{ \delta d\right \} ^{T}\left [ N\right ] ^{T} \tag {2} \end {align}

Since \(\left \{ \varepsilon \right \} =\partial \left \{ \mathbf {u}\right \} \) then \(\left \{ \varepsilon \right \} =\partial \left [ N\right ] \left \{ d\right \} =\) \(\left [ B\right ] \left \{ d\right \} \) where \(B\) is the strain displacement matrix \(\left [ B\right ] =\partial \left [ N\right ] \), hence\begin {align} \left \{ \varepsilon \right \} & =\left [ B\right ] \left \{ d\right \} \nonumber \\ \left \{ \delta \varepsilon \right \} & =\left [ B\right ] \left \{ \delta d\right \} \nonumber \\ \left \{ \delta \varepsilon \right \} ^{T} & =\left \{ \delta d\right \} ^{T}\left [ B\right ] ^{T} \tag {3} \end {align}

Now from the stress-strain relation \(\left \{ \sigma \right \} =\left [ E\right ] \left \{ \varepsilon \right \} \), hence\begin {equation} \left \{ \sigma \right \} =\left [ E\right ] \left [ B\right ] \left \{ d\right \} \tag {4} \end {equation} Substituting Eqs. (2,3,4) into (1) results in\[ \left \{ \delta d\right \} ^{T}\left \{ p\right \} -\int _{V}\left \{ \delta d\right \} ^{T}\left [ B\right ] ^{T}\left [ E\right ] \left [ B\right ] \left \{ d\right \} dV=0 \] Since \(\left \{ \delta d\right \} \) and \(\left \{ d\right \} \) do not depend on the integration variables they can be moved outside the integral, giving\[ \left \{ \delta d\right \} ^{T}\left ( \left \{ p\right \} -\left \{ d\right \} \int _{V}\left [ B\right ] ^{T}\left [ E\right ] \left [ B\right ] dV\right ) =0 \] Since the above is true for any admissible \(\delta d\) then the only condition is that \[ \left \{ p\right \} =\left \{ d\right \} \int _{V}\left [ B\right ] ^{T}\left [ E\right ] \left [ B\right ] dV \] This is in the form \(P=K\Delta \), therefore \[ \left [ K\right ] =\int _{V}\left [ B\right ] ^{T}\left [ E\right ] \left [ B\right ] dV \] knowing \(\left [ B\right ] \) allows finding \(\left [ k\right ] \) by integrating over the volume. For the beam element though, \(\mathbf {u=}v\left ( x\right ) \) the transverse displacement. This means \(\left [ B\right ] =\frac {d^{2}}{dx^{2}}\left [ N\right ] \). Recalling from the above that for the beam element,\[ \left [ N\right ] =\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} \] Hence \[ \left [ B\right ] =\frac {d^{2}}{dx^{2}}\left [ N\right ] =\begin {pmatrix} \frac {1}{L^{3}}\left ( -6L+12x\right ) \ & \frac {1}{L^{2}}\left ( -4L+6x\right ) & \ \frac {1}{L^{3}}\left ( 6L-12x\right ) & \ \frac {1}{L^{2}}\left ( -2L+6x\right ) \end {pmatrix} \] Hence\begin {align*} \left [ K\right ] & =\int _{V}\left [ B\right ] ^{T}\left [ E\right ] \left [ B\right ] dV\\ & =\int _{V}\begin {Bmatrix} \frac {1}{L^{3}}\left ( -6L+12x\right ) \\ \frac {1}{L^{2}}\left ( -4L+6x\right ) \\ \ \frac {1}{L^{3}}\left ( 6L-12x\right ) \\ \frac {1}{L^{2}}\left ( -2L+6x\right ) \end {Bmatrix} E\begin {pmatrix} \frac {1}{L^{3}}\left ( -6L+12x\right ) \ & \frac {1}{L^{2}}\left ( -4L+6x\right ) & \ \frac {1}{L^{3}}\left ( 6L-12x\right ) & \ \frac {1}{L^{2}}\left ( -2L+6x\right ) \end {pmatrix} dV\\ & =EI\int _{0}^{L}\begin {Bmatrix} \frac {1}{L^{3}}\left ( -6L+12x\right ) \\ \frac {1}{L^{2}}\left ( -4L+6x\right ) \\ \ \frac {1}{L^{3}}\left ( 6L-12x\right ) \\ \frac {1}{L^{2}}\left ( -2L+6x\right ) \end {Bmatrix}\begin {pmatrix} \frac {1}{L^{3}}\left ( -6L+12x\right ) \ & \frac {1}{L^{2}}\left ( -4L+6x\right ) & \ \frac {1}{L^{3}}\left ( 6L-12x\right ) & \ \frac {1}{L^{2}}\left ( -2L+6x\right ) \end {pmatrix} dx\\ & =EI\int _{0}^{L}\begin {pmatrix} \frac {1}{L^{6}}\left ( 6L-12x\right ) ^{2} & \frac {1}{L^{5}}\left ( 4L-6x\right ) \left ( 6L-12x\right ) & -\frac {1}{L^{6}}\left ( 6L-12x\right ) ^{2} & \frac {1}{L^{5}}\left ( 2L-6x\right ) \left ( 6L-12x\right ) \\ \frac {1}{L^{5}}\left ( 4L-6x\right ) \left ( 6L-12x\right ) & \frac {1}{L^{4}}\left ( 4L-6x\right ) ^{2} & -\frac {1}{L^{5}}\left ( 4L-6x\right ) \left ( 6L-12x\right ) & \frac {1}{L^{4}}\left ( 2L-6x\right ) \left ( 4L-6x\right ) \\ -\frac {1}{L^{6}}\left ( 6L-12x\right ) ^{2} & -\frac {1}{L^{5}}\left ( 4L-6x\right ) \left ( 6L-12x\right ) & \frac {1}{L^{6}}\left ( 6L-12x\right ) ^{2} & -\frac {1}{L^{5}}\left ( 2L-6x\right ) \left ( 6L-12x\right ) \\ \frac {1}{L^{5}}\left ( 2L-6x\right ) \left ( 6L-12x\right ) & \frac {1}{L^{4}}\left ( 2L-6x\right ) \left ( 4L-6x\right ) & -\frac {1}{L^{5}}\left ( 2L-6x\right ) \left ( 6L-12x\right ) & \frac {1}{L^{4}}\left ( 2L-6x\right ) ^{2}\end {pmatrix} \allowbreak dx\\ & =EI\begin {pmatrix} \frac {12}{L^{3}} & \frac {6}{L^{2}} & -\frac {12}{L^{3}} & \frac {6}{L^{2}}\\ \frac {6}{L^{2}} & \frac {4}{L} & -\frac {6}{L^{2}} & \frac {2}{L}\\ -\frac {12}{L^{3}} & -\frac {6}{L^{2}} & \frac {12}{L^{3}} & -\frac {6}{L^{2}}\\ \frac {6}{L^{2}} & \frac {2}{L} & -\frac {6}{L^{2}} & \frac {4}{L}\end {pmatrix} \\ & =\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} \end {align*}

Which is the stiffness matrix found earlier.

5.2 Potential energy (minimize a functional) method to derive the stiffness matrix

This method is very similar to the first method actually. It all comes down to finding a functional, which is the potential energy of the system, and minimizing this with respect to the nodal displacements. The result gives the stiffness matrix.

Let the system total potential energy by called \(\Pi \) and let the total internal energy in the system be \(U\) and let the work done by external loads acting on the nodes be \(\Omega \), then\[ \Pi =U-\Omega \] Work done by external loads have a negative sign since they are an external agent to the system and work is being done onto the system. The internal strain energy is given by \(\frac {1}{2}\int _{V}\left \{ \sigma \right \} ^{T}\left \{ \varepsilon \right \} dV\) and the work done by external loads is \(\left \{ d\right \} ^{T}\left \{ p\right \} \), hence\begin {equation} \Pi =\frac {1}{2}\int _{V}\left \{ \sigma \right \} ^{T}\left \{ \varepsilon \right \} dV-\left \{ d\right \} ^{T}\left \{ p\right \} \tag {1} \end {equation} Now the rest follows as before. Assuming that displacement can be written as a function of the nodal displacements \(\left \{ d\right \} \), hence \[ \mathbf {u=}\left [ N\right ] \left \{ d\right \} \] Since \(\left \{ \varepsilon \right \} =\partial \left \{ \mathbf {u}\right \} \) then \(\left \{ \varepsilon \right \} =\partial \left [ N\right ] \left \{ d\right \} =\) \(\left [ B\right ] \left \{ d\right \} \) where \(B\) is the strain displacement matrix \(\left [ B\right ] =\partial \left [ N\right ] \), hence\[ \left \{ \varepsilon \right \} =\left [ B\right ] \left \{ d\right \} \] Now from the stress-strain relation \(\left \{ \sigma \right \} =\left [ E\right ] \left \{ \varepsilon \right \} \), hence\begin {equation} \left \{ \sigma \right \} ^{T}=\left \{ d\right \} ^{T}\left [ B\right ] ^{T}\left [ E\right ] \tag {4} \end {equation} Substituting Eqs. (2,3,4) into (1) results in\[ \Pi =\frac {1}{2}\int _{V}\left \{ d\right \} ^{T}\left [ B\right ] ^{T}\left [ E\right ] \left [ B\right ] \left \{ d\right \} dV-\left \{ d\right \} ^{T}\left \{ p\right \} \] Setting \(\frac {\partial \Pi }{\partial \left \{ d\right \} }=0\) gives\[ 0=\left \{ d\right \} \int _{V}\left [ B\right ] ^{T}\left [ E\right ] \left [ B\right ] dV-\left \{ p\right \} \] Which is on the form \(P=\left [ K\right ] D\) which means that\[ \left [ K\right ] =\int _{V}\left [ B\right ] ^{T}\left [ E\right ] \left [ B\right ] dV \] As was found by the virtual work method.

Chapter 6
References

  1. A first course in Finite element method, 3rd edition, by Daryl L. Logan
  2. Matrix analysis of framed structures, 2nd edition, by William Weaver, James Gere

1Instead 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