4. Finite element simulation for crash deformation

Finite element (FE) simulation is the computational backbone of vehicle crashworthiness engineering. It allows engineers to reproduce a full-scale barrier impact, a pole side strike, or a pedestrian head impact in software before a single physical prototype is built. For deformation prediction specifically, the FE model must answer three questions accurately: how the structure deforms over time, how the material behaves under high strain rate, and how sensitive the predicted deformation is to the discretization of the mesh. This chapter addresses all three through the lens of explicit dynamics, constitutive material models, and mesh sensitivity analysis.

Why Explicit Dynamics Dominates Crash Simulation

Crash events unfold over 50–150 ms, with peak decelerations reaching 30–80 g. The governing equation is the dynamic equilibrium:

M·a(t) + C·v(t) + F_int(u(t)) = F_ext(t)

where M is the mass matrix, C the damping matrix, F_int the internal nodal force vector (from element stresses), and F_ext the external load vector. In implicit integration, each time step requires solving a large nonlinear system of equations, and convergence is difficult when contact conditions change rapidly or when elements distort severely. Explicit integration sidesteps this entirely: the central difference scheme computes acceleration directly from the equation of motion at the current step, without assembling or inverting a global stiffness matrix.

The explicit time step is governed by the Courant–Friedrichs–Lewy (CFL) condition:

Δt ≤ L_min / c,   c = sqrt(E / ρ)

Here L_min is the smallest element characteristic length, c the speed of sound in the material, E the Young's modulus, and ρ the density. For a typical steel with E = 210 GPa and ρ = 7.85 × 10⁻⁹ t/mm³, the wave speed is approximately 5,170 m/s. A 5 mm element therefore yields a stable time step of roughly 1 μs. A 100 ms simulation thus requires about 100,000 time steps — computationally intensive but highly parallelizable, which is why explicit solvers scale well on GPU and multi-core clusters.

In LS-DYNA, the keyword *CONTROL_TIMESTEP with DTINIT set to 0 lets the solver compute the initial step automatically. Adding *CONTROL_TIMESTEP with TSSFAC = 0.9 provides a 10% safety margin against numerical instability.

Material Models for Crash Deformation

Vehicle crash deformation involves strain rates from 10⁰ to 10³ s⁻¹, large plastic strains exceeding 30%, and complex failure modes including necking, shear localization, and ductile fracture. Linear elastic models are inadequate. Three material model families are standard in crash FE analysis:

Model Typical Keyword (LS-DYNA) Key Parameters Application
Piecewise linear plasticity *MAT_024 True stress–plastic strain curve, strain-rate scaling (LCSS) Mild steel, HSLA steel body panels
Johnson–Cook *MAT_015 A, B, n, C, m, melting temperature High-strength steel, aluminum alloys
Crushable foam *MAT_063 Stress–volumetric strain curve, tensile cutoff Energy-absorbing foam, honeycomb

The Johnson–Cook model is particularly important because it explicitly separates the effects of plastic strain, strain rate, and temperature:

σ = (A + B·ε_p^n) · (1 + C·ln(ε̇/ε̇₀)) · (1 − T*^m)

where A is the yield stress, B the hardening modulus, n the hardening exponent, C the strain-rate sensitivity, and T* the homologous temperature. For a typical DP780 dual-phase steel, representative calibrated values are A = 450 MPa, B = 800 MPa, n = 0.45, C = 0.015, and m = 0.8. These are example values for illustration; actual calibration requires Split Hopkinson Pressure Bar (SHPB) tests at multiple strain rates.

For the deformation prediction task, the critical output from the material model is the equivalent plastic strain (EPS) field. Elements exceeding a defined failure strain (often 0.15–0.30 for automotive steels) are eroded or deleted, simulating material fracture. The choice of failure criterion directly affects the predicted crash mode — too low a failure strain causes premature element deletion and underestimates energy absorption; too high a value allows unphysical stretching.

Mesh Sensitivity: The Hidden Variable

Mesh sensitivity is arguably the most underappreciated source of error in crash FE simulation. The same geometry, material model, and boundary conditions can produce deformation results that differ by 20–40% depending on element size. This occurs because:

  • Element erosion is mesh-dependent. A coarse mesh deletes a larger volume of material per eroded element, artificially reducing the structure's energy absorption capacity.
  • Strain localization is not captured. Shear bands and necking zones require fine discretization (1–3 mm) to resolve; coarse meshes smear these features and overestimate ductility.
  • Contact stiffness scales with element size. Penalty-based contact forces depend on the master segment area, which changes with mesh density.

The standard mitigation strategy is a mesh convergence study. The procedure is straightforward: build the same model at three or four element sizes (e.g., 10 mm, 7 mm, 5 mm, 3 mm), run each to completion, and plot a scalar output — typically maximum intrusion or peak force — against element size. The result is considered converged when further refinement changes the output by less than 5%.

Do not assume convergence from a single refinement step. A model may appear stable between 10 mm and 7 mm but diverge significantly at 5 mm when a new deformation mode (e.g., local buckling) becomes resolvable.

For thin-walled structures, shell elements with 5–7 integration points through thickness are standard. Solid elements are reserved for thick castings and foam. A practical rule: the element size in the deformation zone should not exceed the wall thickness. For a 1.5 mm thick bumper beam, this means 1.5 mm shell elements in the crush zone, transitioning to 5 mm in non-critical areas.

Explicit Crash Simulation Workflow — Deformation Prediction Geometry & Mesh CAD cleanup Shell: 1.5–5 mm Solid: 3–8 mm 5–7 IP through thickness Material Model MAT_024 / MAT_015 Strain-rate curves Failure strain: 0.15–0.30 SHPB calibration Boundary & Contact Rigid wall velocity Self-contact penalty Friction: μ = 0.1–0.3 Initial velocity: 50–64 km/h Explicit Solver Δt ≈ 1 μs 100k steps TSSFAC = 0.9 Deformation Output Intrusion depth (mm) Equivalent plastic strain Energy absorption (kJ) Mesh Convergence — Intrusion vs. Element Size Element size (mm) Intrusion (mm) 10 7 5 3 80 120 160 200 Converged (<5% change) Coarse mesh overestimates intrusion Intrusion depth Convergence band

Practical Calibration and Verification

After the model is built, calibration against physical test data is mandatory. The standard procedure uses a component-level test — for example, a three-point bending test on a bumper beam or a drop tower test on a rocker panel. The simulation is run with the same boundary conditions, and the force–displacement curve is compared. Key calibration targets are:

  1. Peak force magnitude — sensitive to yield stress and strain-rate parameters.
  2. Post-peak force plateau — sensitive to hardening modulus and failure strain.
  3. Total energy absorption — the integral of the force–displacement curve; sensitive to mesh density and erosion criteria.
  4. Deformation mode — visual comparison of fold patterns, buckling locations, and fracture paths.

A common failure symptom is a simulation that predicts a clean, symmetric fold pattern while the physical test shows a mixed mode with fracture. This usually indicates the failure strain is set too high or the mesh is too coarse to resolve the shear band. Reducing element size to 2 mm in the hinge zone and lowering the failure strain by 0.05 typically resolves the discrepancy.

Key takeaway: Explicit dynamics provides the computational framework, material models provide the constitutive physics, and mesh sensitivity analysis provides the confidence that the predicted deformation is not an artifact of discretization. All three must be validated together — a converged mesh with a poorly calibrated material model still produces wrong deformation.

In the next stage of the workflow, these FE-predicted deformation fields become the training data for deep learning surrogates, which is the subject of the following chapters.