The direct stiffness method

How a finite element solver turns a model into K u = f and back — stiffness, transformation, assembly, constraints, solution and recovery.

Theory notes 6 min read

Every result in Elementarium comes out of one procedure. This page walks through it once, with enough algebra to follow but not to implement, so that when a number surprises you the surprise has somewhere to go.

The idea

A structure is a set of elements joined at nodes. Each element resists the relative motion of its nodes with forces proportional to that motion (that is what "linear elastic" means), and the constant of proportionality is a matrix. Add up every element's contribution at every node and you have the structure's stiffness matrix $K$ , which relates every nodal displacement $u$ to every nodal force $f$ :

Ku = f

Given the loads, solve for the displacements. Given the displacements, every internal force follows. Everything else is bookkeeping.

1. Element stiffness in local axes

Each element type has a stiffness matrix written in its own axes, derived from its formulation. For a bar of length $L$ , area $A$ and modulus $E$ , it is 2 × 2:

k = \frac{E A}{L} [1 - 1 - 1 1]

Push one end towards the other by $δ$ and the bar pushes back with $E A δ / L$ . A beam's matrix is 12 × 12 (six degrees of freedom at each of two nodes), and its bending terms come from cubic interpolation, which is why a prismatic beam under uniform load is exact without subdivision. A shell element's is 24 × 24 and is integrated numerically at Gauss points.

The familiar beam terms are in there. The moment needed to rotate one end of a beam by $θ$ while the far end is fixed is $4 E I θ / L$ , and the moment it induces at the far end is $2 E I θ / L$ . With both ends held against translation, the end rotations and moments are related by

[M_{1} M_{2}] = \frac{E I}{L} [4224] [θ_{1} θ_{2}]

These are the $4 E I / L$ and $2 E I / L$ of moment distribution, which is the same method by hand.

2. Transformation to global axes

The element's axes point along the element; the structure's do not. Each element matrix is rotated into the global system with a transformation matrix $T$ built from the direction cosines of the member:

\overset{ˉ}{k} = T^{T} k T

This is where the local axis convention enters: a member drawn from A to B and one drawn from B to A have different $T$ , report their internal forces with opposite signs, and contribute identically to $K$ .

3. Assembly

Each global element matrix is added into $K$ at the rows and columns belonging to its nodes' degrees of freedom:

K = e \sum A_{e}^{T} \overset{ˉ}{k}_{e} A_{e}

where $A_{e}$ only picks out element $e$ 's degrees of freedom from the global list. A node shared by three members receives three contributions in the same rows. That is what "connected" means numerically, and it is why two members whose ends are 3 mm apart are not connected: they have no rows in common.

Assembly for a portal frame. Each member adds its stiffness to the four node blocks it couples; where two members meet at a node their contributions sum on the diagonal. Node pairs no member joins stay zero, which is the sparsity the solver exploits. The fixed bases' rows and columns are removed before solving.

$K$ is sparse: a node is only coupled to the nodes it shares an element with, so nearly every entry is zero. It is also symmetric, by Maxwell's reciprocal theorem. The solvers exploit both; see Choosing a linear solver.

Loads are assembled into $f$ the same way. A nodal load goes straight in. A distributed load on a member is converted to equivalent nodal forces: for a uniform $w$ on a beam, $w L /2$ at each end and fixed-end moments $w L^{2} /12$ . The element's own end-force treatment later corrects for these, so the moment along the member is exact.

4. Constraints

Supports remove degrees of freedom. A restrained one has a prescribed displacement (zero, or a settlement), so its row and column leave the system and its reaction is recovered afterwards from the equilibrium of that row. Splitting the degrees of freedom into free ( $f$ ) and supported ( $s$ ):

[K_{f f} K_{s f} K_{f s} K_{ss}] [u_{f} u_{s}] = [f_{f} f_{s} + r]

The solve uses only the top row, with the prescribed $u_{s}$ moved to the right; the reactions $r$ come from the bottom row afterwards:

K_{f f} u_{f} r = f_{f} - K_{f s} u_{s} = K_{s f} u_{f} + K_{ss} u_{s} - f_{s}

A spring does not remove the degree of freedom; it adds its stiffness to the diagonal.

Releases modify an element before assembly. A full release decouples that degree of freedom from the node, so the member carries no moment across that end; a partial release puts a finite stiffness between the two, which is the semi-rigid connection.

Rigid arms are constraints that tie one node's degrees of freedom to another's by a kinematic relation, applied at assembly instead of as a very stiff element. That is why they do not damage conditioning the way a stiff element does.

After this step, $K$ is singular if and only if the structure can move without resistance: a mechanism. Elementarium checks for that before solving, and names the free motion.

5. Solution

$Ku = f$ is solved once per load state. With a direct solver $K$ is factorised, as $L L^{T}$ or $LD L^{T}$ , and each load case or combination is a cheap resolve against the factors. That is why fifty combinations cost barely more than one. With an iterative solver the system is converged separately for each right-hand side. See Choosing a linear solver.

The result is $u$ : every translation and rotation at every node, for that load state.

6. Recovery

Internal forces come from the displacements, element by element: transform the element's nodal displacements back to its local axes, multiply by its local stiffness, and add back the fixed-end forces $f_{0}$ from any load along its length:

f_{e} = k T u_{e} + f_{0}

For a beam that gives $N$ , $V$ , $T$ and $M$ at the ends. Between the ends they follow by statics from the end forces and the load on the member, which is why the moment under a uniform load comes out parabolic although the shape function is cubic. For a shell, stress resultants are evaluated at the Gauss points and carried to the corners by the recovery method the analysis setup chooses.

Reactions are the forces at the restrained degrees of freedom, and their sum must equal the applied load: the check described in Your first analysis.

Where the other analyses depart

Analysis	What changes
Modal	A mass matrix $M$ is assembled the same way; $(K - ω^{2} M) ϕ = 0$ is an eigenvalue problem, not a linear system
Linear buckling	A geometric stiffness $K_{g}$ is assembled from the axial forces of a static solve; $(K + α K_{g}) ϕ = 0$
Second-order	$K_{g} (u)$ depends on the answer, so steps 5 and 6 repeat until $u$ stops changing

All three assemble their extra matrix from the same elements, with the same transformation and the same connectivity. That is why a modelling error in step 3, a member that does not connect, appears in every analysis type at once.

What the method assumes

Linearity. $K$ does not depend on $u$ ; the material stays elastic and the geometry stays where it was drawn. The exceptions are the three rows above.
Small displacements. Equilibrium is written on the undeformed shape, except in second-order analysis.
Compatibility at nodes only. Elements are joined at their nodes and nowhere else. Two surfaces that overlap without sharing nodes pass through each other.
The element's shape functions. A beam is exact; a shell is an approximation that converges with refinement. Meshing is the consequence.

The method is the one taught in every structural analysis course and implemented in every finite element package.

Something wrong or unclear on this page? Tell us, and it gets corrected.