# Cauchy Six-Dimensional Formalism for Lamb Waves in Multilayered Plates

Authors: Sergey V. Kuznetsov
Journal: ISRN Mechanical Engineering
 (2013)

---

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


---

## Body

## 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
(1)cP∇divu-cSrotrotu=∂2u∂t2,
where u is the displacement field and cP and cS are velocities of the longitudinal and transverse bulk waves, respectively:
(2)cP=λ+2μρ,cS=μρ.In (2), λ and μ are Lamé constants and ρ is the material density. Then, the displacement field was represented in terms of scalar (Φ) and vector (Ψ) potentials
(3)u=∇Φ+rotΨ.The potentials were assumed to be harmonic in time:(4)Φ(x,t)=Φ′(x)eiωt,Ψ(x,t)=Ψ′(x)eiωt.Substituting representation (4) into (1) yields two independent Helmholtz equations:
(5)(Δ+ω2cP2)Φ′=0,(Δ+ω2cS2)Ψ′=0.To define the spatial periodicity and to simplify the analysis, the splitting spatial argument is needed:(6)x=(x·n)n+(x·ν)ν+(x·w)w,
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 argumentx·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(7)Φ′(x)=φ(x′′)ex′,Ψ′(x)=ψ(x′′)ex′,
where the dimensionless complex coordinates x′ and x′′ are
(8)x′=irx·n,x′′=irx·ν.In (8), i=-1 and r is the wave number related to the wavelength l by
(9)r=2πl.Substituting representations (7) into (5) results in the decoupled ordinary differential equations
(10)d2φdx′′2+(1-c2cP2)φ=0,d2ψdx′′2+(1-c2cS2)ψ=0,
where the phase speed c relates to the frequency and the wave number by the following relation:
(11)c=ωr.
The general solution of (10) can be written in the form
(12)φ(x′′)=C1sinh(γ1x′′)+C2cosh(γ1x′′),ψ(x′′)=C3sinh(γ2x′′)+C4cosh(γ2x′′),
where
(13)γ1=(1-c2cP2)1/2,γ2=(1-c2cS2)1/2.The unknown coefficients in (12) are defined (up to a multiplier) from the following boundary conditions on the free surfaces:
(14)tν≡(λtr(∇u)I+μ(∇u+∇ut))·ν=0,x·ν=±h,
where 2h is the depth of the plate. Substitution representation (3) into boundary conditions (14) yields boundary conditions written in terms of potentials φ′′ and ψ′′:
(15)(λΔΦ′I+2μ(∇∇Φ′+12(+(∇rotΨ′)t∇rotΨ′+(∇rotΨ′)t)(∇∇Φ′+12)(∇∇Φ′+12)·ν=0,x·ν=±h.
Substituting solutions (12) into (15), in view of Remark 1, yields the desired dispersion equation, found in [1, 3]:
(16)tanh(γ2rh)tanh(γ1rh)-(4γ1γ2(1+γ22)2)±1=0.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 c 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 atrh→∞ in (16) yields
(17)1=(4γ1γ2(1+γ22)2).
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
(18)c1,lim=cR.
Analysis of (16) at rh→0 (the long wave limit) yields the following equation [5]:
(19)γ2γ1=(4γ1γ2(1+γ22)2),
from where the values for two limiting velocities can be obtained:
(20)c2,lim=2cS1-cS2cp2.(21)c3,lim=cS.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
(22)cgroup=dωdr.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 F±(r,ω), and assuming that ω is a function of r, the derivative of (16) with respect to r takes the form
(23)∂F±(r,ω)∂r+∂F±(r,ω)∂ωdωdr=0,
from where, in view of (22), the secular equation for the group speed takes the form
(24)cgroup=-(∂F±(r,ω))/∂r(∂F±(r,ω))/∂ω.
In (24), r 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:
(25)u(x,t)=(∑k=16Ckmkeγkx′′︸uk(x′′))ei(rn·x-ωt),
where Ck are arbitrary complex coefficients determined up to a multiplier by satisfying conditions at the plate boundaries, uk is the displacement field of the kth partial wave, mk is a vectorial, generally, complex amplitude, determined by the Christoffel equation (this equation will be introduced later on), and γk is a root of the Christoffel equation. Note that, according to (8), the coordinate x′′ 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
(26)A(∂x,∂t)u≡divxC··∇xu-ρu¨=0,
where C is the fourth-order elasticity tensor assumed to be positive definite, yields the Christoffel equation
(27)[(γkν+n)·C·(n+γkν)-ρc2I]·mk=0,
where I is the unit diagonal matrix. Equation (27) admits the equivalent form
(28)det[(γkν+n)·C·(n+γkν)-ρc2I]=0.
The left-hand side of (28) represents a polynomial of degree six with respect to the Christoffel parameter γk. Equations (27), (28) show that roots γk and the corresponding eigenvectors mk can be considered as functions of the phase speed c.Remark 2.
(a) For Rayleigh waves, the rootsγk in representation (25) should be complex with Im(γk)<0; this ensures attenuation of Rayleigh wave in the “lower” half-space (ν·x)<0. The condition Im(γk)<0 confines the admissible speed interval and decreases the number of summation terms in (25). If Re(γk)=0 for all partial waves composing Rayleigh wave, then such a wave is called the genuine Rayleigh wave; if Re(γk)≠0 for some k, then such a wave is called the generalized Rayleigh wave [28]. For Lamb waves, the cases Re(γk)=0 and Re(γk)≠0 are usually not distinguished.
(b) Within the discussed formalism, the case of appearing multiple rootsγk and the coincident kernel eigenvectors mk 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(29)tν|x·ν=±h≡±ν·C··∇xu|x·ν=±h=0,
where 2h denotes the overall thickness of the plate. Substituting representation (25) into the boundary conditions (29) yields the dispersion equation in the form
(30)∑k=16Ck(γkν·C·ν+ν·C·n)·mke±irγkh=0.Finally, (30) can be rewritten in a form better suited for numerical computations:(31)det((γ1ν·C·ν+ν·C·n)·m1e+irγ1h⋯(γ6ν·C·ν+ν·C·n)·m6e+irγ6h(γ1ν·C·ν+ν·C·n)·m1e-irγ1h⋯(γ6ν·C·ν+ν·C·n)·m6e-irγ6h)=0.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)u(x,t)=f(x′′)ei(rn·x-ωt)
with the unknown amplitude f regarded as a function of the imaginary coordinate x′′ defined by (8). Substituting representation (32) into equation of motion (26) yields
(33)(A1∂x′′2+(A2+A2t)∂x′′+A3)·f(x′′)=0,
where
(34)A1=ν·C·ν,A2=ν·C·n,A3=n·C·n-ρc2I.
The matrix A1 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:
(35)A1=(λ+2μ000μ000μ),A2=(0λ0μ00000),A3=(μ-ρc2000λ+2μ-ρc2000μ-ρc2).Up to the harmonic multiplierei(rn·x-ωt), 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)tν(x′′)≡(A1∂x′′+A2)·f(x′′).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 A1-1 yields the following expression for the derivative ∂x′′f(x′′):
(37)∂x′′f(x′′)=A1-1·tν(x′′)-A1-1·A2·f(x′′).Combining now (33)–(37) produces the desired equation of motion written in terms of vectors f and tν:
(38)∂x′′(ftν)=N·(ftν),
where N is the fundamental matrix [58]
(39)N=(-A1-1·A2A1-1A2t·A1-1·A2-A3-A2t·A1-1).
Now, the general solution of (38) can be represented in the form
(40)(ftν)=exp(x′′N)·C→6,
where C→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(41a)(f(+irh)0)=exp(+irhN)·C→6,(41b)(f(-irh)0)=exp(-irhN)·C→6.Excluding vector C→6 from (41a) and substituting the resultant expression into (41b) gives the following dispersion equation for an anisotropic plate:
(42)det((0,I)·exp(-2irhN))·(I0)=0.
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 multilayeredisotropic 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 Ck 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 containingn-homogeneous anisotropic layers, the matrix equation related to the global matrix method can be written in the following form:
(43)0→6n=((0,I)·D1+D1--D2+D2--D3+⋯⋯(0,I)·Dn-)·(C→1C→2C→3⋮C→n),
where C→k∈ℂ6,k=1,…,n, are six-dimensional generally complex vectors containing the unknown coefficients in representation (25); Dk± are 6×6-matrices defining displacements by (25) and surface tractions by (30) on both surfaces of a layer: (44)Dk±=(mk1e±irγk1hk⋯mk6e±irγk6hk(γk1ν·Ck·ν+ν·Ck·n)·mk1e±irγk1hk⋯(γk6ν·Ck·ν+ν·Ck·n)·mk6e±irγk6hk).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.


## 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
(1)cP∇divu-cSrotrotu=∂2u∂t2,
where u is the displacement field and cP and cS are velocities of the longitudinal and transverse bulk waves, respectively:
(2)cP=λ+2μρ,cS=μρ.In (2), λ and μ are Lamé constants and ρ is the material density. Then, the displacement field was represented in terms of scalar (Φ) and vector (Ψ) potentials
(3)u=∇Φ+rotΨ.The potentials were assumed to be harmonic in time:(4)Φ(x,t)=Φ′(x)eiωt,Ψ(x,t)=Ψ′(x)eiωt.Substituting representation (4) into (1) yields two independent Helmholtz equations:
(5)(Δ+ω2cP2)Φ′=0,(Δ+ω2cS2)Ψ′=0.To define the spatial periodicity and to simplify the analysis, the splitting spatial argument is needed:(6)x=(x·n)n+(x·ν)ν+(x·w)w,
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 argumentx·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(7)Φ′(x)=φ(x′′)ex′,Ψ′(x)=ψ(x′′)ex′,
where the dimensionless complex coordinates x′ and x′′ are
(8)x′=irx·n,x′′=irx·ν.In (8), i=-1 and r is the wave number related to the wavelength l by
(9)r=2πl.Substituting representations (7) into (5) results in the decoupled ordinary differential equations
(10)d2φdx′′2+(1-c2cP2)φ=0,d2ψdx′′2+(1-c2cS2)ψ=0,
where the phase speed c relates to the frequency and the wave number by the following relation:
(11)c=ωr.
The general solution of (10) can be written in the form
(12)φ(x′′)=C1sinh(γ1x′′)+C2cosh(γ1x′′),ψ(x′′)=C3sinh(γ2x′′)+C4cosh(γ2x′′),
where
(13)γ1=(1-c2cP2)1/2,γ2=(1-c2cS2)1/2.The unknown coefficients in (12) are defined (up to a multiplier) from the following boundary conditions on the free surfaces:
(14)tν≡(λtr(∇u)I+μ(∇u+∇ut))·ν=0,x·ν=±h,
where 2h is the depth of the plate. Substitution representation (3) into boundary conditions (14) yields boundary conditions written in terms of potentials φ′′ and ψ′′:
(15)(λΔΦ′I+2μ(∇∇Φ′+12(+(∇rotΨ′)t∇rotΨ′+(∇rotΨ′)t)(∇∇Φ′+12)(∇∇Φ′+12)·ν=0,x·ν=±h.
Substituting solutions (12) into (15), in view of Remark 1, yields the desired dispersion equation, found in [1, 3]:
(16)tanh(γ2rh)tanh(γ1rh)-(4γ1γ2(1+γ22)2)±1=0.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 c 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 atrh→∞ in (16) yields
(17)1=(4γ1γ2(1+γ22)2).
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
(18)c1,lim=cR.
Analysis of (16) at rh→0 (the long wave limit) yields the following equation [5]:
(19)γ2γ1=(4γ1γ2(1+γ22)2),
from where the values for two limiting velocities can be obtained:
(20)c2,lim=2cS1-cS2cp2.(21)c3,lim=cS.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
(22)cgroup=dωdr.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 F±(r,ω), and assuming that ω is a function of r, the derivative of (16) with respect to r takes the form
(23)∂F±(r,ω)∂r+∂F±(r,ω)∂ωdωdr=0,
from where, in view of (22), the secular equation for the group speed takes the form
(24)cgroup=-(∂F±(r,ω))/∂r(∂F±(r,ω))/∂ω.
In (24), r 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:
(25)u(x,t)=(∑k=16Ckmkeγkx′′︸uk(x′′))ei(rn·x-ωt),
where Ck are arbitrary complex coefficients determined up to a multiplier by satisfying conditions at the plate boundaries, uk is the displacement field of the kth partial wave, mk is a vectorial, generally, complex amplitude, determined by the Christoffel equation (this equation will be introduced later on), and γk is a root of the Christoffel equation. Note that, according to (8), the coordinate x′′ 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
(26)A(∂x,∂t)u≡divxC··∇xu-ρu¨=0,
where C is the fourth-order elasticity tensor assumed to be positive definite, yields the Christoffel equation
(27)[(γkν+n)·C·(n+γkν)-ρc2I]·mk=0,
where I is the unit diagonal matrix. Equation (27) admits the equivalent form
(28)det[(γkν+n)·C·(n+γkν)-ρc2I]=0.
The left-hand side of (28) represents a polynomial of degree six with respect to the Christoffel parameter γk. Equations (27), (28) show that roots γk and the corresponding eigenvectors mk can be considered as functions of the phase speed c.Remark 2.
(a) For Rayleigh waves, the rootsγk in representation (25) should be complex with Im(γk)<0; this ensures attenuation of Rayleigh wave in the “lower” half-space (ν·x)<0. The condition Im(γk)<0 confines the admissible speed interval and decreases the number of summation terms in (25). If Re(γk)=0 for all partial waves composing Rayleigh wave, then such a wave is called the genuine Rayleigh wave; if Re(γk)≠0 for some k, then such a wave is called the generalized Rayleigh wave [28]. For Lamb waves, the cases Re(γk)=0 and Re(γk)≠0 are usually not distinguished.
(b) Within the discussed formalism, the case of appearing multiple rootsγk and the coincident kernel eigenvectors mk 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(29)tν|x·ν=±h≡±ν·C··∇xu|x·ν=±h=0,
where 2h denotes the overall thickness of the plate. Substituting representation (25) into the boundary conditions (29) yields the dispersion equation in the form
(30)∑k=16Ck(γkν·C·ν+ν·C·n)·mke±irγkh=0.Finally, (30) can be rewritten in a form better suited for numerical computations:(31)det((γ1ν·C·ν+ν·C·n)·m1e+irγ1h⋯(γ6ν·C·ν+ν·C·n)·m6e+irγ6h(γ1ν·C·ν+ν·C·n)·m1e-irγ1h⋯(γ6ν·C·ν+ν·C·n)·m6e-irγ6h)=0.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)u(x,t)=f(x′′)ei(rn·x-ωt)
with the unknown amplitude f regarded as a function of the imaginary coordinate x′′ defined by (8). Substituting representation (32) into equation of motion (26) yields
(33)(A1∂x′′2+(A2+A2t)∂x′′+A3)·f(x′′)=0,
where
(34)A1=ν·C·ν,A2=ν·C·n,A3=n·C·n-ρc2I.
The matrix A1 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:
(35)A1=(λ+2μ000μ000μ),A2=(0λ0μ00000),A3=(μ-ρc2000λ+2μ-ρc2000μ-ρc2).Up to the harmonic multiplierei(rn·x-ωt), 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)tν(x′′)≡(A1∂x′′+A2)·f(x′′).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 A1-1 yields the following expression for the derivative ∂x′′f(x′′):
(37)∂x′′f(x′′)=A1-1·tν(x′′)-A1-1·A2·f(x′′).Combining now (33)–(37) produces the desired equation of motion written in terms of vectors f and tν:
(38)∂x′′(ftν)=N·(ftν),
where N is the fundamental matrix [58]
(39)N=(-A1-1·A2A1-1A2t·A1-1·A2-A3-A2t·A1-1).
Now, the general solution of (38) can be represented in the form
(40)(ftν)=exp(x′′N)·C→6,
where C→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(41a)(f(+irh)0)=exp(+irhN)·C→6,(41b)(f(-irh)0)=exp(-irhN)·C→6.Excluding vector C→6 from (41a) and substituting the resultant expression into (41b) gives the following dispersion equation for an anisotropic plate:
(42)det((0,I)·exp(-2irhN))·(I0)=0.
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 multilayeredisotropic 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 Ck 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 containingn-homogeneous anisotropic layers, the matrix equation related to the global matrix method can be written in the following form:
(43)0→6n=((0,I)·D1+D1--D2+D2--D3+⋯⋯(0,I)·Dn-)·(C→1C→2C→3⋮C→n),
where C→k∈ℂ6,k=1,…,n, are six-dimensional generally complex vectors containing the unknown coefficients in representation (25); Dk± are 6×6-matrices defining displacements by (25) and surface tractions by (30) on both surfaces of a layer: (44)Dk±=(mk1e±irγk1hk⋯mk6e±irγk6hk(γk1ν·Ck·ν+ν·Ck·n)·mk1e±irγk1hk⋯(γk6ν·Ck·ν+ν·Ck·n)·mk6e±irγk6hk).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)v(x′′)=∂x′′f(x′′),
the equation of motion (26) can be reduced to the Cauchy normal form
(46)∂x′′(fv)=G·(fv),
where
(47)G=(0I-A1·A3-A1·(A2+A2t))
and matrices A1, A2, and A3 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
(48)(f(x′′)v(x′′))=exp(x′′G)·C→6,
where C→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
(49)(f(±irh)tν(±irh))=Z·exp(±irhG)·C→6,
where
(50)Z=(I0A2A1).
Note that matrix Z is nondegenerate due to positive definiteness of the elasticity tensor.Equation (49) allows us to exclude C→6 by expressing the displacements and surface tractions on one side of the plate through the corresponding vectors on the other side; that yields
(51)(f(-irh)tν(-irh))=T·(f(+irh)tν(+irh)),
where
(52)T=Z·exp(-2irhG)·Z-1.
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.


### 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 c2,lim at r→0.


#### 2.2.1. Traction-Free Plate

If both sides of a plate are traction free, then(53)tν(±irh)=0.
Condition (53) means that a 3×3-mapping
(54)R3→(B1,B2)·T·(B2B1)R3
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)B1=0,B2=I.MatricesB1 and B2 are introduced for consistency with other types of boundary conditions that will be considered later in this section.The condition (54) is equivalent to
(56)det((B1,B2)·T·(B2B1))=0.
Thus, (56) is the desired dispersion equation for a plate with free boundaries.Taking in (54) limit at r→0 does not lead to a meaningful condition, because of vanishing determinant in (56) at r=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 c2,lim:
(57)det((B1,B2)·Z·G·Z-1·(B2B1))=0.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(58)f(±irh)=0.
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 B1 and B2:
(59)B1=I,B2=0.
The equation for the limiting speed c2,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(60)ν·tν(±irh)=0,(I-ν⊗ν)·f(±irh)=0.
For boundary conditions (60), matrices B1 and B2 become
(61)B1=I-ν⊗ν,B2=ν·tν.
Similar to expressions (60), boundary conditions (ii) can be written in the form
(62)(I-ν⊗ν)·tν(±irh)=0,ν·f(±irh)=0,
and B1, B2 become
(63)B1=ν·tν,B2=I-ν⊗ν.


### 2.3. Dispersion Relations for a Multilayered Plate

Consider ann-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
(64)(f(-irhn)tν(-irhn))=(∏k=1nTk)·(f(+irh1)tν(+irh1)),
where the lower index indicates the number of a layer and the transfer matrices Tk are defined by (52).Now, multiplying both sides of (64) by matrices B1, B2, similar to (54)–(56) for the traction-free plate, we arrive at the dispersion equation for the multilayered traction-free plate:
(65)det((B1,B2)·(∏k=1nTk)·(B2B1))=0.
Similar to (54), (65) reflects degeneracy of the mapping
(66)R3→(B1,B2)·(∏k=1nTk)·(B2B1)R3.
The dispersion equations for other types of boundary conditions have the same form as (65) but with obvious replacement of matrices B1 and B2.As it was with (57), taking limit at r→0 of (65) does not lead to a meaningful condition, because of vanishing the determinant at r=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:
(67)0→6n=((B1,B2)·D1+D1--D2+D2--D3+⋯⋯(B1,B2)·Dn-)·(C→1C→2C→3⋮C→n).
Within the considered Cauchy formalism, matrices Dk± change:
(68)Dk±=Zk·exp(±irhkGk),
where index k indicates the corresponding layer.


## 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)v(x′′)=∂x′′f(x′′),
the equation of motion (26) can be reduced to the Cauchy normal form
(46)∂x′′(fv)=G·(fv),
where
(47)G=(0I-A1·A3-A1·(A2+A2t))
and matrices A1, A2, and A3 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
(48)(f(x′′)v(x′′))=exp(x′′G)·C→6,
where C→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
(49)(f(±irh)tν(±irh))=Z·exp(±irhG)·C→6,
where
(50)Z=(I0A2A1).
Note that matrix Z is nondegenerate due to positive definiteness of the elasticity tensor.Equation (49) allows us to exclude C→6 by expressing the displacements and surface tractions on one side of the plate through the corresponding vectors on the other side; that yields
(51)(f(-irh)tν(-irh))=T·(f(+irh)tν(+irh)),
where
(52)T=Z·exp(-2irhG)·Z-1.
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.


## 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 c2,lim at r→0.


### 2.2.1. Traction-Free Plate

If both sides of a plate are traction free, then(53)tν(±irh)=0.
Condition (53) means that a 3×3-mapping
(54)R3→(B1,B2)·T·(B2B1)R3
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)B1=0,B2=I.MatricesB1 and B2 are introduced for consistency with other types of boundary conditions that will be considered later in this section.The condition (54) is equivalent to
(56)det((B1,B2)·T·(B2B1))=0.
Thus, (56) is the desired dispersion equation for a plate with free boundaries.Taking in (54) limit at r→0 does not lead to a meaningful condition, because of vanishing determinant in (56) at r=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 c2,lim:
(57)det((B1,B2)·Z·G·Z-1·(B2B1))=0.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(58)f(±irh)=0.
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 B1 and B2:
(59)B1=I,B2=0.
The equation for the limiting speed c2,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(60)ν·tν(±irh)=0,(I-ν⊗ν)·f(±irh)=0.
For boundary conditions (60), matrices B1 and B2 become
(61)B1=I-ν⊗ν,B2=ν·tν.
Similar to expressions (60), boundary conditions (ii) can be written in the form
(62)(I-ν⊗ν)·tν(±irh)=0,ν·f(±irh)=0,
and B1, B2 become
(63)B1=ν·tν,B2=I-ν⊗ν.


## 2.2.1. Traction-Free Plate

If both sides of a plate are traction free, then(53)tν(±irh)=0.
Condition (53) means that a 3×3-mapping
(54)R3→(B1,B2)·T·(B2B1)R3
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)B1=0,B2=I.MatricesB1 and B2 are introduced for consistency with other types of boundary conditions that will be considered later in this section.The condition (54) is equivalent to
(56)det((B1,B2)·T·(B2B1))=0.
Thus, (56) is the desired dispersion equation for a plate with free boundaries.Taking in (54) limit at r→0 does not lead to a meaningful condition, because of vanishing determinant in (56) at r=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 c2,lim:
(57)det((B1,B2)·Z·G·Z-1·(B2B1))=0.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(58)f(±irh)=0.
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 B1 and B2:
(59)B1=I,B2=0.
The equation for the limiting speed c2,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(60)ν·tν(±irh)=0,(I-ν⊗ν)·f(±irh)=0.
For boundary conditions (60), matrices B1 and B2 become
(61)B1=I-ν⊗ν,B2=ν·tν.
Similar to expressions (60), boundary conditions (ii) can be written in the form
(62)(I-ν⊗ν)·tν(±irh)=0,ν·f(±irh)=0,
and B1, B2 become
(63)B1=ν·tν,B2=I-ν⊗ν.


## 2.3. Dispersion Relations for a Multilayered Plate

Consider ann-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
(64)(f(-irhn)tν(-irhn))=(∏k=1nTk)·(f(+irh1)tν(+irh1)),
where the lower index indicates the number of a layer and the transfer matrices Tk are defined by (52).Now, multiplying both sides of (64) by matrices B1, B2, similar to (54)–(56) for the traction-free plate, we arrive at the dispersion equation for the multilayered traction-free plate:
(65)det((B1,B2)·(∏k=1nTk)·(B2B1))=0.
Similar to (54), (65) reflects degeneracy of the mapping
(66)R3→(B1,B2)·(∏k=1nTk)·(B2B1)R3.
The dispersion equations for other types of boundary conditions have the same form as (65) but with obvious replacement of matrices B1 and B2.As it was with (57), taking limit at r→0 of (65) does not lead to a meaningful condition, because of vanishing the determinant at r=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:
(67)0→6n=((B1,B2)·D1+D1--D2+D2--D3+⋯⋯(B1,B2)·Dn-)·(C→1C→2C→3⋮C→n).
Within the considered Cauchy formalism, matrices Dk± change:
(68)Dk±=Zk·exp(±irhkGk),
where index k indicates the corresponding layer.


## 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
(64)(f(-irhn)tν(-irhn))=(∏k=1nTk)·(f(+irh1)tν(+irh1)),
where the lower index indicates the number of a layer and the transfer matrices Tk are defined by (52).Now, multiplying both sides of (64) by matrices B1, B2, similar to (54)–(56) for the traction-free plate, we arrive at the dispersion equation for the multilayered traction-free plate:
(65)det((B1,B2)·(∏k=1nTk)·(B2B1))=0.
Similar to (54), (65) reflects degeneracy of the mapping
(66)R3→(B1,B2)·(∏k=1nTk)·(B2B1)R3.
The dispersion equations for other types of boundary conditions have the same form as (65) but with obvious replacement of matrices B1 and B2.As it was with (57), taking limit at r→0 of (65) does not lead to a meaningful condition, because of vanishing the determinant at r=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:
(67)0→6n=((B1,B2)·D1+D1--D2+D2--D3+⋯⋯(B1,B2)·Dn-)·(C→1C→2C→3⋮C→n).
Within the considered Cauchy formalism, matrices Dk± change:
(68)Dk±=Zk·exp(±irhkGk),
where index k 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 matricesGk,k=1,…,n be decomposed into their Jordan normal forms
(69)Gk=Wk-1·Jk·Wk,
where Wk is a matrix composed of eigenvectors and possibly the generalized eigenvectors, while Jk is a matrix containing eigenvalues and possibly Jordan blocks if Gk is not a semisimple matrix.With decomposition (69), the exponential matrices in (52) can be represented in the following form [79, 80]:
(70)exp(-2irhkGk)=Wk-1·exp(-2irhkJk)·Wk.Let(λk(m))m=1,…,6 denote (not necessary aliquant) eigenvalues of Gk, and
(71)a=max︸k,mRe(-2irhkλk(m)).
Multiplying all the exponent matrices in (65) by e-a 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:(72)D=diag(ε,ε-,α),
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:(73)D=(1+ε(t)-1-11-ε(t)),
where |ε(t)|<10-15 is a small function of the argument t. It is obvious that computing the determinant with double precision accuracy gives zero, inspite of the applied numerical algorithm. However, analytical evaluation yields
(74)detD=ε2(t),
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:(75)δ≡|C1111C1122+2C1212-1|.In (75), C1122 and C1212 would correspond to Lamé parameters λ and μ, and C1111 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(76)C1111=1,C1122=0.1,C1212=0.5.
For the considered cubic crystal, the measuring parameter is relatively small: δ≈0.09.The second cubic crystal has the following elastic parameters:(77)C1111=1,C1122=0,C1212=0.25.
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]:
(78)C1111=1,C1122=-0.1,C1212=0.5.
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):
(79)ν=(1;0;0),n=(0;cosψ;sinψ).
Variation of the limiting velocity c2,lim with respect to angle ψ is shown in Figure 1.Figure 1
Variation of limiting velocityc2,lim versus angle ψ.The plots in Figure1 reveal that the limiting velocity c2,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 c2,lim values with respect to ψ.Remark 5.
In the NDT, the limiting Lamb wave propagating with the velocityc2,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 ish=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/m3 and the following elastic constants (in GPa) [82]:
(80)C1111=503,C1122=109,C1133=52C2222=503,C2233=52,C3333=552C2323=C1313=160,C1212=C1111-C11222;
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 Si3N4 (γ-phase) with ρ=2329 kg/m3 and elastic constants (in GPa) [83]
(81)C1111=514,C1122=182,C1133=182C2222=514,C2233=182,C3333=514C1212=C1313=C2323=327.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 normaln are
(82)SiC:cPn=12498cSnin-plane=7822cSnoutofplane=7049Si3N4:cPn=14856cSn=11849.
For the considered multilayered plate, the dispersion curves for Lamb waves are plotted in Figure 2.Dispersion curves for 3-layerd (a) and 11-layered (b) traction-free plates with alternating SiC-Si3N4 layers.
(a)
(b)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 thatc2,lim<cPnSi3N4. 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 cRnSiC/cPnSi3N4≈0.497, where cRnSiC denotes velocity of Rayleigh wave in SiC for the considered propagation direction.


## 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 matricesGk,k=1,…,n be decomposed into their Jordan normal forms
(69)Gk=Wk-1·Jk·Wk,
where Wk is a matrix composed of eigenvectors and possibly the generalized eigenvectors, while Jk is a matrix containing eigenvalues and possibly Jordan blocks if Gk is not a semisimple matrix.With decomposition (69), the exponential matrices in (52) can be represented in the following form [79, 80]:
(70)exp(-2irhkGk)=Wk-1·exp(-2irhkJk)·Wk.Let(λk(m))m=1,…,6 denote (not necessary aliquant) eigenvalues of Gk, and
(71)a=max︸k,mRe(-2irhkλk(m)).
Multiplying all the exponent matrices in (65) by e-a 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:(72)D=diag(ε,ε-,α),
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:(73)D=(1+ε(t)-1-11-ε(t)),
where |ε(t)|<10-15 is a small function of the argument t. It is obvious that computing the determinant with double precision accuracy gives zero, inspite of the applied numerical algorithm. However, analytical evaluation yields
(74)detD=ε2(t),
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.1.1. The Overflow Problem

This problem can be relatively easily solved by a suitable normalization. Let the fundamental matricesGk,k=1,…,n be decomposed into their Jordan normal forms
(69)Gk=Wk-1·Jk·Wk,
where Wk is a matrix composed of eigenvectors and possibly the generalized eigenvectors, while Jk is a matrix containing eigenvalues and possibly Jordan blocks if Gk is not a semisimple matrix.With decomposition (69), the exponential matrices in (52) can be represented in the following form [79, 80]:
(70)exp(-2irhkGk)=Wk-1·exp(-2irhkJk)·Wk.Let(λk(m))m=1,…,6 denote (not necessary aliquant) eigenvalues of Gk, and
(71)a=max︸k,mRe(-2irhkλk(m)).
Multiplying all the exponent matrices in (65) by e-a 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:(72)D=diag(ε,ε-,α),
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:(73)D=(1+ε(t)-1-11-ε(t)),
where |ε(t)|<10-15 is a small function of the argument t. It is obvious that computing the determinant with double precision accuracy gives zero, inspite of the applied numerical algorithm. However, analytical evaluation yields
(74)detD=ε2(t),
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:(75)δ≡|C1111C1122+2C1212-1|.In (75), C1122 and C1212 would correspond to Lamé parameters λ and μ, and C1111 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(76)C1111=1,C1122=0.1,C1212=0.5.
For the considered cubic crystal, the measuring parameter is relatively small: δ≈0.09.The second cubic crystal has the following elastic parameters:(77)C1111=1,C1122=0,C1212=0.25.
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]:
(78)C1111=1,C1122=-0.1,C1212=0.5.
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):
(79)ν=(1;0;0),n=(0;cosψ;sinψ).
Variation of the limiting velocity c2,lim with respect to angle ψ is shown in Figure 1.Figure 1
Variation of limiting velocityc2,lim versus angle ψ.The plots in Figure1 reveal that the limiting velocity c2,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 c2,lim values with respect to ψ.Remark 5.
In the NDT, the limiting Lamb wave propagating with the velocityc2,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 ish=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/m3 and the following elastic constants (in GPa) [82]:
(80)C1111=503,C1122=109,C1133=52C2222=503,C2233=52,C3333=552C2323=C1313=160,C1212=C1111-C11222;
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 Si3N4 (γ-phase) with ρ=2329 kg/m3 and elastic constants (in GPa) [83]
(81)C1111=514,C1122=182,C1133=182C2222=514,C2233=182,C3333=514C1212=C1313=C2323=327.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 normaln are
(82)SiC:cPn=12498cSnin-plane=7822cSnoutofplane=7049Si3N4:cPn=14856cSn=11849.
For the considered multilayered plate, the dispersion curves for Lamb waves are plotted in Figure 2.Dispersion curves for 3-layerd (a) and 11-layered (b) traction-free plates with alternating SiC-Si3N4 layers.
(a)
(b)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 thatc2,lim<cPnSi3N4. 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 cRnSiC/cPnSi3N4≈0.497, where cRnSiC denotes velocity of Rayleigh wave in SiC for the considered propagation direction.


## 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:(75)δ≡|C1111C1122+2C1212-1|.In (75), C1122 and C1212 would correspond to Lamé parameters λ and μ, and C1111 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(76)C1111=1,C1122=0.1,C1212=0.5.
For the considered cubic crystal, the measuring parameter is relatively small: δ≈0.09.The second cubic crystal has the following elastic parameters:(77)C1111=1,C1122=0,C1212=0.25.
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]:
(78)C1111=1,C1122=-0.1,C1212=0.5.
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):
(79)ν=(1;0;0),n=(0;cosψ;sinψ).
Variation of the limiting velocity c2,lim with respect to angle ψ is shown in Figure 1.Figure 1
Variation of limiting velocityc2,lim versus angle ψ.The plots in Figure1 reveal that the limiting velocity c2,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 c2,lim values with respect to ψ.Remark 5.
In the NDT, the limiting Lamb wave propagating with the velocityc2,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 ish=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/m3 and the following elastic constants (in GPa) [82]:
(80)C1111=503,C1122=109,C1133=52C2222=503,C2233=52,C3333=552C2323=C1313=160,C1212=C1111-C11222;
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 Si3N4 (γ-phase) with ρ=2329 kg/m3 and elastic constants (in GPa) [83]
(81)C1111=514,C1122=182,C1133=182C2222=514,C2233=182,C3333=514C1212=C1313=C2323=327.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 normaln are
(82)SiC:cPn=12498cSnin-plane=7822cSnoutofplane=7049Si3N4:cPn=14856cSn=11849.
For the considered multilayered plate, the dispersion curves for Lamb waves are plotted in Figure 2.Dispersion curves for 3-layerd (a) and 11-layered (b) traction-free plates with alternating SiC-Si3N4 layers.
(a)
(b)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 thatc2,lim<cPnSi3N4. 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 cRnSiC/cPnSi3N4≈0.497, where cRnSiC 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.



---

*Source: 698706-2013-09-16.xml*