Remote vibrometry recognition of nonlinear eigen-states for object coverage of randomly large size

For objects of"large"vibration size such as waves on the sea surface, the choice of measurement method can create different understandings of system behavior. In one case, laser vibrometry measurements of a vibrating bar in a controlled laboratory setting, variation in probe spot size can omit or uncover crucial structural vibration mode coupling data. In another case, a finite element simulation of laser vibrometry measures a nonlinearly clattering armor plate system of a ground vehicle. The simulation shows that sensing the system dynamics simultaneously over the entire structure reveals more vibration data than point measurements using a small diameter laser beam spot, regardless of the variation of footprint (coverage) boundaries. Furthermore, a simulation method described herein allows calculation of transition probabilities between modes (change-of-state). Wideband results of both cases demonstrate the 1/$f$ trend explained within -- that the energy of discrete structural vibration modes tends to decrease with increasing mode number (and frequency), and why. These results quantify the use of less expensive non-imaging classification systems for vehicle identification using the remote sensing of surface vibrations \emph{while mitigating spectral response distortion due to coverage variation on the order of the structural wavelength} (spectral reduction or elimination).


I. INTRODUCTION
Remote vibrometry has many forms. The use of optics to sense vibration is a mature technology dates back to methods of observations with knife-edge imaging invented by Jean Bernard Léon Foucault in 1858 that can be used to measure mirror shape. The method was also used by Töpler [1] in 1866 to measure phase and later used to remove unwanted optical phase effects: ". . . the Schlieren method, where all the spectra on one side of the central order are excluded." [2] The Schlieren method, a German word for streak, uses Foucault's knife-edge in the "Fourier plane" at a distance of a focal length from a collecting lens. A related method, the streak tube, can optically process an entire slice of underwater scenery nearly instantaneously [3], [4], [5]. More recent physics discoveries use vibration to measure atomic size effects. Microscopy resolution smaller than the optical diffraction limit can be accomplished with the Nobel prize winning [6]) scanning tunneling microscope. Other microscopes image S-shells of large atoms and, in some atoms, their P-shells. The microscopic range, remote vibrometry variant, the atomic force microscope (AFM) is able to detect separate molecules by measuring the difference in van der Waals force as a tiny vibrating cantilever sensitive to the atomic outer shell locations, scans over the sample. These are some examples of remote vibrometry. This article describes how the nature of observations can provide an different pictures of the vibrational modes of common vehicles and other integrated structures built from plates, shells, and membranes (the "body-in-white").

SUPPORTING RESULTS:
Three types of analysis provide support for this investigation of observation methods for structural surface waves. (i) Experimental results of a driven clamped bar with full fixity at both ends provides a laser vibrometry spectrum for both small and large probe beam spot sizes. (ii) The author's master's thesis [7] contains calculation results of imaging a vehicle hull at 4 kilometers. Identification of the target vehicles based on spectral fingerprints formed from laser vibrometry. The target armor clattering against the vehicle hull provides a surface vibration that modulates the probe beam. (iii) Several analytical calculations using explicitly nonlinear mechanics validate the distinct clattering spectrum. These experimental, simulation, and analytical results provide cross-supporting reinforcement for the three observational characteristics described in this work.   1. The lowest frequency modes have high resonance quality Q compared to high frequency response, unless systems have large damping on these "lower" modes or frequency-dependent control systems. For adequate signal to noise ratios, these frequencies identify what are effectively discrete modes, as represented in the stem chart to the right.

OBSERVATIONS:
This work defines three characteristics of composite structures such as vehicles. (1) The lower frequency modes are nearly discrete. 1 (2) Energy from driven modes flows into other modes. And finally, (3) due to theory and observations defined herein, vibration strain-energy is higher in the fundamental mode, with energy decreasing as modal frequencies increase. The latter two characteristics can be considered [8] to be consequences [9] of the second law of thermodynamics [10].
The probe beam spot size affects vibration modes that are observed with optics. Some vibration modes may not be obvious from return due to a probe beam with a small spot size, as discussed with Flight Lieutenant Ngoya Pepela's thesis [11] in Figures 5 and 6 for a vibrating bar. In the simulation of clattering armor plates, more mode information is acquired when viewing the entire optical image and processing its spectrum at multiple locations. In spite of some interference, the non-imaging 2 (spatially integrated radiant flux) spectrum of the probe return from the entire target's surface provides more information than the small spot size in that phase information over time is embedded in the sensed signal [7]. On-average spectral elimination (or spectral reduction). 3 [End of the table of contents. This will be removed in the published paper.]

A. Three observational characteristics
Based on the three previously mentioned supporting results, the experiments, simulation, and analytical results, three observations become apparent: (1) the lower frequency modes have sharp resonances that are effectively discrete responses. Figure 1 is a sketch of a notional representation of this case. These modes can be modeled on center frequencies when damping is limited to structural damping [12], as later discussed in Figures 5 and 6. (2) energy transfers from one structural component and its vibration modes into the lower modes of nearby components through joints, including welds, which have nonlinear load-deflection curves. Each mode is comprised of a 3-D mode shape (deflection vectors) of the entire structure ( Figure 2). The strain-energy (modal energy) of each mode couples to lower energy modes as the system transforms from transient to steady state. Finally, (3) modal energy flows into the lowest frequencies first, saturating them ( Figure 3), depleting the energy flow as the mode number increases from f 0 , f 1 , f 2 , · · · and up through the closely-spaced high frequency modes that become practically infinitely dense. For discrete modes the finite element analysis (FEA) results are limited to the number of degrees of freedom (DOF) in the model.
The energy at each of the modes shown in the notional sketch ( Figure 1) transfers to lower energy modes if physical load paths exist that are conducive to mode coupling. These sketches summarize decades of trade secret testing. At lease one example appears in the literature for automotive modal analysis, ride and handling [13]. Helicopter data is available from MIL-HDBK-810 [14]. A transfer of energy from a strongly driven mode to other modes appears in Figure 2. That driver might be in a component located far from the receiving modes [15], but in 'frequency space' these modes might be coupled. Vibrations might be driven by the engine at one frequency, yet sensed in the back seat at another frequency. Nonlinearities in response can create harmonics of resonances that allow energy to couple with lower modes. The low energy tails of wide resonances (damped broadening) also "touch" other modes of lower frequencies. Another energy transfer mechanism occurs in clattering plates where spatial distribution of mode amplitudes causes contact that excites different modes. This transfer of 1 In practice, modal engineers see effectively discrete modes for the fundamental and other low frequency modes. Here fundamental means the lowest frequency mode, rather than the music definition as the largest common integer-based factor. For 300 Hz beat with 500 Hz our modal fundamental mode would be at 300 Hz, not 100 Hz. The experience with effectively discrete modes occurs within trade secret production engineering analysis tasks such as automotive ride and handling -one of the most guarded secrets. These frequencies and mode shapes are rarely published in the literature. 2 Imaging systems require more hardware to calculate pixel values. Non-imaging systems integrate power incident on the entire optical aperture. This power is termed the radiant flux, measured in watts. It is the flux of the Poyting vector S through an enclosing surface Φ E = S ·ndA. 3 The laser vibrometry industry coins the term spectral elimination to represent reductions to represent destruction or distortion of spectral ID features. Other industries use the same term for entirely different concepts. Fig. 2. Vibrational strain-energy from different components transfers from a driven mode, the indicated most powerful response, to lower energy modes at different frequencies. For example, the nonlinearity of all common, practical fasteners (bolted and riveted joints) allows energy transfer. In this notional case, a high-Q response at the driving frequency bleeds energy into other modes through friction or contact, a nonlinear process. If the response at the driving frequency, fdriven, has a wider full width at half height, then energy would transfer faster for nonlinear systems. Some energy would then also transfer for linear systems as well, more so modes with a smaller frequency difference ∆f = |f i − fdriven|. As described by the Zienkiewicz quote in Section III, high frequency modes eschew energy that then tends to flow into lower frequency modes. These lower modes tend to "fill up" first and they retain more energy than higher frequency modes. There are at least two exceptions: the lack of vibration coherence can block the energy transfer, and interfaces that act like active systems can manipulate the energy or act like passive structural filters. energy is much similar to plucking a musical string at different locations to change the overall sound. Energy transfer between states depends on the frequency of the mode, which is the square root of the eigenvalue in an analysis of normal modes. Figure 3 is a sketch where the final steady state of the vibrational system has energy ordered from strong resonances at the fundamental frequency to a paucity of energy at higher frequencies. A quote from Zienkiewicz, displayed later in Section III, provides an explanation for this common 1/f dependence of modal energy. The exceptions noted also include strong damping or added mass. 4 The models, structures, and observational characteristics that comprise this article are summarized in the Conclusion in Table  VI. B. Physical model: symmetries, non-imaging, and coverage By using a "pencil-thin" beam, conventional laser vibrometry systems can sometimes better identify modes in a simple academic structure such as the solitary vibrating bar in the following photos. However, a large spot size can adequately measure those modes and more. This can be shown by simulation and analytical calculation. The simulation uses FEA for vibration deformation data that is handed-off to a set of MATLAB™functions that perform Fresnel propagation for the optical sensing. The FEA shows the mode shapes and calculates the energy per mode which are the modal participation factors (MPFs) [16], [17]. The FEA surface vibration results are input for the MATLAB image propagation calculation (Fresnel propagation Fig. 4. Coherent radiance undergoes a spatially harmonic phase modulation as a function of the probe beam location due to target vehicle surface skin structural vibration. The Object modulates radiance from the outbound probe beam at a distance of up to several kilometers while the receiver system keeps track of the phase of this exitance using the reference beam from articulating beam-splitter A and return mirror B, and appropriate delays (delay fiber or other systems not shown) where appropriate. [18]). The physical schematic for the simulation appears in Figure 4. Laser vibrometers use many different methods that can have an effect on the precision of the measurements. 5 The laboratory measurements of the vibrating bar show the signal to noise ratio (SNR) is quite high for structural components and, unless the system is purposely damped, the resonances are sharp. These sharp resonances lead to the discrete mode picture (stem plot at right in Figure 1) that the field of modal engineering evolved to exploit. The literature provides the theoretical reasons that strain-energy tends to eschew higher frequency modes and thus populate the lower modes. Years of experience of vehicle design and analysis support these conclusions, as does the FEA of clattering armor plate reported here and in detail in a thesis [7]. The FEA results show that unsymmetrical modes 6 of the clattering plates have higher frequencies than symmetrical (in-phase) vibration of the parallel plates, and that these low frequency modes are discrete below 500 Hz, as is the case with typical automobiles. These modes are usually easy for laser vibrometers to detect. The optical image propagation model shows that even a non-imaging spectrum of the surface of the clattering plates contains energy-ordered modes whose vibration strain-energy drops off with frequency, as described in the Zienkiewicz quotation provided here in Section III.
Commercial laser vibrometers tend to use pencil-thin spot sizes on the object being measured. This spot's cross-section can be as small as a one pixel response. Large spot size imaging systems have far more information, much of which can help spectral identification (ID), such as the method developed herein. Further, non-imaging signals average responses from all pixels for each time step. An optimum automated target recognition (ATR) system could use a combination of small spot size and full size illumination. But a large spot size tends to help tracking in the paint-the-target phase, while still enhancing spectral ID.

II. STRUCTURAL NONLINEARITIES
Several different types of laser vibrometers can measure dynamical properties of vibrating surfaces. The type of laser vibrometer is less important than the phenomenology of what can be measured remotely. A priori project restraints such as requirements to calculate a transfer function, or to use pencil-thin probe beams, can limit the sensed spectra and thus omit or change perceived vehicle behavior(s) in a manner unrelated to errors in measurement or analysis. In many fields such as acoustics, the modal results can be less accurate than a kinematic approach would provide. Nevertheless, especially for complicated nonlinear systems, it sometimes pays to ignore detailed kinematic analysis, which are often invalid for nonlinear systems, in favor of modal engineering to tabulate the modes (deformed shapes for particular resonance frequencies), their frequencies, and their energy or modal participation factors [17].
For over half a century modal engineers routinely used modal analysis for nonlinear vehicle structures. This work uses a structural model of a plate bolted to another plate on a base structure. The work also uses simplified and even "linear" component studies. These studies cannot deliver the behavioral metrics that the full system model provides. Component studies do verify limited aspects of the full nonlinear model, which is a crucial part of the overall analysis. In the full analysis, the structure is a free-floating system using D'Alembert reactions, an analysis technique common to vehicle engineering. See, for example, the numerous 'quarter-panel models' in the literature [19], [20], which are mostly meant to show how to model some features for reasons stated previously. 7 The sensor model is an optical model without the usual air turbulence analysis, except to check for Fante's wavefront coherence breakup range [21] and a few other imaging issues [7] that are out of the scope of this paper.
The physical model for the clattering armor assumes a linear superposition of vibration modes on a nonlinear structure. This is not a linear time-invariant (LTI) system. The details of removal of these typical LTI assumptions and its relationship to stability, stabilization, and stabilizability, are in the thesis [7], available through DTIC dot mil. The analysis includes a treatment of how changes in fixity causes changes in the modes. Similar changes can model the aging of parts. Therefore, the modulation of the sensed optical radiant flux is related to the analytical nonlinear expressions in the Appendix that help explain the contact nonlinearity in the FEA simulation. The modeling complications related to fixity, nonlinearities, and modal participation factors [16] are dependent on control of the nonlinear solution per time step and over time using load increments, iterations, and other parameters. The analyst also needs to make sure the solution follows the millions of load deflection curves [22], [7] in the full system model.
The finite element analysis (FEA) output represents the three-dimensional structural model of the clattering armor plate system and provides the time history of a structural vibration that modulates the optical return recorded by the sensor. The Appendix contains a set of single and two degree of freedom analyses that result in analytical expressions. The objective of the analyses is to use the optical images, or non-imaging radiant flux (watts), to determine the frequency of structural spectral modes for target identification (ID).

A. Modal Participation Factors
Recent conversations with Navy scientists and engineers tend to evolve into studies of the role that the nature of system observation plays in the perception of the modal system, and how the low frequency strain-energy modes of the structure are effectively quantized. The discrete nature of the modes has been an underlying assumption of modal analysis and in its application to ride-and-handling for half a century. The modal participation factors [17] describe the dynamics of energy flow by using vectors of MPF's that span the pertinent modes [16]. Each MPF is an element in vectors Φ[f i ] that are arrays ordered by mode number i, used in the FEA of vibrating structures such as vehicles, antennas, and spacecraft. For example, NASTRAN™may require 'DMAP alters' (macros) to read out some of the components of Φ[f i ].
An example of the unintended consequences of design decisions for the 1990 era Corvette follows this detailed restatement of the three modal characteristics introduced earlier, which are essential to the correct system ID. They will be used in the example. (1) The lower frequency modal states (eigenvectors) form discrete modes, even for the complicated system synthesis models and the vehicles they represent. The work of vehicle design for these issues focuses on (2) transition probabilities 8 and (3) the flow of strain-energy that preferentially fills the lower modes with energy, as discussed in the quoted passage in the next section.
B. Observational characteristics: Discrete, Transition, and Ordering 1. Discrete modes are countable sets of distinct and high signal to noise ratio (SNR) response of bandwidth narrow enough to be obviously a single mode. The strain-energy levels within each modal state are tabulated by the structural design team in order to determine critical components. These critical regions absorb most of the ride-and-handling and structural engineering effort during vehicle design, resulting in a final list of lower frequency energy levels 9 The lower modes relate to modes shapes that are usually characteristic of the structure, such as bending, twisting, and stretching. A real-life example appears later in this section.
2. Transitions between modal states, from one eigenvector of the structure to another, are often sigmoid in nature 10 where the region of rapid mode coupling depends on meeting strain-energy thresholds that allow more efficient coupling between modes. Most of the energy tends to flow to components with lower fundamental modes according to characteristic (3) below. Nearly 7 As of 21 November 2019 only the first two hits for "quarter-panel dynamic response" were related to vehicles. No experimental results showed up in the top (most pertinent) 10 pages of 10 hits each of 1,050 hits. 8 Transition probabilities are related to structural coherence spectra. Modal engineers use them to reduce undesired vibration modes (resonances) in the structure. [23], [24]. 9 Acceleration spectra are measured in units-per-hertz such as strain (strain-energy amplitudes) or acceleration. In practice the "power" is displayed, signal processing "power" being the square, element-by-element, in g 2 /Hz. However, the optically sensed units are dB/Hz, which for non-imaging is re radiant flux Φ e,ref in W/Hz). 10 While there are systems where joints change behavior on a short-term basis during operation (e.g., magneto-rheological and similar semi-active systems), the sigmoid transition is usually driven by two factors: (A) energy thresholds (e.g., bolted joints have a bi-linear load curve [25]), and (B) design changes during the drawing release phase of vehicle development. all joints are structurally nonlinear 11 for the transfer of forces and reactions, a situation that provides substantial coupling nonlinearity.
3. Ordering of modal participation: The low frequency end of the spectrum tends to hoard the majority of the strain-energy, as the quote from Zienkiewicz explains. These characteristics are summarized in the table VI, the Conclusion. Analysis of MPF's (developed by characteristic 2) can determine state transition probabilities. However, modal engineers tend to be focused on solving critical failure issues. Engineers might use similar terms for entirely different kinds of 'modes' such as 'failure mode effects analysis' (FMEA). However, a large part of modal engineering involves creation of a record of the order of the vibration modes, followed by an analysis of the tracks of the modes as the design variable change.
A structural vibration example described below shows how a structural design change can de-couple modes in order to improve ride-and-handling. However, this change has consequences other than aesthetic design constraints. The changes can increase crashworthiness risk and other seemingly unrelated systems, and they can degrade previously adequate vibration and modal engineering balances in the design. For example, stiffening one joint or component usually provides a much more efficient flow of energy through that stiffened element of the structure. Strain-energy is drawn from the high frequency modes to the low modes in a surprisingly efficient manner once such a structural "channel" is open (components have structural coherence). Then the MPFs start to re-balancing as the entire integrated structure starts to equilibrate. So a stiffener can solve one problem and cause another. The engineering group that deals with braking and crash loads that transfer from the front to rear bumper can have its margins go from adequate to negative due to a design change by the ride-and-handling team that solves a mode-coupling problem. This is similar to an old solution to convertible mode coupling:

C. Inter-organizational effects of design changes
The 1980 era Chevrolet Corvette convertibles an initial design had a fundamental bending mode just below 30 Hz, which put it very close to the suspension mode. This is a overall vehicle bending mode that curves about a lateral axis, where the bumpers move vertically and together, both in direct opposition to the vertical motion of the seats. It is usually the lowest frequency vehicle mode for a convertible (as if the vehicle is trying to do sit-ups). The Corvette designers needed to stiffen this mode to reduce coupling to the nearby suspension modes, which were also near 30 Hz. The models in the late 1980's had very tall rockers that satisfied this requirement. This was the most efficient structural fix because the stiffness of a beam to bending vertically is proportional to the cube of its section height [27]. Apart from the aesthetic issue with having to step high to get into the Corvette, these tall rockers also changed other load paths including a for-aft transmission of load, important to crashworthiness, as well as making the vehicle slightly more heavy. Automotive is one of the few industries that has a price on the engineering sufficient to remove a kilogram of mass, which was $50,000 in 1990. These and other costs were accepted in order to move the fundamental bending mode away from the suspension mode to effectively eliminate coupling. In later years, other technologies helped solve this issue.
With the exception of the Pininfarina chassis test results [13], few modal analyses such as these (chassis or body-in-white) are in the literature. To a scientist, the ride-and-handling issues and how they relate to the three observational characteristics described above might seem to be basic enough to warrant several scholarly articles. 12 However, this is clearly the arena of trade secrecy at such a high level of value that corporate lawyers would be unsurprised at the paucity of measured data available to the public.

D. Measured data: Vibrating clamped-clamped bar
In his thesis Fl. LT. Ngoya Pepela (he was an Australian Flight Lieutenant in 2006) showed that the modes are not just maxima of a spectrum, but spikes with huge SNRs in the spectral response [11]. These are the "discrete" modes described above in 'Observation 1.' His Figures 5 and 6 each show a photo of the vibrating bar illuminated by a laser vibrometer, with the corresponding spectral density plot are beneath the photos, plotted in the lower pane of each figure. The first ( Figure 5) shows all the major modes whereas the second ( Figure 6) is missing the modes at 1460 Hz. In Figure 6 the vibrating bar is seen to be illuminated by a laser beam with a large round spot size centered on the 1460 Hz node of the vibrating bar. In Figure 5 the right half of the laser beam of Figure 6 is blocked, illuminating only the left half of the symmetric 1460 Hz mode shape. The halved beam size ( Figure 5) allows the laser vibrometry spectral analysis system to show a stronger 1460 Hz mode return due to the elimination of spatial averaging. That averaging is introduce by phase-related destructive interference from both sides of the 1460 Hz mode shape in the full spot size result of Figure 5.
Initially spectral elimination might appear to be a problem for the laser vibrometry industry for use with non-imaging sensors if the spot size is large. However, Ngoya Pepela's thesis [11] and another article by the author [28], with support from simulation 11 Even spot welds have at least geometric nonlinearities in the Green's strain tensor. This is due to large strain related to dissimilar stiffness and offset load moments ('kick's' in aerospace jargon). Nonlinearities are usually due to frictional hold of bolted or riveted joints that "must be ignored" in most structural analyses [26]. Therefore, most system FEA cannot model MPFs unless the model is nonlinear and designed to develop MPFs. 12 Physicists in Ann Arbor. mostly graduate students in 1998, were baffled, surprised that a wheelbase change of less than 2 cm would require over a year of mechanisms analysis to re-balance just the toe, caster, and camber -in negotiations with competing requirements from groups that engineer tires, turn radius, vibration, and crashworthiness. Most vehicle designs barely consolidate the competing requirements, which is more difficult than it would seem.  [11], [28]. and analytical calculations, show numerically how unlikely spectral elimination is for commonly manufactured items, even if the spot size encompasses the entire vehicle. Furthermore it is necessary, but not sufficient that super-symmetric structures produce spectral elimination in laser vibrometry [7]. The vibrating bar is the simplest form of structural super-symmetry that results in spectral elimination, but it is rarely the only structural component a laser vibrometer can use for identification of a manufactured structural system (vehicle), unless the target is the vibrating bar (a structural oscillator) that requires identification. But even for such oscillators, the typical transducer is a piezoelectric system that is not super-symmetric, and thus not prone to spectral elimination; the source of oscillator energy is separate from the driven bar. Therefore, a vehicle that is purposefully fitted with vibrating bars can still be identified with its spectral "fingerprint" in spite of possible spectral elimination from oscillators. Pepela's result [11] was that use of small spot sizes reduces spectral elimination from the laser vibrometry spectral result ( Figure 6 vice 5).

E. The McKinley observation
In early 2018 Michael McKinley of Arlington, Texas, observed that Figure 5 contains modes (at approximately 2740 and 2900 hertz) that are not visible in the small spot size collection [29] in Figure 6. While the small spot size is still large enough to average out higher frequency (smaller vibration shape wavelength) modes, another spectral elimination that small spot size vibrometers are susceptible to is the effect of the probe spot being on a Chladni line of nodes [30] of the vibration shape for particular frequencies [31]. Detection of these modes at 2300, 2740, and 2900 Hz is a positive observational difference provided by using a large spot size, in this case.
Therefore, spectral reduction and elimination (SR and SE) happens for all spot sizes, large and small, for super-symmetric structures such as 1-D bars. The prior paragraph discusses why this is not a problem for manufactured vehicles; nature provides

F. Observational differences
The difference between small spot size ( Figure 5), and larger spot size ( Figure 6) comprises one dimension of observation variation that can change the observer's view of reality. Physical systems such as the Corvette design example can be used to explain the existence of the three characteristics of modal systems listed above. Contact nonlinearity is the main physical model of interest in this paper because contact is a relatively simple and ubiquitous form of nonlinearity, found in nearly all structural joints.

III. ENERGY ORDERING
Vibration strain-energy transfers from one component to another, usually through joints that are necessarily nonlinear (even spot welds), depending on the amount of 'fixity' of that interface or joint. An extensive treatment of how fixity changes the modes and therefore the sensed optical radiant flux is in the author's thesis [7, B.2.3, p. 187], including nonlinear examples that compare to the contact nonlinearity of the FE model. Fixity is an engineering variable for a particular DOF that determines the ratio from zero to one (0% to 100%) of force, moment, or torsion [27] that will transfer from one component to another through the joint member for whom the fixity is defined [26].
Even "simple" vehicle structures have complicated load paths that are nonlinear (mostly through joints) where forces from one region of the vehicle affect parts elsewhere on the vehicle. These are the pathways for vibration energy flow that tend to dump energy into low frequency structural modes. Modal engineers use isolators and suppressors to stymie this natural tendency, but mostly just for critical components. Some components tend to have non-negligible vibration spectral energy down at these "fundamental" and near fundamental modes 13 even without the nonlinear effects. Damping helps drive the energy into the lowest modes because the resonance spreads (widens) to encompass a larger bandwidth.
Zienkiewicz showed energy transfer using only viscous damping [22,[340][341] to form a ratio of damping to its critical value, 14 c i = 2ω i c i depending on the modal frequency ω/2π, as quoted below. For vibration systems that are driven by a forcing function f, Zienkiewicz describes the differential equation (DE) Mä + Cȧ + Ka + f = 0 using a mass matrix (M), a damping matrix (C), and a stiffness matrix (K), acting on deflection and force vectors a and f . α and β are damping parameters.
Zienkiewicz quote: " ... we have indicated that the damping matrix is often assumed as C = αM + βK. Indeed a form of this type is necessary for the use of modal decomposition, although other generalizations are possible [references given in the book]. From the definition of c i , the critical damping ratio [described above], we see that this can now be written as as is the case with most structural damping, c i grows with ω i and at high frequency an over-damped condition will arise. This is indeed fortunate as, in general, an infinite number of high frequencies exist which are not modeled by any finite element discretizations." [22,340,341] Therefore, vibrational strain-energy tends to naturally migrate from higher to lower frequency modes. Some of the strainenergy can find its way down to the lower modes because over-damping causes more overlap of modal response resonances. Nonlinearities add to the methods of energy transport as discussed later.
This results in a set of modes where the mode number increases monotonically with frequency and has an energy that decreases monotonically with mode number. Possible exceptions to this natural effect include artificial structures 15 or temporary transients such as might be introduced with magneto-rheological fluid or other systems controlled by magnetic or electric field variations. The 1/f phenomenon is found in many different fields of engineering. 16 The local mode structure in Ngoya Pepela's laboratory test measurement in Figures 5 and 6, which do not follow the 1/f phenomenon within the displayed 3 kHz, appears to be a contrary example. This an example of not being able to 'see the trees for the forest' within its part of the response spectrum of the very stiff clamped-clamped structure. For those results, the vibrating bar only had a half-dozen distinct modes across 0-3000 Hz, one of which was diminished in Figure 6. The utility of the structurally super-symmetric bar is that it shows discrete modes, and that the energy in some modes are spectrally reduced or eliminated by spatial "averaging" -integration over the area of a surface of variable phase that introduces varying levels of continuous interference depending in large part on the limits of spatial integration [28].
That the vibrating bar of Figures 5 and 6 is consistent with modal ordering becomes apparent by analysis of the physics of a simple model of a point mass at the end of a spring of stiffness k. This explanation considers multi-variate non-uniform parametric change in the system. The maximum potential energy is approximately U max ≈ kx 2 /2 where x is the largest extension of a simple spring. Assume that the frequency f 2 = k/(4π 2 m) is increased, then the energy grows as the square of the frequency, parabolically as U ≈ 2π 2 mf 2 x 2 . In practice, the deflection x decreases with frequency, as it must. However, there is a practical limit to reduction in deflection. When that limit is reached, while the frequency continues to increase, the potential energy must grow. This means that higher frequency modes require more energy, even for smaller deflectionespecially when deflection is already small, when x 2 U/(2π 2 mf 2 ). The kinetic energy side of this simple calculation is Rayleigh's method [32].
It may help to see this from a force point of view using the gradient of U . If the restoring force is conservative, it is F = ∂U/∂x = 4π 2 mf 2 x, and is thus a quadratic function of the frequency, f . It might be tempting to consider a "conjugate force" for frequency, where the gradient ∂U/∂f = 4π 2 mx 2 f (units of action : Energy × time, or momentum × distance) [33]. However, the force (F = ∂U/∂x), cannot increase ad infinitum with f . At some point assumptions of linearity and structural integrity start breaking down. Combining the observations of this upper limit on the force, and the understanding that the spring extension is limited to x > 0, the energy required to have a particular mode increases approximately quadratically with f . This increase strain energy based on the modal displacement u i ≈ dxdyw(x, y) is "on average" energy in this sensethat the other variables are limited, and using experience that dynamical systems naturally tend to avoid populating the higher modes with energy. Many children experience this effect when they try to forcefully excite a higher frequency mode in a rope, only to fail unless they exert considerable effort to the point of overexertion.
Using the observations above, a conjecture can be stated for the system restoring force of simple harmonic motion, based on the gradient of potential energy, along with the FEA simulation included herein (Figure 7), including industrial experience, and the test data from the childhood of most of us. While system energy decays with frequency "on average," a plot of the 13 The 'fundamental' mode is the mode that has the lowest frequency, contrary to the definition for a musical fundamental frequency. 14 Critically damped systems remove all oscillation. 15 Examples of artificial exceptions to the 1/f phenomenon could be exotic. Perhaps a feedback system might "manually" excite high harmonics without exciting lower harmonics, such tapping very near the base on a taut chord would tend to excite the higher frequency mode. 16 One example is the spectrum of a surrogate missile plume. The power spectral density at low optical frequencies is large, decaying exponentially with higher frequency in an approximately straight line for a log-log plot. spectrum within its lowest half dozen modes might not show the 1/f phenomenon, but that the full pattern does. Sometimes the spectrum plot is 'in the trees' and the "forest" of the 1/f phenomenon over the entire large bandwidth spectrum is not seen within the zoomed-in window of smaller width. For example, consider the 'in the trees' ranges of 150-250 Hz and 400-500 Hz for the much lower fundamental frequency system of clattering armor in the FEA results of Figure 7.
This figure plots several spectra of the simulated structure, 1 × 2 meter clattering plate (homogeneous rolled armor, discussed later). There is one curve for each level of damping (see legend) and they show a different method of inter-modal energy transport. The 1/f pattern shows up in its structural vibration spectrum. A collection of the deformed shapes for transient analysis in the CSC thesis [7] shows how the transient result at each time step is a superposition of normal modes (eigenvectors) populated with energy according to the MPFs Φ(f ). Those normal modes plots show how one mode at one frequency for one plate excites another mode at a different frequency on the other plate because of deflections that line up along the surface to contact the other plate at anti-nodes of a different mode. Again, the lower modes tend to be receiving modes.
Using appropriate structural damping, the frequency response function in Figure 7 is the result of an impulse load that rings all the modes of the FE model. It appears to be natural that MPFs decrease with the mode number of the mode they characterize; in a time-averaged sense, the vector is an array of monotonically decreasing participation values usually measured as the strain-energy for each mode. This is part of the reason for the ubiquitous 1/f phenomenon. At least for structural vibration [34], we can thank the Zienkiewicz section quoted above [22,340,341] provides an explanation of at least some of this 1/f fall-off.

IV. SIMPLE HARMONIC MOTION, NOT SO SIMPLE
A spring-damper-mass single DOF system (SDOF) provides a one dimensional model of the clattering armor plate which can provide analytical solutions to the vibration DEs. Assume 'small deflection,' where strain-energy density is low enough for linear elasticity assumptions (lack of permanent set in the structure being modeled). Also assume that modes are nearly monochromatic functions of sines and cosines.
In a more formal development of the transfer of energy from one mode to another, Lord Rayleigh points out that pure sine or cosine waves do not exist in reality. He derived the differential equations for Newton's theory for isothermal compressive-vacuum vibration 17 resulting in this second order form [35]. The footnote describes these generic variables in Equation 1 for a wave in y with respect to x and time t. These concepts are covered in engineering vibration textbooks [36, p. 30].
Lord Rayleigh then comments on the ability of simple harmonic motion to maintain shape in nature. The extent to which the mode shapes are not harmonic introduces a possible corruption of the pure energy-per-mode concept assumed by MPFs. In the scientific literature a 'mode' is derived from the parameters used for the study of statistics which include the median, mean, and the mode. In a field where Gaussian distributions abound thanks to the physics described by the Central Limit Theorem, 18 the mode is merely the location on the abscissa, in terms of the frequency for a spectrum, at the maximum of the response. In this sense the typical emphasis on energy-in-mode used by modal engineers is appropriate [23] to identify the deformed shape of the eigenvector for that approximate center frequency. However, transfer functions are inappropriate for these nonlinear systems. If transfer functions are modified as Bendat does for simple solitary nonlinear systems [15], which is out of the scope of this work, there may be some consistency between the methods that would be useful to develop, separately.

A. Contact nonlinear response
The author's CSC thesis [7, App F] describes and plots the control law as a sigmoid for the stiffness, which is the slope of the load-deflection curve. The control law is a mathematical representation of the nonlinear stiffness associated with the contact state between armor and hull. This classifier is an arctangent function, a smooth version of the typical bilinear gap element load-deflection table. As the armor and hull make and break contact, the stiffness switches between high and low values of stiffness, respectively. When the armor and hull are in contact, negative deflection (compressive penetration) results and there is a large stiffness due to the large slope of the deflection curve. When the armor has separated from the hull, the nonlinear contact stiffness becomes a deflection curve of small slope. This small value of extra stiffness provides computational stability. The author's thesis [7] contains a description of the stability and application of the control law for this case, demonstrating that the damping DOF and frequency DOFs are no longer related. The following analysis expands on one of Dr. Winthrop's control laws listed in his 2004 AFIT dissertation [42]. His is a different control law, but the analysis methods are complementary.
Phase space (state space) plots [7, Fig 12] This model of a "welded-together" panel system 19 uses the parameters in Equation 3. A "manual," analytical calculation [7, App C] helped validate the FEA results for the simulated vibration created before application of the laser vibrometry model. Ngoya Pepela's optically measured modes provided qualitative support for the application of the laser vibrometry model that used the FEA results to modulate the probe beam. This section and the Appendix describes the former, a structural SDOF model.
In the arctangent gap model [7, Fig 11] (not plotted here), the stiffness will transition at contact, x = 0. Contact surfaces are imperfect due to microscopic protuberances that comprise the roughness of the surfaces. As the two plates come into contact (x ≤ 0) the roughness of the surfaces deform, compressing the protuberances, and a transition from low to high stiffness occurs rapidly in order to match the surrounding material. For a "low" damping system, the damping is µ = 0.02 kg/s and the open to closed dimensionless stiffness ratio is k open /k closed = 0.01. The Appendix describes the dimensioned and dimensionless models. A plot of speed versus gap opening in the thesis [7, Fig 12] reveals that the high stiffness during compression of 17 "Since the relation between the pressure [p] and the density [ρ] of actual gases is not that expressed in [ p = const − (u 2 o ρ 2 o /ρ)], we conclude that a self-maintaining stationary aerial wave is an impossibility, whatever may be the velocity uo of the general current, or in other words that a wave cannot be propagated relatively to the undisturbed parts of the gas without undergoing an alteration of type. Nevertheless when the changes of density concerned are small, [ p = const − (u 2 o ρ 2 o /ρ)] may be satisfied approximately; and we can see from [ dp/dρ = (u 2 o ρ 2 o /ρ 2 ) ] that the velocity of stream necessary to keep the wave stationary is given by [ uo = dp/dρ ] which is the same as the velocity of the wave estimated relatively to the fluid." [35] 18 The extent to which the literature supports the Central Limit Theorem can be found in papers from 1937 through 2003 with Uspensky [37], Landon and Norton [38], Khinchin [39], North [40], and Le Cam [41]. 19 The term welded assumes the plate deflections at both points are the same using an infinitely stiff connection. The closed gap stiffness used in this work is more like hard rubber [7, App B.2.5], but the effect is the same. the base is a shallow orbit in x versusẋ phase space. Changes in damping lead to changes in the range of the orbits during stabilization. These control law formulations were implemented using a MATLAB system [43].
Lyapunov function analysis shows the contact control law is asymptotically stable. The Lyapunov function in this case is the total energy with simple damping loss. In rare cases vibration energy gain to the hull upon which the armor plate clatters matches the cycle's damping energy losses and is in synch with plate contact. Persistent energy loss is due to structural damping, material heat losses due to bending, and fastener frictional losses. The energy gain per cycle is less than the energy loss due to damping [7, App F]. If the losses where small enough to just equal the gains, the system would be 'stable in the sense of Lyapunov.' This analysis above summarizes a clattering system. Details of the physics of boundary conditions, initial conditions, and operating states are found in the thesis [7, App F].

B. Laser vibrometry return from the probe beam
The noise floor did not rise in Figure 6 but rather, the surrounding spatial areas had phase differences within the large spot size that, due to spatial averaging, underwent destructive interference in a bandwidth that removed most of the 1460 hertz mode. 20 This interference is a combination of optical phase shifts along-range caused by reflection from deflection amplitudes that are related to the structural wavelength along the bar for each mode.
The purpose of the vibrating bar measurement was to show spectral elimination, but it also shows the discrete nature of high quality (low damping) structural modes. The FEA of this system [7] complemented this laboratory measurement [11] by showing that, while spectral elimination can occur with structures that might be constructed in the lab (one dimensional modes shapes are the simplest super-symmetric structures), it is impractical to build commercial structures that are substantially super-symmetric.
The main result of the CSC simulations was that academic models can produce reductions of vibration sensitivity that theoretically verify part of Ngoya Pepela's spectral elimination thesis [11], within meaningful assumptions. The single and two DOF models described here show why this is the case, theoretically. A separate 2014 article [28] describes the lack of spectral elimination (or spectral reduction) for vehicle and other types of structures that are manufactured economically.
The lab vibration measurement shared some spectral reduction features seen in the nonlinear cross-spectral covariance simulation [7]. The transient vibration results for the surface of the 3-D armor plate were the input for a MATLAB system that forms an image similar to one that a laser vibrometer would detect. Commercial laser vibrometers typically only provide images processed with their proprietary system. For the purpose of simulating the image of the vibrating plate, the optical model for this work used conventional Fresnel diffraction in a MATLAB system of functions. The Fresnel propagation method introduced by Goodman [18] provides a computationally adequate method to propagate an image of the return from the target back to the detector. In this sense, Pepela's laboratory measurement provided a validation for the FEA and vice versa.
The laboratory measurements in Figures 5 and 6 use a clamped-clamped bar that has a high frequency fundamental mode (compared to vehicle fundamental modes that are below 50 Hz). The spectrum has not dropped off with frequency by 3 kHz because this is still the "lower frequency region" for the structure, as discussed in the section III. It has a stiff fixity meant to test the laser vibrometer for structural wavelengths that would undergo a difference in phase across the beam spot size.
The spectral elimination effect, seen mostly in (Figure 6), provided laser vibrometer manufacturers information on what spot size to suggest or program for their probe beams. They want to avoid such 'spectral elimination.' The simulation and analysis in the CSC thesis [7] shows users and manufacturers that, except for 1-D structures (bars), even specially manufactured structures are difficult to create in the super-symmetric form, and thus do not exhibit spectral elimination.
A detailed study is available that has comparisons of both structural and optically sensed CSCs, to imaging versus nonimaging returns [7]. The latter choice of observation type, non-imaging, analyzes a scalar signal composed of the spatial average of the entire image per time step. The simulation analyzes the surface of homogeneous rolled armor (HRA) using only one scalar metric, radiant flux that is the spatial integration of the irradiance reflected from the HRA. These coupled FEA-optics results, using non-imaging return, identify the modes on the target (HRA clattering on a hull), nearly as well as the imaging version of the CSC.

C. Finite element modeling and analysis
The FEA produces a surface of vibration that is curved in 3-D. The vibrational deflections are the input for the scripts and functions run within MATLAB that produce an optical simulation of the target's image. Unlike the small and large spot sizes (Figures 5 and 6) in Ngoya Pepela's measurement [11], the remote sensing of the armor-plate clattering uses a complete coverage large spot size such that all points on the armor modulate the probe beam. Assume the surrounding clutter to be time-gated or otherwise removed. The input that the MATLAB script uses is the set of displacements to simulate images of the transient dynamical response of the surface of the clattering armor for hundreds of time steps in a duration of a few seconds. One such 'vibration deformed shape' is seen in Figure 8. The FEA representation of the vibration shape was applied to the MATLAB model of the probe beam as a phase modulation. Using Fresnel diffraction [18], the return was imaged onto the detector 4 km away. 21 A 10 µm probe wavelength satisfied the numerical requirements for results shown in Figure 9.
A spectrum of the time sequence of sensed images, such as the two shown in Figure 9, appears in Figure 7. More detailed plots are available [7].
V. NONLINEAR EIGEN-STATES Creating a simple model for a clattering structure may seem easy, until the fact that the eigenvalues must be nonlinear complex functions becomes apparent. A two degree of freedom (DoF) model developed in the Appendix herein was meant to validate the FEA. This FEA validates a major dynamical 2 DOF analysis result: Antisymmetrical and unsymmetrical modes increase in frequency for increases in contact stiffness. 22 Investigation of these closed form simple damped sprung mass systems 23 provides not only insight, but qualitative numerical validation. Whether an eigenvalue increases for anti-symmetric and un-symmetrical modes when gap-closed stiffness increases. This simple model result provides another analytical basis for results seen in the FEA. Additionally, the behavior of the modes indicates that symmetrical modes are better target identification features than un-symmetrical modes. Hence, there is a need to analyze the energy balance and stability [7, App F] to validate the simple one and two DOF models as is summarized in the Appendix herein, including the time variation of amplitude and phase.
Using the SDOF model as a basis, the Appendix uses Equation 4 to provide the 2DOF solution. Equation 4 may appear deceptively simple because the time variation of the amplitude and phase is not apparent until explicitly formed. Eigenvalues 21 For the computers at that time, even when running MATLAB simultaneously on several machines on the Air Force Institute of Technology cluster, there was a delicate balance between adequate structural grid densities and optical grid densities for the spatial Fourier transform that performs the Fresnel propagation of the return. (It is more of a numerical modeling issue than a matter of compute power.) 22 The thesis [7] defines several forms of symmetry. These symmetries include structural, clattering (as opposed to synchronized hull and armor), and modes similar to the deformed shape of the components to inertial loads. The latter mode is similar to a mode shape (eigenvector) where the hull and armor move synchronized together, w hull (x, y) = wMRA(x, y) from which comes the 1 degree of freedom (DOF) and 2 DOF models where the plates are averaged over x and y into two point masses. The 1 DOF model grounds the "hull" mass point to zero displacement, so it is removed in the 1 DOF model. 23 'Sprung mass' systems are an automotive term for a systems where vibration is being isolated or suppressed. At an academic or 'free body diagram' level, they can sometimes be approximated by spring-mass-damper systems. λ i are complicated functions of natural frequencies of individual modes which require 'mode tracking' [44] with respect to stiffnesses [45], due to their transient nature in reality [46]. Figure 10 plots the real part of the dimensionless eigenvalues of Equation 4 as a function of dimensionless frequency versus dimensionless closed stiffness k closed , scaled by the hull stiffness k 1 for a closed stiffness ≡ 1/ = k open /k closed [7]. Two models appear in Figure 10. The hull system is eight times stiffer than the "armor" system.
From a chaos point of view, nonlinear attractors exist for the system defined by Equation 4. We know this because the systems exist in reality and their vibrations modes do fluctuate, albeit not monotonically. The nonlinear attractors are related to the underlying hull forcing functions (D'Alambert's forces due to vehicle inertial loads), the timing and shape of the clattering, the plate stiffnesses, the structural and added damping (e.g., washers or armor-hull batting), and the fixity of the joints.
Initially the two eigenvalues for the 2DOF system are equal. As the closed stiffness increases, a critical value of stiffness, k crit , is reached where the armor can no longer follow the hull. It starts to clatter in an unsymmetrical mode. A 3-D contact surface cannot be anti-symmetric unless it is axisymmetric (1-D), as are the SDOF and 2DOF models. When the stiffness exceeds k crit the two curves become distinct, branching into two separate paths. This can be see for the plots of two models where the base stiffness is 128 and 1024 N/mm. In the region before clattering, < k crit , the slope of eigenvalue curves, for the 8×k closed higher base stiffness curve, decreases eight-fold compared to the smaller base stiffness. At this higher 1024 N/mm base stiffness there is also an eight-fold increase in k crit . The intercept also increases 8-fold. Derivations of these concepts and further results, including how the imaginary parts of the frequencies (not shown, Figure 10 is the real part) contribute to the transfer of energy between the modes are developed in the Appendix.

VI. ANALYSIS AND CONCLUSION
A system model of the structural vibration of a contact-plate system was shown to have surface waves with low frequency modes that are effectively discrete modes of the PSD [17]. This clattering armor system provided an example Fig. 10. Real parts of the eigenvalues for the dimensionless DE show why higher-energy antisymmetric modes are less likely to be excited. Above a k crit the symmetric (lower frequency branch) and antisymmetric (higher frequency branch) modes [36,167] for this two DOF problem break out into modes of well separated energy.  Fig 7 FRF, coh(f ), MPFs Typ. Figure 7 (super-sym. target is rare) of nonlinear response. Both SDOF and 2DOF models derived in the Appendix have qualitatively verified the symmetrical versus unsymmetrical modes, and the higher frequency of the unsymmetrical clattering modes. The full finite element model transient results also showed the tendency for vibration strain-energy to collect in the lower frequency modes, as most model engineers have seen [13], and as was calculated analytically in the classical literature [22,340,341].
(1) The lab tests on the doubly clamped bar show experimentally that high SNR modes are effectively discrete. In the absence of substantial damping they are not just the maxima of a spectrum, but spikes in the spectral response [11]. The system in Figure 5 has a signal-to-noise-ratio that is huge. These massive SNRs show that the lower modes of structure, that are fairly high in frequency for this case due to the double clamped nature of the bar, can easily be considered discrete. Hence, quod erat demonstrandum (Q.E.D.), the observational characteristic (1) is demonstrated -lower frequency modes are essentially discrete.
This analysis of nonlinear eigenvalues provides useful results for theory and simulation for common nonlinear structures. Symmetrically moving parallel plates have more strain (and thus more strain-energy) than clattering plates where the unsymmetrical motion interrupts their nearly sinusoidal in time out-of-plane-motion to spew strain-energy into acoustics and even permanent set (deformation of the material). Through simulation 24 and by experience the observational characteristic (1) appears true, that lower modes are effectively discrete for high quality systems (low damping), and that, when there are pathways (nonlinear joints) that allow energy transfer, (2) and (3) are in effect; damping and clattering help the energy flow into the lower modes from higher frequency modes -although overall energy is lower as damping increases (see Figure 7). Most vehicle components have bolted or riveted joints that allow such energy transfer [15] as can be seen with plots of coherence spectra, coh(f ) [47].
(2) The tools in industry that provide MPFs include mature engineering methods [17] from the 1980's such as the FEA tools put in place, for example, by MacNeal-Schwendler Corporation engineers in NASTRAN [16]. These MPFs stored in the FEA vector Φ can measure the energy transmitted between states (between eigenvectors of the system) where the state (mode) i changes it energy with mode j.
For surface waves on the ocean [48] this observation would include wave energy of large ocean swells driving up capillary waves and eventually creating foam as seen in large storms. Energy transfers from one type of wave to another, Rossby-Kelvin, gravity-capillary, internal-surface, ..., and direction variations thereof, may prove useful.
Modal engineering for industry, especially for vehicles (i.e., ride-and-handling) is a trade secret endeavor. Vibration spectral plots for vehicles found in the literature, such as the Pininfarina paper [13] are rare. It makes use of analysis tools in order to calculate modal participation factors [16]. These MPFs show the 'participation' (energy per mode [17]) using modes developed from normal modes that are FEA-produced eigenvectors. The participation flows from energetic modes to lower energy modes similar to how large swells on the ocean in a storm are accompanied by unsettled surfaces, rather than smooth large waves that a tsunami has before it nears the shore [48]. In the latter case, not enough time has passed to transfer the energy to other modes until the wave crashes on shore where sufficient coupling to other modes exists due to the constraints of the shore structures.
Therefore, due to (3) the energy ordering of modal states for complicated system synthesis models and the ordinary vehicles they represent, the work of vehicle design for these issues focuses on (2) transition probabilities 25 and (1) energy levels Φ(f i ).
There are usually other types of oscillator interactions that produce an ordering of energy levels similar to that described by Zienkiewicz [22,[340][341]. These concepts are summarized in Table VI. Pencil-thin probe beams like that used in the lab measurement have observational characteristics quite different from the large spot size used for the FEA-MATLAB model of fully illuminated clattering armor. The former can be less susceptible to spectral elimination. However, Mr. McKinley's observations discussed in section II-D on page 6 show different modes appear and disappear for either change in spot size. The latter low fidelity beam method is adequate for spectral identification of economically manufactured vehicles. Such full coverage probe beams are less likely to be subject to spatial coherence issues or illuminate solely a node of the Chladni zone [31]. Pencil-thin beam returns fail to convey vibration modulation in this manner [30]. For the large (full coverage) spot size, some of the beam will nearly always get through.
While these results are generic for plate structures, application to waves on the surface of a volume appear to fulfill similar behavior of (1) discrete modes, (2) modal participation dissipation, and (3) higher energy in lower frequency modes. A tsunami has a much lower spatial frequency and temporal frequency than typical 2-4 meter gravity waves. However, the larger Rossby and Kelvin waves are also candidates for the study of the extreme low frequency application of the principles observed in this article.

VII. APPENDIX -ONE AND TWO DOF CONTACT EIGENVALUES
This Appendix summarizes the nonlinear dynamics of a simplified one-dimensional (1-D) and 2-D forms of a structural contact system, its relationship to the full 3-D FEA, and in the end, test measurements. The 2-D form provides for a difference in foundation stiffness, and transition to 3-D.

A. Closed-form SDOF nonlinear contact response
This subsection displays the Mathcad™output for the closed form solution to the damped SDOF oscillator meant to represent a lumped mass model of the HRA-hull contact vibration system. The mathematical derivations are based on the nonlinear solutions [42] with a simplification of the control law, u(x), that models simple contact along the x-axis.
The SDOF dimensionless system describes the effect of nonlinear contact stiffness from a solution composed of symmetric and antisymmetric one dimensional modes. This SDOF solution applies to the symmetrical and unsymmetrical modes of the 3-D FE model for the HRA-hull system, respectively. In 3-D the unsymmetrical modes are non-uniform, occurring when some parts of the armor is moving opposite to the hull (i.e. clattering). For clarity 1-D modes are distinguished as symmetric or 24 Figure 7 shows that decreasing baffle stiffness between the plates (see its legend) drives unsymmetrical resonances lower. 25 Transition probabilities are related to CSCs [7] and the coherence spectrum [47]. antisymmetric, while 3-D systems are labeled symmetrical 26 and unsymmetrical. 2 Application of derivatives to the dimensionless system is first made under the assumption that all variables are nonlinear functions of time. Then after starting with a restricted case, the variables are brought into explicit nonlinear use, one at a time, to refine the calculation.
The state-space representation [49] is shown here in its phase space form (location x(t) and derivatives). 27 With the definition cos φ ψβ ≡ cos (ψt + ψ o )t + (βt + β o ) , 28 Even the first derivative of Equation 5, is non-trivial. Mathcad allows factoring in several ways. Collecting on cos φ and then on sin φ, Equation 6 provides a compact expression for the undamped speed in Equation 6.
With more assumptions restricting the nonlinearity of the solution where appropriate, the acceleration is also collected on cos φ and then sin φ in Equation 7 for this simplified expression: The first order nonlinear solution in Equation 7 uses a constant amplitude, a = constant. While not explicitly a function of time, a has a constant time rate of change,ȧ. A further order of nonlinearity to allow the amplitude change rate to be a function of time would follow the above assumption with the use ofψ, the time rate of change of a dimensionless frequency. The following subsections describe how Figure 10 shows that while the stabilizing stiffness increases above the critical liftoff frequency, symmetric modes remain at a constant frequency while antisymmetric modes increase in frequency. This 1-D behavior implies the same effect for 3-D modes, symmetrical and unsymmetrical, as is seen in the FEA results [7].

B. Closed-form Mathcad 2 DOF contact
Mathcad symbolic equations in the following calculations comprise the closed form solutions for a damped two DOF (2DOF) contact vibration model. Stiffness k x is the "rate" that defines the sharpness (hardness) of the contact.

C. Two DOF DE's and solutions
Equation 8 is a dimensioned form of the two DOF damped oscillator DE where P , m, k, d, , ξ, and t are applied force, outboard mass, foundation stiffness, damping, control law stiffness, surface displacement, and time. 26 The symmetrical 3-D mode shows up in frequency response curves [7] but it is not perfectly synchronous across the surface, thus it is not a perfectly symmetric mode like the 1-D model. It exists in 3-D the same way quantum mechanical modes exits -because of their energies -spikes in the spectra at the correct frequency, not because we saw them optically. Although with FEA plenty of studies on modal decomposition are in the literature to reinforce the combination of modes that modal engineers use to analyze vibration. The symmetric (1-D) and symmetrical (3-D) modes do not clatter while the other modes do. 2 27 The state-space definition of the system combines the state x(t) with its output y(t) and changes in state (∂x/∂t). These matrix relationships [50] of the state change vector (ẋ) and system output and its state and input (u(t)) are one of many sets of A, B, C, and D matrices, easily confused with other field's ABCD systems such as optical ABCD ray matrices for laser resonator orientation [51], [52].
x = Ax + Bu using the state and input matrices y = Cx + Du via output and feed-through matrices All matrices and variables can be functions of time. Some definitions distinguish state space as discrete compared to continuous phase space. However the field of linear systems often uses continuous output and even state variables and the resulting Kalman filter is nearly ubiquitous [49]. 28 ψ and β are canonical rate terms for frequency and phase. Their subscripted 'o' terms are the initial frequency and phase values. The 'naught' terms have units of hertz and radians while the un-subscripted terms are time derivatives of frequency and phase.
The dimensionless ratio k = (k 1 + k 2 )/k 1 is the stiffness from the second oscillating point mass (the 1-D representation of the armor) to the base of the system as a whole, which includes the hull for this SDOF system. The hull is freed to oscillate in the 2DOF system. The dimensioned variable ξ is the gap opening in the dimensionless x direction. The dimensioned control law for structural contact for this two DOF problem is the displacement u( − → ξ ) of Equation 9. The dimensionless stiffness from the first oscillator to ground is unity.
ξ indicates the dimensioned form of gap opening in the dimensionless x direction. The dimensioned control law for structural contact for this two DOF problem is the displacement u(ξ) of Equation 9.
The use of Winthrop's method [42] on Equation 8 produces the dimensionless DE in Equation 10. The dimensionless applied load is F applied and the dimensionless stiffness only applies to DOF one: Continuing with Winthrop's assumptions [42] assumes a straightforward solution. Equation 11 shows all the variables that vary with time.
x i = a i (t)e µi(t)t cos(ψ i (t)t + β i (t)) Substituting assumed solutions from Equation 11 into Equation 10 results in a Special Eigenvalue Problem (SEVP). The dimensionless system frequencies ψ i are functions of the individual frequencies . Two assumptions help make the solutions tractable. First, the removal of the driving load provides the homogeneous solution. Second, variables dependent on time and space (uniformity with respect to location) vary differently. To first order the time variation of the dimensionless frequency, ψ, is the derivative in this nonlinear system that has the largest effect on the solution [42]. Equation 12 summaries some of the simplifications: In words, these simplifications are: • Constant amplitude, a, cancels out of the DE.
• But uniform phase, β, is realistic. x = − aµe −µt cos φ + aψe −µt sin φ Maintain the coefficients as variables of time in the subsequent derivative: x = a(µ 2 − ψ 2 )e −µt cos φ + a(µψ + µψ)e −µt sin φ (14) Using Equations 13 and 14, time derivatives formed below come from application of these differential operators: These relations use a different dimensionless stiffness, ≡ 1 = k closed /k open where k closed = k 2 and k open k 2 for numerical stability. Equation 16 is of the SEVP format as shown below in Equation 17.
Equation 18 shows the format of the system matrix A for submission to an eigen-solver. Some of the terms in A 1,1 and A 2,2 were kind enough to cancel.
Solution of the characteristic polynomial, |[A] − λ[I]| = 0, provides the eigenvalues of the system where both DOF have equal displacements. Using a Mathcad '97 symbolic solver 29 two solutions are available: (1) Free classical vibration, and (2) a "welded" system vibration solution. The latter welded components (the two masses) is equivalent to being epoxied with hard rubber [7]. As defined in Equation 9, the welded model reduces to free vibration for slight contact, u = 0, just before compression when u turns negative. The system dynamics are revealed via the welded solution as the control law is varied in the range 0 ≤ u ≤ 1. The switching control law, u(x), determines which frequency is active, switching from one frequency value for the open gap state to another for closed gap state. The nonlinear FEA results for 3-D show that both modes are active at the same time at different locations on the surface. Spectral energy flows into and out of both 'open gap' and 'closed gap' modes when averaged over many cycles of time. The contact gap system is a time composite system with dynamics that are more easily modeled using modal analysis techniques. The noise and vibration industry developed these techniques with nonlinearities such as contact in mind [24], [53], [54], [55].
States evolve in an oscillatory manner. This energy flow between alternately created and destroyed states might be a useful model for other systems where the physics appears to forbid continuous energy level ordering. For vehicle structures, experience indicates that the oscillation period is usually less than an hour, on the order of minutes.

D. Two DOF nonlinear contact eigenvalues
Eigenvalues in Equation 19 are a dual set for the control law, u(x), for an open and closed gaps. In the dimensionless system for the open gap state the stiffness for DOF 1 is normalized to unity and DOF 2 has a stiffness of k = 1 + (k 2 /k 1 ). When the system goes into contact u(x) adds to DOF 1 which becomes 1 + , as seen in Equation 10 and 16. This addition uses the FEA maxim that "stiffnesses add" [56]. The total closed stiffness is 1 + + k. The derivation of the complete dimensionless "eigenvalues" λ i in Equation 19 are Mathematica results for the solution of the eigenvalue problem described in Equations 17 and 18 above [7].
For subsequent nonlinear calculations assume u = 1 so that contact is active; there is no gap between the HRA and the hull. Assume > k (dimensionless Equation 16) and for stability 1 k 1 > 2 k 2 = k base (dimensioned 2 DOF Equation 8). Note that all epsilons are dimensionless but they also are a part of the control law that modulates the dimensioned system, Equation 8. Figure 10 on page 15 plots the real part of values of the eigenvalues. DOF 1 and 2 represent the hull and armor, respectively, in the 2DOF system. After a critical , the higher frequency mode (the clattering, antisymmetric mode) separates from the symmetric mode. When the hull has a stiffness that sufficiently exceeds the critical stiffness, the symmetric mode settles into a constant eigenvalue with respect to regardless of the dimensionless stiffness between the two masses. (∂λ 1 /∂ → 0 and ∂λ 1 /∂k → 0 but not for λ 2 .) From the 3-D FEA results in Figure 7 on page 10, the higher frequency clattering mode is the lower energy mode. This is consistent with experience and the Zienkiewicz quote on page 9 [22,. The excitation an antisymmetric mode requires more energy for the same deflection amplitude as an equivalent symmetric mode. If there is an avenue for vibration strain and strain energy to flow into a lower frequency mode the energy "conduit" can be a joint or structural connection (nonlinear structural components). Plucking a taut cord close to the held end, to vibrate it at higher frequency will similarly excite the lowest frequency mode given enough time.
Application of the SDOF model over all DOF's using a theory of linear structural response gives a relationship for a linear transfer function [47] (the Fourier transform of the impulse response [49]). Each of the many DOFs in a linear time-invariant (LTI) system has a relationship that has the same form as shown in a transfer function [7, p. 81]. But these functions are often inapplicable to such nonlinear systems. Common vehicles are filled with nonlinear joints, all of which can be adequately modeled in FEA, with some effort to avoid misapplication 30 of St. Venant's principle [57], [58], [59]. For the ensuing 2 DOF system, the results of the prior subsections and the 2 DOF extension of the SDOF model are formulated with Mathematica.

E. Eigenvalues: Damped 2 DOF sprung mass
Repeated roots occur for the eigenvalues λ( , k) plotted in Figure 10 on page 15 when the dimensionless closed stiffness is less than a critical stiffness, < k crit . Repeated roots indicate the symmetric and antisymmetric modes have the same frequency. However, since there are real and imaginary parts to the eigenvalues [36,171], that there is a growth in energy for one mode and a decrease in energy for the other. For most joints that undergo cyclic loads the contact frictional footprint area also oscillates, which causes the stiffness to oscillate. Therefore, the imaginary parts of the eigen-frequencies act to re-balance the strain-energy according to the control law u(t). Over time energy will move from the higher energy antisymmetric mode to the lower frequency symmetric mode, unless the lower mode is suppressed.
The repeated roots also show that when k base is small enough, the two masses that represent the HRA-hull system vibrate as if they were attached and in free space. In this case the two masses vibrate together, toward and away from the base, with negligible relative motion. For this small k base trivial solution the center of mass is oscillates according to f 2 CG = k base /4π 2 [m hull +m MRA ]. This is the symmetric mode. Alternately, they can vibrate apart and together with their center of mass remaining stationary. For such antisymmetric modes the increase in the eigenvalue that occurs when gap stiffness increases provide a basis for similar results seen in the unsymmetrical FEA modes for stiffer bolts and batting material. The hull and armor 2DOF model will allow DOF 1 and 2 to have different unsymmetrical frequencies sin in general m hull = m MRA . Where the deflection at higher frequencies does not modulate of the probe beam as well as lower frequencies, this nonlinear eigenvalue behavior has the result that symmetrical modes are better target identification features than the unsymmetrical modes. 3) The combination of low but sufficient base stiffnesses ("open" stiffness) and high "closed" stiffness is stable. (4) For the analytical results, a seemingly severe assumption sets dimensionless frequencies to be equal, ψ 1 = ψ 2 . The FEA more properly analyzes damping and simulates these results via transient nonlinear time-integrated results using Newmark-Beta methods for algorithmic stability [22], [56], for part of the solutions for the proper 3-D dimensioned DEs for all modes. Therefore, the millions of 3-D DOF in the FEA remove this last assumption.
In 1-D each of the two DOFs obey the assumed nonlinear solution in Equation 20 where the amplitude is x max which is scaled by a characteristic length L * providing a dimensionless solution x(t). L * is best chosen by using the Buckingham π theorem [56], [60] could range between the average of the microscopic roughness height (microns) to the gap averaged over the surface area in the 3-D model where some part of the HRA-hull system is barely touching. The latter is most easily be accomplished using the normal modes along with a few pertinent static deformed shapes.
The dimensionless frequencies, ψ, in this system matrix occur in the same form as an eigenvalue λ from [A 1 ] = |[A] − λ[I]|. This form where ψ 2 1 = ψ 2 2 = λ i are the eigenvalues, only applies to the above assumptions of (2) equal damping and (4) equal frequencies for the two masses which are necessary for a practical nonlinear solution. There is a physical superposition of both symmetric and antisymmetric modes at the same time, resulting in hull -HRA resonation at both ψ 1 = ± √ λ 1 and ψ 2 = ± √ λ 2 , in a linear combination of modes. Each of the two Λ i modes will have both ψ 1 and ψ 2 active for that one mode. Therefore, ψ 1 and ψ 2 represent the same energy for λ 1 (symmetric vibration), but ψ 1 and ψ 2 produce different energy for λ 2 (unsymmetric vibration).
The FEA eigenvalues (i.e., values of the diagonal NASTRAN Λ matrix) are related to their eigenvectors − → u i T , as plotted in the thesis [7,. The FEA modal frequencies are the square roots of the eigenvalues, summarized therein [7, Table 7, Fig 20]. Those frequencies and plots of − → u i T represent eigenvalues λ i within [Λ] and eigenvectors Γ i within [Σ], in analogy to the 2DOF model, where the columns of [Σ] are the eigenvectors. As discussed above, the similar-frequency argument makes physical sense, but the similar damping assumption is an artifice used to obtain a practical, simplified solution that is useful for comparisons. The simplification allows the extra ψ terms to cancel. Experience and the FEA results both validate that this technique is adequately appropriate for use in this particular case.
The eigenvalues in matrix form are [Λ] = [λ 1 , 0; 0, λ 2 ] but it is more convenient to display them in the vector form (Equation 22), where the parameter κ was defined in Equation 19.
The 2 DOF eigenvectors in Equation 23 are one dimensional mode shapes for displacement of the hull, u, and the armor, v, for modes 1 and 2:

G. Analysis of the 2 DOF SEVP DE
The eigenvalues λ i = ω 2 i for each column Γ of the eigenvector matrix (Σ) are the result of solving Equation 21. Equation  24 represents an extreme fixity of this Mathematica™solution, a low frequency "DC" limit ψ = 0. lim ψ=0 − → λ = 1 2 The parameter κ, defined in 19, in this first order correction to the linear eigenvalues (defined for Equation 22) is only a function of the stiffnesses, the square of the base and hull stiffnesses. Compared to prior Mathematica™results, this ψ = 0 model provides a "DC" eigenvalue solution that is otherwise not available.
This relation is only valid for the assumptions described earlier: (2) µ 1 = µ 2 and (4) ψ 1 = ψ 2 . Otherwise more nonlinear terms remain and the system is not susceptible to SEVP solution for modes in Equation 25. . This is only approximate because of the many combinations of nonlinear and approximately linear variables (e.g. "linear" ψ) selected in [42] and used the "Mathematical Preliminaries" section of [7]. .
Equation 26 represents a 2 DOF system with springs of stiffnesses k and , and uniform damping of the same magnitude for both DOFs. The eigenvectors and eigenvalues are most easily recognized in relation to the standard EVP in Equation 27 assuming (λ i = ψ 2 ∀ i ∈ [1, 2]).

I. Synthesis of unmatched DE
For clarification and to bring us back to the full nonlinear DE, stipulate that the eigenvalues for each DOF are the same for each mode λ 1 ≡ λ, that the physical argument leading to Equation 20 holds the frequencies equal. Therefore, the frequency subscript matches the subscript for the eigenvector, rather than matching the damping modes as was done in the prior subsection; the combinations of frequency and damping are (ψ k , µ i ) ∀ i = k.
The control law, u(ξ 1 ) of Equation 28 is the nonlinear stiffness [7, Fig 11], based on the relative displacement, ξ 1 = L * x 1 . The arctangent switch changes the stiffness between the hull lumped point mass and the HRA plate lumped point mass, k open ↔ k closed .
The phase φ k = ψ k t + β(t) ≈ ψ k t relates to the eigenvector whose temporal dynamics λ i describes both DOFs, the hull and HRA point masses. Here the frequencies ψ k are not matched to the damping µ i . Redefining Π i and Ω i for this mixed-index format: Equation 29 provides an independent identically distributed (iid) estimate of the dynamics modeled as N iid oscillators. Since oscillators for each DOF in a FE model and in continuous media (reality) are not independent nor of a random distribution in a structure, this relationship merely guides independent behavior that quickly effects other DOF's. More importantly, the "nonlinear" eigenvalue behavior in Figure 10 on page 15 shows how the symmetrical modes are affected by stiffness change [7, p. 148].
A Lyapunov function of energy balance defined for this nonlinear contact system in the thesis [7, App F] validates the stability of these composite SDOF and 2DOF relations.