PDF (letter size)

Engineering report

Derivation of stiﬀness matrix for a beam

July 7, 2016   Compiled on May 25, 2022 at 7:46pm

## Chapter 1Introduction

A short review for solving the beam problem in 2D is given. The deﬂection curve, bending moment and shear force diagrams are calculated for a beam subject to bending moment and shear force using direct stiﬀness method and then using ﬁnite elements method by adding more elements. The problem is solved ﬁrst by ﬁnding the stiﬀness matrix using the direct method and then using the virtual work method.

## Chapter 2Direct method

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 ﬁeld (the unknown being solved for) will be a function of one independent variable which is the $$x$$ coordinate. The displacement ﬁeld 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 ﬁeld $$v\left ( x\right )$$ In the direct method of ﬁnding the stiﬀness matrix, the forces at the ends of the beam are found directly by the use of beam theory. In beam theory the signs are diﬀerent from what is given in the ﬁrst diagram above. Therefore, the moment and shear forces obtained using beam theory ($$M_{B}$$ and $$V_{B}$$ in the diagram below) will have diﬀerent signs when compared to the external forces. The signs are then adjusted to reﬂect 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 stiﬀness 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 deﬂection $$\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 deﬂection $$\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 ﬁeld 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 ﬁeld 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 veriﬁed. The stiﬀness 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 stiﬀness 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 stiﬀness matrix means knowing the nodal displacements $$\left \{ d\right \}$$ when given the forces at the nodes. The power of the ﬁnite 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 ﬁnite diﬀerence method. Eq. 5A above can now be used to ﬁnd 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 stiﬀness matrix

2.1.1 Example 1
2.1.2 Example 2
2.1.3 Example 3

The beam stiﬀness 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 ﬁrst 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 ﬁxed). Hence $$v_{1}=0$$ and $$\theta _{1}=0$$. Hence the only unknowns are $$v_{2}$$ and $$\theta _{2}$$. Therefore the ﬁrst 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
];



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 ﬁeld 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 deﬂection 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 stiﬀness matrix. This method might be easier for programming than the ﬁrst 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 stiﬀness 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 stiﬀness 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 ﬁrst two rows and columns, the above system simpliﬁes 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 ﬁrst 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
];



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 ﬁrst 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 ﬁeld is found and the deﬂection 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 ﬁxed 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 ﬁnding 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 stiﬀness method to work. Equivalent load for the above point load 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 stiﬀness 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 ﬁrst, 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 ﬁeld 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 deﬂection 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 3Finite elements or adding more 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 ﬁeld $$v\left ( x\right )$$ becomes more accurate (By comparing the the result with the exact solution based on using the beam 4th order diﬀerential equation. It is found to be almost the same with only 2 elements)

### 3.1 Example 3 redone with 2 elements

The ﬁrst 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 stiﬀness matrix for the whole beam (including both elements) will be 6 by 6. For each element however, the same stiﬀness matrix will be used as above and that will remain as before 4 by 4. The stiﬀness matrix for each element is found and then the global stiﬀness matrix is assembled. Then $$\left \{ p_{global}\right \} =\left [ K_{global}\right ] \left \{ d_{global}\right \}$$ is solved as before.  The ﬁrst 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 ﬁrst 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 stiﬀness 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 ﬁrst node is ﬁxed, and $$v_{3}=0$$. And putting 1 on the diagonal of the stiﬀness 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 ﬁrst element and for the second element are fully determined. The following code displays the deﬂection 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
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]



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
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 deﬂection 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 deﬂection curve using the above 2 equations When the above plot is compared to the case with one element, the deﬂection is seen to be larger now. Comparing the above to the analytical solution shows that the deﬂection 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 deﬂection curve of problem three above when using one element and two elements on the same plot to help illustrate the diﬀerence in the result more clearly. The analytical deﬂection for the beam in problem three above (ﬁxed 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 deﬂection 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 deﬂection 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 deﬂection with the two elements deﬂection calculated using Finite elements above. In the above, the blue dashed curve is the analytical solution, and the red curve is the ﬁnite elements solution using 2 elements. It can be seen that the ﬁnite element solution for the deﬂection is now in a very good agreement with the analytical solution.

## Chapter 4Generating shear and bending moments diagrams

After solving the problem using ﬁnite elements and obtaining the ﬁeld 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 ﬁnite 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 ﬁnite elements method is a third degree polynomial and after diﬀerentiating 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 ﬁnite elements, more elements will be needed.

## Chapter 5Finding the stiﬀness matrix using methods other than direct method

There are three main methods to obtain the stiﬀness matrix

1. Variational method (minimizing a functional). This functional is the potential energy of the structure and loads.
2. Weighted residual. Requires the diﬀerential 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 stiﬀness 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 ﬁeld 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 ﬁnding $$\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 stiﬀness matrix found earlier.

### 5.2 Potential energy (minimize a functional) method to derive the stiﬀness matrix

This method is very similar to the ﬁrst method actually. It all comes down to ﬁnding a functional, which is the potential energy of the system, and minimizing this with respect to the nodal displacements. The result gives the stiﬀness 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.

1. A ﬁrst 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 stiﬀness matrix for that boundary conditions. I will do this example again using this method