# A Newmark space-time formulation in structural dynamics

Franz Bamer, Nima Shirafkan, Xiaodan Cao, Abdelbacet Oueslati, Marcus

Stoffel, Gery de Saxce, Bernd Markert

► To cite this version:

Franz Bamer, Nima Shirafkan, Xiaodan Cao, Abdelbacet Oueslati, Marcus Stoffel, et al.. A Newmark space-time formulation in structural dynamics. Computational Mechanics, 2021, -, 10.1007/s00466-021-01989-4 . hal-04508736

HAL Id: hal-04508736

https://hal.science/hal-04508736v1

Submitted on 12 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Newmark space-time formulation in structural dynamics

Franz Bamer<fnref n="1" /> Ⓔ · Nima Shirafkan<fnref n="1" /> Ⓔ · Xiaodan Cao<fnref n="2" /> · Abdelbacet Oueslati<fnref n="2" /> · Marcus Stoffel<fnref n="1" /> Ⓔ · Géry de Saxcé<fnref n="2" /> Ⓔ · Bernd Markert<fnref n="1" /> Ⓔ

Received: 7 May 2020 / Accepted: 5 February 2021 / Published online: 19 March 2021
© The Author(s) 2021

## Abstract

In this contribution, we present a space-time formulation of the Newmark integration scheme for linear damped structures under both harmonic and transient excitations. The incremental set of equations of motion and the Newmark approximations are transformed into their corresponding space-time equivalents. The dynamic system is then represented by one algebraic space-time equation only. This equation is projected into a coupled pair of space-time equations, which is solved via the fixed point algorithm. The solution is iteratively assembled by enrichments, each of which is decomposed by a dyadic product of spatial and temporal enrichment vectors. The evolution of the spatial enrichment vectors is investigated during convergence and interpreted by comparing them to the set of linear modes of vibration. The new method is demonstrated by means of four numerical examples, presenting not only the excellent convergence behavior and the numerical efficiency but also the limits of the proposed approach.

**Keywords** Newmark integration scheme · Proper generalized decomposition · Space-time integration · Modal analysis

## 1 Introduction

Evaluating the response time history of engineering structures subjected to transient excitation has become a key issue in structural dynamics. The wide range of conventional methods includes analytic strategies for linear systems [1] and step-by-step time integration schemes [2,3]. Especially with emphasis on computational efficiency, several approximation schemes and model order reduction strategies have been proposed in the literature [4–7].

As a statistical pattern technique in incremental structural dynamics, the proper orthogonal decomposition is used to find the dominant motion patterns for an optimal low-order description of linear and nonlinear systems [8]. Here, the transformation into low-order subspaces depends on the set of proper orthogonal modes that are obtained from parts of

the structural response, the so called observation time period [9]. The applicability of the proper orthogonal decomposition to transient excitations mainly depends on the choice of this observation time period [10]. If nonlinearities are not captured in this observation time period, they cannot be reproduced in the low-order model. However, if the nonlinear response behavior is captured, a surprisingly small number of proper orthogonal modes is sufficient to accurately describe the dynamic response behavior of the structure [11,12]. Another approach is to divide the structure of interest into various substructures. The main goal is to provide low-order representations of the inside areas of the substructures while keeping the number of degrees of freedom of the boundaries unaltered. This strategy is pursued both in the Craig-Bampton method [13] and in the dual Craig-Bampton method [14]. Although these two methods are mainly applicable to linear systems, a hierarchical formulation of the Craig-Bampton method has lately been presented for vibro-impact problems with linear contact formulations [15,16] and highly nonlinear contact laws [17]. Recently, also a combination of the proper orthogonal decomposition and substructure techniques has been proposed in the literature [18]. However, the success and efficiency of these strategies is mainly based on the chosen modes or deformation patterns that must be a priori evaluated.

<footnote>✉ Franz Bamer
Bamer@iam.rwth-aachen.de</footnote>

<footnote>Nima Shirafkan
Shirafkan@iam.rwth-aachen.de</footnote>

<footnote n="1">¹ Institute of General Mechanics, RWTH Aachen University,
Eilfschornsteinstraße 18, 52062 Aachen, Germany</footnote>

<footnote n="2">² Univ. Lille, CNRS, Centrale Lille, UMR 9013 - LaMcube -
Laboratoire de Mécanique, Multiphysique, Multiéchelle,
F-59000 Lille, France</footnote>
Another type of approach is to apply non-incremental structural analysis. In this context, one methodology is to apply finite elements simultaneously in space and time [19]. Using this strategy, different space-time element shapes have been introduced, investigating the numerical stability [20,21]. Furthermore, the application of space-time finite elements in elastodynamics and structural dynamics has been introduced [22–24]. It has also been proposed that the space-time finite element approach could be highly effective in countering discontinuous dynamic problems like fracture [25]. Furthermore, space-time finite elements have been applied to contact mechanics [26–28], multiscale modeling [29,30], landslide dynamics [31], problems including fluid-structure interaction [32] and modeling of viscoelastic materials [33–35]. Lately, an efficient solution algorithm for space-time finite element problems has been suggested that is based on the summation of Kronecker products of temporal and spatial submatrices [36].

The proper generalized decomposition (PGD) is an another strategy that does not depend on a priori response identification [37–39], as required for proper orthogonal decomposition based methods [11]. It allows one to consider the whole problem simultaneously in space and time, while using conventional space discretization strategies such as the finite element method. In this way, the method uses an evolving summation of functionals, iteratively evaluated at once within the whole time domain [40]. This solution procedure was firstly introduced by Ladevéze [41]. Subsequently, the PGD was introduced as model order reduction method for multi dimensional problems [42,43]. The strategy has then been applied in many fields such as chemistry [44], thermo-mechanics [45], fluid dynamics [46–48], geophysics [49], material science [50], solid and fracture mechanics [51–54], as well as error estimations and optimization strategies in mechanics [55–57]. Especially in structural mechanics, first approaches have been proposed to use the proper generalized decomposition ansatz to solve elastoplastic and viscoplastic structures subjected to quasistatic loading [58,59]. Boucinha et al. gave a first introduction to how to deal with linear undamped dynamic systems using non-incremental formulations [60]. Shirafkan et al. physically interpreted the evolution of spatial and temporal functions of a quasistatic problem, which they refer to as proper generalized decomposition modes, and analyzed these modes during convergence [61]. Additionally, the solution of parametric dynamic structures in the frequency domain has been investigated [62,63]. However, space-time formulations in structural dynamics are far from well understood. In particular, a detailed description of numeric space-time integration approaches including a dynamical interpretation during convergence is, to date, missing from the literature.

In this paper, we present a new space-time solution strategy in structural dynamics. We provide the fundamental basis

of the continuous and discrete space-time decomposition, based on which we present the space-time equivalents of the set of equations of motion and the incremental Newmark equations. This leads to a coupled space-time matrix equation pair. In order to solve the coupled space-time problem, we propose an iterative approach based on the fixed point algorithm and investigate the evolution of spatial and temporal modes after each enrichment. This paper aims to provide the basis for future space-time approaches applicable to nonlinear transient structural dynamics.

In Sect. 2, we introduce the problem of the conventional incremental elastodynamics in a nutshell, including the strong form, the weak form, the space discretization, the matrix set of equations of motion and the Newmark equations. In Sect. 3, we present the new Newmark space-time integration strategy, introducing the continuous space-time formulation, space-time discretization, space-time decomposition of the equations of motion and the Newmark equations as well as the implementation of the new strategy. In Sect. 4, the Newmark space-time algorithm is demonstrated on four numerical examples, showing the accuracy, efficiency and limits of the method. Finally, in Sect. 5, the conclusions are drawn.

## 2 Conventional incremental elastodynamics in a nutshell

### 2.1 Problem formulation

We define the initial boundary value problem starting with the linear momentum balance at the material point in a certain body $\Omega$ in the time domain of interest $[0, T]$ [64]:

$$ \mathrm{div}(\boldsymbol{\sigma}) + \boldsymbol{f}_v = \rho \ddot{\boldsymbol{u}} \quad \text{in } \Omega \times [0, T], \qquad (1) $$

where $\boldsymbol{\sigma}(\mathbf{x}, t)$ is the continuous second order Cauchy stress tensor field, $\boldsymbol{f}_v$ denotes the vector field of the volume forces and $\mathbf{u}(\mathbf{x}, t)$ is the continuous displacement field. The field $\ddot{\boldsymbol{u}}$ denotes the second derivative of $\mathbf{u}(\mathbf{x}, t)$ with respect to time. The scalar field $\rho(\mathbf{x})$ describes the density distribution of the material within the body $\Omega$. Dirichlet boundary conditions are defined on the surface $\Gamma_u$ as:

$$ \mathbf{u} = \bar{\mathbf{u}} \quad \text{on} \quad \Gamma_u \times [0, T], \qquad (2) $$

while Neumann boundary conditions are defined on the surface $\Gamma_\sigma$ as:

$$ \boldsymbol{\sigma} \mathbf{n} = \bar{\boldsymbol{\mathbf{f}}} \quad \text{on} \quad \Gamma_\sigma \times [0, T]. \qquad (3) $$
Furthermore, two necessary conditions hold on the whole boundary $\partial\Omega$ of the body $\Omega$:

$$ \Gamma_u \cup \Gamma_\sigma = \partial\Omega \quad \text{and} \quad \Gamma_u \cap \Gamma_\sigma = \emptyset. \qquad (4) $$

The initial conditions are defined by the initial displacement and initial velocity fields on the whole body at the time instant $t = 0$:

$$ u|_{t=0} = u_0 \quad \text{in} \quad \Omega, \qquad (5) $$

$$ \dot{u}|_{t=0} = v_0 \quad \text{in} \quad \Omega. \qquad (6) $$

Multiplying Eq. (1) by a weight function that obeys the kinematic boundary conditions and integrating over the body $\Omega$, leads to the weak form of the initial boundary value problem [64]. Interpreting the weight functions as virtual displacements and considering a symmetric stress tensor field leads to the principle of virtual work [64]:

$$ \begin{aligned} & \int_{\Omega} \rho \delta u \ddot{u} \, d\Omega + \int_{\Omega} \delta \boldsymbol{\varepsilon} \boldsymbol{\sigma} \, d\Omega \\ & = \int_{\Omega} \delta u \, \bar{f}_v \, d\Omega + \int_{\Gamma_{\sigma}} \delta u \, \bar{f} \, d\Gamma_{\sigma}. \end{aligned} \qquad (7) $$

## 2.2 Space discretization

The body $\Omega$ is discretized by finite elements in this paper. Using interpolation functions $N_j^{(e)}(\mathbf{x})$ for every element $e$, the nodal element displacement vector $\mathbf{u}^{(e)}$ and its second derivative with respect to time as well as the virtual nodal element displacement vector are introduced [65]:

$$ u^{(e)} = \sum_{j=1}^{n^{(e)}} N_j^{(e)} \mathbf{u}_j^{(e)} = \mathbf{N}^{(e)} \mathbf{u}^{(e)}, \qquad (8) $$

$$ \ddot{u}^{(e)} = \sum_{j=1}^{n^{(e)}} N_j^{(e)} \dot{\mathbf{u}}_j^{(e)} = \mathbf{N}^{(e)} \ddot{u}^{(e)}, \qquad (9) $$

$$ \delta u^{(e)} = \sum_{j=1}^{n^{(e)}} N_j^{(e)} \delta \mathbf{u}_j^{(e)(t)} = \mathbf{N}^{(e)} \delta \mathbf{u}^{(e)}. \qquad (10) $$

Equivalently, the strain field $\boldsymbol{\varepsilon}$ and the virtual strain field $\delta\boldsymbol{\varepsilon}$ are expressed by the nodal element displacements and the strain-displacement relationship $B_j^{(e)}$ [65]:

$$ \boldsymbol{\varepsilon}^{(e)} = \sum_{j=1}^{n^{(e)}} B_j^{(e)} \mathbf{u}_j^{(e)} = \mathbf{B}^{(e)} \mathbf{u}^{(e)} \qquad (11) $$

$$ \delta \boldsymbol{\varepsilon}^{(e)} = \sum_{j=1}^{n^{(e)}} B_j^{(e)} \delta \mathbf{u}_j^{(e)} = \mathbf{B}^{(e)} \delta \mathbf{u}^{(e)} \qquad (12) $$

Inserting relations (8)–(12) into equation (7) and summing up all element volumes leads to the discrete version of the principle of virtual work. The assembling procedure over all elements $\mathbf{A}$ leads from the local element node coordinates $\mathbf{u}^{(e)}$ to the global nodal coordinate $\mathbf{u}$. The spatially discretized version of equation (7) is written in the global nodal coordinate as:

$$ \begin{aligned} & \delta \mathbf{u} \left\{ \underbrace{\left[ \sum_{j=1}^{n^{(e)}} \int_{\Omega^{(e)}} \mathbf{N}^{(e)}{}^T \rho \mathbf{N}^{(e)} \, d\Omega^{(e)} \right]}_{\mathbf{M}^{(e)}} \ddot{\mathbf{u}} \\ & + \underbrace{\left[ \sum_{j=1}^{n^{(e)}} \int_{\Omega^{(e)}} \mathbf{B}^{(e)}{}^T \boldsymbol{\sigma} \, d\Omega^{(e)} \right]}_{\mathbf{f}_{\text{int}}^{(e)}} \\ & + \underbrace{\left[ \sum_{j=1}^{n^{(e)}} \int_{\Gamma_{\sigma}^{(e)}} \mathbf{N}^{(e)}{}^T \bar{\tau} \, d\Gamma_{\sigma}^{(e)} \right]}_{\mathbf{f}_{\text{ext}}^{(e)}} \\ & + \underbrace{\left[ \sum_{j=1}^{n^{(e)}} \int_{\Omega^{(e)}} \mathbf{N}^{(e)}{}^T \mathbf{f}_v^{(e)} \, d\Omega^{(e)} \right]}_{\mathbf{f}_v^{(e)}} \right\} = 0. \end{aligned} \qquad (13) $$

The assembly procedure leads from the element mass matrix $\mathbf{M}^{(e)}$ to the global mass matrix $\mathbf{M} \in \mathbb{R}^{n_s \times n_s}$, from the internal restoring element force $\mathbf{f}_{\text{int}}^{(e)}$ to the global internal force vector $\mathbf{f}_{\text{int}} \in \mathbb{R}^{n_s}$, from the element volume force $\mathbf{f}_v^{(e)}$ to the global volume force vector $\mathbf{f}_v \in \mathbb{R}^{n_s}$ and from the external element force vector $\mathbf{f}_{\text{ext}}^{(e)}$ to the external global force vector $\mathbf{f}_{\text{ext}} \in \mathbb{R}^{n_s}$. In this paper, we present the new strategy applying a linear elastic material law. Thus, the internal restoring element force vector is described by the relationship:

$$ \mathbf{f}_{\text{int}}^{(e)} = \left[ \int_{\Omega^{(e)}} \mathbf{B}^{(e)}{}^T \bar{\mathbf{C}} \mathbf{B}^{(e)} \, d\Omega^{(e)} \right] \mathbf{u}^{(e)} = \mathbf{K}^{(e)} \mathbf{u}^{(e)}, \qquad (14) $$

where $\bar{\mathbf{C}}$ is the elasticity matrix following Voigt's notation and $\mathbf{K}^{(e)}$ is the element stiffness matrix. Assembly of all element stiffness matrices leads then to the global stiffness matrix:

$$ \mathbf{K} = \sum_{j=1}^{n^{(e)}} \mathbf{K}^{(e)}, \quad \mathbf{K} \in \mathbb{R}^{n_s \times n_s}. \qquad (15) $$

## 2.3 The undamped set of equations of motion and modal analysis

Considering Eqs. (13) and (15) leads to the matrix set of equations of motion:

$$ \mathbf{M}\ddot{\mathbf{u}}(t) + \mathbf{K}\mathbf{u}(t) = \mathbf{f}(t), \qquad (16) $$
where we combined the excitation on the right hand side by adding up the external and the volume force vector $\mathbf{f}(t) = \mathbf{f}_v + \mathbf{f}_{\text{ext}}$. In order to obtain the linear modes of vibration, we are looking at the homogeneous problem only. Thus, setting the force at the right hand side of Eq. (16) to the zero vector and defining the solution ansatz function $\mathbf{u}(t) = \mathbf{\Phi}e^{i\omega t}$, leads to the eigenvalue problem:

$$[\mathbf{K} - \omega_i^2 \mathbf{M}] \mathbf{\Phi}^{(i)} = \mathbf{0}. \qquad (17)$$

As the matrix $\mathbf{M}^{-1}\mathbf{K}$ is positive definite for the problem in this paper, one obtains $n_s$ real and positive eigenvalues (squared eigenfrequencies) $\omega_i^2$ ($i = 1...n_s$). Magnitude-wise arrangement of the eigenvalues from low to high values leads to the corresponding sorted set of the modes of vibration $\mathbf{\Phi}^{(i)}$ ($i = 1...n_s$) [1].

## 2.4 Structural damping

In order to model dissipative effects, linear viscous damping is applied in this paper. A structural damping matrix $\mathbf{C}$ is included, using a linear combination of stiffness and mass matrix:

$$\mathbf{C} = a_0 \mathbf{M} + a_1 \mathbf{K}, \quad \mathbf{C} \in \mathbb{R}^{n_s \times n_s}, \qquad (18)$$

where $a_0$ and $a_1$ are the Rayleigh damping coefficients to be evaluated as elaborately discussed in [1].

The linear viscously damped matrix set of equations of motion is then written as [1]:

$$\mathbf{M}\ddot{\mathbf{u}}(t) + \mathbf{C}\dot{\mathbf{u}}(t) + \mathbf{K}\mathbf{u}(t) = \mathbf{f}(t). \qquad (19)$$

## 2.5 Time discretization using the Newmark scheme

Firstly, we assume that Eq. (19), i.e. dynamic equilibrium, is fulfilled at fixed time instants $t_i$, $i = 1...n_t$. Every time interval $\Delta t = t_i - t_{i-1}$ is chosen to be constant in this paper. The set of equations of motion at fixed time instants $t_i$ reads:

$$\mathbf{M}\ddot{\mathbf{u}}_i + \mathbf{C}\dot{\mathbf{u}}_i + \mathbf{K}\mathbf{u}_i = \mathbf{f}_i, \qquad (20)$$

where $\mathbf{u}_i = \mathbf{u}(t_i)$, $\dot{\mathbf{u}}_i = \dot{\mathbf{u}}(t_i)$ and $\ddot{\mathbf{u}}_i = \ddot{\mathbf{u}}(t_i)$. The conventional strategy is an incremental procedure in time, evaluating the solution $\mathbf{u}_i$ based on the already known solution of all preceding time steps in a step-by-step manner, as shown in Fig. 1.

In between the fixed time instants $t_i$, polynomial approximations are chosen. Within each time interval $[t_{i-1}, t_i]$, a second order approximation of the solution function is chosen using the implicit Newmark ansatz functions, cf. [66,67]:

$$\mathbf{u}_i = \mathbf{u}_{i-1} + \Delta t \dot{\mathbf{u}}_{i-1} + (0.5 - \beta) \Delta t^2 \ddot{\mathbf{u}}_{i-1}$$

<figure><img src="image_2.png" /><figcaption>Fig. 1 Conventional step-by-step integration procedure</figcaption></figure>

$$+ \beta \Delta t^2 \ddot{\mathbf{u}}_i, \qquad (21)$$

$$\dot{\mathbf{u}}_i = \dot{\mathbf{u}}_{i-1} + (1 - \gamma) \Delta t \ddot{\mathbf{u}}_{i-1} + \gamma \Delta t \dot{\mathbf{u}}_i. \qquad (22)$$

The parameters $\beta$ and $\gamma$ are chosen as $\frac{1}{4}$ and $\frac{1}{2}$ in this paper. These settings result in constant acceleration within each interval $[t_{i-1}, t_i]$ and unconditional stability with respect to the size of the chosen time step length $\Delta t$ [1,68,69].

The displacement increment and its derivatives at time instant $t_i$ are introduced as:

$$\Delta \mathbf{u}_i = \mathbf{u}_i - \mathbf{u}_{i-1}, \qquad (23)$$

$$\Delta \dot{\mathbf{u}}_i = \dot{\mathbf{u}}_i - \dot{\mathbf{u}}_{i-1}, \qquad (24)$$

$$\Delta \ddot{\mathbf{u}}_i = \ddot{\mathbf{u}}_i - \ddot{\mathbf{u}}_{i-1}, \qquad (25)$$

and the incremental force at time instant $t_i$ is evaluated as:

$$\Delta \mathbf{f}_i = \mathbf{f}_i - \mathbf{f}_{i-1}. \qquad (26)$$

The set of discretized equations of motion (20) at time instant $t_i$ is rewritten in incremental form as:

$$\mathbf{M} \Delta \ddot{\mathbf{u}}_i + \mathbf{C} \Delta \dot{\mathbf{u}}_i + \mathbf{K} \Delta \mathbf{u}_i = \Delta \mathbf{f}_i. \qquad (27)$$

Equivalently, the Newmark formulations (21) and (22) are rewritten in incremental form as:

$$\Delta \ddot{\mathbf{u}}_i = \frac{1}{\beta \Delta t^2} \Delta \mathbf{u}_i - \frac{1}{\beta \Delta t} \dot{\mathbf{u}}_{i-1} - \frac{1}{2\beta} \ddot{\mathbf{u}}_{i-1}, \qquad (28)$$

$$\Delta \dot{\mathbf{u}}_i = \frac{\gamma}{\beta \Delta t} \Delta \mathbf{u}_i - \frac{\gamma}{\beta} \dot{\mathbf{u}}_{i-1} + \Delta t \left(1 - \frac{\gamma}{2\beta}\right) \ddot{\mathbf{u}}_{i-1}. \qquad (29)$$

The three unknowns $\Delta \mathbf{u}_i$, $\Delta \dot{\mathbf{u}}_i$ and $\Delta \ddot{\mathbf{u}}_i$ are solved using the three equations (27), (28) and (29).

# 3 A non-incremental approach in elasto-dynamics

We start with an illustration, shown in Fig. 2. On the left hand side of this figure, a cantilever beam as an example of our body $\Omega$ with a time dependent punctual force function at its end, is depicted. On the right hand side, the state of this cantilever is shown at the time instant $t = 0$. Due to the excitation time history of the punctual force, the body $\Omega$ undergoes a certain dynamic response along the time axis
<figure><img src="image_3.png" /><figcaption>Fig. 2 Illustration of a cantilever beam in the space-time framework $\Omega \times [0, T]$; space-time displacement vector field $U(\mathbf{x}, t)$; spatial and temporal functions $u_s(\mathbf{x})$ and $u_t(t)$ highlighted in red and blue color, respectively</figcaption></figure>

$t$, as indicated by the continuous space-time response surface $U(\mathbf{x}, t)$. The goal is to provide this space-time response surface introducing a methodology that does not proceed in a step-by-step manner but that solves the whole space-time problem at once. In doing so, a non-incremental Newmark formulation in elasto-dynamics is proposed.

## 3.1 Continuous space-time decomposition

Let us consider a continuous system of which the evolution is described by a space-time displacement vector field $(\mathbf{x}, t) \mapsto U(\mathbf{x}, t)$. Let $\mathbf{x} \mapsto u_s(\mathbf{x})$ and $t \mapsto u_t(t)$ be spatial and temporal functions, respectively. Firstly, we define the tensor product $u_s \otimes u_t$ such that

$$ (u_s \otimes u_t)(\mathbf{x}, t) = u_s(\mathbf{x}) \circ u_t(t), \quad (30) $$

where the right hand side product is being performed component-wise indicated by the symbol $\circ$.

Observing Fig. 2, one can easily see that, in a first attempt, the space-time displacement vector field $U(\mathbf{x}, t)$ can be approximated by the tensor product of two functions $u_s(\mathbf{x})$ and $u_t(t)$, indicated by the red and blue lines in Fig. 2. As the displacement vector field $U(\mathbf{x}, t)$ can be arbitrary complex, an approximation of only one tensorial product could result in a bad representation of the solution. Thus, we claim that one can define the exact response time history by an infinite sum of displacement enrichments:

$$ U(\mathbf{x}, t) = \sum_{j=1}^{\infty} U^{(j)}(\mathbf{x}, t), \quad (31) $$

where each enrichment term is a tensor product $U^{(j)}(\mathbf{x}, t) = (u_s \otimes u_t)^{(j)}$ such that

$$ (u_s \otimes u_t)^{(j)}(\mathbf{x}, t) = (u_s)^{(j)}(\mathbf{x}) \circ (u_t)^{(j)}(t). \quad (32) $$

The response patterns in elastodynamics strongly depend on the type of the structure as well as the excitation time history. Following this line of thought, it is reasonable to hypothesize that a finite or even small number of enrichments will lead to a satisfying approximation of the displacement vector field $U(\mathbf{x}, t)$. This inevitably leads to the idea of model order reduction.

The separation ansatz allows one to easily define the derivatives of the displacement fields with respect to time, i.e. velocity and acceleration fields:

$$ \dot{U}(\mathbf{x}, t) = \sum_{j=1}^{\infty} (u_s)^{(j)}(\mathbf{x}) \circ (\dot{u}_t)^{(j)}(t), \quad (33) $$

$$ \ddot{U}(\mathbf{x}, t) = \sum_{j=1}^{\infty} (u_s)^{(j)}(\mathbf{x}) \circ (\ddot{u}_t)^{(j)}(t). \quad (34) $$

Equivalently to the displacement velocity and acceleration fields, we define the decomposition of the force vector field $(\mathbf{x}, t) \mapsto \mathbf{f}(\mathbf{x}, t)$. Let $\mathbf{x} \mapsto \mathbf{f}_s(\mathbf{x})$ and $t \mapsto \mathbf{f}_t(t)$ be the spatial and temporal force functions. Applying the definition of the tensor product from Eq. (30), the force vector field is also defined as a sum of enrichments:

$$ \mathbf{f}(\mathbf{x}, t) = \sum_{j=1}^{\infty} (\mathbf{f}_s)^{(j)}(\mathbf{x}) \circ (\mathbf{f}_t)^{(j)}(t). \quad (35) $$

The case of the stiffness is less straightforward because it is a functional linear operator mapping the space of displacement fields onto the space of the force fields. Considering a purely static problem, it provides the resulting force at a position $\mathbf{x}$ from a displacement acting at position $\mathbf{x}'$. It may be seen as the kernel $(\mathbf{x}, \mathbf{x}', t) \mapsto K(\mathbf{x}, \mathbf{x}', t)$ (or Green's function) of a functional operator $u \mapsto \mathbf{f}$ such that:

$$ \mathbf{f}(\mathbf{x}, t) = \int_{\Omega} K(\mathbf{x}, \mathbf{x}', t) u(\mathbf{x}', t) d\Omega, \quad (36) $$

where $\Omega$ is the body defined in equation (1), on which the contracted product of the two rank tensor field $K$ and the vector field $u$ occurs. Let $(\mathbf{x}, \mathbf{x}') \mapsto K_s(\mathbf{x}, \mathbf{x}')$ and $t \mapsto k_t(t)$ be spatial and temporal functions, respectively. In general, the tensor product $K_s \otimes k_t$ is such that:

$$ (K_s \otimes k_t)(\mathbf{x}, \mathbf{x}', t) = K_s(\mathbf{x}, \mathbf{x}') \circ K_t(t), \quad (37) $$

where the product in the right hand side being performed component-wise. The evolution of the stiffness is then also
performed in terms of an infinite sum of enrichments

$$K(\mathbf{x}, \mathbf{x}', t) = \sum_{j=0}^{\infty} K^{(j)}(\mathbf{x}, \mathbf{x}', t), \quad (38)$$

each of which is decomposed by $K^{(j)} = (K_s \otimes K_t)^{(j)}$ so that:

$$(K_s \otimes k_t)^{(j)}(\mathbf{x}, \mathbf{x}', t) = (K_s)^{(j)}(\mathbf{x}, \mathbf{x}') \circ (K_t)^{(j)}(t). \quad (39)$$

In this context, we note that the decomposition of the mass is equivalent to the decomposition of the stiffness.

## 3.2 Discrete space-time decomposition

In its discrete version, there are $n_s$ degrees of freedom known at $n_t$ instants of time. The discrete representative of the continuous field $U$ is a $n_s \times n_t$ matrix $U$:

$$U = [\mathbf{u}_1, \dots, \mathbf{u}_{n_t}] \in \mathbb{R}^{n_s \times n_t}, \quad (40)$$

where $\mathbf{u}_i = \mathbf{u}(t_i)$ is the vector containing the displacement data at the consecutive time instants $t_i$ of the discrete model. The discrete representatives of the spatial and temporal functions $u_s$ and $u_t$ are the column vectors:

$$\mathbf{u}_s \in \mathbb{R}^{n_s}, \quad \mathbf{u}_t \in \mathbb{R}^{n_t}. \quad (41)$$

The discrete counterpart of the tensor product $U^{(j)}(\mathbf{x}, t) = (u_s \otimes u_t)^{(j)}$ of functions is the tensor product of vectors:

$$U^{(j)} = \mathbf{u}_s^{(j)} \otimes \mathbf{u}_t^{(j)} \in \mathbb{R}^{n_s \times n_t}. \quad (42)$$

The object $U^{(j)}$ is one $n_s \times n_t$ matrix of which the element $u_{kl}$ is the product of $(u_s)_k$ and $(u_t)_l$. The discrete representation of equations (31) and (32) is then the summation of $M$ tensor products:

$$U = \sum_{j=1}^{M} \mathbf{u}_s^{(j)} \otimes \mathbf{u}_t^{(j)}. \quad (43)$$

Each summation term can be seen as an enrichment by a discrete space-time solution field. An upper boundary for the number $M$, necessary for an exact representation of $U$ by the sum of dyadic vector product enrichments (43), is the number of the spatial degrees of freedom times the number of the time steps. Thus, the condition holds:

$$M \le n_s \cdot n_t. \quad (44)$$

However, it is the objective of this paper to find a significantly smaller number $m \ll M$ that enables for an approximation of the space-time displacement matrix $U$, while preserving the

required level of accuracy for problems in structural dynamics. In Sect. 4, we will show that indeed a surprisingly small number of enrichments is sufficient for a reliable response identification of the dynamical system.

Equivalently, the discrete representations of the continuous velocity and acceleration fields, (33) and (34), are written as:

$$\dot{\mathbf{U}} = \sum_{j=1}^{M} \mathbf{u}_s^{(j)} \otimes \dot{\mathbf{u}}_t^{(j)}, \quad (45)$$

$$\ddot{\mathbf{U}} = \sum_{j=1}^{M} \mathbf{u}_s^{(j)} \otimes \ddot{\mathbf{u}}_t^{(j)}. \quad (46)$$

The discrete version of the continuous operator $K(\mathbf{x}, \mathbf{x}', t)$ is the usual stiffness matrix $K$, as presented in Eq. (15). Resuming the framework developed for the displacement and taking into account that the stiffness operator depends on two spatial variables and one time variable, we claim that the discrete representatives of the functions $(K_s)^{(j)}$ and $(K_t)^{(j)}$ are:

$$\mathbf{K}_s^{(j)} \in \mathbb{R}^{n_s \times n_s}, \quad \mathbf{k}_t^{(j)} \in \mathbb{R}^{n_t}. \quad (47)$$

Then, the discrete counterpart of the tensor product $(K_s)^{(j)} \otimes (k_t)^{(j)}$ of functions is:

$$\mathbf{K}_s^{(j)} \otimes \mathbf{k}_t^{(j)} \in \mathbb{R}^{n_s \times n_s \times n_t}. \quad (48)$$

Thus, we claim that the stiffness in space-time is a tensor of third order evaluated as a sum of decomposed enrichments:

$$\mathbf{K} = \sum_{j=1}^{M} \mathbf{K}_s^{(j)} \otimes \mathbf{k}_t^{(j)}. \quad (49)$$

Equivalently, the space-time mass is defined as a tensor of third order:

$$\mathbf{M} = \sum_{j=1}^{M} \mathbf{M}_s^{(j)} \otimes \mathbf{m}_t^{(j)}. \quad (50)$$

As discussed in Appendix A, the following rule can be applied:

$$(K_s \otimes k_t)(\mathbf{u}_s \otimes \mathbf{u}_t) = (K_s \mathbf{u}_s) \otimes (k_t \circ \mathbf{u}_t). \quad (51)$$

This is a general rule and, in this paper, it is applicable to all third order tensors and the respective tensor product of vectors. It will be applied in the upcoming sections of the paper.
In the special case of linear elastodynamics, stiffness and mass are both invariant with respect to time. The third order tensors are defined by one tensor product only:

$$ \mathbf{K}^3 = \mathbf{K}_s \otimes \mathbf{k}_t, \quad \mathbf{K}_s = \mathbf{K}, \quad \mathbf{k}_t = \mathbf{1}_t, \qquad (52) $$

$$ \mathbf{M}^3 = \mathbf{M}_s \otimes \mathbf{m}_t, \quad \mathbf{M}_s = \mathbf{M}, \quad \mathbf{m}_t = \mathbf{1}_t. \qquad (53) $$

Here, $\mathbf{1}_t \in \mathbb{R}^{n_t}$ demonstrates a vector with every entry equal to one. Taking Rayleigh damping into account, as defined in Eq. (18), also the damping tensor is defined as:

$$ \mathbf{C}^3 = \mathbf{C}_s \otimes \mathbf{c}_t, \quad \mathbf{C}_s = \mathbf{C}, \quad \mathbf{c}_t = \mathbf{1}_t. \qquad (54) $$

### 3.3 Space-time set of equations of motion

Taking into account the stiffness, mass and damping terms, presented in Eqs. (52)–(54), as well as the discrete representation of the decomposed displacement, velocity and acceleration fields, presented in Equations (43)–(46), the set of equations of motion in tensorial space-time formalism is written as:

$$ \begin{aligned} & (\mathbf{M}_s \otimes \mathbf{m}_t) \left( \sum_{j=1}^{M} \mathbf{u}_s^{(j)} \otimes \ddot{\mathbf{u}}_t^{(j)} \right) + \\ & (\mathbf{C}_s \otimes \mathbf{c}_t) \left( \sum_{j=1}^{M} \mathbf{u}_s^{(j)} \otimes \dot{\mathbf{u}}_t^{(j)} \right) + \\ & (\mathbf{K}_s \otimes \mathbf{k}_t) \left( \sum_{j=1}^{M} \mathbf{u}_s^{(j)} \otimes \mathbf{u}_t^{(j)} \right) = \mathbf{F}(t), \end{aligned} \qquad (55) $$

where $\mathbf{F}(t) = [\mathbf{f}(t_1), \dots, \mathbf{f}(t_{n_t})]$ denotes the discretized excitation history matrix. As the mass, the damping and the stiffness terms are invariant with respect to time, they can be moved into the sum of the respective time-acceleration-, time-velocity- and time-displacement decomposition terms. Applying the rule for the product of tensors, defined in equation (51), leads to a simplified version of Eq. (55) that is only represented by tensors of rank two:

$$ \mathbf{M}_s \ddot{\mathbf{U}} + \mathbf{C}_s \dot{\mathbf{U}} + \mathbf{K}_s \mathbf{U} = \mathbf{F}. \qquad (56) $$

### 3.4 Newmark space-time formulations

Firstly, the displacement increment $\Delta \mathbf{u}$ is transformed into its space-time equivalent. In matrix form, it is written as:

$$ \Delta \mathbf{U} = [\mathbf{u}(t_1) - \mathbf{u}(t_0), \dots, \mathbf{u}(t_{n_t}) - \mathbf{u}(t_{n_t-1})]. \qquad (57) $$

The transformation from the incremental formulation into the space-time equivalent is then defined as:

$$ \Delta \mathbf{u} = \mathbf{u}_i - \mathbf{u}_{i-1} \rightarrow \mathbf{I}_s \mathbf{u}_s^{(j)} \otimes \mathbf{A}_t \mathbf{u}_t^{(j)} - \mathbf{u}_0 \otimes \mathbf{a}_0, \qquad (58) $$

where $\mathbf{I}_s$ denotes the space identity matrix $n_s \times n_s$, $\mathbf{A}_t$ is a time coefficient matrix with the dimension $n_t \times n_t$ and $\mathbf{a}_0$ is an $n_t$-dimensional coefficient vector considering the initial displacement $\mathbf{u}_0$:

$$ \mathbf{A}_t = \begin{bmatrix} 1 & 0 & 0 & 0 & \cdots & 0 \\ -1 & 1 & 0 & 0 & \cdots & 0 \\ 0 & -1 & 1 & 0 & \cdots & 0 \\ \vdots & \ddots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & 0 & -1 & 1 & 0 \\ 0 & \cdots & 0 & 0 & -1 & 1 \end{bmatrix}, \quad \mathbf{a}_0 = \begin{bmatrix} 1 \\ 0 \\ 0 \\ \vdots \\ 0 \\ 0 \end{bmatrix} \qquad (59) $$

For clarification, the Newmark interpolations (28) and (29) are represented by five terms, $\mathbf{t}_1-\mathbf{t}_5$, each of which are transformed into their space-time equivalent, adapting the above proposed approach. In doing so, the Newmark formulations, taken from Eqs. (21) and (22), are rewritten as:

$$ \underbrace{\frac{\beta}{\Delta t} (\mathbf{u}_i - \mathbf{u}_{i-1})}_{\mathbf{t}_1} - \underbrace{\frac{\beta \dot{\mathbf{u}}_{i-1}}{\Delta t} +}_{\mathbf{t}_2} \beta \Delta t \left( \left( \beta - \frac{1}{2} \right) \ddot{\mathbf{u}}_{i-1} - \beta \ddot{\mathbf{u}}_i \right) = 0, \qquad (60) $$

$$ \underbrace{\frac{\beta}{\Delta t} (\dot{\mathbf{u}}_i - \dot{\mathbf{u}}_{i-1})}_{\mathbf{t}_4} + \underbrace{\beta ((\gamma - 1) \ddot{\mathbf{u}}_{i-1} - \gamma \ddot{\mathbf{u}}_i))}_{\mathbf{t}_5} = 0. \qquad (61) $$

Equivalent to the transformation (58), the space-time formulations of the defined terms in Eqs. (60) and (61) are defined as:

$$ \mathbf{t}_1 \rightarrow \mathbf{I}_s \mathbf{u}_s^{(j)} \otimes \mathbf{A}_1 \mathbf{u}_t^{(j)} - \mathbf{u}_0 \otimes \mathbf{a}_5, \qquad (62) $$

$$ \mathbf{t}_2 \rightarrow \mathbf{I}_s \mathbf{u}_s^{(j)} \otimes \mathbf{A}_2 \dot{\mathbf{u}}_t^{(j)} - \dot{\mathbf{u}}_0 \otimes \mathbf{a}_6, \qquad (63) $$

$$ \mathbf{t}_3 \rightarrow \mathbf{I}_s \mathbf{u}_s^{(j)} \otimes \mathbf{A}_3 \ddot{\mathbf{u}}_t^{(j)} - \ddot{\mathbf{u}}_0 \otimes \mathbf{a}_7, \qquad (64) $$

$$ \mathbf{t}_4 \rightarrow \mathbf{I}_s \mathbf{u}_s^{(j)} \otimes \mathbf{A}_1 \dot{\mathbf{u}}_t^{(j)} - \dot{\mathbf{u}}_0 \otimes \mathbf{a}_5, \qquad (65) $$

$$ \mathbf{t}_5 \rightarrow \mathbf{I}_s \mathbf{u}_s^{(j)} \otimes \mathbf{A}_4 \ddot{\mathbf{u}}_t^{(j)} - \ddot{\mathbf{u}}_0 \otimes \mathbf{a}_8. \qquad (66) $$

The matrices $\mathbf{A}_1, \mathbf{A}_2, \mathbf{A}_3$ and $\mathbf{A}_4$ are space-time Newmark coefficient matrices, containing $\beta$ and $\gamma$, as well as the time step size $\Delta t$. They are all of dimension $n_t \times n_t$ and listed in Appendix B. Furthermore, the space-time equivalents of the initial conditions $\mathbf{u}_0, \dot{\mathbf{u}}_0, \ddot{\mathbf{u}}_0$ are defined by the vectors $\mathbf{a}_5 - \mathbf{a}_8$. Regarding all four vectors only the first entry is not equal to zero. They are listed in Appendix C. We note that $\mathbf{A}_1-\mathbf{A}_4$ are sparse matrices as most of their components are equal to
zero. This property can be considered when implementing the proposed strategy in order to further increase the level of efficiency.

Inserting the space-time formulations in equations (62–66) into the Newmark equations leads to the space-time equivalent of the Newmark equations:

$$\mathbf{U}\mathbf{A}_1^T + \dot{\mathbf{U}}\mathbf{A}_2^T + \ddot{\mathbf{U}}\mathbf{A}_3^T = \mathbf{u}_0\mathbf{a}_5^T + \dot{\mathbf{u}}_0\mathbf{a}_6^T + \ddot{\mathbf{u}}_0\mathbf{a}_7^T, \quad (67)$$

$$\dot{\mathbf{U}}\mathbf{A}_1^T + \ddot{\mathbf{U}}\mathbf{A}_4^T = \dot{\mathbf{u}}_0\mathbf{a}_5^T + \ddot{\mathbf{u}}_0\mathbf{a}_8^T. \quad (68)$$

We note again that, in these equations, $\mathbf{U}$, $\dot{\mathbf{U}}$ and $\ddot{\mathbf{U}}$ are identified as rank two tensors, represented as $n_s \times n_t$ matrices.

The discretized space time equation We identify three space time Eqs. (56), (67) and (68) with three unknowns, i.e. space-time displacement, $\mathbf{U}$, space-time velocity $\dot{\mathbf{U}}$ and space-time acceleration $\ddot{\mathbf{U}}$. However, we take advantage of the fact that, velocity and acceleration can be written as a function of displacement. Therefore, velocity and acceleration are expressed from the two space-time Newmark interpolations, (67) and (68), with respect to displacement and inserted into the space-time equation of motion (55). In doing so, the linear dynamic problem is expressed by one algebraic space-time equation only:

$$\mathbf{K}_s \mathbf{U} - \mathbf{C}_s \mathbf{U} \mathbf{W} + \mathbf{M}_s \mathbf{U} \mathbf{Y} = \mathbf{L}. \quad (69)$$

In this equation, matrices $\mathbf{W}$ and $\mathbf{Y}$ are evaluated as:

$$\mathbf{W} = \mathbf{A}_1^T \mathbf{D}^{-1} \quad \text{and} \quad \mathbf{Y} = \mathbf{A}_1^T \mathbf{D}^{-1} \mathbf{H}, \quad (70)$$

where

$$\mathbf{D} = \mathbf{A}_2^T - \mathbf{A}_1^T (\mathbf{A}_4^T)^{-1} \mathbf{A}_3^T, \quad (71)$$

$$\mathbf{H} = \mathbf{A}_1^T (\mathbf{A}_4^T)^{-1}. \quad (72)$$

The $n_s \times n_t$ matrix $\mathbf{L}$ is identified as a space-time excitation term and is evaluated as:

$$\mathbf{L} = \mathbf{F} + \sum_{d=1}^{8} \mathbf{L}_d, \quad (73)$$

with

$$\mathbf{L}_1 = \mathbf{M}_s \mathbf{u}_0 \mathbf{a}_5^T \mathbf{D}^{-1} \mathbf{H}, \quad \mathbf{L}_2 = -\mathbf{C}_s \mathbf{u}_0 \mathbf{a}_5^T \mathbf{D}^{-1}, \quad (74)$$

$$\mathbf{L}_3 = \mathbf{M}_s \dot{\mathbf{u}}_0 \mathbf{e} \mathbf{D}^{-1} \mathbf{H}, \quad \mathbf{L}_4 = -\mathbf{M}_s \dot{\mathbf{u}}_0 \mathbf{b}, \quad (75)$$

$$\mathbf{L}_5 = -\mathbf{C}_s \dot{\mathbf{u}}_0 \mathbf{e} \mathbf{D}^{-1}, \quad \mathbf{L}_6 = -\mathbf{M}_s \dot{\mathbf{u}}_0 \mathbf{n}, \quad (76)$$

$$\mathbf{L}_7 = -\mathbf{C}_s \ddot{\mathbf{u}}_0 \mathbf{g} \mathbf{D}^{-1}, \quad \mathbf{L}_8 = \mathbf{M}_s \ddot{\mathbf{u}}_0 \mathbf{g} \mathbf{D}^{-1} \mathbf{H}. \quad (77)$$

The vectors $\mathbf{e}$, $\mathbf{g}$, $\mathbf{n}$ and $\mathbf{b}$ are dependent on the coefficient matrices $\mathbf{A}_1-\mathbf{A}_4$ and the coefficient vectors $\mathbf{a}_5-\mathbf{a}_8$:

$$\mathbf{e} = \mathbf{a}_6^T - \mathbf{a}_5^T (\mathbf{A}_4^T)^{-1} \mathbf{A}_3^T, \quad (78)$$

$$\mathbf{g} = \mathbf{a}_7^T - \mathbf{a}_8^T (\mathbf{A}_4^T)^{-1} \mathbf{A}_3^T, \quad (79)$$

$$\mathbf{n} = \mathbf{a}_8^T (\mathbf{A}_4^T)^{-1}, \quad (80)$$

$$\mathbf{b} = \mathbf{a}_5^T (\mathbf{A}_4^T)^{-1}. \quad (81)$$

We note that the right hand side of equation (69) is still not decomposed using the space-time separation ansatz. Therefore, the proper orthogonal decomposition [15,70,71] is applied. In doing so, the space-time excitation term on the right hand side of the space-time set of equations of motion (69) is decomposed into a sum of dyadic products of uncorrelated space and time vectors. This is realized using the singular value decomposition [72]:

$$\mathbf{L} = \mathbf{L}_s \mathbf{S} \mathbf{L}_t \approx \sum_{q=1}^{n_{tr}} s^{(q)} \mathbf{I}_s^{(q)} \otimes \mathbf{I}_t^{(q)}. \quad (82)$$

Here, $\mathbf{I}_s^{(q)}$ are orthonormal left singular vectors, taken from the matrix $\mathbf{L}_s$ with dimension $n_s \times n_s$ and $\mathbf{I}_t^{(q)}$ are the right singular vectors taken from the matrix $\mathbf{L}_t$ with dimension $n_t \times n_t$. The matrix $\mathbf{S}$ is an $n_s \times n_t$ rectangular diagonal matrix with the singular values $s^{(q)}$ in descending order. In this paper, the number of singular basis vectors $n_{tr}$ is chosen so that 99.9% of the total energy, evaluated by the sum of all singular values $\sum_{q=1}^{n_{tr}} s^{(q)}$, is captured [10]. In structural dynamics it is, generally, the case that an already surprisingly small number of basis vectors is sufficient to obtain such a high level of accuracy. [11,12,61].

## 3.5 The space-time solution strategy

In this paper, the fixed point algorithm is applied to solve the problem in the space-time framework. The solution at the $m^{th}$ enrichment, $\mathbf{U}^{(m)}$, is divided into already evaluated $m-1$ known enrichments and the current $m^{th}$ enrichment solution to be calculated:

$$\mathbf{U}^{(m)} = \sum_{j=1}^{m-1} \mathbf{u}_s^{(j)} \otimes \mathbf{u}_t^{(j)} + \mathbf{u}_s^{(m)} \otimes \mathbf{u}_t^{(m)}. \quad (83)$$

Figure 3 visualizes the iterative enrichment procedure as a step-wise adaption of the solution in the whole time domain.

The solution ansatz (83) as well as the decomposition of the space-time excitation term (82) are inserted into the space-time equation (69):

$$\mathbf{K}_s \left[ \sum_{j=1}^{m-1} \mathbf{u}_s^{(j)} \otimes \mathbf{u}_t^{(j)} + \mathbf{u}_s^{(m)} \otimes \mathbf{u}_t^{(m)} \right] - \mathbf{C}_s \left[ \sum_{j=1}^{m-1} \mathbf{u}_s^{(j)} \otimes \mathbf{u}_t^{(j)} + \mathbf{u}_s^{(m)} \otimes \mathbf{u}_t^{(m)} \right] \mathbf{W}$$
<figure><img src="image_4.png" /><figcaption>Fig. 3 Iterative space time procedure</figcaption></figure>

$$ 
\begin{aligned}
& + \mathbf{M}_s \left[ \sum_{j=1}^{m-1} \mathbf{u}_s^{(j)} \otimes \mathbf{u}_t^{(j)} + \mathbf{u}_s^{(m)} \otimes \mathbf{u}_t^{(m)} \right] \mathbf{Y} \\
& = \sum_{q=1}^{n_{tr}} s^{(q)} \mathbf{l}_s^{(q)} \otimes \mathbf{l}_t^{(q)}.
\end{aligned} 
\quad (84) $$

Equation (84) has two unknowns, the spatial enrichment vector $\mathbf{u}_s^{(m)}$ and the temporal enrichment vector $\mathbf{u}_t^{(m)}$. Multiplication with $(\mathbf{u}_s^{(m)})^T$ leads to a projection into time:

$$ 
\mathbf{u}_t^{(m)} = (a\mathbf{I}_t - b\mathbf{W}^T + c\mathbf{Y}^T)^{-1} (\hat{\mathbf{I}}_t - \mathbf{v}_{s1} + \mathbf{v}_{s2} - \mathbf{v}_{s3}). \quad (85) $$

Analogously, multiplication with $(\mathbf{u}_t^{(m)})^T$ leads to a projection into space:

$$ 
\mathbf{u}_s^{(m)} = (e\mathbf{K}_s - f\mathbf{C}_s + g\mathbf{M}_s)^{-1} (\hat{\mathbf{I}}_s - \mathbf{v}_{t1} + \mathbf{v}_{t2} - \mathbf{v}_{t3}). \quad (86) $$

Two coupled Eqs. (85) and (86) with two unknowns $\mathbf{u}_s^{(m)}$ and $\mathbf{u}_t^{(m)}$ are obtained. In order to efficiently solve this coupled space-time equation pair, the fixed-point algorithm is applied. A one-dimensional representation of the algorithm is depicted in Fig. 4. As also shown in this figure, one starts by randomly choosing all entries of $\mathbf{u}_s$. The vector $\mathbf{u}_t^{(m)}$ is then evaluated by solving Eq. (85) while keeping $\mathbf{u}_s^{(m)}$ fixed. This step is indicated by the first vertical arrow in Fig. 4. In the next step, $\mathbf{u}_s^{(m)}$ is evaluated by solving equation (86) while keeping $\mathbf{u}_t^{(m)}$ fixed. This procedure is indicated by the first horizontal arrow in Fig. 4. These two steps are alternately repeated until convergence is achieved, indicated by the intersection point $(u_s^*, u_t^*)$ of the red and blue curves in Fig. 4.

In Eqs. (85) and (86) the scalar quantities a, b, c, e, f and g are evaluated by:

$$ a = (\mathbf{u}_s^{(m)})^T \mathbf{K}_s \mathbf{u}_s^{(m)}, \quad b = (\mathbf{u}_s^{(m)})^T \mathbf{C}_s \mathbf{u}_s^{(m)}, \quad (87) $$

<figure><img src="image_5.png" /><figcaption>Fig. 4 One-dimensional illustration of the fixed point algorithm to solve the coupled equation pair (85) and (86)</figcaption></figure>

$$ c = (\mathbf{u}_s^{(m)})^T \mathbf{M}_s \mathbf{u}_s^{(m)}, \quad e = (\mathbf{u}_t^{(m)})^T \mathbf{I}_t \mathbf{u}_t^{(m)}, \quad (88) $$

$$ f = (\mathbf{u}_t^{(m)})^T \mathbf{W} \mathbf{u}_t^{(m)}, \quad g = (\mathbf{u}_t^{(m)})^T \mathbf{Y} \mathbf{u}_t^{(m)}. \quad (89) $$

The excitation vectors, projected into space and time for the mth enrichment, are written as:

$$ \hat{\mathbf{I}}_s = \sum_{q=1}^{Tr} s^{(q)} (\mathbf{l}_s^{(q)})^T \mathbf{u}_t^{(m)} \mathbf{l}_s^{(q)}, \quad (90) $$

$$ \hat{\mathbf{I}}_t = \sum_{q=1}^{Tr} s^{(q)} (\mathbf{l}_s^{(q)})^T \mathbf{u}_s^{(m)} \mathbf{l}_t^{(q)}. \quad (91) $$

The vectors $\mathbf{v}_{s1} - \mathbf{v}_{s3}$ in equation (85) and $\mathbf{v}_{t1} - \mathbf{v}_{t3}$ in Eq. (86) are written as:

$$ \mathbf{v}_{s1} = \sum_{j=1}^{m-1} (\mathbf{u}_s^{(j)})^T \mathbf{K}_s \mathbf{u}_s^{(m)} \mathbf{I}_t \mathbf{u}_t^{(j)}, \quad (92) $$

$$ \mathbf{v}_{t1} = \sum_{j=1}^{m-1} (\mathbf{u}_t^{(j)})^T \mathbf{I}_t^T \mathbf{u}_t^{(m)} \mathbf{K}_s \mathbf{u}_s^{(j)}, \quad (93) $$

$$ \mathbf{v}_{s2} = \sum_{j=1}^{m-1} (\mathbf{u}_s^{(j)})^T \mathbf{C}_s \mathbf{u}_s^{(m)} \mathbf{W}^T \mathbf{u}_t^{(j)}, \quad (94) $$

$$ \mathbf{v}_{t2} = \sum_{j=1}^{m-1} (\mathbf{u}_t^{(j)})^T \mathbf{W} \mathbf{u}_t^{(m)} \mathbf{C}_s \mathbf{u}_s^{(j)}, \quad (95) $$

$$ \mathbf{v}_{s3} = \sum_{j=1}^{m-1} (\mathbf{u}_s^{(j)})^T \mathbf{M}_s \mathbf{u}_s^{(m)} \mathbf{Y}^T \mathbf{u}_t^{(j)}, \quad (96) $$
$$v_{t3} = \sum_{j=1}^{m-1} (\mathbf{u}_t^{(j)})^T \mathbf{Y} \mathbf{u}_t^{(m)} \mathbf{M}_s \mathbf{u}_s^{(j)}. \quad (97)$$

Convergence is achieved if the enrichment vectors $\mathbf{u}_t^{(m)}$ and $\mathbf{u}_s^{(m)}$ fulfill the Frobenius norm of the space-time equation of motion with sufficient accuracy. This is defined through a scalar parameter $\varepsilon_{fp}$:

$$\left\| \mathbf{K}_s \mathbf{u}_s^{(m)} \otimes \mathbf{u}_t^{(m)} - \mathbf{C}_s \mathbf{u}_s^{(m)} \otimes \mathbf{W}^T \mathbf{u}_t^{(m)} + \mathbf{M}_s \mathbf{u}_s^{(m)} \otimes \mathbf{Y}^T \mathbf{u}_t^{(m)} - \sum_{q=1}^{Tr} s^{(q)} \mathbf{l}_s^{(q)} \otimes \mathbf{l}_t^{(q)} \right\| < \varepsilon_{fp}. \quad (98)$$

After each enrichment step, the enrichment vectors $\mathbf{u}_t^{(m)}$ and $\mathbf{u}_s^{(m)}$ are inserted into Eq. (83) in order to update the total displacement function. Global convergence is achieved if the Frobenius norm of Eq. (84), denoted by the scalar quantity $\varepsilon_{enr}^{(FN)}$, is fulfilled with the required level of accuracy:

$$\underbrace{\left\| \mathbf{M}_s \ddot{\mathbf{U}}^{(m)} + \mathbf{C}_s \dot{\mathbf{U}}^{(m)} + \mathbf{K}_s \mathbf{U}^{(m)} - \mathbf{F} \right\|}_{\varepsilon_{enr}^{(FN)}} < \varepsilon_{enr}, \quad (99)$$

where $\varepsilon_{enr}$ is the global convergence limit to be defined.

Additionally, the proposed space-time algorithm is summarized in Algorithm 1.

## 4 Numerical examples

For the purpose of demonstration, we introduce a two-dimensional frame structure subjected to four different excitation functions. Regarding all the numerical examples, the proposed Newmark space-time approach is presented and compared to the conventional Newmark integration scheme. In doing so, the focus is not only on the development of the spatial and temporal enrichment vectors during convergence but also on the limits of the proposed strategy.

Figure 5a presents an illustration of the frame structure subjected to the different chosen excitations. In particular, three different horizontal ground excitations $\ddot{x}_g$ and one impact force function $f(t)$ are considered.

Introducing ground excitation, the external force on the right hand side of the set of equations of motion (19) is substituted by the term of inertia induced by the seismic acceleration of the ground [1]:

$$f(t) = -Mf_H \ddot{x}_g. \quad (100)$$

**Algorithm 1:** Implementation of the Newmark space-time integration scheme.

**Result:** Response time history U initialization;

calc. $\mathbf{M}_s, \mathbf{K}_s, \mathbf{C}_s$;

set initial conditions $\mathbf{u}_0$ and $\dot{\mathbf{u}}_0 \to$ calc. $\ddot{\mathbf{u}}_0$;

set coefficient matrices and vectors $\mathbf{A}_1 - \mathbf{A}_4$;

set coefficient vectors $\mathbf{a}_5 - \mathbf{a}_8$;

set enrichment counter $j = 1$;

assign memory location for $\mathbf{U}^{(m)}, \mathbf{u}_s^{(m)}$ and $\mathbf{u}_t^{(m)}$;

set first arbitrarily valued vector $\mathbf{u}_s^{(1)}, \mathbf{u}_s^{(1)} \neq 0$;

while $\varepsilon^{(m)} > \varepsilon_{enr}$ do

calc. projected excitation $\hat{\mathbf{u}}_s$ and $\hat{\mathbf{u}}_t$ using the POD approach;

set coefficients $a, b, c, e, f, g$;

set coefficients vectors $\mathbf{v}_{s1} - \mathbf{v}_{s3}$ and $\mathbf{v}_{t1} - \mathbf{v}_{t3}$;

set fixed point counter $\kappa = 1$;

while $\varepsilon^{(\kappa)} > \varepsilon_{fp}$ do

update $\mathbf{u}_t^{(m)} \to$ solve Eq. (85);

update $\mathbf{u}_s^{(m)} \to$ solve Eq. (86);

update $\varepsilon^{(\kappa)}$ using Eq. (98);

$\kappa + 1$;

end

update $\mathbf{U}^{(m)}$;

update $\varepsilon^{(m)}$ using Eq. (99);

$m + 1$;

end

<figure><img src="image_6.png" /><figcaption>Fig. 5 Overview of the numerical demonstrations considered in this paper; a representative plane frame structure; b chosen types of the excitation time histories: sinusoidal ground excitation (top), ground excitation recorded in Kobe, Japan [73] (middle), impact force function (bottom)</figcaption></figure>
In this paper, $\mathbf{f}_H$ is an influence vector, defining all horizontal degrees of freedom that are affected by the ground acceleration time history $\ddot{x}_g(t)$.

As first ground excitation example, a sinusoidal function $\ddot{x}_g^{(h)}(t) = a_0 \sin \nu t$ is chosen with an amplitude of $a_0 = 0.1 \text{ m/s}^2$ and an excitation frequency of $\nu = 1.61 \text{ rad/s}$. As second ground excitation example, also a sinusoidal function is chosen with the same amplitude of $a_0 = 0.1 \text{ m/s}^2$ but with a significantly higher excitation frequency compared to the first example. In this context, the excitation frequency for the second example is chosen as $5.0 \text{ rad/s}$. As third ground excitation example, a transient time history is chosen that is taken from an earthquake recorded in Kobe, Japan [73], as shown by the plot in the middle of Fig. 5b. The time duration of the record is 22 seconds and the peak ground acceleration (PGA) is $6.8 \text{ m/s}^2$. Finally, as fourth example, an impact force function $f(t)$, as shown in the bottom plot of Fig. 5b, is presented. In particular, a sinusoidal function is considered during the first 0.1 seconds followed by 0.9 seconds of a force function equal to zero. All excitations are also summarized in Table 1.

The frame structure is discretized by two-node truss-beam elements, applying the Euler-Bernoulli beam theory [74]. Both columns and beams of each story are discretized by five elements, which leads to a total number of 128 degrees of freedom. The time step $\Delta t$ is chosen as $10^{-2}$ seconds and the parameters of the numerical example are chosen as: stiffness $E = 210 \text{ GPa}$, density $\rho = 7850 \text{ kg/m}^3$. All cross sections of the columns and beams are modeled by quadratic hollow rectangular shapes with the dimension $0.3 \text{ m} \times 0.1 \text{ m}$ and a thickness of $t = 10 \text{ mm}$.

The output is compared to the conventional step-by-step Newmark solution. Regarding the first three examples, the output quantity considered for comparison is the horizontal displacement history in the left corner of the frame structure, as indicated by the magenta cross and the magenta arrow in Fig. 5a. Regarding the fourth example, the output quantity of interest considered for comparison is the degree of freedom $\nu(t)$ indicated by the red arrow in Fig. 5a. The convergence limit $\varepsilon_{enr}^{(\text{FN})}$ is calculated as the average value over all time steps. Thus, if the quantity $\varepsilon_{enr}^{(\text{FN})}/n_t$ falls below the value of $1.0e-4$, convergence is achieved.

Modal analysis is applied to the structure by solving the eigenvalue problem in Eq. (17). The solution is the set of eigenvalues (squared eigenfrequencies) in ascending order and the corresponding eigenvectors. In Fig. 6, the first four eigenvectors $\mathbf{\Phi}^{(1)} - \mathbf{\Phi}^{(4)}$ are depicted. The corresponding circular eigenfrequencies are evaluated as: $\omega^{(1)} = 1.61 \text{ rad/s}$, $\omega^{(2)} = 7.30 \text{ rad/s}$, $\omega^{(3)} = 11.85 \text{ rad/s}$, $\omega^{(4)} = 28.81 \text{ rad/s}$.

<figure><img src="image_7.png" /><figcaption>**Fig. 6** First four linear modes of vibration of the frame structure, $\mathbf{\Phi}^{(1)} - \mathbf{\Phi}^{(4)}$; the corresponding circular eigenfrequencies are: $\omega^{(1)} = 1.61 \text{ rad/s}$, $\omega^{(2)} = 7.30 \text{ rad/s}$, $\omega^{(3)} = 11.85 \text{ rad/s}$, $\omega^{(4)} = 28.81 \text{ rad/s}$</figcaption></figure>

## 4.1 Example 1: harmonic excitation-resonance

Firstly, the structure is excited by a sinusoidal ground excitation at a frequency of $\nu = 1.61 \text{ rad/s}$. This frequency is equal to the first eigenfrequency $\omega^{(1)}$ that corresponds to the first eigenmode, depicted in Fig. 6. For this numerical example, damping is neglected. The results are depicted in Fig. 7.

The Newmark space-time algorithm is presented by the dashed blue line in the bottom plot of Fig. 7, while the Newmark benchmark solution is presented by the black line. The space-time algorithm converges already after the first enrichment $m = 1$. This is caused by the fact that, in the special case of resonance with respect to the first eigenfrequency, the structure responds in the first mode shape only. In this special situation, the first normalized spatial enrichment vector coincides with the first normalized mode of vibration $\frac{\mathbf{u}^{(1)}}{\|\mathbf{u}^{(1)}\|} = \frac{\mathbf{\Phi}^{(1)}}{\|\mathbf{\Phi}^{(1)}\|}$. The corresponding temporal enrichment vector is depicted in the right top subplot in Fig. 7, clearly following the response characteristics of the resonance effect excited by a sinusoidal excitation. The full solution is then exactly evaluated by the dyadic product of the first spatial and temporal enrichment vectors $\mathbf{u}_s^{(1)} \otimes \mathbf{u}_t^{(1)}$.
<table><caption>**Table 1** Overview of the four applied excitations</caption><thead><tr><th>Ex.</th><th>Function</th><th>Type</th><th>Parameters</th></tr></thead><tbody><tr><td>1</td><td>x̃<sub>g</sub><sup>(h)</sup> = a<sub>0</sub> sin vt</td><td>Harmonic</td><td>a<sub>0</sub> = 0.1 m/s<sup>2</sup>, v = 1.61 rad/s</td></tr><tr><td>2</td><td>x̃<sub>g</sub><sup>(r)</sup> (record)</td><td>Harmonic</td><td>a<sub>0</sub> = 0.1 m/s<sup>2</sup>, v = 5.0 rad/s</td></tr><tr><td>3</td><td>f(t)</td><td>Earthquake</td><td>PGA ... 6.8 m/s<sup>2</sup></td></tr><tr><td>4</td><td>f(t)</td><td>Impact</td><td>max {f(t)} = 2.0e4 N</td></tr></tbody></table>

<figure><img src="image_8.png" /></figure>

## 4.2 Example 2: harmonic excitation-no resonance

Secondly, the structure is excited by a sinusoidal ground excitation, but, in this case, with a frequency that is higher than the first eigenfrequency. It is chosen as $v = 5.0$ rad/s. Also, regarding this example, damping is neglected. The results are depicted in Fig. 8.

Convergence is already achieved after two enrichments. The spatial enrichments modes are presented in red, while the corresponding temporal enrichments vectors are presented in blue. Both spatial modes, $u_s^{(1)}$ and $u_s^{(2)}$, are linearly independent and they show a high similarity with the first two linear

<caption>**Fig. 8** Results applying a sinusoidal excitation function $\ddot{x}_g(t) = a_0 \sin vt$ with $v = 5.0$ rad/s; spatial enrichment vectors in red, temporal enrichment vectors in blue; bottom plot: response function of the output degree of freedom $h(t)$ applying the Newmark space-time approach (blue dashed line) and the conventional Newmark algorithm (black solid line)</caption>
modes of vibration, shown in Fig. 6. This phenomenon is not surprising as, naturally, the response is mainly composed of the first two eigenmodes if the excitation frequency lies between the first and the second eigenfrequency. However, looking at the corresponding temporal enrichment functions shows that the proposed strategy converges differently than the conventional modal truncation method would, as both the first and the second spatial enrichment function contain more than only one frequency. The horizontal response time history $h(t)$ is shown in the bottom plot. The space-time solution is depicted by the dashed blue line, while the Newmark benchmark solution is represented by the black solid line. Also regarding this example, excellent agreement is observed when comparing the new strategy with the benchmark solution after the two enrichments. The whole convergence behavior is also presented in the bottom plot of Fig. 8, where the time step average of the Frobenius norm $\epsilon_{enr}^{(FN)}/n_t$ versus the number of considered enrichments is depicted.

### 4.3 Example 3: transient excitation

Thirdly, the structure is excited by an earthquake, recorded in Kobe, Japan [73], as shown in the middle plot of Fig. 5b. In this numerical example, Rayleigh damping is introduced using a damping ratio of 4% for the second and sixth mode of vibration [1]. The proposed Newmark space-time strategy already converges after the first two enrichments. However, we also demonstrate the results of the third and fourth enrichment for the purpose of illustration. The results are depicted in Fig. 9.

The horizontal response time history $h(t)$ after two enrichments is shown in the bottom plot of Fig. 9. The blue dashed line shows the Newmark space-time solution, while the black solid line shows the conventional Newmark solution. Also regarding this example, all spatial modes are linearly independent and a high similarity with the first four linear modes of vibration, shown in Fig. 6, is observed.

### 4.4 Example 4: impact force

Finally, the structure is excited by a nodal force time history $f(t)$, as shown in Fig. 5a. This excitation can be interpreted as an impact load with a duration of 0.1 seconds. While the three other examples induce global response patterns, exciting all horizontal degrees of freedom equally introducing a mass dependent term of inertia, the last example is characterized by a rather localized effect. The results regarding the fourth example are presented in Fig. 10.

Also here, the spatial enrichment modes are depicted in red color, while the temporal enrichment modes are depicted in blue color. All spatial modes are linearly independent. While, regarding the other three examples, the first enrich-

<figure><img src="image_9.png" /><figcaption>**Fig. 9** Results applying the transient excitation: spatial enrichment vectors in red, temporal enrichment vectors in blue; bottom plot: horizontal displacement response function $h(t)$</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Fig. 10 Results of the impact force function: first four spatial and temporal enrichment modes and the horizontal displacement response function $v(t)$</figcaption></figure>

<figure><img src="image_11.png" /></figure>

ment modes show a high similarity with the first few lowest modes of vibration, the enrichment modes of the present example clearly show a distinct localized effect due to the point transmission of the applied force time history. The response time histories of the Newmark space-time procedure and the conventional Newmark benchmark solution of the output degree of freedom $v(t)$ are shown in the bottom plot of Fig. 10. While the solution functions of the other three examples are converging after the first or second enrichment step, the space-time algorithm converges after 20 enrichment steps regarding this example. The downwards trend of the time step average of the Frobenius norm $\epsilon_{enr}^{(FN)}/n_t$ is considerably less distinct compared to the convergence trend of the other examples. However, the converged solution after the 20th enrichment is also here in excellent agreement with the step-by-step Newmark benchmark response, which indicates the robustness of the proposed algorithm. The classical method of modal truncation would, regarding the present example, not lead to a satisfying approximation of the response function as the local deformation patterns could kinematically not be represented by a spectral basis, evaluated by considering a small number of low-frequency modes from equation (17). Thus, an investigation of the mode shapes would a priori be necessary in order to extract precisely those modes that are responsible for these localized deformation pattern around the impact point. However, using the Newmark space-time approach, the localized deformation patterns are automatically included within the first spatial enrichment vectors as shown in Fig. 10 and no a priori investigation of the dynamic system is required. In this context, this example also shows the clear advantage of the Newmark space-time approach to the classical method of modal truncation.

## 4.5 Computational efficiency

The computational cost of the new strategy depends on the number of enrichments necessary for convergence and, therefore, on the convergence interval $\epsilon_{enr}^{(FN)}$ chosen. As a scalar error quantity $\epsilon_U$, we introduce the trace of the covariance matrix of the displacement matrix after the mth enrichment, $U^{(m)}$, which is equivalent to the Frobenius norm of the global error matrix:

$$ \epsilon_U = \frac{1}{n_t} \left\{ (U - U^{(m)}) (U - U^{(m)})^T \right\}^{\frac{1}{2}}, \quad (101) $$

where $U$ is the response time history obtained using the conventional step-by-step Newmark method. The computational efficiency is investigated, using the numerical example 3, presented in Sect. 4.3, by varying the number of enrichments considered and comparing the corresponding calculation time. The result is depicted in Fig. 11.
<figure><img src="image_12.png" /><figcaption>**Fig. 11** Computational efficiency versus accuracy of the proposed space-time procedure</figcaption></figure>

<table><caption>**Table 2** Computational efficiency of the Newmark space-time algorithm compared to the conventional Newmark method</caption><thead><tr><th>Num. enr.</th><th>ε<sub>U</sub></th><th>Rel. calc. t.</th><th>Speed-up</th></tr></thead><tbody><tr><td>1</td><td>2.81072 × 10<sup>-4</sup></td><td><strong>0.23</strong></td><td><strong>4.35</strong></td></tr><tr><td>2</td><td>1.76001 × 10<sup>-5</sup></td><td><strong>0.48</strong></td><td><strong>2.08</strong></td></tr><tr><td>3</td><td>5.81316 × 10<sup>-6</sup></td><td><strong>0.81</strong></td><td><strong>1.23</strong></td></tr><tr><td>4</td><td>7.37326 × 10<sup>-7</sup></td><td><strong>0.99</strong></td><td><strong>1.01</strong></td></tr><tr><td>5</td><td>7.34054 × 10<sup>-7</sup></td><td>1.33</td><td>0.75</td></tr><tr><td>6</td><td>7.29546 × 10<sup>-7</sup></td><td>1.63</td><td>0.61</td></tr><tr><td>7</td><td>7.29500 × 10<sup>-7</sup></td><td>2.02</td><td>0.50</td></tr></tbody></table>

The blue line presents the error quantity $ε_U$ versus the number of enrichments considered. A clear strong downward trend is observed within the first few enrichments. We note that the convergence criterion chosen for all examples ($ε_{enr}^{(FN)} ≤ 1.0e-4$) is here already satisfied after the second enrichment. However, we increase the number of enrichments for the sake of illustration and do not stop the algorithm after the second enrichment. The magenta line denotes the calculation time versus the number of enrichments. The black dashed line presents the calculation time that is needed when using the conventional Newmark algorithm. For comparison, all other calculation times are normalized with respect to that value. Using four or less enrichments, the proposed method is faster than the classical Newmark method. Increasing the number of enrichments, the calculation time grows linearly with the number of enrichments. Regarding the example discussed in Sect. 4.3, a speed-up by a factor of 2.08 is observed. Additionally, the results regarding computational efficiency are presented in Table 2.

## 5 Conclusion

In this paper, a Newmark space-time formulation for linear transient structural analysis is proposed.

The proper generalized decomposition is applied to decompose the incremental formulations into their space-time equivalent. Instead of proceeding incrementally in time, the system is solved in the whole time domain in terms of enrichment iterations. Within every enrichment iteration a system of two coupled space-time equations is solved using the fixed point algorithm. During convergence, the spatial and temporal enrichment functions are investigated.

The new method is demonstrated, using four different numerical examples. A frame structure is excited by a sinusoidal ground excitation at the first eigenfrequency, a second sinusoidal ground excitation that is larger than the first eigenfrequency, a transient ground excitation and an impact force function that is applied locally.

Exciting the structure by a sinusoidal ground acceleration with the excitation frequency equal to the first eigenfrequency leads to exactly one enrichment step until convergence is achieved. In the special case of resonance, the normalized spatial enrichment vector is equivalent to the normalized first linear mode of vibration and the shape of the temporal enrichment vector follows the response time history of the reference solution.

The structure is then excited by a sinusoidal ground excitation with an excitation frequency that lies between the first and the second eigenfrequency as well as by a transient ground excitation. In both cases the strategy leads to convergence already after two enrichments and the spatial enrichment functions show a high similarity with the corresponding low-frequency eigenmodes.

Finally, the structure is excited by a local impact-type force function. In this case, the first spatial enrichment modes do not reveal a high correlation with the first few modes of vibration. Instead deformation patterns, containing distinct localized shapes, are observed regarding the spatial modes describing both the local and global structural response behavior caused by the impact force function. This indicates that the proposed strategy is superior to classical strategies like the method of modal truncation, as it reveals an automatic adaption of the response patterns in form of an appropriate choice of spatial and temporal enrichment modes. In this regard, the method also outperforms the proper orthogonal decomposition method, as the performance of the latter depends on an a priori system identification for the generation of an appropriate reduced order basis [10,11].

Generally, the computational efficiency of the proposed strategy depends on the number of enrichments necessary to obtain a required level of accuracy. The numerical demonstrations in this paper show that a computational speed-up factor up to a value of four compared to the conventional Newmark method can be achieved with a high level of accuracy regarding problems in structural dynamics. While the conventional Newmark integration scheme proceeds in a step-by-step manner in time and is, therefore, only realiz-
able by a considerably high number of calculation steps, it is shown in this paper that the Newmark space-time strategy needs only a surprisingly small number of enrichment iterations until convergence is achieved.

The proposed strategy will serve as a guidance for future space-time approaches in nonlinear structural dynamics. The next step should be an extension of the proposed approach to nonlinear structures. It is essential to extend the framework to geometrically nonlinear problems as, e.g., the consideration of large deformation strain-displacement relationships or contact constraints, as well as problems considering materially nonlinear, history-dependent material behavior. Finally, it would be of high interest to establish the space-time formulations of other integration schemes and answer crucial research questions regarding applicability and efficiency. In particular, the elaboration of a space-time equivalent of the popular central difference integration scheme [2] that reveals only conditional stability in its classical formulation could provide a new perception of explicit integrators in structural dynamics.

**Acknowledgements** This research was financially supported by Deutsche Forschungsgemeinschaft (DFG) and the Agence Nationale de la Recherche (ANR) with project number MA2233/10-1.

**Funding** Open Access funding enabled and organized by Projekt DEAL.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

## A Multiplication rule of decomposed rank-three tensors tensors with decomposed rank-two tensors

In terms of components for a fix enrichment number $j$ we get:

$$K_{kli} = K_{kl}(t_i) = (K_s)_{kl}(k_t)_i. \qquad (102)$$

In the following, we need to marry the tensor products for the displacement and the stiffness with the contracted product

between stiffness and displacement. In index notation we get:

$$\sum_{l=1}^{n} K_{kl}(t_i) u_l(t_i) = \sum_{l=1}^{n} (K_s)_{kl}(k_t)_i (u_s)_l(u_t)_i \qquad (103)$$

$$\sum_{l=1}^{n} K_{kl}(t_i) u_l(t_i) = \sum_{l=1}^{n} (K_s)_{kl}(u_s)_l (k_t)_i (u_t)_i \qquad (104)$$

$$\sum_{l=1}^{n} K_{kl}(t_i) u_l(t_i) = (K_s u_s)_k (k_t u_t)_i \qquad (105)$$

or in abstract notations:

$$\mathbf{K} \mathbf{U} = (\mathbf{K}_s \mathbf{u}_s) \otimes (\mathbf{k}_t \circ \mathbf{u}_t). \qquad (106)$$

Combining the relations (42) and (48) allows to define the rule:

$$(\mathbf{K}_s \otimes \mathbf{k}_t)(\mathbf{u}_s \otimes \mathbf{u}_t) = (\mathbf{K}_s \mathbf{u}_s) \otimes (\mathbf{k}_t \circ \mathbf{u}_t). \qquad (107)$$

## B Newmark coefficients in space-time

The four Newmark space-time coefficients are $n_t \times n_t$ matrices. They are defined as:

$$A_1 = \frac{\beta}{\Delta t} \begin{bmatrix} 1 & 0 & 0 & 0 & \cdots & 0 \\ -1 & 1 & 0 & 0 & \cdots & 0 \\ 0 & -1 & 1 & 0 & \cdots & 0 \\ \vdots & \ddots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & 0 & -1 & 1 & 0 \\ 0 & \cdots & 0 & 0 & -1 & 1 \end{bmatrix}, \qquad (108)$$

$$A_2 = \beta \begin{bmatrix} 0 & 0 & 0 & 0 & \cdots & 0 \\ -1 & 0 & 0 & 0 & \cdots & 0 \\ 0 & -1 & 0 & 0 & \cdots & 0 \\ \vdots & \ddots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & 0 & -1 & 0 & 0 \\ 0 & \cdots & 0 & 0 & -1 & 0 \end{bmatrix}, \qquad (109)$$

$$A_3 = \begin{bmatrix} -\beta & 0 & 0 & 0 & \cdots & 0 \\ \beta - \frac{1}{2} & -\beta & 0 & 0 & \cdots & 0 \\ 0 & \beta - \frac{1}{2} & -\beta & 0 & \cdots & 0 \\ \vdots & \ddots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & 0 & \beta - \frac{1}{2} & -\beta & 0 \\ 0 & \cdots & 0 & 0 & \beta - \frac{1}{2} & -\beta \end{bmatrix}, \qquad (110)$$

$$A_4 = \begin{bmatrix} -\gamma & 0 & 0 & 0 & \cdots & 0 \\ \gamma - 1 & -\gamma & 0 & 0 & \cdots & 0 \\ 0 & \gamma - 1 & -\gamma & 0 & \cdots & 0 \\ \vdots & \ddots & \ddots & \ddots & \ddots & \vdots \\ 0 & \cdots & 0 & \gamma - 1 & -\gamma & 0 \\ 0 & \cdots & 0 & 0 & \gamma - 1 & -\gamma \end{bmatrix}. \qquad (111)$$
C Space-time initial conditions

The space-time equivalents of the initial conditions are
defined by four vectors **a**<sub>5</sub> - **a**<sub>8</sub> that are all of the dimen-
sion *n*<sub>*t*</sub>. Only the first component of each of the four vectors
is not equal to zero:

$$
\mathbf{a}_5 = \begin{bmatrix} \frac{\beta}{\Delta t} \\ 0 \\ 0 \\ \vdots \\ 0 \\ 0 \end{bmatrix}, \quad \mathbf{a}_6 = \begin{bmatrix} \beta \\ 0 \\ 0 \\ \vdots \\ 0 \\ 0 \end{bmatrix}, \\
\mathbf{a}_7 = \begin{bmatrix} \frac{1}{2} - \beta \\ 0 \\ 0 \\ \vdots \\ 0 \\ 0 \end{bmatrix}, \quad \mathbf{a}_8 = \begin{bmatrix} \beta(1-\gamma) \\ 0 \\ 0 \\ \vdots \\ 0 \\ 0 \end{bmatrix}. \qquad (112)
$$

References

1. Chopra AK (2007) Dynamics of structures. Pearson Education, New Jersey

2. Dokainish M, Subbaraj K (1989) A survey of direct time-integration methods in computational structural dynamics-I. Explicit Methods Comput Struct 32:1371–1386

3. Subbaraj K, Dokainish M (1989) A survey of direct time-integration methods in computational structural dynamics-II. Implicit Methods Comput Struct 32:1387–1401

4. Qu ZQ (2004) Model order reduction techniques with application in finite element analysis. Springer, London

5. Ding C, Deokar R, Cui X, Li G, Cai Y, Tamma K (2019) Proper orthogonal decomposition and Monte Carlo based iso-geometric stochastic method for material, geometric and force multi-dimensional uncertainties. Comput Mech 63:521–533

6. Im S, Kim E, Cho M (2019) Reduction process based on proper orthogonal decomposition for dual formulation of dynamic sub-structures. Comput Mech 64:1237–1257

7. Ostrowski Z, Białecki R, Kassab A (2005) Estimation of constant thermal conductivity by use of proper orthogonal decomposition. Comput Mech 37:52–59

8. Kerschen G, Golivani J (2002) Physical interpretation of proper orthogonal modes using the singular value decomposition. J Sound Vib 249:849–865

9. Rega G, Troga H (2005) Dimension reduction of dynamical systems: methods, models applications. Nonlinear Dyn 41:1–15

10. Bamer F, Bucher C (2012) Application of the proper orthogonal decomposition for linear and nonlinear structures under transient excitation. Acta Mech 223:2549–2563

11. Bamer F, Amiri AK, Bucher C (2017) A new order reduction strategy adapted to nonlinear problems in earthquake engineering. Earthq Eng Struct Dyn 46:537–559

12. Bamer F, Markert B (2017) An efficient response identification strategy for nonlinear structures subject to non-stationary generated seismic excitations. Mech Des Struct Mach 45:313–330

13. Craig D, Bampton M (1968) Coupling of substructures for dynamics analysis. AIAA J 6:1313–1319

14. Rixen D (2004) A dual Craig-Bampton method for dynamic sub-
structuring. J Comput Appl Math 168:383–391

15. Bamer F, Shi J, Markert B (2018) Efficient solution of the multiple seismic pounding problem using hierarchical substructure techniques. Comput Mech 62:761–782

16. Shi J, Bamer F, Markert B (2018) A structural pounding formulation using systematic modal truncation. J Shock Vib 54:1–15

17. Shi J, Bamer F, Markert F (2019) A substructure formulation for the earthquake-induced nonlinear structural pounding problem. Earthq Struct 17:101–113

18. Zhou L, Simon JW, Reese S (2018) Proper orthogonal decomposition for substructures in nonlinear finite element analysis: coupling by means of tied contact. Arch Appl Mech 88:1975–2001

19. Sanz-Serna JM (1983) On finite elements simultaneously in space and time. Int J Numer Methods Eng 19:623–624

20. Bajer CL (1986) Triangular and tetrahedral space-time finite elements in vibration analysis. Int J Numer Methods Eng 23:2031–2048

21. Bajer CI (1987) Notes on the stability of non-rectangular space-time finite elements. Int J Numer Methods Eng 24:1721–1739

22. Hughes TJR, Hulbert GM (1988) Space-time finite element methods for elastodynamics: formulation and error estimates. Comput Methods Appl Mech Eng 66:339–363

23. Li XD, Wiberg NE (1998) Implementation and adaptivity of a space-time finite element method for structural dynamics. Comput Methods Appl Mech Eng 156:211–229

24. Idesman AV (2007) Solution of linear elastodynamics problems with space-time finite elements on structured and unstructured meshes. Comput Methods Appl Mech Eng 196:1787–1815

25. Huang H, Constanzo F (2002) On the use of space-time finite elements in the solution of elasto-dynamic problems with strain discontinuities. Comput Methods Appl Mech Eng 156:211–229

26. Pattillo PD II, Tortorelli DA (2004) A contact algorithm for the Signorini problem using space-time finite elements. Int J Numer Methods Eng 60:1197–1213

27. Blum H, Jansen T, Rademacher A, Wienert K (2008) A contact algorithm for the Signorini problem using space-time finite elements. Int J Numer Methods Eng 76:1632–1644

28. Gimperlein H, Stocek J (2014) Space-time adaptive finite elements for nonlocal parabolic variational inequalities. Comput Methods Appl Mech Eng 352:137–171

29. Yang Y, Chirputkar S, Alpert DN, Eason T, Spottswood S, Qian D (2012) Enriched space-time finite element method: a new paradigm for multiscaling from elastodynamics to molecular dynamics. Int J Numer Methods Eng 92:115–140

30. Zhang R, Wen L, Naboulsi S, Eason T, Vasudevan VK, Qian D (2016) Accelerated multiscale space-time finite element simulation and application to high cycle fatigue life prediction. Comput Mech 58:329–349

31. Pasenow F, Zilian A, Dinkler D (2013) Extended space-time finite elements for landslide dynamics. Int J Numer Methods Eng 92:115–140

32. Takizawa K, Moormann C, Wright S, Christopher J, Tezduyar TE (2010) Wall shear stress calculations in space-time finite element computation of arterial fluid-structure interactions. Comput Mech 46:31–41

33. Giorla AB, Scrivener KL, Dunant CF (2014) Finite elements in space and time for the analysis of generalised disco-elastic materials. Int J Numer Methods Eng 92:115–140

34. Buch M, Idesman A, Niekamp R, Stein E (1999) Finite elements in space and time for parallel computing of viscoelastic deformation. Comput Mech 24:386–395
35. Idesman A, Niekamp R, Stein E (2001) Finite elements in space and time for generalized viscoelastic Maxwell model. Comput Mech 27:49–60

36. Zhang R, Wen L, Xiao J, Qian D (2019) An efficient solution algorithm for space-time finite element method. Comput Mech 63:455–470

37. Chinesta F, Ladevèze P, Cueto E (2010) A short review on model order reduction based on proper generalized decomposition. Arch Comput Methods Eng 18:395–404

38. Ammar A (2010) The proper generalized decomposition: a powerful tool for model reduction. Int J Mater Form 199:1872–1880

39. Chinesta F, Ladevèze P (2014) Separated representations and PGD-based model reduction, fundamentals and applications. Springer, London

40. Ladevèze P, Chamoin L (2011) On the verification of model reduction methods based on the proper generalized decomposition. Comput Methods Appl Mech Eng 200:2033–2047

41. Ladevèze P (1999) Nonlinear computational structural mechanics. New approaches and non-incremental methods of calculation. Springer, Berlin

42. Pruliere E, Chinesta F, Ammar A (2010) On the deterministic solution of multidimensional parametric models using the proper generalized decomposition. Math Comput Simul 81:791–810

43. Ammar A, Chinesta F, Cueto E, Doblaré M (2011) Proper generalized decomposition of time multiscale models. Int J Numer Methods Eng 90:569–596

44. Chinesta F, Ammar A, Cueto E (2012) On the use of proper generalized decompositions for solving the multidimensional chemical master equation. Eur Comput Mech 19:53–64

45. Aguardo J, Huerta A, Chinesta F, Cueto E (2014) Real-time monitoring of thermal processes by reduced-order modeling. Int J Numer Methods Eng 102:991–1017

46. Ammar A, Mokdad B, Chinesta F, Keunings R (2006) A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids: part I. J Non-Newtonian Fluid Mech 139:153–176

47. Ammar A, Mokdad B, Chinesta F, Keunings R (2007) A new family of solvers for some classes of multidimensional partial differential equations encountered in kinetic theory modeling of complex fluids: part II: transient simulation using space-time separated representations. J Non-Newtonian Fluid Mech 144:98–121

48. Dumon A, Allery C, Ammar A (2015) Proper generalized decomposition method for incompressible Navier-Stokes equations with a spectral discretization. Appl Math Comput 219:8145–8162

49. Signorini M, Zlotnik S, Díez P (2017) Proper generalized decomposition solution of the parameterized Helmholtz problem: application to inverse geophysical problems. Int J Numer Methods Eng 109:1085–1102

50. Sibileau A, García-González A, Auricchio F, Morganti S, Díez P (2018) Explicit parametric solutions of lattice structures with proper generalized decomposition (PGD). Comput Mech 62:871–891

51. González D, Cueto E, Chinesta F (2014) Real time direct integration of reduced solid dynamics equations. Int J Numer Methods Eng 99:633–653

52. Chinesta F, Leygue A, Bordeu F, Aguado V, Cueto E, Gonzalez D, Alfaro A, Ammar A, Huerta A (2013) PGD based computational vademecum for efficient design, optimization and control. Arch Comput Methods Eng 20:31–59

53. Ammar A, Chinesta F, Díez P, Huerta A (2010) An error estimator for separated representations of highly multidimensional models. Comput Methods Appl Mech Eng 3:88–102

54. Garikapati H, Zlotnik S, Díez P, Verhoosel C, van Brummelen H (2020) A proper generalized decomposition (PGD) approach to crack propagation in brittle materials: with application to random field material properties. Comput Mech 65:451–473

55. Blanco R, Borzacchiello D, Chinesta F, Díez P (2016) Monitoring a PGD solver for parametric power flow problems with goal oriented error assessment. Int J Numer Methods Eng 109:1–28

56. Heyberger C, Boucard P, Néron D (2012) Multiparametric analysis within the proper generalized decomposition framework. Comput Mech 49:277–289

57. Nadal E, Leygue A, Chinesta F, Beringhier F, Ródenas J, Fuenmayor F (2015) A separated representation of an error indicator for the mesh refinement process under the proper generalized decomposition framework. Comput Mech 55:251–266

58. Bergheau M, Zuchiatti S, Roux J, Feulvarch S, Tissot S, Perrin G (2016) The proper generalized decomposition as a space-time integrator for elastoplastic problems. C R Mecanique 344:759–768

59. Nasri M, Robert C, Ammar A, Arem S, Morel F (2018) Proper generalized decomposition for the numerical simulation of polycrystalline aggregates under cyclic loading. C R Mecanique 346:132–151

60. Boucinha L, Gravouil A, Ammar A (2013) Space-time proper generalized decomposition for the resolution of transient elastodynamic models. Comput Methods Appl Mech Eng 255:67–88

61. Shirafkan N, Bamer F, Stoffel M, Markert B (2020) Quasistatic analysis of elastoplastic structures by the proper generalized decomposition in a space-time approach. Mech Res Commun 10(1016/2020):103500

62. Chevreuil M, Nouy A (2011) Model order reduction based on proper generalized decomposition for the propagation of uncertainties in structural dynamics. Int J Numer Methods Eng 89:241–268

63. Quaranta G, Martin C, Ibanez R, Duval J, Cueto E, Chinesta F (2019) From linear to nonlinear PGD-based parametric structural dynamics. C R Mecanique 347:445–454

64. Wriggers P, Zavarise G (2004) Computational contact mechanics. In: Stein E, Borst R, Hughes JTR (eds) Encyclopedia of computational mechanics

65. Bathe KJ (2014) Finite element procedures. Watertown, MA

66. Newmark JL (1979) A method of computation for structural dynamics. J Eng Mech Div 85:67–94

67. Daniel WJT (1997) The subcycled Newmark algorithm. Comput Mech 20:272–281

68. Zupan E, Saje M, Zupan D (2013) Dynamics of spatial beams in quaternion description based on the Newmark integration scheme. Comput Mech 51:47–64

69. Pasetto M, Waisman H, Chen JS (2019) A waveform relaxation Newmark method for structural dynamics problems. Comput Mech 63:1223–1242

70. Bamer F, Markert B (2018) A nonlinear visco-elastoplastic model for structural pounding. Earthq Eng Struct Dyn 47:2490–2495

71. Bamer F (2018) A Hertz-pounding formulation with a nonlinear damping and a dry friction element. Acta Mech 229:4485–4494

72. Chatterjee A (2000) An introduction to the proper orthogonal decomposition. Curr Sci 78:808–817

73. Kobe Takarazuka earthquake measurement 1995-01-16 (UTC: 20:46:52), ath.KOBE.TAZ000.AT2. (2016) Strong Motion Virtual Data Center (VDC). www.strongmotioncenter.org

74. Falsone G, Settineri D (2011) An Euler–Bernoulli-like finite element method for Timoshenko beams. Mech Res Commun 38:12–16

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.