3. Explicit time integration fundamentals
Explicit dynamics solves the semi-discrete equation of motion M a + C v + F_int(u) = F_ext by advancing nodal kinematics one small step at a time, without assembling or factorizing a global stiffness matrix. The dominant algorithm in production codes such as LS-DYNA, Abaqus/Explicit, and Radioss is the central difference method, a second-order accurate, conditionally stable scheme. This chapter develops the update equations, explains why the stable step size is bounded, and connects that bound to the Courant–Friedrichs–Lewy (CFL) condition.
Central Difference Update
Let un, vn, and an denote displacement, velocity, and acceleration at time tn. The central difference approximations of velocity and acceleration are
v_n = (u_{n+1} - u_{n-1}) / (2 Δt)
a_n = (u_{n+1} - 2 u_n + u_{n-1}) / Δt²Substituting into the equation of motion and solving for un+1 gives the displacement update. In practice, codes use the equivalent half-step velocity (leapfrog) form, which is algebraically identical but numerically cleaner:
v_{n+1/2} = v_{n-1/2} + Δt · M⁻¹ (F_ext,n - F_int,n - C v_{n-1/2})
u_{n+1} = u_n + Δt · v_{n+1/2}Two properties matter for impact work. First, the scheme is explicit: un+1 depends only on known quantities, so no linear solve is required. Second, if M is lumped (diagonal), the inverse M⁻¹ is a scalar division per degree of freedom — the entire step is O(N) and vectorizes trivially. This is why explicit codes scale to tens of millions of elements.
Note the staggered storage: displacements live at integer steps, velocities at half steps. This avoids the singularity that a naive same-time central difference would produce for the velocity term when damping is present.
Critical Time Step and Stability
Conditional stability means the scheme blows up if Δt exceeds a limit. For a linear undamped system, modal analysis of the difference operator yields the requirement
Δt ≤ Δt_crit = 2 / ω_maxwhere ωmax is the highest natural frequency of the discretized mesh. Because ωmax is governed by the stiffest, smallest element, the critical step is set by the mesh, not by the loading. For a rod or truss element of length L, density ρ, and Young's modulus E, the element frequency is ωe = 2c/L with wave speed c = √(E/ρ), so
Δt_crit = L / cThis is the physical statement of the CFL condition: a stress wave must not traverse more than one element per step. If it does, the numerical domain of dependence no longer contains the physical one, and information propagates faster than the scheme can resolve it — the classic instability.
E = 210 GPa, ρ = 7850 kg/m³, giving c ≈ 5172 m/s. If the smallest solid element has characteristic length L = 2 mm, then Δt_crit ≈ 0.002 / 5172 ≈ 3.9 × 10⁻⁷ s. A safety factor of 0.9 gives a working step of about 3.5 × 10⁻⁷ s. Simulating 20 ms of response therefore requires roughly 57,000 steps.
CFL Condition in Practice
Real meshes are not uniform, so codes compute a per-element step and take the minimum:
Δt_e = L_e / c_e
Δt_crit = min over elements ( Δt_e )
Δt = k · Δt_crit , typically k = 0.9The characteristic length Le depends on element type: for solids it is volume divided by the largest face area; for shells it is the shortest edge or a diagonal-based measure; for beams it is the element length. The safety factor k absorbs nonlinearity, contact stiffness, and the fact that ωmax grows as elements distort.
Three practical consequences follow. First, a single tiny element — a sliver from a poor mesh — can dominate the global step and make the run prohibitively expensive. Second, mass scaling deliberately increases ρ in small elements to raise Δtcrit, at the cost of added inertia; it must be limited so that kinetic energy stays a small fraction of internal energy. Third, contact penalties add stiffness that can lower the stable step, so contact stiffness is often scaled to remain consistent with Δt.
Stability Limits and Failure Symptoms
Exceeding Δtcrit does not produce a graceful error; it produces exponential growth of high-frequency modes. The tell-tale signs are nodal velocities or energies diverging within a few dozen steps, often localized to the smallest or most distorted element. Distinguishing numerical instability from a physical result is essential:
| Symptom | Likely cause | Action |
|---|---|---|
| Energy grows without bound, no external work | Δt above critical step | Reduce Δt or apply mass scaling |
| Single element shows extreme distortion | Sliver element setting Δt_crit | Remesh or scale its mass |
| Divergence after contact initiation | Penalty stiffness too high | Lower contact stiffness or Δt |
| Gradual energy drift, stable otherwise | Excessive mass scaling | Reduce added mass, check energy ratio |
The energy balance is the primary verification tool: total energy should equal initial energy plus external work minus dissipated energy, within a few percent. A sudden jump in kinetic energy with no corresponding load is the signature of instability.
Choosing and Controlling the Step
Most codes offer fixed and automatic time stepping. Fixed stepping uses a constant Δt computed at initialization; automatic stepping recomputes Δtcrit as elements deform, which is safer for large-deformation impact but introduces variable-step bookkeeping. For impact problems, automatic stepping with a conservative safety factor is the usual choice. The element time step should be monitored in the output; a steadily decreasing step indicates progressive element distortion and warns of impending trouble.
The figure summarizes the two coupled ideas: the staggered central difference update on the time axis, and the CFL bound that ties the allowable step to element length and wave speed. Together they define the cost and the correctness of every explicit impact simulation.