Elastic wave simulation based on modal excitation in 3 D medium

Paper describes a simulation of an elastic wave in homogenous isotropic waveguide with generic cross-section using semi-analytical finite element (SAFE) formulation. The wave is considered a travelling displacement field in the waveguide as a result to forcing excitation and is expressed via modes’ superposition. The solutions for modes are obtained solving SAFE governing eigen problem. Attenuation of the medium in SAFE framework differently from prior researchers is simulated via Rayleigh damping. It is shown, that severe damping is not properly supported by the SAFE formulation and revision for properly accepting linear viscosity is needed.


Introduction
The guided elastic waves are widely used in non-destructive long-range inspection for the elongated structures (waveguides).Simulation of propagating waves in the rods of the carbon fiber reinforced plastics (CFRP) waveguide enabled to detect the delamination regions in the samples [1], to test the structural integrity of pipes and rails [2], to measure the density of fluids [3], etc.
For characterization of guided waves the dispersion curves are used, which display the relationships of wave velocities against frequencies.Simultaneously, each mode shape of the propagating wave has a characteristic distribution of oscillatory displacements over the cross-section of a waveguide.The dispersion curves also provide information about phase and group velocities of the waves at particular frequency.This information provides the basis for the design of angle-and comb-type transducers used for ultrasonic inspection [4].
The first theoretical approach to guided waves dates back to Pochhammer and Chree who were the first researchers, who obtained the solution in the form of dispersion curves of the rod waveguides [5,6].The dispersion curves of homogenous and isotropic plates can be analytically derived from Rayleigh -Lamb equation [7].Laminated composite structures of complex geometry require numerical calculation techniques such as the finite element method (FEM).Many authors, among them Gavric [8], Hladky-Hennion [9], Damljanovic and Weaver [10], Zienkiewicz [11] have researched on FEM in application for the propagating waves in waveguides.An advantageous computational technique for investigating waves in the infinite-length waveguides is the semi analytical finite element method (SAFEM).SAFE computational costs are much lower when compared with the FEM implementation of the similar computational model.SAFE formulation involves two-dimensional discretization of the cross-section of the waveguide.Complex exponential functions ensure the harmonic propagating wave solutions along the waveguide.The SAFE approach was first used by Lagasse and Aalami [12,13] and further developed by Viola et al. [14], Hayashi et al. [4] and Ahmad [15].Along with obtaining the dispersion curves, the extensions of SAFE approach may simulate the forced response of the waveguide excited by the piezoelectric transducer at one end [2].
As a rule, in most investigations SAFE provided the dispersive solutions of waves in undamped waveguides of arbitrary cross-sections.The approach works for weekly linearly damped waves as well.The propagating waves in viscoelastic isotropic and orthotropic plates were extensively studied by Scalea [16].Two different models of linear viscoelasticity, hysteretic and Kelvin-Voigt, were incorporated in SAFE in order to extract viscoelastic dispersion solutions.Time dependent response to force excitation in isotropic viscoelastic waveguide was presented by Marzani and Bartoli [17].
In this work, we analyzed the free and forced propagating wave response of the waveguide subjected to Rayleigh damping by employing SAFE method.The focus was to analyze the behavior of the solutions as the damping is increased, as well as, to verify the obtained solutions against the ones simulated by full 3D FEM model.The forced response was presented as the superposition of the modes obtained SAFEM.

Governing equation of simulation
The detailed mathematical framework of SAFE is presented by following Gavric and Viola [8,14].Consider a 3D waveguide, which is uniform and infinite in 0 direction with stress-free surface.The vectors of displacements , Voigt's stresses and strains at each point of the cross-section are presented as: The generalized Hook's law is used for determining the 3D volumetric stress-strain relation as: where is the stiffness tensor.
The small strain against displacement relation reads as: where: In SAFE formulation the displacement field along the waveguide is assumed as harmonic propagating wave described analytically by complex exponential function (  ) , where [rad/s], [rad/m] are the angular frequency of the wave and the wavenumber correspondingly.The cross-section of the waveguide is discretized into 2D first-order Serendipity finite elements (FE).
Approximate the displacement field within the FE as: where is the nodal displacements vector, containing displacements of each node in the directions of the three axes, ( , ) is the matrix of approximation functions; is the node number within the FE; , , are the nodal displacements of th node in 0 , 0 , 0 directions correspondingly, and stands for imaginary units.
By substituting Eq. ( 4) into Eq.( 2), the strain-displacement relation reads as: where: Applying Hamilton's virtual work principle for a single FE yields: denotes matrix conjugate transpose and stands for a virtual quantity.The terms in Eq. ( 7) can be expanded to: The complex wavenumber = ( ) + ( ) ⋅ describes the spatial configuration of the wave propagating along the waveguide.( ) means the exponent of wave decay in space.In practical applications the attenuation is preferably measured in Np/m (Nepers per meter, 1Np = 8.6859 dB), that is, 1/ ( ) is a distance over which the amplitude of travelling wave decreases by 1/ .As Eq. ( 8) and ( 9) are complex quantities, the conjugate transposes of matrices are used, where the transposition converts the complex numbers to their conjugates, as: where * is the complex conjugate of .By substituting Eq. ( 10) and (11) into Eq.( 7) we obtain the finite element equation of the FE of the waveguide cross-section as: The FE matrices read as: where is the mass density, is the mass matrix, is stiffness matrix related with planar deformations of the cross-section, is symmetric stiffness matrix related with out-of-plane deformations, is the skew symmetric stiffness matrix, which couples planar and out-of-plane effects.
The integrals in Eq. ( 12) are computed by using the Gaussian quadrature rule.Since term ( ) ⁄ can be arbitrary (infinite length in direction), thus SAFE eigenvalue equation for FE becomes: As external excitation in terms of cross-section tractions is given, SAFE equation reads as: where is a nodal vector of external forces acting upon FE throughout all the length of the waveguide given as a propagating excitation wave.
As SAFE matrices and vectors are assembled to SAFE structural matrix and vector, the structural SAFE equation is obtained as: where is global vector of unknown nodal displacements over the cross-section.
In case = 0 the homogeneous version of Eq. ( 15) is obtained as which is treated as eigenproblem with or treated as the unknowns.

Dispersion curves
The solution of complex eigenvalue problem (Eq.(15a)) at numerous values within the selected range provides a number of propagating wave modes.The modes are presented in terms of , pairs of the wavenumber and oscillatory displacements shape over the cross section.The number of obtained modes at each value depends on the refinement of the FE mesh over the cross-section of the waveguide.Each modal displacement shape defines the type of the wave.The relationship of the wavenumbers of the same type of the wave against the frequency values is referred to as dispersion relationship (dispersion curve) of the particular type of the wave.Dispersion relationships of real waveguides are all non-linear.
Wavenumbers may be real, complex or pure imaginary values.Real values represent propagating waves.In such case * = , therefore Eq. (15a) can be simplified as: Complex wavenumbers represent damped waves.Complex wavenumbers are obtained from Eq. ( 16), however, they are mathematically correct only in case ( ) ≫ ( ) , until relation * ≈ may be assumed as approximately valid.In a general case the matrix expression in Eq. ( 16) has to be adjusted by taking into account real and imaginary parts of * , thus generating a non-linear eigenvalue problem.
Pure imaginary values ( * = − ), or ( ) ≪ ( ), represent end-mode solutions as "non-propagating waves".In such case Eq. ( 16) reads as: In this paper we are using eigenvalue problem formulation Eq. ( 17) for obtaining solutions, which are able to approximately the weakly-evanescent guided waves.
Equation Eq. (15a) is a quadratic eigenvalue problem with respect to .It may be transformed into the first order problem as: where: The homogeneous version of equation Eq. ( 19) reads as: where for each given real value 2×Total Degrees of freedom = 2 eigenvalues may be obtained.
Simultaneously, we formulate the left eigenvalue problem as: For representation purposes the set numerous obtained EWN was classified in order to distinguish the propagating, evanescent and end-mode waves form each other.As the EWN are always calculated with certain numerical error, the symmetry of the values over the complex plane implies that arithmetic averages of real and imaginary parts should be calculated.Several similarity measures are used in order to determine the type of symmetry of the obtained EWN.

Phase and group velocities
Along with the dispersion relations, phase velocity = ( ) ⁄ is another way for expressing relations between the wavelength, which is the inverse of real part of , and frequency.as the frequency increases.Velocity coincides with the velocity of mode at zero frequency.Phase velocities of higher order modes velocities approach at higher frequencies, Fig. 2(b).Phase velocity of the fundamental torsional mode , which in case of plates and cylindrical bars coincides with and is nondispersive, reaches value of ≈ 0.92 [4].However, over higher frequency range it does not approach the horizontal line of any more, thus implying low dispersion, Fig. 2(b).The phase velocities of , tend towards due to their numerical errors.The group velocity represents the velocity at which the mechanical energy travels in lossless waveguide and is defined as = ⁄ .Tracking separate mode over frequency axis might be a challenging task, since many modes with close wavenumbers exist, it is evident in Fig. 1(a).Therefore, usage of direct numerical differentiation + Δ + Δ ⁄ might produce misleading results by losing observed mode.
For given the group velocity of th mode could be derived from homogenous equation Eq. (15a) by taking the derivative with respect to wavenumber as: Pre-multiplying Eq. ( 19) by the left eigenvector gives: Eventually: Fig. 3 shows the results obtained by applying Eq. ( 21) for the waveguide with 4×4 FEs mesh.It may be observed that that modes with negative group velocities exist, which indicate that the mechanical energy of travels in direction opposite to that of phase velocity at particular frequency ranges.Unlike in plates and cylindrical bars, the group velocity of mode also has low dispersion.

Handling force response
The th mode and satisfies the equation: where the right eigenvector has the clear physical interpretation as the displacement profile of the propagating wave over the cross-section of the waveguide.The eigensolution and is obtained as a solution of the left eigenproblem as: Pre-multiplying of Eq. ( 22) by and post-multiplying Eq. ( 23) by yields following relations: This expands to: Which is equivalent to left and right eigenvector orthogonality relations: Any vector of length can be casted in basis of independent vectors.In case of external loading presented in Eq. ( 18) as , the response is presented as the weighted superposition of eigenvectors: And Eq. ( 18) can now be presented as: Pre-multiplying of Eq. ( 28) by left eigenvector yields: Because of orthogonality of the eigenvectors, the scalar weight of th eigenvector reads as: And finally, the nodal displacement vector is: The force is considered to be applied at particular point on the 0 axis over the cross-section of the waveguide as nodal amplitudes vector .Weight of the th mode contribution to overall solution can now reads as: Substitution of Eq. (32) to Eq. ( 31) yields: Since ( ) is a Fourier transform with respect to wavenumber , the application of Cauchy residue theorem ( [18]) leads to: Displacements expressed by relation Eq. (34) at points > are expressed via superposition of modes propagating forwards.At points < in relation Eq. (34) the backwards propagating modes are engaged.Relation Eq. (34) can also be considered as a response to general excitation, while both types of modes (propagating and evanescent) are engaged.Finally, displacement vector, as response to a harmonic load with a particular angular frequency , is obtained as: It is worth noting that calculation of left eigenvectors for lossless waveguides is not necessary.They can be substituted by the right ones as both and are Hermitian matrices ( = and = ) because of the symmetry of matrices , , and skew-symmetry of matrix .Under such conditions Eq. ( 23) may be re-casted as: For propagating modes corresponding to real valued EWN left eigenvectors can be obtained as = , since * = .In case of complex EWN ( ) + ( ) ⋅ and 0 + ( ) ⋅ the left eigenvectors are conjugate transposes of the right ones as ( ) − ( ) ⋅ and 0 − ( ) ⋅ respectively.Graphical point of feasibility for mentioned substitutions is clearly depicted in Fig. 1(b).In the presence of damping the technique for overcoming the necessity of left eigenvector calculation has been proposed by Gavric [8] and Hladky -Hennion [9] and by Viola [14].They demonstrated that in SAFE formulations the left eigenvectors can be always obtained from the right ones as: In case of the general time law, the excitation can be presented as a superposition of numerous Fourier components.In such case the response to each single harmonic component is obtained as described in this section.The overall response equals the sum of time relationships corresponding to all single-frequency responses.

Experimental exploration
Consider aluminum waveguide of mass density = 2700 kg/m 3 , Young's modulus = 7•10 9 Pa and Poisson's ratio = 0.33.The cross-section of the waveguide was 0.001×0.0011m, which slightly different from a square in order to avoid multiple EWN was considered in simulation.Consider the external loading as shearing force applied along a source (Fig. 4(a)).The time law is given as Hanning-windowed 5 cycle sinus burst centered at 250 kHz (Fig. 4(b)).
Fig. 4(c) displays the frequency spectrum of the modulated loading impulse.It is obtained by applying the fast Fourier transform (FFT) to the time relationship shown in Fig. 4(b), which is followed by zero padding 50 times the duration of impulse.The threshold for the meaningful harmonic components is assumed as 0.05 %.FFT sampling frequency was 250 GHz.The sum of displacement vectors of all harmonic components of the response in accordance with relation Eq. (35) provides the overall response to the loading impulse.
Fig. 5(a) presents the resulting displacements in direction registered at the monitored point at 0.1 m from the source line.Certain wrap-around effects may be observed, such as the displacements in direction at the time moment just before the arrival of the impulse.Most probably they appear because of the artificial periodicity of the impulse due to finite duration of the time relationship processed by FFT.Similar discrepancies are visible near to the incoming slope of the impulse in direction, Fig. 6.Such effects can be reduced by increasing the FFT sampling frequency and the length of zero padding.Results of explicit simulation by means of the 3D FE model subjected to the same excitation scenario are shown in Fig. 5(b).The 3D mesh density was 4×4×100 elements up to the monitored point in , and directions.Reasonable agreement of the results provided by both models was demonstrated.
Fig. 6 shows nodal responses obtained by using all modes contributions (Fig. 6(a)) and propagative-only modes' contributions (Fig. 6(b)) to the response at points near to the source line.The results imply that evanescent and end-modes' contribution to response decays rapidly and could be neglected in the far-field response analysis.Similar findings were reported by Weaver and Damljanovic [10].

Rayleigh damping
The Rayleigh damping is one of the models of internal losses appearing during the wave propagation.It enables to analyze the effects of viscosity of the waveguide material.As the practical value of the analysis, the characterization of the material properties on the basis of wave propagation features, should be mentioned.
In case of Rayleigh damping the structural dynamic equation is supplemented by the damping term as: where = + is the damping matrix presented as weighted sum of mass and stiffness matrices.In SAFE formulation, the damping term appears as: where: ).
In this way complex mass and stiffness matrices are introduced similarly to Kelvin-Voigt attenuation model, which involves the frequency-dependent complex stiffness matrix.In practice, the Rayleigh damping coefficients for a particular structure are determined experimentally [19].By the computational experiments presented here we investigate how the dispersion curves of different wave modes depend on values and .The results are presented in Fig. 7, where the cross-section mesh of the waveguide was 4×4.In the damped case all the EWN are complex, so the obtained modes are all evanescent.The phase velocities are obtained by using the real parts of EWN as = ( ) ⁄ .The presence of damping severely modifies the appearance of phase velocities relationships.Due to non-zero and the branches of fundamental modes tend to bend-away from the Rayleigh's wave's velocities relationship Fig. 7(a).Simultaneously, the attenuation of higher modes is less affected as the frequency increases, Fig. 7(b).
Fig. 8 shows how the phase velocities of fundamental modes are affected by severe damping.The scenario of heaviest simulated damping clearly demonstrates the behavior, which may be hardly explicable.This may serve as indication that assumptions on which Eq. ( 17) was based probably lack a solid background, though they were used by many researchers before for analysis of slightly damped structures.Thus implies the need for further research in order to revise the solutions of Eq. (15a) by fully considering the non-linear eigenvalue problem.  of Rayleigh damping at selected frequency = 1 MHz.As damping increases, the EWN do not retain symmetries with respect to and axes any more.Fig. 10 presents the response frequency spectra at the monitored point in case of damped and undamped waveguides with 4×4 FE cross-section.As the damping increases, the number of dominant modes directly influencing the displacement field at the observed point decreases.aluminum waveguide with 3×3 mesh on cross-section at = 20 μs

Conclusions
The analysis of propagating waves in waveguides by SAFE method enables to obtain the travelling-wave responses in infinite or almost-infinite waveguide structures.The approach enabled to obtain the wave modes, as well as, the forced responses of the waveguides to harmonic and general excitations at much lower cost than the direct FEM simulation of a waveguide as solid 3D structure.The new aspect of this work is that the waves have been analyzed at significant level of Rayleigh damping.The influences of damping on dispersion and phase-velocity curves have been investigated.
An important accent is that the linear eigenvalue problem as the model for wave modes investigation is only very approximate in case of severe damping.The presented SAFE framework includes assumptions that allow to simplify the governing equation, which are conventionally accepted by many researchers.However, we have found evidences, that the governing equation does not properly support severe to moderate damping.The reason is that even in case the physical and geometrical behaviors of the waveguide are linear, the eigenvalue problem equation has to be considered as dependent on the eigenvalue and therefore as non-linear one.The detailed analysis remains for future investigations in order to extend the SAFE framework in order to match the general linear damping model.
Moreover, an important application may find the SAFE models, which include the possibilities of simulating the impact of complex waveguide-surrounding media on the modal properties of the waveguide.
though vector × has no immediate physical meaning.Further we use notations for first halves of vectors , as × and × , where = and = [ ].The eigen-wavenumber (EWN) values are obtained from Eq. (16) in pairs.Solutions with positive and negative ( ) parts represent propagating wave modes in forward and backward 0 directions correspondingly.Complex solutions ±( ( ) + ( ) ⋅ ) represent evanescent waves with the amplitude decaying along 0 .The solutions ±( ( ) − ( ) ⋅ ) have no physical meaning [14].Pairs of pure imaginary solutions with positive and negative ( ) values represent non-propagating end-modes.The immediate practical value has the propagating and evanescent wave solutions.For model verification purposes, dispersions curves of a steel waveguide with rectangular cross-section 0.00508×0.00508m were calculated.The material properties were: mass density = 7850 kg/m 3 , Young's modulus = 2•10 11 Pa and Poisson's ratio = 0.3.The FE mesh was 4×4 over the cross-section of the waveguide.Fig. 1(a) show ( ) solutions presented as dispersion curves for forward direction travelling waves, where = /(2 ) is frequency in cycles per second.The results agreed reasonably well with those published by Hayashi and were backed up by experimental data in [4].Fig. 1(b) shows the classified EWN for forward and backward travelling waves at = 0.2 MHz in the complex plane, where symmetry in-between wavenumbers are evident.Four forward -propagating modes were obtained as two identical pairs because of the square shape of the waveguide.a) b) Fig. 1. a) ( ) solutions with 0.2 MHz dotted line mark; b) Eigen wavenumbers in complex number plane at = 0.2 MHz

Fig. 2 (
a) and (b) show the phase velocities of propagative modes of the same waveguide presented by 10×10 and 6×6 FEs over the cross-section.Theoretical values of longitudinal , shear and bar = ⁄ (speed of sound for pressure wave) velocities are presented by dashed lines.Phase velocities of fundamental flexural and longitudinal modes approach Rayleigh wave velocity ≈ . .

Fig. 2 .
Phase velocity of square bar with mesh of a) 10×10 and b) 6×6 on cross-section

Fig. 4 .Fig. 5 . 6 .
a) Schematic representation of excitation; b) excitation's modulation; c) excitation's frequency spectrum a) b) Displacements of monitored point ( = 0.1 m distance from source in the waveguide) in 0 direction obtained using a) SAFEM and b) FEM a) General excitation in lossless waveguide b) Propagative modes excitation in lossless waveguide Fig. Nodal response in direction to a) general and b) propagative modes excitations at distances of 0.0001 m 0.0002 m and 0.0004 m from the source in axis

Fig. 7 .
and values impact to phase velocity a) and attenuation b) of waveguide with 4×4 mesh on the cross-section.

Fig. 8 .
Fig. 8. Damping impact to phase velocity of fundamental modesin the waveguide with 4×4 mesh on the cross-section Fig.9shows the appearance of EWN values on the complex plane in case of different severity

Fig. 11 Fig. 9 .Fig. 10 .Fig. 11 .Fig. 12 .
demonstrates the dispersion effect, which is caused by different phase velocities of the modes excited by the loading in lossless and damped ( = 500 and = 10 -8 ) aluminum waveguide.Wavenumbers in complex plane at = 1 MHz with damping coefficients: a) = 100 and = 10 -9 ; b) = 500 and = 10 -8 ; c) = 1000 and = 10 -Frequency spectrums of the time transient response in 0 direction of the observed point at the distance of = 0.1 m from the forcing source in a) undamped and b) damped, c) Nodal response to force in lossless and damped waveguides at distances a) = 0.07 m, b) = 0.14 m and c) =0.21 m in 0 direction from sourceThe 3D view of the waveguide structure deformed by the wave is presented in Fig.12.The displacement field is shown at time = 20 μs in case of lossless and damped ( = 500 and = 10 -7 ) aluminum waveguide.Separate layers over the cross-section have been drawn.High sampling frequency has been used in order to avoid aliasing in animation.a) b) 3D animation representing excitation in a) lossless and b) damped ( = 500 and = 10 -7 )