## Abstract

The problem of Lamb wave propagation in waveguides with varying height is treated by a multimodal approach. The technique is based on a rearrangement of the equations of elasticity that provides a new system of coupled mode equations preserving energy conservation. These coupled mode equations avoid the usual problem at the cut-offs with zero wavenumber. Thereafter, we define an impedance matrix that is governed by a Riccati equation yielding a stable numerical computation of the solution. Incidentally, the versatility of the multimodal method is exemplified by treating analytically the case of slowly varying guide and by showing how to get easily the Green tensor in any geometry. The method is applied for a waveguide whose height is described by a Gaussian function and the energy conservation in verified numerically. We determine the Green tensor in this geometry.

## 1. Introduction

The interest in waves propagating in elastic waveguide comes, at least, from two applications. First in non-destructive testing, guided waves are believed to have a potential for improving inspection efficiency and sensitivity, compared with bulk-wave technique. Second in geophysics, they are the basic picture for seismic surface waves propagating in the crust and upper mantle.

To tackle the problem of waveguide thickness variation, different techniques have been proposed, such as hybrid boundary-element methods (Cho & Rose 1996; Cho 2000; Galan & Abascal 2003), finite-element methods (Koshiba *et al*. 1984; Galan & Abascal 2002; Wu *et al*. 2003) or modal methods (Abram 1974; Kennett 1984; Maupin 1988; Tromp 1994; Galanenko 1998; Folguera & Harris 1999). The modal approach offers the advantage of discretizing the transverse direction with transverse modes which implies, for instance, an exact solution for straight waveguide, and this method permits the reduction of the problem to an ordinary differential equation that governs the modal components resulting from the projection onto the modal basis.

The problem of wave propagation in a two-dimensional waveguide can be formally written as an evolution problem , where *x* is the waveguide axis, (*u*, *v* are the displacement components and *s*≡_{xx}, *t*≡_{xy}, with the stress tensor) and a differential operator. After has been diagonalized, this vector ** z** can be projected on the usual Lamb modes as , where satisfies the eigenvalue problem . Then, using the biorthogonality relation satisfied by the eigenfunctions

*z*_{n}, the original evolution problem is thus reduced to a first-order ordinary differential equation on the components

*c*

_{n}(e.g. Maupin 1988).

There are two problems in this coupled mode method. The first problem concerns cut-offs that occur at isolated longitudinal coordinates *x*. At cut-offs, the eigenfunctions *z*_{n} cannot form a base anymore. This problem is related to the normalization coefficient of the biorthogonality relation that vanishes at the cut-offs. It is also related to the impossibility to decompose ** z** into forward and backward modes at cut-offs. Relatively few works have been done to overcome this difficulty. Galanenko (1998) proposes to treat cut-offs by using the regular singularity theory of differential equation. The same kind of techniques has also been recently proposed by Perel

*et al*. (2005). The second problem concerns the numerical implementation of the coupled mode equation which is not obvious. Indeed, the coupled mode equation corresponds to a boundary-value problem and the presence of evanescent modes makes the integration of the coupled mode equation numerically unstable if directly performed. To avoid this difficulty, Kennett (1984) used the technique of invariant embedding to obtain coupled Riccati equations on the reflection and transmission matrices.

In this paper, we propose a new formulation of the coupled mode method that is numerically stable and that partially solves the problem of cut-off (namely, the cut-offs with vanishing wavenumber *k*_{n}). Our formulation uses the symmetry properties of the Lamb modes between forward and backward modes. It permits one to decompose the two two-vectors ** X**=(

*u*,

*t*)

^{T},

**=(−**

*Y**s*,

*v*)

^{T}into two biorthogonal bases , , which correspond to the forward modes only, instead of decomposing the four vector

**into both backward and forward modes. It is then possible to renormalize separately and , so that the renormalized**

*z*

*X*_{n}and

*Y*_{n}remain bases at cut-offs with zero wavenumber. This gives two coupled evolution equations on

*a*

_{n}and

*b*

_{n}, the components of

**and**

*X***on the two renormalized bases. To avoid numerical instability, we define an impedance matrix , which links the components through**

*Y***=**

*b***, and which is governed by a Riccati equation. This kind of technique has been used in the scalar case, where the impedance matrix is equivalent to a Dirichlet to Neumann operator (Pagneux**

*a**et al*. 1996). We could say that the matrix for Lamb waves is an extension of the concept of the Dirichlet to Neumann operator to vector waves.

The rest of this paper is organized as follows. In §2, we present the derivation of the coupled equations. In §3, we define the impedance matrix and the Riccati equation. In §4, the energy conservation is discussed. In §5, some analytical solutions for slowly varying guides are presented. In §6, the derivation of the Green tensor in elastic waveguide is presented. Section 7 presents the numerical implementation that will be used to get the results of §8.

## 2. Coupled mode equations

In this section, the problem of guided wave is set using a modal expansion. The modal expansion we choose has been previously used in Pagneux & Maurel (2002, 2004). It is different from the coupled mode expansion classically used; indeed, one could think to write the problem as ** z**′=

**(for instance, as in Maupin 1988), where**

*z***=(**

*z**u*,

*v*,

*s*,

*t*)

^{T}is a four-vector and to project on the usual Lamb modes to obtain an evolution problem for the coefficient

**of**

*c***. Instead, we choose to split the four-vector**

*z***into two two-vectors**

*z***=(**

*X**u*,

*t*)

^{T},

**=(−**

*Y**s*,

*v*)

^{T}and to perform the projection onto

*X*_{n}=(

*U*

_{n},

*T*

_{n})

^{T},

*Y*_{n}=(−

*S*

_{n},

*V*

_{n})

^{T}, where

*X*_{n},

*Y*_{n}are renormalized and build from the usual Lamb modes ,

A first motivation to do that comes from Fraser's biorthogonality relation (1976) that simply translates into when the form is defined. Also, the main motivation is to take care of the mode cut-offs (see for instance Perel *et al*. 2005), problems that are clearly translated to in our formalism. In other words, the whole task in our formalism is to define the bases *X*_{n} and *Y*_{n}, from and , in a way such that does not vanish at cut-offs. It is shown in the present work that this is possible for cut-offs with zero wavenumber. Finally, the problem is reduced to an evolution problem for the projection coefficients ** a** and

**of**

*b***and**

*X***on the renormalized Lamb modes.**

*Y*### (a) Position of the problem

We are interested in the propagation of Lamb wave through a two-dimensional waveguide with a varying height, described by the function *h*(*x*) (see figure 1), with free boundaries, and for which displacements are in the (*x*, *y*) plane (in-plane motion). For the sake of clarity, the waveguide is considered to be symmetric with respect to the horizontal axis, but the method can be easily extended to non-symmetric geometry .

The time dependence is e^{−iωτ} and will be omitted in the following. The equation of motion is(2.1)where *ρ* is the density,is the displacement vector andis the corresponding stress tensor, with(2.2)where (*λ*, *μ*) are the Lamé constants. The boundary condition at the faces *y*=±*h*(*x*) are free of traction, corresponding to boundary conditions:(2.3)(2.4)where *h*′(*x*) is d*h*/d*x*.

### (b) Modal expansion

It will be convenient to work on two quantities ** X** and

**presented below. That formalism allows us to easily tackle the projection on the Lamb modes. The idea is to write the equations as an evolution equation (with respect to the coordinate**

*Y**x*of the waveguide) on

**and**

*X***, which leads to a canonical eigenvalue problem in the transverse direction when transverse modes are sought. This formulation is similar to the one presented recently in Folguera & Harris (1999), in that it describes the evolution of a stress–displacement four-vector, but here that four-vector is suitably split into**

*Y*It is shown in appendix A that the elasticity equations (2.1)–(2.4) can be written as (see also Pagneux & Maurel 2004)(2.5)where and are the operator matrices:(2.6)with and . Concerning the boundary conditions (2.4), since the problem is considered as an evolution equation (2.5) on ** X** and

**, we always use the expressions of**

*Y**r*written as a linear function of

**, , which implies that the boundary conditions are entirely expressed in terms of**

*Y***and**

*X***.**

*Y*Then, in order to project the evolution equation (2.5), we introduce the modes solutions of the eigenvalue problem associated to (2.5),(2.7)with boundary conditions and at *y*=±*h*. Equation (2.7) is a canonical eigenvalue problem, the solutions of which are the Lamb modes (see Viktorov 1967; Achenbach 1987). Assuming the completeness of the set of the Lamb modes (Kirrmann 1995; Besserer & Malishewsky 2004), the vector (** X**,

**)**

*Y*^{T}is expended as(2.8)Here

*n*>0 refers to right-going mode and

*n*<0 refers to left-going modes. The eigenvalues

*k*

_{n}for right-going modes are sorted in ascending order of their imaginary part and descending order of their real part, and if

*k*

_{n}corresponds to a right-going mode, −

*k*

_{n}corresponds to a left-going mode. As can be seen from equation (2.7), the symmetry properties of these bases impose , , and in the sequel we arbitrarily choose , , as in Viktorov (1967). Owing to these symmetry properties, it is possible to write the decomposition (2.8) in term of the right-going modes only,(2.9)with and .

For the projections, we use the inner product between two component functions that is defined byNote that this inner product is not positive definite, since the vectors ** X** and

**we are going to use are complex.**

*Y*One advantage of our formalism (equation (2.5)) is that the operators and have the following remarkable properties (see appendix A):(2.10)That means that and are formally self-adjoint. In the particular case of the modes which have *r*_{n}(±*h*)=0 and *t*_{n}(±*h*)=0 and thus make the boundary terms vanish, the operators and are self-adjoint. Using the property of and (see (A 5)), we get , that is . This is an easy way to derive the biorthogonality relation obtained by Fraser (1976): .

As could be anticipated, it appears that vanishes for each mode cut-off that corresponds to coalescence of modes (Kirrmann 1995). To avoid this latter problem for zero cut-offs (*k*_{n}=0), we choose to work with renormalized bases *X*_{n} for ** X** and

*Y*_{n}for

**built from and in a way that is detailed in appendix B and where (**

*Y*

*X*_{n},

*Z*

_{c,n}

*Y*_{n})

^{T}is proportional to . Here,

*Z*

_{c,n}appear as renormalization coefficients and correspond to the diagonal entries of the characteristic impedance matrix

_{c}(see §3). These new bases satisfy(2.11)The modal decomposition is now done as(2.12)(the expression of

*J*

_{n}is given in appendix C) and(2.13)One important point is that

*X*_{n}and

*Y*_{n}remain bases for cut-offs with zero wavenumber (

*k*

_{n}=0).

### (c) Evolution equation

The task is now to derive an evolution equation for the modal components ** a**(

*x*) and

**(**

*b**x*). To do that, we project the system (2.5) on the bases

*Z*

_{c,n}

*Y*_{n}and

*X*_{n}:(2.14)Each term is then calculated using (2.12) and properties (A 4):Similarly,This leads to a system of first-order differential equations governing

**and**

*a***:(2.15)where matrices**

*b*_{1}to

_{4}are given by(2.16)Expressions of matrices

_{1}and

_{4}, easier for numerical calculation, are given in appendix D.

## 3. The impedance matrix

In principle, equations (2.15) might be directly integrated, but there are two reasons that prevent one from doing that. First, the problem will be posed as a boundary-value problem with a radiation condition at the outlet of the waveguide and a source at the inlet, and the integration of (2.15) as an initial value problem would be very awkward. Second, the numerical integration of (2.15) appears to be unstable because of the presence of the evanescent modes that induce exponential growth (Pagneux *et al*. 1996). In order to circumvent these two problems, it has been proven to be useful to define the impedance matrix (Pagneux *et al*. 1996) as the linear operator that links together vectors ** a**(

*x*) and

**(**

*b**x*) at a given

*x*position,(3.1)Using the definition of the impedance matrix (

*x*), it is easily obtained from (2.15) that (

*x*) obeys a Riccati matrix differential equation,(3.2)This allows us to solve the problem of guided wave in any configuration by simply solving the Riccati differential equation on , starting from a radiation condition at the outlet, e.g. at the guide terminations, with either

*a*^{+}=0 or

*a*^{−}=0, we get =±

_{c}, where

_{c}is the so-called characteristic impedance matrix which is diagonal with entries

*Z*

_{c,n}.

Once has been calculated, the whole fields can be obtained by integrating(3.3)obtained from equation (2.15) and where the source imposes *a*^{+} at the inlet. The successive integrations of equations (3.2) and (3.3) are numerically stable, even in the presence of evanescent modes.

### (a) Reflection and transmission matrices

If the interest is in the determination of the elastic fields between the inlet and the outlet, one has to integrate successively the above equations (3.2) and (3.3), but this necessitates storing the matrix at each point *x* during the first integration (3.2) to perform the second integration (3.3). On the other hand, if the interest is only in the scattering matrix, i.e. the knowledge of the reflection matrix and the transmission matrix , this storage is not necessary, as shown below.

To determine , which is defined at the inlet abscissa *x*_{ini} by(3.4)it is sufficient to know (*x*_{ini}), since is related to (*x*_{ini}) through(3.5)easily obtained from equations (2.13).

To determine , we define the matrix as the operator that links together the coefficients ** a** between the outlet abscissa

*x*

_{f}and any

*x*≤

*x*

_{f}through(3.6)Note that corresponds to the propagator of equation (3.3). Then, by differentiating equation (3.6) with respect to

*x*, we get the following differential equation governing :(3.7)that has to be integrated from

*x*=

*x*

_{f}to

*x*=

*x*

_{ini}with the initial condition . The transmission matrix is defined by(3.8)and is related to the matrix through(3.9)

Eventually, to obtain matrices and , it is sufficient to integrate the coupled equations (3.2) and (3.7) from *x*=*x*_{f} to *x*=*x*_{ini} and this without storing any matrix.

## 4. Energy flux conservation

In this section, we want to verify that the coupled mode equations conserve the energy flux. At each *x* position, the energy flux *W* is defined by , where the average is taken over time and where is the Poynting vector (overdot indicates the time derivative). In the harmonic regime, we get(4.1)where the overbar means complex conjugation. Introducing the matrix , the energy flux can be expressed as . It is sufficient to use the properties: gives and *k*_{n}→−*k*_{n} gives , , to build the matrix , such as and : for real and imaginary *k*_{n} values , we get *X*_{n} and *Y*_{n} real, so locally equals identity; for *k*_{n} complex, say *n*=2, we consider , so , , and we get locally equalsWith , we thus get = ( equals its transpose, with ^{2}=) and the final expression of the flux,(4.2)

### (a) Energy conservation

Let us calculate *W*′=d*W*/d*x*,(4.3)using . We now write the system (2.15) as(4.4)where (*a*) denotes the diagonal matrix with *a*_{n} as *n*th diagonal elements and , . By differentiating the biorthogonality relation in (2.12), it is easy to check that(4.5)Also, , so and, as well(4.6)

We get, using equations (4.4) and (4.5),(4.7)The matrix equal zero because of property (4.6) with ^{2}=. The two first quantities of the right-hand side term also equal zero: for *k*_{n} real or imaginary, both i*k*_{n}*J*_{n}/*Z*_{c,n} and i*k*_{n}*J*_{n}*Z*_{c,n} are real; for *k*_{n} complex, say *n*=2, it is sufficient to consider (and , ) to check that is equal to(cc denotes the complex conjugate), with zero imaginary part.

We finally obtain *W*′=0, which shows that the system is written in a form that implies energy conservation.

### (b) Reflection and transmission coefficients

To define the fraction of reflected/transmitted energy, the energy flux is written as(4.8)The structure of matrix is as follows: for *k*_{n} real, it locally equals . For imaginary *k*_{n} value, it locally equals . For two complex conjugate, say , we get locally,Consequently, both terms of the right-hand side term in equation (4.8) can be written asThe first term accounts for the energy flux carried by the propagating modes, while the second terms account for the energy flux carried by evanescent modes.

Let us consider the case where *N* propagating modes are sent at the waveguide inlet *x*=*x*_{i}, that is and a radiation condition at the waveguide outlet *x*=*x*_{o}, that is *a*^{−}(*x*_{o})=**0**. We get(4.9)where *N*′ denotes the number of propagating modes at *x*=*x*_{o} (*N*′≠*N* *a priori* if *h*(*x*_{o}) differs from *h*(*x*_{i})). The energy conservation implies *W*(*x*_{i})=*W*(*x*_{o}) and we define the coefficients of reflected energy *F*_{R} and transmitted energy *F*_{T}=1−*F*_{R},(4.10)

## 5. Analytical solutions for slowly varying guides

We present here the WKB (Wentzel–Kramer–Brillouin) approximation of equation (2.5) by using the biorthogonality relation (2.12). This approximation is valid for a slowly varying waveguide. The height *h* is supposed to vary slowly and we denote *ζ*=*ϵx* the slow variable, where *ϵ* measures the slowness. The following WKB ansatz is proposed for the four-vector (as in Folguera & Harris 1999),(5.1)In the following, we denoteand . We insert (5.1) in equation (2.5). At zero order in *ϵ*, we obtain(5.2)whereand with the boundary conditions *r*^{0}=*t*^{0}=0. We thus deduce that (*X*^{0}, *Y*^{0})^{T} and *ϕ*′ are, respectively, an eigenvector and an eigenvalue of the eigenvalue problem (2.7) for a given mode number *n*,(5.3)and(5.4)where the proportionality factor *α*_{n} will be determined by the equation at first order. First order in *ϵ* leads to(5.5)with boundary conditions *t*^{1}=*h*′*s*^{0} and *r*^{1}=0. Equation (5.5) and the associated boundary conditions determine the evolution law for *α*_{n}. Equation (5.5) is projected on the four-vector

(5.6)

Properties (A 4) give and . We thus obtainand using equation (5.4)(5.7)With , we get(5.8)Assuming *α*_{n} non-zero, equation (5.8) is equivalent to(5.9)Eventually, the WKB solution is(5.10)with *φ* and *α*_{n} determined by equations (5.3) and (5.9). Obviously, this WKB solution shows no coupling between modes and it conserves the energy.

## 6. Green tensor

To get the usual Green tensor *G*_{ij} (*i*=1, 2), giving the displacements from a delta function force, we start from(6.1)where ** r**=(

*x*,

*y*) is the vector position. If

**=(**

*u**u*

_{1}=

*u*,

*u*

_{2}=

*v*) denotes the displacement fields satisfying(6.2)where is a delta function force located in

*r*_{0}, we have(6.3)As a consequence, the four components of the Green tensor can be simply deduced through the displacement fields(6.4)

In our formalism, equation (6.2) becomes, from equation (2.5),(6.5)Integrating equation (6.5) w.r.t. *x* between *x*_{0}−*ϵ* and *x*_{0}+*ϵ*, it is easy to see that , and (where ).

We then project on *Y*_{n}, withand on *X*_{n}, withto get(6.6)

For the sake of simplicity and without loss of generality, we choose *x*_{0}=0 in the following. The numerical resolution can be now solved using *Z*(*L*^{+})=*Z*_{c} at the waveguide exit and calculating *Z*(*x*) from *L*^{+} to *x*=0^{+}, as described in §8*b*. Similarly, *Z*(*x*<0) is calculated integrating from *Z*(−*L*^{−})=−*Z*_{c} to *x*=0^{−}. The source is taken into account at *x*=0 with(6.7)where and are given in equations (6.6). From these initial conditions on ** a** and

**at**

*b**x*=0

^{±}, we can calculate and until (and thus the corresponding displacement fields).

### (a) Green tensor of a straight waveguide

In the case of a straight waveguide (*h* constant), the calculation can be done explicitly, since we simply have and , for *x*>0, and , for *x*<0. With , we get(6.8)Using equations (6.4) and (6.6)–(6.8), we obtain(6.9)These expressions are identical to those found by Karhitonov (1978). Note also that at *x*=0, the quantities and defining and vanish, as does any series with terms or for a given function *f*. Here, this means, for instance, that a force applied along the *y*-direction does not produce any displacement along the *x*-direction (this could also be deduced simply by symmetry argument, for instance the symmetry for ).

## 7. Numerical resolution

To solve a typical problem in a waveguide, namely with a radiation condition and a source, one has to solve first the Riccati equation (3.2) and, second, equation (3.3). By integrating equation (3.3), and are known in the whole space and, thus, the stress and displacement fields also.

We propose two numerical resolutions of equations (3.2) and (3.3). One resolution method uses a Magnus method for both equations and is detailed below. This method has three advantages: (i) it gives an exact solution for straight waveguide, (ii) the step size is not imposed by the wavelength, but rather by the typical variation length of the waveguide, (iii) it is not sensitive to the quasi-resonances that may be displayed by the behaviour of the impedance (Schiff & Shnider 1999). However, this Magnus method is not adapted to pass through cut-offs with non-zero wavenumber.

The other resolution method is used when cut-off has to be taken into account. In that case, we add a small dissipation to transform the singularity into quasi-singularity at cut-off. Nevertheless, passing through this quasi-singularity requires a smaller step size and is more time-consuming. This resolution, that uses two classical integration schemes, is not detailed below. To integrate the equation (3.2), we use a classical Runge–Kutta scheme with adaptative step size. The details of the scheme are not developed here and can be found in Press *et al*. (1993) for instance. Then, is then simply calculated solving equation (3.3) using a classical Crank–Nicholson scheme, well adapted to preserve the energy conservation.

### (a) Magnus method

Our scheme is inspired by the techniques proposed by Schiff & Shnider (1999) and Iserles *et al*. (1999). The radiation condition gives at and the source is imposed at (figure 2). Then the interval is discretized with d*x* step, so and a second set is defined.

If we start from equation (2.15), the Magnus method gives(7.1)with matrixand matrices _{1} to _{4} defined by(calculated at midpoint between and ), and that eventually permits one to obtain the following scheme for :(7.2)with . Note that integration is performed from right to left.

Then, equation (7.1) is used once again to get(7.3)where the calculation is done from left to right, starting from .

### (b) Calculation in the inlet/outlet portion with constant cross-section

When the waveguide begin or ends with a portion of constant cross-section, the displacement field can be determined analytically. Suppose that the portion with varying cross-section corresponds to (figure 2). In this portion, the fields are numerically calculated, using either a Runge–Kutta scheme or matrix exponential. To obtain the field between and (waveguide inlet) and the field between and *x*_{f} (waveguide outlet), we use:

For the waveguide outlet, , for both symmetric and antisymmetric modes, with no left-going modes. is known from the numerical calculation between and .

For the waveguide inlet, , we have , with that accounts for the incident wave at . Again, is known from the numerical calculation between and .

## 8. Results

We report in this section results obtained with our method. The spectrum for Lamb modes is determined using the spectral method described in Pagneux & Maurel (2001), with a relative accuracy of 10^{−9}. The material constituting the waveguide has the following properties: Poisson ratio , , and *ρ*=1.

### (a) Reflection in a non-uniform guide

We consider here the calculation of the coefficient of energy reflection and transmission. The geometry consists of a waveguide whose cross-section is described by a Gaussian function,(8.1)with , and .

The incident wave at contains only the first antisymmetric mode and the geometry being symmetric with respect to *y*=0, only antisymmetric modes are considered. The range for *ω* is such that for , only is propagating in the whole waveguide. For , the mode is evanescent for and propagating for . Finally, for , both modes and are propagating in the whole waveguide. The calculation is performed using matrix exponential with and 170 steps in the whole range of considered frequencies.

Figure 3 shows the variation of the coefficient of energy reflection as a function of the frequency *ω*. The energy conservation relation is satisfied in the whole range of frequency with an accuracy of around 10^{−5}.

The curve in plain line corresponds to the energy ratio transported by the mode , always propagating. At the cut-off frequency , the mode becomes propagating and, therefore, transports a part of the energy (curve in dotted line). Finally, at frequency , reaches a maximum very close the value 1. This behaviour indicates that this frequency corresponds to a quasi-trapped mode, whose shape is indicated in figure 4*b*.

### (b) Green tensor

We focus here on the Green tensor in a waveguide whose cross-section is described by the Gaussian function, as in §8*a*, with , , *L*=1. The frequency is *ω*=5. For , there are four symmetric and five antisymmetric propagating modes and for , there are six symmetric modes and six antisymmetric propagating modes.

Since there are cut-offs with non-zero wavenumber in this configuration, the numerical calculation is performed using the Runge–Kutta scheme with a relative tolerance of 10^{−6}. The source point being located at , the calculations are divided into two parts, between 0 and and between 0 and , using the initial conditions of equations (6.7), as described in §6

With and , the Runge–Kutta calculation needs *N*=1500 and 750 steps (respectively for antisymmetric and symmetric modes) for the calculation in the right part and *N*=300 and 350 steps for the calculations in the left part.

The displacement fields obtained are shown in figure 5. Note that we observe wiggles at the vertical of the source point, characteristic of the modal decomposition of the Green tensor. This calculation has been performed with an imaginary part *ϵ* of the frequency equal to 10^{−2} in order to avoid the cut-offs with non-zero wavenumber. The value of *ϵ* is small enough not to influence the final result, as shown in figure 6.

## 9. Closing remarks

The method developed in this paper is a multimodal method for waveguide with height variation. The two main advantages of this method are: (i) it avoids singularities at cut-offs with zero wavenumber and (ii) it can be implemented without numerical instability owing to the introduced impedance matrix. The way to avoid singularities at cut-offs with non-zero wavenumber remains an open question.

We think that it is a valuable alternative to purely numerical methods, e.g. the finite-element method or the boundary-element method, because it allows both the brute force numerical calculation and the asymptotic approximation analysis.

## Footnotes

- Received June 24, 2005.
- Accepted November 8, 2005.

- © 2006 The Royal Society