Mathematics.

numerical analysis

Finite Element Method

Numerical Analysis50 minDifficulty9 out of 10

You should know: numerical methods, partial differential equations, linear transformation

Overview

The finite element method (FEM) is the dominant numerical technique for solving the partial differential equations of continuum mechanics — stress and deflection in structures, heat flow, seepage, and vibration. Its central idea is discretization: a complicated body is subdivided into many small, simply-shaped pieces (the 'finite elements' — bars, triangles, quadrilaterals, tetrahedra) connected at nodes. On each element the unknown field is approximated by simple interpolation (shape) functions, the governing PDE is enforced in an averaged 'weak' sense, and each element contributes a small stiffness matrix. Assembling all of them produces one large linear system K·d = F, whose solution gives the nodal displacements (or temperatures, or heads) of the whole structure. FEM is what powers essentially every modern structural-analysis and CAE package (SAP2000, ETABS, ABAQUS, ANSYS).

Intuition

You cannot write a single formula for the stress in a bridge with a complicated shape and loading — but you can for a single straight bar. FEM exploits exactly that: cut the hard object into many easy pieces, describe each piece with an elementary model (a spring/stiffness relation), and then reconnect them by insisting that shared nodes move together and that forces balance at every node. The 'weak form' is the mathematical trick that makes this rigorous: instead of demanding the PDE hold exactly at every point (which the simple piecewise approximations can't do), you require it to hold on average against a set of test functions. The payoff is that a differential equation over a continuous body becomes a matrix equation K·d = F — algebra a computer can crush.

Formal Definition

Definition

The heart of FEM is the element stiffness relation and the assembled global system. For a linear-elastic axial bar element of length L, area A, modulus E, the 2-node stiffness matrix and the global system are:

ke=AEL[1−1−11]\mathbf{k}^e = \dfrac{AE}{L}\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}
Axial bar element stiffness matrix
K d=F\mathbf{K}\,\mathbf{d} = \mathbf{F}
Assembled global system (stiffness × displacement = force)
ue(x)=N1(x) d1+N2(x) d2u^e(x) = N_1(x)\,d_1 + N_2(x)\,d_2
Field interpolated by shape functions Nᵢ within an element
N1(x)=1−xL,N2(x)=xLN_1(x) = 1-\dfrac{x}{L},\qquad N_2(x)=\dfrac{x}{L}
Linear (2-node) shape functions

Notation

NotationMeaning
ke\mathbf{k}^eElement stiffness matrix — relates one element's nodal forces to its nodal displacements
K\mathbf{K}Global (assembled) stiffness matrix for the whole structure
d\mathbf{d}Vector of unknown nodal displacements (degrees of freedom)
F\mathbf{F}Vector of nodal forces (applied loads and reactions)
NiN_iShape / interpolation function — 1 at its own node, 0 at the others

Derivation

A 2-node axial bar illustrates where kᵉ comes from. The bar behaves as a linear spring of stiffness AE/L; equilibrium of its two nodal forces in terms of the two nodal displacements gives the 2×2 element matrix, which then embeds into the global K by adding contributions at shared degrees of freedom:

k=AELk = \dfrac{AE}{L}

A uniform elastic bar acts as a spring with stiffness AE/L (force per unit stretch).

f1=k(d1−d2),f2=k(d2−d1)f_1 = k(d_1 - d_2),\qquad f_2 = k(d_2 - d_1)

Newton's third law at the two nodes: internal force depends on the relative displacement (the stretch).

[f1f2]=AEL[1−1−11][d1d2]\begin{bmatrix} f_1 \\ f_2 \end{bmatrix} = \dfrac{AE}{L}\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}\begin{bmatrix} d_1 \\ d_2 \end{bmatrix}

Write the two equations in matrix form — this is kᵉ.

K=∑eke  (assembled at shared DOFs)\mathbf{K} = \sum_e \mathbf{k}^e \;(\text{assembled at shared DOFs})

Overlap element matrices at common nodes to build the global stiffness matrix; apply supports, then solve K d = F.

Properties

Symmetry of K

K=KT\mathbf{K} = \mathbf{K}^{T}

Condition: Linear elastic (conservative) system

Example: Reciprocity (Maxwell–Betti): the deflection at A from a load at B equals the deflection at B from the same load at A.

Positive definiteness (after supports)

dTK d>0 for d≠0\mathbf{d}^{T}\mathbf{K}\,\mathbf{d} > 0 \text{ for } \mathbf{d} \ne \mathbf{0}

Condition: Enough restraints to remove rigid-body motion

Example: Guarantees a unique displacement solution; the quadratic form equals twice the stored strain energy.

Convergence with refinement

∥u−uh∥→0 as element size h→0\|u - u_h\| \to 0 \text{ as element size } h \to 0

Condition: Compatible, complete shape functions (patch test passed)

Example: A finer mesh gives a more accurate stress field, approaching the exact PDE solution.

Singularity of the free structure

det⁡(Kunrestrained)=0\det(\mathbf{K}_{\text{unrestrained}}) = 0

Condition: Before boundary conditions are applied

Example: An unsupported structure can float (rigid-body modes); you must fix supports before K can be inverted.

Applications

Structural analysis of frames, trusses, plates, and shells: the direct stiffness method — assembling element matrices into K d = F — is the computational engine of SAP2000, ETABS, STAAD, and RISA for buildings and bridges.

Worked Examples

  1. Compute the axial stiffness AE/L. Use consistent units: A = 500 mm², E = 200 kN/mm² (since 200 GPa = 200 kN/mm²), L = 2000 mm.

    AEL=500×2002000=50 kN/mm\dfrac{AE}{L} = \dfrac{500 \times 200}{2000} = 50\ \text{kN/mm}

  2. Insert AE/L into the standard bar template.

    ke=50[1−1−11] kN/mm\mathbf{k}^e = 50\begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix}\ \text{kN/mm}

  3. So the element matrix has +50 on the diagonal and −50 off-diagonal.

    ke=[50−50−5050] kN/mm\mathbf{k}^e = \begin{bmatrix} 50 & -50 \\ -50 & 50 \end{bmatrix}\ \text{kN/mm}

Answer: kᵉ = [[50, −50], [−50, 50]] kN/mm (AE/L = 50 kN/mm).

Practice Problems

Difficulty 5/10

A truss member has E = 70 GPa (aluminium), A = 300 mm², L = 1.5 m. Compute its axial stiffness AE/L in kN/mm.

Difficulty 5/10

A 2-D plane truss has 8 joints (nodes). Before applying supports, how many degrees of freedom does the global stiffness matrix have, and what is its size?

Difficulty 6/10

Why is the global stiffness matrix K singular (non-invertible) BEFORE boundary conditions are applied?

Difficulty 6/10

A linear shape function N₁(x) for a 2-node bar equals 1 at node 1. What is its value at node 2?

Difficulty 7/10

State one physical consequence of the stiffness matrix being symmetric (K = Kᵀ) for an elastic structure.

Difficulty 8/10

In seismic design, FEM assembles stiffness K and mass M matrices. What mathematical problem yields the natural frequencies, and what do the solutions represent?

Common Mistakes

Common Mistake

Trying to invert the global stiffness matrix before applying supports.

The unrestrained K is singular (rigid-body modes). You must impose boundary conditions — delete or penalize the restrained DOFs — before solving K d = F.

Common Mistake

Mixing unit systems when computing AE/L or assembling K.

Keep units consistent throughout: e.g. E in kN/mm² (= GPa), A in mm², L in mm gives stiffness in kN/mm. A single unit slip corrupts the whole matrix.

Common Mistake

Assuming a finer mesh is always necessary everywhere.

Refine where gradients are steep (near holes, corners, point loads); a uniform ultra-fine mesh wastes computation. Convergence should be checked, not assumed.

Common Mistake

Forgetting that FEM gives an approximate solution.

FEM solves the weak form on a finite mesh; stresses in particular converge more slowly than displacements, so results must be checked against equilibrium and mesh-refinement studies.

Quiz

What does the finite element method fundamentally convert a continuum PDE into?
The axial stiffness of a 2-node bar element is:
Why must supports (boundary conditions) be applied before solving K·d = F?
In modal/seismic FEM analysis, natural frequencies come from solving:

Flashcards

1 / 5

Summary

  • FEM solves continuum PDEs by discretizing a body into simple elements connected at nodes.
  • On each element the field is approximated with shape functions and the PDE is enforced in a weak (averaged) sense.
  • Each element contributes a stiffness matrix; assembling them gives the global system K·d = F.
  • The global stiffness matrix is symmetric and, after supports are applied, positive-definite and invertible.
  • FEM underpins commercial structural, geotechnical, thermal, and seismic analysis software; refining the mesh converges to the exact solution.

References

  1. BookHughes, T. J. R. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Dover, 2000.
  2. BookLogan, D. L. A First Course in the Finite Element Method, 6th ed. Cengage, 2016.
  3. BookBathe, K.-J. Finite Element Procedures, 2nd ed. Prentice Hall, 2014.
  4. BookCook, R. D. et al. Concepts and Applications of Finite Element Analysis, 4th ed. Wiley, 2001.
  5. WebsiteWikipedia — Finite element method