Data-driven model order reduction for structures with piecewise linear nonlinearity using dynamic mode decomposition

Data-dr iv en model order reduction f or s tr uctures with piece wise linear nonlinear ity using dynamic mode decomposition ∗ Akira Saito † Masato T anaka ‡ Abstract Piece wise-linear nonlinear sy stems appear in man y engineering disciplines. Prediction of the dynamic beha vior of such sys tems is of great impor tance from practical and theoretical viewpoint. In this paper, a data-driven model order reduction method f or piecewise-linear sy stems is proposed, which is based on dynamic mode decomposition (DMD). The o v erview of the concept of DMD is provided, and its application to model order reduction for nonlinear sys tems based on Galerkin projection is explained. The proposed approach uses impulse responses of the sys tem to obtain snapshots of the state v ar iables. The snapshots are then used to extract the dynamic modes that are used to f or m the projection basis vectors. The dynamics descr ibed b y the equations of motion of the or iginal full-order sy stem are then projected onto the subspace spanned by the basis vectors. This produces a system with much smaller number of degrees of freedom (DOFs). The proposed method is applied to tw o representative ex amples of piece wise linear systems: a cantile v ered beam subjected to an elastic stop at its end, and a bonded plates assembl y with par tial debonding. The reduced order models (ROMs) of these sys tems are constructed b y using the Galerkin projection of the equation of motion with DMD modes alone, or DMD modes with a set of classical constraint modes to be able to handle the contact nonlinearity efficiently . The obtained R OMs are used for the nonlinear forced response analy sis of the sy stems under har monic loading. It is sho wn that the R OMs constructed b y the proposed method produce accurate f orced response results. Ke yw ords: Dynamic mode decomposition, Reduced order model, Piece wise-linear sys tems 1 Introduction Piece wise-linear (PWL) systems appear in many str uctural dynamics sys tems that are subjected to dynamic loadings. In par ticular , damaged structures such as with breathing cracks [ 1 – 3 ] or delaminated composites [ 4 ] are kno wn to be PWL sys tems, because they contain repetitive opening and closing of the contacting boundar ies. Theref ore, prediction of dynamic beha vior of such nonlinear sy stems is impor tant. How ev er , the dynamics of PWL sys tems are kno wn to show strong nonlinearities because of sudden chang es in their stiffness, which results in r ich dynamic phenomena including per iod doubling bifurcation and chaos [ 5 – 7 ]. This simple y et strong nonlinear ity hinders the application of standard linear anal ysis tools such as modal analy sis. A general nonlinear model-order reduction (MOR) methods exis t, which can reduce the size of the degrees of freedom (DOFs) of the systems [ 8 , 9 ]. Ho w ev er , its applicability to non-smooth PWL sy stems has nev er been studied. Theref ore, there hav e been many attempts to de v elop analy sis methods to predict dynamics of PWL sys tems in a computationally efficient manner , such as the ones based on nonlinear nor mal modes [ 10 ] and generalized bilinear amplitude approximation [ 11 , 12 ]. How ev er , due to the computational cost that stems from the increasing comple xity in the underl ying sys tems with the PWL nonlinearity , application of such anal ysis methods is difficult. This is because of the computational cost that comes ∗ This manuscr ipt is the accepted v ersion of the article: A. Saito and M. T anaka, “Data-driven model order reduction for str uctures with piecewise linear nonlinear ity using dynamic mode decomposition, ” Nonlinear Dynamics, 111, pp. 20597–20616 (2023). https: //doi.org/10.1007/s11071- 023- 08958- x © 2023 The A uthor(s). Published by Spr inger N ature. † Department of Mechanical Engineer ing, Meiji U niversity , Ka wasaki, Kanaga wa 214-8571, Japan. E-mail: asaito@meiji.ac.jp, OR CID: 0000-0001-8102-3754 ‡ T oy ota Research Institute of North America T oy ota Motor North Amer ica, Inc. Ann Arbor, Michigan, 48105, US A, OR CID: 0000-0003-0787-0335 1 from the ev aluation of the PWL nonlinearity . Thus, to circumv ent this difficulty , the application of a data-driven Galerkin-projection based MOR to structural dynamics problems in v ol ving PWL nonlinear ity is considered in this paper . With the increasing demand f or the analy sis of larg e-scale nonlinear sy stems, constr ucting reduced order models (R OMs) is becoming impor tant in v arious engineer ing disciplines. That is, reduction of the number of DOFs in the sys tem equations is of great impor tance to produce accurate results with reasonable amount of time, without compromising the computational accuracy . Theref ore, many model-order reduction (MOR) methodologies hav e been dev eloped. This paper proposes an MOR methodology based on Galerkin projection of the sys tem equations onto a set of basis v ectors obtained from dynamic mode decomposition (DMD), which is a data-dr iv en MOR method. A no v el approach f or compensating the dynamics of modal coordinates by using classical cons traint modes is also proposed to handle the contact f orce in the reduced order space effectiv ely . This combination of data-dr iv en DMD modes and the constraint modes that are obtained based on the sy stem model makes the proposed method a gray -box modeling approach [ 13 ], instead of purely data-dr iv en black -bo x modeling. The DMD was first de v eloped b y Schmid [ 14 ] as a method to e xamine coherent structures in the fluid flo w field. DMD has then been applied to many applications inv olving nonlinear systems, such as acoustic mode identification in a three dimensional chamber [ 15 ], a swirling flow problem [ 16 ], combustion [ 17 ], pressure sensitive paint data [ 18 ], and man y others. One of the greatest characteristics of the DMD o ver the proper orthogonal decomposition (POD) [ 19 , 20 ], which is also a data-dr iv en dimensionality reduction method, is that the DMD is able to e xtract not onl y spatial mode shapes but also their temporal inf or mation, i.e., frequency and deca y-rate. If the sy stem is linear , by using the mode shapes and the cor responding temporal information of the sys tem f or a given initial condition, DMD can be used to capture the behavior of the system b y a data-driven spectral decomposition [ 21 ], which is equation-free. Ho w ev er , if the sy stem is nonlinear , the modes obtained by the DMD cannot be used as such an equation-free method, unless the DMD is applied to find a modal representation of the Koopman operator acting on carefull y chosen obser vables of state variables such that the go v er ning nonlinear equations becomes linear . Otherwise, the obtained DMD modes are linear appro ximations of the nonlinear manifold and the linear superposition of the DMD modes does not necessar il y give accurate prediction of the or iginal nonlinear dynamical systems. Instead of equation-free approach, there hav e been attempts to utilize the DMD modes to construct the basis v ectors f or the R OMs based on Galerkin projection . The Galerkin projection model reduction has been successfully applied to various models such as the ones obtained by discretizing the advection diffusion equation b y finite difference method [ 22 ]. It has also been applied to a nonlinear reaction-diffusion equation [ 23 ]. In this paper , the Galerkin projection with DMD modes is applied for str uctural dynamics problems. There are onl y a f ew attempts to date to apply the DMD to structural dynamics problems [ 24 ]. F or instance, Simha et al [ 25 ] hav e applied the DMD to estimate natural frequencies of a linear elastic structure. T o the know ledge of the authors, the Galerkin projection has nev er been applied with DMD modes f or structural dynamics problems in v olving nonlinearity and their applicability as the basis v ectors has nev er been discussed. The remainder of this paper is structured as f ollow s. In Section 2 , general projection based model order reduction procedure f or equations of motion with displacement dependent nonlinear ity is discussed and the theor y of DMD is briefly re vie w ed. The DMD is then combined with the MOR procedure of the sys tem. In Section 3 , the application of the proposed method to the f orced response problem of a cantilev ered beam with an elastic stop is presented. In Section 4 , the proposed method is applied to a bonded shell assembly with par tial debonding, and a h ybr id projection basis of DMD modes and constraint modes is ex amined. Finally in Section 5 , the concluding remarks are pro vided. 2 Method 2.1 Model Order Reduction Assuming that there is an elastic body that is subject to dynamic loading and PWL nonlinear ity , and that the resulting oscillation is infinitesimall y small. Also assume that the go verning par tial differential equations are discretized b y a numer ical method, such as finite element (FE) method, which results in a discrete dynamical system with piece wise-linear nonlinear ity of size 𝑚 . Denoting the nodal displacement v ector as u ( 𝑡 ) , the gov er ning equations of 2 the sys tem are then wr itten as: M ¥ u ( 𝑡 ) + C ¤ u ( 𝑡 ) + Ku ( 𝑡 ) + f ( u ) = b ( 𝑡 ) (1) where M , C , and K are mass, damping, and stiffness matr ices, f ( u ) is the nonlinear force caused by the piece wise- linear nonlinear ity , b ( 𝑡 ) is a time-dependent e xter nal f orcing v ector , and M , C , K ∈ R 𝑚 × 𝑚 , and f ( u ) , b ( 𝑡 ) ∈ R 𝑚 . In many cases, it is desirable if the sys tem size can be reduced by applying MOR on the sys tem Eq. ( 1 ), i.e., the Galerkin projection with a set of appropriately chosen basis v ectors. That is, a set of v ectors in R 𝑚 , which are linearl y independent, or preferabl y or thogonal set of v ectors of dimension 𝑚 is selected, i.e., 𝚽 = [ ϕ 1 , . . . , ϕ 𝑝 ] (2) where 𝑝 is the number of chosen basis vectors, ϕ 𝑖 ∈ R 𝑚 , 𝑖 = 1 , . . . , 𝑝 , and in many cases, 𝑝 ≪ 𝑚 . Then, a coordinate transf or mation cor responding to Eq. ( 2 ) is introduced, i.e., u ( 𝑡 ) ≈ 𝑝  𝑖 = 1 𝜂 𝑖 ( 𝑡 ) ϕ 𝑖 = 𝚽 η ( 𝑡 ) , (3) where 𝜂 𝑖 ( 𝑡 ) is the modal coordinate cor responding to ϕ 𝑖 that acts as the variable to be sol ved in the reduced order space, and η ( 𝑡 ) = [ 𝜂 1 ( 𝑡 ) , . . . , 𝜂 𝑝 ( 𝑡 ) ] T . Then, the go verning equations of reduced size can be obtained with respect to η ( 𝑡 ) , as f ollo ws: ˜ M ¥ η ( 𝑡 ) + ˜ C ¤ η ( 𝑡 ) + ˜ K η ( 𝑡 ) + ˜ f ( η ) = ˜ b ( 𝑡 ) , (4) where ˜ M = 𝚽 T M 𝚽 , ˜ C = 𝚽 T C 𝚽 , ˜ K = 𝚽 T K 𝚽 , ˜ b = 𝚽 T b , and ˜ M , ˜ K , ˜ C ∈ R 𝑝 × 𝑝 , ˜ b ∈ R 𝑝 . In general, the nonlinear f orcing vector ˜ f ( η ) needs to be ev aluated in full-order space with u especiall y f or the PWL sys tems of interest, i.e., ˜ f ( η ) = 𝚽 T f ( 𝚽 η ) , because the ev aluation of the piece wise linear ter m needs to be done in ph y sical coordinate sy stems. Theref ore, one needs to expand the modal coordinate v ector back to full-order coordinate sys tem at ev er y single time step if time integration scheme is applied. This is unav oidable if one uses Galerkin projection, unless special treatment is given to the choice of the modes, as it will be shown in Section 4.1 . It is impor tant to choose appropr iate set of ϕ 𝑖 ’ s for creating accurate R OMs especially f or large-scale nonlinear problems because the nonlinear problem of interest needs to be solv ed in the subspace spanned b y ϕ 𝑖 , f or 𝑖 . . . 𝑝 . Theref ore, the subspace should be sufficientl y larg e such that it encompasses the dynamics of interest, y et it is small enough to conduct nonlinear calculations efficiently with reasonable accuracy . In this paper , the basis v ectors der ived from DMD are considered and the impact of the choice of ϕ 𝑖 ’ s on the resulting dynamics of Eq. ( 1 ) is discussed. Other impor tant bases, such as linear normal modes (LNM), which are obtained as the eigen vectors of M and K and denoted as L = [ L 1 , . . . , L 𝑝 ] , L 𝑖 ∈ R 𝑚 , and eig en v ectors obtained from the POD, which will be discussed later , are also taken into account f or compar ison. 2.2 Dynamic Mode Decomposition In this section, mathematical background of the DMD is br iefly revie w ed. First, let us assume that x ( 𝑡 ) ∈ R 𝑚 denotes a time-dependent v ector of generalized coordinates, or a state v ector of an 𝑚 -dimensional dynamical sys tem. No w assume that a snapshots of the generalized coordinates are taken at discrete time instants 𝑡 𝑗 f or 𝑗 = 1 , . . . , 𝑛 , and are stored in a single matr ix so called snapshot matrix , X = [ x 1 , . . . , x 𝑛 ] f or x 𝑗 ≜ x ( 𝑡 𝑗 ) and 𝑡 𝑗 + 1 = 𝑡 𝑗 + Δ 𝑡 where the size of Δ 𝑡 is fix ed. Ne xt, in DMD formulation, another snapshot matr ix Y is considered, where the columns of Y are the time-shifted versions of X , i.e., Y ≜ [ x 2 , . . . , x 𝑛 + 1 ] . The ke y assumption w e make here is that there e xists a linear map represented by a matr ix A between x 𝑗 and x 𝑗 + 1 , i.e., x 𝑗 + 1 = Ax 𝑗 . (5) This yields, Y = AX . (6) It means that the f ollowing relationship holds: A = YX † , (7) 3 and X † denotes the Moore-P enrose pseudo in verse matrix of X . The DMD modes and eigen values are defined as the eig en v ectors and eigen values of A [ 14 ]. By applying singular value decomposition (SVD) to X , X = U 𝚺 V ∗ where U ∈ C 𝑚 × 𝑟 , 𝚺 ∈ R 𝑟 × 𝑟 , and V ∈ C 𝑛 × 𝑟 , 𝑟 is the rank of X , and ∗ denotes the Hermitian transpose. N ote that U contains the eigen v ectors of the POD [ 20 ], or the POD modes , and the y are designated as [ U 1 , . . . , U 𝑟 ] . The diagonal terms in 𝚺 are denoted as 𝜎 𝑖 , 𝑖 = 1 , . . . , 𝑟 . The pseudo in v erse X † can then be computed b y using V , U , and 𝚺 , as follo w s: X † = V 𝚺 − 1 U ∗ . (8) The DMD eigen value is no w obtained b y computing the eig env alues of the projected matr ix A onto the POD modes [ 26 ]. Namel y , defining the projected matr ix as ˜ A ≜ U ∗ A U , = U ∗ YX † U , = U ∗ YV 𝚺 − 1 U ∗ U , = U ∗ YV 𝚺 − 1 . (9) No w the DMD eigen values of ˜ A can be obtained b y sol ving the f ollo wing eigen value problem, i.e., ˜ Aw 𝑖 = 𝜇 𝑖 w 𝑖 , (10) where 𝜇 𝑖 is the 𝑖 th DMD eigen value. The DMD mode φ 𝑖 cor responding to 𝜇 𝑖 of dimension 𝑚 is then defined b y , as in the definition of the exact DMD in R ef. [ 26 ]:, φ 𝑖 = ( YV 𝚺 − 1 ) w 𝑖 . (11) This holds because φ 𝑖 is the eigen vector of A corresponding to 𝜇 𝑖 , i.e., A φ 𝑖 = YV 𝚺 − 1 U ∗ YV 𝚺 − 1 w 𝑖 , = YV 𝚺 − 1 ˜ Aw 𝑖 , = YV 𝚺 − 1 𝜇 𝑖 w 𝑖 , = 𝜇 𝑖 φ 𝑖 . (12) It is noted that most DMD eigen values appear with their comple x-conjugate pairs, so do the cor responding DMD modes. Theref ore, both real and imaginar y par ts are used f or constructing the projection basis. It is also kno wn that the DMD yields an approximate eigen-decomposition of the best-fit linear operator relating tw o data matrices X and Y [ 27 ], which in this case is the matrix A . One of the advantag es of the DMD o ver POD is that the DMD eigen values can be associated with their frequency and decay rate, or damping ratios. Namel y , defining 𝑠 𝑖 , 𝑖 = 1 , . . . 𝑟 , as 𝑠 𝑖 = log ( 𝜇 𝑖 ) / Δ 𝑡 , (13) which is also wr itten as 𝑠 𝑖 = − 𝜁 𝑖 𝜔 𝑖 ± j 𝜔 𝑖  1 − 𝜁 2 𝑖 where j = √ − 1 with the assumption that all DMD modes decay and 0 ⩽ 𝜁 𝑖 ⩽ 1 . Then, the f ollowing relationships hold. 𝑓 𝑖 = 𝜔 𝑖 / 2 𝜋 = | 𝑠 𝑖 | / 2 𝜋 , (14) 𝜁 𝑖 = − Re ( 𝑠 𝑖 ) / | 𝑠 𝑖 | , (15) where 𝜔 𝑖 , 𝑓 𝑖 and 𝜁 𝑖 are the angular frequency, frequency and damping ratio cor responding to the 𝑖 -th DMD mode, respectiv ely . Although the y are er ror prone, it was sho wn that 𝑓 𝑖 and 𝜁 𝑖 coincide with the undamped natural frequency and modal damping ratio if the matrices X and Y come from the response of a linear sys tem [ 28 ]. If the sys tem is 4 nonlinear , the y are not necessarily the natural frequencies and modal damping ratios. F ur thermore, the magnitude of 𝜇 𝑖 can be used to rank the DMD modes based on its dominance in the snapshot matr ices. Namel y , 𝜇 𝑖 = e xp { R e ( 𝑠 𝑖 ) Δ 𝑡 ± jIm ( 𝑠 𝑖 ) Δ 𝑡 } , = e xp { − 𝜁 𝑖 𝜔 𝑖 Δ 𝑡 } e xp  ± j 𝜔 𝑖  1 − 𝜁 2 𝑖 Δ 𝑡  . (16) Theref ore, | 𝜇 𝑖 | = | e − 𝜁 𝑖 𝜔 𝑖 Δ 𝑡 | = 1 / | e 𝜁 𝑖 𝜔 𝑖 Δ 𝑡 | . (17) This means that if | 𝜇 𝑖 | is larg e, then 𝜁 𝑖 𝜔 𝑖 is small f or a fixed Δ 𝑡 , i.e., the decay rate of the oscillation of the cor responding DMD mode component is small. In other w ords, it deca ys slow ly . If | 𝜇 𝑖 | is small, then 𝜁 𝑖 𝜔 𝑖 is larg e, i.e., the decay rate of the oscillation of the cor responding component is larg e. It means that it deca y s fas t. Theref ore, when constructing a R OM using φ 𝑖 b y the Galerkin projection, it is desirable to include modes with low decay rate, or , with larg e | 𝜇 𝑖 | as man y as necessary , because the y are e xpected to dominate the dynamics. 2.3 Analy sis Procedure In this paper, f orced response of structures with PWL nonlinearity is considered. The process of the anal y sis based on DMD-based MOR is threef old: (1) computation of the snapshots, (2) R OM construction with the DMD modes, and (3) the nonlinear f orced response calculation with the R OM. Firs t, in order to obtain the DMD-based R OMs, snapshots [ x ( 𝑡 1 ) , x ( 𝑡 2 ) , . . . , x ( 𝑡 𝑛 + 1 ) ] need to be obtained numerically or experimentally . Suppose that it is obtained numer ically , the difficulty of the procedure of the computation of the snapshots is problem dependent. Also, if the snapshot is obtained by sol ving initial value problems, it is impor tant to set initial conditions or initial loading such that they ex cite sufficiently many DMD modes that are possibly in v ol v ed in the dynamics to be predicted with the R OM. Second, with the snapshot matr ices, DMD needs to be conducted. The challeng e here is that the DMD spectr um ma y contain DMD eig en values that cor respond to f ast-deca ying modes or non-e xisting modes that do not represent an y phy sical dynamic phenomena. Theref ore, such DMD modes need to be remov ed. In this paper, DMD mode selection is conducted based on singular value rejection [ 29 ], and successiv e application of DMD with increasing sampling frequency [ 28 ]. N ote that the sampling frequency of the or iginal time history cannot be chang ed once it is computed. Theref ore, it is re-sampled with larg er sampling per iods, and the DMD is conducted to see whether the DMD eigen values consistentl y appear regardless of the sampling frequency . This is inspired b y the concept of the stability diag ram [ 30 ] that is used in the frequency-domain pole es timator used in the area of e xper imental modal anal ysis. Moreov er , the computed eigen values are ranked based on their magnitudes to further promote the sparcity of the DMD eigen values, which results in the reduction of the sys tem size in the R OM. Third, the obtained DMD mode shapes φ 𝑖 ’ s are used to f or m a transf or mation matrix 𝚽 = [ φ 1 , . . . , φ 𝑝 ] . The reduced sys tem matrices and reduced f orcing v ectors are then obtained. There are possibl y two methods to handle nonlinear force that comes from the PWL nonlinear ity . The first approach is to compute the nonlinear force in full-order space and project it back to reduced space. Namel y , it is schematicall y wr itten as follo w s. f = f ( 𝚽 η ( 𝑡 ) ) → ˜ f = 𝚽 T f . (18) This process requires the e xpansion of the modal coordinates η ( 𝑡 ) back to larg er set of phy sical DOFs u ( 𝑡 ) at ev er y time step, because all modes in 𝚽 contribute to the nonlinear f orce. The second approach is to compute the nonlinear f orce in reduced-order space directly , without e xpanding the modal coordinates back to phy sical DOFs. In this paper , it is sho wn that this can be achiev ed by adding constr aint modes [ 31 ] in 𝚽 so that the DOFs subjected to the nonlinearity are accessible in the reduced-order subspace. The process is schematicall y wr itten as f ollo ws. ˜ f = ˜ f ( η ( 𝑡 ) ) . (19) This process is computationall y more efficient than the first approach represented as Eq. ( 18 ) , because the multiplications of 𝚽 T η and 𝚽 T f are not necessar y . The ov erall MOR algorithm e xplained abov e is shown in Algorithm 1 . 5 Generation of the f ull-order system (e.g., FEM) (D) Nonlinear f orced response calcula tions (a) Full-order ana lysis (A) Snapshot extraction (initial value problem) (B) DMD calculations (C) Galerkin projec tion using DMD modes and matrix-based modes (D) Nonlinear f orced response calcula tions (c) Hybrid data- driven DMD-ROM (b) Classical model-driven MOR (C) Galerkin projec tion using matrix-based modes (D) Nonlinear f orced response calcula tions (d) Data-dr iven DMD-ROM (A) Snapshot extraction (initial value problem) (B) DMD calculations (D) Nonlinear f orced response calcula tions (C) Galerkin projec tion using DMD modes Figure 1: Ov er vie w of the f orced response analy sis T o better visualize the structure of the proposed methodology , analy sis flo w char ts of typical f orced response calculations are sho wn in Fig. 1 . First, one needs to generate the full-order sys tem b y discretizing the spatial v ar iables of the go verning equations by numer ical methods, such as FEM. A naive approach for obtaining the f orced response solution is to r un the f orced response with full-order sys tem, as sho wn in Fig. 1 (a). One may w ant to appl y classical model-based MOR to the system that are based on super position of modes that can be obtained by modal and static analy ses of the sy stem matrices with special boundar y conditions, such as Craig-Bampton matr ix condensation [ 31 ], as shown in Fig. 1 (b). This process consists of the Galerkin projection of the equation onto the set of modes, and the nonlinear forced response calculations of the reduced sy stem. This has been a good approach f or gaining computational efficiency without degrading the accuracy . How ev er , the modes obtained purel y from static or modal analy ses are from linear analy ses. Theref ore, the y may not necessarily span the subspace where the nonlinear response lies. Data-driv en MOR, which was proposed in Ref. [ 22 ] and is e xtended to PWL sys tem in this paper , on the other hand, is e xpected to be able to use a good basis v ectors f or nonlinear f orced response calculations, because the y come from DMD of nonlinear response of an initial v alue problem of the system to be studied in the f orced response calculations. This process is sho wn in Fig. 1 , which consists of the e xtraction of snapshots and DMD calculations from time ser ies data. The Galerkin projection is then applied to the sytem equations to obtain the R OM. Nonlinear f orced response is then computed. This is considered in Section 3 . F ur thermore, when one f or ms R OMs for PWL sy stems, handling the contact f orce at the contacting boundar ies can be e xpensiv e ev en with R OMs either b y classical model-based MOR or data-dr iv en MOR, because the ev aluation of the contact f orce needs to be done in phy sical domain, which results in repetitive computation of modal coordinates back to ph y sical coordinates at ev er y time instant. Given the number of DOFs that are subject to PWL nonlinear ity is not prohibitiv ely large, one is able to efficiently compute the nonlinear force in the reduced subspace b y keeping the phy sical DOFs in the R OM, b y using static modes called constraint modes. This process is sho wn in Fig. 1 (c), which is denoted here as Hybr id data-driven MOR, and this approach is considered in Section 4 . There ha ve been a couple of methods to apply projection-based methods to nonlinear sy stems. For instance, Discrete Empir ical Inter polation Method (DEIM), which is based on POD, has been proposed to reduce the number of DOFs where nonlinear ter ms are ev aluated [ 32 ]. The method has been e xtended to be used with DMD as w ell [ 22 , 33 ]. Hyper -reduction has also been proposed f or enhancing the numer ical efficiency of the R OM b y using a set of sampling gr id points in the computational domain to specify reduced integration domains where the problem is sol v ed [ 34 ]. The problems studied in this paper are assumed to be fundamentall y localized , i.e., the size of the 6 Algorithm 1 Model order reduction based on Dynamic Mode Decomposition 1: Obtain [ x 1 , x 2 , . . . , x 𝑛 + 1 ] f or 𝑡 𝑗 + 1 = 𝑡 𝑗 + Δ 𝑡 by solving initial v alue problem for a given initial condition 2: f or 𝑘 = 1 , . . . 𝑘 𝑚 𝑎 𝑥 do 3: Δ 𝑡 𝑘 ← ( 𝑘 𝑚 𝑎 𝑥 / 𝑘 ) Δ 𝑡 4: R e-sample the snapshots with Δ 𝑡 𝑘 and obtain [ x 1 , x 2 , . . . , x 𝑛 + 1 ] f or 𝑡 𝑗 + 1 = 𝑡 𝑗 + Δ 𝑡 𝑘 5: X ← [ x 1 , x 2 , . . . , x 𝑛 ] and Y ← [ x 2 , x 3 , . . . , x 𝑛 + 1 ] 6: Compute SVD: U 𝚺 V ∗ ← X 7: ˜ A ← U ∗ YV 𝚺 − 1 8: Sol v e ˜ Aw 𝑖 = 𝜇 𝑖 w 𝑖 9: 𝑠 𝑖 ← log ( 𝜇 𝑖 ) / Δ 𝑡 𝑘 , 𝑓 𝑖 ← | 𝑠 𝑖 | / 2 𝜋 10: φ 𝑖 =  YV 𝚺 − 1  w 𝑖 11: end f or 12: Form 𝚽 = [ φ 1 , . . . , φ 𝑝 ] 13: Compute ˜ M = 𝚽 T M 𝚽 , ˜ K = 𝚽 T K 𝚽 , ˜ C = 𝚽 T C 𝚽 , ˜ b = 𝚽 T b ( 𝑡 ) Figure 2: Cantile vered beam subject to har monic f orcing and elastic stop region where the PWL nonlinear ity appears is relativ ely small in compar ison with the entire computational domain, y et it is subjected to the strong PWL nonlinear ity . Hence, these methods that are specifically tailored f or reducing the number of DOFs that are subject to nonlinearities do not hav e to be applied to the systems of interest in this study . The f ollowing sections provide numer ical e xamples of the application of the proposed approach to the f orced response problems of mechanical sys tems that are subject to piecewise-linear nonlinear ities. 3 F orced response of a beam with an elastic stop with data-driv en model-order reduction The approach is applied to a cantilev ered beam with an elastic stop at its end, and is subject to an external f orcing. The schematics of the problem is shown in F ig. 2 . The go verning equation of motion of the beam is obtained with the Euler -Ber noulli beam theor y , in strong f or m as follo ws, 𝜌 𝐴 𝜕 2 𝑤 𝜕 𝑡 2 + 𝐸 𝐼 𝜕 4 𝑤 𝜕 𝑥 4 = 0 , (20) where 𝐸 is the Y oung’ s modulus, 𝐼 is the second moment of area, 𝜌 is the density , 𝐴 is the cross-sectional area of the beam. The boundar y conditions at 𝑥 = 0 and 𝑥 = ℓ are, 𝑤 ( 0 , 𝑡 ) = 0 , (21) 𝐸 𝐼 𝜕 3 𝑤 𝜕 𝑥 3 ( ℓ , 𝑡 ) = 𝑏 ( 𝑡 ) + 𝑓 𝑛𝑙 ( 𝑤 ( ℓ , 𝑡 ) ) , (22) 7 0 0.5 1 1.5 2 2.5 Time (s) 10 -4 0 20 40 60 80 100 Force (N) (a) Impulsiv e force applied at the tip 0 0.05 0.1 0.15 Time ( s ) -8 -6 -4 -2 0 2 4 6 8 Displacement (m) 10 -3 Linear (without elastic stop) Nonlinear (with elastic stop) (b) Response of the beam at the tip Figure 3: Applied force and the response of the beam f or both linear and nonlinear cases where 𝑏 ( 𝑡 ) is the external har monic f orcing, and 𝑓 𝑛𝑙 ( 𝑤 ( ℓ , 𝑡 ) ) is the PWL nonlinear restoring f orce acted upon by the elastic stop. Namel y , the spring generates the restoring force only when 𝑤 ( ℓ , 𝑡 ) ⩽ 0 , i.e., 𝑓 𝑛𝑙 ( 𝑤 ( ℓ , 𝑡 ) ) = − 𝑘 𝑐 max ( − 𝑤 ( ℓ , 𝑡 ) , 0 ) , (23) where max ( 𝑥 , 𝑦 ) returns 𝑥 if 𝑥 > 𝑦 , 𝑦 if 𝑥 < 𝑦 , or 𝑥 ( = 𝑦 ) if 𝑥 = 𝑦 , and 𝑘 𝑐 is the spr ing constant of the elastic stop. With these gov er ning equations, the go verning equation w as discretized by Euler -Ber noulli beam elements, which resulted in 32 elements of 64 DOFs where each node has a transv erse and a rotational DOFs. The fixed boundar y condition of Eq. ( 21 ) has been achie v ed b y removing the DOFs at the root. The resulting equations of motion hav e the f or m: M ¥ w ( 𝑡 ) + C ¤ w ( 𝑡 ) + Kw ( 𝑡 ) = b ( 𝑡 ) + f ( w ) (24) where w is the displacement v ector containing all transverse and rotational DOFs of all nodes. 3.1 Snapshot extraction and dynamic mode decomposition As discussed in Section 2 , snapshots should be obtained to compute DMD or e v en POD, which is the characteristics of these data-dr iv en approach. The choice of snapshots is arbitrar y , but they need to contain r ich information about the dynamics to be captured b y the R OM using DMD modes. Theref ore, for this numer ical e xample, snapshots w ere obtained b y numer icall y solving Eq. ( 24 ) for an impulsiv e f orce, i.e., impulse response , which is e xpected to ex cite the motion of the beam for a wide frequency rang e. The impulse response of the beam can be obtained b y applying impulsive f orce at the tip of the beam and the time response of the system was computed b y a time integration method. The applied impulsiv e force was defined as a sinusoidal w a v e of half-per iod define as, 𝑏 ( 𝑡 ) =  𝑏 0 sin ( 2 𝜋 𝑡 / 𝑇 ) , 0 ⩽ 𝑡 ⩽ 𝑇 / 2 0 , 𝑡 > 𝑇 / 2 (25) where 𝑇 is the per iod of the sinusoid with 𝑇 = 1 . 0 × 10 − 4 s and 𝑏 0 = 100 N. The value of 𝑇 w as chosen such that it is smaller than the per iod of the f astest dynamics to be considered in this problem, which in this case was judiciously chosen to be 4000 Hz in frequency or 2 . 5 × 10 − 4 s in time duration. The time history of the applied f orce is sho wn in Fig. 3 (a). The equation of motion Eq. ( 24 ) with w ( 0 ) = 0 and ¤ w ( 0 ) = 0 w as solv ed b y ode45 sol v er implemented in Matlab ® . The value of the spring constant f or the elastic stop was chosen to be 𝑘 𝑐 = 1000 N/m. The obtained impulse response of the beam subject to the elastic stop measured at the tip is shown in Fig. 3 (b). T o see the effects of the elastic stop at the tip, the response without the elastic stop, or the case where 𝑘 𝑐 = 0 is labeled as Linear in 8 0 1000 2000 3000 4000 Frequency (Hz) 10 -6 10 -5 10 -4 10 -3 10 -2 Displacement (m) 2 4 6 8 Sampling Frequency (Hz) 10 4 Lowest Highest DMD eigenvalue rank (a) Linear (no elastic stop) 0 1000 2000 3000 4000 Frequency (Hz) 10 -6 10 -5 10 -4 10 -3 10 -2 Displacement (m) 2 4 6 8 Sampling Frequency (Hz) 10 4 Lowest Highest DMD eigenvalue rank (b) Nonlinear (with elastic stop) Figure 4: Pseudo-s tability diagram of the DMD eigen values, along with the FFT results of the displacments. — : FFT spectra of all DOFs, — : av erage spectr um, × : DMD eig en v alues. T able 1: DMD eigen values in Hz f or linear and nonlinear cases Mode number Linear Nonlinear De viation (%) 1 50.9 65.3 28.29 2 320 316 -1.250 3 893 886 -0.7839 4 1745 1735 -0.5731 5 2886 2874 -0.4158 Fig. 3 (b). As can be seen, both responses are identical until the tip str ik es the elastic stop where 𝑤 = 0 . After the first strike, the tip displacement chang es its direction due to the elastic stop, which cannot be seen in the linear response. This process is repeated until this bouncing beha vior stops when the kinetic energy disappears. It is noted that higher frequency components are e x cited in the nonlinear response than in the linear response. Ne xt, using the obtained snapshots of the impulse response, DMD has been applied. Based on the algor ithm sho wn in Algor ithm. 1 , DMD modes were extracted for both linear and nonlinear cases. The results are shown as the pseudo-stability diagram of the DMD [ 28 ]. In the pseudo-stability diag ram, DMD eigen values with increasing sampling frequencies and the spectra of the displacement computed by the F ast Fourier T ransf or m (FFT). The color of the DMD eigen value represents the ranking of the eig en value among the ones computed f or a specific sampling frequency . The larg er the magnitude of the DMD eigen value is, the higher its ranking is. From this diag ram, w e can see if a DMD eigen value and the cor responding DMD mode shape should be taken into account in the R OM, i.e., if the DMD eigen value appears with high ranking regardless of the sampling frequency , then the eigen value is stable and hence is impor tant. It means that it is likel y that it needs to be taken into account in the R OM. Based on the diag rams shown in Fig. 4 , we can see that there are fiv e dominant DMD eigen values for both linear and nonlinear cases. The DMD eigen values cor respond to the peak frequencies of the FFT spectra, especially f or the linear case. In fact, it is known that DMD is equiv alent to the Ibrahim’ s time domain modal parameter e xtraction method [ 28 , 35 ]. Theref ore, if the DMD is applied to the linear response, the obtained DMD modes match ex actly the linear normal modes obtained from the mass and stiffness matr ices, and the corresponding DMD eig env alues contain the undamped natural frequencies. Fur thermore, the modal damping ratio obtained from the DMD eigen values ex actly match the modal damping ratios obtained from the mass, stiffness, and damping matr ices 9 0 0.2 0.4 0.6 0.8 1 -1 -0.5 0 0.5 1 Normalized displacement Mode 1 Mode 2 Mode 3 Mode 4 Mode 5 (a) LNM 0 0.2 0.4 0.6 0.8 1 -1 -0.5 0 0.5 1 Normalized displacement Mode 1 Mode 2 Mode 3 Mode 4 Mode 5 (b) POD modes 0 0.2 0.4 0.6 0.8 1 -1 -0.5 0 0.5 1 Normalized displacement Mode 1 Mode 2 Mode 3 Mode 4 Mode 5 (c) DMD modes Figure 5: Comparison between LNM, POD modes and DMD modes of the system if the damping is modeled as propor tional damping where the damping term can also be diagonalized b y the linear nor mal modes. On the other hand, as can be seen in Fig. 4 (b) f or the nonlinear case, the FFT spectra contain many frequency components and clear dominant peaks cannot be deter mined from the FFT alone, especiall y f or the lo w frequency range. Ho we ver , the DMD eig en values coincide with the peak frequencies of the FFT spectra. The obtained frequencies cor responding to the DMD eigen values are shown in T able 1 . Note that the frequencies e xtracted from the linear response are identical to the natural frequencies that can be obtained from the mass and stiffness matr ices of the sys tem. On the contrar y , the DMD eigen values extracted from the nonlinear response do not match the natural frequencies as expected. The most significant shift in the eigen value is 28.29% f or the first mode in comparison with that of the linear case. This appears to be the stiffening effect due to the elastic stop at the tip. The de viations f or all the other modes are much smaller than that f or the first mode. 1 2 3 4 5 6 7 8 9 10 Number of modes 0 0.5 1 1.5 2 Angle (rad) Figure 6: Pr incipal angles between the subspaces spanned by LNMs, POD modes, and DMD modes f or different number of modes 10 0 50 100 150 Frequency (Hz) 10 -3 10 -2 Displacment amplitude (m) Full (64DOF) DMD(p=1) POD (p=1) DMD (p=5) POD (p=5) LNM(p=1) LNM(p=5) 80 81 82 83 84 Frequency (Hz) 7 8 9 Displacment amplitude (m) 10 -3 Figure 7: F orced response results of full-order FEM, POD, and DMD Ne xt, the obtained mode shapes are discussed. First fiv e modes are shown in Fig. 5 f or LNMs, POD modes, and the DMD modes. For the POD and DMD modes, snapshots obtained from the nonlinear case w ere used where the elastic stop at the tip e xists. Each mode shape is normalized with respect to its infinity nor m. As can be seen from Fig. 5 (a) and (c), despite the slight differences in their amplitude v alues, the DMD modes resemble the LNMs f or the first five modes, ex cept the second mode. On the other hand, the POD modes do not resemble neither LNMs nor DMD modes f or the first fiv e modes, as seen in Fig. 5 (b). Ne xt, to see the nature of the subspaces spanned b y the modes, the pr incipal angles [ 36 ] between the subspaces are e xamined. Let L , P , and D denote the subspaces spanned b y the set of LNMs, POD modes, and DMD modes, respectiv ely , or L = Span ( L 1 , . . . , L 𝑝 ) , P = Span ( U 1 , . . . , U 𝑝 ) , D = Span ( φ 1 , . . . , φ 𝑝 ) . Then, the pr incipal angles between these subspaces w ere computed and sho wn in Fig. 6 , where Θ ( X , Y ) designates the pr incipal angle betw een tw o subspaces X and Y . From Fig. 6 , w e can see that the principal angle betw een D and L is small when the number of modes is one, which means the y are almost linearl y dependent. This ag rees with the observation that the mode shapes of the first mode f or both LNM and DMD modes w ere almost identical to each other . How ev er , as the number of modes increases, the pr incipal angle reaches almost 𝜋 / 2 when the number of modes is three, which means that the y become almost or thogonal to each other, then slightly decreases down to 1.0. The pr icipal angle betw een P and D gradually increases and become close to 𝜋 / 2 when the number of modes is nine. It means that the y are almost or thogonal with each other . On the other hand, the pr incipal angle betw een P and L is between appro ximately 0.5 and 1.5, which means their relationship lies between linearl y dependent and or thogonal. From these results, w e can see that the subspaces spanned by these sets of v ectors are not parallel to each other , i.e., LNMs, PODs, and DMDs do not f or m the same subspace. Theref ore, their capability to express the nonlinear dynamics should also be different from each other , and needs to be fur ther e xamined with f orced response calculations, as f ollo ws. 3.2 F orced response calculation using reduced order models Forced response calculations were conducted using full FE models, reduced order models with LNMs, POD modes, and DMD modes, and the results are ex amined. F or the reduced order models with DMD modes, the MOR procedure f ollow s the Algor ithm 1 . F or the other modes, the process of the Galerkin projection is rather straightf or ward, i.e., with 𝚽 = U = [ U 1 , . . . , U 𝑝 ] f or POD modes, or 𝚽 = L = [ L 1 , . . . , L 𝑝 ] f or LNMs and ˜ M = 𝚽 T M 𝚽 , ˜ K = 𝚽 T K 𝚽 , 11 -2 0 2 Displacement (m) 10 -3 -1 -0.5 0 0.5 1 Velocity (m/s) 0 0.5 1 Time (s) 0 5 10 15 Displacement (m) 10 -4 DMD(p=1) DMD(p=3) DMD(p=5) Full (64DOF) Figure 8: Phase por trait of the tip mo v ement at the resonance ˜ C = 𝚽 T C 𝚽 . A har monic f orcing was applied at the beam tip, i.e., 𝑏 ( 𝑡 ) = 𝐴 sin ( 𝜔 𝑡 ) (26) where 𝐴 = 0 . 01 N, and the steady -state response of the system was sought f or the frequency range of 𝜔 / ( 2 𝜋 ) = 0 to 120 Hz to e x cite a resonance that was e xpected to occur near the first linear natural frequency of 50.9 Hz. The cor responding external f orcing v ector is ev aluated as ˜ b = 𝚽 T b ( 𝑡 ) . The nonlinear f orce cor responding to Eq. ( 23 ) is ev aluated with phy sical DOFs that are transf or med back from the modal coordinates. The results of the f orced response calculations with full-order FEM model, the R OMs constructed by LNMs, POD modes, and DMD modes f or 𝑝 = 1 and 5 are sho wn in Fig. 7 . The displacement was ev aluated at the tip of the beam. As can be seen, there is a resonant peak at 82 Hz f or the FEM model. All R OMs with 𝑝 = 1 do not capture neither the maximum amplitude v alue nor the frequency of the resonant frequency , as sho wn in Fig. 7 . The R OMs with the LNM and the DMD mode better capture the resonance than the R OM with the POD modes, but the y are still lo wer than the one predicted b y the full model by 10%. On the other hand, the R OMs with 𝑝 = 5 capture this resonance behavior quite well f or all R OM types. This means that all three projection basis v ectors are equiv alently good to capture the resonant beha vior if sufficiently larg e number of modes are kept in the R OM. In par ticular , the resonance predicted by the DMD modes is the best among the ones considered, as seen in the enlarg ed plot sho wn in Fig. 7 . The results of full order and the R OM are considered to match quantitativ ely in this case, consider ing that the er rors between the nonlinear resonant frequencies obtained by the full order sy stem and the R OM is approximatel y 1%, and the er rors in the peak amplitude are appro ximately 6%. With these, we belive that the accuracy of the R OM is satisfactory . Moreo ver , to see the behaviors of the sys tems at the resonance, the time histories of the tip displacement and the cor responding phase por traits are ex amined f or the R OMs with the DMD alone. F igure 8 show s the time history of the tip displacement at the resonance frequency . As expected, R OMs with 𝑝 = 1 and 𝑝 = 3 do not capture the tra jector y of the full-order system well. On the other hand, the R OM with 𝑝 = 5 captures the temporal behavior of the tip displacement computed b y the full-order FEM model. From these results sho wn abov e, w e can see that DMD is able to capture the dynamics of a simple beam structure that is subjected to an elastic stop, which causes the sys tem to hav e piecewise-linearity . The accuracy of the R OM is slightl y better than those of LNM and POD modes, which hav e already been kno wn to be good basis v ectors f or MOR of linear and nonlinear structural dynamics problems. 4 F orced response of a bonded shell assembly with partial debonding with h ybrid data-driv en model-order reduction This section provides a numer ical ex ample of the f orced response problem of a bonded shell assembl y with par tial debonding. The FE model of the assembl y is sho wn in Fig. 9 . N amel y , tw o plates are bonded tog ether at one of their ends, but 75% of the bonded area is assumed to be debonded. This causes the plates to contact with each other repeatedl y if they are subject to per iodic loading, such as har monic f orcing, which results in nonlinear dynamic 12 External force Debonded nodes Fixed Fixed Figure 9: Bonded shell assembl y with par tial debonding beha vior . The shell assembly has been discretized b y MITC3 isotropic tr iangular shell elements [ 37 ], where each node of the element has three translational and two rotational DOFs. The total number of DOFs of the sy stem is 500. Both ends are assumed to be fixed, and an e xter nal f orce is assumed to be applied in the middle of the assembly . The equation of motion of the sy stem is written as a set of second-order ODEs, i.e., M ¥ v ( 𝑡 ) + C ¤ v ( 𝑡 ) + Kv ( 𝑡 ) = b ( 𝑡 ) + f ( v ) , (27) where v denotes the nodal displacement vector that contains three translational and tw o rotaitonal DOFs per node. The nonlinear f orce f ( v ) is computed based on relative nodal displacements at the debonded area. The force is assumed to be computed b y a penalty method, i.e., 𝑓 𝑖 ( v ) = − 𝑘 𝑝 max { 0 , 𝑔 ( 𝑥 𝑢 , 𝑥 ℓ ) } (28) where 𝑘 𝑝 denotes a pre-defined penalty parameter and 𝑔 represents a g ap function defined here as 𝑔 ( 𝑥 𝑢 , 𝑥 ℓ ) ≜ 𝑥 𝑢 − 𝑥 ℓ (29) where 𝑥 𝑢 and 𝑥 ℓ denote the displacements of upper and lo w er nodes on the sur faces of a specific location of the debonded areas, which are per pendicular to the debonded sur f aces. 4.1 Galerkin projection based on DMD and constraint modes With FE models containing multiple contact DOFs, it is time consuming to e v aluate the contact nonlinearity because contact f orces need to be ev aluated in phy sical coordinate sy stem. Theref ore, in order to handle the contact problems efficientl y , w e propose that the DMD modes be used with the so called constr aint modes , which can be obtained by sol ving a set of static problems where unit displacements are applied to activ e DOFs that are to be kept in the R OM, with all the other DOFs being fix ed [ 31 ]. These modes ha v e been used in the conte xt of Component Mode Synthesis to apply displacement compatibility conditions to connect multiple bodies. This time, we use them to appl y unilateral contact f orce to the body in the reduced-order space. Namel y , the dynamics of the entire sys tem is projected onto the subspace spanned by the DMD modes and the cons traint modes. Eq uation ( 3 ) is re-written as f ollow s x ( 𝑡 ) ≈ 𝑛 𝑚  𝑖 = 1 𝜂 𝑖 ( 𝑡 ) φ 𝑖 + 𝑛 𝑎  𝑖 = 1 𝑥 𝑖 ( 𝑡 ) ψ 𝑖 (30) where φ 𝑖 denotes the 𝑖 th DMD mode, 𝑛 𝑚 denotes the number of DMD modes to be kept in the R OM, ψ 𝑖 denotes the 𝑖 th constraint mode, 𝑛 𝑎 denotes the number of constraint modes, which is equal to the number of active DOFs to be kept in the R OM, 𝜂 𝑖 ( 𝑡 ) and 𝑥 𝑖 ( 𝑡 ) are the modal coordinates cor responding to the DMD modes and the constraint modes, respectivel y . The modal coordinates 𝑥 𝑖 ( 𝑡 ) beha v e like phy sical displacements in the reduced order space. This enables ones to handle contact nonlinear ity in the reduced order space without e xpanding the dynamics of the modal coordinates back to the full-order sys tem at ev er y time instant because the DOFs that are subject to the contact nonlinearity are accessible in the reduced-order space. 13 Figure 10: Linear nor mal modes of the assembly T able 2: N atural frequencies f or (A) per f ectly bonded and (B) par tiall y debonded assemblies. Mode Number (A) (Hz) (B) (Hz) De viation (%) 1 19.01 18.93 -0.40 2 57.02 56.75 -0.48 3 114.74 113.28 -1.27 4 129.60 126.87 -2.11 5 195.78 194.11 -0.85 6 311.31 304.49 -2.19 7 338.92 317.08 -6.45 8 449.00 444.39 -1.03 9 522.41 478.19 -8.46 4.2 Effects of debonding on the vibration modes T o illustrate the effects of debonding on the dynamic characteristics of the assembly , linear modal analy ses ha v e been conducted on the perfectl y bonded shell assembl y , and that with par tial debonding. Note that with the partial debonding, the contact nonlinearity at the debonded areas has not been enf orced, i.e., the sys tem is linear . Also the external f orce has been ignored f or the modal analy ses. T able 2 show s the first nine natural frequencies f or both cases. As can be seen, natural frequencies all decrease because of the debonding, which makes sense because debonding results in softening the stiffness of the sy stem. The effect of debonding is small f or the first six and the eighth modes, considering that the deviations in the natural frequencies are less than three percent. On the other hand, the natural frequencies for the se v enth and the ninth modes are greatly affected by the debonding consider ing their deviations are greater than five percent. T o fur ther e xamine these variations, the mode shapes cor responding to the sev enth and the ninth modes are sho wn in Fig. 10 f or both per f ectly bonded and par tially debonded cases. As can be seen in Fig. 10 (a), the mode sev en sho ws the twisting shape, whereas the ninth mode also sho ws the twisting, but higher -order mode shape, as shown in Fig. 10 (b). As can be seen in the mode shapes of the assembly with par tial debonding f or both modes, a gap at the debonded area appears because the amount of twist f or the upper shell is 14 10 -6 0 1.0 2.0 3.0 4.0 D i sp l acement ( m ) Figure 11: Initial deformation of the plate assembl y -5 0 5 g (m) 10 -6 -0.02 0 0.02 dg/dt (m/s) 0 0.02 0.04 0.06 Time (s) -5 0 5 Displacement (m) 10 -6 Upper node Lower node (a) Linear case allo wing interpenetration of the sur faces -5 0 5 g (m) 10 -6 -0.02 0 0.02 dg/dt (m/s) 0 0.02 0.04 0.06 Time (s) -5 0 5 Displacement (m) 10 -6 Upper node Lower node (b) Nonlinear case with contact f orce applied by the penalty method Figure 12: Phase por traits of the gap functions and time histories of the cor responding nodes f or 0 ⩽ 𝑡 ⩽ 0 . 0625 s smaller than that f or the low er shell at the debonded area. This is because the stiffness with respect to the twisting motion decreased due to the debonding. 4.3 Dynamic mode decomposition T o obtain DMD modes, snapshots of the displacements need to be obtained. As in the pre vious numer ical e xample, the snapshots are obtained by solving initial value problems. This time, instead of applying impulsiv e f orcing, initial displacements were applied to the assembly and the resulting free responses w ere computed. The initial displacements were obtained by appl ying static f orces to open-up debonded areas and solving the cor responding static problem for the def or mation of the plate assembly . The shape of the applied initial def or mation is sho wn in Fig. 11 . This initial condition is emplo yed to ensure the localized motion of the debonded areas is captured in the snapshot. With the initial condition, time integration has been conducted b y Ne wmark - 𝛽 method f or 0 ⩽ 𝑡 ⩽ 1 . 0 s. The time integ ration has been conducted with a fix ed time increment of Δ 𝑡 = 3 . 0519 × 10 − 5 s, which means that the dynamic response of frequency up to 16.384 kHz are e xpected to be captured, which is enough f or computing DMDs with much lo w frequencies. The obtained time histor ies and the cor responding phase por traits of the gap functions are shown in Fig. 12 , where the gap function is defined as Eq. ( 29 ) . Figure 12 (a) show s a linear response where no contact f orce is imposed on the potentially contacting nodes ev en when the y are supposed to be in contact. As can be seen in the time histories of the nodes, both nodes oscillate with high frequency with lo wer node’ s amplitude being smaller than that f or the upper node. The cor responding gap function on the left sho ws an elliptical orbit with shr inking major and minor ax es as time proceeds. This is a typical behavior of damped linear sys tems. Figure 12 (b) sho ws the response of the nodes with PWL contact f orce being enf orced. As can be seen, the upper node sta ys abo v e the low er node due to the contact f orce that occurs when they are in contact, which hinders the inter penetration of the nodes. The trajectory of the gap function in the phase por trait is asymmetr ic with respect to 𝑔 = 0 , because gap function is f orced to be 𝑔 ⩾ 0 with slight penetration due to the penalty coefficent of finite magnitude. 15 0 100 200 300 400 500 Frequency (Hz) 10 -10 10 -5 Displacement (m) Linear Nonlinear Figure 13: R esults of FFT of the displacement at the f orcing point 0 100 200 300 400 500 Frequency (Hz) 10 -10 10 -8 10 -6 Displacement (m) 500 1000 1500 2000 Sampling Frequency (Hz) Lowest Highest DMD eigenvalue rank (a) Linear (no contact) 0 100 200 300 400 500 Frequency (Hz) 10 -10 10 -8 10 -6 Displacement (m) 500 1000 1500 2000 Sampling Frequency (Hz) Lowest Highest DMD eigenvalue rank (b) Nonlinear (contact at debonded area) Figure 14: Pseudo-stability diag ram of the DMD eig en values, along with the a v erage Fourier spectrum of the displacement of all DOFs. — : FFT spectra of all DOFs, — : a v erag e spectr um, × : DMD eig en values. The results of the FFT of the displacement at the f orcing point f or both linear and nonlinear cases are sho wn in Fig. 13 . From these results, w e can see that the peak frequencies belo w 200 Hz do not chang e much regardless of the state of the debonded areas. This ag rees with the fact that natural frequencies of the first fiv e modes did not g reatl y de viate due to the debonding, as sho wn in T able 2 . Also, we can see that the peak frequency at 304 Hz f or linear case increases to 306 Hz f or the nonlinear case, which produced 0.6% de viation. The peak frequency at 317 Hz f or the linear case, which cor responds to the sev enth mode, increases up to 327 Hz f or the nonlinear case, with its amplitude decreased b y 13%. It is because, f or the sev enth mode, the inter penetration of the debonded areas is so significant without the nonlinear boundary condition that the enf orcement of the nonlinear boundar y condition at the debonded areas greatly influences the motion of the debonded areas, which results in the suppression of the e xcitation of such modes. With the computed snapshot, the DMD based on the Algorithm 1 has been applied to both linear and nonlinear cases f or increasing sampling frequency to not only capture enough dominant dynamic modes but also to see their ph y sical significance. The results of DMD are shown in Fig. 14 . A gain the DMD spectra are sho wn in conjunction with the FFT spectra of all DOFs and their av erag es. As can be seen in Fig. 14 (a), the DMD spectra ag ree with the peak frequencies in the a v erage FFT spectr um. Moreov er , as discussed in Section 3 , they coincide with the sy stem’ s 16 T able 3: DMD frequencies extracted from the free response of the shell assembl y by initial displacement No. Linear (Hz) Nonlinear (Hz) De viation (%) 1 19.09 18.95 -0.74 2 56.73 56.88 0.27 3 113.25 113.64 0.34 4 126.82 127.56 0.58 5 194.05 195.07 0.52 6 304.22 306.28 0.67 7 316.69 331.19 4.57 8 443.53 445.28 0.39 9 477.09 491.20 2.95 Figure 15: Mode shapes of DMD modes ordered based on frequencies of the cor responding DMD eig en values natural frequencies when the DMD is applied to the linear response. Figure 14 (b) show s the DMD spectra computed from the nonlinear case where the contact force at the debonded areas was taken into account. As seen, the dominant DMD frequencies coincide ag ain with the peak frequencies obtained from the FFT . This time, ho we ver , they are not necessarily natural frequencies because the response comes from the underl ying piecewise linear sys tem. In addition, man y DMD eigen values appear between significant peaks with lo w rank among the ones computed. This again sho ws the characteristics of the DMD spectra f or nonlinear responses. T able 3 sho ws the frequencies obtained from the DMD eigen values, which are the first nine dominant ones ordered based on the cor responding frequencies. As can be seen, most frequencies are not different from each other betw een the linear and the nonlinear cases, e x cept the mode sev en, whose mode shape f eatures localized vibration near the debonded area. 4.4 Model order reduction b y DMD With the obtained DMD modes, Galerkin projection has been applied to Eq. ( 27 ) , based on Eq. ( 30 ) . Mode shapes used f or the projection are shown in Fig. 15 . In addition, f or the sake of comparison, the projection has been applied using LNMs and POD modes. The cor responding shapes of those modes are sho wn in Figs. 16 and 17 . As can be seen, the mode shapes of DMD modes resemble those of LNMs, with slight differences especially in the vicinity of debonded regions. On the other hand, the POD modes do not necessarily resemble neither the LNMs nor the DMDs. Instead, it emphasizes the localized motion at the debonded area. For this numerical e xample, 19 constraint modes cor responding to 18DOFs of the nodes on the debonded sur faces and a DOF cor responding to the ex citation point w ere used f or all cases, i.e., 𝑛 𝑎 = 19 . Firs t, the obtained R OMs are compared in ter ms of the er rors in their natural frequencies. Namel y , the eig en values 17 Figure 16: Mode shapes of linear nor mal modes ordered based on frequencies Figure 17: Mode shapes of POD modes ordered based on the magnitude of the cor responding singular values 𝜎 𝑖 of the R OMs w ere computed b y varying the number of modes in the R OMs. N ote that since 𝑛 𝑎 = 19 f or all cases, the number of modes in Eq. ( 30 ) , or 𝑛 𝑑 is varied to see its effects on the R OM accuracy . Also, R OMs with the fix ed-interface normal modes w ere considered, which is widely -used accurate model-based MOR and is kno wn to be the Craig-Bampton (CB) method. The results are shown in Fig. 18 . As can be seen, LNM conv erg es the fas test with 28 modes. The R OMs with the POD modes also sho w con v erg ence with 28 modes, but the a verag e er ror is the larg est among the ones compared. The er rors computed with the DMD R OMs become appro ximately 10 times smaller than those obtained with the POD R OMs. The CB-R OM also giv es good results. Consider ing that both CB and LNM R OMs are model-based, i.e., all projection basis v ectors are computed directl y from the sy stem matr ices M and K , their capability to capture the eigenspace of M and K w ell is not surpr ising. On the contrar y , both POD-R OMs and DMD-R OMs are data-driven, although the y are compensated by the model-based constraint modes. Theref ore, they are not guaranteed to be capable of representing the eig enspace of the sys tem matr ices. Ho w ev er , based on these results, we can see that the y can produce accurate eigen values of the or iginal sy stems if w e include enough number of modes in the projection basis. 4.5 F orced response anal y sis Using the FEM model and the R OMs, f orced response analy sis under har monic loading has been conducted where 𝑏 ( 𝑡 ) = 𝐴 sin ( 2 𝜋 𝑓 𝑡 ) . Time integration calculations ha v e been conducted on the equations f or 295 Hz ⩽ 𝑓 ⩽ 320 Hz to e x cite the sixth and se v enth mode until the responses reach the steady -state response. 18 20 25 30 35 40 Number of modes in the ROM 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 Average error in natural frequencies (%) Craig-Bampton LNM DMD POD Figure 18: Comparison of er rors in the natural frequencies f or the R OMs compared The amplitude of the displacement response measured at the e x citation point is plotted in Fig. 19 f or FEM, LNM, CB, DMD, and POD R OMs. Also, responses of the perfectl y bonded assembly and of par tially debonded but no contact f orce is imposed are shown in Fig. 19 . From the plots, w e can see that the response of the per f ectly bonded assembl y sho ws the larg est resonant frequency at 311 Hz. On the other hand, par tially debonded plate assembly sho ws the smallest resonant frequency at 304 Hz, if no contact force is applied. How ev er , if the contact f orce that stems from the contact at the debonded region is applied, i.e., if the sys tem is PWL, the responses show the resonance at around 307 Hz. Contrary to the e xpectation from the FFT spectra sho wn in Fig. 14 and from the linear response shown in Fig. 19 , the resonant peak cor responding to se v enth mode did not appear . Overall, the R OMs predicted the resonance accuratel y . In ter ms of the resonance height, CB-ROM produced the best results. Both the LNM and the DMD R OMs produced accurate but slightl y larg er height than that predicted by the CB-R OM. The resonance height computed by the POD-R OM is the highest and the least accurate. Figure 20 show s the time histories of the gap functions at the resonance computed b y FEM, LNM, CB, DMD, and POD-R OMs. As can be seen, all R OMs accurately predicted the response in compar ison with that computed b y FEM, but produced slightl y larg er amplitude than that computed by the FEM. The responses comptued b y DMD and LNM produced smooth response that resemble the response computed by the FEM. The responses computed b y POD and CB-R OMs contains high-frequency components that do not e xist in the one obtained b y FEM. Further more, to visualize the capability of the R OM to represent high frequency components contained in the nonlinear response especially with the DMD modes, spectrogram of the displacement at the ex citation point and the gap function w as computed by taking FFT on the responses and shown in Fig. 21 and Fig. 22 . The color represents displacement sho wn in decibel where 1 𝜇 m equals 0dB. First, the response at the e x citation point is discussed. The one predicted by the FEM is shown in Figure 21 (a). It sho ws a strong line cor responding to the first-order of the e xcitation frequency . The second and the third order lines are also obser v ed, though the y are much weak er than the first-order line. Also, there are hor izontal lines sho wing relativ ely strong response at around 20 Hz, and 120 Hz, which cor respond to the first and the f our th modes, respectiv ely , as shown in T able 3 . N ote that unlike the order lines, the y are not dependent on the e xcitation frequency . When the e xcitation frequency is at around 307 Hz, there is a resonance, as it w as shown in Fig. 19 . Ov erall, Figure 21 (a) rev eals that man y linear nor mal modes and motions with frequencies of integer multiples of e x citation frequency occur . Figure 22 (a) sho ws the spectrogram of the displacement ev aluated at the e x citation point, predicted b y the DMD-R OM. Ov erall, it agrees quite well with the one predicted b y the FEM. W e can see the strong first-order line and weak er second and third order lines in Fig. 22 (b). Moreo ver , the horizontal lines cor responding to the frequencies of the first and the fourth modes can also be obser v ed. 19 295 300 305 310 315 320 325 330 335 Fre q uenc y (Hz) 0 1 2 3 4 5 6 7 Displacement amplitude (m) 10 -5 Linear (no contact) Perfectly bonded LNM ( n m =20) FEM CB ( n m =20) DMD ( n m =20) POD ( n m =20) Figure 19: R esults of f orced response anal ysis 0 0.2 0.4 0.6 0.8 1 Normalized time 0 2 4 6 Displacement (m) 10 -5 LNM ( n m =20) FEM Figure 20: T ime histories of gap functions at the resonance (307 Hz) The R OM also captures the nonlinear nature of the response at the resonance that contains many linear nor mal mode components and motions cor responding to the integer -multiples of the e x citation frequency . Second, the response of the gap function is discussed. The one predicted by the FEM, ev aluated at a location on the debonded region is shown in Fig. 21 (b). This time, w e can see strong, frequency-dependent, first, second, and third order lines in the response. U nlike the response at the e xcitation point, frequency components cor responding to the linear nor mal modes are obser v ed onl y when the e x citation frequency is in the vicinity of the resonant frequency . Fig. 22 (b) sho ws the one computed b y the DMD-R OM. As can be seen, it agrees quite well with the one computed b y the FEM. From these, w e can see that the response at the debonded area is s trongly nonlinear that entails motions with integer -multiples of the e xcitation frequencies. This is caused by repetitiv e opening and closing of debonded surfaces. Lastl y , the benefit of the proposed MOR methodology in ter ms of the computational time is discussed. The f orced response calculations were conducted b y the Ne wmark 𝛽 time integration method with the R OMs and the FEM f or 𝑓 = 307 Hz and 0 ⩽ 𝑡 ⩽ 0 . 208 s, where the number of time steps was fix ed to 65,536 f or all cases, resulting in the fixed time step size of Δ 𝑡 = 3 . 17 × 10 − 6 s . Also, 𝑛 𝑚 = 20 f or all the R OMs considered. The f orced response calculations w ere repeated five times with the same conditions and the a v erag e computational time f or 20 300 310 320 330 Excitation fre q uenc y (Hz) 0 200 400 600 800 1000 Frequency (Hz) -80 -60 -40 -20 0 20 Displacement (dB re 1.0 m) (a) Ex citation point 300 310 320 330 Excitation fre q uenc y (Hz) 0 200 400 600 800 1000 Frequency (Hz) -80 -60 -40 -20 0 20 Displacement (dB re 1.0 m) (b) Gap function Figure 21: Spectrogram of displacement at the e xcitation point and gap function computed with FEM 300 310 320 330 Excitation fre q uenc y (Hz) 0 200 400 600 800 1000 Frequency (Hz) -80 -60 -40 -20 0 20 Displacement (dB re 1.0 m) (a) Ex citation point 300 310 320 330 Excitation fre q uenc y (Hz) 0 200 400 600 800 1000 Frequency (Hz) -80 -60 -40 -20 0 20 Displacement (dB re 1.0 m) (b) Gap function Figure 22: Spectrogram of displacement at the e xcitation point and gap function computed with DMD-R OM each R OM per frequency w as computed. The measurements were conducted on a w orkstation with the CPU of AMD Threadr ipper PR O 3975WX (4.2 GHz) and 128GB of RAM. During the measurements, the application of the numerical simulation used only a single core of the CPU f or all cases. The results are shown in T able 4 where CPU time f or each computational step is sho wn f or FEM, LNM, CB, DMD and POD. The total CPU time sho wn in the bottom row is sho wn with the assumption that the CPU time f or f orced response calculation per frequency is constant throughout the frequency range. As can be seen in the table, although the number of modes was chosen to be the same, there are slight differences betw een the computational times f or the R OMs. In par ticular , LNM and CB sho w smaller computational times than DMD and POD. This is because the y are classical model-based R OM that does not require the computation of the snapshots. On the other hand, both DMD and POD take longer CPU time than the model-based R OMs but they are comparable. This means that with the proposed method, w e are able to g enerate R OMs with a reasonable amount of time and accuracy . Further more, we can see that the computations with the R OMs w ere much fas ter than that with the FEM. Indeed, the computations with DMD-R OM w ere appro ximately 100 times f aster than those with the FEM model. From these, we can see the efficiency of using the R OMs in compar ison with the full-order FEM model when predicting the f orced response of the sys tem with PWL nonlinear ity . 21 T able 4: Compar ison of computational time required to complete a set of f orced response calculations at the resonance (second) FEM LNM CB DMD POD (A) Snapshot extraction N/A N/A N/A 6 . 7 × 10 3 6 . 7 × 10 3 (B) DMD (POD) calculations N/A N/A N/A 5 . 9 × 10 1 1 . 2 × 10 − 1 (C) Galerkin projection N/A 7 . 2 × 10 − 3 1 . 2 × 10 − 2 5 . 8 × 10 − 3 4 . 8 × 10 − 3 (D) Forced response calculation (per frequency) 7 . 6 × 10 3 2 . 7 × 10 1 1 . 8 × 10 1 2 . 2 × 10 1 2 . 7 × 10 1 T otal (A)+(B)+(C)+(D) × 𝑛 𝑓 9 . 7 × 10 5 3 . 5 × 10 3 2 . 3 × 10 3 9 . 6 × 10 3 1 . 0 × 10 4 5 Conclusion In this paper, a model order reduction method f or piecewise-linear sys tems based on data-driven dynamic mode decomposition has been proposed. The k e y idea of the proposed method is that the snapshots are obtained based on impulse response or response due to initial def or mation of the system, which results in r ich spectral contents in the response due to the piecewise-linear nonlinear ity . This giv es one the capability to capture the complicated dynamics due to PWL nonlinearity , which nor mall y contains man y frequency components. For the sys tems with man y DOFs with piecewise-linear nonlinear ity , it was proposed that DMD modes be used with the constraint modes to be able to handle the nonlinear ity efficiently . The proposed methodology has been applied to a f orced response problems of a cantile v ered beam with an elastic stop at the tip, and a bonded shell assembly with par tial debonding subjected to harmonic loading. The f orced response problems were sol v ed with the proposed method, and the results ag reed well with the ones obtained b y the full-order FEM models. Statements and Declarations Funding This project w as suppor ted in par t b y Japan Society f or the Promotion of Science (JSPS), Grant-in-Aid f or Scientific R esearch(C), grant number JP20K11855. Compe ting interests The authors hav e no relev ant financial or non-financial interests to disclose. A uthor contributions Both authors contr ibuted to the study conception, design and analy sis. Simulations w ere implemented and perf or med mainl y by A. Saito. The first draft of the manuscr ipt was wr itten b y A. Saito and M. T anaka commented on the draft of the manuscr ipt. All authors read and appro v ed the final manuscript. Data a vailability The dataset generated dur ing and/or analy sed dur ing the cur rent study are a vailable from the cor responding author upon reasonable request. R ef erences [1] Akira Saito, Matthew P . Castanier , and Chr istophe Pier re. Estimation and veering analy sis of nonlinear resonant frequencies of crack ed plates. Journal of Sound and Vibr ation , 326(3–5):725 – 739, 2009. [2] Paolo Casini and F abr izio V estroni. Character ization of bifurcating non-linear normal modes in piece wise linear mechanical sy stems. International Journal of N on-Linear Mechanics , 46(1):142 – 150, 2011. 22 [3] Mohammad A AL-Shudeifat and Eric A Butcher . On the dynamics of a beam with switching crack and damag ed boundar ies. Jour nal of Vibr ation and Control , 19(1):30–46, 2013. [4] V .N. Burla yenk o and T . Sadow ski. Finite element nonlinear dynamic analy sis of sandwich plates with partially detached facesheet and core. Finit e Elements in Analysis and Design , 62:49 – 64, 2012. [5] S. W . Sha w and P . J. Holmes. A periodically forced piece wise linear oscillator . Jour nal of Sound and Vibr ation , 90(1):129–155, 1983. [6] Akira Saito. Nonlinear resonances of chains of thin elastic beams with inter mittent contact. Jour nal of Computational and Nonlinear Dynamics , 13(8):081005–081005–10, 2018. [7] K ohei Noguchi, Akira Saito, Meng-Hsuan Tien, and Kiran D’Souza. Bilinear Sys tems With Initial Gaps In v olving Inelastic Collision: F orced Response Exper iments and Simulations. Jour nal of Vibr ation and Acous tics , 144(2):021001 (9 pages), 2022. [8] S. F . Masr i, R. K. Miller , H. Sassi, and T . K. Caughe y . A method f or reducing the order of nonlinear dynamic sy stems. Journal of Applied Mec hanics , 51(2):391–398, 1984. [9] S. F . Masr i, J. P . Caffre y , T . K. Caughey , A. W . Sm yth, and A. G. Chassiak os. A general data-based approach f or dev eloping reduced-order models of nonlinear mdof sy stems. Nonlinear Dynamics , 39(1):95–112, 2005. [10] D. Jiang, C. Pierre, and S.W . Sha w . Lar ge-amplitude non-linear nor mal modes of piecewise linear sys tems. Journal of Sound and Vibr ation , 272(3–5):869 – 891, 2004. [11] Meng-Hsuan Tien and Kiran D’Souza. A g eneralized bilinear amplitude and frequency appro ximation f or piece wise-linear nonlinear sys tems with gaps or prestress. Nonlinear Dynamics , 88(4):2403–2416, 2017. [12] Meng-Hsuan Tien and Kiran D’Souza. T ransient dynamic anal ysis of crac ked structures with multiple contact pairs using generalized hsnc. Nonlinear Dynamics , 96(2):1115–1131, 2019. [13] Giuseppe Quaranta, W alter Lacarbonara, and Sami F . Masri. A revie w on computational intelligence f or identification of nonlinear dynamical sy stems. Nonlinear Dynamics , 99(2):1709–1761, 2020. [14] P eter J. Schmid. Dynamic mode decomposition of numerical and experimental data. Jour nal of Fluid Mec hanics , 656:5–28, 8 2010. [15] Guillaume Jourdain, Lars-Er ik Er iksson, Su Ho Kim, and Chae Hoon Sohn. Application of dynamic mode decomposition to acoustic-modes identification and damping in a 3-dimensional chamber with baffled injectors. Journal of Sound and Vibr ation , 332(18):4308 – 4323, 2013. [16] Diana A. Bistrian and Ionel M. N av on. The method of dynamic mode decomposition in shallo w water and a swirling flow problem. International Journal for N umerical Methods in Fluids , 83(1):73–89, 2017. [17] Franck Richecoeur , Lay al Hakim, Antoine Renaud, and Laurent Zimmer . Dmd algor ithms f or e xperimental data processing in combustion. In Proceeding of the 2012 Summer Pr ogr am , pag es 459–468. Center f or T urbulence R esearch, Stanf ord U niversity , December 2012. [18] Mohd Y Ali, Anshuman Pande y , and James W Gregor y . Dynamic mode decomposition of fast pressure sensitiv e paint data. Sensors , 16(6):862, 2016. [19] G Berk ooz, P Holmes, and J L Lumley . The proper orthogonal decomposition in the anal ysis of turbulent flo w s. Annual Review of Fluid Mechanics , 25(1):539–575, 1993. [20] Gaetan Kersc hen, Jean-Claude Golin val, Alex ander F . V akakis, and La wrence A. Bergman. The method of proper or thogonal decomposition f or dynamical characterization and order reduction of mechanical sy stems: An o v erview . Nonlinear Dynamics , 41(1):147–169, 2005. 23 [21] S. L. Br unton and J. N . Kutz. Data-driv en science and engineering . Cambridge Univ ersity Press, 2019. [22] Alessandro Alla and J. Nathan Kutz. Nonlinear model order reduction via dynamic mode decomposition. SIAM Journal on Scientific Computing , 39(5):B778–B796, 2017. [23] Riasat Khan and Kw ong T . Ng. DMD-Galerkin model order reduction f or cardiac propag ation modeling. A CES Journal , 33(10):1096–1099, 2018. [24] Barbara Cunha, Chr istophe Droz, A bdelmalek Zine, Stéphane Foulard, and Mohamed Ichc hou. A revie w of machine lear ning methods applied to structural dynamics and vibroacoustic, 2022. [25] C. Hari Manoj Simha and Mohammad Biglarbegian. A note on the use of dynamic mode decomposition in mechanics. Mechanics Researc h Communications , 120:103848, 2022. [26] Jonathan H. T u, Clarence W . R ow ley , Dirk M. Luchtenburg, Stev en L. Br unton, and J. Nathan Kutz. On dynamic mode decomposition: Theory and applications. Journal of Computational Dynamics , 1(2):391–421, 2014. [27] Jonathan H. T u. Dynamic mode decomposition: theor y and applications . PhD thesis, Princeton U niversity , 2013. [28] Akira Saito and T omohiro Kuno. Data-dr iv en e xper imental modal analy sis by dynamic mode decomposition. Journal of Sound and Vibr ation , 481:115434, 2020. [29] A. N. Thite and D. J. Thompson. The quantification of structure-borne transmission paths by inv erse methods. par t 1: impro v ed singular value rejection methods. Jour nal of Sound and Vibr ation , 264:411–431, 2003. [30] Bar t Peeters, Her man V an der A uw eraer , Patrick Guillaume, and Jan Leur idan. The polymax frequency -domain method: A new standard f or modal parameter estimation? Shoc k and Vibration , 11(3-4):395–409, 2004. [31] R. R. Craig and M. C. C. Bampton. Coupling of substructures f or dynamic analy ses. AIAA Jour nal , 6(7):1313–1319, 1968. [32] Saif on Chaturantabut and Danny C. Sorensen. Nonlinear model reduction via discrete empir ical inter polation. SIAM Journal on Scientific Computing , 32(5):2737–2764, 2010. [33] Jan S. Hestha ven, Cecilia Pagliantini, and Nicolò Ripamonti. A daptiv e symplectic model order reduction of parametric particle-based v lasov -poisson equation, 2022. [34] Felix Fr itzen, Ber nard Haasdonk, Da vid R yc kel ynck, and Sebastian Schöps. An algorithmic comparison of the h yper -reduction and the discrete empir ical inter polation method f or a nonlinear ther mal problem. Mathematical and Computational Applications , 23(1), 2018. [35] D. J. Ewins. Modal T esting: Theor y, Pr actice and Application Second Edition . Wile y , 2009. [36] G. H. Golub and C. F . van Loan. Matrix Computations Third Edition . The Johns Hopkins U niversity Press, 1996. [37] Phill-Seung Lee and Klaus-Jürg en Bathe. Dev elopment of MITC isotropic tr iangular shell finite elements. Computer s & Structures , 82(11):945–962, 2004. 24