Lesson 4 of 13 25 min

Element Stiffness Matrix

The element stiffness matrix is the heart of FEA. It relates nodal forces to nodal displacements within a single element — essentially a multi-dimensional version of the spring equation $F = kx$.

In this lesson, we'll derive the stiffness matrix from first principles and build physical intuition for what each entry means.

The Spring Analogy

Before diving into matrices, let's recall the simple spring:

$$F = k \cdot \delta$$

Where:

  • $F$ = force applied to spring
  • $k$ = spring stiffness (N/m)
  • $\delta$ = displacement (extension)

A spring with higher $k$ requires more force to achieve the same displacement — it's "stiffer."

See how individual spring stiffnesses combine into a system stiffness matrix. Click nodes to apply displacements.

Two Springs in Series

For two springs connected end-to-end with a shared node:

NodeForces
Node 1$F_1 = k_1(u_1 - u_2)$
Node 2$F_2 = k_1(u_2 - u_1) + k_2(u_2 - u_3)$
Node 3$F_3 = k_2(u_3 - u_2)$

In matrix form:

$$\begin{bmatrix} F_1 \\ F_2 \\ F_3 \end{bmatrix} = \begin{bmatrix} k_1 & -k_1 & 0 \\ -k_1 & k_1+k_2 & -k_2 \\ 0 & -k_2 & k_2 \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \\ u_3 \end{bmatrix}$$

This is exactly what FEA does — but for continuous structures instead of discrete springs!

From Weak Form to Stiffness Matrix

Recall the weak form for a 1D bar:

$$\int_0^L EA\frac{du}{dx}\frac{dv}{dx} \, dx = \int_0^L v \cdot f \, dx + v(L) \cdot F_L$$

Now we apply the Galerkin method:

  • Approximate $u(x) = \sum_j N_j(x) u_j$
  • Choose $v(x) = N_i(x)$

Substituting:

$$\int_0^L EA \frac{d}{dx}\left(\sum_j N_j u_j\right) \frac{dN_i}{dx} \, dx = \int_0^L N_i \cdot f \, dx + N_i(L) \cdot F_L$$

Since $u_j$ are constants (nodal values), we can pull them out:

$$\sum_j \left[\int_0^L EA \frac{dN_j}{dx} \frac{dN_i}{dx} \, dx\right] u_j = \int_0^L N_i \cdot f \, dx + N_i(L) \cdot F_L$$

This is the matrix equation:

$$[K]\{u\} = \{F\}$$

Where:

$$K_{ij} = \int_0^L EA \frac{dN_i}{dx} \frac{dN_j}{dx} \, dx$$

Deriving [K] for a 2-Node Bar Element

Consider a bar element with:

  • Length $L$
  • Constant $EA$ (uniform cross-section and material)
  • Nodes at $x = 0$ and $x = L$

Shape Functions

$$N_1(x) = 1 - \frac{x}{L}, \quad N_2(x) = \frac{x}{L}$$

Shape Function Derivatives

$$\frac{dN_1}{dx} = -\frac{1}{L}, \quad \frac{dN_2}{dx} = \frac{1}{L}$$

Computing Stiffness Entries

$K_{11}$:

$$K_{11} = \int_0^L EA \left(-\frac{1}{L}\right)\left(-\frac{1}{L}\right) dx = \frac{EA}{L^2} \int_0^L dx = \frac{EA}{L^2} \cdot L = \frac{EA}{L}$$

$K_{12} = K_{21}$:

$$K_{12} = \int_0^L EA \left(-\frac{1}{L}\right)\left(\frac{1}{L}\right) dx = -\frac{EA}{L^2} \cdot L = -\frac{EA}{L}$$

$K_{22}$:

$$K_{22} = \int_0^L EA \left(\frac{1}{L}\right)\left(\frac{1}{L}\right) dx = \frac{EA}{L^2} \cdot L = \frac{EA}{L}$$

The Element Stiffness Matrix

$$[K_e] = \frac{EA}{L} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}$$

This is the element stiffness matrix for a 2-node bar element!

Adjust E, A, and L using the sliders. Apply unit displacement to see the resulting forces.

Physical Interpretation

Each entry $K_{ij}$ has a clear physical meaning:

$K_{ij}$ = Force at node $i$ due to unit displacement at node $j$ (with all other displacements = 0)

For the 2-node bar:

EntryMeaning
$K_{11} = +EA/L$Push node 1 right by 1 → Reaction force at node 1 pointing left
$K_{12} = -EA/L$Push node 2 right by 1 → Force at node 1 pointing right (bar pulls)
$K_{21} = -EA/L$Push node 1 right by 1 → Force at node 2 pointing right (bar pushes)
$K_{22} = +EA/L$Push node 2 right by 1 → Reaction force at node 2 pointing left

Key Properties

  • Symmetric: $K_{ij} = K_{ji}$ — Maxwell-Betti reciprocal theorem
  • Positive diagonal: $K_{ii} > 0$ — Positive stiffness in direction of DOF
  • Row/column sum = 0: Rigid body motion produces no internal forces

The General Formula

For any element, the stiffness matrix is:

$$[K_e] = \int_V [B]^T [D] [B] \, dV$$

Where:

  • $[B]$ = Strain-displacement matrix (relates nodal displacements to strain)
  • $[D]$ = Material constitutive matrix (relates strain to stress)
  • $V$ = Element volume

For 1D:

$$[B] = \begin{bmatrix} \frac{dN_1}{dx} & \frac{dN_2}{dx} \end{bmatrix} = \begin{bmatrix} -\frac{1}{L} & \frac{1}{L} \end{bmatrix}$$

$$[D] = E \quad \text{(Young's modulus)}$$

$$[K_e] = \int_0^L [B]^T E [B] \cdot A \, dx$$

Working this out:

$$[K_e] = \int_0^L \begin{bmatrix} -1/L \\ 1/L \end{bmatrix} E \begin{bmatrix} -1/L & 1/L \end{bmatrix} A \, dx$$

$$= \frac{EA}{L^2} \int_0^L \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} dx = \frac{EA}{L} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}$$

Same result!

Extending to 2D and 3D

The same formula applies, but $[B]$ and $[D]$ become larger:

2D Plane Stress (3-node triangle)

$$[B] = \frac{1}{2A_e} \begin{bmatrix}

b_1 & 0 & b_2 & 0 & b_3 & 0 \\

0 & c_1 & 0 & c_2 & 0 & c_3 \\

c_1 & b_1 & c_2 & b_2 & c_3 & b_3

\end{bmatrix}$$

Where $b_i$, $c_i$ are geometric constants from node coordinates.

$$[D] = \frac{E}{1-\nu^2} \begin{bmatrix}

1 & \nu & 0 \\

\nu & 1 & 0 \\

0 & 0 & \frac{1-\nu}{2}

\end{bmatrix}$$

The resulting $[K_e]$ is $6 \times 6$ (3 nodes × 2 DOFs each).

Effect of Element Parameters

How does changing $E$, $A$, or $L$ affect stiffness?

ParameterEffect on $K$Physical Meaning
$E$ ↑$K$ ↑Stiffer material → less deformation
$A$ ↑$K$ ↑More cross-section → more material to resist
$L$ ↑$K$ ↓Longer bar → more flexible (like a longer spring)
Stiffness is like axial rigidity divided by length: $K = \frac{EA}{L}$

This matches intuition:

  • Steel ($E = 200$ GPa) is stiffer than aluminum ($E = 70$ GPa)
  • Thicker bars are stiffer
  • Longer bars are more flexible

Force Vector

The right-hand side of $[K]\{u\} = \{F\}$ comes from:

$$F_i = \int_0^L N_i \cdot f \, dx + N_i(L) \cdot F_L$$

For a uniformly distributed load $f$ (force per unit length):

$$F_1 = \int_0^L \left(1 - \frac{x}{L}\right) f \, dx = f \cdot \frac{L}{2}$$

$$F_2 = \int_0^L \frac{x}{L} \cdot f \, dx = f \cdot \frac{L}{2}$$

The distributed load is lumped equally to the two nodes!

For a point load $P$ at the right end:

$$F_1 = N_1(L) \cdot P = 0, \quad F_2 = N_2(L) \cdot P = P$$

Example: Computing Stiffness

Problem: Steel bar, $E = 200$ GPa, $A = 100$ mm², $L = 500$ mm Solution:

$$K = \frac{EA}{L} = \frac{200 \times 10^9 \times 100 \times 10^{-6}}{500 \times 10^{-3}}$$

$$K = \frac{200 \times 10^3 \times 100 \times 10^{-6}}{0.5} = \frac{20}{0.5} = 40 \text{ MN/m} = 40 \text{ kN/mm}$$

$$[K_e] = 40 \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} \text{ kN/mm}$$

If we apply a force of 10 kN at the free end (node 2 fixed):

  • $u_1 = F/K = 10/40 = 0.25$ mm

Key Takeaways

  • Element stiffness matrix relates nodal forces to displacements: $[K_e]\{u\} = \{F\}$
  • Derivation: Apply Galerkin method to weak form → $K_{ij} = \int [B]^T [D] [B] dV$
  • For 2-node bar: $[K_e] = \frac{EA}{L} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}$
  • Physical meaning: $K_{ij}$ = force at $i$ due to unit displacement at $j$
  • Properties: Symmetric, positive diagonal, singular (rigid body modes)
  • Stiffness scales as: $EA/L$ — stiffer material, larger area, shorter length → stiffer element

What's Next

Now we have element stiffness matrices. But a real structure has many elements. In the next lesson, we'll learn how to assemble element matrices into a global system that represents the entire structure.