AN INTERPOLATION METHOD FOR DETERMINING THE FREQUENCIES OF PARAMETERIZED LARGE-SCALE STRUCTURES

Parametric Model Order Reduction (pMOR) is an emerging category of models developed with the aim of describing reduced first and second-order dynamical systems. The use of a pROM turns out useful in a variety of applications ranging from the analysis of Micro-Electro-Mechanical Systems (MEMS) to the optimization of complex mechanical systems because they allow predicting the dynamical behavior to be predicted at any values of the quantities of interest within the design space, e.g. material properties, geometric features or loading conditions. The process underlying the construction of a pROM using an SVD-based method accounts for three basic phases: a) construction of several local ROMs (Reduced Order Models); b) projection of the state-space vector onto a common subspace spanned by several transformation matrices derived in the first step; c) use of an interpolation method capable of capturing the values of the quantity of interest for one or more parameters. One of the major difficulties encountered in this process has been identified at the level of the interpolation method and can be encapsulated in the following contradiction: if the number of detailed finite element analyses is high then an interpolation method can better describe the system for a given choice of a parameter, but the computation time is higher. In this paper, a method is proposed for removing this contradiction by introducing a new interpolation method (RSDM). This method allows us to restore and make available to the interpolation tool certain natural components belonging to the matrices of the full FE model that are related, on one hand to the process of reduction and, on the other to the characteristics of a solid in the FE theory. This approach shows higher accuracy than methods used for assessing the eigenbehavior of the system. A Hexapod will be analyzed to confirm the usefulness of the RSDM.


Parametric Model Order Reduction and its Interpolation
Within the engineering community increasing efforts have been focusing on Parametric Models Order Reduction (pMOR) as a way of diminishing the computation time for simulating large systems.By interpolating among the reduced matrices obtained at specific values of 1 or more parameters -where a parameters refers, for example, to the length of a plate, the magnitude of an externally applied force, a material property, or specific boundary or initial conditions -the analyst is able to obtain very fast and accurate models that can then be used for a number of purposes.
One indicator showing how increasing attention has been spreading throughout the scientific community is shown by the application of the parametric order reduction for solving electrochemical models for simulating cyclic voltammograms in the field of Polarography [8]; for simulating 3D Micro-Electro-Mechanical-Systems (MEMS) [19]; for design optimization [9], and also for control [18].In these studies, ordinary first and second order differential equations are solved to estimate the frequency response of the system within a range of frequencies of interest.To reduce these equations, most works employ MOR techniques spanning from the Truncated Balanced Realization method [2][3][4] to the Krylov subspace method [5][6][7]10].

Problems related to the Interpolation Method
The construction of a pROM is accompanied by the use of an interpolation method for predicting the values of the quantity of interest for one or more parameters at any point in the design space.Computing of parametric stiffness and mass matrices requires the application of a suitable interpolation scheme to deliver accurate results.An interpolation scheme for the computation of parametric matrices is suitable when it: (1.) preserves the property of positive definiteness of the parametric ROMs; (2.) preserves the quantity of interest within the design space with the least error; (3.) requires the lowest computational cost with the simplest or minimal mathematical apparatus.
The first property is essential because if it is violated then the computation of, for example, the eigenfrequencies of the system under investigation yields complex numbers.An intuitive approach used for interpolating between two or more precomputed ROMs would consist of lettin a first, second or higher order interpolating polynomial passing through the values of the K or M entries for different values of the parameter of interest.Therefore, if the parameter that is used is, for example, the length of a cantilever, and three points have been chosen to explore the design space spanning from the first length, say L 1 = 200 mm to the third length L 3 = 300 mm through L 2 = 250 mm, a quadratic polynomial could be used to interpolate between points: (L 1 , K , where the superscript indicates at what length value the matrix has been computed.Although this scheme seems attractive due to the use of reduced matrices and available functions for computing interpolating polynomials, it produces ineffective results.The reason is that homologous entries, e.g., of a stiffness matrix calculated at different length values exhibits a highly nonlinear behavior.Therefore, one consequence of applying this method is that the property (1) above is violated.There is a possible way out; it would entail calculating the matrices for length values very close to each other.In this case, the parametric matrix would be readily available for any length value.The problem is that this procedure would require the interval between two adjacent matrices be so narrow that it would discourage the use of a pROM, and would therefore find no application at all for engineering purposes.
A simple way to obtain a parametric ROM for a given set of parameters entails performing a linear approximation of a ROM that lies between two precomputed ROMs obtained for different values of the parameters.These precomputed ROMs are weighted by appropriately chosen functions and the final, parametric matrix is generated by superimposing system matrices, as shown in [16].In order to improve the effectiveness of the interpolation tool, and due to the difficulty arising from describing the parametric effects for nonlinear systems a spline-based interpolation method is introduced in [17] and also in [15], where element-wise interpolation in the tangent space of the matrix manifold is employed.[15] also addresses the mode-veering problem and shows how the method can detect the narrowing of the frequency gaps between two adjacent mode shapes.
In [18], the interpolation method that is used is based on the Direct Matrix Interpolation Methodin this paper indicated as DMM -which consists of the weighted sum of the reduced matrices calculated in correspondence with specific values of the chosen parameter(s).The typical form of, for example, a parametric reduced stiffness matrix is given by where k j=1 α j (p 1 , p 2 , . . . ) = 1, α j (p i ) = δ ij for i, j = 1, . . ., k, p ∈ Π.Here, α j (p 1 , p 2 , . . . ) are the coefficients depending on the selected parameters while K j are the reduced local stiffness matrices.In [18], the DMM is employed for computating of the frequency response of first-order Linear Time Invariant dynamical systems, and two local models -i.e., two pairs of reduced K and M -are used for 1 or 2 parameters.Unfortunately, the application of this superposition of weighted matrices relies on the use of only 1 parameter for a linear interpolation, which implies the use of only two precomputed ROMs.This can strongly affect the precision with which the final parametric matrices can be employed with success in assessing the eigenbehavior of the system.In [15], a wing with store configuration is presented, and first 8 frequencies are calculated, but there is error frequency of up 10 %.
Furthermore, in the aforementioned works it have not yet pointed out the importance of recovering the components of the FE matrices associated to the process of reduction and the deformation of the solid.In order to take these components into account, it is necessary to use a different viewpoint, i.e. to proceed backwards from the ROMs to the mathematical structure that characterizes the computation of the matrices of an FE assemblage.

Aims
The first objective of this paper is to introduce a new interpolation method -named RSDM -for adapting the framework of the DMM to a different structure of stiffness and mass matrices.Instead of performing element-wise interpolation of precomputed ROMs, the proposed method transforms the interpolating matrices into a form that allows interpolation among these deformation components, and changes in the coordinates of the system, which are made visible by appropriate factorization.This process seems essential, because certain components describing the matrices of an FE assemblage have remained encapsulated after the process of reduction, and have blent with the components of the transformation matrix; while other components, e.g.those pertaining to the rigid rotation, which do not participate in the deformation of the solid, are not removed.As result, the application of a Matrix Interpolation based Method without this step remains, so to say, "blind" and incapable of discerning between the useful components.In the first step of the method, two pairs of matrices -the stiffness and mass matrices -are factorized in such a way as to rule out the useless components and restore those components that have been blended after the process of reduction.The overall computational complexity and the numerical properties take advantage of this procedure.In fact, in order to obtain very accurate pROM, especially within the broad design space of the chosen parameters, it might be necessary to run several detailed finite element analyses.Consequently, it might be time-consuming to run several simulations in order to keep the accuracy within a given tolerance.This problem can be condensed in the following contradiction: if the number of detailed finite element analyses is high (3 or more), then an interpolation method can better describe the system because it can more effectively predict the frequencies pertaining to intermediate values of the chosen parameter(s), but the computation time becomes higher.As it will be shown, instead of searching for a trade-off between these two conflicting requirements -computation time and computation accuracy -it is possible to improve the effectiveness of the method and to remove the above contradiction by making minimal changes to the Direct Matrix Method for the system.Since this solution to the above problem involves minimal changes, it is more ideal because it provides an easier way to introduce fewer resources (only a carefully manipulated SVD), while reducing the expenses (the computational time and the mathematical apparatus required by the RSDM) as much as possible.The intuitive promise of the proposed method is illustrated in Figure 1.
The second objective of the paper is to show the performance of RSDM in terms of accuracy.A numerical application to a Hexapod discretized with almost 200 thousand degrees of freedom will be presented.

Foundations and derivation of the proposed Interpolation method. An analogy with continuum mechanics
Solid mechanics is a branch of continuum mechanics that has been developed to predict the behavior of solids -e.g.changes in shape, internal forces -subjected to the action of mechanical and thermal loads.
A general deformation of a solid -using the idealization often considered in engineering dynamics that a body is rigid -can be defined as a combination of a rigid rotation and a stretch (or contraction) [13].
Once the deformation is known, the stresses induced by external mechanical and/or thermal loads can be calculated.However, these stresses are generated by the deformational component of the stretch and not of the rigid rotation.Therefore, in order to define the deformation gradient for a solid, the stretch component is separated from the rigid rotation component by means of a mathematical factorization known as Polar Decomposition [10].
In order to predict the deformation and the stress fields, most of the engineering design calculations make use of finite element techniques that are employed to solve the governing equations of a solid, both for static analysis and for dynamic analysis,.With reference to figure 2 and using the Principle of Virtual Work the equations of motion for linear elastic solids are described as follows ( [11], p.492): This can be in turns be written in the compact form where u a i , is the displacement vector at each nodal point, while b i and t i are respectively the vectors of body and traction forces.This set of linear equations can also be used to study the eigenbehavior of the system.In order to calculate the frequencies and the mode shapes of this structure, it is necessary then to solve the eigenproblem of the type that stems from solving the characteristic equation obtained after introducing the trial solution, within the system of equations ( 3).The calculation of the eigensolution to this set of equations is therefore dependent on K aibk and M aibk , which are the stiffness and mass matrices of the finite element assemblage illustrated in Figure 2. Referring to (2), these matrices and the vector of the externally applied loads have the following mathematical structure: where R and ∂ 2 R are the shape of the solid in its unloaded condition and its boundary, respectively.
Analyzing these relations, it can be observed that both the finite element stiffness and mass matrices depend on the material properties (the density and the elastic modulus tensor C ijkl ) and on the element interpolation functions or shape functions, N (x 1 , x 2 ) with x 1 , x 2 being the coordinates of a generic node with respect to a reference system O e 1 e 2 (Figure 2); while, indices a, b, i, k, l indicate the directions pertaining to the Cartesian basis, {e 1 , e 2 }.
In turns, the shape functions depend on the type of finite element and on the reciprocal displacement of the nodes of the discretized structure.That is, the entries of the stiffness and mass matrices of a finite element assemblage contain information related to the material and the deformation of the solid.
The general deformation of a solid can be defined as combination of a rigid rotation and a stretch (or contraction) [13].The process of general deformation undergone by a finite element is illustrated in part (b) of Figure 2 and it is described as a stretch followed by a rigid rotation or, which is the same thing, by a rigid rotation followed by a stretch.Once the deformation is known, the stresses can be derived.However, these stresses are generated only by the deformation component of the stretch and not of the rigid rotation; that is, a rigid rotation does not participates actively in the stretch or contraction of the solid.As a result, the matrices K aibk and M aibk are dependent only on the two components of rigid rotation and stretch.Now, without loss of generality, let us suppose that the finite element assemblage of Figure 2(a) has a large number of elements, and that we desire to reduce the size of K aibk and M aibk by means of a transformation matrix T which, by definition, is a rotation matrix.
Therefore, the full stiffness and mass matrices, in addition to containing information related to the rigid rotation and the stretch of the solid, also retain information related to the change in the coordinates system undergone after the process of reduction.It then appears evident then that, whenever a large-scale system is reduced and then analyzed for modal analysis by means of K r = T T KT , M r = T T M T , three types of components are blended to one another.
Consequently, these components are not actively exploited when applying the standard Matrix interpolation method, because the weights that are used can never distinguish between them.A direct consequence of this situation is that the resulting ROMs cannot guarantee satisfactory accuracy.
The question that emerges is therefore: How can these three components be taken out of K r and M r in order to make them visible to the interpolation tool?
The next section introduces the important concept of Polar Decomposition, which allows the useless component of rigid rotation to be ruled out, and brings out the useful component associated to the deformation.

Polar decomposition
The Polar Decomposition can be defined as the analogous, for matrices, of a complex number z = re iϑ , r ≥ 0.
If rank A = n, then S is positive definite and R is uniquely determined. 1onsidering the partition (8) where the factor V , known as the right-hand singular vector of S [12], is an orthogonal matrix.In continuum mechanics, the gradient of deformation is expressed in terms of R and S, which describe the rotation and the stretch of a solid, respectively.This is the feature that is exploited by the proposed interpolation method.The two matrices K r and M r are decomposed in terms of R and S, and since the matrix R does not participate in the deformation of the solid, it can be ruled out.Matrices K r and M r can therefore be expressed in terms of the S-component only.As several numerical experiments have shown, ruling out the rotational component of deformation improves the accuracy of results.
Therefore, the term RSDM2 stems from the process that separates the rigid Rotation from Stretch by means of a Decomposition Method in order to bring out the components of transformation and deformation.
Figure 3 depicts the situation just described.At the top, the reduction process of the stiffness matrix is shown: T is an orthogonal matrix, which allows the state variables to be projected onto the subspace spanned by its columns.Matrix K accounts for the deformation due to the reciprocal displacement of the nodes.From continuum mechanics, however, rigid rotation is ruled out by applying a Decomposition to describe the deformation through the action of the stretch which, unlike the rigid rotation, is responsible for the stress.Consequently (Figure 3 bottom), the rotation and the stretch, which, after the process of reduction are originally blended, can each be made visible by the factors of a Singular Value Decomposition (applied in this example to the stiffness matrix): an orthogonal matrix V and the matrix of singular values that define the part responsible for the deformation.Polar Decomposition can therefore be considered as the conceptual tool that, by analogy with continuum mechanics, allows us to identify which component of the rigid rotation to rule out; SVD practically allows us to compute the components of the reduced matrix pertaining to the change of coordinates and stretch.
To compute the S-component, it will therefore be necessary to apply an SVD method to the matrices of interest.The following section shows this step.

Application of the proposed Interpolation method
In order to obtain the S matrix, it is therefore necessary to carry out a Singular Value Decomposition, and to exploit only two of its factors.It is important to highlight this, because it might lead to the idea that SVD is applied to the matrices in its entirety, whereas only a part of its generated factors is in fact used.Since the SVD method is numerically stable and is implemented in different computational software programs, its use for calculating S is simple and fast.Suppose that two finite element analyses have been run, and that the matrices of mass M and stiffness K have been derived.Two couples of matrices are now available: (K 1 , K 2 ) and (M 1 , M 2 ).The first step consists of performing an SVD and extracting the matrices V and D for each of these 4 matrices.
Let us use (V k1 , V k2 , D k1 , D k2 ) and (V m1 , V m2 , D m1 , D m2 ) to refer to the 8 matrices pertaining to the SVD applied to each of the 4 matrices (K 1 , K 2 ) and (M 1 , M 2 ).
The coefficients are an essential part of this model.More specifically, for those values of the parameter(s)  at which (K 1 , M 1 ) are obtained, α i = β i = 1 while in the correspondence of those values at which (K 2 , M 2 ) are computed, α i = β i = 0, i = 1, 2, 3.The advantages of this method rely on three factors: (1.) Unlike the Direct Matrix Method, this method accounts for 6 coefficients, thus allowing the model to be tuned up efficiently.
( .)The decomposition discussed above is guaranteed to act on more components that, without decomposition, would remain encapsulated into the reduced stiffness and mass matrices; namely the component pertaining to the transformation matrix and the component associated to the shape functions of the matrices of the FE assemblage. (3.)The useless component related to the rigid rotation has been ruled out, so it does not participate actively and harmfully in the process of interpolation.

Results
We now present as a numerical example a study of a Hexapod structure with piezoactuators that was created at the Faculty of Mechanical Engineering of the Czech Technical University in Prague.The aim is to investigate the suppression of vibrations of compliant mechanical structures.This structure will show how the implementation of RSDM copes with the use of a large-scale system.The matrices of the original finite element model have dimensions 198015 × 198015 and they are reduced to 70 × 70 (0.035 %).It is interesting to assess the degree of fidelity with which RSDM is able to generate the ROM at a desired test point.
Only one parameter was considered for this structure, namely the radius of the legs.The initial value of the parameter was increased by 100 % and only two finite element analyses were run.To assess the accuracy with which the ROM interpolated by RSDM was computed, a full model was compared with the computed ROM at a chosen test point lying within the design space identified by the radius of the legs.
The reduction process of reduction for obtaining the two initial ROMs employed compund SEREP-CMS MOR, i.e. first the SEREP method was used to preserve the lowest 20 frequencies and subsequently a component-mode synthesis method (Craig-Bampton) was applied to the reduced matrices obtained by SEREP.
The accuracy of the ROM was assessed by applying two criteria: computation of the Normalized Relative Frequency Difference or NRFD, and the MAC number.The error estimate therefore relies on the following: NRFD = 1 − exact frequency(i) approximatingfrequency(i) • 100 %, where Φ and Ψ are the modal matrix pertaining to the direct ROM (DR) and the interpolated ROM (pR) obtained at the test point (TP), respectively.Figure 5 shows the NRFD using a SEREP-CMS MOR, and how successful the reduction process was successful in capturing the dynamics of the original system despite the considerable reduction.

Performance of RSDM for frequency prediction of a Hexapod
Despite the high accuracy of the MOR throughout the range of the 20 frequencies, Figure 6 shows how the dramatic reduction along with a 100 % change of the parameter has influenced the preservation of the whole range of frequencies at which the initial SEREP method was calibrated.The picture highlights the following features.First, the error committed by the ROM for calculating the fundamental frequency is about 0.0006 %, and it is less than 0.7 % for the first three frequencies.Second, the maximum error committed for the first 5 frequencies is below 1.5 %.Third, the maximum error is observed for the mode number 6, at which the error is slightly smaller than 3.9 %.The initial 5 frequencies are the 3 Number of frequencies of interest frequencies of main interest for an engineer.Consequently, RSDM is able to preserve the most important range.For the present study-case, NRFD is also supported also by the calculation of the Modal Assurance Criterion (MAC), the results for which are plotted in Figure 6(b).
For the first 6 frequencies, the values of MAC along the main diagonal range between a minimum of 98.7 % and a maximum of 99.9 %.These results provide evidence that the ROM generated at the test point by RSDM is accurate for the first 6 frequencies, even though a consistent reduction was performed by SEREP-CMS (0.035 % of the full model).

Conclusions
The main challenge that had to be faced in improving the accuracy of parameterized large-scale systems was to optimize the Direct Matrix interpolation method, i.e. to minimize the number of detailed finite element analyses and to guarantee high accuracy, while introducing as few as possible mathematical concepts.In an attempt to achieve this outcome, RSDM has been proposed as promising interpolation method for obtaining accurate parametric ROMs.The application of this method relies on a carefully manipulated SVD of the reduced stiffness and mass matrices, and on the use of only two finite element analyses.However, there are several reasons why this optimization may be compromised.First, if the number of the parameters and the design space are large, only two finite element analyses might not be sufficient.Second, if the size of the ROM is reduced to less than 0.1 % of the size of the original discretized structure, then the model might not be able to capture consistently capture the eigenbehavior of the system for the whole range of frequencies of interest.Nevertheless, this increase in ideality is manageable, and our paper has shown how to address the related difficulties by introducing of a new method of interpolation, named RSDM, with the aim of restoring and making available to the interpolation tool certain natural components belonging to the matrices of the full model that are related, on the one hand, to the process of reduction and, on the other hand, to the characteristics of a solid in the FE theory.This is a point that has been neglected in previously published works, yet it seems to be an indispensable step to adopt.Cautious use of the SVD method for achieving this purpose is effective for 2 reasons: first, it is numerically stable [12] and its algorithm is embedded in different computational software programs; second, the use of reduced matrices limits the computation time that it needs.When the proposed method is applied according to algorithm 2 the following performances have been observed: (1.)For the large-scale Hexapod structure the error committed lies between 0.0006 % and less than 1.5 % for the first 5 frequencies despite the consistent, initial reduction.
( .)There is very accurate correlation, as shown by the MAC number.
Engineering systems can be modeled by effective ROMs displaying various degrees of accuracy.On the basis of the results obtained here, it is evident that despite the increase in ideality attained by the RSDM, this degree of fidelity of ROM is in some way undermined by a combination of the dramatic reduction of the original system and the large number of parameters.One challenge that must be pointed out here is the need to further optimize the interpolation tool, in order to deal ideally with opposing requirements of computation time and accuracy.This of course entails first investigating how the error depends on the number of finite element analyses, on the interpolation method, and on the size of the final ROM.

Figure 1 .
Figure 1.Normalized Relative Frequency Difference (NRFD) vs Mode Number for a comparison between the proposed RSDM and the Direct Matrix Method (DMM).The proposed method preserves the first three frequencies better than DMM.

Figure 2 .
Figure 2.An example of a finite element mesh (a); a general deformation can be decomposed into a rotation and a stretch (b).The gradient of deformation is defined in continuum mechanics as, F = RS.The analogy with continuum mechanics is used to show the importance of recovering those components of stretch that, more than a rigid rotation, contribute to the change in the shape of the solid.

Figure 3 .
Figure 3. Correspondence between the process of reducing the stiffness matrix and the factors calculated by SVD, a Singular Value Decomposition Method.These factors have been identified after using Polar Decomposition to rule out the part associated with rigid rotation.

Figure 5 .
Figure 5. Error committed with the use of SEREP-CMS for calculating the first 20 frequencies of interest of a Hexapod.

Figure 6 .
Figure 6.Computation of the first 6 frequencies by the parametric ROM at the test point (a); Modal Assurance Criterion (MAC) [1].The computed ROM at the test point exhibits a good correlation, despite the relevant reduction of the full system (b) i = 1, 2, . . ., Nfoi 3 , while computing the MAC number we used MAC(i, j)