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."
Two Springs in Series
For two springs connected end-to-end with a shared node:
| Node | Forces |
|---|---|
| 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!
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:
| Entry | Meaning |
|---|---|
| $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?
| Parameter | Effect 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) |
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.