New Asynchronous Time Integrator Preserves Second-Order Accuracy in Multi-Scale Dynamic Simulations

Modern structural mechanics and multiphysics engineering simulations frequently confront systems whose physical subcomponents evolve across vastly disparate spatial resolutions and characteristic temporal rates. When modeling intricate phenomena such as wave propagation through soil-structure interfaces, high-velocity impact responses, or transient fluid-structure interactions within nuclear reactor vessels, resolving the entire computational domain with a uniformly tiny time increment imposes an immense and often prohibitive computational burden. To circumvent this inefficiency, computational mechanicians rely on domain decomposition techniques paired with Heterogeneous Asynchronous Time Integrators. These algorithmic frameworks divide the primary domain into non-overlapping subdomains, granting analysts the freedom to assign distinct spatial discretizations, numerical integration schemes, and local time steps tailored specifically to the physical behavior of each distinct zone. However, retaining both algorithmic stability and high temporal accuracy across the boundaries where these disparate time domains intersect has remained a long-standing numerical challenge.

Derived from the dual Schur formulation of the Finite Element Tearing and Interconnecting method, the GC algorithm employs Lagrange multipliers to glue non-overlapping subdomains by enforcing kinematic velocity continuity along the shared interface at the fine, or micro, time scale. Prescribing continuity at the micro temporal resolution provided exceptional algorithmic flexibility, allowing engineers to handle localized nonlinearities explicitly without demanding costly tangent matrix assemblies, Newton-Raphson equilibrium iterations, or convergence checks at every micro interval. Despite its robustness and unconditional energy-method stability, the classic GC formulation suffers from a critical theoretical drawback: whenever differing time steps are assigned to adjacent subdomains, the global convergence rate drops from second-order to first-order accuracy, accompanied by spurious numerical energy dissipation along the interface that accumulates continuously as simulation duration grows.

To overcome this convergence degradation without sacrificing the operational advantages of micro-scale kinematic enforcement, researchers Yvan Le Nôtre, Michael Brun, Clement Grenat, and Anthony Gravouil introduced an advanced multi-time-step formulation designated as the BLG method. Detailed in a study published in Advanced Modeling and Simulation in Engineering Sciences, the BLG framework revises the interface coupling conditions across each macro time step. Rather than enforcing velocity continuity across all consecutive micro intervals, the BLG method prescribes velocity continuity strictly for intermediate micro steps while deliberately switching to an acceleration continuity constraint at the final micro time step that coincides with the macro step boundary. This hybrid continuity strategy injects the elevated accuracy characteristic of acceleration constraints precisely at the macro milestone while retaining the proven stability and non-iterative simplicity that made micro-scale velocity gluing so attractive for nonlinear transient dynamics.

The foundational mathematical formulation begins with the semi-discrete Lagrangian describing linear elastic bodies under small displacement assumptions, where continuous balance equations are transformed into nodal displacements, velocities, and consistent mass and stiffness matrices through finite element space discretization. When partitioning the global body into two non-overlapping subdomains, designated subdomain A and subdomain B, the variational action integral is augmented with Lagrange multipliers that enforce kinematic compatibility along the interface. The resulting coupled system is segregated into an unconstrained free problem and an interface-correcting link problem. The free problem computes provisional structural accelerations and velocities driven strictly by external applied loads and internal elastic restoring forces, whereas the link problem computes corrective kinematic adjustments produced by the interface Lagrange multiplier forces, which act to eliminate boundary mismatches.

The temporal discretization of the BLG framework establishes an integer time step ratio denoted as m, defined by the macro time step of subdomain A divided by the micro time step of subdomain B. Within each macro interval, intermediate micro steps running from index one through m minus one enforce velocity continuity between the two subdomains through an interface operator matrix derived from the effective dynamic mass matrices and Newmark integration parameters. However, when the micro discretization reaches the final sub-interval where the local counter equals m, the coupling condition transitions directly into an acceleration continuity equation.

To rigorously verify the stability and precision of this modified interface coupling, the investigators conducted a comprehensive spectral analysis based on a Single Degree of Freedom split oscillator. By formulating a seven-component discrete state vector encompassing nodal displacements, velocities, free velocities, and accelerations across both subdomains, the authors derived the global amplification matrix governing state transitions over a macro interval. Assessing the spectral radius of this amplification matrix across varying values of the reduced micro angular frequency enabled an exact mathematical evaluation of algorithmic stability limits without relying on simplified continuum analogies.

The spectral investigations examined three distinct integration configurations: coupling two implicit Constant Average Acceleration schemes, coupling an implicit Constant Average Acceleration scheme with an explicit Central Difference scheme, and coupling two explicit Central Difference schemes across different time step ratios. In the mixed implicit-explicit coupling configuration, which represents the most computationally advantageous setup for transient impact and wave interaction problems, the explicit subdomain must satisfy the standard Courant-Friedrichs-Lewy stability criterion. The spectral radius evaluations revealed that while an alternative pure acceleration scheme drastically curtailed the stable operating zone, the BLG method remained remarkably stable. At a time step ratio of twenty, the BLG stability boundary settled at a reduced frequency of 1.93, positioned just below the theoretical Courant limit of two. As the temporal ratio expanded to one hundred, the critical stability threshold converged smoothly toward 1.99, demonstrating that the BLG formulation asymptotically recovers the full stability domain of uncoupled explicit solvers.

Beyond stability, analyzing the amplification matrix yielded crucial mathematical insights into local truncation error, period elongation, and algorithmic damping ratios. Evaluating the difference between exact analytical state transitions and discrete matrix operations confirmed that the BLG formulation dramatically suppresses high-order truncation errors. In accordance with the classical Lax equivalence theorem, an algorithm whose amplification matrix satisfies stability constraints while exhibiting a local truncation error of order k plus one achieves a global convergence rate of order k. Consequently, for all practical multi-scale ratios greater than two, the BLG method delivers genuine second-order convergence, gaining an entire order of accuracy over the progenitor GC algorithm.

The analytical evaluation of algorithmic damping and period distortion further corroborated the technical superiority of the hybrid interface approach. Expressing the principal complex conjugate eigenvalues in exponential form allowed direct quantification of the numerical damping ratio and period elongation as functions of the macro angular frequency. While the GC method exhibits an algorithmic damping ratio that scales linearly with the reduced frequency, the BLG method demonstrates an algorithmic damping order proportional to the cube of the reduced frequency. This cubic relationship indicates that low-frequency operational modes experience negligible artificial attenuation, preventing unphysical energy loss during prolonged structural vibrations. Concurrently, both the GC and BLG methods maintain an identical, second-order accurate period elongation profile, ensuring that wave propagation phase velocities remain undistorted by the introduction of the acceleration-switching interface boundary.

Direct numerical simulations of the free vibration response of the split oscillator validated these spectral derivations against analytical harmonic benchmarks. Calculating the relative errors in displacement, velocity, and acceleration at the completion of a free vibration interval encompassing three full periods proved that the BLG method maintains an unbroken slope of two in logarithmic error plots across a broad spectrum of macro time steps. In stark contrast, identical simulations conducted with the GC scheme collapsed to a slope of one, affirming the degradation of the older algorithm to first-order accuracy. Furthermore, computing discrete interface energy balance norms highlighted a dramatic curtailment of spurious interface dissipation under the BLG formulation. When tested with a tiny micro time step of one-tenth of a microsecond over an extended simulation duration spanning thirty vibration periods, the normalized relative interface energy dissipation reached negative 0.96 under the GC method, signifying massive energy drainage, whereas the BLG method maintained a modest relative dissipation of only negative 0.07.

To confirm that these outstanding operational characteristics extend to continuum engineering mechanics, the research team implemented the BLG formulation in a two-dimensional finite element model of a clamped-free structural beam. Measuring two and a half meters in length with a twenty-five centimeter square cross-section, the beam was composed of concrete characterized by an elastic modulus of twenty gigapascals, a Poisson ratio of 0.2, and a density of twenty-five hundred kilograms per cubic meter. Discretized under plane stress assumptions using four-node linear quadrilateral elements, the structure was divided at its mid-span into two distinct subdomains. The left partition, clamped securely at its base, was integrated using a fine micro time step of one microsecond, while the right partition, carrying a rapid transient vertical tip load of negative twenty-seven hundred newtons, operated with a coarse macro time step of ten microseconds, establishing a localized time step ratio of ten.

The resulting dynamic tip displacement profiles computed via the BLG multi-time-step explicit coupling closely matched the trajectory of an unpartitioned reference model integrated with a uniform time step of one-tenth of a microsecond. The maximum dynamic deflection of approximately four millimeters accurately mirrored twice the static deflection, oscillating with a fundamental period of 53.4 milliseconds in exact alignment with structural modal predictions. Successive grid convergence evaluations performed at the termination time of 0.125 seconds revealed that as the coarse time step systematically decreased, the relative displacement error of the BLG framework decreased along a clean second-order asymptotic convergence trajectory, whereas the GC method remained trapped in first-order convergence. By successfully combining micro-scale non-iterative explicit simplicity with rigorous second-order temporal precision and minimal interface damping, the BLG method establishes a robust computational foundation for simulating complex nonlinear impacts, earthquake soil-structure interactions, and dynamic multiphysics systems.

Subject of Research: Heterogeneous asynchronous time integration algorithms for multi-scale structural dynamics and multiphysics modeling

Article Title: A new micro HATI second order preserving for transient dynamics

Article References: Le Nôtre, Y., Brun, M., Grenat, C., & Gravouil, A. (2025). A new micro HATI second order preserving for transient dynamics. Advanced Modeling and Simulation in Engineering Sciences, 12(1), Article 9. https://doi.org/10.1186/s40323-025-00286-6

Image Credits: AI Generated

DOI: 10.1186/s40323-025-00286-6

Keywords: heterogeneous asynchronous time integrator, BLG method, structural dynamics, time integration, finite element method, domain decomposition, Newmark schemes, central difference scheme, constant average acceleration, amplification matrix, spectral stability, local truncation error, numerical damping, interface energy, spurious dissipation, split oscillator, Lagrange multipliers