Chapter 3
Finite elements (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 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.