We advocate for the use of Dynamic Mode Decomposition (DMD) in mechanics. Snapshots of the temper- ature or displacement ﬁelds are obtained; and the DMD algorithm, which, remarkably, requires a few lines of code, is applied to matrices obtained by time shifting the snapshots. Only the higher-order modes are retained to obtain a lower-rank matrix and used for analysis. We show that DMD can be used to: predict temperature ﬁelds during transient heat conduction, estimate natural frequency of a linear elastic structure; identify the onset of plasticity in an elastic-plastic plate with a hole; and obtain a reduced-order model that describes the large deformation of a non-linear elastic structure. Keywords: Data-driven; Dynamic Mode Decomposition; Solid Mechanics; Reduced-Order Model 1. Introduction 1 Eulerian measurements or results from Compu- 2 tational Fluid Dynamics (CFD) computations of 3 velocity, temperature, and pressure in ﬂow ﬁelds, 4 can be post-processed using a popular data-driven 5 algorithm called Dynamic Mode Decomposition 6 (DMD). Developed by Schmid [10] and Rowley 7 et al. [7], and applied to ﬂow in channels and jets by 8 Rowleyetal.[7]andSchmidetal.[9], anditisnote- 9 worthy that within a decade of its introduction it 10 hasbeenappliedinthermo-ﬂuids, neuroscience, im- 11 age processing, and epidemiological studies; Brun- 12 ton and Kutz [3]. The most recent version of the al- 13 gorithm, theoretical background, limitations of the 14 method,itsconnectionwiththeKoopmanoperator, 15 openquestions, andassociatedbibliographymaybe 16 found in the monograph by Kutz et al. [6] and text 17 by Brunton and Kutz [3]. 18 At its essence, it takes as input a sequence of 19 temporal snapshots of ﬁeld variables (velocity tem- 20 perature, displacement measurements), and sets up 21 a data matrix by time shifting the snapshots. The 22resultant data matrix may be excessively large, es- 23 pecially in CFD computations. However, the DMD 24 algorithm carries out a Singular Value Decomposi- 25 tion (SVD), which allows selection and retention of 26 only the dominant modes in the system. Retaining 27 onlythedominantmodes, leadstoalower-rankma- 28 trix,which is signiﬁcantly smaller than the larger 29 data matrix and whose eigenvalues and eigenvec- 30 tors describe the dynamics of the system. For lin- 31 ear systems, the DMD algorithm obtains the lead- 32 ing eigenvalues and eigenvectors; and for periodic 33 data, the decomposition is equivalent to a discrete 34 Fourier transform; Rowley et al. [7]. 35 Despite its popularity, DMD has not been ap- 36 plied by the solid mechanics community save two 37 instances. Saito and Kuno [8] applied it to mea- 38 surements on a vibrating cantilever and obtained 39 its natural frequency and damping coeﬃcient. Das 40 et al. [4] used DMD as well as Proper Orthogonal 41 Decomposition (POD), with states obtained from 42 ﬁnite element computation, to model the deformed 43 state of a pneumatic ﬁnger; they, however, used 44 Preprint submitted to Mechanics Research Communications March 1, 2023 displacements and velocities in their representation 45 of states. In solid mechanics, with the prevalence 46 of Digital Image Correlation (DIC) systems, snap- 47 shots of displacements are readily obtained, but not 48 velocities. 49 This letter demonstrates the use of a data-driven 50 approach in typical mechanics problems: transient 51 heat conduction; vibration of a circular plate; iden- 52 tifying onset of plasticity in a elastic-plastic struc- 53 ture; and reduced-order model for a non-linear elas- 54 tic structure undergoing large deformation. Sys- 55 tems in the ﬁrst two instances are linear. Non- 56 linearity in the third is owing to the material re- 57 sponse, and in the last because of the material re- 58 sponse andlarge deformation. 59 2. Background on DMD 60 Notation and treatment in this section is similar to that in Brunton and Kutz [3]. Assume that the data (computational or experimental) is collected from a dynamical system dx=dt=f(x;t;); where xis a state vector at time t,fdenotes the dynam- ics, andrepresents parameters of the system. The DMD approach assumes local linearity of the above and takes dx dt=Ax; (1) which has the solution x(t) = exp( At)x(0), where 61 x(0)is the initial condition. The modal ma- 62 trix ofAcomprising its eigenvectors is M= 63 [1;  2::: n], and it is a standard result that 64 x(t) =Mexp( t)M1x(0), where is the diago- 65 nal matrix containing eigenvalues of A; see section 66 10 of Brogan [2], for instance. 67 The DMD approach solves Eq. (1) by obtain- ing snapshots xof the system f(). Assuming that there arem+1snapshots (vectors) each containing nmeasurement that describe the state of the sys- tem, time-shift them and insert as columns in Xof sizenmandX0of sizenmto get X= [x1x2xm] and (2)X0= [x2x3xm+1]: (3) Laterweshowamodiﬁedstrategyforassembling X for systems with standing waves; see Tu et al. [12] in this regard. Here, we assume that the snapshots are sampled at constant intervals t. The essence of DMD is determination of a best-ﬁt operator A that advances the initial state to the next as X0=AX: (4) Using the Moore-Penrose inverse, denoted by y, leads to A=X0Xy: (5) For connections of Awith the Koopman operator, see Brunton and Kutz [3]. The states are propa- gated through an operator Athat transitions the states from snapshot to snapshot. Thus, one can write the equation for any general state, xk+1, as xk+1Axk: (6) Comparing with the solution for Eq. (1) we can view the above as a discrete-time system sampled every t, andA= exp(At); Kutzetal.[6]. Using Eq. (6) and noting that Acan be used in an iter- ative fashion, with a knowledge of the initial state x1, to predict the terminal state xmusing X= x1Ax 1Am1x1 (7) The ﬁrst step in the DMD algorithm is to compute an SVD ofXas XeUeeV; (8) eU2Cnr,e2Crr, andeV2Cmr;andrm is either the exact or approximate rank of the data matrixX. Here,denotes the complex conjugate transpose. Columns of eUare the proper orthogo- nal modes of X(see section 3.4 of monograph by Holmes et al. [5]) and satisfy eUeU=eUeU=I. Similarly,eVeV=eVeV=I. Taking the Moore- 2 Penrose inverse of Xin Eq. (8) yields Xy= eVe1eU;where1denotes conventional inverse, and combining with Eq. (5) yields Aas A=X0eVe1eU: Throughasimilaritytransformation, theexpression forAis projected on to the rPOD modes, Schmid [10], as eA=eUAeU=eUX0eVe1; (9) which has the same non-zero eigenvectors as A. We chooserby trial-and-error, and this subjectivity is a characteristic of all dimensionality reduction schemes. Matrix eU, can be used to reconstruct he full statexk, Brunton and Kutz [3], or Chap. 13 of Smith [11] as xk=eUexk: (10) Eigenvectors and eigenvalues of eAcan be obtained througheAW =W, where andWare the eigenvalues and eigenvectors of eA, respectively. Schmid [10] proposed the projected DMD modes of sizenrusing Eq. (10) as =eUW; how- ever, these are not guaranteed to be exact and Tu et al. [12] proposed the eigenvectors of the high- dimensional system as: =X0eVe1W; (11) and showed that A=; see also page 238 of Brunton and Kutz [3]. That is, the DMD modes are the eigenvectors of the higher-dimensional ma- trix corresponding to the eigenvalues in . Now deﬁne a transformation xk=bk, which implies bk+1=bk; using the latter two expressions and the eigenvalue equation in Eq. (6) gives xk+1= Abk=bk=kb1. In the last expression, we reference the initial state; and b1=yx1andwe can write the expression as xk+1=rX j=1 jk jb1; (12) where jandjare the eigenvectors and eigenval- ues, respectively. By analogy with the solution of Eq. (1), and taking !j=log(j)=t, and passing to continuous limit, we get x(t) =rX j=1 je!jtbj=exp( t)b1;(13) see also page 240 of Brunton and Kutz [3] and in 68 passing we note that ex1=Wb1. 69 3. Transient Heat Conduction 70 We consider the application of DMD to the one- dimensional transient heat-conduction problem for which an analytical solution exits: @u @t=k@2u @x2; t> 0;0<x<L; whereuis temperature, ttime,xspace, and kthe 71 thermal conductivity. With the boundary condi- 72 tionsu(0;t) =u(L;t) = 0fort > 0; and ini- 73 tial condition u(x;0) = 100 for0< x < L: 74 Takingk= 0:01,L= 1;obtains a solution of 75 u(x;t) =P1 n=0Bnsinnx Lexp(k(n=L )2t);where 76 Bn= 200(1cosn)=n: 77 Assuming t= 2=10s and a spatial frequency 78 ofx= 1=100m, we get 10 snapshots with each 79 containing 100 measurements. These are post- 80 processed using the algorithm described in Section 81 2. MatrixXhas a size of 1009 andX0has a size 82 of 1009, and it follows that Ahas a size of 100 83 100. Using r= 2andeAof sizerr, we get an 84 approximate reconstruction shown in Fig. 1 along 85 with the error of the reconstruction juDMDuj=u 86 , which is of the order of 10%. There is a spatial pe- 87 riod to the error, which is on account of the spatial 88 modes of the decomposition. To see this we plot 89 3 10.5 xr= 2 00t1200 0100 2uDMD 10.5 x00t10.1 0.05 0 2relative error x0 0.5 1φj -0.200.210.5 xr= 4 00t10200 100 2uDMD 10.5 x00t1×10-3 2 1 0 2relative error x0 0.5 1φj -0.200.2Figure 1: Left: r= 2, reconstructed temporal snapshots of temperature, associated error, and spatial mode shapes jof the last snapshot . Right: r= 4, reconstruction, error and spatial mode shapes. Note periodicity in error owing to spatial frequency of the mode shapes. the components of matrix , that has a size of 100 90 r. Only the mode shapes of the last snapshot are 91 shown. 92 When the dimension of eAis increased to r= 4, 93 by retaining more modes, the accuracy of the sur- 94 face representation improves, and the error displays 95 more spatial periodicity as shown in the plots in 96 Fig. 1 (right side). 97 4. Natural Frequency of Circular Plate 98 Consider vibrations of a circular plate, shown in 99 Fig. 2, of radius 200 mm and thickness 10 mm ﬁxed 100 along its perimeter. A quarter-symmetry model 101 wasconstructedintheAbaqusﬁniteelementsolver, 102 and meshed using linear, four-noded tetrahedra, 103 and the Young’s modulus and Poisson’s ratio were 104 taken to be 200 GPa and 0.33, respectively. A mesh 105 size of 6 mm was used and found to guarantee con- 106 vergence. 107 On account of needing time dependent displace- 108 ments, we excited the plate with a sinusoidal veloc- 109ity at its center; and the excitation frequency was 110 taken to be 622 Hz which is close to the ﬁrst natu- 111 ral frequency of the plate. The latter allowed us to 112 carryoutanimplicittransientsolution. Inpractice, 113 a hammer strike will be used to excite several fre- 114 quencies in the plate and drive the resonant modes, 115 and though explicit time integration can be used to 116 model the latter, the calculations are time consum- 117 ing. The present approach serves the purpose for 118 our demonstration. 119 Python scripts were used to output snapshots of nodal coordinate histories zwith t= 0:5ms. For linear problems with standing waves, the DMD method fails and the snapshots have to be stacked and time shifted (also called stack shifting). This increases the rank of the data matrix allowing the computation of complex conjugate eigenvectors; see 4 Fixed boundary 2.5 ms 4.5 msFigure 2: Circular plate clamped at the edges and loaded with a sinusoidal velocity condition. Contour plots of zdis- placement at two instants. page 244 of Brunton and Kutz [3]. X=" x1; :::;xm x2; :::;xm+1# and X0=" x2; :::;xm+1 x3; :::;xm+2# : A total of 18,999 nodes were used to mesh the 120 plate and 10 snapshots were recorded leading to a 121 matrixA, on account of the stacking, of size 75,996 122 75,996; in contrast, because r= 4,eAhas a size 123 of 44. We obtain the !s of theeAas described 124 in Section 2 and the ﬁrst natural frequency of 625 125 Hz, which is close to the natural frequency of the 126 system. 127 5. Deformation of Elastic Perfectly-Plastic 128 Plate with hole 129 Consider an elastic perfectly-plate with a hole, 130 one-quarter of which is shown in the Fig. 3, and 131 take the Young’s Modulus, Poisson’s ratio, and 132 yield strength as 200 GPa, 0.33, and 50 MPa, 133 respectively. We use the Abaqus ﬁnite element 134 software [1] to obtain snapshots of displacements. 135 Whentheimplicitsolverisusedtosolvestaticprob- 136 lems, time is a ﬁctitious parameter and serves as 137 a ﬁducial for the loading. However, snapshots of 138 nodal coordinates at equal intervals of time can be 139 obtained. 140 y-symmetryx-symemtry Displacement control35 mm rad.= 5 mm 25 mmFigure 3: One-quarter of a plate with hole. Symmetry and displacement boundary conditions are shown. Contours of eﬀective plastic strain are shown. The plate was meshed with linear, constant 141 strain, plane stress triangles with a ﬁne mesh 142 around the perimeter of the hole. Appropriate 143 symmetry boundary conditions were imposed and 144 the plate was subjected to displacement controlled 145 loading along the x-direction and the computation 146 was carried up to the point of plastic collapse. The 147 plate was meshed with 423 nodes, and 20 snapshots 148 were acquired leading to an Amatrix of size 846 149 846. Using r= 5, leads to a low-dimensional eA: 150 Python scripts were used to obtain a total of m 151 snapshots, doublets (x;y)in this case. We then 152 post-processed the data and obtained the state of 153 the plate for snapshot musing the sequence which 154 can be used to obtain a prediction of the state of 155 the plate for the terminal snapshot. 156 xDMD m = [x1;eUeAeUx1;eUeA2eUx1;::: eUeAm1eUx1]: (14) Accordingly, we compute a prediction error Er= 157p (xpxt)2+ (ypyt)2, which is the distance be- 158 tween the predicted location of the node, and the 159 location computed using the ﬁnite element solver. 160 The surface plot of the error is shown in Fig. 4. 161 As can be seen the prediction error is very small 162 throughout the plate indicating that the reduced 163 order model has captured the displacements accu- 164 5 4030 2010 001020×10-4 1.5 1 0.5 0 30Er(mm) x yFigure 4: Surface plot of error Erfor the elastic perfectly- plastic plate with hole. rately. 165 As snapshots are obtained, change in the state of 166 theplatecanalsobeidentiﬁedusingtheeigenvalues 167 ofeA. For snapshot i, we can obtain a correspond- 168 ingeAand its eigenvalues. As each new snapshot 169 is acquired, the eigenvalues are plotted on the com- 170 plex plane as shown in Fig. 5. With the snapshot 171 associated with the start of plasticity owing to the 172 stress concentration, the corresponding eigenvalue 173 lies outside the unit circle and subsequent eigen- 174 values will also be outside the circle. This can be 175 interpreted as a growth mode in the physical sys- 176 tem and that corresponds to the spreading plastic 177 zone. 178 6. Reduced Order Model for Large Defor- 179 mation of a Non-Linear Elastic Material 180 Consider a Mooney-Rivlin cantilever of length 250 mm and 25.4-mm cross section with strain en- ergy described by U=C10I13 +C01I23 +(1=D1) Jel12; where I1andI2are the ﬁrst and second deviatoric 181 strain invariants and Jelis the elastic volume, [1] 182 We assumed C10=25 MPa,C01=0.7 MPa,D1183 = 0.00157, and a density of 1.7 kg=m3. A total 184 of 2129 nodes and 9328 four-noded tetrahedral ele- 185 Re(λ)-1.5 -1 -0.5 0 0.5 1 1.5Im(λ) -1.5-1-0.500.511.5Figure 5: Plot of unit circle with eigenvalues of eA. When it falls out of the unit circle (arrow), it is forced to a circle of radius 1.5 for clarity. ments were used to mesh the cantilever. Tip of the 186 cantilever was displaced by (80, 80, 80) mm. 187 As before, implicit solver was used to carry 188 out the simulation, and 20 snapshots of coordi- 189 nates were output at equal time intervals. These 190 were post-processed to obtain a reduced order 191 eAwithr= 4:Using Eq. (14) the deforma- 192 tion of the terminal snapshot was predicted and 193 the prediction error was computed using Er= 194p (xpxt)2+ (ypyt)2+ (zpzt)2. Figure 6 195 plotsthecantileverinitsoriginalanddeformedcon- 196 ﬁguration, and the error contours are also shown. 197 It bears emphasis that the size of the Awas 6387 198 6387, but we are able to describe the system with 199 eAofsize88andeUofsize63878. Toemphasize 200 theimpactofincludingfewermodes, Fig. 7displays 201 the deformed conﬁguration of the cantilever when 202 r=2 andr= 8:It is apparent, that the error is 10 203 times higher and the shape diﬀers signiﬁcantly. 204 7. Summary 205 With the prevalence of DIC systems, structural 206 displacement ﬁelds are readily available and the 207 foregoing demonstrates the potential for the use of 208 6 Figure 6: Square-cross section Mooney-Rivlin cantilever. Tip is displaced by (80, 80, 80) mm. Contours of Erare shown. DMDinstructuralhealthmonitoringandsurrogate 209 models for non-linear structures. 210 8. Acknowledgments 211 Fundingfrom theNaturalSciencesand Engineer- 212 ing Research Council of Canada through the Dis- 213 covery Grant is gratefully acknowledged. 214