At its mathematical foundation, FEA approximates the behavior of a continuous physical body by dividing it into a finite number of discrete spatial subdomains called elements. These elements are interconnected at shared geometric coordinates designated as nodes.
By defining polynomial shape functions that interpolate displacement across each element, the software translates the differential equations of elasticity into a global system of algebraic equations:
[K]{u} = {F}
Where [K] represents the global structural stiffness matrix, {u} is the unknown nodal displacement vector, and {F} represents the externally applied force vector. Once displacements are calculated, strains are derived via spatial differentiation, and stresses are computed through the material’s constitutive stiffness tensor (σ = [D]ε).
To evaluate FEA software objectively, an engineering team must understand the three mechanical domains that govern solver architecture.
1. Element Formulations and Numerical Pathology
The accuracy of an FEA calculation depends heavily on the element formulation chosen to discretize the continuum:
- First-Order (Linear) vs. Second-Order (Quadratic) Elements: Linear tetrahedral elements (4-node tets) are computationally inexpensive, but they exhibit artificial numerical stiffness (“shear locking”) when subjected to bending loads. Using linear tets in critical stress regions invariably underestimates true deformation and underreports peak stresses. High-fidelity structural analysis demands either second-order elements (10-node tets or 20-node hexes) or specialized formulations with interior integration schemes.
- Hexahedral (“Brick”) vs. Tetrahedral Meshing: Hexahedral elements provide superior convergence, minimal numerical diffusion, and clean boundary-layer stress resolution with significantly fewer total degrees of freedom than tetrahedral grids. However, generating an all-hex mesh on complex, organic cast or forged geometries requires substantial manual geometry partitioning.
- Volumetric Locking and Incompressibility: When analyzing materials with Poisson’s ratios approaching 0.5 (such as natural rubbers, elastomers, or metals undergoing massive plastic flow), standard element formulations lock up, falsely reporting near-zero deformation under heavy pressure. Advanced solvers incorporate selective reduced integration or mixed u-P (displacement-pressure) hybrid element formulations to prevent volumetric locking.
2. Time Integration: Implicit Static vs. Explicit Dynamic Solvers
The mathematical strategy used to advance a solution across load increments or time steps fundamentally divides the FEA landscape:
- Implicit Solvers (Newton-Raphson Equilibrium): Standard implicit engines solve for static or quasi-static equilibrium at every discrete load step by computing the tangent stiffness matrix and inverting it iteratively until the residual force vector falls below a convergence tolerance. Implicit integration is unconditionally stable across large time increments, making it the industry standard for linear elastic frameworks, steady-state thermal analysis, structural fatigue, and low-frequency vibration. Its Achilles’ heel is convergence: if complex contact friction slips, localized material cracking initiates, or geometric buckling occurs, the tangent stiffness matrix becomes ill-conditioned, the Newton-Raphson iterations diverge, and the solver aborts.
- Explicit Solvers (Wave Propagation): Explicit dynamic integration bypasses global matrix inversion entirely. Instead, it computes acceleration directly from the net nodal forces and advances velocities and displacements forward using central difference integration across minute time steps bounded by the Courant-Friedrichs-Lewy (CFL) condition:
Δt ≤ L_min / c_d
Here, Δt is the maximum stable time increment, L_min represents the smallest characteristic element length in the mesh, and c_d is the dilatational wave speed through the material. Because explicit solvers do not require global convergence loops, they handle severe non-linear events without aborting: automotive crash testing, consumer electronics drop simulations, explosive ballistic impact, and sheet-metal stamping.
3. The Three Non-Linearities: Geometric, Material, and Boundary
Linear analysis assumes three conditions: displacements are infinitesimal, material stress-strain relationships follow a straight line (E = constant), and boundary constraints never change during loading.
Real engineering structures routinely violate these assumptions:
- Geometric Non-Linearity (Large Deflection / Rotation): When a slender column deflects, its internal moment arms change; when a thin diaphragm inflates, membrane stiffening alters load resistance. Solvers must continuously update element coordinate systems and formulate strain tensors (such as Green-Lagrange or logarithmic strains) to capture structural reality.
- Material Non-Linearity (Plasticity, Viscoelasticity, Hyperelasticity): Metals deform past yield into strain-hardening regimes; polymers creep under prolonged mechanical stress; elastomeric seals exhibit non-linear S-shaped hyperelastic stress-strain curves. Solvers must support advanced constitutive plasticity algorithms (von Mises, Johnson-Cook, Drucker-Prager) and strain-energy density functions (Mooney-Rivlin, Ogden).
- Boundary Non-Linearity (Contact Mechanics): Contact is inherently non-linear because contact status (open or closed) changes abruptly. Calculating contact pressure across sliding, deforming interfaces requires advanced numerical formulations—such as the penalty method or augmented Lagrange multipliers—to prevent unphysical penetration without destabilizing the global solver matrix.
