Newmark-β Method for Transient Structural Analysis

Transient structural analysis computes how a structure responds dynamically to time-varying loads — earthquakes, impacts, blast pressure waves, machine start-up vibrations. Rather than solving a steady-state problem, we need to march the displacement, velocity, and acceleration forward in time. Newmark-β is an implicit time integration scheme proposed by Nathan Newmark in 1959. "Implicit" means the unknowns at the new time step $n+1$ appear on both sides of the equations — you must solve a linear system at every step. The reward: unconditional stability for the right parameter choice, meaning you can use large time steps without the solution blowing up. This makes it ideal for long-duration seismic events lasting tens of seconds.

Explicit methods (central difference scheme) compute the next step's acceleration directly from current forces — no linear system solve needed. This makes each step extremely cheap. The catch: stability requires the time step to satisfy the CFL condition, $\Delta t \le T_{min}/\pi$, where $T_{min}$ is the smallest natural period in the mesh. For a fine crash mesh with tiny elements, $T_{min}$ might be $10^{-7}$ s, forcing hundreds of thousands of time steps per millisecond. That's actually fine for crash (events last 100–150 ms) but completely impractical for seismic events lasting 20–60 seconds. Implicit Newmark flips this: each step costs more (a linear solve), but you can take time steps 100–1000× larger, making long-duration events tractable.

The multi-degree-of-freedom (MDOF) equations of motion from FEM discretization:

$[M]$ is the assembled mass matrix, $[C]$ is the damping matrix, $[K]$ is the stiffness matrix, and $\{F(t)\}$ is the time-dependent external force vector. For seismic analysis, $\{F(t)\} = -[M]\{r\}\ddot{u}_g(t)$ where $\{r\}$ is the influence vector and $\ddot{u}_g(t)$ is the ground acceleration record. Displacement $\{u\}$, velocity $\{\dot{u}\}$, and acceleration $\{\ddot{u}\}$ must all be marched forward simultaneously in time from known initial conditions.

The Newmark update equations for displacement and velocity at time $t_{n+1} = t_n + \Delta t$:

The parameters $\beta$ and $\gamma$ control how the acceleration is interpolated within the time step. Since $\{\ddot{u}\}_{n+1}$ is unknown, we substitute the displacement equation into the equations of motion evaluated at $t_{n+1}$ and solve for $\{\ddot{u}\}_{n+1}$. This is the "corrector" step that makes the method implicit.

After substituting the Newmark equations into $[M]\{\ddot{u}\}_{n+1} + [C]\{\dot{u}\}_{n+1} + [K]\{u\}_{n+1} = \{F\}_{n+1}$, the system can be rearranged as:

where the effective stiffness matrix is:

and $\{\hat{F}\}_{n+1}$ is an effective load vector that incorporates the known quantities from time step $n$. The key insight: $[\hat{K}]$ only changes when $\Delta t$ changes (or when material nonlinearity updates $[K]$). For linear problems with constant $\Delta t$, $[\hat{K}]$ is factorized once at the start using LU decomposition, and only back-substitutions are needed at each subsequent step — making the incremental cost per step much lower than it might appear.

Here are the key parameter combinations and their properties:

The standard choice for linear problems is $\beta = 1/4$, $\gamma = 1/2$ (constant average acceleration). It's unconditionally stable and second-order accurate — but introduces zero numerical damping, so high-frequency modes can sustain artificial oscillations indefinitely. For nonlinear analyses or when high-frequency noise is a concern, HHT-α is far preferable.

The Hilber-Hughes-Taylor (HHT) method modifies the equilibrium equation used at each step by taking a weighted average between time steps $n$ and $n+1$ for the internal forces:

The parameters relate to Newmark's as:

$\alpha = 0$ recovers the standard constant average acceleration Newmark. $\alpha = -1/3$ gives maximum high-frequency numerical damping while maintaining second-order accuracy — something impossible with plain Newmark. The practical effect: high-frequency spurious oscillations are damped out while the structural response in the frequency range of interest remains accurate. Abaqus defaults to $\alpha = -0.05$; Ansys Mechanical to $\alpha = -0.1$.

Related but not identical. The Generalized-α method (Chung and Hulbert, 1993) introduces separate interpolation parameters for inertial, damping, and stiffness forces, giving more freedom to optimize numerical properties. The method is parameterized by $\rho_\infty \in [0, 1]$ — the spectral radius at infinite frequency:

$\rho_\infty = 1$ gives zero numerical damping (equivalent to Newmark). $\rho_\infty = 0$ gives maximum damping. The user simply specifies $\rho_\infty$ and the solver computes all other parameters. Generalized-α is the de facto default in COMSOL and many modern research codes because $\rho_\infty$ has a direct physical interpretation and is easy to tune.

Rayleigh damping constructs the damping matrix as a linear combination: $[C] = \alpha_M [M] + \beta_K [K]$. The coefficients that give damping ratio $h$ at two frequencies $\omega_1$ and $\omega_2$ are:

The damping ratio at any frequency $\omega$ is $h(\omega) = \frac{1}{2}\!\left(\frac{\alpha_M}{\omega} + \beta_K\omega\right)$, which is exactly $h$ at $\omega_1$ and $\omega_2$, higher between them, and lower outside. Strategy: choose $\omega_1$ as the fundamental frequency (or lowest dominant frequency) and $\omega_2$ as the highest frequency of interest. For seismic analysis with ground motion dominated by 1–10 Hz, set $\omega_1 = 2\pi \times 1$ rad/s and $\omega_2 = 2\pi \times 10$ rad/s. Caution: making $\beta_K$ large introduces artificial stiffness proportional damping that can cause numerical issues in nonlinear analyses — keep it as small as feasible.

Here's the complete workflow for seismic time-history analysis:

1. Select ground motion records: Use 3–7 scaled and matched acceleration time histories consistent with the target design spectrum (code-specified or site-specific PSHA). Scale records to the target spectral acceleration at the fundamental period.
2. Set time step: $\Delta t \le T_1/10$ where $T_1$ is the fundamental period. For a building with $T_1 = 1.5$ s, use $\Delta t \le 0.15$ s. If the record sampling interval is finer (e.g., 0.01 s), use that as the time step.
3. Damping setup: For RC structures, target 5% damping per code. For steel, 2–3%. Set Rayleigh coefficients covering the frequency range of the first three or four contributing modes.
4. Apply input motion: Option 1 — apply as support point acceleration (inertia loading $-[M]\{r\}\ddot{u}_g$). Option 2 — impose prescribed base accelerations with large penalty stiffness. Support point acceleration is numerically cleaner.
5. Post-process: Extract inter-story drift ratios (code limit typically 2% for life safety), floor accelerations for equipment qualification, and base shear time histories. For 7 records, use the average of results; for 3 records, use the envelope.
6. Validation: Run an equivalent linear modal superposition (response spectrum) analysis and check that the Newmark peak base shear is within ±15% for the dominant mode contributions.

3000 rpm = 50 Hz. The accuracy rule for the Newmark method is $\Delta t \le T_{response}/10$. For 50 Hz response, $T = 0.02$ s, so $\Delta t \le 0.002$ s. But if you're also interested in harmonics up to the 5th order (250 Hz), you need $\Delta t \le 0.0004$ s. For steady-state forced response at a known frequency, consider using the frequency domain (harmonic analysis) or modal superposition instead — these are far more efficient than time-stepping for steady-state problems. Time-domain Newmark is best when the transient startup behavior or the response to a complex non-periodic excitation must be captured.

Here's a comparison of major tools:

*DYNAMIC (in Abaqus/Standard) is the implicit HHT-α integration — unconditionally stable, each step requires solving a linear system. Use this for seismic, impact on compliant structures, machine vibration, and other events where you want large time steps and automatic time stepping. *DYNAMIC, EXPLICIT (Abaqus/Explicit) is the central difference explicit scheme — conditionally stable, no linear solve, tiny time steps, ideal for crash, high-speed impact, explosive loading, and sheet metal forming. A critical detail: Abaqus/Standard and Abaqus/Explicit use completely separate solver kernels. You can run a co-simulation (Abaqus Cosimulation Interface), but for most problems you choose one or the other based on the physics and duration of the event.

Several important extensions have emerged since 1959:

• Generalized-α Method (Chung & Hulbert, 1993): Already covered above — the practical de facto standard in modern solvers, using $\rho_\infty$ as the single user parameter.
• GSSSS Framework (Tamma, Zhou, Sha): A comprehensive parametric framework that includes Newmark, HHT, and Generalized-α as special cases, offering maximum flexibility for algorithm design.
• IMEX (Implicit-Explicit) Splitting: Treat stiff DOFs (e.g., contact forces, fast-wave structural components) implicitly and flexible DOFs explicitly. Enables coupled FEM-DEM simulations and reduces the cost for systems with mixed stiffness scales.
• Energy-Conserving Algorithms (Simo & Tarnow, 1992): For highly nonlinear elastodynamics (flexible multibody dynamics, thin shells), standard Newmark can permit energy drift. Energy-momentum conserving algorithms ensure discrete energy and momentum balance, critical for long-time integration of conservative mechanical systems.
• Isogeometric Transient Analysis: Combining Newmark time integration with isogeometric spatial discretization (NURBS-based FEM) enables higher-order spatial accuracy and direct CAD-to-analysis workflows for structural dynamics.

Here are the characteristic failure modes and their remedies: