Numerical Methods for CAE
FEM, FVM, FDM Theory and Solver Selection

Every CAE software — Abaqus, Ansys, OpenFOAM, COMSOL — is ultimately a large numerical computation engine. Understanding the mathematical methods beneath the GUI doesn't just satisfy academic curiosity: it determines whether you interpret results correctly, recognize when a simulation is going wrong, and know how to fix it efficiently.

1. PDE Classification and Numerical Methods Overview

1.1 PDE Classification for a 2nd-Order PDE

For $A u_{xx} + 2B u_{xy} + C u_{yy} + \ldots = 0$, the discriminant $\Delta = B^2 - AC$ determines:

1.2 Methods Overview

2. Finite Difference Method (FDM)

2.1 Finite Difference Approximations

From Taylor series expansion, the first derivative:

Second derivative:

2.2 Stability Analysis: Von Neumann Method

For the 1D diffusion equation $u_t = \alpha u_{xx}$ with explicit FDM (FTCS):

Von Neumann stability analysis: assume $u_i^n = g^n e^{i\xi x_i}$. The amplification factor must satisfy $|g| \leq 1$:

This CFL-like stability condition is the fundamental constraint of explicit time integration: halving the grid spacing requires quartering the time step — a severe restriction for 3D fine meshes.

2.3 1D Advection Equation — Upwind Schemes

For $u_t + c\,u_x = 0$ (wave speed c > 0), central differences give an unstable scheme. The first-order upwind scheme is stable:

The CFL number ν must be ≤ 1 for stability — information cannot travel more than one cell per time step.

3. Finite Element Method (FEM) — Mathematical Foundations

3.1 The Variational (Ritz-Galerkin) Foundation

For the model problem $-\nabla\cdot(k\nabla u) = f$ with $u = 0$ on boundary:

Minimization form: The FEM solution minimizes the potential energy functional:

Galerkin (weak) form: Find $u^h \in V^h$ such that $\forall v^h \in V^h$:

3.2 Galerkin Orthogonality and Best Approximation

The Galerkin error $e = u - u^h$ satisfies:

This "Galerkin orthogonality" means the error is orthogonal to the approximation space in the energy inner product. Cea's lemma then gives the best approximation property:

3.3 Convergence Rates

For polynomial order p elements and smooth solutions, the h-convergence rate:

where h is the mesh size and s is the smoothness of the exact solution. For smooth problems with linear elements (p=1): energy error converges as h¹, L₂ error as h².

3.4 Element Types and Polynomial Orders

3.5 Locking Phenomena and Remedies

Low-order elements can exhibit artificial stiffness ("locking") in special cases:

4. Finite Volume Method (FVM)

4.1 FVM Discretization

Integrating a conservation law $\partial_t(\rho\phi) + \nabla\cdot(\mathbf{F}) = S$ over a control volume $V_P$:

Approximating as:

The challenge: computing the face flux $\mathbf{F}_f$ from cell-center values requires interpolation — the choice of interpolation scheme determines accuracy, stability, and conservation.

4.2 Gradient Reconstruction

Cell-center gradients are needed for the viscous fluxes and for second-order interpolation to faces. Least-squares gradient reconstruction:

where sum is over all face neighbors N, and $\mathbf{w}_{PN}$ are weights. This is more accurate than the Green-Gauss theorem approach on skewed meshes.

4.3 Flux Limiters for Convection

High-order (2nd+ order) convection schemes can produce oscillations (Gibbs phenomenon) near sharp gradients. TVD (Total Variation Diminishing) limiters enforce monotonicity:

Common limiters: Minmod ($\psi = \max(0, \min(1, r))$), Van Leer ($\psi = (r+|r|)/(1+|r|)$), Superbee. Each offers different tradeoffs between accuracy and monotonicity.

5. Linear System Solvers (Direct vs. Iterative)

5.1 Direct Solvers: Cholesky Factorization

For symmetric positive definite (SPD) systems (elastic FEM without rigid body motion): $\mathbf{K} = \mathbf{L}\mathbf{L}^T$.

Steps: (1) Analyze fill-in pattern → optimal reordering (AMD, nested dissection). (2) Symbolic factorization → allocate sparse storage. (3) Numerical factorization O(n_fill). (4) Forward/back substitution for each RHS.

5.2 Iterative Solvers: Krylov Methods

Krylov methods build an approximation from the Krylov subspace $\mathcal{K}_k = \text{span}\{r_0, \mathbf{A}r_0, \mathbf{A}^2r_0, \ldots, \mathbf{A}^{k-1}r_0\}$:

• CG (Conjugate Gradient): For SPD systems. Converges in at most n steps; with good preconditioning, much fewer. Computational cost: 1 matrix-vector product per iteration.
• GMRES: For general (non-symmetric) systems. Minimizes residual over Krylov subspace. More memory intensive (stores full Krylov basis), often restarted as GMRES(m).
• BiCGSTAB: For non-symmetric systems; fixed memory cost per iteration; less stable than GMRES but cheaper.

5.3 Preconditioning

Convergence rate depends on the condition number $\kappa(\mathbf{A}) = \lambda_{\max}/\lambda_{\min}$. Preconditioning replaces $\mathbf{A}\mathbf{x} = \mathbf{b}$ with:

Algebraic Multigrid (AMG) is the state-of-the-art preconditioner for FEM structural problems — it achieves nearly O(n) computational complexity. Ansys uses MEBA (Mechanical AMG), OpenFOAM uses GAMG for pressure equations.

6. Eigenvalue Solvers

6.1 The Lanczos Algorithm

For $(\mathbf{K} - \lambda\mathbf{M})\boldsymbol{\phi} = 0$, the shift-and-invert form avoids instability for clustered or zero eigenvalues:

After m steps, the original eigenvalues near shift σ are approximated by eigenvalues of the m×m tridiagonal matrix. Convergence of the extreme Ritz values is rapid (super-linear for well-separated eigenvalues).

6.2 AMLS and Component Mode Synthesis

For models with millions of DOFs where even Lanczos is too expensive, AMLS (Automated Multi-Level Substructuring) partitions the model into substructures, computes local eigenmodes, and assembles a reduced-order model. This can reduce a 10M DOF problem to a 50,000 DOF eigenvalue problem, with less than 1% error for the lowest modes. Abaqus, Nastran, and Ansys all offer this option.

7. Numerical Accuracy and Error Management

7.1 Richardson Extrapolation and GCI

Given results on three successively refined meshes with refinement ratio r:

Safety factor $F_s = 1.25$ (recommended by ASME V&V standards). GCI provides a reliable error bound for the discretization error.

7.2 Sources of Error in CAE

7.3 Condition Number and Ill-Conditioned Systems

The condition number of the stiffness matrix:

For structural FEM with linear elements: $\kappa \sim (L/h)^2$ — scales with the square of the mesh fineness. For double precision (64-bit), significant digits lost = $\log_{10}(\kappa)$. With κ = 10⁸ and 15 significant digits, you have 7 digits of accuracy — usually fine. With κ = 10¹² (very fine mesh or near-singular problem), only 3 digits remain — dangerous.

8. Solver Selection Guide

8.1 Structural Analysis Solver Selection

8.2 CFD Solver Selection

Related Interactive Tools

• Mesh Convergence (GCI) Calculator — Richardson extrapolation and grid convergence index
• Truncation Error Visualizer — Compare FD scheme accuracy orders
• 1D Heat Diffusion (FTCS) — Explore stability limit with CFL number
• Uncertainty Propagation — Monte Carlo error propagation for CAE inputs