Stiffness Matrix (K) — Direct Stiffness Method & FEA Fundamentals

The stiffness matrix (denoted K) is the mathematical heart of modern structural analysis. It encodes the force-displacement relationship of a structure or structural element into a compact matrix form. Combined with the direct stiffness method of assembly, the stiffness matrix transforms the intractable problem of analyzing a statically indeterminate structure into a systematic, computer-friendly linear algebra exercise.

At its core, the stiffness method states: {F} = [K]{u} — for a given vector of nodal forces {F}, the structure deforms by displacements {u}, and the stiffness matrix [K] defines the relationship between them.

PRELIMINARY — NOT FOR CONSTRUCTION. All content is for educational and reference use only. Must be independently verified by a licensed Professional Engineer (PE) or Structural Engineer (SE) before use in any project.

Physical Interpretation of Kij

An individual entry Kij in a stiffness matrix has a clear physical meaning:

Kij = the force (or moment) required at degree of freedom i to produce a unit displacement at degree of freedom j, with all other degrees of freedom held at zero.

This interpretation follows directly from matrix multiplication: when we multiply [K] by a displacement vector with uj = 1 and all other uk = 0, the resulting force vector is the j-th column of [K].

The Direct Stiffness Method — Step by Step

The direct stiffness method follows a six-step procedure:

Step 1: Discretize

Divide the structure into elements connected at nodes. Each node has degrees of freedom (DOFs) — translations and rotations that are unknown.

Step 2: Element Stiffness Matrices

For each element, compute its stiffness matrix in local coordinates. For a 2D beam element with nodes i and j, the 6×6 matrix is:

        [  EA/L       0          0      −EA/L       0          0    ]
        [    0    12EI/L³    6EI/L²       0    −12EI/L³    6EI/L² ]
k_e =   [    0     6EI/L²    4EI/L        0     −6EI/L²    2EI/L  ]
        [ −EA/L       0          0       EA/L       0          0    ]
        [    0    −12EI/L³  −6EI/L²       0     12EI/L³   −6EI/L² ]
        [    0     6EI/L²    2EI/L        0     −6EI/L²    4EI/L  ]

Where EA/L = axial stiffness and EI/L³, EI/L², EI/L = flexural stiffness terms.

Key properties of element stiffness matrices:

Step 3: Transformation to Global Coordinates

Elements are rarely aligned with the global coordinate system. A rotation matrix R transforms local displacements to global coordinates:

u_local = R × u_global
k_global = Rᵀ × k_local × R

For a 2D frame element, R is a 6×6 matrix containing direction cosines.

Step 4: Assembly

The global stiffness matrix is assembled by adding element contributions at shared nodes. For a structure with n nodes, each with d DOFs, the global K is (n×d) × (n×d). The assembly process is:

K_global = Σ (Lᵢᵀ × k_global,i × Lᵢ)

Where Lᵢ is the Boolean connectivity matrix that maps element DOFs to global DOFs.

Step 5: Apply Boundary Conditions

Without boundary conditions, K_global is singular — the structure can undergo rigid-body motion. Boundary conditions (supports, constraints) are applied by:

Penalty method: Add a large number (e.g., 10¹⁰ × max|Kii|) to the diagonal term for each constrained DOF and set the corresponding force to (large number × prescribed displacement).

Condensation method (preferred): Partition K into free and constrained DOFs:

[ Kff  Kfc ]  { uf }   { Ff }
[ Kcf  Kcc ]  { uc } = { Fc }

Then solve: Kff × uf = Ff − Kfc × uc. If uc = 0 (fixed supports), simply Kff × uf = Ff.

Step 6: Solve for Displacements and Back-Calculate Forces

Solve the linear system: {uf} = [Kff]⁻¹ × {Ff}. With displacements known, element forces are recovered:

{fe} = [k_e] × {ue}

Types of Element Stiffness Matrices

2D Truss Element (4×4)

        [  1   0  −1   0  ]
k_e = EA [  0   0   0   0  ]   (in local coordinates, aligned with element axis)
       L  [ −1   0   1   0  ]
        [  0   0   0   0  ]

Truss elements carry only axial force (no moment or shear). The off-diagonal zeros reflect zero coupling between axial and transverse displacements.

2D Beam Element (6×6)

Includes axial stiffness (EA/L), flexural stiffness (EI/L³, EI/L², EI/L), and the coupling between shear and moment at the ends.

2D Frame Element (6×6)

Same as the beam element above — the term "frame element" is used when members are connected at angles and carry combined axial + flexural loads.

3D Beam Element (12×12)

Each node has 6 DOFs: translations (ux, uy, uz) and rotations (θx, θy, θz). The matrix includes:

Spring Element (2×2)

        [  k   −k  ]
k_e =   [ −k    k  ]

The simplest element. Used for elastic supports, soil springs, and connection flexibility modeling.

Matrix Sparsity and Bandwidth

For large structures, the global stiffness matrix is sparse — most entries are zero. Each row has non-zero entries only for DOFs that share an element with the DOF corresponding to that row. For a typical frame structure, the matrix bandwidth (half-bandwidth) is approximately:

b ≈ (max node number difference in any element) × (DOFs per node)

Node numbering matters. Poor numbering produces a wide bandwidth; optimal numbering minimizes it. Modern FEA software uses algorithms like Reverse Cuthill-McKee (RCM) or Gibbs-Poole-Stockmeyer to renumber nodes automatically.

Stiffness vs. Flexibility

The stiffness method is the dual of the flexibility (force) method:

Stiffness Method Flexibility Method
Primary unknowns Displacements (u) Forces (redundants)
Matrix K (stiffness) F (flexibility)
Relationship {F} = [K]{u} {u} = [F]{F}
Best for Computer implementation (systematic assembly) Hand calculation (few redundants)
Boundary conditions Easy (condense constrained DOFs) Requires selecting a determinate base system

The stiffness method dominates modern practice because assembly is algorithmic, boundary conditions are straightforward, and the matrix is well-conditioned for most structures.

Stiffness Matrix Assembly Example

Consider a simple 4-DOF system: two springs in series:

Spring 1: stiffness k₁, connects node 1 to node 2
Spring 2: stiffness k₂, connects node 2 to node 3

[ k₁    −k₁    0   ] { u₁ }   { F₁ }
[ −k₁  k₁+k₂ −k₂  ] { u₂ } = { F₂ }
[ 0     −k₂   k₂  ] { u₃ }   { F₃ }

After applying boundary condition u₁ = 0 (fixed support), condensing row/column 1:

[ k₁+k₂  −k₂  ] { u₂ }   { F₂ }
[ −k₂     k₂  ] { u₃ } = { F₃ }

Solving for u₂, u₃ → spring forces: f₁ = k₁(u₂ − u₁) = k₁u₂, f₂ = k₂(u₃ − u₂).

Computational Considerations

Storage

A full N×N matrix for N = 100,000 DOFs would require 80 GB (double precision). In practice, sparse storage formats (CSR, skyline/profile) reduce this to a few hundred MB — only non-zero entries near the diagonal are stored.

Solution Methods

Ill-Conditioning

If some elements are much stiffer than others (e.g., a beam with very short span adjacent to a long span), the stiffness matrix can become ill-conditioned. The condition number κ = |λmax|/|λmin| measures this — large κ leads to loss of precision. Remedies include scaling, higher-precision arithmetic, or reformulating the element.

Frequently Asked Questions

Why is the stiffness matrix symmetric?

Symmetry follows from Maxwell-Betti's reciprocal theorem: the work done by forces of type i due to displacements of type j equals the work done by forces of type j due to displacements of type i. Mathematically, Kij = Kji. This holds for linear elastic structures with conservative loading.

Why is the stiffness matrix singular without boundary conditions?

Without boundary conditions, the structure has rigid-body modes — the whole structure can translate or rotate without straining any element. These rigid-body motions correspond to zero eigenvalues of K. Applying sufficient boundary conditions removes these modes, making K non-singular and solvable.

How does the direct stiffness method handle member releases (pins)?

Member-end releases (pinned connections, moment releases) are handled by condensing the released DOF out of the element stiffness matrix before assembly. For a beam with a pinned end at node j (zero moment), the rotational DOF θj is condensed out, reducing the 6×6 matrix to 5×5. This is equivalent to static condensation at the element level.

What is the geometric stiffness matrix (Kg)?

The geometric stiffness matrix (also called stress stiffness matrix) accounts for the effect of axial force on flexural stiffness — the P-δ effect. Compression reduces the effective flexural stiffness (K_total = K_elastic + K_geometric, where K_geometric is negative for compression). This is the basis for linearized buckling analysis: solve det(K_elastic + λK_geometric) = 0 for the buckling load factor λ.

Related Terms and Pages

Educational reference only. Structural analysis results from any stiffness-based method must be independently verified by a licensed Professional Engineer using an alternative method or software for critical structures.

Disclaimer: This content is for educational purposes only. Results must be verified by a licensed professional engineer. Steel Calculator provides preliminary design tools — NOT a substitute for professional engineering judgment.