Non-incremental response evaluation in geometrically nonlinear structural dynamics using a space-time stiffness operator

This contribution presents a proper generalized decomposition-based nonlinear solver for an efficient solution of geometrically nonlinear dynamic problems. The solution is built as a sum of dyadic products of space and time modes, and this sum of so-called enrichments is truncated when the required accuracy is achieved. In the proposed algorithm, we apply a consistent linearization of the residual vectors around the currently known solution over the whole space-time domain. At first, the set of vectorized tangent stiffness matrices is separated in space and time using the singular value decomposition. Then, the left and right singular vectors are reshaped into matrices to separate the space-time stiffness operator. The latter can be incorporated into the alternating fixed-point algorithm to compute couples of space and time modes. Numerical examples of a two-dimensional geometrically exact beam model demonstrate the accuracy, efficiency, and limits of the method.


Introduction
An extensive range of numerical analysis tools can predict the dynamical behavior of complex structures and mechanical systems. The associated governing equations of motion are commonly discretized in space with the finite element method [1] and in time using time integration schemes [2,3]. In order to solve these algebraic equations, the classical, incremental, step-by-step approach is widely used in linear and nonlinear numerical analysis. Geometrical nonlinearity is required to model structures undergoing large displacements and rotations, which is a relevant feature in many engineering applications. However, solving the corresponding nonlinear equations of motion is computationally demanding, especially when the transient response of the B Franz Bamer bamer@iam.rwth-aachen.de Tahar Arjoune arjoune@iam.rwth-aachen.de Bernd Markert markert@iam.rwth-aachen.de 1 Institute of General Mechanics, RWTH Aachen University, Eilfschornsteinstraße 18, Aachen 52062, North Rhine-Westphalia, Germany system over a large time interval is needed since an iterative Newton-Raphson procedure [4] with several residual and Jacobian evaluations is required within every calculation step. In order to find sufficiently accurate solutions with significantly lower computational costs, several approximation methods and model order reduction strategies have been proposed in the literature. They can be divided into methods that require offline computations to achieve online speedups and methods that do not need any offline step. Projection-based model order reduction requires the execution of offline computations to define a reduction basis before performing more efficient online computations. It aims to decrease the computational effort by considering a reduced number of unknowns, which can describe almost all relevant configurations of the structure for the problem at hand. This smaller subspace, where the solution is assumed to exist, is spanned by a set of basis vectors, which define possible displacement shapes and form a reduction basis. In linear structural dynamics, an appropriate reduction basis can be identified with the modal truncation method, which goes back to lord Rayleigh [5]. It approximates the dynamical system by a superposition of a selected set of vibration modes. The Krylov Subspace method [6] exploits the information about the location of external excitation points to build the reduction basis from static displacement fields and higher-order approximations. Lately, Sonneville et al. [7] presented an extensive review of component mode synthesis methods a comprehensive review of modal reduction procedures for flexible multibody dynamics introducing a general approach for the modal reduction of geometrically nonlinear problems. The component mode synthesis and substructuring techniques [8][9][10][11][12][13] enable the definition of individual reduction bases in different subdomains of a partitioned structure. These methods enable the assembly of the reduced substructures by keeping the degrees of freedom of the boundaries unaltered.
Although the methods mentioned above are very well applicable to linear systems, they perform very poorly for geometrically nonlinear structures. Idelson and Cardona [14] propose the augmentation of the linear reduction basis with modal derivatives in order to capture the nonlinear behavior of the system. The computation of these additional basis vectors is based on the perturbation of the eigenvalue problem around the equilibrium position. This method is extensively used for the reduction of geometrically nonlinear structures [15][16][17][18][19][20] since it exploits intrinsic physical properties of the system and is less sensitive to changes in model parameters, boundary conditions, or excitations in comparison with data-driven approaches [21]. The most prominent data-driven approach, the proper orthogonal decomposition, also known as Karhunen-Loève decomposition, identifies the optimal subspace from a set of displacement snapshots, which are computed using training simulations of the full, high dimensional system [22,23]. The basis vectors can be extracted from the set of snapshots by applying the singular value decomposition. The proper orthogonal decomposition is widespread in nonlinear structural dynamics [24][25][26][27][28] and fluid mechanics [29,30] as well as in many other fields, as it is independent of the considered system and can be applied to any model as long as representative training simulations are available.
The projection of the nonlinear equations of motion onto a smaller subspace spanned by a meticulously chosen reduction basis reduces the required computational effort by reducing the linear systems of equations that need to be solved in every Newton-Raphson iteration. However, the evaluation of the nonlinear internal force vectors and tangent stiffness matrices of the reduced system is still related to the size of the high dimensional system because the internal forces in geometrically nonlinear finite elements are formulated on the element level. In order to overcome this computational bottleneck and fully make use of the dimensionality reduction, hyperreduction techniques were proposed in the literature. The polynomial tensors method [31,32] exploits the polynomial structure of the internal forces of geometrically nonlinear structures. Computational savings are related to the cheaper evaluation of a multidimensional polynomial with fewer variables and precomputed coefficients. The energy-conserving mesh sampling and weighting method [33] reduces the computational costs by evaluating only a subset of finite elements and weighting their contributions by suitably chosen weights during the assembly process. The discrete empirical interpolation method [34] approximates the high dimensional nonlinear force vector with an empirical force basis using the proper orthogonal decomposition in a similar manner to the approximation of displacements.
The a priori identification of a suitable subspace for the application of projection-based reduction is, in general, a tedious and problem-specific task, in particular in the nonlinear regime. It requires the execution of offline computations and a certain familiarity with the applied reduction methods. Note that the smaller set of nonlinear equations obtained with this method is solved with the same incremental Newton-Raphson algorithm used for the full model. On the contrary, the proper generalized decomposition (PGD) method is a non-incremental solution procedure, which does not require the execution of any problem-specific offline computations before being able to achieve speedups in online simulations [35][36][37][38]. Using an iterative process, it approximates the entire solution over the space-time domain as a sum of dyadic products of space and time mode couples. Therefore, it leads to a space-time separated representation of the solution. A so-called enrichment couple composed of a space mode and its corresponding time mode is computed in each iteration and added to the sum of the previous ones until the required accuracy is achieved. Different strategies for the computation of these space and time modes are presented by Nouy [39]. The main advantage of PGD compared to projection-based methods is that the relevant subspace for the problem is identified during the online computation. This non-incremental solution procedure was first introduced by Ladevéze [40,41] for efficiently solving nonlinear problems in structural mechanics. The PGD was then implemented for the reduction of multidimensional problems [42,43] and multiparametric models [74]. In addition, It was used for solving crack propagation problems in brittle materials [75] and lattice structures problems [76]. It was also applied in many other fields such as fluid dynamics [44][45][46][47][48], heat transfer [49,50] and chemistry [51].
The equations of motion derived from structural dynamics problems contain damping and inertia, dependent on velocities and accelerations. A discretization scheme is hence necessary for the approximation of these time derivatives. Implicit time integration schemes, such as the Newmark [52] scheme, are very popular in the field of structural dynamics. They are usually employed for incremental solution procedures so that the displacements, velocities, and accelerations at the current time step can be computed using the solution of the previous time step. Boucinha et al. proposed a tensorial formulation for the integration of the Newmark scheme into the non-incremental multi-field PGD solution procedure [53], where both displacements and velocities are defined as primary variables. Bamer et al. extended the application of this tensorial formulation of the Newmark time discretization scheme to the single-field PGD solution procedure, where the primary variables of the equations of motion are only the displacements [54].
The proper generalized decomposition was efficiently implemented for solving several nonlinear dynamic and static problems, where the nonlinear terms were handled with various strategies. The most straightforward and primarily used alternative is the fixed point Picard iteration [38]. It was implemented within the PGD framework for solving nonlinear problems in computational rheology [55], multidimensional parametric modeling [42], magneto-thermal coupling [56] and heat transfer [49,57]. The fixed point Picard iteration approximates the nonlinear terms by partially substituting the unknown variables with the solution computed in the previous PGD iteration. This approach can only be applied to nonlinear equations with a specific structure [58]. The nonlinear algebraic equations describing geometrically nonlinear models do not generally satisfy this condition and hence can not be solved with the fixed point Picard technique. The asymptotic numerical method proposed in [59] was implemented within the PGD framework in order to approximate the nonlinear expressions in terms of an expansion parameter around the equilibrium configuration [60,61]. Although this efficient method leads to a solution of a sequence of linear problems with the same constant tangent stiffness matrix, it only delivers a good approximation of the solution in the neighborhood of the equilibrium configuration [62]. In the context of the quasistatic analysis of problems with material nonlinearities, the linearization of the internal forces using a unique constant tangent stiffness matrix at the equilibrium configuration in an analogous manner to the modified Newton method was successfully implemented in a PGD framework [63,64]. However, while quasistatic materially nonlinear problems have successfully and efficiently been solved using the PGD approach it has yet not been possible to provide appropriate PGD algorithms that ensure the convergence in the geometrically nonlinear case.
Thus, this paper presents a new space-time formulation for geometrically nonlinear static and dynamic problems in structural dynamics using the PGD method. In doing so, this paper proposes a consistent linearization of the nonlinear algebraic equations over the whole space-time domain within a progressive PGD framework. For this purpose, the residual vector, evaluated at every time step, is linearized around its corresponding currently known solution by evaluating the first two terms of a Taylor expansion analogously to the Newton-Raphson procedure. The evaluated set of tangent stiffness matrices is then separated in space and time using the singular value decomposition to be incorporated into the proposed PGD based nonlinear solver. This valuable Jaco-bian information ensures the convergence of the space-time solution of a broader range of nonlinear problems, amongst others geometrical nonlinearity.
We organize the paper as follows: In Sect. 2, the semidiscretized equation of motion, describing the dynamical behavior of a geometrically nonlinear hyperelastic structure, is introduced. Subsequently, the discretization in time and the solution procedure of the obtained algebraic system of equations by the incremental Newton-Raphson method are presented. In Sect. 3, the fundamentals of the proper generalized decomposition method and the embedding of the Newmark time integration scheme are described for the linear case. Building on this, the proposed solution procedure for the nonlinear case is introduced. In Sect. 4, the suggested algorithm is demonstrated on a static and a dynamic numerical example of a two-dimensional geometrically exact beam model, showing the accuracy, efficiency, and limits of the method. Finally, conclusions are drawn, and further discussions are presented in Sect. 5.

Preliminaries
The spatial discretization of the equation of motion of a geometrically nonlinear hyperelastic structure with the finite element method is written as: with the vector of nodal displacements u(t) ∈ R n s , the mass matrix M ∈ R n s ×n s , the internal restoring force vector f int (u(t)) ∈ R n s and the external force vector f ext (t) ∈ R n s , where n s is the number of degrees of freedom. The explicit dependence of u and f ext on time will be henceforth omitted for the reason of abbreviation. The internal restoring force is a nonlinear function of the nodal displacements, since a nonlinear strain measure is required for the description of large displacements and finite rotations. A more detailed description of the two-dimensional geometrically exact beam element used in this paper is presented in appendix A. We introduce viscous damping using the Rayleigh assumption: This assumption is a weighted sum of the mass matrix and the initial tangent stiffness matrix evaluated at zero displacement, i.e., u = 0 [65]. Although this purely phenomenological damping approximation is based on the linearized model, it is often used in the nonlinear regime leading to the following damped nonlinear equations of motion:

Time discretization
The semi-discretized equations of motion (1) and their damped counterpart (3) are ordinary differential equations stemming from the spatial discretization by the finite element method. These equations are solved based on specified initial displacements u(t = 0) and velocitiesu(t = 0) in the time interval t ∈ [t 0 , t end ]. A temporal discretization with a time integration scheme is required in order to transform the ordinary differential equations to a set of algebraic equations, which can be solved with established algorithms. This means, that the solution u is defined on a set of discrete time instants P = {t 0 , t 1 , . . . , t n t = t end } and in between these time instants one considers polynomial approximations. The implicit one-step time integration scheme proposed by Newmark [52] is one the most popular schemes in the realm of structural dynamics since it makes use of the second-order structure of the equations of motion (3). The approximation of the velocitiesu j+1 and displacements u j+1 at the time step t j+1 is expressed in terms of the unknown accelerationsü j+1 and all quantities of the previous time step as: where t is the time step size and the parameters γ ∈ [0, 1] and β ∈ 0, 1 2 represent the weighting factors of the previous and the current time steps, which determine the stability and accuracy of the time integration scheme. The rearrangement of the Eqs. (4) and (5) yields the approximation of the velocitiesu j+1 and accelerationsü j+1 of the time step t j+1 in terms of the unknown displacements u j+1 and all quantities of the previous time step: Inserting Eqs. (6) and (7) into the semi-discretized equations (3) leads to the nonlinear algebraic balance equations: where j = 0, . . . , n t − 1. They are solved for the unknown displacements u j+1 for a number n t of time steps.

Incremental solution procedure
The nonlinear system of equations (8) are usually solved incrementally, where the solution u j+1 is computed iteratively using a Newton-Raphson solver based on the converged solution u j of the previous time step [65]. This iterative process within each time step calculation tries to find the root of the residual function of the system of equations (8): In each iteration i of the Newton-Raphson loop, a linearization of the residual around the current solution i u j+1 is evaluated and the linear system of equations: is solved in order to compute the correction i u j+1 yielding the updated solution for the next iteration step: The steps (10) and (11) are repeated until convergence is achieved. A commonly used convergence criterion is the decrease of the L 2 -norm of the residual vector below a certain tolerance, which must be chosen carefully. The linearization of the residual (9) around a certain displacement i u j+1 is given as: where K ( i u j+1 ) is the tangent stiffness matrix (A23), evaluated for the corresponding configuration of the structure. In the case of static analysis, the same aforementioned procedure is adopted to find a solution satisfying the equilibrium: where only internal and external forces play a role. Although there is no time discretization involved in this static problem, the external forces are incrementally increased, and the solution is computed in a step-wise manner to ensure the convergence of the Newton-Raphson nonlinear solver.

The proper generalized decomposition
In this section, the non-incremental solution procedure, using the proper generalized decomposition (PGD), is introduced. For the reason of convenience, we apply a tensorial formulation by means of numerical tensors of order 1, 2, 3 and 4. Using the PGD, the solution matrix U = u 1 · · · u n t ∈ R n s ×n t containing the n t solution vectors of the n t time steps of the discretized equations of motion (9) and (13) are approximated by the following rank-m separated representation: The approximation U (m) of the second-order tensor U is constructed by the sum of a number m of rank-one tensors, also called enrichments [53,66]. Each enrichment constitutes a tensor product of a space mode u (i) s ∈ R n s and a time mode u (i) t ∈ R n t . The size n s of a space mode corresponds to the number of degrees of freedom of the considered mechanical system, and the size n t of a time mode represents the number of time steps associated with the time discretization. Using the non-incremental PGD procedure, we are searching for a solution U (m) , for which the Frobenius norm of the matrix composed of the residual vectors of the equations of motion of all time steps is below a certain tolerance res : where If we assume that the decomposition U (m−1) of rank m −1 has been already computed in a previous step, then a new optimal enrichment couple (u where f p is a carefully chosen tolerance. Note, that the fixed point problem can also start with an arbitrary space mode 0 u (m) s as an initial guess and that a maximum number of iterations can be specified with the parameters k max . The normalization performed in line 6 of algorithm 1 in not absolutely necessary, however, it avoids very large or very small entries in the computed space and time modes [39,69].
The decomposition (14) can have different definitions and can be computed with several methods. An overview of these different possibilities is presented in [39]. In this paper, the progressive Galerkin-based PGD presented in algorithm 1 is used since it offers the best trade-off between robustness and computational costs and is, therefore, of high practical interest [39,69]. It projects the problem onto a subspace orthogonal to the residual in each iteration k, which leads to a time-independent space problem, when the new time mode is fixed, and to a scalar time problem, when the new space mode k u (m) s is fixed. Note that there is a more robust but computationally less efficient alternative to the progressive Galerkin-based PGD, the progressive, minimal residual PGD method. In this formulation, the new optimal enrichment couple (u (m) is computed in each enrichment iteration m by minimizing the residual norm: The main advantage of the minimal residual formulation in comparison with the Galerkin-based strategy is that its convergence in the residual norm is proven to be monotonic [39,69]. The so-called optimal PGD is also proposed in the literature as an alternative to the progressive construction of the space and time modes [39,53,69]. It leads to an optimal decomposition (14) of rank m in the sense of the Galerkin projection. This means that the residual (16) is simultaneously orthogonal to the subspace spanned by the space modes, s ]) ⊂ R n s and to the subspace spanned by the time modes, The optimal PGD is not used in this paper, as it turns out to be computationally inefficient for the problems discussed in this paper.

Solving linear dynamic problems with the PGD
The PGD solution procedure was extensively discussed in the literature for solving partial differential equations associated with different linear dynamic problems, such as the diffusion problem in [69] and the advection-diffusion-reaction problem in [39]. In the field of structural dynamics, PGD was implemented in [53] and [54] in order to solve the linear elastodynamic problem. When the governing partial differential equations are linear with respect to the variables, e.g., the displacements u, then the resulting systems of equations shown in lines 5 and 7 of algorithm 1 have to be solved in each iteration k of the fixed point loop. They are linear and can, therefore, be solved with standard linear solvers [53].
The residual vector of the discretized equations of motion associated with the linear elastodynamic problem is given as: where j = 0, . . . , n t − 1 and the tangent stiffness matrix around the initial configuration is given as follows: Since the matrices M, C and K 0 are constant and independent from the configuration of the structure, we can express the matrix composed of the residual vectors at all time steps as: where and In the following, we introduce a tensorial formulation proposed in [53] and implemented in [54], which enables to express the matrices of velocities and accelerations,Ü andU, as functions of the displacements U and the initial conditions, u 0 andu 0 , applying the Newmark scheme, presented in equations (6) and (7). First, equations (4) and (5) of the Newmark time integration scheme are rearranged as follows: The key idea of this tensorial formulation is to build a matrix of the following linear combination of all couples of two successive displacement vectors: with j = 0, . . . , n t − 1, as: where the coefficients c 1 , c 2 ∈ R depend on the parameters of the time integration scheme and U = u 1 · · · u n t ∈ R n s ×n t is the matrix of the solutions at all time steps. The coefficient matrix is defined as: and the coefficient vector as: The tensorial formulation (27) can be also applied to the velocity and acceleration vectors in order to express the equations (24) and (25) at all time steps as: which can be also written as: where the definitions (28) and (29) are used as follows: The initial acceleration is computed by solving the system of equations of motion for a given initial displacement and velocity: Inserting the time integration equations (32) and (33) in the residual matrix of the equations of motion (21), leads to a set of algebraic equations with only the displacements U as variables: In these equations, the matrices Y , W and L are defined as: and where The solution of the space problem, as shown in line 5 of algorithm 1 and associated with the discretized equations of motion (36), is computed by solving the linear system of equations E s k u (m) s = b s with the corresponding coefficient matrix: and the right hand side is written as: The solution of the corresponding time problem, as shown line 7 of algorithm 1, is computed by solving the linear system of equations E t k u (m) t = b t with the corresponding coefficient matrix: where I t ∈ R n t ×n t is the identity matrix, and the right hand side is written as:

Solving nonlinear static and dynamic problems with the PGD
In this section, we introduce an extension of the PGD solution procedure, presented in algorithm 1, to the nonlinear case by linearizing the equations of motion of each time step j at the corresponding solution u (m−1) j , which represents a column vector of the space-time solution U (m−1) , computed in the previous enrichment step. In this way, we can apply the fixed point iteration in each enrichment step m in a similar manner to the linear case, discussed in Sect. 3.2. The only difference is that in the nonlinear case, the linearization is not constant anymore and only evaluated once at the initial undeformed configuration, but it is dependent on the state of the structure.
In this paper, the PGD solution procedure for nonlinear structures is presented with a hyperelastic geometrically nonlinear finite element model, where the nonlinear term in the system of equations of motion (8) is the vector of internal forces f int (u). In this case, the matrix, composed of the residual vectors at all time steps, is written for the static problem as: and for the dynamic problem as: where The time discretization of the nonlinear equations of motion (45) is performed in the same manner as for the linear case (36) and is written using the tensorial formulation as: In order to compute a new enrichment couple (u , the nonlinear terms of the residual vectors, arranged in the matrices (44) and (47), have to be linearized at the corresponding configuration of the structure, which is a column vector of the space-time solution ] computed in the previous enrichment step. The matrix of linearized internal force vectors is given as: where the column vectors of the new enrichment matrix are written as: Inserting (48) in (47) yields the matrix of linearized residual vectors: which has to be projected alternately on a fixed space mode and on a fixed time mode in order to find a solution (49) for the fixed point problem. For this purpose, a space-time separated representation of the second term of the matrix of linearized internal force vectors (48) is formulated as follows: In this equation, the term n k l=1 K (l) s ⊗ K (l) t ∈ R n s ×n s ×n t ×n t is a separated representation of rank n k of the bilinear stiffness operator, which constitutes a tensor of fourth-order and defines a linear mapping of a second-order tensor into another one [53,66]. The main idea of the proposed PGD solution procedure in algorithm 2 for solving nonlinear problems, defined in (47) and (44) is the introduction of a consistent linearization of the residual vectors, as shown in equations (48) to (52). Hence, the space problem solved in line 6 of algorithm 2 and associated with the nonlinear dynamic case (47) yields the linear system of equations with the coefficient matrix: The right hand side is written as: Algorithm 2 PGD solution procedure for nonlinear problems The corresponding time problem, as shown in line 8 of algorithm 2, is equivalent to the linear system of equations with the coefficient matrix: The right hand side is written as: The space-time separated representation of the bilinear stiffness operator, as shown in equation (52) for the solution U (m−1) , can be computed for a certain space-time solution U = u 1 · · · u n t with the singular value decomposition based on the stiffness matrices of the structure at all time steps K (u 1 ) , . . . , K u n t . For this purpose, we first need to vectorize each stiffness matrix K u j = [k j,1 · · · k j,n s ] ∈ R n s ×n s as follows: with vec : R n s ×n s → R n 2 s and vec −1 : R n 2 s → R n s ×n s . The operator vec(·) creates one column vector from a matrix K by stacking its column vectors below one another. Then, the singular value decomposition of the matrix Z(U) ∈ R n 2 s ×n t , composed of the vectorized stiffness matrices of all time steps, yields: where the matrixÛ ∈ R n 2 s ×n t contains the left singular vectors, the diagonal matrix = diag(σ 1 , . . . , σ n t ) contains the decreasingly ordered singular values, and the orthogonal matrixV ∈ R n t ×n t contains the right singular vectors. Using the singular value decomposition (58), an optimal low rank approximation Z app with rank n k of the matrix of vectorized stiffness matrices Z(U) can be built as a sum of n k rank one matrices: where the vectorsû l andv l are the first n k column vectors of the matricesÛ andV , respectively, and are associated with the largest singular values σ l . The left and right singular vectors and the corresponding singular values used in the approximation (59) enable to build a low-rank approximation of the bilinear stiffness operator as follows: which is required for the implementation of algorithm 2. Note, that the inverse operator vec −1 (·) of the vectorization operator vec(·) transforms back a column vectorû l ∈ R n 2 s into a matrix K (l) s ∈ R n s ×n s and the operator diag(·) builds a diagonal matrix K (l) t ∈ R n t ×n t using the entries of a vector σ lvl ∈ R n t . The steps given in equations (57) to (60) for computing a space-time separated representation of the stiffness operator for a certain space-time solution U are graphically illustrated in Fig. 1. The computation time and memory requirements of these operations, especially the singular value decomposition, become prohibitive for models with a large number of degrees of freedom n s , because the number of rows of the matrix Z(U) increases quadratically with n s . Since the stiffness matrices of finite element models are not fully populated sparse matrices with the same sparsity pattern, we propose to only consider the nonzero entries of the stiffness matrices in the vectorization step (57). The number n nz of nonzero entries of a stiffness matrix K (u) ∈ R n s ×n s is, in general, significantly smaller than n 2 s , especially for large Fig. 2 Finite element model of a cantilever beam composed of 10 three-node geometrically exact beam elements. The translational and rotational degrees of freedom of the left-most node are fixed and a single moment is applied at the right-most node finite element models, which yields a considerably smaller matrix Z(U) ∈ R n nz ×n t and, therefore, a more efficient but equally valuable singular value decomposition step (58).

Numerical examples
In this section, we investigate the ability of the proposed PGD solver presented in algorithm 2 to approximate the solution of a nonlinear problem with a space-time separated representation. For this purpose, a two-dimensional finite element model of a cantilever beam is built using a three-node geometrically exact beam element, presented in Appendix A.
The cantilever beam has a length l = 1 m and a square cross-section with a side length of 0.02 m. Its finite element model is composed of 10 elements and 21 nodes, as illustrated in Fig. 2. Each node is associated with two translational and one rotational degree of freedom. The fixation of the translational and rotational degrees of freedom of the left-most node is performed by eliminating these variables and the equilibrium equations associated with them from the global system of equations so that the mechanical system contains 60 degrees of freedom in total. The single moment at the right-most node is applied by prescribing the value of the moment at the corresponding entry of the global vector of external forces f ext . The considered hyperelastic Saint Venant-Kirchhoff material model is characterized by a Young's modulus E = 210 GPa, a Poisson's ratio ν = 0.3 and a density ρ = 7800 kg/m 3 . In order to compare between different approximations U app of the solution, the following relative error measure with respect to a reference solution U re f is introduced: where the operator · defines the Frobenius norm of a matrix. The averaging with the number of time steps n t in the formula (61) of R E is necessary when simulations with a different number of time steps are compared since errors accumulate. In addition, the observation of the displacements of the right-most node of the model over time is provided to investigate the accuracy of the PGD solutions with respect to a reference solution.
Solving the equations of motion of the static problem in the nonlinear regime requires the same procedure as for the dynamic problem since the external forces should preferably change in sufficiently small increments to achieve convergence, which makes it also a time-dependent problem. The only difference is that inertia and damping forces drop out from the equilibrium equation in the static case, and the residual is composed only of internal and external forces. Therefore, we will first consider a geometrically nonlinear static problem in Sect. 4.1 and then we will investigate a geometrically nonlinear dynamic problem in Sect. 4.2. In both numerical demonstrations, we present the cantilever beam model, as illustrated in Fig. 2.

Geometrically nonlinear static problem
For the static problem, represented by the cantilever beam model shown in Fig. 2, a system of nonlinear equations of the form (13) has to be solved for a number of time steps with a stepwise increasing external moment. Within the PGD solver in algorithm 2, the norm of the matrix composed of the residual vectors of all time steps (44) has to be minimized by adding new enrichments to the space-time solution U. The magnitude of the external moment applied at the right-most node is increased linearly with respect to time according to the following increment: where M s defines the rate of increase of the external moment. The solution of the static problem, illustrated in Fig. 3 at different time points, is computed with an incremental Newton-Raphson solver for a time interval T = [0, 15] s with a time step t = 0.01 s and M s = 10 3 N · m · s −1 . It will be used as a reference solution in the assessment of the accuracy of the solutions computed with the PGD solver.

Convergence behavior and performance assessment
The numerical simulations discussed in this subsection are performed by the numerical setup presented in Table 2. The PGD solver presented in algorithm 2 is implemented to solve a static problem with n s = 60 degrees of freedom and n t = 150 time steps. The initial solution U (0) is assumed to be the undeformed configuration u = 0 for all time steps. The termination of the enrichment loop is only controlled by setting a maximum number of enrichments m max and not with a tolerance res since we are interested in analyzing the convergence behavior when new enrichments are added to the current solution. The fixed point loop is terminated either when the maximum number of iterations k max is reached, or when an absolute tolerance f p,abs and a relative tolerance     The configurations of the cantilever beam corresponding to the PGD solution with different numbers of enrichments and the Newton-Raphson reference solution are shown in Fig. 4a after the last external moment increment (n t = 150). The configuration PGD-1 in red of the cantilever beam represents the last column of the rank one PGD solution U (1) , which means that the solution includes only the first enrichment u (1) s ⊗ u (1) t . This first solution is nothing else than the solution of the linearized system around the initial undeformed configuration since the vectors of internal forces of all time steps are linearized around the initial solution U (0) = 0 in the first iteration m = 1. Note that, in this case, the spacetime separated representation of the stiffness operator is only composed of one summand K (u = 0)⊗ I n t , where I n t is the identity matrix of size n t . The first space and time modes are illustrated in the top of Fig.5a, b, respectively. The first time mode shows that the response, represented by the first enrichment, is proportional to the external moment and is, therefore, a linear response. The configuration PGD-2 in orange represents the last column of the rank two PGD solution U (2) and includes already the geometrically nonlinear behavior of the structure by considering the linearization of the internal force vector of each time step around the corresponding configuration within the previous solution matrix U (1) . The solution U (2) is the sum of the outer products of the two first space modes with their respective time modes, illustrated in Fig. 5. The same procedure is then repeated until a convergence criterion is achieved. Figure 4a shows that the PDG solution converges towards the Newton-Raphson reference solution presented by the black line when a larger number of enrichments is considered. The plot in Fig. 4b illustrates the time evolution of the displacement in x-direction of the rightmost node, which undergoes the largest displacements in the whole cantilever beam. It demonstrates that the nonlinear response of the structure in x-direction is already accurately represented after the second enrichment. Figure 5 represents the normalized space and time modes. It is shown that the contribution of new enrichments to the solution becomes less significant with increasing enrichment, while the temporal evolution becomes more complex.
The convergence behavior of the solution of the static problem computed with the proposed nonlinear PGD solver can also be shown quantitatively by observing the evolution of the relative error measure in Fig. 6a or the Frobenius norm of the residual matrix in Fig. 6b with respect to the number of enrichments. The relative error is computed in accordance with formula (61), where the PGD solution U (m) after the m-th enrichment is inserted as U app , and the highly accurate Newton-Raphson solution is denoted as U re f . The convergence of the residual norm of the solution, computed by algorithm 2, is, in general, non-monotonic because it is based on a Galerkin formulation, as discussed in Sect. 3.1 [39,69].
The total computation time needed by the PGD solver increases linearly with respect to the number of enrichments, as shown in the plot of Fig. 6a, although the number of fixedpoint iterations within each enrichment differs significantly. This phenomenon is caused by the fact that the evaluation of the residual vectors and the stiffness matrices at U (m−1) during the linearization step, which is shown in line 4 of algorithm 2 and depends only on the size of the finite element model and the number of time steps, take the major part of the computational effort (97.2%). In comparison, the space-time decomposition of the stiffness operator and the space and time problems solved within the fixed point iterations take a negligible part of the total computational time. For the numerical setup in Table 2, the plot 6a shows that the computation of the PGD solution with two enrichments (m max = 2) is approximately two times faster than the computation of the reference Newton-Raphson solution, which is represented by the horizontal dashed green line, while it exhibits a lower level of accuracy. The PGD solution based Fig. 7 Linear relation between the problem size (number of dofs) and the required time for the computation of five enrichments on three enrichments requires as much computation time as the Newton-Raphson solution, and starting from the fourth enrichment, it becomes slower. Therefore, the enhancement of the linearization step can reduce the slope of the green linear curve and pave the way for more efficient computations with the proposed nonlinear PGD solver.

Influence of the number of degrees of freedom
Since the number of degrees of freedom of our benchmark model is relatively small we investigate the efficiency of the proposed algorithm with an increasing number of degrees of freedom. The parameter setting of this study is presented in Table 2. As shown in this table, we increase the number of degrees of freedom of the cantilever beam problem from 60 to 1800. The computational time versus the number of degrees of freedom is depicted in Fig. 7. As shown in this figure, one observes a linear increase in computation time when increasing the number of degrees of freedom. We note that this linear relation can be achieved since we apply a vectorization operator on the sparse tangential stiffness matrices in Eq. (60) that considers only the non-zero entries. This approach reduces the memory requirements of the spacetime stiffness matrix Z(U) considerably and allows for a significantly less expensive decomposition of the space-time stiffness matrix Z(U) using the singular value decomposition, as presented in Eq. (58).

Influence of the number of time steps and the size of external force increments
In contrast to the incremental Newton-Raphson solver, the convergence behavior of a PGD solver depends on the number of time steps considered in the computation. In order to illustrate this dependency, we performed different computations of static problems with different numbers of time steps n t for a particular external moment increment size. The same computation was then repeated with external moment incre-  60 10 −2 20 ments M of different sizes. We present the numerical setup of the performed numerical experiments in Table 3.
The termination criterion of the PGD solver is only based on reaching a certain residual norm R < res , where the tolerance res is proportional to the number of time steps and, therefore, to the size of the matrix R. The termination criterion does, therefore, not depend on the considered number of time steps. The plots in Fig. 8 show that the computation time of PGD solutions increases exponentially with respect to the number of considered time steps, while it increases linearly for an incremental Newton-Raphson solver. This phenomenon is because the convergence of the PGD solver worsens with an increasing number of time steps and is even lost at a certain point, whereas a higher number of time steps does not influence the convergence behavior of an incremental Newton-Raphson solver. It should also be mentioned that when the size of the external moment increment is larger, the required computation time increases faster with an increasing number of time steps. The numerical results, depicted in Fig. 8a, indicate that the convergence is lost in all considered three cases when the sum of all external moment increments of the static problem reaches nearly the same value of M = 1.8 · 10 3 N · m. For this reason, the PGD solver in algorithm 2 can not be implemented for solving the static problem illustrated in Fig. 3 within the whole time domain at once using only one space-time decomposition with n t = 1500. Thus, in the next section, we propose a piecewise formulation of the PGD solver that operates in a sequence of smaller time intervals. Within each time interval, the time modes have a smaller size, and convergence can be ensured.

A piecewise PGD formulation
In this section, we present a piecewise PGD formulation for solving problems over large time intervals, such as the static (a) (b) Fig. 8 a Relative computation time of PGD solutions with different numbers of time steps with respect to the computation time required for ten time steps (n t = 10). The same numerical experiment is performed for different sizes of external moment increments. The curves are interrupted, since convergence was not achieved for a higher number of time steps with the corresponding external moment increment size; b Relative computation time of incremental Newton-Raphson solutions with different numbers of time steps with respect to the computation time required for ten time steps (n t = 10). The same numerical experiment is performed for different sizes of external moment increments Fig. 9 Computation time of the whole static problem (n t = 1500) using the PGD piecewise formulation with different subinterval sizes. In each computation the size of all subintervals is equal to n ts time steps problem presented in Fig. 3. The motivation of this formulation is to ensure the convergence of the solution and minimize computation time by solving the problem sequentially on a number of smaller subintervals. We propose to divide the Table 4 Numerical setup of the piecewise PGD formulation for solving the static problem with a large number of time steps (n t = 1500) illustrated in Fig. 3 Fixed parameters n s n t t [s] 60 1500 10 −2 overall time interval of n t = 1500 time steps into a number of n p subintervals with n ts = n t /n p time steps, so that a spacetime solution h U ∈ R n s ×n ts is computed in each subinterval h = 1, · · · , n p . The first and last columns of the matrix h U correspond to the solutions at time steps j = n ts (h − 1) + 1 and j = n ts h, respectively. The plot in Fig. 9 shows the evolution of the required computation time to solve the whole static problem of n t = 1500 time steps when the size n ts of the subintervals is varied. The numerical setup of the computations is presented in Table 4. The termination criterion of the enrichment loop is only based on reaching a certain residual matrix norm, which is proportional to its size, so that all the considered computations exhibit approximately the same accuracy. The initialization of all columns of the space-time solution 1 U (0) of the first subinterval is defined as the zero vector associated with the initial undeformed configuration. All the columns of the space-time solution h U (0) of the subsequent subintervals are initialized with the solution h−1 u n ts at the last time step of the corresponding previous subinterval. We conclude from the plot in Fig. 9 that the piecewise PGD formulation is more efficient for a larger number n p of subintervals, which have a smaller number of time steps n ts . This conclusion is in accordance with the behavior observed in Fig. 8. Although the overall time interval has to be divided into subintervals, this method enables solving problems with an arbitrarily large time interval using the PGD solver in algorithm 2.

Geometrically nonlinear dynamic problem
In this section, we present the numerical demonstration of the PGD approach using the cantilever beam geometry, as shown in Fig. 2, subjected to dynamic excitations. In doing so, we present the applicability of the new strategy considering two different excitation types, the harmonic excitation problem in Sect. 4.2.1 and the white noise excitation problem in Sect. 4.2.2.

Harmonic excitation problem
In order to solve the dynamic problem associated with the cantilever beam model in Fig. 2, the equations of motion (9) have to be satisfied at all time steps. Using the PGD solver in algorithm 2, the norm of the matrix composed of the residual vectors of all time steps (47) is minimized by adding new enrichments to the space-time solution U. The values of the external moment, applied at the right-most node, are defined with the following sinusoidal function M(t) = M 0 sin ωt. The equations of motion are discretized in time using the Newmark integration scheme, defined by the equations (4) and (5), choosing β = 0.25 and γ = 0.5 for the Newmark parameters. The integration of the time discretization using the PGD formulation is discussed in Sect. 3.2. The damping matrix C is defined according to the Rayleigh damping approximation (2). The weighting factors α d and β d are adjusted, so that the first two eigenfrequencies of the linearized problem around the undeformed configuration have a damping ratio of ζ = 0.005.
The numerical simulations discussed in this section demonstrate the ability of the proposed PGD solver to solve geometrically nonlinear dynamic problems. They are performed in accordance with the numerical setup presented in Table 5. For the considered dynamic problem, the motion of the cantilever beam, consisting of n s = 60 degrees of free- and the incremental Newton-Raphson solver. A PGD-m max solution includes a number of m max enrichments. The illustrated displacement is one of the 60 degrees of freedom of the considered finite element model. Its evolution over time is one part of the total solution and is extracted from the corresponding row in the PGD solution matrix U. The red curve PGD-1 is extracted from the rank one solution matrix U (1) , which is the outer product of the first space mode u (1) s and its corresponding time mode u (1) t , shown in the top of Fig. 11. The orange curve PGD-2 and the green curve PGD-3 in Fig. 10 are extracted from the PGD solution matrices U (2) of rank two and U (3) of rank three, which are defined as a sum of the two first enrichments and the three first enrichments, respectively. Convergence of the PGD solution towards the Newton-Raphson reference solution in black is observed when a larger number of enrichments is considered. It should be mentioned that a faster convergence of the solution is achieved in the periodic steady-state at the end of the time interval than during the transient behavior at the beginning. The space and time modes corresponding to the first five enrichments are illustrated in Fig. 11. Note that the space modes are normalized, so that the extent of the contribution of a certain enrichment to the total solution is represented by the time modes. With some exceptions, the contribution of a new computed enrichment has, in general, a smaller magnitude than the previous enrichments.
The convergence behavior of the solution of the dynamic problem computed with the proposed nonlinear PGD solver can also be shown quantitatively by observing the evolution of the relative error measure in Fig. 12a or the Frobenius norm of the residual matrix in Fig. 12b with respect to the number of enrichments. The convergence in residual norm is nonmonotonic since the proposed nonlinear PGD solver is based on the Galerkin formulation. However, the convergence of the relative error measure shown in blue in figure 12a exhibits monotonic behavior. The relative error measure is computed in accordance with the formulation equation (61), where the PGD solution U (m) after the m-th enrichment is inserted as U app and a highly accurate Newton-Raphson solution is used as U re f . The evolution of the computation time, illustrated by the green line in Fig. 12a, is linear with respect to the number of enrichments and exhibits similar behavior as in the static problem since the largest part of the computational effort is required for the evaluation of the stiffness matrices and internal force vectors. The computational time of these evaluations is approximately the same in all enrichment steps because they are, in general, independent from the considered current solution U (m−1) . For the numerical setup in Table 5, Fig. 12a shows that the computation of the PGD solution with two enrichments (m max = 2) is more than twice as fast as the computation of the reference Newton-Raphson solution, which is represented by the horizontal dashed green line, while it exhibits a lower level of accuracy. The PGD solution based on three enrichments is slightly faster than the Newton-Raphson solution, and starting from the fourth enrichment, it   becomes slower. The achieved speedups in solving nonlinear dynamic problems with the proposed nonlinear PGD solver can be increased by making the evaluations of the stiffness matrices and the internal force vectors more efficient.

White noise excitation problem
The cantilever beam is excited by a concentrated load, as depicted in Fig. 13a. The time history is realized as a sta-tionary Gaussian random white noise [71], as depicted in Fig. 13b. The required parameters for the calculation are presented in Table 6. The solution of the transient excitation problem using the proposed PGD approach is depicted in Fig. 14 and compared to the classical step-by-step Newton-Raphson procedure. After one enrichment the PGD solution (red line) converges to a time history that is very similar to the geometrically linear response (gray dashed line). At this stage of the algorithm, the solution is spanned by the dyadic product of the first space and time enrichment modes, presented in the first row of Fig. 15. One observes a high correlation of the space mode with the first mode of vibration while the time mode reveals the time history of the linear response function. The linear solution, presented by the gray dashed line in Fig. 14, is, in this example, essentially governed by the response of the lowest vibration mode, which is, therefore, responsible for the highest amount of the vibration energy of the geometrically linear system. Since the PGD solution almost coincides with the geometrically linear solution after the first enrichment one observes a high similarity of the response with the modal truncated representation considering only the lowest linear mode of vibration. We further present the convergence behavior of the PGD algorithm in Fig. 14, plotting the solution after two, three, and 20 enrichments, in orange, green and blue color, respectively. One observes very good agreement with the Newton-Raphson solution after 20 enrichments.   Table 5. They constitute the enrichment couples computed in five different iterations m = [1,5,10,15,20] (a) (b) Additionally, we present the space and time modes after the fifth, tenth, 15th, and 20th enrichment in the second, third, fourth, and fifth row of Fig. 15, respectively. Looking at the time modes in this figure, one clearly observes that the participation of the enrichments to the full solution decreases significantly with increasing enrichment, although we again emphasize here that the convergence behavior of the PGD algorithm is not strictly monotonic decreasing with respect to the Frobenius norm of the residual matrix, as presented in Fig. 12b.

Conclusion
In this paper, a proper generalized decomposition based solver for nonlinear dynamic and static structural problems is proposed. The proper generalized decomposition based solver builds the solution of an initial boundary value problem in the whole time domain as sum of rank one space-time enrichments in an iterative process. Within every iteration a space mode and its corresponding time mode are computed by solving two coupled nonlinear problems, one in space and the other in time, in an alternating manner using the fixed point algorithm until convergence is achieved.
The performed numerical examples show that the proposed PGD solver achieves a limited speedup in comparison with the incremental Newton-Raphson solver when a small number of enrichments is considered. One key feature of a successful numerical speedup is the consideration of the sparsity of the tangential stiffness matrix during the vectorization operation of the stiffness matrix in Equation (57), which is required to obtain a space-time stiffness formulation. The computational bottleneck causing the limitation of the proposed approach to date is mainly the evaluation of the stiffness matrices and the internal force vectors during the linearization step, which is also a key feature in classical model order reduction approaches that follow the step-by-step scheme. However, the proposed method follows an algorithmic approach that differs significantly when compared to classical step-by-step methods and the possibilities regarding numerical speed-up and applicability to various types of mechanical problems are yet to discover.
Thus, the next step should be the integration of strategies into the proposed algorithm, which decreases the computational effort by reducing the number of required costly stiffness matrix evaluations while preserving a good convergence behavior. The additional application of hyperreduction strategies for a more efficient evaluation of internal force vectors and stiffness matrices is also an interesting topic for future research. Furthermore, future research should deal with the extension of the space-time solutions to three dimensional beam structures and buckling of beams as well as problems including external loads that depend on the displacements. It is of high interest if the space-time algorithm can also be applied if such challenging types of nonlinearities are present or if additional adaptions may be required to ensure stable convergence. This ambitious enterprise should not only include extensions or adaptions of the proposed space-time algorithm but also the implementation of a numerically efficient framework for practical applications.
Finally, we note that the new method should be compared to existing model reduction strategies for geometrically nonlinear problems, as presented by Sonneville et al. [7]. This objective demands rigorous numerical investigations that are beyond the scope of this paper, however, they are on the agenda for future research projects.
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://creativecomm ons.org/licenses/by/4.0/.

Appendix A: Two-dimensional geometrically exact beam element
In the following, we present the theory of the two-dimensional geometrically exact beam proposed in [72] and its finite element discretization, discussed in [73], for initially straight beam axis. This theory is based on the assumption that plane cross sections remain plane and takes into account the transverse shear deformation, which means that the cross sections do not necessarily remain normal to the axis of the beam. It results in nonlinear strain displacement relations and is valid for finite deflections, rotations, and strains. In many applications, it is sufficient to consider small strains, even when deflections and rotations are large, so that the Saint Venant-Kirchhoff constitutive equations can be used. This hyperelastic constitutive law defines linear stress-strain relations and can be considered as an extension of Hooke's law of linear theory of elasticity to the geometrically nonlinear regime. The nonlinear kinematics of the beam are shown in Fig. 16. In the undeformed reference configu- Fig. 16 Nonlinear kinematics of the geometrically exact beam ration B 0 the local axis X coincides with the initially straight beam axis. The beam can exhibit an arbitrary initial inclination with respect to the global axis X g .
In the current configuration B, the position of each material point is expressed with respect to the local coordinate system by the following mapping function: where the three independent fieldsū(X ),w(X ) andψ(X ) are the displacement in axial direction, the deflection and the rotation of the cross section, respectively. The fields are arranged in a vectord = [ū,w,ψ] T , which can be transformed to the global coordinate system using a rotation matrix as: The nonlinear strain measures derived in [72] by using the principal of virtual work are the axial strain, , the shear strain, γ and the curvature, κ. They can be written in matrix notation as: where (·) denotes the partial derivative with respect to the coordinate X and the following matrices and vectors are introduced: The linear hyperelastic constitutive law, Saint Venant-Kirchhoff model, yields the following cross-sectional force vector: The entries of the force vector S are the axial force N , the shear force Q and the bending moment M. The first diagonal entry of the matrix D is the product of the Young's modulus E and the cross section area A, which describes the axial stiffness of the beam. The second diagonal entry is the product of the shear modulus G and the shear-corrected cross section area A s , which defines the shear stiffness of the beam. A shear-corrected cross section area is introduced in order to correct the violation of the boundary condition at Y = ±h/2, where the shear stresses should vanish. The violation of this boundary condition is caused by the assumption that plane cross-sections remain plane. The third diagonal entry, the bending stiffness, is computed as the product of the Young's modulus E and the second moment of area I .
The weak form of equilibrium, which is equivalent to the principle of virtual work, is written for the considered geometrically nonlinear beam model as: where the virtual displacements: and the virtual strains: are introduced. The vector g is the vector of external forces and moments. The spatial discretization is performed by the finite element method using quadratic Lagrangian shape functions [73]. The three quadratic shape functions associated with the three nodes of the beam element are expressed by the local coordinate ξ ∈ [−1, 1] in the parameter space as: The mapping between the physical space and the parameter space is illustrated in figure 17. The geometry X and the displacement fieldsd are both discretized using the shape functions (A13) as: and where X e i is the initial position of node i along the X axis of the local coordinate system of an element e and d e i ∈ R 3 is the corresponding vector of local nodal quantities of the three displacement fields, which are the primary variables in this formulation. The discretized displacement fields can also be written as: Fig. 17 The mapping between the parameter space and the physical space of a beam element e with three nodes