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

1.1. 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. Fig. 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 [11], as later discussed in Fig. 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 (Fig. 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 (Fig. 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 infinite in density. 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 (Fig. 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 least one example appears in the literature for automotive modal analysis, ride and handling [12]. Helicopter data is available from MIL-HDBK-810 [13]. A transfer of energy from a strongly driven mode to other modes appears in Fig. 2. That driver might be in a component located far from the receiving modes [14], 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 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.

Fig. 1The 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

Fig. 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 3, provides an explanation for this common 1/$f$ dependence of modal energy. The exceptions noted also include strong damping or added mass.

Fig. 2Vibrational 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, f driven , 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=|fi-f driven |

The models, structures, and observational characteristics that comprise this article are summarized in the Conclusion in Table 1.

Fig. 3As described by the Zienkiewicz quote in Section 3, 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

1.2. 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) [15, 16]. The FEA surface vibration results are input for the MATLAB image propagation calculation (Fresnel propagation [17]). The physical schematic for the simulation appears in Fig. 4. Laser vibrometers use many different methods that can have an effect on the precision of the measurements.

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 Fig. 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 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 3.

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.

Fig. 4Coherent 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

2.1. 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 [16] describe the dynamics of energy flow by using vectors of MPF’s that span the pertinent modes [15]. Each MPF is an element in vectors $\mathbf{\Phi }\left[f_i\right]$ 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 $\mathbf{\Phi }\left[f_i\right]$.

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 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.

2.2. 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. The lower modes relate to mode 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 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 all joints are structurally nonlinear 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 Table 1, 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 variables 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 MPF’s start to re-balance 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:

2.3. Inter-organizational effects of design changes

The 1980 era Chevrolet Corvette convertibles initial design had a fundamental bending mode just below 30 Hz, which put it very close to the suspension mode. This is an 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 [26]. 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 [12], 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. 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.

2.4. 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 [10]. These are the “discrete” modes described above in Observation 1. His Figs. 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 (Fig. 5) shows all the major modes whereas the second (Fig. 6) is missing the modes at 1460 Hz. In Fig. 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 Fig. 5 the right half of the laser beam of Fig. 6 is blocked, illuminating only the left half of the symmetric 1460 Hz mode shape. The halved beam size (Fig. 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 Fig. 5.

Fig. 5The image from a small spot size makes a spectrum with a f3 = 1460 Hz mode. Permission N. Pepela [10, 27]

Fig. 6A larger laser spot size creates a laser vibrometry spectrum where the f3 = 1460 Hz mode phase-averages away

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 [10] and another article by the author [27], with support from simulation 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 [10] was that use of small spot sizes reduces spectral elimination from the laser vibrometry spectral result (Fig. 6 vice 5).

2.5. The McKinley observation

In early 2018 Michael McKinley of Arlington, Texas, observed that Fig. 6 contains modes (at approximately 2740 and 2900 hertz) that are not visible in the small spot size collection [28] in Fig. 5. 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 [29] of the vibration shape for particular frequencies [30]. 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 a mitigation (vibrometry is a useful classifier system for realizable structures) of the perceived problem that the McKinley observation clarifies (small spot sizes do not mitigate SE).

2.6. Observational differences

The difference between small spot size (Fig. 5), and larger spot size (Fig. 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.

4.1. 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 [39]. His is a different control law, but the analysis methods are complementary.

Phase space (state space) plots [7, Fig. 12] validate this model’s behavior. The arctangent function that does the switching in this simple closed form model appears in state space formulation as shown below in Eq. (2), which has the constants listed as Eq. (3). The state space description for the single mass speed, $v$, and acceleration, $\dot{v}$, shown in Eq. (2) is a function of the open and closed stiffnesses, $k_{open}$ and $k_{closed}$, for a damping coefficient of $\zeta =d/m$ in units of hertz, where $m$ is a point mass:

This model of a “welded-together” panel system uses the parameters in Eq. (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\le 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 $\mu =$ 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 the base is a shallow orbit in $x$ versus $\dot{x}$ 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 [40].

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 were 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].

4.2. Laser vibrometry return from the probe beam

The noise floor did not rise in Fig. 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. 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 [10] 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 [10], within meaningful assumptions. The single and two DOF models described here show why this is the case, theoretically. A separate 2014 article [27] 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 [17] 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 Figs. 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 3. 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 (Fig. 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 non-imaging 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.

4.3. 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 (Figs. 5 and 6) in Ngoya Pepela’s measurement [10], 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 Fig. 8.

Fig. 8FEA of the deformed shape (green) for the clattering armor-on-hull model. Exaggerated micron deflections would not be visible compared to the large physical geometry of the plates. There is a 14.24 ms time difference (many time steps) between these two output time frames. The views are into the fore-shortened long direction (1 m) of the 0.5 m wide plate. The un-deformed shape is a light tan color

The FEA representation of the vibration shape was applied to the MATLAB model of the probe beam as a phase modulation. Using Fresnel diffraction [17], the return was imaged onto the detector 4 km away. A 10 μm probe wavelength satisfied the numerical requirements for results shown in Fig. 9.

A spectrum of the time sequence of sensed images, such as the two shown in Fig. 9, appears in Fig. 8. More detailed plots are available [7].

Fig. 9Simulated magnitude of radiant flux image of the return after modulation by the vibrating plate in Fig. 8. The image horizontal width appears rotated 90 degrees from the long fore-shortened axis in the prior image in order to better display the 20×40 mm detector grid shape. The maximum deflection “hill” along the 1 m length shows up as a vertical maximum magnitude in the lower (later) pane of both figures

7.1. 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 [39] with a simplification of the control law, $u\left(x\right)$, 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 antisymmetric, while 3-D systems are labeled symmetrical and unsymmetrical. 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 [47] is shown here in its phase space form (location $x\left(t\right)$ and derivatives). With the definition $\mathrm{c}\mathrm{o}\mathrm{s}{\varphi }_{\psi \beta }\equiv \mathrm{c}\mathrm{o}\mathrm{s}\left[\left(\dot{\psi }t+{\psi }_o\right)t+\left(\dot{\beta }t+{\beta }_o\right)\right]$:

Even the first derivative of Eq. (5), is non-trivial. Mathcad allows factoring in several ways. Collecting on $\mathrm{c}\mathrm{o}\mathrm{s}\varphi$ and then on $\mathrm{s}\mathrm{i}\mathrm{n}\varphi$, Eq. (6) provides a compact expression for the undamped speed:

With more assumptions restricting the nonlinearity of the solution where appropriate, the acceleration is also collected on $\mathrm{c}\mathrm{o}\mathrm{s}\varphi$ and then $\mathrm{s}\mathrm{i}\mathrm{n}\varphi$ in Eq. (7) for this simplified expression:

The first order nonlinear solution in Eq. (7) uses a constant amplitude, $a=$ constant. While not explicitly a function of time, $a$ has a constant time rate of change, $\dot{a}$. 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 $\dot{\psi }$, the time rate of change of a dimensionless frequency. The following subsections describe how Fig. 10 shows that while the stabilizing stiffness increases above the critical lift-off 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].

7.2. 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 represents the sharpness (hardness) of the contact.

7.3. Two DOF DE’s and solutions

Eq. (8) is a dimensioned form of the two DOF damped oscillator DE where $P$, $m$, $k$, $d$, $\epsilon$, $\xi$, and $t$ are applied force, outboard mass, foundation stiffness, damping, control law stiffness, surface displacement, and time:

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 $\xi$ 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\left(\overrightarrow{\xi }\right)$ of Eq. (9). The dimensionless stiffness from the first oscillator to ground is unity.

$\xi$ 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\left(\xi \right)$ of Eq. (9):

The use of Winthrop’s method [39] on Eq. (8) produces the dimensionless DE in Eq. (10). The dimensionless applied load is $F_{\mathrm{}applied\mathrm{}}$ and the dimensionless stiffness only applies to DOF one: $\overrightarrow{\epsilon }=[{\epsilon }_1;0]$:

Eq. (11) uses the vector $\overrightarrow{x}=[x_1\mathrm{}\mathrm{}{x}_2{]}^T$, the dimensionless location of the two masses where $x_2-x_1$ is the gap opening. Continuing with Winthrop's assumptions [39] assumes a straightforward solution. Eq. (11) shows all the variables that vary with time:

Substituting assumed solutions from Eq. (11) into Eq. (10) results in a Special Eigenvalue Problem (SEVP). The dimensionless system frequencies ${\psi }_i$ are functions of the individual frequencies $f_i={\omega }_i/2\pi$ from ${\omega }_i^2=(k_i/m_i)-2{\zeta }_i^2$ where ${\zeta }_i=d_i/2m_i$ from Eq. (8). 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, $\psi$, is the derivative in this nonlinear system that has the largest effect on the solution [39]. Eq. (12) summaries some of the simplifications:

In words, these simplifications are:

• Constant amplitude, $a$, cancels out of the DE.

• Uniform damping, $\mu$, is for simplicity.

• But uniform phase, $\beta$, is realistic.

• Solve for free vibration first (homogeneous solution)

Using a new variable for dimensionless phase, $\varphi =\psi t+\beta$, to simplify the state space gap opening (Eq. (5) simplified in Eq. (11)), the solution starts with a definition of the speed Eq. (13) and acceleration Eq. (14), where coefficients are maintained as variables of time in the derivatives:

Maintain the coefficients as variables of time in the subsequent derivative:

Using Eqs. (13, 14), time derivatives formed below come from application of these differential operators:

From Eq. (10) the resulting dimensionless DE in Eq. (16) below starts to take the form of an SEVP Eq. (17), easily solved by eigenvalue methods. The time derivative operators $\mathrm{\Pi }$ and $\mathrm{\Omega }$ (defined above) allow the vector $\overrightarrow{x}=[x_1\mathrm{}\mathrm{}{x}_2{]}^T$ to factor out. The damping constant, $\mu$, and dimensionless frequency, $\psi$, remain in the calculation until the “nonlinear” eigenvalues and eigenvectors are developed in the formal SEVP solution. A purely imaginary frequency $\mathrm{e}\mathrm{x}\mathrm{p}\left(if\mathrm{\text{'}}t\right)=\mathrm{e}\mathrm{x}\mathrm{p}\left(i\right[i\mu \left]t\right)$ is in the exponential term in $\mathrm{e}\mathrm{x}\mathrm{p}(-\mu t)$ where the sign of $\mu$ indicates damping. The ensuing 1-D analysis provides insight into the dynamics and validates the transient “nonlinear” modes output by the 3-D FEA, modes that exist in reality:

These relations use a different dimensionless stiffness, $\epsilon \equiv {\epsilon }_1=k_{closed}/k_{open}$ where $k_{closed}=k_2$ and $k_{open}\ll k_2$ for numerical stability. Eq. (16) is of the SEVP format as shown below in Eq. (17):

Eq. (18) shows the format of the system matrix $A$ for submission to an eigen-solver. Some of the terms in $A_{\mathrm{1,1}}$ and $A_{\mathrm{2,2}}$ were kind enough to cancel:

Solution of the characteristic polynomial, $\left|\right[A]-\lambda [I\left]\right|=0$, provides the eigenvalues of the system where both DOF have equal displacements. Using a Mathcad’97 symbolic solver 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 Eq. (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\le u\le 1$. The switching control law, $u\left(x\right)$, 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 [23, 50-52].

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.

7.4. Two DOF nonlinear contact eigenvalues

Eigenvalues in Eq. (19) are a dual set for the control law, $u\left(x\right)$, for a gap that opens and closes. 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\left(x\right)$ adds $\epsilon$ to DOF 1 which becomes $1+\epsilon$, as seen in Eqs. (10, 16). This addition uses the FEA maxim that “stiffnesses add” [53]. The total closed stiffness is $1+\epsilon +k$. The derivation of the complete dimensionless “eigenvalues” ${\lambda }_i$ in Eq. (19) are Mathematica results for the solution of the eigenvalue problem described in Eqs. (17, 18) above [7]:

For subsequent nonlinear calculations assume $u=\text{1}$ so that contact is active; there is no gap between the HRA and the hull. Assume $\epsilon >k$ (dimensionless Eq. (16)) and for stability ${\epsilon }_1{k}_1>{\epsilon }_2{k}_2=k_{base\mathrm{}}$ (dimensioned 2 DOF Eq. (8)). Note that all epsilons are dimensionless but they also are a part of the control law that modulates the dimensioned system, Eq. (8).

Fig. 10 plots the real part of the eigenvalues. DOF 1 and 2 represent the hull and armor, respectively, in the 2DOF system. After a critical $\epsilon$, 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 $\epsilon$ regardless of the dimensionless stiffness between the two masses. ($\partial {\lambda }_1/\partial \epsilon \to 0$ and $\partial {\lambda }_1/\partial k\to 0$ but not for ${\lambda }_2$.) From the 3-D FEA results in Fig. 7, the higher frequency clattering mode is the lower energy mode. This is consistent with experience and the Zienkiewicz quote [21, p. 340-341]. The excitation of 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 [44] (the Fourier transform of the impulse response [46]). 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 of St. Venant’s principle [54-56]. 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.

7.5. Eigenvalues: damped 2 DOF sprung mass

Repeated roots occur for the eigenvalues $\lambda (\epsilon ,k)$ plotted in Fig. 10 when the dimensionless closed stiffness is less than a critical stiffness, $\epsilon <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 [29, p. 171], 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\left(t\right)$. 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 oscillates according to $f_{\mathrm{}CG\mathrm{}}^2=k_{\mathrm{}base\mathrm{}}/4{\pi }^2[m_{\mathrm{}hull\mathrm{}}+m_{\mathrm{}MRA\mathrm{}}]$. 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 provides 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 in general $m_{\mathrm{}hull\mathrm{}}\ne m_{\mathrm{}MRA\mathrm{}}$. While deflection at higher frequencies does not modulate the probe beam as much as lower frequencies, this nonlinear eigenvalue behavior has the result that symmetrical modes are better target identification features than the unsymmetrical modes.

7.6. 1-D eigenvectors compared to 3-D FEA

Eqs. (16, 17) introduced the SDOF DE and its 'special' eigenvalue problem (SEVP). In this SEVP the system matrix $\left[A_1\right]$ models a simplified 1-D nonlinear contact system. The derivation of the system matrix $\left[A_1\right]$, from $\left[A\right]$ is in Eqs. (16-18). The solution to the SEVP $\left[A_1\right]\mathrm{\Lambda }=0$ provides two eigenvectors $\overrightarrow{{\mathrm{\Lambda }}_i}$ from $\mathrm{\Lambda }$.

This is the solution of the dynamical DE in Eq. (10). It uses four main assumptions [7] to model the contact point for 1-D as a surface of slideline contacts in NASTRAN: (1) The damping is uniformly constant per unit area and equivalent to the 1-D model, (2) for both the foundation base to the hull as with the hull to the HRA armor. (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, ${\psi }_1={\psi }_2$. The FEA more properly analyzes damping and simulates these results via transient nonlinear time-integrated results using Newmark-Beta methods for algorithmic stability [21, 53], for part of the solutions for the proper 3-D dimensioned DEs for all modes. Therefore, the millions of 3-D DOFs in the FEA removes this last assumption.

In 1-D each of the two DOFs obey the assumed nonlinear solution in Eq. (20) where the amplitude is $x_{\mathrm{}max\mathrm{}}$, which is scaled by a characteristic length $L^{\mathrm{*}}$ providing a dimensionless solution $x\left(t\right)$. $L^{\mathrm{*}}$ is best chosen by using the Buckingham $\pi$ theorem [53, 57], which 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 can be most easily be accomplished using the normal modes along with a few pertinent static deformed shapes:

The single DOF solution in Eq. (20) applies to both DOF's in matrix form, with cross-terms, in the dimensionless system matrix from Eqs. (16-18) as shown in Eq. (21) below:

The dimensionless frequencies, $\psi$, in this system matrix occur in the same form as an eigenvalue $\lambda$ from $\left[A_1\right]=\left|\right[A]-\lambda [I\left]\right|$. In this form the eigenvalues are ${\psi }_1^2={\psi }_2^2={\lambda }_i$, only apply to the above assumptions of (2) equal damping, and (4) equal frequencies for the two masses, which are necessary assumptions 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 ${\psi }_1=\pm \sqrt{{\lambda }_1}$ and ${\psi }_2=\pm \sqrt{{\lambda }_2}$, in a linear combination of modes. Each of the two ${\mathrm{\Lambda }}_i$ modes will have both ${\psi }_1$ and ${\psi }_2$ active for that one mode. Therefore, ${\psi }_1$ and ${\psi }_2$ represent the same energy for ${\lambda }_1$ (symmetric vibration), but ${\psi }_1$ and ${\psi }_2$ produce different energy for ${\lambda }_2$ (unsymmetric vibration).

The FEA eigenvalues (i.e., values of the diagonal NASTRAN $\mathrm{\Lambda }$ matrix) are related to their eigenvectors ${\overrightarrow{u_i}}^T$, as plotted in the thesis [7, Figs. 21-23]. The FEA modal frequencies are the square roots of the eigenvalues, summarized therein [7, Table 7, Fig. 20]. Those frequencies and plots of ${\overrightarrow{u_i}}^T$ represent eigenvalues ${\lambda }_i$ within $\left[\mathrm{\Lambda }\right]$ and eigenvectors ${\mathrm{\Gamma }}_i$ within $\left[\mathrm{\Sigma }\right]$, in analogy to the 2DOF model, where the columns of $\left[\mathrm{\Sigma }\right]$ 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 $\psi$ 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 $\left[\mathrm{\Lambda }\right]=[{\lambda }_1,0;0,{\lambda }_2]$ but it is more convenient to display them in the vector form Eq. (22), where the parameter $\kappa$ was defined in Eq. (19):

The 2 DOF eigenvectors in Eq. (23) are one dimensional mode shapes for displacement of the hull, $u$, and the armor, $v$, for modes 1 and 2:

7.7. Analysis of the 2 DOF SEVP DE

The eigenvalues ${\lambda }_i={\omega }_i^2$ for each column $\mathrm{\Gamma }$ of the eigenvector matrix ($\mathrm{\Sigma }$) are the result of solving Eq. (21). Eq. (24) represents an extreme fixity of this Mathematica solution, a low frequency “DC” limit $\psi =0$:

The parameter $\kappa$, defined in Eq. (19), in this first order correction to the linear eigenvalues (defined for Eq. (22)) is only a function of the stiffnesses, the square of the base and hull stiffnesses. Compared to prior Mathematica results, this $\psi =0$ model provides a “DC” eigenvalue solution that is otherwise not available:

This relation is only valid for the assumptions described earlier: (2) ${\mu }_1={\mu }_2$ and (4) ${\psi }_1={\psi }_2$. Otherwise more nonlinear terms remain and the system is not susceptible to SEVP solution for modes in Eq. (25).

7.8. Synthesis of the 2 DOF SEVP DE

This subsection describes the use of Mathematica to synthesize the eigenvalues and eigenvectors from the system matrix [$A$]. Since the kernel of the SEVP is $\left[A_1\right]=\left|\right[A]-\lambda [I\left]\right|$, and the eigenvalue solutions are the dimensionless frequencies squared ${\lambda }_i={\psi }_i^2$, the approximate system matrix [$A$] derives from $\left[A_1\right]$, as defined for Eq. (22). This is only approximate because of the many combinations of nonlinear and approximately linear variables (e.g. “linear” $\psi$) selected in [39] and used in the “Mathematical Preliminaries” section of [7]. Adding a diagonal of ${\lambda }_i\left[I\right]={\psi }_i^2\left[I\right]$ to the matrix $\left[A_1\right]$ allows the extraction of [$A$] from $\left[A_1\right]=\left|\right[A]-\lambda [I\left]\right|$ as shown in Eq. (26):

Eq. (26) represents a 2 DOF system with springs of stiffnesses $k$ and $\epsilon$, 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 Eq. (27) assuming (${\lambda }_i={\psi }^2\mathrm{}\forall \mathrm{}i\in \left[\mathrm{1,2}\right]$):

7.9. 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 ${\lambda }_1\equiv \lambda$, that the physical argument leading to Eq. (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 $\left({\psi }_k,{\mu }_i\right)\mathrm{}\forall \mathrm{}i\ne k$.

The control law, $u\left({\xi }_1\right)$ of Eq. (28) is the nonlinear stiffness [7, Fig. 11], based on the relative displacement, ${\xi }_1=L^{\mathrm{*}}{x}_1$. The arctangent switch changes the stiffness between the hull lumped point mass and the HRA plate lumped point mass, $k_{open}\leftrightarrow k_{closed}$:

The phase ${\varphi }_k={\psi }_{k}t+\beta \left(t\right)\approx {\psi }_{k}t$ relates to the eigenvector whose temporal dynamics ${\lambda }_i$ describes both DOFs, the hull and HRA point masses. Here the frequencies ${\psi }_k$ are not matched to the damping ${\mu }_i$. Redefining ${\mathrm{\Pi }}_i$ and ${\mathrm{\Omega }}_i$ for this mixed-index format:

Eq. (29) provides an independent identically distributed (iid) estimate of the dynamics modeled as $N$ iid oscillators. Since oscillators for each DOF in an 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 Fig. 9 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.