Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 56, 2, pp. 511-521, Warsaw 2018 DOI: 10.15632/jtam-pl.56.2.511 MODEL ORDER REDUCTION OF MIKOTA’S VIBRATION CHAIN INCLUDING DAMPING EFFECTS BY MEANS OF PROPER ORTHOGONAL DECOMPOSITION Wolfgang E. Weber Chair of Structural Analysis, Helmut-Schmidt-University/University of the Federal Armed Forces Hamburg, Hamburg, Germany; e-mail: wolfgang.weber@hsu-hh.de BERND W. ZASTRAU Institute of Mechanics and Shell Structures, Technische Universität Dresden, Dresden, Germany e-mail: bernd.zastrau@tu-dresden.de In engineering disciplines, both in scientific and practical applications, systems with a tre- mendous number of degrees of freedom occur. Hence, there is a need for reducing the com- putational effort in investigating these systems. If the systembehaviour has to be calculated for many time instances and/or load scenarios, the need for efficient calculations further increases. Model order reduction is a common procedure in order to cope with such large systems.The aimofmodel order reduction is to reduce the (computational) effort in solving the given task while still keeping main features of the respective system. One approach of model order reduction uses the proper orthogonal decomposition. This approach is applied toMikota’s vibration chain, a linear vibration chain with remarkable properties, where two cases of an undamped and a damped structure are investigated. Keywords: Mikota’s vibration chain, POD, damping, multibody system,mode shape 1. Introduction Model order reduction (MOR) is a common method in engineering disciplines allowing for ef- ficient calculations of e.g. dynamic behaviour of structures. The basic concept is to reduce the order of the system, for example by reducing the number of degrees of freedom. However, the resulting loss of information should not exceed a certain level. Several techniques forMORexist. Someof themare physical subspacemethods such asGuyan reduction,modal subspacemethods or Krylov subspace method (Freund, 2003; Guyan, 1965). While these methods are physically motivated, there are other approaches which do not take physical meaning into consideration. However, these methods still require extraction of main features of the underlying system. One method for this feature extraction is the proper orthogonal decomposition (POD). Once the arbitrary system has been characterised by POD, it is the task of MOR to only take those cha- racteristics into consideration which are needed to adequately describe the (dynamic) behaviour of the system. In this contribution, the model order reduction is applied to a special vibration chain, namely Mikota’s vibration chain. Herein, POD is used. Using Mikota’s vibration chain has the advantage that its dynamic characteristics have been investigated quite well, such that there exist analytical solutions which serve as reference solutions for the reduced system. Addi- tionally, damping effects are taken into account as some damage phenomena can be modelled bymeans of damping. To copewith these tasks,Mikota’s vibration chain is introduced in Section 2.Abrief descrip- tion of the approach of model order reduction involving POD is given in Section 3. In Section 4, POD is applied toMikota’s vibration chain. As a next step, a discrete damping element is added 512 W.E.Weber, B.W. Zastrau to Mikota’s vibration chain. In doing so, the linearity of the system is kept. POD is then used withinMORto approximate the dynamic behaviour of the damped vibration chain in Section 5. Finally, conclusions and an outlook are given with Section 6. 2. Mikota’s vibration chain Fig. 1. Undamped vibration chain (Weber et al., 2015) A linear undamped vibration system with n degrees of freedom (DOF) can be described by Mẍ(t)+Kx(t)=0 (2.1) with x= (x1, . . . ,xi, . . . ,xn) T and ẍ= (ẍ1, . . . , ẍi, . . . , ẍn) T representing the column matrix of displacements and accelerations, respectively. Herein,M andK denote the (diagonal) mass and (tri-diagonal) stiffnessmatrix, respectively. For a linear vibration chain according to Fig. 1 these matrices are M= diag(mi) K= diag {(k1+k2,−k2), . . . ,(−ki,ki+ki+1,−ki+1), . . . ,(−kn,kn)} (2.2) Mikota set the masses and stiffness coefficients to mi = 1 i m and ki =(n− i+1)k where i=1,2, . . . ,n, i∈ N (2.3) cf. (Mikota, 2001). Herein, m is the first mass and k denotes the stiffness of the last spring. Mikota conjectured that this specific vibration chain has eigenfrequencies Ωl = lΩ= l √ k m with again l=1,2, . . . ,n, l∈ N (2.4) where Ω = √ k/m is the first eigenfrequency. As can readily be seen, enlarging the system from n DOF to n+1 DOF changes the mass matrix in such a way that the element mn+1 is appended at the lower right corner leading to Mn,n →Mn+1,n+1, while all other entries of the mass matrix remain the same. In contrast, the corresponding matrix Kn+1,n+1 is obtained by adding one rowandone columnat theupper left corner to the formerKn,n. For discussionof this oppositebehaviour in the set-upof thematrices and the resultingdifficulties for provingMikota’s conjecture to be right, one is referred to e.g., (Müller andHou, 2007;Müller andGürgöze, 2006). However, two proofs were proposed byWeber et al. (2015), Müller and Hou (2007). In order to fully describeMikota’s vibration chain also themode shapes have to be looked at, whichwas not the focus of Mikota’s work. An approach based on polynomial coefficients is given with (Müller and Hou, 2007), a modification of the well-known Laguerre polynomials allowing for evaluating the mode shapes of Mikota’s vibration chain is presented in (Weber et al., 2013). But these approaches are quite laborious and do not reveal a structure in order to obtain general formulae for the mode shapes ul. This gap is closed with (Weber et al., 2017), some results of the latter contribution will be used in what follows. Model order reduction of Mikota’s vibration chain including damping effects... 513 According to Kochendörffer (1963), the coordinates of the eigenvectors of a matrix – and thus themode shapes dealt with here – can be expressed by polynomials pl(i) in the coordinate i leading to ul =(ul,i=1,ul,i=2, . . . ,ul,i=n) T =(pl(i=1),pl(i=2), . . . ,pl(i=n)) T (2.5) For an arbitrary nDOF and l¬n, the first three mode shapes of Mikota’s vibration chain are represented by the following polynomials pl=1(i) = i pl=2(i) = i 2− 2n+1 3 i pl=3(i) = i 3− 3 5 (2n+1)i2+ 1 5 [3 2 n(n+1)+1 ] i (2.6) see (Weber et al., 2017). There, an approach for determining all n polynomials, i.e. all n mode shapes, in a successive manner is suggested, too. For this special vibration chain it could be proved that the polynomial degree of pl is l. Another neat property of the mode shapes of Mikota’s vibration chain is the tri-diagonality of thematrix productUTU, whereU denotes the modalmatrix. It should be noted that, in general, themode shapes to different eigenfrequencies are not perpendicular to each other with respect to the (standard) scalar product, as only ulMu T k { =0 for l 6= k 6=0 for l= k ulKu T k { =0 for l 6= k 6=0 for l= k (2.7) holds. Exemplary, a graphical representation of the eigenvalues λl = Ω 2 l and the first five mode shapes of Mikota’s vibration chain with n=10 DOF is given in Fig. 2. Fig. 2. Eigensolutions of Mikota’s vibration chain for n=10DOF, where only the mode shapes ul=1, . . . ,ul=5 are shown. The displacements between the coordinates i are interpolated linearly, i=0 is at the fixed support according to Fig. 1 3. A brief introduction to MOR by means of POD The solution of typical large systems of equations, which occur in engineering disciplines, requ- ires a huge computational effort. Thus, strategies are sought which allow for reduction of this 514 W.E.Weber, B.W. Zastrau computational cost. One strategy is the model order reduction (MOR). In this strategy, the dimensionality of the underlying mechanical system is reduced while keeping the loss of infor- mation within acceptable bounds. For mechanical systems following the system of equations Mẍ+Dẋ+Kx= f (3.1) withDdenoting the dampingmatrix, ẋ the vector of velocities and thevector of applied forces f, this order reduction can be performed by projecting the involved vectors from the full space Rn to a lower dimensional space Rnred using a projection matrixΦ of dimension (n,nred) x≈Φxred ẋ≈Φẋred ẍ≈Φẍred (3.2) Inserting these quantities into Eq. (3.1) leads to MΦẍred+DΦẋred+KΦxred = f (3.3) Multiplying the transpose of the projectionmatrix from the left finally yields a reduced problem Φ T MΦẍred+Φ T DΦẋred+Φ T KΦxred =Φ T f Mredẍred+Dredẋred+Kredxred = fred (3.4) The question arises how this projection matrix Φ can be obtained. For reducing the dimen- sionality, the main features of the system have to be extracted. Within this contribution, the proper orthogonal decomposition (POD) is used to extract the main features and thus shall be introduced in what follows. In the first step, a suitable amount m of observations is necessary. These observations may result frommeasurements or from analytical or numerical calculations. In the present case, ana- lytical expressions of the displacements will be evaluated numerically and used as observations. The displacement history obtained by the so-called pre-computations is saved in an observation matrix Q= [x1,x2, . . . ,xm] (3.5) In general, the numberm of observations differs from the number n of DOF and, consequently, Q is a rectangularmatrix. This observationmatrix is decomposedbymeans of the singular value decomposition Q=PΣVT (3.6) according to Golub and Kahan (1965). Herein, P denotes the matrix of the left-singular vec- tors φk, which will be called proper orthogonal modes (POMs) in what follows. For special cases, these POMs are equivalent to themode shapes of the respective system.This issuewill be addressed in Section 4. ThematrixΣ is a pseudo-diagonalmatrix containing singular values σk, with k¬min(n,m), in a descending order at its main diagonal, whereas all other entries of the rectangular matrix are zero. ThematrixV of the right-singular vectors will not be used within this contribution. An energymeasure for the matrix Q consisting of a vector sequence is the Frobenius norm, which itself equals the sum of the squared singular values Epseudo(Q)= ‖Q‖ 2 F = 1 m n ∑ i=1 m ∑ j=1 Q2i,j = min(n,m) ∑ k=1 σ2k (3.7) see e.g., (Kerschen andGolinval, 2002). Thus, the value of σk is related to the so-called pseudo- -energy associated to the k-th POM. As a consequence, one may choose the first nred POMs, Model order reduction of Mikota’s vibration chain including damping effects... 515 with nred ≪m, which capture a certain amount of the system total pseudo-energy. From these POMs, the sought projection matrix Φ= [ φ1,φ2, . . . ,φnred ] nred ≪m (3.8) is constructed. With this projection matrix, the reduced system according to Eq. (3.4)2) is solved and the results are transferred back to the full system by using Eqs. (3.2). For a detailed review of the PODand some applications the reader is referred to e.g. (Bamer et al., 2017; Fangye et al., 2016; Radermacher and Reese, 2013; Kerschen et al., 2005; Chatterjee, 2000). 4. Applying POD to Mikota’s vibration chain As the first step,Mikota’s vibration chain is exposed to such an initial displacement which only excites its first mode shape. Afterwards, Mikota’s vibration chain performs free vibration. By doing so, the POD should identify the first mode shape only. This is due to the fact that the whole (vibrational) energy of the system is kept in this mode shape. Solving the differential equation with the following parameters m= k=1 and thus Ω= √ k m =1 (4.1) n=10 t0 =0 ∆t=0.01 tN =20 N =2001 (4.2) where the units have been omitted, yields the displacements xi for each time step tj ofN. These displacements are written into the observation matrix Q, which is then analysed by means of the POD according to Section 3. All calculations have been performed withMatlab R○. As can be seen in Fig. 3, there is indeed only one (dominant) singular value. The remaining singular values are not of practical relevance, as the quotient σk>1/σ1 < 10 −15 is in the order of numerical accuracy. Consequently, only the first POM is plotted whereas the remaining POMs are not taken into account. As expected, the POM 1 equals the first mode shape as given by Eq. (2.6)1 and Fig. 2. Fig. 3. Singular values and POMofMikota’s vibration chain, n=10DOF and initial excitation of the 1st mode shape, according to the magnitude of σk>1/σ1 < 10 −15 only the first POM is plotted 516 W.E.Weber, B.W. Zastrau Inwhat follows,Mikota’s vibration chain is exposed toa loadwith avery short timeduration, which simulates an impulse or an impact load. In detail, the force f1(t)= { 1 for 0¬ t¬ 0.05= δf 0 for δf =0.05 0 and startingwith the 1st DOF, energy is successively transferred to the remainingDOF.HigherDOFsare characterised byweaker springs and lower masses, see also Eqs. (2.3). The respective displacements increase with the increasing DOF, too. At time instances equal to integer multiples of π, all DOFs have zero displacements. At time instances equal to integermultiples of 2π, all DOFs but the 1st have zero velocities, too. With Eq. (4.4), the impulse loadmay be regarded as a Dirac-type loading. According toMüller and Schiehlen (1985), such an applied force in fact leads to non-vanishing initial velocities while Model order reduction of Mikota’s vibration chain including damping effects... 517 maintaining effectively zero initial displacements. In the present case, the 1stDOFgets an initial velocity (more precisely, the initial velocity is at the time instance δf). At this time instance, the whole system effectively still is at rest, see Fig. 4a. For this reason, the system starts a free oscillation from a state of zero displacement and, consequently, has to return into this state periodically.This happensat the aforementioned time instances at integermultiples of 2π, where the time duration δf has been omitted for clarity. At these time instances, all DOFs but the 1st must have both zero displacements and velocities. The velocity of the 1st DOF equals its initial velocity which is due to the Dirac-type loading, see also Fig. 4b. Due to the special structure of the eigenfrequencies, which is given by Eq. (2.4), the temporal factor at integer multiples of 2π of all n mode shapes is identical and equal to 1. Although Fig. 4a may lead to the assumption that there exist pronounced time intervals within which one or more masses mi are at rest, it should be emphasised that this is not the case. On the contrary, there are only distinct time instances – at integer multiples of π as discussed above – at which all massesmj have vanishing displacements. All masses but the first are at rest only at integer multiples of 2π. The displacements for all time instances according to Eq. (4.2), which are plotted in Fig. 4a, arewritten in theobservationmatrixQ, andthe latter is investigatedbymeansofPODaccording to Section 3. The resulting singular values and some POMs are given in Fig. 5. Fig. 5. Singular values and some POMs ofMikota’s vibration chain, n=10DOF and impulse load according to Eq. (4.3) applied to the 1st DOF, only POMs 1, . . . ,5 are shown Compared to Fig. 3, it can readily be seen that in the present load case the so-called pseudo- -energy is distributed over all POMs.However, the respective values, i.e. the singular values, are not equal. Hence, the single POMs each have a different contribution to the oscillation pattern of themechanical system. The descent of the singular values gives an indication of which POMs may be omitted within themodel order reductionwhile not exceeding the given error tolerance. This aspect will be investigated in the next Section, where additionally damping is taken into account. 5. Applying MOR to modified Mikota’s vibration chain including damping In this Section, a single absolute damping element is added to Mikota’s vibration chain. By doing so, both the eigenfrequencies and mode shapes of the system are changed. For a general introduction to the topic of (optimal) damping, the reader is referred to e.g. (Gürgöze and 518 W.E.Weber, B.W. Zastrau Müller, 1992). For a brief and exemplary discussion concerning the optimal position of absolute and relative damping elements inMikota’s vibration chain, seeWeber et al. (2008). In order to apply a useful model order reduction, the number of DOFs of the system is increased to n = 300. The absolute damping element is fixed at the 7th DOF. Only one non- -vanishing initial excitation x3 =1 is prescribed at the 3rd DOF. The other parameters are m= k= d=1 t0 =0 ∆t=0.01 tN =20 N =2001 (5.1) where again the units have been omitted. It should be noted that d = 1 does not lead to weak damping anymore. However, within this contribution, a parameter study is presented. The resulting displacement history for some DOFs is shown in Fig. 6. Fig. 6. Displacements for dampedMikota’s vibration chain with n=300DOF, initial excitation of the 3rd DOF, absolute damping element at the 7th DOF, only DOF 1,3,7,100,250 are shown In what follows, the model order reduction using POD is performed. As the initial step, the observation matrix Q is set by considering the displacements of all n = 300 DOFs within the time span t= 0 to t= 10. The POD is then applied to this observation matrix leading to the singular values and POMs as given in Fig. 7. Fig. 7. Singular values and POMs of dampedMikota’s vibration chain, n=300DOF, initial excitation at the 3rd DOF, absolute damping element at the 7th DOF, only POMs 1,10,50,70,100 are shown Model order reduction of Mikota’s vibration chain including damping effects... 519 The so-called pseudo energyEpseudo of the system can be calculated using Eq. (3.7) Epseudo = min(n,m) ∑ k=1 σ2k =977.5528 2 (5.2) A common approach to themodel order reduction is to consider such an amount nred of POMs that for the corresponding pseudo energy Epseudo,red ≈ 0.99Epseudo (5.3) holds (Bamer et al., 2017, Feeny andKappagantu, 1998).Kerschen et al. (2005) even recommend Epseudo,red ≈ 0.9999Epseudo. For the present case, the former criterion is fulfilled for nred = 76 while the latter criterion gives nred =98. Thus, in the first step, only 76 POMs are considered. As in the present case σ77/σmax = 0.0184, this means that all POMs, for which σk / 0.018σmax holds, are omitted. Although the observation matrix Q only contains data up to t = 10, the calculations in the reduced system have been performed until t=20. The resulting reduced system is solved and then transferred back to the full system. Both the displacement and velocity history within the time interval 3.25π ¬ t ¬ 3.75π for the arbitrarily chosen 150th DOF is given in Fig. 8. Additionally, the diagramcontains thedisplacementandvelocity history resulting fromcalculating the full system. As can be seen, there is a good agreement between these two results for both the displacement and velocity history. Fig. 8. Displacement and velocity history of the 150thDOF for the full (“full”, n=300) and the reduced (“MOR”, nred =76) system For comparison, an additional calculation is performed with nred = 98, thus neglecting all POMs for which σk / 2.4 ·10 −4σmax. The respective results can be taken from Fig. 9 and do not show any (observable) differences between the results obtainedwith the full and the reduced system. Besides Eq. (5.3), an additional relation is introduced tomeasure the deviation between the results obtained with the full system and the results obtained with the reduced system ∆Ephase = √ √ √ √ √ √ √ n ∑ i=1 [(xi−xred,i) 2+(ẋi− ẋred,i) 2] n ∑ i=1 (x2i + ẋ 2 i) (5.4) This relation givesmore reliable results as compared to the relationwhich only takes the pseudo energy into account. This is due to the fact that the latter relation does not reveal pronounced 520 W.E.Weber, B.W. Zastrau Fig. 9. Displacement and velocity history of the 150thDOF for the full (“full”, n=300) and the reduced (“MOR”, nred =98) system phasedifferenceswithin thedisplacements andvelocities, see also (Kappagantu andFeeny, 1999). In the present case, ∆Ephase =0.144 (compared to 0.5149 for nred =76) thus indicating a very good correlation between the results of the full and the reduced (nred =98) system. Hence, the model order reduction has been applied successfully. 6. Conclusions An approach to the model order reduction has been successively applied to Mikota’s vibration chain, a special vibration chain having remarkable properties.The chosen approach to themodel order reduction involves the proper orthogonal decomposition which, therefore, has been shortly introduced. Some basic insights into the proper orthogonal decomposition were given using the standard (that is, undamped) Mikota’s vibration chain with n = 10 DOFs. Finally, more advanced calculations including the model order reduction were performed after introducing a single absolute damping element into the vibration chain with n=300 DOFs. It was observed that the dynamic characteristics of Mikota’s vibration chain could be kept if the underlying mechanical system was reduced in such a way, that > 99% of the so-called pseudo energy was considered. In the present case, an excellent correlation between the results obtained with the full and the reduced system is obtained for nred =98 DOF or – equivalently – for reduction in the dimensionality by≈ 70%. In this contribution, a linear vibration system has been investigated. However, non-linear systems play an important role in engineering applications, too. Thus, more scientific work has to be done in the field of model order reduction of non-linear systems. This non-linearity may additionally be caused by damage of the respective structure. Such effects must be taken into account if the model order reduction is used in, e.g., structural health monitoring. References 1. Bamer F., Amiri A.K., Bucher C., 2017, A new model order reduction strategy adapted to nonlinear problems in earthquake engineering, Earthquake Engineering and Structural Dynamics, 46, 4, 537-559 2. Chatterjee A., 2000,An introduction to the proper orthogonal decomposition,Current Science, 78, 7, 808-818 Model order reduction of Mikota’s vibration chain including damping effects... 521 3. FangyeY.F.,WeberW.E., ZastrauB.W., BalzaniD., 2016, Some basic ideas for the simu- lation of wave propagation inmicrostructures using proper orthogonal decomposition,Proceedings in Applied Mathematics and Mechanics, 16, 1, 333-334 4. Feeny B.F., Kappagantu R., 1998, On the physical interpretation of proper orthogonalmodes in vibrations, Journal of Sound and Vibration, 211, 4, 607-616 5. Freund R.W., 2003, Model reduction methods based on Krylov subspaces, Acta Numerica, 12, 2, 267-319 6. Golub G., Kahan W., 1965, Calculating the singular values and pseudo-inverse of a matrix, Journal of the Society for Industrial and Applied Mathematics, Series B: Numerical Analysis, 2, 2, 205-224 7. Gürgöze M., Müller P.C., 1992, Optimal positioning of dampers in multi-body systems, Jo- urnal of Sound and Vibration, 158, 3, 517-530 8. Guyan R.J., 1965, Reduction of stiffness andmass matrices,AIAA Journal, 3, 2, 380-380 9. Kappagantu R., Feeny B.F., 1999, An “optimal” modal reduction of a system with frictional excitation, Journal of Sound and Vibration, 224, 5, 863-877 10. Kerschen G., Golinval J.-C., 2002, Physical interpretation of the proper orthogonal modes using the singular value decomposition, Journal of Sound and Vibration, 249, 5, 849-865 11. KerschenG.,Golinval J.-C.,VakakisA.F.,BergmanL.A., 2005,Themethodof proper or- thogonal decomposition for dynamical characterization and order reduction ofmechanical systems: an overview,Nonlinear Dynamics, 41, 1-3, 147-169 12. Kochendörffer R.K., 1963, Determinanten und Matrizen, B.G. Teubner Verlagsgesellschaft, Leipzig 13. Mikota J., 2001, Frequency tuning of chain structure multibody oscillators to place the natural frequencies atΩ1 andN−1 integermultiplesΩ2, . . . ,ΩN, Zeitschrift für AngewandteMathematik und Mechanik, 81, S2, 201-202 14. Müller P.C., Gürgöze M., 2006, Natural frequencies of a multi-degree-of-freedom vibration system,Proceedings in Applied Mathematics and Mechanics, 6, 1, 319-320 15. MüllerP,C.,HouM., 2007,Onnatural frequencies and eigenmodes of a linear vibration system, Zeitschrift für Angewandte Mathematik und Mechanik, 87, 5, 348-351 16. Müller P.C., Schiehlen W.O., 1985,Linear Vibrations, M. Nijhoff Publishers, Dordrecht 17. Radermacher A., Reese S., 2013, A comparison of projection-basedmodel reduction concepts in the context of nonlinear biomechanics,Archive of Applied Mechanics, 83, 8, 1193-1213 18. Weber W., Anders B., Müller P.C., 2015, A proof on eigenfrequencies of a special linear vibration system, Zeitschrift für Angewandte Mathematik und Mechanik, 95, 5, 519-526 19. Weber W., Anders B., Zastrau B.W., 2008, Some damping characteristics of a chain struc- tured vibration system,Proceedings in Applied Mathematics and Mechanics, 8, 1, 10391-10392 20. Weber W., Anders B., Zastrau B.W., 2013, Calculating the right-eigenvectors of a special vibration chain by means of modified Laguerre polynomials, Journal of Theoretical and Applied Mechanics, Sofia, 43, 4, 17-28 21. Weber W.E., Müller P.C., Anders B., 2017, The remarkable structure of the mode shapes and eigenforces of a special multibody oscillator, Archive of Applied Mechanics, to appear, DOI: 10.1007/s00419-017-1327-9 Manuscript received February 22, 2018; accepted for print March 1, 2018