2. LS-DYNA explicit solver essentials for crash simulation

Frontal crash simulation in LS-DYNA is built on an explicit central-difference time integration scheme. Unlike implicit solvers that assemble and factorize a global stiffness matrix at every step, the explicit solver advances the solution by evaluating nodal accelerations from internal and external forces, then integrating velocity and displacement. This makes each increment cheap and highly parallelizable, which is exactly what a 1–2 million element vehicle model needs. The trade-off is conditional stability: the time step must stay below a critical value determined by the smallest element and the local wave speed.

Time Integration: The Central-Difference Loop

At each step n, LS-DYNA solves the semi-discrete equation of motion for nodal accelerations:

M · a_n = F_ext,n − F_int,n − C · v_n
a_n     = M⁻¹ · (F_ext,n − F_int,n − C · v_n)
v_(n+1/2) = v_(n−1/2) + a_n · Δt
x_(n+1)   = x_n + v_(n+1/2) · Δt

Here M is the lumped (diagonal) mass matrix, so no matrix inversion is required — the acceleration update is a vector division. Internal forces come from element stresses integrated over the current geometry, which is why the formulation is naturally suited to large deformation, contact-dominated crash events. The scheme is second-order accurate and conditionally stable; it does not require convergence iterations, but it does require that Δt never exceed the stability limit.

Stable Time Step and the Courant Condition

The critical step for a solid element is governed by the Courant–Friedrichs–Lewy (CFL) condition:

Δt_crit = L_e / c
c = sqrt(E / ρ)          (elastic bar wave speed)
Δt = 0.9 · min(Δt_crit)  (default scale factor, TSSFAC = 0.9)

L_e is a characteristic element length (for a hexahedron, roughly volume divided by the largest face area), E is Young's modulus, and ρ is density. For steel with E = 210 GPa and ρ = 7.85e-9 t/mm³, the wave speed is about 5.17 mm/ms. A 5 mm element therefore yields Δt_crit ≈ 0.97 µs, and with the 0.9 scale factor the solver runs at roughly 0.87 µs. A typical 100 ms frontal impact then needs about 115,000 cycles. Because the step is set by the smallest element in the model, a handful of tiny elements in a spot-weld or a fillet can dominate the entire runtime.

Check the initial time step in the LS-DYNA output (d3hsp) right after initialization. The line reporting "smallest timestep" and the element ID that controls it tells you immediately where to focus mesh cleanup or mass scaling.

Mass Scaling: Trading Mass for Runtime

Mass scaling increases the density of selected elements so that their stable step rises to a target value. Since Δt ∝ sqrt(ρ), raising density by a factor of 4 doubles the step. LS-DYNA offers two modes:

  • Fixed mass scaling (*CONTROL_TIMESTEP with DT2MS > 0): every element whose step is below DT2MS has its density scaled up to meet that target. Simple and predictable, but it can add significant artificial mass if the target is aggressive.
  • Selective mass scaling (DT2MS < 0): only elements below the target are scaled, and the added mass is reported per part. This is the preferred approach for production crash models because it limits the mass penalty to the offending regions.

A practical rule is to keep added mass below 1–2% of the total model mass for a frontal crash, and to verify that the added mass is not concentrated in energy-absorbing structures such as the front rail or crash box, where it would artificially stiffen the response. The mass scaling summary in d3hsp and the glstat file both report the percentage of added mass over time.

Energy Balance: The Primary Sanity Check

Because the explicit solver does not enforce equilibrium iteratively, energy balance is the most reliable indicator of solution quality. LS-DYNA tracks the following in glstat:

QuantitySymbolExpected behavior
Kinetic energyKERises from initial velocity, decays as structure deforms
Internal energyIEMonotonically increases as deformation absorbs energy
Hourglass energyHE< 5% of peak IE for reduced integration elements
Sliding interface energySIESmall; large values indicate contact problems
Added mass energy—Reported as % added mass; keep < 2%
Total energyTEShould remain nearly constant (drift < 5%)

The energy balance equation is TE = KE + IE + HE + SIE + .... If total energy drifts upward by more than a few percent, the model is gaining energy numerically — usually from hourglassing, contact penetration, or an unstable time step. If it drifts downward, energy is being lost, often through excessive damping or a contact that is not transmitting force correctly.

A common failure mode is a sudden spike in hourglass energy in a specific part. This usually means the element formulation is too soft for the loading mode (e.g., fully integrated elements are not used where bending dominates) or the mesh is too coarse. Switch to ELFORM = -1 or -2 (fully integrated) for the affected part, or refine the mesh, then re-check glstat.

Putting It Together: A Minimal Frontal Crash Deck

The following snippet shows the essential control cards for a frontal impact run. It is a teaching example, not a validated production deck.

*CONTROL_TERMINATION
$  ENDTIM
   100.0

*CONTROL_TIMESTEP
$  DTINIT   TSSFAC   DT2MS    LCTM     ERODE   MS1ST
   0.0      0.9     -1.0e-3   0        0       0

*CONTROL_ENERGY
$  HGEN     RWEN     SLNTEN   RYLEN
   2        2        2        2

*CONTROL_HOURGLASS
$  IHQ      QH
   4        0.05

*DATABASE_GLSTAT
$  DT
   0.1

*DATABASE_MATSUM
$  DT
   0.1

With DT2MS = -1.0e-3 (i.e., −1 µs), any element whose stable step is below 1 µs is selectively mass-scaled to that target. IHQ = 4 activates the Flanagan–Belytschko hourglass control with a coefficient of 0.05, which is a reasonable starting point for shell-dominated body structures. HGEN = 2 and RWEN = 2 request hourglass and rigid-wall energy tracking in glstat.

Reading the Results: What to Check and When

After the run, open glstat and plot total energy, kinetic energy, internal energy, and hourglass energy against time. A healthy frontal crash shows KE dropping from its initial value to near zero as the vehicle stops, IE rising correspondingly, and TE staying flat within a few percent. If hourglass energy exceeds 5% of internal energy at any point, investigate the responsible parts. If added mass exceeds 2%, reduce the mass-scaling target or refine the offending mesh. If the time step collapses mid-run, look for element inversion or negative volume — often caused by contact penetration or an overly aggressive material failure criterion.

Explicit Crash Solver: Time Step, Mass Scaling, and Energy Flow Central-Difference Loop aₙ = M⁻¹ (F_ext − F_int − C·vₙ) vₙ₊½ = vₙ₋½ + aₙ · Δt xₙ₊₁ = xₙ + vₙ₊½ · Δt next step Stable Time Step (CFL) c = √(E / ρ) Δt_crit = L_e / c Δt = 0.9 · min(Δt_crit) smallest element controls runtime Mass Scaling Δt ∝ √ρ → scale density DT2MS > 0 : fixed DT2MS < 0 : selective keep added mass < 2% Energy Balance Check (glstat) KE decays to ~0 IE monotonic rise HE < 5% of IE SIE small, watch spikes TE = KE + IE + HE + SIE drift < 5% TE drift > 5% → check hourglass, contact, Δt re-run after fix

The workflow above is iterative: run a short initialization check, inspect the reported time step and added mass, then run the full event and validate energy balance. Only after these checks pass should you interpret occupant or structural results. Skipping them is the most common reason a frontal crash model produces plausible-looking but physically wrong predictions.