Skip to main content
Engineering LibreTexts

30.3: Direct Stiffness Method and the Global Stiffness Matrix

  • Page ID
    32862
  • Although there are several finite element methods, we analyse the Direct Stiffness Method here, since it is a good starting point for understanding the finite element formulation. We consider first the simplest possible element – a 1-dimensional elastic spring which can accommodate only tensile and compressive forces. For the spring system shown, we accept the following conditions:

    • Condition of Compatibility – connected ends (nodes) of adjacent springs have the same displacements
    • Condition of Static Equilibrium – the resultant force at each node is zero
    • Constitutive Relation – that describes how the material (spring) responds to the applied loads

    alt
    Model spring system

    The constitutive relation can be obtained from the governing equation for an elastic bar loaded axially along its length:

    \[ \frac{d}{du} (AE \frac{\Delta l}{l_0}) + k = 0 \]

    \[ \frac{\Delta l}{l_0} = \varepsilon \]

    \[ \frac{d}{du} (AE \varepsilon) + k = 0 \]

    \[ \frac{d}{du} (A \sigma) + k = 0 \]

    \[ \frac{dF}{du} + k = 0 \]

    \[ \frac{dF}{du} = -k \]

    \[ dF = -kdu \]

    The spring stiffness equation relates the nodal displacements to the applied forces via the spring (element) stiffness. The minus sign denotes that the force is a restoring one, but from here on in we use the scalar version of Eqn.7.

    Derivation of the Stiffness Matrix for a Single Spring Element

    single spring element

    From inspection, we can see that there are two degrees of freedom in this model, ui and uj. We can write the force equilibrium equations:

    \[ k^{(e)}u_i - k^{(e)}u_j = F^{(e)}_{i} \]

    \[ -k^{(e)}u_i + k^{(e)}u_j = F^{(e)}_{j} \]

    In matrix form

    \[ \begin{bmatrix}
    k^{e} & -k^{e} \\
    -k^{e} & k^{e}
    \end{bmatrix}\begin{Bmatrix}
    u_i\\
    u_j
    \end{Bmatrix} =
    \begin{Bmatrix}
    F^{(e)}_i\\
    F^{(e)}_j
    \end{Bmatrix} \]

    The order of the matrix is [2×2] because there are 2 degrees of freedom. Note also that the matrix is symmetrical. The ‘element’ stiffness relation is:

    \[ [K^{(e)}] \begin{bmatrix} u^{(e)} \end{bmatrix} = \begin{bmatrix} F^{(e)} \end{bmatrix} \]

    Where Κ(e) is the element stiffness matrix, u(e) the nodal displacement vector and F(e) the nodal force vector. (The element stiffness relation is important because it can be used as a building block for more complex systems. An example of this is provided later.)

    Derivation of a Global Stiffness Matrix

    For a more complex spring system, a ‘global’ stiffness matrix is required – i.e. one that describes the behaviour of the complete system, and not just the individual springs.

    Complex spring system

    From inspection, we can see that there are two springs (elements) and three degrees of freedom in this model, u1, u2 and u3. As with the single spring model above, we can write the force equilibrium equations:

    \[ k^1u_1 - k^1u_2 = F_1 \]

    \[ -k^1u_1 + (k^1 + k^2)u_2 - k^2u_3 = F_2 \]

    \[ k^2u_3 - k^2u_2 = F_3 \]

    In matrix form

    \[ \begin{bmatrix}
    k^1 & -k^1 & 0\\
    -k^1 & k^1 + k^2 & -k^2\\
    & -k^2 & k^2
    \end{bmatrix}
    \begin{Bmatrix}
    u_1\\
    u_2\\
    u_3
    \end{Bmatrix}
    =
    \begin{Bmatrix}
    F_1\\
    F_2\\
    F_3
    \end{Bmatrix} \]

    The ‘global’ stiffness relation is written in Eqn.16, which we distinguish from the ‘element’ stiffness relation in Eqn.11.

    \[ \begin{bmatrix}
    K
    \end{bmatrix}
    \begin{Bmatrix}
    u
    \end{Bmatrix}
    =
    \begin{Bmatrix}
    F
    \end{Bmatrix} \]

    Note the shared k1 and k2 at k22 because of the compatibility condition at u2. We return to this important feature later on.

    Assembling the Global Stiffness Matrix from the Element Stiffness Matrices

    Although it isn’t apparent for the simple two-spring model above, generating the global stiffness matrix (directly) for a complex system of springs is impractical. A more efficient method involves the assembly of the individual element stiffness matrices. For instance, if you take the 2-element spring system shown,

    Complex spring system

    split it into its component parts in the following way

    Complex spring part a Complex spring part b

    and derive the force equilibrium equations

    \[ k^1u_1 - k^1u_2 = F_1 \]

    \[ k^1u_2 - k^1u_1 = k^2u_2 - k^2u_3 = F_2 \]

    \[ k^2u_3 - k^2u_2 = F_3 \]

    then the individual element stiffness matrices are:

    \[ \begin{bmatrix}
    k^1 & -k^1 \\ k^1 & k^1 \end{bmatrix}
    \begin{Bmatrix} u_1\\ u_2 \end{Bmatrix}
    =
    \begin{Bmatrix} F_1\\ F_2 \end{Bmatrix} \]

    and

    \[ \begin{bmatrix} k^2 & -k^2 \\ k^2 & k^2 \end{bmatrix}

    \begin{Bmatrix} u_2\\ u_3 \end{Bmatrix}

    =

    \begin{Bmatrix} F_2\\ F_3 \end{Bmatrix} \]

    such that the global stiffness matrix is the same as that derived directly in Eqn.15:

    Global matrix

    (Note that, to create the global stiffness matrix by assembling the element stiffness matrices, k22 is given by the sum of the direct stiffnesses acting on node 2 – which is the compatibility criterion. Note also that the indirect cells kij are either zero (no load transfer between nodes i and j), or negative to indicate a reaction force.)

    For this simple case the benefits of assembling the element stiffness matrices (as opposed to deriving the global stiffness matrix directly) aren’t immediately obvious. We consider therefore the following (more complex) system which contains 5 springs (elements) and 5 degrees of freedom (problems of practical interest can have tens or hundreds of thousands of degrees of freedom (and more!)). Since there are 5 degrees of freedom we know the matrix order is 5×5. We also know that it’s symmetrical, so it takes the form shown below:

    Spring system with 5 degrees of freedom

    We want to populate the cells to generate the global stiffness matrix. From our observation of simpler systems, e.g. the two spring system above, the following rules emerge:

    • The term in location ii consists of the sum of the direct stiffnesses of all the elements meeting at node i
    • The term in location ij consists of the sum of the indirect stiffnesses relating to nodes i and j of all the elements joining node i to j
    • Add a negative for reaction terms (–kij)
    • Add a zero for node combinations that don’t interact

    By following these rules, we can generate the global stiffness matrix:

    Global stiffness matrix

    This type of assembly process is handled automatically by commercial FEM codes

    Drag the springs into position and click 'Build matrix', then apply a force to node 5. You will then see the force equilibrium equations, the equivalent spring stiffness and the displacement at node 5.

    Solving for (u)

    The unknowns (degrees of freedom) in the spring systems presented are the displacements uij. Our global system of equations takes the following form:

    Global system of equations

    To find {u} solve

    \[ u = F[K]^{-1} \]

    Recall that

    \[ [k][k]^{-1} = I = Identity Matrix = \begin{bmatrix} 1 & 0\\ 0 & 1\end{bmatrix}\]

    Recall also that, in order for a matrix to have an inverse, its determinant must be non-zero. If the determinant is zero, the matrix is said to be singular and no unique solution for Eqn.22 exists. For instance, consider once more the following spring system:

    Complex spring system

    We know that the global stiffness matrix takes the following form

    \[ \begin{bmatrix}
    k^1 & -k^1 & 0\\
    -k^1 & k^1+k^2 & -k^2\\
    0 & -k^2 & k^2
    \end{bmatrix}
    \begin{Bmatrix}
    u_1\\
    u_2\\
    u_3
    \end{Bmatrix}
    =
    \begin{Bmatrix}
    F_1\\
    F_2\\
    F_3
    \end{Bmatrix} \]

    The determinant of [K] can be found from:

    \[ det
    \begin{bmatrix}
    a & b & c\\
    d & e & f\\
    g & h & i
    \end{bmatrix}
    = (aei + bfg + cdh) - (ceg + bdi +afh) \]

    Such that:

    \[ (k^1(k^1+k^2)k^2 + 0 + 0) - (0 + (-k^1-k^1k^2) + (k^1 - k^2 - k^3)) \]

    \[ det[K] = ({k^1}^2k^2 + k^1{k^2}^2) - ({k^1}^2k^2 + k^1{k^2}^2) = 0 \]

    Since the determinant of [K] is zero it is not invertible, but singular. There are no unique solutions and {u} cannot be found. If this is the case in your own model, then you are likely to receive an error message!