Evaluation of complex spectrum of the symmetrical Lamb waves for three-layered plates at low frequency. Part I: Foundations

Composite Structures 230 (2019) 111429

Contents lists available at ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

Evaluation of complex spectrum of the symmetrical Lamb waves for threelayered plates at low frequency. Part I: Foundations

T

D.D. Zakharova, , A.V. Nikonovb ⁎

a b

Department of Mathematical Analysis, Moscow State University of Railway Engineering, Moscow, Russia Faculty of Mechanical Engineering, University of Ljubljana, Ljubljana, Slovenia

ARTICLE INFO

ABSTRACT

Keywords: Asymptotic expansion Complex spectrum Layered plates Dispersion equations Lamb waves

A new method combining the asymptotic and iteration approaches to study the total wavenumber spectrum for a plate made of isotropic layers is suggested. In this part we use the propagator matrices to deduce the dispersion equation and its static limit in a closed form with the use of symbolic calculations. The roots in statics which are analyzed at the first step play the key role. Then, the low frequency approximations for the dispersion curves are derived. The method is aimed to evaluate a dispersion curve of any order n and to describe it in a closed form by + which causes difficulty for the pure numerical an asymptotic formula. Especially it concerns the case n approaches.

1. Introduction The implementation of eigenwaves for solving different problems for structural members becomes widely used and embraces more and more complex composite structures. Especially it concerns the structural health monitoring and non-destructive evaluation of solids which can be considered as elastic and viscoelastic waveguides. As known, at each given frequency the elastic waveguide with a finite size of its cross section may have a finite number of real wavenumbers and an infinite set of complex wavenumbers. In the case of viscoelasticity all the wavenumbers are complex. However, the problem of finding the total complex spectrum is now resolved for just a few situations, representing such simplest isotropic and homogeneous waveguides as a plate or a circular cylindrical solid. Our paper is mostly focused on Lamb waves in layered plates. In the classical papers [1,2] the appropriate tools have been introduced first, and then their results exhibited the succession with respect to earlier results for statics [3,4]. As far as the complex dispersion curves are concerned, the larger is the index n of the curve, the wider is the frequency interval where this curve is almost constant and the wavenumber k coincides with its static value [5–7]. The static values + (or can be efficiently approximated by their asymptotics as n |k| + ) and then calculated exactly for any index n . The complex wavenumbers of a circular rod behave similarly [8,9,4]. For other waveguides which may possess a layered configuration and/or anisotropy of materials, similar results are not obtained yet.

⁎

Currently, in order to obtain the real, imaginary and some complex branches of dispersion curves for anisotropic and/or layered plates in a reasonable frequency range, the numerical methods are used. Their common deficiency is the principal inability to describe the total in. For many engineering problems the finite spectrum as Rek, Imk information about some fist wavenumbers is sufficient, but not for all problems, not to mention the lack of the fundamental knowledge in this area. For example, the stress field expansion into mode series using the complete complex spectrum remains a natural representation in context of the wave interaction with the edge, the wave diffraction by an inhomogeneity (see, e.g., [10–16]) or when applying the non-reflecting boundary condition to reduce the size of the area considered [17]. In addition, the complete spectrum permits us the efficient use of the orthogonality relations of modes for both 2D waves [18–22] and 3D waves, introduced in [23–25] and partially repeated in a recent paper [26]. It should be noted that 3D waves in an infinite plate can be described using spectra for 2D waves [23–25] which extends the applicability of the respective results. The needs of the structural health monitoring have been reported in [27]: “As the research into analytical modeling of the interaction between GWs and damage in composite materials progresses, the importance of understanding the behavior of imaginary and complex wavenumbers in composite materials grows. However, the study of the imaginary and complex wavenumbers for GWs in composite materials is not as well documented in the literature as it is for isotropic metallic

Corresponding author. E-mail address: [email protected] (D.D. Zakharov).

https://doi.org/10.1016/j.compstruct.2019.111429 Received 21 February 2017; Received in revised form 1 July 2019; Accepted 11 September 2019 Available online 13 September 2019 0263-8223/ © 2019 Elsevier Ltd. All rights reserved.

Composite Structures 230 (2019) 111429

D.D. Zakharov and A.V. Nikonov

plates” (p.3). Speaking about the existing methods, we may refer to the combinations of the direct numerical approach with variational techniques (see the detailed review in [28,29]), and the ultimate use of the modal expansion and collocation method [30,31]. For multilayers, the most appropriate analytical approach provided by Thompson and Haskel [32,33], is based on the direct and recurrent formulas for the propagator matrices relating the magnitudes of waves in the neighboring layers and permitting one to obtain the reflection or transmission coefficients and dispersion relations. The alternative – and numerically more stable – method of the determinant decomposition is introduced and then modified for the inelastic media by Knopoff and Schwab [34,35]. The methods of exact and approximate analysis of the propagator matrices developed by Lowe [36] are now often used in combination with the direct discretization of the waveguide across the thickness (e.g., SAFE method, see reviews in [37,38], or scaled boundary finite element method [39,40]). A generalization of the spectral collocation method used in the DISPERSE [36] for the case of viscoelastic and anisotropic layered waveguides can be found in [41–43]. A stable recursive algorithm for multilayered media (recursive asymptotic stiffness matrix) with some generalization to the case of piezoelectric and functionally graded materials has been reported in [44–47]. An original method avoiding the pre-history of finding wavenumbers in the limited frequency range is presented in [48]. This method projects the system of differential equations governing the Lamb waves on a special basis of orthogonal functions. In contrast to the standard approach, it does not solve the transcendental dispersion equation but calculates the solutions of a classical eigenvalue problem for a “projection” matrix resulting from the previous step. The method allows approximation of the complex branches of dispersion curves and its accuracy and matrix dimension depend on the number of complex curves. The larger number of curves, the more basis functions are retained. The results on the symmetrical modes are presented for one isotropic layer. The generalization of the “classical” analytical and asymptotical tool for a linearly viscoelastic isotropic layer is performed in [49]. Our paper has four goals: (a) to suggest a complementary analytical tool to existing methods with the aim to obtain the total spectrum of the wavenumbers of solids, made of isotropic layers, and to clarify its regular behavior; (b) to demonstrate how this tool works with a relatively simple model problem; (c) to evaluate the dispersion curve of any order n and to describe it by a relatively simple asymptotic formula (especially as n ) with an estimation of the error and the interval of validity; (d) to reveal the most important features of this method for its potential application in the case of other structural members. In this Part I, we present a new method allowing us to achieve the goals (a) and (b) for a layered plate in the low frequency range. As a model problem we consider the propagation of the symmetrical modes in a three-layered plate. In general, the symmetrical layup is not a limitation and is taken as an example. The case of bending was considered in [50] but the symmetrical waves were not addressed. The limit behavior of the desired wavenumbers in statics plays key role. To this end, we obtain the explicit form of the dispersion equation and its static limit using propagator matrices. After that, the asymptotic behavior of the wavenumbers is considered in statics, and then, in the low frequency dynamics. In the Part II, we pursue our goals (c), (d) and show how the total complex spectrum is organized with the suggestion how to improve the iterative method of calculations. Note, that we don’t consider the interaction between the dispersion curves which happens at larger frequency but consider the initial segments of the dispersion curves. The length of this “initial” segment is increased with the curve order as we see in our analysis. So, the advantage of this result is that the larger is the modulus |k| the less correction in the asymptotic solution of the

dispersion equation is required. At the end, we discuss the main steps of the suggested approach in context of such perspectives as taking into account many layers, possible nonsymmetry of the stacking sequence and the geometry of structural members. 2. Mathematical statement of the problem Consider a multilayer plate, in which the jth layer of thickness < x, y hj = z j + 1 z j occupies a region < + , z j z zj + 1, (j = 1, 2, …, N ) , where x , y , z are Cartesian coordinates, and the z axis is directed across the layer. In an evident case the layer number j is omitted. The Lame constants and the mass density of the layer are denoted by j , µj , and j , respectively. We assume a time harmonic dependence exp( i t ), where ω is an angular frequency. This term is also omitted in the evident case in what follows. Let us introduce in each layer the displacement vector u = (u x , u y , uz ) and the symmetrical stress tensor, satisfying Hooke's law

=

xx

xy

xz

yx

yy

yz

zx

zy

zz

= divu + 2µ x u x = µ ( y u x + x u y ) , (x xz = µ ( z u x + x u z )

xx

,

xy

y

z) (1)

The desired field satisfies the equation of motion

div

+

2u

(2)

=0

and the conditions of the displacement and stress continuity on each interface

[

xz ]

=[

yz ]

=[

zz ]

(3)

= 0, [u x ] = [u y] = [u z ] = 0

The upper and lower faces are supposed to be stress free xz

=

yz

=

zz

(4)

=0

We study the spectrum of wavenumbers {k } for waves propagating as exp(ikx i t ) and satisfying the equation of motion (2) and boundary conditions (3) and (4), i.e., the Lamb waves in a multilayer plate [51]. Our main attention is focused on the infinite set of complex dispersion curves k ( ) with Rek ( ), Imk ( ) 0 . 3. Matrix formulas In this section we derive the explicit form of the propagation matrix for an elastic layer of thickness h . Let the coordinates of the layer faces be denoted as z+ = z2 and z = z1, and the thickness of the layer as h = z2 z1. For the time harmonic process, introduce the 6 × 1 motion-stress vector [52] = ( ik 1ux , k 1u z , k 2 zz , ik 2 zx , ik 1u y , ik 2 zy )T , containing the scaled displacements and stress components, and the following diagonal matrix

= diag[ EP EP+ ES ES+ ES ES+ ] where EP = e

P kh ,

S kh ,

ES = e

P, S

kP2, S k 2)1/2,

= (1

= 2k 2kS 2 = 2c 2cS2,

1 1 = are dimensionless quantities, c = k 1 is a phase speed and kP , S = cP,1S are wavenumbers, corresponding to the speed of bulk compressional and shear waves cP = [( + 2µ)/ ]1/2 and cS = (µ / )1/2 , respectively. Performing additional simplifications, such as,

c 2 = 2k 2kS 2 = 2µ , =

2k 2

=2

= 2µ,

1 (2

k 2 (2µkS2) 1

=

1

kS2 k 2) (4µ )

c2 = ( =

1

1

1) c 2 = 2µ (2

c 2) 1

=

µkS2 k k2 (2

2 2) 1

1

we arrive at the explicit form of the propagator matrix P (h) with 2

Composite Structures 230 (2019) 111429

D.D. Zakharov and A.V. Nikonov

the final basic relation on the layer faces

P (h ) = M ×

where

1

×M

(z2) = P (h) × (z1) with the following explicit expressions for M and M

(5)

E+ = diag[0 FP+ 0 FS+ 0 FS+], E = diag[ FP 0 FS 0 FS 0], e

(6)

Thus, the propagator matrix P (h) in the formula (6) can be represented in the following form:

= e kh, FP , S = e

1

( P, S 1) kh

P = eP+ + e 1P ,

P =M ×E ×M

1

(8)

where 1 ( 2

+ P FP

µ(

FP+

(

+ P FP

µ

1 2

FS+)

1

(

1 2

P+ =

FP+

)

+ 1 FS S

µ

)

+ 1 FS S

+ 1 FP

+

P

1 ( 2

FS+)

1

(

+ 1 FP

(

+ 1 FP P

µ(

)

FS+

S

+

FS+)

+

S

4µ

)

FS+)

P

(

4µ

)

FS+

FP+ +

S

P

2

)

FS+

FP+ +

S

(

FP+

( 2

(FP+

)

FS+

0

0

0

0

FS+

0

0

FS+)

0

0

S

P

4µ

FP+ + FS+)

(

2

2

( FP+ + FS+)

(

FS+

FP+ +

1

4µ

FS+)

FP+

)

S

P

(FP+

0

0

0

0

1 + F 2 S

0

0

0

0

1 µ F+ 2 S S

FS+ 2µ S 1 + F 2 S

and 1 ( 2 1 2

P = µ

FP

1

(

P FP

µ(

FP

(

P FP

1

M

1

=

FS ) S

2µ 2µ 0 0 2

1

2

1

1 (2 S ) 1 (2 S )

1 FP

FP +

1

(

µ(

FS

S

P

4µ

FS )

1 FP

FS

S

P

FP +

1

)

(

4µ

)

2

FS )

2

( FP + FS ) P

(

(

)

FS

FP

4µ

S

FP + FS ) P

FS

FP

S

FP

( 2

+

P

4µ 2

)

(

(FP

S

FS

)

FS )

FP

+

P

(FP

S

FS

FS )

)

0

0

0

0

0

0

0

0

0

0

0

0

1 F 2 S

0

0

0

0

1 µ F 2 S S

P

P

1 ( 2

µ

)

1 FS

1 P

M =

S

1

+

)

1 FS

+

(

1 2

FS )

2µ 2µ P 0 0

S

1 2µ 2µ 0 0 1 (2 P ) 1 (2 P ) 2 1

1 1

2

1

0 0 0 0

S

S

1 1

1 2µ S 2µ 0 1 0 µS (4µ)

1

(4µ)

1

(4µ S )

1

(9)

(z n + 1) = P × (z1), P = P (hN ) × P (hN 1). ..×P (h1)

+ ej P j , where the propagator matrix for j th layer is P (hj ) = kh hj = zj + 1 z j , ej = e j and the propagator P for any number of layers can be presented in the following convenient form with the exponential terms collected ej P+j

1 µS 1

0

(4µ P )

0

0

(4µ)

1

0

0

(4µ)

1

0

0

0

0

0

2

1

0

0

0

0

2

1

(2µ S ) (2µ S )

{

e k (sN hN +…+ s1 h1 ) PNsN (hN ) × ...×P1s1 (h1), sN = +

P =

0

1

0

s1, … , sN

1

(10)

The matrix products in (10) can be easily analyzed step by step. In particular, for a two-layered plate it yields

1

P (h2 ) P (h1) = e2 e1 A + e2 1 e1 B + e2 e1 1 C + e2 1 e1 1 D

1

(11)

Therefore, the properties of the subsequent dispersion equations are determined by the properties of the matrices

The obtained matrices correspond to those presented in [52,53]. The difference appears due to the scaling introduced for the components of vector . The order of these components corresponds to the Pand SV-polarized waves (first 4 positions), and to SH-waves (last 2 positions), and it results in the above 4 × 4 and 2 × 2 nonzero submatrices of M . Some other consequences are considered in Sections 4 and 9. For further analysis we need more convenient form of matrices and to this end we single out the exponential term exp(kh) . Our representations yield the following explicit formula for the matrix (h)

(h) = eE+ + e 1E

1 F 2 S

For a multilayer we consequently obtain

0 0 0 0

(4µ P ) 1

(4µ S )

FS 2µ S

A = P+2 × P+1, B = P 2 × P+1, C = P+2 × P 1, D = P 2 × P 1 4. Asymptotic formulas for the components of the propagator matrix in case of the in-plane and out-of-plane problems and the static limits Let us now express the leading terms in the asymptotic expansion of + . The coefthe submatrices of P and P+ for a single layer as |k| ficients of such expansion are matrices which are composed of the

(7) 3

Composite Structures 230 (2019) 111429

D.D. Zakharov and A.V. Nikonov

respective coefficients for each element of a submatrix. Consider the submatrices with rows and columns 1,2,3,4 describing the most interesting wave motion with P- and SV-polarization. Here and in what follows in the submatrices the upper figures mean rows and the lower figures mean columns. We obtain

( )

1234 P = k1P (1) + P (0) + k 1P ( 1234

P

1)

2P ( 2)

+k

P

cS2 cP2

h 1 2

J

1 J 2µ

2µJ

J

kS2 h2 1 8 (1

cS4 cP4

1 J 2µ

2µJ

J

1 (1 2µ

I)

2µ (1

J

I)

I)

+…

(1

I)

1 (1 2µ

I

1 2

+

2µ (1

I)

1

1

1

1

,I= 1 0 ,1 0 1

Note, that P (

s)

2 ),

= O(

I)

I

0

+

cS2 2cP2

+

1 2

(1

cS2 cP2

1 (1 2µ

I)

2µ (1

I)

(1

2µ (1

I)

0 yields

J

1 J 2µ

2µJ

J

I)

I)

I

.

For the matrix P+ we similarly obtain

P+(1) =

cS2 cP2

h 1 2

JT

1 T J 2µ

2µJ T

JT

, JT = 1 1

zz

1 1

zy

=

(1

s)

I)

= O(

JT

1 T J 2µ

2µJ T

JT

1 (1 2µ

I)

2µ (1 with P+(

cS4 cP4

I)

(1 2 ),

lim P+ ( , k ) = lim 0

0

s

I)

+

+

P+(0) )

+

1 (1 2µ

I 1 2 2µ (1

I)

I 1 2 2µ (1

1 (1 2µ

(1 1 (1 2µ

I)

cS2 cP2

kh = 1 2

(1 I ) cS2 2cP2 2µ (1 I )

1

,

2 (0)

.

= 0,

zx

(14)

=0

(15)

=0

For the symmetrical waves the functions uz , zx , zy are odd, the functions u x , zz , u y are even with respect to the transversal coordinate measured from the middle and the vice versa in the case of the antisymmetrical waves. Thus, the boundary conditions (14), (15) on the lower face can be replaced by the following conditions on the middle plane

cS2 2cP2 I)

I

1 The corresponding limit of P+ as

(kP+(1)

µ

and the SH-waves satisfy the equation

P+(0) kS2 h2 1 8

s

1 µ

1

Let us use the propagator matrix obtained in the previous section to derive and analyze the dispersion equation for a three layered plate. The plate consists of a core and two identical coating layers. Medium 1 is considered as a core whereas upper and lower coatings are marked as medium 2 (see Fig. 1). In order to derive the dispersion equation, the equation of motion are to be solved under stress free boundary conditions on the upper and lower faces of the plate. Thus, the upper face of the coating in Fig. 1 in any case is stress free, but the same conditions on the lower face can be transformed into conditions on the middle plane. Namely, the waves with P- and SV- polarization satisfy the two equations on the faces

I)

1 (1 2µ

I

2 ),

kS2 h 4

5. Dispersion equations for the case of a three layered plate

1, and the limit of P as

kh + P (0) ) = 1 2

lim P ( , k ) = lim (kP (1) 0

I= 0 1 1 0

s

= O(

=

Let us also mention that the same limit results can be obtained directly from statics. Note, that according to our results, the limit value of the determi0 is a nant of any submatrix of the propagator matrix (10) as function which can be expressed as a linear combination of finite number of terms like k me klm (h1, h2, … , hn) , where m runs some integer values, and lm (h1, h2, …, hn ) are just linear combinations of the thicknesses. This function is single valued with respect to k . The problem to find its zeroes is reduced to the problem to search zeroes of the entire function of complex variable k . When such function cannot be represented as the product of a polynomial and another entire function with no zeroes, we are facing a general case of an infinite countable set of complex zeroes with an accumulation point at infinity |k| + . This means that for any finite distance R 0 > 0 the number of roots satisfying the inequality |k| > R0 is infinite and the number of roots inside the circle |k| R0 is finite. More details about these mathematical aspects can be found in [54].

(12)

where

J=

( s)

1

( 1)

0

cS2 2cP2

+

µ

,P

with respective limit lim P ( , k ) = P

P (0) =

P

1 µ

1

1 2

=

and

For the matrix P the respective first and zeroeth matrix terms acquire the form

P (1) =

(0)

0 yields

JT

1 T J 2µ

2µJ T

JT

uz = 0,

zx

= 0 (PSV symmetrical waves)

(16)

u x = 0,

zz

= 0, (bending waves)

(17)

I) I)

I)

I

.

Similarly, we obtain the submatrices of the rows and columns 5,6 which describe the wave motion with SH-polarization

P

P

(5656) = P

(0)

+ k 1P

( 1)

+ k 2P

( 2)

+…

(13) Fig. 1. Geometry of the coated plate.

For P ** the respective initial matrix terms acquire the form 4

Composite Structures 230 (2019) 111429

D.D. Zakharov and A.V. Nikonov

= 0, (SH symmetrical waves)

(18)

G 0 = k 3H30 + kG10 + O (µ1 h13 k S41 k 1) + O (µ1 h15 k S61 k 1)

u y = 0 (SH antisymmetrical waves)

(19)

where

zy

As one can see, the conditions on the upper face correspond to the rows 3,4,6 of the propagator matrix. The conditions on the middle plane correspond to the non-zero components of the motion-stress vector, i.e., the rows 1,3,5 correspond to the symmetrical modes, and the rows 2,4,6 correspond to bending. The dispersion equation is obtained from the equivalence to zero of the determinant of the submatrix in accordance with the homogeneous boundary conditions. Thus, we may write the respective dispersion equations as follows 346 135

=

34 6 13 · 5

= 0 (symmetrical waves )

(20)

346 246

=

34 6 24 · 6

= 0(bending ),

(21)

H11 =

H01 =

( ) = det P ( 6 ) , 5

346 = det P , 135

34 13

( )

34 = det P , 13

(13, 5

k S21 h1 h22 4µ 1 1 4µ1

(µ 2

(

b22

2µ1 µ 2 a12 a22 1 1 2

(

b12 a12

b12 (µ a12 2

a22

µ1) + µ 2 + µ1

2

b22 a22

µ1 ) 1 b 2b 2

2

b22 a22

µ1 ) 1

+ µ22 1 +

µ2 1 +

b12 a12

2

+ µ1 1 +

b14 a14

) + µ (1 ) (1 + ) b12

2 1

) (1 )

b22 2

a12

b24 a24

+

a22

+ O (µ1 h12 k S21 ) + O (µ1 h14 k S41 ),

with the determinants of the submatrices of the propagator matrix P 346 135

h22 (µ 2 µ1 2

H21 =

H12 =

6 5

h1 (µ 2µ1 2

H30 =

24, 6)

b12 a12

µ1 ) 1

2h1 h22 (µ 2 µ1

µ1 ) 2 1

1

b22 a22

b12 a12

1

b22 (µ a22 2 b22 a22

µ1) + µ 2 + µ1

2

0

G10 = H1 + O (µ1 h 2 k S2 )

6. Dispersion equation for symmetrical modes of a three layered plate

0

H1 =

The corresponding problem is reduced to solving the dispersion equation 34 13

(k ( ), )

34 34 34 34 = (e1 e2 )2A13 + (e1/ e2 )2B13 + (e2/ e1 )2C13 + (e1 e2) 2D13 + e12 G+1 + e22

(23)

G+2 + e1 2 G 1 + e2 2 G 2 + (e1 e2 )0G 0.

A0 + O (µ1 h k S2 k 1),

34 C13 = B 0 + O (µ1 h k S2 k 1),

34 B13 =

where A0 and B 0 .are the leading order terms for given by A0 =

1

B0 =

1

(b2/ a2 )2 8µ1

b 2b 2 2µ1 µ2 1 + 12 22 a1 a2

(b2/ a2 )2 2µ1 µ2 8µ1

+ µ12 1

b 2b 2 1 + 12 22 + µ12 1 a1 a2

b12 a12 b12 a12

b2 1 + 22 a2

34 A13 ,

34 B13 ,

b2 + µ22 1 + 12 a1

b2 b2 1 + 22 + µ 22 1 + 12 a2 a1

1

B0 =

µ2 2

Assuming

k2H21

kH11

G =

k 2H21 +

kH11

a22

(1 ) > 0. b22

+

+

+

G =

+

(24)

0 1 kH1 = 0 2

(25)

= = = 0, =0 In this specific case Let us consider the limit behavior of the stretching equation. H12

H30

B0

1) First, assume that media 1 and 2 are identical. Then, the static form 34 µ of the dispersion equation yields 13 2 2 (1 a b )[sinh2k (h1 + h2) + 2k (h1 + h2)] = 0 . +0 when the thickness of the coating 2) Next, consider the case h2 vanishes. In this case, we obtain the equation

O (µ1 h1 k S21 k 1)

O (µ1 h15 k S41 k 1)

kH12

H01sinh2kh1 +

H21

G+2 = kH12 + O (µ1 h12 k S21 ) 2

0 , we rewrite the Eq. (22) in the form

A0 sinh2k (h1 + h2)

a22

O (µ1 h13 k S41 k 1)

a22

b22 a22

A particular form of this equation for µ1 = µ 2 yields

H01 + O (µ1 h13 k S41 k 1) + O (µ1 h1 k S21 k 1)

H01

b22

B 0sinh2k (h1 h2) [k 2H21 + H01]sinh2 0 1 1 kh1 + kH12cosh2kh2 + k 3H30 + kH1 = 0 2 2

+ O (µ1 h15 k S61 k 1) 1

a22

A0 sinh2k (h1 + h2)

b22

The second group of the coefficients in (23) looks as follows

G+1 =

1 2 2

b22

7.1. Static limit of the dispersion equation

Here μα are the shear moduli, and the more compact notations b = c S , a = c p are introduced for the wave speeds in media 1 and 2. Note also that A0

b24 a24

a12

7. Analysis of the dispersion equation

34 and D13 ,

b22 a22

1

b12 a12

b12

0

= 1, 2)

34 C13

1 2 1

a22

A0 > 0 , and H30 < 0 , H1 < 0 . For each of the coefficients, B 0 and H01, there are two possibilities, e.g., B 0 > 0 or B 0 < 0 , and H01 > 0 or H01 < 0 with 2 1 1 signH1 = signH1 = signH2 = sign(µ 2 µ1) .

B 0 + O (µ1 h k S2 k 1)

34 D13 = A0 + O (µ1 h k S2 k 1), ( ,

a12

Remark 1. For the further search of complex wavenumbers we may note that the asymptotically leading terms have the following signs:

The explicit form of the asymptotics of the coefficients in (23) is obtained using the MAPLE code of symbolic computations. The first group of these coefficients acquires the form 34 A13 =

b22 2

b12

One may see that the static form of the dispersion Eq. (22) can be obtained by neglecting the terms of the asymptotic orders O (h 2 k S2 ), O (µ1 h k S2 k 1), O (µ1 h 2 k S2 ) or higher orders in the coefficients of (23). All the terms above have been obtained using symbolic computations.

where the exact expression in the left hand side of (22) follows from the explicit formula (11): 34 13

) (1 ) + 2µ µ h (1 ) (1 )+ + µ h (1 ) (1 + ) + 2µ µ h (1 ) . (

µ22 h1 1 2 1 1

(22)

=0

1 µ1

34 13

O (µ1 h12 k S21 ) 5

µ1 2

(1 ) (sinh2kh + 2kh ) = 0 b12 a12

1

1

Composite Structures 230 (2019) 111429

D.D. Zakharov and A.V. Nikonov

When the thickness of the central layer vanishes (h1 tain similarly

µ2

34 13

2

(

b22

1

a22

)

+ 0) we ob-

9. Symmetry of the wavenumbers

(sinh2kh2 + 2kh2) = 0

It is interesting to note that the expressions in the left hand side of k if we the dispersion Eq. (22) remain unchanged when replacing k 2 2 1/2 recall that EP, S = e h (k kP, S ) = e p,S kh . Indeed, the matrices M and M 1 explicitly depend on k 2 instead of k and this property immediately follows from the relations (5) and (7). Additionally, in the case of pure elasticity the dispersion equation possesses symmetry with respect to _ the replacement of k by its complex conjugation k . Note, that we specially single out the exponential terms in the relations (7) and (8), but they can be merged in the hyperbolic functions in the final expression for P (h) . To be brief, we may refer to a similar expression in [52], chapter 9, page 397 which differs from our form by the normalization, the order of terms in the motion-stress vector and by the use of the ordinary trigonometric (non hyperbolic) functions cosine and sine. Note also that our hyperbolic functions appear in the form cosh(kh P , S ), sinh(kh P, S )/ P, S , or P, S sinh(kh P, S ) and they can be represented by respective power series with final even powers of radicals. Since all the coefficients in the matrices P and P (h) , and the coefficients in the power series of the mentioned functions are real, we immediately arrive at the conjugation symmetry. Hence, the wavenumbers are placed symmetrically by the quadrants on the complex plane and it is sufficient to consider the first quadrant Rek 0 , Imk 0 .

Therefore, each of the equations coincides with the equation for a single layer of a thickness h in statics sinh2kh + 2kh = 0 , within a constant factor (see, e.g., [3,4,5]). 7.2. General case of h2

0 and µ1

µ2

Introducing a new variable = 2k (h1 + h2) and the parameters

h2 ,1 h1 + h2

=

h1 ,1 h1 + h2

=

(dimensionless

2 =

wavenumber)

h1 h2 h1 + h2

(26a)

1 1 2 H21 2H01 H12 2B0 , h2 = , h0 = , h1 = 2 A0 2A0 (h1 + h2 ) A0 A0 (h1 + h2 )

B0 =

(26b)

0

0

h3 =

0 H30 H1 , h1 = 8A0 (h1 + h2 )3 2A0 (h1 + h2 )

(26c)

the Eq. (24) can be rewritten as follows

B0 sinh (1 2

sinh +

3h

cosh

2h

2 )+

0 3

1 2

1

2

+ h0 sinh (1 2

h1 2

)

1

h0 =0 2

2

Remark 2. In general, for the case of linearly viscoelastic materials (e.g., the Kelvin-Voight or Maxwell media) the dispersion equation at 0 keeps the symmetry k k but loses the symmetry of complex k¯ . However, the static limit of the dispersion equation conjugation k still possesses both kinds of symmetry.

(27)

or in the equivalent form below (28)

e =f( ) f( ) 3

=

0

0

|h 3 |

sinh(

|h1 | + e

1

1

1) + [ 2h 2 + h0 ]

+ B0 sinh(2

2

10. Low frequency asymptotics of the dispersion curves

(29)

1) + h1 cosh

In this section we remind some duality of the function (k, ) in the left hand side of the dispersion Eq. (22) which can be rewritten in the form

or in the logarithmic form

= Lnf ( ) = ip + ln f ( )

(p = 2

(30)

n)

In what follows, let us call the integer number n as the index of the current complex root n . As seen, the root n is determined by its index n and by the behavior of the function f ( ) in (29) and (30).

(k ( ), ) =

kn (0) = kn (

+

(e1/ e2 ) B56

At the limit

A56 =

D56 =

+

(e2/ e1 ) C56

+ (e1 e2)

1D 6 5

=0

kn2 (

(31)

4

,

B56 =

C56 =

µ1

µ2 4

kn (

and the Eq. (31) can be rewritten in the form

(µ1 + µ 2 ) sinh k (h1 + h2 ) + (µ1

µ 2) sinh k (h1

h2) = 0

e (1 ) ], M 0 (µ1 where g ( ) 1 + M 0 [e g ( ). the form (30) when replacing f ( )

=0

=

n/2(h1

+ h2)

2)

2 2 d (k n ) d ( 2)

= kn2 (0) +

+ =0

d 2 (kn2 ) 2! d ( 2) 2 4

2:

+…

(35)

=0

2)

= kn (0) 1 +

d (kn2 ) kn2 (0) d ( 2) 2

+

2

×o

=0

1/2

2

kn2 (0)

and the next terms in the series (35) provide the smaller order of frequency ω. Omitting the cumbersome formulas, let us formulate the main results. It can be shown that

(32)

+0 , or h2 +0 , the Eq. At the evident limit cases µ1 µ 2 , or h1 (32) is converted to one of the equations (2µ ) sinh k (h1 + h2) = 0 , or (2µ 2 ) sinh kh2 = 0 , or (2µ1) sinh kh1 = 0 , respectively. Again, we arrive at the coincidence of our results with the known equation for a single layer sinh kh = 0 . In the previous dimensionless notations, the Eq. (32) is reduced to the equivalent form

e = g( )

2 )|

Thus, if the values of derivatives in (35) are restricted as |kn (0)| + 0 then

0 the obtained coefficients satisfy the relations

µ1 + µ 2

(34)

2)

expand the function into the Maclaurent power series of

For the case of SH-polarization of the displacement vector, the dispersion equation acquires the form

(e1 e2 ) A56

2 ),

Assuming an arbitrary index n of the complex dispersion curve kn ( ) (in fact, kn ( 2) ) with its initial value in statics

8. Dispersion equation and its static limit for SH-modes of a three layered plate

6 5

(k 2 (

/ (kn2) = F1 (kn) + O (

2 ),

/ (

2)

= F2 (kn ) + O (

2 ),

(

0) (36)

d (kn2 ) d ( 2)

(33)

= =0

/ ( 2) / (kn2 )

= =0

F2 (kn (0)) F1 (kn (0))

(37)

Retaining two terms in the series (35), we obtain the following lowfrequency asymptotics for the wavenumber kn knas ( )

µ 2 )(µ1 + µ 2) 1, or to 6

Composite Structures 230 (2019) 111429

D.D. Zakharov and A.V. Nikonov

kn2 (0)

12. Discussion and conclusive remarks

1/2

/ ( 2) / (kn2 )

2

knas = kn (0) 1

(38)

=0

An algorithm to treat the dispersion equation on the base of the propagator approach is suggested. The form of the propagator in use is especially convenient to analyze the case of small frequency and large wavenumber. The implementation of the algorithm is demonstrated for the case of symmetrical modes of a three-layered plate and it can be generalized for the case of more layers and unsymmetrical layup. At the first step, the dispersion equation is obtained in a closed form with the aid of codes of symbolic computation though the explicit terms are too lengthy and are not presented. At the second step, the static limit of the dispersion equation is obtained in a closed form. Then we analyze its complex dimensionless roots n as . At the third step, we focus our attention on the Re n, Im n asymptotic behavior of n and on the iterative procedure of the root calculations. Forth, the leading low frequency approximation of the dispersion curves is obtained using symbolic computations. The results of the Part I are to be used in Part II for finding the largewavenumber static asymptotics and low-frequency dispersion curves. Thus, the total infinite set of wavenumbers can be built up at low frequency.

+ , the If the derivative (37) decreases or is restricted as n complex dispersion curve should have a long flat segment at low frequency. The respective estimation is shown in the Part II. For the completeness, let us say a few words about the real wavenumbers which lay on the fundamental curve S0 in the low frequency region. As shown in [55–57], even for the case of an anisotropic multilayered plate with nonsymmetrical stacking sequence it is possible to obtain asymptotically accurate low frequency equation of motion and respective dispersion relations. For the S0 -mode and for the case of isotropic layers, it yields the familiar relation 2

µ j (z j + 1

zj ) = 2k 2

j (z j + 1 j

zj )

1

j

(39)

j

It can be used together with the approximate formula for the fundamental bending curve A0 2

j (z j + 1

µj (z 3j + 1

2 3

zj ) =

j

1

j

z 3j )

k4

(40)

j

under the additional condition for the position of the origin

µj (z j2+ 1

z j2 )

1

j

=0

Acknowledgments This work was partially supported by the charity fund “Osnovanie” and partially by the Slovenian Research Agency (ARRS) which are gratefully acknowledged. The authors are very thankful to the Reviewers for the thorough analysis of the manuscript and for the help to improve the paper.

(41)

j

Here j are Poisson’s ratios. This condition is automatically satisfied for the symmetrically-layered plate and becomes a linear equation with respect to z1 at the given thicknesses z j + 1 zj = hj .

References

11. Asymptotic behavior of the static roots

[1] Mindlin RD (Onoe M, Medick MA). Mathematical theory of vibrations of elastic plates. In: Proc. 11th Annual Symp. on Frequency Control, Fort Monmouth, Army Signal Engineering Laboratories (USA); 1957. 1–40. [2] Mindlin RD. An introduction to the mathematical theory of vibrations of elastic plates. Singapore: World Scientific; 2006. [3] Dougall J. An analytical theory of the equilibrium of an isotropic elastic plate. Trans R Soc Edinburgh 1904;41:129–228. [4] Zlatin AN. Roots of several transcendental equations encountered in elasticity theory. Sov Appl Mech 1980;16(12):69–74. (in Russian). [5] Merkulov LG, Rokhlin SI, Zobnin OP. Calculation of the spectrum of wave numbers for Lamb waves in a plate. Sov J Nondestruct Test 1970;6:369–73. [6] Viktorov IA. Rayleigh and Lamb waves. New York: Plenum Press; 1967. [7] Auld BA. Acoustic Fields and Waves in Solids 2. Krieger Inc; 1990. [8] Dougall J. An analytical theory of the equilibrium of an isotropic elastic rod of circular section. Trans R Soc Edinburgh 1914;49:895–978. [9] Zemanek J. An experimental and theoretical investigation of elastic wave propagation in a cylinder. J Acoust Soc Am 1972;51:265–83. [10] Rokhlin S. Diffraction of Lamb waves by a finite crack in an elastic layer. J Acoust Soc Am 1980;67:1157–65. [11] Rokhlin S. Resonance phenomena of Lamb waves scattering by a finite crack in a solid layer. J Acoust Soc Am 1981;69:922–8. [12] Gregory RD, Gladwell I. Reflection of a symmetric Rayleigh-Lamb wave at a fixed or free edge of a plate. J Elast 1983;13:185–206. [13] Crane LJ, Gilchrist MD, Miller JJH. Analysis of Rayleigh-Lamb wave scattering by a crack in an elastic plate. Comput Mech 1997;19:533–7. [14] Shkerdin G, Glorieux C. Lamb mode conversion in a plate with a delamination. J Acoust Soc Am 2004;116:2089–100. [15] Shkerdin G, Glorieux C. Lamb mode conversion in an absorptive bi-layer with a delamination. J Acoust Soc Am 2005;118:2253–64. [16] Flores-Lopez MA, Gregory RD. Scattering of Rayleigh-Lamb waves by a surface breaking crack in an elastic plate. J Acoust Soc Am 2006;119(4):2041–9. [17] Zakharov DD. Dirichlet-Neumann conditions and the orthogonality of 3D guided waves in layered solids. Comput Mathemat Math Phys 2010;50(9):1522–35. [18] Auld A, Kino GS. Normal mode theory for acoustic waves and its application to interdigital transducer. IEEE Trans Electron Dev 1971;ED-18(10):898–908. [19] Bobrovnitskii YI. Orthogonality relations for Lamb waves. Soviet Acoust Phys 1973;18(4):432–3. [20] Fraser WB. Orthogonality relations for Rayleigh-Lamb modes of vibration of a plate. J Acoust Soc Am 1976;59:215–6. [21] Prakash BG. Generalized orthogonality relations for rectangular strips in elastodynamics. Mech Res Commun 1978;5(4):251–5.

Let us discuss the behavior of the input value n for the asymptotics + as (38). The logarithmic Eq. (30) contains a large parameter p n + . The function f ( ) has mixed polynomial-exponential form (29) and the coefficients of n in its exponents are real and less than 1. So, we may expect the asymptotic behavior (42)

Im = O (p), Re = O (ln p)

+ . In particular, the dimensionless roots of the static equaas p tions for a single isotropic layer and for a circular rod have exactly the same behavior [3–5,8,9]. The iteration scheme [49], which was applied + for the case of a single layer as (0) n

( + 1) n

= ip, = i [an(

+ 1)

= ip + ln f (

p + bn(

results in the limit

+ 1)

( ) n )

] + cn(

(as ) n

+ 1)

ln p + dn(

= lim +

( ) n

(43)

+ 1)

satisfying the static equation as

+ where the error n behaves as O (p ) for any 0 < <1. The exact roots n were then found using n(as) with a bit tricky form of the logarithmic Eq. (30) n

= ip + (

as n

ip) + ln[e

(as ) n

as n + ipf

( n)] =

as n

+ ln[e

as n f

( n)]

(44)

and the following refined recurrent procedure ( + 1) n

=

as n

+ ln[e

as n f

(

( ) n )],

(0) n

=

as n

(45)

Of course, in general this algorithm cannot directly repeat all the simplest derivations for such a case and needs additional analysis and improvement. However, we may follow the same main idea as in [49] and consider the sequence (43) with some optimization steps in the Part II. 7

Composite Structures 230 (2019) 111429

D.D. Zakharov and A.V. Nikonov [22] Slepyan LI. Betti theorem and orthogonality relations for eigenfunctions. Mech Solids 1979;1:83–7. [23] Zakharov DD. Generalised orthogonality relations for eigenmodes in 3D dynamic problem for elastic layer. Mech of Solids [MTT] 1988;6:62–8. [24] Zakharov DD. Orthogonality of 3D guided waves in viscoelastic laminates and far field evaluation to a local acoustic source. Int J Solids Struct 2008;45(6):1788–803. [25] Zakharov DD. The property of orthogonality and energy transfer by 3D eigenwaves in transversely isotropic laminated plates with and without contact with a fluid. J Appl Math Mech [PMM] 2013;77(1):39–50. [26] Li L, Haider MF, Mei H, Giurgiutiu V, XiaY.. Theoretical calculation of circularcrested Lamb wave field in single- and multi-layer isotropic plates using the normal mode expansion method. Struct Health Monitor 2019:1–16. https://doi.org/10. 1177/1475921719848149. [27] Giurgiutiu V, Haider MF. Propagating, evanescent, and complex wavenumber guided waves in high-performance composites. Materials 2019;12(269):1–33. [28] Miamoto T, Yasuura K. Numerical analysis on isotropic elastic waveguides by mode-matching method. IEEE Trans Sonics Ultrason 1977;SU- 24:359–75. [29] Hayashi T, Tamayama C, Murase M. Wave structure analysis of guided waves in a bar with an arbitrary cross-section. Ultrasonics 2006;44:17–24. [30] Hayashi T, Song WJ, Rose JL. Guided wave dispersion curves for a bar with an arbitrary cross section, a rod and rail example. Ultrasonics 2003;41:175–83. [31] Bondarenko AA. Normal modes of the elastic waveguide with rectangular crosssection. Akustichnyi visnik 2007;10(4):12–27. (in Russian). [32] Thomson WT. Transmission of elastic waves through a stratified solid medium. J Appl Phys 1950;21:89–93. [33] Haskel NA. The dispersion of surface waves on multilayered media. Bull Seismol Soc Am 1953;43:17–34. [34] Knopoff AL. A matrix method for elastic wave problem. Bull Seismol Soc Am 1964;54(1):431–8. [35] Schwab F, Knopoff AL. Surface waves in multilayered anelastic media. Bull Seismol Soc Am 1971;61(4):893–912. [36] Lowe MJS. Matrix techniques for modeling ultrasonic waves in multilayered media. IEEE Trans Ultrason Ferroelectr Freq Control 1995;42:525–41. [37] Predoi MV, Castaings M, Moreau L. Influence of material viscoelasticity on the scattering of guided waves by defects. J Acoust Soc Am 2008;124(5):2883-2894. [38] Mazzotti M, Marzani A, Bartoli I, Viola E. Guided waves dispersion analysis for prestressed viscoelastic waveguides by means of the SAFE method. Int J Solids Struct 2012;49:2359–72. [39] Gravenkamp H, Birk C, Song C. Simulation of elastic guided waves interacting with defects in arbitrarily long structures using Scaled Boundary Finite Element Method. J Comp Phys 2015;295:438–55. [40] Krome F, Gravenkamp H. Analyzing modal behavior of guided waves using high order eigenvalue derivatives. Ultrasonics 2016;71:75–85. [41] Hernando Quintanilla F, Fan Z, Lowe MJS, Craster RV. Guided waves’ dispersion

[42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56] [57]

8

curves in anisotropic viscoelastic single- and multi-layered media. Proc R Soc A 2015;471(2183):1–23. Hernando Quintanilla F, Fan Z, Lowe MJS, Craster RV. Modeling guided elastic waves in generally anisotropic media using a spectral collocation method. J Acoust Soc Am 2015;137:1180–94. Hernando Quintanilla F, Fan Z, Lowe MJS, Craster RV. Full 3D dispersion curve solutions for guided waves in generally anisotropic media. J Sound Vibr 2016;363:545–59. Rokhlin SI, Wang L. Stable recursive algorithm for elastic wave propagation in layered anisotropic media: Stiffness matrix method. J Acoust Soc Am 2002;112:822–34. Rokhlin SI, Wang L. Recursive geometric integrators for wave propagation in a functionally-graded multilayered elastic medium. J Mech Phys Solids 2004;52:2473–506. Rokhlin SI, Wang L. Modeling of wave propagation in layered piezoelectric media by a recursive asymptotic method. IEEE Trans Ultrason Ferroelectr Freq Control 2004;51:1060–71. Wang L, Rokhlin SI. Recursive asymptotic stiffness matrix method for analysis of surface acoustic wave devices on layered piezoelectric media. Appl Phys Lett 2002;81:4049–51. Pagneux V, Maurel A. Determination of Lamb mode eigenvalues. J Acoust Soc Am 2001;110(3):1307–14. Zakharov DD, Castaings M, Singh D. Numerical and asymptotic approach for evaluating complex wave-numbers of guided modes in viscoelastic plates. J Acoust Soc Am 2011;130(2):764–71. Zakharov DD. Parametric analysis of complex dispersion curves for flexural Lamb waves in layered plates in the low-frequency range. Acoust Phys 2018;64(4):387–401. Achenbach JD. Wave propagation in elastic solids. Amsterdam: Elsevier; 1984. Aki K, Richards PG. Quantitative seismology, theory and methods Vol. 1. San Francisco: Freeman; 1980. Zhang B. Elastic waves excited by a plane source on the surface of a multilayered medium. J Acoust Soc Am 2007;121(3):1440–8. Yu B. Levin/in collaborations with Lubyanskii Yu, Sorokin M, Tkachenko V/ Lectures on entire functions. Providence: American Mathematical Society; 1996. Zakharov DD. Asymptotic analysis of three dimensional dynamic elastic equations for a thin multilayered anisotropic plate of arbitrary structure. J Appl Math Mech 1992;56(5):637–44. Zakharov DD. Asymptotic integration of 3D dynamic equations for thin multilayered anisotropic plates. Comptes Rendus de l’Academie des Scences II 1992;315(8):915–20. Zakharov DD, Becker W. 2-D Problems of thin asymmetric laminates. Zeut fur Angew Matematik und Physik 2000;51(4):49–66.