2. Governing equations and discretization errors in explicit crash dynamics

Every explicit crash simulation is an approximation of a continuous initial-boundary-value problem. The solver never integrates the real equations of motion; it integrates a discrete algebraic surrogate. Acceptance decisions for crash evidence therefore depend on understanding exactly which continuous equations are being approximated, how the spatial and temporal discretizations introduce error, and which error sources dominate in the quantities used as evidence (intrusion, acceleration pulses, energy absorption, restraint forces). This section establishes that chain from the strong form of the governing equations to the discrete update performed at each time step, and isolates the spatial and temporal error mechanisms that later acceptance criteria must bound.

Strong form and semi-discrete form

For a deformable body occupying domain Ω with boundary Γ, the momentum balance in the updated Lagrangian frame is

ρ ẍ_i = ∂σ_ij/∂x_j + ρ b_i        in Ω
σ_ij n_j = t_i                     on Γ_t
u_i = ū_i                          on Γ_u

where ρ is the current density, σ the Cauchy stress, b the body force per unit mass, t the prescribed traction, and n the outward normal. The material response closes the system through a constitutive law (elastic-plastic, viscoplastic, damage, or foam crush) and, for crash problems, a contact condition that prevents interpenetration between bodies. The contact constraint is a non-smooth inequality, and it is the single most important reason crash dynamics behaves differently from smooth structural dynamics.

Finite element discretization replaces the continuous displacement field with nodal unknowns u(t) and the weak form with a system of ordinary differential equations:

M ü(t) = f_ext(t) − f_int(u, u̇) − f_contact(u, u̇)
u(0) = u₀,  u̇(0) = v₀

This is the semi-discrete form: space is discretized, time is still continuous. The mass matrix M is usually lumped in explicit crash codes, which makes it diagonal and the update trivially parallel. The internal force f_int is assembled from element stresses evaluated at integration points, and f_contact is evaluated from penetration and relative velocity at active contact segments.

Explicit central-difference update

The explicit solver advances the semi-discrete system with the central difference scheme. With lumped mass, the nodal update is

üⁿ = M⁻¹ (f_extⁿ − f_intⁿ − f_contactⁿ)
u̇ⁿ⁺¹ᐟ² = u̇ⁿ⁻¹ᐟ² + Δt · üⁿ
uⁿ⁺¹   = uⁿ + Δt · u̇ⁿ⁺¹ᐟ²

Two properties follow directly. First, no system of equations is solved, so cost per step scales linearly with element count. Second, the scheme is conditionally stable: it requires the time step to satisfy the Courant–Friedrichs–Lewy (CFL) limit

Δt ≤ Δt_crit = 2 / ω_max

where ω_max is the highest natural frequency of the discretized system. For a linear element of characteristic length L_e and dilatational wave speed c, this is commonly estimated as Δt_crit ≈ L_e / c, with a safety factor of 0.9 or lower applied in practice. The wave speed itself depends on the material:

c = √(E / ρ)          (bar, 1-D)
c = √((K + 4G/3) / ρ) (solid, dilatational)

This coupling is the central engineering fact of explicit crash dynamics: the stable time step is set by the smallest element and the stiffest material in the model, not by the physics of interest. A single 1 mm steel element in a 5 mm mesh model will force the whole model to run at the 1 mm step.

Sources of spatial discretization error

Spatial error arises wherever the discrete mesh cannot represent the continuous field. In crash problems the dominant contributors are:

• Wave dispersion and numerical diffusion. A coarse mesh cannot propagate short-wavelength components at the correct speed. Stress waves smear, and the peak acceleration transmitted through a rail or rocker is under-predicted. This directly biases acceleration-based injury criteria.
• Locking and hourglassing. Fully integrated low-order elements lock in bending and shear; reduced-integration elements are cheaper but admit zero-energy hourglass modes. Hourglass control adds stiffness that is not physical and can stiffen the response.
• Contact discretization. Contact is enforced at nodes or segments, not on the true surface. Penetration depth, contact stiffness, and the choice of penalty versus kinematic enforcement all introduce error in the contact force history, which is often the quantity of interest in occupant and barrier interactions.
• Mesh distortion and element deletion. Severe distortion degrades element accuracy and shrinks the stable time step. Element deletion removes mass and energy from the model, so the energy balance must be monitored as an error indicator.
• Boundary and loading representation. A prescribed velocity applied to a single node is not the same as a distributed crush initiator; the difference is a modeling error, not a numerical one, but it propagates through the same acceptance metrics.

Sources of temporal discretization error

Temporal error has three distinct origins. The first is the truncation error of the central difference scheme, which is second order in Δt for displacement but only first order for velocity at the half step. The second is the stability limit itself: as Δt approaches Δt_crit, high-frequency content becomes marginally stable and can grow. The third, and most consequential in crash work, is mass scaling. When the stable step is unacceptably small, engineers add non-physical mass to small or stiff elements to raise Δt_crit. This changes the inertia of the model and therefore the dynamics. The added mass must be tracked as a percentage of total mass and as a shift in the model's kinetic energy; a common acceptance guard is to keep added mass below 1–2% of total mass and to verify that the global center-of-mass velocity and the total kinetic energy are not materially altered.

Time step choice also interacts with contact. Penalty contact introduces a contact stiffness that can be stiffer than the material, and the stable step must account for it. If the step is too large relative to the contact stiffness, contact chatter and spurious energy appear; if it is too small, the run becomes impractical and mass scaling is invoked, reintroducing the inertia error above.

Error budget and what to monitor

A practical way to keep these error sources under control is to maintain an explicit error budget for each run. The table below lists the quantities that should be extracted from the solver output and the direction in which each error typically biases the result.

The energy balance is the single most useful global diagnostic. For an explicit crash run, the sum of internal, kinetic, contact, and hourglass energy should track the external work within a few percent, and the ratio of hourglass to internal energy should stay small (commonly below 5%). A drifting energy ratio is a direct signal that one of the error sources above is active.

Practical consequence for acceptance

The error structure above explains why convergence acceptance cannot be reduced to a single mesh-size check. A run can be mesh-converged in displacement while still under-predicting the acceleration peak because of dispersion; it can satisfy the CFL limit while carrying 4% added mass that shifts the pulse timing; it can show a clean energy ratio while contact penalty is over-stiffening the barrier force. The acceptance strategy developed in later sections therefore combines a mesh refinement study on the evidence quantities, a time-step and mass-scaling audit, an energy and hourglass budget, and a contact-quality check. Each of these targets a specific term in the error chain shown above, and together they bound the total discretization error to a level that is defensible as crash evidence.

A useful discipline at this stage is to record, for every run, the minimum element size, the resulting stable time step, the applied mass scaling factor and added-mass percentage, the hourglass-to-internal energy ratio, and the total energy ratio at the end of the event. These five numbers, taken together, tell an reviewer which error sources were active and how strongly they acted. They form the quantitative bridge between the governing equations and the acceptance decision, and they are the first thing to inspect when a simulated pulse disagrees with a physical test.