4. Implicit time integration fundamentals

Explicit integration schemes are conditionally stable: the critical time step is governed by the highest natural frequency of the mesh, so stiff structures or fine meshes force impractically small increments. Implicit schemes remove that restriction by solving a system of equations at each step, which couples the unknown accelerations, velocities, and displacements. This chapter develops the three families that dominate structural dynamics codes — Newmark-beta, HHT-alpha, and the generalized-alpha method — and explains how numerical damping and unconditional stability are engineered into them.

The semi-discrete equation of motion

After spatial discretization, the equations of motion for a linear system are

M a(t) + C v(t) + K d(t) = f_ext(t)

where M, C, and K are the mass, damping, and stiffness matrices; d, v, and a are the nodal displacement, velocity, and acceleration vectors; and f_ext is the external load. Time integration advances the solution from step n to n+1. The essential difference between explicit and implicit methods lies in when the equilibrium equation is enforced. Explicit methods evaluate the residual at time t_n using known quantities. Implicit methods enforce equilibrium at t_{n+1}, where d_{n+1} is unknown, producing an algebraic system that must be solved.

Newmark-beta: the parent family

The Newmark-beta method approximates displacement and velocity with two parameters, beta and gamma:

d_{n+1} = d_n + dt*v_n + dt^2*[(0.5 - beta)*a_n + beta*a_{n+1}]
v_{n+1} = v_n + dt*[(1 - gamma)*a_n + gamma*a_{n+1}]

Substituting these into the equilibrium at t_{n+1} and collecting the unknown d_{n+1} yields the effective stiffness and effective load:

K_eff = K + (gamma/(beta*dt))*C + (1/(beta*dt^2))*M
f_eff = f_ext_{n+1} + M*[ (1/(beta*dt^2))*d_n + (1/(beta*dt))*v_n + (1/(2*beta) - 1)*a_n ]
      + C*[ (gamma/(beta*dt))*d_n + (gamma/beta - 1)*v_n + dt*(gamma/(2*beta) - 1)*a_n ]

Solving K_eff * d_{n+1} = f_eff gives the new displacement, after which the velocity and acceleration are recovered from the update formulas. The two parameters control accuracy and stability:

• Average acceleration (beta = 0.25, gamma = 0.5): second-order accurate, unconditionally stable, zero numerical damping. This is the default in most structural codes.
• Linear acceleration (beta = 1/6, gamma = 0.5): second-order accurate but only conditionally stable, with a critical step roughly dt_crit = 0.551*T_min.
• Fox-Goodwin (beta = 1/12, gamma = 0.5): also conditionally stable.

Unconditional stability for the Newmark family requires gamma >= 0.5 and beta >= 0.25*(gamma + 0.5)^2. When gamma > 0.5, the method introduces numerical damping, but it is high-frequency damping that also contaminates low-frequency modes, which is why plain Newmark is rarely used for its damping alone.

HHT-alpha: shifting equilibrium in time

The Hilber-Hughes-Taylor (HHT-alpha) method modifies the equilibrium equation by evaluating stiffness and damping at a weighted time between t_n and t_{n+1}:

M a_{n+1} + (1 + alpha)*C v_{n+1} - alpha*C v_n
          + (1 + alpha)*K d_{n+1} - alpha*K d_n = f_ext(t_{n+1})

The Newmark update formulas are retained, with the parameter relations

beta = (1 - alpha)^2 / 4
gamma = 0.5 - alpha

For alpha in [-1/3, 0], the method is unconditionally stable and second-order accurate. The parameter alpha controls numerical dissipation: alpha = 0 recovers average-acceleration Newmark with no damping, while alpha = -0.1 or -0.2 provides moderate high-frequency dissipation. The key advantage over Newmark with gamma > 0.5 is that HHT-alpha damps the high-frequency modes far more aggressively than the low-frequency modes, preserving the accuracy of the physically important response.

Generalized-alpha: damping control decoupled from accuracy

The generalized-alpha method (Chung and Hulbert) extends HHT-alpha by also weighting the inertial term, introducing a second parameter alpha_m:

M[(1 - alpha_m)*a_{n+1} + alpha_m*a_n]
  + (1 - alpha_f)*C v_{n+1} + alpha_f*C v_n
  + (1 - alpha_f)*K d_{n+1} + alpha_f*K d_n = f_ext(t_{n+1})

The Newmark relations become

beta  = (1 - alpha_m + alpha_f)^2 / 4
gamma = 0.5 - alpha_m + alpha_f

Second-order accuracy and unconditional stability hold when

alpha_m = (2*rho_inf - 1) / (rho_inf + 1)
alpha_f = rho_inf / (rho_inf + 1)

where rho_inf is the spectral radius at the infinite-frequency limit, 0 <= rho_inf <= 1. This single parameter gives direct control over high-frequency dissipation: rho_inf = 1 means no damping (energy preserved), rho_inf = 0 means asymptotic annihilation of the highest modes, and intermediate values like rho_inf = 0.8 or 0.9 are common in production analyses. Because the low-frequency amplification stays close to unity, the physical response is barely affected.

Numerical damping and why it matters

Finite element discretization produces a spectrum of natural frequencies. The lowest modes represent the physical response of interest; the highest modes are numerical artifacts of the mesh, often with periods comparable to the time step. In impact and wave propagation problems, these spurious high-frequency modes can be excited by the loading and then pollute the solution with non-physical oscillations. Numerical damping selectively attenuates these modes without distorting the low-frequency response.

The generalized-alpha method is the preferred choice in modern explicit-implicit hybrid codes because rho_inf provides a single, intuitive knob. A typical workflow for a stiff structural impact problem might use rho_inf = 0.8, which gives about 20% amplitude reduction per step at the highest frequency while keeping the low-frequency amplification within 0.1% of unity. Setting rho_inf = 1.0 is appropriate when energy conservation is critical, such as in long-duration vibration or wave propagation studies where artificial damping would corrupt the physics.

Implementation notes and verification

An implicit step requires solving K_eff * delta_d = residual. For linear problems, K_eff is constant and can be factorized once; for nonlinear problems, it must be reformed and factorized each iteration, which is the dominant cost. The effective stiffness for generalized-alpha is

K_eff = (1 - alpha_f)*K + ((1 - alpha_f)*gamma/(beta*dt))*C
        + ((1 - alpha_m)/(beta*dt^2))*M

To verify an implementation, run a single-degree-of-freedom oscillator with known period T and compare the computed period elongation and amplitude decay against the analytical amplification matrix. For rho_inf = 0.8 and dt/T = 0.1, the period error should be below 0.5% and the amplitude decay per step below 0.1%. A common failure symptom is a solution that grows without bound despite unconditional stability — this usually indicates an error in the beta/gamma relations or a sign mistake in the effective load vector.

The choice among these methods is not merely academic. In a crash simulation with shell elements, the highest frequencies can be two to three orders of magnitude above the frequencies of interest. Without numerical damping, those modes would require an explicit step size thousands of times smaller than the physical time scale. Generalized-alpha with rho_inf = 0.8 allows a step size set by accuracy requirements rather than stability, which is the entire point of implicit integration for impact and structural response.