Finite element analysis basics

In short

Finite element analysis replaces a real component — whose governing differential equations cannot be solved by hand — with a mesh of simple pieces (elements) joined at nodes, writes an approximate algebraic equation for each piece, assembles them into one large system [K]{u} = {F}, and solves it for displacements, from which strains and stresses follow. It is the standard way Indian product companies check strength, stiffness, vibration, heat and crash performance before cutting metal, because a simulation costs hours of CPU time while a prototype costs lakhs of rupees. It is examined as a full university subject (AICTE model curriculum carries Finite Element Methods/Analysis as a core or departmental elective in B.E./B.Tech Mechanical) with numericals on element stiffness matrices, shape functions and two- or three-element assemblies; note that FEM is not an explicitly listed topic in the current GATE ME syllabus, so weight it mainly for university exams, placement technical rounds and simulation-job interviews. Employers hiring on it include Tata Motors, Mahindra, Ashok Leyland, Bajaj, L&T, BHEL, HAL, ISRO, DRDO, Bosch, Cummins and the ER&D services firms (L&T Technology Services, Tata Technologies, Cyient, KPIT, Quest Global, Tata Elxsi).

Explained

What FEA is, and the eight steps of a solution

• Idea: a continuum has infinite degrees of freedom; FEA *discretises* it into a finite number of nodal unknowns and assumes a simple (usually polynomial) variation of the field inside each element.
• Mathematical basis: either the principle of minimum potential energy (Π = ½{u}ᵀ[K]{u} − {u}ᵀ{F}, minimised ⇒ [K]{u} = {F}) or the Galerkin weighted-residual / weak form of the governing equation. Rayleigh–Ritz is the ancestor; FEM differs by using *piecewise* trial functions instead of one global function.
• Standard procedure (write this in exams):
  1. Discretise the domain into elements and number the nodes.
  2. Choose the element type and interpolation (shape) functions.
  3. Derive the element stiffness matrix [k]ᵉ and load vector {f}ᵉ.
  4. Assemble all elements into the global [K] and {F} by adding each term into the row/column of its global DOF.
  5. Apply boundary conditions (supports, prescribed displacements) — until this is done [K] is singular, because the body can move as a rigid body.
  6. Solve the reduced system for unknown nodal displacements {u}.
  7. Recover element strains and stresses ({ε} = [B]{u}, {σ} = [D]{ε}) and reactions.
  8. Interpret and validate — convergence check, hand calculation, test data.
• Pre-processor / solver / post-processor is the software split: geometry cleanup and meshing; equation solving (direct Cholesky/sparse for static, iterative or eigen solvers for large/modal); contour plots, animations, probes.

Element types, DOF and shape functions

• 1D: spring/bar (truss) element — 1 axial DOF per node, carries only axial force; beam element — transverse displacement *v* and rotation *θ* per node (2D) or 6 DOF per node (3D frame), carries bending, shear and torsion.
• 2D: CST (constant strain triangle), 3 nodes × 2 DOF; LST 6-node triangle; Q4 4-node quadrilateral; Q8 8-node quad. Used for plane stress (thin plate loaded in its plane — σ_z = 0, e.g. a bracket, gusset plate), plane strain (long body, ε_z = 0, e.g. a dam, long pipe cross-section, extrusion die) and axisymmetric elements (a 2D mesh of the half-section representing a solid of revolution — pressure vessel, flywheel, shaft).
• 3D: tetrahedron (tet4 linear, tet10 quadratic), hexahedron (hex8, hex20), wedge/pyramid — 3 translational DOF per node, no rotational DOF.
• Plate/shell: 4- or 3-noded with 6 DOF per node; used when thickness ≪ other dimensions (sheet-metal body panels, pressure-vessel walls); thickness is a property, not geometry.
• Shape (interpolation) functions Nᵢ: satisfy Nᵢ = 1 at node *i*, 0 at other nodes, and ΣNᵢ = 1 (partition of unity, required to represent rigid-body motion and constant strain). Linear bar: N₁ = 1 − x/L, N₂ = x/L. Quadratic 1D in natural coordinate ξ ∈ [−1, 1]: N₁ = ξ(ξ−1)/2, N₂ = 1 − ξ², N₃ = ξ(ξ+1)/2.
• Isoparametric formulation: the same shape functions interpolate geometry and displacement, so curved edges are possible; a Jacobian maps natural to physical coordinates and the element integrals are evaluated numerically by Gauss quadrature (2-point rule, ξ = ±1/√3, w = 1, integrates up to cubic exactly).
• Linear vs quadratic: linear elements have constant strain/stress inside and need a fine mesh in gradient regions; quadratic elements capture linear stress variation and curved boundaries, so far fewer are needed — a linear tet4 mesh is notoriously stiff and should not be used for stress work.
• Convergence requirements: elements must be complete (contain all lower-order terms) and compatible/conforming (field continuous across edges — C⁰ for bars/solids, C¹ slope continuity for beams); the patch test checks an element formulation.

Assembly, boundary conditions and loads

• Each element's stiffness terms are added into the global matrix at positions given by its node numbers (direct stiffness method); [K] is symmetric, sparse, banded and positive-definite only after restraints.
• Boundary conditions: eliminate the restrained rows/columns, or use the penalty/large-stiffness method, or Lagrange multipliers for prescribed non-zero displacement and MPCs (rigid links, couplings, bolt connectors).
• Loads: nodal forces/moments; pressure, body force and thermal loads are converted into consistent equivalent nodal loads by {f}ᵉ = ∫[N]ᵀ q dV, not guessed. For a beam under UDL *w*: {wL/2, wL²/12, wL/2, −wL²/12}.
• Symmetry halves or quarters the model: apply zero normal displacement on the symmetry plane. Valid only if geometry, material, loads and the expected response are symmetric — never for buckling or asymmetric modes.
• Saint-Venant's principle justifies simplifying load and restraint details far from the region of interest; submodelling then re-analyses a small region with a fine mesh driven by the coarse-model displacements.

Analysis types

• Linear static: small displacement, linear elastic, loads applied slowly; gives stress, deflection, reactions. Valid while stress < yield and deflection ≪ thickness/length.
• Modal (eigenvalue) analysis: ([K] − ω²[M]){φ} = 0 gives natural frequencies and mode shapes; needs no load. Used to keep resonances away from engine/road/blade-passing frequencies. Consistent mass matrix (from ∫ρ[N]ᵀ[N]dV) is more accurate; lumped mass is cheaper and used in explicit codes.
• Harmonic / random vibration / PSD: steady sinusoidal or statistical base excitation; needs damping.
• Linear buckling: ([K] + λ[K_g]){φ} = 0 gives the critical load factor λ; it overestimates real strength because imperfections are ignored.
• Thermal: steady-state or transient conduction with convection/radiation boundary conditions; conductance element kA/L is mathematically identical to a bar's stiffness. Followed by a thermal-stress run where ΔT becomes a load.
• Nonlinear: three sources — geometric (large deflection, snap-through, follower loads), material (plasticity, creep, hyperelastic rubber, composites), contact (gap, friction, pre-loaded bolted joints). Solved incrementally by Newton–Raphson with convergence tolerances.
• Explicit dynamics (LS-DYNA, Radioss, Abaqus/Explicit): crash, drop, forming, ballistics; no stiffness inversion, but the stable time step is limited by the smallest element (Δt ≤ L_min/c, c = √(E/ρ)) — one bad tiny element can double the run time.
• Fatigue / durability: FE stress or strain history fed into S-N or ε-N life calculation (nCode, FEMFAT) with Goodman/Gerber mean-stress correction.

Meshing, element quality and convergence

• Mesh density must be fine where stress gradients are steep (fillets, holes, keyways, weld toes) and may be coarse elsewhere; aim for at least 3–4 quadratic elements through a fillet arc.
• Quality metrics: aspect ratio, skew, included angle, taper, warpage, Jacobian ratio, element size. Typical OEM acceptance limits (they vary by company and code): aspect ratio < 5, quad angles 45°–135°, tria angles 20°–120°, Jacobian > 0.7, warpage < 15°, ≤ 5 % trias in a shell mesh.
• Mesh convergence (h-refinement): refine the mesh in steps and plot the quantity of interest; accept the mesh when the change is below a few per cent. p-refinement raises the polynomial order instead. Displacement-based FEA converges from below — a coarse mesh is too stiff, so it under-predicts deflection and usually under-predicts peak stress.
• Stress singularities: at a sharp re-entrant corner, a point load or a knife-edge restraint, stress rises without limit as the mesh is refined — it never converges. Fix the model (add the real fillet, spread the load, use submodelling), do not report the number.
• Numerical pathologies: shear locking (fully integrated linear elements too stiff in bending), hourglassing (zero-energy modes with reduced integration), volumetric locking (near-incompressible rubber/plasticity — use hybrid or selectively reduced integration).
• Reading results honestly: stresses are computed at Gauss points and extrapolated; the difference between averaged and unaveraged nodal stress is an error indicator — a big jump across element boundaries means the mesh is too coarse. Deformation plots are scaled, not real.

Errors, validation, and Indian codes and practice

• Three error sources: modelling/idealisation error (wrong geometry, material, boundary condition — the largest and most common), discretisation error (mesh), numerical error (round-off, solver tolerance, incomplete convergence of a nonlinear run).
• Verification vs validation: verification asks "did I solve the equations right" (mesh convergence, benchmarks such as the NAFEMS test set, hand calculations); validation asks "are these the right equations" (compare with strain-gauge, modal-hammer or rig-test data).
• Design by analysis in India: pressure vessels are usually designed by rule to IS 2825 (unfired pressure vessels — formula-based, an old standard) or ASME BPVC Section VIII Division 1; where FEA governs the design, Indian fabricators and consultants follow ASME Section VIII Division 2, Part 5 — elastic stress linearisation and categorisation (membrane P_m, membrane + bending P_m + P_b, peak), limit-load and elastic–plastic analysis, and fatigue screening. Piping follows ASME B31.3 with pipe-stress FE software (Caesar II).
• Other Indian touchpoints: IS 800:2007 limit-state steel design and IS 1893 seismic provisions for structural/plant FE work; IS 4000 for high-strength bolted joints; automotive crash and structural requirements under the Automotive Industry Standards (AIS) and CMVR, administered by ARAI/ICAT and tested under Bharat NCAP (effective from 2023); AERB safety codes govern analysis for nuclear components; welded structures draw on IS/ISO weld-quality classes for fatigue classification. Pharmaceutical, nuclear and aerospace work also demands documented, traceable analysis reports under the organisation's quality system.
• Deliverable: an FE report states geometry source and simplifications, material data and its source, mesh statistics and quality, loads and BCs with justification, solver settings, convergence evidence, results against the allowable, and the conclusion — unsupported colour plots are not an engineering report.

Key points

• FEA - numerical solution of field problems by dividing a body into elements connected at nodes.
• Governing equation - [K]{u} = {F}: stiffness × displacement = force; displacement is the primary unknown in structural FEA.
• Node vs element - nodes hold the DOF and results; elements hold the interpolation, material and integration.
• Shape function - interpolates the field from nodal values; Nᵢ = 1 at its own node, 0 elsewhere, ΣNᵢ = 1.
• [B] matrix - strain–displacement matrix, {ε} = [B]{u}; constant for a CST, so strain and stress are constant inside it.
• Element stiffness - [k] = ∫[B]ᵀ[D][B] dV; [D] is the constitutive (material) matrix.
• Singular [K] - the global matrix has no inverse until restraints remove rigid-body motion.
• Plane stress vs plane strain - thin in-plane loaded part (σ_z = 0) vs long prismatic body (ε_z = 0).
• Axisymmetric element - a 2D mesh of the half-section represents a full solid of revolution.
• Isoparametric element - same Nᵢ for geometry and displacement; integrated by Gauss quadrature.
• Consistent load vector - distributed loads converted by ∫[N]ᵀq dV, not split equally by eye.
• Convergence from below - refining the mesh increases deflection towards the true answer.
• Stress singularity - sharp corners and point loads never converge; model the real fillet instead.
• Modal analysis - ([K] − ω²[M]){φ} = 0; frequencies and mode shapes, no load needed.
• Nonlinearity - geometric, material and contact; solved incrementally by Newton–Raphson.
• Implicit vs explicit - implicit inverts stiffness for static/slow problems; explicit marches tiny time steps for crash and impact.
• Quality metrics - aspect ratio, skew, Jacobian, warpage; bad elements corrupt results locally and slow explicit runs.
• Garbage in, garbage out - the boundary conditions and material data decide accuracy far more than the element count.

Formulas and examples

Example 1 — stepped bar, two bar elements (classic exam question). A steel bar (E = 200 GPa = 200 000 MPa) is fixed at node 1. Element 1: A = 400 mm², L = 500 mm. Element 2: A = 200 mm², L = 500 mm. A 20 kN axial pull acts at node 3.

• k₁ = AE/L = 400 × 200 000 / 500 = 160 000 N/mm; k₂ = 200 × 200 000 / 500 = 80 000 N/mm.
• Reduced system (nodes 2, 3): [[240 000, −80 000], [−80 000, 80 000]]{u₂, u₃} = {0, 20 000}.
• Row 2: u₃ − u₂ = 20 000/80 000 = 0.25 mm. Row 1: 240 000u₂ = 80 000u₃ ⇒ u₃ = 3u₂ ⇒ u₂ = 0.125 mm, u₃ = 0.375 mm.
• σ₁ = E(u₂ − u₁)/L₁ = 200 000 × 0.125/500 = 50 MPa; σ₂ = 200 000 × 0.25/500 = 100 MPa (= P/A in each step, so the two-element FE answer is exact here).
• Reaction at node 1 = −k₁u₂ = −20 000 N = 20 kN (balances the load).

import numpy as np
E = 200e3                 # MPa = N/mm^2
A = [400.0, 200.0]        # mm^2
L = [500.0, 500.0]        # mm
K = np.zeros((3, 3))
for e in range(2):        # direct stiffness assembly
    k = (A[e]*E/L[e]) * np.array([[1.0, -1.0], [-1.0, 1.0]])
    K[e:e+2, e:e+2] += k
F = np.array([0.0, 0.0, 20e3])        # N, load at node 3
free = [1, 2]                         # node 0 (node 1) is fixed
u = np.zeros(3)
u[free] = np.linalg.solve(K[np.ix_(free, free)], F[free])
# u = [0.0, 0.125, 0.375] mm

Example 2 — interpolation inside an element. A 100 mm linear bar element has u₁ = 0.020 mm, u₂ = 0.060 mm, E = 200 GPa.

• At x = 25 mm: N₁ = 0.75, N₂ = 0.25 ⇒ u = 0.75 × 0.020 + 0.25 × 0.060 = 0.030 mm.
• ε = (0.060 − 0.020)/100 = 4 × 10⁻⁴ (constant) ⇒ σ = 200 000 × 4 × 10⁻⁴ = 80 MPa, the same everywhere in the element — the reason a single linear element cannot show a stress gradient.

Example 3 — one beam element for a cantilever. L = 1 000 mm, E = 200 000 MPa, I = 1 × 10⁶ mm⁴, tip load P = 1 kN downward. Node 1 fixed, so only v₂ and θ₂ remain (rows/columns 3–4 of [k]):

• (EI/L³)[[12, −6L], [−6L, 4L²]]{v₂, θ₂} = {−1 000, 0}; EI/L³ = 2 × 10¹¹/10⁹ = 200 N/mm.
• From the second equation θ₂ = 3v₂/(2L); substituting gives v₂ = −PL³/(3EI) = −1 × 10¹²/(6 × 10¹¹) = −1.667 mm, θ₂ = −PL²/(2EI) = −1 × 10⁹/(4 × 10¹¹) = −2.5 × 10⁻³ rad (slope downward).
• Identical to the textbook PL³/3EI: a single cubic beam element is *exact* for a tip point load, and with consistent nodal loads the nodal values stay exact even under UDL — but intermediate deflection and the bending-moment variation then need more elements.

Example 4 — mesh convergence at a hole. Thin plate, hole small compared with the width, remote tension 50 MPa; theory gives K_t = 3.0, so peak σ = 150 MPa. Successive Q8 meshes (quarter model, symmetry):

The sequence rises monotonically towards the exact value — convergence from below. Accepting mesh 2 would have under-sized the part by a factor of more than 1.5 on stress.

Example 5 — first natural frequency, FE vs theory. Steel strip cantilever: L = 0.3 m, b = 25 mm, t = 3 mm, E = 200 GPa, ρ = 7 850 kg/m³.

• I = bt³/12 = 25 × 27/12 = 56.25 mm⁴ = 56.25 × 10⁻¹² m⁴ ⇒ EI = 11.25 N·m².
• ρA = 7 850 × 75 × 10⁻⁶ = 0.58875 kg/m; L⁴ = 8.1 × 10⁻³ m⁴.
• EI/(ρAL⁴) = 11.25/(0.58875 × 8.1 × 10⁻³) = 2 359 s⁻²; √ = 48.57.
• ω₁ = 3.516 × 48.57 = 170.8 rad/s ⇒ f₁ = 27.2 Hz. An FE modal run with 10 beam elements should land within about 1 % of this — the standard sanity check before trusting a complex modal model.

Where it is used

• Automotive: body-in-white stiffness and NVH modal targets, crash and BNCAP ratings in LS-DYNA/Radioss, engine-block thermal-structural runs, suspension durability — Maruti Suzuki, Tata Motors, Mahindra, Ashok Leyland, Bajaj, TVS and their Tier-1 suppliers.
• Power, process and heavy engineering: turbine blades and rotors, boiler and pressure-vessel nozzles by ASME Sec VIII Div 2 Part 5, pipe stress to B31.3, lifting and structural steel to IS 800 — BHEL, L&T, Thermax, Godrej & Boyce, EIL.
• Aerospace and defence: launch-vehicle structures, random-vibration and shock qualification, composite laminates — ISRO/VSSC, HAL, DRDO labs, Safran/Collins/Boeing India engineering centres.
• Consumer and electronics: drop-test simulation of appliances and handsets, PCB thermal and vibration analysis, moulded-plastic snap-fit nonlinear contact.
• ER&D services roles: titles such as CAE Engineer, Structural/Durability Analyst, NVH Engineer, Crash Analyst, Thermal Engineer at L&T Technology Services, Tata Technologies, Cyient, KPIT, Quest Global, Tata Elxsi; typical toolchain Hypermesh/ANSA for meshing plus Ansys Mechanical, Abaqus, Nastran, LS-DYNA, Optistruct, nCode.
• Design-office daily use: embedded solvers (SolidWorks Simulation, Creo Simulate, Fusion, Inventor Nastran) for first-cut stiffness, weight optimisation and topology optimisation of castings and brackets.

Interview questions

• Q: Why is the global stiffness matrix singular before boundary conditions are applied? A: Without restraints the structure can translate or rotate as a rigid body at zero strain energy, so [K] has zero-energy modes and no unique solution.
• Q: What properties must a shape function satisfy? A: Unity at its own node and zero at all others, sum of all Nᵢ = 1 at every point, and enough continuity across element edges (C⁰ for solids, C¹ slope continuity for beams) so the element passes the convergence requirements.
• Q: Why does a CST element give poor stress results? A: Its [B] matrix is constant, so strain and stress are uniform inside each element; representing a gradient needs very many elements.
• Q: When would you use plane strain rather than plane stress? A: Plane strain for a long prismatic body restrained along its length (dam, long pipe, extrusion die), plane stress for a thin part loaded in its own plane with free faces.
• Q: Your peak stress keeps rising every time you refine the mesh. What is happening? A: You have a stress singularity at a sharp corner, point load or point restraint; model the actual fillet, distribute the load, or use submodelling — the number itself is meaningless.
• Q: Does a finer mesh always give a larger displacement? A: In standard displacement-based FEA, yes — the discretisation is too stiff, so refinement increases deflection monotonically towards the exact value.
• Q: Difference between consistent and lumped mass matrices? A: Consistent mass comes from ∫ρ[N]ᵀ[N]dV and couples DOF, giving more accurate frequencies; lumped mass puts mass at nodes as a diagonal matrix and is cheaper, standard in explicit codes.
• Q: What are the three types of nonlinearity? A: Geometric (large deflections/rotations, stress stiffening), material (plasticity, creep, hyperelasticity) and contact/boundary nonlinearity.
• Q: Implicit or explicit solver for a car crash, and why? A: Explicit — the event is milliseconds long with severe contact and material nonlinearity; explicit marches very small stable time steps without inverting a stiffness matrix, where implicit iterations would struggle to converge.
• Q: How do you decide the mesh is good enough? A: Run a convergence study on the quantity of interest, check element quality metrics, check that averaged and unaveraged nodal stresses agree closely in the region of interest, and compare with a hand calculation or benchmark.
• Q: Is linear buckling analysis conservative? A: No — it ignores imperfections, residual stress and plasticity, so it over-predicts the collapse load; thin shells need a nonlinear post-buckling run with imperfections.
• Q: Why does a linear tetrahedron mesh give low stresses in bending? A: The tet4 has constant strain and no bending capability within the element, making it artificially stiff; use tet10 or hex elements for stress work.
• Q: You get 420 MPa von Mises in a mild-steel bracket. What do you do next? A: Check that the mesh has converged and that no singularity drives the peak, confirm the material yield and the load case, then compare against yield with the required factor of safety and redesign (fillet, thickness, rib) if the stress is real.

Common mistakes

• Treating a colourful contour plot as proof — never checking convergence, units, or a simple hand calculation such as P/A or PL³/3EI.
• Mixing unit systems: entering E in GPa with lengths in mm and loads in N, or using kg with mm-N without converting mass to tonnes; FE solvers are unit-agnostic and will happily return answers off by 10³.
• Over-restraining the model (fixing all DOF of a bolted flange face) which kills thermal expansion and Poisson contraction and shows false high stresses at the supports.
• Reporting the singular peak stress at a sharp corner or at a point load instead of fixing the idealisation.
• Using linear tet4 elements for a stress or fatigue assessment, or a coarse one-element-thick mesh through a fillet.
• In exam numericals: forgetting to delete the restrained rows and columns before inverting, or placing an element's stiffness terms in the wrong global rows during assembly.
• Assuming symmetry boundary conditions for buckling, vibration or asymmetric loads, which hides exactly the modes being looked for.

Related topics

Mechanics of materials: stress, strain and deflection · CAD and 3D modelling · Vibrations and modal analysis · Machine design: failure theories and factor of safety · Heat transfer: conduction, convection and radiation · Computational fluid dynamics basics