Finite Element Method in Structural Mechanics: Trusses, Beams & Frames

1. The Direct Stiffness Method

Structural FEM predates general continuum FEM by decades — it grew out of matrix analysis of frameworks developed for aircraft design in the 1950s. The core idea is simple: model a structure as a network of one-dimensional members (bars, beams, columns) connected at discrete nodes. Each member has a known relationship between the forces at its ends and the displacements at its ends — its element stiffness matrix. Combine every member's stiffness matrix into one large global stiffness matrix, apply loads and supports, and solve a linear system for the unknown displacements.

This is exactly the same K·u = F system solved by continuum FEM, but here k_e comes from an exact analytical solution of the governing ODE for a bar or beam rather than from numerically integrating shape functions over an arbitrary element shape. That is why hand calculation of small trusses was practical even before computers — and why the direct stiffness method is still the fastest and most transparent way to teach FEM fundamentals.

2. Truss (Bar) Elements

A truss member carries only axial force — tension or compression along its length — and is pinned at both ends, so it transmits no bending moment. In its own local coordinate system (x′ along the bar), a 2-node bar with axial stiffness EA/L has the well-known 1D stiffness matrix:

In 2D or 3D, each node has 2 or 3 translational degrees of freedom, so the local axial stiffness must be rotated into global x-y(-z) coordinates using the direction cosines of the bar. A truss model of a bridge or roof might have a few hundred bars and a few hundred nodes — small enough to solve by hand in the pre-computer era, which is exactly why the method was developed for trusses first.

3. Euler-Bernoulli Beam Elements

Beams resist bending in addition to axial load, so each node needs a transverse displacement and a rotation DOF. The classical Euler-Bernoulli beam element assumes plane sections remain plane and perpendicular to the neutral axis (no shear deformation), giving a cubic (Hermite) displacement field between nodes:

The Hermite cubic shape functions used here interpolate both displacement and slope, guaranteeing continuity of both across element boundaries (C¹ continuity) — necessary because the governing equation EI·v'''' = q is fourth order. Timoshenko beam elements extend this to include transverse shear deformation, important for short, deep beams where shear flexibility is significant (span-to-depth ratio below about 10).

4. Frame Elements & Coordinate Transformation

Real structures — building frames, roof trusses on sloped members, offshore jackets — have members oriented at arbitrary angles. Every element's local stiffness matrix must be transformed into the shared global coordinate system before assembly:

In 3D, frame elements also carry torsion (GJ/L, where G is shear modulus and J the torsion constant) and bending about two orthogonal axes, giving a full 12×12 element stiffness matrix (6 DOF per node: 3 translations, 3 rotations). This is the workhorse element behind almost every structural engineering package — SAP2000, ETABS, STAAD.Pro and Robot Structural Analysis all build their global stiffness matrix from exactly this kind of 12×12 frame element.

5. Assembly & Solving the Global System

Assembly maps each element's local DOFs to their position in the global DOF numbering and adds the transformed element matrix into the global matrix at those locations:

K is sparse and banded — a node only connects to a handful of neighbouring members, so most entries are zero. Renumbering DOFs to minimise bandwidth (Cuthill-McKee, nested dissection) dramatically speeds up direct solvers. Once D_f is known, each element's local displacements are recovered, and its internal axial force, shear and bending moment diagrams follow directly from k_local · d_e.

6. Static Condensation & Substructuring

Large structures are often broken into substructures (a floor, a bridge segment, a repeated truss module) that are each reduced to an equivalent stiffness at just their boundary nodes — the internal DOFs are eliminated algebraically without approximation:

This is the same Guyan/static reduction used to shrink huge FE models for dynamic analysis, and it's how CAD-integrated frame software assembles precomputed "super-elements" for standard components (a precast beam, a bridge pier) without re-deriving their internal behaviour every time.

7. Modal Analysis & Natural Frequencies

Structures don't just carry static loads — they vibrate under wind, seismic excitation, footfall, or machinery. Adding a consistent mass matrix M turns the static problem into a generalized eigenvalue problem:

Solving for the lowest eigenvalues gives the fundamental natural frequencies and mode shapes that dominate a structure's dynamic response. Resonance occurs when an external excitation frequency approaches ω — this is why tall buildings are checked against wind-induced vortex shedding frequencies, and why footbridges (famously London's Millennium Bridge in 2000) are checked against pedestrian pacing frequencies around 2 Hz.

8. Linear Buckling Analysis

Slender compression members can fail by buckling — sudden lateral instability — at loads well below their material yield strength. Linear (eigenvalue) buckling analysis introduces a geometric stiffness matrix K_g that captures how axial load reduces effective bending stiffness:

FEM buckling analysis generalizes Euler's column formula to arbitrary frames, where instability can involve sway of an entire storey rather than bending of a single member. The lowest eigenvalue λ₁ gives the critical load factor; design codes typically require λ₁ well above 1.0 with a safety margin, and nonlinear (large-displacement) analysis is used when post-buckling behaviour or imperfection sensitivity matters.

9. Implementing a 2D Frame Solver

A minimal 2D frame stiffness solver needs only three pieces: element stiffness generation, assembly, and boundary-condition handling. Here is the core of an element stiffness routine in JavaScript, matching the 6×6 matrix derived in Section 4:

function frameElementStiffness(E, A, I, L, angleRad) {
  // local 6x6 stiffness: [u1,v1,th1,u2,v2,th2]
  const k = zeros(6, 6);
  const EA_L  = E * A / L;
  const EI_L3 = E * I / (L ** 3);
  const EI_L2 = E * I / (L ** 2);
  const EI_L  = E * I / L;

  k[0][0] =  EA_L;  k[3][3] =  EA_L;  k[0][3] = -EA_L;  k[3][0] = -EA_L;

  k[1][1] =  12 * EI_L3;  k[4][4] =  12 * EI_L3;
  k[1][4] = -12 * EI_L3;  k[4][1] = -12 * EI_L3;

  k[1][2] = k[2][1] =  6 * EI_L2;
  k[1][5] = k[5][1] =  6 * EI_L2;
  k[2][4] = k[4][2] = -6 * EI_L2;
  k[4][5] = k[5][4] = -6 * EI_L2;

  k[2][2] = 4 * EI_L;  k[5][5] = 4 * EI_L;
  k[2][5] = k[5][2] = 2 * EI_L;

  const c = Math.cos(angleRad), s = Math.sin(angleRad);
  const T = buildTransform(c, s); // 6x6 block-diagonal rotation
  return multiply(transpose(T), multiply(k, T)); // Tᵀ k T
}

function assembleGlobal(elements, nNodes) {
  const ndof = nNodes * 3;
  const K = zeros(ndof, ndof);
  for (const el of elements) {
    const ke = frameElementStiffness(el.E, el.A, el.I, el.L, el.angle);
    const map = [el.n1*3, el.n1*3+1, el.n1*3+2,
                el.n2*3, el.n2*3+1, el.n2*3+2];
    for (let i = 0; i < 6; i++)
      for (let j = 0; j < 6; j++)
        K[map[i]][map[j]] += ke[i][j];
  }
  return K; // then partition into free/supported DOFs and solve K_ff Df = Ff
}

From here, a browser-based visualizer needs three more steps: (1) solve K_ff·D_f = F_f with a sparse or dense linear solver (Gaussian elimination is fine below a few thousand DOFs), (2) recover member end forces from k_local·d_e for each element, and (3) render deflected shapes and bending-moment diagrams by scaling and interpolating the Hermite cubic within each beam element — exactly the interactive workflow used in the structural simulations on this site.