Crack development assessment using modal analysis in peridynamic theory

If structural damage remains undetected and is allowed to grow, structure's load-bearing capacity deteriorates, which can lead to costly repairs or in extreme cases its collapse. Modal analysis is widely used to detect structural damage because, when damage, such as cracks, is introduced, structure's geometrical and/or mechanical properties change, and these changes can be used for damage detection. Peridynamics is a non-local alternative to the continuum mechanics theory that represents forces and displacements using integral equations, which are defined even with discontinuous displacement fields, thus making this theory an attractive option for damage modeling. In this paper, authors verify peridynamic (PD) modal analysis against finite-element (FE) results, and validate it against experimental modal analysis results. The modal solver was implemented in the open-source program Peridigm and four different damage configurations were considered for verification and validation. The results show close agreement between the PD and the FE results, and the PD and the experimental results. Moreover, PD modal frequencies are shown to have similar accuracy to experimental data as the FE results. It is also shown that the frequency shifts are comparable between all three types of modal analysis. The PD mode shapes agreed well with both the FE and the experimental mode shapes at all considered damage configurations. Furthermore, the change in mode shapes from the introduced damage is similar in all three analyses.


Introduction
Structural damage is caused by design faults, construction quality shortcomings, or external effects such as improper use, overloading, natural disasters, environmental factors, etc. If the damage is not discovered and is allowed to grow, a structure's load-bearing capacity deteriorates, which can lead to costly repairs or in extreme cases its collapse. Notable examples include Ponte Morandi bridge collapse in Genoa, Italy, in 2018 (Italy bridge: Dozens feared dead in Genoa as motorway collapses, 2018); Florida International University pedestrian bridge collapse in Florida, USA, in 2018 (State: Voicemail about cracks 2 days before bridge fell, 2018); and a department store collapse in Riga, Latvia, in 2013 (Baisā traģēdija "Maxima" -lielākā kopš neatkarības atjaunošanas, 2013). Moreover, damage accumulates as a part of a structure's natural aging process. In situations as, for example, in Latvia where 37% of all bridges are reported as either in poor or very poor condition (Latvian State Roads, 2018) or in Europe where 35% of its roughly half a million rail bridges are over 100 years old (Mainline, 2013), accurate damage detection techniques can extend the life of these structures and provide measurable economic benefit (Bolognani et al., 2018;Thöns & Faber, 2014;Zonta, Glisic, & Adriaenssens, 2014).
On top of natural aging, an increasing number of extreme loading conditions due to climate change also contribute to faster deterioration of the infrastructure. In Europe, expenses of about 0.9 billion €/year are associated with extreme weather events alone and some scenarios have shown that more frequent extreme precipitations and floods in Europe in the period 2 Crack development assessment using modal analysis in peridynamic theory between 2040 and 2100 could result in an extra 50-192 million €/year cost for road transport infrastructures (Nemry & Demirel, 2012).
Modal analysis is widely used to detect structural damage, to control quality in manufacturing, to validate numerical models, etc. Modal properties -modal frequencies and mode shapes -depend on the object's mechanical and geometrical properties. When damage, such as cracks, is introduced in a structure, its geometrical and/or mechanical properties change. In principle, engineers should be able to detect introduced damage from shifts in modal frequencies and changes in mode shapes. In practice, damage detection by modal analysis is hard due to several limitations: r generally, the damage is localized, but computed or experimentally determined modes represent global displacements; therefore, local phenomena must be detected from a change in global data; r frequency shifts due to damage are often small compared to the influence of boundary conditions (BCs) or test conditions, e.g. wind for bridges; r lack of computational approaches that allow both crack growth modeling and modal analysis.
Additionally, there is a clear need for new and open computational modeling and simulation software as noted by the European Materials Modeling Council's roadmap (The European Materials Modeling Council, 2018). Moreover, the need for new modeling software particularly for engineering applications, such as uncertainty quantification, risk analysis, and decision in engineering, is recognized by the Council of the European Union in its decision that established the Horizon 2020 (H2020) framework programme for research and innovation (European Commission, 2013).
Modeling systems with evolving discontinuities, e.g. cracks, still present a challenge, because continuum mechanics theory uses partial differential equations, which are undefined when discontinuities are present in the displacement field. Peridynamic (PD) theory is a non-local reformulation of the classical continuum mechanics originally proposed by Stewart Silling to address these limitations (Silling, 2000;Silling, Epton, Weckner, Xu, & Askari, 2007). Peridynamics represents forces and displacements using integral equations, which are defined even with discontinuous displacement fields. Therefore, discontinuities are a natural part of a PD solution rather than a burden, thus making this theory an attractive option for damage modeling.
A peridynamic body in its continuum definition consists of an infinite number of nodes. Each node is connected through bonds to all other nodes within a range called the horizon. This name was chosen because interaction vanishes for nodes that are farther than one horizon away from each other. Essentially, a node does not "see" past its horizon. The peridynamic theory is non-local because a node can be connected to other nodes that are not its nearest neighbors as long as the distance between them is shorter than the horizon. Moreover, PD models converge to continuum mechanics models in the limit when the horizon approaches zero for smooth deformations (Emmrich & Weckner, 2007;Silling & Lehoucq, 2008) and non-smooth solutions of non-local models (Mengesha & Du, 2014). PD models also recover molecular dynamics results for atomistic resolution; therefore, they can be considered an upscaling of molecular dynamics models (Seleson, Parks, Gunzburger, & Lehoucq, 2009). A brief introduction to peridynamics is presented in Section 2.1, but for an extended overview, the authors recommend Bobaru, Foster, Geubelle, and Silling (2016), and Madenci and Oterkus (2014).
In this study, the authors implement a modal solver and study the well-known modal analysis in a new mechanics theory -peridynamics -in which it has not been explored in detail before. The difference from previous works is that 3D rather than 2D problems are considered, which means that more than just bending modes are present, and the results are verified and validated at several damage configurations using results from finite-element (FE) analysis and experimental modal analysis.

Peridynamic theory
A peridynamic body consists of some number of nodes that in the reference configuration are described by their position vector x i , volume V i , and density ρ i ; see a discretized 2D peridynamic body in Fig. 1. A node x i interacts with other nodes x j through bonds ξ i j where These interactions disappear behind a range called the horizon δ and nodes x j within this range are called the family of Peridynamic theory is split in its bond-based form (Silling, 2000), and its state-based form . In the bondbased form, the Poisson's ratio value is limited to 0.25 for 3D bodies and it is impossible to distinguish between deviatoric and volumetric deformations; therefore, it does not allow plastic incompressibility. The state-based form was introduced in 2007. It allows interdependencies between bonds, thus removing previously mentioned limitations. Only the state-based form is used in this study.
A body undergoes some deformation u i and nodes x i move to their deformed position y i = x i + u i . The bond ξ i j in the deformed configuration is y j − y i . A deformation u i creates a bond force density vector t i j , which depends on the collective deformation of all nodes in H i and an opposite bond force density vector t ji that depends on the collective deformation of H j is created. Since in PD each node describes some volume rather than some area as in the classical continuum mechanics, the bond forces are force densities (force per volume) and not stresses (force per area).
The bond deformation vectors for each bond in a family H i are stored in an infinite-dimensional array called the deformation state Y i : The deformation of a bond ξ i j = x j − x i in a deformation state Y i is by convention referred to using · brackets: (3) The same notation is used also for other states when referring to the value for an individual bond. Similarly, the force state T i at particle x i associates with each bond in a family H i a force density vector t i j : and The bond force density vectors are computed using bond deformations: where T (Y i ) is a state-valued function of a state -a material model. Force states are analogous to stress tensors in the classical continuum mechanics theory. Then, the peridynamic equation of motion in the integral form at some time t is where ρ i is density,ü i is acceleration, and b i is external force density.
The PD equation of motion does not contain any spatial derivatives; therefore, a solution does not require BCs, but they are commonly required to solve practical problems. In PD, BCs are applied to a layer of nodes under the surface of a body. This means that forces acting on the body's surface must also be applied to some volume under the surface as external force density b, which has the dimensions of force per volume. Similarly, other BCs such as nodal displacement or velocity must also be applied to some boundary layer. The boundary layer should be taken equal to the horizon δ (Macek & Silling, 2007).
An influence function ω i , introduced in Silling et al. (2007), is used to evaluate the contribution of a bond to the force density at a node x i . Their role has been further explored in Queiruga and Moridis (2017) and Seleson and Parks (2011). The value of an influence function can depend on different parameters, such as bond length or direction, etc. It can also be used to introduce damage by setting the influence function to 0 after reaching some damage criterion, thus removing the interaction between the nodes, i.e. breaking the bond.
One of the simplest damage criteria is the critical stretch. A bond is broken after it has stretched past some critical value s c : where s c is critical bond stretch, and s i j is current bond stretch: The natural damage growth is an advantage in PD over the classical mechanics theory. A damage model is prescribed to some nodes and, as the bonds between them break, the load is distributed to other yet unbroken bonds. The applied forces are distributed over a smaller and smaller number of bonds and, when this number becomes too low, a catastrophic failure occurs. No special elements or an a priori defined damage path is necessary for realistic damage simulations. Other damage models -such as critical energy (Foster, Silling, & Chen, 2011) and fatigue damage model (Silling & Askari, 2014) -have also been introduced.
The damage at a node can be defined as a ratio between the broken and the initial number of bonds [50], for an example please see Fig. 2: In peridynamics damage is not located only on one side of a node as shown in Fig. 2 but located on both sides of a node, thus a crack is a zone of damaged nodes similar to the cohesive zone modeling in the finite-element method.

Computational PD implementation
The modal solver was implemented in open-source peridynamics code Peridigm (Littlewood, Parks, Mitchell, & Silling, 2013;Parks, Littlewood, Mitchell, & Silling, 2012), which was chosen because it provides an easily extensible, robust framework for peridynamic computations. Moreover, as Peridigm already implements quasi-static and explicit solvers, several parts -material models, stiffness matrix creation routines, discretization management, and BCs -could be reused. Peridigm is built using packages from the Trilinos library (Sala, Heroux, Day, & Willenbring, 2010). The modal solver was created using Belos (Bavier, Hoemmen, Rajamanickam, & Thornquist, 2012), Ifpack (Sala, 2005), and Anasazi (Baker, Hetmaniuk, Lehoucq, & Thornquist, 2009) packages. The tangent stiffness matrix in Peridigm is built using a central difference scheme and the exact algorithm is presented in Parks et al. (2012). Peridynamic stiffness matrices are often nonsymmetrical, because even though two nodes x i and x j are connected by the same in magnitude but opposite bonds x j − x i and x i − x j , the force densities in them are influenced by different neighborhoods H i and H j . Due to the non-symmetrical stiffness matrix, solvers for non-Hermitian problems must be used.
The mass matrix describes the mass that acts on each separate degree of freedom. In PD, a node describes some amount of volume around it; therefore, the lumped mass approach is an intuitive way to build the mass matrix: where M ii is an entry in a mass matrix M, V i is the volume described by a node x i , and ρ i is the density of a node x i . The peridynamic mass matrix is built by looping over all nodes, multiplying their volumes with their densities, and storing them in a separate matrix. This matrix is positive and Hermitian. The shift-invert technique was used to speed up the modal analysis. A generalized eigenvalue problem where K is stiffness matrix, M is mass matrix, x is eigenvector, and λ is eigenvalue, is transformed into: where C is a substitution and θ is a substitution and σ is a shift value. Then, (12) is a simple eigenvalue problem that was solved using block Krylov-Schur eigensolver (Stewart, 2002), but the substitution in (13) was solved using flexible block GMRES linear solver (Saad, 1993) with ILU preconditioner (Saad, 2003). The linear peridynamic solid (LPS) material model was used and the material properties were r elastic modulus -5.35 GPa, r Poisson's ratio -0.33, r density -1200.00 kg/m 3 .
No BCs were applied, thus simulating free-free test conditions. Damage in a model was created by specifying a crack plane and breaking any bond that crosses this plane. Damage configurations and locations are presented in Section 2.4.
In PD, three types of convergence can be defined. They will not be discussed in detail here, but they were introduced in Bobaru et al. (2009). The δm-convergence, where δ → 0 and the mesh density m increases as δ decreases, was studied to find an acceptable discretization. Four different model sizes, 40 000, 135 000, 320 000, and 625 000 nodes, and two different horizon lengths, δ = 0.02 m and δ = 0.0016 m, were used. Figure 4 shows convergence for the 1 st bending mode and the 1 st torsional mode. Both plots show that the computed modal frequencies converge to a single value for a single horizon length as the mesh gets refined from 40 000 to 625 000 nodes. Also, the asymptotic convergence range is reached when the model size is 625 000 nodes. The rest of the presented results, therefore, were computed and are presented using the 625 000 node model.
PD simulations were run on the computing cluster Vasara at Riga Technical University. Different numbers of Dell EMC Pow-erEdge R640 nodes each with 2 x Intel(R) Xeon(R) Gold 6154 3.00 GHz CPUs with 36 cores were used. Each node had 348 GB maximum available RAM and a 240 GB SSD. They were connected through an InfiniBand EDR 100 Gb/s connection. The cluster ran Centos 7.5 operating system and used the Torque 6.1.1.1 resource manager.

FE implementation
The FE analysis results were used as the benchmark because it is the most widely used computational method for modal analysis. It is well understood and its results are trusted by both researchers and engineers. Therefore, it was a natural choice for the verification of PD results. A rectangular plate model with the same dimensions as the PD model -0.10 × 0.05 × 0.008 m -was created using Ansys FE software. An elastic material model was used, and the material properties were the same as for the PD simulations. The FE model was meshed using SOLID185 8-node cubic elements and contained 664 146 nodes.
This kind of problem could have been modeled using shell elements, which would have been computationally less expensive; however, solid elements were chosen for two reasons: r the authors wanted FE and PD meshes to have a similar number of nodes, and r it was possible to create cracks through the depth of the FE model, which would not have been possible if the model was created using shell elements.
Cracks in the FE model were created by not connecting the nodes of solid elements on a crack plane. Such an approach is not ideal; however, it is good enough for verification purposes because as will be shown in Section 3 the FE results agree with the experimental results well and are, therefore, reliable. Damage configurations and locations are presented in Section 2.4.
The Block Lanczos solver was used to solve the underlying eigenvalue problem. Nodal mass was approximated us-ing the lumped mass method, in which mass matrix is computed assuming that the mass at a node is the mass of the object's part that is closer to that node than any other node.

Test specimens and experimental setup
Plate-shaped specimens with dimensions of 0.1 × 0.05 × 0.008 m were manufactured from several large polymethyl methacrylate (PMMA) sheets using laser cutting.
Five different damage configurations were considered -one with no damage called Healthy configuration and four configurations with increasing damage called configurations Tb through Td. Moreover, five specimens were tested at each damage configuration and only the average values are presented hereafter. The damage was introduced by cutting a crack through the thickness of a specimen using a laser. Since a laser melts the material under the beam, a piece of paper was slid through the crack to ensure that it is clear from any leftover material. The crack nominal length -which was used in the simulations -was 12.5 and 25.00 mm. Their location and actual lengths are shown in Fig. 5.
The experimental modal analysis test setup is shown in Fig. 6. Test specimens were suspended from an aluminum frame in two cotton thread loops to ensure free-free BCs and they were excited using a loudspeaker.
2D Polytec PSV-400 laser vibrometer measured specimen speed at several hundred points on their surface. Measurement range was between 1000 and 8000 Hz with measurement step f = 2.5 Hz. Each point was measured for 400 ms and velocity decoder VD-07 1 mm/s/V was used.
Since PMMA is translucent, a single layer of paper tape was applied to prevent the laser beam from shining through. Tape increased specimen mass and thickness by less than 1%, so the effect on measured modal properties can be considered negligible. Shallow cuts were made at the bottom of the specimens at 25 mm distance from either end, to have the two cotton loops always positioned at the same place to ensure that BCs between all tests were as similar as possible.

Verification and validation of PD modal frequencies
The first 12 modes were computed using PD and FE modal analyses. The first six were rigid-body-motion modes, but the other six can be used for verification.
Not all modes were computed at all damage configurations. Only the 1 st bending, the 2 nd bending, the 1 st torsional, and the 2 nd torsional mode were computed with PD and FE analysis at all five configurations. Other modes, such as the 1 st transversal bending mode, were computed only at Healthy and Tb configurations, but not at the other ones. This happened because the mode order changed as the crack length increased at different damage configurations.
The % difference between the PD modal frequencies and the FE modal frequencies at Healthy and damage configurations are presented in Fig. 7, but exact modal frequency values are included in the supplementary data.
PD results with shorter horizon δ = √ 2 h agreed better with the FE results than simulations with longer horizon δ = 4h. Such behavior was expected because the FE analysis gives a local solution, which should be better approximated by a shorter horizon, i.e. δ = √ 2 h. Moreover, it has been shown that under sufficiently smooth motion, constitutive model, and non-homogeneities, the non-local PD solution approaches the local solution when the horizon shrinks to zero (Silling & Lehoucq, 2008).
Some differences between the results were expected because the same problem was solved using two different mechanics theories, but the results show a very good agreement. The largest difference between PD and FE modal frequencies when δ = √ 2 h at Healthy configuration is −0.51% for the 1 st bending mode, at Tb damage configuration −0.66% for the 2 nd bending mode, at Tc damage configuration −1.57% for the 1 st bending mode, at Td damage configuration −1.45% for the 2 nd torsional mode, and at Te damage configuration −2.37% for the 2 nd torsional mode. In simulations with δ = 4h, the largest difference at Healthy configuration is −2.02%, at Tb configuration −2.10%, at Tc configuration −2.91% all for the 1 st bending mode, at Td damage configuration −2.95% for the 2 nd torsional mode, and at Te damage configuration −3.92% for the 2 nd torsional mode. PD modal frequencies were lower than FE modal frequencies at all damage configurations. The PD and the FE mass matrices were created using the same lumped-mass approach, material density was the same and the total volume of the model was the same. The biggest difference between the matrices was that the FE mass matrix was larger because there were more nodes in the FE model. Therefore, the differences in computed modal frequencies most likely occurred due to different stiffness matrices. Generally, the PD stiffness matrix is denser than the FE stiffness matrix for an equivalent problem; also its bandwidth is wider. This happens because a PD node is connected to all other nodes within its neighborhood rather than only its nearest neighbors as is the case for a FE node. Since the PD modal frequencies were lower, the PD computational stiffness must have been lower than the FE computational stiffness.
Moreover, PD simulations with the horizon δ = 4h resulted in lower modal frequencies than PD simulations with δ = √ 2 h. This shows that the computational stiffness is influenced by the length of the horizon and that it is lower when longer horizon values are used.
The reduced stiffness in models with a longer horizon is caused by the so-called peridynamic "surface effect." PD material models are derived for a node in the bulk of a model, i.e. it is assumed that the neighborhood of a node is full of other nodes. When a node is located near an edge of a model or a material surface, its neighborhood is "missing" some nodes and therefore in ordinary state-based PD theory the model may be stiffer or softer on the surface than in bulk along with being slightly anisotropic even when trying to model isotropic material, under non-homogeneous deformations (Le & Bobaru, 2017). In a 2D peridynamic object, shown in Fig. 8, a node with δ = 3h (Fig. 8a) next to a material surface has its neighborhood filled by 60.71%, but, if its horizon is reduced to √ 2h (Fig. 8b), then its neighborhood is filled 62.50%. Since the material model was derived under the assumption of 100% full neighborhood, the stiffness of nodes with a longer horizon is reduced more than the stiffness of nodes with a shorter horizon.
Several studies have tried to remedy this effect to greater or lesser success. The most well-known methods include: the volume method , the force density method Oterkus, 2008), the energy method Oterkus, 2008), the force normalization method (Macek & Silling, 2007), position-aware linear solid constitutive model (Mitchell, Silling, & Littlewood, 2015), and the fictitious nodes method (Gerstle, Sau, & Silling, 2005;Oterkus, Madenci, & Agwai, 2014). For a comparison between different methods, please see Le and Bobaru (2017). Most of the methods mentioned previously are adhoc solutions and their effect depends on the model's geometry, applied BCs, etc. Recently, a promising method that eliminates the need for ad-hoc solutions has been presented in Madenci, Dorduncu, Barut, and Phan (2018).
Five modes were measured at Healthy and Tb damage configurations, but six modes were measured at Tc, Td, and Te configurations. More modes were measured for these three configurations because as damage increased an additional mode shifted in the measured frequency range.
The % difference between PD, experimental, and FE modal frequencies at Healthy and T damage configurations are presented in Fig. 9. PD results at Healthy configuration were within a range of −1.31% to 2.42%, at Tb configuration within −1.68% to 2.59%, at Tc configuration within −1.44% to 2.55%, at Td configuration within −3.20% to 1.65%, and at Te damage configuration within −2.84% to 1.98% of the experimental modal frequencies.
The differences between the PD results and the experimental results are low, which shows that PD modal analysis closely reflects the modal behavior of the test specimens.
PD results are as similarly close to the experimental results as the FE results. All PD modal frequencies were within −3.20% to 2.59% from the experimental modal frequencies, but the FE results were in a range from −1.20% to 3.57%. These ranges are reasonably close, therefore, PD modal analysis is at least just as accurate as FE modal analysis for these test conditions. When separate damage configurations and modes are considered, then the FE results are closer to the experimental results than the PD results only for the 2 nd bending mode at Td and Te damage configurations (Fig. 9d).
Computational models are an idealized representation of test specimens as they do not contain surface cracks, microvoids, and other imperfections. So, the authors expected that their stiffness would be higher, hence, computed frequencies should be higher. This was not true for all PD frequencies. The computational mass and the specimens' mass should have been almost equal because the computational model and test specimens had roughly the same size (nominal rather than real dimensions were used for computational models) and the density used in simulations was calculated from the specimen mass. Therefore, the computational stiffness in some cases must have been lower than the real specimen stiffness. Moreover, the difference between the experimental results and the PD results for the same mode but at different damage configurations was not constant, which means the stiffness reduction in computational models and the specimens must have been different. Most likely it is related to the differences in crack length. Computational models used the nominal crack length, which was different from the real crack length in the test specimens.

Frequency shift
The frequency shift is the change in modal frequencies when an object's stiffness changes due to the introduction of damage or some other influence. The shift occurs, because an object's material properties or geometric properties have changed. In the presented cases material properties stayed the same, but the introduced cracks changed the geometrical properties of the computational models or the test specimens.
Computed modal frequencies varied between 1000 and 7000 Hz. Their absolute values and the difference between them would depend on the mode in question, therefore, relative frequency shifts will be presented. Modal frequencies were normalized against their respective modal frequency at Healthy configuration and the difference between the normalized modal frequencies at Healthy configuration and at the other damage con- where ν δ= √ 2 h Tc, n is the relative modal frequency shift at Tc damage configuration, ν δ= √ 2 h Tc is the modal frequency computed with δ = √ 2 h for Tc damage configuration, and ν δ= √ 2 h H is the modal frequency computed with δ = √ 2 h for Healthy configuration. Relative frequency shifts at Healthy and T damage configurations are presented in Fig. 10, but exact shift values are included in the supplementary data. The PD results show a clear frequency shift. At the most severe -Te -damage configuration the largest shift is for the 2 nd torsional mode −40.46%, but the smallest shift is for the 1 st bending mode −25.61%. The results show a large 14.85 percentage points difference, which means that crack influence on the stiffness of the PD models differs considerably depending on the mode in question.
The length of the horizon had a minuscule influence on the PD frequency shift. For the 1 st bending mode, the frequency shift was −25.73% and −25.61%, for the 1 st torsional mode −26.50% and −26.55%, for the 2 nd torsional mode −40.39% and −40.46%, for the 2 nd bending mode −32.38% and −32.33% when δ = √ 2 h and δ = 4h respectively.
The agreement between the maximum PD and the FE frequency shift and between PD and experimental frequency shifts was excellent. The smallest and the largest differences between PD and FE frequency shift were −1.38 and −0.04 percentage points for the 2 nd torsional and the 1 st bending mode. The difference between the PD and the experimental relative fre-quency shifts ranged between −1.47 (the 2 nd bending mode at Td damage configuration) and +0.34 (the 1 st bending mode at Tb damage configuration) percentage points. Relative frequency shift differences within such a small range show that PD frequency shifts closely mimic the FE and the experimental results.

Verification and validation of PD mode shapes
Mode shapes for Healthy, Tc, and Te damage configurations will be presented. Authors chose not to present mode shapes at all damage configurations because little useful information would be added and they would take up too much space. These three damage configurations were chosen because they represent the undamaged, the most damaged, and one damage configuration between them.
Three different computer programs -Paraview, Ansys, and Polytec Scan Viewer -were used to visualize the PD, FE, and experimental mode shapes. The programs used similar, but not completely identical, color schemes. The color distribution in the mode shapes, therefore, can appear slightly different. Figures 12-14 present the mode shapes as a sum of displacements in X, Y, and Z directions for a particular shape. Red color shows the extreme displacements, but blue color shows displacements close to zero; the whole color scheme is shown in Fig. 11. A color legend with numerical values cannot be provided, because a mode shape is a dimensionless representation of a structure vibrating at a modal frequency.     Figures 12-14 show the PD and the FE mode shapes at Healthy, Tc, and Te configurations. PD modal analysis and FE modal analysis yielded the same mode shapes at each single damage configuration. At Healthy configuration the mode order was: the 1 st bending, the 1 st torsional, the 2 nd torsional, the 2 nd bending, the 1 st transversal bending, and the 1 st in-plane bending mode. The mode shape displacement field is regular and the zero displacement areas are roughly parallel or orthogonal to one of the three coordinate axes. At Tc damage configuration, Fig. 13, the mode order was: the 1 st bending, the 1 st torsional, the 1 st in-plane bending, the 2 nd torsional, the 2 nd bending, the 3 rd torsional mode. At Te damage configuration the mode order was: the 1 st bending, the 1 st torsional, the 2 nd torsional, the 1 st in-plane bending, the 2 nd bending, the 2 nd in-plane bending mode.
The introduced cracks changed the model's stiffness, but stiffness acting against different kinds of deformations was changed differently. The in-plane bending stiffness was reduced more than other stiffness, therefore, in-plane bending mode's modal frequencies experienced a larger frequency shift and the 1 st in-plane bending mode shifted from being the sixth mode at Healthy configuration to the fourth position at Te configuration. Similarly, the 2 nd in-plane bending mode, which was not between the six simulated modes at Healthy damage configuration, appeared at Te damage configuration. Therefore, the mode order is different between Healthy and Te damage configurations. But the order is the same at every single damage configu- ration, e.g. the same mode order at Te damage configuration in both PD and FE results.
The introduced cracks created a sudden discontinuity in the displacement field and the displacement field shifted in response. As an example, please see Figs 12b, 13b, and 14b. The maximum displacements of the 1 st torsional mode at Healthy configuration ( Fig. 12b) are located at the four corners of the model. When the first crack is created (Tc damage configuration in Fig. 13b) the maximum displacements shift to the top right and the bottom left corner and to the right side from the created crack. Lastly, at Te configuration ( Fig. 14b) maximum displacements again shift to the two top corners, and also to the whole bottom side to the left from the left crack and to the right from the right crack. All mode shapes show such shifts in the displacement field and in all considered cases they were the same in PD and FE analysis results.
Only 2D modes could be obtained with the experimental analysis because a 2D laser vibrometer, which cannot measure in the in-plane direction, was used. The authors chose to also present the PD modes in 2D so that they can be compared with the experimental results more easily. Furthermore, due to this limitation, only four modes -the 1 st and the 2 nd bending and the 1 st and the 2 nd torsional mode -were both computed and measured at all damage configurations. Therefore, these four will be compared.
The color scheme used in Figs 16-18 represents displacement in the out-of-plane direction. Therefore, the color scheme used here is different. Numerical values cannot be used for the reasons discussed previously. The color scheme is presented in Fig. 15. Peridynamic and experimental mode shapes in Fig. 16b and c and Fig. 17a-c are out of phase by π. The authors are unable to draw all mode shapes in the same phase due to limitation of the software used for illustration.   The PD and the experimental mode shapes are similar for all considered modes and damage configurations. The created cracks can be seen clearly in all PD mode shapes, but in the experimentally measured mode shapes the created cracks can be seen clearly only in the torsional mode shapes, but not in either of the bending mode shapes. This was expected because, in the bending mode shapes, the crack location coincided with the location of the extreme displacements. The crack would not appear as a discontinuity in the displacement field because the measured displacements on both sides of a crack would have been roughly equal.
The mode shapes shifted when damage was introduced. For example, please see Figs 16a, 17a, and 18a, in which the 1 st bending mode at Healthy, Tc, and Te damage configurations are shown. Figure 16a shows the 1 st bending mode of an undamaged specimen, the extreme displacements are in the middle of a specimen and at both ends. When a single crack is introduced, the mode shape shifts clockwise (see Fig. 17a), and it is no longer symmetric. Now the extreme displacements are located in the top left corner, the bottom right corner, and in the center of a specimen. Lastly, in Fig. 18a, when two cracks are introduced, the extreme displacement is again located in the middle and along the ends of the specimen. Since cracks are located on the bottom side, the extreme displacements have shifted slightly downwards. These changes in PD mode shapes were the same as in the experimental mode shapes. This shows that the modal analysis of a PD model with damage well describes the modal properties of a real object with similar damage and that the change in modal properties due to cracks is comparable in both PD simulations and real measurements.

Conclusions
This study considered a novel approach to modal analysis in the peridynamic theory. A peridynamic modal solver was implemented in an open-source peridynamics code Peridigm. It was used to verify and validate peridynamic modal analysis against finite-element and experimental modal analysis results. Four different damage configurations with increasingly damaged cross-section were considered.
The agreement between the peridynamic results and the FE results is excellent at all damage configurations. The differences in the computed modal frequencies ranged between 0.00% and −4.00%. PD modal frequencies were lower than FE frequencies at all damage configurations, which means that PD computational stiffness was lower. As expected, PD modal frequencies agreed with the FE modal frequencies better when the horizon was shorter because the analysis with a shorter horizon is "more local." Additionally, peridynamic modal frequencies are lower in simulations with a longer horizon. This is explained by the increasing peridynamic "surface effect." PD results also agree well with the experimental results. The PD modal frequencies were within ± 3.2% of the experimental results, while the FE results were within a range of −1.20% to 3.57% of the experimental frequencies. This shows that PD and FE modal analysis results are similarly accurate when compared to the experimental modal analysis results.
Moreover, the PD frequency shift is similar to the frequency shifts in FE and experimental analyses. The largest difference between the PD and the FE frequency shifts was −1.38 percentage points but between the PD and the experimental shifts it was −1.47 percentage points.
The PD and the FE mode shapes agree well and are in the same order at each damage configuration. The agreement between the PD and the experimental mode shapes is good; furthermore, the change in the mode shapes from the introduced damage is similar in both analyses.
Potential coupling of peridynamic damage simulations and modal analysis for damage detection remains an area for future research.

Supplementary Data
Supplementary data are available at JCDENG online.