Why Does the Bridge Sag Like That? The Direct Stiffness Method for Trusses#

In 1956, Boeing engineer M. J. Turner hit a wall computing swept-wing stresses by hand. You cannot chase equilibrium through a structure of hundreds of members one equation at a time. The answer he and his colleagues published was the direct stiffness method. Build the stiffness matrix of a single element, add contributions at shared degrees of freedom into one big matrix, and solve $Ku = f$ once.

This post takes a 2D truss and codes the whole chain in Python: derive one bar element's stiffness matrix, rotate it from local to global coordinates, assemble the global matrix, apply boundary conditions, and solve for displacements. Once you see that the skeleton of the finite element method is really one line of linear algebra, the structural part of an FSI solver and the assembly loop of a commercial FEA code read the same way.

Stiffness Starts With a Single Spring#

A bar element is a spring that only carries force along its axis. For a bar of length $L$, cross-section $A$, and modulus $E$, the axial stiffness is $k = EA/L$ (force per unit stretch). With end-node axial displacements $u_1, u_2$, the nodal forces are:

Here $f_1, f_2$ are the end forces and the bracket is the local stiffness matrix $k_e$. Each row sums to zero, meaning a rigid-body motion (both nodes sliding together) costs no force. That singularity is exactly why boundary conditions become mandatory later.

From Local to Global: The Rotation#

Truss members point every which way. To assemble them, the local axial displacement has to become global $x, y$ components. Let $\theta$ be the angle the member makes with the $x$-axis, with $c = \cos\theta$ and $s = \sin\theta$. The axial displacement is a projection of the global one: $u_{\text{axial}} = c\,u_x + s\,u_y$. Multiplying this transformation onto both sides of the stiffness matrix gives the 4×4 global element stiffness.

Each entry is the stiffness linking node 1's $x, y$ to node 2's $x, y$ degrees of freedom. Rotate the bar yourself in the simulation below.

At $\theta = 0$ we have $s = 0$, so the vertical-DOF rows and columns go entirely to zero. The bar is stiff only along its own axis. That is why each 2×2 block has rank one.

Adding Overlapping DOFs: Global Assembly#

Assembly sounds grand, but it is just adding stiffness at shared degrees of freedom. Node $i$'s DOFs live at global indices $2i, 2i+1$ (its $x, y$). Scatter each element's four local DOFs to these global indices, then accumulate the element stiffness there.

$\mathbf{L}_e$ is the selection matrix that sends element DOFs to global DOFs, and $K$ is the global stiffness matrix. Real code never forms $\mathbf{L}_e$; it adds directly through an index array. Where several members meet at one node, their contributions stack in that diagonal block.

Boundary Conditions: Erase the Supports#

The freshly assembled $K$ is singular. The structure still floats freely in space, so rigid-body motion is unconstrained. You have to pin displacements to zero at the supports before it can be solved. The cleanest way is to split the DOFs into a free set $f$ and a constrained set $c$.

Since $u_c = 0$ (fixed), keeping only the top row gives $K_{ff}\,u_f = F_f$. This reduced system is non-singular and solves. The reactions $R_c$ are recovered afterward as $R_c = K_{cf}\,u_f$. Change the load and the stiffness in the truss below.

Raise EA and the same load bends the truss far less. The stiffness matrix is the load-to-displacement map.

Python: Solving a 12-Member Truss#

Assemble and solve a wall-mounted, three-panel cantilever truss (8 nodes, 12 members) with numpy. The inputs are node coordinates, member connectivity, and loads; the outputs are nodal displacements and member axial forces.

import numpy as np
 
nodes = np.array([[0,0],[0,1],[1,0],[1,1],[2,0],[2,1],[3,0],[3,1]], float)
members = [(0,2),(2,4),(4,6),(1,3),(3,5),(5,7),
           (2,3),(4,5),(6,7),(1,2),(3,4),(5,6)]
EA = 8.0e6           # axial rigidity EA [N]
fixed = [0, 1]       # nodes pinned to the wall
 
def bar_stiffness(p1, p2, EA):
    d = p2 - p1
    L = np.hypot(*d)
    c, s = d / L
    k = EA / L * np.array([[ c*c,  c*s, -c*c, -c*s],
                           [ c*s,  s*s, -c*s, -s*s],
                           [-c*c, -c*s,  c*c,  c*s],
                           [-c*s, -s*s,  c*s,  s*s]])
    return k, L, (c, s)
 
ndof = nodes.shape[0] * 2
K = np.zeros((ndof, ndof))
geom = []
for a, b in members:                       # assemble global stiffness
    k, L, cs = bar_stiffness(nodes[a], nodes[b], EA)
    dof = [2*a, 2*a+1, 2*b, 2*b+1]
    K[np.ix_(dof, dof)] += k
    geom.append((L, cs))
 
F = np.zeros(ndof)                          # 24 kN load at free end (nodes 6,7)
for n in (6, 7):
    F[2*n+1] -= 12.0e3
 
fixed_dof = [d for n in fixed for d in (2*n, 2*n+1)]
free_dof = [d for d in range(ndof) if d not in fixed_dof]
 
u = np.zeros(ndof)                          # reduced system K_ff u_f = F_f
u[free_dof] = np.linalg.solve(K[np.ix_(free_dof, free_dof)],
                              F[free_dof])
 
for m, (a, b) in enumerate(members):        # axial force N = (EA/L)[-c,-s,c,s]·u
    L, (c, s) = geom[m]
    dof = [2*a, 2*a+1, 2*b, 2*b+1]
    N = EA / L * np.array([-c, -s, c, s]) @ u[dof]
    print(f"member {a}-{b}: N = {N/1e3:+7.2f} kN "
          f"({'tension' if N > 0 else 'compression'})")
 
print(f"free-end drop = {u[2*6+1]*1e3:.3f} mm")

np.ix_ builds the index grid that adds each element stiffness at the right global slots. These twenty lines are structurally identical to the heart of a commercial FEA solver.

When the Stiffness Matrix Goes Singular#

The most common failure in the field is a "singular matrix" error. The cause is usually one of three things.

First, missing boundary conditions. Too few supports to block rigid-body motion, and $K_{ff}$ stays singular. In 2D you must constrain at least 3 DOFs, in 3D at least 6.

Second, a mechanism. A quad panel that was never triangulated has too few members and flops. Always fill a truss with triangles.

Third, zero-length or duplicate nodes. Two nodes at the same coordinate make $L = 0$ blow up the division. Filter coincident coordinates before merging a mesh.

The Takeaway#

• The skeleton of the finite element method is: build element stiffness, rotate to global, add at shared DOFs, apply boundary conditions, and solve $Ku = f$.
• Before assembly with constraints, $K$ is always singular. The supports that block rigid-body motion are what make it invertible.
• A bar element's 2×2 block has rank one: it is stiff only along its axis. The angle transformation scatters that stiffness into global coordinates.