Semi-Analytical Finite Element Analysis of the Influence of Axial Loads on Elastic Waveguides

Guided wave ultrasound is a promising technology for non-destructive inspection and monitoring of long slender structures such as pipes and rails. These structures are effectively one-dimensional waveguides and therefore a large length of the structure can be inspected from a single transducer location. A good introduction to guided wave inspection is available in (Rose, 2002). During the design of guided wave inspection systems it is advantageous to understand the modes of propagation, to predict how these modes interact with the damage to be detected and to be able to transmit and receive these modes independently (Lowe et al., 1998). Analytical solutions for wave propagation in circular cylinders are well known (Gazis, 1959) and may be used when developing pipe monitoring systems. When the cross-sectional geometry is more complex, such as in the case of rails, it becomes necessary to employ numerical solutions.


Introduction
Guided wave ultrasound is a promising technology for non-destructive inspection and monitoring of long slender structures such as pipes and rails.These structures are effectively one-dimensional waveguides and therefore a large length of the structure can be inspected from a single transducer location.A good introduction to guided wave inspection is available in (Rose, 2002).During the design of guided wave inspection systems it is advantageous to understand the modes of propagation, to predict how these modes interact with the damage to be detected and to be able to transmit and receive these modes independently (Lowe et al., 1998).Analytical solutions for wave propagation in circular cylinders are well known (Gazis, 1959) and may be used when developing pipe monitoring systems.When the cross-sectional geometry is more complex, such as in the case of rails, it becomes necessary to employ numerical solutions.
It is possible to compute the wavenumber -frequency relations of propagating modes using conventional three -dimensional continuum finite element models, with appropriate boundary conditions, as was demonstrated by Thompson (1997).The advantage of this method is that it can be implemented using commercially available finite element codes, but requires significant user post processing, and the solution of many different model lengths.Alternately, semi-analytical finite elements (SAFE) can be specially formulated to efficiently analyse these problems.The formulation of these elements includes complex exponential functions to describe the wave propagation along the waveguide and conventional finite element interpolation functions across the cross-section.Therefore, only a two -dimensional mesh of the cross-section of the waveguide is required, which results in a significant reduction in the amount of computation required.These elements are sometimes called waveguide finite elements and have been implemented by a number of research groups (Aalami, 1973;Lagasse, 1973;Gavrić, 1995;Hayashi et al., 2003;Damljanović & Weaver, 2004a;Bartoli et al., 2006;Predoi Mihai et al., 2007;Castaings & Lowe, 2008 and Ryue et al., www.intechopen.comFinite Element Analysis -From Biomedical Applications to Industrial Developments 440 2008).In addition to computing dispersion properties of complex geometry waveguides, the SAFE approach has been applied to a number of other problems, including forced response computations (Damljanović & Weaver, 2004b), modelling of piezoelectric transducers attached to waveguides (Loveday, 2008), determining the power input to waveguides (Nilsson & Finnveden, 2007), and analysing the scattering of waves at discontinuities (Baronian et al., 2009).Authors use slightly different formulations and two popular formulations will be compared in Section 2.
It is well known that the application of an axial load to a beam can influence the natural frequencies of the beam.Similarly, the application of an axial load to a waveguide can influence the wave propagation characteristics of the waveguide.In the case of continuously welded rail, the axial load in the rail is an important parameter.These rails are installed to generally be in tension but the amount of tension depends on the temperature.If the rail goes into compression there is a danger of buckling, which can lead to derailments while if the tension is excessive, the probability of fatigue cracks is increased and this can lead to rail breaks and derailments.A guided wave system was developed to continuously monitor rails for breaks (Loveday, 2000) and it would be a valuable addition if the same system could monitor the rail for compression in the rail before buckling.There has been considerable research into the possibility of using guided waves to measure the axial load in rails.Low frequency flexural waves were investigated by Damljanović & Weaver (2005) who proposed to use a scanning laser vibrometer to measure displacements of points along the rail and then to extract the wavenumber of the flexural wave at 200Hz (Damljanović & Weaver, 2004b).This method requires that the rail be released from the sleepers for a considerable length.Chen & Wilcox (2006) investigated the use of higher frequency guided waves for measuring loads in rods.Simulation results clearly demonstrated that the phase velocity and group velocity are sensitive to changes in load.
The influence of the axial load on wave propagation characteristics may be analysed by three -dimensional finite element models, although the process is tedious and difficult at higher frequencies where numerous waves propagate (Chen & Wilcox, 2007).Recently the SAFE method was extended to include axial loads (Loveday, 2009) and it was demonstrated that the required modification to an existing SAFE code is trivial.This extension is described in Section 3 of this chapter.The use of SAFE to analyse the influence of axial loads on wave propagation in rails offers the advantage of computational efficiency.In addition because of the analytical nature of the method, it is possible to directly compute certain sensitivities that would otherwise have to be computed by a finite difference method (Loveday & Wilcox, 2010).These possibilities are described in Section 4 while numerical results are presented in Section 5.

Semi-analytical finite element formulation for elastic waveguides
The formulation presented here follows that of Gavrić (1995) who choose the complex exponential functions along the waveguide to have the axial terms phase shifted by π/2 relative to the in-plane terms.This choice leads directly to symmetric matrices and is described in Section 2.1.Many of the other SAFE formulations used do not result in a symmetric stiffness matrix, which is inconvenient for eigensolvers.These formulations will not be considered here.A different approach was adopted by Damljanović & Weaver (2004a) who used a transformation to obtain symmetric matrices.These matrices are the same as those obtained by Gavrić except that the sign of one matrix is reversed.A comparison of the two approaches is outlined in Section 2.2.

Formulation by Gavrić
The displacement field in an elastic waveguide, extending in the z direction, may be written as a complex exponential along the waveguide and a finite element approximation over the cross-section.The displacement fields (u, v, w) employed by are Gavrić (1995) were chosen to take the form: where, z is the coordinate in the direction along the waveguide, κ the wavenumber and ω the frequency.u(x,y), v(x,y) and w(x,y) are the interpolated displacements in the x, y and z directions respectively.The strain energy of an infinitesimal element of the waveguide is, where, * denotes complex transpose, 0 k , 1 k and 2 k are defined below,  and c are the strain and elasticity matrices defined in standard form as: where ν is Poisson's ratio and E is elastic modulus.The strains due to the displacement field in (1) are separated into a function over the cross-section, multiplied by the complex exponential.Furthermore, for convenience the terms that are dependent on the wavenumber are separated from those independent of wavenumber.The result is: The terms in the strain energy, which are independent, linearly and quadratically dependent can now be expressed as follows: Four-noded isoparametric elements were used in this research although quadratic elements are more efficient (Andhavarapu et al., 2010).The standard interpolation matrix N, which is a function of local coordinate (,)   , relates the nodal degrees of freedom to the coordinate and displacement distributions as follows: Where u is the vector of nodal displacements.The strains may now be written,

 
by the Jacobian , J , as usual.Integration is performed over the area of the elements to give the elemental matrices.
The element matrices are assembled into the system equations of motion.For free vibration the equation of motion is: where capitals indicate the assembled counterpart of the lowercase elemental matrices.It should be noted that the mass ( M ) and stiffness ( 012 ,, KKK ) matrices are all real and symmetric.

Comparison with formulation by Damljanović and Weaver
Damljanović & Weaver (2004a) employed the following functions to describe the displacement field: After applying a procedure similar to that in Section 2.1 they obtain the system of equations of motion: In this case the matrix 1 K  is skew symmetric and they use the transformation matrix T to transform this matrix to a symmetric matrix.The matrix T has the form, and can be thought of as a transformation from the physical coordinates to a transformed coordinate system, The transformation has no effect on the symmetric matrices, but the skew symmetric matrix is transformed to a symmetric matrix.The equations of motion in the transformed coordinates are then, In order to compare the two formulations it is noted that the displacement functions used by Gavrić (1) may be written as, If one were to start with the functions: and follow Garvić's method, the transformed matrices of Damljanović & Weaver would be obtained directly.

Extension of SAFE to include axial load
The presence of an initial stress or load introduces additional terms in the strain energy which therefore lead to additions to the stiffness matrix.The SAFE method is extended to include axial loads in one-dimensional waveguides in this section.
When analysing small amplitude elastic waves we generally use the linear infinitesimal strain -displacement relations.However, when initial finite strains are present in the structure, it is necessary to make use of the full strain-displacement relation (Rose, 1999).
The linear (  ) and full (E) strain -displacement relationships may be written as, where u, v and w are displacements in the x, y and z directions respectively, as before.
The potential energy per unit volume, k, may be written as the sum of the strain energy associated with the small amplitude elastic wave (  ) and the work performed by the initial stress ( (0) If we restrict our analysis to consider only an axial load on the waveguide, which extends in the z direction, we need only retain the initial stress (0) zz  .We can assume that the stresses associated with the small amplitude elastic wave are at least an order of magnitude smaller than the axial stress and therefore the product of these stresses and the non-linear strain terms are negligible.
The strain energy can then be written as, (0) All of these terms except the term containing the initial load are already included in the linear strain energy used in the SAFE method.The term containing the initial load can be expanded as follows, The term (0) disappears when the variation of the Hamiltonian of the waveguide is taken as in Gavrić (1995) or alternatively when the Lagrange equations are applied as in Damljanović & Weaver (2004a).This term represents the interaction of the initial stress with the linear strain, which does not feature in the equations of motion for linear systems.
Therefore the only term we need to add to the linear strain energy expression, previously used in the SAFE, is, Substituting the displacement interpolation functions and integrating as before produces the additional strain energy term with the form, The form of this term is identical to that of the kinetic energy.Therefore the additional stiffness matrix is proportional to the mass matrix and the equations of motion can be written simply as, where, As no new matrices have to be created it is trivial to extend existing software to analyse the influence of axial load on the wave propagation characteristics.

Computation of dispersion characteristics
If we consider free harmonic vibration ( 24) and ( 25) provide the eigenvalue problem, including the initial axial stress, To obtain the relationship between wavenumber and frequency it is necessary to specify one of these and solve the eigenvalue problem to obtain the other.If one is interested in the behaviour at a frequency or a range of frequencies these may be computed by complementing ( 26) with an identity as suggested by Hayashi et al. (2003) and then solving the following complex eigenvalue problem: The wavenumbers that are obtained, by solving this problem, can be real, imaginary or c o m p l e x a n d o c c u r i n p a i r s w i t h o p p o s i t e s i g n c o r r e s p o n d i n g t o w a v e s t r a v e l l i n g i n opposite directions.If the number of nodes in the model is denoted N, the eigensolution results in 6N eigenvalue-eigenvector pairs.
If only the propagating modes are required the real wavenumber may be specified and the eigenvalue problem ( 26) may be solved.The propagating modes have real frequency,  and real mode shape  .The dispersion characteristics, of the propagating modes, may be obtained by solving this eigenvalue problem for a range of different real wavenumbers and collecting the real frequencies that are produced.This approach is used here.At each specified wavenumber a set of frequency points are obtained but there is no relationship between the frequencies at one wavenumber and those at the next wavenumber.If we want to plot the dispersion curves it is necessary to track the modes from one wavenumber to the next.A technique was developed to track the modes, utilizing the orthogonality property of the mode shapes, as expressed in (28).Equation 28expresses the mass orthogonality of two arbitrary modes r and s but the stiffness orthogonality could have been used instead. ( It is reasonable to assume that if small wavenumber steps are taken, the mode shapes will not change significantly between steps, and that a mode shape at one wavenumber would almost be mass-orthogonal to those at the next step.The mass-orthogonality of the mode shapes at wavenumber step k, to those at wavenumber step k+1, is computed, i.e., If the wavenumber versus frequency curves have not crossed in this wavenumber interval, the largest terms in the matrix Θ will be the diagonal terms.The presence of an off-diagonal term that is larger than the corresponding diagonal term, indicates that the curves have crossed.For example, if the terms Θ i,i+1 and Θ i+1,i are larger than the terms Θ i,i and Θ i+1,i+1 this indicates that the i th mode at wavenumber step k becomes the i+1 th mode at wavenumber step k+1 and that the i+1 th mode at wavenumber step k becomes the i th mode at wavenumber step k+1.This is then taken into account in the numbering of the modes so that a curve for each mode can be plotted.If the complex eigenvalue problem ( 27) is solved at a range of frequencies a similar approach may be used to plot the dispersion curves of the propagating modes (Loveday, 2008).
The group velocities can be computed using an analytical expression similar to the one presented by Hayashi et al. (2003) for their element formulation.The analytical expression for the group velocity is obtained by differentiating the solution of the eigenvalue problem with respect to frequency,

www.intechopen.com
Finite Element Analysis -From Biomedical Applications to Industrial Developments 448 leading to an expression which can be rearranged to obtain the group velocity as the subject.
The analytical nature of the SAFE method also makes it possible to directly compute the sensitivity of the wavenumber to changes in axial load at a particular frequency as shown in (32).
In a similar way, the sensitivity of the wavenumber to changes in elastic modulus can be obtained It is interesting to note the similar form and the role of the group velocity in these sensitivities.This feature was used to compare the change in wavenumber that would be expected to result from a change in temperature via these two mechanisms (Loveday & Wilcox, 2010).

Results & discussion
The influence of axial load on the wave propagation in an aluminium rod was computed and compared to an analytical solution.The analytical solution is based on Euler -Bernoulli beam theory and was presented by Chen and Wilcox (2007).The phase velocity v ph , at frequency, ω, for a beam with Young's modulus, E, second moment of area, I, mass per unit length, m, and subjected to a tensile axial load, T, is: A 1 mm diameter aluminium rod was modelled with SAFE method without axial load and then with a tensile axial load corresponding to 0.1% axial strain.Frequency Thickness (Hz m) Phase Velocity / Shear Velocity 0% strain 0.1% strain Fig. 1.Influence of strain on 1mm diameter aluminium rod computed by Euler-Bernoulli beam theory (top) and SAFE (bottom).
The ability of the SAFE method to compute the influence of axial loads on complex shapes such as rails was demonstrated by computing the dispersion curves for a UIC60 rail profile.The model used a density of 7700 kg/m3, a Young's modulus of 215 GPa and a Poisson's ratio of 0.3.The axial load applied corresponds to 0.1% strain.This analysis was performed by setting the wavenumber and computing the frequencies of the propagating modes at this wavenumber.The wavenumber was increased in a number of steps and the orthogonality conditions were used to track the evolution of the modes.The wavenumber -frequency curves are shown in fig.2, which shows numerous propagating modes.The case with no axial load is shown with dashed lines while the solid lines represent the case with axial load.The first four propagating modes, which propagate at all frequencies, are identified in the legend.It is difficult to see the influence of the axial load when the curves are plotted over 50 kHz and 120 rad/m.The second plot shows a zoomed in view and the small influence can be observed in the fundamental horizontal and vertical bending modes in the frequency range between 5 kHz and 10 kHz.The influence of axial load on the group velocities was computed using (31) and is shown in fig. 3. Again, the influence is small and the zoomed in view is necessary to see the influence.
The group velocity of the two fundamental flexural modes appears to be most sensitive in the frequency range from 3 kHz to 6 kHz.In this range the group velocity changes rapidly with frequency and significant dispersion would be observed in experiments.Finally, the sensitivity of the wavenumbers of the propagating modes to the axial load was computed using (32).The sensitivities are plotted in fig. 4 along with the relative sensitivities, which are the sensitivities normalised by the wavenumber.These sensitivities indicate modes and frequencies where the axial load has an influence on the wavenumber or phase velocity more clearly than the curves in fig. 2. Again, the two fundamental flexural modes are the most sensitive.The relative wavenumber sensitivity provides a measure of how much the wavelength changes relative to the wavelength and if one considers transmitting and then receiving at a number of wavelengths further along the rail one could calculate the amount of phase change that would occur due to the axial load.The relative wavenumber sensitivity is highest at very low frequencies.At slightly higher frequencies the two fundamental flexural modes show maxima in the relative wavenumber sensitivities at approximately 5kHz and 6 kHz.It is interesting to note that these frequencies appear to correspond to minima in the group velocities of these two modes.This information would be useful when trying to design a system to exploit the wavenumber sensitivity to measure the axial load in a rail.

Conclusion
The SAFE method has become popular for analysing guided wave propagation in structures with complex geometries.The method, applied to one-dimensional waveguides, was extended to include the presence of an axial load.It was shown that the software modifications required for this extension are trivial.One of the advantages of the SAFE method is that it allows group velocities to be computed analytically.The analytical nature of the method can further be exploited to compute sensitivities analytically and the sensitivity of the wavenumber to axial load was computed in this manner.Results for propagating modes up to 50 kHz in UIC60 rail were computed.

Acknowledgment
This work was funded in part by the CSIR SRP Project -TB_2009_19, "Ultrasonic Transducers for Long-Range Rail Monitoring".
U are the displacements in the physical coordinates and t U are the displacements in the transformed coordinates.The transformation has the properties of the different transformations is that, in the transformed coordinates, the matrix 1 K in Garvić's method is equal to 1 K  in Damljanović & Weaver's formulation.The other matrices are the same.
The results obtained are compared to the Euler-Bernoulli beam solution in fig 1.It is clear that the results of the SAFE analysis are practically identical to the Euler-Bernoulli beam model results.This result confirms the accuracy of the formulation and the numerical implementation of the method.www.intechopen.comSemi-AnalyticalFinite Element Analysis of the Influence of Axial Loads on Elastic

Fig. 2 .
Fig. 2. Influence of strain on dispersion in UIC60 rail over a large frequency range (top) and zoomed in (bottom).

Fig. 3 .
Fig. 3. Influence of strain on group velocity in UIC60 rail over a large frequency range (top) and zoomed in (bottom).
Finite Element Analysis -From Biomedical Applications to Industrial Developments Semi-Analytical Finite Element Analysis of the Influence of Axial Loads on Elastic Waveguides www.intechopen.com