Surface adapted partial waves for the description of elastic vibrations in bilayered plates

Wave Motion 92 (2020) 102430

Contents lists available at ScienceDirect

Wave Motion journal homepage: www.elsevier.com/locate/wamot

Surface adapted partial waves for the description of elastic vibrations in bilayered plates S. Cojocaru Horia Hulubei National Institute For Physics And Nuclear Engineering, RO-077125, Magurele, Romania

article

info

Article history: Received 8 October 2018 Received in revised form 5 April 2019 Accepted 19 September 2019 Available online 21 September 2019

a b s t r a c t An approach is proposed in which the form of the partial waves in the decomposition of the displacement field is chosen in such a way that boundary conditions on the outer surfaces of a layered structure are matched exactly. So found basis functions are further used to solve for the remaining boundary conditions. Explicit analytic expressions for the vibration spectrum and normalized amplitudes of a free plate in the long wavelength regime are derived and their behavior in the full space of material parameters is discussed. Thus, it is found that properties of the fundamental flexural mode change in a strongly non-monotonous way. In particular, by increasing the thickness of added layer one may observe the same propagation velocity for up to three different values of layer thicknesses ratio. Existence of an important criterion based on the relative signs taken by the amplitudes on the outer surfaces in this regime is revealed. It allows to establish a clear correspondence between symmetric and antisymmetric solutions for a single layer and those for a bilayered plate. It is suggested that the criterion is not limited to case of two layers. Branches of the spectrum are classified accordingly and their evolution from long to short wavelength is discussed. Analysis of a specific example which includes volume and surface modes (Rayleigh and Stoneley) is presented. The results demonstrate the efficiency of the proposed scheme and possible extensions are suggested. © 2019 Elsevier B.V. All rights reserved.

1. Introduction Propagation of waves in confined systems is of a quite general relevance in physics and its description in terms of continuum medium is applicable to a broad range of problems from macro- to nanoscale. For mechanical vibrations elasticity theory has developed a rich variety of models and methods, e.g., [1–8], describing wave propagation. Thus, as was shown by Rayleigh, for a half-infinite medium the lowest energy excitations are surface guided acoustic waves where amplitude of particle motion exponentially decays on a wavelength distance from the surface. Such guided waves are much less dissipated with propagation distance compared to bulk waves [2] and have important applications in areas where the relevant length-scale varies from that of geophysics [9] to nanophysics [10,11] or from those related to classical material properties [12] to those studying quantum information [13,14]. Modern technology utilizes modulated structures, e.g., phononic crystals, where propagation of acoustic waves through a plate of finite thickness results in a rather complex dispersion including multiple band gaps [15,16]. In many cases, micro and nanoelectronic devices like ultrasensitive detectors or electronic micro-coolers [17–19] require low temperature regimes where lower energy and longer wavelength part of the excitation spectra plays the leading role. At higher temperatures, when a large number E-mail address: [email protected]. https://doi.org/10.1016/j.wavemoti.2019.102430 0165-2125/© 2019 Elsevier B.V. All rights reserved.

2

S. Cojocaru / Wave Motion 92 (2020) 102430

Fig. 1. Bilayered free standing slab of thickness d with rigidly bonded interface.

of higher energy modes are excited, the statistical weight of the lowest ones may still be dominating [20]. Thus, in graphene at room temperatures the three fundamental acoustic phonon modes are responsible for its superior thermal conductivity [21,22]. When Rayleigh wavelength reaches the confinement scale the wave becomes dispersive and the finite size and surface related effects come to the forefront. The reference example are Lamb’s modes in a single layer of an isotropic material when multiple reflections from the two boundaries are relevant. Particles of the medium move within the sagittal plane so that the polarization pattern is a mixture of longitudinal and shear vertical displacements and in the short wavelength limit one recovers the Rayleigh equation [1,5]. Due to the symmetry with respect to the middle plane there exist two classes of Lamb’s solutions, symmetric and antisymmetric, and their frequencies become degenerate in this limit. However, the knowledge of the frequency dispersion is often insufficient and one also needs a detailed knowledge of the wave pattern. For instance, in micro- and nanodevices, the proper quantization and normalization of elastic waves is required to study their interaction with other quantum excitations, e.g., electronic or magnetic [10,23,24]. The complexity of the physical picture increases for inhomogeneous anisotropic or layered structures so that new analytical methods are being developed to meet the challenges, [7,8,25–28]. In many cases of practical interest the structure, e.g., multilayered or sandwiched, still preserves the symmetry with respect to the middle plane [29–33]. The simplest extension of Lamb’s model where the intrinsic inhomogeneity of the interface (z = 0 in Fig. 1) breaks the symmetry of the structure is that of a free standing plate of finite thickness d composed of two rigidly bonded layers, da + db = d. The index γ = {a, b} identifies the layer: the upper and lower boundaries correspond to z = da and z = −db , respectively. This model includes the key characteristics of many actual devices (sensors, spintronic microchips etc. [34]) where rigid bonding is possible to realize even for a very large lattice mismatch [35]. The boundary between two half-spaces can serve as a waveguide for the propagation of Stoneley or Love waves [5,36–39] and, similar to Rayleigh modes, in a confined layered system these would become dispersive when the wavelength compares to d. However, this regime has been studied mostly by numerical methods for specific values of parameters. Below it is shown that a different treatment of boundary conditions and generalization of the symmetry properties of Lamb’s modes allows to significantly simplify the problem. The explicit solution for the long wavelength part of the spectrum obtained in this way allows to grasp the general behavior in the full parametric space and recovers the Lamb solution in the respective limit. field vector U(r, t) at a point r at time t is linearly related to the strain tensor εkm (r, t) = ) ( The displacement ∂ Uk ∂ rm

∂ Um ∂ rk

/2. For an isotropic medium it is linearly related to the stress Σij by means of the stiffness tensor λijkm = ( ) λδij δkm + µ δik δjm + δim δjk , where λ and µ are the two Lamé parameters and δij is the Kronecker δ−symbol with {i, j, k, m} = x, y, z: +

Σij = λδij εkk + 2µεij .

(1)

This expression is then substituted in the equations of elastodynamics

ρ

∂ 2 Ui (r, t) ∂ Σij (r, t) = . 2 ∂t ∂ rj

(2)

It is further assumed that material parameters (mass densities ρa , ρb etc.) take constant values in each layer. 2. The bilayered plate For the geometry in Fig. 1 the displacement can be represented by a superposition of normal modes, n, with eigenfrequencies ωn : U(r, t) =

∑∫

un z , q∥ , ωn exp iq∥ r∥ − iωn t

(

)

(

n

) dq∥ , (2π )2

(3)

and similarly for the stress tensor

Σij (r, t) =

∑∫ n

σij (z ) exp (iqx − iωn t )

dq 2π

.

(4)

S. Cojocaru / Wave Motion 92 (2020) 102430

3

Here the in-plane wavevector, q∥ , is used to fix the direction of the x - axis, so that q∥ r∥ = qx x + qy y = qx, and the n index will be dropped henceforth. The equations decouple according to the polarization pattern into the shear horizontal and mixed, shear vertical and longitudinal [40]. We will further consider only the latter ones u (z ) = (ux , 0, uz ) ,

(5)

and introduce the following representation for the amplitudes ui (z) = θ (z ) ui,a (z) + (1 − θ (z )) ui,b (z), σxz (z ) = θ (z ) σxz ,a (z ) + (1 − θ (z )) σxz ,b (z ) ,

σzz (z ) = θ (z ) σzz ,a (z ) + (1 − θ (z )) σzz ,b (z ) ,

(6)

where the second index identifies the layer and θ (z ) is the Heaviside θ− function. From (2)–(4) we obtain the relations between displacement and stress amplitudes

∂σxz , ∂z ∂σzz . − ρ uz ω2 = iqσxz + ∂z

− ρ ux ω2 = iqσxx +

(7)

The contribution of the derivatives in (7) contains terms with Dirac delta-function

[ ] δ (z ) σiz ,a (z ) − σiz ,b (z ) ,

(8)

which vanish by continuity condition on the interface. The remaining terms have either θ (z ) or (1 − θ (z )) as a prefactor and refer to one of the two layers where the material parameters are constant. Also note that σxx (z ) as well as σyy (z ) are generally not continuous functions of z: the respective tensor components represent forces acting along the normals to the surfaces, respectively, x = const and y = const. Since these forces are parallel to the interface, they are counterbalanced by equal and opposite forces within each layer and may therefore differ across the interface. The above conclusions are valid for anisotropic media as well. From (1) and (4) we obtain iσzz = λ

σxz = µ

(

) ( ) ∂ (iuz ) ∂ (iuz ) − qux + 2µ , ∂z ∂z

(

) ∂ ux + q (iuz ) . ∂z

(9)

These expressions are consistent with the representation (6) by virtue of the property (8) with displacement amplitudes in the square bracket. After taking the remaining derivatives in (2) we recover the commonly used form of the equations describing Rayleigh–Lamb waves in stratified isotropic media with rigid bonding, e.g., [6,41]:

∂ 2 ux,γ ∂ iuz ,γ − (λγ + 2µγ )q2 wγ2 ux,γ + (λγ + µγ )q = 0, (10) ∂ z2 ∂z ∂ 2 iuz ,γ ∂ ux,γ − µγ q2 vγ2 iuz ,γ − (λγ + µγ )q = 0, (λγ + 2µγ ) ∂ z2 ∂z where v and w are √ √ ( )2 ( )2 c c wγ = 1 − , vγ = 1 − . (11) ℓγ sγ √ √ Here c = ω/q is the phase velocity and ℓ = (λ + 2µ) /ρ , s = µ/ρ are the bulk longitudinal and transverse sound µγ

velocities in the respective materials. There are overall 8 equations associated with the boundary conditions on the outer surfaces,

) ( ) ( σxz ,γ z = dγ = 0, σzz ,γ z = dγ = 0, ) γ = a; dγ = −db , γ = b. (

da ,

(12)

(13)

and on the solid–solid interface ux,a (z = 0) = ux,b (z = 0) , uz ,a (z = 0) = uz ,b (z = 0) ,

σzx,a (z = 0) = σzx,b (z = 0) , σzz ,a (z = 0) = σzz ,b (z = 0) .

(14)

4

S. Cojocaru / Wave Motion 92 (2020) 102430

The dimensionless quantities v and w satisfy the identity

λ + µ(1 + v 2 ) = (λ + 2µ)w2 ,

(15)

and play a key role in the behavior of the amplitudes u (z ). For instance, if c < sγ , ℓγ , v and w take real values and the wave has therefore an attenuated character and may be described by hyperbolic functions. Note that phase velocity c = c (q) diverges in the long wavelength limit for the gapped (ω (q −→ 0) ̸ = 0, “pseudo optical” or “cutoff”) modes [2]. As a consequence, the above inequality is reversed and all the solutions for these plate resonance modes are volume waves described by harmonic functions. An equivalent representation

{

κℓ = iqw,

}

κs = iqv,

(16)

is also often used in the literature. Since the Rayleigh–Lamb waves are composed of the two polarizations, longitudinal and shear vertical, these quantities correspond to the wavevector angles of the two polarization components of the wave for a given in-plane wavenumber q : cot θs = κs /q and cot θℓ = κℓ /q, e.g., [40]. For instance, when a longitudinal wave is incident on the outer surface at an angle θℓ , it is converted into a shear vertical one, outgoing at an angle θs , and a specularly reflected longitudinal wave. Therefore one may interpret the identity (15) or its equivalent

) ( ) ( ω2 = s2 q2 + κs2 = ℓ2 q2 + κℓ2 ,

(17)

as a form of Snell law for the Lamb waves. The set of discrete solution values of v and w for a given q corresponds to the branches of the vibration spectrum. 3. Surface adapted partial waves The standard approach to the problem of layered structures is to use the partial wave decomposition, which can be applied either directly to the displacement field or to its Helmholtz representation, i.e., to the curl-free and divergencefree potentials, U(r, t) = ∇φ + ∇ × ψ [1,39]. Description is usually formulated in terms of incoming and outgoing partial plane or harmonic waves [42] which are a convenient basis set when considering the system embedded in an infinite medium [41,43]. In our case for the Eqs. (10)–(14) the minimal basis contains 8 partial waves: direct and reflected longitudinal components with ±wγ (or ±qℓ,γ according to (16)) and, similarly, the shear wave components with ±vγ . This basis allows to map the problem into one of a linear algebra and results in a 8 × 8 matrix of the secular equation, e.g., [44,45]. Although the complexity of the problem keeps growing when the model is extended to account for anisotropy or larger number of layers, advanced methods of solving such matrices have been developed during the years which allow to greatly reduce the computational effort, e.g., [46,47]. The proposed approach adopts a different strategy: instead of operating with plane wave basis one first finds such a form of partial waves which satisfies the conditions on the outer boundaries “automatically”, i.e., exactly, even if the solution of the problem is found with approximation. The problem is then reduced to solving only for the conditions on the interfaces. Aside from leading directly to a reduced matrix problem, the approach reveals some exact properties of the wave amplitudes on the outer boundaries which allow to establish an important connection to the reference Lamb model. The scheme of the approach contains a possibility of extension to layered structures with a curved surface, so that the present example of a plate geometry may be considered as its simplest realization. Thus, the amplitudes are expressed in the form where the equation of the plate surface enters as the argument of the partial wave ux,γ = ξw,γ sinh((z − dγ )qwγ ) + ξv,γ sinh((z − dγ )qvγ ) + θw,γ cosh((z − dγ )qwγ ) + θv,γ cosh((z − dγ )qvγ ), iuz ,γ = ζw,γ sinh((z − dγ )qwγ ) + ζv,γ sinh((z − dγ )qvγ ) + ηw,γ cosh((z − dγ )qwγ ) + ηv,γ cosh((z − dγ )qvγ ).

(18)

It is clear that in this way the explicit dependence on dγ is eliminated from (12) for the components of the stress tensor (9) since arguments of functions in (18) define the respective surfaces. At first glance the approach seems counterproductive since we have now twice the number of unknown coefficients {ξ , θ, . . .}, i.e. 16, and, in contrast to the method of potentials, the two polarizations remain coupled, Eqs. (10). However, this “excessive” freedom can now be used to adapt our partial waves as required. It can then be easily seen that this choice produces massive simplifications upon substitution into the wave equations (10)

ηw = wξw , ηv = ξv /v,

ζw = wθw , ζv = θv /v,

and (12) 1 + v2

( ξw = −

2wv

) ξv , θw = − (

2θv

1 + v2

).

(19)

For obvious reasons these relations are the same for both layers and we have therefore dropped the γ -index. We are now left with only 4 unknown coefficients which are henceforth relabeled as Xγ = ξv,γ , Yγ = θv,γ .

S. Cojocaru / Wave Motion 92 (2020) 102430

5

This procedure can be extended for the case of a larger number of layers, each additional layer contributing with 4 additional unknown coefficients since the new boundaries are of the interface type which contribute with 4 conditions as in Eqs. (14). Expressions for the components of the amplitudes ux , uz and stresses σxz , σzz in the bilayer system become:

(

(

ux,γ (q, z ) = Xγ (q) sinh

z − dγ qvγ −

((

)

)

1 + vγ2

2

+ Yγ (q) cosh z − dγ qvγ − (

)

sinh

2wγ vγ

( ((

)

)

(

iuz ,γ (q, z ) vγ = Xγ (q) cosh

z − dγ qvγ −

((

)

)

1 + vγ2

2wγ vγ

+ Yγ (q) sinh z − dγ qvγ − (

((

)

)

) )

) cosh z − dγ qwγ

1 + vγ2

)

) )

,

)

) cosh

2

(

z − dγ qwγ

)

((

1 + vγ2

(

((

((

z − dγ qwγ

)

) sinh z − dγ qwγ ((

)

)

) )

,

)) ) ) (( ) ( (( vγ σxz ,γ (q, z ) ( ) = Xγ (q) cosh z − dγ qvγ − cosh z − dγ qwγ 1 + vγ2 µγ q ( ) ) ) (( ) ) (( 4wγ vγ + Yγ (q) sinh z − dγ qvγ − ( )2 sinh z − dγ qwγ , 1 + vγ2 iσzz ,γ (q, z ) 2µγ q

(

1 + vγ2

(

= Xγ (q) sinh z − dγ qvγ −

)

((

)

)2

4wγ vγ

(20)

) sinh

((

z − dγ qwγ

)

)

( (( ) ) (( ) )) + Yγ (q) cosh z − dγ qvγ − cosh z − dγ qwγ . Note that the amplitudes on the outer surfaces simplify to the following expressions ux,γ q, z = dγ

)

= −Yγ (q)

iuz ,γ q, z = dγ

)

= Xγ (q)

(

(

1 − vγ2 1 + vγ2

1 − vγ2 2vγ

,

.

(21)

The normal mode spectrum is now determined solely by the four continuity conditions on the interface: 1. ux,a (z = 0) = ux,b (z = 0) , 2. uz ,a (z = 0) = uz ,b (z = 0) , 3. σxz ,a (z = 0) = σxz ,b (z = 0) , 4. σzz ,a (z = 0) = σzz ,b (z = 0) .

(22)

The explicit { form of }these equations is given in the Appendix (note that it differs from the equations obtained in [44]). Solutions ω, Xγ , Yγ are functions of q and of the parameters for the a and b layers (we will occasionally indicate this explicitly as, e.g., Xa = Xa (a, b)). The phase velocity c }and three ratios of the coefficients can now be found from (22), { while to completely determine the coefficients Xγ , Yγ the normalization condition is used,

∫

da

|u (q, z )|2 dz = 1.

(23)

−db

It can now be easily seen that the chosen form of the amplitudes indeed “automatically” satisfies the boundary conditions in Eq. (12), even when solutions of (22) are obtained in an approximation. On the other side, the error of an approximation is easily quantifiable by the mismatch of the amplitudes on the interface. 4. Identification of the modes Each natural frequency or phase velocity, solution Eqs. (22)–(23) for a given wavenumber q, belongs to one of the branches of the spectrum and determines the wave pattern specific to the branch. The branch is then easily identifiable if the pattern has some symmetry property which leaves the frequency spectrum invariant. In our case symmetry with respect to the middle plane is absent, the set of equations cannot be reduced to a simpler form and classification of excitation branches as symmetric (S) or antisymmetric (A) is not applicable. Indeed, the shapes of the amplitudes resulting

6

S. Cojocaru / Wave Motion 92 (2020) 102430

from numerical solution demonstrate an apparent lack of symmetry. For instance, it can be shown that the four equations describing the Stoneley modes localized near the interface correspond to the short wavelength limit of the Eqs. (22). It is known that there exist parametric restrictions for the real valued solutions of these equations to exist, e.g., [5]. However, for a finite value of q it is highly non-trivial to identify such a mode because its amplitude is spread over the hole sample. On the other side, when the symmetry is restored in the single-layer limit (“a = b”) one should recover the pair of Lamb equations

(

1 + v2

)2

4wv

( =

tanh (qw d/2) tanh (qv d/2)

)±1

,

(24)

were the upper sign is for the S −modes with an easily identifiable pattern ux (z ) = ux (−z ) , uz (z ) = −uz (−z ) ,

(25)

and the lower sign is for the A−modes, respectively ux (z ) = −ux (−z ) , uz (z ) = uz (−z ) .

(26)

If there exist some invariance property of the equations then it could simplify their solution and allow to classify the modes and clarify the structure of the spectrum, in particular, to establish a continuity relation with the solutions for a single layer. Such symmetry transformations for a composite plate “a ̸ = b” which generalize the above relations for a uniform plate are defined as follows: 1. For the D waves we have,

(

ux,a (q, z ) = ux,b (q, −z , a ⇄ b) , ux,b (q, z ) = ux,a (q, −z , a ⇄ b) ,

)

uz ,a (q, z ) = −uz ,b (q, −z , a ⇄ b) , uz ,b (q, z ) = −uz ,a (q, −z , a ⇄ b) .

(27)

2. For the F waves we have,

(

ux,a (q, z ) = −ux,b (q, −z , b ⇄ a) , ux,b (q, z ) = −ux,a (q, −z , a ⇄ b) , uz ,a (q, z ) = uz ,b (q, −z , b ⇄ a) ,

uz ,b (q, z ) = uz ,a (q, −z , a ⇄ b) .

) (28)

The meaning of the operations, e.g., in the first line is that for the D waves the value of the ux in the upper (lower) stratum for some given set of values {da , db , µa , µb , etc.} is the same as the value of ux in the lower (upper) stratum with inverted coordinate z −→ −z and interchanged (b ⇄ a) values for the respective material parameters, {db , da , µb , µa , etc.}. The notation “a ⇄ b” emphasizes the continuity with the single layer limit, “a = b”: D −→ S, F −→ A. Validity of these relations is verified by the Eqs. (10) and the boundary conditions. We will further use more compact notations:

γ −→ −γ (i.e., a −→ b, b −→ a) : µa −→ µb , µb −→ µa ; etc .

(29)

The transformation rules for the coefficients of the partial waves follow from (27) and (28) and ensure invariance of the equations in (22) and (23) (c (γ ) = c (−γ )): 1. D−waves, Xγ (q) = −X−γ (q) , Yγ (q) = Y−γ (q) .

(30)

2. F −waves, Xγ (q) = X−γ (q) , Yγ (q) = −Y−γ (q) .

(31)

Since the amplitudes are defined up to the choice of a common unitary prefactor or sign, care should be taken when numerical solution is considered. Indeed, by flipping the sign of the amplitude relation (30) is converted to (31) and vice versa. Flipping of the sign may occur when computations for a single point of the spectrum is carried out without taking into account the continuity between solutions belonging to the same branch. As it will be shown below, at long wavelengths one can make a clear distinction between the two types of modes so that correspondence to the single layer be established. When explicit expressions are known (e.g., analytic solutions at long wavelengths) the “sign problem” does not pose difficulties, similarly as for the Lamb’s solutions (which are the single layer limit of the present case). The sign degeneracy can also be lifted, e.g., by applying a weak coupling to an external field, change of boundary conditions or considering weak non-linearity. In this case, by analogy to Lamb case, it is to be expected that the two types of modes would respond in a qualitatively different way. Note also that according to (21) the coefficients X , Y determine the values of respective components of the amplitude on outer surfaces. When “a” and “b” media become equivalent we recover Lamb’s equations and the known symmetry relations. This can be verified directly by taking the “a = b” limit in the governing equations: {X = Xa = −Xb , Y = Ya = Yb } for the D −→ S waves and {X = Xa = Xb , Y = Ya = −Yb } for the F −→ A waves. Respectively, one gets X = Y tanh (v qd/2) for the symmetric Lamb modes and Y = X tanh (v qd/2) for the antisymmetric ones. The simplest way to see this is for a particular choice da = db = d/2 in (22) when the symmetric components of the amplitudes give trivial identities and expressions in (20) simplify to the familiar forms. Since at long wavelengths v takes finite values for the two fundamental modes A0 and S0 , the above exact relations set the relative scale for the wavenumber dependence of the unknown parameters X , Y and are useful in deriving the respective long wavelength solutions in the form of Taylor series. As will be shown below the same scale of relations remain qualitatively valid for the X γ and Yγ solutions of our Eqs. (22), corroborating the continuity conjecture F0 −→ A0 , D0 −→ S0 in the parametric space.

S. Cojocaru / Wave Motion 92 (2020) 102430

7

5. Long wavelength behavior Let us first consider the low energy sector of the spectrum, more specifically, the gapless acoustic branches of the normal modes (i.e., c (q −→ 0) ̸ = ∞). The form of the Taylor series for the long wavelength expansion in Eqs. (22) and (23) depends on the behavior of the parameters v and w in (11) which, in their turn, depend on the type of wave. Thus, for the D−wave consistency of the equations is provided by the following expressions cD = c0 + c2 q2 + · · · ,

(32)

and ∞ ∑

Xγ (q) =

x2n+1,γ q2n+1 ,

Yγ (q) =

n=0

∞ ∑

y2n,γ q2n .

(33)

n=0

We see that the scale of relations between X and Y is the same as for the single layer plate. By substituting (33) into the linear system of four equations (22) one can see that it is decoupled into matrices of the second rank that give the secular equations for the respective terms of the series for the phase velocity and the three ratios of the coefficients in (33). The value of the fourth coefficient is then fixed by the normalization condition (23) which provides the non-homogeneous closure equation. The procedure is iterated for the higher order terms in q so that on every iteration step the highest rank of the matrices for the coefficients in the Taylor expansion is 2. Below are the resulting expressions for the leading terms of the expansion. For the phase velocity

√ c0 = 2

da µa (1 − Ja ) + db µb (1 − Jb ) db ρb + da ρa

,

(34)

where Jγ = s2γ /ℓ2γ ,

(35)

The amplitudes y0,a =

c02 − 2s2a

√

dc02

c 2 − 2s2 , y0,b = 0√ 2 b ,

(36)

dc0

x1,a = −x1,b (a ⇄ b) ,

(37)

verify the symmetry relation (30). Expression for the velocity of the dilatational mode in (34) is well known [29,44,45]. It is also interesting to note that essentially the same expression [32,33,48] is obtained for the symmetric “sandwich” structure composed of two materials (one at the core and the other for the outer layers) while expressions describing the flexural wave do not fall in the same pattern, as will be seen from the result below. Unlike the D0 − mode, which shares some generic features with bulk acoustic waves, such as a linear frequency dispersion, the fundamental flexural mode F0 differs qualitatively and is more difficult to describe analytically even for a single layer plate. Nevertheless, following the scheme outlined above it can be shown that the following form of the Taylor series is consistent with the Eqs. (22) and (23): cF = c1 q + c3 q3 + · · · .

(38)

and Xγ =

∞ ∑

x2n,γ q2n , Yγ =

n=0

∞ ∑

y2n+1,γ q2n+1 ,

(39)

n=0

From the Eq. (22.3) we find x0,a x0,b

s2a

=

s2b

.

(40)

The two expansion coefficients, y1,a and y1,b , are then obtained from the Eqs. (22.1) and (22.2): y1,a =

x0,a 2

×

d2 µb (1 − Jb ) + d2a (µa (1 − Ja ) − µb (1 − Jb )) da µa (1 − Ja ) + db µb (1 − Jb )

,

(41)

and the relation y1,b = −y1,a (a ⇄ b) verifies the property (31). Substitution of the expressions (40) and (41) into Eq. (22.4) gives the required result for the leading term of the Taylor series in (39) for the phase velocity:

√ c1 =

)2

4µa µb da db (1 − Ja ) (1 − Jb ) d2 + µa d2a (1 − Ja ) − µb d2b (1 − Jb )

(

3 (da µa (1 − Ja ) + db µb (1 − Jb )) (ρa da + ρb db )

,

(42)

8

S. Cojocaru / Wave Motion 92 (2020) 102430

Fig. 2. Calculated amplitudes ux (z ) and σxz (z ) (u is normalized and σxz is in arbitrary units). Aluminum/Tungsten bilayer with dAl /dW = 1/4 for the next-to-fundamental flexural branch at qd = 1.87 and phase velocity c = 6408 m/s. The interface is at z = 0.

Finally, the normalization equation (23) leads to the following expressions x0,a = 2

s2a 2 2 c1 q

√ ,

x0,b = 2

d

s2b 2 2 c1 q

√ .

(43)

d

The singularity of the above two coefficients and vanishing of the phase velocity (11) in the limit q −→ 0 means that the leading term in the series for v and w in Eq. (11) is equal to 1. This ensures that the amplitudes in (20) remain finite and properly normalized: da

(∫ lim

q−→0

|ua (q, z )| dz + 2

0

∫

|ub (q, z )| dz 2

)

= 1.

−db

0

The higher order terms in (38) and (39) can be obtained by the iterative procedure as already mentioned. In contrast to the smooth behavior of the uniform plate, a distinctive feature of the layered system is that amplitudes have a kink on the interface, Fig. 2, which becomes more pronounced for shorter wavelengths. Parameters for the bulk materials, aluminum (a) and tungsten (b) [37,49], are as follows: ℓAl = 6419.23 m/s, sAl = 3118.67 m/s, ρa = 2702 kg/m3 , ℓW = 5214.41 m/s, sW = 2881.89 m/s, ρb = 19250 kg/m3 . The origin of the kink is in the interface boundary conditions discussed earlier (see, e.g., Eqs. (6), (9) and (22)) which imply that derivatives of the amplitudes are piecewise continuous functions of z and should have a discontinuity on the interface. The parametric space of the bilayered system is overwhelmingly multidimensional and revealing some general features by direct observation seems impossible. As follows from the analytic expressions obtained above, a better way to grasp the general behavior is to choose material “a” as a reference single-layer plate of the same total thickness d = da and consider the three scaled characteristics of the composite plate:

µ=

µb 1 − Jb db ρb × ,δ = ,ρ = . µa 1 − Ja da ρa

(44)

The reference layer corresponds to Lamb’s fundamental S0 and A0 modes. Then the respective phase velocities can be presented in the form containing only scaled variables:

(

cD

)2 =

cS0

(

cF

)2

cA0

1 + µδ 1 + ρδ

( =

,

(45)

1 − µδ 2

)2

+ 4µδ (1 + δ)2

(1 + δ)2 (1 + µδ) (1 + ρδ)

.

(46)

The D0 wave, Fig. 3, confirms the expectation that increasing the thickness of the added layer modifies the properties in a monotonous way: the velocity either increases or decreases with thickness db depending on the density and the elastic properties of the material “b”. The condition s2b

( 1−

s2b

ℓ2b

) =

which follows from

s2a

∂ ∂δ

(

(

1− cD cS

0

)2

s2a

ℓ2a

)

,

= 0, identifies the threshold µ = ρ (the thick horizontal line on the surface in the figure)

which separates the two regimes with increasing and decreasing velocity of the composite plate. On this line one has cD = cS0 irrespective of the ratio db /da . The pattern remains the same for any value of ρ away from the line: increasing the thickness of the “b” — layer produces a strictly monotonous variation of cD /cS0 .

S. Cojocaru / Wave Motion 92 (2020) 102430

9

Fig. 3. Phase velocity cD /cS0 , Eq. (45), of the lowest dilatational mode in a composite plate scaled with velocity of the lowest symmetric Lamb mode in a uniform plate of material a with thickness da = d. The scaled parameters are defined in (35) and (44) and mass density ratio is ρ = 3.

Fig. 4. Phase velocity of the lowest flexural mode of a composite plate scaled with the velocity of the respective antisymmetric Lamb mode cF /cA0 , Eq. (46), for the same value ρ = 3 as in Fig. 3.

Fig. 5. Cross section of the surface in Fig. 4 illustrates the non-monotonous change of the phase velocity on the layer thickness ratio for the flexural mode (elastic constants ratio is µ = 3).

In contrast, velocity of the F0 − wave ( cF /cA0 ) depends on the thickness of the added layer in a very non-monotonous way, Fig. 4 . Thus, in Fig. 5 one can see that for a choice of parameters above the line ρ = 1 it is first sharply increasing and then drops at db /da ≃ 0.1 to start growing again at larger values. Respectively, the frequency of the flexural mode would change in the same way when a coating layer “b” is applied. In Fig. 6 the non-monotonous dependence is illustrated for a choice of parameters below the line ρ = 1. 6. Thickness resonances The other class of solutions in the long wavelength part of the spectrum may be referred to as gapped modes. There exist a number of studies discussing the high-frequency long-wave vibrations in a more general parametric framework,

10

S. Cojocaru / Wave Motion 92 (2020) 102430

Fig. 6. Cross section of the flexural mode velocity surface for a different choice of mass and elastic constants ratios: ρ = 0.3 and µ = 0.2.

see, e.g., [8,50]. Here we present the solution for the bilayered isotropic model in the context of the proposed approach to make the connection with the finite wavelength behavior discussed in the next section. In particular, it can be rigorously proven using Eq. (21) and analytic expressions for the amplitudes in the previous Section that in a bilayered plate for any combination of material parameters there exist a pattern of the fundamental flexural mode F0 which is also a characteristic of Lamb’s fundamental mode A0 : the signs taken by the components of amplitude on the outer surfaces sign (ux (z = da ) × ux (z = −db )) = −1, sign (uz (z = da ) × uz (z = −db )) = +1. Similarly, it can be proven that for the S0 and D0 pair the pattern is: sign (ux (z = da ) × ux (z = −db )) = +1. In this Section it will be shown that a similar A − F and S − D correspondence can be established for the thickness plate resonances. The real values of the parameters v and w in (20) indicate that the gapless modes considered above belong to the class of surface acoustic waves. In contrast, the plate resonance or cutoff modes, e.g., [2,5], correspond to the gapped solutions in the q −→ 0 limit and are volume waves represented by harmonic functions. This follows from the definition in (11) and the diverging dependence ivγ , iwγ ∼

1 q

.

(47)

Since the amplitudes (20) should remain finite, there exist only two possibilities for the {X , Y } coefficients to be compatible with the divergency in (47): I . ux ∼ q, uz ∼ 1 :

Yγ ∼ Xγ ∼ q;

(48)

II . ux ∼ 1, uz ∼ q :

Yγ ∼ 1, Xγ ∼ q2 .

(49)

Relations (48) and (49) demonstrate the decoupling of the two polarizations in the q −→ 0 limit which can be inferred directly from the wave equations due to vanishing of the last terms in the l.h.s. of Eqs. (10). Thus, also in a composite plate the resonance vibrations are either of type I, “thickness-stretch”, or type II, “thickness-shear”, similar to a single layer case [2], Ch. 8. For these standing waves the secular equation is reduced to the product of second rank matrices. The non-vanishing components of displacement amplitudes in (20) for the “thickness-stretch” modes take the following form uz ,γ (z ) = Xγ vγ cos

((

z − dγ ω/ℓγ /2.

)

)

(50)

From the Eqs. (22.2) and (22.4) we then obtain Xa /Xb =

sa cos (ωdb /ℓb ) sb cos (ωda /ℓa )

,

(51)

ρa ℓa sin (ωda /ℓa ) cos (ωdb /ℓb ) + ρb ℓb sin (ωdb /ℓb ) cos (ωda /ℓa ) = 0.

(52)

Waves of the type II are polarized in the x− direction ux,γ (z ) = Yγ cos

((

z − dγ ω/sγ ,

)

)

(53)

and from the Eqs. (22.1) and (22.3) we find Ya /Yb =

cos (ωdb /sb ) cos (ωda /sa )

,

ρa sa sin (ωda /sa ) cos (ωdb /sb ) + ρb sb sin (ωdb /sb ) cos (ωda /sa ) = 0.

(54)

(55)

Eqs. (52) and (55) are well known and can be arrived at as a limiting case from different approaches even when the general equations differ from (22), e.g., Refs. [44,51–53]. This may occur not because of the difference in methods, but

S. Cojocaru / Wave Motion 92 (2020) 102430

11

Fig. 7. Solution of the Eqs. (52) and (55) for the scaled dimensionless variables. In addition to the trivial single layer limit, e.g., db = 0, the free composite plate has an equidistant type I or type II resonance spectrum when the special conditions are satisfied, ρa ℓa = ρb ℓb or ρa sa = ρb sb .

due to mistakes and typos. For instance, up to the choice of notations, our main set of Eqs. (22) differ by a factor of 2 of some terms in two of the Eqs. (34) in Ref. [44] and a factor δ1 in the last equation; respectively, our result in (42) differs from Eqs. (19), (20) in Ref. [44]. After normalization we obtain the final expressions for the amplitudes of the resonance modes. Thickness-stretch modes: uz ,γ (z ) = G (ℓ) cos ωd−γ /ℓ−γ cos ω z − dγ /ℓγ .

)

(56)

ux,γ (z ) = G (s) cos ωd−γ /s−γ cos ω z − dγ /sγ .

(57)

)

(

)

( (

Thickness-shear modes:

)

(

)

( (

)

Here the normalization constant G is a function of the respective bulk velocities V

{

[ G (V ) =

( 2

2dγ cos

ω−γ

)

sin 2ωγ

(

1+

2ωγ

)} + 2d−γ

{ ( ) }]−1/2 ) ( sin 2 ω −γ 1+ , cos2 ωγ 2ω−γ

where we have used the convention (29) and the notations V = s, ℓ;

ωγ = ωdγ /Vγ ,

ω−γ = ωd−γ /V−γ .

It can be verified that the standard solutions for the “a = b” limit are indeed reproduced by these expressions. For instance, the sin (ωz /ℓ) dependence of the uz (z ) amplitudes, Eq. 8.1.101 of Ref. [5], for the thickness-stretch modes follows from (56) for the solutions ωd/ℓ = π (2n + 1) of the (52) in this limit (Eq. 8.1.99) where da = db = d/2, while the cos (ωz /ℓ) dependence results from the solution ωd/ℓ = 2π n. These resonance conditions, when an integer number of halfwavelength matches the width of the plate, are generally not valid for a composite plate. However, as one can see from the Eqs. (56) and (57), the generic feature of the resonances is that their magnitudes reach maximum at the free surfaces. Note that for the special condition

ρa Va = ρb Vb , the Eqs. (52) and (55) reduce to

ω (da /sa + db /sb ) = π n.

(58)

In this case also the spectrum of composite plate becomes equidistant. Fig. 7 gives a general picture of the resonance spectrum of the bilayered plate in terms of scaled variables for the thickness-shear, V = s, and thickness-stretch, V = ℓ, types of waves. The separate sheets for each type of resonances correspond to the frequencies of the flexural and dilatational modes which follow in an alternating order. As can be seen, also the gaps between different branches change in a non-monotonous way with the variation of material parameters. Fig. 8 illustrates this behavior for the W /Al free standing plate.

12

S. Cojocaru / Wave Motion 92 (2020) 102430

Fig. 8. Scaled resonance frequencies for a b/a = Tungsten/Aluminum plate ρb sb /ρa sa = 6.58 (thick dashed) versus scaled thickness ratio dW /dAl as compared to the equidistant spectrum for materials satisfying the special condition (58) (thin continuous).

Fig. 9. Amplitudes of the first lowest resonance modes in a W /Al free plate with thickness ratio dW /dAl = 1/3.

Considering a 100 nm thick plate with the ratio db /da = 1/3 as an example, we see that the first five resonances are dominated by the shear modes: the shear F1 , stretch D1 , shear D2 , shear F2 and shear D3 modes with the corresponding frequencies ω/2π ≃ 11.7 GHz (a), 24 GHz (b), 31 GHz (c), 49.8 GHz (d), 59 GHz respectively. A discussion of a threelayered structure which includes a situation when the lowest shear thickness resonance frequency approaches zero is presented in Ref. [33]. The non-zero components of the amplitudes are shown in Fig. 9. Comparison of the amplitude in Fig. 9 with Fig. 2 illustrates its evolution with the wavenumber q. Note that F and D solutions do not follow in a regular order with increasing resonance frequency, however, they do follow in alternating order within the same type, stretch or shear. Let us now consider the behavior of amplitudes (56) and (57) on the outer surfaces, which is determined by the coefficients Xγ and Yγ (see also Eqs. (21)). The ratios (51) and (54) of these coefficients are given by the solutions for the resonance frequencies of the corresponding equation, (52) or (55). The l.h.s. of each equation is a continuous function of ω and can be rewritten in the form where the ratio appears explicitly. It can then be easily proven that the sign of the ratio changes from one solution for the frequency to the next. In the considered example of Al/W plate the lowest resonance (shear) reproduces the sign pattern of an A mode at q = 0 (sign (ux (z = da ) × ux (z = −db )) = −1, sign (uz (z = da ) × uz (z = −db )) = +1) and is therefore assigned to F1 . The next shear resonance reproduces the sign pattern of an S mode, etc. Similarly, the stretch resonances start with an S wave pattern which corresponds to the D1 mode, etc. Thus, “sign rule” reveals a characteristic property of the single layered plate which is preserved despite of breaking the symmetry of the structure by adding a layer and which allows to identify the parametric continuity of the two symmetry classes of Lamb waves (F −→ A and D −→ S). Its validity is, however, limited to long wavelengths since the physical interpretation is that the composite system may indeed be considered

S. Cojocaru / Wave Motion 92 (2020) 102430

13

Fig. 10. Phase velocity versus dimensionless wavenumber qd for a few lowest branches of the spectrum in W /Al plate, thickness ratio is dW /dAl = 1/4. Regions of “avoided crossing” A–C are discussed in the text.

as roughly uniform in this case. As follows from numerical calculations, deviations from the rule may appear at larger wavenumbers, however for some branches of the spectrum the sign pattern remains unchanged. In the next Section the above results are used for the analysis of finite wavelengths spectrum. 7. Finite wavelength behavior The analytical expressions have been verified to be in agreement with numerical approaches, such as Disperse software [41,54], and experimental data on, e.g., nickel/bronze coated plates [55]. In this section we consider the lowest branches of the phase velocity spectrum in Al/W plate by numerically solving the Eqs. (22)–(23) in order to analyze the behavior of the amplitudes. Since the wave pattern varies continuously along each branch, it is possible to connect the solutions to the ones described in the previous Sections. It can be shown that in the short wavelength limit, when considering normal modes localized at the respective surfaces, the equations reduce to those of the Rayleigh surface waves (R) or Stoneley interface guided wave (St) [37,42]. Thus, for the aluminum/tungsten plate we obtain the following phase velocities cR (W ) = 2664 m/s, cR (Al) = 2914 m/s and cSt (Al/W ) = 2787 m/s (parameters vW and vAl take real values for the respective materials) in agreement with known results for the half-space boundary conditions [49]. To obtain the volume waves it is necessary that one of the vγ , wγ vanishes in this limit while taking imaginary values at large but finite wavenumbers qd. Note that in the considered example we have sAl > cR (Al) > sW so that vW takes imaginary values close to the Rayleigh wave limit in aluminum. In this case for the solution converging to the Rayleigh wave in aluminum the amplitude in tungsten is exponentially small. In Fig. 10 a few lowest branches are shown for the choice of layer thicknesses dW /dAl = 1/4. The slowest of the two Rayleigh modes (W ) originates at small wavenumbers from the fundamental flexural mode and belongs to the F0 branch. Similarly, the Stoneley mode (both vW and vAl take real values when qd ≫ 1) originates from the fundamental dilatational wave and belongs to the branch D0 . The fastest of the Rayleigh modes is localized near the outer surface of the Al layer and seems to belong to the branch F1 which originates from the lowest thickness shear resonance [52], this Rayleigh mode is also known as Sezawa wave, e.g., [2,5]. The fourth branch, D1 , approaches the bulk shear wave in Al at short wavelengths (both vW and vAl take imaginary values) and originates from the lowest thickness stretch resonance discussed in the previous section. The “terracing” behavior of the dispersion curves c (q) with the regions of “avoided crossing”, (A, B, C ) in Fig. 10, is well known [44,56]. Thus, it can be seen that c (q) is almost constant for the two branches at q > qB and q > qC and respective velocities are close to sAl and cR (Al). It thus appears that we have found the connection between long and short wavelength properties of the spectrum. Let us, however, analyze the behavior of amplitudes in more detail. Fig. 11 gives a closer view on dispersion in the C region. From Fig. 12 it can be seen that, aside from an “exchange of slopes” between the incident (I and III) and outgoing (II and IV) lines in Fig. 11, there is also an “exchange” of the wave profiles between the F1 and D0 branches. The exchange is obviously not exact, in particular, one can see the difference in sign of the amplitudes on the outer surfaces, which, incidentally, remains a valid characteristic of the D and F classification scheme even at such short wavelengths. As a consequence, the amplitudes undergo a large but continuous shape modification on a relatively narrow interval of q along each of the branches, so that the typical pattern of the Rayleigh wave localized at the Al surface for the D0

14

S. Cojocaru / Wave Motion 92 (2020) 102430

Fig. 11. Phase velocity dispersion in the region C Fig. 10 near the turning point qd = 36.5. Amplitudes for the respective segments are shown in Fig. 12.

Fig. 12. The four graphs show the component ux (q, z) of the amplitude for the upper, D1 , and lower, F1 , curves in Fig. 11 in the immediate vicinity of the turning point. Note the signs of the amplitude on the outer surfaces for the branches I–II (D1) and III–IV (F 1). The uz component (not shown) follows the same sign pattern as described in the text.

branch at qd < 36.5 is “transferred” to the F1 branch at qd > 36.5. The shape taken by the amplitude on the D0 branch becomes more localized near the interface and at larger values of qd finally converges to the typical profile of Stoneley wave while cD0 −→ cSt . However, by comparing the above results to the value of the lowest bulk wave velocity, sW , it can be concluded that the picture in Fig. 10 is incomplete because no branch appears to converge to this value. This “discrepancy” is resolved due to the “shape swapping” phenomenon taking place in the other regions of “avoided crossing” which are separated by a rather extended interval on the qd axis. Thus, in the region B the amplitude of the F1 branch receives the localized shape, inset I in Fig. 12, by transfer from D1 , which in turn continues as a volume wave of mixed polarization till about qd = 65. From this region the D1 branch emerges as a different volume wave of mixed polarization, as shown in Figs. 13 and 14. Its amplitude is of harmonic character in tungsten and is exponentially suppressed in aluminum. At an even larger value at about qd = 141.6 close to the region c (q) ≃ cR (Al) the two branches, D1 and F1 , “nearly collide” again to exchange their profiles for the last time. The outgoing F1 branch disperses down in velocity and converges to the true shear bulk wave (ux −→ 0) in tungsten with cF1 −→ sW , while the D1 branch takes over in converging to the Rayleigh surface wave in aluminum. The above analysis allows to establish a connection between long and short wavelengths modes (bulk and surface guided) as belonging to the same dispersion branch of the vibration spectrum, the branches being resolved by their long wave origins described earlier. The correspondence for the first few branches in the considered example of an

S. Cojocaru / Wave Motion 92 (2020) 102430

15

Fig. 13. Detailed view of the downturn of the D1 branch, see Fig. 10, in the region near qd = 65.

Fig. 14. The ux (q, z) component of the amplitude for the two branches in Fig. 13 close to qd = 65.

aluminum/tungsten plate is as follows: F0 (lowest velocity Rayleigh wave in the softer material, W ), D0 (Stoneley interface wave), F1 (shear bulk wave in W ), D1 (Rayleigh wave in the stiffer material, Al), D2 (shear bulk wave in Al). One can conclude that for the concidered few branches the correspondence follows the “natural ordering” rule: the lower velocity mode at long waves evolves along the dispersion curve into the lower velocity mode at short wavelengths. 8. Summary The proposed approach addresses the problem of wave propagation in layered structures from a fresh perspective when one first finds a partial wave basis in such a form which satisfies the boundary conditions on the outer surfaces exactly due to the choice of functional form and not by the values of the unknown parameters. The exact match of boundary conditions is achieved by choosing the arguments of the partial waves in the form which defines the respective surfaces. The remaining interface boundary conditions are then solved in terms of the given basis. Efficiency of the method has been demonstrated by application to the description of a bilayered plate. Known analytic results for the long wavelength part of the spectrum are recovered and new ones are obtained within the proposed framework. This has allowed to find a simple and exact characteristic of the wave pattern unifying the solutions for single and bilayered plates in terms of the signs taken by the amplitudes on the outer surfaces. It has further allowed to clarify the connection between long wave modes and their short wavelength counterparts as belonging to the respective branches of the spectrum. It has been shown that general features of the dependence on material parameters can be represented in terms of scaled quantities.

16

S. Cojocaru / Wave Motion 92 (2020) 102430

The analysis has revealed a non-monotonous dependence on parameters of the fundamental flexural mode which would show up in other physical properties, such as thermal conductivity, electron–phonon interaction etc. The approach can be extended to multilayered systems, where powerful matrix methods using the plane wave basis have been developed in recent years. Another possibility is to consider curved plates or buckled composite configurations where the surface adapted partial waves contain the surface equation in the argument of the respective basis function. It would also be interesting to consider the validity of the “sign rule” correspondence with Lamb solutions for other traction-free layered or curved geometry systems. Acknowledgment This work has been financially supported by MCI Romania, Project No. PN 19060101/2019 Appendix The explicit form of the main equations follows the numbering in Eq. (22).

(( −Xa

1 + va2

sinh (da qwa ) − sinh (da qva )

2wa va

( +Ya

)

2 1 + va2

(

) cosh (da qwa ) − cosh (da qva )

((

1 + vb2

= Xb + Yb

Xa

1 + vb2

1 + va2

2wa 1 + va2

(

cosh (da qwa ) − cosh (da qva )

) sinh (da qwa ) −

1 + vb2

vb

2wb 1 + vb2

(

1 + va2

(

va (

−µa Ya

4wa

1 + vb2

+ µb Yb

((

sinh (da qva )

) cosh (db qwb ) − cosh (db qvb )

) sinh (db qwb ) −

1

vb

) sinh (db qvb ) .

(A.2)

(cosh (da qwa ) − cosh (da qva )) 1 + va2

(

1 + va2

= Xb µb

va

)

)

( (

1

)

2

+ Yb

(A.1)

)

((

1

−µa Xa

) cosh (db qwb ) − cosh (db qvb ) .

)

2

(

Xa µ a

)

(

va (

= Xb

sinh (db qwb ) − sinh (db qvb )

2

((

−Ya

)

)

2wb vb

(

1

)

)

) sinh (da qwa ) −

va

)

)

sinh (da qva )

)

(cosh (db qwb ) − cosh (db qvb )) ) ( ) 1 + vb2 4wb ( ) sinh (db qwb ) − sinh (db qvb ) . vb 1 + vb2

v (b

1 + va2

)2

2wa va

) sinh (da qwa ) − 2 sinh (da qva )

+2µa Ya (cosh (da qwa ) − cosh (da qva ))

(A.3)

S. Cojocaru / Wave Motion 92 (2020) 102430

(( = µb Xb

1 + vb2

)2

2wb vb

17

) sinh (db qwb ) − 2 sinh (db qvb )

+ 2µb Yb (cosh (db qwb ) − cosh (db qvb )) .

(A.4)

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56]

B.A. Auld, Acoustic Fields and Waves in Solids, second ed., Krieger, Malabar, FL, 1990. K.F. Graff, Wave Motion in Elastic Solids, Dover Publications, NY, 1991. L.M. Brekhovskikh, O.A. Godin, Acoustics of Layered Media I, Springer, New York, 1990. N.A. Haskell, The dispersion of surface waves on multilayered media, Bull. Seismol. Soc. Am. 43 (1953) 17. J.D.N. Cheeke, Fundamentals and Applications of Ultrasonic Waves, CRC Press LLC, USA, 2012. J.L. Rose, Ultrasonic Guided Waves in Solid Media, Cambridge University Press, NY, 2014. K.C. Le, Vibrations of Shells and Rods, Springer, Berlin, 2013. J.D. Kaplunov, Dynamics of Thin Walled Elastic Bodies, Academic Press, USA, 2012. A.A. Kaufman, A.L. Levshin, Acoustic and Elastic Wave Fields in Geophysics, Elsevier, Amsterdam, 2005. M.A. Stroscio, M. Dutta, Phonons in Nanostructures, CUP, United Kingdom, 2004. A.N. Cleland, Foundations of Nanomechanics, Springer-Verlag Berlin Heidelberg, 2003. Pham Chi Vinh, A. Aoudi, Vu Thi Ngoc Anh, Wave Motion 83 (2018) 202. M.J.A. Schuetz, E.M. Kessler, G. Giedke, L.M.K. Vandersypen, M.D. Lukin, J.I. Cirac, Phys. Rev. X 5 (2015) 031031. M.V. Gustafsson, P.V. Santos, G. Johansson, P. Delsing, Nat. Phys. 8 (2012) 338–343. A. Balandin, D. Nika, Phononics in low-dimensional materials, Mater. Today 15 (2012) 266. Xiyue An, Hualin Fan, Chuanzeng Zhang, Wave Motion 80 (2018) 69. H.Q. Nguyen, M. Meschke, H. Courtois, J.P. Pekola, Phys. Rev. Appl. 2 (2014) 054001. K.L. Viisanen, J.P. Pekola, Phys. Rev. B 97 (2018) 115422. J.T. Muhonen, M. Meschke, J.P. Pekola, Rep. Progr. Phys. 75 (2012) 046501. I. Maasilta, A.J. Minnich, Phys. Today 67 (2014) 27. Jin-Wu Jiang, Bing-Shen Wang, Jian-Sheng Wang, H.S. Park, J. Phys.: Condens. Matter 27 (2015) 083001. D.L. Nika, A.A. Balandin, Two-dimensional phonon transport in graphene, J. Phys.: Condens. Matter 24 (2012) 233203. D.-V. Anghel, S. Cojocaru, Eur. Phys. J. B 90 (2017) 260. M. Weiler, H. Huebl, F.S. Goerg, F.D. Czeschka, R. Gross, S.T.B. Goennenwein, Phys. Rev. Lett. 108 (2012) 176601. W.J. Parnell, I.D. Abrahams, Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 464 (2008) 1461. C. Schaal, A. Mal, Wave Motion 82 (2016) 86. K. Kanda, T. Sugiura, Wave Motion 66 (2018) 177. D.L. Nika, Phonon Engineering in Graphene and Semiconductor Nanostructures, USM, Chisinau, 2015. G.A. Rogerson, K.J. Sandiford, Int. J. Solids Struct. 37 (2000) 2059. E.P. Pokatilov, D.L. Nika, A.A. Balandin, Phys. Rev. B 72 (2005) 113311. S. Uno, J. Hattori, K. Nakazato, N. Mori, J. Comput. Electron. 10 (2011) 104–120. M. Yu. Ryazantseva, F.K. Antonov, Internat. J. Engrg. Sci. 59 (2012) 184. J.Kaplunov D.A. Prikazchikov, L.A. Prikazchikova, Int. J. Solids Struct. 113–114 (2017) 169. S. Andrieu, et al., Phys. Rev. Mater. 2 (2018) 064410. V. Panella, G. Carlotti, G. Socino, L. GioVannini, M. Eddrief, K. Amimer, C. Sebenne, J. Phys.: Condens. Matter 9 (1997) 5575. W.M. Ewing, Elastic Waves in Layered Media, McGraw-Hill, NY, 1957. G.W. Farnell, E.L. Adler, Elastic Wave Propagation in Thin Layers, Physical Acoustics, Vol. IX, Academic Press, New York, 1972, chap. 2. M. Destrade, Y.B. Fu, Internat. J. Engrg. Sci. 44 (2006) 26–36. S.K. Datta, A.H. Shah, Elastic Waves in Composite Media and Structures, CRC Press, USA, 2009. H. Ezawa, Ann. Phys., NY 67 (1971) 438. M.J.S. Lowe, Matrix techniques for modeling ultrasonic waves in multilayered media, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 42 (1995) 525. Haidang Phan, Tinh Quoc Bui, Hoai T.-L. Nguyen, Chi Vinh Pham, Wave Motion 79 (2018) 10. W.Q. Chen, H.M. Wang, R.H. Bao, Compos. Struct. 81 (2007) 233. J.P. Jones, J. Appl. Mech. 31 (2) (1964) 213–222. Wang Xiao-Min, Lian Guo-Xuan, Li Ming-Xuan, Chin. Phys. Lett. 20 (2003) 1084. S.I. Rokhlin, L. Wang, J. Acoust. Soc. Am. 112 (2002) 822. E.L. Tan, J. Acoust. Soc. Am. 119 (2006) 45. P. Lee, N. Chang, J. Elasticity 9 (1979) 51. Y. Chevalier, M. Louzar, G.A. Maugin, J. Acoust. Soc. Am. 90 (1991) 3218. L.Y. Kossovitch, R.R. Moukhomodiarov, G.A. Rogerson, Acta Mech. 153 (2002) 89. R. Jones, E.N. Thrower, J. Sound Vib. 2 (1965) 167. R. Jones, E.N. Thrower, J. Sound Vib. 2 (1965) 328. P.C.Yang C. H. Norris, Y. Stavsky, Int. J. Solids Struct. 2 (1966) 665. Zheng Fan, M.J.S. Lowe, Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci. 465 (2009) 2053. Yung-Chun Lee, Sheng-Wen Cheng, IEEE Trans. Ultrason. Ferroelectr. Freq. Control 48 (2001) 830. R.D. Mindlin, Structural Mechanics, in: Proc. 1st Symp. on Naval Struct. Mech., Pergamon Press, New York, 1960, pp. 199–232.