Non-linear dynamics of a stiffened composite laminated panel with debonds

.


Introduction
Composite laminated panels with stiffeners are often applied in engineering, due to their large resistance to buckling, high specific stiffness and strength [1][2][3][4].Flaws, which can be produced during the manufacturing process or the service life [5], compromise the integrity of these structures.In the framework of the Optimised Design for Inspection (ODIN) action [6] of the European Cooperation in Science and Technology (COST) [7], a structure that represents an airplane wing component was designed, manufactured and tested [8].Numerical and experimental analyses of this structure indicate that the bonded connections between flat composite panels and aluminium stiffeners are critical regions [8,9].One of these laminated composite panel reinforced with partially debonded aluminium stiffeners is numerically analysed in the present paper.It is intended to understand if the non-linear dynamic behaviour due to this type of flaw is useful for vibration-based structural health monitoring (VSHM).
VSHM technologies have the potential to detect flaws and enable online inspection [10][11][12].Some VSHM methods are based on the fact that the local reduction of stiffness caused by damage modifies modal parameters [13][14][15][16].However, the natural frequencies of vibration may lack sensitivity to damage [11,[17][18][19], since damage starts as a local phenomenon that does not significantly affect global stiffness.According to a review by Tsyfansky and Beresnevich [20], a damaged area of 15-20% on the cross section of a component may lead to a reduction of the natural frequencies smaller than 1.5%.The mode shapes of vibration are more sensitive to damage than natural frequencies, but are more difficult to obtain due to the need for many sensors.Progress has been made to solve this problem by using contactless instrumentation, such as a scanning laser Doppler vibrometer [21].Alternative methods, that do not rely directly on modal data, have been developed, as in [22], where the mutual information between two signals measured on a vibrating structure is suggested as a delamination detection method.
Response-based methods are generally based on the comparison of measured vibration signals from the undamaged and the damaged structures.If the undamaged system is linear, the undamaged results may not even be needed.In comparison with modal-based methods, response-based methods have the additional benefit of being more sensitive to small flaws and less affected by noise, operational and environmental conditions [15,24].Tools to analyse motions [26] that can be used to detect damage are addressed in the next paragraphs.
Time histories of the response of undamaged and damaged structures can be compared.Alternatively, the response can be compared with the input signal, using the fact that flaw induced non-linearities lead to the appearance of harmonics in the response to harmonic inputs [14].Displacements are frequently shown in numerical analysis, but accelerations, which are more sensitive than displacements to the presence of harmonics, because the displacement amplitude of each harmonic is multiplied by the square of the respective frequency [25], are often measured in experimental testing and can be advantageously employed.
Trajectory projections in a phase plane, typically in a displacement versus velocity graphic, can also be used.Phase space trajectories can be complex but provide relevant information.Reaching steady state, the trajectories converge to an invariant subspace, the attractor, which is easier to characterise than the full response.A comparison between the results of the undamaged and the damaged structure allows conclusions to be drawn [24]; invariants of the attractors, as Lyapunov exponents and state space dimensions can be employed.In statistical approaches, the statistical distribution of the attractor's points is used to extract damage sensitive features.For instance, the skewness and the variance have shown sensitivity to damage [27,28].
A Poincaré section is a projection of the phase space response in a lower dimensional space.Often, the Poincaré section contains a sequence of points at moments in time t, equally spaced, with t a multiple of the forcing periodic function's period T. If the response is characterised by a period nT, the section contains n different points, if it is quasi-periodic, it defines a line and if it is chaotic an infinite number of points [14,29].Poincaré sections can be quite sensitive to changes in the system and contain far fewer data than phase plane plots or time histories.To profit from this, Manoach and co-workers developed damage indexes based on Poincaré sections [21,24,30].
Finite element models provide a means to test the applicability of structural health monitoring methods to specific structures.Cho and Kim [31] applied a higher-order zig-zag theory to analyse the effect of delamination on natural frequencies and buckling loads.Hu et al. [32], employed a higher-order shear deformation theory to investigate how modal properties are affected by delamination.Kim et al. [33] resorted to a layerwise theory for a similar purpose.Burlayenko and Sadowski [34][35][36][37][38][39] analysed in detail linear and non-linear vibrations of sandwich plates with debonded layers using Abaqus®.Natural frequencies, frequency response curves, phase plane plots and Poincaré sections were extracted from the time history signals.In [40] and [41] non-linearities due to large displacements and to contact were taken into account.
In addition to [8,9], several experimental and numerical studies have found that stiffened panels are prone to debonding between the stiffener and the panel, e.g.[15,[42][43][44][45][46][47].In [42] a fracture mechanics analysis of skin-stiffener debonding is carried out using shell finite elements.In [44] numerical and experimental investigations are performed on the detection of defects in composite T-stiffened panels using vibration modal analysis.The mode shape curvature was found to be a more useful damage indicator than the actual mode shapes.An experimental study on the damage evolution under fatigue load and on the shear-afterimpact-fatigue behaviour of stiffened composite panels is presented in [45].One of the main types of failure found was debonding between the stiffeners and the skin.In [48], a semi-analytical approach to analyse the consequences of skin-stiffener bonding flaws on the linear buckling and vibration behaviours of T-stiffened is presented.Large debonding defects significantly affect the critical linear buckling load and the first natural frequency, but the same is not true for debondings of smaller areas.The progressive separation of the skin and stringer in stiffened CFRP panels under compression is analysed in [46].Tan et al. [47] also found skin/stiffener interface debonding in their investigation of Tstiffened composite panels subjected to impacts.Low-velocity impact tests and compression after impact tests were conducted on a composite stiffened panel in [49]; debonding between the stiffener and the skin panel was one of the failure mechanisms found.Carminati and Ricci [15] analyse an aluminium rectangular plate with five T-shaped stiffeners reinforcing the longitudinal direction.Three damaged scenarios are recreated by the absence of connections between one stiffener and the base plate.
To the best of the authors' knowledge, a published analysis of the non-linear dynamics of stiffened plate with partially debonded stiffeners the goal of this workdoes not exist.
The structure analysed in this paper is presented in Section 2. The theoretical assumptions adopted in the implementation of the finite element model in Abaqus® are justified in Section 3.Then, in Section 4, modal and response-based analysis are applied, namely resorting to natural frequencies, mode shapes, time histories, phase plane plots, Poincaré sections and frequency spectra.Section 5 concludes the paper with a summary of the findings that can be useful to detect damage in stiffened composite panels.

Structure to analyse
The panel analysed is a rectangular plate, reinforced with T shaped stiffeners as represented in Fig. 1.The plate has 18 stacked plies (-45 ) s of transversely isotropic carbon fibre-epoxy with reference USN150B [50]; all plies have the same thickness.The five stiffeners are in aluminium 6063A T6.Araldite® 420 A/B is used to bond the stiffeners to the panels.
A linear elastic constitutive material relation is assumed.The elastic properties and density of the aluminium alloy and of Araldite® 420 A/B  are summarised in Table 1.E stands for Young's modulus, ν for Poisson's ratio and ρ for mass density.The properties of aluminium are from [51] and the ones of Araldite® 420 A/B are from [52] and [53].Both materials are considered isotropic.Not all the necessary properties of USN150B were found in manufacturer datasheets; however, a detailed datasheet of a very similar composite, USN150A, was.Knowing that USN150A and USN150B only differ in the fibre content, the missing mechanical properties of USN150B were inferred using micromechanics-based models.Details on these deductions are provided in [9]; the values of the diverse properties are given in Table 2. E ii and G ij represent, respectively, the Young's modulus and shear modulus related to directions i and j; these are: the fibre direction, 1; the direction perpendicular to the fibre, in-plane (2); the direction perpendicular to the fibre, out-of-plane (3).These properties were compared with the existing literature and with tests carried out in the framework of COST action 18,203 [6,9].
Numerical and experimental static analysis were carried out on the whole representative wing structure and on the composite panel with stiffeners, with the goal of identifying regions of the panel where damage is likely to occur [8,9].These analyses indicate that the adhesive connection between the composite plate and the stiffener is a critical region, due to debonding at the Araldite® 420 A/B layer.The fragility of this region is due to the large mismatch bending that results from the large difference in the stiffness of each part.Thus, although other types of damage may appear, depending on loads and boundary conditions, we will here focus on debonding between the composite plate and the stiffeners.

Finite element model and solution of equations of motion
Abaqus®/Standard code [54] was used to represent the stiffened panel.The composite plate was discretised by reduced integrated 8-node continuum shell elements SC8R, the stiffeners by 8-node hexahedral continuum solid elements with incompatible modes C3D8I.The latter elements have the advantage of eliminating shear locking and alleviating volumetric locking.The meshes are conformal.Reference [54] provides more information on these finite elements, as well as on the methods employed to consider contact and to solve the equations of motion, which are briefly described in the following paragraphs.
Since we are interested in understanding the differences introduced by debonding in the dynamic behaviour of stiffness panels, it is important to consider contact between debonded surfaces.A brief description of the contact kinematics and contact conditions of interest to the present analysis is provided here.A comprehensive theoretical background is given in [55][56][57].Damage progress is not addressed because it is not part of the objectives of this analysis.
A typical general formulation is employed, where contact interactions between two deformable bodies, Ω s , Ω m ∈ R 3 , Fig. 2, are assumed.Index m indicates the so-called master body, whilst index s corresponds to the slave body.
The boundaries of the two bodies, ∂Ω i , i = m,s, are divided into three disjoint subsets where Γ u,i stands for the Dirichlet boundary where the displacements u are imposed, Γ σ,i for the Neumann boundary where the tractions t are applied and Γ c,i for the potential contact surface.A characteristic of contact problems is that the active contact surface Γ a,i ∈ Γ c,i is unknown, may change over time and thus must be determined as a part of the solution process.
To measure the proximity, potential contact, and penetration of the two bodies, a geometric measure named gap function is introduced.It can be defined as In it, xm is the contact point on the master surface associated with each point x s of the slave surface and n c is its contact normal unit vector.The identification of the contact point on the master surface may be performed by a closest point projection The argminf(x) function provides the input x that results in the minimum of f(x).So, in this context, the previous equation returns the point x m , from the potential contact master surface, closest to point x s of the slave surface.The Euclidean norm is employed in (3).
Neglecting any tangential interactions in the contact (frictionless contact), the contact traction on the slave surface t c,s can be written by where p n is the normal contact pressure.Due to the action-reaction principle on the contact interface, the contact traction on the master surface is t c,m = − t c,s .
For the normal contact constraints, the Karush-Kuhn-Tucker (KKT) conditions are used These define a non-adhesive and frictionless contact.Equation ( 5) represents the geometric constraint of non-penetration and Equation (6) implies that only compressive stresses are allowed in the contact zone.In conjunction with the former equations, Equation (7) forces the gap to be closed when contact pressure is not zero (contact occurs) and the contact pressure to be zero when there is no contact (gap is open); it also allows for the transition case where both the pressure and the gap are null.In the context of contact dynamics, the following persistency condition must be considered where ġn ( x 0 , t ) is the material time derivative of the gap function.In the present model, the constraints are enforced using the Lagrange multiplier method [57], with the Lagrange multipliers forming vector λ.
Using the interpolations of displacements [9] and assembling the contributions from all elements, the final equation of motion can be written as where M is the global mass matrix, F int is the global vector of internal elastic forces, F ext is the global vector of external forces and F c is the global vector of the contact forces.d and d are global vectors of nodal values of accelerations and displacements (generalised coordinates), respectively.
To take energy dissipation into account, a global damping matrix C is introduced in Equation ( 9) where ḋ is the global vector with nodal values of velocities.As the other matrices, the global damping matrix can be obtained by assembling the contribution from each elemental damping matrix C e .Using a Rayleightype damping, these are a linear combination of the elemental mass matrix M e and the elemental initial tangent stiffness matrix [58] for each element where α e and β e are proportional damping constants of each element.
Equation (10) denotes a system of n nod × n e non-linear ordinary differential equations.These must be subjected to initial conditions, to boundary conditions on Γ u,h and to impenetrability conditions, equations ( 5)-( 7), on Γ c,h , with subscript h standing for 'finite element approximation'.
The discretisation of the contact interface is accomplished by the contact algorithm, responsible for the discretisation of the virtual work and for the enforcement of the contact constraints.In this paper, a surface-to-surface contact algorithm was applied.Compared to node-tosurface discretisation, surface-to-surface discretisation provides more accurate stress and pressure results and is less sensitive to master and slave surface designations.However, it has a larger computational cost [54].
The Hilber-Hughes-Taylor (HHT) temporal integrator [59] is employed so solve the equations of motion.This implicit scheme is second order accurate, unconditionally stable and with controllable numerical damping.The α -dissipation is introduced in the integration by shifting the evaluation of the internal, damping and contact forces to a generalised midpoint t n+αHHT , with α HHT ≤ 0, so that Equation ( 10) is replaced by [37] with the initial conditions d 0 = d and ḋ0 = ḋ, given displacement boundary conditions and boundary conditions updated due to developing contact.

Introduction
In this section, dynamic analyses are performed on the healthy and damaged reinforced composite panel.First, section 4.2, the modes of vibration are obtained.Then, damping parameters are established in section 4.3.Lastly, section 4.4, the influences of the flaw size, the measurement location and the excitation frequency on non-linear dynamics of the structure are evaluated.

Modal analysis
The flaw is inserted in the central region of the structure, between the composite panel and the stiffener.This is implemented by not tying this region as presented in Fig. 3; contact is not considered in this section, so this analysis only provides an indication to the effect of debonding in the modes of vibration.
Flaws are quantified by The boundary conditions are as in Fig. 4, where the structure is sustained on two roller supports along its width and has pinned points in the middle of edges parallel to the x and y directions to avoid a rigid body motion in the y and x direction, respectively.Lanczos method [54,60] is applied to compute the eigenvalues.
Table 3 lists the first six natural frequencies for simulations with

Table 3
First six natural frequencies obtained for different flaw sizes.
The natural frequencies of the undamaged (γ = 0%) structure are the largest, because flaws cause a loss of the structure's stiffness.However, even the largest of the flaws considered has an extremely small influence on most natural frequencies of vibration.The difference between the values of the first natural frequency with a flaw defined by γ = 8% and without flaw is only 0.031%.With lower sized flaws, the values of some natural frequencies of vibration are, at least up to the 5th digit, equal to the ones of the undamaged structure.In engineering practice, it would be nearly impossible to detect flaws from the majority of the natural frequencies on Table 3.
Only the third natural frequency of the structure with γ = 8% differs more than 1% from the respective frequency of the undamaged structure.The reason for this distinction is that larger displacement amplitudes are attained in the damaged region when the stiffened panel vibrates in the third mode of vibration, Fig. 5 (this figure only contains the mode shapes of vibration of the undamaged structure; the mode shapes obtained for the damaged structure are similar).This indicates that one can use the combined information from natural frequencies and natural mode shapes of vibration to estimate damage location.Still, the third natural frequencies of less damaged structures (from γ = 0.4% to γ = 4%) differ very little or not at all (within five digits) from the ones of the undamaged panel, being the difference lower than 0.15%.

Damping parameters
In the absence of experimental modal analysis on the structure, values determined from experiments on similar materials were used to establish a first approach for proportional damping constants α and β of equation (11).The initial values of α and β for the USN150B composite plate come from experimental modal analysis on a carbon-epoxy composite plate [61].The ones for the stiffeners in aluminium 6063AT6 are calculated using experimental data from an aluminium beam [62].We are here interested in analysing steady-state oscillations and increasing damping decreases the computational time required to achieve these.Therefore, after verifying in a few simulations that the essential outcomes of this work would not be affected by such decision, the values of the proportional damping constants were doubled.The parameters employed are listed in Table 4.

. Load, boundary conditions and solver parameters
Now, a harmonic load is applied in the central region of the composite panel -Fig.6.The amplitude is 50 N, the frequency is defined in each of the following sub-sections.The boundary conditions are the ones of the modal analysis (Fig. 4).
The parameters used in the HHT time integrator, Table 5, are the

Table 5
Parameters for the HHT integrator used in the dynamic analysis of the structure.ones suggested by Abaqus® for contact situations, benefitting from "numerical dissipation to damp out any spurious participation of the higher modes" [54,59].Geometrical non-linearities are neglected, because large displacements are not achieved.
Convergence tests were carried out to estimate adequate time increments, Δt.Different numbers of points per excitation wave were considered and applied in the following equation: Then, a Δt that appeared to be smaller than necessary was employed, to assure accuracy.In the undamaged structure case, 40 points per excitation cycle provided an accurate and smooth description of the response, with a reasonable simulation time.Contact between surfaces  of the damaged reinforced panel introduces higher-order harmonics Hence, the number of points per cycle used in each time simulation increases with the flaw size: 120 points for γ = 0.4%, 160 points for γ = 2%, 200 points for γ = 4% and 240 points for γ = 8%.Furthermore, it was verified that the duration of each simulation is enough for the dynamic response to reach the steady-state stage.

Numerical results of the dynamic analysis on the structure -damage size influence
Multiple simulations with different flaws, characterised by parameter γ, Equation ( 13), were run: γ = 0.4%, γ = 2%, γ = 4% and γ = 8%.Contact interactions were considered between the composite panel and the central stiffener where a flaw exists -Fig.3. The simulation hence mimics interaction between the detached surfaces, without interpenetration.Particular interest will be paid to the displacements and velocities in two points in the centre of the detached area.The latter two points are represented by black dots in Fig. 3.
Figs. 7-10 display the time histories, phase plane plots, frequency spectra and Poincaré sections in the top central point of the composite panel.The damaged structure plots are represented together with the ones from the undamaged structure.
Several conclusions may be drawn from the time histories (Fig. 7): irrespective of the flaw size, the response remains periodic with the period of the excitation; damage introduces an asymmetry in the response, in relation to the static equilibrium configuration; this occurs because when displacements are positive the plate and the stiffener separate and the displacement of the plate achieves large magnitudes; in the part of the cycle where displacements are negative, the behaviours of the damaged and undamaged structure are quite similar; the plate and the stiffener appear -we will see in other analyses that it is not always so -to remain in contact, as if the flaw was not there; damage leads to responses that are not harmonic, as more evidently happens when γ = 8%; the two previous findings are due to the non-linearity of the damaged structure.The greater the flaw, the more important the non-linearity, increasing the difference between the responses.Hence, with smaller damaged areas it is only possible to perceive the asymmetry of the response, but for larger damaged areas it is also noticeable that the wave loses its sinusoidal shape.
The phase plane trajectory (Fig. 8) corresponding to the undamaged structure is a perfect ellipse, centred at the origin: a linear response.With the introduction of damage, the biggest differences in the phase plane plot occur for positive displacements due to the discrepant stiffnesses of the components, which detach in part of the cycle.In some phase plots, differences between the undamaged and damaged structure in the negative displacement part of the cycle are also visible.In comparison to the displacement time histories, the introduction of velocity in the analysis amplifies the importance of higher harmonics, turning flaws more conspicuous.
Regarding the importance of the flaw size in the phase plane plots, for a small damaged area (here represented by the γ = 0.4% case) the difference to the undamaged results is visible but discreet.For larger damaged areas, the phase plane plot becomes more complex and the asymmetry with respect to the axis defined by null velocity increases.The asymmetry occurs because when velocities are positive the two parts of the structure detach, whilst when the velocity is negative, the two parts collide and then reattach or very nearly so.The impacts lead to abrupt variations in velocity, observable in all phase plane plots but which become more patent as the flaw size increases.
To obtain the frequency spectra of Fig. 9, the time histories calculated are processed by using the Fast Fourier Transform (FFT) within the MATLAB® software environment.The frequency axis is normalised, by dividing frequencies by the excitation frequency ω = 500rad.s− 1 .
The response of the undamaged structure is harmonic.In the simulations with damage, the fundamental harmonic is still the most important, but higher-order harmonics that are multiples of the excitation frequency exist.As the flawed area increases, so does the contribution of higher harmonics.These have two physical explanations: (1) the variation of stiffness along a cycle of vibration, as detachment/ reattachment occurs; (2) the collisions between the detached surfaces.
In the γ = 0.4% simulation the amplitudes of higher-order harmonics are very small.In the γ = 2% and γ = 4% simulations, higher-order harmonics start to stand out, especially the 2nd one.Finally, in the γ = 8% simulation, the contributions of the 2nd and 4th harmonics are important, but many other harmonics are visible.The appearance of even harmonics in the response is consistent with theflaw inducedabsence of symmetry in this dynamic system.
Lastly, the Poincaré sections are constructed, Fig. 10, by registering the displacement and velocity near the maximum positive displacement, separated by period T [26] Since the periodicity of all the responses is one (i.e., the period is the one of the excitation) only one point is visible in the images.The location of points in each Poincaré section increasingly differs as damage increases.This can be used to establish damage indexes [21,30].
Fig. 11 displays the displacement of the two points represented by black dots in Fig. 3, for a damage of γ = 8%.The contact constraints are properly implemented, since there is no interpenetration.The detached parts come into contact once every excitation period and remain in contact for about half a cycle, when the displacement is negative.The thin composite plate is more flexible than the aluminium stiffener: when the displacement is negative, the rigid stiffener limits the displacement of the composite panel, so the two parts tend to reattach.

Response location influence
In this section, the previous analysis with a flaw of γ = 8% is replicated, but outputs are registered at twelve distinct points, shown in Fig. 12.It is intended to inspect how response signals change with location.Six points are in the line of symmetry and the other six lie 75 mm to the side.Points number 2, 4, 6, 8, 10 and 12 are located above the stiffeners; the others are between them.
In Fig. 13 and Fig. 14, the time history of each of the twelve points is represented.From these figures, some important observations can be highlighted: the fact that the structure is damaged is easily detected in all points, because the oscillations are not harmonic; in the undamaged structure, the displacements of points further away from the loading point attain lower magnitudes; in the damaged structure an exception to this apparent rule occurs in point 1, with a stronger excitation of at least one higher harmonic; higher-order harmonics are dominant in the responses of points 1, 2,  damage affects the response of point 12 much less than the response of other points, a reminder that measurements should always be performed in more than one point.
To better illustrate the harmonic content of these responses, Fig. 15 contains the frequency spectra of the responses of points 1 and 4. The spectra extend the findings of Fig. 9 to the other locations of the panel: the interaction between the detached parts leads to higher-order harmonics.It is worthy of note that higher-order harmonics stand out in the spectrum of point 1 (as it does in the spectra of points 2, 3, 7, 8 and 9, not shown), because the contribution of the fundamental harmonic is small.
Poincaré sections of the response in these points were also defined [9], but are not shown here, for conciseness sake.As transpires from the  phase plane plots and frequency spectra, the responses are always periodic with the period of the excitation; therefore, a single point was obtained in each Poincaré section.Debonding led to a change in the location of points in Poincaré sections, qualitatively similar to what can be seen in Fig. 9.

Excitation frequency influence
Simulations were carried out at two additional excitation frequencies: ω = 1000 rad⋅s − 1 and ω = 650 rad⋅s − 1 .The former is closer to the first resonance frequency, where the response amplitude increases.
The latter is about half the first resonance frequency.This case is interesting because, as verified in the previous sections, the damageinduced non-linearity introduces the second harmonic in the response and the second harmonic of 650 rad⋅s − 1 is close to the first natural frequency (Table 3).An excitation frequency equal to the first natural frequency was also attempted, but this test case was not concluded, because it was computationally challenging to achieve steady state by numerical integration so near resonance, in this lightly damped, multidegree-of-freedom, non-linear system.
Convergence with the number of discretisation points per wave was again verified [9].
Excitation frequency ω = 1000 rad⋅s − 1 .Fig. 16 displays the time histories of the response without, γ = 0%, and with damage, γ = 8%, excitation frequency ω = 1000 rad⋅s − 1 .In the flaw case, the time histories of the detached parts, i.e. top and bottom, are shown.Here and in the following figure, Fig. 17, results are from the central point of the composite panel (Fig. 3).
The non-linearities introduced by the flaw clearly show in the response.As with the lower excitation frequency, the intact structure and the detached parts oscillate with the period of the excitation.
Because the response amplitude is larger for ω = 1000 rad • s − 1 than for ω = 500 rad • s − 1 , damage leads to a greater variation in the response when ω = 1000 rad • s − 1 .Furthermore, now, even when the displacement is negative, the response is not the one of the undamaged structure.
The frequency spectrum of the displacement when ω = 1000 rad • s − 1 confirms that the most important harmonics are the first, the second and the fourth [9].Fig. 17 represents the phase plane plot corresponding to ω = 1000 rad • s − 1 .The difference between the response with and without damage is plain.Not only the plot of the damaged response is not symmetric with respect to any of the reference axes, but also it is very far from elliptical, confirming that many harmonics are excited.Both when the two parts of the damaged structure detach (positive velocity) and reattach (negative velocity), spikes appear in the phase plot.These are very pronounced in the reattachment case, due to the impacts.
It is concluded that, due to the larger response magnitudes attained, an excitation frequency closer to resonance results in a more noticeable effect of flaws in the response than a frequency that has no relation with resonance frequencies.
Excitation frequency ω = 650 rad⋅s − 1 .Fig. 18 shows time histories of the responses at diverse points of the structure with, γ = 8%, and without damage, when the excitation frequency is approximately half the first resonance frequency.
The difference between the dynamics of the sound and damaged stiffened panels is remarkable, turning this the most interesting of the excitation frequencies attempted.Higher harmonics are important in the displacement of all points.As occurred with the other excitation frequencies, but more pronounced now, debonding turns the magnitude of the response larger.In points 1 and 7, where the first mode shape attains large amplitudes, oscillations barely occur without damage, but they became large with debonding; this impact of debonding was not so notorious with the previously attempted excitation frequencies.
The phase plots -Fig.19 -confirm the remarkable (and larger than with the other excitation frequencies) impact of the debond when the structure is excited at about half its first resonance frequency.In addition to the extraordinary differences in magnitude, the higher harmonics turn some plots rather complicated.Loops before a plot closes on itself are plain; nonetheless, the response is still periodic

Conclusion
A composite panel reinforced with aluminium stiffeners, connected by an adhesive layer of Araldite® 420 A/B was analysed.Debonds were considered in the latter connection, which had been identified as sensitive in previous works.
A simplified modelcontact was not accounted forwas employed to estimate the modes of vibration.It was found that the latter are only slightly affected by flaws, both in what concerns the frequencies and the natural mode shapes.The differences between the modes of the damaged and sound structures are imperceptible, except for modes of vibration with large deflection in an area where the flaw is.The latter feature may be useful to estimate the damage location.
It was verified that the steady-state response to harmonic forces can be employed to detect damage.Debonding between the aluminium stiffener and the composite plate led to non-harmonic and nonsymmetric responses, which affected time histories, phase plane plots, frequency spectra and Poincaré sections.The greater the flaw, the greater the effect on the dynamics; very small flaws, in this particular study represented by a ratio between the debonded area and the total glued area equal to 0.4%, barely affect the response.Variations in velocities and accelerations are more significant and more likely to be noticed than alterations in displacements.The change in the response due to the flaw also varies with the point where measurements are taken.Higher-order harmonicsa sign that debonding occurredare more pronounced in the oscillations of some points (e.g., in the extremities of the composite panel) than in other, depending also on the excitation frequency.
Exciting the structure close to resonance increases the amplitude of the oscillations and, hence, the difference between the damaged and undamaged responses.But a more promising strategy is to employ an excitation frequency that is a sub-harmonic of the fundamental natural frequency (natural and resonance frequencies are close, due to low damping).The appearance of a higher-order harmonic that coincides with the fundamental natural frequency leads to very noticeable differences between the response of the damaged (non-linear) and the undamaged (linear) structure.Excitations at half the resonance frequency appear to be in particular promising, because damaged led to symmetry breaking.
In summary, the time histories, phase plane plots and frequency spectra are very useful in practical non-destructive VSHM on structures that behave linearly without defects.Those tools do not need a comparative response and allow to detect the presence of damage simply by exhibiting non-linearities.Poincaré sections also show the presence of damage, butin the cases here analysedrequire knowledge on the response of the undamaged structure.Therefore, methods based on Poincaré sections appear to be less practical.However, they should not be despised, because they use a smaller amount of data and permit defining quantitative indexes.A promising procedure, to employ in practical non-destructive VSHM on structures that behave linearly without defects, is to apply a single-frequency excitation at half a resonance frequency and obtain response spectra, preferably of accelerations.A second harmonic with significant amplitude is likely to be conspicuous in the response if a defect exists.This procedure should be repeated at several half-resonance frequencies for two reasons.One is to avoid that a defect goes unnoticed for being near a nodal line.A second reason is to obtain an indication of the position of the flaw: the response will be dominated by one mode in each resonance, if the flaw is in a position that undergoes large amplitude oscillations when this specific mode is excited, then its effects are more prone to appear in the response.Measurements should be taken at several points, again because mode shapes experience larger amplitudes in some points.

Data Availability Statement
The raw/processed data required to reproduce these findings cannot be shared at this time due to technical or time limitations.

B
. da Silva Henriques et al.

Fig. 3 .
Fig. 3. Tie interaction between the composite panel and the stiffeners.In pink the slave surface and in red the master one; the central damaged region, length l in , is not tied.The black dots mark points of interest for the analysis of section 4.4.

Fig. 4 .
Fig. 4. Boundary conditions in the modal analysis of the structure.

Fig. 6 .
Fig. 6.Harmonic loading applied on the central region in the structure.

Table 1
Material properties of aluminium alloy and of Araldite® 420 A/B.