FEM for structures — assembling global stiffness

Overview

The global stiffness matrix assembly is the heart of the finite element method. We take local element stiffness matrices (already computed in the element's own coordinate system) and assemble them into one giant system matrix that represents the entire structure. This transforms a complex continuous structure into a solvable discrete system of linear equations: $\mathbf{K}\mathbf{u} = \mathbf{F}$ .

Why this matters for spacecraft: Every structural analysis—from launch loads to thermal distortion—relies on this assembly step. A solar panel bracket, a satellite bus, or a rocket interstage is meshed into thousands of elements, and we must correctly assemble them to predict deformations, stresses, and failure modes.

[!intuition] The Core Idea

Think of a structure as a network of springs. Each finite element (beam, truss bar, shell) is like a spring with its own stiffness. When you connect them at nodes, forces at shared nodes must balance, and displacements must be compatible (no gaps or overlaps).

The assembly process:

Each element "votes" for the global stiffness at its connected nodes
We add (superimpose) all element contributions at shared nodes
The result: one big matrix where row/column $i$ corresponds to degree-of-freedom (DOF) $i$ globally

Key insight: Elements share nodes → their stiffness matrices overlap at those DOFs → we sum the overlapping entries. This is the direct stiffness method.

[!definition] Key Concepts

Global Stiffness Matrix $\mathbf{K}$ : A square matrix of size $N \times N$ , where $N$ is the total number of DOFs in the mesh. Entry $K_{ij}$ is the force at DOF $i$ due to unit displacement at DOF $j$ (all other DOFs fixed).

Element Stiffness Matrix $\mathbf{k}^{(e)}$ : The stiffness of element $e$ in its local coordinate system. For a 2-node bar element with2 DOFs per node (4 DOFs total), $\mathbf{k}^{(e)}$ is $4 \times 4$ .

Degrees of Freedom (DOF): The independent displacement/rotation variables. A3D node has 6 DOFs: $u_x, u_y, u_z$ (translations) and $\theta_x, \theta_y, \theta_z$ (rotations). We number all DOFs sequentially across the mesh.

Connectivity (Topology): The mapping from element-local node numbers to global node numbers. If element $e$ connects global nodes 5 and 12, we know exactly where to add $\mathbf{k}^{(e)}$ into $\mathbf{K}$ .

[!formula] Derivation from First Principles

Step 1: Energy Principle (Why Assembly Works)

The total strain energy of the structure is the sum of strain energies in all elements:

$U = \sum_{e=1}^{ne} U^{(e)} = \sum_{e=1}^{n_e} \frac{1}{2} \mathbf{u}^{(e)T} \mathbf{k}^{(e)} \mathbf{u}^{(e)}$

where $\mathbf{u}^{(e)}$ is the element displacement vector (extracted from the global displacement vector $\mathbf{u}$ ).

Why this step? Energy must be additive. Each element contributes independently to the total.

Using the Boolean localization matrix $\mathbf{L}^{(e)}$ that extracts element DOFs from global DOFs:

$\mathbf{u}^{(e)} = \mathbf{L}^{(e)} \mathbf{u}$

Substitute:

$U = \frac{1}{2} \sum_{e=1}^{n_e} (\mathbf{L}^{(e)} \mathbf{u})^T \mathbf{k}^{(e)} (\mathbf{L}^{(e)} \mathbf{u}) = \frac{1}{2} \sum_{e=1}^{n_e} \mathbf{u}^T \mathbf{L}^{(e)T} \mathbf{k}^{(e)} \mathbf{L}^{(e)} \mathbf{u}$

Why this step? We factor out the global displacement vector $\mathbf{u}$ to identify the global stiffness.

Since $\mathbf{u}$ is global and common, the total energy becomes:

$U = \frac{1}{2} \mathbf{u}^T \left( \sum_{e=1}^{n_e} \mathbf{L}^{(e)T} \mathbf{k}^{(e)} \mathbf{L}^{(e)} \right) \mathbf{u}$

Identify the global stiffness:

$\boxed{\mathbf{K} = \sum_{e=1}^{n_e} \mathbf{L}^{(e)T} \mathbf{k}^{(e)} \mathbf{L}^{(e)}}$

Why this step? The quadratic form $\frac{1}{2}\mathbf{u}^T \mathbf{K} \mathbf{u}$ defines the global stiffness matrix. The sum of transformed element stiffnesses gives us $\mathbf{K}$ .

Step 2: Practical Assembly Algorithm

The matrix multiplication $\mathbf{L}^{(e)T} \mathbf{k}^{(e)} \mathbf{L}^{(e)}$ is expensive. In practice, we use the scatter operation:

For each element $e$ :

Compute $\mathbf{k}^{(e)}$ (in local coords, transform to global if needed)
Look up the global DOF numbers for element $e$ 's nodes: $[g_1, g_2, \ldots, g_m]$
Add each entry $k^{(e)}_{ij}$ to $K_{g_i, g_j}$

Algorithm in pseudo-code:

Initialize K = zeros(N, N)
for each element e:
    k_e = element_stiffness(e)
    dofs = global_dof_indices(e)  # e.g., [3, 4, 9 10] for a 2-node bar
    for in range(len(dofs)):
        for j in range(len(dofs)):
            K[dofs[i], dofs[j]] += k_e[i, j]  # Scatter-add

Why scatter-add? Because $\mathbf{L}^{(e)T} \mathbf{k}^{(e)} \mathbf{L}^{(e)}$ is mostly zeros—it only populates entries at the DOFs that element $e$ touches. Direct addition is faster than matrix multiplication.

Step 3: Properties of $\mathbf{K}$

Symmetry: If each $\mathbf{k}^{(e)}$ is symmetric (Maxwell-Betti reciprocal theorem), then $\mathbf{K}$ is symmetric.

Sparsity: $\mathbf{K}$ is extremely sparse. Entry $K_{ij} \neq 0$ only if DOF $i$ and DOF $j$ belong to the same element. For a mesh with $N$ DOFs, typical bandwidth is $\mathcal{O}(\sqrt{N})$ in2D, $\mathcal{O}(N^{2/3})$ in 3D.

Positive semi-definite (before BCs): $\mathbf{u}^T \mathbf{K} \mathbf{u} = 2U \geq 0$ . Rigid-body modes (translations/rotations with zero strain) give eigenvalues. After applying boundary conditions (fixing DOFs), $\mathbf{K}$ becomes positive definite → invertible.

[!example] Worked Example 1: Two-Bar Truss

Setup:

2 bar elements in series along $x$ -axis
3 nodes: Node 1 (fixed), Node 2 (middle), Node 3 (free end, force $F$ )
Each bar: length $L$ , area $A$ , Young's modulus $E$
Each node has 1 DOF (horizontal displacement $u$ )

Element stiffness (1D bar):

$\mathbf{k}^{(e)} = \frac{EA}{L} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}$

Why this formula? From Hooke's law: force $= \frac{EA}{L}(\Delta L) = \frac{EA}{L}(u_j - u_i)$ . Enforcing equilibrium at both nodes gives the $2 \times 2$ stiffness.

Global DOF numbering:

Node 1 → DOF 1
Node 2 → DOF 2
Node 3 → DOF 3

Element 1 (Nodes 1-2): $\mathbf{k}^{(1)} = \frac{EA}{L} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}, \quad \text{DOFs: } [1, 2]$

Scatter into $\mathbf{K}$ (3×3): $\mathbf{K} \leftarrow \mathbf{K} + \frac{EA}{L} \begin{bmatrix} 1 & -1 & 0 \\ -1 & 1 & 0 \\ 0 & 0 \end{bmatrix}$

Why this step? Element 1 connects DOFs 1 and 2. We place its $2\times 2$ block in rows/columns 1-2 of the global $3\times 3$ matrix.

Element 2 (Nodes 2-3): $\mathbf{k}^{(2)} = \frac{EA}{L} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}, \quad \text{DOFs: } [2, 3]$

Scatter: $\mathbf{K} \leftarrow \mathbf{K} + \frac{EA}{L} \begin{bmatrix} 0 & 0 & 0 \\ 0 & 1 & -1 \\ 0 & -1 & 1 \end{bmatrix}$

Assembled global stiffness: $\mathbf{K} = \frac{EA}{L} \begin{bmatrix} 1 & -1 & 0 \\ -1 & 2 & -1 \\ 0 & -1 & 1 \end{bmatrix}$

Why does $K_{22} = 2$ ? Node 2 is shared by both elements. Both contribute $\frac{EA}{L} \cdot 1$ to $K_{22}$ , so we sum: $1 + 1 = 2$ .

Apply boundary condition: Fix Node 1 → $u_1 = 0$ → remove row/column 1:

$\mathbf{K}_{\text{reduced}} = \frac{EA}{L} \begin{bmatrix} 2 & -1 \ -1 & 1 \end{bmatrix}, \quad \mathbf{F}_{\text{reduced}} = \begin{bmatrix} 0 \ F \end{bmatrix}$

Solve: $\begin{bmatrix} 2 & -1 \\ -1 & 1 \end{bmatrix} \begin{bmatrix} u_2 \\ u_3 \end{bmatrix} = \frac{L}{EA} \begin{bmatrix} 0 \\ F \end{bmatrix}$

From the second equation: $-u_2 + u_3 = \frac{FL}{EA}$ . From the first: $2u_2 - u_3 = 0 \Rightarrow u_3 = 2u_2$ .

Substitute: $-u_2 + 2u_2 = \frac{FL}{EA} \Rightarrow u_2 = \frac{FL}{EA}$ . Then $u_3 = \frac{2FL}{EA}$ .

Physical check: Total elongation is $u_3 = \frac{2FL}{EA}$ , which matches two bars in series: $\Delta L_{\text{total}} = \frac{FL}{EA} + \frac{FL}{EA} = \frac{2FL}{EA}$ . ✓

[!example] Worked Example 2: D Truss with Coordinate Transformation

Setup:

Two bar elements forming a 90° corner
Node 1 at origin (fixed), Node 2 at $(L, 0)$ , Node 3 at $(L, L)$
Element 1: horizontal (1→2), Element 2: vertical (2→3)
Each node has 2 DOFs: $(u_x, u_y)$
Global DOF numbering: Node $i$ has DOFs $2i-1$ (x) and $2i$ (y)

Element stiffness in local coordinates (aligned with bar axis $\xi$ ):

For a bar along its local axis: $\mathbf{k}_{\text{local}} = \frac{EA}{L} \begin{bmatrix} 1 & 0 & -1 & 0 \\ 0 & 0 & 0 \\ -1 & 0 & 1 & 0 \\ 0 & 0 & 0 \end{bmatrix}$

Why these zeros? A bar only resists axial deformation. Transverse displacements (perpendicular to the bar) produce zero strain, hence zero stiffness.

Transformation matrix $\mathbf{T}$ : Converts global displacements $(u_x, u_y)$ to local $(u_\xi, u_\eta)$ :

$\mathbf{T} = \begin{bmatrix} \cos\theta & \sin\theta & 0 & 0 \\ -\sin\theta & \cos\theta & 0 & 0 \ 0 & 0 & \cos\theta & \sin\theta \\ 0 & 0 & -\sin\theta & \cos\theta \end{bmatrix}$

where $\theta$ is the angle of the bar axis from the global $x$ -axis.

Global element stiffness: $\mathbf{k}_{\text{global}} = \mathbf{T}^T \mathbf{k}_{\text{local}} \mathbf{T}$

Element 1 (horizontal, $\theta = 0°$ ): $\mathbf{T}^{(1)} = \mathbf{I} \quad (\text{already aligned})$ $\mathbf{k}^{(1)}_{\text{global}} = \frac{EA}{L} \begin{bmatrix} 1 & 0 & -1 & 0 \\ 0 & 0 & 0 \\ -1 & 0 & 1 & 0 \\ 0 & 0 & 0 \end{bmatrix}$

DOFs: $[1, 2, 3, 4]$ (Node 1: DOFs 1,2; Node 2: DOFs 3,4).

Element 2 (vertical, $\theta = 90°$ ): $\mathbf{T}^{(2)} = \begin{bmatrix} 0 & 1 & 0 & 0 \\ -1 & 0 & 0 \\ 0 & 0 & 1 \\ 0 & -1 & 0 \end{bmatrix}$

After transformation: $\mathbf{k}^{(2)}_{\text{global}} = \frac{EA}{L} \begin{bmatrix} 0 & 0 \\ 0 & 1 & 0 & -1 \\ 0 & 0 \\ 0 & -1 & 0 & 1 \end{bmatrix}$

DOFs: $[3, 4, 5, 6]$ (Node 2: DOFs 3,4; Node 3: DOFs 5,6).

Assembly: Initialize $\mathbf{K}$ as $6 \times 6$ zeros.

Scatter $\mathbf{k}^{(1)}$ into rows/cols $[1,2,3,4]$ : $\mathbf{K}_{1:4,1:4} \leftarrow \mathbf{k}^{(1)}$

Scatter $\mathbf{k}^{(2)}$ into rows/cols $[3,4,5,6]$ : $\mathbf{K}_{3:6,3:6} \leftarrow \mathbf{K}_{3:6,3:6} + \mathbf{k}^{(2)}$

Result (before BCs):