Numerical solution of Rayleigh-Lamb frequency equation for real , imaginary and complex wavenumbers

Guided waves, especially Lamb waves or shear-horizontal waves, are widely used types of waves for ultrasonic inspection of large structures. Well known property of guided waves is their dispersive character, which means that the propagation velocity of the particular wave mode is not only a function of physical properties of the material, in which the wave propagates or the wave ́s frequency, but also depends on the geometry of the structure in itself. Dispersion curves provide us the information related to the dependency between the wavenumber and the frequency of the particular mode and can be obtained by a numerical solution of Rayleigh-Lamb frequency equation. A solution of Rayleigh-Lamb frequency equation forms for a given frequency and plate thickness a set of a finite number of real and pure imaginary wavenumbers and an infinite number of complex wavenumbers. Proposed paper presents a complete procedure of how to obtain all three kinds of wavenumbers for a given geometry and frequency interval. The main emphasis is placed on the effectiveness of the procedures, which are used for finding the roots of dispersion equation for all three kinds of wavenumbers.


Introduction
Ultrasonic testing is, besides eddy current testing method, one of most commonly used nondestructive techniques in the case of inspection of plate-like structures [1].Guided waves, especially Lamb waves, are ideal solution for inspection of large structures such as metal sheets or composite structures [2].The ability of the effective crack detection by use of various modes of Lamb wave has been demonstrated in several publications [2][3][4].The structure of the particular mode of guided Lamb wave in itself is relatively complex compared to the longitudinal or transversal bulk waves.This fact appears to be an advantage, since it is possible to select a particular mode shape in order to increase the sensitivity of such mode to specific types of flaws [3].This advantage, is however, directly related with increased complexity of the results interpretation.In most cases, one cannot perform further analyses without the numerical predictions, since the more modes are present in the structure, the more information about the flaw is available with rapid increase of the interpretation difficultness at the same time.
The main feature of guided Lamb waves is their dispersive character, defined by dispersion curves, which describe us the dependency between the wavenumber and the frequency of particular wave mode for given geometry.The dispersion curves in case of single or multi-layered media can be obtained numerically and, as well as in particular frequency range, also experimentally [3][4][5][6].Various methods of numerical approaches for dispersion curves computations has been published so far [3,7].The main scope of the article is to present fast and robust procedure for calculation a set of real, imaginary and complex wavenumbers, which are the roots of dispersion equation, for wide range of frequency thickness product.

Theory of Lamb Waves
The governing equation will be derived under plane strain condition with the assumption of traction-free surfaces (Fig. 1).

Fig. 1. Geometry of the free plate
The unknown displacement vector u will be defined with the use of the Helmholtz decomposition theorem [8,9]: where x and z are spatial coordinates, and is potential and vector function respectively and is differential operator.The equation of motion can be expressed in terms of displacements, resulting in following form: where and denotes Lame constants, is density of the material and ̈ denotes second time t derivate of displacement vector u.Substitution of Eq. 1 into Eq. 3 will result in two independent governing equations for longitudinal and transversal waves: where c L is the speed of longitudinal (dilatational) waves, c T is the speed of transversal (shear) waves.The solutions of Eqs. 4 and 5 can be assumed in following forms: (6) where k represents wavenumber, denotes angular frequency, functions and represent standing waves in z direction.Functions and (Eq.6) have to be further substituted into Eq. 4 and 5 in order to obtain final solutions for and : where and .
, , , denotes constants dependent of boundary conditions, and represents complex components along z axis, and represents wavenumbers of longitudinal and transversal wave and k represents wavenumber along direction of wave propagation.The displacements and stresses can be with help of Lame equations, Eq. 1 and Eqs.6-9 expressed in following form: Finally, Eqs. 10 and 11 can be split into two sets of modes according to odd or even functions about z = 0 axis [3].
The dispersion equation, or in other words Rayleigh-Lamb frequency equation, for both modes can be determined by applying the traction-free boundary condition (See Fig. 1), which results in: where +1 applies for modes of symmetric Lamb wave and -1 applies for modes of antisymmetric Lamb wave.

Solution of Rayleigh-Lamb Frequency Equation
As mentioned in chapter 1, the solution of Eq. 12 for given frequency will produce a set of finite number of propagating modes with real wavenumber, finite number of nonpropagating modes with imaginary wavenumber and an infinite number of inhomogeneous modes with complex wavenumber [1].In the following chapters are discussed individual approaches of solving the Rayleigh-Lamb frequency equation for individual types of wavenumbers.It has to be noted, that the entire code has been written in Matlab ® software.

Numerical solution for real wavenumbers
For solving the real wavenumbers, it is preferable to rewrite Eq. 12 in following form [3]: for symmetric modes and for antisymmetric modes.The procedure of finding the real roots of the dispersion equation is the same for both modes (Eqs.13 and 14).The robustness of the procedure is based on the vectorization, which is incorporated in Matlab ® software.The wavenumber vector, through which were computed the values of the left-hand side of dispersion equations ranged from zero to , where: (15) where is maximum angular frequency within the computed frequency interval.The incremental step of wavenumber interval is equal to: (16) the larger incremental step is suitable in case of lower values of frequency .thickness product (approximately up to 3 MHz .mm), which delivers sufficient resolution for finding the roots.For higher values of frequency .thickness product, it is necessary to increase the resolution of the wavenumber vector for proper root localization.The roots of the dispersion equation for given frequency are located at the zero crossings when the left-hand side changes its functional value from positive to negative sign.-see Fig 2. Resulting dispersion curves in terms of phase velocity as a function of frequency and wavenumber as a function of frequency respectively are displayed in Fig. 3.It has to be noted, that the dispersion curves were constructed for aluminium plate of c L = 6300 m/sec, c T = 3100 m/sec and the thickness of the plate was equal to 8 mm.The frequency interval has been divided into 200 steps with a solving time equal to 5.5 seconds (CPU i7 3.4 GHz, 32 GB RAM).For proper estimation of the cut-off frequencies, the number of frequency interval division has to be refined.

Numerical solution for imaginary wavenumbers
The pure imaginary wavenumbers, which represent evanescent waves, are supposed to be in following form: (17) Eqs. 13 and 14 can be therefore rewritten as follows: With Matlab´s incorporated Meshgrid function, which is also supported by, above mentioned, vectorization property, it is possible to relatively quickly and efficiently obtain the zero crossings for the real and imaginary parts of left-hand side of the dispersion relations for symmetric and antisymmetric modes.Figure 5 displays us the evaluated real part of left-hand side of dispersion equation for symmetric modes for particular value of frequency .thickness product.The resulting real and imaginary parts of complex wavenumbers are given by the intersections between the solution of the real and imaginary parts of left-hand side of the dispersion relation for particular mode and frequency .thickness product (See Fig. 6).As mentioned before, the solution of the dispersion equation will give an infinite set of complex wavenumbers.In our case, we will present the results for aluminium plate of c L = 6300 m/sec, c T = 3100 m/sec, thickness of 8 mm and the frequency interval was ⟩ , divided into 100 steps with a solving time equal to 9.2 seconds (CPU i7 3.4 GHz, 32 GB RAM).The ranges for real and imaginary parts of complex wavenumbers were within .The incremental step of wavenumber vector, for both the real and imaginary part, was equal to max(k im/real ) . 10 -3 m -1 .Fig. 7 displays the real and the corresponding imaginary parts of complex wavenumbers for above-mentioned input parameters.

Conclusion
Dispersion curves are an indispensable tool in case of wave scattering problems [1,9,10] or problematics of wave propagation.The main aim of this paper was to present an overall overview of the procedures of finding the real, imaginary and complex wavenumbers for given material parameters and frequency .thickness interval.The procedures of finding the individual types of wavenumbers have been written in Matlab software with incorporation of functions, which support the vectorised operations.Thanks to this fact, there was no extra need to put emphasis on the proper selection of the wavenumber intervals in order to decrease the number of steps, related to the interval range, with the aim of reducing the computational time.The entire set consisting from all three types of wavenumbers can be for given material parameters and frequency range with reasonable division obtained within 20 seconds.Currently, authors are working on development of more time-efficient algorithm in order to apply the results with cooperation of Lamb wave structural health monitoring for various experiments [11][12][13][14][15][16] and also for further application in case of composite structures [17].This work was supported by Specific Research (SP2017/106 and SP2017/136) and by The Ministry of Education, Youth and Sports from the National Programme of Sustainability II (LQ1602).The support is acknowledged.

Fig. 2 .
Fig. 2. Functional values of the left-hand side of dispersion relation for given wavenumber interval (left graph) and the detail of roots (right graph) for aluminium plate of c L = 6300 m/sec, c T = 3100 m/sec and frequency .thickness = 4 MHz .mm

Fig. 3 .
Fig. 3. Dispersion curves in terms of c phase as a function of frequency (left) and wavenumber as a function of frequency (right).The symmetric Lamb modes are represented by red dots, the antisymmetric modes by blue dots

Fig. 5 .
Fig. 5. Real part of left-hand side of dispersion equation (symmetric mode) in case of aluminium plate of c L = 6300 m/sec, c T = 3100 m/sec, thickness of 8 mm and frequency equal to 0.2 MHz

Fig. 6 .
Fig. 6.Real and imaginary parts of complex wavenumbers in case of aluminium plate of c L = 6300 m/sec, c T = 3100 m/sec, thickness of 8 mm and frequency equal to 0.2 MHz

Fig. 7 .
Fig. 7. Real (left) and imaginary (right) parts of complex wavenumbers in case of aluminium plate of c L = 6300 m/sec, c T = 3100 m/sec, thickness of 8 mm and frequency interval of ⟩ .Capital letters denotes corresponding real and imaginary parts