Structural Mechanics for CAE Engineers
Stress, Strain, FEM & Nonlinear Analysis
Structural mechanics is the mathematical backbone of every FEM simulation. Whether you're checking a bracket for yield or simulating a full automotive crash, the same fundamental equations govern everything your software calculates. This article builds that foundation from first principles — with real engineering context at every step.
1. Stress and Strain — The Two Pillars of Structural Analysis
When I run a simulation in Ansys or Abaqus, I always see these colorful contour plots labeled "stress" and "strain." I know red means bad, but what's the actual difference between the two? Are they just two ways of saying the same thing?
They're closely related but fundamentally different. Think of pulling a rubber band: stress is how hard the material is being pulled internally — force per unit cross-sectional area. Strain is how much it actually deforms — the fractional change in length. You need both because a steel rod and a rubber band can have the same strain, but completely different stresses.
1.1 Normal Stress and Shear Stress
When an internal force $F$ acts on a cross-section of area $A$, the normal stress $\sigma$ (perpendicular to the cut) and shear stress $\tau$ (parallel to the cut) are defined as:
Units: Pa = N/m² in SI. In structural engineering, MPa = N/mm² is far more common. This is a critical unit system point: when your model geometry is in mm and loads in N, stress comes out in MPa automatically — a good setup. If geometry is in m but loads are in kN, be careful.
1.2 Normal Strain and Shear Strain
For a bar of original length $L_0$ that stretches by $\Delta L$:
In large-deformation analyses (crash, metal forming), true strain (logarithmic strain) is used instead:
For small strains (below ~5%), the difference is negligible. For 50% elongation, engineering strain gives 0.5 while true strain gives 0.405 — a 19% difference that matters a lot in forming simulations.
Wait, so the contour plots in Ansys — when it says "Equivalent Stress (von Mises)" — that's not the same as the normal stress in one direction?
Exactly right — and that's one of the most common sources of confusion. The software computes a full 3D stress state (6 components), then collapses it into a single "equivalent" number for display. Von Mises stress is that number. It tells you how close the material is to yielding, regardless of the direction. We'll cover that in Section 3. For now, just remember: the real stress state is a 3×3 tensor, not a single number.
1.3 Unit System Pitfalls in CAE
| Length | Force | Stress | Density | Notes |
|---|---|---|---|---|
| m | N | Pa | kg/m³ | SI — standard in most solvers |
| mm | N | MPa | tonne/mm³ | Most common in structural CAE |
| mm | kN | GPa | tonne/mm³ | Rarely used, large civil structures |
| in | lbf | psi | lbf·s²/in⁴ | US aerospace legacy |
Common mistake: Importing a CAD model in mm, setting E = 200,000 (thinking GPa), but the solver interprets it as MPa because the unit system is mm-N-MPa. Result: a steel structure that deflects like rubber. Always double-check: in the mm-N-MPa system, steel has E = 210,000 MPa, and density = 7.85×10⁻⁹ tonne/mm³.
2. Hooke's Law and Elastic Constants
Hooke's Law — I know it's $F = kx$ for a spring. But in a 3D solid, there must be more to it than one equation, right?
Right. For a 3D isotropic material, Hooke's law becomes a matrix equation relating 6 stress components to 6 strain components. But the key insight is that for an isotropic material, you only need two numbers to describe it completely: Young's modulus $E$ and Poisson's ratio $\nu$. Everything else — shear modulus, bulk modulus — follows from those two.
2.1 The Three Elastic Constants
Young's modulus E — stiffness in tension/compression. For steel, E ≈ 210 GPa. For aluminum, E ≈ 70 GPa. For concrete, E ≈ 30 GPa. Carbon fiber composites can reach 500+ GPa in the fiber direction.
Poisson's ratio ν — lateral contraction under axial load. When you stretch a bar, it gets thinner:
For most metals, ν ≈ 0.25–0.35. For rubber, ν ≈ 0.499 (nearly incompressible). Theoretically, $-1 \leq \nu \leq 0.5$. Auxetic materials (negative ν) exist but are rare.
Shear modulus G — resistance to shear deformation. Not independent:
For steel (E = 210 GPa, ν = 0.3): G = 210/(2×1.3) ≈ 80.8 GPa.
2.2 The 3D Isotropic Constitutive Law
The generalized Hooke's law for an isotropic material in matrix form (Voigt notation):
Inverted, this gives the elasticity matrix D (stiffness form):
This D matrix appears explicitly in FEM stiffness calculations as $\mathbf{K} = \int \mathbf{B}^T \mathbf{D} \mathbf{B}\, dV$.
2.3 Material Properties Reference Table
| Material | E (GPa) | ν | G (GPa) | ρ (kg/m³) | σ_y (MPa) |
|---|---|---|---|---|---|
| Structural steel (S235) | 210 | 0.30 | 80.8 | 7,850 | 235 |
| Stainless steel 304 | 193 | 0.29 | 74.8 | 8,000 | 215 |
| Aluminum 6061-T6 | 68.9 | 0.33 | 25.9 | 2,700 | 276 |
| Titanium Ti-6Al-4V | 113.8 | 0.34 | 42.5 | 4,430 | 880 |
| Cast iron (gray) | 100–170 | 0.26 | — | 7,200 | —* |
| Concrete | 25–35 | 0.20 | — | 2,400 | —* |
| Carbon fiber (unidirectional) | 135 (fiber) | 0.27 | — | 1,600 | 1,500 (fiber) |
| HDPE plastic | 0.8–1.4 | 0.44 | — | 950 | 18–30 |
* Cast iron and concrete are brittle materials — yield stress is not a meaningful design criterion; use ultimate tensile/compressive strength instead.
3. Multiaxial Stress and Yield Criteria
In a tensile test, yield is easy — the material yields at the yield stress. But in a real part, you have stress going in multiple directions at once. How do you even know when it's going to yield?
Great question. The trick is to define an "equivalent" scalar stress that captures the combined effect of all stress components. The two main criteria are von Mises (most common for metals) and Tresca. Both reduce the full 6-component stress state to a single number you can compare against the uniaxial yield stress from a standard tensile test.
3.1 The Stress Tensor
At any point in a 3D body, stress is completely described by a 3×3 symmetric tensor:
Symmetry ($\tau_{ij} = \tau_{ji}$) follows from moment equilibrium, leaving 6 independent components.
3.2 Principal Stresses
We can always rotate the coordinate system to find the principal directions, where all shear stresses vanish. The resulting normal stresses are the principal stresses $\sigma_1 \geq \sigma_2 \geq \sigma_3$, found by solving the eigenvalue problem:
This gives a cubic in $\sigma$ whose roots are $\sigma_1, \sigma_2, \sigma_3$.
3.3 Von Mises Yield Criterion
The von Mises criterion states that yielding occurs when the distortional strain energy reaches a critical value. In terms of principal stresses:
In terms of stress tensor components (as the solver computes it):
Physical interpretation: von Mises stress is zero for hydrostatic stress states (equal compression from all sides), which makes physical sense — you can squeeze a metal hydrostatically to enormous pressures without yielding it. Yielding is driven purely by shape distortion.
I've also heard of the Tresca criterion. When would I use that instead of von Mises?
In practice, von Mises is the default for ductile metals — it fits experimental data better and is smoother mathematically. Tresca is used in pressure vessel codes (ASME Section III) because it's slightly more conservative and was historically easier to apply. For most CAE work, just use von Mises unless your design standard explicitly requires Tresca.
3.4 Tresca Yield Criterion
Yielding when the maximum shear stress reaches the shear yield stress:
Tresca is always more conservative (lower predicted yield load) than von Mises. The maximum difference is 15.5%, occurring in pure shear.
3.5 Mohr's Circle
For plane stress (σ_zz = τ_yz = τ_xz = 0), Mohr's circle is a graphical method to find principal stresses and maximum shear stress:
The circle has center at $\left(\frac{\sigma_{xx}+\sigma_{yy}}{2},\ 0\right)$ and radius equal to $\tau_{\max}$.
4. Beam Theory
I know that beams are a simplified model of slender structures. But I've seen both "Euler-Bernoulli beam" and "Timoshenko beam" in the solver settings. Which one should I use?
Simple rule of thumb: if your beam's length-to-depth ratio is greater than about 10, Euler-Bernoulli is fine. If it's shorter and stockier — like a thick composite sandwich panel or a short machine component — use Timoshenko, which accounts for shear deformation that Euler-Bernoulli ignores. In most shell and 3D element models, this distinction doesn't apply; it's mainly relevant when you're using 1D beam elements.
4.1 Euler-Bernoulli Beam
Assumes plane sections remain plane and perpendicular to the neutral axis (no shear deformation). The governing equation:
where $w$ is transverse deflection, $EI$ is flexural rigidity, and $q(x)$ is distributed load. The bending moment and shear force are:
The bending stress in a beam cross-section:
where $y$ is the distance from the neutral axis and $I$ is the second moment of area.
4.2 Common Deflection Formulas
| Configuration | Max Deflection $w_{\max}$ | Location |
|---|---|---|
| Simply supported, central point load P | $\dfrac{PL^3}{48EI}$ | Midspan |
| Simply supported, UDL q (total W=qL) | $\dfrac{5qL^4}{384EI} = \dfrac{5WL^3}{384EI}$ | Midspan |
| Cantilever, tip load P | $\dfrac{PL^3}{3EI}$ | Free end |
| Cantilever, UDL q | $\dfrac{qL^4}{8EI}$ | Free end |
| Fixed-fixed, central load P | $\dfrac{PL^3}{192EI}$ | Midspan |
4.3 Timoshenko Beam
Includes shear deformation by relaxing the perpendicularity constraint. The rotation of the cross-section $\psi$ is now independent of $dw/dx$:
where $\kappa$ is the shear correction factor (≈ 5/6 for rectangular sections, ≈ 9/10 for circular sections).
5. FEM Mathematical Foundations
I use FEM software every day, but I've never really understood what's happening mathematically. Where does the stiffness matrix $\mathbf{K}$ actually come from?
It comes from the virtual work principle — essentially a weighted average of the equilibrium equation. The key idea is: instead of solving the PDE exactly (which is usually impossible for complex shapes), we approximate the displacement field using polynomial shape functions and minimize residual energy. The result is a linear system $\mathbf{K}\mathbf{u} = \mathbf{f}$. Let me walk you through it step by step.
5.1 Strong Form and Weak Form
The strong form of the equilibrium equation for linear elasticity is:
with boundary conditions $\mathbf{u} = \bar{\mathbf{u}}$ on $\Gamma_u$ (Dirichlet) and $\boldsymbol{\sigma}\cdot\mathbf{n} = \bar{\mathbf{t}}$ on $\Gamma_t$ (Neumann).
The weak form is derived by multiplying by a test function $\mathbf{v}$ (virtual displacement) and integrating over the domain:
This form requires only first derivatives (reduces continuity requirements) and naturally incorporates Neumann boundary conditions — two major advantages over the strong form.
5.2 Shape Functions and Isoparametric Elements
Within each element $e$, the displacement field is approximated by:
where $N_a$ are shape functions and $\mathbf{u}_a$ are nodal displacements.
For a linear 4-node quadrilateral element in natural coordinates $(\xi, \eta) \in [-1,1]^2$:
In an isoparametric element, the same shape functions map both geometry and displacement from reference to physical space. The Jacobian of this mapping:
5.3 The B-Matrix and Element Stiffness
The strain-displacement matrix $\mathbf{B}$ relates nodal displacements to strains:
The element stiffness matrix follows directly from the weak form:
5.4 Gaussian Quadrature
The integral above is evaluated numerically. For a 2D element:
For a 4-node element, 2×2 Gauss points (4 total) integrate the polynomial exactly. An 8-node serendipity element needs 3×3 = 9 Gauss points. Reduced integration (one point fewer per direction) is often used to avoid shear locking in bending-dominated problems — but can introduce spurious zero-energy modes ("hourglass modes") that must be controlled.
So once we have all the element stiffness matrices, how do they come together into the global system?
Assembly — the global stiffness matrix is just the sum of all element contributions, placed at the right rows and columns corresponding to shared degrees of freedom. It's like adding up springs in a network. After applying boundary conditions (which usually means removing rows/columns for constrained DOFs), you solve $\mathbf{K}\mathbf{u} = \mathbf{f}$ with a direct or iterative solver. For a million-DOF model, this system is sparse and symmetric, which solvers exploit heavily.
5.5 Global Stiffness Assembly
The global stiffness matrix is assembled by summing element contributions using the connectivity (incidence) array:
where $\mathbf{L}_e$ is the Boolean connectivity matrix mapping local to global DOFs. After assembly and application of essential boundary conditions, the system $\mathbf{K}\mathbf{u} = \mathbf{f}$ is solved.
6. Nonlinear Analysis Overview
When I try to run a nonlinear analysis in Abaqus, there are all these options: NLGEOM, plasticity, contact. What makes a problem "nonlinear" and why is it so much harder to solve?
Linear analysis assumes: small deformations, elastic material, and no contact changes. Nonlinearity means at least one of those assumptions breaks down. The hard part is that the stiffness matrix $\mathbf{K}$ now depends on the current state — so you can't just solve once. You have to iterate: apply load in increments, solve, update the geometry/material state, solve again, repeat until convergence. That's Newton-Raphson, and it can be expensive and tricky to converge.
6.1 Three Sources of Nonlinearity
- Geometric nonlinearity (NLGEOM): Large deformations cause the stiffness to change — a bent sheet is stiffer than a flat one. The deformed geometry must be tracked. Relevant when displacements are large compared to structure dimensions (typically > 5–10% of characteristic length). Automotive crash, flexible structures, snap-through buckling.
- Material nonlinearity: Stress-strain relationship is nonlinear. Plastic deformation (yielding), rubber elasticity (hyperelastic), rate-dependent creep. The material tangent modulus $E_T$ replaces the elastic $E$ in the stiffness calculation.
- Contact nonlinearity: Bodies may or may not be in contact at any given load step — a fundamentally discontinuous constraint that changes the topology of the stiffness matrix. Gear meshing, bolted joints, press-fit assemblies.
6.2 Newton-Raphson Iteration
The residual (out-of-balance force) is:
Linearizing about the current state:
where $\mathbf{K}_T$ is the tangent stiffness matrix. This converges quadratically when the initial guess is good, but can diverge if the load step is too large or the material/geometry changes rapidly. In practice, load is applied in increments, and convergence is checked via:
My Abaqus job keeps failing with "too many attempts" in a nonlinear step. What's usually causing that?
The most common causes: (1) Load step is too large and the solver can't converge — let Abaqus use automatic time stepping and reduce the minimum step size; (2) Material instability — checking if your plasticity parameters are physically reasonable; (3) Rigid body motion — something is insufficiently constrained and can float freely; (4) Contact issues — try softened contact or adjust penalty stiffness. I'd start by looking at the message file — Abaqus tells you which nodes have the largest residuals, which usually points straight at the problem.
7. Dynamic Analysis Basics
7.1 Equation of Motion
The governing equation for structural dynamics is:
where $\mathbf{M}$ is the mass matrix, $\mathbf{C}$ is the damping matrix, $\mathbf{K}$ is the stiffness matrix, and $\mathbf{f}(t)$ is the time-varying external force vector.
7.2 Eigenvalue Problem and Natural Frequencies
For free undamped vibration ($\mathbf{C} = 0$, $\mathbf{f} = 0$), assuming $\mathbf{u} = \boldsymbol{\phi} e^{i\omega t}$:
The eigenvalues $\omega_i^2$ give the natural frequencies $f_i = \omega_i / (2\pi)$, and the eigenvectors $\boldsymbol{\phi}_i$ are the mode shapes. In practice, only the lowest few dozen modes are computed (Lanczos or subspace iteration methods), as higher modes have little energy participation in typical loading scenarios.
7.3 Rayleigh Damping
Since physical damping matrices are almost never known directly, Rayleigh damping is the standard approximation:
where $\alpha$ and $\beta$ are mass-proportional and stiffness-proportional damping coefficients. For two target modes at frequencies $\omega_i$ and $\omega_j$ with damping ratio $\zeta$:
The modal damping ratio at frequency $\omega_k$ is:
This means Rayleigh damping over-damps low-frequency modes (if $\alpha$ is large) or high-frequency modes (if $\beta$ is large). Choose $\omega_i$ and $\omega_j$ to bracket the frequency range of interest.
7.4 Newmark-β Time Integration
For transient analysis, the equation of motion is integrated step-by-step. The Newmark-β family:
| Parameters | Scheme | Accuracy | Stability |
|---|---|---|---|
| β=0, γ=½ | Central difference (explicit) | 2nd order | Conditionally stable: Δt < Δt_crit |
| β=¼, γ=½ | Average acceleration (implicit) | 2nd order | Unconditionally stable |
| β=1/6, γ=½ | Linear acceleration | 2nd order | Conditionally stable |
Explicit methods (β=0) require no matrix factorization — each step is cheap but the time step is tiny (typically microseconds for crash simulations). Implicit methods (β=¼) allow larger steps but require solving a linear system at each step. For structural dynamics below ~1 kHz, implicit is usually more efficient. For wave propagation and explicit crash (LS-DYNA), explicit wins.
8. Practical Notes and Common Mistakes
What are the most common mistakes that junior engineers make when setting up structural simulations?
In my experience, these come up again and again: (1) Wrong units — always verify by doing a sanity check on deflection or natural frequency analytically; (2) Missing or wrong boundary conditions — check that your model has no rigid body modes by running a modal analysis first; (3) Using too-coarse mesh and trusting the first result — always do a mesh convergence study on the quantity of interest; (4) Not checking contact penetration in nonlinear analyses — large penetration means your contact stiffness is too low; (5) Confusing the coordinate system of material properties in composites.
8.1 Quick Sanity Check Checklist
- Units: Do a hand calculation of deflection or frequency for a simple sub-structure. If the FEM result is 1000× off, it's a unit error.
- Boundary conditions: Run a free-free modal analysis. Should get 6 zero-frequency modes for 3D, 3 for 2D plane, none for axisymmetric.
- Mesh quality: Check aspect ratio (<5:1 preferred), Jacobian (>0.7), and angle limits. Many solvers report these as warnings.
- Results reasonableness: Maximum stress >> yield stress at every node means something is wrong (singularity at boundary condition, element distortion, etc.).
- Load magnitude: Total reaction force at fixed supports should equal applied load exactly (±numerical tolerance).
8.2 Stress Singularities
A classic trap: applying a point load or a perfectly sharp corner to an elastic model will give stress that grows without bound as you refine the mesh — a stress singularity. This is not a physical result; it's a consequence of the mathematical idealization. Solutions:
- Distribute the load over an area (more physical)
- Use sub-modeling: coarse global model → refined local model with realistic load distribution
- Apply fillets to sharp internal corners (the part should have them anyway for manufacturability)
- Accept that peak stress at the singularity is not meaningful; report stress a characteristic distance away
- Stress = force/area (MPa), strain = deformation ratio (dimensionless). Related by Hooke's law: $\boldsymbol{\sigma} = \mathbf{D}\boldsymbol{\varepsilon}$
- Von Mises stress collapses the 6-component stress tensor to a single yield indicator
- FEM stiffness matrix comes from the weak form: $\mathbf{K}^e = \int \mathbf{B}^T\mathbf{D}\mathbf{B}\,dV$
- Nonlinear analysis requires incremental loading and Newton-Raphson iteration
- Dynamic analysis adds mass and damping: $\mathbf{M}\ddot{\mathbf{u}} + \mathbf{C}\dot{\mathbf{u}} + \mathbf{K}\mathbf{u} = \mathbf{f}(t)$
- Always verify with a hand calculation before trusting FEM results
Related Interactive Tools
- Beam Deflection Calculator — Compute deflection, bending moment, shear force diagrams interactively
- Mohr's Circle Visualizer — Enter σx, σy, τxy → principal stresses and max shear
- SDOF Frequency Response — Explore damping ratio effects on resonance amplification
- Goodman Diagram — Check fatigue safety factors under mean + alternating stress