15.2 Computational Methods: FEA & Numerical Analysis

15.2 Computational Methods: FEA & Numerical Analysis

Modern machine design relies on computer-aided engineering (CAE), but the exam tests whether you understand what the software is doing: how a finite element model is built, when its answers can be trusted, and which numerical methods sit behind it. Hand-checkable numerical techniques also appear directly in problems.

---

Finite Element Analysis (FEA) Fundamentals

Finite Element Analysis discretizes continuous structures into finite subdomains (elements) to solve boundary-value stress, thermal, and dynamic problems.

Theoretical Formulation: The Stiffness Matrix Equation

In structural FEA, equilibrium is governed by the global stiffness equation:

$$[K] \{u\} = \{F\}$$

where $[K]$ is the global structural stiffness matrix ($N \times N$, symmetric and banded), $\{u\}$ is the vector of unknown nodal displacements, and $\{F\}$ is the vector of applied nodal forces and boundary reactions.

The element stiffness matrix $[k]^e$ is derived using the principle of virtual work:

$$[k]^e = \int_{V^e} [B]^T [D] [B] \, dV$$

where $[B]$ is the strain-displacement matrix (spatial derivatives of shape functions $N_i$) and $[D]$ is the material constitutive elasticity matrix (containing Young's modulus $E$ and Poisson's ratio $\nu$).

Element Types & Formulations

1. 1D Elements:
   • Truss / Bar: 2 nodes, transmits pure axial force ($T_x$).
   • Beam (Euler-Bernoulli / Timoshenko): Transmits axial force, shear force, bending moments, and torsion. Timoshenko formulation incorporates shear deformation for deep beams.
2. 2D Continuum Elements:
   • Plane Stress: Assumes zero out-of-plane normal and shear stresses ($\sigma_z = \tau_{xz} = \tau_{yz} = 0$). Applicable to thin plates loaded strictly in their plane.
   • Plane Strain: Assumes zero out-of-plane normal and shear strains ($\epsilon_z = \gamma_{xz} = \gamma_{yz} = 0$). Generates out-of-plane stress $\sigma_z = \nu(\sigma_x + \sigma_y)$. Applicable to long, thick bodies of uniform cross-section (e.g., dams, tunnels, high-pressure pipes, retaining walls).
   • Axisymmetric: Models 3D bodies of revolution subjected to axisymmetric loading using a 2D planar section ($r, z$).
3. 3D Solid Elements:
   • Tetrahedral (Tets): Used for automated meshing of complex organic CAD geometry. First-order (4-node linear) tets must be avoided in stress analysis because they suffer from severe artificial shear locking (overly stiff). Second-order (10-node quadratic) tets possess mid-side nodes and provide excellent accuracy.
   • Hexahedral (Bricks): 8-node linear or 20-node quadratic elements. Exhibit superior computational convergence, avoid shear locking, and are preferred for critical fatigue and contact stress zones.

Mesh Convergence & The Stress Singularity Trap

  Peak Stress
       ^
       |                          /--- Stress Singularity (Sharp Corner r=0)
       |                         /     (Stress -> Infinity as Mesh -> 0)
       |                        /
       |       +---------------+------ Physical Stress Concentration (Fillet r>0)
       |      /                       (Converges to True Peak Kt * sigma_nom)
       |     / 
       |    /
       +---+-------------------------> Number of Elements (Refinement)

• Mesh Convergence Study: The process of systematically refining element size ($h$-refinement) or increasing element polynomial order ($p$-refinement) until the output of interest (peak von Mises stress, displacement, or strain energy) changes by less than a specified threshold (typically $< 2\%$ to $5\%$).
• Physical Stress Concentration vs. Stress Singularity:
  • Physical Stress Concentration: Occurs at a geometric transition with a finite physical radius (e.g., a shaft fillet radius $r = 3\,\text{mm}$ or an elliptical hole). As the mesh is refined, the calculated stress converges smoothly to the true physical stress: $\sigma_{\max} = K_t \sigma_{\text{nom}}$.
  • Stress Singularity: Occurs at idealized mathematical anomalies—such as a point load, a fixed boundary corner constraint, or a sharp re-entrant internal corner with zero radius ($r = 0$). Classical elasticity theory dictates that stress at a sharp crack or zero-radius corner is theoretically infinite ($\sigma \propto r^{-\lambda} \rightarrow \infty$).

Verification vs. Validation (V&V in Computational Mechanics)

• Verification: "Solving the equations right." The mathematical verification that the software correctly solves the underlying discrete model equations (e.g., checking mesh convergence, verifying reaction force equilibrium $\sum F = 0$, confirming energy balance, comparing against closed-form beam solutions).
• Validation: "Solving the right equations." The physical confirmation that the computational simulation accurately reflects real-world physics (e.g., comparing FEA strain and displacement predictions against physical strain gauge telemetry, photoelastic experiments, or full-scale destructive load tests).

Building a Trustworthy Model

• Workflow: geometry (CAD) is simplified (defeatured), meshed, given material properties, loads and boundary conditions, solved, and post-processed.
• Boundary conditions cause most errors. Over-constraining (fixing every degree of freedom on a face) adds false stiffness and creates artificial stress peaks; under-constraining leaves rigid-body motion and a singular stiffness matrix.
• Symmetry: modeling one half or one quarter with symmetry constraints reduces size, but only if geometry, loads and supports are all symmetric.
• Sanity checks: reaction forces must balance the applied loads, deformed shapes must make physical sense, units must be consistent (a model built in millimeters with stresses in MPa needs forces in newtons), and key results should be compared with a hand calculation.

---

Numerical Methods Behind Engineering Software

Root Finding

Newton-Raphson iteration solves $f(x) = 0$ by repeatedly following the tangent:

$$x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)}$$

It converges quadratically near a simple root (the number of correct digits roughly doubles each step) but can diverge from a poor starting guess. Bisection is slower but always converges once the root is bracketed by a sign change.

Example: for $f(x) = x^3 - 2x - 5$ with $x_0 = 2$: $x_1 = 2 - (-1)/10 = 2.1$; $x_2 = 2.1 - 0.061/11.23 = 2.09457$; $x_3 = 2.094551$, already correct to six digits. Nonlinear FEA (contact, plasticity, large deflection) uses the same idea, iterating each load step until the force residual is small.

Numerical Integration

For equally spaced data with step $h$:

$$\text{Trapezoidal rule: } \int y\,dx \approx h\left[\frac{y_0 + y_n}{2} + y_1 + y_2 + \dots + y_{n-1}\right]$$ $$\text{Simpson's 1/3 rule (even number of intervals): } \int y\,dx \approx \frac{h}{3}\left[y_0 + 4y_1 + 2y_2 + 4y_3 + \dots + 4y_{n-1} + y_n\right]$$

Example: a spring test records force $0$, $2.4$, $4.2$, $5.4$ and $6.0\,\text{kN}$ at $0$, $5$, $10$, $15$ and $20\,\text{mm}$. Simpson's rule gives the work $W = (5/3)[0 + 4(2.4) + 2(4.2) + 4(5.4) + 6.0] = (5/3)(45.6) = 76.0\,\text{J}$ (since $\text{kN}\cdot\text{mm} = \text{J}$); the trapezoidal rule gives $5[(0 + 6.0)/2 + 2.4 + 4.2 + 5.4] = 75.0\,\text{J}$.

Differential Equations and Linear Systems

• Euler's method advances $y' = f(t, y)$ by $y_{n+1} = y_n + h\,f(t_n, y_n)$; it is simple but first-order accurate and can go unstable with large steps. The fourth-order Runge-Kutta method averages four slope estimates per step and is the usual general-purpose choice for dynamics and vibration simulations.
• FEA ultimately solves the linear system $[K]\{u\} = \{F\}$ by Gaussian elimination or a related factorization (or by iterative solvers for very large models). Explicit dynamic codes instead march in very small time steps limited by the time a stress wave takes to cross the smallest element.

Errors come from discretization (element size, time step), truncation of the numerical method and round-off; refining the mesh or step and confirming that results stop changing is the universal check.