Higher-Harmonic Generation Analysis in Complex Waveguides via a Nonlinear Semianalytical Finite Element Algorithm

Propagation of nonlinear guided waves is a very attracting phenomenon for structural health monitoring applications that has received a lot of attention in the last decades. They exhibit very large sensitivity to structural conditions when compared to traditional approaches based on linear wave features. On the other hand, the applicability of this technology is still limited because of the lack of a solid understanding of the complex phenomena involved when dealing with real structures. In fact the mathematical framework governing the nonlinear guided wave propagation becomes extremely challenging in the case of waveguides that are complex in either materials damping, anisotropy, heterogeneous, etc. or geometry multilayers, geometric periodicity, etc. . The present work focuses on the analysis of nonlinear second-harmonic generation in complex waveguides by extending the classical Semianalytical Finite Element formulation to the nonlinear regime, and implementing it into a powerful commercial Finite Element package. Results are presented for the following cases: a railroad track and a viscoelastic plate. For these casestudies optimum combinations of primary wave modes and resonant double-harmonic nonlinear wave modes are identified. Knowledge of such combinations is critical to the implementation of structural monitoring systems for these structures based on higher-harmonic wave generation.


Introduction
Traditional techniques in nondestructive evaluation and structural health monitoring applications rely on measuring "linear" parameters of the waves amplitude, speed, and phase shifts to infer salient features of the inspected structure.Several studies, however, have shown that "nonlinear" parameters are, in general, more sensitive to structural condition than linear parameters 1 .Furthermore, the use of nonlinear guided waves is extremely attractive because guided waves combine the mentioned high sensitivity typical of nonlinear parameters with large inspection ranges 2-9 .
From a mathematical standpoint, the framework behind nonlinear guided waves propagation is relatively challenging since the Navier elastodynamic equations are further complicated by stress-free conditions at the waveguide's cross-sectional boundaries.For this reason, most of the previous works on elastic waves in waveguide solids considered the propagation to be in the linear elastic regime with the assumption of infinitesimal deformations coincidence between deformed and initial configurations .However, as the amplitude of the wave increases or the structure starts experiencing finite deformations i.e., nonlinear elasticity or another cause of nonlinear effects is present, the nonlinearity in the structural response becomes relevant and must be introduced in the analysis.Hence cubic and eventually higher-order terms in the particle displacements gradients must be included in the elastic strain energy density expression 10, 11 .
Among the manifestations of the nonlinear behavior, higher-harmonic generation is considered in detail in the present work.In this scenario, supposing to excite an ultrasonic wave into the waveguide structure at a fixed frequency, ω Fundamental Frequency , the nonlinearity manifests itself in the generation of multiple harmonics of ω, for example, 2ω second harmonic , 3ω third harmonic , and so on.For a practical use, this nonlinearity can be quantified via an ultrasonic nonlinear parameter, β, well documented in literature 2 .
In the last thirty years, several successful applications of nonlinear guided waves have been discussed, spanning from assessing the fatigue damage of metals 12-14 and concrete 15 , to the efficient location of internal cracks and dislocations 16-20 .The authors of the present paper recently exploited the features of nonlinear guided wave propagation in sevenwire steel strands and proposed a methodology to measure the stress level acting on these structural elements based on the theory of contact acoustic nonlinearity 21 .While several investigations pertaining to nonlinear effect in solids and second harmonic generation were reported in the past 22, 23 , most of them were limited in their applicability to structures with simple geometries plates, rods, and shells where analytical solutions for the primary linear wave field are available in literature.Very few studies tried to analyze the nonlinear wave propagation phenomena in geometrically complex waveguides using specialized SAFE codes 24 .
In the present work, the propagation of waves in nonlinear solid waveguides with complex geometrical and material properties is investigated theoretically and numerically.For the solution of the nonlinear boundary value problem, perturbation theory and modal expansion are used 22 .The main novelty consists in the development of a powerful numerical algorithm, able to efficiently predict and explore the nonlinear wave propagation phenomena in several types of structural waveguides.It is based on the implementation of a nonlinear semianalytical finite element formulation into a commercial multipurpose finite element package.Compared to the classical finite element formulation, the proposed solution is computationally more efficient since it simply requires the finite element discretization of the cross-section of the waveguide and assumes harmonic motion along the wave propagation direction.Furthermore, compared to traditional spectral or waveguide element method approaches, no new elements need to be developed, the full power of ready-touse high-order shape functions crucial for the development of the present theory can be easily exploited though friendly GUI, and immediate and extensive postprocessing for all the required quantities can be developed.
The applicability of the proposed analysis is quite wide, since it can efficiently handle general prismatic structures, viscoelastic waveguides with damping effects, multilayered composite laminate panels, and heterogeneous systems, all cases where theoretical wave solutions are either nonexistent or extremely difficult to determine.In addition, the proposed approach requires simple modifications to the original commercial FEM code so that the nonlinear semianalytical formulation can be taken into account and translated to match the required formalism.After a brief discussion on the background of the present work and the proposed algorithm, two case studies have been analyzed in detail: a railroad track and a viscoelastic plate.They were considered to show the potential of the algorithm in handling complex geometry as well as viscoelastic material properties.The proposed code was successful in identifying optimal combinations of resonant primary and secondary modes.The knowledge of these nonlinear resonance conditions is of paramount importance for the actual implementation of conditions assessment systems for these structures that are based on the measurement of nonlinear ultrasonic guided waves.

Nonlinear Guided Waves Propagation
In the present section, a brief overview of the generalized nonlinear theory of elasticity for wave propagation involving finite deformations is presented 25 .Following 22 , assuming that the body is homogeneous, isotropic, and hyperelastic, it possesses a strain energy density ε which is an analytic function of the Green-Lagrange strain tensor E ij via its invariants; in this scenario, the Second Piola-Kirchoff stress tensor S ij can be expressed as: where ρ 0 is the initial density of the body.
According to finite strain theory, in 2.1 we have assumed the following: The strain energy density expression becomes where I 1 , I 2 , and I 3 are the first three invariants of the Green-Lagrange strain tensor defined as 2.5 Characterizing the system of 2.5 to the "guided" wave propagation case stress-free boundary condition , the governing equations can be recast in vector notation as: where u is the particle displacement vector, ρ 0 , λ and μ are defined above, f is the nonlinear term acting as a body force, n r is the unit vector normal to the surface of the waveguide Γ, and S L and S are the linear and nonlinear parts of the second Piola-Kirchoff stress tensor, respectively.The nonlinear waveguide system is illustrated in Figure 1.
Considering higher harmonics up to the second order, the nonlinear boundary value problem presented in 2.6 is solved using perturbation theory.The solution of the primary wave field can be obtained analytically for simple geometries plates, rods, shells, etc. and numerically using the classical SAFE formulation for waveguides with generic crosssection 27 .Following 22, 28 , if ω is the primary frequency that is excited into the system, the first-order nonlinear solution is calculated through modal expansion using the existing propagating guided modes 2ω as: where x, y are the cross-sectional coordinates of the waveguide, z is the lengthwise coordinate of the waveguide, c.c. denotes complex conjugates, v m is the particle velocity vector referred to the mth mode at 2ω, and A m is the higher-order modal amplitude given by: where k represents the wavenumber.The amplitude A m z quantifies how strong is the excitation of the mth secondary mode in the modal expansion.In 2.8 , the amplitude of the secondary modes is expressed in two different forms depending on the existence or not of the phase-matching condition synchronism .The latter occurs between two modes having the same phase velocity.The expressions are where P mn is the complex power flow along the direction of wave propagation and f vol n and f surf n are identified as the complex external power due to surface sources and volume force, respectively.
It is possible to notice how the nonlinearity of the waveguide transforms a monochromatic single frequency wave input into a distorted output where primary wave and second harmonic coexist Figure 1 .Furthermore the modal amplitude of the generic mth secondary mode oscillates in value if the solution is asynchronous, while it increases with propagation distance z if the solution is synchronous.The latter is the known cumulative behavior occurring for nonlinear resonant modes.Further details concerning the terms appearing in 2.9 -2.10 can be found in 22 .The internal resonance mechanism relies on the simultaneous occurrence of two conditions, namely: 2 Nonzero power transfer from primary to secondary wave field: Recent investigations performed by Deng et al. have analyzed the influence of an additional requirement for the occurrence of internal resonance, namely, the group velocity matching 29 .In this study, the authors showed analytically and experimentally that, as long as the two aforementioned conditions phase-matching and nonzero power transfer are satisfied, the cumulative effect of the secondary resonant mode takes place even when the group velocity matching condition is not satisfied.They concluded that group velocity matching does not represent a necessary requirement for cumulative second-harmonic generation.For this reason in the present work, phase-matching and power transfer only are considered in detail.
In nonlinear structural monitoring, the key consists of the identification of an optimal combination of synchronous primary and secondary modes.The rest of this paper presents a numerical tool that enables to identify these resonant conditions for various complex waveguides, that would be extremely difficult to study by other means, and that include cases of periodic structures, damped structures, multilayered geometries and heterogeneous structures.

Nonlinear Semianalytical Finite Element Algorithm
Linear SAFE formulation has shown in the past its great potential in calculating the dispersion characteristics of complex waveguides where the analytical solution is not available in a very efficient way 27, 30 .The knowledge of these curves is the starting point for the development of any application based on the use of guided waves.The present work focuses on the extension of this approach to the nonlinear regime and its implementation, into a highly flexible COMSOL commercial code, of a nonlinear SAFE formulation to solve complex waveguides CO.NO.SAFE Algorithm .
The implementation combines the full power of existing libraries and routines of the commercial code with its ease of use and extremely capable postprocessing functions; hence internal resonance conditions of structural waveguides with different level of complexity can be conveniently analyzed via user-friendly interfaces.Furthermore, since all the nonlinear parameters involve gradients of the displacement field up to the third order, high-order finite elements at least cubic need to be used in order to obtain meaningful results; this task is not trivial to implement in general SAFE algorithms.
Starting from the nonlinear boundary value problem stated in 2.6 , the displacement field is approximated in the cross-section of the waveguide x, y and is enforced to be harmonic in time and along the direction of wave propagation z in accordance with the classical SAFE formulation.For the generic eth element, this condition reads u e x, y, z, t N e x, y U e e i kz−ωt , 3.1 where N e x, y is the matrix of shape functions, t is time, k is the wavenumber, and U e is the nodal displacement vector for the eth element.The enforcement of this particular displacement field in 2.6 constitutes the main modification that needs to be applied in the original cross-sectional FEM formulation.Hence, after the original quadratic eigenvalue problem in wavenumbers has been reformulated in a linear fashion by doubling the space dimension 27 , the nonlinear boundary value problem can be implemented in COMSOL using the general PDE solver engine 31 .COMSOL formalism for the boundary value where n is the outward unit normal vector on the surface of the waveguide, c is the diffusion coefficient, α is the conservative flux convection coefficient, d a is a damping coefficient, β is the convection coefficient, a is the absorption coefficient, γ is the conservative flux source term, f is the source term, q is the boundary absorption term, λ is the eigenvalue and U represents the set of dependent variables to be determined.All these coefficients generally admit complex values appropriate for viscoelastic materials 32 .The formalism introduced in 3.2 -3.3 is very general and can be used for a broad range of physical problems governed by a system of PDEs, once every coefficient has been conveniently characterized to the particular physics governing the considered problem.Once all the parameters have been defined, dispersion curves for the selected waveguide can be promptly calculated.Next, after a particular frequency has been selected as primary excitation, second harmonic generation and internal resonance occurrence can be analyzed.
In the next section, the proposed algorithm is benchmarked with two case studies of interest in structural engineering.

Railroad Track
A136RE railroad track was considered first for this study.Due to the complex geometry of the cross section, solutions for the dispersion curves and, consequently, for the higher harmonic generation analysis cannot be calculated analytically.After a preliminary study involving the selection and the analysis of internal resonance conditions for several primarysecondary wave field combinations, two exemplary cases were selected as representative.In the first case, phase matching between primary and secondary modes is verified.However, due to the characteristic energy distribution over the rail cross-section, no power transfer is present between the modes and, consequently, internal resonance does not occur; hence, the secondary modal amplitude is bound in value and oscillates with distance along the direction of wave propagation 3.1 .In the second case, instead, both required conditions are verified and internal resonance takes place, leading to a resonant secondary wave field growing linearly with wave propagation distance.
The material properties considered are given in Table 1.The Landau-Lifshitz thirdorder elastic constants are detailed in 33 .The geometry of the railroad track cross-section, with the FE mesh used for the analysis, is shown in Figure 2 a .To correctly explore the displacement field and all the derived quantities essential for the calculation of all the terms during the nonlinear postprocessing , 618 cubic Lagrangian triangular isoparametric finite elements were employed 34 .In Figure 2 b , the resultant phase-velocity dispersion curve in the 0 ÷ 200 kHz frequency range is represented.As detailed in the following, the same figure also pinpoints the two selected combinations of primary and secondary modes as exemplary cases for the internal resonance analysis.
The complexity of the guided wave propagation for this particular waveguide is evident considering the abundance of propagative modes present and their dispersion characteristics especially at higher frequencies .Selecting a primary excitation frequency of 80 kHz, the eigenvalue problem has been solved, and 500 propagative modes real eigenvalues have been extracted at ω 80 kHz and at 2ω 160 kHz .Next, Figure 3 shows some propagative modes found in this specific frequency range.It can be noted how differently the energy is concentrated within the waveguide.

Nonresonant Combination
A flexural horizontal primary mode was selected as primary excitation input for the CO.NO.SAFE algorithm .The nonlinear analysis revealed the presence of a synchronous secondary mode at 2ω similar flexural horizontal displacement distribution .However, the second required condition concerning the power transfer is not met for this particular combination, leading to an oscillating secondary modal amplitude value and absence of internal resonance.At the same time, a conspicuous power transfer is present between the selected primary mode and some asynchronous secondary modes; here again, because of the lack of one of the two essential requirements phase matching internal resonance does not take place.This fact translates into the very small value associated of modal amplitude associated with the only synchronous mode and the relatively higher values associated to the asynchronous secondary modes.
The following Figures 4 a and 4 b illustrate the selected primary and secondary modes, respectively.Figure 4 c plots the modal amplitude results as calculated from 3.2 for the propagative secondary modes present at 160 kHz.

Resonant Combination
In this case a flexural vertical mode was selected as primary excitation.The results of the nonlinear SAFE analysis disclosed the presence of some synchronous secondary modes with one in particular slightly different flexural vertical type able to verify both requirements producing internal resonance and a nonlinear double harmonic growing linearly with distance.As in the previous case, Figures 5 a -5 b display the selected modes, while Figure 5 c spotlights the very high value of modal amplitude related to the secondary resonant mode; small amplitude values associated to the other synchronous modes, for which power transfer is absent, are also shown in the same figure.
The previous results point up an optimal combination of primary and secondary wave fields able to maximize the nonlinear response of the waveguide.Furthermore, it is worthwhile to notice how the selected primary mode is not only able to produce a resonant condition, but also very attractive in terms of practical excitability by a piezoelectric transducer.

Viscoelastic Isotropic Plate
A viscoelastic isotropic high-performance polyethylene HPPE plate was investigated next to extend the applicability to dissipative waveguides.This system is of primary importance in aerospace and mechanical engineering and has been studied quite extensively in the past assuming linear elastic regime to obtain dispersion curves and associated waveguide modes 27, 35, 36 .In the present work, these results are confirmed and extended to the nonlinear regime; an efficient combination of resonant primary and secondary modes is identified and discussed in detail.Material and geometrical properties for the plate are illustrated in Table 2 35, 36 , where ρ is the density, h is the thickness, c L is the longitudinal bulk wave velocity, c T is the shear bulk wave velocity, k L is the longitudinal bulk wave attenuation, and k T is the shear bulk wave attenuation.The dissipative behavior of the plate was modeled via the Hysteretic formulation 27 .Hence, the resultant stiffness matrix is frequency-independent and was calculated just once at the beginning of the analysis once the complex Lame's constants were evaluated.The results for the present case are

4.1
In 4.1 the complex bulk wave velocities longitudinal and transverse are calculated as follows:

4.2
The resultant viscoelastic stiffness matrix, with terms expressed in GPa, is given by: First, the plate system was solved in the linear regime in order to calculate the dispersion curves and obtain the propagative modes, necessary for the nonlinear analysis.For this purpose, an extension of the linear SAFE algorithm 32 was employed.It allows the study of the guided wave propagation along structures exhibiting material/geometrical periodicity along their width which is normal to the direction of propagation and to the thickness and considered infinite by applying the so-called periodic boundary conditions PBCs .With this powerful tool, a generally complex periodic structure grooved panel, reinforced concrete elements, just to mention a couple can be modeled simply by considering a very small cell and applying PBCs on its sides.Mathematically, they represent a particular case of Neumann boundary conditions: the variables and their derivatives up to the element order are forced to take identical values on the pair of boundaries of the structure where the PBCs are applied.This tool is very attractive since it opens new possibilities to study the guided wave propagation linear and nonlinear for a general class of periodic structures by developing the analysis just on a small portion periodic cell .
According to this approach, the present plate system was modeled using a mesh of just 60 quadrilateral cubic Lagrangian elements mapped and deployed in a 3.17 × 12.7 mm periodic cell Figure 6 a .The resulting Lamb wave solutions are displayed in Figures 6 b -6 c in the 0 ÷ 500 kHz frequency range.They are found to be in perfect agreement with well-known results previously published in literature.Primary and secondary modes for the nonlinear analysis are highlighted with white circles in the same figures.
Due to the lack of studies in literature concerning specifically the HPPE material, the third-order Landau-Lifshitz elastic constants of a very similar plastic polymer Polystyrene were adopted for the nonlinear analysis 37 .The assumed values are A −10.8 GPa, B −7.85 GPa, and C −9.81 GPa.
The nonlinear analysis was developed between 250 kHz primary mode and 500 kHz secondary mode .The waveguide being dissipative, all the eigenvalues and eigenvectors are complex.Propagative modes were separated from evanescent and nonpropagative solutions by using a threshold of 10% between imaginary and real parts of each eigenvalue.After a preliminary analysis on different potential combinations among the propagative modes, one particular mode was selected as input primary mode for the nonlinear postprocessing.It is  The application of the CO.NO.SAFE algorithm in this case is simplified because of the assumption of 2D strain regime the plate is considered infinite in the width direction .For this reason all the terms used in the nonlinear postprocessing discussed before are evaluated on a line segment running through the thickness.This approach is sometimes referred as 1D SAFE 32 , and was first introduced almost four decades ago 38, 39 .
The results of the analysis pinpointed the presence of a resonant secondary mode.As mentioned before, while the contribution of all other modes is oscillatory and bounded 2.9 , this secondary mode shows a cumulative behavior and represents the dominant term in the expansion equation 2.7 with a contribution that linearly increases with distance.In fact, after all the secondary modal amplitudes were calculated from 2.10 for the synchronous case, the identified resonant secondary mode exhibits a value which is orders of magnitude larger than those associated to the asynchronous modes Figure 7 .The same figure also illustrates primary and secondary modes as contour plots height and color gradients are proportional to the out-of-plane displacement component along the propagation direction and 3D rendered views global modeshape considering a length of 1 cm.The amplitudes of the displacement fields are not normalized and, consequently, they supply exact information about the mode shapes.At the same time, the values are therefore not comparable from one mode to another.Figure 7 shows that the selected primary mode is a complex axial symmetric mode.The mode at the double harmonic shows also features typical of axial modes.This resonant secondary mode at 500 kHz looks very promising in a possible structural monitoring system because it keeps the majority of the energy in the central area of the cross-section and minimizes wave leakage into the surrounding medium.Furthermore, Figure 6 c shows that both primary and secondary modes have very small attenuation values especially the secondary mode at 500 kHz ; this fact makes the studied combination even more attractive because of the large inspection range that can potentially be achieved.

Conclusions
Nondestructive evaluation and structural health monitoring communities are showing an increasing interest in nonlinear guided waves because of their significant potential in several applications.However, proper application of nonlinear features requires a complete understanding of the higher-harmonic generation phenomenon that can be expected for the test waveguide.This paper discussed the extension of the classical SAFE algorithm to the nonlinear regime and its implementation in a powerful multipurpose commercial FEM code COMSOL .The result is an innovative tool that opens new possibilities for the analysis of dispersion characteristics and, most importantly here, nonlinear internal resonance conditions, for a variety of complex structural waveguides that do not lend themselves to alternative analyses such as purely analytical solutions.The specific cases that were examined in this paper include complex geometry railroad track and viscoelastic waveguides with damping effects HPPE plate .In all these cases, the proposed algorithm successfully identified optimal combinations of resonant primary and secondary wave modes that exhibit the desired conditions of synchronicity and large cross-energy transfer.These properties can be exploited in an actual system aimed at monitoring the structural condition of the waveguide by nonlinear waves detect defects, measure quasi-static loads or instability conditions, etc. .

Figure 1 :
Figure 1: Generic nonlinear waveguide finite element mesh just on the cross section with second harmonic generation highlighted.

Figure 2 :
Figure2: a Geometry and finite element mesh adopted for the railroad track nonlinear analysis.b Phase-velocity dispersion curve in the 0 ÷ 200 kHz frequency range with selected combinations of primary and secondary modes pinpointed.

Figure 3 :
Figure 3: Propagative modes in the 80 ÷ 160 kHz frequency range.a Flexural vertical mode energy mainly concentrated in the rail's head .b Flexural horizontal mode energy exclusively, confined in the rail's web .c Axial mode.d Complex mode involving a mixture of axial, torsional, and flexural displacements color online .

Figure 4 :
Figure 4: a Selected primary mode at 80 kHz.b Phase-matched synchronous but nonresonant secondary mode at 160 kHz.c Modal amplitude plot for propagative secondary modes.

Figure 5 :
Figure 5: a Selected primary mode at 80 kHz.b Resonant secondary mode at 160 kHz.c Modal amplitude plot for secondary propagative modes.

Figure 6 :
Figure6: a Geometry with associated mesh for the 2D periodic cell representative of the 12.7 mm thick HPPE plate dimensions in mm .b Phase-velocity dispersion curve in the 0 ÷ 500 kHz frequency range with primary and secondary modes for nonlinear analysis highlighted white circles .c Attenuation curve expressed in dB/m in the 0 ÷ 500 kHz frequency range with primary and secondary modes for nonlinear analysis highlighted white circles .

Figure 7 :
Figure 7: Modal amplitude plot for secondary propagative modes along with contour plots and 3D views of the selected primary and secondary modes for the viscoelastic HPPE plate color online .

Table 1 :
Material properties assumed for the railroad track analysis.