Highly Sensitive Nonlinear Identification to Track Early Fatigue Signs in Flexible Structures

This study describes a physics-based and data-driven nonlinear system identi ﬁ cation (NSI) approach for detecting early fatigue damage due to vibratory loads. The approach also allows for tracking the evolution of damage in real-time. Nonlinear parameters such as geo- metric stiffness, cubic damping, and phase angle shift can be estimated as a function of fatigue cycles, which are demonstrated experimentally using ﬂ exible aluminum 7075-T6 structures exposed to vibration. NSI is utilized to create and update nonlinear frequency response functions, backbone curves and phase traces to visualize and estimate the structural health. Findings show that the dynamic phase is more sensitive to the evolution of early fatigue damage than nonlinear parameters such as the geometric stiffness and cubic damping parameters. A modi ﬁ ed Carrella – Ewins method is introduced to calculate the backbone from nonlinear signal response, which is in good agreement with the numer- ical and harmonic balance results. The phase tracing method is presented, which appears to detect damage after approximately 40% of fatigue life, while the geometric stiffness and cubic damping parameters are capable of detecting fatigue damage after approximately 50% of the life-cycle. [DOI: 10.1115/1.4052420]


Introduction
It is quite remarkable to witness just in the past decade the technological achievements made to improve fuel efficiency, speed and travel range in modern aerospace and automotive vehicles, and robots. Some of the key enablers of these achievements are recent advancements in light-weight materials [1], morphing structures [2], structural optimization [3], and fluid-solid interaction modeling [4]. While those achievements have been pivotal in generating a new class of advanced structures optimized for performance, these modern structures are becoming highly flexible and thus susceptible to vibration, shock, or impact damage [5]. Consequently, three major engineering challenges have emerged from light structure designs: (i) amplification of dynamic nonlinearities due to high elastic deformation [6]; (ii) limited physics-based identification tools for connecting nonlinearities due to local material degradation to collective global dynamics [7]; thus (iii) detecting early fatigue damage via current diagnostics and health monitoring systems is often formidable. To ameliorate those challenges, significant structural health diagnostic, prognostic, and monitoring technologies are beginning to leverage new advancements in system identification (SI) tools for improving detection sensitivity to early signs of fatigue damage [8,9]. In fact, structural health monitoring [10] and prognostics and health management [11] utilize a variety of SI algorithms for identifying and tracking changes in the baseline structural dynamic features intrinsic to the component of interest. When using SI for damage detection and health monitoring systems, the main objective is to relate anomalies in the baseline structural dynamic features to early signs of damage before reaching a critical failure [12], thus ensuring the safety and extending the service life of high-value engineering systems.
In this section, a brief background on linear and nonlinear system identification (NSI) is provided in Sec. 1.1 with focus on damage detection and monitoring. For completeness, some of the latest NSI approaches for improving detection sensitivity to early fatigue damage are highlighted in Sec. 1.2. The proposed physicsbased NSI approach is described in Sec. 1.3. Sections 2 and 3 detail the analytical development and experimental approach for the proposed fatigue detection approach, respectively. The results and findings of this study are reported and discussed in Sec. 4.

Linear and Nonlinear System
Identification. SI can be defined as an engineering process for modeling dynamic systems mostly based on their input and output signals [13,14]. In aerospace, automotive and robotic applications and structures are often exposed to harsh dynamic environments [15]. Such structures can be thought of as systems with their own intrinsic dynamic parameters. When an excitation is applied to a structure, it is recorded as an input signal and its resulting dynamics as an output signal. The dynamic parameters are then extracted from the structural features embedded in the system output signals. Any changes in the parameters from the established baseline dynamics of the system are leveraged to identify potential damage or failures [16,17]. The estimated dynamic parameters from an output signal can be utilized for detecting instabilities [18], faults [19], or disturbances in the system components [20], and subsequently for estimating their impact on the overall performance [21,22].
In general, the distinction between linear and nonlinear systems can be made by examining how the output signals are mapped mathematically to the inputs. A system is considered to be linear when it obeys the principle of superposition: where a and b are constants, and x and y are input signals. The superposition principle implies that the additive and homogeneity conditions are satisfied, otherwise, the system is considered to be nonlinear [23]. In general, Hooke's law is a good approximation for structures exhibiting linear elastic response due to the application of external loads in the form of restoring forces proportional to the structural displacement, x, as shown in Fig. 1(a). Flexible structures can endure large reversible deformations before their inherent stress-strain characteristics show any significant departure from their linear regime. In fact, kinematic nonlinearity generates dynamic hardening at the first flextural mode, as shown in Fig. 1(b), which is a restoring force proportional to the cubic structural displacement, x 3 [24]. For a beam element, flextural displacement >1.0% than the beam length is considered nonlinear. The total restoring force becomes a function of the linear structural stiffness, k 1 , and nonlinear geometric stiffness, k 3 , as follows: One of the early SI techniques applied for detecting damage is linear dynamic parameters extraction, also referred to as modal analysis [25]. This technique provides insights into dynamic features such as vibration modes [26], internal resonances [27], damping [28], and structural stiffness [21]. Modal analysis can be employed to relate fatigue life to the dynamic parameters of structures subjected to shock [29], harmonic [30], rotational [24], or random excitation [31]. When conducting modal analysis, the effect of nonlinearities are often neglected, due to the assumption of invariant dynamic parameters [32]. However, modal analysis cannot accurately describe the dynamic behaviors of most engineering systems since the parameters are dependent on the magnitude of the applied excitation [33,34]. In reality, most engineering structures exhibit various nonlinear behaviors [35,36]. Engineers are often confronted with geometric and damping nonlinearities in a variety of systems with flexible structural components such as aircraft wings [37], blades [38], composite structures [29,39], electronics [40], robots [41], and micro devices [42,43].
Through NSI, the detection sensitivity to microscopic damage features can be amplified by extracting the nonlinear parameters from the outputs such as geometric stiffness and gyroscopic inertia [32]. Changes in the nonlinear geometric stiffness were related experimentally to material fatigue precursors such as the occurrence of microplasticity [44,45], crystal reorientation [46], and microscopic cracks [3] using techniques such as atomic force microscopy [47], Eddy current [48], and ultrasound [49]. The backbone curves can be leveraged to track the softening or hardening behaviors due to the development of fatigue crack, Fig. 1(b).

System Identification Tools.
It is important to highlight commonly used methods to overcome challenges related to the computational requirement for solving nonlinear dynamic models and high signal-to-noise ratio demand for data processing [20]. The equations of motion can be solved algebraically using techniques such as the harmonic balance method (HBM) [50], incremental HBM [51], method of averaging [52], Volterra series [53], normal form [54], nonlinear normal modes [55], or multiple scales [56,57]. In this study, HBM was applied to analytically solve the equation of motion, where the results were in good agreement with the numerical solution. From a data processing perspective, methods for extracting nonlinear dynamic features are, but not limited to, adaptive and deep learning [58][59][60], data and deep fusion [13,61], cascaded optimization [62], neural networks [63], multi-agent-based modeling [64], and Bayesian dynamic network [65]. Various filtering techniques such as modified Kalman and adaptive particle filters were utilized successfully to improve fatigue predictions [65,66]. The Hilbert transform and wavelet transforms [67][68][69], and zero-crossing [70] were leveraged for constructing the backbone curves and estimating nonlinear damping coefficient.

Physics-Based System Identification Approach.
Our proposed physics-based and data-driven NSI approach for detecting and tracking early signs of fatigue is summarized in Fig. 2. The linear and nonlinear identification processes are shown in gray and orange colors, respectively. The fatigue step, which includes real-time phase and resonance shifts, is displayed in blue (Fig. 2). This approach was demonstrated experimentally and analytically using beam-like structures fabricated from aluminum 7075-T6 aerospace-grade alloy, which was exposed to vibration fatigue. Prior to performing the fatigue test, the baseline linear and nonlinear dynamic features must be established first for each pristine specimen (Fig. 2). The system identification consists of three major steps: (i) identifying the damping coefficients from free vibration, (ii) extracting the stiffness coefficients from slow sine-sweep, and (iii) in situ tracing of the response phase and resonance shifts during vibration fatigue. Free-vibration tests are utilized to calculate the viscous and cubic damping coefficients, c 1 and c 3 , respectively, from the time decay signal [71]. Low and high amplitude slow sine-sweeps are leveraged to extract k 1 and k 3 , the linear and nonlinear stiffness coefficients, respectively [72]. For both freevibration and sine-sweep steps, linear and nonlinear behaviors are defined as tip displacements <1% and >2% of the total beam length, respectively (Fig. 1). The identification steps are repeated after completing each specified number of fatigue cycles (Fig. 2). In this study, the fatigue life was divided into four 25% segments, approximately 3 × 10 5 cycles per segment. The test was terminated when the total resonance shift reached 6 Hz. Details are provided in Sec. 3.
In this study, the sensitivity of each dynamic feature to vibration fatigue was examined as a function of structures' life-cycle. The phase tracing appeared to be the most sensitive feature to early fatigue damage followed by the nonlinear stiffness and damping features. Early fatigue damage is defined as a structure having a remaining useful life >50%. Changes in the linear dynamic features due to structural degradation were insignificant before exhausting 50% of the remaining useful life.

Mathematical Development
The development of the analytical model for approximating the displacement response of flexible cantilevered beams due to vibration fatigue is detailed in this section. Each beam has a uniform cross-sectional area, A, with width b and thickness h along its entire length L, as shown in Fig. 3. The beam distributive mass and the rotary inertia are m s and I r , respectively. The beam is clamped to a rigid support boundary. The vibratory loads are applied to the rigid support in the transverse direction using harmonic excitation with an amplitude Y, as shown in Fig. 3.
The beam displacement is approximated using nonlinear Euler-Bernoulli theory since L/h is ≥ 10 [24]. It is also reasonable to assume that the beam first mode undergoes purely planar flextural vibrations since the structure length to width ratio is ≤ 10 [44]. Therefore, the effects of torsion and shear deformation can be ignored.
Two coordinate systems are used to describe the deformed and undeformed geometries of the beam. The global Cartesian system describing the undeformed geometry is xy. The local ξη orthogonal curvilinear coordinates capture the nonlinear elastic deformation (Fig. 3). Each differential beam element has infinitesimal thickness ds. Exciting the rigid support causes each point on the undeformed cross section of the beam to experience an elastic displacement. The deformation with respect to the x and y axes along the beam's undeformed arclength from the fixed-end to a reference point, s, at time, t, are expressed in terms of a rotational angle, as well as the axial and transverse displacements, θ(s, t), u(s, t), and v(s, t), respectively. The displacement vector can be expressed as Substituting the time derivative of R into the kinetic energy, T, yields the following [24]: where the overdots denote the temporal partial derivatives. The strain energy due to flextural deformation is where E and I are the material elastic constant and second moment of inertia, respectively. The normalized flextural curvature associated with pure bending is ρ. The dissipated energy can be expressed as follows [24]: The total damping force is [37] The linear and nonlinear damping coefficients are c 1 and c 3 , respectively.
To express the nonlinear displacement and its partial derivatives in convenient forms, Taylor series expansion is performed, where only terms up to the cubic order are kept. The axial displacement u is expressed as a function of v since u and v are small but finite [24]. The spatiotemporal energy equations are converted into manageable temporal forms through the application of the Rayleigh-Ritz method using assumed functions for the beam mode shapes [46]: The generalized coordinate q n represents the time modulation of the nth mode expressed by an eigenfunction Ψ n . Only the first flextural mode is considered since pure bending was the beam dominating mode in this study. The assumed function, Eq. (8), is then substituted into Eqs. (4), (5), and (6) to express the temporal forms of the kinetic, potential, and dissipation energies, respectively [24]: Applying Lagrangian mechanics and normalized by the modal mass yield, the final equation of motion [24] is Therefore, the equivalent damping is β eq = β 1 + β 3 q 2 and equivalent stiffness is k eq = k 1 + k 3 q 2 (12)  Transactions of the ASME where d is the longitudinal distance of the material damaged region, starting from the clamped end of the beam (Fig. 3). During the early development of fatigue damage, it is reasonable to assume that the moment of inertia, I, will remain constant. As the damage buildups near the beam root, k 1 and k 3 are assumed to remain constant in the healthy section of the beam; i.e., L − d region, in accordance to the fractional energy theory [7]. E d in a damaged region can be expressed as a function of the damage ratio, D, as follows [7]: To improve the model accuracy, E d is included in the stiffness coefficients to account for the spatially variation in material modulus due to damage accumulation. The steady-state solution of the nonlinear ordinary differential equation (10) is obtained numerically using Matlab ODE solver, ODE45 function. The HBM is also utilized to perform algebraic calculation of the frequency response function (FRF), thus significantly reducing the computation time [34]. The harmonic trial solution is q = Asin(ωt), with an amplitude, A. The excitation function in Eq. (10) can be expressed asŸ = Y/m b sin(ωt − ϕ), where Y is the excitation amplitude. Upon substituting the trial solution and excitation function in Eq. (10), the FRF becomes [34] where the equivalent stiffness is The equivalent mass and natural frequency are m eq and ω n , respectively. The backbone curve (Fig. 1) is Example of the backbone curve for nonlinear beam experiencing a dynamic hardening due to the geometric stiffness effect is shown in the figure.

Experimental Analysis
As stated previously, complete failure was defined as a 6 Hz net decrease in the structure resonance frequency, which turned out to be ∼1.2 × 10 5 cycles. The fatigue life was divided into three cases: 25%, 50%, and 75% of the total number of cycles to complete fatigue failure, i.e., each case consisted of ∼300 × 10 3 cycles. At the end of each fatigue case, the system identification steps shown in Fig. 2 were performed to extract the dynamic parameters and compare them to the beam baseline values. The experimental steps consisted of (i) free vibration, (ii) slow sine-sweep, and (iii) vibration fatigue.

Experimental Setup.
The beam samples utilized in this study were aerospace-grade aluminum alloy, Al7075-T6. The material density and elastic modulus were approximately 2810 kg m −3 and 71 GPa, respectively. The dimensions of each specimen were 150 × 50 × 1 mm 3 . The beam geometries were chosen to provide a generous separation between the first, ω 1 , and second, ω 2 , structural vibration modes to minimize potential modal coupling, where ω 2 must be ≥3ω 1 [32]. In this study, ω 2 ≅ 5ω 1 , thus each beam was always under pure flextural displacement only. Table 1 reports the first three modal frequencies, which were obtained from finite element analysis using ABAQUS [30].
After clamping each beam sample to a rigid steel fixture, the system was attached to an electrodynamic shaker, as shown in Fig. 4. To maintain consistent boundary condition and uniform clamping pressure distribution, the fixture fasteners were torqued  to 10 N m, and checked after completing each fatigue case. The base excitation provided by the shaker was controlled using an accelerometer sensor, PCB Piezotronics Model 333M0. The beam flextural displacement was measured at 30 mm from the clamped end using Polytec OFV-3001 laser vibrometer controller and OFV-303 sensor head (Fig. 4). The reason for locating the laser sensor at the 30 mm location from the fixture was to collect reliably velocity measurements within the laser sensor range limits. Moving the laser beam beyond the 30 mm distance will cause data clipping.

Experimental Procedure.
To generate correlation between the materials degradation and changes in the dynamic parameters, each specimen was tested in two clamping configurations, A and B, as shown in Fig. 5 and highlighted in blue and orange, respectively. The clamping configuration of side A was performed for the 25%, 50%, and 75% fatigue cases. The beams were flipped to side B to perform fatigue testing until complete failure is reached. The example provided in Fig. 6 shows a beam in configuration A fatigued until 6 × 10 5 cycles were reached, then flipped to configuration B to fatigue until complete failure. Side A configuration allowed for performing material degradation analysis at 25% (3 × 10 5 cycles), 50% (6 × 10 5 cycles), and 75% (9 × 10 5 cycles) prior to final fatigue. For side B, material characterization was performed only after complete failure.
The stepping frequency rate for both low and high amplitude sine-sweeps was set to 0.1 Hz. The linear and nonlinear parameters were calculated from the beam measured displacements for low and high amplitude sine-sweeps at 100 mV and 500 mV, respectively. Examples of the linear and nonlinear frequency response and backbone curves prior to vibration fatigue are shown in Fig. 6. Once a set of baseline parameters was established for an undamaged case, the vibration fatigue was performed every 3 × 10 5 cycles (Fig. 2).
The fatigue tests were carried out by setting the excitation frequency equal to the structural resonance obtained from the low amplitude sine-sweep tests. For all fatigue tests, the beam displacement amplitude was kept constant at 1.4 mm, which was measured by the laser sensor and controlled using proportional-integralderivative. With fixed excitation frequency and response amplitude, the phase was traced, as detailed in Sec. 3.3. To ensure constant beam displacement, adjustments were made to the excitation frequency when the response phase exceeded a threshold limit, which was the beam resonance. Deviation from the threshold limit will require supplying the shaker with high electrical power and inducing high stresses on its components.
The example provided in Fig. 6 shows the frequency shift for configurations A and B. In Fig. 6 inset plots, the linear and nonlinear frequency response and backbone curves are shown in black and red marks, respectively. Clear distortions in the frequency response and backbone curves can be noticed when the fatigue exceeds 6 × 10 5 cycles, which is indicative of crack propagation. Additional details are provided in Sec. 4.

Modified Carrella-Ewins Method.
Extracting the backbone curves can be achieved by using the Carrella-Ewins method [73]. Briefly, when applying the Carrella-Ewins method, the dynamic features are extracted by equating two receptance points at the same amplitude at either side of the FRF maximum peak. Therefore, the resonance and damping ratio can be calculated at that specific vibration response level. Subsequently, a discrete function describing the backbone and damping curves can be constructed by repeating this process from the lowest to highest response amplitudes. The shortfall of the Carrella-Ewins method, however, its assumption that both receptance branches are available, is not true for manifold points; notice the red curve in the Nyquist circle in Fig. 7. The Nyquist circles of the same dynamical system with its linear and nonlinear parameters are shown in black and red markers, respectively (Fig. 7). To overcome this limitation, our modified Carrella-Ewins method, instead, takes three mobility points (Fig. 7): two low amplitudes at either side of the FRF peak, shown in green markers, and a third point anywhere on the Nyquist plot, shown in a blue marker. Therefore, three is the minimum number of points for defining an FRF circle. The mobility points can be obtained as follows: The three points are utilized to calculate the resonance frequency and linear damping. Keeping the same two low amplitudes points (nonlinearity is not active) and sweeping the third point, several triplets can be created each of which identifies a new linear Nyquist circle relative to the amplitude of the sweeping frequency point. Taking advantage of the Nyquist circularity, the linear damping is calculated from the diameter inverse. In Fig. 7, for example, two linear and nonlinear Nyquist circles are shown in black and red traces, respectively. When examining the linear response, it can be observed that every three frequency points obtained from the Nyquist data set always resolve the same circle. For linear response, all the frequency points are equally spaced within the circle's center. When the response is nonlinear, every three mobility points describe a new circle. Therefore, the Nyquist circle parameters are dependent on the magnitude of each selected point.
While there are several approaches for performing a circle fit or inverse FRF method, the Dobson approach was selected in this paper [74]. Briefly, the Dobson approach uses FRF points to build a series of straight lines each of which is curve-fitted. Both the angular coefficient and intercept are collected for each straight line to create two new straight lines. The Dobson approach is   3.4 Phase Tracking Method. The phase was defined as the lag between the input force and response singles over the number of excitation cycles. The evolution in the phase was a more sensitive fatigue indicator than shift in the resonance (Fig. 8). Similar findings were reported for delamination onset in composites [75,76]. Utilizing the finite element method, Voudouris et al. showed that vibration response phase and delamination propagation were linearly correlated [77].
In this study, we borrow experimental strategy for composites reported in Refs. [75,76] to track the phase for aluminum samples exposed to vibration fatigue. In Fig. 9, the phase plotted as segmented slopes due to adjustments in the excitation frequency and details are provided in Sec. 3.2. For enhanced visualization of the response phase and its sensitivity to early fatigue damage, only phase traces between 4.8 × 10 5 and 7.2 × 10 5 cycles (Fig. 9).

Results and Discussion
This section is divided into two subsections. In Sec. 4.1, the main focus is on reporting and discussing the dynamic response of the tests samples before, during, and after exposure to high cycle fatigue. Both experimental and modeling results, in Sec. 4.1, consist of the FRF and backbone topology of pristine and fatigued cases. Analysis of in situ frequency and phase traces during fatigue testing is included. The comparison study between the experimental and analytical results are discussed in detail. Finally, Sec. 4.2 provides detailed insights into the evolution of the dynamic parameters as a function of fatigue cycles, as well as their sensitivity to material degradation over time.

In Situ Identification Using the Phase Tracking
Method. It can be observed that the evolution trends in the phase traces for configurations A and B are similar (Fig. 8). This is an indication that flipping the beam from side A to B did not have any significant effect on the dynamic behavior of the beam. Overall, each phase trace is shown as a series of segmented lines plotted at their associated excitation frequencies, which were not updated until approximately 5 × 10 5 cycles.
While the emerging segmented phase slopes were indicative of early crack development, it was not clear why there was a slight jump in the phase between approximately 5 × 10 5 and 6 × 10 5 cycles, prior to updating the excitation frequency (Fig. 8). The phase increased after 6 × 10 5 cycles, followed by a drop after 9 × 10 5 cycles (Fig. 9). A potential explanation is the phase that appears to indicate material hardening followed by softening due to fatigue accumulation-a phenomenon reported in a previous study for Al 7075-T5 [47]. Insights into the relationship between the material hardening-softening phenomenon and phase shift is beyond the scope of this investigation and the subject of future studies.
Finally, calculating the phase change over the number of cycles Δϕ/ΔN, described in Sec. 3.4, one can estimate the crack growth rate by means of phase evolution instead of the crack tip length. In Fig. 10, Δϕ/ΔN is reported for all three specimens, which illustrates the onset damage rate over the number of cycles. The rate of damage growth indicated by Δϕ/ΔN appears to be consistent for all three samples (Fig. 10), which means the method is repeatable. The data scatter was due to the manual collection of Δϕ/ΔN. For the sake of completeness, crack growth in each specimen within each fatigue was verified, qualitative, using optical microscopy.

Evolution in the Dynamic Parameters Due to Fatigue.
An example of the dynamic parameter identification for a beam experiencing fatigue is provided in Fig. 11, where the FRF and backbone curves are shown prior to fatigue testing, in Figs. 11(a) and 11(b), after 300 × 10 3 fatigue cycles in Figs. 11(c) and 11(d), 600 × 10 3 cycles in Figs. 11(e) and 11( f ), and 900 × 10 3 cycles in Figs. 11(g) and 11(h), respectively. The experimental and modeling results are shown in markers and sold lines, respectively. The numerical solutions for the equation of motion, Eq. (10), and HBM, Eq. (16), are the blue and red curves, respectively, superimposed on top of the experimental results. The backbones were calculated from the experimental data using the modified Carrella-Ewins method, Eq. (17), and compared to the HBM (Figs. 11(b), 11(d ), 11(f ), and 11(h)). The amplitudes in Fig. 11 are reported as the maximum velocity normalized by the input acceleration.
The input frequencies were normalized by the natural frequency of the pristine specimens to show that the linear and nonlinear backbone curves overlapped with excellent agreements up to 600 × 10 3 cycles. Near 300 × 10 3 cycles, the fatigue damage led to distortion in the nonlinear response, where the maximum amplitude increased from approximately 0.82 m/s/g to 0.83 m/s/g (Fig. 11(c)). The insignificant drop in the natural frequency was approximately 0.2%. After completing ∼ 600 × 10 3 fatigue cycles, the FRF distortion remained unchanged, but the natural frequency dropped by approximately 0.8% (Figs. 11(e) and 11( f )). For the 900 × 10 3 cycles case, the linear FRF curve shows dynamic softening manifested by material compliance due to crack development [44,48]. At this point of the fatigue process, breathing-crack was developing [6].
The dynamic softening seen in the response signal was also captured analytically by updating the parameters k 1 , k 3 , c 1 , and c 3 , as discussed in Sec. 2. Using the sine-sweep tests, the linear stiffness, k 1 , was calculated from the backbone curves. The decrease in k 1 as a function of cycles is shown in Fig. 12. The linear and nonlinear damping coefficients, c 1 and c 3 , respectively, were obtained from the free-vibration experiments and reported in Fig. 13 as a function of fatigue cycles. After obtaining k 1 , c 1 , and c 3 , the nonlinear geometric stiffness, k 3 , in Eqs. (10) and (16) was to match the experimental results (Fig. 11).
Updating parameters k 1 and k 3 was achieved by adjusting the damage ratio, D, Eq. (13), up to 600 × 10 3 cycles. Due to the FRF distortion severity in the 900 × 10 3 cycle case, D was adjusted only for k 1 . However, k 3 was adjusted numerically to curve fit the linear and nonlinear FRFs. The decay in k 1 and k 3 parameters as a function of fatigue cycles is reported in Fig. 12. Both k 1 and k 3 were normalized by their respective initial values. Both k 1 and k 3 decreased logarithmically due to the increase in number of fatigue cycles. The logarithmic degradation in k 3 , however, was more severe than that for k 1 after 600 × 10 3 cycles, which was consistent with fatigue studies by Andreaus and Casini [67] and Newman [78]. Figure 13 shows the beam displayed logarithmic decrease in c 1 . However, c 3 increased exponentially as a function of fatigue cycles (Fig. 13). It is important to point out that the standard deviations of the calculated nonlinear parameters were higher than those for linear cases.

Conclusion
The paper illustrated that NSI holds promise as a viable tool for early fatigue damage diagnostics and monitoring. NSI was utilized for detecting and tracking fatigue in flexible aluminum structures exposed to vibratory loads. The nonlinear geometric stiffness, cubic damping, and phase shift were tracked as a function of fatigue cycles. Extracting the stiffness and damping coefficients required interrupting the fatigue tests. However, in situ phase tracing was accomplished successfully without interruptions using the proposed modified Carrella-Ewins method. The manifestation of material degradation due to fatigue accumulation can be identified from the (i) logarithmic decrease in the rate of change in phase and (ii) increase in the cubic damping accompanied with a decrease in the nonlinear geometric stiffness. It is important to point out that linear parameter identification was an inadequate technique in detecting early signs of fatigue.
Tracking the rate of changes in phase as a function of fatigue cycles appeared to be even more sensitive in detecting early signs of fatigue than any of the nonlinear parameters. Therefore, further studies are needed to gain deeper understanding of the physics required to exploit those parameters in a variety of materials and geometries to improve our capabilities for localizing, quantifying, and classifying early fatigue damage.