A Non-Intrusive Space-Time Interpolation from Compact Stiefel Manifolds of Parametrized Rigid-Viscoplastic FEM Problems
Orestis Friderikos, Marc Olive, Emmanuel Baranger, Dimitrios Sagris, Constantine David
AA NON-INTRUSIVE SPACE-TIME INTERPOLATION FROMCOMPACT STIEFEL MANIFOLDS OF PARAMETRIZEDRIGID-VISCOPLASTIC FEM PROBLEMS
O. FRIDERIKOS, M. OLIVE, E. BARANGER, D. SAGRIS, AND C. DAVID
Contents
1. Introduction 22. Space–Time POD, Grassmann and compact Stiefel manifolds 63. Space-Time Interpolation on compact Stiefel manifolds 144. Rigid-Viscoplastic FEM Formulation 205. Numerical Investigations 246. Conclusions 29References 42
Abstract.
This work aims to interpolate parametrized Reduced Order Model(ROM) basis constructed via the Proper Orthogonal Decomposition (POD) toderive a robust ROM of the system’s dynamics for an unseen target parametervalue. A novel non-intrusive Space-Time (ST) POD basis interpolation scheme isproposed, for which we define ROM spatial and temporal basis curves on compactStiefel manifolds . An interpolation is finally defined on a mixed part encoded ina square matrix directly deduced using the space part, the singular values andthe temporal part, to obtain an interpolated snapshot matrix, keeping track ofaccurate space and temporal eigenvectors. Moreover, in order to establish a well-defined curve on the compact Stiefel manifold, we introduce a new procedure,the so-called oriented SVD. Such an oriented SVD produces unique right and lefteigenvectors for generic matrices, for which all singular values are distinct. It isimportant to notice that the ST POD basis interpolation does not require theconstruction and the subsequent solution of a reduced-order FEM model as clas-sically is done. Hence it is avoiding the bottleneck of standard POD interpolationwhich is associated with the evaluation of the nonlinear terms of the Galerkin pro-jection on the governing equations. As a proof of concept, the proposed methodis demonstrated with the adaptation of rigid-thermoviscoplastic finite elementROMs applied to a typical nonlinear open forging metal forming process. Strongcorrelations of the ST POD models with respect to their associated high-fidelityFEM counterpart simulations are reported, highlighting its potential use for nearreal-time parametric simulations using off-line computed ROM POD databases. a r X i v : . [ m a t h . DG ] F e b O. FRIDERIKOS, M. OLIVE, E. BARANGER, D. SAGRIS, AND C. DAVID
Notations
Mat n,p ( R ) Set of n × p matrices in R I p Identity matrix in Mat p,p ( R )[ y , . . . , y p ] Matrix in Mat n,p ( R ) Matrix with column vectors y i ∈ R n O( p ) Orthogonal group on R p (cid:8) Q ∈ Mat p,p ( R ) , Q T Q = I p (cid:9) G ( p, n ) Grassmann manifold Set of p linear subspaces in R n π − ( m ) Fiber at m ∈ G ( p, n ) If y , . . . , y p is an orthonormal basis of m π − ( m ) = { Y Q, Q ∈ O( p ) , Y = [ y , . . . , y p ] }T m := T m G ( p, n ) Tangent space of G ( p, n ) at m For Y ∈ π − ( m ), one model of T m is (cid:8) Z ∈ Mat n,p ( R ) , Z T Y = 0 (cid:9) S t ( p, n ) Stiefel manifold Set of ordered p -tuples independentvectors in R n S t c ( p, n ) Compact Stiefel manifold Set of ordered p -tuples of orthonormalvectors in R n S t c ( p, n ) = (cid:8) Y ∈ Mat n,p ( R ) , Y T Y = I p (cid:9) Hor Y Horizontal space at Y Hor Y := (cid:8) Z ∈ Mat n,p ( R ) , Z T Y = (cid:9) v ∈ T m Velocity vector on the tan-gent plane T m Represented by a horizontal lift Z ∈ Hor Y , with Y ∈ π − ( m ) S ( i ) Snapshot matrix S ( i ) ∈ Mat n,m ( R ) corresponding to pa-rameter value λ i Introduction
Computational metal forming has been widely used in academic laboratories andthe manufacturing industry over the last decades, becoming nowadays a mature, wellestablished technology. Nevertheless, new challenging fields are emerging, amongothers, uncertainty quantification, optimization of processes and parameter identifi-cation in design analysis [1, 2]. One of the key challenging topics mentioned in [1]is the introduction of Model Order Reduction (MOR) methods to combat the highcomputational cost, which is also of paramount interest in the above-mentionedfields. Moreover, due to the multiple sources of strong non-linearities inherent inmanufacturing problems, design optimization and multi-parametric studies of largescale models turns out to be prohibitively expensive. Indeed, simulation of complexconfigurations can be intractable since the computational times can highly increase.To this end, meta-model techniques are often used to tackle the computationalburden. These rely on a manifold learning stage during which we need to capturethe original space where the solution of the model problem lies. This data collectionconsists of solving the full-scale model for an ensemble of training data over theparametric range and is commonly referred to as the offline stage. Even though
T INTERPOLATION OF PARAMETRIZED RIGID-VISCOPLASTIC FEM PROBLEMS 3 meta-models can speed up the simulation time, nevertheless their construction withstandard computations based on full-order models is expensive.Closely related to the concept of metamodeling, Reduced Order Models (ROMs)have been chosen to reduce the problem’s dimensionality while at the same timemaintaining solution accuracy. ROMs can decrease the computational complexityof large-scale systems, solving parametrized problems and offering the potential fornear real-time analysis. The methods for building ROMs can be classified into twogeneral families: a priori and a posteriori ones. The well known a priori MORincludes methods such as the Proper Generalized Decomposition (PGD) [3], andthe a priori reduction method (APR) [3, 4]. The main characteristic of all thesemethods is that they do not require any precomputed ROMs. In the second class ofmethods, the reduced basis is built, a posteriori , from the state variable snapshotsin the parametric space. One popular method is the POD [5, 6, 7], also knownas Kharhunen-Lo`eve Decomposition (KLD) [8, 9], Singular Value Decomposition(SVD) [10] or Principal Component Analysis (PCA) [11, 12, 13, 14].For nonlinear systems, even though a Galerkin projection reduces the number ofunknowns, however, the computational burden for obtaining the solution could stillbe high due to the computational costs involved in the evaluation of nonlinear terms.Hence, the nonlinear Galerkin projection principally leads to a ROM, but its evalua-tion could be more expensive than the corresponding one of the original problem. Tothis effect, to make the resulting ROMs computationally efficient, a sparse samplingmethod is used, also called hyper reduction, to mention among others, the missingpoint estimation (MPE) [15], the empirical interpolation method (EIM) [16], thediscrete empirical interpolation method (DEIM) [17], the Gappy POD method [18],and the Gauss-Newton with approximated tensors (GNAT) method [19]. Thus, allthese methods imply the solution of a new ROM FEM problem.In the case of a parametric analysis using POD basis interpolation on Grassmannmanifolds [20, 21], the method starts with a training stage during which the problemis solved for several training points. Then, using the FEM solutions, the full-orderfield ‘snapshots’ are compressed using the POD to generate a ROM that is expectedto reproduce the most characteristic dynamics of its high-fidelity counterpart solu-tion. However, the relevant information is contained in the vector spaces generatedby the (left or right) singular vectors of the snapshot matrices. Now, for a newparameter value, interpolation methods have to be defined from such relevant sub-spaces spanned by the POD basis vectors [20]. Other approaches obviously couldbe considered, such as interpolations computed on the space of matrices of a fixedrank, whereby the mechanical origin of the problem imposes to consider the vectorsubspaces, and not the matrices themselves [21]. Nevertheless, such methods as the
O. FRIDERIKOS, M. OLIVE, E. BARANGER, D. SAGRIS, AND C. DAVID one of interpolation between two positive semidefinite matrices of fixed rank [22],may not capture the important elements obtained from the mechanical equations.To interpolate between different vector spaces of the same dimension (encodedinto the mode p of the POD), a Grassmann manifold [23] must be used, which isthe set of p -dimensional subspaces of R n . Such a manifold is in fact a Riemann-ian manifold [24], so we can construct geodesics between two points, and use suchgeodesics to define a logarithm map to linearize , and conversely using the exponen-tial map to return back to the Grassmann manifold. While an interpolation cannotbe done directly on Grassmann manifolds, linearization allows computing such aninterpolation, at least locally once a reference point has been selected [20, 21]. Toany new parameter value, thus we get a new subspace obtained from interpolationbetween all subspaces related to the spatial eigenvectors of the snapshot matrices.Another approach using inverse distance weighting was initiated in [21, 25], but italso relies on several choices (as one of the weights). Furthermore, an extension ofNeville-Aitken’s algorithm to Grassmann manifolds which computes the Lagrangeinterpolation polynomial in a recursive way from the interpolation of two points wasrecently presented [26].In the standard POD interpolation mentioned above [20], the spatial ROM basiscorresponding to the target point is used to generate a ROM FEM, which is ex-pected to have a lower computational cost compared to the high-fidelity problem.The key idea in the Space-Time (ST) POD basis interpolation proposed by [27, 28],is that the reduced spatial and temporal basis are considered separately, both defin-ing points on two different Grassmann manifolds. However, such points are stronglyrelated: a spatial vector directly corresponds to a temporal vector, and vice versa .From this, firstly we need to consider the p -tuples of spatial (and temporal) vectors,instead of the p -dimensional subspace, which defines points on an associated compactStiefel manifold , strongly connected to Grassmann manifolds. Contrary to what issuggested in [28], we propose a different interpolation scheme, as we do not performinterpolation of the singular values, followed by spatial and temporal calibration.Instead, we exploit the dependence between the spatial and temporal parts. Indeed,using an interpolation algorithm defined on a Grassmann manifold, we derive curveson a compact Stiefel manifold , which are no longer interpolating, but which never-theless allow us to obtain new singular vectors for the spatial part, and separatelyfor the temporal part. Such space and temporal singular vectors finally are taken todefine a mixed part on which a classical interpolation can be computed. In the end,we get in this way a ROM matrix corresponding to a new parameter value. Notethat in order to obtain a well-defined curve on compact Stiefel manifolds, we haveto introduce a new procedure, the so-called oriented SVD . Such an oriented SVD T INTERPOLATION OF PARAMETRIZED RIGID-VISCOPLASTIC FEM PROBLEMS 5 produces unique right and left eigenvectors for snapshot matrices, supposed to be generic matrices , for which all non–zero singular values are distinct.The off-line stage in the ST approach consists of solving FEM problems whichare corresponding to the training points of the given parameter. The on-line stageconcerns the use of a curve defined on a compact Stiefel manifold to determinethe spatial and temporal ROM basis for the target point, in order to construct therelated ROM snapshot matrix. In fact, the ST interpolation offers the advantage ofreconstructing a snapshot matrix without relaunching ROM FEM computations. Tothis end, it results in near-real-time solutions due to direct matrix multiplications inthe on-line stage.We could also mention some other ST approaches [29, 30, 31, 32], where neitherGrassmann nor compact Stiefel manifolds are considered. For instance, an approx-imation of the spatial and temporal basis functions by linear interpolation of theirmodes is proposed in [29] to study the flow past a cylinder at low Reynolds numbers.A non-intrusive ROM approach for nonlinear parametrized time-dependent PDEsbased on a two-level POD method by using Radial Basis Functions interpolation ispresented in [30, 33].The method proposed in this work is applied to a coupled thermomechanical rigidvisco-plastic (RVP) FEM analysis based on an incremental implicit approach [34,35, 36, 37]. Note that the RVP formulation specifically is tailored for metal formingsimulations, where the plastic flow is unconstrained and usually of finite magni-tude, involving large strain-rates and high temperatures. In the present study, allsimulations are performed by using an in-house Matlab code which consists of twoindependent FEM solvers. A mechanical solver for the viscoplastic deformation anal-ysis [38] and a thermal solver for the heat transfer analysis. A staggered procedureis used to solve the system of coupled equations.The paper is organized as follows: in section 2, the Proper Orthogonal Decom-position is presented, followed by an introduction to some basic notions about thegeometry of the Grassmann and Stiefel manifolds to make the article reasonable self-contained. POD basis interpolation on Grassmannian manifolds is introduced con-sidering the underlying formulation of the logarithm and the exponential map. Thecore of this paper is illustrated in section 3, where the computational framework forthe ROM adaptation based on a novel non-intrusive Space-Time POD basis interpo-lation on compact Stiefel manifolds is developed. The following section 4 covers therigid visco-plastic formulation, the general framework of the thermal field equations,and the thermomechanical coupling. In section 5, the interpolation performance ap-plied to a metal forming process is shown, as well as further computational aspectsare discussed. Finally, section 6 highlights the main results and some importantoutcomes.
O. FRIDERIKOS, M. OLIVE, E. BARANGER, D. SAGRIS, AND C. DAVID Space–Time POD, Grassmann and compact Stiefel manifolds
Let us recall here the important link between Proper Orthogonal Decompositionand Grassmann manifold [20, 39, 40, 21, 25].Assume S ∈ Mat n,m ( R ) to be any real matrix of size n × m (with n ≥ m ), takenhere to be a snapshot matrix with n = 3 N S obtained from the spatial discretization N s , and m = N t obtained from the time one. Any spatial POD of mode p leads toa p -dimensional vector space V p ⊂ R m such that the Frobenius norm (cid:107) S − Π p S (cid:107) is minimal, where matrix Π p corresponds to the orthogonal projection on V p (see[21] for more details). Such a matrix Π p is directly obtained from a Singular ValueDecomposition (SVD) of S . Indeed, let us write a SVD S = ΦΣΨ T with Φ = [ φ , . . . , φ r ] and Ψ = [ ψ , . . . , ψ r ], where the columns φ k ∈ R n and ψ k ∈ R m form a set of orthonormal vectors, and Σ ∈ Mat r,r ( R ) is a diagonal matrix, where r denotes the rank of S . Then, we can define Φ p := [ φ , . . . , φ p ] ∈ Mat n,p ( R ) and weobtain Π p = Φ p Φ Tp .In this classical approach, the relevant object is not the reduced matrix S p := Π p S ,supposed to be of maximal rank, but the p -dimensional vector space V p spanned byvectors φ , . . . , φ p , and thus the image of the matrix Φ p . From this, interpolationhas to be considered on the set of all p -dimensional vector spaces, that is on theso–called Grassmann manifold G ( p, n ): G ( p, n ) := {V p ⊂ R n , dim( V p ) = p } . Note here that the point m := V p ∈ G ( p, n ) defines a vector space spanned by theset φ , . . . , φ p represented by matrix Φ p , however this matrix representation is notunique (see Example 2.2).Take now a set { λ , . . . , λ N } of parameter values leading to snapshot matrices S (1) , . . . , S ( N ) with SVD S ( k ) = Φ ( k ) Σ ( k ) Ψ ( k ) , Φ ( k ) = [ φ ( k )1 , . . . , φ ( k ) r ] , Ψ ( k ) = [ ψ ( k )1 , . . . , ψ ( k ) r ] , where φ ( k ) i are orthonormal vectors in R n and ψ ( k ) j are orthonormal vectors in R m .The classical approach [20, 21] then considers the spatial POD of the snapshotmatrices S (1) p , . . . , S ( N ) p of mode p , so that we obtain points m i ( i = 1 , . . . , N ) on G ( p, n ), respectively represented by the matrices Φ ( k ) p := [ φ ( k )1 , . . . , φ ( k ) p ] ∈ Mat n,p ( R ) , (cid:16) Φ ( k ) p (cid:17) T Φ ( k ) p = I p . T INTERPOLATION OF PARAMETRIZED RIGID-VISCOPLASTIC FEM PROBLEMS 7
To any new parameter value (cid:101) λ , it is possible to make an interpolation consideringthe spatial part based on the points m i ∈ G ( p, n ), using a local chart given bynormal coordinates [20, 21, 25], in order to obtain a point (cid:101) m ∈ G ( p, n ) representedby a matrix (cid:101) Φ . From such a point (cid:101) m ∈ G ( p, n ), we deduce a p -dimensional vectorspace on which some POD-Galerkin approach [21] can lead to a new ROM model.On the contrary, we propose another approach as we consider a Space–Time in-terpolation , using both the spatial vector spaces represented by matrices Φ ( k ) p andthe temporal vector spaces represented by matrices Ψ ( k ) p := [ ψ ( k )1 , . . . , ψ ( k ) p ] ∈ Mat m,p ( R ) , (cid:16) Ψ ( k ) p (cid:17) T Ψ ( k ) p = I p . An important observation now is that matrices Φ ( k ) p (resp. Ψ ( k ) p ) directly definean ordered p -tuple of orthonormal vectors in R n (resp. R m ), that is a point on the compact Stiefel manifold S t c ( p, n ) := { Ordered orthonormal p -tuple of vectors in R n } . To obtain a Space-Time POD interpolation (instead of a spatial POD interpolationfollowed by Galerkin approach), we finally adopted the following strategy, whendealing with a parameter value (cid:101) λ :(1) Define a curve on the compact Stiefel manifold corresponding to the spatialpart λ (cid:55)→ Φ ( λ ) ∈ S t c ( p, n )obtained using the already known interpolation algorithm on Grassmannmanifold.(2) In the same way, define a curve on the compact Stiefel manifold correspond-ing to the temporal part λ (cid:55)→ Ψ ( λ ) ∈ S t c ( p, m ) . (3) Construct an interpolated curve λ (cid:55)→ S ( λ ) passing through the POD of mode p snapshot matrices S ( k ) p , in order to obtain an interpolation of a ROM matrix (cid:101) S := S ( (cid:101) λ ).In the next subsections, we give all important details to obtain such an interpolatedcurve λ (cid:55)→ S ( λ ). First, in subsection 2.1 we explain how to compute on Grass-mann manifolds using their Riemannian structure to obtain explicit formulae for the geodesics defining normal coordinates . From this explicit formulae, we can deducein subsection 2.2 a target algorithm in order to define the curves λ (cid:55)→ Φ ( λ ) ∈ S t c ( p, n ) , λ (cid:55)→ Ψ ( λ ) ∈ S t c ( p, m ) O. FRIDERIKOS, M. OLIVE, E. BARANGER, D. SAGRIS, AND C. DAVID on compact Stiefel manifolds. The question on how to define an interpolated curvefor matrices S ( k ) p will then be addressed in section 3.2.1. Riemannian geometry on Grassmann manifolds.
We will summarize nowsome essential results about Grassmann manifolds. Such manifolds are in fact com-plete Riemannian manifolds [24], meaning for instance that we can define the length of a curve. Moreover, we can always construct a curve of the shortest length betweentwo points, which is called a geodesic , and it will be the starting point to define nor-mal coordinates via the exponential and logarithm map (Definition 2.6 and 2.7). Aswe cannot do direct computations on Riemann manifolds, normal coordinates enableus to obtain formulae of curves, such as the Lagrangian polynomials. Note finallythat a rigorous mathematical background of all of this is given in [41].After we give a definition of the Grassmann manifold and how to represent itspoints with matrices, we propose to define the tangent plane using matrix represen-tative, to have formulae for a scalar product, given by (3). From this, we deduce aclassical expression for geodesics (Theorem 2.3).Let p ≤ n be two non-zero integers and G ( p, n ) the Grassmann manifold of p -dimensional subspaces in R n . In fact, Grassmann manifolds are special cases of quotient manifolds , meaning that a point on such a manifold can have many repre-sentatives . Let us consider indeed a p -dimensional linear subspace V of R n . Such asubspace can be defined using any ordered set of p independent vectors v , . . . , v p in R n , encoded into a full rank matrix M := [ v , . . . , v p ] ∈ Mat n,p ( R ) . Any other basis v (cid:48) , . . . , v (cid:48) p of V will then lead to another full rank matrix M (cid:48) := [ v (cid:48) , . . . , v (cid:48) p ] ∈ Mat n,p ( R ) , and we necessary have M (cid:48) = M Pwhere P ∈ GL( p ) is some invertible matrix in Mat p,p ( R ). From all this, we deducethat the point m := V ∈ G ( p, n ) is represented by the infinite set of matrices { M P , P ∈ GL( p ) } . Now, the ordered set of p independent vectors in R n and thus the set of full rankmatrices in Mat n,p ( R ) define the Stiefel manifold (see Figure 1) S t ( p, n ) := { M = [ v , . . . , v p ] ∈ Mat n,p ( R ) , rg( M ) = p } so that we obtain a natural map from such Stiefel manifold and the Grassmannmanifold G ( p, n ) (see Figure 2): M = [ v , . . . , v p ] ∈ S t ( p, n ) (cid:55)→ m = { M P , P ∈ GL( p ) } . T INTERPOLATION OF PARAMETRIZED RIGID-VISCOPLASTIC FEM PROBLEMS 9
In our situation, nevertheless, we will only focus on orthonormal bases of p -dimensionalsubspaces. Doing so, we thus consider matrices defined by orthonormal vectors, lead-ing to the so-called compact Stiefel manifold S t c ( p, n ) := (cid:8) Y ∈ Mat n,p ( R ) , Y T Y = I p (cid:9) (1)and any point m ∈ G ( p, n ) will then be represented by the infinite set { Y Q , Q ∈ O( p ) } where Y = [ y , . . . , y p ] is defined using an orthonormal basis y , . . . , y p of m . Thisdefines a surjective map π : Y ∈ S t c ( p, n ) (cid:55)→ m = π ( Y ) = { Y Q , Q ∈ O( p ) } ∈ G ( p, n )and the set of all matrices representing the same point m ∈ G ( p, n ) is called the fiber of π at m (see Figure 3 for an illustration of a fiber): π − ( m ) = { Y Q, Q ∈ O( p ) } . Remark . An important point here is that, from now on, any computationon G ( p, n ) will be done using a choice in the fibers . Nevertheless, for any point m ∈ G ( p, n ), there is no canonical way to choose an element Y ∈ π − ( m ), so anycomputation has to be independent of that choice.We need now to define the geodesics of Grassmann manifold, which can be doneonce we have defined the tangent plane at each point m ∈ G ( p, n ) and a Riemaniannmetric. Take any point m ∈ G ( p, n ) represented by a matrix Y = [ y , . . . , y p ] oforthonormal vectors, the tangent plane T m := T m G ( p, n ) is then represented by the p ( n − p ) dimensional vector spaceHor Y := (cid:8) Z ∈ Mat n,p ( R ) , Z T Y = 0 (cid:9) , (2)called the horizontal space , where Y ∈ π − ( m ). From all this, a vector v ∈ T m willbe called a velocity vector , which can be represented by a matrix Z ∈ Mat n,p ( R ) suchthat Z T Y = 0, and Z is called a horizontal lift of v . Example . Take here p = 2 and n = 5, so that G (2 ,
5) is the set of planes in a fivedimensional space. The matrices Y = − √ √ , Y (cid:48) = √ − √ −√
24 2+ √ √
24 2 −√ √ − √ are in the compact Stiefel manifold S t c (2 , m ∈G (2 , Y defined by (2) is a 6-dimensional vector space ofmatrices Z , for instance given by Z = u v u v u v − u − v − u + u − u − v + v − v , u i , v i ∈ R . Taking now velocity vectors v , v ∈ T m with respective horizontal lifts Z , Z ∈ Hor Y we define the point–wise scalar product [42, 40]: (cid:104) v , v (cid:105) m := tr (cid:0) Z T Z (cid:1) . (3)Such a Riemannian metric leads to explicit geodesics given by [40, 39]: Theorem 2.3.
Let m ∈ G ( p, n ) represented by Y ∈ S t c ( p, n ) . For any v ∈ T m withhorizontal lift given by Z in Hor Y , let Z = UΣV T be a thin SVD of Z . Then α v : t ∈ R (cid:55)→ α v ( t ) := π (cid:2) ( YV cos( t Σ ) + U sin( t Σ )) V T (cid:3) ∈ G ( p, n ) (4) is the unique maximal geodesic such that α v (0) = m and initial velocity ˙ α v (0) := ∂α v ( t ) ∂t | t =0 = v. Remark . Up to our knowledge, there is no proof that Y ( t ) := ( YV cos( t Σ ) + U sin( t Σ )) V T ∈ S t c ( p, n ) . (5)In fact, this follows by direct computation. Indeed, Z = UΣV T being a thin SVD,we have V ∈ O( p ) and Z T Y = VΣU T Y = 0 = ⇒ ΣU T Y = 0so that sin( t Σ ) U T Y = 0 and Y T U sin( t Σ ) = 0 . Finally, we have: Y T ( t ) Y ( t ) = V cos ( t Σ ) + sin( t Σ ) U T Y (cid:124) (cid:123)(cid:122) (cid:125) =0 V cos( t Σ )+cos( t Σ ) V T Y T U sin( t Σ ) (cid:124) (cid:123)(cid:122) (cid:125) =0 + sin ( t Σ ) V T which concludes the proof. T INTERPOLATION OF PARAMETRIZED RIGID-VISCOPLASTIC FEM PROBLEMS 11
Remark . In many cases, formulas of the geodesic do not use the right multi-plication by V T , as for instance in [40, 21]. Of course, as V being in O( p ) bothmatrices ( YV cos( t Σ ) + U sin( t Σ )) V T and YV cos( t Σ ) + U sin( t Σ )define the same point on G ( p, n ). Now, the choice of such right multiplication in (5)is related to the choice of the horizontal lift Z = UΣV T . Indeed, taking back thepath given by (5), we have˙ Y ( t ) = ( − YVΣ sin( t Σ ) + UΣ cos( t Σ )) V T = ⇒ ˙ Y (0) = UΣV T = Z which corresponds to the choice of the horizontal lift for velocity vector v ∈ T m .A consequence of Theorem 2.3 is an explicit formula for the exponential map [40,21] (see Figure 4): Definition 2.6.
Let m ∈ G ( p, n ) be represented by Y ∈ S t c ( p, n ). For any velocityvector v ∈ T m with horizontal lift Z ∈ Hor Y , take Z = UΣV T to be a thin SVD of Z . Then we define the exponential mapExp m : T m −→ G ( p, n ) ,v (cid:55)→ Exp m ( v ) := π (cid:2) ( YV cos( Σ ) + U sin( Σ )) V T (cid:3) = α v (1) . Now, it is possible to define directly some inverse map of the exponential map,called the logarithm map [40], but only locally . For any m and Y in its fiber, let usfirst define the open spaceU m := { m ∈ G ( p, n ) , Y T Y is invertible , Y ∈ π − ( m ) } . (6)Then we have: Definition 2.7 (Logarithm map on Grassmannian manifold) . Let m ∈ G ( p, n ) berepresented by a matrix Y ∈ S t c ( p, n ). For any point m in the open space U m represented by a matrix Y ∈ S t c ( p, n ), define a thin SVD Y (cid:0) Y T Y (cid:1) − − Y = UΣV T . Then the logarithm Log m ( m ) ∈ T m is the velocity vector in T m with horizontal lift Z = U arctan( Σ ) V T ∈ Hor Y . Remark . The logarithm map is only defined on some open set U m . This meansthat for any point m / ∈ U m , the associated matrix Y T Y is not invertible, so thatthe computation of Y (cid:0) Y T Y (cid:1) − − Y can not be done. Note finally that such an open set is strongly related to the cut-locus of a Grassmann manifold [43]. Target Algorithm on compact Stiefel manifolds.
All the mathematicalbackground summarized in subsection 2.1 can be used to obtain an interpolationcurve between points m , . . . , m N on Grassmann manifold G ( p, n ) [20, 21], whereeach point m i corresponds to a parameter value λ i . Indeed, once a reference point m i ∈ { m , . . . , m N } is chosen (see Figure 4): • We use the logarithm map Log m i to linearize , i.e., meaning we define ve-locity vectors v i := Log m i ( m i ) on the vector space T m i . • We obtain an interpolation curve λ (cid:55)→ v ( λ ) between vectors v i , using forinstance Lagrangian polynomial, and thus v ( λ i ) = v i , ∀ i = 1 , . . . , N. • Taking the exponential map Exp m i , we obtain back an interpolation curve λ (cid:55)→ m ( λ ) := Exp m i ( v ( λ ))between the points m , . . . , m N on G ( p, n ), so that m ( λ i ) = m i , ∀ i = 1 , . . . , N. We propose here to define curves on the compact Stiefel manifold S t ( p, n ) insteadof the ones defined on the Grassmann manifold G ( p, n ). The starting point is a setof matrices Y , . . . , Y N in the compact Stiefel manifold S t ( p, n ), corresponding toparameter values λ , . . . , λ N . Once a reference parameter value λ ı has been chosen,we obtain a curve λ (cid:55)→ Y ( λ )where in general, Y ( λ i ) (cid:54) = Y i . As a consequence, such a curve will not be an interpolation curve between thematrices Y , . . . , Y N (see Remark 2.11). Before doing so, and to obtain well-definedcurves, we need to make a specific definition: Definition 2.9 (Genericity) . A matrix is said to be generic if all its non–zerosingular values are distinct. The set of generic matrices in Mat n,p ( R ) is denotedMat n,p ( R ).For any generic matrix M ∈ Mat n,p ( R ), we know that its thin SVD M = UΣV T is well defined. Indeed, taking σ > . . . > σ p to be its ordered singular values, wecan write M = p (cid:88) i =1 σ i u i v Ti (7)where u i (resp. v i ) is a left singular vector associated to σ i (resp. a right singularvector). All singular values being distinct, the only other possibility is to consider T INTERPOLATION OF PARAMETRIZED RIGID-VISCOPLASTIC FEM PROBLEMS 13 singular vectors (cid:15) i u i and (cid:15) i v i , with (cid:15) i = ±
1, so that the decomposition (7) remainsthe same. We thus deduce that the target Algorithm below is well defined:
Algorithm 2.10 (Target algorithm) . • Inputs : – Matrices Y , . . . , Y N in S t c ( p, n ) , corresponding to parameter values λ < . . . < λ N . – A reference parameter value λ i with i ∈ { , . . . , N } . – A parameter value λ . • Output : A matrix Y ( λ ) ∈ S t c ( p, n ) . (1) Define Z i := and for each k ∈ { , . . . , N } with k (cid:54) = i compute a thinSVD of the generic matrix Y k (cid:0) Y Ti Y k (cid:1) − − Y i = U k Σ k V Tk and define Z k := U k arctan( Σ k ) V Tk , with assumption Z k ∈ Mat n,p ( R ) . (2) Define an n × p matrix and compute a thin SVD Z ( λ ) := N (cid:88) i =1 (cid:89) i (cid:54) = j λ − λ j λ i − λ j Z i = U ( λ ) Σ ( λ ) V ( λ ) T , with assumption Z ( λ ) ∈ Mat n,p ( R ) . (3) Define the n × p matrix in S t c ( p, n ) (see Remark 2.5): Y ( λ ) := [ Y i V ( λ ) cos ( Σ ( λ )) + U ( λ ) sin ( Σ ( λ ))] V ( λ ) T . (8) Note: cos and sin act only on diagonal entries.
In this algorithm, as already noticed and following the assumptions of genericity,the matrices Z k , Z ( λ ) and Y ( λ ) do not depend on the choice of matrices in theassociated thin SVD. Remark . Using this target Algorithm to parameter value λ := λ k leads tosome matrix Y ( λ k ) generally different from Y k (except for k = i ). Thus, such analgorithm computed on compact Stiefel manifold do not produce an interpolation onthe points Y , . . . , Y N (see Figure 3). Indeed, to represent an interpolation curvebetween these points means that if we consider the parameter value λ = λ k (with k ∈ { , . . . , N } ) as input in the algorithm, one should expect to return as output Y ( λ k ) (given by (8)) the initial matrix Y k , which is not the case in general.Nevertheless, matrices Y ( λ k ) and Y k define the same point on the Grassmannmanifold G ( p, n ), meaning that they both define an orthonormal basis of the samesubspace m k (see Remark 2.8). As a consequence, a projection matrix onto thesubspace m k is given by Y ( λ k ) T Y ( λ k ) or equivalently by Y Tk Y k . Example . Take for instance the compact Stiefel manifold S t c (2 , Y := , Y := √ √ √ √ − √ √ , Y := √ − √ √ √ √ √ √ √ which correspond respectively to λ = 15, λ = 22 and λ = 27. Choosing thereference parameter value to be λ and following the target Algorithm 2.10 we obtain Y ( Y T Y ) − − Y = − − , Y ( Y T Y ) − − Y = Taking λ = λ and λ = λ as inputs in the algorithm, we finally obtain the matrices(with computation done using 5 digits): Y ( λ ) = . . . . . − . − . . (cid:54) = Y , Y ( λ ) = . − . − . . . . − . . . . (cid:54) = Y . Space-Time Interpolation on compact Stiefel manifolds
As already noticed, POD is extracting the optimal space structures and the as-sociated temporal modes. An important property is that the spatial and temporalorthogonal modes are coupled : each space component is associated with a temporalcomponent partner and there is a one-to-one correspondence between both spaces.Taking advance of this decomposition into orthogonal modes, it is natural to trya
Space-Time interpolation on compact Stiefel manifolds based on the target Algo-rithm 2.10, instead of an interpolation of the space part alone, followed by a Galerkinapproach as is classically done [20, 21].As a starting point, take a set of snapshot matrices S (1) , . . . , S ( N ) , where eachmatrix S ( k ) ∈ Mat n,m ( R ) corresponds to a given parameter value λ k ∈ R , with λ < . . . < λ N and n = 3 N s corresponding to the spatial part, while m = N t corresponds to the temporal part. For a given mode p ≤ N t , our goal is to T INTERPOLATION OF PARAMETRIZED RIGID-VISCOPLASTIC FEM PROBLEMS 15 (1) Extract in a unique way a POD of mode p of each matrix S ( k ) , so that wehave a well defined map S ( k ) ∈ Mat n,m ( R ) (cid:55)→ S ( k ) p ∈ Mat n,m ( R ) . (2) Obtain for each S ( k ) p ∈ Mat n,m ( R ) a unique matrix Φ ( k ) p ∈ S t c ( p, n ) for thespatial part and another unique matrix Ψ ( k ) p ∈ S t c ( p, m ) for the temporalpart.(3) Use the target Algorithm 2.10 on matrices Φ ( k ) p first, and then on matrices Ψ ( k ) p , in order to obtain two curves λ (cid:55)→ Φ ( λ ) , λ (cid:55)→ Ψ ( λ ) (9) which are not interpolated curves , as in general Φ ( λ k ) (cid:54) = Φ ( k ) p and Ψ ( λ k ) (cid:54) = Ψ ( k ) p (see Remark 2.11).(4) Define an interpolation curve λ (cid:55)→ S ( λ ) between matrices S (1) p , . . . , S ( N ) p ,using curves obtained by (9).We now detail two key points: the first concerns a new type of SVD, called orientedSVD, which allows defining the matrices Φ ( k ) p and Ψ ( k ) p in a unique way. Finally, wewill explain in subsection 3.2 how to construct the curve λ (cid:55)→ S ( λ ), which requiresthe introduction of a mixed part.3.1. Oriented SVD on generic matrices.
As already noticed in section 2, anycomputation of a POD of mode p of a matrix S ∈ Mat n,m ( R ) can be obtainedfrom a SVD. Suppose now that S is of rank r ≥ p . Any SVD of S with singularvalues σ > . . . > σ r leads to spatial orthonormal vectors φ , . . . , φ r in R n (theleft singular vectors) and temporal orthonormal vectors ψ , . . . , ψ r in R m (the rightsingular vectors). A POD of mode p then writes S p = Φ p Σ p Ψ Tp , Φ p := [ φ , . . . , φ p ] , Σ p := diag( σ , . . . , σ p ) , Ψ p := [ ψ , . . . , ψ p ] . (10)Now, because of sign indeterminacy of the spatial vectors φ i and temporal vectors ψ i , the matrices Φ p , Ψ p are not uniquely defined.To overcome this difficulty, we need to introduce a new SVD so that, under theassumption of genericity (see Definition 2.9), the matrices Φ p and Ψ p given by (10)can be well-defined.The main idea of the new SVD introduced here is to make an intrinsic choice onthe orientation for each space and temporal vector. Indeed, for each spatial vector φ , only two choices can occur: φ or − φ (thus inducing a choice on the associatedtemporal vector). A choice of orientation is then made as follows. Taking the column vectors S = [ s , . . . , s m ] and s to be the first column vector such that thescalar product (cid:104) s , φ (cid:105) is non zero, we impose the sign taking (cid:104) s , φ (cid:105) > S to choose orientation: Lemma 3.1.
Let us consider s , . . . , s m ∈ R n to be the column vectors of S ∈ Mat n,m ( R ) and take φ ∈ R n to be a unit spatial vector of S , associated with a non–zero singular value σ . Then, there exists i ∈ { , . . . , m } such that (cid:104) s i , φ (cid:105) = s Ti φ (cid:54) = 0 . From this, for any unit spatial vector φ ∈ R n of S , let us define s ( φ ) to be thefirst column vector s i in S = [ s , . . . , s n ] such that (cid:104) φ, s i (cid:105) (cid:54) = 0: s ( φ ) := s i , i := min { j, (cid:104) s j , φ (cid:105) (cid:54) = 0 } . (11)Any spatial eigenvector can therefore have a specific orientation: Definition 3.2 (Oriented eigenvectors) . Let S ∈ Mat n,m ( R ) and φ ∈ R n a unitspatial vector associated to a non–zero singular value σ . Then φ is said to be oriented if (cid:104) s ( φ ) , φ (cid:105) > Lemma 3.3 (Oriented SVD) . Let S ∈ Mat n,m ( R ) ( m ≤ n ) of rank r such that allits non-zero singular values are distinct. Then, there exists one and only one coupleof matrices Φ = [ φ , . . . , φ r ] ∈ Mat n,r ( R ) , Ψ = [ ψ , . . . , ψ r ] ∈ Mat m,r ( R ) (12) such that (cid:104) φ i , φ j (cid:105) = (cid:104) ψ i , ψ j (cid:105) = δ ij , S = ΦΣΨ T , Σ := Diag ( σ , . . . , σ r ) ∈ Mat r,r ( R ) (13) and φ i are oriented spatial unit eigenvectors: (cid:104) s ( φ i ) , φ i (cid:105) > with s ( φ i ) defined by (11) . Such a decomposition is called an oriented SVD .Proof. First, any couple ( φ, ψ ) of spatial–temporal unit eigenvector for S is definedmodulo ±
1, and ψ is obtained in a unique way from φ .Let us suppose now we do not have uniqueness, so that there exist two unit spatialvectors φ and φ (cid:48) associated to σ such that (cid:104) s ( φ ) , φ (cid:105) > (cid:104) s ( φ (cid:48) ) , φ (cid:48) (cid:105) > . We necessary have φ (cid:48) = − φ and s ( φ ) = s ( φ (cid:48) ) so we deduce that (cid:104) s ( φ (cid:48) ) , φ (cid:48) (cid:105) > = −(cid:104) s ( φ ) , φ (cid:105) > (cid:3) T INTERPOLATION OF PARAMETRIZED RIGID-VISCOPLASTIC FEM PROBLEMS 17
We give now an algorithm to obtain such an oriented SVD:
Algorithm 3.4 (Oriented SVD) . • Inputs : m ≤ n and S ∈ Mat n,m ( R ) ofrank r . • Output : Unique matrices Φ and Ψ for an oriented SVD of S . (1) Compute a SVD of S so that to obtain spatial unit vectors φ , . . . , φ r andtemporal unit vectors ψ , . . . , ψ r . (2) Consider the column vectors s , . . . , s m of S . (3) For i = 1 , . . . , r define ε i := (cid:104) s ( φ i ) , φ i (cid:105)|(cid:104) s ( φ i ) , φ i (cid:105)| where s ( φ i ) is the first column vector s of S such that (cid:104) φ i , s (cid:105) (cid:54) = 0 , see (11) . (4) For i = 1 , . . . , r , make sign replacement φ i ← ε i φ i , ψ i ← ε i ψ i . Example . Assume the rank 3 matrix S = − − − = [ s , s , s ]where a unit spatial vector corresponding to the largest singular value is given by(with 5 digits) φ = − . . . − . − . and we can check that s ( φ ) = s with (cid:104) φ , s (cid:105) < − φ insteadof φ , and so on.3.2. Space–Time interpolation algorithm.
In this subsection, we define a Space–Time interpolation on any family of POD of mode p taken from generic snapshotmatrices (see an overview in Figure 5). Such interpolation captures both the spatialand temporal part of such matrices, which is necessary from the point of view ofmechanical equations, but we will also need to define a specific mixed part of eachPOD (see lemma 3.6). Take back parameter values λ < . . . < λ N , corresponding to snapshot matrices S (1) , . . . , S ( N ) in Mat n,m ( R ), with n = 3 N s and m = N t . To make use of the orientedSVD, let us suppose: Genericity assumption:
All snapshot matrices S (1) , . . . , S ( N ) have distinct non–zero singular values.Take now p to be some integer (less or equal than the minimum rank of all matrices S ( k ) ). Using the oriented SVD given by Algorithm 3.4, we can consider a POD ofmode p on each matrix S ( k ) : S ( k ) p := Φ ( k ) p Σ ( k ) p Ψ ( k ) p T ∈ Mat n,m ( R ) (15)where Σ k corresponds to singular values, and Φ ( k ) p as well as Ψ ( k ) p uniquely definepoints in a compact Stiefel manifold: Φ ( k ) p := [ φ ( k )1 , . . . , φ ( k ) p ] ∈ S t c ( p, n ) , Ψ ( k ) p := [ ψ ( k )1 , . . . , ψ ( k ) p ] ∈ S t c ( p, m ) , (16)Recall that in previous equation, φ ( k )1 , . . . , φ ( k ) p (resp. ψ ( k )1 , . . . , ψ ( k ) p ) correspond tospatial oriented eigenvectors (resp. temporal ones) of S ( k ) .Using the target Algorithm 2.10 first for the spatial matrices Φ ( k ) p and then forthe temporal matrices Ψ ( k ) p , we obtain two curves λ (cid:55)→ Φ ( λ ) , λ (cid:55)→ Ψ ( λ ) . Now, our goal is to produce an interpolation curve between the matrices S ( k ) p , takinginto account both spatial and temporal curves defined above. Such a curve is givenby λ (cid:55)→ S ( λ ) := Φ ( λ ) M ( λ ) Ψ ( λ ) T with S ( λ k ) = S ( k ) p . (17) Lemma 3.6.
For a curve defined (17) to be an interpolation curve between thematrices S ( k ) p , we necessary have M ( λ k ) = Φ ( λ k ) T S ( k ) p Ψ ( λ k ) . (18) Proof.
To satisfy (17), we necessary have S ( λ k ) = Φ ( λ k ) M ( λ k ) Ψ ( λ k ) T = S ( k ) p , (19)where, by construction, we have Φ ( λ k ) ∈ S t c ( p, n ) and Ψ ( λ k ) ∈ S t c ( p, m ) (see targetAlgorithm 2.10). We thus have Φ ( λ k ) T Φ ( λ k ) = Ψ ( λ k ) T Ψ ( λ k ) = I p , so that by the left and right multiplication of (19) we obtain formula (18) for themixed part. (cid:3) Now we have:
T INTERPOLATION OF PARAMETRIZED RIGID-VISCOPLASTIC FEM PROBLEMS 19
Lemma 3.7.
Let λ (cid:55)→ M ( λ ) ∈ Mat p,p ( R ) be any interpolated curve between themixed part matrices M k := Φ ( λ k ) T S ( k ) p Ψ ( λ k ) ∈ Mat p,p ( R ) so that M ( λ k ) = M k for k = 1 , . . . , N . Then, using the curve λ (cid:55)→ Φ ( λ ) (resp. λ (cid:55)→ Ψ ( λ ) ) defined by the target Algorithm 2.10 applied on the matrices Φ ( k ) p (resp. Ψ ( k ) p ), the curve κ : λ (cid:55)→ Φ ( λ ) M ( λ ) Ψ ( λ ) T (20) is an interpolated curve between the matrices S (1) p , . . . , S ( N ) p , so that κ ( λ k ) = S ( k ) p foreach k = 1 , . . . , N .Proof. We need to check that κ ( λ k ) = S ( k ) p for each k = 1 , . . . , N . Now: κ ( λ k ) = Φ ( λ k ) M ( λ k ) Ψ ( λ k ) T = Φ ( λ k ) M k Ψ ( λ k ) T = Φ ( λ k ) Φ ( λ k ) T S ( k ) p Ψ ( λ k ) Ψ ( λ k ) T = Φ ( λ k ) Φ ( λ k ) T Φ ( k ) p Σ ( k ) p ( Ψ ( k ) p ) T (cid:124) (cid:123)(cid:122) (cid:125) S ( k ) p Ψ ( λ k ) Ψ ( λ k ) T where Φ ( λ k ) Φ T ( λ k ) corresponds to the projection matrix on the subspace m k := π ( Φ ( λ k )) = π ( Φ k ) (see Remark 2.11) so that Φ ( λ k ) Φ T ( λ k ) Φ ( k ) p = Φ ( k ) p (21)and the same being true for the temporal part, we obtain the proof of the lemma. (cid:3) The Space–Time interpolation algorithm is now given by:
Algorithm 3.8 (Space–Time interpolation) . • Inputs : – Generic matrices S (1) , . . . , S ( N ) in Mat n,m ( R ) ( m ≤ n ), correspondingto parameter values λ < . . . < λ N . – A reference parameter value λ i with i ∈ { , . . . , N } . – A mode p ≤ m . – A parameter value (cid:101) λ . • Output : A matrix (cid:101) S ∈ Mat n,m ( R ) . (1) Compute an oriented SVD on each matrix S ( k ) and write a POD of mode p S ( k ) p := Φ p Σ ( k ) p ( Ψ ( k ) p ) T ∈ Mat n,m ( R ) with Φ ( k ) p ∈ S t c ( p, n ) and Ψ ( k ) p ∈ S t c ( p, m ) uniquely defined. (2) Consider the target Algorithm 2.10 applied to the spatial parts Φ (1) p , . . . , Φ ( N ) p ,reference parameter value λ i and each of the N +1 parameter values λ , . . . , λ N , (cid:101) λ ,so from (8) we can define matrices in S t c ( p, n ) : Φ ( λ k ) := Y ( λ k ) , Φ ( (cid:101) λ ) := Y ( (cid:101) λ ) . (22)(3) Consider the target Algorithm 2.10 applied to the temporal parts Ψ (1) p , . . . , Ψ ( N ) p reference parameter value λ i and each of the N +1 parameter values λ , . . . , λ N , (cid:101) λ ,so from (8) we can define matrices in S t c ( p, m ) : Ψ ( λ k ) := Y ( λ k ) , Ψ ( (cid:101) λ ) := Y ( (cid:101) λ ) . (23)(4) For each k = 1 , . . . , N , define the square matrix of the mixed part M k := Φ ( λ k ) T S ( k ) p Ψ ( λ k ) ∈ Mat p,p ( R ) . (24)(5) Use a standard interpolation on square matrices M , . . . M N , for instance: M ( (cid:101) λ ) := N (cid:88) i =1 (cid:89) i (cid:54) = j (cid:101) λ − λ j λ i − λ j M i (25)(6) Using the spatial part Φ ( (cid:101) λ ) ∈ S t c ( p, n ) from (22) , the temporal part Ψ ( (cid:101) λ ) ∈S t c ( p, m ) from (23) , and the mixed part M ( (cid:101) λ ) ∈ Mat p,p ( R ) from (25) , theinterpolated snapshot matrix corresponding to (cid:101) λ is finally given by (cid:101) S := Φ ( (cid:101) λ ) M ( (cid:101) λ ) Ψ ( (cid:101) λ ) T ∈ Mat n,m ( R ) . Rigid-Viscoplastic FEM Formulation
The main defining characteristic of the RVP formulation is that it neglects theelasticity effects. This idealization is based on the fact that elastic componentsof strain remain small as compared with irreversible strains. This means that theadditive decomposition of the total strain-rate tensor ˙ ε ij = ˙ ε eij + ˙ ε pij simplifies to˙ ε ij = ˙ ε pij , where ˙ ε eij is the elastic component of the strain-rate tensor, ˙ ε pij is theplastic component and ˙ ε ij is the total strain-rate tensor. Therefore, the RVP formu-lation turns out to be very similar to fluid flow problems, and it is also called flowformulation [44]. Although it is not possible to calculate the residual stresses andthe spring-back effect, the flow formulation presents several advantages. Unlike theelastoplastic FEM, the RVP formulation, even though more approximate, is morestable, simpler to be implemented in computer codes, and can use relatively largertime increments, thus improving the computational efficiency. A thorough overviewof the foundation of the theory can be found in [34, 1]. T INTERPOLATION OF PARAMETRIZED RIGID-VISCOPLASTIC FEM PROBLEMS 21
Governing Field Equations.
Classical rigid viscoplastic problems considerthe plastic deformation of an isotropic body occupying a domain Ω ⊂ R . The do-main Ω and its boundary ∂ Ω represent the current configuration of a body accordingto the
Updated Lagrangian formulation. The governing equations that have to besatisfied are:(a) Equilibrium condition: σ ij,j = 0(b) Compatibility conditions: ˙ ε ij = 12 ( v i,j + v j,i )(c) Yield criterion: ¯ σ := (cid:18) σ (cid:48) ij σ (cid:48) ij (cid:19) = ¯ σ (¯ ε, ˙¯ ε, T )(d) Constitutive equations: σ (cid:48) ij = 23 ¯ σ ˙¯ ε ˙ ε ij , ˙¯ ε = (cid:18)
23 ˙ ε ij ˙ ε ij (cid:19) (26)(e) Incompressibility condition: ˙ ε v := ˙ ε kk = 0(f) Boundary conditions: v = ˆ v on ∂ Ω v F = ˆ F on ∂ Ω F friction and contact on ∂ Ω c In the above equations σσσ = ( σ ij ) is the stress tensor, ˙ εεε = ( ˙ ε ij ) is the strain ratetensor, v i are velocity components, ¯ σ is the effective stress, ˙¯ ε is the second invariantof ˙ εεε called effective strain rate, and σσσ (cid:48) = ( σ (cid:48) ij ) is the deviatoric stress tensor definedby σ (cid:48) ij = σ ij − δ ij σ kk / ∂ Ω consistsof three distinct parts: over ∂ Ω v velocity conditions are prescribed (essential bound-ary conditions), ∂ Ω F is the part where the traction conditions are imposed in theform of nodal point forces (natural boundary conditions), while the boundary con-ditions along ∂ Ω c are mixed, and neither the velocity nor the force can be described.Therefore, we have the disjoint union: ∂ Ω = ∂ Ω v ∪ ∂ Ω F ∪ ∂ Ω c (27) Variational form.
In a variational formulation, the functional Π (energy rate)is defined by an integral form in accordance with the virtual work-rate principleΠ( v ) := (cid:90) Ω ¯ σ ˙¯ εdV − (cid:90) ∂ Ω F F i v i dS (28)where the first term in (28) represents the internal deformation work-rate, whereasthe second term represents the work-rate done by the external forces. F i denotesprescribed surface tractions on the boundary surface ∂ Ω F . Recalling the MarKov-Hill [45, 46] variational principle, among all virtual (admissible) continuous andcontinuously differentiable velocity fields v i satisfying the conditions of compati-bility and incompressibility, as well as the velocity boundary conditions, the realvelocity field gives to the functional Π a stationary value, i.e., the first-order varia-tion vanishes. Moreover, in order to relax the incompressibility constraint condition˙ ε v = ˙ ε kk = 0 on an admissible velocity field, a classical penalized form is used δ Π := (cid:90) Ω ¯ σδ ˙¯ εdV + 12 (cid:90) Ω K ˙ ε v δ ˙ ε v dV − (cid:90) ∂ Ω F F i δv i dS = 0 (29)where K is a large positive constant which penalizes the dilatational strain-ratecomponent. It can be shown that the mean stress is σ m = K ˙ ε kk . Remark . A limitation of the Updated Lagrangian method for large deforma-tion problems is the excessive element distortion. To this end, remeshing processesare necessary to simulate unconstrained plastic flows. A mesh generation processis activated in case of zero or negative determinant of the Jacobian matrix, or dueto various element quality criteria. Then, a new mesh is calculated conforming tothe current state of the geometry followed by an interpolation of the state variablesbetween the old and the newly generated mesh. Thus, the information of the remap-ping process has to adequately be transferred to the ROM basis obtained using thePOD snapshot method. We remark that at this first attempt, we avoid remeshingsof the workpiece during the course of the simulation. This topic will be addressedin a future investigation.4.3.
Discretization and iteration.
The discretization of the functional followsthe standard procedure of the finite element method. Eq. (29) is expressed in termsof nodal point velocities v i and their variations δv i . Using the variational principle δ Π = M (cid:88) m =1 ∂ Π ( m ) ∂v i δv i = 0 , i = 1 , , ..., N s , (30) T INTERPOLATION OF PARAMETRIZED RIGID-VISCOPLASTIC FEM PROBLEMS 23 where δv i are arbitrary except that they must be zero to satisfy the correspondingessential boundary conditions, and M denotes the number of elements. From thearbitrariness of δv i , a set of algebraic equations (stiffness equations) are obtained ∂ Π ∂v i = M (cid:88) m =1 ∂ Π ( m ) ∂v i = 0 . (31)As the resulting algebraic equations are highly nonlinear, they linearized by theTaylor expansion near an assumed velocity field v = v as ∂ Π ∂v i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) v = v + ∂ Π ∂v i ∂v j (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) v = v ∆ v j = 0 (32)where the first factor of the second term is also known as the Jacobian of the sys-tem (Hessian matrix), and ∆ v j is a first-order correction of the velocity component v j . Solving (32) with respect to ∆ v j , the assumed velocity field is updated by theform (written in vector notation) v ( i ) = v ( i − + α (∆ v ) ( i ) (33)where 0 ≤ α ≤ i is the iteration step. The solution is obtained by the Directiteration method [34, 47] and/or by Newton-Raphson type methods. The iterationprocess is repeated until the following described convergence criteria are satisfiedsimultaneously (cid:107) ∆ v (cid:107) L (cid:107) v (cid:107) L ≤ e , (cid:13)(cid:13)(cid:13)(cid:13) ∂ Π ∂ v (cid:13)(cid:13)(cid:13)(cid:13) L ≤ e (34)namely, the velocity error norm and the norm of the residual equations, where e and e are sufficiently small specified tolerance numbers.4.4. Heat Transfer Analysis.
In the present model, a thermodynamically soundderivation is adopted using the conservation of energy − ρc ∂T∂t + k ∇ T + ξ ¯ σ ˙¯ ε = 0 (35)where ρc is the volume-specific heat of the material, ξ ¯ σ ˙¯ ε represents the work heatrate per unit volume due to plastic deformation, k is the thermal conductivity, T isthe temperature and ξ is a coefficient that presents the fraction of the deformationenergy dissipated into heat also known as the Taylor-Quinney coefficient.In a weak form, and using the divergence theorem − (cid:90) Ω ξ ¯ σ ˙¯ εδT dV + (cid:90) Ω k ∇ T δ ( ∇ T ) dV + (cid:90) Ω ρc ϑTϑt δT dV − (cid:90) ∂ Ω q n δT dS = 0 (36) where q n := k ∂T∂n (37)is the heat flux across the boundary ∂ Ω and n denotes the unit normal vector tothe boundary surface ∂ Ω.In standard finite element books, e.g. [48], it can be seen that the heat balanceequations such as (36), upon finite element discretization are reduced to the form: C ˙ T + KT = Q (38)where C is the heat capacity matrix, K denotes the heat conduction matrix, Q is the heat flux vector, T is the vector of nodal point temperatures, and ˙ T is therate of temperature increase vector of nodal points.The theory necessary to integrate (38) can be found in numerical analysis books [49,50]. It suffices to say that one-step time integration is used. The convergence of ascheme requires consistency and stability. Consistency is satisfied by a general timeintegration scheme t +∆ t T = t T + ∆ t [(1 − θ ) t ˙ T + θ t +∆ t ˙ T ] (39)where θ is a parameter varying between 0 and 1 ( θ = 0: Forward difference, θ = 1 /
2: Crank-Nicholson, θ = 2 /
3: Galerkin, θ = 1: Backward difference). Remark . Unconditional stability is obtained for θ ≥ .
5. This is important,because it is desirable to take time steps as large as the deformation formulationallows, since this is the most expensive part of the process.4.5.
Computational Procedure for Thermo-Mechanical Analysis.
For solv-ing coupled thermomechanical problems, two different approaches can be used. Inthe traditional monolithic approach, a single solver is in charge of the solution of theentire system of equations. In an alternative approach, the mechanical and thermalsolvers deal respectively with the viscoplastic flow and the thermal field equations.Thus, in the so-called staggered solution procedure used here, the state of the systemis advanced by sequentially executing and exchange information between these twosolvers [51]. The equations for the mechanical analysis and the temperature calcu-lation are strongly coupled , thereby making necessary the simultaneous solution ofthe finite element counterparts [34, 52, 53].5.
Numerical Investigations
The purpose of this section is to evaluate the performance of the ST POD interpo-lation using the velocity and temperature fields during the course of the simulation
T INTERPOLATION OF PARAMETRIZED RIGID-VISCOPLASTIC FEM PROBLEMS 25 of the forming process. As a benchmark test case, a rectangular cross-section baris compressed between two parallel flat dies under the condition of a constant shearfriction factor m at the die-workpiece interface. The initial workpiece has dimensions h = 20 mm (height) and w = 20 mm (width). Plane strain conditions are consid-ered. Due to the symmetry of the problem, only one quarter of the cross-section isanalyzed. The velocity of the upper and the lower die is set to v = 1 mm/s. Theinitial temperature of the die and the workpiece is set to T = 25 ◦ C. The bar iscompressed until a 35% reduction in height is achieved. The final simulation stateis accomplished in 7-time steps with a constant time increment ∆ t = 0 . σ ( ˙¯ ε ) = 1000 ˙¯ ε . (MPa) (40)The solution convergence is assumed when the velocity error norm and the forceerror norm (34) becomes less than 10 − . The type of element used is the linearisoparametric rectangular element with four-point integration. However, one pointintegration is used for the dilatation term, the second integral of the functionalin (29). This is known as the reduced integration scheme which imposes the volumeconstancy averaged over the linear rectangular element. The computational gridcomposed of 100 elements interconnected at N s = 121 nodes with 2 degrees of free-dom, resulting in a global stiffness matrix of size 242 × ε is chosen to be 0.01 and the penalty constant (orbulk modulus) K is set to 10 .Among the various models of friction, the one proposed in [54] is adapted to modelthe sliding contact at the tool-workpiece interface. This model allows the variationof the tangential traction with the relative velocity at the tool-workpiece interface t f = − mk v s | v s | (cid:39) − mk (cid:40) π arctan (cid:32) | v s | v (cid:33)(cid:41) v s | v s | where v s is the relative velocity in the tangential direction between the tool andthe workpiece, and v is a positive constant several orders of magnitude smallerthan v s ; m is the friction factor (0 < m <
1) and k is the material shear yield stress k = ¯ σ/ √
3. For the compression tests considered here, the relative tangential velocityat the tool-workpiece interface at the beginning of deformation is zero. The presentanalysis assumes that the friction factor remains constant throughout compression.Investigations on frictional shear stress measurements over the interface between acylindrical workpiece and a die during plastic compression are reported in [55]. Thebasic characteristics of algorithms used in the RVP FEM analysis are summarizedin Table 1.
Basic characteristics of algorithms in RVP FEM
Type of problem Two dimensional, plane strain, rigid viscoplastic ma-terial flow, isotropic, homogeneousThermomechanical prob-lem solution Loose coupling (staggered) - Backward Euler differ-ence ( θ = 1)Type of elements 4-node quadrilateral isoparametric elements, bilinearshape functionsFlow stress equation Power law: ¯ σ ( ˙¯ ε ) = c ˙¯ ε p , c, p constantsIteration method Direct, BFGS with line searchRemeshing N/ABoundary conditions Sliding friction on S c Table 1.
Numerical algorithms.
Remark . Note that during the course of the simulation we avoid remeshing of the work-piece. As discussed in [56], remeshing techniques can be taken into account provided thatmesh transfer operations are applied to the reduced-basis.5.1.
Mechanical field.
The first case for numerical illustration of the method considersthe velocity field during the simulation of the forming process using the shear friction factor m as the investigated parameter. From now on, let the parametric points corresponding tothe shear friction factor m denoted with λ for convenience with the previous sections. Forthe numerical study, the following training points are selected λ ∈ Λ t = { . , . , . } . Thechoice made here, is to use a minimum number of sampling points equi-distributed over theparametric range. The target point is set to (cid:101) λ = 0 .
3. See the FEM solutions for the trainingand target points at the final state of the computation in Figure 6.For each parametric simulation, a sequence of snapshots uniformly distributed over timeusing an increment of ∆ t = 0 . S ( i ) ∈ Mat N s ,N t ( R ) with 2 N s = 242 and N t = 7, corresponding toparameter values λ i , are associated with the nodal velocity field in x and y directions.For the parametric Space-Time interpolation, the snapshot matrix (cid:101) S of mode p cor-responding to the target point (cid:101) λ is computed via the target Algorithm 3.8. The targetAlgorithm 2.10 is applied to the spatial Φ (1) p , . . . , Φ ( N ) p and temporal parts Ψ (1) p , . . . , Ψ ( N ) p ,with reference parameter value λ i = 0 .
5. In order to assess the interpolation acuracy, thesnapshot matrix (cid:101) S is compared against the high-fidelity FEM solution by introducing thefollowing a posteriori errors. Using the interpolated and the HF-FEM snapshot matrices (cid:101) S and S FEM , respectively, the relative L -error measure is defined as e L ( (cid:101) s i ) := (cid:107) (cid:101) s i − s FEM i (cid:107) L (cid:107) s FEM i (cid:107) L , i = 1 , . . . , p ≤ N t . (41)Additionally, the relative Frobenius error norm of (cid:101) S and S FEM is defined as e F ( (cid:101) S ) := (cid:107) (cid:101) S − S FEM (cid:107) F / (cid:107) S FEM (cid:107) F . (42) T INTERPOLATION OF PARAMETRIZED RIGID-VISCOPLASTIC FEM PROBLEMS 27
The eigenvalue spectrum of snapshot matrices S ( i ) corresponding to training points λ i ∈ Λ t is exhibited in a semi-log scale in Figure 7. We can observe that the distance between thefirst and the last eigenvalue is from 5 up to 6 orders of magnitude. Moreover, the percentageof energy E ( k ) = (cid:80) ki =1 σ i / (cid:80) N t i =1 σ i captured from the POD modes is shown in Figure 8.It is evident that most of the 99 .
9% of the total energy is contained by the first two PODmodes.The relative L -error norm e L ( (cid:101) s i ) (see (41)) between the interpolated and the HF-FEMsolution for various POD modes is displayed in Figure 9. In general, the relative error for allPOD modes lie within a range of 0.0175 up to 0.038. It can be observed that the interpolatedST POD solution delivers good accuracy and is reliable enough to predict the velocity fieldfor the investigated target point. Remark . In the case of using p = 7 POD modes for the temporal basis interpolation,the Grassmannian manifold G ( p, p ) reduces to one point , so it is not relevant to use thetarget Algorithm 2.10: any new parameter value will give rise to the same matrix Ψ i in theassociated compact Stiefel manifold, corresponding to the reference point.Additionally, the position vector error e L ( (cid:101) x ( t )) = (cid:107) (cid:101) x ( t ) − x FEM ( t ) (cid:107) L at the nodal pointsis computed for p =2,3,5 and 7 POD modes, where (cid:101) x ( t ) and x FEM ( t ) denotes the positionvector of the ST POD and the high-fidelity FEM solutions, respectively, at the time incre-ments during the deformation. Figure 10 presents the local error e L ( (cid:101) x ( t )) superimposed atthe final loading state t = 0 .
35 s obtained from the high-fidelity FEM solution. Differentpatterns of the spatial error distribution can be observed concerning the number of PODmodes p . It is interesting to observe that in both cases, the maximum error is located nearthe upper-right location of the deforming workpiece.The evolution of the deformation process can be also represented using the time-displacementhistories of some selected nodes of the workpiece (Figure 11). The ST POD predictions arecompared against the high-fidelity FEM counterpart solution using p = 2 POD modes.Again, it can be observed that the interpolated ST POD solution is accurate and reliable topredict the evolution of the displacement field for the investigated target point during theforming process.For the preceding numerical investigations, the ST POD efficiency is demonstrated usinga single target point, i.e., (cid:101) λ = 0 .
3. To further assess the interpolation performance, anew target point is now considered, (cid:101) λ = 0 .
8. Interpolation is performed using the sameset of training points λ ∈ Λ t = { . , . , . } , with reference parameter value λ i = 0 . L -error norm e L ( (cid:101) s i ) for various POD modes p corresponding to target point (cid:101) λ = 0 . Temperature field.
To further investigate the performance of the proposed ST PODinterpolation, the temperature field obtained from the coupled thermomechanical simulationof the forming process is considered. Again, for the temperature field, we consider the shearfriction factor m as the investigated system parameter. The training points selected forthe mechanical field analysis are also used in this study, i.e., λ ∈ Λ t = { . , . , . } . The target point is set to (cid:101) λ = 0 .
3. For each parametric problem, snapshots are uniformlydistributed over time using an increment step size ∆ t = 0 . t = 0 .
35 s. The space-time snapshot matrices S ( i ) ∈ Mat N s ,N t ( R ) of size121 ×
7, corresponding to λ i , are associated with nodal temperatures. We will now comparethe Space-Time interpolation (see Algorithm 3.8) against the high-fidelity FEM solution.Again, for the target Algorithm 2.10 applied to the spatial Φ (1) p , . . . , Φ ( N ) p and temporalparts Ψ (1) p , . . . , Ψ ( N ) p , the reference parameter value λ i = 0 . m (represented by parameter λ ). The temperaturerises due to plastic work conversion to heat assuming a constant value for the Taylor-Quinneycoefficient ξ = 0 .
9. In all cases, the maximum temperature is located at the center of theworkpiece with values ranging from T = 89 . ◦ C up to T = 98 ◦ C.The eigenvalue spectrum of snapshot matrices S ( i ) corresponding to training points λ i ∈ Λ t is shown in a semi-log scale in Figure 14. We can observe that the distance between thefirst and the last eigenvalue of the curves is of the order of 5 up to 6 orders of magnitude.Moreover, the system energy E ( k ) = (cid:80) ki =1 σ i / (cid:80) N t i =1 σ i captured from the POD modes isshown in Figure 15. Most of the 99 .
9% of the total energy is contained by the first two PODmodes.The relative L -error norm e L ( (cid:101) s i ) (41) between the interpolated and the HF-FEM snap-shot matrices (cid:101) S and S FEM , respectively, for various modes p is shown in Figure 16. Addi-tionally, the Frobenius relative error norm (42) for the POD modes is presented in Figure 17.In general, the obtained results are found to have less than 1% relative error for POD modes p > p = 7 modes. The predictions are compared against the high-fidelitycounterpart solution, and it is difficult to distinguish differences among these plots. Itis revealed that the interpolated ST POD solution delivers good accuracy for all selectednodes.5.3. Computational complexity.
The computational cost of the ST POD interpolationscales with the computational complexity of SVD and the matrix operations in the target STAlgorithm 3.8. It is evident, that the cost of ST POD interpolation will be lower comparedto the standard POD Galerkin nonlinear approaches and even lower than the full order FEMsolution. The coupled thermomechanical FEM simulation for the target point takes 35.123seconds in wall-clock time. On the other hand, the ST interpolation for the mechanicalproblem using a ROM POD basis of mode p = 4 results in 0.147 seconds in wall-clock time.The ST interpolation for the thermal problem using a ROM POD basis of mode p = 4results in 0.153 seconds in wall-clock time. Therefore, the total ST interpolation takes 0.3seconds in wall-clock time corresponding to a time speed-up of 116.96. All experiments inthis section were implemented in Matlab and run on a 4th Generation Intel(R) Core(TM)i7-4600U CPU @ 2.10GHz, 8GB RAM, 250 GB SSD, Debian 9 x64. T INTERPOLATION OF PARAMETRIZED RIGID-VISCOPLASTIC FEM PROBLEMS 29 Conclusions
A novel non-intrusive Space-Time POD basis interpolation scheme on compact Stiefelmanifolds is developed and applied to parametric high nonlinear metal forming problems.Apart from the separate interpolation of POD spatial and temporal basis on associatedGrassmannian manifolds, an interpolation function is defined on a set of parametric snapshotmatrices. This function results from curves, which are defined on compact Stiefel manifoldsboth for space and the temporal part, and also the use of some mixed part encoded bya square matrix. This latter matrix provides a link between the interpolated space andtemporal basis for the construction of the target ROM snapshot matrix. To prove theefficiency of the method it has been used a coupled thermomechanical rigid-viscoplastic FEMformulation which is integrated into the manufacturing industry in a variety of applications.The performed numerical investigations have considered the reconstruction of the ROMsnapshot matrices both of the velocity and the temperature fields. Moreover, the errornorms of the Space-Time POD interpolated ROM models concerning the associated high-fidelity FEM counterpart solutions are validating the accuracy of the proposed interpolationscheme. In conclusion, the overall results demonstrate the potential use of the proposed STPOD interpolation scheme for near real-time parametric simulations using off-line computedROM POD databases, supporting thus manufacturing industries to accelerate design-to-production timespans, and thereby reducing costs while ensuring the design of superiorprocesses. • Y = [ u , u ] • Y = [ v , v ] u u v v Set of 2 independent vectors in R Points on Stiefel manifold S t (2 , Figure 1.
Points on Stiefel manifold. The linearly independent vec-tors in R spanning the red and blue planes correspond to points in S t c (2 , Vect( u , u ) = Vect( u , u )Vect( v , v ) • Y = [ v , v ] • Y = [ u , u ] • Y = [ u , u ] • • m = Vect( v , v ) m = Vect( u , u ) u u u u v v R Stiefel manifold S t (2 , G (2 , Figure 2.
Points on Stiefel S t (2 ,
3) and Grassmann manifold G (2 , T INTERPOLATION OF PARAMETRIZED RIGID-VISCOPLASTIC FEM PROBLEMS 31 π − ( m ) • Y • Y P • Y • Y G ( p, n ) S t c ( p, n ) • Y ( λ ) • Y ( λ ) • Y (˜ λ ) m • • m • m • ˜m Figure 3.
There is a natural projection π : S t c ( p, n ) −→ G ( p, n )from the compact Stiefel manifold S t c ( p, n ) to the Grassmannian G ( p, n ) of p -dimensional subspaces in R n which sends a p -frame tothe subspace spanned by that frame. The fiber over a given point m on G ( p, n ) is the set of all orthonormal p -frames spanning the sub-space m . Computations on S t c ( p, n ) using the target Algorithm 2.10for λ := λ k , lead to some matrix Y ( λ k ) generally different from Y k (except for the reference point), and thus do not produce an interpo-lation on the points Y , . . . , Y N . • m • Exp m ( v ) G ( n, p ) • m v Exp m ( · )Log m ( · ) Log m ( m ) T m G ( n, p ) Figure 4.
The exponential Exp m and the logarithm Log m map onthe Grassmann manifold G ( p, n ). Input • Snapshot matrices S (1) , . . . , S ( N ) corresponding to parameter values λ < . . . < λ N • A new parameter value e λ ∈ [ λ , λ N ] • A mode p S (1) . . . . . . Φ (1) Σ (1) ( Ψ (1) ) T S (1) p = Φ (1) p Σ (1) p ( Ψ (1) p ) T Φ (1) p Ψ (1) p S ( N ) Φ ( N ) Σ ( N ) ( Ψ ( N ) ) T S ( N ) p = Φ ( N ) p Σ ( N ) p ( Ψ ( N ) p ) T Φ ( N ) p Ψ ( N ) p Φ ( λ ) Ψ ( λ ) Φ ( λ N ) Ψ ( λ N ) M := Φ ( λ ) T S (1) p Ψ ( λ ) M N := Φ ( λ N ) T S ( N ) p Ψ ( λ N )1 Oriented SVD2 POD mode p e λ e S := Φ ( e λ ) M ( e λ ) Ψ ( e λ ) T Target Alg. on spatial part: Φ ( e λ )Target Alg. on temporal part: Ψ ( e λ )Interpolation on mixed part: M ( e λ ) Figure 5.
The Space-Time Algorithm.
T INTERPOLATION OF PARAMETRIZED RIGID-VISCOPLASTIC FEM PROBLEMS 33 (a)
For λ = 0 . (b) For λ = 0 . (c) For λ = 0 . (d) For λ = 0 . Figure 6.
Deformation patterns of the benchmark metal formingexample using different values for the shear friction factor m repre-sented by the parameter λ . -6 -4 -2 Figure 7.
The eigenvalue spectrum of snapshot matrices S ( i ) corre-sponding to training points λ ∈ Λ t = { . , . , . } . Figure 8.
Energy captured by the singular values of snapshot ma-trices S ( i ) corresponding to training points λ ∈ Λ t = { . , . , . } . T INTERPOLATION OF PARAMETRIZED RIGID-VISCOPLASTIC FEM PROBLEMS 35
Figure 9.
Performance of ST POD interpolation using the relative L -error norm e L ( (cid:101) s i ) for various modes p ; training points λ ∈ Λ t = { . , . , . } ; reference parameter value λ i = 0 .
5; target point (cid:101) λ =0 . Figure 10.
The position vector error e L ( (cid:101) x ( t )) = (cid:107) (cid:101) x ( t ) − x FEM ( t ) (cid:107) L of the nodal points at the final deformation state t = 0 .
35 s superim-posed on the high-fidelity FEM solution; POD modes p = { , , , } ;training points λ ∈ Λ t = { . , . , . } ; reference parameter value λ i = 0 .
5; target point (cid:101) λ = 0 . T INTERPOLATION OF PARAMETRIZED RIGID-VISCOPLASTIC FEM PROBLEMS 37
Figure 11.
Comparison of the total displacement of selected nodesagainst the high-fidelity FEM solution; training points λ ∈ Λ t = { . , . , . } ; reference parameter value λ i = 0 .
5; target point (cid:101) λ =0 .
3; POD modes p = 2. Figure 12.
Performance of ST POD interpolation using the relative L -error norm e L ( (cid:101) s i ) for various POD modes p ; training points λ ∈ Λ t = { . , . , . } ; reference parameter value λ i = 0 .
5; target point (cid:101) λ = 0 . Temperature field (C) (a)
For m = 0 . Temperature field (C) (b)
For m = 0 . Temperature field (C) (c)
For m = 0 . Temperature field (C) (d)
For m = 0 . Figure 13.
Temperature profiles at the final compression state t =0 .
35 s obtained using different values of the shear friction factor m represented by parameter λ . T INTERPOLATION OF PARAMETRIZED RIGID-VISCOPLASTIC FEM PROBLEMS 39 -4 -2 Figure 14.
The eigenvalue spectrum of snapshot matrices S ( i ) cor-responding to training points λ ∈ Λ t = { . , . , . } . Figure 15.
Energy captured by the singular values of snapshot ma-trices S ( i ) corresponding to training points λ ∈ Λ t = { . , . , . } . -3 Figure 16.
Performance of ST POD interpolation using the relative L -error norm e L ( (cid:101) s i ) for various POD modes p ; training points λ ∈ Λ t = { . , . , . } ; reference parameter value λ i = 0 .
5; target point (cid:101) λ = 0 . Figure 17.
Performance of the POD interpolation using the relativeFrobenius error norm e F ( (cid:101) S ) against the number of POD modes p ;training points λ ∈ Λ t = { . , . , . } ; reference parameter value λ i = 0 .
5; target point (cid:101) λ = 0 . T INTERPOLATION OF PARAMETRIZED RIGID-VISCOPLASTIC FEM PROBLEMS 41
20 40 60 80 1000.050.10.150.20.250.30.35
POD modes = 7
Figure 18.
Temperature evolution of selected nodal points validatedagainst the high-fidelity FEM solution; ST POD and HF-FEM solu-tions virtually coincide; training points λ ∈ Λ t = { . , . , . } ; ref-erence parameter value λ i = 0 .
5; target point (cid:101) λ = 0 .
3; POD modes p = 7. References [1] J.-L. Chenot. Recent contributions to the finite element modelling of metal forming processes.
Journal of Materials Processing Technology , 34(1-4):9–18, sep 1992.[2] Z. Gronostajski, Z. Pater, L. Madej, A. Gontarz, L. Lisiecki, A. (cid:32)Lukaszek-So(cid:32)lek, J. (cid:32)Luksza,S. Mr´oz, Z. Muskalski, W. Muzykiewicz, M. Pietrzyk, R.E. ´Sliwa, J. Tomczak, S. Wiewi´orowska,G. Winiarski, J. Zasadzi´nski, and S. Zi´o(cid:32)lkiewicz. Recent development trends in metal forming.
Archives of Civil and Mechanical Engineering , 19(3):898–941, may 2019.[3] Francisco Chinesta, Pierre Ladeveze, and El´ıas Cueto. A short review on model order reductionbased on proper generalized decomposition.
Archives of Computational Methods in Engineering ,18(4):395–404, oct 2011.[4] C. Allery, A. Hamdouni, D. Ryckelynck, and N. Verdon. A priori reduction method for solvingthe two-dimensional burgers’ equations.
Applied Mathematics and Computation , 217(15):6671–6679, apr 2011.[5] Philip Holmes, John L Lumley, Gahl Berkooz, and Clarence W Rowley.
Turbulence, coherentstructures, dynamical systems and symmetry . Cambridge university press, 2012.[6] Thibault Henri and Jean-Pierre Yvon. Convergence estimates of POD-galerkin methods forparabolic problems. In
IFIP International Federation for Information Processing , pages 295–306. Kluwer Academic Publishers, 2005.[7] Nadine Aubry. On the hidden beauty of the proper orthogonal decomposition.
Theoretical andComputational Fluid Dynamics , 2(5-6):339–352, aug 1991.[8] Kari Karhunen. Zur spektraltheorie stochastischer prozesse.
Ann. Acad. Sci. Fennicae, AI , 34,1946.[9] M. Lo`eve. Elementary probability theory. In
Probability Theory I , pages 1–52. Springer NewYork, 1977.[10] Gene H Golub and CFV Loan.
Matrix Computations, 3rd edn., vol. 1 . JHU Press, 1996.[11] Ian T Jolliffe. Springer series in statistics.
Principal component analysis , 29, 2002.[12] Herv´e Abdi and Lynne J Williams. Principal component analysis.
Wiley interdisciplinary re-views: computational statistics , 2(4):433–459, 2010.[13] J. Edward Jackson. Principal components and factor analysis: Part i principal components.
Journal of Quality Technology , 12(4):201–213, oct 1980.[14] J. Edward Jackson. Principal components and factor analysis: Part II—additional topics relatedto principal components.
Journal of Quality Technology , 13(1):46–58, jan 1981.[15] Patricia Astrid, Siep Weiland, Karen Willcox, and Ton Backx. Missing point estimation in mod-els described by proper orthogonal decomposition.
IEEE Transactions on Automatic Control ,53(10):2237–2251, nov 2008.[16] Annika Radermacher and Stefanie Reese. POD-based model reduction with empirical inter-polation applied to nonlinear elasticity.
International Journal for Numerical Methods in Engi-neering , 107(6):477–495, dec 2015.[17] Saifon Chaturantabut and Danny C. Sorensen. Nonlinear model reduction via discrete empiricalinterpolation.
SIAM Journal on Scientific Computing , 32(5):2737–2764, jan 2010.[18] R. Everson and L. Sirovich. Karhunen–lo`eve procedure for gappy data.
Journal of the OpticalSociety of America A , 12(8):1657, aug 1995.[19] Kevin Carlberg, Charbel Farhat, Julien Cortial, and David Amsallem. The GNAT method fornonlinear model reduction: Effective implementation and application to computational fluiddynamics and turbulent flows.
Journal of Computational Physics , 242:623–647, jun 2013.
T INTERPOLATION OF PARAMETRIZED RIGID-VISCOPLASTIC FEM PROBLEMS 43 [20] David Amsallem, Julien Cortial, Kevin Carlberg, and Charbel Farhat. A method for interpolat-ing on manifolds structural dynamics reduced-order models.
International journal for numericalmethods in engineering , 80(9):1241–1258, 2009.[21] Rolando Mosquera Meza.
Interpolation sur les vari´et´es grassmanniennes et applications `a lar´eduction de mod`eles en m´ecanique . PhD thesis, La Rochelle, 2018.[22] Silv`ere Bonnabel and Rodolphe Sepulchre. Riemannian metric and geometric mean for posi-tive semidefinite matrices of fixed rank.
SIAM Journal on Matrix Analysis and Applications ,31(3):1055–1070, jan 2010.[23] John M Lee. Smooth manifolds. In
Introduction to Smooth Manifolds , pages 1–31. Springer,2013.[24] Sylvestre Gallot, Dominique Hulin, and Jacques Lafontaine.
Riemannian geometry , volume 2.Springer, 1990.[25] Rolando Mosquera, , Aziz Hamdouni, Abdallah El Hamidi, and Cyrille Allery. POD basisinterpolation via inverse distance weighting on grassmann manifolds.
Discrete & ContinuousDynamical Systems - S , 12(6):1743–1759, 2019.[26] G. Muhlbach. The general neville-aitken-algorithm and some applications.
Numerische Math-ematik , 31(1):97–110, mar 1978.[27] Y. Lu, N. Blal, and A. Gravouil. Space–time POD based computational vademecums for para-metric studies: application to thermo-mechanical problems.
Advanced Modeling and Simulationin Engineering Sciences , 5(1), feb 2018.[28] M. Oulghelou and C. Allery. Non intrusive method for parametric model order reduction usinga bi-calibrated interpolation on the grassmann manifold.
Journal of Computational Physics ,426:109924, feb 2021.[29] Vilas Shinde, Elisabeth Longatte, Franck Baj, Yannick Hoarau, and Marianna Braza. Agalerkin-free model reduction approach for the navier–stokes equations.
Journal of Compu-tational Physics , 309:148–163, mar 2016.[30] C. Audouze, F. De Vuyst, and P. B. Nair. Reduced-order modeling of parameterized PDEsusing time-space-parameter principal component analysis.
International Journal for NumericalMethods in Engineering , 80(8):1025–1057, nov 2009.[31] Youngsoo Choi and Kevin Carlberg. Space–time least-squares petrov–galerkin projection fornonlinear model reduction.
SIAM Journal on Scientific Computing , 41(1):A26–A58, jan 2019.[32] Youngsoo Choi, Peter Brown, William Arrighi, Robert Anderson, and Kevin Huynh.Space–time reduced order model for large-scale linear dynamical systems with application toboltzmann transport problems.
Journal of Computational Physics , 424:109845, jan 2021.[33] Christophe Audouze, Florian De Vuyst, and Prasanth B. Nair. Nonintrusive reduced-ordermodeling of parametrized time-dependent partial differential equations.
Numerical Methods forPartial Differential Equations , 29(5):1587–1628, feb 2013.[34] Shiro Kobayashi, Shir¯o Kobayashi, Soo-Ik Oh, and Taylan Altan.
Metal forming and the finite-element method , volume 4. Oxford University Press on Demand, 1989.[35] C. H. Lee and S. Kobayashi. New solutions to rigid-plastic deformation problems using a matrixmethod.
Journal of Engineering for Industry , 95(3):865–873, aug 1973.[36] Shiro Kobayashi. Rigid-plastic finite element analysis of axisymmetric metal forming processes.
Numerical Modeling of Manuf. Process (ASME, New York, 1977) , pages 49–65, 1977.[37] ZQ Feng and G De Saxc´e. Rigid-plastic implicit integration scheme for analysis of metal form-ing.
European Journal of Mechanics, A/Solids , 15(1):51–66, 1996. [38] O. Friderikos. Two-dimensional rigid-plastic fem simulation of metal forming processes in mat-lab. Proceedings of the 4th International Conference on Manufacturing and Materials Engi-neering (ICMMEN), 3-5 October, Thessaloniki, Greece, 2011.[39] Alan Edelman, Tom´as A Arias, and Steven T Smith. The geometry of algorithms with or-thogonality constraints.
SIAM journal on Matrix Analysis and Applications , 20(2):303–353,1998.[40] P.-A. Absil, R. Mahony, and R. Sepulchre. Riemannian geometry of grassmann manifolds witha view on algorithmic computation.
Acta Applicandae Mathematicae , 80(2):199–220, jan 2004.[41] S. E. Kozlov. Geometry of the real grassmannian manifolds. parts i, ii.
Zapiski NauchnykhSeminarov POMI , 246:84–107, 1997.[42] Yung-Chow Wong. Differential geometry of grassmann manifolds.
Proceedings of the NationalAcademy of Sciences of the United States of America , 57(3):589, 1967.[43] S. Berceanu. On the geometry of complex grassmann manifold, its noncompact dual and co-herent states.
Bulletin of the Belgian Mathematical Society - Simon Stevin , 4(2):205–243, 1997.[44] Rodney Hill.
The mathematical theory of plasticity , volume 11. Oxford university press, 1998.[45] A.A. Markov.
On variational principles in the theory of plasticity . Division of Applied Mathe-matics, Brown University, 1948.[46] R Hill. A variational principle of maximum plastic work in classical plasticity.
The QuarterlyJournal of Mechanics and Applied Mathematics , 1(1):18–28, 1948.[47] S.I. Oh. Finite element analysis of metal forming processes with arbitrarily shaped dies.
Inter-national Journal of Mechanical Sciences , 24(8):479–493, jan 1982.[48] O. C. Zienkiewicz and P. N. Godbole. Flow of plastic and visco-plastic solids with specialreference to extrusion and forming processes.
International Journal for Numerical Methods inEngineering , 8(1):1–16, 1974.[49] Anthony Ralston and Philip Rabinowitz.
A first course in numerical analysis . Courier Corpo-ration, 2001.[50] Germund Dahlquist and ˚Ake Bj¨orck.
Numerical Methods in Scientific Computing, Volume I .Society for Industrial and Applied Mathematics, jan 2008.[51] C.A. Felippa and K.C. Park. Staggered transient analysis procedures for coupled mechanicalsystems: Formulation.
Computer Methods in Applied Mechanics and Engineering , 24(1):61–111,oct 1980.[52] N. Rebelo and S. Kobayashi. A coupled analysis of viscoplastic deformation and heat trans-fer—i.
International Journal of Mechanical Sciences , 22(11):699–705, jan 1980.[53] N. Rebelo and S. Kobayashi. A coupled analysis of viscoplastic deformation and heat trans-fer—II.
International Journal of Mechanical Sciences , 22(11):707–718, jan 1980.[54] CC Chen. Rigid-plastic finite-element analysis of ring compression.
Applications of numericalmethods of forming processes , 1978.[55] G.T. van Rooyen and W.A. Backofen. A study of interface friction in plastic compression.
International Journal of Mechanical Sciences , 1(1):1–27, jan 1960.[56] D. Ryckelynck. Hyper-reduction of mechanical models involving internal variables.
Interna-tional Journal for Numerical Methods in Engineering , 77(1):75–89, jan 2009.
T INTERPOLATION OF PARAMETRIZED RIGID-VISCOPLASTIC FEM PROBLEMS 45 (Orestis Friderikos)
Universit´e Paris-Saclay, ENS Paris-Saclay, CNRS, LMT - Labora-toire de M´ecanique et Technologie, 91190, Gif-sur-Yvette, France
Email address : [email protected] (Marc Olive) Universit´e Paris-Saclay, ENS Paris-Saclay, CNRS, LMT - Laboratoirede M´ecanique et Technologie, 91190, Gif-sur-Yvette, France
Email address : [email protected] (Emmanuel Baranger) Universit´e Paris-Saclay, ENS Paris-Saclay, CNRS, LMT - Labo-ratoire de M´ecanique et Technologie, 91190, Gif-sur-Yvette, France
Email address : [email protected] (Dimitris Sagris) Mechanical Engineering Department, Laboratory of ManufacturingTechnology & Machine Tools, International Hellenic University, GR-62124 SerresCampus, Greece.
Email address : [email protected] (Constantine David) Mechanical Engineering Department, Laboratory of Manufac-turing Technology & Machine Tools, International Hellenic University, GR-62124Serres Campus, Greece.
Email address ::