Finite Element Method (FEM) in Practice

Solving a Simple Beam Problem by FEM

An Interactive Example

Beams are components which are subjected to bending. Analysis of beams is a fundamental part of solid mechanics (strength of materials). Therefore, a beam in transverse loading can be a reasonable starting point to apply the FEM.

1. Mesh

   Meshing is the process of dividing the geometry (domain) of the problem into a finite number of elements with simple shapes. The vertices of the elements are called nodes which are crucial in the FEM formulation.

   Here, meshing is quite minimal as there are only two elements as shown above, and they are enough for this problem. The nodes are markes as red dots in the figure above.

2. Shape Function

   Having the geometry reduced to a collection of elements by meshing, we now look for a relation to relate the variables of interest throughout the elements, to the values of those variable at the nodes, thus, roughly speaking, reducing the element to its nodes. Such a mathematical equation is referred to as the shape function for that element.

   Beams do deflect and bend when loaded transversely, therefore the shape function here needs to relates deflections (i.e. transverse displacements and rotations) at points throughout the element to the nodal deflections and rotations:

   The differential equation for a beam element is: $\frac{d^{4}y}{d{x}^4}=0$. (1)

   This is due to the fact that for transversely loaded beams, if $w\left(x\right)$ is the transverse load per unit length and $y\left(x\right)$ the transverse deflection: $\frac{d^{4}y}{d{x}^4}=$ $w\left(x\right)$. And for an element, loads are only applied at the nodes, i.e. there is no load inside the element.

   therefore a 3rd degree polynomial satisfies the equation: $y\left(x\right) = a_0+a_{1}x +a_2{x}^2+a_3{x}^3$. (2)

   If such a polynomial is used to derive the shape function, the solution will be exact.

   The relation between rotation and deflection is : $\theta \left(x\right)=$ $\frac{dy}{dx}$. Thus, here: $\theta \left(x\right)=a_1+2a_{2}x+3a_3{x}^2$. (3)

   Now using the above equations for deflection and rotation, coefficients $a_0$ to $a_3$ can be written in terms of transverse deflections and rotations at the two nodes and the length of the element.

   Let $y_0$ and ${\theta }_0$, be the transverse deflection and rotation at $x=0$, and $y_1$ and ${\theta }_1$ the transverse deflection and rotation at $x=L$, respectively. Then for $x=0$ and $x=L$ equations (2) and (3) will give:

   $a_0=y_0$ ,     $a_1={\theta }_0$ ,         $y_1=a_0+a_{1}L+a_2{L}^2+a_3{L}^3$ ,     and     ${\theta }_1=a_1+2a_{2}L+3a_3{L}^3$.

   Trying to write $a_0$, $a_1$ $a_2$, $a_3$ in terms of $y_0$, $y_1$, ${\theta }_0$, ${\theta }_1$ will give:

   $a_0=y_0$ , $a_1={\theta }_0$ , $a_2=\frac{3}{L^2}(y_1-y_0)-\frac{1}{L}(2{\theta }_0+{\theta }_1)$ , $a_3=\frac{2}{L^3}(y_1-y_0)-\frac{1}{L^2}({\theta }_0+{\theta }_1)$

   Substituting theses into equation (2), and rearranging gives:

   $y\left(x\right)=$ $\left[1-3{\left(\frac{x}{L}\right)}^2+2{\left(\frac{x}{L}\right)}^3\right]y_0+\left[x-2\frac{x^2}{L}+\frac{x^3}{L^2}\right]{\theta }_0+\left[3{\left(\frac{x}{L}\right)}^2-2{\left(\frac{x}{L}\right)}^3\right]y_1+\left[-\left(\frac{x^2}{L}\right)+\frac{x^3}{L^2}\right]{\theta }_1$.

   The last equation relates the quantities of deflection and rotation for every point throughout the element to the deflection and rotation at the nodes. Therefore it defines the shape function.

   Let the coefficients of nodal displacements and rotations be named as:

   $N_{y_0}=\left[1-3{\left(\frac{x}{L}\right)}^2+2{\left(\frac{x}{L}\right)}^3\right]$ , $N_{{\theta }_0}=\left[x-2\frac{x^2}{L}+\frac{x^3}{L^2}\right]$ , $N_{y_1}=\left[3{\left(\frac{x}{L}\right)}^2-2{\left(\frac{x}{L}\right)}^3\right]$ , $N_{{\theta }_1}=\left[-\left(\frac{x^2}{L}\right)+\frac{x^3}{L^2}\right]$

   The shape function can be written in matrix form:         $y\left(x\right)=$ $\left\{N_{y_0} N_{{\theta }_0} N_{y_1} N_{{\theta }_1}\right\}\left\{\begin{array}{l}{y}_0\\ {\theta }_0\\ y_1\\ {\theta }_1\end{array}\right\}$ .(4)

3. Element Stiffness Matrix

   Now, we are looking for an equation to describe the load-deflection behavior of individual elements. A linear equation like the spring formula $F=Kx$ , which in FEM will appear as a matrix equation relating loads applied at the nodes of an element to deflections and rotations of the nodes: $\left\{{\mathbf{F}}_e\right\}=\left[{\mathbf{K}}_e\right]\left\{{\mathbf{u}}_e\right\}$ , where $\left\{{\mathbf{F}}_e\right\}$ is the column matrix (vector) of the applied nodal loads, $\left\{{\mathbf{u}}_e\right\}$ the nodal deflections (and rotations) vector, and $\left[{\mathbf{K}}_e\right]$ is the element stiffness matrix.

   For deriving the element stiffness matrix, total potential energy methods can be used, which is used here. Using the spring analogy, the approach is to reach at an equation similar to U=12Kx2, which here will be in the form of: $U_e=\frac{1}{2}{\left\{{\mathbf{u}}_e\right\}}^T\left[{\mathbf{K}}_e\right]\left\{{\mathbf{u}}_e\right\}$. Here, $\left\{{\mathbf{u}}_e\right\}$ is the nodal deflections vector and as said before, $\left[{\mathbf{K}}_e\right]$ is the element stiffness matrix.

   For a beam in bending, Ub, the potential energy due to bending moment M(x) along the beam is:

   $U_b={\int }_0^L\frac{M^2}{2EI}dx$. (5.a)

   For beams, it is also known that: $M\left(x\right)=$ $EI\frac{d^{2}y}{d{x}^2}$ . (5.b)

   For an element $EI$ is constant, therefore substituting equation (5.b) into equation (5.a) gives:

   $U_b=\frac{EI}{2}{\int }_0^L(\frac{d^{2}y}{d{x}^2}{)}^{2}dx$ . (5.c)

   From equation (4), we can have:

   $\frac{d^2}{d{x}^2}$ $\left(N_{y_0}\right)=$ $-\frac{6}{L^2}+\frac{12}{L^3}x$ ,   $\frac{d^2}{d{x}^2}$$\left(N_{{\theta }_0}\right)=$$-\frac{4}{L}+\frac{6}{L^2}x$ ,   $\frac{d^2}{d{x}^2}$$\left(N_{y_1}\right)=$$\frac{6}{L^2}-\frac{12}{L^3}x$ ,   $\frac{d^2}{d{x}^2}$$\left(N_{{\theta }_1}\right)=$$-\frac{2}{L}+\frac{6}{L^2}x$ .

    

   Let us say:     $\left\{\frac{d^2{N}_{y_0}}{d{x}^2}\frac{d^2{N}_{{\theta }_0}}{d{x}^2}\frac{d^2{N}_{y_1}}{d{x}^2}\frac{d^2{N}_{{\theta }_1}}{d{x}^2}\right\}=\left\{B_1{B}_2{B}_3{B}_4\right\}=\left\{\mathbf{B}\right\}$ . (5.d)

   If $\left\{y_0 {\theta }_0 y_1 {\theta }_1\right\}={\left\{\mathbf{u}\right\}}^T$ (superscript $T$ is the transpose operator) , we will have: $\frac{d^{2}y}{d{x}^2}=\left\{\mathbf{B}\right\}\left\{\mathbf{u}\right\}$. (5.e)

    

   For every $n\times 1$, i.e. column vector $\left\{v\right\}$:     $v^2={\left|v\right|}^2={\left\{v\right\}}^T\left\{v\right\}$ . Therefore substituting (5.e) in (5.c) and rearrangement using matrix algebra will give:

    

       $U_b=$ $\frac{EI}{2}{\int }_0^L{\left\{\mathbf{u}\right\}}^T{\left\{\mathbf{B}\right\}}^T\left\{\mathbf{B}\right\}\left\{\mathbf{u}\right\}dx$ . (5.f)

    

   Equation (4) shows that $\left. {\mathbf{u}\right\}$ and therefore also ${\left\{\mathbf{u}\right\}}^T$ are independent of the variable of the integral in equation (5.f),

   thus:     $U_b=\frac{1}{2}{\left\{\mathbf{u}\right\}}^T\left[\left(EI\right){\int }_0^L{\left\{\mathbf{B}\right\}}^T\left\{\mathbf{B}\right\}dx\right]\left\{\mathbf{u}\right\}$ . (5.g)

   If we say:     Ke = EI ∫0L BT B d x , (5.h)

   equation (5.g) can be written as:     $U_b=\frac{1}{2}{\left\{\mathbf{u}\right\}}^T\left[K_e\right]\left\{\mathbf{u}\right\}$ . (5.i)

   Equation (5.d) means that:     ${\left\{\mathbf{B}\right\}}^T\left\{\mathbf{B}\right\}=$ $\left[\begin{array}{cccc}{B}_1^2 & B_1{B}_2 & B_1{B}_3& B_1{B}_4\\ B_1{B}_2& B_2^2& B_2{B}_3 & B_2{B}_4\\ B_1{B}_3& B_2{B}_3& B_3^2& B_3{B}_4\\ B_1{B}_4& B_2{B}_4& B_3{B}_4& B_4^2\end{array}\right]$ . (5.j)

    

   Equation (5.h) shows that integrating the elements of the right-hand side of equation (5.j) will give the element stiffness matrix for a beam in pure bending:

   For a beam element with Young's moulus of $E$, cross section moment of inertia of $I$ and length of $L$, the following matrix equation holds:

   $\frac{EI}{L^3}\left[\begin{array}{cccc}12& 6L& -12& 6L\\ 6L& 4L^2& -6L& 2L^2\\ -12& -6L& 12& -6L\\ 6L& 2L^2& -6L& 4L^2\end{array}\right]$ (6) $\left\{\begin{array}{l}{y}_0\\ {\theta }_0\\ y_1\\ {\theta }_1\end{array}\right\}$ $=\left\{\begin{array}{l}{P}_0\\ M_0\\ P_1\\ M_1\end{array}\right\}$ ,           which can be written as:   $\left[K_e\right]\left\{u_e\right\}=\left\{F_e\right\}$ .

   In equation (6), $y_0$, ${\theta }_0$, $P_0$ and $M_0$ are the transverse deflection, rotation, transverse load and moment at $x=0$, and $y_1$, ${\theta }_1$, $P_1$ and $M_1$ are the transverse deflection, rotation, transverse load and moment at $x=L$ , respectively.

   The Above equation, can therefore serve as the beam element stiffness matrix for a beam element of length L.

   for a beam with a rectangular cross section, $I=\frac{b{h}^3}{12}$, (7)

   where $b$ and $h$ are the (flange) width and depth of the beam, respectively.

   Element (1):

   Young's modulus, $E$ [GPa] :   |   length, $L$ [m]

   width, $b$ [cm]   |   depth, $h$ [cm] , then element stiffness matrix is:

   $\left[K_e\right]=\left[\begin{array}{}\\ \\ \\ \end{array}\right.$

   12		-12
   -12		12

   $\left.\begin{array}{}\\ \\ \\ \end{array}\right]$

   (8)

   
   

   Element (2):

   Young's modulus, $E$ [GPa] :   |   length, $L$ [m]

   width, $b$ [cm]   |   depth, $h$ [cm] , then element stiffness matrix is:

   $\left[K_e\right]=\left[\begin{array}{}\\ \\ \\ \end{array}\right.$

   12		-12
   -12		12

   $\left.\begin{array}{}\\ \\ \\ \end{array}\right]$

   (9)

   
   

4. Assembling Global Stiffness Matrix

   The equation (10) shows how the stiffness matrices of the individual elements are combined to make the global stiffness matrix.

   $\left[\begin{array}{cccccc}{K}_{11}^1& K_{12}^1& K_{13}^1& K_{14}^1& 0& 0\\ K_{21}^1& K_{22}^1& K_{23}^1& K_{24}^1& 0& 0\\ K_{31}^1& K_{32}^1& (K_{33}^1+K_{11}^2)& (K_{34}^1+K_{12}^2)& K_{13}^2& K_{14}^2\\ K_{41}^1& K_{42}^1& (K_{43}^1+K_{21}^2)& (K_{44}^1+K_{22}^2)& K_{23}^2& K_{24}^2\\ 0& 0& K_{31}^2& K_{32}^2& K_{33}^2& K_{34}^2\\ 0& 0& K_{41}^2& K_{42}^2& K_{43}^2& K_{44}^2\end{array}\right]$ (10) $\left\{\begin{array}{l}{y}_0\\ {\theta }_0\\ y_1\\ {\theta }_1\\ y_2\\ {\theta }_2\end{array}\right\}$ $=\left\{\begin{array}{l}{P}_0\\ M_0\\ P_1\\ M_1\\ P_2\\ M_2\end{array}\right\}$ ,   or:     $\left[K_g\right]\left\{u_g\right\}=\left\{F_g\right\}$

   
   
   

5. Applying Boundary Conditions

   Equation (11) shows how applying BCs makes the problem to be simplified as corresponding rows and columns are deleted for each BC.

   BCs:

   $y_0=0$

   ${\theta }_0=0$

   $y_1=0$

   ${\theta }_1=0$

   $y_2=0$

   ${\theta }_2=0$

   Loads:

   $P_0=$ [N]

   $M_0=$ 0 [Nm]

   $P_1=$ [N]

   $M_1=$ 0 [Nm]

   $P_2=$ [N]

   $M_2=$ 0 [Nm]

   (11)

   $\left[\begin{array}{}\\ \\ \\ \\ \\ \end{array}\right.$

   				0	0
   				0	0
   0	0
   0	0

   $\left.\begin{array}{}\\ \\ \\ \\ \\ \end{array}\right]$ $\left\{\begin{array}{}\\ \\ \\ \\ \\ \end{array}\right.$

   $y_0$

   ${\theta }_0$

   $y_1$

   ${\theta }_1$

   $y_2$

   ${\theta }_2$

   $\left.\begin{array}{l}\\ \\ \\ \\ \\ \end{array}\right\}$ $=\left\{\begin{array}{}\\ \\ \\ \\ \\ \end{array}\right.$

   $P_0$

   $M_0$

   $P_1$

   $M_1$

   $P_2$

   $M_2$

   $\left.\begin{array}{l}\\ \\ \\ \\ \\ \end{array}\right\}$

6. Solving the Reduced Matrix Equation and Results

   The reduced global stiffness matrix (after application of the BCs) is shown below, as well as the final results. The results for displacements are in meters and rotations in radians.

   (12)