Sources and References

Primary references — D.L. Logan, A First Course in the Finite Element Method, Cengage Learning; D.V. Hutton, Fundamentals of Finite Element Analysis, McGraw-Hill.

Online resources — MIT OCW 2.094 (Finite Element Analysis); Klaus-Jürgen Bathe, Finite Element Procedures, Prentice Hall (2nd ed.); O.C. Zienkiewicz & R.L. Taylor, The Finite Element Method (standard reference, 3 vols.); T.J.R. Hughes, The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover.

Chapter 1: Introduction and the FEM Philosophy

1.1 What Is the Finite Element Method?

The finite element method (FEM) is a general numerical technique for obtaining approximate solutions to boundary-value and initial-value problems in engineering and science. Its core idea is elegant: a complex continuous domain — a loaded truss, a heat-conducting fin, a vibrating turbine blade — is subdivided into a finite number of simple subregions called elements, connected at discrete points called nodes. The governing differential equations, which cannot in general be solved analytically over the original domain, are reformulated as a system of algebraic equations over these elements. Assembling the element equations and imposing boundary conditions yields a global system that can be solved numerically to yield approximate values of the unknown field (displacement, temperature, pressure) at every node.

What sets FEM apart from other discretisation methods — finite differences, boundary element methods, spectral methods — is its natural handling of complex geometries and material heterogeneity, its systematic assembly procedure, and its solid variational foundation. The method was developed independently by structural engineers in the 1950s and applied mathematicians in the 1960s; the name “finite element method” was coined by Clough in 1960. Today it underlies virtually every commercial simulation package (ANSYS, Abaqus, COMSOL, LS-DYNA) and is central to product design across aerospace, civil, biomedical, and mechanical engineering.

1.2 The Three Fundamental Problem Classes

Virtually every physical problem amenable to FEM treatment falls into one of three classes:

Equilibrium (static) problems. The system does not evolve in time; boundary conditions are fixed. The governing equation takes the form \( \mathbf{K}\mathbf{u} = \mathbf{f} \), where \(\mathbf{K}\) is the stiffness (or conductivity) matrix, \(\mathbf{u}\) the nodal unknowns, and \(\mathbf{f}\) the load vector. Examples: static deflection of a loaded beam, steady-state temperature distribution in a fin, stress analysis of a pressurised vessel.

Propagation (transient) problems. The field evolves in time. The semi-discrete FEM equation is

where \(\mathbf{M}\) is the mass matrix, \(\mathbf{C}\) the damping matrix, and \(\mathbf{f}(t)\) a time-varying load. Examples: transient heat conduction, structural dynamics under impact, wave propagation.

Eigenvalue problems. One seeks natural frequencies or critical loads without specifying a particular forcing function. For free vibration:

The eigenvalues \(\omega^2\) are squared natural frequencies and the eigenvectors \(\boldsymbol{\phi}\) are mode shapes. For buckling, \(\mathbf{K}_G\) (geometric stiffness) replaces \(\omega^2 \mathbf{M}\).

1.3 Review of Relevant Mechanics

1.3.1 Stress and Strain Notation

In three dimensions the stress state at a point is described by six independent components:

The engineering strain vector is similarly:

where \(\gamma_{ij} = 2\varepsilon_{ij}\) is the engineering shear strain. For a linear elastic material (generalised Hooke’s law):

where \(\mathbf{D}\) is the \(6 \times 6\) constitutive (elasticity) matrix. For an isotropic material it depends only on Young’s modulus \(E\) and Poisson’s ratio \(\nu\).

1.3.2 Strain–Displacement Relations

For small deformations, the strain–displacement relations are:

In matrix form \(\boldsymbol{\varepsilon} = \mathbf{B}\mathbf{u}\), where \(\mathbf{B}\) is the strain–displacement matrix derived by differentiating the shape functions.

1.3.3 Virtual Work and Equilibrium

The principle of virtual work states that equilibrium is satisfied if, for every kinematically admissible virtual displacement \(\delta\mathbf{u}\):

where \(\mathbf{b}\) is the body force per unit volume and \(\mathbf{t}\) is the surface traction. This is the weak form from which FEM derives its element equations.

1.4 Overview of the FEM Pipeline

The complete FEM workflow, which this text develops chapter by chapter, consists of:

1. Pre-processing: Define geometry, material properties, element type, mesh.
2. Element formulation: Derive shape functions; form element stiffness \(\mathbf{k}^{(e)}\), mass \(\mathbf{m}^{(e)}\), and load \(\mathbf{f}^{(e)}\).
3. Assembly: Assemble global \(\mathbf{K}\), \(\mathbf{M}\), \(\mathbf{F}\) from element contributions.
4. Boundary conditions: Apply essential (displacement) BCs by elimination or penalty.
5. Solution: Solve \(\mathbf{K}\mathbf{u} = \mathbf{F}\) (or the appropriate eigenvalue/transient system).
6. Post-processing: Recover stresses, strains, reactions; visualise results.

Chapter 2: Spring Elements and the Direct Stiffness Method

2.1 The Spring Element

The simplest finite element is the axial spring. A single spring element has two nodes — labelled \(i\) and \(j\) — each with one degree of freedom (DOF): the axial displacement, \(u_i\) and \(u_j\). The spring has stiffness constant \(k\) (force per unit displacement).

The element equilibrium equation is therefore:

Note that \(\mathbf{k}^{(e)}\) is singular because a rigid body displacement (both nodes moving equally) produces no force. The matrix becomes non-singular only after boundary conditions are applied.

2.2 The Direct Stiffness Method

The direct stiffness method assembles the global stiffness matrix by placing each element stiffness coefficient into the appropriate position determined by the element’s node numbering. If element \((e)\) connects global node \(p\) to global node \(q\), then:

2.3 Applying Boundary Conditions

After assembly, the global system \(\mathbf{K}\mathbf{U} = \mathbf{F}\) is singular. Boundary conditions (BCs) are imposed to remove rigid-body modes:

Elimination method: Rows and columns corresponding to known (zero or specified) displacements are removed from the system, reducing it to a non-singular reduced system.

Penalty method: Add a large number \(\alpha \gg K_{max}\) to the diagonal entry corresponding to the constrained DOF, and set the corresponding load to \(\alpha u_{prescribed}\). This is less elegant numerically but convenient for implementation.

2.4 Multiple Degrees of Freedom and Symmetry

The global stiffness matrix \(\mathbf{K}\) is always:

• Symmetric: \(K_{ij} = K_{ji}\) (follows from the symmetry of each \(\mathbf{k}^{(e)}\)).
• Positive semi-definite before BCs; positive definite after sufficient BCs.
• Banded (sparse): Non-zero entries are confined to a band around the diagonal. The half-bandwidth equals the maximum difference in node numbers sharing an element, plus one. Efficient solvers exploit this sparsity.

Chapter 3: Bar and Truss Elements

3.1 The Bar Element

A bar (rod) element resists only axial force. Unlike the spring element (which models a concentrated stiffness), the bar element has distributed stiffness governed by its cross-sectional area \(A\), Young’s modulus \(E\), and length \(L\).

3.1.1 Shape Functions for the Bar Element

In matrix form: \(u(x) = \mathbf{N}(x)\,\mathbf{u}^{(e)}\) where \(\mathbf{N} = \left[N_i,\; N_j\right]\).

3.1.2 Strain–Displacement for the Bar

The axial strain is:

For the linear bar:

The axial stress follows from Hooke’s law: \(\sigma = E\varepsilon = E\mathbf{B}\mathbf{u}^{(e)}\).

3.1.3 Bar Element Stiffness Matrix from Minimum Potential Energy

3.2 Truss Elements

A truss element is a bar element oriented at an arbitrary angle in 2D or 3D space. Each node has two DOFs in 2D (\(u\), \(v\)) or three in 3D (\(u\), \(v\), \(w\)).

3.2.1 Local and Global Coordinate Systems

Let the element make angle \(\theta\) with the global \(x\)-axis. The local axial DOFs \(\bar{u}_i\), \(\bar{u}_j\) are related to global DOFs by the transformation:

In matrix form: \(\bar{\mathbf{u}}^{(e)} = \mathbf{T}\mathbf{u}^{(e)}\) where the transformation (rotation) matrix is:

3.2.2 Element Stiffness in Global Coordinates

Since \(\mathbf{k}^{(e)} = \frac{EA}{L}\left[\begin{smallmatrix}1&-1\\-1&1\end{smallmatrix}\right]\) in local coordinates, the global element stiffness is:

This \(4 \times 4\) matrix is the standard plane truss element stiffness in global coordinates.

3.2.3 Elemental Stress and Strain Recovery

After solving for global nodal displacements \(\mathbf{U}\), element displacements in local coordinates are recovered:

The axial strain and stress in element \((e)\) are:

Chapter 4: Beam and Frame Elements

4.1 Euler–Bernoulli Beam Theory

The Euler–Bernoulli beam theory assumes that plane sections remain plane and perpendicular to the neutral axis after bending (no shear deformation). For a beam with bending stiffness \(EI\), the governing differential equation is:

where \(v(x)\) is the transverse displacement and \(q(x)\) is the distributed transverse load. The bending moment is \(M = EI \, d^2v/dx^2\) and the shear force is \(V = EI \, d^3v/dx^3\).

4.2 Beam Element DOFs and Hermite Shape Functions

A two-node beam element has four DOFs: transverse displacement \(v_i\), rotation \(\theta_i\) at node \(i\), and \(v_j\), \(\theta_j\) at node \(j\). Let \(\theta = dv/dx\). The element DOF vector is \(\mathbf{d}^{(e)} = \left[v_i,\; \theta_i,\; v_j,\; \theta_j\right]^T\).

4.3 Beam Element Stiffness Matrix

The curvature–displacement relation is \(\kappa = d^2v/dx^2 = \mathbf{B}_{beam}\mathbf{d}^{(e)}\) where:

The element stiffness matrix is:

This is the Euler–Bernoulli beam stiffness matrix, one of the most important results in structural FEM.

4.4 Consistent Nodal Loads for Distributed Loads

For a uniform distributed load \(q_0\) over the full length \(L\):

The end moments \(q_0 L^2/12\) are the “fixed-end moments” familiar from structural analysis.

4.5 Frame Elements

A plane frame element combines bar (axial) DOFs with beam (bending) DOFs. Each node has three DOFs: \(u\) (axial), \(v\) (transverse), and \(\theta\) (rotation). In local coordinates the element stiffness is block diagonal:

where the \(2\times2\) bar stiffness and \(4\times4\) beam stiffness are as derived above. In global coordinates, a \(3\times3\) rotation matrix transforms each node’s DOFs, giving a \(6\times6\) global frame element stiffness matrix.

The rotation matrix for each node is:

and the full transformation is \(\mathbf{T} = \text{diag}(\mathbf{R}, \mathbf{R})\), so \(\mathbf{K}^{(e)} = \mathbf{T}^T \mathbf{k}_{local}^{(e)} \mathbf{T}\).

Chapter 5: Galerkin Weighted Residual Method

5.1 Weighted Residual Methods: Overview

The direct stiffness method and energy approaches are intuitive for structural problems where an energy functional is available. For problems without a clear energy principle — such as convection-dominated heat transfer, fluid flow, or electromagnetic fields — the method of weighted residuals provides a more general framework.

Given a differential equation \(\mathcal{L}(u) = f\) on domain \(\Omega\) with BCs on \(\partial\Omega\), an approximate solution \(\tilde{u} = \sum_j c_j \phi_j\) is sought. The residual is \(R = \mathcal{L}(\tilde{u}) - f \neq 0\). The weighted residual condition demands:

Different choices of weight functions \(w_i\) give different methods: Galerkin (\(w_i = \phi_i\)), Petrov-Galerkin (different test and trial spaces), collocation (Dirac deltas), least squares, and subdomain methods.

5.2 The Galerkin Approach

The Galerkin method sets the weight functions equal to the approximation basis functions: \(w_i = N_i\), the same shape functions used to interpolate the unknown. This choice is not arbitrary — it produces symmetric stiffness matrices for self-adjoint operators and has strong variational justification.

5.3 Galerkin Weak Form for the Bar Element

Consider the bar element governed by:

with essential BC \(u(0) = 0\) and natural BC \(EA \, du/dx|_{x=L} = P\) (applied force).

5.4 Galerkin Method Applied to 1D Heat Conduction

The steady-state heat equation with internal generation \(Q\) is:

Chapter 6: Heat Transfer by FEM

6.1 1D Heat Conduction: Element Formulation

For a two-node element \([x_i, x_j]\) with linear shape functions, the element conductivity matrix is:

where \(k\) is thermal conductivity (W/m·K) and \(A\) is cross-sectional area. Note the structural identity with the bar stiffness matrix, with \(EA/L \to kA/L\).

The heat source (load) vector for uniform generation \(Q\) (W/m³):

6.1.1 Boundary Conditions in Heat Transfer FEM

Three types of thermal boundary conditions arise:

1. Essential (Dirichlet): Prescribed temperature \(T = T_0\). Applied by the elimination or penalty method, exactly as for displacement BCs.

2. Natural (Neumann): Prescribed heat flux \(q_n = -k \, dT/dn\). Appears in the natural boundary condition of the weak form; simply add \(q_n A\) to the nodal load.

3. Convection (Robin): \(q_n = h(T - T_\infty)\) on the boundary. Adds to both the stiffness matrix (convective conductance) and the load vector (heat input from ambient).

For a convection boundary at node \(j\), the additional terms are:

where \(A_s\) is the surface area at the node.

6.2 Transient Heat Transfer

For transient problems \((\rho c \, \partial T/\partial t = \nabla \cdot (k \nabla T) + Q)\), the FEM discretisation yields:

where \(\mathbf{C}\) is the heat capacity matrix with element entries \(C_{ij}^{(e)} = \int \rho c \, N_i N_j \, dV\), giving for linear elements:

Time integration is performed by finite differences (backward Euler, Crank–Nicolson, or higher-order methods).

% Assemble conductivity matrix
K = zeros(n_nodes);
F = zeros(n_nodes, 1);
for e = 1:n_elements
    nodes = connectivity(e, :);  % [i, j]
    k_e = kA/L_e * [1 -1; -1 1];
    K(nodes, nodes) = K(nodes, nodes) + k_e;
    F(nodes) = F(nodes) + Q*A*L_e/2 * [1; 1];
end
% Apply convection at right boundary
K(end,end) = K(end,end) + h*As;
F(end) = F(end) + h*As*T_inf;
% Apply essential BC: T(1) = T_prescribed
K_red = K(2:end, 2:end);
F_red = F(2:end) - K(2:end, 1)*T_prescribed;
T_free = K_red \ F_red;

Chapter 7: Isoparametric Formulation and Higher-Order Elements

7.1 Motivation for Isoparametric Elements

Real structures have curved boundaries, non-uniform element sizes, and distorted shapes. The isoparametric formulation elegantly handles these complications by using the same shape functions to interpolate both the geometry and the field variable.

7.2 Natural Coordinates for 1D Elements

For a 1D bar element, the natural coordinate \(\xi \in [-1, 1]\) maps to physical coordinate:

The Jacobian is \(J = (x_j - x_i)/2 = L/2\).

Shape functions in natural coordinates:

Strain–displacement matrix:

— identical to before, confirming consistency.

7.3 Higher-Order 1D Elements

7.4 Gaussian Quadrature

Exact integration of complex expressions (products of shape function derivatives and the Jacobian) is rarely possible analytically. Gaussian quadrature provides exact integration for polynomials of degree \(2n-1\) using only \(n\) integration points.

7.5 2D Quadrilateral Elements

For a bilinear four-node quadrilateral (Q4) element with nodes in natural coordinates \((\xi, \eta) \in [-1,1]^2\):

The Jacobian matrix is:

The strain–displacement matrix \(\mathbf{B}\) in physical coordinates is obtained by inverting \(\mathbf{J}\):

The stiffness integral is evaluated numerically:

where \(t\) is the element thickness.

7.6 Triangular Elements: Area Coordinates

For triangular elements the natural coordinates are area (barycentric) coordinates \(L_1, L_2, L_3\) satisfying \(L_1 + L_2 + L_3 = 1\).

Chapter 8: 2D Plane Stress and Plane Strain

8.1 2D Stress State Assumptions

Many engineering structures — thin plates, long dams, pressure vessels in cross-section — can be analysed accurately using 2D simplifications of the full 3D problem.

8.1.1 Constitutive Matrices for Plane Problems

Plane stress constitutive matrix:

Plane strain constitutive matrix:

Note: for plane strain, replace \(E \to E/(1-\nu^2)\) and \(\nu \to \nu/(1-\nu)\) in the plane stress matrix, or equivalently use the expression above.

8.2 The Constant Strain Triangle (CST)

The CST is the simplest 2D solid element. It has three nodes, each with two DOFs (\(u\), \(v\)), giving six DOFs per element.

8.2.1 CST Shape Functions from Area Coordinates

8.2.2 CST Strain–Displacement Matrix

Since the shape functions are linear in \(x\) and \(y\), their derivatives are constant — hence the element name “Constant Strain Triangle.”

The strain vector is:

The \(3 \times 6\) strain–displacement matrix is:

where \(b_i = y_j - y_k\), \(c_i = x_k - x_j\) (with \(i,j,k\) cyclic permutations of \(1,2,3\)).

8.2.3 CST Element Stiffness Matrix

Since \(\mathbf{B}\) is constant over the element, the stiffness matrix integrates trivially:

This is a \(6 \times 6\) symmetric matrix — the complete element stiffness for the CST.

8.3 Distributed Surface Loads on Triangular Elements

8.4 2D Plane Stress FEM: Complete Example

8.5 Strain Energy and von Mises Criterion

The elastic strain energy density for a plane stress element is:

The von Mises (distortion energy) failure criterion is:

In commercial FEM codes, post-processing computes \(\sigma_{VM}\) at each Gauss point or node and contours it over the mesh. Regions where \(\sigma_{VM} > \sigma_y\) indicate potential yielding.

Chapter 9: 3D Solid and Axisymmetric Elements

9.1 Tetrahedral Elements

The 4-node tetrahedron (Tet4) is the 3D analogue of the CST. Each of four nodes has three DOFs (\(u\), \(v\), \(w\)), giving 12 DOFs per element. Linear shape functions in volume coordinates \(L_1, L_2, L_3, L_4\) (\(L_i \geq 0\), \(\sum L_i = 1\)) give constant strains throughout the element.

The displacement field:

where \(N_i = L_i = (a_i + b_i x + c_i y + d_i z)/(6V)\), with \(V\) the element volume and the coefficients derived from node coordinates.

The \(6 \times 12\) strain–displacement matrix \(\mathbf{B}\) is constant (linear field → constant derivatives), and the element stiffness is:

where \(\mathbf{D}_{3D}\) is the full \(6 \times 6\) isotropic elastic constitutive matrix:

9.1.1 Higher-Order Tetrahedral Elements (Tet10)

The 10-node quadratic tetrahedron (Tet10) adds mid-side nodes to each of the six edges. Shape functions are quadratic in volume coordinates. The strain field is now linear, providing much better accuracy for stress gradients. Tet10 is the default solid element in many commercial codes for unstructured meshes.

9.2 Hexahedral Elements

The 8-node hexahedron (Hex8 or “brick element”) has trilinear shape functions in natural coordinates \((\xi, \eta, \zeta) \in [-1,1]^3\):

where \((\xi_i, \eta_i, \zeta_i) = \pm 1\) defines the node positions. The Hex8 element has 24 DOFs, a linearly varying displacement field, and uses \(2 \times 2 \times 2\) Gauss quadrature (8 integration points) for full integration.

The 20-node serendipity hexahedron (Hex20) adds mid-edge nodes and uses quadratic shape functions — significantly more accurate for curved geometries and large stress gradients.

9.3 Axisymmetric Elements

For structures with rotational symmetry (pressure vessels, flywheels, circular foundations), 3D analysis reduces to a 2D problem in the \(r\)–\(z\) plane.

The strain vector for an axisymmetric problem has four components:

where the hoop strain \(\varepsilon_{\theta\theta} = u_r / r\) (not a derivative) requires special treatment.

The axisymmetric strain–displacement matrix for a 3-node triangular element in the \(r\)–\(z\) plane is:

The stiffness integral includes the circumference factor \(2\pi r\):

Because \(\mathbf{B}_{axi}\) is not constant (\(N_i/r\) varies), numerical integration (Gauss quadrature in the \(r\)–\(z\) plane) is required.

9.4 Rectangular Shell Elements

For thin shell structures, the shell element combines membrane (in-plane stretching) and bending (out-of-plane) behaviour. The simplest approach is to superpose a plane stress membrane element with a plate bending element — this gives a “flat shell” element.

Each node of a flat shell element has five or six DOFs: \(u, v, w, \theta_x, \theta_y\) (and sometimes a drilling rotation \(\theta_z\) for numerical stability). The stiffness matrix is assembled from membrane and bending contributions.

For a rectangular flat shell element of dimensions \(a \times b\) with Kirchhoff plate theory:

The membrane stiffness uses bilinear shape functions (Q4 plane stress). The bending stiffness uses the ACM (Adini–Clough–Melosh) plate element or Mindlin–Reissner formulation (preferred for thick shells).

Chapter 10: Mesh Design, Convergence, and Verification

10.1 Convergence Theory

The FEM approximate solution \(u_h\) converges to the exact solution \(u\) as the mesh is refined (element size \(h \to 0\)). Convergence requires:

1. Completeness: The displacement field must be able to represent rigid body modes and constant strain states exactly. This is guaranteed if the shape functions form a complete polynomial of degree at least \(p\).

2. Compatibility (conforming elements): Displacements must be continuous across element boundaries. For C0 elements (bar, CST, Q4), continuity of \(u\) is required; for C1 elements (Euler–Bernoulli beam), continuity of \(u\) and \(du/dx\).

The convergence rate in energy norm is:

where \(p\) is the polynomial degree of the shape functions and \(h\) is the characteristic element size. For the energy (strain energy):

The strain energy converges from below (the FEM solution is stiffer than reality) for conforming elements.

10.2 h-Refinement, p-Refinement, and hp-Methods

h-refinement: Reduce element size \(h\) while keeping polynomial order \(p\) fixed. Convergence rate \(O(h^p)\). Most common approach in practice.

p-refinement: Increase polynomial order \(p\) while keeping mesh fixed. Exponential convergence for smooth solutions (Babuska and Szabo). Useful for problems where the solution is globally smooth.

hp-methods: Combine both. Optimal convergence for problems with singularities (crack tips, re-entrant corners) when \(h\) is reduced near the singularity and \(p\) is increased elsewhere.

10.3 Patch Test

The patch test (due to Irons and Razzaque, 1965) is a necessary and sufficient condition for convergence of an element. It requires that when the exact solution is a polynomial of degree \(p\) (corresponding to the element’s completeness requirement), the FEM solution must reproduce it exactly on an arbitrary patch of elements.

Procedure: Assemble a patch of elements with arbitrary nodal coordinates. Apply nodal displacements corresponding to a state of constant strain. The element stresses should equal the exact constant stress everywhere. If the patch test is passed, the element converges.

10.4 Stress Recovery and Superconvergence

The raw FEM stress field (computed as \(\boldsymbol{\sigma} = \mathbf{D}\mathbf{B}\mathbf{d}\)) is discontinuous across element boundaries because \(\mathbf{B}\) is element-wise constant (for CST) or piecewise polynomial. Several stress recovery techniques improve accuracy:

Stress averaging: Average stresses from all elements sharing a node. Simple but first-order.

Superconvergent patch recovery (SPR): Due to Zienkiewicz and Zhu (1987). Fit a polynomial to stresses at superconvergent (Barlow) points within a patch of elements around each node. For Q8 elements, the optimal sampling points are the \(2 \times 2\) Gauss points. SPR gives a recovered stress field one order higher in accuracy than the raw field.

Error estimation: The difference between recovered and raw stresses gives a local error indicator for adaptive mesh refinement.

10.5 Sources of Error in FEM

FEM solutions contain three main sources of error:

1. Discretisation error: From approximating the geometry and field with a finite mesh. Controlled by \(h\)- or \(p\)-refinement.

2. Integration error: From Gaussian quadrature with insufficient integration points. For linear elements (Q4), full \(2 \times 2\) integration is exact for the bilinear stiffness integrand. Under-integration (\(1 \times 1\) for Q4) introduces “spurious zero-energy modes” (hourglass modes) that pollute the solution unless stabilised.

3. Equation solving error: Floating-point round-off in the linear solver. The condition number of \(\mathbf{K}\) grows as \(O(h^{-2})\) for 2D elasticity, eventually limiting accuracy. Preconditioning is used in large-scale problems.

10.6 Locking Phenomena

Shear locking: In Euler–Bernoulli beam or thin shell elements formulated with displacement-based C0 elements (rather than the exact Hermite C1 formulation), artificial shear strains arise as thickness \(t \to 0\). The element becomes over-stiff. Remedies: reduced/selective integration; assumed natural strain (ANS) formulation; MITC (Mixed Interpolation of Tensorial Components) elements.

Volumetric locking: In nearly incompressible materials (\(\nu \to 0.5\)), the pressure field is poorly represented by polynomial interpolation, causing the mesh to “lock” (spuriously high stiffness). Remedies: mixed \(u/p\) formulations; reduced integration; B-bar method.

10.7 Practical Mesh Guidelines

The following heuristics are drawn from industrial practice and standard FEM textbook guidance (Logan, Hutton, Bathe):

• Element aspect ratio: Keep \(a/b \leq 4\) for quads, \(\leq 3\) for triangles. High aspect ratios degrade accuracy severely near stress concentrations.
• Element distortion: Avoid highly skewed elements (interior angles outside 45°–135°). The Jacobian determinant \(\det\mathbf{J}\) must remain positive everywhere; near-zero \(\det\mathbf{J}\) signals a degenerate element.
• Mesh density at stress concentrations: Use at least 4–6 elements per fillet radius. A Saint-Venant decay length of \(\sim 3 \times\) the characteristic dimension guides transition from fine to coarse mesh.
• First-order vs. second-order elements: Second-order elements (Tet10, Hex20, Q8) are generally preferred in modern practice. They represent curved boundaries exactly and reduce the number of elements needed for the same accuracy.
• Symmetry: Exploit symmetry planes to reduce problem size. Symmetric geometry + symmetric loading → symmetric solution: apply zero-normal-displacement and zero-tangential-traction BCs on symmetry planes.

Appendix: Complete FEM Reference Tables

A.1 Summary of Element Stiffness Matrices

A.2 Gauss Point Locations and Weights

A.3 Common Boundary Condition Types

A.4 Key FEM Equations Summary

Element stiffness (general):

Consistent mass matrix:

Consistent nodal load vector:

Global equilibrium:

Stress recovery:

Von Mises stress (plane stress):

Supplementary Examples

Key Formula Reference

Bar Element

Beam Element (Euler–Bernoulli)

Plane Truss Element (Global)

CST Element (Plane Stress/Strain)

1D Heat Conduction Element