Finite Element Analysis

Intuition

FEA replaces an infinite-dimensional field u(x) with a finite-dimensional approximation. Shape functions Nᵢ define compact support. Nodal unknowns uᵢ fill the ansatz. The solver finds uᵢ values that satisfy the PDE in a weak (integral-averaged) sense. Everything else is bookkeeping.[fn:: Best here means in a Galerkin sense. The residual is made orthogonal to the test space. For symmetric positive-definite operators, the test space coincides with the trial space. Galerkin is not the only choice. Least-squares, collocation, and Petrov-Galerkin methods exist. Your commercial code uses Galerkin by default unless you ask otherwise.]

FEA handles geometries with no closed-form solution. It scales to millions of elements. This leverage makes FEA dangerous. The method converges to the PDE you specify, not to reality. Wrong constitutive laws, missing contact, or missing nonlinearities produce a high-precision portrait of the wrong problem. A coarse mesh gives a low-precision portrait of the right problem. A stress result that varies by a factor of two between meshes is not a truth bound. It is a flag for insufficient discretization.[fn:: A common managerial failure mode is reporting the finer of two meshes as the answer because it looks more authoritative. Convergence is a property of the sequence, not of any single member. Richardson extrapolation exists precisely to discipline this behavior.]

FEA projects a hard problem onto an easy subspace. Always ask whether the projection discarded the feature you care about. Stresses at re-entrant corners often need local refinement or singularity removal. Displacements away from singularities usually converge well.

First Principles

The weak form and variational principle

Start from the strong form of linear elastostatics on a body omega. The divergence of sigma plus body force b equals zero. Sigma equals D times epsilon. Dirichlet conditions set u equal to u bar on gamma u. Neumann (traction) conditions set sigma n equal to t bar on gamma t. Multiply by a test function v that vanishes on gamma u. Integrate by parts. You obtain the weak form:

\[ \int_\Omega \boldsymbol{\varepsilon}(\mathbf{v})\!:\!\mathbf{D}\!:\!\boldsymbol{\varepsilon}(\mathbf{u})\,dV = \int_\Omega \mathbf{v}\cdot\mathbf{b}\,dV + \int_{\Gamma_t} \mathbf{v}\cdot\bar{\mathbf{t}}\,dA. \]

This is also the stationarity condition of the total potential energy Pi equals U minus W external. This is the principle of minimum potential energy for the conservative case. The equivalence between weak form and stationary potential gives FEA its dual interpretation as either a variational method or a weighted-residual method. Pick whichever framing makes the proof you need easier.[fn:: For non-conservative problems such as plasticity, contact with friction, or large-deformation path-dependence, the potential-energy interpretation breaks down because energy is no longer a state function. The weak form still holds. Textbooks lead with the weak form and mention the variational view as a special case. Treat any commercial-code strain-energy value reported for a path-dependent problem as a bookkeeping number, not a thermodynamic state function.]

Discretization: element stiffness and assembly

Restrict u and v to the span of shape functions. After partitioning omega into elements omega e and substituting the ansatz, the element stiffness matrix emerges as:

\[ \mathbf{K}_e = \int_{\Omega_e} \mathbf{B}^T \mathbf{D}\,\mathbf{B}\,dV, \]

where B is the strain-displacement matrix. B contains derivatives of shape functions mapped from parent to physical coordinates via the Jacobian. D is the constitutive matrix. The global system K u equals f assembles by scattering the Ke contributions to the appropriate global degrees of freedom according to the element connectivity. Assembly is sparse, symmetric, and embarrassingly parallel. The engineering cost is dominated by the linear solve, not by forming K.

Element types

Tetrahedra mesh automatically. Hexahedra give better accuracy per degree of freedom but resist automatic meshing. Element stiffness is Kₑ = ∫ BᵀDB dV. B maps shape-function derivatives to physical coordinates via the Jacobian. D is the constitutive matrix. The global system Ku=f assembles by scattering Kₑ contributions to global degrees of freedom. Assembly is sparse, symmetric, and parallelizable. The linear solve dominates the cost, not K formation.

Hexahedra mesh manually but deliver accuracy per degree of freedom. Tetrahedra mesh automatically. Hex offers five to ten times better energy-norm accuracy per degree of freedom than tet at equal count. In practice, a finished tet mesh beats an unfinished hex mesh.[fn:: This is an instance of a general pattern in numerical methods: the theoretically-superior discretization often loses to the practically-available one. Compare the persistence of low-order finite volumes in CFD against spectral element methods.]

Reduced integration and the hourglass disease

Using fewer Gauss points than the element polynomial degree requires for exact integration saves cost and relieves locking in nearly-incompressible and thin-bending regimes. The spurious stiff modes that full integration over-counts are exactly the ones that pollute the under-integrated quadrature. The price is hourglassing. These are zero-energy deformation modes. The element can wrinkle without storing strain energy. These modes propagate through the mesh as a checkerboard pattern. Modern codes control hourglassing with hourglass stiffness, assumed-strain formulations, or enhanced-strain methods. None of these solutions is free. Each injects a length scale or a user-settable parameter. You should know these parameters exist before trusting the output.

Shape functions, meshing, refinement

• h-refinement: shrink elements and keep polynomial order fixed. This approach is simple, robust, and automatic for tetrahedra. It is what most engineers mean by refine the mesh.
• p-refinement: keep the mesh and raise the polynomial order. This gives spectral convergence for smooth solutions. It is excellent for stress intensity where the singularity is algebraic. p-methods capture it with geometric convergence.
• hp-refinement: adaptively choose both element size and polynomial order. This is the theoretically optimal strategy according to Babuska and Suri. It is the practical choice when your code supports it. Most codes do not support it well.

Meshing is where the art lives. The rule of thumb is to refine where the gradient of the solution is large, not where the geometry is complicated. These often coincide. Stress concentrations, contact zones, crack tips, and boundary layers all demand local refinement. A globally fine mesh is a confession that you did not think about where the action is.

Boundary conditions: Dirichlet vs Neumann

Dirichlet (prescribed displacement) conditions are enforced by constraint. You can modify rows and columns of K directly. Modern codes use Lagrange multipliers or penalty methods instead. Neumann (prescribed traction) conditions enter naturally as contributions to the load vector f on the right-hand side. The asymmetry matters. Neumann is cheap and stable. Dirichlet is where constraint-handling bugs live. A common silent error is over-constraining a face that should slide. This error builds in a spurious reaction couple that dominates the very stress you are computing.[fn:: In 3D, a fully-clamped face is rarely what you actually want physically. Most real supports have some compliance. Modeling every support as a rigid clamp systematically biases stress upward and natural frequencies downward. If the answer is sensitive to the support compliance, that sensitivity is itself an important result.]

Convergence: Richardson, and the converged does not mean correct caveat

The proper convergence study is not run it once more with a finer mesh and see if it changed. The proper study produces a sequence uh, uh over 2, uh over 4 with a geometric refinement ratio of r equals 2. You assume uh equals u exact plus C hp. You then solve for p and the extrapolated u exact using Richardson extrapolation. This gives you an observed order p that you can sanity-check against the known theoretical order of your element. If the observed p is much lower than the theoretical one, you are either in a singular region, in a pre-asymptotic regime where the mesh is still too coarse, or you have a bug. If you cannot afford uh over 4, you cannot do Richardson extrapolation. You can only claim monotonicity, which is much weaker.

Validation against analytical limits

Before trusting FEA on the real geometry, run it on a problem you can solve by hand. Try a cantilever beam, a pressurized thick cylinder, a plate with a hole, or a contact Hertz problem. Agreement to less than one percent in a quantity you can verify analytically is the minimum bar before believing results on parts you cannot verify. This is not paranoia. It is the same epistemic posture as unit testing. A model that fails the analytical check has a setup bug. A model that passes it may still have a different setup bug on the real geometry. At least the platform is sane.

Locking: shear and volumetric

• Shear (transverse) locking: low-order elements resist bending by spurious shear strain when they are slender. Linear shape functions cannot represent the quadratic bending displacement profile. The element becomes orders of magnitude too stiff. Reduce integration, use assumed-strain methods such as MITC, or use quadratic elements to cure this.
• Volumetric locking: in near-incompressible regimes where Poisson nu approaches 0.5 or in plastic flow at constant volume, the pressure field cannot be represented by the displacement interpolation and the element locks. Use mixed u-p formulations such as B-bar, F-bar, or Q1P0. Selective reduced integration on the volumetric part or enriched elements also cure this.

Both are discretization artifacts, not material behavior. If your nearly-incompressible rubber model is stiffer than reality, suspect locking before suspecting the material data.

Worked Example — Cantilever Convergence Study

Take a steel cantilever. L equals 100 mm. Rectangular section b times h equals 10 by 10 mm. E equals 200 GPa. Nu equals 0.3. End load P equals 100 N. The Euler-Bernoulli tip deflection is:

\[ \delta = \frac{P L^3}{3 E I} = \frac{100 \cdot 100^3}{3 \cdot 200000 \cdot (10\cdot10^3/12)} = 0.200\ \text{mm}. \]

Mesh the beam with plane-stress Q4 (bilinear quad) elements. Use length L/n per element. Use full 2-by-2 integration. Clamp the left end. Distribute the end load over the rightmost nodes. Track dh at the tip:

Richardson extrapolation on n equals 16 and 32 gives delta infinity approximately 0.2001 mm. This agrees with the Euler-Bernoulli result to 0.05 percent. Three observations:

1. The coarse mesh (4 elements along the span) is off by 13 percent. That mesh would be entirely defensible in a design review that did not ask did you converge. This is the leverage of asking the cheap question. 2. The observed order p approximately equals 2 matches the theoretical rate for Q4 in the bending regime with full integration. If you saw p approximately equals 1, suspect locking or a singular load. 3. Repeating with a single concentrated node load at the tip gives a stress field that diverges as you refine. The analytical load is a Dirac delta. FEA produces a singular stress. Distribute the load over the end face or accept that the peak stress at the load point is mesh-dependent forever.

Units sanity check. E in MPa (200000), forces in N, lengths in mm gives stresses in MPa and deflections in mm. Mixing units is the number one silent bug in FEA scripts. Pick one consistent unit system. The mm-N-MPa system is the standard for mechanical parts. Never deviate.

Pitfalls

These are the mistakes that cost a week. They are not subtle. They are silent. The solver does not error. The result is just wrong.

1. Point load on a single node. A nodal force is a Dirac load. The local stress is formally unbounded and diverges as h to the minus one power under refinement. Cures include distributing the force over the nodes of the loaded face via consistent nodal forces (work-equivalent), or modeling the actual contact patch. Never report peak stress at the load point. It is a discretization artifact.

2. Reduced-integration hourglassing. A mesh that looks fine but shows a checkerboard of zero-energy modes in the displacement or stress plot. Always check the deformed shape visually. If it has spurious wiggle where the physics should be smooth, your hourglass control is inadequate or off.

3. Not checking reaction force sums. Sum the reaction forces at the supports. They must equal the applied loads to solver tolerance, typically about 10 to the minus 6th power. If they don't, your constraints are wrong, your loads are applied to the wrong place, or the solve did not actually converge. This 10-second check catches a large fraction of setup errors. Do it every time.

4. Coarse mesh at stress concentrations. A fillet root, a notch, a hole edge: these are where the answer matters and where the gradient is steepest. A uniformly coarse mesh will systematically under-predict the peak stress. A designer who reads von Mises equals 0.8 sigma y from a coarse mesh is being told a comfortable lie. The truth is at best the coarse-mesh lower bound on stress is 0.8 sigma y.

5. Treating a converged answer as a correct one. Convergence means the discretization error is small. It says nothing about whether the PDE you solved is the right one, whether the boundary conditions are right, whether the material model is calibrated, or whether the loading envelope is realistic. A converged-but-wrong answer is more dangerous than an unconverged one because it carries false confidence. The discipline of FEA is largely the discipline of separating these error sources and not letting convergence in one disguise errors in another.

6. Mesh dependency of plasticity and other local results. Plasticity localizes. A softening material can collapse into a band one element wide. The dissipated energy and the peak strain both scale with h. Without a length scale in the constitutive law such as nonlocal, gradient, or cohesive-zone models, the answer is a function of the mesh you chose. This is not a minor caveat. It is a foundational limitation of classical local plasticity in FEA. If your problem involves fracture or localization and your code does not have a regularization, treat all peak-strain and energy numbers as mesh-dependent no matter how fine you go.[fn:: The same applies to damage mechanics, to Eulerian advection of discontinuities, and to many mesh-adaptive contact algorithms. The rule is that any quantity that scales with h as h goes to zero is not a material property. It is a discretization parameter wearing a material-property costume.]

References

• K.-J. Bathe, /Finite Element Procedures/, 2nd ed., 2014. — the reference for nonlinear and the source of much of the right-way thinking.
• T. J. R. Hughes, /The Finite Element Method: Linear Static and Dynamic Finite Element Analysis/, Dover, 2000. — the rigorous mathematical FEM text. If you want to know why Galerkin works, start here.
• O. C. Zienkiewicz, R. L. Taylor, J. Z. Zhu, /The Finite Element Method: Its Basis and Fundamentals/, 7th ed., 2013. — the encyclopedic reference. Consult by topic.
• R. D. Cook, D. S. Malkus, M. E. Plesha, R. J. Witt, /Concepts and Applications of Finite Element Analysis/, 4th ed., 2002. — the practical undergraduate-to-grad text. Best first book for an engineer.
• I. Babuska, B. Szabo, "On the rates of convergent approximations in the finite element method", /Int. J. Numer. Methods Eng./, 1987. — the hp-convergence theory.
• T. Belytschko, W. K. Liu, B. Moran, K. Elkhodary, /Nonlinear Finite Elements for Continua and Structures/, 2nd ed., 2014. — the nonlinear bible. Essential for plasticity, contact, large deformation.
• G. Strang, G. Fix, /An Analysis of the Finite Element Method/, 1973 (2nd ed. 2008). — the cleanest mathematical introduction, still.

Related

• Mechanics: Statics — the underlying continuum mechanics. FEA is the numerical arm of this.
• Mechanics: Dynamics (Modal FEA) — eigenvalue extraction, modal analysis, the dynamics cousin of static FEA.
• Failure Criteria — what to do with the stress field FEA produces. von Mises, Tresca, Mohr-Coulomb.
• Fatigue & Fracture — where FEA peak stresses feed S-N and delta-K life calculations. Beware the singularity caveat.
• OpenFOAM (CFD cousin) — finite volume method for fluids. The same weak/strong-form dichotomy and mesh-convergence discipline applies.
• OpenFOAM Meshing — boundary-layer meshing, y+ requirements. The CFD analog of refine where the gradient is large.