## Abstract

The scattering of an acoustic wave propagating in a non-uniform waveguide is inspected by revisiting improved multimodal methods in which the introduction of additional modes, so-called boundary modes, allows to better satisfy the Neumann boundary conditions at the varying walls. In this paper, we show that the additional modes can be identified as evanescent modes. Although non-physical, these modes are able to tackle the evanescent part of the field omitted by the truncation and are able to restore the right boundary condition at the walls. In the low-frequency regime, the system can be solved analytically, and the solution for an incident plane wave including one or two boundary modes is shown to be an improvement of the usual Webster equation.

## 1. Introduction

Since the pioneering work of Stevenson [1,2], the multimodal propagation method has been widely used to describe propagation in non-uniform waveguides in acoustics [3–7] and in elasticity [8,9]. In the two-dimensional acoustic case, the pressure *p* satisfies the Helmholtz equation:
1.1in a waveguide of height *h*(*x*) (located in *y*∈[0,*h*(*x*)]) with boundary conditions ∂_{y}*p*(*x*,0)=0 and
1.2In a uniform waveguide of constant height *h*, the solution *p*(*x*,*y*) is expanded in a series , using the transverse modes of the waveguide (where *A*_{n} are normalization coefficients). The function *ψ*_{n} satisfies the boundary conditions *ψ*_{n}′(0)=0 and *ψ*_{n}′(*h*)=0 (as *p*). For a non-uniform waveguide, *p*(*x*,*y*) is still expanded as a series using the local transverse modes (*ψ*_{n} depends now on *x*, and we still have ∂_{y}*ψ*_{n}_{|y=h(x)}=0). Although modal approaches have been shown to be efficient when sufficient evanescent modes are taken into account [4], this representation of *p*(*x*,*y*) is incompatible with the boundary condition (1.2) if *h*′(*x*)≠0. Further, for each *x*-value, the convergence of the series is poor, and the derivative of the series does not converge uniformly in the interval [0,*h*(*x*)]. In a series of papers, Athanassoulis & Belibassakis [10,11] propose an improved representation in which the function *q*(*x*,*y*) is introduced
1.3where *ξ* is a function satisfying ∂_{y}*ξ*_{|y=h(x)}=1 and ∂_{y}*ξ*_{|y=0}=0. This new function satisfies ∂_{y}*q*(*x*,0)=0=∂_{y}*q*(*x*,*h*) (as *ψ*_{n}), so that the series , and its derivative with respect to *y* converge uniformly in the interval *y*∈[0,*h*(*x*)]. However, the value of ∂_{y}*p*(*x*,*h*) is, in general, unknown *a priori,* so that it is not possible to write *p* as a function of *q*.

An exception occurs in the case of an impedance boundary condition of the form ∂_{y}*p*(*x*,*h*)=*Y* _{0}*p*(*x*,*h*) as considered in Bi *et al.* [6]. Choosing *ξ*(*h*;*x*)=0, we deduce from (1.3) that *q*(*x*,*h*)=*p*(*x*,*h*), leading to *p*(*x*,*y*)=*q*(*x*,*y*)+*Y* _{0}*q*(*x*,*h*)*ξ*(*y*;*x*) and thus the solution can be written as
In this case, the convergence of the truncated series *p*(*x*,*y*) is the same as the convergence of the series and the boundary condition ∂_{y}*p*(*x*,*h*)−*Y* _{0}*p*(*x*,*h*)=0 is satisfied strictly for any truncation number *N*.

When considering the boundary condition (1.2), ∂_{y}*p*(*x*,*h*) is unknown, and the approach must be modified. Athanassoulis & Belibassakis [10,11] propose to keep the additional transverse mode *ξ* and to introduce an additional unknown modal component *q*_{N+1}(*x*) in the expansion:
1.4where *q*_{N+1}(*x*) is expected *in fine* to behave as ∂_{y}*p*(*x*,*h*). However, this cannot be guaranteed for any finite truncation. Using this decomposition in (1.1) and (1.2), a system of coupled ordinary differential equations is obtained
1.5for *n*=0,…*N*+1. The system is found to remain coupled in the straight part of the waveguide which implies that the classical radiation condition *q*_{n}′=±i*k*_{n}*q*_{n} cannot be applied directly at the inlet/outlet of the scattering region (where *h*′≠0), with the wavenumbers . In Athanassoulis & Belibassakis [10,11] and Hazard & Lunéville [12], the system (1.5) was solved using finite difference methods by imposing *q*_{N+1}=0 at the inlet/outlet and the usual radiation conditions on the *q*_{n≤N}. The convergence of the improved representation has been shown [12].

In this paper, we revisit the coupled mode equation (1.5) in order to derive an improved system adapted to define radiation conditions. To do so, we need to define the modal components:
in an orthogonal basis of transverse functions (*ψ*_{n}). This is possible, defining from (1.4), the supplementary transverse function with *α*_{n}≡(*ξ*,*ψ*_{n}), *n*≤*N*. We obtain
1.6with
1.7The resulting coupled mode equations are a partially decoupled system for the (*p*_{n}) components, *n*=0 to *N*+1, which is the key point of our new formulation:
1.8where *D*, *E* vanish in the straight parts of the waveguide. Now, approximate solutions can be found using successive Born approximations, as in Maurel & Mercier [13]. The leading-order solution is found for an incident wave *p*^{inc}(*x*,*y*) and
is solution of our system (1.8) for *h*′=0, namely . Then, this solution can be used as the first step of an iterative process to obtain the solution *p*^{(1)}
1.9Applications of the first Born approximation are presented in §4, leading to improved approximate equations compared with the usual Webster equation.

The main advantage of system (1.8) is that, in the straight part of the waveguide, we obtain , for *n*≤*N*, as expected. More surprisingly, the equation for *n*=*N*+1 is
and introduces a new wavenumber *K* (we prefer the notation *K* rather than *k*_{N+1} to avoid confusion with the usual wavenumber). Although this wavenumber is non-physical, it turns out that *K*^{2}<0, indicating that the boundary mode can be interpreted as an evanescent mode. Instead of imposing vanishing value of *p*_{N+1} at the inlet/outlet of the scattering region, as in Athanassoulis & Belibassakis [10,11] and Hazard & Lunéville [12], we can associate to the boundary mode a radiation condition *p*′_{N+1}=±i*Kp*_{N+1}. This makes possible the implementation of efficient numerical multimodal methods, as proposed in Pagneux *et al.* [4], Felix & Pagneux [5] and Pagneux & Maurel [8].

The paper is organized as follows. Section 2 presents in detail the derivation of the improved modal system written using the modal components (*p*_{n}) for the case of a waveguide with one varying wall. The advantages of this formulation are discussed when compared with the classical formulation (without boundary mode). It is shown that, thanks to the boundary mode, the system can be practically truncated at the number of propagating modes. The boundary condition on the varying wall is also inspected as a function of the truncation number, leading to the conclusion that the desired behaviour of the truncated solution *p*^{N}(*x*,*y*), namely ∂_{y}*p*^{N}(*x*,*h*)=*h*′(*x*)∂_{x}*p*^{N}(*x*,*h*), is reached asymptotically, for . The case of two varying walls, where two boundary modes are considered, is presented in §3. Although more involved, the conclusions are the same in this case. Finally, for low frequencies, the systems can be solved in the first Born approximation. The contribution of the boundary mode(s) to the plane wave is shown to improve the Webster equation. This is considered in the §4 and illustrated with examples. Technical calculations and discussion on the derivation of the modal systems are collected in the appendices.

## 2. The case of one varying wall

### (a) Derivation of the improved coupled representation

The initial geometry of the waveguide is set in the (*x*,*y*) plane with *y*∈[0,*h*(*x*)] and . Let us introduce a change of variables to obtain a problem set in a straight guide of unitary height. Considering the transformation *X*=*x*,*Y* =*y*/*h*(*x*) and with *p*(*x*,*y*) satisfying the Helmholtz equation (1.1), the field *p*(*X*,*Y*) satisfies
where *J* is the Jacobian of the transformation (*x*,*y*)→(*X*,*Y*):
We deduce the modified Helmholtz equation:
2.1for , *Y* ∈]0,1[. The boundary conditions ∂_{y}*p*(*x*,0)=0 and (1.2) become
2.2Equation (2.1) is projected onto a basis *φ*_{n}(*Y*), and we obtain
2.3where we have used the boundary conditions equation (2.2). The system of coupled mode equations resulting from this projection depends on the expansion chosen for *p*. First, we will derive the system of equations using the usual expansion without boundary mode. Then, we will show that using the expansion as used in Athanassoulis & Belibassakis [10] and in Hazard & Lunéville [12] results in a system of equations that remain coupled in the straight parts of the waveguide (where *h*′(*X*)=0). As already mentioned, this is not suitable to use the Born approximation or to define radiation conditions. Finally, we use our new formulation based on a similar expansion, leading to a system that decouples the modes in the straight parts of the waveguide.

#### (i) Classical representation

In the classical modal approach (without additional boundary mode):
where *φ*_{n}(*Y*) are the transverse modes of a straight unitary waveguide:
Because of the orthogonality of the *φ*_{n}-functions, (*φ*_{n},*φ*_{m})=*δ*_{mn}, the functions *p*_{n} are the modal components . We obtain, from equation (2.3) and for 0≤*n*≤*N*:
2.4where we have used (*φ*′_{n},*φ*_{m}′)=(*nπ*)^{2}*δ*_{nm} to obtain *k*_{n}(*x*)^{2}≡*k*^{2}−(*nπ*/*h*(*x*))^{2} and where we have defined:
2.5

The above system differs from the derivation proposed in Pagneux *et al.* [4] (this derivation corrects an error in the original derivation of Stevenson 1951 [2]). This is discussed in appendix A. As expected in the straight parts of the waveguide, we recover the usual one-dimensional Helmholtz equations .

#### (ii) Coupled mode equations for the *q*_{n}-components

We consider now the expansion of *p* as carried out in Athanassoulis & Belibassakis [10]:
and we choose *χ* such as *χ*′(1)≠0. For instance, a convenient choice is
Note that (*φ*_{0},…*φ*_{N}) and *χ* are linearly independent as soon as *N* is finite. The system of coupled mode equations can be found in Hazard & Lunéville [12], see eqn. (4.11). Here, we just inspect the form of this system in the straight parts of the waveguide. We obtain, for *h*′(*X*)=0:
2.6with
Surprisingly, the *q*_{n} components, *n*≤*N*, appear to be coupled with the *q*_{N+1} component, although we expect *q*_{N+1} to be useless in the straight parts of the guide. In Athanassoulis & Belibassakis [10], Hazard & Lunéville [12], *q*_{N+1}=0 is imposed outside a bounded calculation domain. Then, (2.6) implies for *n*≤*N*, and the usual radiation conditions *q*_{n}′=±i*k*_{n}*q*_{n} (0≤*n*≤*N*) for can be applied. Obviously, this means that the solutions (*q*_{n})_{n≤N+1} depend on the size of the bounded calculation domain. In the following section, we show that it is possible to decouple the mode equations outside the scattering region independently of the size of the calculation domain. In the process, the mode *q*_{N+1} is identified as an evanescent mode, associated with a new wavenumber *K*. Our coupled mode equation does not restrict to a bounded calculation domain, because a natural radiation condition *p*_{N+1}=±i*Kp*_{N+1} is obtained.

#### (iii) New coupled mode equations on the *p*_{n}-components

To prevent the modal components from remaining coupled in the straight parts of the waveguide, we use a reformulation of the expansion in terms of new unknowns *p*_{n} such that
2.7with
2.8The *p*_{n} and the *q*_{n} components are linked by the relations:

We obtain, for 0≤*n*,*m*≤*N*+1,
2.9where *a*_{nm} and *d*_{nm} are defined in equation (2.5) (and applied here for 0≤*n*,*m*≤*N*+1) and with
2.10The above properties are important. The functions (*φ*_{n})_{n≤N+1} are orthogonal, by construction and the functions, (*φ*′_{n})_{n≤N+1} are also orthogonal. It follows that our reformulation of the modal expansion succeeds in decoupling the modal components in the straight parts of the waveguide: for the *N*+1 first modes, 0≤*n*≤*N*, we recover the expected propagation equations , as in the classical projection. For *n*=*N*+1, we obtain
2.11We have (*β*_{n})_{n≤N}=1, (*γ*_{n})_{n≤N}=(*nπ*)^{2} and
2.12This shows that the additional mode is associated with a wavenumber *K* such that
2.13Let us prove now that *p*_{N+1} is an evanescent mode. For each *x*-value, we note *n*_{p}(*x*) the number of propagating modes (*n*_{p} is the integer part of *kh*(*x*)/*π* plus one). If we introduce where the sup is taken on all the *x*-values, then the mode *φ*_{Np} is the first mode evanescent in the whole waveguide. Assuming that the truncation does not eliminate the propagating modes (*N*≥*N*_{p}−1), we have *K*^{2}<0 which corresponds to an evanescent mode. Indeed, using *nπα*_{n}≥(*N*+1)*πα*_{n} for all *n*≥*N*+1, we obtain
and for all *x*-values by definition of *N*_{p}.

We return to our system of coupled mode equations. In the part of the waveguide with varying cross section, we have 2.14

A reasonable question is what happens for . The equations for *p*_{n}, *n*≤*N*+1 involve (*a*_{n,N+1}, *a*_{N+1,n}) and (*d*_{N+1,n},*d*_{n,N+1}). We have for *n*≤*N*
2.15and
2.16The expressions of the coefficients *a*_{nm} and *d*_{nm} are given in the appendix B. When , it is easy to check the following behaviours:
2.17(an example of derivation of this equivalent is given in appendix A). Let us consider the equation on *p*_{n}, *n*≤*N* in equation (2.14). It is easy to check that the equation tends to the classical equation (2.4) when (without boundary mode). For instance, *a*_{nm}∼*m*^{2}/(*n*^{2}−*m*^{2}), so that *a*_{nm}→−1 for large *m*, while *a*_{n,N+1}∼1/*N*→0. This is expected, because the additional degree of freedom *p*_{N+1} becomes unnecessary. The equation on *p*_{N+1} has the same structure, but with an non-physical wavenumber *K* that depends on the truncation. Asymptotically, equation (2.14) for *n*=*N*+1 leads to constant *p*_{N+1} (the dominant term is (*h*′/*h*)*δ*_{m,N+1}*p*′_{m} in the right-hand side term, which is *O*(*N*^{3})) and uncoupled to the other modes *p*_{m}.

### (b) Convergence and boundary condition in the improved representation

The convergence and the errors of the improved method compared with the classical method are illustrated in figures 1 and 2. In the calculation, the non-uniform part of the waveguide is given by for *x*∈[−0.5,0.5] (geometry A). A plane wave is sent from the left at a frequency *kh*=2, for which *N*_{p}=2, two propagating modes exist in the largest part of the waveguide. The presented results have been obtained using a numerical scheme based on the use of the admittance matrix, as described in Pagneux *et al.* [4], (A. Maurel, J.-F. Mercier & V. Pagneux 2013, unpublished data), implementing the evanescent boundary mode such as other evanescent modes.

Figure 1 shows the real part of the acoustic field, in the improved method for a truncation just at the two propagating modes (*N*=1, the field is described by *N*_{d.f.}=3 degrees of freedom (d.f.), the two propagating modes and the boundary mode) and in the classical method for truncations including the two propagating modes plus 1–41 additional evanescent modes. It can be seen that the upper wall boundary condition is reasonably verified in the improved method, whereas the classical method suffers from oscillations, and it needs about 40 evanescent modes to behave satisfactorily at the walls. More quantitatively, figure 2*a* shows the errors on the total field as a function of the truncation *N*. The reference field has been calculated considering *N*=100, for which the fields obtained with or without boundary modes are equal up to less than 1 per cent. The error is defined in *L*^{2} norm: ∥*p*^{N}−*p*^{100}∥_{L2}/∥*p*^{100}∥_{L2}. The general trends of the errors are in agreement with the prediction of Hazard & Lunéville [12], with an error varying such as *N*^{−3/2} in the classical method and *N*^{−7/2} in the improved method (the dependence of the ∥*q*_{n}∥_{L2} and ∥*p*_{n}∥_{L2} modal components is also reported in figure 2*b* as a function of *n*; ∥*q*_{n}∥_{L2}∝*n*^{−4} in the improved method, whereas ∥*p*_{n}∥_{L2}∝*n*^{−2} in the classical modal method). More remarkable is the case where *N*_{d.f.}=3 d.f. (three terms taken into account in the modal expansion) are considered, already seen qualitatively in figure 1. Adding to the two propagating modes, the third d.f. is the first evanescent mode in the classical method and it is the boundary mode in the improved method. The resulting error on the total field is about 30 per cent in the classical method, whereas it is only 5 per cent in the improved method, and one has to consider more than 10 evanescent modes in the classical method to reach the same small error (note that, although the error on the total field is small for *N* around 10, one can see in figure 1 that the field still suffers from oscillations).

We now focus on the boundary condition on the non-uniform wall at *y*=*h*(*x*). The mode *p*_{N+1} has been shown to be an evanescent mode. This mode is excited in the perturbed region *h*′≠0 but is not equal to zero outside this region, it is only exponentially decreasing in the straight part of the waveguide. This seems to be in contradiction with the desired condition ∂_{y}*p*(*x*,*h*)=*h*′∂_{x}*p*(*x*,*h*)=0 in the straight part. To clarify this point, let us consider the truncated condition as a function of the truncation number. The change of variable implies that
2.18and we will prove now that this quantity is of order *O*(1/*N*). In this aim, the equation for *p*_{N+1} in (2.14) is rearranged in the following way:
2.19Owing to the behaviour of the coefficients in equations (2.12) and (2.17), it appears that the left-hand side in (2.19) is of order 1/*N* (let us recall that *p*_{n} decreases with *n*, |*p*_{n}|∼1/*n*^{2}), whereas the right-hand side is of order 1/*N*^{2}. More precisely, we obtain at dominant order:
Using *φ*_{0}(1)=1, , 0<*n*≤*N* and , the above equation can be written as
It is also easy to see, from *p*_{n}=*q*_{n}+*α*_{n}*p*_{N+1}, that
(using ). It follows that the truncated boundary condition (2.18) satisfies
and the right boundary condition is verified for .

## 3. The case of two varying walls

In the case of two varying walls of shapes *h*_{1}(*x*) and *h*_{2}(*x*), two extra degrees of freedom *p*_{N+1}(*x*) and *p*_{N+2}(*x*) are considered in the expansion of *p*(*x*,*y*), and it is expected that each degree of freedom will tackle one of the two boundary conditions on the non-uniform walls.

The procedure to derive a family of one-dimensional wave equations is similar to the case of one varying wall. The transformation (*x*,*y*)→(*X*,*Y*) is now defined by *x*=*X* and *y*=*h*_{1}+*Y* *h*, where *h*=*h*_{2}−*h*_{1}. The Jacobian of this transformation is
with *f*(*Y*)≡*h*_{1}′+*Y* *h*′, from which we deduce the modified Helmholtz equation:
with boundary conditions:
3.1where *p*(*X*,*Y*) is expanded onto the usual basis (*φ*_{n})_{n≥0,} and on two additional modes *φ*_{N+1} and *φ*_{N+2}:
and we choose *φ*_{N+1} and *φ*_{N+2} as:
with
and
3.2The normalizations are and . The projection of the wave equation (3.1) is
for 0≤*n*,*m*≤*N*+2, with *a*_{nm} and *d*_{nm} defined in equation (2.5) and

The key point is that we have chosen *χ*_{1} and *χ*_{2} such that (*φ*_{N+1},*φ*_{N+2})=(*φ*′_{N+1},*φ*′_{N+2})=0 (otherwise, for *n*≤*N*, (*φ*_{N+1},*φ*_{n})=(*φ*_{N+2},*φ*_{n})=0 and (*φ*_{N+1}′,*φ*_{n}′)=(*φ*_{N+2}′,*φ*_{n}′)=0 by construction). This is because we have chosen (*χ*_{1},*χ*_{2})=(*χ*_{1}′,*χ*_{2}′)=0 and such that for any *n*-value (see equation (B8)), from which we deduce
3.3Again, the equations on *p*_{n} are decoupled in the straight parts of the waveguide, with (and ) for the usual modes *n*≤*N* and
3.4and
3.5In the part of the waveguide with varying cross section, we have, for *n*=0,…(*N*+2)
3.6The coefficients (*a*_{nm},*b*_{nm},*c*_{nm},*d*_{nm}) and (*γ*_{n},*β*_{n}) are given in the appendix B. As for one varying wall, we find that the modal components *p*_{n≤N} are the usual modal components associated with the wavenumbers , whereas the two boundary modes appear to be evanescent with purely imaginary wavenumbers , *j*=1,2.

Typical results are shown in figure 3 for two varying walls with for *x*∈[0,1] and *h*_{2}(*x*)=*h*_{1}(*x*)+1. We consider *α*=0.25*π* (geometry B) and *α*=0.37*π* (geometry C). At a frequency *kh*=2, only the plane mode is propagating. It can be seen that the boundary conditions on each walls are reasonably verified as soon as *N*=2 in the improved method (*N*_{d.f.}=4), whereas the classical method suffers from oscillations after truncation at a few modes and it needs more than 10 evanescent modes to approach reasonably the boundary conditions. The conclusions (not reported) on the convergence versus *N* and the boundary conditions are the same as in the case of one varying wall.

## 4. Low-frequency limit: revisiting Webster equation

In this section, the low-frequency regime is inspected and the possibility to improve the usual Webster equation, namely
4.1where *p*_{0}(*x*)=e^{ikx} at leading order is investigated. At a given frequency *kh* (*h* is the height outside the perturbed area), the Webster equation inspects the first-order correction in *O*(*h*′) owing to the inhomogeneity in the cross section. A derivation of the Webster equation can be found in Rienstra [14] (see also a brief discussion in the appendix C).

We will show that (2.14) for one boundary mode and (3.6) for two boundary modes degenerate in coupled equations of the general form
4.2where refer to the wavenumber and modal component of one boundary mode, *F*,*G* are constants that will be given later. The coupled system involves now so that the plane mode *p*_{0} can be determined up to *O*(*h*^{′2}). In the following, we refer to the improved Webster approximation as IWA (note that a different improvement is proposed for axisymmetric three-dimensional geometry in reference [15]).

Although Webster’s equation (4.1) and the improved representation (4.2) can be solved numerically, we consider here analytical solutions by using the first Born approximation (namely *p*_{0}(*x*)=e^{ikx} at leading order). We inspect two cases leading to simplifications. In the first case, *h*(*x*) varies on a typical scale *L* much larger than *h* (slowly varying section), and in the second case, *L* is much smaller than any scale of the problem (sudden variation).

### (a) General expressions

In this section, we propose an analytical solution of the IWA equation (4.2), considering the case of one boundary mode. The case of two boundary modes follows simply by linearity. The second equation of (4.2) is solved in the first Born approximation (*p*_{0}(*x*)≃e^{ikx}):
4.3where *g*_{K}(*x*)≡e^{iK|x|}/(2i*K*) is the Green function for the one-dimensional wave equation. Reporting the solution for in the first equation of (4.2) leads to:
4.4where refers to the solution of the Webster equation (4.1) and where *g*′_{k}(*x*−*y*) denotes the derivative with respect to *x*. Note that we have used *h*(*x*)≃*h* (remember that |*h*′|≪1) in the second equation.

#### (i) Slowly varying waveguide

We consider a cross section varying of *Δh*≪*h* over a length *L* with *h*/*L*≪1. As previously said, *K* is imaginary, and here, *K*^{2}=*k*^{2}−*γ*_{N+1}/*β*_{N+1}*h*^{2}≃−*σ*^{2}/*h*^{2}, with *σ*^{2}≡*γ*_{N+1}/*β*_{N+1}. We use for *f* varying over a typical length *L*≫*l*. Applying this result with *σL*/*h*≫1, we obtain from equation (4.3), at dominant order and in the perturbed area (where *h*′≠0):
4.5Equation (4.4) simply becomes
4.6where *C*≡*FG*/(2*σ*^{2}).

#### (ii) Sudden variation in the section

We consider now a sudden variation located at *x*=*x*_{0}, in which *h* varies of *Δh* on a typical scale *L* much smaller than any scale of the problem. Then, *h*′(*x*)=*Δhδ*(*x*−*x*_{0}) and equation (4.3) reduces to:
4.7Then, the solution of the equation (4.4) simply follows (note that here the solution takes a very simple form, because its second order in (*Δh*/*h*)^{2} vanishes):
4.8with *D*=*FG*/4*σ*.

### (b) The case of one varying boundary

Here, the case of one boundary mode is treated, and a discussion on the use of two boundary modes is proposed in the following section. With *a*_{00}=*d*_{00}=0, and, for *N*=0, we have from equations (B3)–(B4): *a*_{0+1,0}=*d*_{0,0+1}=0 and (and *γ*_{0+1}=(*π*/2)^{2}, *β*_{0+1}=1−8/*π*^{2}, from equations (B1)–(B2)), the system (2.14) reduces to the following system at dominant order:
which can be identified to the system (4.2) with *F*=*a*_{0,0+1}/*β*_{0+1} and *G*=*a*_{0,0+1}. Then, the expressions (4.6) and (4.8) can be directly applied, taking and . The comparisons of IWA solutions to the Webster approximation (WA) are illustrated in figures 4 and 5 for a slowly varying waveguide with cosine modulation and for a waveguide with the upper boundary varying suddenly. In both cases, a significant improvement is observed although a quite high frequency *kh*=2 is considered. In the case of the sudden variation of the upper boundary (figure 5), the IWA solution is able to capture the discontinuity of the plane mode at the expansion location (from the Webster equation (4.1) written in the conservative form (*hp*_{0}′)′+*k*^{2}*hp*_{0}=0 it is possible to prove that its solution is continuous, even for *h* discontinuous).

### (c) The case of two varying walls

We use *a*_{00}=*b*_{00}=*c*_{00}=*d*_{00}=0, to obtain, from equations (B10)–(B11): for *j*=1,2, , , and ; we also have , , , , the system of equation (3.6) leads, at dominant order, to the simplified system for :
By identifying the above system with equation (4.2), it is straightforward to use equations (4.6) and (4.8). Indeed, tackles the variation of cross section *h* with *F*_{1}≡*ξ*_{1}(1−4/*π*)/*β*_{1} and *G*_{1}=*ξ*_{1}(1−4/*π*), whereas tackles the variation of the centreline of the waveguide (here 2*h*_{1}+*h*=*h*_{1}+*h*_{2} plays the role of *h*′ in (4.2)) with *F*_{2}=−*ξ*_{2}/*β*_{2} and *G*_{2}=−*ξ*_{2}. We also have, from equations (B8) and (B9), *β*_{1}=1−(4*ξ*_{1}/*π*)^{2}, and *β*_{2}=1, .

More precisely, equations (4.5)–(4.6) become, for slowly varying *h* and *h*_{1}:
4.9where and . Note that in the case of one varying wall, *h*_{1}′=0, the solution for *p*_{0}(*x*) in (4.9) is identical to the solution of equation (4.6) using one boundary mode by considering the change *C*→*C*_{1}+*C*_{2} and, as expected, we have *C*_{1}+*C*_{2}∼0.1655, very close to the value *C*∼0.1643 when using one boundary mode.

An illustration of a slowly varying waveguide is given in figure 6 for a constant section guide but varying *h*_{1} (the deviation is sudden but *h*_{1} is continuous). We considered the simplified expression obtained for a portion of waveguide between *a* and *b* with constant height *h* and forming an angle *α* with *x* (leading to *h*′=0, ). Equation (4.9) simplifies in
4.10The result is shown in figure 6 for *α*=0.1 and (*a*=2.5,*b*=7.5) at a frequency *kh*=2. Here, because *h*′=0, the Webster’s equation does not predict any effect of the bending on the plane mode, and the effect captured by the boundary mode in IWA follows a law.

If the variations are sudden for both walls, then we consider *h*_{1}(*x*)=*Δh*_{1}*δ*(*x*−*x*_{0}) in addition to *h*(*x*)=*Δhδ*(*x*−*x*_{0}), which leads to, using equations (4.7) and (4.8):
4.11where and (where we used , *j*=1,2). Note again that if only the upper boundary experiences a sudden change (figure 5 with *Δh*_{1}=0), then the solution for *p*_{0}(*x*) in (4.11) is identical to the solution of equation (4.8) using one boundary mode with *D*→*D*_{1}+*D*_{2}. Although we have here *D*_{1}+*D*_{2}∼0.3648, quite different from *D*∼0.2964 when using one boundary mode, no significant difference is observed in the IWA profile *p*_{0}(*x*) (the result is not reported).

The result for a symmetric sudden expansion at frequency *kh*=2 is shown in figure 7. The agreement between the reference solution and the IWA solution is excellent, better than the agreement found for a sudden change in the upper boundary only (figure 5). This is probably due to the symmetry of the sudden expansion: because odd modes are not allowed, the first evanescent mode that can be excited is the mode 2, associated with the cut-off frequency *k*_{2}*h*=2*π*. This frequency is far from the incident wave frequency *kh*=2, which means that the mode 2 is very evanescent, well captured by the added boundary modes. On the contrary, for the non-symmetric expansion, the mode 1 of cut-off frequency *k*_{1}*h*=*π* seems to be too close to a propagating mode to be satisfactorily tackled by the boundary mode.

## 5. Conclusion

In this paper, we have revisited an efficient method developed earlier which consists of adding an extra non-physical mode to the usual modal expansion to obtain a better convergence of the modal series. By performing a change of unknowns we are able to partially decouple the modal components, which improves the boundary mode method and leads to at least two interesting consequences. (i) It allows to identify the nature of the boundary mode and its relation with the usual modes. This defines radiation conditions and thus facilitates the use of efficient numerical methods such as the admittance matrix method. The numerical tests have shown that our method is very efficient in reducing the number of degrees of freedom: adding to the boundary modes, it is sufficient to take only the propagating modes to obtain very good results. Works are in progress to investigate in detail the strength of such approach in multimodal numerical schemes (A. Maurel, J.-F. Mercier & V. Pagneux 2013, unpublished data). (ii) In the low-frequency regime, which is the main goal of the present paper, the boundary mode is used to derive new approximate equations improving the Webster equation. Extensions to three-dimensional axisymmetric waveguides and to bent waveguides with varying cross section are under progress.

## Acknowledgements

We acknowledge the support from the Agence Nationale de la Recherche, ANR-10-INTB-0914 ProComedia. The authors thank Christophe Hazard and Eric Lunéville for useful discussions.

## Appendix A. Projection of the wave equation

The derivation of the wave equation (2.9) (with *h*_{1}′=0) is here recovered, and compared with classical projection. For the classical projection, we refer to Pagneux’s derivation [4], who avoided an error in Stevenson 1956 that is commented below. The direct projection of the wave equation onto the eigenfunctions *ψ*_{n} is classic Stevenson (1951), Pagneux *et al.* [4], with
and *ψ*_{0}(*y*;*x*)=1, , satisfying *hr*_{n}=(*hp*_{n})′−*h*′*p*_{m}(*f*_{mn}−*a*_{mn}), with *f*_{mn}≡*ψ*_{m}(*h*;*x*)*ψ*_{n}(*h*;*x*). The projection of is, using the boundary condition at *y*=*h*:
The projection of starts with
A1and the difference in the obtained representation comes from the treatment of this derivative.

In the following, we need the following relations:

— Starting from the definition of , we deduce

*y*∂_{y}*ψ*_{n}(*y*;*x*)=*a*_{ln}*ψ*_{l}(*y*;*x*). Then, it follows*a*_{ln}*f*_{lm}=0 for any (*n*,*m*), where we used ∂_{y}*ψ*(*h*;*x*)=0.— We also have .

— Integrating by part , we have

*a*_{mn}+*a*_{nm}=−*δ*_{mn}+*f*_{mn}.

**(a) Classical derivation by direct projection of the wave equation**

In Pagneux *et al.* [4], the second derivative is, from equation (A1),
using the boundary condition for *p* at *y*=*h*(*x*) to obtain (d/d*x*)*p*(*x*,*h*)=(1+*h*^{′2})∂_{x}*p*(*x*,*h*) and our *ψ*_{n}(*h*;*x*) is independent of *x*. Note that it is this derivation that is done in Stevenson by deriving term by term, namely Stevenson considers . Using rather [−2*a*_{mn}+*f*_{mn}(1+*h*^{′2})]*hr*_{m}=[−2*a*_{mn}+*f*_{mn}(1+*h*^{′2})]*p*′_{m}+2(*h*′/*h*)*d*_{mn}*p*_{m}, we obtain the system of coupled differential equations on the *p*_{n}:
A2

**(b) Alternative projection**

Alternatively, we can consider the derivation of equation (A1):
Using *a*_{mn}*hr*_{m}=*a*_{mn}*hp*_{m}′−*d*_{mn}*h*′*p*_{m}, we obtain
that corresponds to our equation (2.4).

## Appendix B. Expressions of (*α*,*β*), (*a*,*b*,*c*,*d*)

For *m*≤*N* and *n*≤*N*, we have obviously *β*_{n}=1 and *γ*_{n}=(*nπ*)^{2}. Then
and

**(a) For one degree of freedom**

We have
B1The coefficients *β*_{N+1} and *γ*_{N+1} are given using (*χ*,*χ*)=1, (*χ*′,*χ*′)=(*π*/2)^{2} and (*φ*_{n},*φ*_{m})=*δ*_{mn}, (*φ*_{n}′,*φ*_{m}′)=(*nπ*)^{2}*δ*_{mn} for (*m*,*n*)≤*N*
B2We also need the coefficients *a*,*b*,*c*,*d* for the boundary mode (although not necessary for one varying wall, we define also the coefficients *b* and *c*, see following section).
B3where *a*_{m}≡(*Y* *φ*_{m},*χ*′), , *b*_{m}≡(*φ*_{m},*χ*′), , *c*_{m}≡(*Y* ^{2}*χ*′,*φ*_{m}′) and *d*_{m}≡(*Y* ^{2}*χ*′,*φ*_{m}′).
B4and
B5The asymptotic forms of *a*_{m,N+1}, *a*_{N+1,m} and *d*_{m,N+1} can be checked numerically, but, because they are rests of series, they can be evaluated explicitly. For instance, we have
B6We can now evaluate and for large *N*. It follows that *d*_{m,N+1}∼8(−1)^{m}*m*^{2}/(3*πN*^{3}).

Finally, we also need , (*χ*,*χ*′)=−1, and to obtain
B7

**(b) For two degrees of freedom**

It is more convenient to express the coefficients as a function of the formers. We have
B8with *α*_{m} given in equation (B1). We also have
B9Then, we define, for *j*=1,2
B10with *c*_{m,N+j} and *d*_{m,N+j} being the symmetrical forms. The coefficients can be expressed as a function of the coefficients in equation (B5).
B11

The coefficients for *n*=*N*+*j* or *m*=*N*+*j*, *j*=1,2 are of the form (we give the example for *a*_{N+1,N+2})
B12and this can be carried out for all coefficients *b*,*c*,*d*. It appears that we need the following integrals , , (*Y* *χ*_{1},*χ*′_{2})=−*ξ*_{1}*ξ*_{2}/(2+*π*), (*Y* *χ*_{2},*χ*′_{1})=−*ξ*_{1}*ξ*_{2}/(2−*π*) (*χ*_{1},*χ*′_{1})=0, (*χ*_{2},*χ*′_{2})=0, (*χ*_{1},*χ*′_{2})=−*πξ*_{2}/(2*ξ*_{1}) (*χ*_{2},*χ*′_{1})=*πξ*_{1}/(2*ξ*_{2}) , , (*Y* *χ*_{1}′,*χ*′_{2})=*ξ*_{1}*ξ*_{2}/2, , , (*Y* ^{2}*χ*_{1}′,*χ*′_{2})=−*ξ*_{1}*ξ*_{2}/2.

## Appendix C. On the Webster equation

The plane wave approximation, where *p*(*x*,*y*)≃*p*_{0}(*x*) is known to produce the Webster equation:
In the low-frequency regime (*ϵ*≡*kh*≪1), we have *p*_{0}=*O*(1) (incident mode). For *n*≠0, from equation (2.4), using we deduce that *p*_{n}=*O*(*h*^{′2},*ϵh*′). For *n*=0, with *a*_{00}=*d*_{00}=0, equation (2.4) leads to:
from which the Webster equation can be deduced if *h*^{′2}≪*ϵ*≪1. A particular case satisfying such condition is used in [14] with *h*′=*ϵ*. Equation (A2) leads to a slightly different equation:
As it was already mentioned in Pagneux *et al.* [4], equation (A2) seems to have an extra term in *h*′^{3}. However, the extra term,
has to be neglected, leading to the usual Webster equation.

- Received April 9, 2013.
- Accepted May 22, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.