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

## I. INTRODUCTION

Ultrasonic guided waves play a key role in the nondestructive inspection (NDI) or structural health monitoring (SHM) of single and multilayered plates and pipes.^{1–8} The significant complication inherent to guided waves is, however, their dispersive nature. Lamb waves, as a prominent example, represent a case of resonance in the waveguide, where some fundamental modes and a theoretically unlimited number of higher order modes can exist. To select the proper mode and frequency for a given inspection task, a dispersion diagram is essential. Calculating such diagrams is the topic of ongoing research. Historically, isotropic plates were the first specimens considered for NDI,^{9,10} and thanks to the fundamental work of Lamb^{11} 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 Thomson^{13} in 1950 and edited by Haskell^{14} in 1953. In TMM, the transfer matrix relates the displacements and stresses at the top of a layer to those at its bottom. Through the multiplication of the individual layer transfer matrices, assuming the continuity of stresses and displacements across the layer boundaries, a global transfer matrix is obtained, relating the displacements and stresses at the top of the whole structure to those at its bottom. Modal solutions are obtained through imposing the boundary condition of stress-free top and bottom surfaces. Nayfeh and Chimenti^{15,16} and Nayfeh^{17,18} extended TMM to generally anisotropic multilayered media. This work is well-documented in the excellent book by Nayfeh.^{19} However, Dunkin^{20} discovered in 1965 that TMM becomes unstable when the product of the layer thickness and frequency is large.

A different approach is the global matrix method (GMM), originally proposed for isotropic media by Knopoff^{21} 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 Tango^{23} and Schmidt and Jensen^{24,25} showed that this method is numerically stable if implemented properly. Mal^{26,27} extended GMM to anisotropic multilayered media. Kausel and Roesset^{28} conducted a reformulation for isotropic cases and Wang and Rajapakse^{29} conducted a reformulation for an orthotropic layer in the plane of symmetry. Wang and Rokhlin^{30,31} used this approach for the computation of the cell stiffness matrix of a cross-ply composite. Another stable formulation for isotropic media was introduced by Kennett^{32} and Kennett and Kerry,^{33} which was subsequently extended to generally anisotropic media by Fryer and Frazer^{34,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 Wang^{39,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 Sause^{49,50} used SMM in 2018 to calculate dispersion diagrams for thick-walled composites containing up to 400 layers.

The most often applied discretizing method is the semi-analytical finite element (SAFE) method. It was introduced by Gavric^{51} 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 fluid^{54} and elastic waves propagating along weld joints between plates^{55} were modeled by Fan *et al.* Bartoli *et al.* and Marzani *et al.* investigated guided waves in viscoelastic waveguides of rectangular, arbitrary,^{56} and axissymmetric^{57} 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 Craster^{65} for guided wave modeling in elastic media in 2004. SCM is similar to the SAFE approach in that it uses a one-dimensional mesh over the system's cross section, but SCM possesses a higher accuracy and speed of computation. Instead of solving a differential equation directly, SCM uses a spectral approximation, which satisfies the differential equation and boundary conditions. Karpfinger *et al.* applied SCM for the modeling of wave dispersion along multilayered isotropic cylindrical structures^{66} and, later, for axisymmetric wave modes in a porous elastic cylinder^{67} 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 Pavlakovic^{78,79} since the early 1990s, and it uses GMM. The first free software, GUIGUW, appeared in 2011. It was developed by Marzani and Bocchini,^{80,81} based on SAFE method. The MATLAB^{®}-based (The MathWorks, Inc., Natick, MA) open source dispersion calculator (DC), released by Huber^{82} 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, Kiefer^{88} 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.

## II. NUMERICAL MODEL

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

### A. Field equations and geometrical setup

*ρ*is the material's mass density. The stiffness tensor, $ c ijkl \u2032$, contains the elastic material constants. Through symmetry properties and strain energy density considerations, the 81 elements in $ c ijkl \u2032$ 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 \u2032$, Eq. (3) can be written as

Consider a multilayered composite such as that sketched in Fig. 1. Composites are stackings of *m* layers with layer thicknesses, *d _{m}*, consisting of a fiber-epoxy combination. For each layer, a local (crystallographic) coordinate system $ x i , m \u2032 = ( x 1 \u2032 , x 2 \u2032 , x 3 \u2032 ) m$ is assigned, residing on the top of the

*m*th layer such that the layers lie in the $ x 1 , m \u2032$- $ x 2 , m \u2032$-planes. The fibers are oriented along $ x 1 , m \u2032$, and $ x 3 , m \u2032$ is normal to the layer. To describe a system of arbitrary layer orientations in a convenient way, the nonprimed global coordinate system, $ x i = ( x 1 , x 2 , x 3 )$, is introduced. It resides on the top layer of the composite, and it is oriented such that (bulk) wave propagation takes place in the

*x*

_{1}-

*x*

_{3}-plane. With respect to the global coordinate system, the local coordinate systems are yielded by a counterclockwise rotation through an angle, $ \Phi m$, between

*x*

_{1}and $ x 1 , m \u2032$ about the

*x*

_{3}-axis. Therefore, the stiffness matrix, $ C \u2032 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

*C*are given in Appendix A. Notice that

_{ij}*C*has the form of monoclinic material symmetry, in general, having 13 independent elastic constants. For $ \Phi = 0 \u2218 , 90 \u2218$, the elements

_{ij}*C*

_{16},

*C*

_{26},

*C*

_{36}, and

*C*

_{45}vanish, returning orthotropic material symmetry. This will be of importance when the decoupling into pure Lamb and SH waves is studied.

^{91}one quasi-longitudinal wave, L, and two quasi-shear waves, SV and SH. These waves propagate downward ( $ \u2212 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 \u2212$, $ 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

*u*, $ \sigma i 3$ correspond to the originating layer, and $ u \xaf i , \sigma \xaf 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 = \u2212 d$) of the plate vanish, i.e.,

_{i}### B. The SMM

*ξ*( $ = \omega / c p$) along

*x*

_{1}, as required by Snell's law. Here, $ c p$ is the phase velocity in the

*x*

_{1}-direction, and

*ω*is the angular frequency ( $ = 2 \pi 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

*x*

_{3}-component of the wavenumber

*ξ*

_{3}directly, instead of searching for the ratio of the bulk waves' wavenumber components in the

*x*

_{3}and

*x*

_{1}-directions,

*α*. This change is necessary since

*ξ*turns complex once attenuation is introduced. The sought polarization vector is represented by the bulk wave amplitudes,

*U*.

_{i}#### 1. Coupled polarization

*x*, Eq. (8) can be written out as

_{i}*M*. By vanishing the determinant of

_{ij}*M*, the sixth-degree bicubic polynomial equation,

_{ij}*A*

_{1},

*A*

_{2}, and

*A*

_{3}given in Appendix B. Equation (11) has six solutions, $ \xi 3 , q , \u2009 q = 1 , 2 , \u2026 , 6$, namely, three pairs, corresponding to the upward and downward propagating bulk waves,

*m*-layered plate, Eqs. (16) and (17) can be written in matrix form, relating the displacement and stress components at the top $ u m , \u2009 \sigma m$ ( $ x 3 , m = 0$) and bottom $ u m + 1 , \u2009 \sigma m + 1$ ( $ x 3 , m = \u2212 d m$) of the

*m*th layer to the wave amplitude vectors, $ U m \xb1$,

^{39,40}is used. The submatrices, $ P \xb1 , \u2009 D \xb1$,

**H**, and the vectors, $ U \xb1$, are expanded in Appendix C. By eliminating $ U m \xb1$ in Eqs. (18) and (19), the transfer matrix representation is obtained, relating the displacement and stress field components at the top of the

*m*th layer to those at its bottom,

**A**

_{m}is the local transfer matrix of the

*m*th layer. The numerical instability of TMM occurs when contributing bulk waves become evanescent at high frequency-thickness products. Then, $ \xi 3 , q$ becomes imaginary, leading to decaying functions, $ e \u2212 \xi 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.

*m*th layer to the displacements at the top, $ u m$, and bottom, $ u m + 1$, via the local stiffness matrix, $ K m$,

*m*-layered plate, the individual layer stiffness matrices, $ K m$, must be calculated first. Considering the stiffness matrices of two neighboring layers,

**K**, of the whole, symmetric laminate is given by

#### 2. Decoupled polarization

In case of wave propagation along principal axes, polarization in the sagittal plane (*x*_{1}-*x*_{3}) decouples from polarization in the shear horizontal direction *x*_{2}. 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 \u2212$, $ 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.

*x*

_{1}-

*x*

_{3}-plane, Eq. (9) reduces to

*M*. By vanishing its determinant, the fourth-degree bicubic polynomial equation,

_{ij}##### b. Shear horizontal waves.

*x*

_{2}-direction such that the displacement and stress fields are given by

^{19}

**A**

_{m}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.

**A**, was obtained by consecutively multiplying all

**A**

_{m}, starting with the lowermost layer. At least in the case of SH waves, however, one can also start with the top layer, namely,

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

_{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

_{4}is efficiently calculated by the power of

*n*,

*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

*n*= 0 in Eq. (57) shows that the fundamental SH mode is nondispersive at the shear wave speed $ c p = C 66 / \rho $. 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 $ \xi = \omega / c p$ does not hold anymore. Therefore, Eq. (57) is written in the form

*α*, are obtained from the real and imaginary wavenumber parts, $ \xi re$ and $ \xi im$, respectively, according to

### C. Viscoelasticity and fluid-loading

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

^{73}given in Appendix D, serves for the numerical computations. The complex stiffness tensor is given by $ C i j \u2032 = C i j \u2032 re \u2212 i C i j \u2032 im$. Defining $ C i j \u2032$ 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 \u2032$, are given by, respectively,

*α*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

^{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 $ \rho s = 1500 \u2009 \u2009 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 , \u2212 d$),

^{19}wave motion in a nonviscous fluid is described in the global coordinate system by specializing the field equations [Eqs. (2) and (3)] to

*u*

_{2}vanishes. Through combination of Eqs. (65) and (66), the expanded forms,

*q*= 1,2,

*W*, and from Eq. (66), one gets the stress amplitude,

_{q}*D*, as

^{94}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

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

*T*is the transmission coefficient. Specializing Eqs. (73) and (74) to $ x 3 = 0 , \u2212 d$ gives

**K**to obtain

*R*, it is more convenient to derive it with the compliance matrix,

**S**,

*W*

_{1}and

*D*, Eq. (83) simplifies to

### D. Dispersion curve tracing

Dispersion curve tracing in the attenuated case requires a two-dimensional search for modal solutions in the complex wavenumber plane rather than a one-dimensional search on the real axis only. The basic procedure used by DC for finding a modal solution is rather conventional but robust. A particular dispersion equation is evaluated on a raster of complex wavenumbers covering a certain range of real and imaginary parts. For example, Fig. 3 plots the dispersion equation amplitude of Eq. (95) calculated at 500 kHz for wave propagation along $ \Phi = 0 \u2218$ in a 2 mm thick viscoelastic composite layup [0/90/45/–45]_{2s} immersed on both sides in water. Equation (95) is evaluated on a 51 × 51 point raster in the complex wavenumber plane and plotted in the phase velocity-attenuation space. The minimum indicates the modal solution of the S_{1} Lamb wave, which is located at $ c p = 5.98516$ m/ms and $ \alpha = 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}

### E. Mode shape calculation

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

^{39}Accordingly, the displacements at the top and bottom of the system are given by

*m*th 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

*m*th layer is calculated by

*d*is the thickness of the

*m*th layer in this case, and

*x*

_{3}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,

*x*

_{3}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, $ \xi 3 , 2 u$ and $ \xi 3 , 2 l$, must be used in the calculations. All equations given in Sec. II E can be also used in decoupled cases in the reduced form. In the case of SH modes, Eqs. (101)–(103) and (107) can be ignored.

## III. CLASSIFICATION OF SOLUTIONS

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

In general, a classification of solutions is possible according to the symmetry and coupling properties of the considered waveguide and modes. The four possible configurations with the corresponding mode families appearing in a fluid-loaded waveguide are listed in Table I. The prerequisite is that one must know which kinds of mode families occur in a given configuration. The decoupling into pure Lamb, Scholte, and SH waves depends on material symmetries, the layup, and propagation direction of the guided waves. In orthotropic and transversely isotropic media such as unidirectional composite plies, decoupling occurs for wave propagation parallel or normal to the fibers ( $ \Phi = 0 \u2218 , 90 \u2218$). This is possible in a single layer or cross-ply laminate [0/90] $ n ( s ) , \u2009 n = 1 , 2 , \u2026$. In cubic materials, decoupling occurs for guided wave propagation along $ \Phi = 0 \u2218 , 45 \u2218 , 90 \u2218$. 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.

Symmetric . | Decoupled . | Mode families . |
---|---|---|

Yes | Yes | S, A, S^{SH}, A^{SH}, S^{Scholte,} A^{Scholte} |

Yes | No | S, A, S^{Scholte}, A^{Scholte} |

No | Yes | B,a B^{SH} (S^{SH}, A^{SH}),b B^{Scholte} |

No | No | B, B^{Scholte} |

Symmetric . | Decoupled . | Mode families . |
---|---|---|

Yes | Yes | S, A, S^{SH}, A^{SH}, S^{Scholte,} A^{Scholte} |

Yes | No | S, A, S^{Scholte}, A^{Scholte} |

No | Yes | B,a B^{SH} (S^{SH}, A^{SH}),b B^{Scholte} |

No | No | B, B^{Scholte} |

^{a}

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

^{b}

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

Once it is clear which mode families must be expected, the appropriate dispersion equations can be applied to trace them separately. Table II lists those for Lamb and Scholte waves in coupled and decoupled cases, depending on the presence or absence of viscoelasticity and fluid-loading, and Table III lists the dispersion equations of SH waves. Only in nonattenuated cases, one can search for a zero in the dispersion equation amplitude, such as $ det \u2009 K = 0$. Actually, because $ det \u2009 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 \u2009 K$ also can switch between positive and negative infinity (but not *a*_{21} and *a*_{22} in the case of SH waves). Therefore, it must be checked that the absolute value is a (local) minimum, i.e., $ | det \u2009 K | = min$. In attenuated cases, one searches only for minima in the complex dispersion equation amplitude.

Family . | Viscoelasticity and fluid-loading . | Dispersion equations . |
---|---|---|

S | — | det $ K S = 0 , \u2009 |$ det $ K S | = min$ |

A | — | det $ K A = 0 , \u2009 |$ det $ K A | = min$ |

B | — | det $ K = 0 , \u2009 |$ det $ K | = min$ |

S | V | $|$ det $ K S | = min$ |

A | V | $|$ det $ K A | = min$ |

B | V | $|$ det $ K | = min$ |

S | F or F + V | $ | R S | = min$ |

A | F or F + V | $ | R A | = min$ |

B | F or F + V | $ | R | = min$ |

S^{Scholte} | F | $ | R S | = min$ |

A^{Scholte} | F | $ | R A | = min$ |

B^{Scholte} | F | $ | R | = min$ |

S^{Scholte} | F + V | $ | R S | = min , \u2212 \xi 3 f$a |

A^{Scholte} | F + V | $ | R A | = min , \u2212 \xi 3 f$a |

B^{Scholte} | F + V | $ | R | = min , \u2212 \xi 3 f$a |

Family . | Viscoelasticity and fluid-loading . | Dispersion equations . |
---|---|---|

S | — | det $ K S = 0 , \u2009 |$ det $ K S | = min$ |

A | — | det $ K A = 0 , \u2009 |$ det $ K A | = min$ |

B | — | det $ K = 0 , \u2009 |$ det $ K | = min$ |

S | V | $|$ det $ K S | = min$ |

A | V | $|$ det $ K A | = min$ |

B | V | $|$ det $ K | = min$ |

S | F or F + V | $ | R S | = min$ |

A | F or F + V | $ | R A | = min$ |

B | F or F + V | $ | R | = min$ |

S^{Scholte} | F | $ | R S | = min$ |

A^{Scholte} | F | $ | R A | = min$ |

B^{Scholte} | F | $ | R | = min$ |

S^{Scholte} | F + V | $ | R S | = min , \u2212 \xi 3 f$a |

A^{Scholte} | F + V | $ | R A | = min , \u2212 \xi 3 f$a |

B^{Scholte} | F + V | $ | R | = min , \u2212 \xi 3 f$a |

^{a}

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

Family . | Viscoelastic . | Layers . | Dispersion equations . |
---|---|---|---|

S^{SH}, B^{SH} | No | Multi | $ a 21 = 0$ |

A^{SH} | No | Multi | $ a 22 = 0$ |

S^{SH} | No | Single | $ c p = C 66 \rho \u2212 C 44 ( n / ( 2 f d ) ) 2$, n even |

A^{SH} | No | Single | $ c p = C 66 \rho \u2212 C 44 ( n / ( 2 f d ) ) 2$, n odd |

S^{SH}, B^{SH} | Yes | Multi | $ | a 21 | = min$ |

A^{SH} | Yes | Multi | $ | a 22 | = min$ |

S^{SH} | Yes | Single | $ \xi = \rho \omega 2 \u2212 C 44 ( n \pi / d ) 2 C 66$, n even |

A^{SH} | Yes | Single | $ \xi = \rho \omega 2 \u2212 C 44 ( n \pi / d ) 2 C 66$, n odd |

Family . | Viscoelastic . | Layers . | Dispersion equations . |
---|---|---|---|

S^{SH}, B^{SH} | No | Multi | $ a 21 = 0$ |

A^{SH} | No | Multi | $ a 22 = 0$ |

S^{SH} | No | Single | $ c p = C 66 \rho \u2212 C 44 ( n / ( 2 f d ) ) 2$, n even |

A^{SH} | No | Single | $ c p = C 66 \rho \u2212 C 44 ( n / ( 2 f d ) ) 2$, n odd |

S^{SH}, B^{SH} | Yes | Multi | $ | a 21 | = min$ |

A^{SH} | Yes | Multi | $ | a 22 | = min$ |

S^{SH} | Yes | Single | $ \xi = \rho \omega 2 \u2212 C 44 ( n \pi / d ) 2 C 66$, n even |

A^{SH} | Yes | Single | $ \xi = \rho \omega 2 \u2212 C 44 ( n \pi / d ) 2 C 66$, n odd |

## IV. NUMERICAL EXAMPLES

In this section, dispersion diagrams of phase velocity and attenuation are presented. For the calculations, the viscoelastic composite from Ref. 73, whose material constants are given in Appendix D, was used. In DC, the material is available in the list of orthotropic materials as *CarbonEpoxy_Hernando_2015_Viscoelastic*. Fluid-loading is present in the form of water-filled ( $ \rho f = 1000 \u2009 \u2009 kg / m 3 , \u2009 \u2009 v f = 1500 \u2009 \u2009 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 $ \Phi = 0 \u2218$. The first two layups contain 16 layers, summing up to a thickness of 2 mm. These were calculated with DC version 2.4 and DISPERSE version 2.0.20f and serve as validation examples. The third layup contains 400 layers with a total thickness of 50 mm, proving that DC can handle thick-walled composite structures.

### A. Sixteen layers composite

Figure 4 displays the phase velocity dispersion diagram for the [0/90/45/–45]_{2s} layup immersed in water. The bold curves were calculated with DISPERSE, and the thin curves were obtained with DC. As the system is symmetric, one can distinguish symmetric (gray lines) and antisymmetric (black lines) modes. Solid lines represent leaky Lamb modes while the dashed-dotted lines indicate Scholte modes. All modes have coupled polarization. Although DC missed one higher order mode, whereas DISPERSE missed A_{0}, it is obvious that the other curves are in good agreement. The corresponding attenuation dispersion diagram is depicted in Fig. 5, where DC and DISPERSE are in good agreement again. Most remarkable in Fig. 5 is how strongly the attenuation of A_{0} rises at low frequency-thicknesses compared to the other fundamental modes. This is because the water weighs heavily against the strong and predominantly out-of-plane vibration possessed by A_{0}. The out-of-plane displacement is generally more pronounced in composites than in isotropic media as high stiffness carbon fibers are usually oriented in the *x*_{1}-*x*_{2}-plane while the *x*_{3}-direction (out-of-plane) is governed by the low stiffness of the polymer matrix. Another feature is the strong rise of S_{2} 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.

As the next example, the [0/90]_{4s} layup is considered. Because wave propagation is in the $ \Phi = 0 \u2218$-direction, decoupling into sagittal plane motion and shear horizontal motion takes place. Pure leaky Lamb waves and pure SH waves, as well as Scholte waves with only sagittal motion, can propagate in this configuration. The phase velocity dispersion diagram is shown in Fig. 6. The SH modes are represented by dashed lines. Whereas DC did not trace the SH modes up to the top phase velocity of 20 m/ms, DISPERSE missed A_{0} again and left gaps in two dispersion curves. However, the good agreement of both softwares is obvious. The attenuation dispersion is displayed in Fig. 7. A_{0} shows the similar behavior as discussed in the above coupled case. S_{0} 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.

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

### B. Four hundred layers composite

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

Through-thickness mode shapes are presented to shed light on the wave structure inside the waveguide and on the wave field radiated into the water-filled half-spaces. The mode shapes were calculated based on the dispersion curves displayed in Figs. 8 and 9. It should be noted that DC uses the same plotting conventions as DISPERSE, namely, power-normalization of the amplitudes (neglecting the power flow in the half-spaces) and phase normalization with respect to the second sample from the top. In Fig. 10, the displacement component amplitudes of the S_{2} mode at the maximum attenuation located at 4.075 MHz·mm are shown. The white area represents the composite, and the light gray areas indicate the water-filled half-spaces to a depth of 50 mm on both sides. The symmetric character of S_{2} is evident through the in-plane displacement component, *u*_{1}, because it is symmetric about the middle axis of the composite. An interesting finding is that the shear horizontal displacement, *u*_{2}, seems to vanish inside the composite although all modes should possess coupled polarization. In fact, *u*_{2} 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, *u*_{3}, at the surfaces of the composite is related to the transmittance of energy into the half-spaces, which, in turn, leads to the strong damping of S_{2} at 4.075 MHz·mm. As a second example, the $ A 0 Scholte$ mode at 2.5 MHz·mm is considered. It propagates at a phase velocity of 1.31 m/ms and has a low damping of 9.581 Np·mm/m, caused by the composite's viscoelasticity. As depicted in Fig. 11, the most striking difference to leaky waves is that Scholte waves are decaying exponentially inside the half-spaces. From *u*_{1}, the antisymmetric character of the mode becomes obvious. Because of the aforementioned homogenization effect, *u*_{2} 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.

## V. DISCUSSION

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

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

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

## VI. CONCLUSIONS

It was shown that guided wave dispersion diagrams for fluid-loaded viscoelastic composites consisting of up to 400 layers can be obtained with DC, although some higher order modes were missed in some cases. DC is an open source software available at GitHub.^{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.

## ACKNOWLEDGMENTS

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

### APPENDIX A: TRANSFORMED STIFFNESS MATRIX

### APPENDIX B: POLYNOMIAL COEFFICIENTS

### APPENDIX C: SUBMATRICES AND VECTORS

**H**, is suppressed.

### APPENDIX D: MATERIAL PARAMETERS

## REFERENCES

*Rayleigh and Lamb Waves: Physical Theory and Applications*

*Review of Progress in Quantitative Nondestructive Evaluation,*

*Wave Propagation in Layered Anisotropic Media with Applications to Composites*

*Review of Progress in Quantitative Nondestructive Evaluation*

*Review of Progress in Quantitative Nondestructive Evaluation*

*Seismic Wave Propagation in Stratified Media*

*n*-layered anisotropic composite laminate

*Physical Ultrasonics of Composites*

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

*Review of Progress in Quantitative Nondestructive Evaluation*

*Review of Progress in Quantitative Nondestructive Evaluation*