Monitoring prestress in plates by sideband peak count-index (SPC-I) and nonlinear higher harmonics techniques

Propagating guided waves in a homogeneous, isotropic, prestressed, hyperelastic plate show nonlinear characteristics that depend on the state of initial prestress. These nonlinear phenomena include higher harmonic generation, occurring when Lamb wave modes of different frequencies (ωa\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _a$$\end{document} and ωb\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _b$$\end{document}) are allowed to mix within the material generating secondary waves at frequencies 2ωa,2ωbandωa±ωb\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2\omega _a,2\omega _b\, \text {and}\,\omega _a\pm \omega _b$$\end{document}. Further, if prescribed internal-resonance conditions are satisfied, the amplitude of secondary waves increases in space, providing a response quantity which is dependent on prestress and easy to be observed. Using the finite element method, in this paper we investigate the time and space evolution of higher harmonics arising in one-way wave mixing. The influence of prestress on the response is elucidated, observing the nonlinear parameter β\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document}. It is further shown that the nonlinear ultrasonic technique called sideband peak count-index (SPC-I) can provide an effective monitoring tool for prestress.


Introduction
Nondestructive evaluation (NDE) and structural health monitoring (SHM) based on elastic wave propagation in waveguides ranges from traditional ultrasonics, which rely on the linear theory, to nonlinear ultrasonics, which exploit some of the nonlinear phenomena observed in the experimental response and requires appropriate models. These models account for geometric and material nonlinearities and describe the occurrence of nonlinearities while waves propagate. In this paper, we numerically investigate the response characteristics of a prestressed plate in view of an understanding of the potential of nonlinear parameters to determine preexisting stress. This presents considerable challenge, but also has a lot of interesting engineering applications, such as the monitoring of stress in pressurized tanks or in other structural elements like truss members or rails.
Ultrasound-based techniques show high potential in different areas of NDE. Nonlinear ultrasonic (NLU) techniques enable the identification and tracking of material degradation at an early stage and provide an estimate of the damage state which is more reliable than that provided by linear ultrasonic (LU)-based NDE techniques [1]. NLU techniques exhibit high sensitivity to microstructural defects, fatigue [2], creep [3], material degradation [4] and stress [5,6]. A large amount of research was conducted about self and mutual interactions of nonlinear waves for the unique sensitivity of nonlinear wave interactions to material and geometric nonlinearities [7].
When two waves propagate in a nonlinear medium, their interaction occurs in a wave mixing zone, where mutual wave interactions result in combinational harmonics at the sum and difference frequencies [25][26][27][28][29][30][31][32][33][34]. If the two primary waves satisfy certain resonance conditions, the mixed wave is also a propagating wave, and its maximum amplitude is proportional to the size of the mixing zone and the distance travelled [19,25]. Note that one primary wave and its resonant second harmonic can be also such kind of mixed waves, generating a secondary harmonic in a condition of self-wave mixing. In the technique called one-way mixing, two primary waves propagating in the same direction [20,32] are employed. Among the researchers who studied this subject, Sun [19] investigated the wave mixing zone and the backward propagation phenomenon, Ishii [26] studied the interaction of guided elastic waves in an isotropic plate based on perturbation approach, and showed that the amplitude of the resonant harmonics increases linearly with the propagation distance, Hughes [27] and Yeung [29] evaluated the material nonlinearity from mixed waves, Ding [31] demonstrated the resonant wave mixing phenomena caused by micro-cracks, Ju [34] used the wave mixing technique to monitor thermal aging of adhesive joints.
The experimental observability of the cumulative effect of secondary and combination harmonics was demonstrated in several works, among which [27,29,31,35,36], respectively, for edge waves, pipes and shear waves in plates, damage and fatigue evaluation. With specific regard to Lamb waves in plates, Pineda Allen [37] showed that the combination internally resonant waves are spatially cumulative, moreover, Hu [38] observed the cumulative effect of the combination of low-frequency S0 modes. In the last paper, an experiment with an empty and filled tank demonstrates that the amplitude of the secondary combination peaks depends on the presence of the fluid, and hence on the stress too. The results reported in the literature point at the effectiveness of the combination harmonics in characterizing material nonlinearities, which could be extended to a problem of stress monitoring. One practical difficulty is related to the multiple sources of nonlinearity, among which also damage plays an important role. In the case of pressurized tanks, a manometer can measure the gas pressure, which will enable to separate the contribution of damage from that of stress.
Recently, a newly developed NLU technique called sideband peak count-index (or SPC-I) has shown promising results in engineering structural health monitoring. The details of SPC and SPC-I analysis process have been reported in the literature [39,40] and will be described briefly in Sect. 4. The SPC technique was first successfully applied in monitoring degradation of glass fiber reinforced cement composites [41]. Later, applicability of the SPC-I technique in monitoring damage evolution in polymer composite plates has been investigated [42]. Nondestructive evaluation of concrete using SPC-I technique has been also reported [43][44][45][46]. Experimental results show that the SPC-I technique is very sensitive to damage in concrete. Other relevant investigations using SPC-I technique involve crack detection in metallic materials such as aluminum plates and aircraft lugs [47]. Conclusions have been drawn from these investigations that for micro-cracks, SPC-I is more sensitive than LU parameters and it outperforms other nonlinear techniques. Recently, SPC-I has been applied to monitor porosities in additively manufactured parts [48] as well as steel tube welded joints [49] and adhesion defects in FRCM reinforcements [50]. Also, a recently modified SPC-I technique called sideband peak intensity (or SPI) has been proposed for monitoring impact damage and shows consistent trends as SPC-I techniques [51], which reinforces the conclusion that nonlinear SPC-I technique is reli-able and has advantages in monitoring early stages of damage.
In this work, we conduct numerical finite-element (FE) simulations to investigate the generation of higher harmonics in one-way mixing, and elucidate the dependence of their amplitudes on the initial prestress as well. The finite element model implements the second-order approximation of the equations of motion obtainable by perturbation approach. Different states of prestress are considered, including uniaxial and biaxial, with wavefronts orthogonal or parallel to the principal directions of prestress. The influence of prestress on the response is elucidated, observing the nonlinear parameter β. It is further shown that the application of the new nonlinear ultrasonic technique SPC-I enables an effective monitoring of the state of prestress.

Nonlinear equations of motion of a prestressed plate
The equations of motion of a prestressed solid can be formulated using the classical approach of acoustoelasticity. Three different configurations of the material points P are distinguished, which are: the natural configuration C 0 , free of stress and strains, the initial configuration C, which is a stressed and strained equilibrium state, and the current configuration C(t) (Fig. 1a). The coordinates of point P in the natural, initial and current configurations are, respectively: a(P), X(P) and x(P, t). We define these coordinates with respect to the same rectangular Cartesian common frame. A motion of the material body is a one-parameter mapping [7]: where u i is the static displacement vector field taking from the natural to the initial state. Physical quantities in the initial state are referred to with the superscript "i". It is assumed that the initial state is attained by a static displacement of the kind u i = αa, where α and A are given diagonal tensors related to the initially applied stress; a and A can be assumed to be diagonal without loss of generality by letting the principal axes coincide with the reference frame. The final con- figuration is reached superimposing to the initial state a dynamic disturbance u(t).
Strains are defined using the Green-Lagrange tensor of finite strains: where the superscript " f " represents the final state and ∇ is the gradient operator with respect to the natural coordinates a. We will consider a hyperelastic material admitting a strain energy function . If the medium is isotropic, its strain energy function is unaffected by any arbitrary rotation of the reference frame, hence can be written as a power series of the invariants I 1 , I 2 , and I 3 of the strain tensor E f . Neglecting infinitesimals of order higher than three in , the strain energy can be expressed as: where I 1 = trE f , I 2 = (trE f ) 2 − tr(E f 2 ), and I 3 = detE f , λ and μ are the Lamé constants (second order) and l, m and n are the Murnaghan constants (third order) [52]. The constitutive relation expressing the first Piola-Kirchhoff stress S in this hyperelastic body is obtained manipulating the strain energy. According to Murnaghan, the relation between the Cauchy stress tensor σ and the strain energy function is written as: from the relation between the Cauchy stress tensor and the first Piola-Kirchoff stress tensor, that is: ( This results into the following expression of the first Piola-Kirchhoff tensor S: where co indicates cofactor. Finally, we can set up the field equations of motion in the natural configuration and in the absence of body forces other than inertia as: where ρ 0 is the natural material density, the operator Div involves derivatives with respect to natural coordinates a. When considering a wave propagating in a plate (Fig. 1b), free-stress equations on the upper and lower surfaces have to be added to Eq. (7), that is: where n 3 is the unit vector normal to the upper and lower surfaces of the plate.

Perturbation approach to the solution of second-order approximation
A second-order approximation of the equations of motion entails retaining first and second powers of α in terms not involving ∇u; first and second powers of α in terms involving ∇u; discarding terms involving powers higher than two of ∇u. In this way, the stress tensor can be divided into two parts: S I collecting the first order terms, and S II which gathers the second order terms: 3 . Its second-order approximation will then be written as: A classical way to solve of the set of Eqs. (7) and (8) is perturbation approach. This implies writing the solution u as the sum of a primary solution u 1 and a secondary solution u 2 : This approach is illustrated in detail in [12] and [5].
Here we will recall only the main results, which will guide interpretation of the numerical results. Substituting Eq. (10) into Eqs. (7) and (8), and separating linear (superscript I ) and quadratic ( II ) terms of S NL containing derivatives of u 1 and u 2 (superscripts 1 and 2 ), yields [7]: The terms S I1 and S I2 contain only the derivatives of u 1 and u 2 , respectively. The term S II1 instead contains the quadratic terms in u 1 . Finally, S II2 contains quadratic terms in u 2 and mixed products of u 2 and u 1 , and is neglected in the analytical approach to the solution of the nonlinear problem [5,12]. The obtained hierarchy of equations of motion consists in a free-vibration linear problem which coincides with the linear approximation: and a forced linear problem where the forcing terms depend on the solution of the first-order problem: where f 1 = DivS II1 and S II1 are the volume and the surface forcing term, respectively. They are both known once the solution to the first-order problem is determined.
Equations (12) and (13) have some analogies with those of parametrically excited structures, in which case, the parametric excitations appear as forcing terms in the equation of motion of order 1 , similarly to the forcing term in Eq. (13) [54][55][56]. Consequently, when the frequency of the parametric excitation coincides with one of the frequencies of the system, the amplitude of the combination terms is enhanced (principal parametric resonances). However, it must be noticed that, in the case treated in this paper, the forcing term of Eq. (13) is tied to the free-response at the zero-order.
The first-order problem Eq. (12) describes the freevibration of infinite plates. For a given angular frequency ω, the infinite admissible wavenumbers k m can be obtained by solving the eigenvalue problem which is obtained enforcing wave-like solutions. Note that the symbols ω and f are used, respectively, for the frequency in rad/s and in Hz. The related eigenfunction U m (a 3 ) is the m-th wave mode shape. The results are usually represented in the dispersion diagram, which includes a set of curves, each one associated with a wave mode, which represent wavenumbers or phase or group velocity as a function of frequency. Dispersion curves depend on the initial state of prestress (see [5]). The m-th wave mode at frequency ω can be written as: These modes satisfy a reciprocity relationship which is the analogous of the orthogonality condition for wave modes [53]. The reciprocity enables writing the forced response in terms of wave mode superposition and also determining the expansion coefficients. Since the freevibration wave modes are the same for the second-order homogeneous problem, the second-order forced solution to Eq. (13) can also be written in terms of wave mode expansion.

Generation of higher harmonics
It is clear from the equations presented in the former subsection that if two waves with different frequency (ω a and ω b , with ω a > ω b ) propagate in a plate, wave mixing occurs. This means that higher harmonics of the fundamental frequency are generated (2ω a and 2ω b ), as well as the sum (ω a + ω b ) and difference (ω a − ω b ), resulting in the phenomenon called wave mixing. Figure 2 shows a schematic frequency spectrum which depicts the interaction of fundamental waves and the combinational secondary harmonic generation due to wave mixing. Let us assume that two waves with angular frequencies ω a and ω b and corresponding wavenumbers k a and k b are excited at a 1 = 0. The primary waves of amplitudes A and B will be expressed by: Substituting (15) into the second-order equations (13) results into volume f 1 and surface S I2 forcing terms which contain sum and difference frequencies: where f 2ω a and S 2ω a (f 2ω b and S 2ω b ) are the amplitudes of the forcing terms due to the self-interaction of the excited mode at frequency ω a (ω b ) and f + , S + , f − , S − are due to the mutual interaction of the modes at sum (ω a + ω b ) and difference frequencies (ω a − ω b ). The second-harmonic generation is considered as a special case of sum frequency generation, in which only one primary single mode is excited.
The secondary solution is written in terms of modal superposition including all the frequencies deriving from self and mutual interaction: is the spatially modulated amplitude of the m-th mode at frequency w, which can be obtained from the solution of the second-order forced problem exploiting the reciprocity relation. It was shown that the amplitude of the secondary solution remains bounded when the wavenumbers k w m at frequencies w differ from the wavenumbers k a , k b or k a ± k b [12]. On the contrary, when some internal resonance conditions are satisfied, that is: k w r = k a , k b or k a ± k b (synchronism or phase velocity matching) where subscript r stands for internally resonant, and there is nonzero volume and surface power flux between primary and secondary modes (see [12] for the definition of volume and power flux), the amplitude of the secondary r -th mode grows linearly in the direction of propagation. This cumulative effect occurs at the expense of the primary lower harmonics, in the absence of further energy fed to the system, as in our case, so that the energy balance is satisfied. In the absence of internal resonance, all modes are needed to represent the secondary solution. However, when one mode is in resonance with the primary wave, this mode is the dominant term in the solution. Henceforth, it is desirable that waves a and b are excited at the same location and that they satisfy self or internal-resonance conditions with waves at sum or difference frequencies, which ensures an effective mixing.
Similarly to that, in nonlinear vibrations, we have an internal resonance when one or more natural frequencies are commensurable. When an internal resonance coincides with a parametric resonance, the combination of the two types gives rise to simultaneous resonances, and the system vibrates at more than one mode at different frequencies, although only one frequency is directly excited by the parametric excitation.
The amplitudes and rate of accumulation of the harmonics can be obtained by measuring the response at some points following the mixing. This is expressed by the nonlinear parameter β, which is the constant ratio between the amplitude of the secondary resonant harmonics and the product of the amplitude of the primary waves and the abscissa where the displacements are observed [25,27]: where the superscripts 2a, 2b, a + b and a − b refer to the frequencies ω a , ω b , ω a + ω b , and ω a − ω b respectively.

Description of the model
A FE model was used to numerically simulate nonlinear guided wave propagation in a plate. A time-step analysis was carried out with COMSOL, using the equations of motion (7) and (8) [5,6], where the stress tensor is defined as the second-order approximation of Eq. (9). Attenuation is omitted. Different from what happens in the two-scale approach to the solution, in this way we will include also the contribution of quadratic terms in u 2 and mixed products of u 2 and u 1 . The equations of motion are enforced in a 2D domain (Fig. 3a), which is an area with thickness h = 10 mm and length l = 4000 mm, representing a plate with upper and lower surfaces free of stress. The length is sufficient not to see the reflection from the right boundary. The 2D setup of the model implies that geometrical spreading is ignored. The material is 7075-T651 Aluminum, whose mechanical properties are reported in Table 1.
A state of plane strain is assumed for the plate, whose motion takes place in the plane a 1 -a 3 (Fig. 3a), with propagation in the a 1 direction. This corresponds to a displacement field with two components u 1 and u 3 , with u 2 = 0. On the right end of the plate, displacements are constrained to zero. Primary waves are gen- erated enforcing a Dirichlet boundary condition on the left end of the plate. This is the so-called side incidence method, which entails forcing the boundary with a time-history of displacements which, in space, correspond to the u 1 and u 3 components of the primary wave mode shape we are willing to excite [5]. A sketch of the displacement field applied at the boundary is reported in Fig. 3a for the S0 mode. In time, we will use harmonic functions at the frequency of the primary wave modes enveloped by a Gaussian curve so as to include a number of cycles equal to 28, which provides sufficiently narrow peaks in the frequency domain. For the sake of brevity, an example concerning the wave mixing is reported in Fig. 3b. Frequency and mode shapes are the Rayleigh-Lamb waves obtainable from Eq. (12). The FE model was set so that the maximum element size is 1/10 of the shortest wavelength of interest and the time step is 1/100 of the largest frequency of interest. The elements used were four-nodes second-order Lagrange elements.
We will conduct numerical experiments with two different primary waves, whose frequencies and phase velocities are reported in Table 2. One case mixes two wave modes S0 at relatively low frequencies (Fig. 4). In Table 2 Frequencies (kHz) and phase velocities c ph (m/s) of the primary and secondary wave modes for wave mixing and self-resonance  Table 2 and Fig. 4. Despite that, we will show that quite strong secondary harmonics are obtained with cumulative growth of amplitude in space. The other case employs the S1 mode as primary wave, generating the secondary self-resonant S2 wave (Fig. 5). In this case we have exact resonance conditions, which not only occur for modes S1 and S2, but also for higher harmonics, in a condition of multiple internal resonance.

Higher harmonic generation in unstressed plates
For the sake of brevity, in this subsection, the results are monitored only for the numerical experiment involving wave mixing of the S0 mode in the unstressed case. Figure 6a and b report the u 1 displacement component respectively in time and space domains. Primary and secondary waves are present in each plot. The secondary waves can be observed in the tail of the response, in fact, they have slightly lower group velocities than the primary, as it is apparent from Fig. 4. Figure 6c reports the contour plot of the u 1 component, which has a symmetric distribution on the height of the plate compatible with the S0 mode. Also, this demonstrates that in the current simulation, all the waves generated by the mixing are S0 modes. Figure 7 shows the Fourier Transform of the response in terms of u 1 component of displacement on the plate surface, and at a distance of 100 mm from the left boundary. This picture shows the occurrence of the primary harmonics at the frequencies f a and f b contained in the disturbance applied at the boundary, together with secondary harmonics f 2a , f 2b , f a−b , and f a+b . Figure 7 also compares the response to a disturbance containing mixed frequencies to that with one singlefrequency: when excited separately, the 53.7 and 71.6 kHz S0 waves generate secondary harmonics at integer multiples of the fundamental frequency.
Frequency(kHz) Group Velocity(km/s) that are the ratios between amplitude of the secondary harmonics and products of the amplitudes of primary waves at different increasing abscissae (subscript a 1 ), divided by the same quantity these ratios assume at a 1 =100 mm. Note that when w = 2ω a , B = A and when w = 2ω b , A = B. In the four curves, the dots represent the values retrieved from the Fourier transforms of the numerical responses. These normalized amplitudes increase linearly, apart from some predictable slight deviation for the secondary harmonic 2ω a , due to its increased mismatch of the resonance conditions (Fig. 4). It can also be seen that the slopes of the curves do not differ much. These results indicate that cumulative harmonic generation mixing the primary S0 mode is achievable at low frequency. Given that in the problem under study attenuation and geometrical spreading are omitted, the amount of nonlinear cumulative effect which could be observed in experiments can be smaller.

States of prestress under investigation
Three different states of initial prestress are investigated and represented in Fig. 9, specifically: uniaxial stress in the direction of the wave propagation (a 1 , case A), uniaxial stress orthogonal to the direction of wave propagation (a 2 , case B), and plane state of stress in the plane (a 1 , a 2 ) with equal principal stresses along the two directions a 1 , a 2 which is called plane-isotropic state (case C). The corresponding strains are reported   Table 3 States of prestrain in Table 3 (ν = 0.33), where α is the applied initial strain, α 1 , α 2 and α 3 are the strains in a 1 , a 2 and a 3 directions, respectively. For each of the A, B, and C cases, both α > 0 (tensile strain) and α < 0 (compression strain) are investigated. It is taken a limit value α = 0.004, which is about 50% of the yielding strain for 7075-T651 Aluminum. It should be noticed that, due to the coupling induced by the nonlinearity, the prestrain applied in direction a 2 affects the motion in the plane (a 1 ,a 3 ) where displacements take place.

Higher harmonic generation in prestressed plates
In this section, the effect of initial stress on the generation of higher harmonics is illustrated. Figure 10 compares the response spectrum in prestressed cases A, B and C, for S0 wave mixing at frequencies f a = 71.60 kHz and f b = 53.70 kHz. The response quantity observed is the displacement u 1 on the plate surface 400 mm away from the left end, with tensile prestrain set to α = 0.001. We can see from the details reported in Fig. 10b-e that the largest increase in the amplitude of the secondary harmonics is obtained for case C, then A. The case B presents the smallest variation. Figures 11, 12 and 13 illustrate the effect of the increase of strain (α from 0.001 to 0.004) in the dif- Fig. 9 Schematic diagram of three states of stress, wave propagation in a 1 direction. Cases A and B: uniaxial prestress. Case C: prestress in the plane of the plate ferent states of stress A, B and C, respectively. The response quantity observed is again the component u 1 of displacement picked up 400 mm away from the left boundary. Similar to Fig. 10, in cases A and C an increase of strain causes a increase of the FFT amplitudes of all the secondary harmonics, while in case B an increase of strain results in a reduction of the secondary-harmonic amplitudes. Quantitatively, the variation of case B is smaller than in other cases.
The parameter β is useful in quantitative analysis of the material and geometric nonlinearity. To understand the effect of prestress, we introduce and observe the parameter β, which is defined as the difference between β α in prestrained conditions and in stress-free initial state (α = 0), i.e., Figure 14 reports the quantity β a−b a 1 for the different prestress states as a function of α, for a 1 = 400 mm.
The picture reports both the discrete points for those values of α at which the response was observed, and their linear fit. It can be seen that a linear increase in β with increasing α is observed in cases A and C, which corresponds to a cumulative second harmonic generation for the all harmonics. Also, in case B, the magnitude of β linearly increases for increasing α, but it is an out-of-phase contribution with respect to the wave propagating in the unstressed medium. Figure 14 shows that the nonlinear harmonics have the highest sensitivity to stress when they propagate in the same direction as the nonzero component of the uniaxial stress, and the lowest when they propagate orthogonally to it.
It can be concluded that the observed response quantity β has the potential to be employed in stress monitoring. Moreover, the presented results show that the interpretation of the information derived from the   Fig. 13 Details of the secondary harmonic peaks as a function of the initial state of strain-case C measured response is specific to the direction of the wavefront in relation to the direction of the principal directions of the preexisting stress (Fig. 9). Although this is in general a limitation, there are cases in engineering when the principal directions of the stress are known with reasonable approximation (i.e. truss members, rails, pressure vessels).

Nonlinear SPC-I analysis method
The SPC-I analysis method [40] is based on the analysis of the peaks of the secondary harmonics due to internal resonance and wave interaction. In the previous sections it was shown that geometric and material nonlinearities give rise to secondary harmonics whose amplitude depends on the initial prestrain. A sample SPC plot is shown in Fig. 15b. It is generated by moving a threshold line, shown by the horizontal continuous line in Fig. 15a. The threshold line is moved vertically between the lower and upper threshold limits, shown by two dashed lines. All peaks that are above the moving threshold line are counted and plotted against the moving threshold value. It should be noted that both the number of peaks and their strengths affect the SPC plot which is a measure of the degree of material nonlinearity. Thus, the SPC plot (number of peaks as a function of the threshold value) gives a visual representation of the material nonlinearity [58].
The SPC-I is an index value which is the average of SPC values for all threshold positions. This index indicates the degree of material nonlinearity. The higher the material nonlinearity, the greater is this number. In this work, after verifying the occurrence and sensitivity of secondary harmonics to the prestress, we adopt the SPC-I technique for monitoring prestress-induced nonlinear response in the plate structure. In this section, SPC-I technique is applied to the analysis for both S0 wave mixing and S1 self resonance conditions. The three states of stress shown in Fig. 9 with different tensile strain values are considered. We focus the analysis on tensile strains to avoid practical situations which, in thin plates, could be tied to instability.

SPC-I for wave-mixing conditions
To get clearer results, we select the spectral plots which correspond to the case when the distance is 2000 mm from the leading edge of the excitation boundary. As mentioned above, in the SPC-I analysis the number and intensity of peaks above certain threshold is recorded. For case A (the geometry is shown in Fig. 9a), different peaks detected by our software are marked in the logarithmic spectral plots of Fig. 16 for different states of tensile prestrain. When the elastic wave propagates through materials with higher nonlinearity, then, in the spectral plot the number and the strength of the peaks increases. Case B and C can be analyzed in the same manner, but the spectral plots are not shown for the sake of brevity.
Analysis of Fig. 16 shows that additional peaks are generated around the main envelope, and the sideband peak amplitudes are much less than that of the main amplitude. A threshold value of 20% of the maximum peak for each logarithmic spectral plot was used as the lower threshold limit, and 40% of the maximum peak for each logarithmic spectral plot was used as the upper threshold limit. Only the peaks above this moving threshold line were counted. The SPC plots (the number of peaks) above the threshold value as it varies from 20 to 40% of the maximum peak value are shown in Fig. 17 for cases A, B and C for different strain levels.
The SPC-I value is the average of all SPC values (peak counts above the threshold position) for all  Fig. 18. These results clearly show that, in all cases, SPC-I increases with the stress level which is consistent with the variation of the amplitudes of harmonics. It should be noted that although SPC-I is very sensitive to the pre-stress level, for S0 wave mixing it is found to be not so sensitive to the direction of the principal stresses relative to the wave propagation direction. For this reason the three figures in Fig. 18 look very similar.

SPC-I for second harmonic generation conditions
In a previous paper [6], we studied the second harmonic generation of S1-S2 Lamb mode pair and introduced the relative nonlinear parameter to quantify the degree of nonlinearity in the response, as a function of the initial state of prestress. The results showed a linear increase with increasing of tensile stress in cases A and C, and a linear increase in magnitude but opposite phase in case B. Among them, case B has lowest slope while case C is tied to the largest variation (Fig. 19).
In this paper, we carry out the SPC-I analysis as a supplement to the previous study. For cases A, B and C, the lower threshold limits are set as 20% of the maximum peak, and the upper threshold limits are set as 40% of the maximum peak for each logarithmic spectral plot as shown in Fig. 19. The SPC plots are shown in Fig. 20 for cases A, B and C for different strain levels when moving threshold lines are varying between lower and upper threshold limits mentioned above. The response is taken 200 mm away from the left end. The SPC-I results are shown in Fig. 21. They report the results of analyses carried out for a frequency thickness product equal to 3.55 Mhz*mm with tension prestrain values varying from 0.001 to 0.004.
It can be seen that in all cases, SPC-I increases with tensile stress, which corresponds to a higher degree of nonlinearity with larger strain value. At the same time, from Fig. 21, case B presents the lowest variation and case C the highest sensitivity to prestrain, which coincide well with our previous results in [6]. It is also observed that, when self-resonant waves are employed, a remarkable sensitivity to the state of prestress is achieved, which is different from what was found in wave mixing.

Summary and discussion
In this paper, we investigated the frequency mixing of primary wave modes in prestressed plates by using a finite element model. We employed a model accounting for both geometric and material nonlinearities and considered either wave mixing of two S0 wave modes and S1-S2 self-resonant modes. It was shown that the generation of mixed secondary frequencies from two primary S0 modes is achievable. Besides frequencies at natural multiples of the primary frequencies, the wave mixing generates additional harmonics at the sum ( f a + f b ) and difference ( f a − f b ) frequencies, increasing the number of independent harmonics which can be utilized for material nonlinearity measurements. Despite the fact that when using a couple of S0 modes resonance conditions are not strictly satisfied, the amplitude of all of these secondary waves increases while they propagate in space. This increase is quantified by the slope of the linear increase, that is, the nonlinear parameter β, which was employed to quantify the material nonlinearities. Moreover, the use of S0 modes at low frequencies is favorable in experimental practice thanks to their relatively low wave speed.
Different initial states of prestress were investigated, including waves propagating in the same direction as  In conclusion, the sideband peak count-index technique was used to monitor prestress. This technique takes advantage of the number and strength of peaks of the secondary harmonics. In numerical applications SPC-I showed remarkable sensitivity to stress and proved to be a very promising tool for monitoring the amount of stress in plates, both when using S0 mixing waves and when using S1-S2 self-resonant waves. However, self-resonant waves S1-S2 proved to be much more sensitive to the principal directions of the stress than S0 self-resonant waves, which could indicate a preference toward the first kind of excitation when the orientation of the principal directions of the stress are known.
The use of S1-S2 modes and S0 mode has their advantages and disadvantages. Not only S1-S2 proved to be more sensitive to the state prestress, they also expressed higher sensitivity in experimental applications for micro-scale surface damage detection [35] and fatigue damage evaluation [36]. However, their use in practice is challenging for the need of using high frequencies and for some amount of dispersion that can arise. The same applies to the other limited number of mode pairs that meet the internal resonance conditions, such as A2-S4 and S2-S4 [19], and which were considered in the literature as mode pairs for inspection. Some limitations of these mode pairs can be resolved by using the S0 modes at low frequencies, which approximately satisfy phase-velocity matching. However, this low-frequency range is narrow, and the sensitivity to stress is lower. nor the European Commission can be considered responsible for them.
Funding Open access funding provided by Università degli Studi di Roma La Sapienza within the CRUI-CARE Agreement. The authors have not disclosed any funding.

Data availability
The datasets generated and analyzed in the current study are available from the corresponding author on reasonable request.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.