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

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.