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.

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 Lamb11 in 1917, dispersion diagrams of the guided waves named after him can be calculated in the most efficient way since then.12 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 Thomson13 in 1950 and edited by Haskell14 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 Chimenti15,16 and Nayfeh17,18 extended TMM to generally anisotropic multilayered media. This work is well-documented in the excellent book by Nayfeh.19 However, Dunkin20 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 Knopoff21 in 1964. A well-organized and instructive discussion of TMM and GMM is provided by Lowe.22 In GMM, a large single matrix comprises the equations of all layers. Schmidt and Tango23 and Schmidt and Jensen24,25 showed that this method is numerically stable if implemented properly. Mal26,27 extended GMM to anisotropic multilayered media. Kausel and Roesset28 conducted a reformulation for isotropic cases and Wang and Rajapakse29 conducted a reformulation for an orthotropic layer in the plane of symmetry. Wang and Rokhlin30,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 Kennett32 and Kennett and Kerry,33 which was subsequently extended to generally anisotropic media by Fryer and Frazer34,35 and obtained later by Rokhlin and Huang.36 In 2014, Pant et al.37 applied GMM to a multilayered composite laminate, and in 2023, Guo et al.38 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 Wang39,40 in 2001, and well-documented in the book of Rokhlin et al.41 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,42 most notably the hybrid compliance-stiffness matrix method (HCSMM).43 SMM was used by researchers like Kamal and Giurgiutiu,44 Barski and Pajak,45,46 Barski and Muc,47 and Muc et al.48 for the calculation of guided wave dispersion curves for multilayered composites. Huber and Sause49,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 Gavric51 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.52 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.53 Torsional waves propagating along arbitrary cross sectional waveguides immersed in a fluid54 and elastic waves propagating along weld joints between plates55 were modeled by Fan et al. Bartoli et al. and Marzani et al. investigated guided waves in viscoelastic waveguides of rectangular, arbitrary,56 and axissymmetric57 cross sections. Another study of dispersion curve calculation for complex waveguides was conducted by Sorohan et al.58 SAFE has also reached the field of clinical research as Pereira et al.59 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.,60 who proposed a semi-analytical isogeometric analysis (SAIGA), and Treysséde,61 who introduced a reduction strategy.

The local interaction simulation approach (LISA), a method studied by Delsanto et al.62 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 Craster65 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 structures66 and, later, for axisymmetric wave modes in a porous elastic cylinder67 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.70 Zharnikov et al.71 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.74 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.75 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 Pavlakovic78,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 Huber82 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.87 as an improvement to the previous lower order plate theories. Also, in 2022, Kiefer88 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.90 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.

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.

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 = 1 2 ( u l x k + u k x l ) ,
(1)
σ i j x j = ρ 2 u i 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. (3) can be written as
[ σ 11 σ 22 σ 33 σ 23 σ 13 σ 12 ] = [ 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 ] [ ε 11 ε 22 ε 33 2 ε 23 2 ε 13 2 ε 12 ] .
(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, dm, 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 x1-x3-plane. With respect to the global coordinate system, the local coordinate systems are yielded by a counterclockwise rotation through an angle, Φ m, between x1 and x 1 , m about the x3-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.19 The resulting components of Cij are given in  Appendix A. Notice that Cij has the form of monoclinic material symmetry, in general, having 13 independent elastic constants. For Φ = 0 , 90 , the elements C16, C26, C36, and C45 vanish, returning orthotropic material symmetry. This will be of importance when the decoupling into pure Lamb and SH waves is studied.

FIG. 1.

Single composite layer (left) with local (crystallographic) coordinate system x i and global coordinate system xi, i = 1,2,3, and layered composite plate (right) with [0/90/0] orientation with respect to the x i -axis.

FIG. 1.

Single composite layer (left) with local (crystallographic) coordinate system x i and global coordinate system xi, i = 1,2,3, and layered composite plate (right) with [0/90/0] orientation with respect to the x i -axis.

Close modal
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,91 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. (1)–(3), which are given in the global coordinate system by
ρ 2 u i t 2 = c ijkl 2 u l x j x k .
(5)
While the bulk waves in an infinite medium are unconditionally described by Eq. (5), 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 = u ¯ i , σ i 3 = σ ¯ i 3 , i = 1 , 2 , 3 ,
(6)
where ui, σ i 3 correspond to the originating layer, and u ¯ i , σ ¯ 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)
The task is to solve Eq. (5) for the propagation directions and polarization vectors of the bulk waves for which all share the same wavenumber component, ξ ( = ω / c p) along x1, as required by Snell's law. Here, c p is the phase velocity in the x1-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 x3-component of the wavenumber ξ3 directly, instead of searching for the ratio of the bulk waves' wavenumber components in the x3 and x1-directions, α. This change is necessary since ξ turns complex once attenuation is introduced. The sought polarization vector is represented by the bulk wave amplitudes, Ui.

1. Coupled polarization

Considering the general case where the bulk waves contribute to the displacement field with nonzero polarization components in all three directions xi, Eq. (8) 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. (9) into Eq. (5) in global coordinates yields three homogeneous linear equations,
[ 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 ] [ U 1 U 2 U 3 ] = 0 ,
(10)
where the 3 × 3 matrix is the Christoffel matrix, Mij. By vanishing the determinant of Mij, 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 A1, A2, and A3 given in  Appendix B. Equation (11) has six solutions, ξ 3 , q , q = 1 , 2 , , 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. (10) 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 = 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 = 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. (3) as
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 ) .
(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 ) = q = 1 6 ( 1 , V q , W q ) U 1 , q e i ξ 3 , q x 3 ,
(16)
( σ 33 , σ 13 , σ 23 ) = 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. (16) and (17) 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 ±,
[ u m σ m ] = [ P P + H D D + H ] m [ U m U m + ] ,
(18)
[ u m + 1 σ m + 1 ] = [ P H P + D H D + ] m [ U m U m + ] ,
(19)
where the tensor notation introduced by Rokhlin and Wang39,40 is used. The submatrices, P ± , D ±, H, and the vectors, U ±, are expanded in  Appendix C. By eliminating U m ± in Eqs. (18) and (19), 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,
[ u m + 1 σ m + 1 ] = [ P H P + D H D + ] m [ P P + H D D + H ] m 1 [ u m σ m ] = A m [ u m σ m ] ,
(20)
where Am 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. (20) 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. (18) and (19) 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,
[ σ m σ m + 1 ] = [ D D + H D H D + ] m [ P P + H P H P + ] m 1 [ u m u m + 1 ] = K m [ u m u m + 1 ] .
(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,
[ σ 1 σ 2 ] = [ K 11 A K 12 A K 21 A K 22 A ] [ u 1 u 2 ] , [ σ 2 σ 3 ] = [ K 11 B K 12 B K 21 B K 22 B ] [ u 2 u 3 ] ,
(22)
Rokhlin and Wang's recursive algorithm is used to combine them,
[ σ 1 σ 3 ] = [ 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 ] [ u 1 u 3 ] ,
(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. (23) 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. (23) 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
[ σ 1 σ 3 ] = [ 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 ) ] [ u 1 u 3 ] ,
(24)
with Y = ( K 22 A ° I K 22 A ) 1, the Hadamard operator “ °,” and
I = [ 1 1 1 1 1 1 1 1 1 ] .
(25)
Finally, by invoking the stress-free boundary conditions from Eq. (7),
[ σ 0 σ d ] = K [ u 0 u d ] = 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. (27) 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
[ 0 0 0 σ 33 d / 2 σ 13 d / 2 σ 23 d / 2 ] = [ 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 ] half [ u 1 0 u 2 0 u 3 0 u 1 d / 2 u 2 d / 2 u 3 d / 2 ] .
(30)
Then, by using the symmetric boundary conditions, Eq. (28), on Eq. (30), the dispersion equation for symmetric modes is obtained as
det K S = det [ 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 ] half = 0 ,
(31)
and by using the antisymmetric boundary conditions, Eq. (29), on Eq. (30), the dispersion equation for antisymmetric modes reads
det K A = det [ 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 ] half = 0.
(32)
The benefit of Eqs. (31) and (32) in terms of efficiency is that Eq. (24) 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. (27). 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 (31) and (32) are used in DC since version 2.4. In earlier versions, symmetric and antisymmetric modes were obtained using Eq. (27) 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 (x1-x3) decouples from polarization in the shear horizontal direction x2. 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 x1-x3-plane, Eq. (9) reduces to
( u 1 , u 3 ) = ( U 1 , U 3 ) e i ξ 3 x 3 ,
(33)
and Eq. (10) reduces to
[ C 11 ξ 2 + C 55 ξ 3 2 ρ ω 2 ( C 13 + C 55 ) ξ ξ 3 sym C 55 ξ 2 + C 33 ξ 3 2 ρ ω 2 ] [ U 1 U 3 ] = 0 ,
(34)
where the 2 × 2 matrix is the Christoffel matrix Mij. 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,
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 ,
(36)
is obtained. Equation (35) admits four solutions, ξ 3 , q , q = 1 , 2 , 3 , 4, given by
ξ 3 , q = ± A 2 ± 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. (34) delivers the amplitude ratios, W q = U 3 , q / U 1 , q, namely,
W q = m 11 ( ξ 3 , q ) m 12 ( ξ 3 , q ) = ρ ω 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
D 33 , q = i ( C 13 ξ + C 33 W q ξ 3 , q ) , D 13 , q = i C 55 ( ξ 3 , q + W q ξ ) .
(40)
By superposing the bulk waves, the displacement and stress fields can be written as
( u 1 , u 3 ) = q = 1 4 ( 1 , W q ) U 1 , q e i ξ 3 , q x 3 ,
(41)
( σ 33 , σ 13 ) = 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. (21)–(27) in the reduced form with
I = [ 1 1 1 1 ] .
(43)
The dispersion equation for symmetric modes is given by
det K S = det [ k 11 k 12 k 13 k 21 k 22 k 23 k 41 k 42 k 43 ] half = 0 ,
(44)
and that for antisymmetric modes is given by
det K A = det [ k 11 k 12 k 14 k 21 k 22 k 24 k 31 k 32 k 34 ] half = 0.
(45)
b. Shear horizontal waves.
SH waves consist of two shear bulk waves, which possess polarization only in the x2-direction such that the displacement and stress fields are given by
u 2 = q = 1 2 U 2 , q e i ξ 3 , q x 3 ,
(46)
σ 23 = q = 1 2 D 23 , q U 2 , q e i ξ 3 , q x 3 ,
(47)
with
ξ 3 , q = ± ρ ω 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. (20) can be manipulated into the convenient form19 
[ u 2 σ 23 ] m + 1 = [ cos γ m i / D m sin γ m i D m sin γ m cos γ m ] [ u 2 σ 23 ] m = A m [ u 2 σ 23 ] m ,
(49)
where γ m = ξ 3 , 1 , m d m and D m = D 23 , 1 , m. The layer transfer matrix, Am in Eq. (49), contains the matrix product required in Eq. (20) 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 Am, 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 A m ,
(50)
where A relates the displacements and stresses at the top of the whole plate to those at its bottom, according to
[ u 2 σ 23 ] d = A [ u 2 σ 23 ] 0 = [ a 11 a 12 a 21 a 22 ] [ u 2 σ 23 ] 0 .
(51)
Consider two exemplary composite layups with the configurations [0/90]4 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]4 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. (24) has to be imposed on A, i.e., the transfer matrix of the upper half of the symmetric layup given by Eq. (53) is sufficient. Referring to Eq. (51), 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. (54) and (55) do not apply to the symmetric and antisymmetric modes in single layers. Applying stress-free boundary conditions to Eq. (49) leads to the dispersion equation
i D sin γ = 0.
(56)
Because sin γ = 0 for γ = n π with n = 1 , 2 , 3 , , one has ξ 3 d = n π. Substituting this into Eq. (48), and ξ = ω / c p and ω = 2 π f, one reaches at
c p = C 66 ρ C 44 ( n / ( 2 f d ) ) 2 .
(57)
Symmetric modes are obtained for n = 0 , 2 , 4 , and antisymmetric modes are obtained for n = 1 , 3 , 5 , . Using n = 0 in Eq. (57) shows that the fundamental SH mode is nondispersive at the shear wave speed c p = C 66 / ρ. The appreciable advantage of Eq. (57) 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. (57) is only valid for real wavenumbers. In case of attenuation, ξ turns complex, and the relation ξ = ω / c p does not hold anymore. Therefore, Eq. (57) is written in the form
ξ = ρ ω 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 = ω ξ re , α = ξ im .
(59)

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,73 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, α1 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 = π C 11 im C 11 re , c p = C 11 re ρ = ω ξ 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. (8) 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. (61) 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. (27) 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.92 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.93 What happens with increasing ρ is that the fluid begins to modify the boundary conditions. The free wave boundary conditions in Eq. (7) 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 = u ¯ 3 0 , d , σ 33 0 , d = σ ¯ 33 0 , d ,
(63)
and by the vanishing of only the shear stresses,
σ i 3 0 , d = 0 , i = 1 , 2.
(64)
Following Nayfeh,19 wave motion in a nonviscous fluid is described in the global coordinate system by specializing the field equations [Eqs. (2) and (3)] to
σ i j x j = ρ f 2 u i t 2 ,
(65)
σ i j = λ f u k 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 u2 vanishes. Through combination of Eqs. (65) and (66), the expanded forms,
2 u 1 x 1 2 + 2 u 3 x 1 x 3 = 1 v f 2 2 u 1 t 2 ,
(67)
2 u 1 x 1 x 3 + 2 u 3 x 3 2 = 1 v f 2 2 u 3 t 2 ,
(68)
can be derived, where v f = λ f / ρ f is the wave speed in the fluid. Using the general solution [Eq. (8)] with Eqs. (67) and (68) delivers
[ ξ 2 ω 2 / v f 2 ξ ξ 3 sym ξ 3 2 ω 2 / v f 2 ] [ U 1 U 3 ] = 0 ,
(69)
and vanishing the Christoffel matrix determinant yields two solutions, ξ 3 , q, q = 1,2,
ξ 3 , q = ± ω 2 v f 2 ξ 2 ,
(70)
corresponding to downward and upward propagating waves. Substituting ξ 3 , q into Eq. (69) delivers the amplitude ratios, Wq, and from Eq. (66), 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 ) = 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 Cremer94 and investigated more deeply by Schoch,95 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. (63) and (64)], introduced by Rokhlin et al.92 The displacements and stresses in the upper fluid are given by superposition of the incident and reflected waves such that
[ u 1 u 3 σ 33 ] u = [ 1 1 W 1 u W 1 u D u D u ] [ e i ξ 3 , 1 u x 3 R e i ξ 3 , 1 u x 3 ] ,
(73)
where R is the reflection coefficient. In the lower half-space, the displacements and stresses are determined by the transmitted wave,
[ u 1 u 3 σ 33 ] l = [ 1 W 1 l D l ] T e i ξ 3 , 1 l ( x 3 + d ) ,
(74)
where T is the transmission coefficient. Specializing Eqs. (73) and (74) to x 3 = 0 , d gives
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 .
(75)
The global stiffness matrix of the composite has already been obtained and should be repeated here for convenience as
[ σ 0 σ d ] = K [ u 0 u d ] .
(76)
Although one can use K to obtain R, it is more convenient to derive it with the compliance matrix, S,
[ u 0 u d ] = S [ σ 0 σ d ] , S = K 1 .
(77)
In the coupled case, under consideration of the vanishing of the shear stresses according to Eq. (64), the expanded form of Eq. (77) is given by
[ u 1 0 u 2 0 u 3 0 u 1 d u 2 d u 3 d ] = [ 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 ] [ σ 33 0 0 0 σ 33 d 0 0 ] .
(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. (75) into Eqs. (79) and (80), under consideration of Eq. (63), 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. (81) and (82), the reflection coefficient is obtained as
R = ( 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, W1 and D, Eq. (83) simplifies to
R = ( Λ 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. (84) reduces to
R = Λ s 31 Λ + s 31 .
(85)
However, similar to that pointed out in Sec. II B, Eqs. (83)–(85) 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. (28) and (29). 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
[ u 1 0 u 2 0 u 3 0 u 1 d / 2 u 2 d / 2 u 3 d / 2 ] = [ 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 ] half [ σ 33 0 0 0 σ 33 d / 2 σ 13 d / 2 σ 23 d / 2 ] .
(86)
Fortuitously, by using the symmetric boundary conditions, Eq. (28), on Eq. (86), 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 = Λ X S Λ + X S , X S = s 31 s 34 s 61 s 64 .
(89)
The antisymmetric boundary conditions, Eq. (29), used on Eq. (86), 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 = Λ X A Λ + X A ,
(93)
with
X A = s 31 s 35 s 51 s 55 ( s 35 s 56 s 55 s 36 ) 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. (83)–(85) as well as those in Eqs. (89) and (93) 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. (95) [respectively, the denominators of Eqs. (83)–(85), (89), (93)] 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 (95) (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. (70), 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. (95) are replaced by the negative amplitude ratios W 2 u and W 2 l.
FIG. 2.

Fluid-loaded composite plate.

FIG. 2.

Fluid-loaded composite plate.

Close modal
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. (78) vanish, and the compliance matrix reduces to 4 × 4 elements. For this case, Eqs. (79)–(85) and (95) can be used by changing the compliance matrix elements according to s 31 s 21 , s 34 s 23 , s 61 s 41 , s 64 s 43. Likewise, u 2 0 , d / 2 and σ 23 0 , d / 2 vanish in Eq. (86), leading to
X S = s 21 s 23 s 41 s 43 , X A = s 21 s 24 s 31 s 34 .
(96)

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 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 S1 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.82 

FIG. 3.

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 S1 Lamb wave.

FIG. 3.

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 S1 Lamb wave.

Close modal

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 x3. 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.39 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. (97) into
[ σ 0 σ d ] = [ K 11 K 12 K 21 K 22 ] [ u 0 u d ] = 0 ,
(98)
the simpler relation,
[ U r U t ] = [ K 11 P 0 + K 12 P d K 21 P 0 + K 22 P d ] 1 [ K 11 P 0 K 21 P 0 ] U in ,
(99)
is obtained with
U in = [ 1 0 0 ] .
(100)
With the known reflected and transmitted wave amplitudes, U r and U t, Eq. (97) 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
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 ,
(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. (101) into Eq. (79), one finds
U l = 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,
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.
(103)
After having solved either Eq. (97) or Eq. (103), 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
[ U m U m + 1 ] = [ P P + H P H P + ] m 1 [ u m u m + 1 ] .
(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. (16), 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 x3 is a value ranging from zero to d. The first three terms in Eq. (106) 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. (102). 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 x3 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.

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 , 90 ). This is possible in a single layer or cross-ply laminate [0/90] n ( s ) , n = 1 , 2 , . In cubic materials, decoupling occurs for guided wave propagation along Φ = 0 , 45 , 90 . 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, SSH, ASH, SScholte, AScholte 
Yes  No  S, A, SScholte, AScholte 
No  Yes  B,a BSH (SSH, ASH),b BScholte 
No  No  B, BScholte 
Symmetric Decoupled Mode families
Yes  Yes  S, A, SSH, ASH, SScholte, AScholte 
Yes  No  S, A, SScholte, AScholte 
No  Yes  B,a BSH (SSH, ASH),b BScholte 
No  No  B, BScholte 
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, SSH and ASH 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 a21 and a22 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
—  det K S = 0 , | det K S | = min 
—  det K A = 0 , | det K A | = min 
—  det K = 0 , | det K | = min 
| det K S | = min 
| det K A | = min 
| det K | = min 
F or F + V  | R S | = min 
F or F + V  | R A | = min 
F or F + V  | R | = min 
SScholte  | R S | = min 
AScholte  | R A | = min 
BScholte  | R | = min 
SScholte  F + V  | R S | = min , ξ 3 fa 
AScholte  F + V  | R A | = min , ξ 3 fa 
BScholte  F + V  | R | = min , ξ 3 fa 
Family Viscoelasticity and fluid-loading Dispersion equations
—  det K S = 0 , | det K S | = min 
—  det K A = 0 , | det K A | = min 
—  det K = 0 , | det K | = min 
| det K S | = min 
| det K A | = min 
| det K | = min 
F or F + V  | R S | = min 
F or F + V  | R A | = min 
F or F + V  | R | = min 
SScholte  | R S | = min 
AScholte  | R A | = min 
BScholte  | R | = min 
SScholte  F + V  | R S | = min , ξ 3 fa 
AScholte  F + V  | R A | = min , ξ 3 fa 
BScholte  F + V  | R | = min , ξ 3 fa 
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
SSH, BSH  No  Multi  a 21 = 0 
ASH  No  Multi  a 22 = 0 
SSH  No  Single  c p = C 66 ρ C 44 ( n / ( 2 f d ) ) 2, n even 
ASH  No  Single  c p = C 66 ρ C 44 ( n / ( 2 f d ) ) 2, n odd 
SSH, BSH  Yes  Multi  | a 21 | = min 
ASH  Yes  Multi  | a 22 | = min 
SSH  Yes  Single  ξ = ρ ω 2 C 44 ( n π / d ) 2 C 66, n even 
ASH  Yes  Single  ξ = ρ ω 2 C 44 ( n π / d ) 2 C 66, n odd 
Family Viscoelastic Layers Dispersion equations
SSH, BSH  No  Multi  a 21 = 0 
ASH  No  Multi  a 22 = 0 
SSH  No  Single  c p = C 66 ρ C 44 ( n / ( 2 f d ) ) 2, n even 
ASH  No  Single  c p = C 66 ρ C 44 ( n / ( 2 f d ) ) 2, n odd 
SSH, BSH  Yes  Multi  | a 21 | = min 
ASH  Yes  Multi  | a 22 | = min 
SSH  Yes  Single  ξ = ρ ω 2 C 44 ( n π / d ) 2 C 66, n even 
ASH  Yes  Single  ξ = ρ ω 2 C 44 ( n π / d ) 2 C 66, n odd 

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

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 A0, 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 A0 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 A0. 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 x1-x2-plane while the x3-direction (out-of-plane) is governed by the low stiffness of the polymer matrix. Another feature is the strong rise of S2 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.

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

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.

Close modal
FIG. 5.

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

FIG. 5.

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

Close modal

As the next example, the [0/90]4s layup is considered. Because wave propagation is in the Φ = 0 -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 A0 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. A0 shows the similar behavior as discussed in the above coupled case. S0 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.

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

FIG. 6.

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

Close modal
FIG. 7.

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

FIG. 7.

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

Close modal

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 S0 should be called S 0 SH (or SH0) and S1 should be S0. 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 SH0, SH1, SH2, 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.”

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 × 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 S0 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.

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

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.

Close modal
FIG. 9.

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

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

Close modal

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 S2 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 S2 is evident through the in-plane displacement component, u1, because it is symmetric about the middle axis of the composite. An interesting finding is that the shear horizontal displacement, u2, seems to vanish inside the composite although all modes should possess coupled polarization. In fact, u2 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.41 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, u3, 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 S2 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 u1, the antisymmetric character of the mode becomes obvious. Because of the aforementioned homogenization effect, u2 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.

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.

Close modal
FIG. 11.

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.

FIG. 11.

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.

Close modal

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.

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

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.

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 Φ .
(A1)
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 ] / Δ ,
(B1)
with
Δ = C 33 C 44 C 55 C 33 C 45 2 .
(B2)
P = [ 1 1 1 V 1 V 3 V 5 W 1 W 3 W 5 ] , P + = [ 1 1 1 V 1 V 3 V 5 W 1 W 3 W 5 ] , D = [ 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 ] , D + = [ 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 ] , H = [ e i ξ 3 , 1 x 3 0 0 0 e i ξ 3 , 3 x 3 0 0 0 e i ξ 3 , 5 x 3 ] , U = [ U 1 , 1 U 1 , 3 U 1 , 5 ] , U + = [ U 1 , 2 U 1 , 4 U 1 , 6 ] ,
(C1)
where the common term, e i ( ξ x 1 ω t ) in H, is suppressed.
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 = [ 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 ] .
(D2)
Imaginary parts of the stiffness matrix in GPa are
C i j im = [ 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 ] .
(D3)
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
).
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
).
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
).
4.
M.
Castaings
and
B.
Hosten
, “
Ultrasonic guided waves for health monitoring of high-pressure composite tanks
,”
NDT&E Int.
41
,
648
655
(
2008
).
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
).
6.
M. J. S.
Lowe
,
D. N.
Alleyne
, and
P.
Cawley
, “
Defect detection in pipes using guided waves
,”
Ultrasonics
36
,
147
154
(
1998
).
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
).
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
).
9.
M.
Luukkala
,
P.
Heikkila
, and
J.
Surakka
, “
Plate wave resonance—A contactless test method
,”
Ultrasonics
9
(
4
),
201
208
(
1971
).
10.
M.
Luukkala
and
P.
Meriläinen
, “
Metal plate testing using airborne ultrasound
,”
Ultrasonics
11
(
5
),
218
221
(
1973
).
11.
H.
Lamb
, “
On waves in an elastic plate
,”
Proc. R. Soc. London, Ser. A
93
,
114
128
(
1917
).
12.
I. A.
Viktorov
,
Rayleigh and Lamb Waves: Physical Theory and Applications
(
Plenum
,
New York
,
1967
).
13.
W. T.
Thomson
, “
Transmission of elastic waves through a stratified solid medium
,”
J. Appl. Phys.
21
,
89
93
(
1950
).
14.
N. A.
Haskell
, “
The dispersion of surface waves on multilayered media
,”
Bull. Seism. Soc. Am.
43
(
1
),
17
34
(
1953
).
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
).
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
.
17.
A. H.
Nayfeh
, “
The propagation of horizontally polarized shear waves in multilayered anisotropic media
,”
J. Acoust. Soc. Am.
86
(
5
),
2007
2012
(
1989
).
18.
A. H.
Nayfeh
, “
The general problem of elastic wave propagation in multilayered anisotropic media
,”
J. Acoust. Soc. Am.
89
(
4
),
1521
1531
(
1991
).
19.
A. H.
Nayfeh
,
Wave Propagation in Layered Anisotropic Media with Applications to Composites
(
North-Holland
,
Amsterdam
,
1995
).
20.
J. W.
Dunkin
, “
Computation of modal solutions in layered, elastic media at high frequencies
,”
Bull. Seism. Soc. Am.
55
(
2
),
335
358
(
1965
).
21.
L.
Knopoff
, “
A matrix method for elastic wave problems
,”
Bull. Seism. Soc. Am.
54
(
1
),
431
438
(
1964
).
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
).
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
).
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
).
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
).
26.
A. K.
Mal
, “
Wave propagation in layered composite laminates under periodic surface loads
,”
Wave Motion
10
(
3
),
257
266
(
1988
).
27.
A. K.
Mal
, “
Guided waves in layered solids with interface zones
,”
Int. J. Eng. Sci.
26
(
8
),
873
881
(
1988
).
28.
E.
Kausel
and
J. M.
Roesset
, “
Stiffness matrices for layered soils
,”
Bull. Seis. Soc. Am.
71
(
6
),
1743
1761
(
1981
).
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
).
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
.
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.
32.
B. L. N.
Kennett
,
Seismic Wave Propagation in Stratified Media
(
Cambridge University Press
,
New York
,
1983
).
33.
B. L. N.
Kennett
and
N. J.
Kerry
, “
Seismic waves in a stratified half space
,”
Geophys. J. Int.
57
(
3
),
557
583
(
1979
).
34.
G. J.
Fryer
and
L. N.
Frazer
, “
Seismic waves in stratified anisotropic media
,”
Geophys. J. Int.
78
(
3
),
691
710
(
1984
).
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
).
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
).
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
).
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
).
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
).
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
).
41.
S. I.
Rokhlin
,
D. E.
Chimenti
, and
P. B.
Nagy
,
Physical Ultrasonics of Composites
(
Oxford University Press
,
Oxford, UK
,
2011
).
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
).
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
).
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
).
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
).
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
).
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
).
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
).
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
).
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
.
51.
L.
Gavric
, “
Computation and propagative waves in free rail using a finite element technique
,”
J. Sound Vib.
185
(
3
),
531
543
(
1995
).
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
).
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
).
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
).
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
).
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
).
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
).
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
).
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
).
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
).
61.
F.
Treysséde
, “
A model reduction method for fast finite element analysis of continuously symmetric waveguides
,”
J. Sound Vib.
508
,
116204
(
2021
).
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
).
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.
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
).
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
).
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
).
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
).
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
).
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
).
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
).
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
).
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
).
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
).
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
).
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
).
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
).
77.
E.
Georgiades
,
M.
Lowe
, and
R.
Craster
, “
Leaky wave characterisation using spectral methods
,”
J. Acoust. Soc. Am.
152
(
3
),
1487
1497
(
2022
).
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
).
80.
P.
Boccini
,
A.
Marzani
, and
E.
Viola
, “
Graphical user interface for guided acoustic waves
,”
J. Comput. Civil. Eng.
25
(
3
),
202
210
(
2011
).
81.
A.
Marzani
and
P
Bocchini
, “
GUIGUW
,” available at http://www.guiguw.com/ (Last viewed June 26,
2023
).
82.
A. M. A.
Huber
, “
Dispersion Calculator (DC)
,” available at https://github.com/ArminHuber/Dispersion-Calculator (Last viewed June 26,
2023
).
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
).
84.
D. R.
Ramasawmy
,
B. T.
Cox
, and
B. E.
Treeby
, “
ElasticMatrix
,” available at https://github.com/dannyramasawmy/ElasticMatrix (Last viewed June 26,
2023
).
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
).
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
).
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
).
88.
D. A.
Kiefer
, “
GEW dispersion script
,” available at https://doi.org/10.5281/zenodo.7010603 (Last viewed June 26,
2023
).
89.
E.
Ducasse
and
M.
Deschamps
, “
Mode computation of immersed multilayer plates by solving an eigenvalue problem
,”
Wave Motion
112
,
102962
(
2022
).
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
.
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
).
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
).
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
).
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).
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#.
96.
M. J. S.
Lowe
, “
Plate waves for the NDT of diffusion bonded titanium
,” Ph.D. thesis,
University of London
, London,
1992
.
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
).