Matrix stiffness method II: beams

A truss member resists only stretching. A beam element resists bending, so each of its nodes has a transverse displacement and a rotation. This lesson builds the 4 × 4 beam element stiffness matrix, shows how loads applied between nodes are converted into equivalent joint loads, and solves the continuous beams of lessons 4 and 5 by the stiffness method. The result is a general beam program that handles any number of spans, fixed and pinned supports, and support settlements, all with the same few lines of code.

Degrees of freedom for a beam

For a straight beam along the x axis, ignoring axial deformation, each node has two degrees of freedom: a vertical displacement v (positive upward) and a rotation \(\theta\) (positive counterclockwise). Matrix methods conventionally use this right-handed convention rather than the clockwise-positive convention of slope-deflection. The numbers are identical, only the signs differ, so be careful when comparing results: a counterclockwise rotation of \(20.87/EI\) in this lesson is the \(\theta_B = -20.87/EI\) of lesson 5.

Nodes are placed at every support, at every change of section, and at any point where you want displacements reported. Supports restrain degrees of freedom: a pin or roller restrains v; a fixed support restrains both v and \(\theta\); a free end or an interior node restrains neither. An internal hinge is modeled either by giving the element a released end (a modified element matrix) or by two nodes at the same location sharing v but not \(\theta\).

The beam element stiffness matrix

Column j of the element matrix holds the end forces needed to produce a unit value of degree of freedom j with the other three held at zero. Each of these four shapes is a fixed-ended beam with one end displaced or rotated, so the end forces come straight from the slope-deflection equation of lesson 5 and from equilibrium of the element. A unit rotation at i needs a moment \(4EI/L\) there, induces \(2EI/L\) at j, and needs shears of \(\pm 6EI/L^2\) to keep the element in equilibrium. A unit transverse displacement at i, with no rotation, needs end moments of \(6EI/L^2\) and shears of \(\pm 12EI/L^3\). Collecting all four columns:

The matrix is symmetric and singular: it has two rigid-body modes, a uniform translation and a small rigid rotation, which produce no end forces. A beam element is exact for a prismatic member under end loads, because the cubic deflected shape assumed in its derivation is the true elastic curve when there is no load between the nodes. That is a property worth remembering: one element per span gives exact end moments, with no mesh refinement needed.

Loads between nodes

The stiffness method applies loads at nodes. Loads along an element are converted by the fixed-end idea used in slope-deflection: first lock the nodes and compute the fixed-end forces \(\mathbf{q}_F\) that the locked ends would supply; then release the nodes by applying the opposite of these forces as equivalent joint loads.

The equivalent joint loads, \(-\mathbf{q}_F\), are added into the structure load vector. After solving for the displacements, the true end forces of each element are the stiffness response plus the fixed-end forces:

Forgetting to add \(\mathbf{q}_F\) back is one of the most common errors in hand-coded stiffness programs: the joint displacements come out right, but the end moments of loaded members are wrong by exactly the fixed-end moments.

The office floor beam again

The correspondence is not a coincidence. Slope-deflection is the stiffness method with the free degrees of freedom chosen by hand and the support translations eliminated in advance. The matrix version keeps every degree of freedom and lets the partitioning step remove the restrained ones, which is what makes it automatic.

Settlements, springs and hinges

Support settlement is a prescribed displacement. Set the settled value in \(\mathbf{D}_s\) and move its effect to the right-hand side: \(\mathbf{K}_{ff}\mathbf{D}_f = \mathbf{P}_f - \mathbf{K}_{fs}\mathbf{D}_s\). This reproduces the settlement results of lesson 4 without any flexibility coefficients. An elastic support (a footing on soil, or a beam resting on a flexible girder) is a spring: its stiffness k is simply added to the diagonal term of the corresponding degree of freedom, and the support is then treated as free. A rotational spring at a support models a partially fixed connection in the same way, which is how semi-rigid connections are represented in practice.

Programs model internal hinges by releasing the moment at an element end. For a release at end j, the element matrix is condensed to the propped-cantilever form, in which the \(4EI/L\) term at i becomes \(3EI/L\), just as in the modified slope-deflection equation. The released rotation is then excluded from the structure's degrees of freedom at that end. Releasing both ends of every element turns a beam model into a truss model, which is a handy check of the release logic.

Solving it in code

The program below is a complete continuous-beam stiffness solver. It assembles element matrices and equivalent joint loads, partitions out the restrained degrees of freedom, includes prescribed settlements, solves by Gaussian elimination and recovers element end forces. Model 1 repeats the hand example. Model 2 repeats the settlement case of lesson 4.

# Matrix stiffness method for continuous beams (2 DOF per node: v up, rotation counterclockwise)
def gauss(A, b):
    n = len(b); M = [row[:] + [b[i]] for i, row in enumerate(A)]
    for c in range(n):
        p = max(range(c, n), key=lambda r: abs(M[r][c]))
        if abs(M[p][c]) < 1e-12: raise ValueError("singular stiffness matrix: the model is unstable")
        M[c], M[p] = M[p], M[c]
        for r in range(c + 1, n):
            f = M[r][c] / M[c][c]
            for k in range(c, n + 1): M[r][k] -= f * M[c][k]
    x = [0.0] * n
    for r in range(n - 1, -1, -1):
        x[r] = (M[r][n] - sum(M[r][k] * x[k] for k in range(r + 1, n))) / M[r][r]
    return x

def k_beam(EI, L):
    a, b, c, d = 12 * EI / L ** 3, 6 * EI / L ** 2, 4 * EI / L, 2 * EI / L
    return [[a, b, -a, b], [b, c, -b, d], [-a, -b, a, -b], [b, d, -b, c]]

def fixed_end(L, w, P, x):          # member end forces with both ends fixed: [V_i, M_i, V_j, M_j]
    b = L - x
    return [w * L / 2 + P * b * b * (3 * x + b) / L ** 3, w * L * L / 12 + P * x * b * b / L ** 2,
            w * L / 2 + P * x * x * (x + 3 * b) / L ** 3, -w * L * L / 12 - P * x * x * b / L ** 2]

def solve_model(xs, elems, restrained, settle=()):
    """xs: node positions; elems: (i, j, EI, w, P, a); restrained: DOF numbers; settle: (DOF, value)"""
    n = 2 * len(xs); K = [[0.0] * n for _ in range(n)]; F = [0.0] * n
    for i, j, EI, w, P, a in elems:
        L = xs[j] - xs[i]; k = k_beam(EI, L); fe = fixed_end(L, w, P, a)
        dofs = [2 * i, 2 * i + 1, 2 * j, 2 * j + 1]
        for r in range(4):
            F[dofs[r]] -= fe[r]                    # equivalent joint loads
            for c in range(4): K[dofs[r]][dofs[c]] += k[r][c]
    D = [0.0] * n
    for d, val in settle: D[d] = val
    free = [d for d in range(n) if d not in restrained]
    rhs = [F[r] - sum(K[r][c] * D[c] for c in restrained) for r in free]
    for d, val in zip(free, gauss([[K[r][c] for c in free] for r in free], rhs)): D[d] = val
    ends = []
    for i, j, EI, w, P, a in elems:
        L = xs[j] - xs[i]; k = k_beam(EI, L); fe = fixed_end(L, w, P, a)
        dofs = [2 * i, 2 * i + 1, 2 * j, 2 * j + 1]
        ends.append([sum(k[r][c] * D[dofs[c]] for c in range(4)) + fe[r] for r in range(4)])
    return D, ends

z = lambda v: round(v, 4) + 0.0
# Model 1: lesson 5 office floor beam, fixed at A, column at B, pin at C (EI = 1 gives EI x rotation)
D, ends = solve_model([0.0, 24.0, 44.0], [(0, 1, 1.0, 1.25, 0.0, 0.0), (1, 2, 1.0, 0.0, 13.5, 10.0)],
                  restrained=[0, 1, 2, 4])
print("model 1: EI x rotation at B = %.2f, at C = %.2f (counterclockwise +)" % (D[3], D[5]))
for e, (V1, M1, V2, M2) in enumerate(ends):
    print("  member %d: V_i %.2f kips, M_i %.2f kip*ft, V_j %.2f kips, M_j %.2f kip*ft" % (e + 1, V1, z(M1), V2, z(M2)))
# Model 2: lesson 4 corridor beam, EI = 72 000 kip*ft^2, support B settles 0.4 in., no load
D, ends = solve_model([0.0, 20.0, 40.0, 60.0], [(0, 1, 7.2e4, 0, 0, 0), (1, 2, 7.2e4, 0, 0, 0), (2, 3, 7.2e4, 0, 0, 0)],
                  restrained=[0, 2, 4, 6], settle=[(2, -0.4 / 12)])
print("model 2: moments at B = %.2f and at C = %.2f kip*ft (sagging +)" % (ends[0][3], ends[1][3]))

model 1: EI x rotation at B = 29.61, at C = 153.95 (counterclockwise +)
  member 1: V_i 15.31 kips, M_i 62.47 kip*ft, V_j 14.69 kips, M_j -55.07 kip*ft
  member 2: V_i 9.50 kips, M_i 55.07 kip*ft, V_j 4.00 kips, M_j 0.00 kip*ft
model 2: moments at B = 21.60 and at C = -14.40 kip*ft (sagging +)

Model 1 reproduces the rotations and end forces above, and the shears 14.69 + 9.50 = 24.19 kips give the column reaction at B. Model 2 gives \(+21.6\) kip·ft over B and \(-14.4\) kip·ft over C, the force method answer of lesson 4, from a model that needed no redundants at all. A singular-matrix check in the solver catches unstable models: try restraining only the vertical displacement at A in model 1 and removing the column at B.

Displacements between nodes are not computed by the program. To see the deflected shape of a span, either add intermediate nodes (which does not change the end moments, since the element is exact) or evaluate the cubic shape functions with the element end displacements and add the fixed-end deflection of the span load. Commercial programs do the second automatically when they draw deflected shapes.