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

Authors: Pablo Gómez García; José-Paulino Fernández-Álvarez
Journal: Mathematical Problems in Engineering
 (2015)

---

## 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.


---

## Body

## 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 calleddispersion 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ω,k 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(1)f1x=Mxfx,
where f(x) is the solution vector and the matrix M is periodic such that M(x+L)=M(x) for a certain period L. At first sight it might seem that the solution of such a problem would have to be also L-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 F is a fundamental matrix of solutions, then another matrix B can be found such that
(2)Fx+L=BFx.B can be constructed by setting x=0 in (2) such that B=F-1(0)F(L). A simplest case is obtained using F(0)=I so that B=F(L). As there is not a unique choice for the fundamental matrix F and how it is exactly chosen depends on the problem, B 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 ρ 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
(3)f^x+L=ρf^x.
Following the classical nomenclature ρ=exp⁡(kFL) is known as Floquet multiplier, kF being the complex Floquet exponent. Moreover, Floquet found that the solution at any point can also be factored in two terms:
(4)f^x=pxexp⁡kFx.
Here p(x) is a periodic function, playing the role of the eigenvectors if M was a constant matrix and carrying the periodicity L of the coefficients of the problem. The complex exponential distorts the strict periodicity of p(x) incorporating damping or ever growing effects in the amplitudes of the solution depending on the value of |ρ|. 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 |ρ|≤1 [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 x direction), homogeneous, isotropic, and elastic layer with thickness 2b as shown in Figure 1.An infinite homogenous isotropic and elastic layer with thicknessd=2b 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 (v-f) representation, for a given thickness b=0.045m, P and S velocities c1=5800m/s and c2=3200m/s, 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)
(b)
(c)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(5)ρu¨(x,y,t)+∇[C(x,y)∇ε(u(x,y,t))]=0,
where u is the displacement vector, C(x,y) is the elastic constants tensor, and ∇ε(u(x,y,t)) is the strain tensor. If u∈C2, then u=∇ϕ+∇×H. Here, ϕ is a scalar field and H is a vector field. For the plane strain case, the fields P-SV and SH are decoupled. Considering the P-SV field the z component of the displacement field and its partial derivative uz=∂uz=0 vanish. A solution of the form u(x,y,t)=f(y)exp⁡(kx-ωt) appears. It expresses the propagation of a shape f(y) (Figure 1) trapped in the thickness of the layer along the x-direction with wave number k and frequency ω. The dispersion relation for the single infinite layer case is [5]
(6)tanαbtanβb=-4αβk2(k2-β2)2±1
with
(7)α2=ω2c12-k2;β2=ω2c22-k2,
where c1 is the compressional wave (P) velocity and c2 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(v-f) 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 k (resp., ω) 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δv=200m/s, from 0 to 104m/s. When the sign changes, the bisection method is iteratively applied to obtain a root with an accuracy δε=0.3m/s. The process is repeated for different frequency values and the roots are plotted in (v-f) (Figure 1).The elastic parameters input to the R-L equation are thicknessd=0.045m, longitudinal wave (P) velocity c1=5800m/s, and shear wave (S) velocity c2=3200m/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(v-f) 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 (f-k) plane becomes an irregular sampling in the (v-f) 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.


## 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(v-f) 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 k (resp., ω) 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δv=200m/s, from 0 to 104m/s. When the sign changes, the bisection method is iteratively applied to obtain a root with an accuracy δε=0.3m/s. The process is repeated for different frequency values and the roots are plotted in (v-f) (Figure 1).The elastic parameters input to the R-L equation are thicknessd=0.045m, longitudinal wave (P) velocity c1=5800m/s, and shear wave (S) velocity c2=3200m/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(v-f) 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 (f-k) plane becomes an irregular sampling in the (v-f) 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 (v-f) 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 factorF=1 and PML curvature parameter n=1 (Figure 2).Figure 2
Typical scheme for MASW method implementation in a FEM software. The computational domain has been set toL=3m, the PML parameters used F=1 and n=1, the distance between receivers δx=5cm, the PML’s sizes δL=1m, and the thickness d=10cm.


### 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 (Section3), 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 d=10 cm. The length of the simulated profile is Lp=2 m and that of the PML is δL=1 m.A frequency domain study has been performed with a step ofδf=100Hz and sweeping frequencies up to f=50kHz with uniform energy distribution. 40 point accelerometers are separated δx=5cm. 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)
(b)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 Figure3 corresponds to pairs of frequency-wave numbers impossible to be reached with a sensor separation of δx=5cm. For this sampling distance, the highest representable wave number is k=62.82m-1. As vph=2πf/k, the phase velocity stays always over vmin⁡=2πf/62.82 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 Figure4. 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(k-f) 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 Δk=2π/Nδx=3.1, where δx is the spatial sampled rate and N is the number of measured points.
(a)
(b)All the described difficulties with the MASW domain are absent in the Floquet-Bloch technique.


## 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 (Section3), 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 d=10 cm. The length of the simulated profile is Lp=2 m and that of the PML is δL=1 m.A frequency domain study has been performed with a step ofδf=100Hz and sweeping frequencies up to f=50kHz with uniform energy distribution. 40 point accelerometers are separated δx=5cm. 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)
(b)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 Figure3 corresponds to pairs of frequency-wave numbers impossible to be reached with a sensor separation of δx=5cm. For this sampling distance, the highest representable wave number is k=62.82m-1. As vph=2πf/k, the phase velocity stays always over vmin⁡=2πf/62.82 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 Figure4. 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(k-f) 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 Δk=2π/Nδx=3.1, where δx is the spatial sampled rate and N is the number of measured points.
(a)
(b)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 Section3  (5). To solve it, it is necessary to try plane wave type solutions:
(8)ux,y,t=fyexp⁡ikx-ωt.Equation (8) establishes that the displacement solution of the problem is a certain shape f(y), trapped in the thickness of the system, which propagates in x-direction with a certain wave number k and frequency ω 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 Figure5.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 shapef(y).
(a)
(b)
(c)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):
(9)ux,y=fyexp⁡ikx.Therefore, the displacement field at the left side of the unit cell (Figure5(c)) will be related to the displacement at the right side, such that
(10)ux+L,y=fyexp⁡ikxexp⁡⁡ikL=u(x,y)exp⁡(ikL).Note that the functionf(y) 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 f(y).Moreover, a lateral infinite layered system is a trivially periodic medium in the propagationx-direction. Due to it, the F-B theory states that the solution of the problem can be written as
(11)unkFBx,y=unpkFBx,yexp⁡ikFBx.For certain periodic functionunpkFB(x,y), F-B exponent kFB and where the nth-dependency appears due to the relationship between the wave number k and the F-B exponent kFB.From (9) and (11), the solution of the problem will be equivalent if the wave vector k and the F-B exponent kFB are related as
(12)k=2nπL+kFB,
where L is the length of the unit cell in the propagation x-direction (Figures 5(b) and 5(c)) such that substituting (12) into (9) leads to
(13)un,kFBx,y=fyexp⁡2nπxLexp⁡kFBx,
where now the term unpkFB=f(y)exp⁡(2nπx/L) is clearly L-periodic. However, due to the nth-dependent term, a fixed value of the F-B exponent kFB and n values of the wave vector k 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(14)unkBx+L,y=unpkFBx,yexp⁡ikFBL.Due to the relationship betweenkFB and the nth-dependent periodic term, the dispersion relation ω(k) becomes ωn(kFB) which can be calculated over a unit cell of arbitrary length L solving the following eigenvalue problem (Figure 5):
(15)ρωn2kFBun,kFBx,y,ωn+∇C(x,y)×∇ε(un,kFB(x,y,ωn))=0∀(x,y)∈ΩR,∇C(x,y)∇ε(un,kFB(x,y,ωn))y=0∀x,y∈δΩRy,un,kBx,y+L,2b,ωn=un,kB((x,y),ωn)eikFBL∀(x,y)∈δΩRx,
where un,kFB are the displacement field and [∇C(x,y)∇ε(un,kFB(x,y,ωn))]y any component in the y-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 (x-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 functionf(y) (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 kFB and obtaining the corresponding values ωn(kFB).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 then-dependent term. For this reason, given certain F-B wave number kFB, from (12), n values of the wave number kn (and eigenfrequencies ωn) will be obtained.An example is presented in Figure6(a), where the value of the F-B wave vector is fixed as kFB=k0. For this value, all eigenfrequencies ωn corresponding with the propagation wave vector values derived from (12) are obtained (in the example only are presented the fundamental one ω0 and the first one corresponding to n=1 as ω1, which match with the wave vector k1=π/L+k0).Dispersion curves in(f-k) representation obtained with the proposed Floquet-Bloch method for different aspect ratios of the unit cell. Ra=1 (a) and Ra=0.1 (b).
(a)
(b)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Δt=143±7μs.
(a)
(b)
(c)
(d)


### 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 becomesL-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 cellRa=2b/L have been explored Ra=1,0.4,0.2,0.1 to obtain the dispersion curves. Results for the cases Ra=1 and Ra=0.1 are shown in Figure 6. The value of L chosen for the unit cell establishes at π/L 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 caseRa=1 the limit π/L is situated at low wave number values and the reflected branches render the spectrum unclear. However, a good representation is obtained for Ra=0.1 (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.


## 5.1. Aspect Ratio Effects

Using COMSOL Multiphysics software, the eigenvalue problem equation (15) can be solved by sweeping different values of kFB and obtaining the corresponding values ωn(kFB).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 then-dependent term. For this reason, given certain F-B wave number kFB, from (12), n values of the wave number kn (and eigenfrequencies ωn) will be obtained.An example is presented in Figure6(a), where the value of the F-B wave vector is fixed as kFB=k0. For this value, all eigenfrequencies ωn corresponding with the propagation wave vector values derived from (12) are obtained (in the example only are presented the fundamental one ω0 and the first one corresponding to n=1 as ω1, which match with the wave vector k1=π/L+k0).Dispersion curves in(f-k) representation obtained with the proposed Floquet-Bloch method for different aspect ratios of the unit cell. Ra=1 (a) and Ra=0.1 (b).
(a)
(b)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Δt=143±7μs.
(a)
(b)
(c)
(d)


## 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 becomesL-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 cellRa=2b/L have been explored Ra=1,0.4,0.2,0.1 to obtain the dispersion curves. Results for the cases Ra=1 and Ra=0.1 are shown in Figure 6. The value of L chosen for the unit cell establishes at π/L 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 caseRa=1 the limit π/L is situated at low wave number values and the reflected branches render the spectrum unclear. However, a good representation is obtained for Ra=0.1 (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 Figure7. 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 (Figure7). The average impact duration is Δt=143±7μs. Beyond f=15kHz, 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 Figure8(c).Superposition of the dispersion curves obtained with the numerical searching root method (6) and with the proposed F-B method, for aspect ratio Ra=0.1 (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)
(b)
(c)The estimated parameters using the proposed F-B boundary conditions are presented in Table1.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 ± 23The estimated errors are obtained with the propagation law (16), taking into account the temporal and spatial frequency steps δf=166.7 Hz, δk=0.1754 based on the Nyquist criteria and on the fact that v=2πf/k. The computation time in finite element software is under 2 minutes, computing the curves in an i7 PC processor:
(16)δv=∂v∂fδf+∂v∂kδk=2πkδf+-2πfk2δk.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 (Figure7(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. Figure8(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 Figure8, 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.


## 6.1. Results and Discussion

The resulting empirical dispersion curves can be seen in Figure8(c).Superposition of the dispersion curves obtained with the numerical searching root method (6) and with the proposed F-B method, for aspect ratio Ra=0.1 (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)
(b)
(c)The estimated parameters using the proposed F-B boundary conditions are presented in Table1.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 ± 23The estimated errors are obtained with the propagation law (16), taking into account the temporal and spatial frequency steps δf=166.7 Hz, δk=0.1754 based on the Nyquist criteria and on the fact that v=2πf/k. The computation time in finite element software is under 2 minutes, computing the curves in an i7 PC processor:
(16)δv=∂v∂fδf+∂v∂kδk=2πkδf+-2πfk2δk.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 (Figure7(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. Figure8(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 Figure8, 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 thanRa≤0.2 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.



---

*Source: 475364-2015-01-12.xml*