FEM / FEA — Finite Element Method — Engineering Reference

1. What it is

Finite Element Method (FEM) is the numerical solution of partial differential equations (PDEs) and variational problems by discretizing a continuous domain Ω into a finite number of geometrically simple sub-regions (elements) joined at nodes, then approximating the unknown field — displacement, temperature, electric potential, magnetic vector potential, fluid velocity — by piecewise polynomial shape functions N_i(x) that are nonzero only on neighbouring elements. Finite Element Analysis (FEA) is the engineering application of that method: meshing real geometry, applying boundary conditions and loads, solving the assembled linear or nonlinear system, and post-processing the results.

FEM is the dominant numerical method for solid mechanics (stress, strain, vibration, plasticity, contact, fracture), heat transfer (conduction, convection BCs, radiation, transient thermal), electromagnetics (low-frequency magnetostatics and eddy currents, high-frequency wave propagation, antennas), acoustics (vibro-acoustics, mufflers, sound radiation), and most coupled multi-physics (thermo-mechanical, piezoelectric, magneto-mechanical, fluid-structure interaction). For pure fluid flow the finite volume method (FVM) dominates instead — ANSYS Fluent, Star-CCM+, OpenFOAM are FVM — though FEM-based CFD exists (COMSOL, FEniCSx).

The method became industry-standard during the 1970s in aerospace and nuclear, broadened to automotive and civil through the 1980s, and became commodity desktop software around 1995 once Windows NT made workstation-class FEA affordable. Today every modern engineering design with non-trivial geometry — engine block, aircraft wing, MRI gradient coil, dam, prosthetic hip, MEMS sensor, smartphone drop case — is FE-verified before it is manufactured.

Place in the design stack: CAD → CAE pre-processing (meshing + BCs) → FE solve → post-processing → optimization → test correlation → certification.

2. Why it matters

Closed-form solutions exist for a handful of canonical geometries — a uniform bar in tension, an Euler–Bernoulli beam under standard loads, a Lamé thick-walled cylinder, a Kirsch infinite plate with a hole. The rest of engineering reality — engine blocks with cooling galleries, additively-manufactured lattices, composite layups with ply drops, weldments with HAZ gradients, contact between gear teeth, snap-fit hooks, automotive crash structures, scoliosis braces, tibial trays — cannot be solved analytically. Before FEM, the design process was rules-of-thumb plus prototype-and-test: slow, expensive, conservative, and routinely wrong.

With FEM, the loop is closed inside CAD. A solid model goes through a meshing kernel (Ansys Workbench, HyperMesh, ANSA, Cubit) into a solver and out to a visualizer in hours. Engineers iterate dozens of geometric or material variants before any tooling is cut. Topology-optimization workflows (Altair OptiStruct, ANSYS Mechanical Optimization, nTopology) drive the geometry itself from FE-computed sensitivities, producing organic load-path-aligned designs that human intuition would not have proposed.

Codes formalize the dependency. ASME BPVC Section VIII Div. 2 — Design by Analysis explicitly mandates FE-based protection against plastic collapse, local failure, buckling, ratcheting, and fatigue when the simpler Division 1 rules are exceeded. Eurocode 3 Part 1-5 allows FE-based plate-girder stability checks under EN 1993-1-5 Annex C with prescribed imperfection patterns. Aerospace AC 25.571 requires damage-tolerance demonstrations whose load distributions come from FE global/local models. ISO 10303-21 STEP AP242 carries the FE neutral-format exchange between CAD/CAE systems.

The credibility of an FE result is now itself an engineering deliverable. ASME V&V 10-2006 (reaffirmed 2016) lays out the verification-and-validation framework for computational solid mechanics; V&V 20-2009 is the CFD counterpart; V&V 40-2018 extends V&V to medical-device computational modelling for FDA submissions. NAFEMS publishes the canonical benchmark library (linear elastic, nonlinear material, dynamics, heat transfer, contact, vibroacoustic) that every commercial code is tested against and that competent engineers re-run after a major code upgrade.

3. First principles

3.1 Strong form, weak form, Galerkin

A typical solid-mechanics boundary value problem in strong form:

div σ(u) + b = 0          in Ω           (equilibrium)
σ = D : ε(u)               in Ω           (constitutive — linear isotropic D here)
ε(u) = ½ (∇u + ∇u^T)       in Ω           (small-strain kinematics)
u = û                       on Γ_D          (essential / Dirichlet BC)
σ · n = t̂                  on Γ_N          (natural / Neumann BC, traction)

The weighted-residual / weak form multiplies the equilibrium equation by a test function w ∈ V_0 (vanishing on Γ_D), integrates over Ω, and applies the divergence theorem to shift one derivative onto w:

∫_Ω ε(w) : σ(u) dΩ  =  ∫_Ω w · b dΩ  +  ∫_{Γ_N} w · t̂ dΓ          ∀ w ∈ V_0

Choosing the test space equal to the trial space (Galerkin’s method, Boris Galerkin 1915) and discretizing with the same shape functions for both — u_h(x) = Σ N_i(x) u_i and w_h(x) = Σ N_i(x) w_i — produces, after the arbitrary w_i cancel, the assembled linear system:

K · u = F

with global stiffness K and global force vector F.

3.2 Element kinematics and stiffness

For a generic solid element with nodal displacements u_e:

u(x)   = N(x) · u_e
ε(x)   = B(x) · u_e          (B = symmetric gradient of N — strain-displacement matrix)
σ(x)   = D · B(x) · u_e       (D = constitutive matrix; 6×6 in 3-D Voigt notation)

The element stiffness matrix comes from the inner integral of the weak form:

K_e = ∫_{Ω_e} B^T · D · B  dΩ

evaluated by Gauss–Legendre quadrature at a small number of integration points. For a linear 8-node hex with full 2×2×2 = 8 integration points, K_e is 24×24 (3 DOF × 8 nodes); for a quadratic 20-node hex with 3×3×3 = 27 points, K_e is 60×60.

The element force vector combines body force, traction, and any initial strain (thermal, inelastic):

F_e = ∫_{Ω_e} N^T · b dΩ + ∫_{Γ_e ∩ Γ_N} N^T · t̂ dΓ + ∫_{Ω_e} B^T · D · ε_0 dΩ

3.3 Assembly

A DOF-mapping table (LM(e, i) in textbook notation) records, for each element e and its local DOF i, the global DOF index. Assembly is a scatter operation:

for each element e:
  K_e = element_stiffness(e)
  F_e = element_force(e)
  for each local DOF i:
    g_i = LM(e, i)
    F[g_i] += F_e[i]
    for each local DOF j:
      g_j = LM(e, j)
      K[g_i, g_j] += K_e[i, j]

K is sparse, symmetric (positive definite once supports are applied), and large — a modern automotive crash model has 5–20 million DOFs; an aircraft global finite-element model (GFEM) has 1–5 million; a single-component machine-design FEA is typically 50k–500k.

3.4 Boundary conditions

Type	Mathematical form	Physical meaning
Essential (Dirichlet)	u_i = û_i	Prescribed displacement, including u = 0 supports
Natural (Neumann)	σ · n = t̂	Applied traction, point load, pressure
Mixed / Robin	α u + β σ·n = γ	Elastic foundation, convective heat transfer h(T − T∞)
Constraint	Σ c_i u_i = d	MPC, RBE2/RBE3 rigid links, periodic, contact

Essential BCs are enforced by direct elimination (delete the row/column), penalty (add a large diagonal stiffness), or Lagrange multipliers (extend the system with one equation per constraint). The Lagrange method is exact but expands the system; penalty is approximate but cheap.

3.5 Convergence

A finite element method converges to the PDE solution as the mesh is refined if and only if the elements satisfy three properties:

1. Consistency — the element can exactly reproduce a constant-strain (constant-gradient) state. The patch test (Irons 1965; Taylor, Beresford, Wilson 1976) is the operational check: assemble a small patch of elements with prescribed boundary displacements consistent with a uniform strain field; if every interior node satisfies equilibrium and reproduces the strain exactly to machine precision, the element passes.
2. Stability — the element is rank-sufficient (no spurious zero-energy modes), and the discrete inf–sup (LBB) condition holds for mixed formulations.
3. Completeness — the shape function basis can represent rigid-body motion and constant strain.

Given these, the asymptotic error in the energy norm is O(h^p) where h is the element size and p is the complete polynomial order (linear: p = 1; quadratic: p = 2). Stress error converges one order slower than displacement error.

3.6 Alternative variational principles

Pure displacement (Galerkin) is the default for solids, but mixed and hybrid formulations dominate in two regimes:

• Hellinger–Reissner — treats displacement and stress as independent fields. Useful for nearly incompressible materials and for stress-recovery.
• Hu–Washizu — treats displacement, strain, and stress as three independent fields. Underlies many high-performance shell elements (MITC family, Bathe; assumed-strain elements, Belytschko).
• Mixed u-p — displacement and pressure as independent fields. The standard fix for volumetric locking in nearly incompressible elasticity and J2 plasticity.

4. Element types

4.1 Element family table

Family	Dim	Nodes	DOF/node	Interpolation	Typical use	Pitfalls
Truss / bar	1-D	2	1–3 (u)	Linear axial	Trusses, cables, struts	No bending, no shear
Beam (Euler–Bernoulli)	1-D	2	6	Hermite cubic	Slender frames, L/h > 10	Ignores transverse shear
Beam (Timoshenko)	1-D	2–3	6	Linear or quadratic	Deep beams, sandwich	Shear-locking with linear
Pipe	1-D	2	6	Beam + section pressure	Piping stress, ASME B31	Ovalization not captured
Plane stress (CST)	2-D	3	2	Linear	Sheet panels, in-plane only	Constant strain → over-stiff
Plane stress (LST / T6)	2-D	6	2	Quadratic	Better in-plane	Cost vs Q8
Plane stress (Q4)	2-D	4	2	Bilinear	In-plane structures	Shear locking under bending
Plane stress (Q8 serendipity)	2-D	8	2	Serendipity quadratic	High-accuracy 2-D	No interior node — less stable
Plane stress (Q9 Lagrange)	2-D	9	2	Biquadratic Lagrange	High-accuracy	Slight extra DOF cost
Plane strain	2-D	as above	2	as above	Long prismatic bodies, tunnels	Wrong if ends are free
Axisymmetric	2-D	3/4/6/8	2 (u_r, u_z)	as above	Pressure vessels, shafts	Axis-symmetry must be true
Plate (Kirchhoff, DKT/DKQ)	2-D	3/4	3	Discrete-Kirchhoff	Thin plates, t/L < 0.05	No transverse shear
Shell (Mindlin–Reissner)	2-D	4/9	5–6	Bilinear/biquadratic	Thin-walled vessels, fuselage	Shear/membrane locking
Shell (MITC4 / MITC9)	2-D	4/9	5–6	Mixed interpolation of tensorial components	Bathe shells; default in ADINA	Robust against locking
Solid tet (TET4)	3-D	4	3	Linear	Auto-mesh of complex CAD	Over-stiff; locking in plasticity
Solid tet (TET10)	3-D	10	3	Quadratic	Default for complex CAD	More cost; needs node mid-side OK
Solid hex (HEX8)	3-D	8	3	Trilinear	Structured solid mesh	Hourglassing with reduced int.
Solid hex (HEX20 serendipity)	3-D	20	3	Quadratic serendipity	High-accuracy solids	Hard to mesh complex CAD
Solid hex (HEX27 Lagrange)	3-D	27	3	Triquadratic Lagrange	Research / spectral	Rare in commercial
Wedge / Prism	3-D	6 / 15	3	Mixed	Transition between hex and tet	Distortion sensitivity
Pyramid	3-D	5 / 13	3	Special	Hex-to-tet transition	Limited solver support

4.2 Integration: full vs reduced

For a linear 8-node hex element, full integration is 2×2×2 = 8 Gauss points. Reduced integration is 1 point at the centre. Reduced integration:

• Eliminates shear locking in thin-walled regions (the element does not over-resist bending).
• Eliminates volumetric locking in nearly-incompressible regimes (ν → 0.5, or J2 plasticity).
• Reduces cost ~8× per element.
• Introduces hourglass modes — zero-energy deformation patterns that the single integration point cannot detect.

ANSYS LS-DYNA, Abaqus/Explicit, and RADIOSS use single-point reduced integration with hourglass control (Flanagan–Belytschko stiffness or viscous; perturbed-Jacobian) as the default for explicit-dynamic analyses with hex elements. Abaqus/Standard offers B-bar and F-bar formulations to suppress volumetric locking while keeping full integration. ANSYS Mechanical’s SOLID185 has both full-integration and enhanced-strain (B-bar-like) modes selectable via KEYOPT.

4.3 Higher-order: p-FEM and spectral

The h-version refines the mesh; the p-version raises the polynomial order on a fixed coarse mesh. StressCheck (ESRD) is the canonical commercial p-FEM solver, used heavily in aerospace local-detail stress analysis. Spectral elements (Patera 1984) use Lagrange interpolants on Gauss–Lobatto nodes; they converge exponentially for smooth fields and underpin high-order CFD codes like Nek5000 and seismic-wave codes like SPECFEM3D.

5. Common pitfalls (element / mesh)

1. Shear locking — Linear plate / shell / beam elements with full integration over-resist bending because the displacement field cannot satisfy the Kirchhoff constraint (γ = 0) without zero curvature. Symptom: deflection 10×–100× too small as element thickness decreases. Fix: reduced integration, MITC (mixed interpolation of tensorial components), assumed natural strain (ANS), or enhanced assumed strain (EAS). In Abaqus, use S4R (4-node reduced) or S8R (8-node reduced) shells; in ANSYS, SHELL181 with KEYOPT(3) = 0 reduced or 2 full+enhanced strain.

2. Volumetric locking — Standard displacement elements cannot represent isochoric (volume-preserving) deformation when ν → 0.5 (rubber, incompressible plasticity) without spurious pressure oscillation. Symptom: stress field looks checkerboard; reaction force off by 30 %+. Fix: mixed u-p elements (Abaqus C3D8H, ANSYS hybrid), B-bar / F-bar, selective reduced integration, or assumed-strain formulations.

3. Aspect ratio — An element with edge ratio > 5:1 distorts the constitutive map; > 10:1 is usually a meshing failure. Target < 3:1 in stressed regions, < 10:1 in low-gradient regions. NASTRAN, Abaqus, ANSYS all flag elements above ~10 as warnings.

4. Element distortion — Internal angle near 0° or 180°, or a negative Jacobian determinant (J_det < 0), means the element maps inside-out from parametric to physical space. The element matrix is singular or wrong; meshers reject these or flag them. Targets: hex internal angles 45°–135°, tet dihedral angles 20°–160°.

5. Hourglassing — In reduced-integration HEX8 (or QUAD4 in 2-D) there are 12 zero-energy deformation modes that a single integration point cannot detect. The mesh deforms in a tell-tale “hourglass” or “zig-zag” pattern. Fix: enable hourglass control (Flanagan–Belytschko stiffness type 4 or 5 in LS-DYNA, ENHANCED in Abaqus); refine the mesh; or use full-integration elements (C3D8 instead of C3D8R in Abaqus). Always check the hourglass-energy / internal-energy ratio — it should stay below 5 % (10 % maximum) for a credible explicit run.

6. Mesh dependence for softening / damage — Strain-softening constitutive models (concrete cracking, ductile damage with CDM, fibre failure) lose ellipticity at peak load; the FE solution then depends pathologically on the mesh size. As h → 0, the dissipation goes to zero — physically wrong (Bažant’s localization paradox). Fix: fracture-energy regularization (smear crack energy over a characteristic length), nonlocal damage, gradient-enhanced damage, phase-field fracture, or viscous regularization (rate-dependent softening with a fictitious time scale).

7. Stress concentrations — Use h-refinement around holes, fillets, and re-entrant corners — but only down to the physical radius. A “sharp” corner in CAD that is really a 0.1 mm fillet in manufacturing should be modelled at 0.1 mm radius, not as a mathematical singularity. Submodelling (cut-boundary displacement) and p-enrichment are standard.

8. Singular points — At a true mathematical re-entrant corner, a point load, or a bonded bi-material crack tip, the analytical stress is infinite (σ ~ r^(−λ) with 0 < λ < 1). Mesh refinement increases the reported peak without bound. Use J-integral, K_I extraction, or stress-averaged-over-a-region rather than peak nodal stress.

9. Anisotropic material orientation — Composites, rolled plate, single crystals, and wood have material axes that must be assigned per-element. A common error is leaving the material coordinate system at global default in a curved shell — the laminate then has wrong fibre direction everywhere. Always plot the material orientation triad on the deformed shape and verify.

10. Connection idealization — Welded, bolted, riveted, and bonded connections are often modelled as MPCs, RBE2 (rigid), RBE3 (load-distribution), CBUSH (lumped spring), or contact pairs. Each idealization changes local stiffness and stress. RBE2 over-constrains; RBE3 under-constrains. Stress in the connector region is almost always wrong by 20–50 %; use it for global-load distribution, not local-stress allowables.

6. Worked examples

6.1 Example A — Cantilever beam validation

Problem. A steel cantilever beam, length L = 1.0 m, width b = 50 mm, depth h = 10 mm, fixed at one end, tip load P = 1 kN at the free end. E = 200 GPa, ν = 0.30. Compute tip deflection by (i) analytical beam theory, (ii) 10 Euler–Bernoulli beam elements, (iii) 10 × 1 linear 4-node plane-stress quads, (iv) 10 × 1 quadratic 8-node plane-stress quads.

Analytical.

I = b · h³ / 12 = 0.050 · (0.010)³ / 12 = 4.167 × 10⁻⁹ m⁴
δ_tip = P · L³ / (3 · E · I) = 1000 · 1.0³ / (3 · 200×10⁹ · 4.167×10⁻⁹)
      = 1000 / 2500 = 4.00 × 10⁻³ m = 4.00 mm

FE results.

Model	Element	DOFs	δ_tip (mm)	Error
Analytical	—	—	4.000	—
10× Euler–Bernoulli beam (B23 / BEAM188)	1-D Hermite cubic	66	3.996	−0.10 %
10×1 Q4 plane stress (CPS4 / PLANE182)	2-D bilinear, full int.	44	3.60	−10.0 %
10×1 Q8 plane stress (CPS8 / PLANE183)	2-D serendipity quad	138	3.99	−0.25 %
40×4 Q4 refined	bilinear, full int.	410	3.97	−0.75 %
40×4 Q4 with EAS / B-bar	bilinear, enhanced strain	410	4.00	−0.05 %

Lessons. (1) Pick the right element family for the problem — a 1-D beam element captures pure-bending exactly in its native form. (2) Coarse linear quads exhibit shear locking; quadratic shape functions or enhanced-strain formulations fix it without refinement. (3) Always run a convergence study before trusting a peak quantity.

6.2 Example B — Stress concentration at a hole (Kirsch validation)

Problem. A 200 mm wide × 5 mm thick steel plate with a central 20 mm diameter through-hole, under far-field uniaxial tension σ_∞ = 100 MPa. Compute σ_max at the hole edge and compare with the Kirsch (1898) analytical result K_t = 3.00.

Analytical (Kirsch, infinite plate): σ_max = 3 · σ_∞ = 300 MPa. The plate is wide enough relative to the hole (D_hole / W = 20/200 = 0.10) that the finite-width correction (~3 %) keeps K_t close to 3.0 — full Howland correction K_t = 3.03 for D/W = 0.1.

FE setup. Plane stress, exploit 1/4 symmetry, model 100 × 100 mm quadrant with hole. Compare three meshes.

Mesh	Elements around 1/4 hole	Total elements	DOFs	σ_max (MPa)	K_t = σ_max / σ_∞	Error vs Kirsch
Coarse	10	60	158	220	2.20	−27 %
Medium	40	240	562	280	2.80	−7 %
Fine	80	960	2 100	295	2.95	−1.7 %
Very fine + Q8	160	3 840	23 800	302	3.02	+0.7 %

Lessons. (1) K_t is a converged quantity, not the first value you read. (2) A factor of 2 element refinement around the hole gave a 6 % improvement in K_t — well-behaved h-convergence. (3) Going from linear to quadratic shapes is more effective than another factor of 2 h-refinement at the same DOF count. (4) Always report K_t with the mesh that produced it so the reviewer can judge convergence.

6.3 Example C — Nonlinear contact (turbine-blade fir-tree root)

Problem. A HIPed nickel-superalloy turbine blade is retained in a disk slot by a multi-flank “fir-tree” root. Operating centrifugal force at 12 000 rpm pulls the blade radially outward against the disk slot faces. Compute peak contact pressure and worst stress concentration for design fatigue check.

Model setup.

• 3-D Abaqus/Standard model, ~250 k DOFs in the local fir-tree zone (refined HEX8 with B-bar), submodel of the global blade FEM.
• Contact pair: blade-root surface (slave) ↔ disk-slot surface (master), surface-to-surface discretization, finite sliding.
• Friction coefficient μ = 0.30 (typical Ni-superalloy / Ni-superalloy at 600 °C).
• Normal stiffness: hard contact with augmented Lagrange (Abaqus default — penalty with periodic adjustment of penalty to drive penetration < 1 % of element size).
• Loading: 30 quasi-static increments from zero to 12 000 rpm centrifugal body force.
• Solver: Newton–Raphson with line search; automatic time-stepping with initial Δt = 0.05, min 1e-5, max 0.1.

Convergence. Each increment converges in 6–10 Newton iterations. Two restarts needed: one at 60 % load (contact open-close oscillation, fixed with smaller initial increment), one at 95 % load (large friction redistribution, fixed by enabling unsymmetric solver because friction stiffness matrix is non-symmetric for μ > 0).

Results.

Quantity	Peak value	Location
Max contact pressure	1 350 MPa	Lower flank, mid-span, near the fillet corner
Max von Mises stress	980 MPa	At the same fillet corner, ~0.2 mm sub-surface
Slip distance	18 μm	Top flank (less loaded)
Coefficient of utilization (σ_vM / σ_y) at peak	0.84	(σ_y = 1 170 MPa at 600 °C for IN-718-equivalent)

Lessons. (1) Contact pressure peaks at corners; design the flank radii to spread the load — chamfer the load-bearing corner if possible. (2) Friction makes the stiffness matrix non-symmetric; use UNSYMM solver or μ < 0.1 will quietly produce non-convergent residuals. (3) Hourglass-energy check: 2.1 % of total internal energy — acceptable. (4) Fatigue life is governed by sub-surface von Mises and surface contact pressure together (fretting fatigue) — feed into fatigue-analysis and a Findley or Smith–Watson–Topper multi-axial criterion.

7. Nonlinear FEA

Linear FEA assumes small displacements, linear material, no contact, and a single solve K·u = F. Real problems break all four.

7.1 Geometric nonlinearity

When displacements are no longer small relative to characteristic length (cable structures, tension membranes, post-buckling, large-displacement plasticity), the strain measure must be finite-deformation (Green–Lagrange E = ½ (F^T F − I), or logarithmic, or rate-of-deformation). The stress measure becomes the second Piola–Kirchhoff or Cauchy stress, pulled back as needed. The equation set becomes nonlinear: K depends on u.

Three kinematic frameworks:

• Total Lagrangian (TL) — formulate everything on the reference (undeformed) configuration. Strain = E (Green–Lagrange), stress = S (second PK). Default in some Abaqus elements.
• Updated Lagrangian (UL) — update reference configuration each increment. Strain = ε (rate-of-deformation, logarithmic), stress = σ (Cauchy). Default in most explicit codes.
• Arbitrary Lagrangian–Eulerian (ALE) — mesh moves but not pinned to material. Used for free-surface flows, large metal-forming (rolling, extrusion), and erosion.

In commercial codes, geometric nonlinearity is a toggle: NLGEOM, ON in Abaqus, NLGEOM, ON in ANSYS APDL, the LARGE_DISP keyword in MSC Nastran SOL 400.

7.2 Material nonlinearity

Class	Constitutive law	Typical use
J2 plasticity (von Mises)	f = σ_vM − σ_y(ε̄_p) ≤ 0; flow rule dε_p = dλ · ∂f/∂σ	Steels, aluminum, copper
Drucker–Prager	f = √J₂ + α I₁ − k ≤ 0	Soils, rock, concrete
Mohr–Coulomb	τ = c + σ_n tan φ	Soils, slope stability
Hill anisotropic	quadratic in σ with 6 coefficients	Rolled sheet, textured metals
Hyperelastic (Mooney–Rivlin, Ogden, Yeoh, Arruda–Boyce)	W(I_1, I_2, J)	Rubber, elastomers, soft tissue
Viscoelasticity	Prony series G(t), K(t)	Polymers, rubber damping
Viscoplasticity	Perzyna, Chaboche	High-strain-rate metal, creep
Continuum damage (CDM)	σ = (1 − D) D̄ : ε	Concrete cracking, ductile damage
Cohesive zone	traction-separation T(δ)	Delamination, adhesive joints
Phase-field fracture	crack as smeared scalar field	Brittle and ductile fracture

The constitutive integration is implicit-incremental: at each increment, given Δε, solve for σ and the updated internal variables using a return-mapping algorithm (Simo & Hughes 1998 is the canonical reference). The algorithmic tangent (consistent tangent operator) — distinct from the continuum tangent — is essential for quadratic Newton convergence.

7.3 Contact nonlinearity

Three enforcement methods:

• Penalty — add k_n · g (gap) to the residual; cheap but allows penetration ~ F/k_n. Penalty stiffness must be tuned: too soft → penetration, too stiff → conditioning and convergence problems. Typical k_n ≈ 10–100 × neighbour-element stiffness.
• Lagrange multiplier — exact constraint (g = 0) by adding a multiplier λ as DOF. Expands the system; produces a saddle-point problem.
• Augmented Lagrange — penalty inner loop + multiplier update outer loop (Uzawa-type). Default in Abaqus, ANSYS Mechanical for “hard” contact.

Discretization:

• Node-to-surface (NTS) — slave node projected onto master surface. Cheap; can miss; one-pass.
• Surface-to-surface (STS) — both slave and master surfaces integrated over. Smoother pressures; default in modern codes for solids.
• Segment-to-segment (mortar) — full mortar projection of slave onto master; smoothest, most expensive.

Friction: Coulomb (μ_s = μ_k, single coefficient) is default; rate-dependent and stick-slip available. Sliding friction makes the global stiffness matrix non-symmetric — use an unsymmetric solver (Abaqus UNSYMM, ANSYS NROPT, UNSYM) or convergence stalls.

7.4 Buckling

• Linear (eigenvalue) buckling — solve (K + λ K_g) · φ = 0 where K_g is the geometric stiffness from a reference load. Gives buckling load factors λ_i and mode shapes. Fast but optimistic by 1.5×–3× — does not include imperfections, plasticity, or load redistribution.
• Nonlinear (Riks arc-length) — incremental-iterative with arc-length control instead of load or displacement; tracks the equilibrium path through limit and turning points. Required for post-buckling, snap-through, snap-back.

7.5 Newton–Raphson and friends

The unbalanced-force equation R(u) = F_ext − F_int(u) = 0 is solved by Newton:

K_T · Δu = R(u_k)        with K_T = ∂F_int / ∂u  (tangent stiffness)
u_{k+1} = u_k + Δu

iterating until ‖R‖ / ‖F_ext‖ < tol (typically 1e-3 to 1e-6). Variants:

• Full Newton — recompute K_T each iteration; quadratic convergence near the solution.
• Modified Newton — keep K_T fixed for many iterations; cheaper per iteration, more iterations.
• Quasi-Newton (BFGS) — update K_T by rank-one or rank-two correction; widely used in research codes.
• Line search — scale Δu by α ∈ (0, 1] to reduce the residual.
• Arc-length — control a generalized “arc length” parameter to traverse limit points (Riks; Crisfield).

7.6 Dynamics — implicit vs explicit

Implicit time integration (HHT-α, Newmark β, Bossak, generalized-α): solve K_eff · u_{n+1} = F_eff each step, where K_eff = K + a₀ M + a₁ C. Step size limited by accuracy, not stability — typical Δt ≈ T_lowest_mode / 20 = milliseconds for structural transients. Robust, unconditionally stable for linear cases. Used in: low-frequency vibration, seismic response, creep, slow contact problems.

Explicit time integration (central difference): u_{n+1} = 2 u_n − u_{n−1} + Δt² M⁻¹ (F − K u_n); no matrix factorization, only matrix-vector products against a lumped (diagonal) mass matrix. Step size limited by stability — Courant condition Δt < L_min / c, where c = √(E/ρ). Typical Δt ≈ 1 μs for automotive crash (L_min ≈ 5 mm, c ≈ 5000 m/s). Used in: crash, blast, drop impact, sheet-metal forming, ballistics, fragmentation.

Aspect	Implicit	Explicit
Step size	Accuracy-limited (large)	Stability-limited (tiny)
Cost per step	High (matrix solve)	Low (mat-vec only)
Total steps	Few (10²–10⁴)	Many (10⁵–10⁷)
Memory	High (factor K)	Low (no factor)
Stability	Unconditional (most schemes)	Conditional (CFL)
Contact	Hard for many bodies	Trivial
Best for	Structural vibration, slow plasticity, creep	Crash, impact, blast, forming
Commercial	Abaqus/Standard, ANSYS Mech, MSC Nastran SOL 109/112, NX Nastran	LS-DYNA, Abaqus/Explicit, RADIOSS, PAM-CRASH

8. Common analysis types

Analysis	Solver path	Output
Static linear	K · u = F, one direct solve	Displacements, stresses, reactions
Static nonlinear	Incremental + Newton iteration	Load-displacement curve, residual stress
Modal (eigenvalue)	(K − ω² M) · φ = 0	Natural frequencies, mode shapes
Pre-stressed modal	Modal on K + K_g(σ_0)	Frequencies under preload (turbine blade spin, bolted joint)
Frequency response (harmonic)	(−ω² M + iω C + K) · u = F at each ω	Complex displacement vs frequency
Random vibration (PSD)	Cross-PSD propagation through modal basis	RMS stress, fatigue input
Transient (time-history)	Implicit or explicit integration	u(t), σ(t), v(t), a(t)
Linear buckling	(K + λ K_g) · φ = 0	Buckling factors λ, mode shapes
Nonlinear buckling	Arc-length (Riks)	Limit-load, post-buckling path
Steady thermal	K_T · T = Q	Temperature field
Transient thermal	C_T · dT/dt + K_T · T = Q	T(x, t)
Thermal-mechanical coupled	T → ε_th → σ; one-way or two-way	Thermal stress, distortion
Fluid-structure interaction (FSI)	Coupled CFD + FEM, partitioned or monolithic	Aeroelastic flutter, pump impellers, heart valves
Acoustics / vibro-acoustic	Helmholtz / coupled structural-acoustic	Sound pressure, radiation efficiency
Electromagnetic	Maxwell / quasi-static FEM	B, H, current density, eddy losses
Topology optimization	SIMP / level-set / homogenization	Optimal material distribution
Shape optimization	Free-form / morphing	Smoothed shape variant
Sizing optimization	DOE / gradient / GA	Optimal thickness, cross-section

Multi-physics is increasingly the default: COMSOL Multiphysics is the specialty platform (Joule heating, piezoelectric, MEMS, electrochemistry, plasma, microfluidics); ANSYS Workbench couples Mechanical + Fluent + Maxwell + Icepak through System Coupling; Abaqus Co-simulation Engine connects to Star-CCM+, Adams, MATLAB.

9. Verification and validation (V&V)

The terms are precisely defined and mean different things — engineers and managers confuse them constantly.

9.1 Definitions (ASME V&V 10-2006)

Verification — Are we solving the equations right? Code-to-code comparison, mesh-convergence study, comparison to analytic solutions, patch test, single-element tests. Purely mathematical / numerical exercise.

Validation — Are we solving the right equations? Comparison of the model to physical experiments. Includes test-uncertainty, model-form error, and parameter uncertainty.

Uncertainty quantification (UQ) — propagation of input distributions (material scatter, dimensional tolerance, load uncertainty) through the model to output distributions. Polynomial chaos, Monte Carlo, stochastic collocation, surrogate models (Kriging, Gaussian process) are the workhorses.

9.2 Standards and benchmarks

Standard / Suite	Scope
ASME V&V 10-2006	Computational solid mechanics V&V
ASME V&V 20-2009	CFD and heat transfer V&V
ASME V&V 40-2018	Medical-device computational modelling (FDA submission)
NAFEMS benchmark library	Linear, nonlinear, dynamics, contact, thermal — every code is tested against these
NAFEMS R0026	FEA quality-management procedures
NAFEMS QSS001	Quality system for simulation
ISO 9001:2015	General QMS, applied to FE-using engineering process
AIAA G-077-1998	Guide for V&V in CFD

9.3 The patch test

The minimum element-level verification: a small mesh of distorted elements, with boundary nodal displacements set to produce a uniform strain field, must satisfy:

1. Every interior node is in equilibrium (Σ F = 0) to machine precision.
2. Every element reports the prescribed uniform strain and stress.

If the patch test fails, the element is broken — it cannot converge to the correct solution under refinement. Every reputable commercial element is patch-tested by the vendor and re-tested by NAFEMS.

9.4 Mesh-convergence study (h-refinement)

For a target response quantity Q (peak stress, displacement, frequency, J-integral):

1. Solve on mesh 1 (coarse), record Q_1, mesh size h_1.
2. Refine globally or locally; solve on mesh 2, record Q_2, h_2.
3. Refine again; solve on mesh 3, record Q_3, h_3.
4. Apply Richardson extrapolation assuming Q ≈ Q_∞ + C · h^p:

p ≈ log( (Q_1 − Q_2) / (Q_2 − Q_3) ) / log( h_1 / h_2 )
Q_∞ ≈ Q_3 + (Q_3 − Q_2) / (r^p − 1)         (r = h_2 / h_3)

Report Q_∞ ± GCI (Grid Convergence Index, Roache 1994). A converged result with GCI < 5 % is the lowest bar for credibility.

10. Tools / software

10.1 Commercial — general-purpose FEA (Tier 1)

Tool	Vendor	Strengths
ANSYS Mechanical	ANSYS / Synopsys	Industry-standard structural; deep nonlinear, contact, vibration, fatigue, multi-physics via Workbench
Abaqus (Standard + Explicit)	Dassault Systèmes / SIMULIA	Best-in-class nonlinear material, contact, explicit dynamics; aerospace + auto + research
MSC Nastran / Simcenter Nastran	Hexagon / Siemens	Aerospace structural-analysis lingua franca; aeroelasticity; SOL sequences
Altair HyperWorks (OptiStruct)	Altair	Topology optimization; strong dynamics; HyperMesh meshing
Hexagon Apex	Hexagon	Modernized Nastran front-end with direct CAD modelling
Siemens NX CAE	Siemens	CAD-integrated structural, motion, durability
MSC Marc	Hexagon	Heritage nonlinear, contact, large deformation; rubber & polymer industry
ADINA	Bentley	Bathe-school formulation; strong nonlinear and FSI; research-grade

10.2 Commercial — multi-physics

Tool	Vendor	Strengths
COMSOL Multiphysics	COMSOL	Built for coupling; equation-based modelling; MEMS, electrochemistry, plasma
ANSYS Workbench	ANSYS	System Coupling for Mechanical ↔ Fluent ↔ Maxwell ↔ Icepak
Abaqus Co-simulation Engine	Dassault	External coupling to Fluent, Star-CCM+, Adams, MATLAB
CST Studio Suite	Dassault	High-frequency EM; antennas, EMC

10.3 Commercial — explicit dynamics

Tool	Vendor	Strengths
LS-DYNA	ANSYS	Crash, blast, drop, sheet-metal forming; thousands of material models
Abaqus/Explicit	Dassault	Same problem space, tighter integration with Standard
RADIOSS	Altair	European auto-crash standard
PAM-CRASH	ESI	Older auto-crash, still in airbag/occupant niches
AUTODYN	ANSYS	Hydrocode for ballistics, fragmentation, hypervelocity impact

10.4 Commercial — CFD

Tool	Vendor	Strengths
ANSYS Fluent	ANSYS	Industry workhorse, broad turbulence and combustion
Star-CCM+	Siemens	Polyhedral mesher; tightly integrated workflow
OpenFOAM	OpenFOAM Foundation / ESI / OpenCFD	Open-source FVM; widely used in research and consulting
CFX	ANSYS	Turbomachinery focus
SU2	Stanford / open source	Aerodynamic shape optimization, adjoint

10.5 Topology + shape optimization

• Altair OptiStruct — SIMP, MMA, topology + shape + sizing.
• Tosca (Dassault) — non-parametric topology and shape, plugs into Abaqus and ANSYS.
• nTopology — implicit geometry + lattice synthesis + FE.
• ANSYS Mechanical Optimization — built-in level-set and SIMP.

10.6 Mesh generation

• ANSA (Beta CAE) — best-in-class for auto/aero crash and CFD meshing.
• HyperMesh (Altair) — same tier; broad solver support.
• Pointwise (Cadence) — CFD-focused, structured + unstructured.
• Cubit / Trelis (Coreform / Sandia) — hex meshing of complex CAD.
• Gmsh — open-source; tetra/hex/quad, scriptable.

10.7 Open-source

• CalculiX — Abaqus-syntax solver and pre/post-processor, GPL.
• code_aster (EDF) — French nuclear-grade, comprehensive material models.
• Elmer (CSC Finland) — multi-physics, GPL.
• FEBio (UCSD) — biomechanics specialty, soft tissue and contact.
• MFEM (LLNL) — high-order FEM C++ library; spectral/h/p.
• FEniCSx (University of Oxford / Cambridge) — Python/C++ symbolic FEM for research.
• deal.II — C++ adaptive-mesh-refinement FEM library.
• MoFEM (Glasgow) — flexible Galerkin, research.

10.8 Cloud / embedded / education

• Onshape Simulation — browser-native CAD with built-in linear static.
• Autodesk Fusion 360 Simulation — CAD-integrated linear and modal.
• SimScale — browser-hosted Nastran, CalculiX, OpenFOAM.
• Rescale — HPC cloud for large commercial FE jobs.
• MSC Apex — modernized Nastran direct-CAD environment.

10.9 Pre / post-processing

• ParaView (Kitware) — open-source post-processor; the de-facto research standard, used by LLNL, Sandia, NASA.
• Tecplot 360 — commercial CFD-focused post.
• EnSight (ANSYS) — large-data visualisation, originally for LS-DYNA.
• HyperView (Altair) — auto / aero post-processing.

11. Cross-references

• mechanics-of-materials — first-principles stress, strain, Hooke’s law; FEA’s hand-calc sanity check.
• beam-theory — Euler–Bernoulli and Timoshenko beam elements derive from this.
• structural-analysis — matrix stiffness method is the bridge from frame analysis to general FEM.
• materials-steel, materials-aluminum, materials-composites — provide E, σ_y, σ_u, ply stiffness used in FE input decks.
• vibration-dynamics — modal analysis, frequency response, transient dynamics on the FE mesh.
• heat-transfer — thermal FEA conduction-convection-radiation problems.
• fluid-mechanics — FVM CFD heritage; FSI partner of structural FE.
• electromagnetics-engineering — low-frequency magnetostatics FEM, high-frequency wave FEM.
• fracture-mechanics — J-integral, K_I extraction, XFEM (extended FEM with discontinuous enrichment).
• fatigue-analysis — FE stress field is the input to S-N / ε-N / multi-axial fatigue.
• statics-fundamentals — equilibrium prerequisite.
• planned structural-dynamics — seismic and vibration-isolation FE workflows.
• planned cfd-deep — finite volume vs finite element CFD.
• planned scientific — Abaqus .inp, Nastran .bdf, ANSYS .cdb, LS-DYNA .k, CGNS, EXODUS II, VTK file formats.

12. Citations

1. Bathe, K.-J. Finite Element Procedures, 2nd ed. Prentice Hall / KJB, 2014. ISBN 978-0979004957. Canonical rigorous treatment of linear and nonlinear FEM; the textbook used at MIT and ETH for graduate FE courses.
2. Hughes, T. J. R. The Finite Element Method: Linear Static and Dynamic Finite Element Analysis. Prentice-Hall, 1987 (Dover reprint 2000). ISBN 978-0486411811. The other rigorous-modern foundation text; especially strong on time integration and variational principles.
3. Zienkiewicz, O. C.; Taylor, R. L.; Zhu, J. Z. The Finite Element Method, 7th ed., 3 volumes (Vol. 1 Basis, Vol. 2 Solid Mechanics, Vol. 3 Fluid Dynamics). Butterworth-Heinemann, 2013. ISBN 978-1856176330. The encyclopedia of the field; first edition 1967 by Zienkiewicz alone.
4. Cook, R. D.; Malkus, D. S.; Plesha, M. E.; Witt, R. J. Concepts and Applications of Finite Element Analysis, 4th ed. Wiley, 2002. ISBN 978-0471356059. The practical-application companion to the rigor of Bathe and Hughes.
5. Belytschko, T.; Liu, W. K.; Moran, B.; Elkhodary, K. I. Nonlinear Finite Elements for Continua and Structures, 2nd ed. Wiley, 2014. ISBN 978-1118632703. The canonical nonlinear-FE reference; large deformation, contact, plasticity, explicit dynamics.
6. Reddy, J. N. An Introduction to Nonlinear Finite Element Analysis, 2nd ed. Oxford University Press, 2014. ISBN 978-0199641758. Compact treatment of geometric and material nonlinearity.
7. Holzapfel, G. A. Nonlinear Solid Mechanics: A Continuum Approach for Engineering. Wiley, 2000. ISBN 978-0471823193. Theoretical underpinning for finite-deformation FE; rubber, biological tissue, large-strain plasticity.
8. Crisfield, M. A. Non-linear Finite Element Analysis of Solids and Structures, Vol. 1 (Essentials) and Vol. 2 (Advanced Topics). Wiley, 1991, 1997. ISBN 978-0470666449 (combined reissue). The arc-length method and continuation algorithms reference.
9. Simo, J. C.; Hughes, T. J. R. Computational Inelasticity. Springer, 1998. ISBN 978-0387975207. The canonical reference on return-mapping algorithms and algorithmic tangents for plasticity.
10. Galerkin, B. G. “Series solution of some problems in elastic equilibrium of rods and plates.” Vestnik Inzhenerov, 1915. The original Galerkin-method paper (in Russian).
11. Clough, R. W. “The Finite Element Method in Plane Stress Analysis.” Proceedings, 2nd ASCE Conference on Electronic Computation, Pittsburgh, 1960. First published use of the term “finite element.”
12. Argyris, J. H. “Energy theorems and structural analysis.” Aircraft Engineering, 1954–1955 (series). The matrix structural analysis foundation for FEM.
13. Turner, M. J.; Clough, R. W.; Martin, H. C.; Topp, L. J. “Stiffness and Deflection Analysis of Complex Structures.” Journal of the Aeronautical Sciences, vol. 23, no. 9, 1956, pp. 805–823. The Boeing paper widely credited with launching applied FEM.
14. Irons, B. M. “Engineering applications of numerical integration in stiffness methods.” AIAA Journal, vol. 4, no. 11, 1966. Reduced integration and the patch test.
15. Flanagan, D. P.; Belytschko, T. “A uniform strain hexahedron and quadrilateral with orthogonal hourglass control.” International Journal for Numerical Methods in Engineering, vol. 17, 1981, pp. 679–706. The standard hourglass-control formulation.
16. Bathe, K.-J.; Dvorkin, E. N. “A four-node plate bending element based on Mindlin/Reissner plate theory and a mixed interpolation.” IJNME, vol. 21, 1985, pp. 367–383. The MITC4 shell-element paper.
17. Newmark, N. M. “A method of computation for structural dynamics.” Journal of the Engineering Mechanics Division (ASCE), vol. 85, 1959, pp. 67–94. The Newmark β family of time-integration schemes.
18. Roache, P. J. “Perspective: A Method for Uniform Reporting of Grid Refinement Studies.” Journal of Fluids Engineering, vol. 116, 1994, pp. 405–413. The Grid Convergence Index.
19. NAFEMS Benchmark Library and Procedures. International Association for the Engineering Modelling, Analysis and Simulation Community. nafems.org.
20. ASME V&V 10-2006 (R2016) — Guide for Verification and Validation in Computational Solid Mechanics. ASME, 2006.
21. ASME V&V 20-2009 (R2021) — Standard for Verification and Validation in Computational Fluid Dynamics and Heat Transfer. ASME, 2009.
22. ASME V&V 40-2018 — Assessing Credibility of Computational Modeling through Verification and Validation: Application to Medical Devices. ASME, 2018.
23. ASME BPVC Section VIII Division 2 — Rules for Construction of Pressure Vessels: Alternative Rules. American Society of Mechanical Engineers, 2023 edition. The “Design by Analysis” code that explicitly mandates FE.
24. EN 1993-1-5:2024 — Eurocode 3: Design of steel structures — Part 1-5: Plated structural elements, Annex C — Finite element methods of analysis. CEN.
25. ANSYS Mechanical Theory Reference and Element Reference, current release. ANSYS Inc.
26. Abaqus Theory Guide and Abaqus Analysis User’s Guide, current release. Dassault Systèmes / SIMULIA.
27. MSC Nastran Quick Reference Guide and Theoretical Manual, current release. Hexagon.
28. LS-DYNA Theory Manual, current release. ANSYS / LST.
29. CalculiX User’s Manual. Guido Dhondt, dhondt.de.
30. code_aster Reference Manual. EDF R&D. code-aster.org.
31. MFEM: Modular Finite Element Methods Library. Lawrence Livermore National Laboratory. mfem.org.

---

Session log:

node ~/.claude/bin/obsidian-research.mjs log "Built Engineering/fem-fea.md Tier 2 deep note"