Finite Element Method (FEM) in Practice

Solving a Simple Beam Problem by FEM

An Interactive Example

Please note and try:

red   boxes change parameters dynamically.

Yellow   boxes are draggable.

buttons close and open sections

(click for partial and double click for full close and open).

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.

Figure 1. A transversely loaded beam
  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: d4ydx4=0. (1)

    This is due to the fact that for transversely loaded beams, if w(x) is the transverse load per unit length and y(x) the transverse deflection: d4ydx4= w(x). 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(x)  =  a0+a1 x +a2x2+a3x3 . (2)

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

    The relation between rotation and deflection is : θ(x)= dydx . Thus, here: θ(x)=a1+2a2x+3a3x2. (3)

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

    Let y0 and θ0, be the transverse deflection and rotation at x=0, and y1 and θ1 the transverse deflection and rotation at x=L, respectively. Then for x=0 and x=L equations (2) and (3) will give:

    a0=y0 ,     a1=θ0 ,         y1=a0+a1L+a2L2+a3L3 ,     and     θ1=a1+2a2L+3a3L3.

    Trying to write a0, a1 a2, a3 in terms of y0, y1, θ0, θ1 will give:

    a0=y0 , a1=θ0 , a2=3L2(y1-y0)-1L(2θ0+θ1) , a3=2L3(y1-y0)-1L2(θ0+θ1)

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

    y(x)= 1-3xL2+2xL3y0+x-2x2L+x3L2θ0+3xL2-2xL3y1+-x2L+x3L2θ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:

    Ny0=1-3xL2+2xL3 , Nθ0=x-2x2L+x3L2 , Ny1=3xL2-2xL3 , Nθ1=-x2L+x3L2

    The shape function can be written in matrix form:         y(x)= Ny0 Nθ0 Ny1 Nθ1y0θ0y1θ1 .(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: Fe = Ke ue , where Fe is the column matrix (vector) of the applied nodal loads, ue the nodal deflections (and rotations) vector, and Ke 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: Ue = 12 ueT Ke ue . Here, ue is the nodal deflections vector and as said before, Ke is the element stiffness matrix.

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

    Ub=0LM22EIdx. (5.a)

    For beams, it is also known that: M(x)= EId2ydx2 . (5.b)

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

    Ub=EI20L(d2ydx2)2dx . (5.c)

    From equation (4), we can have:

    d2dx2 (Ny0)= -6L2+12L3x ,   d2dx2(Nθ0)=-4L+6L2x ,   d2dx2(Ny1)=6L2-12L3x ,   d2dx2(Nθ1)=-2L+6L2x .


    Let us say:     d2Ny0dx2 d2Nθ0dx2 d2Ny1dx2 d2Nθ1dx2 = B1 B2 B3 B4 = B . (5.d)

    If y0    θ0    y1    θ1 = u T (superscript T is the transpose operator) , we will have: d2ydx2=Bu. (5.e)


    For every n×1, i.e. column vector v:     v2=v2=vTv . Therefore substituting (5.e) in (5.c) and rearrangement using matrix algebra will give:


        Ub= EI20L uT BT B udx . (5.f)


    Equation (4) shows that u and therefore also uT are independent of the variable of the integral in equation (5.f),

    thus:     Ub = 12 uT [ EI 0L BT B dx ] u . (5.g)

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

    equation (5.g) can be written as:     Ub = 12 uT Ke u . (5.i)

    Equation (5.d) means that:     BT B = B12   B1B2  B1B3  B1B4B1B2B22B2B3 B2B4B1B3B2B3B32B3B4B1B4B2B4B3B4B42 . (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:

    EIL3 12 6L -12 6L 6L 4L2 -6L 2L2 -12 -6L 12 -6L 6L 2L2 -6L 4L2 (6) y0 θ0 y1 θ1 = P0 M0 P1 M1 ,           which can be written as:   Ke ue = Fe .

    In equation (6), y0, θ0, P0 and M0 are the transverse deflection, rotation, transverse load and moment at x=0, and y1, θ1, P1 and M1 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=bh312, (7)

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

    Element (1):


      |   , then element stiffness matrix is:

    Ke1 =
    12 -12
    -12 12


    Element (2):


      |   , then element stiffness matrix is:

    Ke2 =
    12 -12
    -12 12


  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.

    K111 K121 K131 K141 00 K211 K221 K231 K241 00 K311 K321 (K331+K112) (K341+K122) K132 K142 K411 K421 (K431+K212) (K441+K222) K232 K242 00 K312 K322 K332 K342 00 K412 K422 K432 K442 (10) y0 θ0 y1 θ1 y2 θ2 = P0 M0 P1 M1 P2 M2 ,   or:     Kg ug = Fg

  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.

    y 0 = 0
    θ 0 = 0
    y 1 = 0
    θ 1 = 0
    y 2 = 0
    θ 2 = 0
    P 0 = [N]
    M 0 = 0 [Nm]
    P 1 = [N]
    M 1 = 0 [Nm]
    P 2 = [N]
    M 2 = 0 [Nm]


    0 0
    0 0
    0 0
    0 0
    y 0
    θ 0
    y 1
    θ 1
    y 2
    θ 2
    P 0
    M 0
    P 1
    M 1
    P 2
    M 2

  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.