Cauchy Six-Dimensional Formalism for Lamb Waves in Multilayered Plates

Abstract

Propagation of Lamb waves in multilayered elastic anisotropic plates is studied in the framework of combination of the six-dimensional Cauchy formalism and the transfer matrix method. The closed form secular equations for dispersion curves of Lamb waves propagating in multilayered plates with arbitrary elastic anisotropy are obtained.

1. Introduction

Herein, a brief introduction to the theory of Lamb waves and a review of some of the most important works on this matter are presented.

1.1. Lamb Waves in a Homogeneous Isotropic Plate

The first works [1, 2] on waves propagating in an infinite isotropic homogeneous plate with the traction-free boundary surfaces were done at the assumption that the wavelength is much longer than the plate thickness.

The complete theory of harmonic Lamb waves free from the long wavelength limit assumption was presented in [3]. The starting point of the Lamb theory is considering the equation of motion in the form 𝑐_𝑃∇div𝐮−𝑐_𝑆𝜕rotrot𝐮=²𝐮𝜕𝑡², (1) where 𝐮 is the displacement field and 𝑐_𝑃 and 𝑐_𝑆 are velocities of the longitudinal and transverse bulk waves, respectively: 𝑐_𝑃=𝜆+2𝜇𝜌,𝑐_𝑆=𝜇𝜌. (2)

In (2), 𝜆 and 𝜇 are Lamé constants and 𝜌 is the material density. Then, the displacement field was represented in terms of scalar (Φ) and vector (Ψ) potentials 𝐮=∇Φ+rotΨ. (3)

The potentials were assumed to be harmonic in time: Φ(𝐱,𝑡)=Φ^′(𝐱)𝑒^𝑖𝜔𝑡,Ψ(𝐱,𝑡)=Ψ^′(𝐱)𝑒^𝑖𝜔𝑡. (4)

Substituting representation (4) into (1) yields two independent Helmholtz equations: 𝜔Δ+²𝑐²_𝑃Φ^′𝜔=0,Δ+²𝑐²_𝑆Ψ^′=0. (5)

To define the spatial periodicity and to simplify the analysis, the splitting spatial argument is needed: 𝐱=(𝐱⋅𝐧)𝐧+(𝐱⋅𝝂)𝝂+(𝐱⋅𝐰)𝐰, (6) where 𝐧 is the unit wave vector, 𝝂 is the unit normal to the median plane of the plate, and 𝐰=𝐧×𝝂.

Remark 1. For the considered waves, it was further assumed that the displacement field does not depend upon argument 𝐱⋅𝐰. That allowed Lamb [3] to consider scalar potentials Ψ and Ψ^′ in (4) instead of vector ones (actually, Lamb considered the vector potential composed of one nonvanishing component (Ψ⋅𝐰)𝐰).

The further assumption relates to the periodicity of the potentials in the direction of propagation Φ^′(𝐱)=𝜑𝑥^′′𝑒^𝑥′,Ψ^′(𝐱)=𝜓𝑥^′′𝑒^𝑥′, (7) where the dimensionless complex coordinates 𝑥^′ and 𝑥^′′ are 𝑥^′=𝑖𝑟𝐱⋅𝐧,𝑥^′′=𝑖𝑟𝐱⋅𝝂. (8)

In (8), 𝑖=√−1 and 𝑟 is the wave number related to the wavelength 𝑙 by 𝑟=2𝜋𝑙. (9)

Substituting representations (7) into (5) results in the decoupled ordinary differential equations 𝑑²𝜑𝑑𝑥^′′2+𝑐1−²𝑐²_𝑃𝑑𝜑=0,²𝜓𝑑𝑥^′′2+𝑐1−²𝑐²_𝑆𝜓=0, (10) where the phase speed 𝑐 relates to the frequency and the wave number by the following relation: 𝜔𝑐=𝑟. (11) The general solution of (10) can be written in the form 𝜑𝑥^′′=𝐶₁sinh𝛾₁𝑥^′′+𝐶₂cosh𝛾₁𝑥^′′,𝜓𝑥^′′=𝐶₃sinh𝛾₂𝑥^′′+𝐶₄cosh𝛾₂𝑥^′′, (12) where 𝛾₁=𝑐1−²𝑐²_𝑃^1/2,𝛾₂=𝑐1−²𝑐²_𝑆^1/2. (13)

The unknown coefficients in (12) are defined (up to a multiplier) from the following boundary conditions on the free surfaces: 𝐭_𝜈≡𝜆tr(∇𝐮)𝐈+𝜇∇𝐮+∇𝐮^𝑡⋅𝝂=0,𝐱⋅𝝂=±ℎ, (14) where 2ℎ is the depth of the plate. Substitution representation (3) into boundary conditions (14) yields boundary conditions written in terms of potentials 𝜑^′′ and 𝜓^′′: 𝜆ΔΦ^′𝐈+2𝜇∇∇Φ^′+12∇rotΨ^′+∇rotΨ^′𝑡⋅𝝂=0,𝐱⋅𝝂=±ℎ. (15) Substituting solutions (12) into (15), in view of Remark 1, yields the desired dispersion equation, found in [1, 3]: tanh𝛾₂𝑟ℎtanh𝛾₁−𝑟ℎ4𝛾₁𝛾₂1+𝛾²₂^2±1=0. (16)

The sign “+” refers to symmetrical and “−” to antisymmetrical modes. In view of expressions (11) and (13), the obtained dispersion equation implicitly determines the phase velocity 𝑐 as a function of frequency. Equations for velocities related to the long wave and short wave limits were found in [3].

Taking the short wave limit at 𝑟ℎ→∞ in (16) yields 1=4𝛾₁𝛾₂1+𝛾²₂². (17) The latter expression coincides with the secular equation for the speed of Rayleigh waves derived in [4], and hence the first limiting speed coincides with the speed of Rayleigh wave 𝑐_1,lim=𝑐_𝑅. (18) Analysis of (16) at 𝑟ℎ→0 (the long wave limit) yields the following equation [5]: 𝛾₂𝛾₁=4𝛾₁𝛾₂1+𝛾²₂², (19) from where the values for two limiting velocities can be obtained: 𝑐_2,lim=2𝑐_𝑆𝑐1−^2𝑆𝑐²_𝑝.𝑐_3,lim=𝑐_𝑆. (20)(21)

Expression (20) coincides with the limiting velocity value found in [1, 2] and differs from the anticipated value for the long wave limiting velocity of the sound waves in rods.

For the antisymmetric fundamental mode, (16) was analytically solved in [6] using the perturbation technique. Earlier studies of the dispersion curves for Lamb waves in a layer contacting with a half-space were mainly associated with the geophysical applications [7–10]. The analysis of the dispersion curves at different Poisson’s ratios (including negative values) was done in [11–13]. The points of intersection of the dispersion curves were studied in [14].

1.2. Group Velocity

Consider now the notion of the group velocity introduced by Stokes [15] for description of the wave package propagation in hydrodynamics and later extrapolated to acoustic waves in the theory of elasticity; see [16–18]. Formally, the group velocity can be defined by 𝑐^group=𝑑𝜔.𝑑𝑟 (22)

Numerical studies [19–25] of the group velocity dispersion, mainly at Poisson’s condition 𝜆=𝜇, confirmed Rayleigh’s anticipation [16, 17] that the negative values of the group velocity can appear at the very small wave numbers. Later on, numerical computations [25] revealed the existence of a broader range of negative group velocities at Poisson’s ratios belonging to the interval 0.31<𝜈<0.45.

Several additional equations for computing dispersion of the group velocity can be constructed from (16). For example, substituting the phase speed defined by (11) into (13), denoting the left-hand side of (16) by 𝐹^±(𝑟,𝜔), and assuming that 𝜔 is a function of 𝑟, the derivative of (16) with respect to 𝑟 takes the form 𝜕𝐹^±(𝑟,𝜔)+𝜕𝑟𝜕𝐹^±(𝑟,𝜔)𝜕𝜔𝑑𝜔𝑑𝑟=0, (23) from where, in view of (22), the secular equation for the group speed takes the form 𝑐^group=−𝜕𝐹^±(𝑟,𝜔)/𝜕𝑟𝜕𝐹^±().𝑟,𝜔/𝜕𝜔 (24) In (24), 𝑟 is considered as a function of 𝜔. Theoretical studies of the phase, group, and ray velocities were done in [26].

1.3. Homogeneous Anisotropic Plates: Three-Dimensional Formalism

In the first works on Lamb waves propagating in anisotropic plates, a three-dimensional formalism was used. Initially, that formalism was developed in [27] for analysis of Rayleigh waves in an anisotropic half-space; later on, this method was applied to half-spaces with different groups of elastic symmetry [28–32]. With necessary modification, the approach [27] was used for analyzing Lamb waves in anisotropic plates [33–43].

All these publications, except [33] where a more complicated case of a cylindrically anisotropic plate was considered, actually exploited the following representation for the displacement field: 𝐮(𝐱,𝑡)=⁶_𝑘=1𝐶_𝑘𝐦_𝑘𝑒^{𝛾𝑘𝑥′′}_{𝐮𝑘𝑥′′}𝑒^{𝑖(𝑟𝐧⋅𝐱−𝜔𝑡)}, (25) where 𝐶_𝑘 are arbitrary complex coefficients determined up to a multiplier by satisfying conditions at the plate boundaries, 𝐮_𝑘 is the displacement field of the 𝑘th partial wave, 𝐦_𝑘 is a vectorial, generally, complex amplitude, determined by the Christoffel equation (this equation will be introduced later on), and 𝛾_𝑘 is a root of the Christoffel equation. Note that, according to (8), the coordinate 𝑥^′′ is imaginary. Six partial waves in (25) correspond to six (not necessary aliquant) roots of the Christoffel equation.

Substituting representation (25) into equation of motion 𝐀𝜕_𝑥,𝜕_𝑡𝐮≡div_𝑥𝐂⋅⋅∇_𝑥𝐮−𝜌𝐮̈=0, (26) where 𝐂 is the fourth-order elasticity tensor assumed to be positive definite, yields the Christoffel equation 𝛾_𝑘𝝂+𝐧⋅𝐂⋅𝐧+𝛾_𝑘𝝂−𝜌𝑐²𝐈⋅𝐦_𝑘=0, (27) where 𝐈 is the unit diagonal matrix. Equation (27) admits the equivalent form det𝛾_𝑘𝝂+𝐧⋅𝐂⋅𝐧+𝛾_𝑘𝝂−𝜌𝑐²𝐈=0. (28) The left-hand side of (28) represents a polynomial of degree six with respect to the Christoffel parameter 𝛾_𝑘. Equations (27), (28) show that roots 𝛾_𝑘 and the corresponding eigenvectors 𝐦_𝑘 can be considered as functions of the phase speed 𝑐.

Remark 2. (a) For Rayleigh waves, the roots 𝛾_𝑘 in representation (25) should be complex with Im(𝛾_𝑘)<0; this ensures attenuation of Rayleigh wave in the “lower” half-space (𝝂⋅𝐱)<0. The condition Im(𝛾_𝑘)<0 confines the admissible speed interval and decreases the number of summation terms in (25). If Re(𝛾_𝑘)=0 for all partial waves composing Rayleigh wave, then such a wave is called the genuine Rayleigh wave; if Re(𝛾_𝑘)≠0 for some 𝑘, then such a wave is called the generalized Rayleigh wave [28]. For Lamb waves, the cases Re(𝛾_𝑘)=0 and Re(𝛾_𝑘)≠0 are usually not distinguished.
(b) Within the discussed formalism, the case of appearing multiple roots 𝛾_𝑘 and the coincident kernel eigenvectors 𝐦_𝑘 was considered in [44] with application to Rayleigh waves and in [45] for obtaining the dispersion equation of subsonic Lamb waves.
(c) For anisotropic plates, the limiting velocities (20), (21) were computed in [46, 47].

To obtain the dispersion equation, consider the traction-free boundary conditions 𝐭_{𝜈𝐱⋅𝝂=±ℎ}≡±𝝂⋅𝐂⋅⋅∇_𝐱𝐮_{𝐱⋅𝝂=±ℎ}=0, (29) where 2ℎ denotes the overall thickness of the plate. Substituting representation (25) into the boundary conditions (29) yields the dispersion equation in the form ⁶_𝑘=1𝐶_𝑘𝛾_𝑘𝝂⋅𝐂⋅𝝂+𝝂⋅𝐂⋅𝐧⋅𝐦_𝑘𝑒^{±𝑖𝑟𝛾𝑘ℎ}=0. (30)

Finally, (30) can be rewritten in a form better suited for numerical computations:𝛾det¹𝝂⋅𝐂⋅𝝂+𝝂⋅𝐂⋅𝐧⋅𝐦₁𝑒^{+𝑖𝑟𝛾1ℎ}⋯𝛾₆𝝂⋅𝐂⋅𝝂+𝝂⋅𝐂⋅𝐧⋅𝐦₆𝑒^{+𝑖𝑟𝛾6ℎ}𝛾₁𝝂⋅𝐂⋅𝝂+𝝂⋅𝐂⋅𝐧⋅𝐦₁𝑒^{−𝑖𝑟𝛾1ℎ}⋯𝛾₆𝝂⋅𝐂⋅𝝂+𝝂⋅𝐂⋅𝐧⋅𝐦₆𝑒^{−𝑖𝑟𝛾6ℎ}=0. (31)The dispersion equation (31) defines the phase speed as a function of the wave number or, in view of (11), as a function of the circular frequency.

1.4. Homogeneous Anisotropic Plates: Stroh Six-Dimensional Formalism

Initially, Stroh formalism [48] was applied to analysis of Rayleigh waves propagating on a free boundary of an anisotropic half-space [49–55]. The case of the nonsemisimple degeneracy of the fundamental matrix was considered in [53] (that case is associated with appearing multiple roots and the coincident kernel eigenvectors in the Christoffel equation). In [56, 57], the Stroh formalism was applied to the description of Lamb waves propagating in the homogeneous anisotropic plates.

Following [49], the displacement field for Rayleigh or Lamb wave is searched in the form 𝐮(𝐱,𝑡)=𝐟𝑥^′′𝑒^{𝑖(𝑟𝐧⋅𝐱−𝜔𝑡)} (32) with the unknown amplitude 𝐟 regarded as a function of the imaginary coordinate 𝑥^′′ defined by (8). Substituting representation (32) into equation of motion (26) yields 𝐀₁𝜕²_𝑥′′+𝐀₂+𝐀^𝑡₂𝜕_𝑥′′+𝐀₃⋅𝐟𝑥^′′=0, (33) where 𝐀₁=𝝂⋅𝐂⋅𝝂,𝐀₂=𝝂⋅𝐂⋅𝐧,𝐀₃=𝐧⋅𝐂⋅𝐧−𝜌𝑐²𝐈. (34) The matrix 𝐀₁ is nondegenerate due to the assumption of positive definiteness of the elasticity tensor.

Remark 3. For an isotropic material, matrices (34) take the following form: 𝐀₁=𝜆+2𝜇000𝜇000𝜇,𝐀₂=,𝐀0𝜆0𝜇00000₃=𝜇−𝜌𝑐²000𝜆+2𝜇−𝜌𝑐²000𝜇−𝜌𝑐². (35)

Up to the harmonic multiplier 𝑒^{𝑖(𝑟𝐧⋅𝐱−𝜔𝑡)}, the surface tractions on the planes parallel to the free boundary (or the median surface of a plate) can be written in terms of matrices (34): 𝐭_𝜈𝑥^′′≡𝐀₁𝜕_𝑥′′+𝐀₂⋅𝐟𝑥^′′. (36)

The main idea of Stroh formalism is in rewriting the equation of motion (26) in terms of the displacements and surface tractions. Multiplying both sides of (36) by matrix 𝐀₁⁻¹ yields the following expression for the derivative 𝜕_𝑥′′𝐟(𝑥^′′): 𝜕_𝑥′′𝐟𝑥^′′=𝐀₁⁻¹⋅𝐭_𝜈𝑥^′′−𝐀₁⁻¹⋅𝐀₂⋅𝐟𝑥^′′. (37)

Combining now (33)–(37) produces the desired equation of motion written in terms of vectors 𝐟 and 𝐭_𝜈: 𝜕_𝑥′′𝐟𝐭_𝜈𝐟𝐭=𝐍⋅_𝜈, (38) where 𝐍 is the fundamental matrix [58] 𝐍=−𝐀₁⁻¹⋅𝐀₂𝐀₁⁻¹𝐀^𝑡₂⋅𝐀₁⁻¹⋅𝐀₂−𝐀₃−𝐀^𝑡₂⋅𝐀₁⁻¹. (39) Now, the general solution of (38) can be represented in the form 𝐟𝐭_𝜈=exp𝑥^′′𝐍⋅𝐶₆, (40) where 𝐶₆ is the six-dimensional generally complex vector of the unknown coefficients; that vector can be defined (up to a multiplier) by the boundary conditions. Substituting solution (40) into boundary conditions (29) yields𝐟()𝟎()⋅𝐶+𝑖𝑟ℎ=exp+𝑖𝑟ℎ𝐍₆𝐟()𝟎()⋅𝐶,−𝑖𝑟ℎ=exp−𝑖𝑟ℎ𝐍₆. (41a)(41b)Excluding vector 𝐶₆ from (41a) and substituting the resultant expression into (41b) gives the following dispersion equation for an anisotropic plate: ()()⋅𝐈𝟎det𝟎,𝐈⋅exp−2𝑖𝑟ℎ𝐍=0. (42) Despite the obvious simplicity of deriving (42), the relevant works on Lamb waves in anisotropic plates [56, 57] use a more complicated procedure.

1.5. Multilayered Plates: The Transfer Matrix Method

The transfer matrix method suggested for analyzing propagation of Lamb and Rayleigh waves multilayered isotropic media was introduced in [59, 60]. Up to now, the transfer matrix method was combined with the three-dimensional formalism. The main idea of the method is to exclude the unknown coefficients 𝐶_𝑘 in representation (25) expressing them in terms of the only partly known surface tractions and displacements on one of the outer surfaces of the plate; such a procedure is similar to (40a), (41b). If the plate contains several layers, the interface conditions can also be expressed in terms of these surface tractions and displacements via specially constructed transfer matrices. Finally, the boundary conditions on the other outer surface are expressed in terms of the surface tractions and displacements from the first outer surface.

Since its introduction, the transfer matrix method was applied to finding the dispersion relations for Lamb waves propagating in both isotropic [60–62] and monoclinic (or with higher symmetries) [63, 64] multilayered plates. Despite the abundance and simplicity, the original variant of the transfer matrix method revealed some numerical instability, especially when high frequencies and large depths of the layers were considered. To overcome this problem, the numerically stable algorithms were suggested [65–70]. Differing in details, these algorithms have the common principle idea of eliminating the terms containing large exponentials. Such a procedure, developed in [65], is known as the 𝛿-matrix method [66, 67]. The 𝛿-matrix method can also be applied for analyzing dispersion of Lamb waves in a single-layered plate; see [65].

In the following sections, the transfer matrix method will be coupled with the Cauchy six-dimensional formalism, and a numerically stable approach resembling [65–70] will be worked out.

1.6. Multilayered Plates: The Global Matrix Method

The global matrix method introduced in [71] appeared more numerically stable than the original version of the transfer matrix method; see [62, 65] for discussions. In [72–77] applications of the global matrix method to analyses of Lamb wave dispersion are presented.

There are several variants of constructing the global matrix. For a traction-free plate containing 𝑛-homogeneous anisotropic layers, the matrix equation related to the global matrix method can be written in the following form: 𝟎_6𝑛=(𝟎,𝐈)⋅𝐃⁺₁𝐃⁻₁−𝐃⁺₂𝐃⁻₂−𝐃⁺₃()⋯⋯𝟎,𝐈⋅𝐃⁻_𝑛⋅𝐶₁𝐶₂𝐶₃⋮𝐶_𝑛, (43) where 𝐶_𝑘∈ℂ⁶,𝑘=1,…,𝑛, are six-dimensional generally complex vectors containing the unknown coefficients in representation (25); 𝐃^±_𝑘 are 6×6-matrices defining displacements by (25) and surface tractions by (30) on both surfaces of a layer: 𝐃^±_𝑘=𝐦_𝑘1𝑒^{±𝑖𝑟𝛾𝑘1ℎ𝑘}⋯𝐦_𝑘6𝑒^{±𝑖𝑟𝛾𝑘6ℎ𝑘}𝛾_𝑘1𝝂⋅𝐂_𝑘⋅𝝂+𝝂⋅𝐂_𝑘⋅𝐧⋅𝐦_𝑘1𝑒^{±𝑖𝑟𝛾𝑘1ℎ𝑘}⋯𝛾_𝑘6𝝂⋅𝐂_𝑘⋅𝝂+𝝂⋅𝐂_𝑘⋅𝐧⋅𝐦_𝑘6𝑒^{±𝑖𝑟𝛾𝑘6ℎ𝑘}. (44)

Despite the reported numerical stability [72–74], the comparative study of the original transfer matrix method [58, 59] and the global matrix method [71] revealed that actually both methods exhibit some intrinsic instability resulting in loss of accuracy at high frequencies and large depths of the layers [78]. To overcome the persistent instability, in [78] a more refined rearrangement of the exponential terms, appearing in the global matrix method, was suggested.

2. Cauchy Six-Dimensional Formalism

In the following, another variant of a six-dimensional formalism for analyzing dispersion of Lamb waves in plates with arbitrary anisotropy is presented. Since the developed formalism is based on reduction of the second-order equation (33) to the Cauchy normal form, it will be called as the Cauchy formalism.

2.1. Basic Notations

For the considered formalism, the representation for harmonic Lamb wave is searched in the form (32). By introducing a new vector function 𝐯𝑥^′′=𝜕_𝑥′′𝐟𝑥^′′, (45) the equation of motion (26) can be reduced to the Cauchy normal form 𝜕_𝑥′′𝐟𝐯𝐟𝐯=𝐆⋅, (46) where 𝐆=0𝐈−𝐀₁⋅𝐀₃−𝐀₁⋅𝐀₂+𝐀^𝑡₂ (47) and matrices 𝐀₁, 𝐀₂, and 𝐀₃ are defined by (34). By the analogy with Stroh formalism, 6×6-matrix 𝐆 will be called the fundamental matrix.

Similar to (40), the general solution of (46) can be represented in Eulerian exponential form 𝐟𝑥^′′𝐯𝑥^′′=exp𝑥^′′𝐆⋅𝐶₆, (48) where 𝐶₆ is the six-dimensional generally complex vector of the unknown coefficients defined (up to a multiplier) by the boundary conditions.

The surface traction vector is defined by (36). Now, combining (48) and (36), the displacements and surface tractions on both sides of the plate take the form 𝐟()𝐭±𝑖𝑟ℎ_𝜈()()⋅𝐶±𝑖𝑟ℎ=𝐙⋅exp±𝑖𝑟ℎ𝐆₆, (49) where 𝐀𝐙=𝐈𝟎₂𝐀₁. (50) Note that matrix 𝐙 is nondegenerate due to positive definiteness of the elasticity tensor.

Equation (49) allows us to exclude 𝐶₆ by expressing the displacements and surface tractions on one side of the plate through the corresponding vectors on the other side; that yields 𝐟()𝐭−𝑖𝑟ℎ_𝜈()𝐟()−𝑖𝑟ℎ𝐭=𝐓⋅+𝑖𝑟ℎ_𝜈(+𝑖𝑟ℎ), (51) where 𝐓=𝐙⋅exp(−2𝑖𝑟ℎ𝐆)⋅𝐙⁻¹. (52) Matrix 𝐓 can be considered as the transfer matrix, as it “transfers” displacements and surface tractions from one side of the plate to another. It should be noted that matrix 𝐓 is independent of boundary conditions.

2.2. Dispersion Relations for a Homogeneous Plate

Equation (51) is a source for constructing the dispersion relations and obtaining an expression for the limiting velocity 𝑐_2,lim at 𝑟→0.

2.2.1. Traction-Free Plate

If both sides of a plate are traction free, then 𝐭_𝜈(±𝑖𝑟ℎ)=0. (53) Condition (53) means that a 3×3-mapping 𝑅^{3𝐁1,𝐁2⋅𝐓⋅}𝐁₂𝐁₁−−−−−−−−−−−−−−−−→𝑅³ (54) from the three-dimensional space of nontrivial displacements and vanishing surface-tractions on the “upper” surface to the three-dimensional space of the surface-tractions on the “bottom” surface is degenerate. In (54): 𝐁₁=𝟎,𝐁₂=𝐈. (55)

Matrices 𝐁₁ and 𝐁₂ are introduced for consistency with other types of boundary conditions that will be considered later in this section.

The condition (54) is equivalent to det𝐁₁,𝐁₂𝐁⋅𝐓⋅²𝐁₁=0. (56) Thus, (56) is the desired dispersion equation for a plate with free boundaries.

Taking in (54) limit at 𝑟→0 does not lead to a meaningful condition, because of vanishing determinant in (56) at 𝑟=0 irrespective of the phase speed. To derive a meaningful equation, consider the first derivative of (54) with respect to the wave number (see [59–61]); this yields the following equation for 𝑐_2,lim: det𝐁₁,𝐁₂⋅𝐙⋅𝐆⋅𝐙⁻¹⋅𝐁₂𝐁₁=0. (57)

Remark 4. Substituting matrices (35) into (57) yields an isotropic plate (20).

2.2.2. Clamped Plate

If both sides of a plate are clamped, then 𝐟(±𝑖𝑟ℎ)=0. (58) Condition (58) means that mapping from the space of vanishing displacements and nontrivial surface tractions on the “upper” surface to the three-dimensional space of the displacements on the “bottom” surface is degenerate. But, this condition can be expressed by the same equations as (56) with the only difference in matrices 𝐁₁ and 𝐁₂: 𝐁₁=𝐈,𝐁₂=𝟎. (59) The equation for the limiting speed 𝑐_2,lim remains the same, as (57).

2.2.3. Plates with Mixed Boundary Conditions

Herein, two types of mixed boundary conditions at a boundary surface are considered: (i) vanishing normal surface-traction and vanishing tangential displacements and (ii) vanishing tangential surface-tractions and vanishing normal displacements. Boundary conditions (i) can be written in the form 𝝂⋅𝐭_𝜈(±𝑖𝑟ℎ)=0,(𝐈−𝝂⊗𝝂)⋅𝐟(±𝑖𝑟ℎ)=0. (60) For boundary conditions (60), matrices 𝐁₁ and 𝐁₂ become 𝐁₁=𝐈−𝝂⊗𝝂,𝐁₂=𝝂⋅𝐭_𝜈. (61) Similar to expressions (60), boundary conditions (ii) can be written in the form (𝐈−𝝂⊗𝝂)⋅𝐭_𝜈(±𝑖𝑟ℎ)=0,𝝂⋅𝐟(±𝑖𝑟ℎ)=0, (62) and 𝐁₁, 𝐁₂ become 𝐁₁=𝝂⋅𝐭_𝜈,𝐁₂=𝐈−𝝂⊗𝝂. (63)

2.3. Dispersion Relations for a Multilayered Plate

Consider an 𝑛-layered plate with the homogeneous anisotropic layers. At the interfaces, ideal mechanical contact is assumed.

2.3.1. The Transfer Matrix Method

Assuming that propagation of the Lamb wave in each layer is determined by (45)–(50) and writing a sequence of (51) that transfers boundary conditions from top to bottom yield 𝐟−𝑖𝑟ℎ_𝑛𝐭_𝜈−𝑖𝑟ℎ_𝑛=^𝑛_𝑘=1𝐓_𝑘⋅𝐟+𝑖𝑟ℎ₁𝐭_𝜈+𝑖𝑟ℎ₁, (64) where the lower index indicates the number of a layer and the transfer matrices 𝐓_𝑘 are defined by (52).

Now, multiplying both sides of (64) by matrices 𝐁₁, 𝐁₂, similar to (54)–(56) for the traction-free plate, we arrive at the dispersion equation for the multilayered traction-free plate: det𝐁₁,𝐁₂⋅^𝑛_𝑘=1𝐓_𝑘⋅𝐁₂𝐁₁=0. (65) Similar to (54), (65) reflects degeneracy of the mapping 𝑅^{3𝐁1,𝐁2⋅𝑛∏𝑘=1𝐓𝑘⋅}𝐁₂𝐁₁−−−−−−−−−−−−−−−−−−−−−→𝑅³. (66) The dispersion equations for other types of boundary conditions have the same form as (65) but with obvious replacement of matrices 𝐁₁ and 𝐁₂.

As it was with (57), taking limit at 𝑟→0 of (65) does not lead to a meaningful condition, because of vanishing the determinant at 𝑟=0 irrespective of the phase speed. To derive a meaningful equation, consider the first derivative of the mapping (66) with respect to the wave number.

2.3.2. The Global Matrix Method

While the global matrix equation (43) mainly preserves its form, a minor modification is needed to account for other types of boundary conditions at the outer boundaries: 𝟎_6𝑛=𝐁₁,𝐁₂⋅𝐃⁺₁𝐃⁻₁−𝐃⁺₂𝐃⁻₂−𝐃⁺₃𝐁⋯⋯₁,𝐁₂⋅𝐃⁻_𝑛⋅𝐶₁𝐶₂𝐶₃⋮𝐶_𝑛. (67) Within the considered Cauchy formalism, matrices 𝐃^±_𝑘 change: 𝐃^±_𝑘=𝐙_𝑘⋅exp±𝑖𝑟ℎ_𝑘𝐆_𝑘, (68) where index 𝑘 indicates the corresponding layer.

3. Numerical Examples

3.1. Improving Numerical Stability of the Transfer Matrix Method

Performing numerical computations with the exponential matrices leads to numerical instability associated with either loss of precision or overflow due to the presence of the exponential terms with large positive powers.

3.1.1. The Overflow Problem

This problem can be relatively easily solved by a suitable normalization. Let the fundamental matrices 𝐆_{𝑘,𝑘=1,…,𝑛} be decomposed into their Jordan normal forms 𝐆_𝑘=𝐖_𝑘⁻¹⋅𝐉_𝑘⋅𝐖_𝑘, (69) where 𝐖_𝑘 is a matrix composed of eigenvectors and possibly the generalized eigenvectors, while 𝐉_𝑘 is a matrix containing eigenvalues and possibly Jordan blocks if 𝐆_𝑘 is not a semisimple matrix.

With decomposition (69), the exponential matrices in (52) can be represented in the following form [79, 80]: exp−2𝑖𝑟ℎ_𝑘𝐆_𝑘=𝐖_𝑘⁻¹⋅exp−2𝑖𝑟ℎ_𝑘𝐉_𝑘⋅𝐖_𝑘. (70)

Let (𝜆_𝑘^(𝑚))_𝑚=1,…,6 denote (not necessary aliquant) eigenvalues of 𝐆_𝑘, and 𝑎=max_𝑘,𝑚Re−2𝑖𝑟ℎ_𝑘𝜆_𝑘^(𝑚). (71) Multiplying all the exponent matrices in (65) by 𝑒^−𝑎 yields the desired normalization. Indeed, after multiplication, all the exponential terms in (70) will have powers with real parts not exceeding unity.

3.1.2. Loss of Precision

As reported in [78], the extension of the 𝛿-matrix approach [65] to anisotropic media [67–70] and the global matrix methods still manifest numerical instability, which is mainly due to loss of numerical precision. To overcome this problem, the following measures are suggested: (i) computing minimal eigenvalue of the resultant matrix (for both transfer and global matrix methods), instead of computing the corresponding determinants, and (ii) using longer mantissas instead of “double precision” computations.

While being more time consuming, the first measure allows one to avoid possible overlook of a root. For example, consider a typical case with the following diagonal matrix: 𝐃=diag(𝜀,𝜀,𝛼), (72) where 𝜀 is a small (generally complex) eigenvalue associated with the root and 𝛼 is a large eigenvalue. In the vicinity of the root the determinant of matrix (72) has no change of sign and, thus, can lead to overlooking of the root. Moreover, if 𝛼 increases at approaching the root, the determinant may even increase, at least being not very close to the root.

Necessity to perform multiprecision computations can be demonstrated by considering the following matrix: (𝑡)(𝑡)𝐃=1+𝜀−1−11−𝜀, (73) where |𝜀(𝑡)|<10⁻¹⁵ is a small function of the argument 𝑡. It is obvious that computing the determinant with double precision accuracy gives zero, inspite of the applied numerical algorithm. However, analytical evaluation yields det𝐃=𝜀²(𝑡), (74) confirming necessity in multiprecision computations.

In numerical examples considered below mantissas from ~25 decimal digits (quadruple precision) to ~1000 decimal digits were used, depending upon the analyzed problem.

3.2. Numerical Examples

3.2.1. Homogeneous Plates with Cubic Symmetry

Three plates made of materials with cubic symmetry are considered, all with unit material density. To estimate deviation of elastic moduli of a cubic crystal from the isotropic elastic constants, the following measure can be introduced: 𝐶𝛿≡¹¹¹¹𝐶₁₁₂₂+2𝐶₁₂₁₂−1. (75)

In (75), 𝐶₁₁₂₂ and 𝐶₁₂₁₂ would correspond to Lamé parameters 𝜆 and 𝜇, and 𝐶₁₁₁₁ would be equal to the denominator in (75) if the material was isotropic.

The first considered crystal has elastic parameters that only slightly differ from the isotropic ones 𝐶₁₁₁₁=1,𝐶₁₁₂₂=0.1,𝐶₁₂₁₂=0.5. (76) For the considered cubic crystal, the measuring parameter is relatively small: 𝛿≈0.09.

The second cubic crystal has the following elastic parameters: 𝐶₁₁₁₁=1,𝐶₁₁₂₂=0,𝐶₁₂₁₂=0.25. (77) These parameters differ from the isotropic case in a much greater extent with 𝛿=1.

The third considered crystal belongs to a set of auxetic materials characterized by negative Poisson’s ratios [81]: 𝐶₁₁₁₁=1,𝐶₁₁₂₂=−0.1,𝐶₁₂₁₂=0.5. (78) For the auxetic crystal (78), 𝛿≈0.09, yielding the same value as for the first crystal.

For these crystals, vectors 𝝂 and 𝐧 have the following coordinates (relative to the principle elasticity directions): 𝝂=(1;0;0),𝐧=0;cos𝜓;sin𝜓. (79) Variation of the limiting velocity 𝑐_2,lim with respect to angle 𝜓 is shown in Figure 1.

Figure 1 

Variation of limiting velocity 𝑐2,lim versus angle 𝜓.

The plots in Figure 1 reveal that the limiting velocity 𝑐_2,lim is fairly informative for identification of the propagation direction of the limiting Lamb wave: even for the first and third crystals, there are evident variations of 𝑐_2,lim values with respect to 𝜓.

Remark 5. In the NDT, the limiting Lamb wave propagating with the velocity 𝑐_2,lim plays a very important role, since such a wave can easily be excited. Indeed, energy needed for excitation is proportional to the square of the frequency, and for the considered limiting wave, the corresponding frequency vanishes.

3.2.2. Multilayered Plates

Herein, two multilayered traction-free plates with alternating anisotropic layers are considered; the depth of each layer is ℎ=14 nm (14×10⁻⁹ m). In the subsequent analysis, the direction of propagation of Lamb waves is 𝐧=(1,0,0) in coordinate axes oriented along principal directions of the elasticity tensors introduced below.

The alternating layers are made of hexagonal crystal SiC with 𝜌=3220 kg/m³ and the following elastic constants (in GPa) [82]: 𝐶₁₁₁₁=503,𝐶₁₁₂₂=109,𝐶₁₁₃₃𝐶=52₂₂₂₂=503,𝐶₂₂₃₃=52,𝐶₃₃₃₃𝐶=552₂₃₂₃=𝐶₁₃₁₃=160,𝐶₁₂₁₂=𝐶₁₁₁₁−𝐶₁₁₂₂2; (80) the unit normals to planes of elastic isotropy are coincident with vector 𝝂 (normal to the median plane of a layer). The second alternating layer is cubic crystal Si₃N₄ (𝛾-phase) with 𝜌=2329 kg/m³ and elastic constants (in GPa) [83] 𝐶₁₁₁₁=514,𝐶₁₁₂₂=182,𝐶₁₁₃₃𝐶=182₂₂₂₂=514,𝐶₂₂₃₃=182,𝐶₃₃₃₃𝐶=514₁₂₁₂=𝐶₁₃₁₃=𝐶₂₃₂₃=327. (81)

Remark 6. It should be noted that the SiC crystal is piezoelectric, so in a more rigorous approach [84, 85] the change of elastic constants reflecting electrically open (constant electric displacement) or closed (constant electric field) circuit conditions should be foreseen.

The computed bulk wave velocities (in m/s) for the plane waves with normal 𝐧 are SiC∶𝑐_𝑃𝐧=12498𝑐_𝑆𝐧^in-plane𝑐=7822_{𝑆𝐧outofplane}=7049Si₃N₄∶𝑐_𝑃𝐧=14856𝑐_𝑆𝐧=11849. (82) For the considered multilayered plate, the dispersion curves for Lamb waves are plotted in Figure 2.

(a)

(b)

Figure 2 

Dispersion curves for 3-layerd (a) and 11-layered (b) traction-free plates with alternating SiC-Si3N4 layers.

These plots were constructed by the described transfer matrix method in combination with the quadruple precision computations (mantissas with ~34 decimal digits).

One of the most interesting observations concerns almost vertical parts of the dispersion curves at low-frequency limits; these are clearly seen in plots for both of the plates. The corresponding relative limiting velocity increases with increase of number of layers; however, the relative limiting velocity does not exceed unity, meaning that 𝑐_2,lim<𝑐^Si3N4_𝑃𝑛. Another observation concerns the presence of the vertical asymptotes at large 𝜔 for lower branches. A more detailed analysis reveals that in both of these plates these asymptotes correspond to the same relative phase velocity 𝑐^SiC_𝑅𝑛/𝑐^Si3N4_𝑃𝑛≈0.497, where 𝑐^SiC_𝑅𝑛 denotes velocity of Rayleigh wave in SiC for the considered propagation direction.

4. Conclusion

The Cauchy six-dimensional formalism is developed for analyzing propagation of Lamb waves in anisotropic multilayered plates. The presented numerical examples reveal its ability to obtain dispersion curves for plates with large number of anisotropic layers.