## Abstract

Guided waves propagating in lossy media are encountered in many problems across different areas of physics such as electromagnetism, elasticity and solid-state physics. They also constitute essential tools in several branches of engineering, aerospace and aircraft engineering, and structural health monitoring for instance. Waveguides also play a central role in many non-destructive evaluation applications. It is of paramount importance to accurately represent the material of the waveguide to obtain reliable and robust information about the guided waves that might be excited in the structure. A reasonable approximation to real solids is the perfectly elastic approach where the frictional losses within the solid are ignored. However, a more realistic approach is to represent the solid as a viscoelastic medium with attenuation for which the dispersion curves of the modes are, in general, different from their elastic counterparts. Existing methods are capable of calculating dispersion curves for attenuated modes but they can be troublesome to find and the solutions are not as reliable as in the perfectly elastic case. In this paper, in order to achieve robust and accurate results for viscoelasticity a spectral collocation method is developed to compute the dispersion curves in generally anisotropic viscoelastic media in flat and cylindrical geometry. Two of the most popular models to account for material damping, Kelvin–Voigt and Hysteretic, are used in various cases of interest. These include orthorhombic and triclinic materials in single- or multi-layered arrays. Also, and due to its importance in industry, a section is devoted to pipes filled with viscous fluids. The results are validated by comparison with those from semi-analytical finite-element simulations.

## 1. Introduction

Robust and reliable computation of dispersion curves is key for the successful development of non-destructive evaluation (NDE) techniques using guided waves. Dispersion curves provide the information required to correctly select propagating modes for the study of the NDE of a particular type of defect, or some property of the materials, in a structure. The usage of dispersion curves for NDE is well established as the copious literature and studies on the subject show: analytical solutions using potentials or the partial wave decomposition for the isotropic plate and cylinder were found by Mindlin [1], Pao [2,3], Gazis [4] and Zemanek [5]. Studies and solutions for anisotropic media in flat, and in cylindrical, geometry are in Solie & Auld [6], Nayfeh & Chimenti [7] and, more recently, Li & Thompson [8] or Towfighi *et al.* [9]. Also, attention has been given to multi-layer systems where fluid layers could also be present, such as fluid-filled pipes or plates surrounded by infinite fluid (or solid). Solutions based on the transfer or global matrix method are in [10,11] or [12] and references therein. For more recent accounts covering these topics and related advances in the field one can look in [13–16]. Finally, some classical texts on the subject of guided waves are [17–19] or [20]. These texts treat a variety of the problems above and their different applications in engineering and NDE.

Most of the calculations cited above assume materials to be perfectly elastic. However, a more general and realistic approach to guided wave problems is to allow the materials to possess some kind of material damping, such as viscoelasticity. By doing so, we allow for dissipation and loss of energy by various microscopic mechanisms within the material. As a consequence of this, waves propagating within the structure generally present attenuation and decay of their amplitude as they propagate. The development of efficient and reliable NDE techniques to inspect these kinds of media, commonly encountered in industry, is based on a thorough understanding of their physical behaviour, as well as on robust and accurate tools to model and plot the dispersion curves of the modes they support. In a recent paper [21], a spectral collocation method (SCM) [22] was generalized as an alternative to the classical partial wave root finding (PWRF) routines for solving elastic (lossless)-guided wave problems. This method presented a number of advantages over the PWRF which made it more robust and reliable. Moreover, due to its generality, the SCM solved cases that were very difficult or not possible to solve with the PWRF. A range of illustrative examples, ranging from single- to multi-layer systems in flat and cylindrical geometry was presented and the different ways of validating the SCM were discussed.

Our focus is on lossy media in the NDE context, but the methodology we develop is relevant in broader settings and it is worth highlighting that attenuated guided waves appear not only in the context of elasticity and NDE but also in other branches of physics and engineering. For this reason, the SCM deployed in this paper is perfectly transferable across a wide range of disciplines. As pointed out in [22], spectral methods were introduced in the 1970s in the field of fluid dynamics by Kreiss & Oliger [23], Orszag [24] and Fornberg [25] and have remained a standard computational tool in the field ever since. Recently, the SCM has been successfully used in the fields of seismology and geophysics to study wave propagation, see for instance [26–28]. Guided wave problems in other contexts such as electromagnetic waveguides can also be tackled by means of the SCM, see for instance [29], where metal–insulator–metal electromagnetic lossless and lossy waveguides are studied. A last example of the wide applicability of the methodology presented in this paper is provided by [30] where the eigenstates of the Schrödinger wave function in quantum rings are studied with the aid of a SCM. In general, the SCM can be applied to solve problems which are posed in the form of an eigenvalue problem.

In addition, laminates consisting of multiple layers of viscoelastic anisotropic materials, such as carbon fibre composites, as studied in this paper, are being widely introduced in many engineering structures particularly in aerospace applications. These materials exhibit damping that varies hugely across the different modes, frequencies and directions of propagation. Also, the popular subject of structural health monitoring (SHM) has strong interest in using guided waves to monitor large areas of composite plate structure by means of permanently attached transducers. The ability to calculate the dispersion curves with damping is essential for that purpose. Therefore, it is intended that the results and methodology presented in this paper will be helpful to applications across a breadth of disciplines.

Regarding the possible mechanisms for attenuation, it must be remembered that decay of elastic-guided waves can be caused by material damping, fluid viscosity or by energy leakage into an infinite medium surrounding the waveguide. The reader must bear in mind that in this paper structures are surrounded by a vacuum so the only mechanisms causing attenuation are material damping and fluid viscosity when fluids are present. Note that attenuation due to leakage of energy happens regardless of the nature of the material, that is, it affects perfectly elastic as well as viscoelastic materials. In the latter case, both mechanisms contribute to the wave's attenuation, which is higher than in the former. In this paper only free viscoelastic structures will be studied, hence no attenuation due to leakage of energy is considered. Being intrinsic properties of the medium, damping and viscosity will cause attenuation of any perturbation within the medium. For propagating modes, attenuation will describe the decay of the amplitude of the fields as the wave travels through the structure. For non-propagating (local) elastic waves, attenuation describes the spatial decay of the fields in the waveguide; the former is the object of study of this paper. However, solutions for attenuated local non-propagating elastic waves can also be found by the same procedure described in the following sections.

One of the methods usually used to approach guided wave problems in viscoelastic materials is the PWRF approach mentioned above, which, *for a given value of the frequency*, searches for the values of the wavenumber *k* satisfying the dispersion relation. It is important to note how this is in clear contrast with perfectly elastic cases, in which the search for roots is performed in *k* is generally complex, and the search must be carried out in

The PWRF has been successfully used by various authors to model guided waves in viscoelastic materials. Nagy & Nayfeh [31] studied the effect of viscosity of the loading fluid on the longitudinal waves propagating in a multi-layer system of cylinders. The model used for the fluid was a hypothetical isotropic solid and the attenuation was described using the Kelvin–Voigt model described in the next section. More recently, members of the NDE group at Imperial College have successfully investigated the propagation of guided waves in viscoelastic composites [32] using DISPERSE [33,34], a partial wave-based software package developed by them. Similar studies in cylindrical geometries can be found in [35,36]. In [32], the authors preferred to use the Hysteretic model although PWRF, and DISPERSE in particular, can accommodate both Hysteretic and Kelvin–Voigt approaches. The solution of the equations for a given value of the frequency, as the PWRF and SCM approaches do, is not limited to those two models but allows for a very wide variety of damping models. More information about them is in [33,34]. Finally, Bernard *et al.* [37] studied how energy velocity was affected by absorbing layers using PWRF.

Another alternative recently and successfully used to solve viscoelastic-guided wave problems [38,39] is the semi-analytical finite-element (SAFE) method. In this approach, it is also possible to implement the Kelvin–Voigt and Hysteretic models as described in [38]. The SAFE methodology is a valuable tool for the study of guided waves, particularly because of its powerful treatment of waveguides of arbitrary cross-section; however, it does present some difficulties such as the overestimation of frequencies due to the higher stiffness of discretized structures which can be dealt with by simply increasing the number of elements as described in [39] or the need to filter spurious modes produced in the simulation used in this paper due to the periodic boundary conditions (BCs) of the SAFE scheme. Although some automation has been used, more usually the filtering has to be done manually. More details about the SAFE scheme used in this paper can be found in §3 and references therein. Besides, the implementation of a SAFE model requires sound knowledge of the finite-element (FE) procedure and use of specialist FE codes.

In this paper, propagating modes in viscoelastic media are studied. The main contribution is the successful extension of the previously implemented SCM approach [22,21] to finding guided waves in generally anisotropic viscoelastic media in flat and cylindrical geometry by means of a companion matrix technique as described in the following paragraphs and in §2.

As mentioned above, the greatest difficulty that arises in the modelling of viscoelastic materials is the search for complex roots. Generally, if the PWRF algorithm is capable of performing searches in both variables, it looks for pairs (*ω*,*k*) that satisfy the dispersion relation, that is it searches in *ω*,*k*) are both real for non-attenuating propagating modes. This extra dimension of the solution space makes the search for roots in the viscoelastic case even harder than in the perfectly elastic case and the probability of missing a root is, therefore, increased. The roots of attenuated modes are also very important for NDE practitioners, engineers or researchers in the field, but they are relatively difficult to find, especially for demanding cases such as those involving multiple layers, anisotropy or cylindrical geometry. Therefore, robust and reliable new methods to compute all roots of the dispersion relation are extremely valuable. In the method presented here, a real frequency *ω* is fixed and the algorithm finds the eigenvalues that are precisely the complex *k*. Again, it must be emphasized that, this approach differs fundamentally from the one used for perfectly elastic materials in a previous paper by Hernando *et al.* [21] in which *k* is fixed and solutions for real *ω* are sought.

The SCM overcomes these difficulties as it solves the equations algebraically by finding the eigenvalues of an analogous matrix problem [22,21]. The SCM is equally applicable to problems with complex or real eigenvalues at no extra cost or effort on the part of the modeller and it has the noteworthy advantage of not missing any modes. As the SCM computes all the eigenvalues rather than looking for zeroes of a function, as the PWRF does, one can be sure to obtain a complete solution. A more detailed description of the SCM's features shall not be pursued here, the interested reader will find an exhaustive discussion of the SCM and its applications to elastic-guided wave problems in [22] or [21]. These references provide detailed explanations about how the number of grid points affects the accuracy of the results. Additionally, the books by Gottlieb & Orszag [40], Boyd [41], Trefethen [42] and Fornberg [43] are established references in the field of Spectral Methods and contain rigorous derivations of several features of the SCM.

The dispersion curves of a viscoelastic material are often very different from those of its perfectly elastic counterparts depending on the value of the damping. Viscoelastic media, no longer makes sense to speak about cut-off frequencies as the solutions are complex. When one sets a limit to the attenuation of the mode for it to be considered a propagating one, the dispersion curves might look incomplete or two different modes appear to merge or cross at a point where no crossing was seen in the perfectly elastic counterpart. These difficulties, added to the ones explained above, sometimes render the PWRF approach misleading or not very robust when modelling viscoelastic media. This highlights the need for a robust and reliable algorithm to find dispersion curves in viscoelastic materials. In addition, a robust and reliable algorithm based on the SCM serves as a solid foundation for the development of more complicated models for the leaky and trapped modes that occur in embedded structures.

The paper is structured as follows. In the second section, the viscoelastic models presented in the paper, as well as the necessary modifications for the SCM scheme to handle them are described, and references are given for descriptions of more basic SCM schemes. The third section summarizes the SAFE models used to validate our results. The fourth section is devoted to systems in flat geometry, single- and multi-layer examples of the most relevant or general cases are shown and compared with the results given by the SAFE simulation. Section five deals with single- and multi-layer systems in cylindrical geometry, the results and how they were validated are presented. In many cases, the validation of the cylindrical cases has followed the same steps as their perfectly elastic counterparts, therefore, when appropriate, references are given to the relevant literature. Owing to its paramount importance in NDE, section six focuses on viscoelastic or perfectly elastic pipes filled with perfect or viscous fluids. The challenges posed to the SCM by this family of problems are described as well as the theoretical framework chosen to model the fluid layers. The paper ends with a summary and a discussion of the results and possible lines for future work.

## 2. Spectral collocation scheme and the companion matrix

We begin with a description of the viscoelastic models used here. The equations of time-harmonic motion for a linear elastic anisotropic homogeneous medium are
*c*_{KL} is the medium's stiffness matrix in reduced index notation, [18], *u*_{j} are the components of the displacement vector field
*X*,*Y*,*Z*} are the crystal axes (blue) which can rotate about the fixed spatial axes {*x*,*y*,*z*} (black) according to the choice of orientation of the material within the waveguide. In the flat case, the {*z*}-axis is the phase direction of the propagating waves (normal to the plane of the wavefront), and in the cylindrical case, it is the axis of the cylinder. For the cylindrical case, {*x*} and {*y*} should be replaced by {*θ*} and {*r*}, respectively. This configuration is assumed throughout the paper unless otherwise stated. Nevertheless, the reader will be reminded of it where appropriate.

For structures in a vacuum, traction-free BCs must be taken into account that require the vanishing of the following three components of the *stress tensor field* defined below, *T*_{ij}. Taking the faces of the plate to be located at *y*=±*h*/2 the BCs are given by
*stress tensor field* in terms of the *strain tensor field* for a *perfectly elastic* material reads
*c*_{ijkl} is the fourth-rank stiffness tensor, related to *c*_{KL} as described in [18]. When material damping is taken into account, the entries of the stiffness tensor of the material are no longer real but also have a complex part. The *strain tensor field*, *S*_{ij}, in terms of the *displacement vector field*, *u*_{j} is

The Kelvin–Voigt (KV) model is well established for describing viscoelastic media and it is briefly described below. In this model, the imaginary part of the stiffness matrix depends on the frequency. Thus, following Auld [18], one has the extended version of (2.4)
*S*_{kl}∝*e*^{−iωt}, etc. assumed henceforth then one is left with
*κ*_{pjkl} has units of sPa, not Pa as *c*_{pjkl}; to have the same units in both entries and avoid possible confusion the following tensor is defined:
*c*_{pjkl} and *η*_{pjkl} have the same units and the prefactor

The second model considered is the Hysteretic (H) model that is a simplification of the previous one. In this model, one takes the viscoelastic stiffness matrix to be frequency independent:

The SCM scheme for viscoelastic materials is similar to that used for their perfectly elastic counterparts, see for instance [22,21]. The novelty in the present case, with respect to the study for elastic materials presented in [21], is that the SCM will be deployed to search for complex values of the wavenumber *k* rather than real ones as in the perfectly elastic case and this has consequences in how the SCM is implemented as described in the next paragraphs. In fact, as has already been emphasized, in the perfectly elastic case, the real value of *k* is fixed and one solves the eigenvalue problem for the real values of frequency. This was more convenient in perfectly elastic cases because the frequency only enters the equations in one term of the PDE so they can be directly recast into a general eigenvalue problem.

In the case of viscoelastic materials, it is better to fix a real frequency (one-dimensional space) and solve for the complex values of *k* (two-dimensional space) that satisfy the dispersion equation. The issue now is that *k* does not enter linearly into the equations of motion. This appears to make it impossible to recast them into a general eigenvalue problem but an algebraic manipulation known as the *Linear Companion Matrix method* allows for this to be rearranged into the more convenient form of a general eigenvalue problem. This is achieved at the cost of doubling the dimension of the matrices involved, see Bridges & Morris [44] and references therein for further mathematical details about this method and their impact on the eigenvalues computation. More recently, other authors have successfully used this approach in guided wave problems, see Pagneux & Maurel [45] or Postnova [46] for instance. This success has motivated the choice of this linearization scheme and as the results have been satisfactory an extension to other schemes has not been pursued. A brief outline of the manipulations is given in the following paragraph. It must be noted that this manipulation has only been performed in the SCM scheme.

For the viscoelastic case, the equations of motion (2.1) and BCs (2.3) are recast into a general eigenvalue problem by the standard SCM procedure: we discretize and substitute the derivatives by differentiation matrices (DMs) computed owing to the Matlab suite provided by [47]. For a single layer in a vacuum, as we have a bounded interval, the appropriate choice is to use Chebyshev DMs, based on a non-uniform Chebyshev grid of *N* points, these are *N*×*N* matrices; the generation of DMs is covered in [42,47]. The *m*th derivative with respect to *y* is approximated by the corresponding *m*th order Chebyshev *N*×*N* DMs:
**U** is the *vector of vectors*: **U**=[**U**_{1},**U**_{2},**U**_{3}]^{T}, these vectors **U**_{i} are the components of the displacement vector field. The matrix *ρ*. The BCs in equation (2.3) are taken into account by appropriately substituting the corresponding rows of (2.13) by those of (2.14) in the following fashion: the 1,*N*,(*N*+1),2*N*,(2*N*+1) and 3*N* rows of the

The *k*^{0},*k*^{1} and *k*^{2}. Let us rearrange the terms in the matrix equation (2.13) and decompose the *k* dependence of the different terms becomes more apparent. One can write
*ω*, this does not have the structure of a general eigenvalue problem in *k*. To achieve this the following definitions are required. Let *companion displacement vector field*:
**I** is the identity matrix and **Z** is a matrix of zeroes. With the above definitions, equation (2.15) is more conveniently expressed as
*k*. Regarding the computational resources, we use the routine `eig` of Matlab (v. R2012b) on an HP desktop computer. Detailed numerical comparisons and studies of the SCM have been already carried out in [22,21], the reader is referred to these papers for more details on convergence and accuracy.

To better process the results obtained and retain only the propagating modes which are the object of study in the present paper, a ratio between the real and imaginary parts of *k*=*α*+i*β* is defined
*k* have been found, *R* is used to select only those with a low value of attenuation within the range of wavenumbers under study, these are the propagating modes. Here, the propagation direction of the harmonic perturbation is taken in the positive direction of the {*z*}-axis, the displacement vector field was taken to be proportional to **U**∼*e*^{ikz} in equation (2.2), therefore, *β*≥0 to ensure the decay of **U** with distance.

## 3. The semi-analytical finite-element models for validation of results

The SAFE method is popular for studying properties of guided waves along waveguides with arbitrary cross-section, such as railway lines [48], beams [49,50], welded [51] or stiffened plates [52], etc. It uses FEs to represent the cross-section of the waveguide, plus a harmonic description along the propagation direction, thus limiting the FE model to two dimensions. SAFE analysis can be deployed using specific programming [38], or advanced use of a flexible commercial code such as COMSOL (2014). In this paper, a method described by Predoi *et al.* [53] is implemented to validate results from the SCM method.

The SAFE method assumes that there is no geometric variation of the cross-section along the axis of the waveguide, so the behaviour in the wave propagation direction can be written in analytical form. Thus, the displacement vector in the waveguide is written as
*k* is the wavenumber, *ω*=2*πf* is the angular frequency, *f* is the frequency, *t* is the time variable and the subscript *j*=1,2,3. The function *U*_{j} represents the behaviour in the cross-section of the waveguide, for which the geometry is irregular, such that it is incorporated in the model by a two-dimensional FE discretization. For general anisotropic media, the equation of dynamic equilibrium is written in the following form of an eigenvalue problem:
*j*=1,2,3 and *q*,*l*=1,2. The coefficients *c*_{ijkl} are the stiffness moduli and *δ*_{ij} is the Kronecker symbol. The equation is reconstructed and solved in the format of a standard eigenvalue problem in the commercial FE code COMSOL (2014), and the full details are provided in [53].

The geometry is meshed by square elements of second order, with side length of 0.05 mm. Periodic BCs [53] were imposed at the lateral boundaries of the domain in the plate models. The number of degree of freedom in our models is less than 11 000, and the typical calculation time for each SAFE model on a standard PC (Intel Core i7, 8 GB memory) was less than 15 s. The convergence has been checked in the models used in our paper, and the results are satisfactory. In fact, this implementation of the SAFE model has already been validated in multi-layer structures [54] as well as solid–fluid structures [55] in which there is a strong impedance contrast. For chosen values of angular frequency *ω*, eigenvalues of complex wavenumber *k* are found, in which the real part describes the harmonic wave propagation while its imaginary part presents the attenuation. Each solution at a chosen frequency reveals the wavenumbers of all possible modes at that frequency; then the full dispersion curve spectrum is constructed by repeating the eigenvalue solutions over the desired range of frequencies. In the following sections, SAFE models will be applied to both flat and cylindrical structures, and the results are used for validation of those obtained by the SCM.

## 4. Flat geometry

In this section, some illustrative examples in flat geometry are presented. Orthorhombic materials are treated firstly as a preparatory example because they have already been studied in the literature (see references). Then, a few novel cases of the most general choice of anisotropic material, triclinic, are presented and the section finishes with a multi-layer example.

Orthorhombic materials are commonly encountered in industry and have already been studied in two references given in the introduction, namely [38,32]. In the SCM context, a code for an orthorhombic medium can also be used for all those materials whose stiffness matrix has a similar block structure, such as hexagonal or isotropic; we begin by presenting an example of a viscoelastic orthorhombic plate in vacuum. The thickness of the plate is 1 mm, the propagation takes place along the {*Z*} crystal axis and the {*Y* } crystal axis is perpendicular to the plane of the plate. More details about the physical properties of this plate are given in the appendix at the end. This example has been done using the Kelvin–Voigt model described in the previous section and the parameter *R*=0.5 for both Lamb and SH modes.

The aforementioned first case is presented in figure 2*a*,*b* that show the phase velocity curves and the attenuation of the Lamb and SH modes, respectively. The comparison between the results given by the SCM (blue circles) and those given by the corresponding SAFE simulation (red asterisks) is excellent. In figure 3, a detail of figure 2*a* is shown, and some spurious solutions found by the SAFE simulation are clearly seen, which are pleasingly not given by the SCM approach.

The next case presented is a 1 mm thick viscoelastic triclinic plate using the Kelvin–Voigt model. The propagation takes place along the {*Z*} crystal axis and the {*Y* } crystal axis is perpendicular to the plane of the plate. Figure 4*a* features the dispersion curves and figure 4*b*, the attenuation for this case. The solution given by the SCM (blue circles) is compared with that given by SAFE (red circles) and the agreement between both solutions is very good. This case is of particular importance because triclinic materials are the most general type of anisotropic material; with 21 independent constants and only one centre of symmetry which imposes no restriction over the stiffness constants this case poses a great challenge for PWRF routines because of its complicated dispersion relation. In [56], the reader can find a brief discussion about these materials. Further physical properties for this case, as well as the parameter *R* used in the SCM code, are given in the appendix.

The last case of single plate systems uses the Hysteretic model to compute the attenuation and dispersion curves showing that the SCM works equally well with this simpler damping model. A 50 mm thick viscoelastic triclinic plate is studied, the physical elastic properties of this plate are the same as the one of the previous examples and are given in the appendix. The propagation takes place along the {*Z*} crystal axis and the {*Y* } crystal axis is perpendicular to the plane of the plate. This last example is plotted in figure 5*a*,*b*, displaying the dispersion curves and attenuation, respectively. The SCM solution (blue circles) is compared with excellent agreement to the solution given by the SAFE simulation (red asterisks).

Three illustrative single plate examples have been presented to illustrate how the SCM can handle equally well two of the most used models to account for attenuation, the Kelvin–Voigt model and the Hysteretic model; the agreement with the SAFE simulation was excellent for all the examples.

This section is closed with a multi-layer example shown in figure 6*a*,*b*. The system is composed of three plates: the top plate is 8 mm thick and it is made of a viscoelastic triclinic material, the middle plate is 5 mm thick and it is made of a viscoelastic orthorhombic material and the bottom plate is a 3 mm thick perfectly elastic triclinic layer. Note that the total thickness has been used for the *x*-axis in both figures. The materials of all three layers are different and the physical parameters of the system are retrieved from the appendix. In all the layers, the propagation takes place along the {*Z*} crystal axis and the {*Y* } crystal axis is perpendicular to the plane of the plate, crystal and spatial axes are aligned and the Hysteretic model was chosen to account for the damping in the materials. Figure 6*a* displays the dispersion curves and figure 6*b* the attenuation. The solutions given by the SCM (blue circles) are once more validated with those given by the SAFE simulation (red asterisks) showing good agreement.

It should be noted that the first case of the section studying orthorhombic materials involved propagation along a principal axis and nothing explicit has been said about propagation at an arbitrary angle with respect to the principal axis which can also be seen as a case of monoclinic material with propagation within its plane of symmetry, see for instance [7]. This case as well as the even more general one of a monoclinic material with propagation direction at an arbitrary angle with respect to the plane of symmetry, which after rotation yields a triclinic stiffness matrix, can be solved by employing the codes for triclinic materials which have been studied and validated here. Therefore, for the sake of brevity, these ‘particular’ cases have not been studied explicitly and attention has been given only to the general triclinic case containing them as limiting cases.

## 5. Cylindrical geometry

For the cylindrical geometry, with guided waves propagating along the {*Z*} crystal axis, the solutions are classified according to their *circumferential harmonic order*: *n*. In general, the solutions have the following form:
*n*=0 yield two decoupled independent modes, we refer the reader to [21] for the procedure to find out whether there is decoupling of modes or not in the different materials. These modes are known in the literature as Torsional and Longitudinal modes, see [17] or [18] for instance. The former are the cylindrical analogues of the plate SH modes, the latter are the analogues of the plate Lamb waves. In the literature (see previous references), the solutions with *n*≠0 are called Flexural modes. The different families are labelled by the *circumferential harmonic order*, *n*. Note that in materials such as monoclinic or triclinic the splitting of the *n*=0 family into Torsional and Longitudinal modes does not take place because of their lower degree of symmetry so we simply call them Flexural modes of order *n*=0. Unless otherwise stated, in the following examples the propagation is along the {*Z*} crystal axis. The crystal axes are also aligned with those of the cylinder as described in §2 and shown in figure 1.

The SCM codes for viscoelastic materials in cylindrical geometries are validated in exactly the same way as their perfectly elastic analogues: by taking the *thin plate limit* making the internal radius of the cylinder very large compared with its thickness, and comparing the results with those of a plate of equal properties. This procedure was explained in detail and illustrated with examples for all types of cylindrical solutions, including circumferential propagation, in [21]. For the cases presented in this paper, the same testing procedure was followed and the reference cases were those presented in the previous section which comprise all types of anisotropic materials and have been successfully validated with the SAFE method.

It is also interesting to see the agreement between the results given by the SCM and SAFE when the *thin plate limit* is not taken. In figure 7*a*,*b*, a comparison between the results given by the SCM and SAFE is presented. It is a viscoelastic hexagonal cylinder with inner radius *r*_{i}=50 mm, thickness *h*=15 mm and axial propagation. The Kelvin–Voigt model has been chosen for modelling the material damping of this case. The figure only displays a short range of frequencies because the number of modes present for higher frequencies increases enormously and this renders the comparison of both solutions unclear. The modes computed by the SCM have been plotted with squares of different colours according to the family to which they belong: *n*=0 in black squares, *n*=1 in blue squares, *n*=2 in red squares and *n*=3 in green squares. The solution computed by SAFE is plotted in asterisks with the same colour scheme as above. Note that the SCM gives the results for each family *separately*, once we have fixed the value of *n* the code computes only the modes belonging to that family.

Because of the central role played by fibre reinforced composites in the industry, this section is closed with an example of a hexagonal (transversely isotropic) viscoelastic material, such as is found in fibre-reinforced composites. An interesting and challenging problem for more conventional approaches such as PWRF routines, though not for the SCM, is the case of a composite where the fibres along the {*Z*}-axis of the crystal were chosen to form an angle of 35° with the axis of the cylinder, this angle does not correspond to any symmetry of the hexagonal system. Thus, the fibres follow a helical path around the axis of the cylinder. The cylinder has inner radius *r*_{i}=20 mm, thickness *h*=7 mm and the propagation is along the axis of the cylinder. The dispersion curves and attenuation for the family *n*=0 using the SCM are shown in figure 8*a*,*b,* respectively. This example was done with the Kelvin–Voigt model for damping and the physical parameters for this example are given in the appendix at the end.

## 6. Pipes containing fluids

Owing to the importance of pipes filled with fluids in various sectors of industry, NDE engineers need to have robust tools to perform studies on them and obtain the necessary information about modes that propagate and of their peculiarities. The examples shown in this section deal with viscous fluids as well as with ideal fluids in viscoelastic or elastic pipes, the different models suitable to describe the fluid are also discussed. They are validated by comparison to relevant examples from the literature and to results obtained by using the PWRF approach. Because of the way the implementation of the SCM is done to handle cases with material damping, the pipes can be chosen to be elastic or viscoelastic; therefore, throughout this section, the reader should bear in mind that whenever we encounter an example of elastic pipe, it is equally possible to study its viscoelastic counterpart using exactly the same code and vice versa. To better study the effects of fluid viscosity and material damping, these two mechanisms are isolated from each other in the study cases. The section closes with a more general and illustrative example of an orthorhombic pipe with viscous fluid inside.

Before presenting the results, a brief discussion about the different approaches to fluid viscosity and their effects on the results is given. In the context of guided wave NDE, the most established model for viscous fluids regards them as hypothetical solids. For constructing this hypothetical solid two main alternatives have been proposed. The first alternative consists of simply regarding the fluid as a *true solid* with the appropriate stiffness constants: the real part of shear elastic modulus *μ* must be taken to be zero, because in the absence of viscosity one must recover the case of ideal fluids that do not support shear waves; also the imaginary part of the longitudinal elastic modulus λ is set to zero because viscosity is assumed to arise solely from shear motion. The stiffness matrix entries are exactly the same as those of an isotropic solid. In the second alternative, λ and *μ* are taken as in the previous one, but the contribution of the viscosity (non-zero imaginary part of *μ*) is split among the stiffness matrix entries differently than usually done in solid mechanics. Viscosity also enters in those entries with only longitudinal contribution, this model is often referred to as the Stokes model. A more detailed discussion is in [57] and the equations for these models can be retrieved in [33]. A final comment about the two models described above: as explained in [57] the numerical results given by each of the models will differ in general. As there is no solid conceptual reason for choosing *a priori* one or the other model, which of them is more suitable in a particular case is determined by the results from experiments or other measurements obtained in similar cases.

The first cases studied in this section are pipes with viscous fluid inside. Without loss of generality, the pipes will be elastic so that the attenuation is solely caused by the viscosity of the fluid. The dynamical viscosity *η* for most fluids of interest ranges from orders of 10^{−1} (very low viscosity fluids) to 10^{+1} (high viscous fluids), and in all cases the coefficients representing this viscosity are several orders of magnitude smaller than those of the elastic solids. When the problem, pipe and fluid, is described by displacement fields, spurious modes are present because of the numerical ill conditioning caused by the big differences among terms in the equations produced in turn by the difference in orders of magnitude of the constants entering the problem. One avoids this by turning to a mixed description of the problem using potentials in the fluid layer and displacements in the solid layer. This enables a homogenization of the terms in the equations of motion of the fluid without any loss of generality regarding the solid layer: recall that a potential description of anisotropic solids is not possible so the displacement description is necessary lest generality is lost. A similar issue was also found in elastic pipes with ideal fluids, see [21].

As the potential description can also be used for isotropic solid materials, the first example computes the longitudinal modes of a hypothetical viscoelastic brass solid rod of radius *r*_{i}=23 mm coated with a layer of elastic steel of thickness *h*=5 mm; the hysteretic model has been chosen to account for material damping within the solid rod and the propagation takes place along the axis of the pipe. This case serves as a testbed for the SCM code that is validated by comparison with results given by the PWRF approach. This case is chosen for its numerical simplicity, as all the elastic constants involved are of similar orders of magnitude and it is not quite as challenging as that of viscous fluids.

The dispersion curves and the attenuation of the modes is shown in figure 9*a*,*b*, respectively. The results given by the SCM are shown in red circles and the solution given by the PWRF approach is shown in black solid lines. The agreement between both approaches is excellent.

The next two examples show the agreement of the results given by the SCM code just validated when used for modelling viscous fluids with those given in [58]. In the following two cases, the Stokes model is assumed throughout to model the viscous fluid. The first example in figure 10 shows a more detailed comparison of the attenuation curves of longitudinal modes propagating along the axis of the pipe obtained with the SCM solution to the PWRF solution given in fig. 11 of [58]. The pipe is made of steel with inner radius *r*_{i}=4.5 mm and thickness *h*=0.5 mm. The fluid is glycerol, the solutions given by the SCM for the viscosities *η*=0.8 Pas, *η*=0.94 Pas and *η*=1.1 Pas are shown in black hollow squares, red diamonds and blue circles, respectively. The corresponding PWRF solutions in solid black lines are indicated by the arrows in the background figure taken from the reference [58], the black solid squares are the experimental results obtained in the experiment described therein. The SCM solution shows excellent agreement not only with the PWRF solution but also with the experimental data shown in solid black squares.

The second example in figure 11 computed by the SCM displays the attenuation curve for one longitudinal mode propagating along the axis of a steel pipe with the same properties as before filled with a Cannon VP8400 fluid of viscosity *η*=18.8 Pa s. The solid black squares in the background corresponding to fig. 14 of [58] display experimental results and the solid black line (lying between the two dashed lines) corresponds to the solution for the same value of viscosity *η*=18.8 Pa s obtained with DISPERSE. The dashed lines correspond to close but different values of viscosity, more details of this figure are in the original reference [58]. Once again, the solution given by the SCM shows excellent agreement with both experiment and PWRF prediction.

The above examples were also run using the *true solid* model described at the beginning, for brevity not shown here. As a consequence of the different splitting of the viscosity among the entries of the stiffness matrix, it was seen that the peaks appearing in the attenuation curves were lower than those yielded by the Stokes model. In the regions of the spectrum where experimental results were available, the solution given by the *true solid* model agrees very well with them, though it was not as good a match as with the Stokes model.

All the cases studied in this section involved a perfectly elastic pipe with either a viscoelastic solid or a viscous fluid filling the inside. Figure 12*a*,*b* show, respectively, the dispersion curves and attenuation of longitudinal modes propagating along the axis of a viscoelastic steel pipe of inner radius *r*_{i}=30 mm and thickness *h*=8 mm filled with water (modelled as an ideal fluid without viscosity). The solution given by the SCM is given in red circles and that computed by the PWRF method in solid black lines. Both solutions show very good agreement and can be compared with the perfectly elastic analogue in [21], thus showing that low material damping has little effect on the dispersion curves. The hysteretic model is used in this case and the values of the viscoelastic constants of the pipe are in the appendix. Note that the solution given by the PWRF is not complete, some of the curves had to be traced manually and even after some work it was very difficult to find the incomplete intervals. This highlights one of the advantages of the SCM, namely completeness of solutions, over conventional PWRF routines.

This section concludes with a brief comment on the ill-conditioning of the SCM when studying problems involving fluids of very low viscosity. It has been found in the course of this investigation that for viscosity values of the order of ∼0.1 Pa s or less, some pieces of the dispersion curves and the corresponding mode shapes obtained in the region close to zero frequency are not as smooth as one normally observes in more simple problems. This phenomenon has been observed in a variety of different combination of materials and geometric parameters which suggests it is an inherent feature of the SCM when used to model fluids which lie at the frontier between ideal and viscous fluids. However, a useful trend has been observed: the higher the viscosity is, the closer one can go towards zero frequency keeping smooth profiles in dispersion curves and mode shapes.

## 7. Discussion and conclusion

In this paper, an extension of the SCM to guided wave problems in viscoelastic media with the possibility of containing viscous fluids has been presented. Cases in flat, as well as in cylindrical, geometry have been studied and validated with numerous examples from the literature, solutions obtained with the conventional PWRF approach and when this approach was not possible, in the case of triclinic media for instance, a SAFE simulation was used to confirm the results. As the wavenumber *k* entered the equations nonlinearly the *Linear Companion Matrix Method* was introduced to rearrange the problem into a generalized eigenvalue problem.

The SCM shows excellent performance solving guided wave problems in viscoelastic media and some of its advantages with respect to the SAFE method or PWRF routines have become apparent in the course of the paper and are summarized below. Even though the aforementioned *Linear Companion Matrix Method* is required to recast the problem into the desired form, the SCM remains conceptually simpler than the SAFE method and the amount of previous knowledge required for its implementation is significantly less. As the examples and comments in the §§3 and 4 have highlighted, the SAFE simulation used here is known to give solutions which do not correspond to plate modes and need to be carefully filtered, mostly by manual intervention. It must be emphasized that other SAFE schemes, using one-dimensional grids for instance [39], can be used for these type of problems in which no non-plate modes in the final results have been reported so far. This phenomenon has not been observed to occur with the SCM even though some very simple post processing is useful to retain the desired propagating mode solutions. This has been done by introducing at the end of §2 the ratio *R* between the real and imaginary parts of the wavenumber. For the cases studied here, the implementation of the method and processing of physical results of the SCM has been found to be invariably easier and more straightforward than that with SAFE.

When complex roots must be found the SCM is easier to code than the conventional PWRF routines because it transforms the set of partial differential equations for the acoustic waves into a purely algebraic problem, see [21,22] for more details on this important point. In addition, the SCM is generally faster than the PWRF method as it does not require any intervention on the part of the practitioner who sometimes has to manually complete or find solutions when using the PWRF approach. This is because the SCM does not miss any eigenvalues and therefore one can be certain that the solution is complete and no modes have been missed. Mode missing is a recurrent problem encountered in PWRF routines specially when solving complicated cases involving complex roots or anisotropic materials. Finally, the SCM is capable of solving cases that are not currently solved by PWRF routines such as propagation in anisotropic media in an arbitrary direction as was pointed out at the end of §4 or as the last case in §4 shows.

In conclusion, the SCM presents itself as a powerful and robust complement to the already established SAFE and PWRF approaches. The SCM scheme, and the study of guided waves in viscoelastic generally anisotropic media, as presented in this paper is expected to be of interest in very different branches of physics dealing with guided wave phenomena and engineering applications working with composite materials and guided waves. Finally, as mentioned in the introduction, the implementation of the SCM described in this paper can also be used to calculate the non-propagating solutions, without any further development of the computational method. For brevity, such examples have not been included here.

## Data accessibility

The datasets supporting this article have been uploaded as part of the electronic supplementary material. And can also be found in the appendix to be published at the end of the paper.

## Author' contributions

F.H.Q. carried out the simulations and coding of the Spectral Collocation Method (SCM) as well as drafting the paper. Z.F. designed and carried out the SAFE model and simulations with which the SCM was validated and wrote the section of the SAFE method. M.J.S.L. and R.V.C. supervised the project as well as helped drafting the manuscript and revising.

## Competing interests

We have no competing interests.

## Funding

F.H.Q. is supported by a PhD studentship from the NDT group at Imperial College London. Z.F. is an Assistant Professor at Nanyang Technological University (Singapore) and is supported by his own institution. Both M.J.S.L. and R.V.C. are Professors at Imperial College London (UK) and are supported by their institution.

## Appendix A. Numerical data

The physical and geometrical information used for the figures presented in the main text are given here. The number of grid points *N* varies from one example to another, but it is always at least double the number of modes plotted in the figure. On a practical level *N* is chosen to achieve the shortest computation time, that is, if one is interested in the first 10 modes, running a code with *N*=100 is unnecessary; a value of *N* between 25 and 30 has consistently been shown to be sufficient.

The parameters for the plate of figure 2*a*,*b* are as follows (with the usual axes orientation shown in figure 1):
*h* stands for the thickness of the plate. The elastic stiffness matrix is given in GPa

The parameters for the plate of figure 4*a*,*b* are as follows (with the usual axes orientation shown in figure 1):
*h* stands for the thickness of the plate and

The parameters for the plate of figure 5*a*,*b* are the same as above, equations (7.4)–(7.6) except for the thickness which now is *h*=50 mm.

The multi-layer system of figure 6*a*,*b* is composed of three layers. The top layer has a thickness of *h*=8 mm and the rest of the properties are given in (7.4)–(7.6). The middle layer has a thickness of *h*=5 mm and the other properties are given in (7.1)–(7.3). The bottom layer parameters are as follows:
*h* stands for the thickness of the plate. The elastic stiffness matrix is given in GPa

The parameters for the cylinder of figure 7*a*,*b* are as follows (with the usual axes orientation shown in 1):
*h* stands for the thickness of the plate and

The parameters for the cylinder of figure 8*a*,*b* are as follows (with the usual axes orientation shown in 1):
*h* stands for the thickness of the plate and *r*}-axis must performed on these matrices so that the fibres follow an helicoidal path around the cylinder.

The parameters for the cylindrical system of figure 9*a*,*b* are as follows. For the outer elastic steel cylinder

For the system in figure 10, the parameters for the steel pipe are

The fluid inside is glycerol and has the following properties:

For the system in figure 11, the parameters for the steel pipe are the same as for figure 10. The fluid inside is Cannon VP8400 and has the following properties

For the system in figure 12*a*,*b* the parameters for the steel pipe are

- Received April 22, 2015.
- Accepted September 22, 2015.

- © 2015 The Author(s)