Cauchy Six-Dimensional Formalism for Lamb Waves in Multilayered Plates

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


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 where u is the displacement field and   and   are velocities of the longitudinal and transverse bulk waves, respectively: In (2),  and  are Lamé constants and  is the material density.Then, the displacement field was represented in terms of scalar (Φ) and vector (Ψ) potentials u = ∇Φ + rot Ψ. ( The potentials were assumed to be harmonic in time: Φ (x, ) = Φ  (x)   , Ψ(x, ) = Ψ  (x)   .
Substituting representation (4) into (1) yields two independent Helmholtz equations: To define the spatial periodicity and to simplify the analysis, the splitting spatial argument is needed: where n is the unit wave vector, ^is the unit normal to the median plane of the plate, and w = n × ^.
Remark 1.For the considered waves, it was further assumed that the displacement field does not depend upon argument x ⋅ w.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 (Ψ ⋅ w)w).
The further assumption relates to the periodicity of the potentials in the direction of propagation where the dimensionless complex coordinates   and   are In (8),  = √ −1 and  is the wave number related to the wavelength  by Substituting representations (7) into (5) results in the decoupled ordinary differential equations where the phase speed  relates to the frequency and the wave number by the following relation: The general solution of (10) can be written in the form where The unknown coefficients in (12) are defined (up to a multiplier) from the following boundary conditions on the free surfaces: where 2ℎ is the depth of the plate.Substitution representation (3) into boundary conditions ( 14) yields boundary conditions written in terms of potentials   and   : Substituting solutions (12) into (15), in view of Remark 1, yields the desired dispersion equation, found in [1,3]: 2 ) 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 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 Analysis of ( 16) at ℎ → 0 (the long wave limit) yields the following equation [5]: from where the values for two limiting velocities can be obtained: 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][8][9][10].The analysis of the dispersion curves at different Poisson's ratios (including negative values) was done in [11][12][13].The points of intersection of the dispersion curves were studied in [14].[15] for description of the wave package propagation in hydrodynamics and later extrapolated to acoustic waves in the theory of elasticity; see [16][17][18].Formally, the group velocity can be defined by

Group Velocity. Consider now the notion of the group velocity introduced by Stokes
Numerical studies [19][20][21][22][23][24][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 from where, in view of (22), the secular equation for the group speed takes the form In (24),  is considered as a function of .Theoretical studies of the phase, group, and ray velocities were done in [26].

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][29][30][31][32].With necessary modification, the approach [27] was used for analyzing Lamb waves in anisotropic plates [33][34][35][36][37][38][39][40][41][42][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: where   are arbitrary complex coefficients determined up to a multiplier by satisfying conditions at the plate boundaries, u  is the displacement field of the th partial wave, m  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
where C is the fourth-order elasticity tensor assumed to be positive definite, yields the Christoffel equation where I is the unit diagonal matrix.Equation ( 27) admits the equivalent form 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 m  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 (^⋅x) < 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 m  was considered in [44] with application to Rayleigh waves and in [45] for obtaining the dispersion equation of subsonic Lamb waves.
To obtain the dispersion equation, consider the tractionfree boundary conditions where 2ℎ denotes the overall thickness of the plate.Substituting representation (25) into the boundary conditions ( 29) yields the dispersion equation in the form Finally, (30) can be rewritten in a form better suited for numerical computations: 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.

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][50][51][52][53][54][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 with the unknown amplitude f regarded as a function of the imaginary coordinate   defined by (8).Substituting representation (32) into equation of motion ( 26) yields where The matrix A 1 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: Up to the harmonic multiplier  (n⋅x−) , 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): 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 A −1 1 yields the following expression for the derivative    f(  ): Combining now ( 33)-( 37) produces the desired equation of motion written in terms of vectors f and t ] : where N is the fundamental matrix [58] N = ( ) .
Now, the general solution of ( 38) can be represented in the form where ⃗  6 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 Excluding vector ⃗  6 from (41a) and substituting the resultant expression into (41b) gives the following dispersion equation for an anisotropic plate: Despite the obvious simplicity of deriving (42), the relevant works on Lamb waves in anisotropic plates [56,57] use a more complicated procedure.

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 (41a), (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][61][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][66][67][68][69][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][66][67][68][69][70] will be worked out.

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][73][74][75][76][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: where ⃗   ∈ C 6 ,  = 1, . . ., , are six-dimensional generally complex vectors containing the unknown coefficients in representation (25); D ±  are 6×6-matrices defining displacements by (25) and surface tractions by (30) on both surfaces of a layer: Despite the reported numerical stability [72][73][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.

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.

Basic Notations.
For the considered formalism, the representation for harmonic Lamb wave is searched in the form (32).By introducing a new vector function the equation of motion (26) can be reduced to the Cauchy normal form where and matrices A 1 , A 2 , and A 3 are defined by (34).By the analogy with Stroh formalism, 6 × 6-matrix G will be called the fundamental matrix.Similar to (40), the general solution of ( 46) can be represented in Eulerian exponential form where ⃗  6 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 where Note that matrix Z is nondegenerate due to positive definiteness of the elasticity tensor.
Equation ( 49) allows us to exclude ⃗  6 by expressing the displacements and surface tractions on one side of the plate through the corresponding vectors on the other side; that yields where Matrix T 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 T is independent of boundary conditions.

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.

Traction-Free Plate.
If both sides of a plate are traction free, then Condition (53) means that a 3 × 3-mapping ) 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): Matrices B 1 and B 2 are introduced for consistency with other types of boundary conditions that will be considered later in this section.
The condition ( 54) is equivalent to 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][60][61]); this yields the following equation for  2,lim : Remark 4. Substituting matrices ( 35) into ( 57) yields an isotropic plate (20).

Clamped Plate.
If both sides of a plate are clamped, then 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 B 1 and B 2 : The equation for the limiting speed  2,lim remains the same, as (57).

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 For boundary conditions (60), matrices B 1 and B 2 become Similar to expressions (60), boundary conditions (ii) can be written in the form and B 1 , B 2 become

Dispersion Relations for a Multilayered Plate.
Consider an -layered plate with the homogeneous anisotropic layers.At the interfaces, ideal mechanical contact is assumed.

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 where the lower index indicates the number of a layer and the transfer matrices T  are defined by (52).Now, multiplying both sides of (64) by matrices B 1 , B 2 , similar to (54)-( 56) for the traction-free plate, we arrive at the dispersion equation for the multilayered traction-free plate: Similar to ( 54), ( 65) reflects degeneracy of the mapping The dispersion equations for other types of boundary conditions have the same form as (65) but with obvious replacement of matrices B 1 and B 2 .
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.

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: Within the considered Cauchy formalism, matrices D ±  change: where index  indicates the corresponding layer.

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.

The Overflow Problem
. This problem can be relatively easily solved by a suitable normalization.Let the fundamental matrices G , =1,..., be decomposed into their Jordan normal forms where W  is a matrix composed of eigenvectors and possibly the generalized eigenvectors, while J  is a matrix containing eigenvalues and possibly Jordan blocks if G  is not a semisimple matrix.
With decomposition (69), the exponential matrices in (52) can be represented in the following form [79,80]: Let ( ()  ) =1,...,6 denote (not necessary aliquant) eigenvalues of G  , and 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.

Loss of Precision.
As reported in [78], the extension of the -matrix approach [65] to anisotropic media [67][68][69][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: 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: where |()| < 10 −15 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 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.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:

Numerical Examples
In (75),  1122 and  1212 would correspond to Lamé parameters  and , and  1111 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 For the considered cubic crystal, the measuring parameter is relatively small:  ≈ 0.09.The second cubic crystal has the following elastic parameters: 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]: For the auxetic crystal (78),  ≈ 0.09, yielding the same value as for the first crystal.For these crystals, vectors ^and n have the following coordinates (relative to the principle elasticity directions): ^= (1; 0; 0) , n = (0; cos ; sin ) .
Variation of the limiting velocity  2,lim with respect to angle  is shown in Figure 1.
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.

Multilayered Plates.
Herein, two multilayered tractionfree plates with alternating anisotropic layers are considered; the depth of each layer is ℎ = 14 nm (14 × 10 −9 m).In the subsequent analysis, the direction of propagation of Lamb waves is n = (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 3 and the following elastic constants (in GPa) [82]: For the considered multilayered plate, the dispersion curves for Lamb waves are plotted in Figure 2.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 <  Si 3 N 4   .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   /

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.

81 ) 6 .
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 3 N 4 (-phase) with  = 2329 kg/m 3 and elastic constants (in GPa)[83]  1111 = 514,  1122 = 182,  1133 = 182  2222 = 514,  2233 = 182,  3333 = 514  1212 =  1313 =  2323 = 327.(Remark 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 n are SiC :   n = 12498   n in-plane = 7822   n out of plane = 7049 Si 3 N 4 :   n = 14856   n = 11849.
3.2.1.HomogeneousPlates with Cubic Symmetry.Three plates made of materials with cubic symmetry are considered, all Figure 1: Variation of limiting velocity  2,lim versus angle .