Stiffness matrix, step by step, for plane trusses and frames

The direct stiffness method turns a whole structure into one matrix equation, K·u = P + R, and solves it. K is the stiffness matrix of the structure, u the displacements of all its nodes, P the loads put on them and R the forces the supports add. The tool works as a person does on paper and prints every step, with its numbers, its symbols and the cells that matter marked. What follows is the same method, in the order the tool lists it.

Degrees of freedom (DOFs) come first (step 1). A truss node has two, a displacement along x and one along y; a frame node has a third, its rotation. They are numbered node by node: node n has DOFs 2n − 1 and 2n in a truss, and 3n − 2, 3n − 1 and 3n in a frame. That order is the order of the rows and columns of every global matrix and of every global vector, so the number of a row says which movement it is about. A DOF that a support holds is restrained, or held; every other one is free.

Each member is first taken in axes of its own (steps 2 and 3), with x̄ running along it from its first node i to its second node j. A bar, which carries axial force alone, has one movement at each end, and its matrix is k = (EA/L)·[1 −1; −1 1], where E is the modulus, A the area and L the length. A column of any stiffness matrix is the set of forces that hold a member in the shape of one unit movement, every other movement being zero. Here the first column says that shortening the bar by one unit takes a push of EA/L on each end.

A frame member is a beam-column with six movements, a displacement along x̄, one along ȳ and a turn at each end, so its matrix is 6 × 6. It is made of five numbers: EA/L for the axial terms and, for bending, 12EI/L³, 6EI/L², 4EI/L and 2EI/L, each cell taking the sign its position gives. Move the first end sideways by one unit and the member answers with a shear of 12EI/L³ and a moment of 6EI/L² at that end. Turn it by one unit and the moment there is 4EI/L, and 2EI/L at the far end.

Members point every way, so each is turned into the global axes before it is added (steps 4 and 5). With c = cosθ and s = sinθ for the member’s direction, pointing away from node i toward node j, the end movements in global axes give the local ones by u_local = T·u_global. The matrix T holds c and s; for a frame member it is [c s 0; −s c 0; 0 0 1] at each end, since a rotation is the same in every set of axes. The member’s matrix in global axes is K_e = T^T·k·T, with T^T the transpose of T. For a bar it is (EA/L)·[c² cs −c² −cs; cs s² −cs −s²; −c² −cs c² cs; −cs −s² cs s²], which the tool prints beside the product as a second route to the same figures.

Assembly adds the members up (step 6). Each member has a code vector: the numbers of the global DOFs that its own DOFs correspond to, in order, as in [1 2 5 6] for a bar from node 1 to node 3. The entry in row a and column b of K_e is added to K at row code(a) and column code(b), K ← K + K_e. Two members that meet at a node both add into the cells of that node’s DOFs, and their sum is the stiffness the joint offers to that movement. The tool shows each member’s contribution, the running K after it and, last, the whole K with the cells that more than one member reached marked.

Loads come next (step 7). A load at a node goes into the load vector P at the DOF it acts along, and a couple at a node goes to its rotation. A load on a member (frame mode) is a uniform load over its length, a point load or an axial load, each along the member’s own axes, and each is exchanged for its fixed-end forces f_FE: what two clamps would exert to keep both ends from moving. A uniform load w has wL/2 and wL²/12 (in size), an axial load q has qL/2 along the member at each end, and a point load P, with a and b the distances to its first and second ends, has Pb²(3a+b)/L³ and Pab²/L² at the first end and Pa²(a+3b)/L³ and Pa²b/L² at the second. The equivalent nodal loads, −T^T·f_FE, are their opposite in global axes and join P; the member’s end forces take f_FE back when they are found.

The DOFs are then split into free (F) and restrained (R), and K, P and u are cut into blocks to match (step 8). The free rows give K_FF·u_F = P_F − K_FR·u_R (step 9). Here u_R holds what the supports allow: zero, unless a support is given a settlement, which then acts as a load, −K_FR·u_R. The restrained rows are kept for later, when they give the reactions.

This reduced equation is solved by Gaussian elimination, without exchanging rows (step 10). At column k the pivot a_kk is the diagonal entry; for each row i below it, l_ik = a_ik/a_kk is the multiple of row k to take away, a_ij ← a_ij − l_ik·a_kj, which leaves zeros under the pivot. When the matrix is triangular the unknowns come from the bottom up, u_k = (b_k − Σ a_kj·u_j)/a_kk. No exchange is needed because K_FF is symmetric and positive definite when the structure is stable, so every pivot is positive.

A pivot under a billionth of the diagonal entry it started as is taken as zero. It means that nothing in the structure resists some movement: the members can swing without stretching or shortening, and that is a mechanism. The tool stops there, shows the steps worked out until then and names the DOF. A smallest pivot under a millionth of its start, but above a billionth, is solved with a warning: the structure is nearly a mechanism, and its displacements change a great deal if any one stiffness changes even a little.

The free displacements and the prescribed ones make up the whole vector u (step 11). The reactions come from the restrained rows, R = K_RF·u_F + K_RR·u_R − P_R, one for each restrained DOF (step 12). A reaction is the force or moment a support applies to the structure.

Each member then reads its end movements u_e from u through its code vector, turns them to local axes, u_local = T·u_e, and multiplies by its matrix, f = k·u_local, adding f_FE if it carries a load (step 13). For a bar, N = (EA/L)(u_j − u_i) in local axes; the tool also gives the stress N/A and the elongation. For a frame member, N, V and M at both ends follow from f with the signs given further down.

Step 14 tests the answer. Loads and reactions are added up along x, along y and as a moment taken about the centre of the box around the nodes; every joint is balanced against the end forces of the members that meet there; the strain energy U = ½·u^T·K·u is compared with the work ½·u^T·P + ½·u_R^T·R; and K·u − P − R is read at every DOF, with the elimination’s own residual. Each is measured against the size of the terms that make it up, with a limit of a hundred-millionth. A structure for which one of them is missed is refused as ill-conditioned, and no results are given for it.

A first worked case is the two-bar truss the page opens with. Nodes 1, 2 and 3 are at (0, 0), (16, 0) and (8, 6) ft. Bar 1 joins nodes 1 and 3, bar 2 joins nodes 2 and 3, and nodes 1 and 2 are pinned. Both bars have E = 29,000 ksi and A = 10 in², and node 3 carries 20 kip to the right and 30 kip down.

Each bar is 10 ft long, so EA = 290,000 kip and EA/L = 29,000 kip/ft. Bar 1 has c = 0.8 and s = 0.6; bar 2 runs from node 2 up to node 3, so c = −0.8 and s = 0.6. The DOFs are 1 and 2 at node 1, 3 and 4 at node 2, and 5 and 6 at node 3, so the code vectors are [1 2 5 6] and [3 4 5 6]. At DOFs 5 and 6, bar 1 adds 29,000 × [0.64 0.48; 0.48 0.36] and bar 2 adds 29,000 × [0.64 −0.48; −0.48 0.36]. The cross terms cancel, and K_FF = [37,120 0; 0 20,880] kip/ft.

With P_F = [20; −30] kip and the supports not moving, nothing lies under the diagonal to eliminate: u_5 = 20/37,120 and u_6 = −30/20,880, which are 0.00647 in to the right and 0.01724 in down. The reactions are R_1 = 10, R_2 = 7.5, R_3 = −30 and R_4 = 22.5 kip. N = (EA/L)(u_j − u_i) is −12.5 kip in bar 1 and −37.5 kip in bar 2, both in compression, and ΣFx = 20 + 10 − 30 = 0 and ΣFy = −30 + 7.5 + 22.5 = 0. By joint equilibrium at node 3, −0.8·N_1 + 0.8·N_2 + 20 = 0 and −0.6·N_1 − 0.6·N_2 − 30 = 0 give the same two forces.

A frame needs bending, so take one member of length L, fixed at node 1 and loaded by P downward at node 2. The load has no axial part: the axial row holds only EA/L, nothing couples it to the others and its movement comes out zero. The other two free rows, for v_2 and θ_2, are K_FF = [12EI/L³ −6EI/L²; −6EI/L² 4EI/L] with P_F = [−P; 0]. The multiplier is l = (−6EI/L²)/(12EI/L³) = −L/2; taking l times the first row from the second leaves EI/L where 4EI/L stood and −PL/2 on the right. Back-substitution gives θ_2 = −PL²/(2EI) and v_2 = −PL³/(3EI), the textbook rotation and deflection at the tip, and the reactions at node 1 are P upward and PL counter-clockwise.

Reading the matrices. Each matrix and vector has its rows and columns named by DOF: 3y is the y movement of node 3 and 2θ the rotation of node 2. Row d of K is the equation of equilibrium of DOF d. Column d is the set of forces, at every DOF, that hold the structure in the shape of a unit movement of DOF d, every other DOF being held still, so the entry in row a and column b is the force at a when b moves one unit.

The unit of an entry follows from its two DOFs: a force per length between two movements (kip/ft, or kN/m), a force per radian between a movement and a rotation (kip/rad, kN/rad), and a moment per radian between two rotations (kip·ft/rad, kN·m/rad).

K is symmetric: the force at a for a unit movement of b is the force at b for a unit movement of a, which is Maxwell and Betti’s reciprocity, and the tool tests it. No diagonal entry is negative. A zero off the diagonal means that moving one DOF asks no force at the other: no member joins their nodes, or the member that does has no stiffness across that direction (a horizontal bar and a movement across it), or terms cancelled, as the cross terms did in the example. Before supports are applied K has no inverse: a movement of the whole structure as a rigid body takes no force at all, so K·u is zero for it, whatever the loads. The supports remove those movements, and K_FF can be solved.

Within a member matrix the axial terms and the bending terms share no cell in local axes, because stretching a member does not bend it; the turn to global axes mixes them. The tool marks the axial and the bending terms of k, the cells where a member lands in the assembly and those that more than one member reached, and in the elimination the pivot, the cells it reduces to zero and those it changes.

The entries are physical stiffnesses. For a column of height h fixed at its foot, 12EI/h³ is the force that sways its top one unit with the top kept from turning. A beam that cannot turn, on two such columns that cannot shorten, makes a stiffness of 24EI/h³, so a horizontal force H sways it by Hh³/(24EI).