Elastic wave simulation based on modal excitation in 3D medium
Rimantas Barauskas^{1} , Audrius Nečiūnas^{2} , Martynas Patašius^{3}
^{1, 2, 3}Kaunas University of Technology, Kaunas, Lithuania
^{1}Corresponding author
Journal of Vibroengineering, Vol. 18, Issue 8, 2016, p. 53215336.
https://doi.org/10.21595/jve.2016.18101
Received 29 September 2016; received in revised form 5 December 2016; accepted 13 December 2016; published 31 December 2016
JVE Conferences
Paper describes a simulation of an elastic wave in homogenous isotropic waveguide with generic crosssection using semianalytical 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.
Keywords: semianalytical finite element, time transient response, Rayleigh damping, complex eigen wavenumbers, evanescent and propagative modes.
1. Introduction
The guided elastic waves are widely used in nondestructive longrange 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 crosssection 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 combtype 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], HladkyHennion [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 infinitelength 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 crosssection 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 crosssections. 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 KelvinVoigt, 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.
2. 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 $0Z$ direction with stressfree surface. The vectors of displacements $\mathbf{u}$, Voigt’s stresses $\sigma $ and strains $\epsilon $ at each point of the crosssection are presented as:
The generalized Hook’s law is used for determining the 3D volumetric stressstrain relation as:
where $\mathbf{D}$ 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 ${e}^{i(kz\omega t)}$, where $\omega $ [rad/s], $k$ [rad/m] are the angular frequency of the wave and the wavenumber correspondingly. The crosssection of the waveguide is discretized into 2D firstorder Serendipity finite elements (FE).
Approximate the displacement field within the FE as:
$={\left[\begin{array}{l}{\sum}_{a=1}^{4}{N}_{a}(x,y){{v}_{xa}}_{}\\ {\sum}_{a=1}^{4}{N}_{a}(x,y){{v}_{ya}}_{}\\ {\sum}_{a=1}^{4}{N}_{a}\left(x,y\right){{v}_{za}}_{}\end{array}\right]}_{e}{\stackrel{}{\mathbf{U}}}_{e}{e}^{i(kz\omega t)}=\mathbf{N}\left(x,y\right){\stackrel{}{\mathbf{U}}}_{e}{e}^{i\left(kz\omega t\right)},$
where ${\stackrel{}{\mathbf{U}}}_{e}$ is the nodal displacements vector, containing displacements of each node in the directions of the three axes, $\mathbf{N}(x,y)$ is the matrix of approximation functions; $a$ is the node number within the FE; ${v}_{x}$, ${v}_{y}$, ${v}_{z}$ are the nodal displacements of $a$th node in $0X$, $0Y$, $0Z$ directions correspondingly, and $i$ stands for imaginary units.
By substituting Eq. (4) into Eq. (2), the straindisplacement relation reads as:
where:
Applying Hamilton’s virtual work principle for a single FE yields:
$T$ denotes matrix conjugate transpose and $\delta $ stands for a virtual quantity. The terms in Eq. (7) can be expanded to:
The complex wavenumber $k=Re\left(k\right)+Im\left(k\right)\cdot i$ describes the spatial configuration of the wave propagating along the waveguide. $Im\left(k\right)$ means the exponent of wave decay in space. In practical applications the attenuation is preferably measured in Np/m (Nepers per meter, $1$Np = 8.6859 dB), that is, $1/Im\left(k\right)$ is a distance over which the amplitude of travelling wave decreases by $1/e$. 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:
$=\underset{S}{\int}\underset{z}{\int}\delta \left({\stackrel{}{\mathbf{U}}}_{e}^{T}\left({\mathbf{B}}_{1}^{T}i{k}^{*}{\mathbf{B}}_{2}^{T}\right){e}^{Im\left(k\right)z}{e}^{i\left(Re\right(k)z\omega t)}\right)$
$\times \mathbf{D}\left({\mathbf{B}}_{1}+ik{\mathbf{B}}_{2}\right){\stackrel{}{\mathbf{U}}}_{e}{e}^{Im\left(k\right)}{e}^{i\left(Re\right(k)z\omega t)}dzdS={\int}_{{V}_{e}}\delta {\mathbf{\epsilon}}_{e}^{T}\mathbf{D}\mathbf{\epsilon}d{V}_{e}$
$={\int}_{S}{\int}_{z}\delta \left({\stackrel{}{\mathbf{U}}}_{e}^{T}\left({\mathbf{B}}_{1}^{T}i{k}^{*}{\mathbf{B}}_{2}^{T}\right){e}^{i({k}^{*}z\omega t)}\right)\mathbf{D}\left({\mathbf{B}}_{1}+ik{\mathbf{B}}_{2}\right){\stackrel{}{\mathbf{U}}}_{e}{e}^{i(kz\omega t)}dzdS$
$=\frac{{e}^{2Im\left(k\right)z}}{2Im\left(k\right)}\delta {\stackrel{}{\mathbf{U}}}_{e}^{T}{\int}_{S}\left[{\mathbf{B}}_{1}^{T}DB{}_{1}{}^{}i{k}^{*}{\mathbf{B}}_{2}^{T}\mathbf{D}{\mathbf{B}}_{1}+ik{\mathbf{B}}_{1}^{T}\mathbf{D}\mathbf{B}{}_{2}{}^{}+{k}^{*}k{\mathbf{B}}_{2}^{T}\mathbf{D}{\mathbf{B}}_{1}\right]dS{\stackrel{}{\mathbf{U}}}_{e},$
$\times {\stackrel{}{\mathbf{U}}}_{e}{e}^{Im\left(k\right)}{e}^{i\left(Re\right(k)z\omega t)}dzd{S}_{e}=\frac{{e}^{2Im\left(k\right)z}}{2Im\left(k\right)}{\omega}^{2}\delta {\stackrel{}{\mathbf{U}}}_{e}^{T}{\int}_{S}{\left(\mathbf{N}\left(x,y\right)\right)}^{T}{\rho}_{e}\mathbf{N}\left(x,y\right)dS{\stackrel{}{\mathbf{U}}}_{e},$
where ${k}^{*}$ is the complex conjugate of $k$. By substituting Eq. (10) and (11) into Eq. (7) we obtain the finite element equation of the FE of the waveguide crosssection as:
The FE matrices read as:
${\mathbf{K}}_{2e}={\int}_{x}{\int}_{y}{\mathbf{B}}_{2}^{T}\mathbf{D}{\mathbf{B}}_{2}dxdy,{\mathbf{K}}_{3e}={\int}_{x}{\int}_{y}{\mathbf{B}}_{2}^{T}\mathbf{D}{\mathbf{B}}_{2}dxdy,$
where ${\rho}_{e}$ is the mass density, ${\mathbf{M}}_{e}$ is the mass matrix, ${\mathbf{K}}_{1e}$ is stiffness matrix related with planar deformations of the crosssection, ${\mathbf{K}}_{3e}$ is symmetric stiffness matrix related with outofplane deformations, ${\mathbf{K}}_{2e}$ is the skew symmetric stiffness matrix, which couples planar and outofplane effects.
The integrals in Eq. (12) are computed by using the Gaussian quadrature rule. Since term ${e}^{2Im\left(k\right)z}\delta {\stackrel{}{\mathbf{U}}}_{e}/2Im\left(k\right)$ can be arbitrary (infinite length in $OZ$ direction), thus SAFE eigenvalue equation for FE becomes:
As external excitation in terms of crosssection tractions is given, SAFE equation reads as:
where ${\stackrel{}{\mathbf{F}}}_{e}$ 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 $\stackrel{}{\mathbf{U}}$ is global vector of unknown nodal displacements over the crosssection.
In case $\stackrel{}{\mathbf{F}}=0$ the homogeneous version of Eq. (15) is obtained as
which is treated as eigenproblem with $\omega $ or $k$ treated as the unknowns.
3. Dispersion solutions for eigenvalue problem
3.1. Dispersion curves
The solution of complex eigenvalue problem (Eq. (15a)) at numerous $\omega $ values within the selected range provides a number of propagating wave modes. The modes are presented in terms of $k$, $U$ pairs of the wavenumber and oscillatory displacements shape over the cross section. The number of obtained modes at each $\omega $ value depends on the refinement of the FE mesh over the crosssection 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 nonlinear.
Wavenumbers $k$ may be real, complex or pure imaginary values. Real values $k$ represent propagating waves. In such case ${k}^{*}=k$, therefore Eq. (15a) can be simplified as:
Complex wavenumbers $k$ represent damped waves. Complex wavenumbers are obtained from Eq. (16), however, they are mathematically correct only in case $Re\left(k\right)\gg Im\left(k\right)$ , until relation ${k}^{*}\approx k$ 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 ${k}^{*}$, thus generating a nonlinear eigenvalue problem.
Pure imaginary $k$ values $\left({k}^{*}=k\right)$, or $Re\left(k\right)\ll Im\left(k\right)$, represent endmode solutions as “nonpropagating 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 weaklyevanescent guided waves.
Equation Eq. (15a) is a quadratic eigenvalue problem with respect to $k$. 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 $\omega $ value 2×Total Degrees of freedom = 2$N$ eigenvalues $k$ may be obtained.
Simultaneously, we formulate the left eigenvalue problem as:
though vector ${\mathbf{W}}_{1\times 2N}$ has no immediate physical meaning.
Further we use notations for first halves of vectors $\mathbf{V}$, $\mathbf{W}$ as ${\mathbf{v}}_{N\times 1}$ and ${\mathbf{w}}_{N\times 1}$, where ${\mathbf{V}}_{j}=\left[\begin{array}{c}{\mathbf{v}}_{j}\\ k{\mathbf{v}}_{j}\end{array}\right]$ and ${\mathbf{W}}_{j}=\left[\begin{array}{ll}{\mathbf{w}}_{j}& {k}_{j}{\mathbf{w}}_{j}\end{array}\right]$.
The eigenwavenumber (EWN) values are obtained from Eq. (16) in pairs. Solutions with positive and negative $Re\left(k\right)$ parts represent propagating wave modes in forward and backward $0Z$ directions correspondingly. Complex solutions $\pm \left(Re\right(k)+Im(k)\cdot i)$ represent evanescent waves with the amplitude decaying along $0Z$. The solutions $\pm \left(Re\right(k)Im(k)\cdot i)$ have no physical meaning [14]. Pairs of pure imaginary solutions with positive and negative $Im\left(k\right)$ values represent nonpropagating endmodes. The immediate practical value has the propagating and evanescent wave solutions.
For model verification purposes, dispersions curves of a steel waveguide with rectangular crosssection 0.00508×0.00508 m were calculated. The material properties were: mass density $\rho =$ 7850 kg/m^{3}, Young’s modulus $E=$ 2·10^{11} Pa and Poisson’s ratio $\nu =$ 0.3. The FE mesh was 4×4 over the crosssection of the waveguide. Fig. 1(a) show $k\left(f\right)$ solutions presented as dispersion curves for forward direction travelling waves, where $f=\omega /\left(2\pi \right)$ 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 $f=$ 0.2 MHz in the complex plane, where symmetry inbetween wavenumbers are evident. Four forward propagating modes were obtained as two identical pairs because of the square shape of the waveguide.
Fig. 1. a) $k\left(f\right)$ solutions with 0.2 MHz dotted line mark; b) Eigen wavenumbers in complex number plane at $f=$ 0.2 MHz
a)
b)
For representation purposes the set numerous obtained EWN was classified in order to distinguish the propagating, evanescent and endmode 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.
3.2. Phase and group velocities
Along with the dispersion relations, phase velocity ${c}_{ph}=\omega /Re\left(k\right)$ is another way for expressing relations between the wavelength, which is the inverse of real part of $k$, and frequency. 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 crosssection. Theoretical values of longitudinal ${c}_{L}$, shear ${c}_{S}$ and bar ${c}_{b}=\sqrt{E/\rho}$ (speed of sound for pressure wave) velocities are presented by dashed lines. Phase velocities of fundamental flexural ${F}_{0}$ and longitudinal ${L}_{0}$ modes approach Rayleigh wave velocity ${c}_{R}\approx \frac{0.862+1.14\nu}{1+\nu}{c}_{S}$ as the frequency increases. Velocity ${c}_{b}$ coincides with the velocity of ${L}_{0}$ mode at zero frequency. Phase velocities of higher order modes velocities approach ${c}_{L}$ at higher frequencies, Fig. 2(b). Phase velocity of the fundamental torsional mode ${T}_{0}$, which in case of plates and cylindrical bars coincides with ${c}_{S}$ and is nondispersive, reaches value of ${c}_{torsional}\approx \text{0.92}{c}_{S}$ [4]. However, over higher frequency range it does not approach the horizontal line of ${c}_{S}$ any more, thus implying low dispersion, Fig. 2(b). The phase velocities of ${L}_{0}$, ${F}_{0}$ tend towards ${c}_{S}$ due to their numerical errors.
Fig. 2. Phase velocity of square bar with mesh of a) 10×10 and b) 6×6 on crosssection
a)
b)
The group velocity represents the velocity at which the mechanical energy travels in lossless waveguide and is defined as ${c}_{gr}=\partial \omega /\partial k$. 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 $\omega +\mathrm{\Delta}\omega /k+\mathrm{\Delta}k$ might produce misleading results by losing observed mode.
For given $\omega $ the group velocity of $j$th mode could be derived from homogenous equation Eq. (15a) by taking the derivative with respect to wavenumber as:
Premultiplying Eq. (19) by the left eigenvector $\mathbf{w}$ 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 ${T}_{0}$ mode also has low dispersion.
Fig. 3. Group velocity
4. Handling force response
The $j$th mode ${k}_{j}$ and ${V}_{j}$ satisfies the equation:
where the right eigenvector ${\mathbf{V}}_{j}$ has the clear physical interpretation as the displacement profile of the propagating wave over the crosssection of the waveguide.
The eigensolution ${k}_{m}$ and ${\mathbf{W}}_{m}$ is obtained as a solution of the left eigenproblem as:
Premultiplying of Eq. (22) by ${\mathbf{W}}_{m}$ and postmultiplying Eq. (23) by ${\mathbf{V}}_{j}$ yields following relations:
This expands to:
Which is equivalent to left and right eigenvector orthogonality relations:
Any vector of length $l$ can be casted in basis of $l$ independent vectors. In case of external loading presented in Eq. (18) as $\mathbf{P}$, the response is presented as the weighted superposition of eigenvectors:
And Eq. (18) can now be presented as:
Premultiplying of Eq. (28) by left eigenvector ${\mathbf{W}}_{m}$ yields:
Because of orthogonality of the eigenvectors, the scalar weight of $m$th eigenvector reads as:
And finally, the nodal displacement vector is:
The force is considered to be applied at particular point ${z}_{0}$ on the $0Z$ axis over the crosssection of the waveguide as nodal amplitudes vector $\stackrel{}{\mathbf{F}}$. Weight of the $m$th mode contribution to overall solution can now reads as:
Substitution of Eq. (32) to Eq. (31) yields:
Since $\mathbf{U}\left(z\right)$ is a Fourier transform with respect to wavenumber $k$, the application of Cauchy residue theorem ([18]) leads to:
Displacements expressed by relation Eq. (34) at points $z>{z}_{0}$ are expressed via superposition of $N$ modes propagating forwards. At points $z<{z}_{0}$ 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 $\omega $, 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 $\mathbf{A}$ and $\mathbf{B}$ are Hermitian matrices (${\mathbf{A}}^{T}=\mathbf{A}$ and ${\mathbf{B}}^{T}=\mathbf{B}$) because of the symmetry of matrices ${\mathbf{K}}_{1}$, ${\mathbf{K}}_{3}$, $\mathbf{M}$ and skewsymmetry of matrix ${\mathbf{K}}_{2}$. Under such conditions Eq. (23) may be recasted as:
For propagating modes corresponding to real valued EWN left eigenvectors can be obtained as $\mathbf{W}={\mathbf{V}}^{T}$, since${k}^{*}=k$. In case of complex EWN $Re\left(k\right)+Im\left(k\right)\cdot i$ and $0+Im\left(k\right)\cdot i$ the left eigenvectors are conjugate transposes of the right ones as $Re\left(k\right)Im\left(k\right)\cdot i$ and $0Im\left(k\right)\cdot i$ 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 singlefrequency responses.
5. Experimental exploration
Consider aluminum waveguide of mass density $\rho =$ 2700 kg/m^{3}, Young’s modulus $E=$7·10^{9} Pa and Poisson’s ratio$\nu =$ 0.33. The crosssection of the waveguide was 0.001×0.0011 m, 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 Hanningwindowed 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 $OZ$ direction registered at the monitored point at 0.1 m from the source line. Certain wraparound effects may be observed, such as the displacements in $OZ$ 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 $OZ$ direction, Fig. 6. Such effects can be reduced by increasing the FFT sampling frequency and the length of zero padding.
Fig. 4. a) Schematic representation of excitation; b) excitation’s modulation; c) excitation’s frequency spectrum
a)
b)
c)
Fig. 5. Displacements of monitored point ($z=$ 0.1 m distance from source in the waveguide) in $0Z$direction obtained using a) SAFEM and b) FEM
a)
b)
Fig. 6. Nodal response in $OZ$ 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 $OZ$ axis
a) General excitation in lossless waveguide
b) Propagative modes excitation in lossless waveguide
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 $OX$, $OY$ and $OZ$ 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 propagativeonly modes’ contributions (Fig. 6(b)) to the response at points near to the source line. The results imply that evanescent and endmodes’ contribution to response decays rapidly and could be neglected in the farfield response analysis. Similar findings were reported by Weaver and Damljanovic [10].
6. 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 $\mathbf{C}=\alpha \mathbf{M}+\beta \mathbf{K}$ is the damping matrix presented as weighted sum of mass and stiffness matrices. In SAFE formulation, the damping term appears as:
$=\stackrel{}{\mathbf{F}}{e}^{i(kz\omega t)},$
where:
In this way complex mass and stiffness matrices are introduced similarly to KelvinVoigt attenuation model, which involves the frequencydependent 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 $\alpha $ and $\beta $ . The results are presented in Fig. 7, where the crosssection 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 ${c}_{ph}=\omega /Re\left(k\right)$.
The presence of damping severely modifies the appearance of phase velocities relationships. Due to nonzero $\alpha $ and $\beta $ the branches of fundamental modes tend to bendaway 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 nonlinear eigenvalue problem.
Fig. 7. $\alpha $ and $\beta $ values impact to phase velocity a) and attenuation b) of waveguide with 4×4 mesh on the crosssection.
a)
b)
Fig. 8. Damping impact to phase velocity of fundamental modes in the waveguide with 4×4 mesh on the crosssection
Fig. 9 shows the appearance of EWN values on the complex plane in case of different severity of Rayleigh damping at selected frequency $f=$ 1 MHz. As damping increases, the EWN do not retain symmetries with respect to $OX$ and $OY$ 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 crosssection. As the damping increases, the number of dominant modes directly influencing the displacement field at the observed point decreases.
Fig. 11 demonstrates the dispersion effect, which is caused by different phase velocities of the modes excited by the loading in lossless and damped ($\alpha =$ 500 and $\beta =$ 10^{8}) aluminum waveguide.
Fig. 9. Wavenumbers in complex plane at $f=1$MHz with damping coefficients: a) $\alpha =$ 100 and $\beta =$ 10^{9}; b) $\alpha =$ 500 and $\beta =$ 10^{8}; c) $\alpha =$ 1000 and $\beta =$ 10^{7}
a)
b)
c)
Fig. 10. Frequency spectrums of the time transient response in $0Z$ direction of the observed point at the distance of $z=$ 0.1 m from the forcing source in a) undamped and b) damped, c) waveguide
a)
b)
c)
Fig. 11. Nodal response to force in lossless and damped waveguides at distances a) $z=$ 0.07 m, b) $z=$ 0.14 m and c) $z=$0.21 m in $0Z$ direction from source
a)
b)
c)
The 3D view of the waveguide structure deformed by the wave is presented in Fig. 12. The displacement field is shown at time $t=$ 20 μs in case of lossless and damped ($\alpha =$ 500 and $\beta =$ 10^{7}) aluminum waveguide. Separate layers over the crosssection have been drawn. High sampling frequency has been used in order to avoid aliasing in animation.
Fig. 12. 3D animation representing excitation in a) lossless and b) damped ($\alpha =$ 500 and $\beta =$10^{7}) aluminum waveguide with 3×3 mesh on crosssection at $t=$ 20 μs
a)
b)
7. Conclusions
The analysis of propagating waves in waveguides by SAFE method enables to obtain the travellingwave responses in infinite or almostinfinite 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 phasevelocity 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 nonlinear 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 waveguidesurrounding media on the modal properties of the waveguide.
References
 Kazys R., et al. Application of ultrasonic guided waves for NDE of composite structures. Proceedings of the National Seminar and Exhibition on NonDestructive Evaluation, 2009, p. 236239. [Search CrossRef]
 Loveday P. W., Long C. S. Time domain simulation of piezoelectric excitation of guided waves in rails using waveguide finite elements. The 14th International Symposium on: Smart Structures and Materials and Nondestructive Evaluation and Health Monitoring, International Society for Optics and Photonics, 2007. p. 65290V. [Publisher]
 Fan Z., et al. Torsional waves propagation along a waveguide of arbitrary cross section immersed in a perfect fluid. The Journal of the Acoustical Society of America, Vol. 124, Issue 4, 2008, p. 20022010. [Publisher]
 HayashiT., Song W. J., Rose J. L. Guided wave dispersion curves for a bar with an arbitrary crosssection, a rod and rail example. Ultrasonics, Vol. 41, Issue 3, 2003, p. 175183. [Publisher]
 Pochhammer Lt. Ueber die fortpflanzungsgeschwindigkeiten kleiner schwingungen in einem unbegrenzten isotropen kreiscylinder. Journal für die Reine und Angewandte Mathematik, Vol. 81, 1876, p. 324336. [Search CrossRef]
 Chree C. The equations of an isotropic elastic solid in polar and cylindrical coordinates their solution and application. Transactions of the Cambridge Philosophical, Society, Vol. 1889, 14, p. 250369. [Search CrossRef]
 Victorov I. A. Rayleigh and Lamb Waves. Plenum Press, New York, 1967. [Publisher]
 Gavric L. Computation of propagative waves in free rail using a finite element technique. Journal of Sound and Vibration, Vol. 187, Issue 3, 1995, p. 531543. [Publisher]
 HladkyHennionA. C. Finite element analysis of the propagation of acoustic waves in waveguides. Journal of Sound and Vibration, Vol. 194, Issue 2, 1996, p. 119136. [Publisher]
 Damljanovic V., Weaver R. L. Forced response of a cylindrical waveguide with simulation of the wavenumber extraction problem. The Journal of the Acoustical Society of America, Vol. 115, Issue 4, 2004, p. 15821591. [Publisher]
 Zienkiewicz O. C., Taylor R. L. The Finite Element Method: Solid Mechanics. ButterworthHeinemann, 2000. [Search CrossRef]
 Lagasse P. E. Higher‐order finite‐element analysis of topographic guides supporting elastic surface waves. The Journal of the Acoustical Society of America, Vol. 53, Issue 4, 1973, p. 11161122. [Publisher]
 Aalami B. Waves in prismatic guides of arbitrary cross section. Journal of Applied Mechanics, Vol. 40, Issue 4, 1973, p. 10671072. [Publisher]
 Viola E., Marzani A., Bartoli I. Semianalytical Formulation for Guided Wave Propagation. Mechanical Vibration: Where do we Stand? Springer, Vienna, 2007, p. 105121. [Publisher]
 Ahmad Z. A. B., Gabbert U. Simulation of Lamb wave reflections at plate edges using the semianalytical finite element method. Ultrasonics, Vol. 52, Issue 7, 2012, p. 815820. [Publisher]
 Bartoli I., et al. Modeling wave propagation in damped waveguides of arbitrary crosssection. Journal of Sound and Vibration, Vol. 295, Issue 3, 2006, p. 685707. [Publisher]
 Marzani A., Bartoli I. High frequency waves propagating in octagonal bars: a low cost computation algorithm. Algorithms, Vol. 2, Issue 1, 2009, p. 227246. [Publisher]
 Renno J. M., Mace B. R. On the forced response of waveguides using the wave and finite element method. Journal of Sound and Vibration, Vol. 329, Issue 26, 2010, p. 54745488. [Publisher]
 Vasudeva S., et al. Estimation of elastic and damping characteristics of viscoelastically constrained carbon strands. Structural Dynamics and Materials Conference (SDM), 2006. [Publisher]
Cited By
Mathematical Modelling and Analysis
Audrius Nečiūnas, Martynas Patašius, Rimantas Barauskas

2018
