The Open Access NDT Database

www.ndt.net/?id=10704

A Viscoelastic plate theory for the fast modelling of Lamb wave solutions in
NDT/SHM applications
Miguel A. TORRES-ARREDONDO 1,2, Claus-Peter FRITZEN 1,2
1
University of Siegen, Centre for Sensor Systems; Siegen, Germany
2
University of Siegen,. Institute of Mechanics and Control Engineering-Mechatronics; Siegen, Germany
Phone: +49 271 7402132, Fax: +49 271 7402769; e-mail: {torres, fritzen}@imr.mb.uni-siegen.de

Abstract
Guided ultrasonic waves have many useful properties that can be exploited for non-destructive testing (NDT)
More info about this article: https://www.ndt.net/?id=10704

and structural health monitoring (SHM) applications. However, much information and analysis regarding the
generation and propagation of these waves is needed before automatic processing and analysis techniques can
provide useful information for reliable fault monitoring. Moreover, the knowledge of dispersion characteristics is
crucial for the optimization of sensor networks in terms of sensor placement and number of sensors. On that
account, the present work introduces a higher order plate theory for modelling disperse solutions in viscoelastic
fibre-reinforced composites. This approach offers a higher computational efficiency and simplicity in
comparison to traditional exact elasticity methods, while providing an adequate description of the structure's
global response in the low frequency range which is the most used in Lamb wave applications. The proposed
method was applied to several examples in order to obtain numerical results in elastic and viscoelastic
anisotropic plates. Some comparisons to experimental data are presented, and the effectives and limitations of
the method are discussed.

Keywords: Lamb Waves, Plate Theory, Viscoelasticity, Fibre Composite, Structural Health Monitoring

1. Introduction
Guided ultrasonic waves are a valuable tool in order to get information regarding the origin
and importance of a discontinuity in a structure for a longer safe life and lower operation costs
[1]. Depending on the material attenuation and geometric beam spreading effects, guided
waves are able to propagate over relative long distances, interact sensitively with and/or being
related to different types of defects like e.g. delaminations, corrosion damage, etc. However, it
is only possible to benefit from these advantages once the complexity of guided wave
propagation are disclosed. Generally, there are two ways to obtain such complex dispersion
characteristics. The first strategy is based on common experimental time-delay measurements
where piezoelectric transducers are attached to a structure and play the role of either actuator
or sensor. A second approach is based on different modelling approaches, where material
constants are measured by means of different identification or material testing techniques and
then fed into an analytic or finite element model to extract the relevant information. The
situation considered in this paper is representative of the latter case.

The modelling of wave propagation in multilayered anisotropic structures has been

extensively studied by several researchers and a considerable amount of literature has been
published on this topic [2]. Analyzing guided waves in these structures is often categorized
into three methods. There are methods based on exact three dimensional elasticity, waveguide
finite element methods and laminated plate theories of different orders. Exact methods are
based on the superposition of bulk waves that include the popular matrix based methods [3].
Waveguide finite element methods have appeared for modelling the guided wave propagation
numerically as an alternative to exact methods by using a finite element discretization of the
cross-section of the waveguide [4]. A different alternative providing simplicity and low
computational cost in comparison to other techniques are laminated plate theories. These
theories expand the displacement fields in terms of the thickness to any desired degree and
reduced the 3-D continuum problem to a 2-D problem. There are many papers dealing with
plate theories and their applications; some notable work in this field is given by Reddy in [5].
The present work proposes a third order plate theory that can approximate five symmetric and
six antisymmetric Lamb wave modes. Here, the two classical models of viscoelastic
attenuation are implemented and discussed. Previous related work to the one proposed here
using a second order expansion of the displacement fields in terms of the thickness was
presented by Calomfirescu in [6]. However, this work was restricted to the analysis of just the
fundamental modes of propagation and the use of one model of attenuation for computational
reasons. The motivation for expanding the displacement field up to the cubic term in the
thickness is to provide better kinematics and accurate interlaminar stress distributions [7].
Moreover, the complex root finding algorithms needed with the exact methods are not
required here since analytical solutions can be obtained and reduce to the solution of simple
polynomials. Since the use of composite materials has extensively increased in the design of
existing engineering structures, what also increases the analysis complexity of such structures,
this poses a necessity for fast modelling tools that can be used for a rapid and reliable
analysis.

2. Material Damping Models for Wave Propagation

The mathematical modelling of viscoelasticity using ideas from elasticity has attracted the
attention of a large number of investigators over the past century. In order to account for
material damping, the stiffness matrix is represented by a complex quantity. The real part C of
this complex term relates to the elastic behaviour of the material and defines the stiffness. The
imaginary component η relates to the material viscous behaviour and defines the energy
dissipative ability of the material. Two models are often used to describe viscoelastic
behaviour. The first model is called the hysteretic model whose complex stiffness matrix is
given by
Cɶ = C + iη . (1)

The hysteretic model assumes no frequency dependence of the viscoelastic constants. The
second model is the Kelvin-Voigt model and assumes a linear dependence of the viscoelastic
coefficients. The complex stiffness matrix is expressed as
ω (2)
Cɶ = C + i η
ωɶ
where ω is the angular frequency and ῶ is the frequency of characterization. The influence in
the attenuation predicted by both methods for the fundamental modes of propagation is
depicted in Figure 1. The characterization frequency is 2MHz.
2500
SH0 Hysteretic Model
S0 Hysteretic Model
2000 A0 Hysteretic Model
SH0 Kelvin-Voigt Model
Attenuation [dB/m]

S0 Kelvin-Voigt Model
1500 A0 Kelvin-Voigt Model

1000

500

0
0 0.5 1 1.5 2 2.5 3
Frequency [MHz]

Figure 1. – Comparison of attenuation between the hysteretic and Kelvin-Voigt models as a function of
frequency.
From the previous picture it can be clearly seen that the attenuation is a linear function of the
frequency in the case of the hysteretic model and a quadratic function of the frequency in the
case of the Kelvin–Voigt model. Additionally, both models are just coincident in the
frequency of characterization and away from this point, the deviation in the prediction of both
models is noteworthy.

3. Mathematical Framework for the Plate Theory

The model considers a linearly viscoelastic, non-piezoelectric layer of material subjected to a
complex stress system in three dimensions. The material is considered to have a monoclinic
symmetry. Figure 2 depicts the definition of stress resultants ( N, M, Q ) in the three
dimensional system for a given propagation direction θ and fiber orientation φ.

Qy M Ny
y
M yx
Mx N yx
M xy z
z y Qx
Nx φ N xy
θ
Qx
N xy x Mx
My Nx
N yx M xy

Qy
Ny M yx

Figure 2. – Complex Stress System in Three Dimensions.

The approximated displacement fields are given by

u = u0 ( x, y, t ) + zψ x ( x, y, t ) + z 2φx ( x, y, t ) + z 3λx ( x, y, t )
v = v0 ( x, y, t ) + zψ y ( x, y, t ) + z 2φ y ( x, y, t ) + z 3λ y ( x, y , t ) (3)
w = w0 ( x, y , t ) + zψ z ( x, y, t ) + z 2φz ( x, y, t )

where u, v, and w are the displacement components in x, y and z directions, ѱx and ѱy

represent rotations having the same meaning as in the first order shear deformation theory
[8]. The additional terms expand the displacement field. The strain energy of each layer can
be represented as [9]
1
U=
2V∫ (C11ε x2 + 2C12ε xε y + 2C13ε xε z + 2C16ε xγ xy + C22ε y2 + 2C23ε y ε z + 2C26ε y γ xy
(4)
+ C33ε z2 + 2C36ε z γ xy + C66γ xy
2
+ C44γ yz
2
+ 2C45γ yz γ xz + C55γ xz
2
)dV .

Additionally,
∂U ∂U (5)
σx = ,...,τ xy = .
∂ε x ∂γ xy

The plate constitutive equations may be derived from the strain energy density in the 3-D
elasticity theory and the linear elastic stress-strain and strain-displacement relations. The
equations of motion may be derived from the dynamic version of the principle of virtual
displacement. Consequently, the system can be expressed in a matrix form, and by imposing
boundary conditions and setting its determinant to zero, a characteristic function relating the
angular frequency to the wavenumber is obtained. For the case of symmetric laminates, the
system of equations can be decoupled into two independent systems of equations for the
symmetric and antisymmetric modes of propagation. A comparison between exact solutions
and a similar elastic third order plate theory is provided by Wang and Yuan in [10]. However,
this paper does not document the relevant results needed for the calculation of the dispersion
relations and uses different shear correction factors to the ones provided here . A complete
description of the numerical strategy for the tracing of the dispersion solutions is presented in
[11]. The complete analytical expressions are given in the appendix of this paper.

4. Results

The proposed viscoelastic plate theory formulation is applied to several examples including
two anisotropic elastic plates and two anisotropic viscoelastic plates.

4.1 Elastic Glass Fibre Reinforced Plastic Plate

In order to validate the modelling approach, a case study has been conducted on a
unidirectional glass-fibre reinforced plastic (GFRP) plate. A single-layered specimen was
selected because of its highly anisotropic character. Figure 3(a) shows the structure that has
the dimensions 800mm×800mm and a thickness of approximately 1.5mm. Nine piezoelectric
transducers are attached to the surface of the structure with equidistant spacing. The piezo
patches have a diameter of 10mm and a thickness of 0.25mm. The elastic properties in the
principal directions of material symmetry provided by the manufacturer are given in Table 1.

E1 (GPa) E2 (GPa) E3 (GPa) G12 (GPa) G13 (GPa) G23 (GPa) ѵ12= ѵ13= ѵ23 ρ (kg/m3)
30.7 15.2 10 4 3.1 2.75 0.3 1700

Table 1. – Material properties of unidirectional glass fibre reinforced plastic.

The experimental group velocities were determined in the defined frequency by means of
time-delay measurements. Numerical results for the group velocities for the fundamental
modes of propagation at a central frequency of 200kHz are depicted in Figure 3(b). The wave
surface for the S0 mode at 200kHz is shown in comparison with some measured values at
discrete angular points (black circles) in order to validate the analytical model with
experimental data. It can be seen that the estimated group velocity matches the theoretical
curve very well, demonstrating the effectiveness of the model.
SHo Mode
So Mode 90 5 [km/s]
Ao Mode 120 60
4
Experimental So Mode So
3
150 30
2
1
Ao
180 SHo 0

210 330

240 300
270

(a) (b)
Figure 3. – Propagation Modes: (a) Experimental Setup and (b) Wave Surface for the 1.5mm thick GFRP in
vacuum.
4.2 Elastic Carbon Fibre Reinforced Plastic Plate

A carbon fibre reinforced plastic (CFRP) plate made of 16 equal layers is analyzed in this
section. The stacking sequence is [0 90 -45 45 0 90 -45 45]s and the total thickness is 4.2mm.
The nominal material parameters of the unidirectional layers provided by the manufacturer
are given in Table 2. Figure 4(a) shows the plate of approximately 500mm x 500mm
instrumented with nine piezo electric transducers. A 3 cycle tone burst signal with a 60kHz
centre frequency was used as the input waveform. The piezoelectric transducer number five
was designated as the actuator. Figure 4(b) depicts the signals received by the sensors PZT4
and PZT2. In contrast to the previous GFRP plate, here it can be observed by checking the
times of arrival and shapes of the received signals that the studied CFRP plate has an evident
lower anisotropic behaviour (see also Figure 4(c)). Furthermore, the fundamental
antisymmetric mode is more strongly excited than that of the fundamental symmetric mode.
The modes of propagation were identified by means of time-frequency analysis.

E1 (GPa) E2 (GPa) E3 (GPa) G12 (GPa) G13 (GPa) G23 (GPa) ѵ12= ѵ13= ѵ23 ρ (kg/m3)
155.0 8.5 8.5 4.0 4.0 4.0 0.33 1600

Table 2. – Material properties of multilayered carbon fibre reinforced plastic.

The calculated group velocities at θ=30° are depicted in Figure 4(d). It can be noticed that the
behaviour of the SH0 and S0 modes is different from the A0 mode in both the low and high
frequency zones. In the relatively low frequency range the higher the frequency of the A0
mode, the faster its group velocity. In an opposite manner for the S0 mode, the higher its
frequency, the slower its group velocity.

20 PZT 4
PZT 2
15 A0
10
Voltage [mVolts]

-5
S0
-10 +
-15 S0reflect Reflections
-20
0 0.1 0.2 0.3 0.4 0.5
Time [ms]
(a) (b)
90 8 km/s
120 60
6
150 4 30

2 So
SHo
180 0
Ao

210 330

240 300
270
(c) (d)

Figure 4. – Propagation Modes: (a) Experimental Setup, (b) Recorded Signals at the PZT Sensors, (c) Wave
Surface at 60kHz, (d) Group Velocity dispersion curve.
4.3 Viscoelastic Orthotropic Carbon-Epoxy Plate

A viscoelastic orthotropic unidirectional lamina of 1mm of thickness is analyzed in this

section. This example was fully studied in references [12-13]. The elastic and viscoelastic
material properties are given in Table 3. The material was characterized at 2.242MHz by the
use of ultrasonic waves transmitted through a plate-shaped sample immersed in water.

C11 C12 C13 C22 C23 C33 C44 C55 C66

132 6.9 12.3 5.9 5.5 12.1 3.32 6.21 6.15
η11 η12 η13 η22 η23 η33 η44 η55 η66
0.4 0.001 0.016 0.037 0.021 0.043 0.009 0.015 0.02

Table 3. – Material properties of unidirectional Carbon-Epoxy 1mm thick plate (units in GPa) [12].

Figure 5 presents the Lamb wave results obtained for both the hysteretic and the Kelvin-Voigt
viscoelastic models. It can be clearly seen from Figure 5(a) and (b) that both models do not
affect the phase velocity results in a substantial manner. The results obtained with the
proposed method are in good agreement with those obtained in references [12-13] for the A0
and SH0 fundamental modes of propagation. Although not depicted here (see [13]), the
approximate solutions of S0 mode have a good agreement with exact solutions at the lower
frequency range, but a big divergence of the S0 mode occurs at the zone of high dispersion
around 1.5MHz. For the higher order modes, solutions can be tracked very well at the
beginning of zones of high dispersion; however, after these zones, results are not longer
accurate and the higher order plate theory fails in providing good estimates for the velocities.
Figure 5(c) and (b) show that the effect in the prediction of attenuation for both models is
appreciable when the working frequency is different from the characterization frequency.
Hysteretic Model Kelvin-Voigt Model
20 20

A2 A3 A2 A3

15 SH1 A1 15 SH1 A1
Phase Velocity [km/s]

Phase Velocity [km/s]

S2 S2

So So
10 10
SH2 S1 SH2 S1

SH3 SH3
5 SHo 5 SHo

A0 A0
0 0
0 1 2 3 4 5 0 1 2 3 4 5
Frequency [MHz] Frequency [MHz]
(a) (b)
90 90 100 dB/m
200 dB/m
120 120 60
60 80
150
60
150 100 30 150 30
40
A0 A0
SH1 50 SH1 20
180 0 180 0

So So
210 330 210 330
SHo SHo
240 300 240 300
(c) 270 (d) 270
Figure 5. – Comparison between the Hysteretic and Kevin-Voigt Models: (a) and (b) Phase Velocity Dispersion
Curves, (c) and (d) Attenuation Polar Plots at 500MHz.
4.4 Viscoelastic Unidirectional Carbon-Epoxy Plate

A viscoelastic unidirectional carbon-epoxy panel of 3.6mm thickness is analyzed here. This

example was chosen since it was fully studied by Neau in [14] by using a method based on
the exact elasticity theory. The elastic and viscoelastic material properties are given in Table
4. The material was characterized at 2MHz using conventional interferometry methods.

C11 C12 C13 C22 C23 C33 C44 C55 C66

86.6 9 6.4 13.5 6.8 14 2.72 4.06 4.7
η11 η12 η13 η22 η23 η33 η44 η55 η66
7.5 0.3 0.6 0.6 0.25 0.28 0.1 0.12 0.28

Table 4. – Material properties of unidirectional Carbon-Epoxy 3.6mm thick plate (units in GPa) [14].

Figure 6 depicts the attenuation of Lamb modes in polar coordinates for the 3.6mm thick
carbon-epoxy plate. The hysteretic model was used for this numerical example. Although not
depicted here, the results obtained with the proposed method (Figure 6(a) and (b)) are
coincident with those obtained by Neau in [14] for the fundamental modes of propagation at
100kHz. However, as one moves along the frequency axis to higher frequencies, i.e. 500kHz ,
the predictions for the symmetric modes of propagation start to highly deviate from the exact
solutions (Figure 6(c) and Ref. [14]). Nevertheless, the good accuracy for the antisymmetric
modes of propagation still holds for this frequency (Figure 6(d)). This can be explained by the
fact that the displacement field for the flexural motion is one order higher than that of the
extensional motion, what assures a better approximation of the flexural modes. It can also be
seen that the mode attenuation is strongly related to the frequency and orientation of
propagation. This is a common characteristic of anisotropic materials where the velocity,
attenuation and energy of propagation of the multiple modes are both frequency and angle
dependent.
90 100 dB/m 90 80 dB/m
120 120 60
80 60
60
Ao
60
150 30 150 40 30
40
SHo 20
20
So
180 0 180 0

210 330 210 330

240 300 240 300

270 270

(a) 90 800 dB/m (b) 90 500 dB/m

120 120
600 60 400 60
S1
400 300
150 30 150 30
SHo 200
200 A2 100
So A1
180 0 180 0
Ao

210 330 210 330

240 300 240 300

(c) 270 (d) 270

Figure 6. – Attenuation of Lamb waves: (a) and (b) Attenuation Polar Plot at 100kHz, (c) and (d) Attenuation
Polar Plot at 500kHz.
5. Discussion and Conclusions
A coupling between viscoelasticity theory and a laminated plate theory has been suggested
which is applicable to viscoelastic fibre reinforced composite materials for the calculation of
wave velocities and attenuation for the different modes of propagation. The proposed theory
poses a compromise between low computational cost and accuracy in the results. It has been
shown that the method has provided good estimates of velocity and attenuation in anisotropic
laminates in the frequency range of Lamb wave applications. Comparisons to experimental
data and results from literature have been presented in order to validate the model. New shear
correction coefficients have been introduced which provide a better matching of the
frequencies from the approximate theory to frequencies obtained from the exact theory.
However, it was also depicted that the model suffers from some limitations that prevent it
from being used to solve the whole spectrum of composite laminate problems, i.e. high
frequency range. It is also well known from literature that higher order theories accuracy
deteriorates as the laminate becomes thicker. Nevertheless, dispersion knowledge gained with
this fast modelling approach is of great importance for NDT and SHM applications since it
plays a critical role in the selection of the optimal inspection frequencies for the improvement
of the sensitivity and for optimization of sensor networks in terms of sensor placement and
number of sensors.

Acknowledgment

The authors would like to express their gratitude to the Research School on Multi-Modal
Sensor Systems for Environmental Exploration (MOSES) and the Center for Sensor Systems
(ZESS) for sponsoring the research presented herein. Furthermore, the authors thank Dr.-Ing.
Rolf T. Schulte from the University of Siegen, Department of Technical Mechanics, for
fruitful discussions concerning the implementation proposed methodology.

References

1. Staszewski, W., C. Boller, and G.R. Tomlinson, Health Monitoring of Aerospace

Structures: Smart Sensor Technologies and Signal Processing. Munich: Wiley, 2004.
2. Kundu, T., Ultrasonic Nondestructive Evaluation: Engineering and Biological Material
Characterization. 1st ed. London: CRC Press, 2003.
3. Pavlakovic, B., et al., Disperse: A General Purpose Program for Creating Dispersion
Curves. Review of Progress in QNDE, 1997. 16: pp. 185-192.
4. Bartoli, I., et al., Modeling guided wave propagation for structural monitoring of damped
waveguides, in Third European Workshop on Structural Health Monitoring. 2006:
Granada, Spain.
5. Reddy, J.N., Mechanics of Laminated Composite Plates and Shells: Theory and Analysis.
2nd ed. New York: CRC Press, 2004.
6. Calomfirescu, M., Lamb Waves for Structural Health Monitoring in Viscoelastic
Composite Materials. PhD Thesis. Faserinstitut Bremen, University of Bremen, 2008.
7. Reddy, J.N., A refined nonlinear theory of plates with transverse shear deformation.
International Journal of Solids and Structures, 1984. 20(9-10): pp. 881-896.
8. Mindlin, R.D., Influence of rotary inertia and shear on flexural motion of isotropic
elastic plates. Journal of Applied Mechanics Transactions of the ASME, 1951. 18: pp.
31-38.
9. Whitney, J.M. and C.T. Sun, A higher order theory for extensional motion of laminated
composites Journal of Sound and Vibration, 1973. 30(1): pp. 85-97.
10. Wang, L. and F.G. Yuan, Lamb wave propagation in composite laminates using a higher-
order plate theory. Vol. 6531. Bellingham, WA, ETATS-UNIS: Society of Photo-Optical
Instrumentation Engineers, 2007.
11. Moll, J., M.A. Torres Arredondo, and C.P. Fritzen, Numerical Aspects of Automatic
Guided Wave Based Damage Localization Algorithms in Flat Anisotropic Structures.
Submitted for publication to "Smart Structures and Systems".
12. Deschamps, M. and B. Hosten, The effects of viscoelasticity on the reflection and
transmission of ultrasonic waves by an orthotropic plate. The Journal of the Acoustical
Society of America, 1992. 91(4): pp. 2007-2015.
13. Bartoli, I., et al., Modeling wave propagation in damped waveguides of arbitrary cross-
section. Journal of Sound and Vibration, 2006. 295(3-5): pp. 685-707
14. Neau, G., Lamb waves in anisotropic viscoelastic plates. Study of the wave fronts and
attenuation. PhD Thesis. L’Universite de Bordeaux, 2003.
Appendix

For the case of symmetric laminates the equations of motion can be decoupled into two
independent matrices. In order to obtain solutions for the symmetric and antisymmetric modes
of propagation, combinations of frequency and wavenumber where the matrices determinants
go to zero must be found. The symmetric modes are governed by

 L11
s L12
s L13
s L14
s L15
s

 12 
 Ls L22
s L23
s L15
s L25
s 
 
L s χ s =  L13
s L23
s L33
s L34
s L35
s  s
χ (A.1)
 L14 L15 L34 L44 L45 
 s s s s s 
 L15
s L25
s L35
s L45
s L55
s 


where χs = {U0 V0 Ѱz Φx Φy} and the terms of Ls are given by

s = A11k x + A66 k y + 2 A16 k x k y − ω I 0

L11 2 2 2

s = A16 k x + A26 k y + ( A12 + A66 ) k x k y

L12 2 2

s = −iκ 3 ( A13 k x + A36 k y )

L13

s = D11k x + D66 k y + 2 D16 k x k y − ω I1

L14 2 2 2

s = D16 k x + D26 k y + ( D12 + D66 ) k x k y

L15 2 2

s = A66 k x + A22 k y + 2 A26 k x k y − ω I 0

L22 2 2 2

s = −iκ 3 ( A36 k x + A23 k y )

L23 (A.2)
L25
s = D66 k x2 + D22 k y2 + 2 D26 k x k y − ω I 22

s = −κ 4 D55 k x − D44κ 5 k y − 2κ 4κ 5 D45 k x k y − κ 3 A33 + ω I 2

L33 2 2 2 2 2 2

s = i (2κ 4 D55 − κ 3 D13 ) k x + i (2κ 4κ 5 D45 − κ 3 D36 ) k y

L34 2

s = i (2κ 4κ 5 D45 − κ 3 D36 ) k x + i (2κ 5 D44 − κ 3 D23 ) k y

L35 2

s = H11k x + H 66 k y + 2 H16 k x k y + 4κ 4 D55 − ω I 4

L44 2 2 2 2

s = H16 k x + H 26 k y + ( H12 + H 66 ) k x k y + 4κ 4κ 5 D45

L45 2 2

s = H 66 k x + H 22 k y + 2 H 26 k x k y + 4κ 5 D44 − ω I 4
L55 2 2 2 2

The antisymmetric modes are governed by

 L11
A L12
A L13
A L14
A L15
A L16
A

 12 
 LA L22
A L23
A L24
A L25
A L26
A 
 13 
LA L23 L33 L34 L35 L36
L A χ A =  14 A A A A A 
χ (A.3)
L L24 L34 L44 L45 L46  A
 A A A A A A 
 L15
A L25
A L35
A L45
A L55
A L56
A

 16 
 LA L26
A L36
A L46
A L56
A L66
A 
where χA = {W0 Ѱx Ѱy Φz Λx Λy} and the terms of LA are given by

A = −(κ1 A55 k x + κ 2 A44 k y + 2κ1κ 2 A45 k x k y − ω I 0 )

L11 2 2 2 2 2

A = iκ1 (κ1 A55 k x + κ 2 A45 k y )

L12

A = iκ 2 (κ1 A45 k x + κ 2 A44 k y )

L13

A = −(κ1κ 7 D55 k x + κ 2κ 8 D44 k y + (κ 2κ 7 + κ1κ 8 ) D45 k x k y − ω I 2 )

L14 2 2 2

A = 3i (κ1κ 7 D55 k x + κ 2κ 7 D45 k y )

L15

A = 3i (κ1κ 8 D45 k x + κ 2κ 8 D44 k y )

L16
A = D11k x + D66 k y + 2 D16 k x k y + κ1 A55 − ω I 2
L22 2 2 2 2

A = D16 k x + D26 k y + ( D12 + D66 ) k x k y + κ1κ 2 A45

L23 2 2

A = i (κ1κ 7 D55 k x + κ1κ 8 D45 k y + 2κ 6 ( D13k x + D36 k y ))

L24

A = H11k x + H 66 k y + 2 H16 k x k y + 3κ1κ 7 D55 − ω I 4

L25 2 2 2

(A.4)
A = H16 k x + H 26 k y + ( H12 + H 66 ) k x k y + 3κ1κ 8 D45
L26 2 2

A = D66 k x + D22 k y + 2 D26 k x k y + κ 2 A44 − ω I 2

L33 2 2 2 2

A = i (κ 2κ 7 D45 k x + κ 2κ 8 D44 k y − 2κ 6 ( D36 k x + D23k y ))

L34

A = H16 k x + H 26 k y + ( H12 + H 66 ) k x k y + 3κ 2κ 7 D45

L35 2 2

A = H 66 k x + H 22 k y + 2 H 26 k x k y + 3κ 2κ 8 D44 − ω I 4
L36 2 2 2

A = −(κ 7 H 55 k x + κ 8 H 44 k y + (κ 7 + κ 8 ) H 45 k x k y + 4κ 6 D33 − ω I 4 )
L44 2 2 2 2 2 2

A = 3i (κ 7 H 55 k x + κ 7κ 8 H 45 k y ) − 2iκ 6 ( H13 k x + H 36 k y )
L45 2

A = 3i (κ 7κ 8 H 45 k x + κ 8 H 44 k y ) − 2iκ 6 ( H 36 k x + H 23k y )
L46 2

A = K11k x + K 66 k y + 2 K16 k x k y + 9κ 7 H 55 − ω I 6
L55 2 2 2 2

A = K16 k x + K 26 k y + ( K12 + K 66 ) k x k y + 9κ 7κ 8 H 45
L56 2 2

A = K 66 k x + K 22 k y + 2 K 26 k x k y + 9κ 8 H 44 − ω I 6
L66 2 2 2 2

where kx= k cos(θ), ky= k sin(θ) and k is the wavenumber. Additionally,

h /2
( Aij , Dij , H ij , Kij ) = ∫ Cɶij (1, z 2 , z 4 , z 6 )dz (A.5)
− h /2
h /2
(I0 , I2 , I4 , I6 ) = ∫ ρ (1, z 2 , z 4 , z 6 )dz (A.6)
− h /2

where i, j= 1,...,6, ρ is the density of the material and h is the thickness of the plate. The shear
correction coefficients κ1 to κ5 where taken according to [9]. The remaining correction factors
were calculated by matching specific frequencies from the approximate theory to frequencies
obtained from the exact elasticity theory as follows: κ6 = π/ 15 and κ7 = κ8 = π/ 17 . The
resulting complex wavenumber k = kRe + ikIm is used to describe the phase velocity of waves
travelling through their real part, kRe, and the amplitude decay through their imaginary part,
kIm.