Floquet-Bloch Theory and Its Application to the Dispersion Curves of Nonperiodic Layered Systems

Abstract

Dispersion curves play a relevant role in nondestructive testing. They provide estimations of the elastic and geometrical parameters from experiments and offer a better perspective to explain the wave field behavior inside bodies. They are obtained by different methods. The Floquet-Bloch theory is presented as an alternative to them. The method is explained in an intuitive manner; it is compared to other frequently employed techniques, like searching root based algorithms or the multichannel analysis of surface waves methodology, and finally applied to fit the results of a real experiment. The Floquet-Bloch strategy computes the solution on a unit cell, whose influence is studied here. It is implemented in commercially finite element software and increasing the number of layers of the system does not bring additional numerical difficulties. The lateral unboundedness of the layers is implicitly taken care of, without having to resort to artificial extensions of the modelling domain designed to produce damping as happens with perfectly matched layers or absorbing regions. The study is performed for the single layer case and the results indicate that for unit cell aspect ratios under 0.2 accurate dispersion curves are obtained. The method is finally used to estimate the elastic parameters of a real steel slab.

Pablo Gomez Garcia 1, 2 and Jose-Paulino Fernandez-Álvarez 1

Academic Editor:Xiao-Qiao He

1, Hydrogeophysics and NDT Modelling Unit, University of Oviedo, C/Gonzalo Gutierrez Quiros s/n, 33600 Mieres, Spain
2, Dynamics Division, Applied Mechanics Department, Chalmers University of Technology, Hörsalsvagen 7, 41296 Gothenburg, Sweden

Received 20 October 2014; Revised 28 November 2014; Accepted 29 November 2014; 12 January 2015

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

1. Introduction

Floquet-Bloch (hereafter F-B) theory provides a strategy to analyze the behavior of systems with a periodic structure. Floquet's seminal paper dealt with the solution of 1D partial differential equations with periodic coefficients [1]. In solid state physics, Bloch generalized Floquet's results to 3D systems and obtained the description of the wave function associated with an electron traveling across a periodic crystal lattice [2]. This wave function is a solution of the Schrödinger equation with a periodic potential and Bloch showed that it was the product of a simple plane wave multiplied by a periodic function with the same periodicity of the lattice. The mathematical description of these ideas, in the context of quantum mechanics, can be found in [3, 4].

In the literature dealing with wave propagation problems in mechanical systems the theory is referred to as Floquet-Bloch theory or, simply, Floquet theory. In layered systems, due to the heterogeneity of the relevant elastic properties, to particular geometric features, or to both, only certain wave modes can physically propagate inside the structure [5]. Each of these modes can be identified by a determined-generally nonlinear-function relating the time frequency and the spatial frequency (or wave number). These relationships are called dispersion curves default and, as they summarize all the oscillatory behavior of the system, their calculation is of paramount importance in NDE applications [6].

Vibrations occur also in objects with periodic structure [7]. These problems usually admit a separation between the time and the spatial dependent parts of the solution. For instance, the Helmholtz equation is a known example of equation describing the spatial behavior [8]. There, the physical periodic structure of the studied object translates into spatial periodicity of its coefficients. Therefore, the F-B theory has been applied to obtain the dispersive properties of different mechanical periodic systems [8-12].

Many relevant structures can be assumed to be layered systems of infinite extent, for example, [13, 14] in civil engineering constructions, [15, 16] in optics, or [17] in electromagnetics. Therefore, theoretical methods and experimental techniques to obtain their dispersion curves have been devised. From the theoretical side, different matrix techniques have been developed to address the calculation. They involve numerical computational methods whose complexity increases with the number of layers in the system [18, 19] or more recently [6].

In laboratory experiments or field work, the dispersion curves can be obtained using, for example, the multichannel analysis of surface waves (MASW) method. The MASW procedure involves collecting equally spaced measures of vibration along a profile on the system surface using, for example, accelerometers. The resulting 2D space-time discrete image is Fourier-transformed to the frequency-wave number [figure omitted; refer to PDF] domain and then processed to build the dispersion curves [20, 21]. The method has some drawbacks inherent to the Fourier transform limitations which will be discussed later. The MASW has been applied successfully in the characterization of pavement systems [13], as a seismic data acquisition technique [20] or for geotechnical characterization [22]. The MASW strategy is here also used to perform a computer numerical simulation of the system, closely mimicking the field setup. The issue of infinite lateral extent is usually tackled by using perfectly matched layers (PML) [23-26] as has been done here or absorbing regions. Both techniques present drawbacks [26].

In this paper, an alternative way to calculate the dispersion curves of layered systems with infinite lateral extent using the F-B theory is presented. The method has never been applied to the dispersion curves calculations of nonperiodic layered systems. Here it is used to obtain the dispersion curves of a single layer case and to estimate the elastic parameters of a real steel slab, for showing the method. However, the novelty in this work is that it can be applied to an arbitrary number of layers, even if the layers are anisotropic or orthotropic, with the same complexity level. The power of the method is that the equations are solved by the finite element software, because the F-B theory only affects the propagation term, which is the same, whatever the nature of the layers. It is not necessary to develop the equations for each specific problem and to generate complex codes to get the dispersion relations.

The F-B theory reduces the problem to calculations performed in the so-called unit cell, subject to certain specific boundary conditions derived from the F-B theory and elastodynamics. The influence of the size of the unit cell is ascertained. The results are first compared with the dispersion curves derived from the Rayleigh equations [5], solved by a searching root numerical method. Comparison is also made with the curves resulting from a FEM computer simulation followed by a 2D Fourier transformation. Finally, a real experiment was performed on a steel slab employing the MASW method and the empirical dispersion curves were compared with the analytical and numerical ones.

The results show that the F-B method compares favorably with other methods and fits accurately the empirical data, providing a good alternative to obtain dispersion curves in layered systems. The F-B technique can be run on a finite element package like COMSOL Multiphysics, can be applied to an arbitrary number of layers in the system, with the same complexity level, and eliminates issues of infinite lateral extent.

2. Floquet-Bloch Theory: Explanation

The Floquet-Bloch theory provides a strategy to obtain a set of solutions of a linear ordinary equations system of the form [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the solution vector and the matrix [figure omitted; refer to PDF] is periodic such that [figure omitted; refer to PDF] for a certain period [figure omitted; refer to PDF] . At first sight it might seem that the solution of such a problem would have to be also [figure omitted; refer to PDF] -periodic. But Floquet showed that this need not be so. There exists, however, a simple relationship between the solution's behaviors inside one period and outside it. If [figure omitted; refer to PDF] is a fundamental matrix of solutions, then another matrix [figure omitted; refer to PDF] can be found such that [figure omitted; refer to PDF] [figure omitted; refer to PDF] can be constructed by setting [figure omitted; refer to PDF] in (2) such that [figure omitted; refer to PDF] . A simplest case is obtained using [figure omitted; refer to PDF] so that [figure omitted; refer to PDF] . As there is not a unique choice for the fundamental matrix [figure omitted; refer to PDF] and how it is exactly chosen depends on the problem, [figure omitted; refer to PDF] is also not unique. But its eigenvalues are intrinsic of the problem and, under the right transformation, can be used as a propagator or evolution factor [figure omitted; refer to PDF] relating the value of the solution at a point inside the period with its value at a point outside of it. Only the solution inside a period is, therefore, needed verifying that [figure omitted; refer to PDF] Following the classical nomenclature [figure omitted; refer to PDF] is known as Floquet multiplier, [figure omitted; refer to PDF] being the complex Floquet exponent. Moreover, Floquet found that the solution at any point can also be factored in two terms: [figure omitted; refer to PDF] Here [figure omitted; refer to PDF] is a periodic function, playing the role of the eigenvectors if [figure omitted; refer to PDF] was a constant matrix and carrying the periodicity [figure omitted; refer to PDF] of the coefficients of the problem. The complex exponential distorts the strict periodicity of [figure omitted; refer to PDF] incorporating damping or ever growing effects in the amplitudes of the solution depending on the value of [figure omitted; refer to PDF] . This is why solutions are, in general, not periodic and also why the Floquet perspective is usually employed to study their stability. A solution will be stable if the Floquet multipliers verify [figure omitted; refer to PDF] [27].

3. Guided Waves in Layers: Analytical Dispersion Curves

The guided wave propagation problem in a homogeneous, isotropic, and infinite single layer has been widely treated in the literature [5, 28]. In this paper, we follow the theory developed in [5]. So consider an infinite (in [figure omitted; refer to PDF] direction), homogeneous, isotropic, and elastic layer with thickness [figure omitted; refer to PDF] as shown in Figure 1.

An infinite homogenous isotropic and elastic layer with thickness [figure omitted; refer to PDF] as part of an infinite 3D layer under plain strain consideration (a). Typical deformation for the symmetric and antisymmetric modes (b) and the dispersion curves in [figure omitted; refer to PDF] representation, for a given thickness [figure omitted; refer to PDF] , P and S velocities [figure omitted; refer to PDF] and [figure omitted; refer to PDF] , obtained with searching root algorithm (c). The parts marked with dashed boxes and the intersection points between modes present difficulties with typical searching root algorithms.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

The conservation of momentum equation plus free traction boundary conditions leads to a system of equations that produces a solution when its determinant vanishes. In the absence of body forces the equation reads [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the displacement vector, [figure omitted; refer to PDF] is the elastic constants tensor, and [figure omitted; refer to PDF] is the strain tensor. If [figure omitted; refer to PDF] , then [figure omitted; refer to PDF] . Here, [figure omitted; refer to PDF] is a scalar field and [figure omitted; refer to PDF] is a vector field. For the plane strain case, the fields P-SV and SH are decoupled. Considering the P-SV field the [figure omitted; refer to PDF] component of the displacement field and its partial derivative [figure omitted; refer to PDF] vanish. A solution of the form [figure omitted; refer to PDF] appears. It expresses the propagation of a shape [figure omitted; refer to PDF] (Figure 1) trapped in the thickness of the layer along the [figure omitted; refer to PDF] -direction with wave number [figure omitted; refer to PDF] and frequency [figure omitted; refer to PDF] . The dispersion relation for the single infinite layer case is [5] [figure omitted; refer to PDF] with [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the compressional wave (P) velocity and [figure omitted; refer to PDF] the shear wave (S) velocity. Expression equation (6) is known as the Rayleigh-Lamb equation (R-L). The plus sign in the (R-L) expression corresponds with symmetric modes and the minus with the antisymmetric ones (Figure 1).

3.1. Calculation of the Dispersion Curves from the Analytical Rayleigh-Lamb Equation

The R-L equation is transcendental, so no closed analytical solution is available. It can, though, be cast into a form amenable to the use of iterative, root finding, local algorithms, but these present various difficulties arising from the nature of the equations. First, due to the tangent functions, the left hand side is discontinuous at certain points where local algorithms for smooth functions will find difficulties [29]. Due to this and to the sampling rate characteristics of the searching algorithm, some roots might be missed.

A visual inspection to the pattern followed by the dispersion curves on a phase velocity versus frequency [figure omitted; refer to PDF] plot (Figure 1) clearly shows the presence of near vertical and horizontal stretches with physical relevance. For instance, the horizontal line towards which the modes A0 and S0 converge contains the information of the so-called Rayleigh waves [30]. The S1 mode, on the other hand, owns a point with vertical tangent, that is, minimum wave number, in an otherwise near vertical portion of the curve with the most important property of having zero group velocity, producing then a useful resonance [31].

For local root searching methods, like the Newton-Raphson method [32], the strategy consists in using one found root as a seed to calculate the next point. Methods performing a 1D search with a fixed value of [figure omitted; refer to PDF] (resp., [figure omitted; refer to PDF] ) will very poorly characterize a near vertical (resp., horizontal) portion of the curve, unless the sampling rate is extremely dense. Other methods start from the frequencies at zero wave number where the roots can be calculated analytically and try to follow up each curve with some sort of linear predictions [33]. Problems arise at the mode intersections (Figure 1). It has been suggested that the search grid should locally mimic the behavior of the dispersion curves [34]. At any rate, the methods are computationally intensive for useful tolerances and their extension to the case of more layers, where the curve pattern is much more intricate, is inefficient.

In this paper, the root loci of the R-L equation have been obtained using the bisection method, which always converges. It is very slow and fails wherever multiple roots exist in the proposed interval. A pair of time frequency, phase velocity values are input in the R-L equation and its sign evaluated. Keeping the frequency fixed, the velocity has been varied with a step of [figure omitted; refer to PDF] , from 0 to [figure omitted; refer to PDF] m/s. When the sign changes, the bisection method is iteratively applied to obtain a root with an accuracy [figure omitted; refer to PDF] m/s. The process is repeated for different frequency values and the roots are plotted in [figure omitted; refer to PDF] (Figure 1).

The elastic parameters input to the R-L equation are thickness [figure omitted; refer to PDF] m, longitudinal wave (P) velocity [figure omitted; refer to PDF] m/s, and shear wave (S) velocity [figure omitted; refer to PDF] m/s because they will be shown to produce the best fit to the empirical tests performed on a steel slab (Section 6).

Some of the problems discussed above can be clearly seen in the [figure omitted; refer to PDF] representation (Figure 1(c)). The upper part of the modes (S1, S2) was impossible to calculate because, due to their almost infinite slope, the number of roots there may be found to be infinite or zero. Starting, instead, from a fixed velocity shifts the difficulty to the near horizontal segments of the curves. Besides, a regular searching grid in the [figure omitted; refer to PDF] plane becomes an irregular sampling in the [figure omitted; refer to PDF] domain, where some stretches with the same amount of complexity independent of the number of layers may not be sufficiently well characterized. The proposed F-B based method will allow an arbitrary degree of accuracy in the calculation of the curves without suffering from these problematic issues. This might also be an advantage for systematic performance of sensitivity analysis, testing a relevant amount of perturbed models around one found solution.

4. Dispersion Curves Calculation Using FEM with a MASW Type Scheme

The multichannel analysis of surface waves (MASW) methodology is a procedure to numerically or empirically calculate dispersion curves. The strategy is to measure (or obtain numerically), on the surface of the system, the wave field in a number of equally spaced points and take readings at a certain temporal sampling rate [35]. This sampled 2D time-space field is then Fourier-transformed into a 2D time frequency-wave number domain. A continuous surface of amplitudes is now obtained in this [figure omitted; refer to PDF] domain using interpolation or fitting methods. The dispersion curves can be obtained as the loci of local maxima of that surface. The MASW method has been discussed in the literature [20, 36, 37].

Given the inverse relationship between the length of the profile on the surface and the wave number sampling interval and, correspondingly, between the range of sensed wave numbers and the distance between adjacent sensors [38], field constraints (i.e., of logistic or economic type) on the feasible profile length or on the number of sensors available for use do influence the representation. Reciprocity [21] alleviates the difficulty allowing keeping one sensor fixed while moving the impactor in the so-called multichannel record with one receiver (MROR) technique.

When the scheme is applied to perform a numerical calculation on guided waves, the simulation domain has to be finite. This brings the problem of unwanted reflections and mode conversions at those boundaries, coming back into the relevant domain and corrupting the signal. The naive option of extending the domain implies increasing the number of nodes and computation times and assumes that boundary reflection events are separated in time from the studied events, which might not be possible. The use of the so-called absorbing regions, where the wave field enters and is computationally absorbed, has been treated in the literature [23-26]. For instance, perfectly matched layers (PML) or absorbing layers using increasing damping (ALID) have been frequently employed. PML are regions attached to the boundaries where the wave enters and decays exponentially [25, 26]. In ALID, the domain is enlarged with layers of the same material but with increasing damping parameters [39]. Although both are implemented in commercial finite element codes, the success of both techniques relies on iteratively finding the optimum design parameters. This can be time consuming and varies with the characteristics lack of the system one is interested in.

In this case, COMSOL Multiphysics has been used for the simulation. The characteristic parameters employed for the PML are PML scaling factor [figure omitted; refer to PDF] and PML curvature parameter [figure omitted; refer to PDF] (Figure 2).

Figure 2: Typical scheme for MASW method implementation in a FEM software. The computational domain has been set to [figure omitted; refer to PDF] m, the PML parameters used [figure omitted; refer to PDF] and [figure omitted; refer to PDF] , the distance between receivers [figure omitted; refer to PDF] cm, the PMLâ[euro](TM)s sizes [figure omitted; refer to PDF] m, and the thickness [figure omitted; refer to PDF] cm.

[figure omitted; refer to PDF]

4.1. Numerical Implementation

In this section the results of the numerical simulations using the FEM software COMSOL Multiphysics, following a MASW procedure and employing PML to take lateral unboundedness into account, are presented. The dispersion curves will be discussed and serve as a reference to be compared with those obtained by the searching root algorithm (Section 3), those calculated using the F-B approach (Section 5), and the field curves presented in Section 6. The elastic parameters are the same as in Section 3 and the thickness now is [figure omitted; refer to PDF] cm. The length of the simulated profile is [figure omitted; refer to PDF] m and that of the PML is [figure omitted; refer to PDF] m.

A frequency domain study has been performed with a step of [figure omitted; refer to PDF] Hz and sweeping frequencies up to [figure omitted; refer to PDF] kHz with uniform energy distribution. 40 point accelerometers are separated [figure omitted; refer to PDF] cm. The results after 2D Fourier transform and interpolation are shown in Figure 3.

Dispersion curves for a single layer in frequency-wave number representation (a) and in phase velocity-frequency representation (b). Parts that are lacking in the phase velocity plot (b) correspond to wave numbers greater than the measured ones.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

Some observations are in order. First, certain portions of some modes are not excited [40-42]. This affects mainly the lower frequency parts of the zero order symmetric mode (S0). Should a half cycle sinusoidal function have been used as impact model [43], the frequency energy input at lower frequencies would have been relatively higher and the A0 low frequency part would be visible whereas the, now visible, higher frequency branch would be absent. This argument does not affect the S0 mode.

The empty triangular space in the bottom right part of Figure 3 corresponds to pairs of frequency-wave numbers impossible to be reached with a sensor separation of [figure omitted; refer to PDF] cm. For this sampling distance, the highest representable wave number is [figure omitted; refer to PDF] . As [figure omitted; refer to PDF] , the phase velocity stays always over [figure omitted; refer to PDF] for any given frequency.

Additional numerical artifacts arise due to the finite length of the profile. As this is mathematically equivalent to multiplication by a boxcar window it generates, in the wave number domain, a convolution with a sinc function. This effect has been zoomed in Figure 4. This spatial convolution shows up, Figure 4(b), as a spurious repetition of some branches. For more complex patterns present in more than one layer system the real position of the dispersion lines becomes more uncertain.

Dispersion curves in [figure omitted; refer to PDF] representation obtained with padded signals (a). Zoom of the box marked in the left plot and the fast Fourier transform of the box function used for padding the signals (b). The distance between zeros of the obtained sinc function is [figure omitted; refer to PDF] , where [figure omitted; refer to PDF] is the spatial sampled rate and [figure omitted; refer to PDF] is the number of measured points.

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

All the described difficulties with the MASW domain are absent in the Floquet-Bloch technique.

5. Floquet-Bloch Theory and Guided Waves in Layers

The equation of movement for the single layer case was presented in Section 3 (5). To solve it, it is necessary to try plane wave type solutions: [figure omitted; refer to PDF]

Equation (8) establishes that the displacement solution of the problem is a certain shape [figure omitted; refer to PDF] , trapped in the thickness of the system, which propagates in [figure omitted; refer to PDF] -direction with a certain wave number [figure omitted; refer to PDF] and frequency [figure omitted; refer to PDF] along the infinite lateral extension (Figure 1).

The layer is considered infinite; however, the solutions can be computed over a finite computational domain (unit cell), subject to certain boundary conditions. Consider an infinite layer as is shown in Figure 5.

Infinite layer and the unit cell elected (a). Computational domain (unit cell) properties and dimensions (b). F-B boundary conditions applied to the laterals of the computational domain (unit cell) for a certain shape [figure omitted; refer to PDF] .

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

According to [43], a layer behavior will appear experimentally if both dimensions of the layer are, at least, ten times the thickness dimension. This goal will be achieved in the experimental study, taking into account the values presented in Figure 7.

Consider, now, the spatial part of (8): [figure omitted; refer to PDF]

Therefore, the displacement field at the left side of the unit cell (Figure 5(c)) will be related to the displacement at the right side, such that [figure omitted; refer to PDF]

Note that the function [figure omitted; refer to PDF] does not change (Figure 5(c)). It defines the form of the considered mode and is only the propagative term (the exponential ones) which defines the propagation of [figure omitted; refer to PDF] .

Moreover, a lateral infinite layered system is a trivially periodic medium in the propagation [figure omitted; refer to PDF] -direction. Due to it, the F-B theory states that the solution of the problem can be written as [figure omitted; refer to PDF]

For certain periodic function [figure omitted; refer to PDF] , F-B exponent [figure omitted; refer to PDF] and where the [figure omitted; refer to PDF] th-dependency appears due to the relationship between the wave number [figure omitted; refer to PDF] and the F-B exponent [figure omitted; refer to PDF] .

From (9) and (11), the solution of the problem will be equivalent if the wave vector [figure omitted; refer to PDF] and the F-B exponent [figure omitted; refer to PDF] are related as [figure omitted; refer to PDF] where [figure omitted; refer to PDF] is the length of the unit cell in the propagation [figure omitted; refer to PDF] -direction (Figures 5(b) and 5(c)) such that substituting (12) into (9) leads to [figure omitted; refer to PDF] where now the term [figure omitted; refer to PDF] is clearly [figure omitted; refer to PDF] -periodic. However, due to the [figure omitted; refer to PDF] th-dependent term, a fixed value of the F-B exponent [figure omitted; refer to PDF] and [figure omitted; refer to PDF] values of the wave vector [figure omitted; refer to PDF] appear.

Now, the relation of the solutions at the left and right sides of the unit cell can be used to define the proposed F-B boundary conditions as [figure omitted; refer to PDF]

Due to the relationship between [figure omitted; refer to PDF] and the [figure omitted; refer to PDF] th-dependent periodic term, the dispersion relation [figure omitted; refer to PDF] becomes [figure omitted; refer to PDF] which can be calculated over a unit cell of arbitrary length [figure omitted; refer to PDF] solving the following eigenvalue problem (Figure 5): [figure omitted; refer to PDF] where [figure omitted; refer to PDF] are the displacement field and [figure omitted; refer to PDF] any component in the [figure omitted; refer to PDF] -direction of the stress tensor (traction free-surface boundaries in the unit cell).

The theory presented in this section applied to the lateral sides of the unit cell can be used for any kind of layered systems, whatever the nature of the layers (isotropic or anisotropic) and the number of them are. The reason is that the proposed F-B boundary conditions affect only the propagative part of the solution ( [figure omitted; refer to PDF] -direction), which appears in all terms of the equations when the problem is treated in a classical analytical way and is simplified disappearing from the equations.

Therefore, the complicated part of the problem, which is to obtain the function [figure omitted; refer to PDF] (given the propagation modes), is solved by the finite element software. This kind of software allows solving a huge kind of complicated systems which are intractable analytically. This is the case, for example, of layered systems with a big number of layers or when the layers are anisotropic composites.

However, the perspective for getting analytical solutions is always the same, to obtain a function (complicated in general) in the thickness direction (the mode) which propagates laterally. In this approach, the propagation part is always extracted in the equations when infinite layered systems are considered.

Because of the above reason, the theory applied to define the F-B boundary conditions can be used for any complicated systems while they are laterally infinite. The boundary conditions only affect the propagation part and allow converting the infinity analytical problem, in a finite problem, which can be solved with commercial Finite Element software. Therefore, it is not necessary to develop equations (generally complicated) for each certain problem and perform complex numerical codes based on searching roots algorithms [44, 45].

5.1. Aspect Ratio Effects

Using COMSOL Multiphysics software, the eigenvalue problem equation (15) can be solved by sweeping different values of [figure omitted; refer to PDF] and obtaining the corresponding values [figure omitted; refer to PDF] .

However, to transform the problem of the infinite layered system to a finite problem through the F-B boundary conditions, introduce artificial aspects due to the periodicity of the [figure omitted; refer to PDF] -dependent term. For this reason, given certain F-B wave number [figure omitted; refer to PDF] , from (12), [figure omitted; refer to PDF] values of the wave number [figure omitted; refer to PDF] (and eigenfrequencies [figure omitted; refer to PDF] ) will be obtained.

An example is presented in Figure 6(a), where the value of the F-B wave vector is fixed as [figure omitted; refer to PDF] . For this value, all eigenfrequencies [figure omitted; refer to PDF] corresponding with the propagation wave vector values derived from (12) are obtained (in the example only are presented the fundamental one [figure omitted; refer to PDF] and the first one corresponding to [figure omitted; refer to PDF] as [figure omitted; refer to PDF] , which match with the wave vector [figure omitted; refer to PDF] ).

Dispersion curves in [figure omitted; refer to PDF] representation obtained with the proposed Floquet-Bloch method for different aspect ratios of the unit cell. [figure omitted; refer to PDF] (a) and [figure omitted; refer to PDF] (b).

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

Outline of the experiment assembly (a), a photograph (b), and the impacts generated in the experiment in the time domain (c) and spectra (d). The average time has been [figure omitted; refer to PDF] .

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

(d) [figure omitted; refer to PDF]

5.2. Results and Discussion Derived from the Proposed Floquet-Bloch Method

Since the calculation is developed in terms of the F-B wave vector, the solution becomes [figure omitted; refer to PDF] -periodic, and the dispersion representation obtained has all eigenfrequencies inside the first periodic zone of the solution (Figure 6). Due to it, the branches of the dispersion curves are reflected back in the limits of the zone, being necessary to take into account the aspect ratio of the computational domain (unit cell) in terms of getting a good dispersion curves representation.

Different aspect ratios of the unit cell [figure omitted; refer to PDF] have been explored [figure omitted; refer to PDF] to obtain the dispersion curves. Results for the cases [figure omitted; refer to PDF] and [figure omitted; refer to PDF] are shown in Figure 6. The value of [figure omitted; refer to PDF] chosen for the unit cell establishes at [figure omitted; refer to PDF] the periodicity in the wave number domain. The reflected lines intersect the canonical modes in the unit cell preventing them from being clearly identified (Figure 6(a)). Aspect ratios smaller than 0.2 are usually enough to obtain a representative number of modes defining clear dispersion curves (Figure 6(b)).

In the case [figure omitted; refer to PDF] the limit [figure omitted; refer to PDF] is situated at low wave number values and the reflected branches render the spectrum unclear. However, a good representation is obtained for [figure omitted; refer to PDF] (emphasized in the blue frame of Figure 6(b)) because the reflected branches arise at higher wave numbers (Figure 6(b)).

A good criterion is to choose the lateral dimension at least five times the thickness in the computational domain.

6. Real Test on a Steel Slab and Results

A real NDT experiment has been conducted on the surface of a steel slab shown in Figure 7. The profile was set up along the symmetry axis and included 57 equally spaced measurement points separated at a distance of 5 cm. An instrumented hammer and accelerometers have been employed following the MROR method described in Section 4. The accelerometers were threaded producing a coupling resonance around 12 KHz. An outline of the experiment configuration, together with the relevant dimensions and photographs of the test setup and the impacts generated, is shown in Figure 7.

The shape of the hammer impacts in the time domain is very consistent (Figure 7). The average impact duration is [figure omitted; refer to PDF] . Beyond [figure omitted; refer to PDF] , the amplitudes of the impact spectra are very weak and fall mostly under the measurement noise level.

6.1. Results and Discussion

The resulting empirical dispersion curves can be seen in Figure 8(c).

Superposition of the dispersion curves obtained with the numerical searching root method (6) and with the proposed F-B method, for aspect ratio [figure omitted; refer to PDF] (a). Overlap of F-B method and the simulated MASW method dispersion curves (b). Overlap of F-B method and the experimental dispersion curves (c).

(a) [figure omitted; refer to PDF]

(b) [figure omitted; refer to PDF]

(c) [figure omitted; refer to PDF]

The estimated parameters using the proposed F-B boundary conditions are presented in Table 1.

Table 1: Estimated elastic parameters of a real steel slab using the experimental signals obtained with MASW method fitted with the proposed F-B method.

Density (kg/m3 )	P wave velocity (m/s)	S wave velocity (m/s)
7850	5800 ± 30	3200 ± 23

The estimated errors are obtained with the propagation law (16), taking into account the temporal and spatial frequency steps [figure omitted; refer to PDF] Hz, [figure omitted; refer to PDF] based on the Nyquist criteria and on the fact that [figure omitted; refer to PDF] . The computation time in finite element software is under 2 minutes, computing the curves in an i7 PC processor: [figure omitted; refer to PDF]

A fit has been achieved where only the A0 mode is clearly seen in the measurements. There are two main reasons for that. On the one hand, the tip of the hammer only significantly inputs frequencies until 15 kHz (Figure 7(d)) so that every mode, including the A0, above this frequency will not be seen. On the other hand, the S0 mode is not there because its excitability is very low [42]. The excitability is a concept related to what parts of the different modes are detectable at the surface. These commercial receivers measure the out-of-plane (perpendicular to the surface) acceleration component, which has a very low excitability value at the low frequency part of the S0 mode. The vertical band of high energy near 1-1.2 kHz is the coupling frequency of the accelerometer.

6.2. Comparison of the Three Methodologies

This section presents the comparison of the different methods to obtain dispersion curves: the searching root method, the numerical MASW method together with the proposed F-B method, and their ability to match the empirical dispersion curves. The results have been calculated for a layer with the same thickness as the slab used for the experiment. The numerical MASW simulation performed in COMSOL Multiphysics has used as impacts the experimental ones. Figure 8(a) presents the dispersion curves obtained with the searching root algorithm and with F-B method and the numerical MASW simulation together with the relevant part of the experiment.

From Figure 8, the proposed F-B method matches perfectly with the analytical solution of the dispersion relation in the layer and with the experimental and simulated dispersion curves obtained with MASW method. As was emphasized in previous sections, the F-B method provides better results in the vertical zones of the modes than the method based on searching roots algorithm (Figure 8(a)). Another feature is that the F-B method as analytical method is not affected by the excitability concept, as happens with the MASW method because the F-B method is not based on the measurements of the displacement field. The results of the simulated experiment (Figure 8(b)) are again affected by the excitability and the S0 mode is not obtained as happens in the real experiment, which might be taken as a proof that the S0 absence is not due to pitfalls in the measuring process. The experimental results are affected by the resonant frequencies of the accomplished method used for the accelerometers (Figure 8(c)) highlighting that the measuring process could be improved proving different coupling systems.

7. Conclusions

The main conclusion of this study is that the F-B theory can be used to compute the theoretical dispersion curves of layered systems with infinite lateral extension over a finite unit cell. The method is applied directly using a finite element commercial software and is free of the drawbacks associated with other numerical procedures used. It is also a tidy method to calculate curves with more than one layer, avoiding ambiguities at crossing points or with particular slopes. Based on the obtained results, aspect ratios with lower values than [figure omitted; refer to PDF] of the computational domain are enough to obtain a good number of modes in a clear plot representation. The method allows obtaining the associated eigenvectors in (15) so the excitability curves could be obtained too in the same computing process and compared with those results obtained with the MASW procedure. The infinite lateral extension is implicitly taken into account without the need to use ad hoc domain extensions.

Acknowledgments

Thanks are due to the Wave Propagation Division of Chalmers University and specially to Professor Anders Boström because of his suggested corrections in the development of this work.

Conflict of Interests

The authors Pablo Gomez Garcia and Jose-Paulino Fernandez-Álvarez declare that there is no conflict of interests regarding the publication of this paper.

References

[1] G. Floquet, "Sur les equations differentielles lineaires a coefficients periodiques," Annales scientifiques de l'École Normale Superieure , vol. 12, pp. 47-88, 1883.

[2] F. Bloch, "Über die Quantenmechanik der Elektronen in Kristallgittern," Zeitschrift für Physik A , vol. 52, no. 7-8, pp. 555-600, 1929.

[3] C. Kittel, H. Y. Fan, "Introduction to solid state physics," American Journal of Physics , vol. 25, no. 5, pp. 330, 1957.

[4] L. Brillouin Wave Propagation in Periodic Structures: Electric Filters and Crystal Lattices , 1946.

[5] K. F. Graff Wave Motion in Elastic Solids , Dover, 1991.

[6] M. Lowe, "Matrix techniques for modeling ultrasonic waves in multilayered media," IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control , vol. 42, no. 4, pp. 525-542, 1995.

[7] M. Collet, M. Ouisse, M. Ruzzene, M. N. Ichchou, "Floquet-Bloch decomposition for the computation of dispersion of two-dimensional periodic, damped mechanical systems," International Journal of Solids and Structures , vol. 48, no. 20, pp. 2837-2848, 2011.

[8] J. Tausch, "Computing Floquet-Bloch modes in biperiodic slabs with boundary elements," Journal of Computational and Applied Mathematics , vol. 254, pp. 192-203, 2013.

[9] J. Butler, G. Evans, L. Pang, P. Congdon, "Analysis of grating-assisted directional couplers using the Floquet-Bloch theory," Journal of Lightwave Technology , vol. 15, no. 12, pp. 2301-2315, 1997.

[10] A. Safaeinili, D. E. Chimenti, B. A. Auld, S. K. Datta, "Floquet analysis of guided waves propagating in periodically layered composites," Composites Engineering , vol. 5, no. 12, pp. 1471-1476, 1995.

[11] S. Dahmen, M. Ben Amor, M. H. Ben Ghozlen, "An inverse procedure for determination of material constants of a periodic multilayer using Floquet wave homogenization," Composite Structures , vol. 92, no. 2, pp. 430-435, 2010.

[12] C. Potel, J.-F. de Belleval, Y. Gargouri, "Floquet waves and classical plane waves in an anisotropic periodically multilayered medium: application to the validity domain of homogenization," The Journal of the Acoustical Society of America , vol. 97, no. 5, pp. 2815, 1995.

[13] N. Ryden, "Surface wave testing of pavements," The Journal of the Acoustical Society of America , vol. 125, no. 4, pp. 2603, 2009.

[14] N. Ryden, M. J. S. Lowe, "Guided wave propagation in threelayer pavement structures," The Journal of the Acoustical Society of America , vol. 116, no. 5, pp. 2902, 2004.

[15] R. C. Thompson, "Optical waves in layered media," Journal of Modern Optics , vol. 37, no. 1, pp. 147-148, 1990.

[16] P. Yeh, A. Yariv, A. Y. Cho, "Optical surface waves in periodic layered media," Applied Physics Letters , vol. 32, no. 2, pp. 104-105, 1978.

[17] H. Chen, C. T. Chan, "Electromagnetic wave manipulation by layered systems using the transformation media concept," Physical Review B-Condensed Matter and Materials Physics , vol. 78, no. 5, 2008.

[18] W. T. Thomson, "Transmission of elastic waves through a stratified solid medium," Journal of Applied Physics , vol. 21, pp. 89-93, 1950.

[19] N. A. Haskell, "The dispersion of surface waves on multi-layered media," Bulletin of the Seismological Society of America , vol. 43, pp. 17-34, 1953.

[20] R. D. Miller, J. Xia, C. B. Park, J. M. Ivanov, "Multichannel analysis of surface waves to map bedrock," Leading Edge , vol. 18, no. 12, pp. 1392-1396, 1999.

[21] N. Ryden, P. Ulriksen, C. B. Park, R. D. Miller, J. Xia, J. Ivanov, "High frequency MASW for non-destructive testing of pavements-accelerometer approach," in Proceedings of the 14th EEGS Symposium on the Application of Geophysics to Engineering and Environmental Problems, vol. 14, of Roadbed Applications, pp. RBA5, 2001.

[22] D. Penumadu, C. B. Park, "Multichannel Analysis of Surface Wave (MASW) method for geotechnical site characterization," Geotechnical Special Publication , pp. 957-966, 2005.

[23] D. Givoli, "Non-reflecting boundary conditions," Journal of Computational Physics , vol. 94, no. 1, pp. 1-29, 1991.

[24] M. Israeli, S. A. Orszag, "Approximation of radiation boundary conditions," Journal of Computational Physics , vol. 41, no. 1, pp. 115-135, 1981.

[25] J.-P. Berenger, "A perfectly matched layer for the absorption of electromagnetic waves," Journal of Computational Physics , vol. 114, no. 2, pp. 185-200, 1994.

[26] M. B. Drozdz Efficient finite element modelling of ultrasound waves in elastic media [Ph.D. dissertation] , Imperial College, 2008.

[27] M. S. P. Eastham http://books.google.es/books?id=LUHvAAAAMAAJ The Spectral Theory of Periodic Differential Equations , of Texts in Mathematics, Scottish Academic Press, Chatto & Windus, London, UK, 1973.

[28] J. D. Achenbach Wave Propagation in Elastic Solids , vol. 16, 1973.

[29] Z. Su, L. Ye, Y. Lu, "Guided Lamb waves for identification of damage in composite structures: a review," Journal of Sound and Vibration , vol. 295, no. 3-5, pp. 753-780, 2006.

[30] I. A. Viktorov Rayleigh and Lamb Waves: Physical Theory and Applications , Plenum Press, New York, NY, USA, 1967.

[31] A. Gibson, J. S. Popovics, "Lamb wave basis for impact-echo method analysis," Journal of Engineering Mechanics , vol. 131, no. 4, pp. 438-443, 2005.

[32] J. L. Rose Ultrasonic Waves in Solid Media , C. U. Press, 1999.

[33] N. Gandhi, J. E. Michaels, "Efficient perturbation analysis of Lamb wave dispersion curves," Review of Progress in Quantitative Nondestructive Evaluation , vol. 29, of AIP Conference Proceedings, pp. 215-222, 2010.

[34] L. de Marchi, A. Marzani, S. Caporale, N. Speciale, "Ultrasonic guided-waves characterization with warped frequency transforms," IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control , vol. 56, no. 10, pp. 2232-2240, 2009.

[35] D. N. Alleyne, P. Cawley, "A 2-dimensional Fourier transform method for the quantitative measurement of Lamb modes," in Proceedings of the IEEE 1990 Ultrasonics Symposium, vol. 2, pp. 1143-1146, December 1990.

[36] C. B. Park, R. D. Miller, J. Xia, "Multichannel analysis of surface waves," Geophysics , vol. 64, no. 3, pp. 800-808, 1999.

[37] R. D. Miller, J. Xia, J. Ivanov, C. B. Park, "Multichannel analysis of surface waves (MASW), active and passive methods," The Leading Edge , vol. 26, pp. 60-64, 2007.

[38] R. G. Lyons Understanding Digital Signal Processing , 2004.

[39] M. Drozdz, L. Moreau, M. Castaings, M. J. S. Lowe, P. Cawley, "Efficient numerical modelling of absorbing regions for boundaries of guided waves problems," AIP Conference Proceedings , vol. 820, no. 1, pp. 126-133, 2006.

[40] K. Aki, P. G. Richards Quantitative Seismology , University Science Books, 2002.

[41] P. Wilcox, "Modeling the excitation of lamb and SH waves by point and line sources," http://link.aip.org/link/?APC/700/206/1&Agg=doi in AIP Conference Proceedings, vol. 700, pp. 206-213, AIP, 2004.

[42] K. Luangvilai Attenuation of ultrasonic Lamb waves with applications to material characterization and condition monitoring [Ph.D. dissertation] , 2007.

[43] M. Sansalone, "Impact-echo: the complete story," ACI Structural Journal , vol. 94, no. 6, pp. 777-786, 1997.

[44] A. Maghsoodi, A. Ohadi, M. Sadighi, "Calculation of wave dispersion curves in multilayered composite-metal plates," Shock and Vibration , vol. 2014, 2014.

[45] A. Kamal, V. Giurgiutiu, "Stiffness transfer matrix method (STMM) for stable dispersion curves solution in anisotropic composites," in Health Monitoring of Structural and Biological Systems, 906410, vol. 9064, of Proceedings of SPIE, San Diego, Calif, USA, March 2014.