The Mathematics of Linear FEA

Understanding Linear FEA: The Fundamental Equation

Finite Element Analysis (FEA) for linear problems is built on a beautifully simple mathematical foundation:

$$\mathbf{K} \mathbf{u} = \mathbf{f}$$

This equation represents the equilibrium of a discretized structure, where:

• $\mathbf{K}$ is the global stiffness matrix
• $\mathbf{u}$ is the displacement vector
• $\mathbf{f}$ is the force (load) vector

---

The Stiffness Matrix (K)

The stiffness matrix $\mathbf{K}$ describes how the structure resists deformation. For structural problems, it relates applied forces to resulting displacements.

Building the Global Stiffness Matrix

The global stiffness matrix is assembled from individual element stiffness matrices:

$$\mathbf{K} = \sum_{e=1}^{n_{elem}} \mathbf{K}^{(e)}$$

Each element stiffness matrix $\mathbf{K}^{(e)}$ is computed as:

$$\mathbf{K}^{(e)} = \int_{V^{(e)}} \mathbf{B}^T \mathbf{D} \mathbf{B} \, dV$$

Where:

• $\mathbf{B}$ is the strain-displacement matrix (relates nodal displacements to strains)
• $\mathbf{D}$ is the material constitutive matrix (relates stresses to strains via material properties)
• $V^{(e)}$ is the volume of element $e$

Material Properties in the Stiffness Matrix

For linear elastic materials, the constitutive matrix $\mathbf{D}$ depends on:

• Young's Modulus ($E$) - material stiffness
• Poisson's Ratio ($\nu$) - lateral strain response

For 3D isotropic elasticity:

$$\mathbf{D} = \frac{E}{(1+\nu)(1-2\nu)} \begin{bmatrix} 1-\nu & \nu & \nu & 0 & 0 & 0 \\ \nu & 1-\nu & \nu & 0 & 0 & 0 \\ \nu & \nu & 1-\nu & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1-2\nu}{2} & 0 & 0 \\ 0 & 0 & 0 & 0 & \frac{1-2\nu}{2} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1-2\nu}{2} \end{bmatrix}$$

---

How Mesh Elements Affect the Stiffness Matrix

Element Shape Functions

Each element type (triangles, quads, tetrahedra, hexahedra) uses shape functions $\mathbf{N}$ to interpolate displacements within the element:

$$\mathbf{u}(x,y,z) = \mathbf{N}(x,y,z) \mathbf{u}^{(e)}$$

The strain-displacement matrix $\mathbf{B}$ is derived from these shape functions:

$$\mathbf{B} = \frac{\partial \mathbf{N}}{\partial x}, \frac{\partial \mathbf{N}}{\partial y}, \frac{\partial \mathbf{N}}{\partial z}$$

Mesh Refinement Impact

Finer meshes (more elements, smaller size):

• Increase the size of $\mathbf{K}$ (more degrees of freedom)
• Improve solution accuracy by better approximating curved geometries and stress gradients
• Capture localized effects like stress concentrations
• Increase computational cost (larger system to solve)

Coarser meshes (fewer elements, larger size):

• Reduce computational cost
• May miss important local behavior
• Can produce inaccurate results in high-gradient regions

Element Quality Matters

Poor element shapes (highly skewed or distorted) lead to:

• Ill-conditioned stiffness matrices
• Numerical errors in integration
• Inaccurate stress and strain calculations

---

The Force Vector (f)

The force vector $\mathbf{f}$ contains all external loads applied to the structure:

$$\mathbf{f} = \mathbf{f}_{point} + \mathbf{f}_{distributed} + \mathbf{f}_{body}$$

Types of Loads

Point loads are directly assigned to nodes:

$$\mathbf{f}_{point} = [0, 0, ..., F_x, F_y, F_z, ..., 0, 0]^T$$

Distributed loads (pressure, surface tractions) are converted to nodal forces via integration:

$$\mathbf{f}_{distributed}^{(e)} = \int_{S^{(e)}} \mathbf{N}^T \mathbf{p} \, dS$$

Where $\mathbf{p}$ is the pressure/traction on surface $S^{(e)}$.

Body forces (gravity, acceleration) are volume loads:

$$\mathbf{f}_{body}^{(e)} = \int_{V^{(e)}} \mathbf{N}^T \rho \mathbf{g} \, dV$$

Where $\rho$ is density and $\mathbf{g}$ is gravitational acceleration.

---

The Displacement Vector (u)

The displacement vector $\mathbf{u}$ contains the unknowns we're solving for:

$$\mathbf{u} = [u_1, v_1, w_1, u_2, v_2, w_2, ..., u_n, v_n, w_n]^T$$

For a 3D problem with $n$ nodes:

• $u_i, v_i, w_i$ are displacements in x, y, z directions at node $i$
• Total degrees of freedom = $3n$ (before applying boundary conditions)

Boundary Conditions

Essential (Dirichlet) boundary conditions fix certain displacements:

• Fixed support: $u = v = w = 0$
• Roller support: one or two components set to zero

These constraints are enforced by modifying the system $\mathbf{K}\mathbf{u} = \mathbf{f}$ before solving.

---

Solving the System

Once $\mathbf{K}$ and $\mathbf{f}$ are assembled and boundary conditions applied, we solve:

$$\mathbf{u} = \mathbf{K}^{-1} \mathbf{f}$$

In practice, direct solvers (Cholesky, LU decomposition) or iterative solvers (Conjugate Gradient) are used because:

• $\mathbf{K}$ is large and sparse (most entries are zero)
• $\mathbf{K}$ is symmetric and positive definite (for stable structures)

Post-Processing

After solving for $\mathbf{u}$, we can compute:

• Strains: $\boldsymbol{\epsilon} = \mathbf{B} \mathbf{u}$
• Stresses: $\boldsymbol{\sigma} = \mathbf{D} \boldsymbol{\epsilon}$
• Reaction forces: $\mathbf{f}_{reactions} = \mathbf{K} \mathbf{u} - \mathbf{f}$

---

Summary

The beauty of FEA is that this same mathematical framework extends to heat transfer, fluid flow, electromagnetics, and beyond — always returning to the same fundamental pattern of discretization, assembly, and solution.