Classification of solutions for guided waves in fluid-loaded viscoelastic composites with large numbers of layers

Huber, Armin M. A.

doi:10.1121/10.0020584

Guided ultrasonic waves are used for the inspection of multilayered composite aerospace structures. Calculating the corresponding dispersion diagrams is challenging for thick-walled composites with more than 100 layers, such as in modern rocket booster pressure vessels. The Dispersion Calculator (DC) is an open source software for calculating such dispersion diagrams and mode shapes of guided waves. Attenuation caused by viscoelasticity and fluid-loading makes the dispersion curve tracing much more difficult than in the nonattenuated case because the modal solutions are sought in the complex wavenumber plane. The tracing problem is mastered by a reliable algorithm. Whereas leaky Lamb and Scholte waves in coupled and decoupled cases are modeled using the stiffness matrix method, shear horizontal (SH) waves are traced using the transfer matrix method without facing the numerical instability. Through implementation of mode family specific dispersion equations in both matrix techniques for nonattenuated and attenuated cases, symmetric, antisymmetric, and nonsymmetric leaky Lamb, Scholte, and SH waves can be traced separately with better efficiency and robustness. The capabilities of DC are demonstrated by calculating dispersion diagrams and mode shapes for a viscoelastic composite with 400 layers immersed in water. Results are compared against DISPERSE (Imperial College London, London, UK) for selected cases.

I. INTRODUCTION

Ultrasonic guided waves play a key role in the nondestructive inspection (NDI) or structural health monitoring (SHM) of single and multilayered plates and pipes.^1–8 The significant complication inherent to guided waves is, however, their dispersive nature. Lamb waves, as a prominent example, represent a case of resonance in the waveguide, where some fundamental modes and a theoretically unlimited number of higher order modes can exist. To select the proper mode and frequency for a given inspection task, a dispersion diagram is essential. Calculating such diagrams is the topic of ongoing research. Historically, isotropic plates were the first specimens considered for NDI,^9,10 and thanks to the fundamental work of Lamb¹¹ in 1917, dispersion diagrams of the guided waves named after him can be calculated in the most efficient way since then.¹² However, evolution of materials and manufacturing led to multilayered anisotropic structures, culminating nowadays in lightweight composite laminates used in aircraft fuselages, rocket booster pressure vessels, or liquid hydrogen tanks. As the complexity of structures increased, new numerical methods for predicting guided wave propagation were required. In general, these methods can be classified in two groups. The first group, to which the Rayleigh-Lamb equations also belong, are the partial wave root-finding (PWRF) methods. Analytic dispersion equations are solved numerically, yielding exact solutions of wavenumber versus frequency. By applying a root-finding algorithm like Newton-Raphson, Muller's method, bisection, etc., the desired accuracy of the solutions can be maintained unconditionally. The second and newer groups are the discretizing methods. As in the well-known finite element method (FEM), the waveguide is discretized by several nodes through the thickness at which the differential equations are evaluated. The difference to the full FEM is that the direction of wave propagation is often solved analytically, which gives the partially discretizing method a better efficiency. The accuracy of the solutions depends on the number of through-thickness nodes. Whereas PWRF methods are limited to flat plates and cylindrical structures, the discretizing methods allow for modeling of waveguides with complex cross sections.

The first PWRF method for predicting guided wave propagation in a layered flat plate was the transfer matrix method (TMM), developed by Thomson¹³ in 1950 and edited by Haskell¹⁴ in 1953. In TMM, the transfer matrix relates the displacements and stresses at the top of a layer to those at its bottom. Through the multiplication of the individual layer transfer matrices, assuming the continuity of stresses and displacements across the layer boundaries, a global transfer matrix is obtained, relating the displacements and stresses at the top of the whole structure to those at its bottom. Modal solutions are obtained through imposing the boundary condition of stress-free top and bottom surfaces. Nayfeh and Chimenti^15,16 and Nayfeh^17,18 extended TMM to generally anisotropic multilayered media. This work is well-documented in the excellent book by Nayfeh.¹⁹ However, Dunkin²⁰ discovered in 1965 that TMM becomes unstable when the product of the layer thickness and frequency is large.

A different approach is the global matrix method (GMM), originally proposed for isotropic media by Knopoff²¹ in 1964. A well-organized and instructive discussion of TMM and GMM is provided by Lowe.²² In GMM, a large single matrix comprises the equations of all layers. Schmidt and Tango²³ and Schmidt and Jensen^24,25 showed that this method is numerically stable if implemented properly. Mal^26,27 extended GMM to anisotropic multilayered media. Kausel and Roesset²⁸ conducted a reformulation for isotropic cases and Wang and Rajapakse²⁹ conducted a reformulation for an orthotropic layer in the plane of symmetry. Wang and Rokhlin^30,31 used this approach for the computation of the cell stiffness matrix of a cross-ply composite. Another stable formulation for isotropic media was introduced by Kennett³² and Kennett and Kerry,³³ which was subsequently extended to generally anisotropic media by Fryer and Frazer^34,35 and obtained later by Rokhlin and Huang.³⁶ In 2014, Pant et al.³⁷ applied GMM to a multilayered composite laminate, and in 2023, Guo et al.³⁸ proposed the damped global matrix method (dGMM) to predict the frequency and propagation direction-dependent attenuation in viscoelastic composite laminates. dGMM extends the conventional damping-free GMM by incorporating viscoelastic damping models. However, the drawback of GMM is that it becomes slow on calculating many layers, resulting from the corresponding large size of the global matrix.

Probably, the most powerful PWRF method for multilayered problems is the stiffness matrix method (SMM), developed by Rokhlin and Wang^39,40 in 2001, and well-documented in the book of Rokhlin et al.⁴¹ Rokhlin and Wang rearranged the transfer matrix formulation in such a way that the numerical instability for large frequency-thickness products is avoided. In contrast to the transfer matrix, the stiffness matrix relates the stresses at the top and bottom of a layer (or the whole structure) to the displacements at the top and bottom of the layer/structure. At the same time, SMM formulation retains the concise form and efficiency of TMM. Some extensions to SMM were made by Tan,⁴² most notably the hybrid compliance-stiffness matrix method (HCSMM).⁴³ SMM was used by researchers like Kamal and Giurgiutiu,⁴⁴ Barski and Pajak,^45,46 Barski and Muc,⁴⁷ and Muc et al.⁴⁸ for the calculation of guided wave dispersion curves for multilayered composites. Huber and Sause^49,50 used SMM in 2018 to calculate dispersion diagrams for thick-walled composites containing up to 400 layers.

The most often applied discretizing method is the semi-analytical finite element (SAFE) method. It was introduced by Gavric⁵¹ in 1995 for the calculation of guided wave dispersion curves for a free rail. The SAFE approach has been used for many different scenarios. Hayashi et al.⁵² calculated guided wave dispersion curves for bar, rod, and rail examples. Leaky waves propagating along waveguides of arbitrary cross sections surrounded by an infinite medium were considered by Castaings and Lowe.⁵³ Torsional waves propagating along arbitrary cross sectional waveguides immersed in a fluid⁵⁴ and elastic waves propagating along weld joints between plates⁵⁵ were modeled by Fan et al. Bartoli et al. and Marzani et al. investigated guided waves in viscoelastic waveguides of rectangular, arbitrary,⁵⁶ and axissymmetric⁵⁷ cross sections. Another study of dispersion curve calculation for complex waveguides was conducted by Sorohan et al.⁵⁸ SAFE has also reached the field of clinical research as Pereira et al.⁵⁹ investigated the diagnosis of osteoporosis through guided wave propagation in cortical bone. More recent studies elaborated on the efficiency of SAFE, such as Seyfaddini et al.,⁶⁰ who proposed a semi-analytical isogeometric analysis (SAIGA), and Treysséde,⁶¹ who introduced a reduction strategy.

The local interaction simulation approach (LISA), a method studied by Delsanto et al.⁶² for a three-dimensional case in 1997, discretizes the structure into a lattice. The advantage of LISA is that discontinuities in the material's properties can be modeled by changing the properties of the lattice at the corresponding locations. LISA was applied for the modeling of isotropic as well as anisotropic specimens.^63,64

The possibly most powerful discretizing method is the spectral collocation method (SCM), introduced by Adamou and Craster⁶⁵ for guided wave modeling in elastic media in 2004. SCM is similar to the SAFE approach in that it uses a one-dimensional mesh over the system's cross section, but SCM possesses a higher accuracy and speed of computation. Instead of solving a differential equation directly, SCM uses a spectral approximation, which satisfies the differential equation and boundary conditions. Karpfinger et al. applied SCM for the modeling of wave dispersion along multilayered isotropic cylindrical structures⁶⁶ and, later, for axisymmetric wave modes in a porous elastic cylinder⁶⁷ as well as for a geophysical application involving boreholes.^68,69 Yu et al. combined SCM with PWRF methods for the modeling of longitudinal guided waves in a multilayered tube.⁷⁰ Zharnikov et al.⁷¹ used SCM to obtain dispersion curves for anisotropic waveguides. In 2015, Quintanilla et al. provided comprehensive studies of guided wave modeling in generally anisotropic media by means of SCM,^72,73 later adding three-dimensional dispersion curve solutions.⁷⁴ These publications claim that SCM is easier to code than PWRF methods, it is faster, and, most importantly, it does not miss modal solutions. In 2017, Quintanilla et al. developed a classification of multilayered anisotropic waveguides according to their symmetry and coupling properties.⁷⁵ Five classes were determined wherein the separation into symmetric and antisymmetric modes, as well as pure Lamb and shear horizontal (SH) waves, was elucidated. A critical benefit is that dispersion curves belonging to the various mode families do not cross or do so less frequently. This helps reduce or completely avoid the well-known jumping mode problem, which is present in PWRF methods in some situations. Furthermore, one can save time by computing only a single mode family if one is not interested in the other mode families. More recently, the method was used also to model leaky waves.^76,77

While the computation of dispersion curves has been studied extensively for several decades, access to computer programs has been difficult or impossible for a long time. The only commercial software specifically designed for calculating dispersion curves is DISPERSE (Imperial College London, London, UK). It has been developed by Lowe and Pavlakovic^78,79 since the early 1990s, and it uses GMM. The first free software, GUIGUW, appeared in 2011. It was developed by Marzani and Bocchini,^80,81 based on SAFE method. The MATLAB^®-based (The MathWorks, Inc., Natick, MA) open source dispersion calculator (DC), released by Huber⁸² in 2018, uses SMM. Another contender, the open source MATLAB^® toolbox ElasticMatrix, was released by Ramasawmy et al.^83,84 in 2020. It uses GMM. Next, the Dispersion Box was released by Orta et al.^85,86 in 2022. Up to six different methods can be used simultaneously, namely, GMM, SMM, HCSMM, SAFE, the Legendre polynomial method (LPM), and the fifth order shear deformation theory (5-SDT). The latter was also developed by Orta et al.⁸⁷ as an improvement to the previous lower order plate theories. Also, in 2022, Kiefer⁸⁸ published GEW dispersion script for MATLAB^®, which uses SCM.

In this contribution, the implementation of SMM in DC is presented. Although discretizing methods are preferred lately,^77,89 the decision made in 2016 to implement SMM still proves justified today. The key requirement imposed on DC was that it should be able to calculate thick-walled composites of up to 400 layers as nowadays they are present in rocket booster pressure vessels.⁹⁰ To facilitate the calculation, it is a common practice to group layers, but this is not possible here because the layups are complicated and irregular in terms of the fiber orientation and layer thicknesses. Therefore, every single layer must be calculated. In the previous study in Ref. 49, it was demonstrated that DC fulfills this requirement. However, the composite was considered perfectly elastic and with stress-free surfaces, thereby circumventing the problem of attenuation. In this paper, the two important damping mechanisms are added, namely, viscoelasticity and energy leakage, due to fluid-loading. This poses a significant complication in the dispersion curve tracing routine because the root search must now be performed in the complex wavenumber plane rather than on the real wavenumber axis only. To improve efficiency and robustness of the dispersion curve tracing, mode family specific dispersion equations and boundary conditions were presented in Ref. 49. These are simplified where possible and extended to the attenuated cases. In particular, symmetric and antisymmetric dispersion equations are derived for Lamb waves as well as leaky Lamb and Scholte waves. These are computationally less expensive and make the previous symmetry check of each solution during the dispersion curve tracing obsolete. In case of wave propagation along principal axes, efficiency and robustness are further enhanced by using the decoupled equations. For SH waves, TMM can be used without encountering the numerical instability because the contributing bulk waves never become evanescent. The global transfer matrix is obtained by a simplified algorithm, and distinct equations for symmetric and antisymmetric modes are given. In summary, by classifying modal solutions according to their mode type as well as symmetry and coupling properties, symmetric, antisymmetric, and nonsymmetric leaky Lamb, Scholte, and SH waves are traced separately with improved speed and robustness.

The paper is organized as follows. In Sec. II, the numerical model used by DC is presented. The mode family specific dispersion equations are derived for the traction-free case and attenuated cases. Section III elaborates which mode families appear depending on the configuration of the waveguide and adjacent half-spaces, and then the appropriate dispersion equations are assigned. Section IV features some numerical examples. Two exemplary water-loaded viscoelastic composite layups containing 16 layers are calculated with DC and DISPERSE for validation. Then, a 400 layer composite is calculated with DC, serving as a showcase of its capabilities.

II. NUMERICAL MODEL

Guided waves form through the propagation and superposition of bulk waves. Therefore, in the context of guided wave modeling, bulk waves are often referred to as partial waves. All bulk waves in every layer must act in accordance with Snell's law to allow for a propagating guided wave, and the waveguide constitutes the geometry for which the case of resonance is sought. In this case, a guided wave propagates along the waveguide, also representing a resonant standing wave in the transverse direction. This phenomenon is often referred to as transverse resonance in the literature. The bulk waves must reflect forth and back at the appropriate wave vectors such that they undergo a phase shift of an integral multiple of 2π between each round trip and consequently interfere constructively. Therefore, it is of fundamental importance to find the propagation directions of the bulk waves. As boundaries and interfaces are involved in a waveguide, guided waves show dispersion, i.e., the phase velocity is frequency dependent. Thus, to obtain dispersion diagrams, the problem must be solved for extended frequency ranges.

A. Field equations and geometrical setup

Wave propagation in linearly elastic solids is described by two basic field equations and one elastic constitutive equation, which in the following is referred to as Hooke's law. Introducing the Cartesian crystallographic coordinate system,

x_{i}^{'} = (x_{1}^{'}, x_{2}^{'}, x_{3}^{'})

⁠, and using subscript notation,

i, j, k, l = 1, 2, 3

⁠, the strain-displacement relation, the equation of motion, and Hooke's law are, respectively, given by

ε_{k l}^{'} = \frac{1}{2} (\frac{\partial u_{l}^{'}}{\partial x_{k}^{'}} + \frac{\partial u_{k}^{'}}{\partial x_{l}^{'}}),

(1)

\frac{\partial σ_{i j}^{'}}{\partial x_{j}^{'}} = ρ \frac{\partial^{2} u_{i}^{'}}{\partial t^{2}},

(2)

σ_{i j}^{'} = c_{ijkl}^{'} ε_{k l}^{'},

(3)

where

σ_{i j}^{'}

and

ε_{k l}^{'}

are the stress and strain tensors, respectively,

u_{i}^{'}

is the displacement vector, and ρ is the material's mass density. The stiffness tensor,

c_{ijkl}^{'}

⁠, contains the elastic material constants. Through symmetry properties and strain energy density considerations, the 81 elements in

c_{ijkl}^{'}

are reduced to 21, describing a generally anisotropic (triclinic) material. In the present study, layups of unidirectional carbon fiber reinforced polymer (CFRP) layers are considered, where the single layer belongs to the orthotropic symmetry class. In this case, only nine independent elastic constants are remaining. By applying Voigt's notation to

c_{ijkl}^{'}

⁠, Eq. can be written as

[\begin{matrix} σ_{11}^{'} \\ σ_{22}^{'} \\ σ_{33}^{'} \\ σ_{23}^{'} \\ σ_{13}^{'} \\ σ_{12}^{'} \end{matrix}] = [\begin{matrix} C_{11}^{'} & C_{12}^{'} & C_{13}^{'} & 0 & 0 & 0 \\ C_{22}^{'} & C_{23}^{'} & 0 & 0 & 0 \\ C_{33}^{'} & 0 & 0 & 0 \\ C_{44}^{'} & 0 & 0 \\ sym & C_{55}^{'} & 0 \\ C_{66}^{'} \end{matrix}] [\begin{matrix} ε_{11}^{'} \\ ε_{22}^{'} \\ ε_{33}^{'} \\ 2 ε_{23}^{'} \\ 2 ε_{13}^{'} \\ 2 ε_{12}^{'} \end{matrix}] .

(4)

By using upper case

C'

's, the contracted matrix form is distinguished from the lower case

c'

's of the expanded tensor form.

Consider a multilayered composite such as that sketched in Fig. 1. Composites are stackings of m layers with layer thicknesses, d_m, consisting of a fiber-epoxy combination. For each layer, a local (crystallographic) coordinate system $x_{i, m}^{'} = {(x_{1}^{'}, x_{2}^{'}, x_{3}^{'})}_{m}$ is assigned, residing on the top of the mth layer such that the layers lie in the $x_{1, m}^{'}$ - $x_{2, m}^{'}$ -planes. The fibers are oriented along $x_{1, m}^{'}$ ⁠, and $x_{3, m}^{'}$ is normal to the layer. To describe a system of arbitrary layer orientations in a convenient way, the nonprimed global coordinate system, $x_{i} = (x_{1}, x_{2}, x_{3})$ ⁠, is introduced. It resides on the top layer of the composite, and it is oriented such that (bulk) wave propagation takes place in the x₁-x₃-plane. With respect to the global coordinate system, the local coordinate systems are yielded by a counterclockwise rotation through an angle, $Φ_{m}$ ⁠, between x₁ and $x_{1, m}^{'}$ about the x₃-axis. Therefore, the stiffness matrix, ${C^{'}}_{i j}$ ⁠, must be transformed from the local into the global coordinate system for each layer to get their respective transformed stiffness matrices, $C_{i j, m}$ ⁠, using the orthogonal transformation method.¹⁹ The resulting components of C_ij are given in Appendix A. Notice that C_ij has the form of monoclinic material symmetry, in general, having 13 independent elastic constants. For $Φ = 0 \circ, 90 \circ$ ⁠, the elements C₁₆, C₂₆, C₃₆, and C₄₅ vanish, returning orthotropic material symmetry. This will be of importance when the decoupling into pure Lamb and SH waves is studied.

FIG. 1.

View large Download slide

Single composite layer (left) with local (crystallographic) coordinate system $x_{i}^{'}$ and global coordinate system x_i, i = 1,2,3, and layered composite plate (right) with [0/90/0] orientation with respect to the $x_{i}^{'}$ -axis.

When an incident wave hits the top interface of the uppermost layer, up to three scattered bulk waves are generated in this layer, namely, using the notation introduced by Li and Thompson,⁹¹ one quasi-longitudinal wave, L, and two quasi-shear waves, SV and SH. These waves propagate downward (⁠

- x_{3}

⁠) at their respective refraction angles until they reach the bottom interface of the first layer. Here, a certain proportion of energy is transmitted into the second layer while the rest gets reflected, propagating upward (⁠

+ x_{3}

⁠) at the same refraction angles. Hence, a maximum of six bulk waves propagates in every layer, namely,

L^{-}

⁠,

L^{+}

⁠, SV⁻, SV⁺, SH⁻, and SH⁺. The prefix quasi indicates that the polarization vector does not point in the wave propagation direction for L-waves and is not normal to the propagation direction in shear waves, except for propagation along principal axes. In general, the polarization vector is oriented under a certain angle with respect to the wave vector. This angle depends on the orientation of the wave vector in an anisotropic solid. The sought propagation directions and polarization vectors are yielded as eigenvalues and eigenvectors of Christoffel's equation. This is a combination of Eqs. , which are given in the global coordinate system by

ρ \frac{\partial^{2} u_{i}}{\partial t^{2}} = c_{ijkl} \frac{\partial^{2} u_{l}}{\partial x_{j} \partial x_{k}} .

(5)

While the bulk waves in an infinite medium are unconditionally described by Eq. , the presence of interfaces leads to restrictions on the guided wave motion, accounted for by boundary conditions. First, rigid bonding between the layers is assumed, implying the continuity of stresses and displacements over the layer boundaries, namely,

u_{i} = {\bar{u}}_{i}, σ_{i 3} = {\bar{σ}}_{i 3}, i = 1, 2, 3,

(6)

where u_i,

σ_{i 3}

correspond to the originating layer, and

{\bar{u}}_{i}, {\bar{σ}}_{i 3}

correspond to the continuing layer. Second, if no fluid-loading is considered, the outer surfaces of the composite are traction-free such that the stress components at the top (⁠

x_{3} = 0

⁠) and bottom (⁠

x_{3} = - d

⁠) of the plate vanish, i.e.,

σ_{i 3}^{0, - d} = 0, i = 1, 2, 3.

(7)

B. The SMM

The task is to solve Eq. for the propagation directions and polarization vectors of the bulk waves for which all share the same wavenumber component, ξ (⁠

= ω / c_{p}

⁠) along x₁, as required by Snell's law. Here,

c_{p}

is the phase velocity in the x₁-direction, and ω is the angular frequency (⁠

= 2 π f

⁠). In contrast to the previous study in Ref. 49, where attenuation was not considered, the general solution of the displacement field of the form

u_{i} = U_{i} e^{i (ξ x_{1} + ξ_{3} x_{3} - ω t)}

(8)

is chosen, where one searches for the x₃-component of the wavenumber ξ₃ directly, instead of searching for the ratio of the bulk waves' wavenumber components in the x₃ and x₁-directions, α. This change is necessary since ξ turns complex once attenuation is introduced. The sought polarization vector is represented by the bulk wave amplitudes, U_i.

1. Coupled polarization

Considering the general case where the bulk waves contribute to the displacement field with nonzero polarization components in all three directions x_i, Eq. can be written out as

(u_{1}, u_{2}, u_{3}) = (U_{1}, U_{2}, U_{3}) e^{i ξ_{3} x_{3}},

(9)

where the common term

e^{i (ξ x_{1} - ω t)}

is implied. Substituting Eq. into Eq. in global coordinates yields three homogeneous linear equations,

[\begin{matrix} C_{11} ξ^{2} + C_{55} ξ_{3}^{2} - ρ ω^{2} & C_{16} ξ^{2} + C_{45} ξ_{3}^{2} & (C_{13} + C_{55}) ξ ξ_{3} \\ C_{66} ξ^{2} + C_{44} ξ_{3}^{2} - ρ ω^{2} & (C_{36} + C_{45}) ξ ξ_{3} \\ sym & C_{55} ξ^{2} + C_{33} ξ_{3}^{2} - ρ ω^{2} \end{matrix}] [\begin{matrix} U_{1} \\ U_{2} \\ U_{3} \end{matrix}] = 0,

(10)

where the 3 × 3 matrix is the Christoffel matrix, M_ij. By vanishing the determinant of M_ij, the sixth-degree bicubic polynomial equation,

ξ_{3}^{6} + A_{1} ξ_{3}^{4} + A_{2} ξ_{3}^{2} + A_{3} = 0,

(11)

is obtained with the coefficients A₁, A₂, and A₃ given in Appendix B. Equation has six solutions,

ξ_{3, q}, q = 1, 2, \dots, 6

⁠, namely, three pairs, corresponding to the upward and downward propagating bulk waves,

ξ_{3, 2} = - ξ_{3, 1}, ξ_{3, 4} = - ξ_{3, 3}, ξ_{3, 6} = - ξ_{3, 5} .

(12)

Substituting

ξ_{3, q}

into Eq. delivers the bulk waves' displacement amplitude ratios,

V_{q} = U_{2, q} / U_{1, q}

and

W_{q} = U_{3, q} / U_{1, q}

⁠, namely,

V_{q} = \frac{m_{11} (ξ_{3, q}) m_{23} (ξ_{3, q}) - m_{13} (ξ_{3, q}) m_{12} (ξ_{3, q})}{m_{13} (ξ_{3, q}) m_{22} (ξ_{3, q}) - m_{12} (ξ_{3, q}) m_{23} (ξ_{3, q})},

(13)

W_{q} = \frac{m_{11} (ξ_{3, q}) m_{22} (ξ_{3, q}) - m_{12} {(ξ_{3, q})}^{2}}{m_{12} (ξ_{3, q}) m_{23} (ξ_{3, q}) - m_{22} (ξ_{3, q}) m_{13} (ξ_{3, q})},

(14)

where

m_{i j} (ξ_{3, q})

are the elements of the Christoffel matrix. The stress amplitudes,

D_{i 3, q}

⁠, are obtained from Eq. as

\begin{matrix} D_{33, q} = i (C_{13} ξ + C_{36} V_{q} ξ + C_{33} W_{q} ξ_{3, q}), \\ D_{13, q} = i (C_{55} (ξ_{3, q} + W_{q} ξ) + C_{45} V_{q} ξ_{3, q}), \\ D_{23, q} = i (C_{45} (ξ_{3, q} + W_{q} ξ) + C_{44} V_{q} ξ_{3, q}) . \end{matrix}

(15)

Now, the displacement and stress field components can be written as superpositions of the six contributing bulk waves,

(u_{1}, u_{2}, u_{3}) = \sum_{q = 1}^{6} (1, V_{q}, W_{q}) U_{1, q} e^{i ξ_{3, q} x_{3}},

(16)

(σ_{33}, σ_{13}, σ_{23}) = \sum_{q = 1}^{6} (D_{33, q}, D_{13, q}, D_{23, q}) U_{1, q} e^{i ξ_{3, q} x_{3}} .

(17)

Considering an m-layered plate, Eqs. and can be written in matrix form, relating the displacement and stress components at the top

u_{m}, σ_{m}

(⁠

x_{3, m} = 0

⁠) and bottom

u_{m + 1}, σ_{m + 1}

(⁠

x_{3, m} = - d_{m}

⁠) of the mth layer to the wave amplitude vectors,

U_{m}^{\pm}

⁠,

[\begin{matrix} u_{m} \\ σ_{m} \end{matrix}] = {[\begin{matrix} P^{-} & P^{+} H \\ D^{-} & D^{+} H \end{matrix}]}_{m} [\begin{matrix} U_{m}^{-} \\ U_{m}^{+} \end{matrix}],

(18)

[\begin{matrix} u_{m + 1} \\ σ_{m + 1} \end{matrix}] = {[\begin{matrix} P^{-} H & P^{+} \\ D^{-} H & D^{+} \end{matrix}]}_{m} [\begin{matrix} U_{m}^{-} \\ U_{m}^{+} \end{matrix}],

(19)

where the tensor notation introduced by Rokhlin and Wang^39,40 is used. The submatrices,

P^{\pm}, D^{\pm}

⁠, H, and the vectors,

U^{\pm}

⁠, are expanded in Appendix C. By eliminating

U_{m}^{\pm}

in Eqs. and , the transfer matrix representation is obtained, relating the displacement and stress field components at the top of the mth layer to those at its bottom,

[\begin{matrix} u_{m + 1} \\ σ_{m + 1} \end{matrix}] = {[\begin{matrix} P^{-} H & P^{+} \\ D^{-} H & D^{+} \end{matrix}]}_{m} {[\begin{matrix} P^{-} & P^{+} H \\ D^{-} & D^{+} H \end{matrix}]}_{m}^{- 1} [\begin{matrix} u_{m} \\ σ_{m} \end{matrix}] = A_{m} [\begin{matrix} u_{m} \\ σ_{m} \end{matrix}],

(20)

where A_m is the local transfer matrix of the mth layer. The numerical instability of TMM occurs when contributing bulk waves become evanescent at high frequency-thickness products. Then,

ξ_{3, q}

becomes imaginary, leading to decaying functions,

e^{- ξ_{3, q} x_{3}}

⁠. These fall below the precision of a computer quickly such that H contains only zeros, eventually. Hence, entire columns in the two matrices in Eq. from which

A_{m}

is calculated turn zero, rendering the matrices singular. The numerical instability comes from the fact that the inverse of a singular matrix is undefined.

To overcome the problem, Rokhlin and Wang rearranged Eqs. and in such a way that the exponential terms are absent from the diagonal elements. The stiffness matrix representation relates the stresses at the top,

σ_{m}

⁠, and bottom,

σ_{m + 1}

⁠, of the mth layer to the displacements at the top,

u_{m}

⁠, and bottom,

u_{m + 1}

⁠, via the local stiffness matrix,

K_{m}

⁠,

[\begin{matrix} σ_{m} \\ σ_{m + 1} \end{matrix}] = {[\begin{matrix} D^{-} & D^{+} H \\ D^{-} H & D^{+} \end{matrix}]}_{m} {[\begin{matrix} P^{-} & P^{+} H \\ P^{-} H & P^{+} \end{matrix}]}_{m}^{- 1} [\begin{matrix} u_{m} \\ u_{m + 1} \end{matrix}] = K_{m} [\begin{matrix} u_{m} \\ u_{m + 1} \end{matrix}] .

(21)

To obtain the global stiffness matrix of an m-layered plate, the individual layer stiffness matrices,

K_{m}

⁠, must be calculated first. Considering the stiffness matrices of two neighboring layers,

[\begin{matrix} σ_{1} \\ σ_{2} \end{matrix}] = [\begin{matrix} K_{11}^{A} & K_{12}^{A} \\ K_{21}^{A} & K_{22}^{A} \end{matrix}] [\begin{matrix} u_{1} \\ u_{2} \end{matrix}], [\begin{matrix} σ_{2} \\ σ_{3} \end{matrix}] = [\begin{matrix} K_{11}^{B} & K_{12}^{B} \\ K_{21}^{B} & K_{22}^{B} \end{matrix}] [\begin{matrix} u_{2} \\ u_{3} \end{matrix}],

(22)

Rokhlin and Wang's recursive algorithm is used to combine them,

[\begin{matrix} σ_{1} \\ σ_{3} \end{matrix}] = [\begin{matrix} K_{11}^{A} + K_{12}^{A} X K_{21}^{A} & - K_{12}^{A} X K_{12}^{B} \\ K_{21}^{B} X K_{21}^{A} & K_{22}^{B} - K_{21}^{B} X K_{12}^{B} \end{matrix}] [\begin{matrix} u_{1} \\ u_{3} \end{matrix}],

(23)

with

X = {(K_{11}^{B} - K_{22}^{A})}^{- 1}

⁠. Calling the obtained matrix

K^{A}

and the stiffness matrix of the third layer

K^{B}

⁠, one can recursively use Eq. to obtain the global stiffness matrix, relating the stresses and displacements at the top and bottom of the whole plate. If the plate consists of periodic repetitions of a unit cell, as the case in composites, one obtains the unit cell stiffness matrix and then uses Eq. with

K^{B} = K^{A}

as many times as repetitions are contained in the composite. Finally, if the layup is symmetric, the so far obtained stiffness matrix, renamed to

K^{A}

⁠, represents one-half of the layup. As an alternative formalism to that presented by Rokhlin and Wang, the global stiffness matrix, K, of the whole, symmetric laminate is given by

[\begin{matrix} σ_{1} \\ σ_{3} \end{matrix}] = [\begin{matrix} K_{11}^{A} + K_{12}^{A} Y K_{21}^{A} & - K_{12}^{A} Y (K_{21}^{A} ° I) \\ (K_{12}^{A} ° I) Y K_{21}^{A} & K_{11}^{A} ° I - (K_{12}^{A} ° I) Y (K_{21}^{A} ° I) \end{matrix}] [\begin{matrix} u_{1} \\ u_{3} \end{matrix}],

(24)

with

Y = {(K_{22}^{A} ° I - K_{22}^{A})}^{- 1}

⁠, the Hadamard operator “

°

⁠,” and

I = [\begin{matrix} 1 & 1 & - 1 \\ - 1 & - 1 & 1 \\ - 1 & - 1 & 1 \end{matrix}] .

(25)

Finally, by invoking the stress-free boundary conditions from Eq. ,

[\begin{matrix} σ^{0} \\ σ^{- d} \end{matrix}] = K [\begin{matrix} u^{0} \\ u^{- d} \end{matrix}] = 0,

(26)

the dispersion equation is found as

det K = 0.

(27)

However, to improve efficiency and robustness of dispersion curve tracing, it is convenient to use mode family specific dispersion equations. In fact, Eq. is used only for the nonsymmetric modes propagating in nonsymmetric layups. By contrast, if the layup is symmetric or a single layer, one can derive distinct dispersion equations for symmetric and antisymmetric modes by applying corresponding boundary conditions. For both mode families, the stress-free upper surface remains as before, but symmetry specific boundary conditions are invoked at the half thickness of the layup. Looking at the mode shapes, it turns out that for symmetric modes, the following components of stress and displacement must vanish at the half thickness:

σ_{13}^{- d / 2}, σ_{23}^{- d / 2}, u_{3}^{- d / 2} = 0,

(28)

and for antisymmetric modes,

σ_{33}^{- d / 2}, u_{1}^{- d / 2}, u_{2}^{- d / 2} = 0.

(29)

Hence, it is sufficient to compute only the stiffness matrix of the upper half of the symmetric layup,

K^{half}

⁠, relating the stresses at the top,

σ^{0}

⁠, and the half thickness,

σ^{- d / 2}

⁠, of the composite to the corresponding displacements,

u^{0}

and

u^{- d / 2}

⁠. With the stress-free upper surface, one can write the expanded form as

[\begin{matrix} 0 \\ 0 \\ 0 \\ σ_{33}^{- d / 2} \\ σ_{13}^{- d / 2} \\ σ_{23}^{- d / 2} \end{matrix}] = {[\begin{matrix} k_{11} & k_{12} & k_{13} & k_{14} & k_{15} & k_{16} \\ k_{21} & k_{22} & k_{23} & k_{24} & k_{25} & k_{26} \\ k_{31} & k_{32} & k_{33} & k_{34} & k_{35} & k_{36} \\ k_{41} & k_{42} & k_{43} & k_{44} & k_{45} & k_{46} \\ k_{51} & k_{52} & k_{53} & k_{54} & k_{55} & k_{56} \\ k_{61} & k_{62} & k_{63} & k_{64} & k_{65} & k_{66} \end{matrix}]}^{half} [\begin{matrix} u_{1}^{0} \\ u_{2}^{0} \\ u_{3}^{0} \\ u_{1}^{- d / 2} \\ u_{2}^{- d / 2} \\ u_{3}^{- d / 2} \end{matrix}] .

(30)

Then, by using the symmetric boundary conditions, Eq. , on Eq. , the dispersion equation for symmetric modes is obtained as

det K_{S} = det {[\begin{matrix} k_{11} & k_{12} & k_{13} & k_{14} & k_{15} \\ k_{21} & k_{22} & k_{23} & k_{24} & k_{25} \\ k_{31} & k_{32} & k_{33} & k_{34} & k_{35} \\ k_{51} & k_{52} & k_{53} & k_{54} & k_{55} \\ k_{61} & k_{62} & k_{63} & k_{64} & k_{65} \end{matrix}]}^{half} = 0,

(31)

and by using the antisymmetric boundary conditions, Eq. , on Eq. , the dispersion equation for antisymmetric modes reads

det K_{A} = det {[\begin{matrix} k_{11} & k_{12} & k_{13} & k_{16} \\ k_{21} & k_{22} & k_{23} & k_{26} \\ k_{31} & k_{32} & k_{33} & k_{36} \\ k_{41} & k_{42} & k_{43} & k_{46} \end{matrix}]}^{half} = 0.

(32)

The benefit of Eqs. and in terms of efficiency is that Eq. can be skipped, and only a 5 × 5 or 4 × 4 matrix determinant must be computed instead of a 6 × 6 matrix determinant, as required in Eq. . The advantage in terms of robustness is that one can trace symmetric and antisymmetric modes separately without the risk of jumping to a mode of the other mode family. Equations and are used in DC since version 2.4. In earlier versions, symmetric and antisymmetric modes were obtained using Eq. and distinguished by calculating the in-plane displacement at the top and bottom surfaces during dispersion curve tracing, as described in Ref. 49. This displacement calculation routine was removed in version 2.4, thereby enhancing speed and robustness of dispersion curve tracing even more.

2. Decoupled polarization

In case of wave propagation along principal axes, polarization in the sagittal plane (x₁-x₃) decouples from polarization in the shear horizontal direction x₂. This situation is possible only in single layers and cross-ply laminates for wave propagation along or normal to the fibers. For the first type, (pure) Lamb waves are formed through the superposition of four bulk waves per layer, namely, $L^{-}$ ⁠, $L^{+}$ ⁠, and either SV⁻ and SV⁺ or SH⁻ and SH⁺. The second type are guided SH waves, evolving through the superposition of only two bulk waves per layer, namely, that pair of shear waves, which does not contribute to the Lamb waves in the respective layer.

a. Lamb waves.

Because the bulk waves contributing to a pure Lamb wave have nonzero polarization components only in the x₁-x₃-plane, Eq. reduces to

(u_{1}, u_{3}) = (U_{1}, U_{3}) e^{i ξ_{3} x_{3}},

(33)

and Eq. reduces to

[\begin{matrix} C_{11} ξ^{2} + C_{55} ξ_{3}^{2} - ρ ω^{2} & (C_{13} + C_{55}) ξ ξ_{3} \\ sym & C_{55} ξ^{2} + C_{33} ξ_{3}^{2} - ρ ω^{2} \end{matrix}] [\begin{matrix} U_{1} \\ U_{3} \end{matrix}] = 0,

(34)

where the 2 × 2 matrix is the Christoffel matrix M_ij. By vanishing its determinant, the fourth-degree bicubic polynomial equation,

A_{1} ξ_{3}^{4} + A_{2} ξ_{3}^{2} + A_{3} = 0,

(35)

with the coefficients,

\begin{matrix} A_{1} = C_{33} C_{55}, \\ A_{2} = (C_{11} C_{33} - 2 C_{13} C_{55} - C_{13}^{2}) ξ^{2} - (C_{33} + C_{55}) ρ ω^{2}, \\ A_{3} = C_{11} C_{55} ξ^{4} - (C_{11} + C_{55}) ξ^{2} ρ ω^{2} + ρ^{2} ω^{4}, \end{matrix}

(36)

is obtained. Equation admits four solutions,

ξ_{3, q}, q = 1, 2, 3, 4

⁠, given by

ξ_{3, q} = \pm \sqrt{\frac{- A_{2} \pm \sqrt{A_{2}^{2} - 4 A_{1} A_{3}}}{2 A_{1}}},

(37)

corresponding to the upward and downward propagating bulk waves,

ξ_{3, 2} = - ξ_{3, 1}, ξ_{3, 4} = - ξ_{3, 3} .

(38)

Substituting

ξ_{3, q}

into Eq. delivers the amplitude ratios,

W_{q} = U_{3, q} / U_{1, q}

⁠, namely,

W_{q} = \frac{- m_{11} (ξ_{3, q})}{m_{12} (ξ_{3, q})} = \frac{ρ ω^{2} - C_{11} ξ^{2} - C_{55} ξ_{3, q}^{2}}{(C_{13} + C_{55}) ξ ξ_{3, q}},

(39)

where

m_{i j} (ξ_{3, q})

are the elements of the Christoffel matrix. The stress amplitudes are

\begin{matrix} D_{33, q} = i (C_{13} ξ + C_{33} W_{q} ξ_{3, q}), \\ D_{13, q} = i C_{55} (ξ_{3, q} + W_{q} ξ) . \end{matrix}

(40)

By superposing the bulk waves, the displacement and stress fields can be written as

(u_{1}, u_{3}) = \sum_{q = 1}^{4} (1, W_{q}) U_{1, q} e^{i ξ_{3, q} x_{3}},

(41)

(σ_{33}, σ_{13}) = \sum_{q = 1}^{4} (D_{33, q}, D_{13, q}) U_{1, q} e^{i ξ_{3, q} x_{3}} .

(42)

The assembly of the layer stiffness matrices, obtaining the global stiffness matrix and modal solutions follow Eqs. in the reduced form with

I = [\begin{matrix} 1 & - 1 \\ - 1 & 1 \end{matrix}] .

(43)

The dispersion equation for symmetric modes is given by

det K_{S} = det {[\begin{matrix} k_{11} & k_{12} & k_{13} \\ k_{21} & k_{22} & k_{23} \\ k_{41} & k_{42} & k_{43} \end{matrix}]}^{half} = 0,

(44)

and that for antisymmetric modes is given by

det K_{A} = det {[\begin{matrix} k_{11} & k_{12} & k_{14} \\ k_{21} & k_{22} & k_{24} \\ k_{31} & k_{32} & k_{34} \end{matrix}]}^{half} = 0.

(45)

b. Shear horizontal waves.

SH waves consist of two shear bulk waves, which possess polarization only in the x₂-direction such that the displacement and stress fields are given by

u_{2} = \sum_{q = 1}^{2} U_{2, q} e^{i ξ_{3, q} x_{3}},

(46)

σ_{23} = \sum_{q = 1}^{2} D_{23, q} U_{2, q} e^{i ξ_{3, q} x_{3}},

(47)

with

ξ_{3, q} = \pm \sqrt{\frac{ρ ω^{2} - C_{66} ξ^{2}}{C_{44}}}, D_{23, q} = i C_{44} ξ_{3, q} .

(48)

As the bulk waves in guided SH waves never become evanescent, the numerical instability does not occur in the transfer matrix. For the relatively simple case of SH waves, Eq. can be manipulated into the convenient form¹⁹

{[\begin{matrix} u_{2} \\ σ_{23} \end{matrix}]}_{m + 1} = [\begin{matrix} cos γ_{m} & i / D_{m} sin γ_{m} \\ i D_{m} sin γ_{m} & cos γ_{m} \end{matrix}] {[\begin{matrix} u_{2} \\ σ_{23} \end{matrix}]}_{m} = A_{m} {[\begin{matrix} u_{2} \\ σ_{23} \end{matrix}]}_{m},

(49)

where

γ_{m} = ξ_{3, 1, m} d_{m}

and

D_{m} = D_{23, 1, m}

⁠. The layer transfer matrix, A_m in Eq. , contains the matrix product required in Eq. inherently and, therefore, allows highly efficient calculation.

In the previous study in Ref. 49, mode family specific dispersion equations had been introduced, allowing to distinguish symmetric, antisymmetric, and nonsymmetric SH waves. However, the algorithm used to obtain the global transfer matrix and dispersion equations depended on many layup parameters, leading to many different, and sometimes complicated, dispersion equations which had to be applied to specific cases. For antisymmetric SH waves in multilayered symmetric laminates, no transfer matrix style solution had been found such that a stiffness matrix style solution had to be used. In the following, a simplified and unified algorithm for calculating the global transfer matrix as well as more convenient dispersion equations are presented. All solutions are obtained in transfer matrix style.

In the original formulation introduced by Nayfeh in Ref. 19, the global transfer matrix, A, was obtained by consecutively multiplying all A_m, starting with the lowermost layer. At least in the case of SH waves, however, one can also start with the top layer, namely,

A = A_{1} A_{2} \dots A_{m},

(50)

where A relates the displacements and stresses at the top of the whole plate to those at its bottom, according to

{[\begin{matrix} u_{2} \\ σ_{23} \end{matrix}]}_{- d} = A {[\begin{matrix} u_{2} \\ σ_{23} \end{matrix}]}_{0} = [\begin{matrix} a_{11} & a_{12} \\ a_{21} & a_{22} \end{matrix}] {[\begin{matrix} u_{2} \\ σ_{23} \end{matrix}]}_{0} .

(51)

Consider two exemplary composite layups with the configurations [0/90]₄ and [0/90]

_{4 s}

⁠, where the unit cell is given by [0/90] and “4” is the number of repetitions. The second layup is symmetric about the middle plane, as indicated by the s. Therefore, the layups contain 8 and 16 layers. In these configurations, besides pure Lamb waves, SH wave propagation is possible in the 0∘ and 90∘ directions. In the first layup, only nonsymmetric modes occur, whereas in the second, symmetric and antisymmetric modes can be distinguished. The unit cell transfer matrix,

A_{uc}

⁠, is obtained by

A_{uc} = A_{1} A_{2},

(52)

and the global transfer matrix for [0/90]₄ is efficiently calculated by the power of n,

A = A_{uc}^{n},

(53)

where n is the number of unit cell repetitions. As in the SMM applied to Lamb waves in coupled and decoupled cases, the useful fact at this point is that for obtaining modal solutions for the symmetric configuration [0/90]

_{4 s}

⁠, no further operation equivalent to Eq. has to be imposed on A, i.e., the transfer matrix of the upper half of the symmetric layup given by Eq. is sufficient. Referring to Eq. , the dispersion equation for modal solutions in nonsymmetric layups and symmetric modes in symmetric layups is given by

a_{21} = 0,

(54)

and for antisymmetric modes in symmetric layups, the dispersion equation is given by

a_{22} = 0.

(55)

Notice that Eqs. and do not apply to the symmetric and antisymmetric modes in single layers. Applying stress-free boundary conditions to Eq. leads to the dispersion equation

i D sin γ = 0.

(56)

Because

sin γ = 0

for

γ = n π

with

n = 1, 2, 3, \dots

⁠, one has

ξ_{3} d = n π

⁠. Substituting this into Eq. , and

ξ = ω / c_{p}

and

ω = 2 π f

⁠, one reaches at

c_{p} = \sqrt{\frac{C_{66}}{ρ - C_{44} {(n / (2 f d))}^{2}}} .

(57)

Symmetric modes are obtained for

n = 0, 2, 4, \dots

and antisymmetric modes are obtained for

n = 1, 3, 5, \dots

⁠. Using n = 0 in Eq. shows that the fundamental SH mode is nondispersive at the shear wave speed

c_{p} = \sqrt{C_{66} / ρ}

⁠. The appreciable advantage of Eq. is that it does not have to be solved numerically. Entering a range of frequencies and looping through n delivers an entire dispersion diagram in an instant of time. Notice, however, that Eq. is only valid for real wavenumbers. In case of attenuation, ξ turns complex, and the relation

ξ = ω / c_{p}

does not hold anymore. Therefore, Eq. is written in the form

ξ = \sqrt{\frac{ρ ω^{2} - C_{44} {(n π / d)}^{2}}{C_{66}}},

(58)

and then phase velocity and attenuation, α, are obtained from the real and imaginary wavenumber parts,

ξ^{re}

and

ξ^{im}

⁠, respectively, according to

c_{p} = \frac{ω}{ξ^{re}}, α = ξ^{im} .

(59)

C. Viscoelasticity and fluid-loading

Viscoelasticity and fluid-loading are the main physical phenomena causing damping of elastic waves. In the first, energy of bulk waves is absorbed with propagated distance and dissipated as heat. Because guided waves consist of bulk waves, the effect appears also in those. Fluid-loading allows bulk waves radiating into the fluid surrounding the plate along which a guided wave propagates. Because of the loss of energy to the surrounding medium with traveled distance, such waves are also termed leaky waves. Although counterintuitive, it is a known fact that the leaky wave amplitude grows exponentially with distance from the plate. A graphic explanation of this is provided in Ref. 77. An additional wave type occurring on a fluid-loaded plate is the Scholte wave. This is generally an interface wave propagating at the wave speed in the fluid along an interface between a solid and a fluid. The displacement field decays exponentially away from the interface with almost all energy carried in the fluid. Whereas Scholte waves are nondispersive when propagating along an interface separating two semi-infinite half-spaces, this can change if the waveguide is a plate. Then, the two separate Scholte waves propagating along the two interfaces can couple and form a single plate wave mode just like Lamb waves, including dispersive behavior, appearance of symmetric and antisymmetric fundamental and possibly higher order modes, and a significant displacement field inside the plate (see Fig. 11). In fact, these Scholte modes are governed by the same equations as are Lamb waves in coupled and decoupled cases. The difference is that Lamb waves are damped through the fluid-loading while Scholte waves are not, provided that the fluid is nonviscous, as in the case in this study. Therefore, Scholte waves are only damped through the viscoelasticity of the plate. The same is true for SH waves. As the fluid does not support shear motion, SH waves do not couple to it.

Viscoelastic material parameters are represented by a complex stiffness tensor, where the real parts describe the elastic stiffness in the usual way, and the imaginary parts describe the damping. In this study, a viscoelastic composite material,⁷³ given in Appendix D, serves for the numerical computations. The complex stiffness tensor is given by

C_{i j}^{'} = C_{i j}^{' re} - i C_{i j}^{'}^{im}

⁠. Defining

C_{i j}^{'}

in this way corresponds to the hysteretic damping model, where a linear increase in attenuation with frequency or a constant loss per wavelength is assumed. For instance, the approximate attenuation, α₁ in Nepers per wavelength, and the phase velocity of a longitudinal bulk wave propagating along the principal axis,

x_{1}^{'}

⁠, are given by, respectively,

α_{1} = - π \frac{{C^{'}}_{11}^{im}}{{C^{'}}_{11}^{re}}, c_{p} = \sqrt{\frac{{C^{'}}_{11}^{re}}{ρ}} = \frac{ω}{ξ^{re}} .

(60)

An attenuated wave is described by a complex wavenumber, consisting of a real part and an imaginary part,

ξ = ξ^{re} + i ξ^{im}

⁠. Therefore, the general solution for the displacement field from Eq. can now be written as

u_{i} = U_{i} e^{i (ξ^{re} x_{1} + ξ_{3} x_{3} - ω t)} e^{- ξ^{im} x_{1}}, α = ξ^{im},

(61)

where the first exponential term describes the harmonic wave propagation, and the second term describes the exponential decay of the wave with propagated distance. The attenuation, α in Eq. in Nepers per meter, is represented by the imaginary part of the wavenumber. If one considers only viscoelastic damping, then all dispersion equations in Sec. II B can be used. Of course, as the dispersion equation amplitudes now have a significant imaginary part with respect to the real part, one cannot search for a zero, respectively, sign change in the real part, but one can search for a minimum in the absolute value. For instance, Eq. becomes

| det K | = min .

(62)

Fluid-loading causes a distortion of the dispersion curves to an extent ruled by the mass density ratio,

ρ = ρ_{f} / ρ_{s}

⁠, where

ρ_{f}

and

ρ_{s}

are the mass densities of the fluid and solid, respectively.⁹² For instance, whereas guided wave propagation in a steel plate immersed in water hardly experiences a notable modification compared to the free wave case, dramatic changes can be observed if the metal plate is replaced by a composite with a typical mass density of

ρ_{s} = 1500 {kg / m}^{3}

⁠.⁹³ What happens with increasing ρ is that the fluid begins to modify the boundary conditions. The free wave boundary conditions in Eq. are replaced by the continuity conditions of the normal displacements and stresses across the solid–fluid interfaces (⁠

x_{3} = 0, - d

⁠),

u_{3}^{0, - d} = {\bar{u}}_{3}^{0, - d}, σ_{33}^{0, - d} = {\bar{σ}}_{33}^{0, - d},

(63)

and by the vanishing of only the shear stresses,

σ_{i 3}^{0, - d} = 0, i = 1, 2.

(64)

Following Nayfeh,¹⁹ wave motion in a nonviscous fluid is described in the global coordinate system by specializing the field equations [Eqs. and ] to

\frac{\partial σ_{i j}}{\partial x_{j}} = ρ_{f} \frac{\partial^{2} u_{i}}{\partial t^{2}},

(65)

σ_{i j} = λ_{f} \frac{\partial u_{k}}{\partial x_{k}} δ_{i j}, i, j = 1, 2, 3,

(66)

where

λ_{f}

is Lamé's first parameter. Motion is confined to the sagittal plane such that u₂ vanishes. Through combination of Eqs. and , the expanded forms,

\frac{\partial^{2} u_{1}}{\partial x_{1}^{2}} + \frac{\partial^{2} u_{3}}{\partial x_{1} \partial x_{3}} = \frac{1}{v_{f}^{2}} \frac{\partial^{2} u_{1}}{\partial t^{2}},

(67)

\frac{\partial^{2} u_{1}}{\partial x_{1} \partial x_{3}} + \frac{\partial^{2} u_{3}}{\partial x_{3}^{2}} = \frac{1}{v_{f}^{2}} \frac{\partial^{2} u_{3}}{\partial t^{2}},

(68)

can be derived, where

v_{f} = \sqrt{λ_{f} / ρ_{f}}

is the wave speed in the fluid. Using the general solution [Eq. ] with Eqs. and delivers

[\begin{matrix} ξ^{2} - ω^{2} / v_{f}^{2} & ξ ξ_{3} \\ sym & ξ_{3}^{2} - ω^{2} / v_{f}^{2} \end{matrix}] [\begin{matrix} U_{1} \\ U_{3} \end{matrix}] = 0,

(69)

and vanishing the Christoffel matrix determinant yields two solutions,

ξ_{3, q}

⁠, q = 1,2,

ξ_{3, q} = \pm \sqrt{\frac{ω^{2}}{v_{f}^{2}} - ξ^{2}},

(70)

corresponding to downward and upward propagating waves. Substituting

ξ_{3, q}

into Eq. delivers the amplitude ratios, W_q, and from Eq. , one gets the stress amplitude, D, as

W_{q} = U_{3, q} / U_{1, q} = ξ_{3, q} / ξ, D = i ρ_{f} ω^{2} / ξ .

(71)

Hence, the complete solution is given by

(u_{1}, u_{3}, σ_{33}) = \sum_{q = 1}^{2} (1, W_{q}, D) U_{1, q} e^{i ξ_{3, q} x_{3}} .

(72)

Consider a composite plate immersed in two different nonviscous fluids, covering the upper and lower half-spaces with the corresponding quantities

ξ_{3, q}^{u}, W_{q}^{u}, D^{u}

⁠, and

ξ_{3, q}^{l}, W_{q}^{l}, D^{l}

⁠, and let a longitudinal plane wave incident from the top half-space on the plate, as sketched in Fig. 2. This generates one reflected wave in the upper half-space and one transmitted wave in the lower half-space. According to Cremer's correspondence principle found by Cremer⁹⁴ and investigated more deeply by Schoch,⁹⁵ the dispersion relation coincides with the onset of total transmission through the plate, i.e., the reflection coefficient must vanish. Although Cremer's principle in its original interpretation is applicable only for small ratios of ρ, one can make use of it also in heavily loaded cases by invoking the fluid-loading boundary conditions [Eqs. and ], introduced by Rokhlin et al.⁹² The displacements and stresses in the upper fluid are given by superposition of the incident and reflected waves such that

{[\begin{matrix} u_{1} \\ u_{3} \\ σ_{33} \end{matrix}]}^{u} = [\begin{matrix} 1 & 1 \\ W_{1}^{u} & - W_{1}^{u} \\ D^{u} & D^{u} \end{matrix}] [\begin{matrix} e^{- i ξ_{3, 1}^{u} x_{3}} \\ R e^{i ξ_{3, 1}^{u} x_{3}} \end{matrix}],

(73)

where R is the reflection coefficient. In the lower half-space, the displacements and stresses are determined by the transmitted wave,

{[\begin{matrix} u_{1} \\ u_{3} \\ σ_{33} \end{matrix}]}^{l} = [\begin{matrix} 1 \\ W_{1}^{l} \\ D^{l} \end{matrix}] T e^{- i ξ_{3, 1}^{l} (x_{3} + d)},

(74)

where T is the transmission coefficient. Specializing Eqs. and to

x_{3} = 0, - d

gives

\begin{matrix} u_{3}^{u, x_{3} = 0} = W_{1}^{u} (1 - R), \\ σ_{33}^{u, x_{3} = 0} = D^{u} (1 + R), \\ u_{3}^{l, x_{3} = - d} = W_{1}^{l} T, \\ σ_{33}^{l, x_{3} = - d} = D^{l} T . \end{matrix}

(75)

The global stiffness matrix of the composite has already been obtained and should be repeated here for convenience as

[\begin{matrix} σ^{0} \\ σ^{- d} \end{matrix}] = K [\begin{matrix} u^{0} \\ u^{- d} \end{matrix}] .

(76)

Although one can use K to obtain R, it is more convenient to derive it with the compliance matrix, S,

[\begin{matrix} u^{0} \\ u^{- d} \end{matrix}] = S [\begin{matrix} σ^{0} \\ σ^{- d} \end{matrix}], S = K^{- 1} .

(77)

In the coupled case, under consideration of the vanishing of the shear stresses according to Eq. , the expanded form of Eq. is given by

[\begin{matrix} u_{1}^{0} \\ u_{2}^{0} \\ u_{3}^{0} \\ u_{1}^{- d} \\ u_{2}^{- d} \\ u_{3}^{- d} \end{matrix}] = [\begin{matrix} s_{11} & s_{12} & s_{13} & s_{14} & s_{15} & s_{16} \\ s_{21} & s_{22} & s_{23} & s_{24} & s_{25} & s_{26} \\ s_{31} & s_{32} & s_{33} & s_{34} & s_{35} & s_{36} \\ s_{41} & s_{42} & s_{43} & s_{44} & s_{45} & s_{46} \\ s_{51} & s_{52} & s_{53} & s_{54} & s_{55} & s_{56} \\ s_{61} & s_{62} & s_{63} & s_{64} & s_{65} & s_{66} \end{matrix}] [\begin{matrix} σ_{33}^{0} \\ 0 \\ 0 \\ σ_{33}^{- d} \\ 0 \\ 0 \end{matrix}] .

(78)

Hence, the normal displacements at the top and bottom interfaces of the composite can be expressed as, respectively,

u_{3}^{0} = s_{31} σ_{33}^{0} + s_{34} σ_{33}^{- d},

(79)

u_{3}^{- d} = s_{61} σ_{33}^{0} + s_{64} σ_{33}^{- d} .

(80)

Substituting Eq. into Eqs. and , under consideration of Eq. , gives

W_{1}^{u} (1 - R) = s_{31} D^{u} (1 + R) + s_{34} D^{l} T,

(81)

W_{1}^{l} T = s_{61} D^{u} (1 + R) + s_{64} D^{l} T .

(82)

From Eqs. and , the reflection coefficient is obtained as

R = \frac{(W_{1}^{l} - s_{64} D^{l}) (W_{1}^{u} - s_{31} D^{u}) - s_{34} s_{61} D^{l} D^{u}}{(W_{1}^{l} - s_{64} D^{l}) (W_{1}^{u} + s_{31} D^{u}) + s_{34} s_{61} D^{l} D^{u}} .

(83)

In case the fluids covering the upper and lower half-spaces share the same properties, W₁ and D, Eq. simplifies to

R = \frac{(Λ - s_{64}) (Λ - s_{31}) - s_{34} s_{61}}{(Λ - s_{64}) (Λ + s_{31}) + s_{34} s_{61}}, Λ = W_{1} / D .

(84)

Suppose the lower half-space is a vacuum such that

σ_{33}^{- d} = 0

⁠, then Eq. reduces to

R = \frac{Λ - s_{31}}{Λ + s_{31}} .

(85)

However, similar to that pointed out in Sec. II B, Eqs. are only used for nonsymmetric modes. Otherwise, if symmetric and antisymmetric modes are allowed, distinct dispersion equations can be derived by invoking the symmetry specific boundary conditions, Eqs. and . For this to happen, the layup must be symmetric, and both half-spaces must be filled with the same fluid. Hence, it is sufficient to compute only the compliance matrix of the upper half of the symmetric layup,

S^{half}

⁠. With the shear stress-free upper surface, the expanded form is given by

[\begin{matrix} u_{1}^{0} \\ u_{2}^{0} \\ u_{3}^{0} \\ u_{1}^{- d / 2} \\ u_{2}^{- d / 2} \\ u_{3}^{- d / 2} \end{matrix}] = {[\begin{matrix} s_{11} & s_{12} & s_{13} & s_{14} & s_{15} & s_{16} \\ s_{21} & s_{22} & s_{23} & s_{24} & s_{25} & s_{26} \\ s_{31} & s_{32} & s_{33} & s_{34} & s_{35} & s_{36} \\ s_{41} & s_{42} & s_{43} & s_{44} & s_{45} & s_{46} \\ s_{51} & s_{52} & s_{53} & s_{54} & s_{55} & s_{56} \\ s_{61} & s_{62} & s_{63} & s_{64} & s_{65} & s_{66} \end{matrix}]}^{half} [\begin{matrix} σ_{33}^{0} \\ 0 \\ 0 \\ σ_{33}^{- d / 2} \\ σ_{13}^{- d / 2} \\ σ_{23}^{- d / 2} \end{matrix}] .

(86)

Fortuitously, by using the symmetric boundary conditions, Eq. , on Eq. , the convenient set of equations,

u_{3}^{0} = s_{31} σ_{33}^{0} + s_{34} σ_{33}^{- d / 2},

(87)

0 = s_{61} σ_{33}^{0} + s_{64} σ_{33}^{- d / 2},

(88)

is obtained, from which the reflection coefficient for symmetric modes,

R_{S}

⁠, results as

R_{S} = \frac{Λ - X_{S}}{Λ + X_{S}}, X_{S} = s_{31} - \frac{s_{34} s_{61}}{s_{64}} .

(89)

The antisymmetric boundary conditions, Eq. , used on Eq. , lead to the more unwieldy system of equations,

u_{3}^{0} = s_{31} σ_{33}^{0} + s_{35} σ_{13}^{- d / 2} + s_{36} σ_{23}^{- d / 2},

(90)

0 = s_{41} σ_{33}^{0} + s_{45} σ_{13}^{- d / 2} + s_{46} σ_{23}^{- d / 2},

(91)

0 = s_{51} σ_{33}^{0} + s_{55} σ_{13}^{- d / 2} + s_{56} σ_{23}^{- d / 2},

(92)

from which the reflection coefficient for antisymmetric modes,

R_{A}

⁠, is derived as

R_{A} = \frac{Λ - X_{A}}{Λ + X_{A}},

(93)

with

X_{A} = s_{31} - \frac{s_{35} s_{51}}{s_{55}} - (\frac{s_{35} s_{56}}{s_{55}} - s_{36}) \frac{s_{41} s_{55} - s_{45} s_{51}}{s_{45} s_{56} - s_{46} s_{55}} .

(94)

As discussed above, modal solutions for leaky waves coincide with the vanishing of the reflection coefficient. Hence, the denominators in Eqs. as well as those in Eqs. and must vanish. Because the reflection coefficient is complex, one searches for a minimum in the absolute value of its denominator. For example, the dispersion equation for leaky waves in the case of two different fluids covering the upper and lower half-spaces is given by

| (W_{1}^{l} - s_{64} D^{l}) (W_{1}^{u} + s_{31} D^{u}) + s_{34} s_{61} D^{l} D^{u} | = min .

(95)

An important note should be made concerning Scholte modes. If the fluids and waveguide are perfectly elastic, Scholte modes are not attenuated. Hence, one can use Eq. [respectively, the denominators of Eqs. , , ] with searching only on the real wavenumber axis. In this paper, the fluids are considered nonviscous. However, if the waveguide is viscoelastic, Scholte waves are damped, and the search must be performed in the complex wavenumber plane. Equation (and the others) can still be used but with an important modification. Instead of using the positive solutions of the out-of-plane wavenumber components in the upper and lower fluids,

ξ_{3, 1}^{u}

and

ξ_{3, 1}^{l}

⁠, as obtained from Eq. , the negative solutions of the out-of-plane wavenumber components in the upper and lower fluids,

ξ_{3, 2}^{u}

and

ξ_{3, 2}^{l}

⁠, must be used. Accordingly, the positive amplitude ratios,

W_{1}^{u}

and

W_{1}^{l}

⁠, in Eq. are replaced by the negative amplitude ratios

W_{2}^{u}

and

W_{2}^{l}

⁠.

FIG. 2.

View large Download slide

Fluid-loaded composite plate.

In the decoupled case, the displacement and stress components in leaky Lamb and Scholte waves,

u_{2}^{0, - d}

and

σ_{23}^{0, - d}

⁠, in Eq. vanish, and the compliance matrix reduces to 4 × 4 elements. For this case, Eqs. and can be used by changing the compliance matrix elements according to

s_{31} \to s_{21}, s_{34} \to s_{23}, s_{61} \to s_{41}, s_{64} \to s_{43}

⁠. Likewise,

u_{2}^{0, - d / 2}

and

σ_{23}^{0, - d / 2}

vanish in Eq. , leading to

X_{S} = s_{21} - \frac{s_{23} s_{41}}{s_{43}}, X_{A} = s_{21} - \frac{s_{24} s_{31}}{s_{34}} .

(96)

D. Dispersion curve tracing

Dispersion curve tracing in the attenuated case requires a two-dimensional search for modal solutions in the complex wavenumber plane rather than a one-dimensional search on the real axis only. The basic procedure used by DC for finding a modal solution is rather conventional but robust. A particular dispersion equation is evaluated on a raster of complex wavenumbers covering a certain range of real and imaginary parts. For example, Fig. 3 plots the dispersion equation amplitude of Eq. (95) calculated at 500 kHz for wave propagation along $Φ = 0 \circ$ in a 2 mm thick viscoelastic composite layup [0/90/45/–45]_2s immersed on both sides in water. Equation (95) is evaluated on a 51 × 51 point raster in the complex wavenumber plane and plotted in the phase velocity-attenuation space. The minimum indicates the modal solution of the S₁ Lamb wave, which is located at $c_{p} = 5.98516$ m/ms and $α = 15.3463$ Np/m. When a minimum has been located by the initial coarse sweep, the second step is to converge on it by using two-dimensional bisection until the desired accuracy has been achieved. Calculating such rasters and following the minimum as it shifts with frequency is how DC traces dispersion curves. A more detailed description of this is given in DC user's manual.⁸²

FIG. 3.

View large Download slide

Dispersion equation amplitude of Eq. (95) calculated at 500 kHz for wave propagation along $Φ = 0 °$ in a 2 mm thick viscoelastic composite layup [0/90/45/–45]_2s immersed on both sides in water. The minimum indicates the modal solution of the S₁ Lamb wave.

E. Mode shape calculation

Once a modal solution at a certain frequency has been obtained, one can calculate the through-thickness mode shape for that frequency-wavenumber pair, i.e., the components of displacement, stress, and strain, as well as derived quantities such as energy density and power flow density versus x₃. The last two quantities are needed to calculate the energy velocity components.

The starting point is to determine the displacements at the top and bottom surfaces of the composite. In case the surfaces are stress-free, one should follow the formalism described by Rokhlin and Wang.³⁹ Accordingly, the displacements at the top and bottom of the system are given by

u^{0} = P_{0}^{-} U_{in} + P_{0}^{+} U_{r}, u^{- d} = P_{- d}^{-} U_{t},

(97)

where

U_{in}, U_{r}

⁠, and

U_{t}

are the incident, reflected, and transmitted wave amplitudes at the top and bottom surfaces, respectively. Now, as an alternative to the procedure applied by Rokhlin and Wang, one can benefit from the stress-free boundary conditions. Therefore, by substituting Eq. into

[\begin{matrix} σ^{0} \\ σ^{- d} \end{matrix}] = [\begin{matrix} K_{11} & K_{12} \\ K_{21} & K_{22} \end{matrix}] [\begin{matrix} u^{0} \\ u^{- d} \end{matrix}] = 0,

(98)

the simpler relation,

[\begin{matrix} U_{r} \\ U_{t} \end{matrix}] = {[\begin{matrix} - K_{11} P_{0}^{+} & - K_{12} P_{- d}^{-} \\ K_{21} P_{0}^{+} & K_{22} P_{- d}^{-} \end{matrix}]}^{- 1} [\begin{matrix} K_{11} P_{0}^{-} \\ - K_{21} P_{0}^{-} \end{matrix}] U_{in},

(99)

is obtained with

U_{in} = [\begin{matrix} 1 \\ 0 \\ 0 \end{matrix}] .

(100)

With the known reflected and transmitted wave amplitudes,

U_{r}

and

U_{t}

⁠, Eq. can be solved. However, if the half-spaces are filled with fluids, a different routine must be used. From the composite, one bulk wave is emitted into the upper half-space and one bulk wave is emitted into the lower half-space. Therefore, at the interfaces

x_{3} = 0, - d

⁠, the normal stresses and displacements are given by

\begin{array}{l} u_{3}^{u, x_{3} = 0} = - W_{1}^{u} U^{u}, \\ σ_{33}^{u, x_{3} = 0} = D^{u} U^{u}, \\ u_{3}^{l, x_{3} = - d} = W_{1}^{l} U^{l}, \\ σ_{33}^{l, x_{3} = - d} = D^{l} U^{l}, \end{array}

(101)

where

U^{u}

is the wave amplitude in the upper half-space, and

U^{l}

is the wave amplitude in the lower half-space. By setting

U^{u}

to unity and substituting Eq. into Eq. , one finds

U^{l} = - \frac{W_{1}^{u} + s_{31} D^{u}}{s_{34} D^{l}} .

(102)

Now, the displacements at the upper and lower surfaces of the composite are given by, respectively,

\begin{matrix} u_{i}^{0} = s_{i 1} D^{u} + s_{i 4} D^{l} U^{l}, i = 1, 2, 3, \\ u_{i}^{- d} = s_{k 1} D^{u} + s_{k 4} D^{l} U^{l}, k = 4, 5, 6. \end{matrix}

(103)

After having solved either Eq. or Eq. , the next step is to get the displacements at the inner layer interfaces. This is achieved by using the backpropagation recursive algorithm proposed by Rokhlin and Wang. The algorithm starts with the lowermost layer. The displacements at the bottom of this layer,

u^{- d}

⁠, have already been obtained. Then,

u_{m}

at the top of this layer is given by

u_{m} = X K_{21}^{*} u^{0} - X K_{12}^{m} u_{m + 1},

(104)

where

X = {(K_{11}^{m} - K_{22}^{*})}^{- 1}, K^{m}

is the local stiffness matrix of the mth layer, and

K^{*}

is the total stiffness matrix of the top m – 1 layers. This procedure has to be repeated until one has the displacements at the bottom of the top layer. Now, the wave amplitudes in each layer are calculated by

[\begin{matrix} U_{m} \\ U_{m + 1} \end{matrix}] = {[\begin{matrix} P^{-} & P^{+} H \\ P^{-} H & P^{+} \end{matrix}]}_{m}^{- 1} [\begin{matrix} u_{m} \\ u_{m + 1} \end{matrix}] .

(105)

Finally, the field components of displacement, stress, and strain are obtained layer wise by superposing the contributing bulk waves. For example, according to Eq. , the normal displacement in the mth layer is calculated by

u_{3} = W_{1} U_{1} e^{i ξ_{3, 1} x_{3}} + W_{3} U_{2} e^{i ξ_{3, 3} x_{3}} + W_{5} U_{3} e^{i ξ_{3, 5} x_{3}} - W_{1} U_{4} e^{i ξ_{3, 1} (d - x_{3})} - W_{3} U_{5} e^{i ξ_{3, 3} (d - x_{3})} - W_{5} U_{6} e^{i ξ_{3, 5} (d - x_{3})},

(106)

where d is the thickness of the mth layer in this case, and x₃ is a value ranging from zero to d. The first three terms in Eq. correspond to the downward propagating bulk waves, and the second three terms correspond to the upward propagating bulk waves. If the guided waves interact with the surrounding fluids, it is interesting to calculate the field components inside the half-spaces also. Only one bulk wave propagates in each, where the wave amplitude in the upper half-space,

U^{u}

⁠, is set to unity and that in the lower half-space,

U^{l}

⁠, is calculated using Eq. . As an example, the normal displacements in the upper and lower fluid-filled half-spaces are obtained using, respectively,

u_{3}^{u} = - W_{1}^{u} e^{i ξ_{3, 1}^{u} x_{3}}, u_{3}^{l} = W_{1}^{l} U^{l} e^{i ξ_{3, 1}^{l} x_{3}},

(107)

where x₃ ranges from zero at the respective composite surface to an arbitrary distance in the corresponding half-space.

Like that mentioned in Sec. II C, if Scholte modes are considered in a viscoelastic scenario, the negative solutions of the out-of-plane wavenumber components, $ξ_{3, 2}^{u}$ and $ξ_{3, 2}^{l}$ ⁠, must be used in the calculations. All equations given in Sec. II E can be also used in decoupled cases in the reduced form. In the case of SH modes, Eqs. (101)–(103) and (107) can be ignored.

III. CLASSIFICATION OF SOLUTIONS

The most adverse effect encountered during dispersion curve tracing is the so-called jumping mode problem. This means that the dispersion curve tracing algorithm jumps from the currently traced mode to another. More generally spoken, it is the problem of assigning a solution to the corresponding mode. This occurs often in situations where two modes come close or cross each other. Clearly, assigning a solution to the wrong mode becomes more likely if different mode families are traced using the same dispersion equation. Therefore, to reduce or even eliminate the problem, each mode family must be obtained by using specific dispersion equations. Then, only jumps to modes of the same family are possible, but this is already less likely because dispersion curves of the same family do not cross each other or do so less frequently, depending on the configuration.

In general, a classification of solutions is possible according to the symmetry and coupling properties of the considered waveguide and modes. The four possible configurations with the corresponding mode families appearing in a fluid-loaded waveguide are listed in Table I. The prerequisite is that one must know which kinds of mode families occur in a given configuration. The decoupling into pure Lamb, Scholte, and SH waves depends on material symmetries, the layup, and propagation direction of the guided waves. In orthotropic and transversely isotropic media such as unidirectional composite plies, decoupling occurs for wave propagation parallel or normal to the fibers (⁠ $Φ = 0 \circ, 90 \circ$ ⁠). This is possible in a single layer or cross-ply laminate [0/90] $_{n (s)}, n = 1, 2, \dots$ ⁠. In cubic materials, decoupling occurs for guided wave propagation along $Φ = 0 \circ, 45 \circ, 90 \circ$ ⁠. Finally, in isotropic materials, decoupling occurs for any propagation direction. In case a layup contains different materials, the system allows decoupling only if all contributing layers support decoupling in the given situation. In symmetric systems, one can distinguish symmetric and antisymmetric Lamb, SH, and, in case of fluid-loading, Scholte waves. A system consisting of a waveguide located between the upper and lower half-spaces is symmetric if the top half mirrors the bottom half in terms of layer (fiber) orientations, layer thicknesses, layer materials, and fluid. If the system is nonsymmetric, so are the modes, and no separation into symmetric and antisymmetric modes is possible. The only exception from this are SH modes. If the waveguide is symmetric, but the upper and lower half-spaces are filled with different nonviscous fluids, the SH modes are still symmetric and antisymmetric because they do not couple to the fluids.

TABLE I.

Possible configurations of mode families propagating in fluid-loaded anisotropic laminates.

Symmetric	Decoupled	Mode families
Yes	Yes	S, A, S^SH, A^SH, S^Scholte, A^Scholte
Yes	No	S, A, S^Scholte, A^Scholte
No	Yes	B,a B^SH (S^SH, A^SH),b B^Scholte
No	No	B, B^Scholte

Symmetric	Decoupled	Mode families
Yes	Yes	S, A, S^SH, A^SH, S^Scholte, A^Scholte
Yes	No	S, A, S^Scholte, A^Scholte
No	Yes	B,a B^SH (S^SH, A^SH),b B^Scholte
No	No	B, B^Scholte

^a

Because there exists no generally accepted notation for nonsymmetric modes, these are denoted here with the letter B.

^b

If the half-spaces are covered by different fluids but the layup is symmetric, S^SH and A^SH can be distinguished.

Once it is clear which mode families must be expected, the appropriate dispersion equations can be applied to trace them separately. Table II lists those for Lamb and Scholte waves in coupled and decoupled cases, depending on the presence or absence of viscoelasticity and fluid-loading, and Table III lists the dispersion equations of SH waves. Only in nonattenuated cases, one can search for a zero in the dispersion equation amplitude, such as $det K = 0$ ⁠. Actually, because $det K$ never becomes zero, one searches for a sign change. Notice, however, that a sign change does not necessarily indicate a modal solution since $det K$ also can switch between positive and negative infinity (but not a₂₁ and a₂₂ in the case of SH waves). Therefore, it must be checked that the absolute value is a (local) minimum, i.e., $| det K | = min$ ⁠. In attenuated cases, one searches only for minima in the complex dispersion equation amplitude.

TABLE II.

Mode family specific dispersion equations of (leaky) Lamb waves and Scholte waves in coupled and decoupled cases.

Family	Viscoelasticity and fluid-loading	Dispersion equations
S	—	det $K_{S} = 0, \|$ det $K_{S} \| = min$
A	—	det $K_{A} = 0, \|$ det $K_{A} \| = min$
B	—	det $K = 0, \|$ det $K \| = min$
S	V	$\|$ det $K_{S} \| = min$
A	V	$\|$ det $K_{A} \| = min$
B	V	$\|$ det $K \| = min$
S	F or F + V	$\| R_{S} \| = min$
A	F or F + V	$\| R_{A} \| = min$
B	F or F + V	$\| R \| = min$
S^Scholte	F	$\| R_{S} \| = min$
A^Scholte	F	$\| R_{A} \| = min$
B^Scholte	F	$\| R \| = min$
S^Scholte	F + V	$\| R_{S} \| = min, - ξ_{3}^{f}$ a
A^Scholte	F + V	$\| R_{A} \| = min, - ξ_{3}^{f}$ a
B^Scholte	F + V	$\| R \| = min, - ξ_{3}^{f}$ a

Family	Viscoelasticity and fluid-loading	Dispersion equations
S	—	det $K_{S} = 0, \|$ det $K_{S} \| = min$
A	—	det $K_{A} = 0, \|$ det $K_{A} \| = min$
B	—	det $K = 0, \|$ det $K \| = min$
S	V	$\|$ det $K_{S} \| = min$
A	V	$\|$ det $K_{A} \| = min$
B	V	$\|$ det $K \| = min$
S	F or F + V	$\| R_{S} \| = min$
A	F or F + V	$\| R_{A} \| = min$
B	F or F + V	$\| R \| = min$
S^Scholte	F	$\| R_{S} \| = min$
A^Scholte	F	$\| R_{A} \| = min$
B^Scholte	F	$\| R \| = min$
S^Scholte	F + V	$\| R_{S} \| = min, - ξ_{3}^{f}$ a
A^Scholte	F + V	$\| R_{A} \| = min, - ξ_{3}^{f}$ a
B^Scholte	F + V	$\| R \| = min, - ξ_{3}^{f}$ a

^a

Use negative out-of-plane wavenumber component solutions in the fluids.

TABLE III.

Mode family specific dispersion equations of SH waves.

Family	Viscoelastic	Layers	Dispersion equations
S^SH, B^SH	No	Multi	$a_{21} = 0$
A^SH	No	Multi	$a_{22} = 0$
S^SH	No	Single	$c_{p} = \sqrt{\frac{C_{66}}{ρ - C_{44} {(n / (2 f d))}^{2}}}$ ⁠, n even
A^SH	No	Single	$c_{p} = \sqrt{\frac{C_{66}}{ρ - C_{44} {(n / (2 f d))}^{2}}}$ ⁠, n odd
S^SH, B^SH	Yes	Multi	$\| a_{21} \| = min$
A^SH	Yes	Multi	$\| a_{22} \| = min$
S^SH	Yes	Single	$ξ = \sqrt{\frac{ρ ω^{2} - C_{44} {(n π / d)}^{2}}{C_{66}}}$ ⁠, n even
A^SH	Yes	Single	$ξ = \sqrt{\frac{ρ ω^{2} - C_{44} {(n π / d)}^{2}}{C_{66}}}$ ⁠, n odd

Family	Viscoelastic	Layers	Dispersion equations
S^SH, B^SH	No	Multi	$a_{21} = 0$
A^SH	No	Multi	$a_{22} = 0$
S^SH	No	Single	$c_{p} = \sqrt{\frac{C_{66}}{ρ - C_{44} {(n / (2 f d))}^{2}}}$ ⁠, n even
A^SH	No	Single	$c_{p} = \sqrt{\frac{C_{66}}{ρ - C_{44} {(n / (2 f d))}^{2}}}$ ⁠, n odd
S^SH, B^SH	Yes	Multi	$\| a_{21} \| = min$
A^SH	Yes	Multi	$\| a_{22} \| = min$
S^SH	Yes	Single	$ξ = \sqrt{\frac{ρ ω^{2} - C_{44} {(n π / d)}^{2}}{C_{66}}}$ ⁠, n even
A^SH	Yes	Single	$ξ = \sqrt{\frac{ρ ω^{2} - C_{44} {(n π / d)}^{2}}{C_{66}}}$ ⁠, n odd

IV. NUMERICAL EXAMPLES

In this section, dispersion diagrams of phase velocity and attenuation are presented. For the calculations, the viscoelastic composite from Ref. 73, whose material constants are given in Appendix D, was used. In DC, the material is available in the list of orthotropic materials as CarbonEpoxy_Hernando_2015_Viscoelastic. Fluid-loading is present in the form of water-filled (⁠ $ρ_{f} = 1000 {kg / m}^{3}, v_{f} = 1500 m / s$ ⁠) half-spaces on both sides of the composite. The same fluid is chosen on both sides to demonstrate the distinction of symmetric and antisymmetric modes. Three different layups are considered, namely, [0/90/45/–45]_2s, [0/90]_4s, and [0/90/45/–45]_50s. The direction of wave propagation is always along $Φ = 0 \circ$ ⁠. The first two layups contain 16 layers, summing up to a thickness of 2 mm. These were calculated with DC version 2.4 and DISPERSE version 2.0.20f and serve as validation examples. The third layup contains 400 layers with a total thickness of 50 mm, proving that DC can handle thick-walled composite structures.

A. Sixteen layers composite

Figure 4 displays the phase velocity dispersion diagram for the [0/90/45/–45]_2s layup immersed in water. The bold curves were calculated with DISPERSE, and the thin curves were obtained with DC. As the system is symmetric, one can distinguish symmetric (gray lines) and antisymmetric (black lines) modes. Solid lines represent leaky Lamb modes while the dashed-dotted lines indicate Scholte modes. All modes have coupled polarization. Although DC missed one higher order mode, whereas DISPERSE missed A₀, it is obvious that the other curves are in good agreement. The corresponding attenuation dispersion diagram is depicted in Fig. 5, where DC and DISPERSE are in good agreement again. Most remarkable in Fig. 5 is how strongly the attenuation of A₀ rises at low frequency-thicknesses compared to the other fundamental modes. This is because the water weighs heavily against the strong and predominantly out-of-plane vibration possessed by A₀. The out-of-plane displacement is generally more pronounced in composites than in isotropic media as high stiffness carbon fibers are usually oriented in the x₁-x₂-plane while the x₃-direction (out-of-plane) is governed by the low stiffness of the polymer matrix. Another feature is the strong rise of S₂ above 3.2 MHz·mm, reaching a maximum of 2363 Np·mm/m at 3.87 MHz·mm. By contrast, the Scholte modes are only weakly attenuated. $A_{0}^{Scholte}$ and $S_{0}^{Scholte}$ reach only 14.81 Np·mm/m at 5 MHz·mm. This is due to viscoelastic damping of the composite, which is a much weaker effect than the heavy water-loading dominating the damping of the leaky modes.

FIG. 4.

View large Download slide

Phase velocity dispersion diagram for wave propagation along $Φ = 0 °$ in a 2 mm thick viscoelastic composite layup [0/90/45/–45]_2s immersed in water. Bold curves were calculated with DISPERSE and thin curves were calculated with DC.

FIG. 5.

View large Download slide

Attenuation dispersion diagram for wave propagation along $Φ = 0 °$ in a 2 mm thick viscoelastic composite layup [0/90/45/–45]_2s immersed in water.

As the next example, the [0/90]_4s layup is considered. Because wave propagation is in the $Φ = 0 \circ$ -direction, decoupling into sagittal plane motion and shear horizontal motion takes place. Pure leaky Lamb waves and pure SH waves, as well as Scholte waves with only sagittal motion, can propagate in this configuration. The phase velocity dispersion diagram is shown in Fig. 6. The SH modes are represented by dashed lines. Whereas DC did not trace the SH modes up to the top phase velocity of 20 m/ms, DISPERSE missed A₀ again and left gaps in two dispersion curves. However, the good agreement of both softwares is obvious. The attenuation dispersion is displayed in Fig. 7. A₀ shows the similar behavior as discussed in the above coupled case. S₀ reaches a maximum attenuation of nearly 2000 Np·mm/m at 4 MHz·mm. This is where the phase velocity increases steeply in Fig. 6. Although SH waves are damped only by the viscoelasticity of the composite, their damping rises quickly above that of the leaky modes as one approaches the cutoff frequencies.

FIG. 6.

View large Download slide

Phase velocity dispersion diagram for wave propagation along $Φ = 0 °$ in a 2 mm thick viscoelastic composite layup [0/90]_4s immersed in water.

FIG. 7.

View large Download slide

Attenuation dispersion diagram for wave propagation along $Φ = 0 °$ in a 2 mm thick viscoelastic composite layup [0/90]_4s immersed in water.

At this point, a few comments should be made about the mode naming applied in DC. Consider, again, the coupled case in Fig. 4. The reader may be concerned that the mode which is denoted S₀ should be called $S_{0}^{SH}$ (or SH₀) and S₁ should be S₀. Notice, however, that all modes are coupled and SH modes, in a strict sense, do not exist. All modes have nonzero polarization components in all three directions, whereas a SH mode is polarized exclusively in the shear horizontal direction. Also notice that $S_{0}^{SH}$ is always nondispersive, as observed in Fig. 6, which is also not the case in Fig. 4. There exists the term of quasi-SH modes, implying that the sagittal motion is small compared to the shear horizontal component. However, considering that the mode shape depends on the frequency, it would be a difficult task to determine for every mode whether it is a quasi-SH mode or not. To avoid any ambiguities and for simplicity, DC denotes all coupled modes of the same family with ascending numbers. For SH modes, DC does not use the traditional notations of SH₀, SH₁, SH₂, etc., but symmetric and antisymmetric modes are counted separately, namely, $S_{0}^{SH}$ ⁠, $S_{1}^{SH}$ ⁠, and $A_{1}^{SH}$ ⁠, $A_{2}^{SH}$ ⁠, and so on to be consistent with the similar notation of Lamb and Scholte waves. Nonsymmetric modes are denoted with the letter “B.”

B. Four hundred layers composite

In this section, the capabilities of DC shall be demonstrated. The idea is to calculate a dispersion diagram for the most demanding but still realistic case. In the previous study in Ref. 49, a perfectly elastic [0/90/45/–45]_50s composite layup with stress-free surfaces had been calculated. This layup contained 400 layers, similar to the thickest laminates in modern rocket booster pressure vessels. However, the calculation did not incorporate attenuation, thereby avoiding a significant complication. Therefore, the same layup is calculated again but with the viscoelastic composite and water-loading on both sides. The computation was performed with DC version 2.4 installed as a stand-alone software on a personal computer with an Intel^® Core™ i7-10850H CPU (central processing unit) at $6 \times 2.7$ GHz (Intel Corporation, Santa Clara, CA) and 16 GB RAM (random access memory). The dispersion diagram was calculated up to a frequency of 100 kHz (5 MHz·mm) with steps of 0.25 kHz, yielding 400 samples in the fundamental modes. To speed up the calculation, dual core computing was enabled in DC such that symmetric and antisymmetric modes were traced simultaneously. The calculation took 38 min. A bottleneck was finding and starting the dispersion curves at the cutoff frequencies. Figure 8 displays the phase velocity dispersion diagram. All modes have coupled polarization. Compared to the 2 mm thick coupled case depicted in Fig. 4, a smaller number of higher order modes was found, and S₀ is nearly nondispersive now. The missing higher order modes are still there, but the corresponding minima in the dispersion equation amplitudes have become too small to be found by DC. The attenuation is plotted in Fig. 9. It has similarities with that of the thin laminate presented in Fig. 5, but it is generally lower.

FIG. 8.

View large Download slide

Phase velocity dispersion diagram for wave propagation along $Φ = 0 °$ in a 50 mm thick viscoelastic composite layup [0/90/45/–45]_50s immersed in water.

FIG. 9.

View large Download slide

Attenuation dispersion diagram for wave propagation along $Φ = 0 °$ in a 50 mm thick viscoelastic composite layup [0/90/45/–45]_50s immersed in water.

Through-thickness mode shapes are presented to shed light on the wave structure inside the waveguide and on the wave field radiated into the water-filled half-spaces. The mode shapes were calculated based on the dispersion curves displayed in Figs. 8 and 9. It should be noted that DC uses the same plotting conventions as DISPERSE, namely, power-normalization of the amplitudes (neglecting the power flow in the half-spaces) and phase normalization with respect to the second sample from the top. In Fig. 10, the displacement component amplitudes of the S₂ mode at the maximum attenuation located at 4.075 MHz·mm are shown. The white area represents the composite, and the light gray areas indicate the water-filled half-spaces to a depth of 50 mm on both sides. The symmetric character of S₂ is evident through the in-plane displacement component, u₁, because it is symmetric about the middle axis of the composite. An interesting finding is that the shear horizontal displacement, u₂, seems to vanish inside the composite although all modes should possess coupled polarization. In fact, u₂ does not vanish inside the composite (unlike in the water), but it reaches a maximum amplitude of only 0.03 nm. This is too small to be observed in Fig. 10. This phenomenon is related to the homogenization of periodically layered laminates. A detailed discussion on this topic is given in the book by Rokhlin et al.⁴¹ Accordingly, if the wavelength of the bulk waves inside the layers is large compared to the layer thickness, the effect of the layering becomes negligible. The laminate behaves like an effective medium with averaged elastic properties, and the coupling of the bulk waves becomes small. The strong out-of-plane component, u₃, at the surfaces of the composite is related to the transmittance of energy into the half-spaces, which, in turn, leads to the strong damping of S₂ at 4.075 MHz·mm. As a second example, the $A_{0}^{Scholte}$ mode at 2.5 MHz·mm is considered. It propagates at a phase velocity of 1.31 m/ms and has a low damping of 9.581 Np·mm/m, caused by the composite's viscoelasticity. As depicted in Fig. 11, the most striking difference to leaky waves is that Scholte waves are decaying exponentially inside the half-spaces. From u₁, the antisymmetric character of the mode becomes obvious. Because of the aforementioned homogenization effect, u₂ becomes small with a maximum amplitude of 0.14 nm. It should be noted that with the decreasing of the frequency-thickness product further, $A_{0}^{Scholte}$ possesses an even higher proportion of displacement inside the composite. This, together with the low damping, makes Scholte modes suitable for the NDI of plates.

FIG. 10.

Through-thickness displacements of the S2 mode at 4.075 MHz·mm for wave propagation along Φ = 0 ° in a 50 mm thick viscoelastic composite layup [0/90/45/–45]50s immersed in water.

View large Download slide

Through-thickness displacements of the S₂ mode at 4.075 MHz·mm for wave propagation along $Φ = 0 °$ in a 50 mm thick viscoelastic composite layup [0/90/45/–45]_50s immersed in water.

FIG. 11.

View large Download slide

Through-thickness displacements of the $A_{0}^{Scholte}$ mode at 2.5 MHz·mm for wave propagation along $Φ = 0 °$ in a 50 mm thick viscoelastic composite layup [0/90/45/–45]_50s immersed in water.

V. DISCUSSION

Although there is no doubt about the advantages of using mode family specific dispersion equations when it comes to the speed and robustness of dispersion curve tracing, the reader may be concerned about the described way of searching and tracing modal solutions by calculating dispersion equation amplitudes on a raster of complex wavenumbers and using bisection to converge on roots. Other root-finding algorithms, such as Newton-Raphson and Muller's method, are more commonly employed because of their better rate of convergence. However, while bisection is the slowest method, it is also the most robust. When trying to accelerate DC using Muller's method, the author found that the initial guess must be sufficiently close to the root to converge on it reliably. Consequently, dispersion curve tracing robustness became much worse, but robustness was considered more important than speed because DC is an open source software, which must be able to handle a wide range of cases successfully. In that sense, it should be noted that DISPERSE uses a similar routine in which a coarse sweep is performed to find minima, followed by two-dimensional iterative bisection.^22,96 Other root-searching algorithms are described in Ref. 85 and references therein as well as in Ref. 97.

Nonetheless, dispersion curve tracing using DC is not perfectly robust as evident through the missing higher order modes in the 400 layers composite in Fig. 8 compared to the 16 layers composite in Fig. 4. The reason for this is that the dispersion equation amplitude minima at the cutoff frequencies of the missing modes, where DC tries to start tracing, had become so weak that the raster of complex wavenumbers was too coarse to detect them. The reason why the minima became so weak, in turn, lies in the homogenization effect, leading to gradual decoupling into quasi-pure Lamb and quasi-SH waves when the layer thickness is small compared to the wavelength. This decoupling is weak in the 16 layers composite but strong in the 400 layers composite. The missing modes in the latter are exactly the quasi-SH waves because their roots cause much smaller minima in the coupled dispersion equation amplitude than those of the quasi-pure Lamb waves. A possible solution to the problem would be to increase the complex wavenumber raster density in DC options and/or to start the tracing at lower phase velocities where the minima are more pronounced.

Finally, the presented numerical model of DC covers only multilayered anisotropic flat plates, which can be viscoelastic and surrounded by nonviscous fluid-filled half-spaces. Viscous fluid-loading as well as liquid inner layers and finite thickness fluid layers surrounding the composite are not included but can be obtained with minor manipulation in the equations. Similarly, guided wave propagation in rods and pipes under consideration of fluid-interaction and viscoelasticity are the next topics to be included using root-finding methods. However, more arbitrary geometries would require a discretizing method.

VI. CONCLUSIONS

It was shown that guided wave dispersion diagrams for fluid-loaded viscoelastic composites consisting of up to 400 layers can be obtained with DC, although some higher order modes were missed in some cases. DC is an open source software available at GitHub.⁸² Validation against DISPERSE was made with two composite layups containing 16 layers.

To improve efficiency and robustness of the dispersion curve tracing, DC determines which kinds of mode families occur in each configuration, and then it applies mode family specific dispersion equations. By classifying modal solutions according to their mode type as well as symmetry and coupling properties, symmetric, antisymmetric, and nonsymmetric leaky Lamb, Scholte, and SH waves can be traced separately. This helps minimize the risk of jumping to a different mode during dispersion curve tracing.

For leaky Lamb and Scholte modes in coupled and decoupled cases, DC uses SMM. Mode family specific dispersion equations were presented for each combination of symmetry as well as presence and absence of fluid-loading and viscoelasticity. For the latter, DC employs the hysteretic damping model. Although viscoelasticity does not affect the stress-free boundary conditions, fluid-loading does affect the stress-free boundary conditions. Therefore, instead of searching for minima in the global stiffness matrix determinant, one searches for minima in the reflection coefficient. The dispersion equations were derived for different and similar fluids covering the upper and lower half-spaces, as well as for the case of only one half-space filled with a fluid. In particular, by applying symmetry specific boundary conditions at the half thickness of the composite, distinct dispersion equations for symmetric and antisymmetric modes were derived for stress-free and fluid-loaded cases. These equations are computationally less expensive than the general equation, and they enable more stable dispersion curve tracing. Another beneficial measure is the decoupling of the 6 × 6 stiffness matrix into a 4 × 4 stiffness matrix for pure leaky Lamb and Scholte waves and a 2 × 2 transfer matrix for SH waves when wave propagation takes place along principal directions.

For SH modes appearing in decoupled cases, DC uses TMM. The global transfer matrix was obtained by a simplified algorithm. Few and concise dispersion equations were presented, covering all combinations of symmetry and viscoelasticity or perfect elasticity in single or multilayered laminates. Again, distinct equations for symmetric and antisymmetric modes were given. Nonviscous fluid-loading does not affect SH waves.

ACKNOWLEDGMENTS

Parts of this work were funded by the German Federal Ministry for Economic Affairs and Climate Action in the frame of the LuFo VI-2 project NEUTRON under the funding indicator 20M2114C. The responsibility for the content of this paper lies with the author. DC supports the ultrasonic inspection of components manufactured in NEUTRON. The author would like to thank Professor Dr. Michael Lowe, Professor Dr. Michel Castaings, Professor Dr. Stanislav Rokhlin, Professor Dr. Marc Deschamps, Dr. Eric Ducasse, Professor Dr. Victor Giurgiutiu, and Professor Dr. Markus Sause for their support during the development of DC.

APPENDIX A: TRANSFORMED STIFFNESS MATRIX

\begin{matrix} C_{11} = C_{11}^{'} {cos}^{4} Φ + C_{22}^{'} {sin}^{4} Φ + 2 (C_{12}^{'} + 2 C_{66}^{'}) {sin}^{2} Φ {cos}^{2} Φ, \\ C_{12} = (C_{11}^{'} + C_{22}^{'} - 2 C_{12}^{'} - 4 C_{66}^{'}) {sin}^{2} {cos}^{2} + C_{12}^{'}, \\ C_{13} = C_{13}^{'} {cos}^{2} Φ + C_{23}^{'} {sin}^{2} Φ, \\ C_{16} = (C_{12}^{'} + 2 C_{66}^{'} - C_{11}^{'}) sin Φ {cos}^{3} Φ + (C_{22}^{'} - C_{12}^{'} - 2 C_{66}^{'}) cos Φ {sin}^{3} Φ, \\ C_{22} = C_{11}^{'} {sin}^{4} Φ + C_{22}^{'} {cos}^{4} Φ + 2 (C_{12}^{'} + 2 C_{66}^{'}) {sin}^{2} Φ {cos}^{2} Φ, \\ C_{23} = C_{23}^{'} {cos}^{2} Φ + C_{13}^{'} {sin}^{2} Φ, \\ C_{26} = (C_{12}^{'} + 2 C_{66}^{'} - C_{11}^{'}) cos Φ {sin}^{3} Φ + (C_{22}^{'} - C_{12}^{'} - 2 C_{66}^{'}) sin Φ {cos}^{3} Φ, \\ C_{33} = C_{33}^{'}, \\ C_{36} = (C_{23}^{'} - C_{13}^{'}) sin Φ cos Φ, \\ C_{44} = {C^{'}}_{44} {cos}^{2} Φ + C'_{55} {sin}^{2} Φ, \\ C_{45} = ({C^{'}}_{44} - {C^{'}}_{55}) sin Φ cos Φ, \\ C_{55} = {C^{'}}_{55} {cos}^{2} Φ + {C^{'}}_{44} {sin}^{2} Φ, \\ C_{66} = C_{66}^{'} + (C_{11}^{'} + C_{22}^{'} - 2 C_{12}^{'} - 4 C_{66}^{'}) {sin}^{2} Φ {cos}^{2} Φ . \end{matrix}

(A1)

APPENDIX B: POLYNOMIAL COEFFICIENTS

\begin{matrix} A_{1} = {[C_{11} C_{33} C_{44} + C_{33} C_{55} C_{66} - C_{36}^{2} C_{55} - C_{13}^{2} C_{44} + 2 (C_{13} C_{36} C_{45} + C_{13} C_{45}^{2} - C_{13} C_{44} C_{55} - C_{16} C_{33} C_{45})] ξ^{2} \\ + (C_{45}^{2} - C_{33} C_{44} - C_{33} C_{55} - C_{44} C_{55}) ρ ω^{2}} / Δ, \\ A_{2} = {[C_{11} C_{33} C_{66} + C_{11} C_{44} C_{55} - C_{11} C_{36}^{2} - C_{11} C_{45}^{2} - C_{13}^{2} C_{66} - C_{16}^{2} C_{33} + 2 (C_{16} C_{36} C_{55} + C_{13} C_{16} C_{36} + C_{13} C_{16} C_{45} - C_{11} C_{36} C_{45} - C_{13} C_{55} C_{66})] ξ^{4} + [C_{13}^{2} + C_{45}^{2} + C_{36}^{2} - C_{11} C_{33} - C_{11} C_{44} - C_{33} C_{66} - C_{55} C_{66} - C_{44} C_{55} \\ + 2 (C_{13} C_{55} + C_{16} C_{45} + C_{36} C_{45})] ξ^{2} ρ ω^{2} + (C_{44} + C_{33} + C_{55}) ρ^{2} ω^{4}} / Δ, \\ A_{3} = [(C_{11} C_{55} C_{66} - C_{16}^{2} C_{55}) ξ^{6} + (C_{16}^{2} - C_{55} C_{66} - C_{11} C_{55} - C_{11} C_{66}) ξ^{4} ρ ω^{2} \\ + (C_{11} + C_{55} + C_{66}) ξ^{2} ρ^{2} ω^{4} - ρ^{3} ω^{6}] / Δ, \end{matrix}

(B1)

with

Δ = C_{33} C_{44} C_{55} - C_{33} C_{45}^{2} .

(B2)

APPENDIX C: SUBMATRICES AND VECTORS

\begin{matrix} P^{-} = [\begin{matrix} 1 & 1 & 1 \\ V_{1} & V_{3} & V_{5} \\ W_{1} & W_{3} & W_{5} \end{matrix}], P^{+} = [\begin{matrix} 1 & 1 & 1 \\ V_{1} & V_{3} & V_{5} \\ - W_{1} & - W_{3} & - W_{5} \end{matrix}], \\ D^{-} = [\begin{matrix} D_{33, 1} & D_{33, 3} & D_{33, 5} \\ D_{13, 1} & D_{13, 3} & D_{13, 5} \\ D_{23, 1} & D_{23, 3} & D_{23, 5} \end{matrix}], D^{+} = [\begin{matrix} D_{33, 1} & D_{33, 3} & D_{33, 5} \\ - D_{13, 1} & - D_{13, 3} & - D_{13, 5} \\ - D_{23, 1} & - D_{23, 3} & - D_{23, 5} \end{matrix}], \\ H = [\begin{matrix} e^{i ξ_{3, 1} x_{3}} & 0 & 0 \\ 0 & e^{i ξ_{3, 3} x_{3}} & 0 \\ 0 & 0 & e^{i ξ_{3, 5} x_{3}} \end{matrix}], U^{-} = [\begin{matrix} U_{1, 1} \\ U_{1, 3} \\ U_{1, 5} \end{matrix}], U^{+} = [\begin{matrix} U_{1, 2} \\ U_{1, 4} \\ U_{1, 6} \end{matrix}], \end{matrix}

(C1)

where the common term,

e^{i (ξ x_{1} - ω t)}

in H, is suppressed.

APPENDIX D: MATERIAL PARAMETERS

Mass density and layer thickness are

ρ = 1500 {kg / m}^{3}, d_{m} = 0.125 mm .

(D1)

Real parts of the stiffness matrix in GPa are

C_{i j}^{' re} = [\begin{matrix} 132 & 6.9 & 5.9 & 0 & 0 & 0 \\ 12.3 & 5.5 & 0 & 0 & 0 \\ 12.1 & 0 & 0 & 0 \\ 3.32 & 0 & 0 \\ sym & 6.21 & 0 \\ 6.15 \end{matrix}] .

(D2)

Imaginary parts of the stiffness matrix in GPa are

C_{i j}^{' im} = [\begin{matrix} 0.4 & 0.001 & 0.016 & 0 & 0 & 0 \\ 0.037 & 0.021 & 0 & 0 & 0 \\ 0.043 & 0 & 0 & 0 \\ 0.009 & 0 & 0 \\ sym & 0.015 & 0 \\ 0.02 \end{matrix}] .

(D3)

REFERENCES

1.

M.

Castaings

and

P.

Cawley

, “

The generation, propagation, and detection of Lamb waves in plates using air-coupled ultrasonic transducers

,”

J. Acoust. Soc. Am.

100

(

5

),

3070

–

3077

(

1996

).

https://doi.org/10.1121/1.417193

Google Scholar

Crossref

2.

M.

Castaings

,

P.

Cawley

,

R.

Farlow

, and

G.

Hayward

, “

Single sided inspection of composite materials using air coupled ultrasound

,”

J. Nondestr. Eval.

17

(

1

),

37

–

45

(

1998

).

https://doi.org/10.1023/A:1022632513303

Google Scholar

Crossref

3.

M.

Castaings

and

B.

Hosten

, “

Lamb and SH waves generated and detected by air-coupled ultrasonic transducers in composite material plates

,”

NDT&E Int.

34

,

249

–

258

(

2001

).

https://doi.org/10.1016/S0963-8695(00)00065-7

Google Scholar

Crossref

4.

M.

Castaings

and

B.

Hosten

, “

Ultrasonic guided waves for health monitoring of high-pressure composite tanks

,”

NDT&E Int.

41

,

648

–

655

(

2008

).

https://doi.org/10.1016/j.ndteint.2008.03.010

Google Scholar

Crossref

5.

P. D.

Wilcox

,

M. J. S.

Lowe

, and

P.

Cawley

, “

Mode and transducer selection for long range Lamb wave inspection

,”

J. Intell. Mater. Syst. Struct.

12

,

553

–

565

(

2001

).

https://doi.org/10.1177/10453890122145348

Google Scholar

Crossref

6.

M. J. S.

Lowe

,

D. N.

Alleyne

, and

P.

Cawley

, “

Defect detection in pipes using guided waves

,”

Ultrasonics

36

,

147

–

154

(

1998

).

https://doi.org/10.1016/S0041-624X(97)00038-3

Google Scholar

Crossref

7.

D. N.

Alleyne

,

B.

Pavlakovic

,

M. J. S.

Lowe

, and

P.

Cawley

, “

Rapid, long range inspection of chemical plant pipework using guided waves

,”

AIP Conf. Proc.

557

,

180

–

196

(

2001

).

https://doi.org/10.1063/1.1373757

Google Scholar

Crossref

8.

R.

Long

,

M. J. S.

Lowe

, and

P.

Cawley

, “

Attenuation characteristics of the fundamental modes that propagate in buried iron water pipes

,”

Ultrasonics

41

(

7

),

509

–

519

(

2003

).

https://doi.org/10.1016/S0041-624X(03)00166-5

Google Scholar

Crossref

PubMed

9.

M.

Luukkala

,

P.

Heikkila

, and

J.

Surakka

, “

Plate wave resonance—A contactless test method

,”

Ultrasonics

9

(

4

),

201

–

208

(

1971

).

https://doi.org/10.1016/0041-624X(71)90387-8

Google Scholar

Crossref

10.

M.

Luukkala

and

P.

Meriläinen

, “

Metal plate testing using airborne ultrasound

,”

Ultrasonics

11

(

5

),

218

–

221

(

1973

).

https://doi.org/10.1016/0041-624X(73)90233-3

Google Scholar

Crossref

11.

H.

Lamb

, “

On waves in an elastic plate

,”

Proc. R. Soc. London, Ser. A

93

,

114

–

128

(

1917

).

https://doi.org/10.1098/rspa.1917.0008

Google Scholar

Crossref

12.

I. A.

Viktorov

,

Rayleigh and Lamb Waves: Physical Theory and Applications

(

Plenum

,

New York

,

1967

).

Google Scholar

Crossref

13.

W. T.

Thomson

, “

Transmission of elastic waves through a stratified solid medium

,”

J. Appl. Phys.

21

,

89

–

93

(

1950

).

https://doi.org/10.1063/1.1699629

Google Scholar

Crossref

14.

N. A.

Haskell

, “

The dispersion of surface waves on multilayered media

,”

Bull. Seism. Soc. Am.

43

(

1

),

17

–

34

(

1953

).

https://doi.org/10.1785/BSSA0430010017

Google Scholar

Crossref

15.

A. H.

Nayfeh

and

D. E.

Chimenti

, “

Ultrasonic wave reflection from liquid-coupled orthotropic plates with application to fibrous composites

,”

J. Appl. Mech.

55

(

4

),

863

–

870

(

1988

).

https://doi.org/10.1115/1.3173734

Google Scholar

Crossref

16.

A. H.

Nayfeh

and

D. E.

Chimenti

, “

Free wave propagation in plates of general anisotropic media

,” in

Review of Progress in Quantitative Nondestructive Evaluation,

edited by D. O. Thompson and D. E. Chimenti (

Plenum

,

New York

1989

), pp.

181

–

188

.

https://doi.org/10.1007/978-1-4613-0817-1_23

Google Scholar

17.

A. H.

Nayfeh

, “

The propagation of horizontally polarized shear waves in multilayered anisotropic media

,”

J. Acoust. Soc. Am.

86

(

5

),

2007

–

2012

(

1989

).

https://doi.org/10.1121/1.398580

Google Scholar

Crossref

18.

A. H.

Nayfeh

, “

The general problem of elastic wave propagation in multilayered anisotropic media

,”

J. Acoust. Soc. Am.

89

(

4

),

1521

–

1531

(

1991

).

https://doi.org/10.1121/1.400988

Google Scholar

Crossref

19.

A. H.

Nayfeh

,

Wave Propagation in Layered Anisotropic Media with Applications to Composites

(

North-Holland

,

Amsterdam

,

1995

).

Google Scholar

20.

J. W.

Dunkin

, “

Computation of modal solutions in layered, elastic media at high frequencies

,”

Bull. Seism. Soc. Am.

55

(

2

),

335

–

358

(

1965

).

https://doi.org/10.1785/BSSA0550020335

Google Scholar

Crossref

21.

L.

Knopoff

, “

A matrix method for elastic wave problems

,”

Bull. Seism. Soc. Am.

54

(

1

),

431

–

438

(

1964

).

https://doi.org/10.1785/BSSA0540010431

Google Scholar

Crossref

22.

M. J. S.

Lowe

, “

Matrix techniques for modeling ultrasonic waves in multilayered media

,”

IEEE Trans. Ultrason. Ferroelec. Freq. Contr.

42

(

4

),

525

–

542

(

1995

).

https://doi.org/10.1109/58.393096

Google Scholar

Crossref

23.

H.

Schmidt

and

G.

Tango

, “

Efficient global matrix approach to the computation of synthetic seismograms

,”

Geophys. J. R. Astron. Soc.

84

,

331

–

359

(

1986

).

https://doi.org/10.1111/j.1365-246X.1986.tb04359.x

Google Scholar

Crossref

24.

H.

Schmidt

and

F. B.

Jensen

, “

Efficient numerical solution technique for wave propagation in horizontally stratified environments

,”

Comput. Math. Appl.

11

(

7-8

),

699

–

715

(

1985

).

https://doi.org/10.1016/0898-1221(85)90166-X

Google Scholar

Crossref

25.

H.

Schmidt

and

F. B.

Jensen

, “

A full wave solution for propagation in multilayered viscoelastic media with application to Gaussian beam reflection at fluid–solid interfaces

,”

J. Acoust. Soc. Am.

77

(

3

),

813

–

825

(

1985

).

https://doi.org/10.1121/1.392050

Google Scholar

Crossref

26.

A. K.

Mal

, “

Wave propagation in layered composite laminates under periodic surface loads

,”

Wave Motion

10

(

3

),

257

–

266

(

1988

).

https://doi.org/10.1016/0165-2125(88)90022-4

Google Scholar

Crossref

27.

A. K.

Mal

, “

Guided waves in layered solids with interface zones

,”

Int. J. Eng. Sci.

26

(

8

),

873

–

881

(

1988

).

https://doi.org/10.1016/0020-7225(88)90038-9

Google Scholar

Crossref

28.

E.

Kausel

and

J. M.

Roesset

, “

Stiffness matrices for layered soils

,”

Bull. Seis. Soc. Am.

71

(

6

),

1743

–

1761

(

1981

).

https://doi.org/10.1785/BSSA0710061743

Google Scholar

Crossref

29.

Y.

Wang

and

R. K. N. D.

Rajapakse

, “

An exact stiffness method for elastodynamics of a layered orthotropic half-plane

,”

J. Appl. Mech.

61

(

2

),

339

–

347

(

1994

).

https://doi.org/10.1115/1.2901450

Google Scholar

Crossref

30.

L.

Wang

and

S. I.

Rokhlin

, “

Time resolved line focus acoustic microscopy of composites

,” in

Review of Progress in Quantitative Nondestructive Evaluation

, edited by D. O. Thompson and D. E. Chimenti (

Plenum

,

New York

, 1999), pp.

1321

–

1328

.

https://doi.org/10.1007/978-1-4615-4791-4_169

Google Scholar

Crossref

31.

L.

Wang

and

S. I.

Rokhlin

, “

Analysis of ultrasonic wave propagation in multiply composites: Homogenization and effective anisotropic media

,” in

Review of Progress in Quantitative Nondestructive Evaluation

, edited by D. O. Thompson and D. E. Chimenti (

Plenum

,

New York

, 2001), pp.

1015

–

1022

, available at https://pubs.aip.org/aip/acp/article/557/1/1015/579859/Analysis-of-ultrasonic-wave-propagation-in.

Crossref

32.

B. L. N.

Kennett

,

Seismic Wave Propagation in Stratified Media

(

Cambridge University Press

,

New York

,

1983

).

Google Scholar

33.

B. L. N.

Kennett

and

N. J.

Kerry

, “

Seismic waves in a stratified half space

,”

Geophys. J. Int.

57

(

3

),

557

–

583

(

1979

).

https://doi.org/10.1111/j.1365-246X.1979.tb06779.x

Google Scholar

Crossref

34.

G. J.

Fryer

and

L. N.

Frazer

, “

Seismic waves in stratified anisotropic media

,”

Geophys. J. Int.

78

(

3

),

691

–

710

(

1984

).

https://doi.org/10.1111/j.1365-246X.1984.tb05065.x

Google Scholar

Crossref

35.

G. J.

Fryer

and

L. N.

Frazer

, “

Seismic waves in stratified anisotropic media—II. Elastodynamic eigensolutions for some anisotropic systems

,”

Geophys. J. Int.

91

(

1

),

73

–

101

(

1987

).

https://doi.org/10.1111/j.1365-246X.1987.tb05214.x

Google Scholar

Crossref

36.

S. I.

Rokhlin

and

W.

Huang

, “

Ultrasonic wave interaction with a thin anisotropic layer between two anisotropic solids: Exact and asymptotic-boundary-condition methods

,”

J. Acoust. Soc. Am.

92

(

3

),

1729

–

1742

(

1992

).

https://doi.org/10.1121/1.403912

Google Scholar

Crossref

37.

S.

Pant

,

J.

Laliberte

,

M.

Martinez

, and

B.

Rocha

, “

Derivation and experimental validation of Lamb wave equations for an n-layered anisotropic composite laminate

,”

Compos. Struct.

111

,

566

–

579

(

2014

).

https://doi.org/10.1016/j.compstruct.2014.01.034

Google Scholar

Crossref

38.

S.

Guo

,

M.

Rébillat

, and

N.

Mechbal

, “

Prediction of frequency and spatially dependent attenuation of guided waves propagating in mounted and unmounted A380 parts made up of anisotropic viscoelastic composite laminates

,”

Struct. Health Monit.

22

(

2

),

1326

–

1352

(

2023

).

https://doi.org/10.1177/14759217221099967

Google Scholar

Crossref

39.

S. I.

Rokhlin

and

L.

Wang

, “

Stable recursive algorithm for elastic wave propagation in layered anisotropic media: Stiffness matrix method

,”

J. Acoust. Soc. Am.

112

(

3

),

822

–

834

(

2002

).

https://doi.org/10.1121/1.1497365

Google Scholar

Crossref

PubMed

40.

L.

Wang

and

S. I.

Rokhlin

, “

Stable reformulation of transfer matrix method for wave propagation in layered anisotropic media

,”

Ultrasonics

39

,

413

–

424

(

2001

).

https://doi.org/10.1016/S0041-624X(01)00082-8

Google Scholar

Crossref

PubMed

41.

S. I.

Rokhlin

,

D. E.

Chimenti

, and

P. B.

Nagy

,

Physical Ultrasonics of Composites

(

Oxford University Press

,

Oxford, UK

,

2011

).

Google Scholar

42.

E. L.

Tan

, “

Stiffness matrix method with improved efficiency for elastic wave propagation in layered anisotropic media

,”

J. Acoust. Soc. Am.

118

(

6

),

3400

–

3403

(

2005

).

https://doi.org/10.1121/1.2118287

Google Scholar

Crossref

43.

E. L.

Tan

, “

Hybrid compliance-stiffness matrix method for stable analysis of elastic wave propagation in multilayered anisotropic media

,”

J. Acoust. Soc. Am.

119

(

1

),

45

–

53

(

2006

).

https://doi.org/10.1121/1.2139617

Google Scholar

Crossref

PubMed

44.

V. G. A.

Kamal

and

V.

Giurgiutiu

, “

Stiffness Transfer Matrix Method (STMM) for stable dispersion curves solution in anisotropic composites

,” in

Proc. SPIE

, San Diego, CA (SPIE, Bellingham, WA,

2014

).

Google Scholar

45.

M.

Barski

and

P.

Pajak

, “

An application of stiffness matrix method to determining of dispersion curves for arbitrary composite materials

,”

J. KONES Powertrain Transp.

23

(

1

),

47

–

54

(

2016

).

https://doi.org/10.5604/12314005.1213520

Google Scholar

Crossref

46.

M.

Barski

and

P.

Pajak

, “

Determination of dispersion curves for composite materials with the use of stiffness matrix method

,”

Acta Mech. Autom.

11

(

2

),

121

–

128

(

2017

).

https://doi.org/10.1515/ama-2017-0019

Google Scholar

Crossref

47.

M.

Barski

and

A.

Muc

, “

The impact of the fiber orientation angle on the dispersion characteristic of guided waves for multilayered GLARE composites

,”

J. Phys. Conf. Ser.

1603

,

012001

(

2020

).

https://doi.org/10.1088/1742-6596/1603/1/012001

Google Scholar

Crossref

48.

A.

Muc

,

M.

Barski

,

A.

Stawiarski

,

M.

Chwał

, and

M.

Augustyn

, “

Dispersion curves and identification of elastic wave modes for fiber metal laminate

,”

Compos. Struct.

255

,

112930

(

2021

).

https://doi.org/10.1016/j.compstruct.2020.112930

Google Scholar

Crossref

49.

A. M. A.

Huber

and

M. G. R.

Sause

, “

Classification of solutions for guided waves in anisotropic composites with large numbers of layers

,”

J. Acoust. Soc. Am.

144

(

1

),

3236

–

3251

(

2018

).

https://doi.org/10.1121/1.5082299

Google Scholar

Crossref

PubMed

50.

A.

Huber

, “

Numerical modeling of guided waves in anisotropic composites with application to air-coupled ultrasonic inspection

,” Ph.D. thesis,

University of Augsburg

,

Augsburg, Germany

,

2020

.

Google Scholar

51.

L.

Gavric

, “

Computation and propagative waves in free rail using a finite element technique

,”

J. Sound Vib.

185

(

3

),

531

–

543

(

1995

).

https://doi.org/10.1006/jsvi.1995.0398

Google Scholar

Crossref

52.

T.

Hayashi

,

W.-J.

Song

, and

J. L.

Rose

, “

Guided wave dispersion curves for a bar with an arbitrary cross-section, a rod and rail example

,”

Ultrasonics

41

(

3

),

175

–

183

(

2003

).

https://doi.org/10.1016/S0041-624X(03)00097-0

Google Scholar

Crossref

PubMed

53.

M.

Castaings

and

M. J. S.

Lowe

, “

Finite element model for waves guided along solid systems of arbitrary section coupled to infinite solid media

,”

J. Acoust. Soc. Am.

123

(

2

),

696

–

708

(

2008

).

https://doi.org/10.1121/1.2821973

Google Scholar

Crossref

PubMed

54.

Z.

Fan

,

M. J. S.

Lowe

,

M.

Castaings

, and

C.

Bacon

, “

Torsional waves propagation along a waveguide of arbitrary cross section immersed in a perfect fluid

,”

J. Acoust. Soc. Am.

124

(

4

),

2002

–

2010

(

2008

).

https://doi.org/10.1121/1.2968677

Google Scholar

Crossref

PubMed

55.

Z.

Fan

and

M. J. S.

Lowe

, “

Elastic waves guided by a welded joint in a plate

,”

Proc. R. Soc. London A

465

,

2053

–

2068

(

2009

).

https://doi.org/10.1098/rspa.2009.0010

Google Scholar

Crossref

56.

I.

Bartoli

,

A.

Marzani

,

F. L.

di Scalea

, and

E.

Viola

, “

Modeling wave propagation in damped waveguides of arbitrary cross-section

,”

J. Sound Vib.

295

(

3-5

),

685

–

707

(

2006

).

https://doi.org/10.1016/j.jsv.2006.01.021

Google Scholar

Crossref

57.

A.

Marzani

,

E.

Viola

,

I.

Bartoli

,

F. L.

di Scalea

, and

P.

Rizzo

, “

A semi-analytical finite element formulation for modeling stress wave propagation in axisymmetric damped waveguides

,”

J. Sound Vib.

318

(

3

),

488

–

505

(

2008

).

https://doi.org/10.1016/j.jsv.2008.04.028

Google Scholar

Crossref

58.

S.

Sorohan

,

N.

Constantin

,

M.

Găvan

, and

V.

Anghel

, “

Extraction of dispersion curves for waves propagating in free complex waveguides by standard finite element codes

,”

Ultrasonics

51

(

4

),

503

–

515

(

2011

).

https://doi.org/10.1016/j.ultras.2010.12.003

Google Scholar

Crossref

PubMed

59.

D.

Pereira

,

J.

Fernandes

, and

P.

Belanger

, “

Ex vivo assessment of cortical bone properties using low-frequency ultrasonic guided waves

,”

IEEE Trans. Ultrason. Ferroelectr. Freq. Control

67

(

5

),

910

–

922

(

2020

).

https://doi.org/10.1109/TUFFC.2019.2958035

Google Scholar

Crossref

PubMed

60.

F.

Seyfaddini

,

H.

Nguyen-Xuan

, and

V.-H.

Nguyen

, “

Wave dispersion analysis of three-dimensional vibroacoustic waveguides with semi-analytical isogeometric method

,”

Comput. Methods Appl. Mech. Eng.

385

,

114043

(

2021

).

https://doi.org/10.1016/j.cma.2021.114043

Google Scholar

Crossref

61.

F.

Treysséde

, “

A model reduction method for fast finite element analysis of continuously symmetric waveguides

,”

J. Sound Vib.

508

,

116204

(

2021

).

https://doi.org/10.1016/j.jsv.2021.116204

Google Scholar

Crossref

62.

P. P.

Delsanto

,

R. S.

Schechter

, and

R. B.

Mignogna

, “

Connection machine simulation of ultrasonic wave propagation in materials III: The three-dimensional case

,”

Wave Motion

26

(

4

),

329

–

339

(

1997

).

https://doi.org/10.1016/S0165-2125(97)00013-9

Google Scholar

Crossref

63.

M.

Ruzzene

,

S. M.

Jeong

,

T. E.

Michaels

,

J. E.

Michaels

, and

B.

Mi

, “

Simulation and measurement of ultrasonic waves in elastic plates using laser vibrometry

,” in

Review of Progress in Quantitative Nondestructive Evaluation

, edited by D. O. Thompson and D. E. Chimenti (

Plenum

,

New York

, 2005), pp.

172

–

179

, available at https://pubs.aip.org/aip/acp/article/760/1/172/945089/Simulation-and-Measurement-of-Ultrasonic-Waves-in.

Crossref

64.

K. S.

Nadella

and

C. E. S.

Cesnik

, “

Local interaction simulation approach for modeling wave propagation in composite structures

,”

CEAS Aeronaut. J.

4

(

1

),

35

–

48

(

2013

).

https://doi.org/10.1007/s13272-012-0061-9

Google Scholar

Crossref

65.

A. T. I.

Adamou

and

R. V.

Craster

, “

Spectral methods for modelling guided waves in elastic media

,”

J. Acoust. Soc. Am.

116

(

3

),

1524

–

1535

(

2004

).

https://doi.org/10.1121/1.1777871

Google Scholar

Crossref

66.

F.

Karpfinger

,

B.

Gurevich

, and

A.

Bakulin

, “

Modeling of wave dispersion along cylindrical structures using the spectral method

,”

J. Acoust. Soc. Am.

124

(

2

),

859

–

865

(

2008

).

https://doi.org/10.1121/1.2940577

Google Scholar

Crossref

PubMed

67.

F.

Karpfinger

,

B.

Gurevich

, and

A.

Bakulin

, “

Modeling of axisymmetric wave modes in a poroelastic cylinder using spectral method

,”

J. Acoust. Soc. Am.

124

(

4

),

EL230

–

EL235

(

2008

).

https://doi.org/10.1121/1.2968303

Google Scholar

Crossref

PubMed

68.

F.

Karpfinger

,

H.-P.

Valero

,

B.

Gurevich

,

A.

Bakulin

, and

B.

Sinha

, “

Spectral-method algorithm for modeling dispersion of acoustic modes in elastic cylindrical structures

,”

Geophysics

75

(

3

),

H19

–

H27

(

2010

).

https://doi.org/10.1190/1.3380590

Google Scholar

Crossref

69.

F.

Karpfinger

,

B.

Gurevich

,

H.-P.

Valero

,

A.

Bakulin

, and

B.

Sinha

, “

Tube wave signatures in cylindrically layered poroelastic media computed with spectral method

,”

Geophys. J. Int.

183

(

2

),

1005

–

1013

(

2010

).

https://doi.org/10.1111/j.1365-246X.2010.04773.x

Google Scholar

Crossref

70.

B.

Yu

,

S.

Yang

,

C.

Gan

, and

H.

Lei

, “

A new procedure for exploring the dispersion characteristics of longitudinal guided waves in a multi-layered tube with a weak interface

,”

J. Nondestruct. Eval.

32

(

3

),

263

–

276

(

2013

).

https://doi.org/10.1007/s10921-013-0179-7

Google Scholar

Crossref

71.

T.

Zharnikov

,

D.

Syresin

, and

C.-J.

Hsu

, “

Application of the spectral method for computation of the spectrum of anisotropic waveguides

,”

Proc. Mtgs. Acoust.

19

,

045065

(

2013

).

https://doi.org/10.1121/1.4799396

Google Scholar

Crossref

72.

F. H.

Quintanilla

,

M. J. S.

Lowe

, and

R. V.

Craster

, “

Modeling guided elastic waves in generally anisotropic media using a spectral collocation method

,”

J. Acoust. Soc. Am.

137

(

3

),

1180

–

1194

(

2015

).

https://doi.org/10.1121/1.4913777

Google Scholar

Crossref

PubMed

73.

F. H.

Quintanilla

,

Z.

Fan

,

M. J. S.

Lowe

, and

R. V.

Craster

, “

Guided waves' dispersion curves in anisotropic viscoelastic single- and multi-layered media

,”

Proc. R. Soc. A

471

(

2183

),

20150268

(

2015

).

https://doi.org/10.1098/rspa.2015.0268

Google Scholar

Crossref

74.

F. H.

Quintanilla

,

M. J. S.

Lowe

, and

R. V.

Craster

, “

Full 3D dispersion curve solutions for guided waves in generally anisotropic media

,”

J. Sound Vib.

363

(

17

),

545

–

559

(

2016

).

https://doi.org/10.1016/j.jsv.2015.10.017

Google Scholar

Crossref

75.

F. H.

Quintanilla

,

M. J. S.

Lowe

, and

R. V.

Craster

, “

The symmetry and coupling properties of solutions in general anisotropic multilayer waveguides

,”

J. Acoust. Soc. Am.

141

(

1

),

406

–

418

(

2017

).

https://doi.org/10.1121/1.4973543

Google Scholar

Crossref

PubMed

76.

D.

Kiefer

,

M.

Ponschab

,

S.

Rupitsch

, and

M.

Mayle

, “

Calculating the full leaky Lamb wave spectrum with exact fluid interaction

,”

J. Acoust. Soc. Am.

145

(

6

),

3341

–

3350

(

2019

).

https://doi.org/10.1121/1.5109399

Google Scholar

Crossref

PubMed

77.

E.

Georgiades

,

M.

Lowe

, and

R.

Craster

, “

Leaky wave characterisation using spectral methods

,”

J. Acoust. Soc. Am.

152

(

3

),

1487

–

1497

(

2022

).

https://doi.org/10.1121/10.0013897

Google Scholar

Crossref

PubMed

78.

B.

Pavlakovic

,

M.

Lowe

,

D.

Alleyne

, and

P.

Cawley

, “Disperse: A general purpose program for creating dispersion curves,” in

Review of Progress in Quantitative Nondestructive Evaluation

, edited by D. O. Thompson and D. E. Chimenti (

Plenum

,

New York

, 1997), pp.

185

–

192

, available at https://link.springer.com/chapter/10.1007/978-1-4615-5947-4_24.

79.

M. J. S.

Lowe

and

B.

Pavlakovic

, “

DISPERSE

,” available at http://www.imperial.ac.uk/non-destructive-evaluation/products-and-services/disperse (Last viewed June 26,

2023

).

Google Scholar

80.

P.

Boccini

,

A.

Marzani

, and

E.

Viola

, “

Graphical user interface for guided acoustic waves

,”

J. Comput. Civil. Eng.

25

(

3

),

202

–

210

(

2011

).

https://doi.org/10.1061/(ASCE)CP.1943-5487.0000081

Google Scholar

Crossref

81.

A.

Marzani

and

P

Bocchini

, “

GUIGUW

,” available at http://www.guiguw.com/ (Last viewed June 26,

2023

).

Google Scholar

82.

A. M. A.

Huber

, “

Dispersion Calculator (DC)

,” available at https://github.com/ArminHuber/Dispersion-Calculator (Last viewed June 26,

2023

).

Google Scholar

83.

D. R.

Ramasawmy

,

B. T.

Cox

, and

B. E.

Treeby

, “

ElasticMatrix: A MATLAB toolbox for anisotropic elastic wave propagation in layered media

,”

SoftwareX

11

(

100397

),

100397

(

2020

).

https://doi.org/10.1016/j.softx.2020.100397

Google Scholar

Crossref

84.

D. R.

Ramasawmy

,

B. T.

Cox

, and

B. E.

Treeby

, “

ElasticMatrix

,” available at https://github.com/dannyramasawmy/ElasticMatrix (Last viewed June 26,

2023

).

Google Scholar

85.

A. H.

Orta

,

M.

Kersemans

, and

K. V. D.

Abeele

, “

A comparative study for calculating dispersion curves in viscoelastic multi-layered plates

,”

Compos. Struct.

294

,

115779

(

2022

).

https://doi.org/10.1016/j.compstruct.2022.115779

Google Scholar

Crossref

86.

A. H.

Orta

,

K.

Van Den Abeele

, and

M.

Kersemans

, “

The Dispersion Box

,” available at https://github.com/adilorta/The-Dispersion-Box (Last viewed June 26,

2023

).

Google Scholar

87.

A. H.

Orta

,

J.

Vandendriessche

,

M.

Kersemans

,

W. V.

Paepegem

,

N. B.

Roozen

, and

K. V. D.

Abeele

, “

Modeling Lamb wave propagation in visco-elastic composite plates using a fifth-order plate theory

,”

Ultrasonics

116

,

106482

(

2021

).

https://doi.org/10.1016/j.ultras.2021.106482

Google Scholar

Crossref

PubMed

88.

D. A.

Kiefer

, “

GEW dispersion script

,” available at https://doi.org/10.5281/zenodo.7010603 (Last viewed June 26,

2023

).

Google Scholar

89.

E.

Ducasse

and

M.

Deschamps

, “

Mode computation of immersed multilayer plates by solving an eigenvalue problem

,”

Wave Motion

112

,

102962

(

2022

).

https://doi.org/10.1016/j.wavemoti.2022.102962

Google Scholar

Crossref

90.

A.

Huber

, “

Non-destructive testing of future rocket boosters using air-coupled ultrasound

,” in

19th World Conference on Non-Destructive Testing

, Munich, Germany (NDT.net, Germany,

2016

), pp.

1

–

9

.

Google Scholar

91.

Y.

Li

and

R.

Thompson

, “

Influence of anisotropy on the dispersion characteristics of guided ultrasonic plate modes

,”

J. Acoust. Soc. Am.

87

(

5

),

1911

–

1931

(

1990

).

https://doi.org/10.1121/1.399318

Google Scholar

Crossref

92.

S. I.

Rokhlin

,

D. E.

Chimenti

, and

A. H.

Nayfeh

, “

On the topology of the complex wave spectrum in a fluid–coupled elastic layer

,”

J. Acoust. Soc. Am.

85

(

3

),

1074

–

1108

(

1989

).

https://doi.org/10.1121/1.397490

Google Scholar

Crossref

93.

D. E.

Chimenti

and

A. H.

Nayfeh

, “

Anomalous ultrasonic dispersion in fluid–coupled, fibrous composite plates

,”

J. Appl. Phys.

49

(

9

),

492

–

493

(

1986

).

https://doi.org/10.1063/1.97127

Google Scholar

Crossref

94.

L.

Cremer

, “Theorie der Schalldämmung dünner Wände bei schrägem Einfall” (“

Theory of sound isolation by thin walls at oblique incidence

”),

Akust. Z.

7

,

81

–

104

(

1942

) (in German).

Google Scholar

95.

A.

Schoch

, “Der Schalldurchgang durch Platten” (“

Sound transmission through plates

”),

Acoustica

2

,

1

–

18

(

1952

) (in German), available at https://www.ingentaconnect.com/content/dav/aaua/1952/00000002/00000001/art00003#.

Google Scholar

96.

M. J. S.

Lowe

, “

Plate waves for the NDT of diffusion bonded titanium

,” Ph.D. thesis,

University of London

, London,

1992

.

Google Scholar

97.

N.

Bochud

,

J.

Laurent

,

F.

Bruno

,

D.

Royer

, and

C.

Prada

, “

Towards real-time assessment of anisotropic plate properties using elastic guided waves

,”

J. Acoust. Soc. Am.

143

(

2

),

1138

–

1147

(

2018

).

https://doi.org/10.1121/1.5024353

Google Scholar

Crossref

PubMed

Classification of solutions for guided waves in fluid-loaded viscoelastic composites with large numbers of layers

I. INTRODUCTION

II. NUMERICAL MODEL

A. Field equations and geometrical setup

B. The SMM

1. Coupled polarization

2. Decoupled polarization

a. Lamb waves.

b. Shear horizontal waves.

C. Viscoelasticity and fluid-loading

D. Dispersion curve tracing

E. Mode shape calculation

III. CLASSIFICATION OF SOLUTIONS

IV. NUMERICAL EXAMPLES

A. Sixteen layers composite

B. Four hundred layers composite

V. DISCUSSION

VI. CONCLUSIONS

ACKNOWLEDGMENTS

APPENDIX A: TRANSFORMED STIFFNESS MATRIX

APPENDIX B: POLYNOMIAL COEFFICIENTS

APPENDIX C: SUBMATRICES AND VECTORS

APPENDIX D: MATERIAL PARAMETERS

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Classification of solutions for guided waves in fluid-loaded viscoelastic composites with large numbers of layers Open Access

I. INTRODUCTION

II. NUMERICAL MODEL

A. Field equations and geometrical setup

B. The SMM

1. Coupled polarization

2. Decoupled polarization

a. Lamb waves.

b. Shear horizontal waves.

C. Viscoelasticity and fluid-loading

D. Dispersion curve tracing

E. Mode shape calculation

III. CLASSIFICATION OF SOLUTIONS

IV. NUMERICAL EXAMPLES

A. Sixteen layers composite

B. Four hundred layers composite

V. DISCUSSION

VI. CONCLUSIONS

ACKNOWLEDGMENTS

APPENDIX A: TRANSFORMED STIFFNESS MATRIX

APPENDIX B: POLYNOMIAL COEFFICIENTS

APPENDIX C: SUBMATRICES AND VECTORS

APPENDIX D: MATERIAL PARAMETERS

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only

Classification of solutions for guided waves in fluid-loaded viscoelastic composites with large numbers of layers