Matrix Stiffness and Direct Stiffness Methods

Introduction and Context

The matrix stiffness method is one of the most important general procedures in modern structural analysis. It provides a systematic relationship between nodal forces and nodal displacements and forms the theoretical foundation of the finite element method used in structural engineering software. Unlike classical force methods, which take member forces as primary unknowns, the stiffness method treats joint displacements as the fundamental unknowns.

For a linear elastic structure, the central relationship is

where K is the assembled global stiffness matrix, d is the vector of unknown nodal displacements and rotations, and F is the vector of applied nodal forces and moments. The equation applies under the assumptions of small displacement, linear elastic material behaviour, compatible deformation, and properly specified boundary conditions.

The direct stiffness method is an implementation-oriented form of the matrix stiffness method. It constructs the structure’s global stiffness matrix directly from individual member stiffness matrices, transforms member quantities into global coordinates, assembles the contributions at common nodes, applies restraints, and solves the resulting simultaneous equations.

The method is particularly valuable for continuous beams, plane and space trusses, rigid frames, grids, plates, shells, and three-dimensional finite element models.

Learning Objectives

By the end of this chapter, the reader should be able to:

1. Explain the displacement-based basis of the matrix stiffness method.
2. Derive and interpret local stiffness matrices for common structural elements.
3. Transform element stiffness matrices between local and global coordinate systems.
4. Assemble a global stiffness matrix using the direct stiffness procedure.
5. Apply support conditions and solve for nodal displacements.
6. Recover member end forces, reactions, and internal actions.
7. Identify singularity, mechanisms, ill-conditioning, and modelling errors.
8. Evaluate the practical use of matrix methods in structural software and finite element analysis.

1. Fundamental Concepts

1.1 Degrees of freedom

A degree of freedom (DOF) is an independent displacement or rotation required to describe the configuration of a structural model. A plane truss joint generally has two translational DOFs, u and v. A plane frame joint generally has three DOFs: horizontal translation u, vertical translation v, and rotation θ.

The nodal displacement vector is written as

The corresponding nodal action vector contains forces and moments:

Translations must be expressed in a length unit and rotations in radians. Forces and moments must be dimensionally compatible with the stiffness coefficients.

1.2 Flexibility and stiffness

Flexibility relates deformation to applied force, whereas stiffness relates force to deformation. For a single axial member,

where N is axial force, E is Young’s modulus, A is cross-sectional area, L is member length, and Δ is axial extension. The axial stiffness EA/L has units of force per length.

For a general structure, the equivalent relationship is F = Kd. The matrix K therefore contains the force required at each degree of freedom when a unit displacement is imposed at one DOF while all other free DOFs are held fixed.

1.3 Energy interpretation

For a linearly elastic structure, the strain energy is

The matrix K is symmetric when the material is elastic and the structural formulation is conservative. Positive definiteness is expected for a stable structure after all rigid-body motions have been restrained. A free structure may have zero-energy rigid-body modes, so its unrestrained stiffness matrix is generally positive semidefinite rather than positive definite.

2. Element Stiffness Matrices

The direct stiffness method begins with an element-level force–displacement relationship. The element matrix is derived using equilibrium, compatibility, constitutive laws, or the principle of virtual work.

2.1 Two-node axial bar in local coordinates

Consider a prismatic axial member with local nodal displacements u₁ and u₂. The displacement field is assumed linear:

The constant axial strain is

Using σ = Eε and N = Aσ, the element force relationship is

The signs depend on the selected positive directions and on whether end actions are defined as forces exerted on or by the element. Consistent sign conventions must be used during assembly and force recovery.

2.2 Plane truss element in global coordinates

Let a two-dimensional truss member connect nodes i and j. Define

where θ is the angle measured from the global x-axis. The global displacement vector is

The extension of the member is

Thus the plane truss stiffness matrix is

The matrix has units of force per length. It is symmetric and has rank one for an ideal axial member because the element resists only extension along its axis. It does not resist transverse deformation or nodal rotation.

2.3 Euler–Bernoulli beam element

For a two-node beam in local coordinates, each node has a transverse displacement and a rotation. The local DOF vector is

For a prismatic Euler–Bernoulli beam with constant EI, the local bending stiffness matrix is

Here E is Young’s modulus, I is the second moment of area, and L is element length. The formulation assumes small deflection, plane sections remaining plane, negligible shear deformation, constant properties, and linear elastic behaviour. For deep beams or short members, a Timoshenko beam formulation may be more appropriate because it includes transverse shear deformation.

2.4 Plane frame element

A two-dimensional frame element combines axial and bending actions. Its local DOF vector is

The local stiffness matrix is

This matrix assumes that axial and bending deformations are uncoupled in the local member coordinate system. It is widely used for plane rigid-jointed frames.

3. Coordinate Transformation

Element matrices are often most easily derived in local coordinates, whereas nodal displacements and loads are stored in global coordinates. A transformation is therefore required.

For a plane frame member, define the transformation matrix T such that

where dₑˡ and dₑᵍ are local and global element displacement vectors. For a member at angle θ,

The global element stiffness matrix is

The transpose appears because virtual work and energy must be invariant under coordinate transformation. The same transformation must be used consistently for displacements, forces, and stiffness matrices.

4. Assembly of the Global Stiffness Matrix

4.1 Direct stiffness procedure

The standard procedure is:

1. Define nodes, elements, coordinates, material properties, sections, and DOF numbering.
2. Construct each element stiffness matrix in local coordinates.
3. Transform each matrix into global coordinates if necessary.
4. Add element contributions into the appropriate locations of the global matrix.
5. Assemble the global load vector, including nodal loads and equivalent nodal loads from distributed actions.
6. Apply support restraints and prescribed displacements.
7. Solve for the unknown nodal displacements.
8. Recover reactions and element end forces.
9. Check equilibrium, compatibility, and, where possible, energy consistency.

4.2 Assembly operator

Let Aₑ be a Boolean or signed connectivity matrix that maps element DOFs to global DOFs. The global stiffness matrix can be expressed compactly as

Similarly, element load vectors may be assembled through

where fₑ is an element equivalent nodal load vector and P is the externally applied nodal load vector. This expression is especially useful in finite element programming because it separates element calculations from global numbering.

4.3 Equivalent nodal loads

A distributed load should not normally be inserted directly into a nodal load vector without conversion. It must be replaced by an equivalent nodal force vector that produces the same virtual work for the selected interpolation functions.

For a two-node beam subjected to a uniform transverse load w over length L, one common local consistent-load vector is

The signs depend on the chosen positive transverse direction, rotation convention, and whether the vector represents applied loads or fixed-end actions. This is an important area where software conventions differ; the load vector must be consistent with the element stiffness and end-force definitions.

5. Boundary Conditions and Solution

Partition the assembled system into free and restrained DOFs:

If restrained displacements are zero, the reduced equations are

If nonzero prescribed displacements exist, then

The reactions at restrained DOFs are obtained from

The exact sign of the reported reaction depends on whether F is defined as the applied external load vector or the total nodal action vector. A complete analysis must clearly state this convention.

Common implementation techniques include elimination, the penalty method, Lagrange multipliers, and multi-point constraint transformations. Elimination is generally simple and numerically reliable for standard support conditions. Penalty methods are convenient but may introduce poor conditioning if the penalty parameter is excessively large or small. Lagrange multipliers enforce constraints exactly but enlarge the system and may produce an indefinite matrix.

6. Worked Example: Two-Element Axial Bar

Problem

A straight axial bar consists of two collinear elements. Node 1 is fixed, node 2 is an internal node, and node 3 carries a tensile load of 30 kN. Each element has E = 200 GPa, A = 600 mm², and length L = 1.5 m. Determine the displacement of node 2, displacement of node 3, and support reaction.

Consistent units

Use N and mm. Then

E = 200,000 N/mm², A = 600 mm², and L = 1,500 mm.

The element stiffness is

Assembly

For two identical axial elements, the global stiffness matrix before applying the support condition is

Node 1 is restrained, so u₁ = 0. The reduced system for nodes 2 and 3 is

Dividing by 80,000 gives

The second equation gives u₃ = u₂ + 0.375. Substitution into the first gives

2u₂ − (u₂ + 0.375) = 0, hence

and

The reaction at node 1 is 30 kN in compression relative to the applied tensile load, so its magnitude is 30 kN. The extension of each element is 0.375 mm, giving an element force of

The result satisfies equilibrium and compatibility. It also illustrates that the node at the interface carries no externally applied load but still has a nonzero displacement because it is part of the continuous structure.

7. Numerical and Computational Considerations

7.1 Symmetry and sparsity

For conventional linear elastic structures, K = Kᵀ. Large structural models also produce sparse matrices because each element connects only a small number of DOFs. Sparse storage and factorization are therefore essential for computational efficiency.

7.2 Singularity and mechanisms

A singular or nearly singular stiffness matrix may indicate:

• an unconstrained rigid-body motion;
• an unstable structural mechanism;
• an incorrectly connected element;
• a missing rotational restraint;
• duplicate or disconnected nodes;
• an incorrect transformation matrix; or
• an inappropriate element formulation.

A singular matrix is not merely a numerical inconvenience. It may represent a physically unstable structure. Artificially restraining a DOF to force a solution can conceal a serious modelling error.

7.3 Conditioning

The condition number of a matrix measures sensitivity of the solution to numerical perturbations. Large differences in material stiffness, element size, section properties, or units can create ill-conditioning. Good practice includes consistent units, appropriate scaling, avoidance of excessively distorted elements, and use of numerically stable solvers.

7.4 Linear versus nonlinear analysis

The basic direct stiffness method is linear when K is constant. In nonlinear analysis, the tangent stiffness matrix changes with displacement, material state, contact condition, or geometry. The equilibrium equations are solved incrementally, commonly using Newton–Raphson iteration:

where Kₜ is the tangent stiffness matrix, R_ext is the external load vector, R_int is the internal resisting-force vector, and Δuᵢ is the iterative displacement correction. Convergence requires specified tolerances on residual forces, displacement increments, or energy error.

8. Applied Case: Analysis of a Small Portal Frame

Background

Consider an illustrative single-storey steel portal frame with rigid beam-to-column joints. The columns are fixed at their bases, and a horizontal wind load acts at the beam level. The frame is modelled using two-dimensional frame elements with six DOFs per element end.

Problem

The design team must estimate lateral drift, column bending moments, and base reactions before selecting preliminary member sizes. A hand calculation based only on independent column cantilevers is considered inadequate because beam–column interaction influences both stiffness and force distribution.

Analysis approach

The frame is idealized using nodes at the two bases and two beam–column joints. Each member is assigned E, A, and I. Local frame matrices are generated, transformed according to member orientation, and assembled into the global matrix. Base translations and rotations are restrained. The horizontal wind load is applied at the upper nodes or represented through equivalent nodal loads if it is distributed.

The reduced system is solved for joint translations and rotations. Member end forces are then calculated using

where qₑˡ is the local member-end action vector and fₑˡ is the consistent local element-load vector. The sign convention must be maintained when interpreting shear, axial force, and bending moment diagrams.

Outcome and engineering interpretation

The solution provides lateral drift, which can be checked against the applicable serviceability limit; base shears and overturning reactions, which support foundation design; and member end moments, which support strength design. A sensitivity study can then be performed by changing I, joint releases, or support conditions. This often reveals that drift is more sensitive to flexural stiffness than to axial area, whereas axial force distribution may be strongly affected by geometry and load path.

Learning points

1. The global model must reflect the actual connectivity and releases.
2. A beam element and a frame element are not interchangeable when axial deformation matters.
3. Member forces must be recovered in local coordinates.
4. Numerical output requires engineering verification through equilibrium and plausible deformation patterns.
5. A software result is not a substitute for selecting an appropriate structural idealization.

9. Professional Practice, Verification, and Limitations

The direct stiffness method is powerful but depends on modelling assumptions. Important limitations include:

• Linear elastic analysis does not predict yielding, cracking, crushing, or large-displacement effects without suitable extensions.
• Euler–Bernoulli elements may overestimate stiffness for deep beams.
• Rigid joints, hinges, offsets, and releases must be explicitly represented.
• Concentrated loads applied at inappropriate locations may create unrealistic local effects.
• Mesh refinement can change results, especially near discontinuities, supports, openings, and singular stress regions.
• Reaction and member-force signs can be misinterpreted if software conventions are not documented.

A professional verification process should include independent checks of total equilibrium, support reactions, symmetry where applicable, limiting cases, and dimensional consistency. Results should be compared with simple hand calculations for selected actions. For large finite element models, mesh convergence studies and sensitivity analyses are essential.

Structural engineers must also comply with the governing building code, material standards, design combinations, serviceability criteria, and quality-assurance requirements. The matrix method itself does not determine whether a structure is safe; it supplies an analysis model whose assumptions and results must be judged by competent professionals.

10. Contemporary Developments

Modern structural analysis systems extend the direct stiffness method through nonlinear material models, geometric stiffness, contact elements, dynamic analysis, probabilistic simulation, and parallel sparse solvers. Automated model generation and AI-assisted error detection can help identify disconnected elements, unusual stiffness values, or inconsistent boundary conditions. However, such tools should be treated as aids to engineering judgement. Explainability, traceability, version control, and independent validation remain necessary, particularly for safety-critical structures.

Chapter Summary

The matrix stiffness method expresses structural equilibrium through the relationship Kd = F, with nodal displacements as the primary unknowns. The direct stiffness method builds the global matrix from element stiffness matrices, coordinate transformations, connectivity information, and boundary conditions. Axial, truss, beam, and frame elements have distinct stiffness formulations and assumptions. Assembly is performed by placing element contributions into the appropriate global DOF locations. Support conditions reduce or constrain the system, after which nodal displacements, reactions, and member forces are recovered. Singularity may indicate a mechanism or modelling error, while poor conditioning can result from incompatible units or extreme stiffness contrasts. Reliable engineering practice requires equilibrium checks, appropriate element selection, code compliance, and critical interpretation of software results.

Key Terms

Degree of freedom: An independent displacement or rotation in a structural model.

Global stiffness matrix: The assembled matrix relating global nodal actions to global nodal displacements.

Element stiffness matrix: The force–displacement matrix of an individual structural element.

Direct stiffness method: A procedure that assembles element matrices directly into a global structural matrix.

Local coordinates: Coordinates attached to an element, usually with one axis along its longitudinal direction.

Global coordinates: A common reference coordinate system for the complete structure.

Transformation matrix: A matrix that converts element quantities between local and global coordinates.

Assembly: The process of adding element contributions to the corresponding global DOF locations.

Equivalent nodal load: Nodal loading that reproduces the virtual work of a distributed element load.

Boundary condition: A prescribed displacement, rotation, force, or support constraint.

Reaction: An action developed at a restrained degree of freedom.

Mechanism: A kinematic mode that permits deformation without corresponding strain energy.

Singular matrix: A matrix that cannot be uniquely inverted, often indicating instability or missing constraints.

Conditioning: The numerical sensitivity of a system of equations to perturbations.

Consistent mass or load formulation: A formulation derived using the same interpolation functions as the element displacement field.

Euler–Bernoulli beam: A beam model neglecting transverse shear deformation.

Timoshenko beam: A beam model including transverse shear deformation and rotational effects.

Tangent stiffness: The incremental stiffness used in nonlinear equilibrium iterations.

Compatibility: The requirement that connected elements have geometrically consistent deformations.

Virtual work: Work-equivalence principle used to derive element equations and equivalent nodal actions.

Review Questions

Short Answer Questions

1. What is the physical meaning of a coefficient in a structural stiffness matrix?
2. Why is the global element matrix obtained from kᵍ = TᵀkˡT?
3. State two assumptions of the Euler–Bernoulli beam element.
4. What are two possible causes of a singular global stiffness matrix?
5. Why must distributed loads be converted into equivalent nodal loads in a displacement-based formulation?

Long Answer Questions

1. Derive the stiffness matrix of a two-node axial bar using the linear displacement field and constitutive relation.
2. Explain the complete direct stiffness algorithm for a plane frame, from structural idealization to member-force recovery.
3. Derive the coordinate transformation concept for a plane truss element and explain the roles of c and s.
4. Compare elimination, penalty, and Lagrange-multiplier techniques for applying support constraints.
5. Discuss the sources and engineering significance of singularity and ill-conditioning in structural analysis.

Analytical / Application Questions

1. A two-element axial bar has EA/L = 50,000 N/mm for each element. Node 1 is fixed and a 20 kN tensile load acts at node 3. Assemble the reduced stiffness matrix and calculate the nodal displacements.
2. Develop a plane-frame model for a single-storey portal frame subjected to lateral load. Specify nodes, DOFs, restraints, element properties, load vector, verification checks, and expected deformation patterns.
3. A beam model gives a lateral drift twice the value predicted by a hand calculation. Identify at least five modelling or interpretation checks that should be completed before accepting either result.

Multiple Choice Questions with Answers

Question 1

In the matrix stiffness method, the primary unknowns are generally:

A. Member bending moments

B. Nodal displacements and rotations

C. Support reactions only

D. Material stresses only

Correct Answer: B. Nodal displacements and rotations

Question 2

For a two-node axial bar, the elemental axial stiffness is:

A. EA/L

B. EI/L³

C. EA/L³

D. EI/L

Correct Answer: A. EA/L

Question 3

The global stiffness matrix of a transformed element is obtained as:

A. TkTᵀ

B. T⁻¹kT

C. TᵀkT

D. kTᵀTᵀ

Correct Answer: C. TᵀkT

Question 4

Which property normally characterizes the stiffness matrix of a conservative linear elastic structure?

A. It is antisymmetric

B. It is symmetric

C. It contains only negative terms

D. It is always diagonal

Correct Answer: B. It is symmetric

Question 5

A singular global stiffness matrix may indicate:

A. A perfectly constrained stable structure

B. A rigid-body mode or structural mechanism

C. An automatically converged nonlinear solution

D. A dimensionally consistent unit system

Correct Answer: B. A rigid-body mode or structural mechanism

Question 6

For the two-element axial-bar example in the chapter, the displacement of node 3 is:

A. 0.125 mm

B. 0.375 mm

C. 0.750 mm

D. 1.500 mm

Correct Answer: C. 0.750 mm

Question 7

The Euler–Bernoulli beam formulation generally neglects:

A. Axial equilibrium

B. Transverse shear deformation

C. Bending stiffness

D. Nodal rotations

Correct Answer: B. Transverse shear deformation

Question 8

The purpose of assembling element matrices is to:

A. Remove all support conditions

B. Place element contributions into the corresponding global DOF locations

C. Convert forces into stresses without displacement calculations

D. Eliminate the need for coordinate systems

Correct Answer: B. Place element contributions into the corresponding global DOF locations

Question 9

In a nonlinear incremental analysis, the matrix used in a Newton–Raphson iteration is commonly called the:

A. Flexibility matrix

B. Tangent stiffness matrix

C. Mass-normalization matrix

D. Connectivity matrix

Correct Answer: B. Tangent stiffness matrix

Question 10

For the axial-bar example, each element force after applying the 30 kN end load is:

A. 10 kN

B. 15 kN

C. 20 kN

D. 30 kN

Correct Answer: D. 30 kN