An integrated simulation method for coupled dynamic systems

The paper addresses the challenge of accurately and efficiently simulating large-scale, complex structural systems that involve multiple materials, components, loading scenarios, and coupled physical fields (e.g., soil–structure and fluid–structure interactions). Conventional experimental and single-code numerical approaches struggle with scalability and comprehensive modeling. Integrated (partitioned) simulation methods enable decomposition of a large problem into smaller dynamic substructures that can be handled by existing tools. However, monolithic coupling is difficult to implement and partitioned methods often require intrusive code modifications and iterative rollbacks, which many proprietary FE codes do not support. The research objective is to develop a non-iterative staggered partitioned method that treats each substructure’s solver as a black box, exchanging only interface displacements and forces, permitting nonhomogeneous time integrators, and providing a practical, verifiable stability criterion (maximum stable time step) for coupled dynamic analyses. The importance lies in enabling robust, reusable integration of disparate FE programs for complex, large-scale dynamic simulations with controlled stability and accuracy.

Prior work distinguishes monolithic (direct) and partitioned (staggered/iterative) domain-decomposition methods for coupled dynamics (Felippa et al., 2001). Monolithic approaches advance all subsystems simultaneously and naturally enforce interface compatibility/equilibrium but are difficult to implement and reuse, especially across different codes. Partitioned methods integrate subsystems separately, exchanging interface variables. Iterative partitioned schemes (e.g., Dirichlet–Neumann iterations) generally offer better stability/accuracy but require rollback and deep code-level access, often impractical in commercial software (Gravouil & Combescure, 2001; Gravouil et al., 2015; Pan et al., 2005, 2006). Staggered (non-iterative) schemes are easier to implement but only conditionally stable, and performance depends on interface predictors (Giles, 1997; Rugonyi & Bathe, 2001). Partitioned approaches have been applied to hybrid simulation and coupled multiphysics, with sensitivity to mass/stiffness ratios across subdomains (Dettmer & Perić, 2013; Kübler & Schiehlen, 2000; Jahromi et al., 2007). The authors previously developed a generalized component-level partitioned framework (Huang & Kwon, 2018); this paper focuses on system-level coupling of dynamic subsystems and proposes a practical stability assessment method.

Formulation: Starting from the system equation of motion M a + C v + R = F, the system is decomposed into two substructures i and j. A generalized time-integration predictor–corrector framework (Newmark family) is used with coefficients (β1, β2, β3, γ1, γ2). Each substructure’s effective equation at time n+1 is (Ai + Aij) ai,n+1 = Bi,n+1 + Dij,n+1, where Ai encodes Mi, Ci, Ki and integration parameters, Aij represents interface mass from the counterpart substructure, and Dij,n+1 is the effective external force that substructure j imposes on i at the interface, assembled from predicted interface displacements/velocities/forces of j. Coupling strategies: Two staggered approaches are proposed to render Dij,n+1 explicit: - Method-1: Enforce interface compatibility (continuity of displacement, velocity, acceleration), expand interface variables of Sub-j using corrector relations, then replace unknown interface quantities of Sub-i by their predictors to obtain an explicit Dij,n+1 and an adjusted Aij. - Method-2: Keep Aij as the interface mass only; explicitly determine Dij,n+1 by replacing Sub-j’s displacement/velocity at n+1 with their predictors from the previous step. This yields the same explicit form of Dij,n+1 as Method-1 but with simpler Aij. Staggered algorithm (per time step): (1) choose Δt; (2) primary substructure j provides interface mass Aij to secondary i (predefined if constant); (3) initialize n=0; (4) j predicts Dij and imposes it on i; (5) solve Sub-i for ai,n+1; (6) update di,n+1, vi,n+1; (7) impose ai,n+1 to j’s interface; (8) solve Sub-j with its own integrator; (9) increment n and repeat from step 4. Only one prediction and one substitution per step, no rollback. Parametric study: The example 3-DOF mass–spring–damper system is used to compare Method-1 and Method-2 across three integrator pairings: NAA-to-NAA (both Newmark average acceleration), ENM-to-ENM (both explicit Newmark), and NAA-to-ENM (mixed), under varying primary/secondary mass and stiffness ratios. Performance is assessed versus a standalone reference model using the same integrator as the coupled secondary. Stability analysis (Newmark family, α=0): An energy method (Richtmyer & Morton) is used since modal analysis does not apply to mixed integrators. By transforming to incremental form and deriving energy inequalities, stability reduces to conditions on positive definiteness of an equivalent matrix Ψ that depends on Δt and interface contributions. Closed-form stability bounds are obtained for NAA-to-NAA, ENM-to-ENM, and NAA-to-ENM cases (e.g., Δt upper bounds from positivity of modified mass-like blocks). Implementation strategy: To minimize code changes, the primary substructure is realized as thin, linear-elastic, undamped interface layers. Then Aij is constant interface mass, and Dij,n+1 can be computed directly from the predefined interface stiffness and predicted interface DOFs without extracting internal damping/restoring force details from black-box FE codes. Secondary substructures (e.g., soil, dam, building segments) are solved in existing FE programs; only interface DOFs’ displacements/accelerations/reaction forces are exchanged. Generalization to α-modified schemes: The method extends to α-modified time integration (Hilber–Hughes–Taylor family) to introduce numerical damping. The effective equations and Dij,n+1 are reformulated with (1+α) weightings. Because the energy method is limited to α=0, a practical stability counterpart procedure is introduced: derive conservative critical Δt by forming stability counterparts of primary and secondary substructures that include equivalent interface mass, damping, and stiffness contributions; then apply known stability criteria of the chosen integrators to these counterparts. The overall Δt must not exceed the minimum critical Δt across substructures. Verification examples: (1) Nonlinear five-story shear building decomposed into three substructures with mixed integrators (NAA and alpha-operator splitting). Stability counterpart yields Δtcr ≈ 0.02 s; integrated results at Δt ≤ 0.02 s match standalone solutions; larger Δt causes instability. (2) Large-scale soil–structure interaction of a concrete gravity dam with two soil domains and the dam as secondary substructures (all modeled in ABAQUS) and thin soil interface layers as primary. Using Δt=0.001 s (below predicted critical), the integrated model reproduces the standalone model’s displacement fields and time histories under an obliquely incident Ricker wavelet.

- Method-2 (explicit prediction of Dij with unchanged Aij) is generally more accurate than Method-1 but less stable. Method-1 is more robust (larger stable Δt) but exhibits larger errors, especially with large primary-to-secondary stiffness ratios. - NAA-to-NAA (both unconditionally stable integrators): Despite individual unconditional stability, the coupled staggered scheme is only conditionally stable. For the benchmark system with mass and stiffness ratios (1,1), analysis indicates stability requires Δt ≤ about 0.038 s; instability is observed at Δt=0.04 s. Method-2 becomes unstable at larger Δt whereas Method-1 remains stable but less accurate. NRMSE examples (Table 3): for (500,1), Method-1 errors increase markedly with Δt (e.g., 13.7%, 60.6%, 511.5%), while Method-2 is accurate at small Δt (e.g., 0.04% at Δt1) but unstable at larger Δt. - ENM-to-ENM (both explicit Newmark): Both methods are conditionally stable with more stringent limits than standalone ENM (Δt ≤ T/π). Additional interface damping term tightens stability (Ψ = M − 0.5 Δt^2 Cji − 0.25 Δt^2 K). Cases show instability for coupled models even when the standalone limit is satisfied (e.g., unstable at Δt=0.001 s for coupled cases where standalone permits it). Method-1 is more stable but less accurate than Method-2. - NAA-to-ENM (mixed): Stability is governed by positive definiteness of a block-diagonal Ψ; for representative parameters, stability requires Δt ≤ about 0.0185 s. Again, Method-2 offers better accuracy but smaller stability margins. - Practical stability assessment: A general stability counterpart method yields conservative Δtcr for arbitrary integrator combinations, including α-modified schemes. If both substructures’ integrators are unconditionally stable, the integrated scheme’s stability reduces to the positive definiteness of the equivalent interface mass matrix; thus stability depends only on the (linear-elastic) primary interface layers. - Implementation feasibility: By modeling thin linear-elastic interface layers as the primary substructure, Aij and Dij can be computed without intrusive access to secondary solvers, enabling black-box coupling across different FE codes. - Verification: (1) Five-story nonlinear building: stability counterpart predicts Δtcr ≈ 0.02 s; integrated results at Δt=0.01–0.02 s match standalone; Δt=0.021 s leads to instability. (2) Gravity dam–soil system: integrated solution at Δt=0.001 s closely matches standalone spatial deformation fields and time histories along sections X–X and Y–Y under an oblique SV-wave input.

The study demonstrates that a non-iterative staggered partitioned approach can effectively couple dynamic substructures analyzed with disparate FE programs by exchanging only interface DOFs’ displacements/forces, treating each code as a black box. The two proposed formulations reveal a trade-off between stability and accuracy: Method-2, which keeps the interface mass unchanged and explicitly predicts interface forces, attains higher accuracy but reduced stability margins, while Method-1 is more stable but less accurate. Stability is not guaranteed by using individually unconditionally stable integrators; the coupling introduces additional constraints that render the overall scheme conditionally stable. The energy-based stability analysis for Newmark schemes and the generalized stability-counterpart procedure (also covering α-modified methods) provide practical means to compute conservative maximum stable time steps without intrusive code modification or modal decomposition. Verification on a nonlinear building and a large-scale dam–soil interaction problem confirms that, with Δt selected via the proposed stability procedure, the integrated solutions closely match standalone results. The approach facilitates code reuse and mixed time integrators, enabling integrated simulations across different physics and scales while maintaining control over stability and accuracy.

The paper introduces a staggered partitioned method for integrating coupled dynamic substructures using only interface displacements and forces, enabling straightforward coupling of different FE programs as black boxes and supporting nonhomogeneous time integrators. Two variants were studied; Method-2 was selected due to superior accuracy, acknowledging reduced stability margins. A practical stability framework was established: (1) an energy-based analysis for Newmark schemes (α=0) yielding explicit Δt bounds, and (2) a general stability-counterpart procedure applicable to α-modified schemes and mixed integrators. Numerical studies show the method is conditionally stable, with stability dependent on time step, mass/stiffness ratios, and interface layer properties. Implementation via thin linear-elastic interface layers allows non-intrusive computation of interface mass and forces. Verification on nonlinear building and large-scale dam–soil interaction problems confirmed accuracy when Δt respects the derived criteria. Future work includes extending the method to allow nonlinear interface layers to relax the layer-size limitation and broaden stability margins, and further demonstrations on additional large-scale, multi-platform simulations (e.g., nuclear containment SSI, multi-scale urban seismic analysis).

- Conditional stability: Even with individually unconditionally stable integrators, the coupled staggered scheme is only conditionally stable and requires careful selection of Δt. - Interface layer assumption: The primary substructure is assumed linear-elastic and undamped; larger layer sizes can improve stability but may introduce modeling bias and are limited by this linear assumption. Nonlinear interface behavior is not yet supported in the presented formulation. - Sensitivity to substructure properties: Stability and accuracy depend on primary-to-secondary mass and stiffness ratios; highly disparate ratios can degrade performance or reduce stable Δt. - Non-iterative coupling: Without corrective iterations, exact interface compatibility/equilibrium at each step is not enforced, leading to approximation errors compared to monolithic or iteratively coupled schemes. - ENM restrictions: When explicit methods are used, stability constraints are stricter than standalone conditions due to additional interface terms. - Publication lacks a precise publication date; examples are limited to two detailed cases (with others referenced but not presented).