A Space-Time POD Basis Interpolation on Grassmann Manifolds for Parametric Simulations of Rigid-Viscoplastic FEM

. Parametric simulations of thermomechanical metal forming processes still remain computational costly and di ﬃ cult due to inherent strong non-linearities. To this end, Reduced Order Models (ROMs) are introduced to decrease the computational time in large scale models and provide near-optimal solutions in acceptable times. ROMs based on the Proper Orthogonal Decomposition (POD) are usually capable of accurately reproducing the dynamics of high-ﬁdelity FEM simulations and o ﬀ er the potential for near real-time analysis. However, ROMs are not robust with respect to parameter changes and must often be rebuilt for each parameter variation. This work aims to interpolate ROM POD basis associated with a limited number of training points on Grassmann manifolds, so as to obtain a robust ROM corresponding to a target parameter. A novel Space-Time (ST) POD basis interpolation, where the reduced spatial and time basis are separately interpolated on Grassmann manifolds, is proposed. Good correlations of the ROM ST models with respect to their associated high-ﬁdelity FEM counterpart simulations are found. Hence, application of the ROM adaptation method for near real-time metal forming simulations using o ﬀ -line computed ROM POD databases can be possible.


Introduction
Computational metal forming has been widely used in a variety of applications in academia laboratories and manufacturing industry over the last decades. Simulations of such problems are characterized by multiple sources of strong non-linearities, thereby are often prohibitively expensive for an exhaustive design analysis of large scale models. Moreover, several core issues are still of significant importance and development, among others, parallel computing, adaptive meshing and remeshing, modeling of micro-structure evolution, multiscale modeling, optimization of processes or parameter identification, stochastic approaches and machine learning (see a thorough review in [1,2]). One of the key feature challenges mentioned in [1], is the model order reduction methods to decrease computational time. Indeed, simulation of complex configurations can be intractable since the computational cost can highly increase for a parametric analysis. ROMs via the POD have been chosen to reduce the problem dimensionality while maintaining solution accuracy. For solving a parametric problem, the method starts by a training stage during which the problem is solved for a limited number of parameters. The full-order field 'snapshots' are then compressed using the POD to generate a ROM that is expected to reproduce the most characteristic dynamics of its high-fidelity counterpart solution. For a new parameter value, an interpolation method is proposed using the underlying spanned subspaces of the POD basis [3]. Since these matrices are obtained via the Singular Value Decomposition (SVD), they describe the optimal projection operators, which reduce the problem from the full space to a low dimensional subspace. POD basis interpolation is performed using local maps on Grassmann manifolds by evaluating the geodesic paths between the subspaces on this manifold, all this being done in the framework of Riemannian geometry, which is a specific matter of differential geometry. In this context, a novel non-intrusive Space-Time POD interpolation is proposed here allowing near-real time predictions since no new ROM FEM models have to be solved, as required by the standard POD approach.
The coupled thermomechanical analysis is based on the Rigid Visco-Plastic (RVP) formulation using an implicit solution strategy. In RVP, the elasticity effects are completely neglected due to the fact that elastic components of strain remain small as compared with irreversible strains. Therefore, the RVP formulation turns out to be very similar of fluid flow problems and it is also called as flow formulation. Although it is not possible to calculate the residual stresses and the spring-back effect, the flow formulation presents outstanding advantages. Unlike the elasto-plastic FEM, the RVP formulation, even though more approximate, is simpler to be implemented in computer codes, more stable and can use relatively larger time increments, thus improving the computational efficiency. A thorough overview of the foundation of the theory can be found in [4]. In this work, RVP simulations are performed using a mechanical FEM solver [5] and FEM solver for the thermal analysis in Matlab. A staggered approach for the thermomechanical coupling is used.

Proper Orthogonal Decomposition and Grassmann Manifolds
Let us consider here some snapshot matrix obtained from a FEM solution of a RVP problem with N s the number of spatial points and N t the number of time steps. Such matrix encodes in fact N t vectors u k := u(·, t k ) ∈ H spatial , and we write As a classical result [6], POD of mode p can be obtained from any singular value decomposition so to define a projection matrix. Taking φ 1 , . . . , φ k (k = min(n, N t )) to be the spatial vectors associated to the singular values σ 1 ≥ σ 2 ≥ . . . ≥ σ k of S, POD of mode p is obtained using . so that the projection matrix Π p is defined using the matrix issued from vectors φ i . A main observation is that a projection matrix can be defined using any orthonormal basis v 1 , . . . , v p of the vector space V p generated by φ 1 , . . . , φ p .
Finally, POD of mode p is thus characterized by the p dimension subspace V p of R n and not by the spatial vectors φ 1 , . . . , φ p . Any interpolation made on POD translate into interpolation made on p dimensional subspaces.
Mathematically, the set of all p dimensional subspaces of R n is known as a Grassmann manifold G(p, n) which is also a Riemannian manifold. Thus, a POD of mode p made on a snapshot matrix S leads to some point m in G(p, n), and any interpolation has to be done on this curved manifold. This means that computations on G(p, n) can only be done using local charts, and the Riemannian structure allows us to use normal coordinates obtained from an exponential map defined using geodesics.
To address the question of interpolation for a target value λ, we have to make interpolation on the Grassmannian manifold G(p, n). A strategy introduced in [3] make use of normal coordinates, and thus the exponential and the logarithm map. This leads to the following, for the spatial part (while the same is done from the time part): • Step 1: Choose a reference point m 0 ∈ {m 1 , . . . , m N } ⊂ G(p, n). For each vector space m i , take Y i to be the matrix associated to one of its orthonormal basis. • Step 2 (the logarithm step): First compute thin SVD Step 4 (the exponential step): Take a thin SVD of Z given by Z = U Σ V T and define A Space-Time (ST) interpolation is finally some ROM snapshot matrix S obtained from a spatial part Φ, a time part Ψ (both obtained from previous algorithm), and a mixed part S * defined to be some square matrix of size p × p, so that

Rigid-Viscoplastic FEM Formulation
Classical RVP problems consider the deformation of an isotropic body occupying a domain Ω ⊂ R 3 and its associated boundary ϑΩ, which represents the current configuration according to the Updated Lagrangian formulation. A weak form is applied to obtain a solution that satisfies the governing equations and the boundary conditions. The essence of the variational formulation is to calculate the total potential of a scalar quantity Π (functional) of the system, and to invoke the stationarity of Π, i.e. δΠ = 0 with respect to arbitrary changes of the state variables. The functional Π is defined by an integral form in accordance with the virtual work-rate principle whereσ is the effective stress,ε is the effective strain rate and F i denotes prescribed surface tractions on the boundary surface Ω F . The first term in eq. (2) represents the internal deformation work-rate, whereas the second term represents the work-rate done by the external forces. The real velocity field gives to the functional Π a stationary value, which means that the first order variation of Π vanishes. In order to meet the incompressibility constraint conditionε v =ε kk = 0 on an admissible velocity field, a penalized form of incompressibility is used where K is a large positive constant which penalizes the dilatational strain-rate. For the thermal problem, a thermodynamically sound derivation is adopted based on the conservation of energy − ρc ϑT ϑt where ρc is the volume specific heat of the material, ξσε represents the work heat rate per unit volume due to plastic deformation, k is the thermal conductivity, T is the temperature and ξ is a coefficient that presents the fraction of the deformation energy dissipated into heat. Using the weak form and after the application of the divergence theorem, the energy balance eq. (4) can be written where q n = kϑT/ϑn is the heat flux across the boundary ϑΩ and n denotes the unit normal vector to the boundary surface ϑΩ. The mechanical and thermal analysis can be loosely coupled in such a way that plastic work and friction are considered as heat source in the thermal analysis while the updated temperature field is used to determine the flow stress behavior during the deformation analysis.

Numerical Investigations
To evaluate the performance of ST POD interpolation, a rectangular cross section bar is compressed between two parallel flat dies under the conditions of constant shear friction factor m at the die/workpiece interface. The initial workpiece has dimensions h = 20 mm in height and w = 20 mm in width. Plane strain conditions are considered. Due to the symmetry of the problem, only one quarter of the cross section is analyzed (Figure 1). The lower die is stationary while the upper die velocity is set to v = 1 mm/s. The workpiece is compressed till a 35% reduction in height is achieved. The final state is accomplished in 7 time steps with a constant time increment ∆t = 0.5 s. A 4-node quadrilateral element with bilinear shape functions and four point integration is used. However, a reduced integration scheme which imposes the volume constancy over the linear element (one point integration) is used for the dilatational term ter, the following training points are selected, denoted by the variable λ ∈ Λ tr = {0.1, 0.5, 0.9}. A simple strategy for sampling the parametric space is used by equi-distribute the points over the associated range. The target point is set toλ = 0.3. For each parametric simulation, a sequence of snapshots uniformly distributed over time using an increment of ∆t = 0.5 s is extracted for all nodes of the workpiece. The space-time snapshot matrices S(λ i ) ∈ R n×N t of size 242 × 7 are associated to nodal velocity fields. Using the POD/SVD, the ROM spatial basis {m i } N i=1 and time basis {m i } N i=1 are constructed for the training points λ i . Then, the set of the low-dimensional POD basis is interpolated on Grassmann manifolds using the algorithm presented in section 3. ST interpolation is compared against the high-fidelity FEM solution. Considering the interpolated and the HF-FEM snapshot matricesS and S FEM , respectively, the following error measure is defined where u i denotes the i-th column of the snapshot matrix S. The eigenvalue spectrum of snapshot matrices S i corresponding to training points λ i ∈ Λ tr is shown in a log scale in Figure 2. We can observe that the distance between the first and the last eigenvalue is of the order of 10 5 -10 6 . The relative L 2 -error norm e L 2 (S) between the interpolated and the HF-FEM solution for various POD modes is shown in Figure 3. The relative error for all POD modes lies within a range of 0.0175 up to 0.038. For this particular case, the interpolation error is minimum with p = 2 and p = 3 POD modes. In general, it can be observed that the interpolated ST POD solution delivers good accuracy and is reliable enough to predict the velocity field for the investigated target point. Due to the fact that the online stage of ST POD interpolation consists solely of algorithmic matrix operations, it outperforms standard POD methods which require a FEM solution for the new ROM. Thus, in ST POD, the computational time can be several orders of magnitude lower than standard POD approaches, and even more against the high-fidelity FEM counterpart solution making it specially tailored for near real-time computations.