Lamb waves are elastodynamic guided waves in plates and are used for non-destructive evaluation, sensors, and material characterization. These applications rely on the knowledge of the dispersion characteristics, i.e., the frequency-dependent wavenumbers. The interaction of a plate with an adjacent fluid leads to a nonlinear differential eigenvalue problem with a square root term describing exchange of energy with the surrounding medium, e.g., via acoustic radiation. In this contribution, a spectral collocation scheme is applied to discretize the differential eigenvalue problem. A change of variable is performed to obtain an equivalent polynomial eigenvalue problem of fourth order, which is linear in state-space and can reliably be solved using modern numerical methods. Traditionally, the leaky Lamb wave problem has been solved by finding the roots of the characteristic equations, a numerically ill-conditioned problem. In contrast to root-finding, the approach described in this paper is inherently able to find all modes and naturally handles complex wavenumbers. The full phase velocity dispersion diagram and attenuation curves are presented and are shown to be in excellent agreement with solutions of the characteristic equation as well as computations made with a perturbation method. The procedure is applicable to anisotropic, viscoelastic, inhomogeneous, and layered plates coupled to an inviscid fluid.

## I. INTRODUCTION

Lamb waves are guided modes confined inside a plate. Understanding their propagation characteristics is important for diverse applications. In the field of non-destructive evaluation, they can be used to achieve long-range inspection of plate-like structures and pipes.^{1} Sensors and actuators based on Lamb wave propagation include flow meters, liquid viscosity sensors, acoustic touch screens, and microfluidic pumping.^{2,3} Similar guided modes in layered media have also been studied extensively in seismology.^{4}

In practice, the plate serving as a waveguide is often in contact with a fluid, either on one or both sides. The mutual interaction between the guided wave in the plate and the pressure wave in the fluid will, henceforth, be referred to as *fluid loading* of the plate. While the energy of waves propagating in a free plate is strictly confined within the structure, fluid loading will usually lead to leakage of energy into the surrounding medium by means of acoustic radiation. Such waves are called leaky Lamb waves.^{5–8} Like most guided waves, they are dispersive, meaning that their wavenumber depends on frequency.^{9,10} Knowledge of the dispersion characteristics is essential for the design of ultrasonic devices for the mentioned applications. In particular, both phase velocity and attenuation due to energy leakage are important parameters.^{1}

The wavenumber-frequency relationship, given implicitly by the characteristic equation of the guided wave problem,^{10,12} is transcendental and cannot be solved analytically. Traditionally, numerical root-finding has been performed to obtain dispersion diagrams.^{13} This, however, is a numerically ill-conditioned problem.^{14} Modern techniques to obtain dispersion characteristics are based on discretization of the guided wave boundary value problem and subsequent solution of the resulting algebraic eigenvalue problem. This discretization-based methodology has significant numerical advantages: it is robust, guarantees to find all eigenvalues, and naturally handles complex wavenumbers.^{15} Different discretization methods have been used for this kind of solution procedure. Spectral collocation with Chebyshev polynomials was used by Adamou and Craster^{15} and later by Hernando Quintanilla *et al.*^{16,17} Pioneering work was done by Lefebvre *et al.*^{18,19} using Legendre and Laguerre spectral collocation. Pagneux and Maurel^{20} presented a spectral numerical method with explicit separation into symmetric and antisymmetric mode families. Finite difference schemes have also been used for this purpose.^{21} Waveguides with arbitrary cross-section are often discretized using finite elements.^{22,23}

The above described efforts mostly focus on the solution of the free waveguide problem. A waveguide adjacent to an infinite medium represents an open-domain problem and is, hence, considerably more difficult to handle. Commonly, fluid loading of Lamb waves has been studied by either perturbation methods^{10} or by finding the roots of the corresponding characteristic equation.^{6,24,25} Perturbation based on reciprocity relations provides a simple way to obtain approximative attenuation values due to fluid loading given the free plate modes. The results are only valid for fluids with low mass density and the method is not able to find any additional modes arising in the fluid loaded case. The second approach, root-finding of the characteristic equations, is potentially able to find the complex wavenumbers of all leaky modes.^{12,13} However, the numerical challenges mentioned for the free plate still apply and are even more severe because most wavenumbers are truly complex. Root-finding of the leaky wave problem is a cumbersome method, which might miss important parts of the dispersion diagram.

In order to avoid these complications, efforts have been made to compute dispersion characteristics of leaky guided waves by solution of the corresponding eigenvalue problem. The main difficulty herein is that fluid loading leads to a nonlinear eigenvalue problem.^{23,26} Numerical algorithms for general nonlinear eigenvalue problems are the subject of current research^{27–29} and several different solvers have been proposed.^{30–36} Polynomial eigenvalue problems represent a special case, as they can be reduced to a generalized linear eigenvalue problem in a higher dimensional state-space. They are, therefore, easier to solve because conventional, highly sophisticated numerical solvers for linear eigenvalue problems might be used.

Hayashi and Inoue^{22} solved the fluid loaded plate problem by discretizing the waveguide with finite elements and assuming plane harmonic waves in the surrounding medium. They further use *a priori* knowledge on the propagation symmetry of leaky Lamb waves to reduce the nonlinear eigenvalue problem to a generalized eigenvalue problem in terms of the transversal wavenumber. This method should be capable of finding all solutions but does not resolve the ambiguity in the quadratic relationship of the wavenumber components. Mazzotti *et al.*^{23,37} calculated the dispersion characteristics of leaky waveguides of arbitrary cross-section. They combine a finite element discretization of the waveguide with a boundary element discretization of the fluid domain. Subsequently, they solve the resulting nonlinear eigenvalue problem by a contour integral method, which requires *a priori* knowledge about the regions in the complex plane where solutions may lie. The assumptions made hereon will not be correct when the guided mode's wavenumber is close to the fluid bulk wavenumber and may lead to inaccurate and incomplete results in this region.

The nonlinear eigenvalue problem can be avoided by additionally discretizing the infinite surrounding medium. Hernando Quintanilla^{38} implemented a perfectly matched layer to obtain a generalized eigenproblem formulation. However, the unnecessary discretization of the fluid domain leads to spurious eigenvalues, which have to be discarded appropriately. Additionally, a bad design of the perfectly matched layer potentially leads to detrimental effects. Another tempting approach would be to combine a Chebyshev spectral collocation of the plate with a Laguerre spectral collocation of the fluid half-space. Nevertheless, this does not work because the radiated acoustic field is known not to be square integrable and can, thus, not be represented by Laguerre polynomials.^{19}

In this contribution, we apply a Chebyshev spectral collocation to discretize the waveguide problem with spectral accuracy. Fluid loading is then incorporated by assuming inhomogeneous plane wave radiation into the medium. The resulting nonlinear eigenvalue problem is transformed to an equivalent polynomial eigenvalue problem of fourth order by an appropriate change of variable. Hence, a formulation is obtained which can be solved robustly by conventional numerical methods. The resulting computational method is easily implemented and able to reliably calculate the full dispersion characteristics of the leaky Lamb wave problem with analytically exact fluid interaction.

## II. PROBLEM FORMULATION

In the following, let us consider a homogeneous, isotropic, linear-elastic plate of infinite extend. An excerpt of the cross-sectional geometry of the problem is shown in Fig. 1. The coordinate system is chosen such that the wave propagation takes place in *x*-direction and the *y* axis is normal to the plate's plane. In this case, the wave motions with displacements in *z*-direction decouple from those with displacements in the *x–y*-plane.^{39} Solely the latter are capable of interacting with the surrounding ideal fluid and are referred to as *Lamb waves*. Hereafter, only this kind of wave propagation shall be considered.

The displacements $u\u2192$ in the *x–y*-plane may be represented by a Helmholtz displacement decomposition^{3,9} as

inside the plate and as

in the fluid domain. Herein, *φ*_{p} and *φ*_{f} represent longitudinal wave displacement potentials and Ψ_{z} designates the *z*-component of the transverse wave potential, respectively. We consider harmonic wave propagation in both time *t* and in the coordinate *x*. This leads to

*k*is the guided mode wavenumber and

_{x}*ω*= 2

*πf*is the angular frequency corresponding to the frequency

*f*.

*φ*(

*y*) and Ψ(

*y*) designate the unknown

*y*-dependent part of the plate potentials. Being time-independent, they represent a through-thickness standing wave in the plate. The potentials in the plate and the fluid share the same harmonic factor $ei(kxx\u2212\omega t)$ in order to satisfy the generalized Snell's law,

^{40}that is, $kx,fluid=kx,lamb=:kx$. Since this factor appears in all equations, it is omitted from now on. The factor $i=\u22121$ in the transversal wave potential is chosen for practical reasons, as the problem then results in a purely real formulation.

Furthermore, the fluid motion is assumed to be an inhomogeneous plane harmonic bulk wave^{40} with unknown complex amplitude *A* and complex wave vector $k\u2192f=[kx,ky]T$. In contrast to guided waves in the plate, which propagate in the *x*-direction, this wave vector has a non-vanishing component normal to the plate that is called the transversal wavenumber *k _{y}*. The magnitude $kf:=|k\u2192f|=\omega /cf$ for an inviscid fluid is real valued and fully determined by the fluid's longitudinal wave speed

*c*

_{f}and the angular frequency

*ω*. Moreover, the components of the wave vector must satisfy the trigonometric relation

With these assumptions and considering that the plate-fluid interface is at *y*_{0} = *h*/2, where *h* is the plate's thickness, the fluid potential is then known to be of the form^{40}

With the above ansatz for the potentials, namely, Eqs. (3a), (3b), and (5), the unknowns are *φ*(*y*), Ψ(*y*), and *A*. For the sake of conciseness, we omit the *y*-dependence of the field variables from now on in the notation unless necessary. The differential equations of *unforced motion*^{11} of the plate may then be formulated in dependence of the unknowns as

*c*

_{l}and the transversal wave speed

*c*

_{t}of the plate's material have been used.

Boundary and interface conditions are needed in addition to the equation of motion [Eq. (6)]. For single-sided fluid loading of the plate, the top boundary in Fig. 1 is traction free,^{10} requiring the normal stress *σ _{yy}* and shear stress

*σ*to vanish at

_{xy}*y*= −

*h*/2. In terms of the unknowns, the conditions are (see Appendix)

*λ*and

*μ*denote the Lamé parameters of the plate's material. The linear displacement-strain relation

^{11}and the linear-elastic constitutive relations

^{11}have been used assuming a plain strain state,

^{41}i.e.,

*ε*=

_{zy}*ε*=

_{zx}*ε*= 0. Additionally, for numerical reasons and conciseness, the above expressions have been normalized to

_{zz}*μ*and i

*μ*/2, respectively.

The interface conditions at the boundary *y* = *h*/2 between the plate and the inviscid fluid are given by the balance of tractions and the continuity of normal displacements.^{24} The first condition is ensured by the continuity of normal stress in the plate *σ _{yy}* and the fluid

*σ*

_{f}

_{yy}, given by

*σ*−

_{yy}*σ*

_{f}

_{yy}= 0, and vanishing of shear stresses, i.e.,

*σ*= 0. The second is guaranteed by demanding that

_{xy}*u*−

_{y}*u*

_{f}

_{y}= 0. Normalizing the two stress requirements as before and the displacement condition with

*h*, these conditions result in (see Appendix)

^{9}

The equations of motion [Eq. (6)] together with the boundary conditions [Eq. (7)] and interface conditions [Eq. (8)] are further referred to as the *leaky Lamb problem* for single sided fluid interaction. In comparison with the free plate problem, it has merely been increased by the scalar degree of freedom *A* and one additional equation given by Eq. (8c). In contrast to other work,^{22,26} we do not reduce this additional degree of freedom. It should be noted that due to the assumption of inhomogeneous plane wave propagation in the fluid, the acoustic field is fully determined by the values at the interface *y* = *h*/2. In this way, the computational domain has been reduced to a simple one-dimensional closed geometry, namely, the interval *y* ∈ [−*h*/2, *h*/2]. This means that no discretization of the fluid domain is required in order to solve the equations. As a consequence, spurious modes due to discretization of the open domain *y* ∈ [*h*/2, *∞*] are avoided.^{38}

The leaky Lamb problem for the *immersed plate*, i.e., loaded with the same fluid on both sides, can be obtained in a similar way by introducing an additional scalar degree of freedom representing the potential amplitude in the second fluid half-space and prescribing interface conditions similar to Eq. (8) at *y* = −*h*/2 instead of the free boundary conditions.

For a specific angular frequency *ω*, one may regard the problem as a differential eigenvalue problem,^{20} where the eigenvalues *k _{x}* and eigenfunctions

*q*(

*y*) = [

*φ*(

*y*), Ψ(

*y*),

*A*]

^{T}need to be found. In the interface condition [Eq. (8c)], the transversal wavenumber

*k*still appears explicitly. Using Eq. (4), it could be rewritten in terms of

_{y}*k*as

_{x}Due to this square root dependence, the differential eigenvalue problem is nonlinear.^{29} We remark that due to the ambiguity of the sign in Eq. (9), technically two problems need to be solved in order to find all solutions. As will be seen in Sec. IV, the proposed change of variable that is performed after discretization resolves this ambiguity.

## III. DISCRETIZATION

The leaky Lamb problem is continuous in the coordinate *y*. In order to solve these equations on a computer, a discrete approximation to the continuous functions of this coordinate as well as their derivatives needs to be found. Independent of the choice of discretization technique, the structure of the resulting discrete eigenvalue problem is similar because it is induced by the continuous formulation. The problem has been formulated on a simple one-dimensional and closed domain. Such problems are often conveniently discretized by spectral collocation.^{42,43} Chebyshev spectral collocation will be used in the present contribution, in accordance with the work on free guided waves by Adamou and Craster^{15} and Hernando Quintanilla *et al.*^{16} In the following, we shortly review the features and advantages of this method based on the monographs by Trefethen^{42} and Fornberg.^{43}

Spectral collocation is a special form of the weighted residuals method. The basic idea is to expand the solution as a sum of *N* weighted and known continuous trial functions. Then, an approximate solution to the continuous problem is found by determining the weights of the expansion terms that yield the minimal residuum. A *collocation* method^{43} determines these weights by requiring the residual of the approximation to vanish at a selected set of *N* points, which are the so called collocation points *y _{i}*. In contrast to many other methods, it is implemented in physical space. The term

*spectral*refers to the fact that smooth global functions are used as test functions,

^{43}i.e., they extend over the whole domain of the problem.

Different test functions can be used to expand the approximate solution. On finite domains, Chebyshev polynomials show optimal behavior for *N* → *∞*.^{43} In order to achieve the best performance and avoid the Gibbs phenomenon, the collocation points must be clustered quadratically at the domain borders.^{42} Chebyshev-Gauß-Lobatto points^{43} fulfill this requirement and are used together with Chebyshev polynomial test functions.

In contrast to finite elements, which exhibit polynomial convergence with increased discretization order *N*, spectral methods provide exponential convergence rates for sufficiently smooth solutions.^{42,43} This property is called *spectral convergence* or *spectral accuracy*. For problems with non-smooth solutions, spectral methods usually still converge polynomially. Spectral convergence is the main reason why these methods are considered superior to other discretization methods on simple domains.^{14,42,43} These features are obtained at cost of dense matrices for the discretized problem. However, the rapid convergence of the method usually compensates for this disadvantage as it leads to small matrices.

Differentiation with spectral methods can be represented explicitly by differentiation matrices.^{42,47} To illustrate this, consider any arbitrary function *g*(*y*). If the function values *g*(*y _{i}*) at the collocation points

*y*are collected in the vector $g\xaf$, then the vector $g\u2032\xaf$ of corresponding derivative values $\u2202g(y)/\u2202y|yi$ can be computed with the differentiation matrix $D\xaf\xafy$ as $g\u2032\xaf=D\xaf\xafy\u2009g\xaf$. This provides a straightforward approach to assemble the discretized form of the leaky Lamb problem. We shortly describe the procedure as presented in references.

_{i}^{15,16,47}As a first step, the unknowns in Eqs. (6), (7) and (8) are collected in

*q*(

*y*) = [

*φ*(

*y*), Ψ(

*y*),

*A*]

^{T}and all equations are rewritten in matrix form depending only on

*q*(

*y*). In a second step, the following formal mappings are performed:

unknowns: $q(y)=[\phi (y)\Psi (y)A]\u21a6q\xaf=[\phi \xaf\Psi \xafA]$,

all constant elements

*ξ*in the equations: $\xi \u2009\u21a6\u2009\xi I\xaf\xaf$,diff. operators: $(\u2202/\u2202y)\u21a6(1/h)D\xaf\xafyn,\u2009(\u22022/\u2202y2)\u21a6(1/h2)D\xaf\xafyyn$,

where $I\xaf\xaf$ denotes the *N* × *N* identity matrix, whereas $D\xaf\xafyn$ and $D\xaf\xafyyn$ refer to the *N* × *N* differentiation matrices of first and second order on the normalized domain *y* ∈ [−1/2, 1/2], respectively. For the boundary and interface conditions, only local equations are needed. For this purpose, $D\xafyn|yi,\u2009D\xafyyn|yi$, and $I\xaf|yi$ shall represent the rows of the corresponding matrices describing the maps for *y* = *y _{i}*. Performing the two steps as explained above and multiplying all equations by

*h*

^{2}yields the corresponding discretized and normalized form of the leaky Lamb problem.

Denoting with “0” the zero matrix with appropriate dimension, the equation of motion [Eq. (6)] results in a system of 2*N* + 1 equations and 2*N* + 1 variables given by

The boundary conditions [Eq. (7)] give a 2 × 2*N* + 1 system for *y* = −*h*/2, namely,

The discrete form of the interface conditions [Eq. (8)] yields a 3 × 2*N* + 1 system for *y* = *h*/2, written as

Lastly, we incorporate the boundary conditions (BCs) and interface conditions (ICs) into the equations of motion to obtain one system of equations.^{47} In the equations of motion [Eq. (10)], the rows 1 and *N* + 1 represent equations at the top boundary point *y* = −*h*/2, while rows *N*, 2*N* and 2*N* + 1 correspond to the point at *y* = *h*/2. These rows are replaced with the boundary conditions [Eq. (11)] and interface conditions [Eq. (12)], respectively. Therefore, for each dependence on the wavenumbers *k _{x}* and

*k*, a matrix is assembled with structure as shown in Fig. 2.

_{y}Naming the resulting matrices $A\xaf\xaf2,\u2009A\xaf\xaf1,\u2009A\xaf\xaf0$, and $B\xaf\xaf$ and introducing the matrix function

the above scheme yields

## IV. CHANGE OF VARIABLE

Several approaches exist to solve nonlinear eigenvalue problems, which are often based on approximation by iterative linearization.^{5,26,36} This is not necessary for polynomial and rational eigenvalue problems because they are linearizable, meaning that an equivalent linear eigenvalue problem in a higher-dimensional state space can be found.^{29} The *gun problem* models a radio-frequency gun cavity and consist of an eigenvalue problem with two square root terms. As shown by Hood,^{44} it can be transformed to a polynomial eigenvalue problem by performing a trigonometric change of variable. As we demonstrate below, the leaky Lamb problem [Eq. (14)] can also be reduced to polynomial form by a similar change of variable.

For this purpose, the variable $\gamma \u2208\u2102\u2216{0}$ is defined such that

holds. From Eq. (9), the corresponding transversal wavenumber is

We now insert the change of variable Eqs. (15) and (16) into the matrix function [Eq. (13)] and multiply by 4*γ*^{2}. Depending on the choice of sign, this leads to two different polynomial matrix functions, given by

It is noted that the matrix polynomials $P\xaf\xaf\xb1(\gamma )$ will be singular for *γ _{n}* if and only if the corresponding

*k*given by Eq. (15) is an eigenvalue of $F\xaf\xaf(kx)$ or

_{xn}*γ*= 0. Hence, the spectrum of $F\xaf\xaf(kx)$ is fully determined by the eigenvalues of $P\xaf\xaf\xb1(\gamma )$. However, additional eigenvalues

_{n}*γ*= 0 are introduced by $P\xaf\xaf\xb1(\gamma )$, which yield wavenumbers

_{n}*k*→

_{xn}*∞*. These are of no further interest and are easily discarded.

It is pointed out that the choice of sign in Eq. (16) and consequently Eq. (18) is irrelevant. By inspecting the coefficient matrices in Eq. (18), we find that $P\xaf\xaf+(\gamma )=\gamma 4P\xaf\xaf\u2212(\gamma \u22121)$. Therefore, the spectrum of $P\xaf\xaf\u2212$ equals the inverted spectrum of $P\xaf\xaf+$, i.e., if *γ _{n}* is an eigenvalue of $P\xaf\xaf+$, then $\gamma n\u22121$ will be an eigenvalue of $P\xaf\xaf\u2212$. Nevertheless, these two solutions yield the same wavenumber

*k*when substituted back into Eq. (15). They will also map to the exact same transversal wavenumber

_{xn}*k*when taking into account the sign chosen in Eq. (16). As a result, the two polynomial eigenvalue problems are equivalent and each of them

_{yn}*fully*and

*uniquely*describes the spectrum of leaky Lamb modes. This feature is due to the fact that the change of variable is not an injective map and the values

*γ*and $\gamma n\u22121$ actually represent the choice of sign in Eq. (16).

_{n}We conclude that in order to obtain the full leaky Lamb spectrum, it is possible to choose without loss of generality $P\xaf\xaf=P\xaf\xaf+$ and then solve the polynomial eigenvalue problem

obtaining the eigenvalues *γ _{n}* and eigenvectors $q\xafn$. Subsequently, the sought wavenumbers

*k*are computed using Eq. (15), while the eigenvectors remain unchanged. Similarly, the corresponding transversal wavenumber

_{xn}*k*may be determined uniquely with Eq. (16) by choosing the sign corresponding to the choice of Eq. (19). This uniqueness is in fact an additional major advantage of the proposed method and is crucial for the solution procedure.

_{yn}To solve the polynomial eigenvalue problem [Eq. (19)], a *companion linearization*^{27,45,46} is usually performed, which is not an approximation but rather an equivalent representation in a higher dimensional state space. Such linearization may be obtained by introducing the state variables

as well as the companion matrices

In doing so, Eq. (19) is equivalent to

which is a generalized linear eigenvalue problem. It is noted that this linearization procedure is not unique, i.e., other equivalent linear forms can be constructed.^{27,45} After linearization, conventional eigenvalue solvers can be employed to obtain the desired solutions, harnessing all advantages these methods provide. In contrast to root-finding methods, the presented eigenvalue solution approach reliably finds all wavenumbers of the leaky Lamb problem and additionally provides the associated eigenvectors free of extra cost. Moreover, the above formulation inherently handles complex wavenumbers without extra effort.

## V. CALCULATION RESULTS

The matrix $A\xaf\xaf0$ in Eq. (18) depends on the frequency *f*. Solving this problem for different frequencies and substituting the obtained eigenvalues *γ* back into Eq. (15) results in the complex wavenumber spectrum *k _{x}*(

*f*). The differentiation matrices are computed with a routine provided by Weideman and Reddy.

^{47}We use the matlab polyeig-function to compute the eigenvalues, which is based on companion linearization in combination with the QZ-algorithm. A discretization order of

*N*= 10 yields converged results. Computing the solutions for 150 frequency values requires approximately one second on a regular personal computer with an Intel Core i5 processor and 8 GB of main memory.

The computations have been performed for a brass plate (*ρ* = 8440 kg/m^{3}, *λ* = 87 GPa, *μ* = 41 GPa) in contact on one side with water (*ρ*_{f} = 1000 kg/m^{3}, *c*_{f} = 1480 m/s). All results are shown for a 1 mm thick plate. We note, however, that similar to Lamb waves in a free plate, the wavenumber-thickness product *k _{x}h* depends only on the frequency-thickness product

*fh*and not explicitly on the thickness

*h*. The resulting wavenumbers

*k*are shown in Fig. 3(a), where Re{•} and Im{•} refer to the real and imaginary parts of •, respectively. Accordingly, the transversal wavenumber

_{x}*k*in the fluid domain has been calculated by Eq. (16) and is plotted in Fig. 3(b). The dispersion diagrams are actually curves in the three-dimensional space given by [Re{•} × Im{•}×

_{y}*f*].

^{10,17}The two graphs in Fig. 3 are projections onto the complex plane $[Re{\u2022}\xd7Im{\u2022}]$. A perspective view in three-dimensional space of the dispersion curves of

*k*and

_{x}*k*is depicted in Figs. 4(a) and 4(b), respectively. Additionally, the spectrum of the free plate is displayed as reference in Fig. 4(a).

_{y}The wavenumber spectrum *k _{x}*(

*f*) exhibits the expected Hamiltonian symmetry, i.e., if

*k*

_{x}_{0}is an eigenvalue, then ${\u2212kx0,kx0*,\u2212kx0*}$ will also be eigenvalues, where •

^{*}denotes the complex conjugate. The axes extends of the plots are chosen such that only the spectra of propagating modes, having $|Im\u2009kx|<0.1\u2009Np/mm$, are shown. In this region, the free plate spectrum is purely real, as these modes are unable to leak energy into the surroundings and, hence, do not experience attenuation. It is important to note that fluid loading leads to

*additional modes*not present in the free plate, a fact that is ignored by perturbation methods. The present method correctly finds these additional modes. One of them is the well known quasi-Scholte (QS) mode,

^{48}which is sometimes also called the Stoneley-Scholte plate wave or A-mode.

^{49}

In the following, we solely consider the open half plane $kx+\u2208{kx|Re\u2009kx>0}$, which represents guided waves traveling in positive *x*-direction. Of these modes, the ones with Im *k _{x}* > 0 radiate a plane wave into the infinite fluid domain, being attenuated thereby. These are further referred to as leaky modes. The reciprocal process, an inhomogeneous plane wave incident on the plate whereby a guided mode is excited, is described by the complex conjugate wavenumbers seen in the quadrant with Im

*k*< 0. On the real half-axis, having Im

_{x}*k*= 0, modes propagating without attenuation are found. The other half-plane with Re

_{x}*k*< 0 is symmetric to the first one because the analogous modes propagate in opposite direction.

_{x}The phase velocity is calculated as

and is shown in Fig. 5(a). The color of the points reflects the corresponding attenuation Im*k _{x}*, which is additionally provided explicitly in Fig. 5(b). For the present case of a plate with high mass density

*ρ*compared to the fluid, the phase velocity dispersion diagram is similar to the one corresponding to waves in a free plate. For this reason, the leaky Lamb modes are termed $Ai\u2032$ or $Si\u2032$ whenever a correspondence to anti-symmetric or symmetric Lamb waves in a free plate of order

*i*can be established. It should be pointed out, however, that the mode structure

*u*(

_{x}*y*) and

*u*(

_{y}*y*) of the single-sided leaky Lamb problem does not strictly exhibit any kind of symmetry because the problem domain itself is not symmetric along the

*y*-coordinate.

^{39}While the phase velocity in the present case is similar to the free plate, attenuation is solely a consequence of fluid loading—the latter often being an important parameter for system design.

The low frequency region of the phase velocity dispersion curves is enlarged in Fig. 6 for better clarity. The frequency where the phase velocity of the A0-mode in the free plate coincides with the fluid wave speed *c*_{f} is called the coincident frequency. In the present case it is at 0.49 MHz. It is often assumed *a priori* that no radiation will occur below this frequency.^{23,26} However, as the attenuation diagram in Fig. 5(b) indicates, energy leakage may take place below the coincident frequency, a fact that has been noted before.^{50} Moreover, we observe splitting of the A0′-mode into two branches with real wavenumbers near 0.33 MHz, named D0 and D1 here. This splitting of the A0′-mode has been analyzed by Dabirikhah and Turner^{50} and was also observed by Bao *et al.*^{51} as well as Freedman.^{52}

In order to validate the presented calculation method, it is checked against solutions obtained by two different and independent models. The first one is root-finding of the characteristic equation, which should yield the same results as the presented model. The equation is provided, e.g., in Überall *et al.*^{49} as well as Grabówska.^{24} It is difficult to reliably find all solutions of these transcendental equations and initial values are needed for the iterative solution procedure. For this purpose, random values between −1% and +1% of *k _{x}* are added to the wavenumbers calculated with the presented spectral collocation method. These perturbed wavenumbers are then used as initial values for the Newton iteration applied to find the closest root. The solutions all converge back to the original values of

*k*, as can be seen in the phase velocity plot in Fig. 6. The relative difference between phase velocities calculated by the two methods is always less than 4 × 10

_{x}^{−5}. Furthermore, the colored background of the diagram represents the local magnitude of the characteristic equation for real wavenumbers. The light regions are close to zero and roots are expected to lie there, which is the case for the obtained solutions.

The second check is performed against a totally different model, namely, perturbation of the free plate solutions. According to Auld,^{10} approximative attenuation values *α* ≈ Im*k _{x}* can be computed given the phase velocity

*c*

_{p0}as well as the displacement structure

*u*

_{y}_{0}(

*y*) of the free plate modes as

where *P*_{0} denotes the power flux inside the plate. These computations were performed for the same plate as before and the results are shown in Fig. 5(b). The obtained attenuation values largely coincide with those computed from the leaky Lamb eigenvalue problem. As expected, the coincidence region cannot be represented by the perturbation method, nor the high-frequency range of the A0′ and S0′ modes.

Last, we compute the modes of the brass plate immersed in water. The results are shown in Fig. 7. An additional, non-dispersive mode is obtained, known as the second quasi-Scholte mode or S-mode.^{53} The attenuation of leaky modes is roughly doubled, as expected from perturbation theory. An exception is the A0″-mode in the very low and high frequency range, where it exhibits substantially different attenuation values compared to single-sided fluid loading.

## VI. CONCLUSIONS

Leaky Lamb waves are modeled by a nonlinear eigenvalue problem. These are, in general, difficult to solve. However, the leaky Lamb problem is linearizable by performing a change of variables. An equivalent polynomial eigenvalue problem of fourth order was obtained, which is linear in state-space and can be solved seamlessly by conventional methods. The found eigenvalues then lead to the corresponding wavenumbers without difficulty.

The main advantage of the presented solution procedure is that it robustly finds all eigenvalues, including in the region of coincidence. As the modeled fluid loading is exact, it is valid for any combination of plate material and fluid. The method is easy to implement and makes use of common numerical algorithms. This allows adaptation for custom requirements, e.g., real time computation of one or more target wavenumbers by using an appropriate eigenvalue solver. The technique is not specific to spectral collocation, but it could in fact be used in combination with most other discretization methods, including finite elements. The method could be extended to viscoelastic, anisotropic, inhomogeneous, and layered plates. However, in the current form, it is restricted to plain geometries of infinite extend because it assumes inhomogeneous plane wave radiation. It might also not be suitable for modeling radiation into viscous fluids, solid media or when two different fluids are in contact with the plate.

## ACKNOWLEDGMENTS

The Chair of Sensor Technology is grateful for financial support by Diehl Metering GmbH.

### APPENDIX: COMPUTATION OF STRAIN AND STRESS

In order to obtain a consistent formulation, stresses and displacements in the boundary and interface conditions need to be written in terms of the three degrees of freedom *φ*, Ψ, and A. The displacements are obtained by inserting the ansatz for the potentials in Eqs. (3) and (5) into the Helmholtz decompositions, Eqs. (1) and (2), yielding

inside the plate and

for the fluid normal displacements. The non-zero strain components are given by

The two relevant stresses are

## References

_{0}and interface Scholte modes in fluid-loaded plates