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

In 1956, Boeing engineer M. J. Turner stalled while computing swept-wing stresses by hand. Hundreds of members cannot be balanced one equation at a time. The direct stiffness method was the way out he and his colleagues published. Form an element stiffness, add it where degrees of freedom overlap, and solve $Ku = f$ once.

Start with the angle. A bar laid flat kills the vertical slots. Drop the angle slider to 0 and the vertical rows and columns of the 4×4 stiffness go to zero together.

Watch which entries empty at $\theta = 0$. The bar is stiff only along its own axis. Raise the angle and those zeros split between $x$ and $y$.

Stiffness Starts With a Single Spring#

A bar element is a spring that carries force only along its axis. Length $L$, area $A$, modulus $E$, and the axial stiffness is $k = EA/L$. Force per unit stretch. With end displacements $u_1, u_2$, the nodal forces are:

$f_1, f_2$ are the end forces, and the bracket is the local stiffness $k_e$. Each row sums to zero. A rigid shift of both nodes costs no force. Leave the supports off and the matrix is singular for that reason.

From Local to Global: The Rotation#

Members point every which way. Assembly needs the local axial displacement written in global $x, y$. Let $\theta$ be the angle with the $x$-axis, $c = \cos\theta$, $s = \sin\theta$. Axial displacement is a projection: $u_{\text{axial}} = c\,u_x + s\,u_y$. Apply that map on both sides and the global element stiffness is 4×4.

Each entry ties node 1's $x, y$ to node 2's $x, y$. At $\theta = 0$, $s = 0$ and $c = 1$, so the matrix times $EA/L$ has a first row of 1, 0, −1, 0 and a second row of zeros. The third row flips the sign of the first. The fourth row is zero too. Those zeros are why each 2×2 block has rank one. Tilt the coordinates and the stiffness tensor turns with them. That turn is spelled out in rotating a constitutive tensor with the metric.

Adding Overlapping DOFs: Global Assembly#

Assembly is addition at shared degrees of freedom. Node $i$ sits at global indices $2i, 2i+1$, $x$ then $y$. Scatter the element's four local DOFs onto those indices and accumulate.

$\mathbf{L}_e$ is the selection matrix from element DOFs to global DOFs, and $K$ is the global stiffness. Code never builds $\mathbf{L}_e$. It adds through the index list. Several members on one node stack inside that diagonal block.

Boundary Conditions: Erase the Supports#

A freshly assembled $K$ is singular. The structure still floats, so rigid motion is free. Displacements at the supports have to be held at zero. Split the DOFs into a free set $f$ and a constrained set $c$.

With $u_c = 0$, only the top block remains: $K_{ff}\,u_f = F_f$. That reduced system is nonsingular. Reactions come back later as $R_c = K_{cf}\,u_f$. $R_c$ is the force the support returns to the structure.

Raise the load slider and EA together and the tip drop changes in front of you.

A larger EA cuts the sag under the same load. The stiffness matrix is the map from load to displacement. Hold the load fixed and double EA. The tip displacement should head toward half.

Python: Solving a 12-Member Truss#

A wall-mounted three-panel cantilever. Eight nodes, twelve members. Coordinates, connectivity, and loads go in. Displacements and axial forces come out.

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_ drops each element stiffness into the right global slots. That loop is the same shape as a commercial FEA assembly.

Node spacing is 1 m. Nine horizontal and vertical members have $L = 1$. With the code's EA of $8.0\times 10^6$ N, $k = EA/L = 8.0\times 10^6$ N/m. One tip node's vertical load is $12\times 10^3$ N, so that bar alone moves $1.5\times 10^{-3}$ m, which is 1.5 mm. The truss tip is not that division. Members share the load, and the matrix that actually gets solved is the reduced $K_{ff}$.

Three members are diagonals. The one from (0, 1) to (1, 0) has length $\sqrt{2}$, so its stiffness is $1/\sqrt{2}$ times a horizontal bar. Same EA, softer diagonal. That member sits at 45°. Then $c^2$ and $s^2$ are $1/2$, and the absolute value of $cs$ is $1/2$ too. Every nonzero slot in the 4×4 is half of $EA/L$. Length $\sqrt{2}$ makes the slot size $EA/(2\sqrt{2})$. A horizontal bar leaves the vertical slots empty. A diagonal fills all four directions at that size. The same EA splits that way from the angle alone. Eight nodes make 16 DOFs. Wall nodes 0 and 1 restrain 4 of them, so the solve is 12×12. A 2D model needs at least 3 restraints. This one uses 4. Miss the $y$ restraint on one wall node and the count sits on that minimum. Miss both and $K_{ff}$ is singular again.

When the Stiffness Matrix Goes Singular#

The error people actually hit is a singular matrix. Three ways it shows up.

Skip a support and a rigid mode remains, so $K_{ff}$ stays singular. Two dimensions need at least 3 DOFs held. Three dimensions need 6.

A quad panel that was never split into triangles is a mechanism. A truss has to be filled with triangles or it flops.

Zero length or stacked nodes make $L = 0$ and the division blows up. Coincident coordinates have to be filtered before a merge.

Miss one wall node and the assembly stops there#

The 1.5 mm figure belongs to a single horizontal bar. Tip sag of the 12-member cantilever comes from solving $K_{ff}\,u_f = F_f$, and that equation is not even defined until the four wall DOFs are removed. One missing support leaves a rigid mode, and the displacement is a null space rather than a number.

Swap the bars for surface elements and the same assembly can freeze for a different reason. How that freeze shows up in a thin shell is in shear locking, and how MITC puts the stiffness back. Check the four supports on the truss first. Open that page only when the model gets thin.